Born-Oppenheimer Approximation reloaded

I'm an undergraduate Chemistry student and I'm currently writing my Final Dissertation.   
 I've written a form of the Born-Oppenheimer approximation with the only initial assumption of parametric dependence on nuclear coordinates. Let me know what do you think about it. Thanks in advance for any help you are able to provide.


For a given molecular system, within the non-relativistic approximation, the molecular Schrodinger equation can be written as $$ \hat{\mathcal{H}}\_{\mathrm{tot}}\Psi\left(\mathbf{R},\,\mathbf{r}\right)=\mathcal{E}\_{\mathrm{tot}}\Psi\left(\mathbf{R},\,\mathbf{r}\right)\quad\left(1\right) $$ where $$\hat{\mathcal{H}}\_{\mathrm{tot}}=\hat{\mathcal{T}}\_{\mathrm{e}}\left(\mathbf{r}\right)+\hat{\mathcal{T}}\_{\mathrm{N}}\left(\mathbf{R}\right)+\hat{\mathcal{V}}\_{\mathrm{ee}}\left(\mathbf{r}\right)+\hat{\mathcal{V}}\_{\mathrm{NN}}\left(\mathbf{R}\right)+\hat{\mathcal{V}}\_{\mathrm{eN}}\left(\mathbf{R},\mathbf{r}\right)=\hat{\mathcal{T}}\_{\mathrm{N}}+\hat{\mathcal{H}}\_{\mathrm{e}}$$


1. Without loss of generality, if there exists a complete basis $ \left\{ \Psi\_{k}\left(\mathbf{R},\,\mathbf{r}\right)\right\} 
$ , the molecular wave function $ \Psi\left(\mathbf{R},\,\mathbf{r}\right) $ can be expanded as $$ \Psi\left(\mathbf{R},\,\mathbf{r}\right)={\displaystyle \sum\_{k}}c\_{k}\Psi\_{k}\left(\mathbf{R},\,\mathbf{r}\right) $$
2. The nuclei in a molecule are heavy enough that, to a good approximation, they behave as classical particles. Thus, it is completely meaningful to talk about nuclear configurations, or in other words to assume a parametric dependence from nuclear coordinates.
$$
\hat{\mathcal{H}}\_{\mathrm{tot}}\Psi\left(\mathrm{R},\,\mathbf{r}\right)=\mathcal{E}\_{\mathrm{tot}}\left(\mathrm{R}\right)\Psi\left(\mathrm{R},\,\mathbf{r}\right)
$$
3. For a given complete basis set $ \left\{\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right\} $, the molecular wavefunction is determined only by the coefficients $ c\_{k} $. Thus, the coefficients $ c\_{k} $
depends on nuclear configuration:
$$
\Psi\left(\mathrm{R},\,\mathbf{r}\right)={\displaystyle \sum\_{k}}c\_{k}\left(\mathrm{R}\right)\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)\quad\left(2\right)
$$
Now by substituting eq.(1) into eq.(2)


$$
\begin{array}{c}
{\displaystyle \sum\_{k}}\left(\hat{\mathcal{T}}\_{\mathrm{N}}+\hat{\mathcal{H}}\_{\mathrm{e}}\right)c\_{k}\left(\mathrm{R}\right)\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)=\mathcal{E}\_{\mathrm{tot}}{\displaystyle \sum\_{k}}c\_{k}\left(\mathrm{R}\right)\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)\\
{\displaystyle \sum\_{k}}\hat{\mathcal{T}}\_{\mathrm{N}}c\_{k}\left(\mathrm{R}\right)\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)+{\displaystyle \sum\_{k}}c\_{k}\left(\mathrm{R}\right)\hat{\mathcal{H}}\_{\mathrm{e}}\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)=\mathcal{E}\_{\mathrm{tot}}{\displaystyle \sum\_{k}}c\_{k}\left(\mathrm{R}\right)\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)
\end{array} $$ 


