Born-Oppenheimer approximation: calculating the expectation values for the molecular hamiltonian Here's what I've understood so far.
The Hamiltonian for a molecule made of $N$ nuclei, and $n$ electrons is:
$${\cal \hat H} = \underbrace{{-\frac {\hbar^2}{2} \sum _{\alpha =1}^N } \frac {\nabla_\alpha ^2}{M_\alpha}}_{\rm {kinetic \ energy \\ of \ all \ nuclei }}
 \underbrace{{-\frac {\hbar^2}{2m} \sum _{i=1}^n} {\nabla_i ^2}}_{\rm {kinetic \ energy \\ of \ all  \ electrons }} 
 \underbrace{+ {\sum _{i<j}} \frac {e^2}{r_{ij}}}_{\rm {\ electric \ repulsion \\ electron-electron }}
 \underbrace{+ {\sum _{\alpha<\beta}} \frac {Z_{\alpha} Z_{\beta}e^2}{R_{\alpha \beta}}}_{\rm {\ electric \ repulsion \\ nucleus-nucleus }}
 \underbrace{- {\sum _{\alpha, i}} \frac {Z_\alpha e^2}{r_{\alpha i}}}_{\rm {\ electric \ attraction \\ electron-nucleus }}$$
Or, in short terms:
$$ {\cal \hat H}= T_N + T_e + V_{ee} + V_{NN} + V_{Ne}$$
The strategy adopted to find a solution for the time-independent Schrödinger equation ${\cal \hat H} \psi = E \psi $ in the Born-Oppenheimer approximation is the following:

*

*Assume a factorized eigenstate of the kind $\psi ({\bf r},{\bf R}) = \phi_N ({\bf R}) \phi_e({\bf r}, {\bf R})$. (Here ${\bf r} =\{ {\vec r_i} \}_{i = 1, \dots, n}$, ${\bf R} =\{ {\vec R_\alpha} \}_{\alpha = 1, \dots, N} $)

*Treat ${\bf R} $ as a parameter and optimize $\langle  \psi | {\cal \hat H}| \psi \rangle =E $ with respect to $\bf R$. (We are hence using a variational approach)

Let's analyze now the matrix element $\langle  \psi | {\cal \hat H}| \psi \rangle =E $ for each term of the Hamiltonian written above.
My troubles start right away with the first term: $T_N$
$$ {T_N}| \psi \rangle = {T_N}| \phi_N \phi_e \rangle = -\frac {\hbar^2}{2} \sum _{\alpha =1}^N  \frac {\nabla_\alpha ^2}{M_\alpha} \phi_N \phi_e $$
applying vector identity: $\nabla ^2 (ab) = b \nabla^2 a + 2 {{\vec \nabla} a \cdot {\vec \nabla} b + a \nabla ^2 b}$
$${T_N}| \phi_N \phi_e \rangle = \phi _e T_N \phi_N + \phi_N T_N \phi_e - \hbar^2 \sum_\alpha {1 \over M_\alpha}{\vec \nabla_\alpha} \cdot \{ \phi_e \vec \nabla_\alpha \phi_N + \phi_N \vec \nabla_\alpha \phi_e \} $$
My professor said the second term is negligible, and if you apply the bra to the third term it vanishes.
uhm... can you show me?
 A: Question: can the term $\phi_N T_N \phi_e$ be neglected?
$${T_N}| \phi_N \phi_e \rangle = \phi _e T_N \phi_N + \phi_N T_N \phi_e - \hbar^2 \sum_\alpha {1 \over M_\alpha}{\vec \nabla_\alpha} \cdot \{ \phi_e \vec \nabla_\alpha \phi_N + \phi_N \vec \nabla_\alpha \phi_e \} \tag{1} $$

The electron and nucleu combined  system has two time scales differ by serveral order of magnitude: the electron motion and the nucleu motion. We can examine Eq.(1) by taking the integral over the fast oscillatory function, multiply eq.(1) by $\phi_e$ and integral over spatial dimention.
$$ \int d^3 \phi_e^* \left\{ {T_N}| \phi_N \phi_e \rangle =\phi _e T_N \phi_N + \phi_N T_N \phi_e - \hbar^2 \sum_\alpha {1 \over M_\alpha}{\vec \nabla_\alpha} \cdot \{ \phi_e \vec \nabla_\alpha \phi_N + \phi_N \vec \nabla_\alpha \phi_e \} \right\} $$
Using $\int d^3\vec r \phi_e^* \phi_e = 1$:
\begin{align}
 \langle \phi_e| T_N| \phi_N \phi_e \rangle = T_N \phi_N + \phi_N \langle \phi_e| T_N| \phi_e \rangle - \hbar^2 \sum_\alpha  {1 \over M_\alpha} \langle \phi_e | { \vec \nabla_\alpha} \cdot \{ |\phi_e \rangle \vec \nabla_\alpha \phi_N + \phi_N \vec \nabla_\alpha |\phi_e \rangle \} 
\end{align}
The spirit of Born-Oppenheimer treatment is taking nucleu in quasi-static mode when looking at the electron motion, therefore the integrals can be approximated:
$$
\langle \phi_e| T_N| \phi_e \rangle \approx T_N \langle \phi_e| \phi_e \rangle = T_N 1 = 0.
$$
A: We remove the electronic coordinate dependency of the equation by projection onto an adiabatic electronic eigenstate and integration over electronic coordinates. The adiabatic electronic eigenstates are orthogonal and satisfy
$$ 
\langle \phi_e^n|\phi_e^m\rangle _r = \delta_{nm}
$$
With this we can simplify the equation,
$$\begin{aligned}
\langle \phi_e^n|T_N|\phi_N\phi_e^m\rangle_r  &= \langle \phi_e^n|\phi _e^m T_N \phi_N\rangle_r + \langle \phi_e^n|\phi_N T_N \phi_e^m\rangle_r \\
&- \hbar^2 \sum_\alpha {2 \over M_\alpha}{\langle \phi_e^n|(\vec \nabla_\alpha}   \phi^m_e) \cdot  (\vec \nabla_\alpha \phi_N) \rangle_r\\
&=\langle \phi_e^n|\phi _e^m \rangle_rT_N \phi_N + \langle \phi_e^n| T_N \phi_e^m\rangle_r \phi_N\\
&- \hbar^2 \sum_\alpha {2 \over M_\alpha}{\langle \phi_e^n|(\vec \nabla_\alpha}   \phi^m_e)  \rangle_r\cdot  (\vec \nabla_\alpha \phi_N)\\
&=\delta_{nm}T_N\phi_N + d^{(2)}_{nm}\phi_N - \hbar ^2 \sum_\alpha\frac{2}{M_\alpha}\vec d^{(1)}_{nm}\cdot \vec \nabla_{\alpha}\phi_N
\end{aligned}$$
where I have defined the first and second order kinetic coupling elements $d^{(1)}_{nm}, d^{(2)}_{nm}$. These coupling elments couple the nuclear wavefunctions of different electronic states. The adiabatic approximation is typically understood as the neglection of these coupling elements. In the adiabatic approximation, the nuclear wavefunctions of different electronic states are uncoupled and we typically work within this framework when we describe molecules at equilibrium.
The approximation is valid as the derivative of the adiabatic electronic states with respect to nuclear coordinates is typically small at molecular equilibrium coordinates.
But the approximation does not hold everywhere and its failure leads to nonadiabatic effects, which are relevant to understand much of photo-induced dynamics of molecular systems, especially when conical intersections are involved.
