I Born Oppenheimer problem

In Born-Oppenheimer approximation there is one step, when you divide your wavefunction into two pieces - first dependent on nuclei coordinates only and second dependent on electron coordinates only (the nuclei coordinates are treated as parameter here). The "global" wavefunction is a product of these two. Why "almost" independent movement of nuclei and electrons determine form of global wavefunction as product of electron and nuclei wavefunctions? What is mathematical proof of that?

If you write the Hamiltonian as ##H_{tot}=T_N+T_e +V(r,R)## with ##T_{N/e}## being the kinetic energy operator of the nuclei and electrons, respectively, and V the Coulomb interaction of the electrons with coordinates r and nuclei with coordinates R, then the electronic hamiltonian is
##H_{el}=T_e+ V(r,R)## with eigenvalues ##\psi_n(r;R)##. These eigenvalues form a complete basis in which also the eigenvalues ##\Psi_m(r,R)## of ##H_{tot}## can be developed, namely ##\Psi_m=\sum_n \psi_n(r;R) \phi_{nm}(R)##.
Born and Oppenheimer now claim that it is - at least sometimes - a good approximation to keep only a single term of the sum. The condition for this to be approximately true is that the action of ## T_N ## on ##\psi_n## can be neglected, which is justified mainly by the dependence of ##T_N## on the small factor ##m/M## where ##M## is the mass of the electron and M a typical mass of the nuclei. A further condition is that the electronic states are energetically well separated - this condition fails for example for Jahn-Teller states, where orbitals become degenerate due to symmetry restrictions.
I invite you to set up an equation for the ##\phi_{nm}## so we can work out the details.