4. **First term of eq.(3)**
$$
\begin{array}{c}
{\displaystyle \sum\_{k}}\hat{\mathcal{T}}\_{\mathrm{N}}c\_{k}\left(\mathrm{R}\right)\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)=-{\displaystyle \sum\_{a,k}}\dfrac{\hbar^{2}}{2M\_{a}}\left\{ \nabla\_{\mathrm{N}\_{a}}\cdot\nabla\_{\mathrm{N}\_{a}}\left[c\_{k}\left(\mathrm{R}\right)\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right]\right\} =\\
=-{\displaystyle \sum\_{a,k}}\dfrac{\hbar^{2}}{2M\_{a}}\left\{ \nabla\_{\mathrm{N}\_{a}}\cdot\left[c\_{k}\left(\mathrm{R}\right)\nabla\_{\mathrm{N}\_{a}}\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)+\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)\nabla\_{\mathrm{N}\_{a}}c\_{k}\left(\mathrm{R}\right)\right]\right\} =\\
=-{\displaystyle \sum\_{a,k}}\dfrac{\hbar^{2}}{2M\_{a}}\left\{ c\_{k}\left(\mathrm{R}\right)\nabla\_{\mathrm{N}\_{a}}^{2}\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)+\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)\nabla\_{\mathrm{N}\_{a}}^{2}c\_{k}\left(\mathrm{R}\right)+2\left[\nabla\_{\mathrm{N}\_{a}}c\_{k}\left(\mathrm{R}\right)\nabla\_{\mathrm{N}\_{a}}\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right]\right\} 
\end{array} $$
Multiplying from the left by $ \Psi\_{j}^{\star}\left(\mathrm{R},\,\mathbf{r}\right) $ and integrating over electron coordinates, remembering that $ \left\langle \Psi\_{j}\left(\mathrm{R},\,\mathbf{r}\right)|\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right\rangle =\delta\_{jk} $
$$
\begin{array}{c}
-{\displaystyle \sum\_{a}}\dfrac{\hbar^{2}}{2M\_{a}}\nabla\_{\mathrm{N}\_{a}}^{2}c\_{k}\left(\mathrm{R}\right)-{\displaystyle \sum\_{a,k}}\dfrac{\hbar^{2}}{2M\_{a}}\left\{ c\_{k}\left(\mathrm{R}\right)\left\langle \Psi\_{j}\left(\mathrm{R},\,\mathbf{r}\right)\left|\nabla\_{\mathrm{N}\_{a}}^{2}\right|\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right\rangle +2\nabla\_{\mathrm{N}\_{a}}c\_{k}\left(\mathrm{R}\right)\left\langle \Psi\_{j}\left(\mathrm{R},\,\mathbf{r}\right)\left|\nabla\_{\mathrm{N}\_{a}}\right|\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right\rangle \right\} =\\
=\hat{\mathcal{T}}\_{\mathrm{N}}c\_{k}\left(\mathrm{R}\right)-{\displaystyle \sum\_{a,k}}\dfrac{\hbar^{2}}{2M\_{a}}\left\{ c\_{k}\left(\mathrm{R}\right)\left\langle \Psi\_{j}\left(\mathrm{R},\,\mathbf{r}\right)\left|\nabla\_{\mathrm{N}\_{a}}^{2}\right|\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right\rangle +2\nabla\_{\mathrm{N}\_{a}}c\_{k}\left(\mathrm{R}\right)\left\langle \Psi\_{j}\left(\mathrm{R},\,\mathbf{r}\right)\left|\nabla\_{\mathrm{N}\_{a}}\right|\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right\rangle \right\} 
\end{array}
$$
5. **Second term of eq.(3)**   
 Multiplying by $\Psi\_{j}\left(\mathrm{R},\,\mathbf{r}\right)$ and integrating over electron coordinates
$$
{\displaystyle \sum\_{k}}c\_{k}\left(\mathrm{R}\right)\left\langle \Psi\_{j}\left(\mathrm{R},\,\mathbf{r}\right)\left|\hat{\mathcal{H}}\_{\mathrm{e}}\right|\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right\rangle =\mathcal{E}\_{k}\left(\mathrm{R}\right)c\_{k}\left(\mathrm{R}\right)
$$
6. **Final form of the molecular electronic function** 
$$
\left(\hat{\mathcal{T}}\_{\mathrm{N}}+\mathcal{E}\_{k}\left(\mathrm{R}\right)\right)c\_{k}\left(\mathrm{R}\right)-{\displaystyle \sum\_{a,k}}\dfrac{\hbar^{2}}{2M\_{a}}\left\{ C\_{2a,jk}c\_{k}\left(\mathrm{R}\right)+2C\_{1a,jk}\nabla\_{\mathrm{N}\_{a}}c\_{k}\left(\mathrm{R}\right)\right\} =\mathcal{E}\_{\mathrm{tot}}c\_{k}\left(\mathrm{R}\right)
$$   
 where
$$
\begin{array}{c}
C\_{1a,jk}=\left\langle \Psi\_{j}\left(\mathrm{R},\,\mathbf{r}\right)\left|\nabla\_{\mathrm{N}\_{a}}^{2}\right|\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right\rangle \\
C\_{2a,jk}=\left\langle \Psi\_{j}\left(\mathrm{R},\,\mathbf{r}\right)\left|\nabla\_{\mathrm{N}\_{a}}\right|\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right\rangle 
\end{array}
$$   

By neglecting the coupling terms
$$
\left(\hat{\mathcal{T}}\_{\mathrm{N}}+\mathcal{E}\_{k}\left(\mathrm{R}\right)\right)c\_{k}\left(\mathrm{R}\right)=\mathcal{E}\_{\mathrm{tot}}c\_{k}\left(\mathrm{R}\right)
$$


**Final conclusion**


Within this picture, the $c\_{k}\left(\mathrm{R}\right)$
 behave, in a natural way, as a nuclear eigenfunction.
$$
\begin{array}{c}
\left(\hat{\mathcal{T}}\_{\mathrm{N}}+\mathcal{E}\_{k}\left(\mathrm{R}\right)\right)\Phi\_{k}\left(\mathrm{R}\right)=\mathcal{E}\_{\mathrm{tot}}\Phi\_{k}\left(\mathrm{R}\right)\\
\hat{\mathcal{H}}\_{\mathrm{e}}\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)=\mathcal{E}\_{k}\left(\mathrm{R}\right)\Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)
\end{array}
$$


**Questions**


1. While computing the **Second term of eq.(3)**, I assumed $ \Psi\_{k}\left(\mathrm{R},\,\mathbf{r}\right)$
as an eigenfunction of the electronic Hamiltonian. How can this assumption be demonstrated?
2. Do you think that the **Assumption 3** is correct?




Coincidentally, I'm also currently writing my dissertation and I do also have section discussing Born-Oppenheimer approximation. :) My description is quite similar to yours, so I do not see any critical issues.


With respect to your question, we do not assume that $\Psi\_k(R,r)$ is an an eigenfunction of the electronic Hamiltonian, we know that for sure. If you look at equation (2) then essentially what is going on here is the following. Any well-behaved function of two independent variables $\Psi(R,r)$ can always be expanded over the complete set of functions $\{ g\_{i}(R) \}$ and $\{ h\_{j}(r) \}$ in the following way
$$
\Psi(R,r) = \sum\limits\_{i} \sum\limits\_{j} c\_{ij} g\_{i}(R) h\_{j}(r) \, ,
$$
or by defining
$$
c\_{j}(R) = \sum\limits\_{i} c\_{ij} g\_{i}(R) \, ,
$$
as
$$
\Psi(R,r) = \sum\limits\_{j} c\_{j}(R) h\_{j}(r) \, .
$$
Now as a complete set $\{ h\_{j}(r) \}$ of functions of electronic coordinates we can choose the set of electronic wave functions $\Psi\_k(R,r)$ which are by definition eigenfunctions of electronic Hamiltonian. The fact that these function carry parametric dependence on nuclear coordinates does not matter here.


And yes, this expansion is correct. For more on that story check [this question at Physics.SE](https://physics.stackexchange.com/questions/64320/expanding-two-variable-function-fx-y-over-the-complete-sets-g-ix) I asked about a year ago. 


 4