Yet another derivation of the Born–Oppenheimer approximation

2018-11-02

There are plenty of existing discussions of the Born–Oppenheimer approximation, but none that I've read so far are entirely satisfying. They tend to use confusing notation, conflate operators with their representations, gloss over some crucial aspects, and so on.

The following is my attempt at a succinct derivation of the approximation that touches on all the important details. Specifically, here are some of the questions that arose when I was first learning about this (and that I try to answer below):

• What is the exact nature of the parametric dependence of the fast wavefunctions $\chi_k(x; X)$ on the slow coordinate $X$? Do the fast wavefunctions form an orthonormal basis in some way?
• How does the wavefunction expansion $\sum_k \varphi_k(X) \chi_k(x; X)$ differ from a standard basis expansion? Is it a Schmidt decomposition?
• How can the kinetic energy operator on the slow space result in a derivative of the fast wavefunctions?
• Why do all the surfaces seem to have the same energy?

Hilbert space

Consider a system with degrees of freedom that we will group into "slow" and "fast". They don't actually need to be slow and fast, but these are the labels we will use. The prototypical example is a molecule with slow nuclei and fast electrons. States of this system live in the tensor product Hilbert space $\mathcal{H} = \mathcal{H}^\mathrm{s} \otimes \mathcal{H}^\mathrm{f}.$

On the slow space, we have the (multivariate) continuous representation $| X \rangle$, and on the fast space we have $| x \rangle$; together, that's $| X \, x \rangle$. This doesn't have to be the position representation, but it almost certainly will be.

Parameterized Hamiltonian

For the time being, we'll keep the Hamiltonian fairly general: $\hat{H} = \hat{K}^\mathrm{s} + \hat{K}^\mathrm{f} + \hat{V}.$ We have kinetic energies for the slow and fast degrees of freedom and a potential energy term that operates on the entire space. One requirement that we'll impose is that the potential energy operator must be diagonal in the continuous representation we've chosen: $\hat{V} | X \, x \rangle = V(X, x) | X \, x \rangle.$ This allows us to express the Hamiltonian as $\langle X \, x | \hat{H} = \langle X | \hat{K}^\mathrm{s} \otimes \langle x | + \langle X | \otimes \langle x | \hat{K}^\mathrm{f} + V(X, x) \langle X \, x |.$ This isn't quite in the position representation, since we haven't given the kinetic energy operators a form yet. If they looked like $\sum_i \frac{\partial^2}{\partial \bigstar_i^2},$ we could express the Hamiltonian properly in the continuous representation: $\langle X \, x | \hat{H} = \left( -\sum_i \frac{\partial^2}{\partial X_i^2} - \sum_i \frac{\partial^2}{\partial x_i^2} + V(X, x) \right) \langle X \, x |.$ We'll come back to this form later, so you should keep it in mind, but we'll stick to being more general for now.

We define a parametrized potential operator with the following eigenvalue equation: $\hat{V}^\mathrm{f}(X) | x \rangle = V(X, x) | x \rangle.$ Using this operator, we construct another Hamiltonian, which we will call the fast Hamiltonian: $\hat{H}^\mathrm{f}(X) = \hat{K}^\mathrm{f} + \hat{V}^\mathrm{f}(X).$ The fast Hamiltonian is parameterized by $X$ and acts only on $\mathcal{H}^\mathrm{f}$. Conceptually, this is the Hamiltonian that describes the remaining (fast) system when we freeze out the slow degrees of freedom (by removing $\hat{K}^\mathrm{s}$) and pin them at a specific position $X$ (by parameterizing $\hat{V}$).

For every $X$, the operator $\hat{H}^\mathrm{f}(X)$ is a perfectly legitimate Hamiltonian for the fast system. That means that we could construct the Hamiltonian $\hat{K}^\mathrm{s} + \hat{H}^\mathrm{f}(X),$ but this is a useless object! It describes a complicated system in the fast degrees of freedom and a collection of free particles in the slow degrees of freedom; there is no coupling whatsoever between the two.

What we'll do instead is note that the above definitions allow us to write $\langle X \, x | \hat{H} = \langle X \, x | \left( \hat{K}^\mathrm{s} + \hat{H}^\mathrm{f}(X) \right).$ This may look like we've simply thrown a $\langle X \, x |$ onto the Hamiltonian that we've only just ridiculed, but there's a vital difference: the $X$ parameter of $\hat{H}^\mathrm{f}(X)$ depends on the $X$ value in the bra. This is what gives rise to the coupling between the slow and fast degrees of freedom, and it's at least a little weird to think about.

Parameterized basis

The infinitely many fast Hamiltonians $\hat{H}^\mathrm{f}(X)$ give rise to infinitely many orthonormal bases for $\mathcal{H}^\mathrm{f}$. For any choice of $X$, the states $| k ; X \rangle$ satisfy the eigenvalue equation $\hat{H}^\mathrm{f}(X) | k ; X \rangle = E_k(X) | k ; X \rangle,$ where the kets are also parameterized by $X$. The wavefunctions for these states are commonly written as $\langle x | k ; X \rangle = \chi_k(x; X).$

To be perfectly clear, we have defined a basis $\{ | k ; X \rangle \}_k$ for each value of $X$. There is a basis $\{ | k ; X' \rangle \}_k$, and another basis $\{ | k ; X'' \rangle \}_k$, and so forth; there is nothing we can say in general about the overlap $\langle k' ; X' | k ; X \rangle$ when $X' \ne X$. Given a wavefunction $| \psi \rangle \in \mathcal{H}^\mathrm{f}$, we can expand it as $| \psi \rangle = \sum_k C^\psi_k(X) | k ; X \rangle,$ where the expansion coefficients are given by $C^\psi_k(X) = \langle k ; X | \psi \rangle,$ and $X$ is arbitrary.

Since we're not mathematicians, we can (and will) take continuity for granted. It's fairly safe to assume that the potential $V(X, x)$ varies continuously as $X$ is changed; after all, an arbitrarily large change in the potential when the configuration undergoes an infinitesimal shift would be unphysical. Hence, the Hamiltonian $\hat{H}^\mathrm{f}(X)$ and its eigenfunctions should also be continuous in the parameter $X$, as should the expansion coefficients $C^\psi_k(X)$ for any state.

One wrinkle that we do expect is that funny things can happen at degeneracies. The adiabatic theorem tells us that if we vary $X$ sufficiently slowly (compared to the gap between $E_k(X)$ and adjacent energies $E_{k'}(X)$), then the ordering of the eigenvalues remains the same and we can treat each $k$ as roughly independent. In that case, we have what looks like multiple hypersurfaces in $(X, E)$ space floating above one another like sheets. However, if the energies become equal, the gap between them vanishes, so these sheets touch and cease to be independent. In making the Born–Oppenheimer approximation, we'll be implicitly assuming that this won't happen, so we won't dwell on this.

The $| X \rangle$ representation is complete for $\mathcal{H}^\mathrm{s}$, and every $\{ | k ; X \rangle \}_k$ basis is complete for $\mathcal{H}^\mathrm{f}$. Thus, we are free to pick a specific $X'$ and use the states $| X \rangle \otimes | k ; X' \rangle$ as a basis for $\mathcal{H}$, but this would be silly, since $| k ; X' \rangle$ is not generally an eigenstate of $\hat{H}^\mathrm{f}(X)$ when $X \ne X'$. Instead, we want to use the states $| X \, k \rangle = | X \rangle \otimes | k ; X \rangle,$ where the same $X$ appears in both kets.

To see that $\{ | X \, k \rangle \}_{X,k}$ also forms a basis for $\mathcal{H}$ (technically some sort of half-basis, half-representation mutant), we show that the transformation matrix $U_{k' k}(X', X; \tilde{X}) = \left( \langle X' | \otimes \langle k' ; \tilde{X} | \right) | X \, k \rangle = \langle X' | X \rangle \langle k' ; \tilde{X} | k ; X \rangle$ is unitary for any fixed $\tilde{X}$. The requirements for this are $\begin{aligned} & \int\! \mathrm{d}X \sum_k U_{k' k}(X', X; \tilde{X}) U^*_{k'' k}(X'', X; \tilde{X}) \\ &= \int\! \mathrm{d}X \, \langle X' | X \rangle \langle X | X'' \rangle \sum_k \langle k' ; \tilde{X} | k ; X \rangle \langle k ; X | k'' ; \tilde{X} \rangle \\ &= \int\! \mathrm{d}X \, \langle X' | X \rangle \langle X | X'' \rangle \langle k' ; \tilde{X} | k'' ; \tilde{X} \rangle \\ &= \langle X' | X'' \rangle \langle k' ; \tilde{X} | k'' ; \tilde{X} \rangle \\ &= \delta(X' - X'') \delta_{k' k''} \end{aligned}$ and $\begin{aligned} & \int\! \mathrm{d}X' \sum_{k'} U^*_{k' k''}(X', X''; \tilde{X}) U_{k' k}(X', X; \tilde{X}) \\ &= \int\! \mathrm{d}X' \, \langle X'' | X' \rangle \langle X' | X \rangle \sum_{k'} \langle k'' ; X'' | k' ; \tilde{X} \rangle \langle k' ; \tilde{X} | k ; X \rangle \\ &= \int\! \mathrm{d}X' \, \langle X'' | X' \rangle \langle X' | X \rangle \langle k'' ; X'' | k ; X \rangle \\ &= \langle X'' | X \rangle \langle k'' ; X'' | k ; X \rangle \\ &= \delta(X'' - X) \delta_{k'' k}. \end{aligned}$ In the last step we used the sampling property of the Dirac delta function outside an integral, with the understanding that it only exists inside an integral anyway.

A consequence of the $| X \, k \rangle$ states forming a complete basis is that a state $| \Psi \rangle \in \mathcal{H}$ has the wavefunction $\langle X \, x | \Psi \rangle = \int\! \mathrm{d}X' \sum_k \langle X \, x | X' \, k \rangle \langle X' \, k | \Psi \rangle = \sum_k \langle x | k ; X \rangle \langle X \, k | \Psi \rangle.$ Alternatively, we could write this as $\Psi(X, x) = \langle X \, x | \Psi \rangle = \sum_k \varphi^\Psi_k(X) \chi_k(x; X),$ where $\varphi^\Psi_k(X) = \langle X \, k | \Psi \rangle.$ The last of these is a strange animal, simultaneously serving both the roles of a basis expansion coefficient and a wavefunction for the slow space.

Since it only has a single index, this expansion looks suspiciously like a Schmidt decomposition, but is it one? To express $| \Psi \rangle$ in Schmidt form, we would need to be able to write $\langle X \, x | \Psi \rangle = \sum_j \sqrt{\lambda_j} \langle X | \varphi_j \rangle \langle x | \chi_j \rangle,$ where $| \varphi_j \rangle$ are orthogonal states on $\mathcal{H}^\mathrm{s}$ and $| \chi_j \rangle$ are orthogonal states on $\mathcal{H}^\mathrm{f}$. While our $| k ; X \rangle$ are orthogonal for a fixed $X$, they have an explicit dependence on $X$, which is not allowed. Additionally, $\int\! \mathrm{d}X \, \varphi^{\Psi*}_{k'}(X) \varphi^\Psi_k(X) = \int\! \mathrm{d}X \, \langle \Psi | X \, k' \rangle \langle X \, k | \Psi \rangle \ne \delta_{k' k},$ so these functions aren't even orthogonal.

Conversely, we also have $\langle X \, k | \Psi \rangle = \int\! \mathrm{d}X' \int\! \mathrm{d}x \, \langle X \, k | X' \, x \rangle \langle X' \, x | \Psi \rangle = \int\! \mathrm{d}x \, \langle k ; X | x \rangle \langle X \, x | \Psi \rangle,$ or $\varphi^\Psi_k(X) = \langle X \, k | \Psi \rangle = \int\! \mathrm{d}x \, \chi_k^*(x; X) \Psi(X, x).$

Adiabatic representation

Now that we believe that the states $| X \, k \rangle$ form a basis, we can try to find the matrix elements of the Hamiltonian: $\begin{aligned} \langle X' \, k' | \hat{H} | X \, k \rangle &= \int\! \mathrm{d}x \int\! \mathrm{d}x' \, \langle k' ; X' | x' \rangle \langle X' \, x' | \hat{H} | X \, x \rangle \langle x | k ; X \rangle \\ &= \int\! \mathrm{d}x \int\! \mathrm{d}x' \, \chi_{k'}^*(x'; X') \langle X' \, x' | \hat{K}^\mathrm{s} | X \, x \rangle \chi_k(x; X) \\ &\qquad + \int\! \mathrm{d}x \int\! \mathrm{d}x' \, \chi_{k'}^*(x'; X') \langle X' \, x' | \hat{H}^\mathrm{f}(X') | X \, x \rangle \chi_k(x; X) \\ &= \langle X' | \hat{K}^\mathrm{s} | X \rangle \int\! \mathrm{d}x \, \chi_{k'}^*(x; X') \chi_k(x; X) \\ &\qquad + E_{k'}(X') \langle X' | X \rangle \int\! \mathrm{d}x \, \chi_{k'}^*(x; X') \chi_k(x; X) \\ &= \langle X' | \hat{K}^\mathrm{s} | X \rangle \langle k' ; X' | k ; X \rangle + E_{k'}(X') \delta(X' - X) \delta_{k' k}. \end{aligned}$ We have a complicated kinetic energy term, but a very diagonal potential energy term.

To proceed, we'll choose the form that was mentioned earlier for the kinetic energy: $\langle X | \hat{K}^\mathrm{s} = -\sum_i \frac{\hbar^2}{2 M_i} \frac{\partial^2}{\partial X_i^2} \langle X |.$ Then, it follows that $\begin{aligned} \langle X' | \hat{K}^\mathrm{s} | X \rangle &= -\sum_i \frac{\hbar^2}{2 M_i} \frac{\partial^2}{\partial X_i^2} \langle X' | X \rangle \\ &= -\sum_i \frac{\hbar^2}{2 M_i} \delta^{(2)}_i(X' - X), \end{aligned}$ where $\delta^{(2)}_i(X' - X)$ is the second (distributional) derivative of the delta function in the $i$th direction, which satisfies $\int\! \mathrm{d}X' \, \delta^{(2)}_i(X - X') f(X') = \frac{\partial^2}{\partial X_i^2} f(X).$

Thus, we can apply the Hamiltonian to a generic state $| \Psi \rangle$ as follows: $\begin{aligned} & \langle X \, k | \hat{H} | \Psi \rangle \\ &= \int\! \mathrm{d}X' \sum_{k'} \langle X \, k | \hat{H} | X' \, k' \rangle \langle X' \, k' | \Psi \rangle \\ &= \int\! \mathrm{d}X' \sum_{k'} \langle X | \hat{K}^\mathrm{s} | X' \rangle \langle k ; X | k' ; X' \rangle \langle X' \, k' | \Psi \rangle \\ &\qquad + \int\! \mathrm{d}X' \sum_{k'} E_k(X) \delta(X - X') \delta_{k k'} \langle X' \, k' | \Psi \rangle \\ &= \int\! \mathrm{d}x \, \chi_k^*(x; X) \sum_{k'} \int\! \mathrm{d}X' \langle X | \hat{K}^\mathrm{s} | X' \rangle \chi_{k'}(x; X') \varphi^\Psi_{k'}(X') \\ &\qquad + E_k(X) \varphi^\Psi_k(X) \\ &= -\int\! \mathrm{d}x \, \chi_k^*(x; X) \sum_{k'} \sum_i \frac{\hbar^2}{2 M_i} \int\! \mathrm{d}X' \, \delta^{(2)}_i(X - X') \chi_{k'}(x; X') \varphi^\Psi_{k'}(X') \\ &\qquad + E_k(X) \varphi^\Psi_k(X) \\ &= -\int\! \mathrm{d}x \, \chi_k^*(x; X) \sum_{k'} \sum_i \frac{\hbar^2}{2 M_i} \frac{\partial^2}{\partial X_i^2} \chi_{k'}(x; X) \varphi^\Psi_{k'}(X) \\ &\qquad + E_k(X) \varphi^\Psi_k(X) \\ &= -\int\! \mathrm{d}x \, \chi_k^*(x; X) \sum_{k'} \sum_i \frac{\hbar^2}{2 M_i} \left[ \frac{\partial^2}{\partial X_i^2} \chi_{k'}(x; X) \right] \varphi^\Psi_{k'}(X) \\ &\qquad - \int\! \mathrm{d}x \, \chi_k^*(x; X) \sum_{k'} \sum_i \frac{\hbar^2}{M_i} \left[ \frac{\partial}{\partial X_i} \chi_{k'}(x; X) \right] \left[ \frac{\partial}{\partial X_i} \varphi^\Psi_{k'}(X) \right] \\ &\qquad - \int\! \mathrm{d}x \, \chi_k^*(x; X) \sum_{k'} \sum_i \frac{\hbar^2}{2 M_i} \chi_{k'}(x; X) \left[ \frac{\partial^2}{\partial X_i^2} \varphi^\Psi_{k'}(X) \right] \\ &\qquad + E_k(X) \varphi^\Psi_k(X) \\ &= -\sum_{k'} \sum_i \frac{\hbar^2}{2 M_i} \left[ \int\! \mathrm{d}x \, \chi_k^*(x; X) \frac{\partial^2}{\partial X_i^2} \chi_{k'}(x; X) \right] \varphi^\Psi_{k'}(X) \\ &\qquad - \sum_{k'} \sum_i \frac{\hbar^2}{M_i} \left[ \int\! \mathrm{d}x \, \chi_k^*(x; X) \frac{\partial}{\partial X_i} \chi_{k'}(x; X) \right] \left[ \frac{\partial}{\partial X_i} \varphi^\Psi_{k'}(X) \right] \\ &\qquad - \sum_{k'} \sum_i \frac{\hbar^2}{2 M_i} \left[ \int\! \mathrm{d}x \, \chi_k^*(x; X) \chi_{k'}(x; X) \right] \left[ \frac{\partial^2}{\partial X_i^2} \varphi^\Psi_{k'}(X) \right] \\ &\qquad + E_k(X) \varphi^\Psi_k(X) \\ &= -\sum_{k'} \sum_i \frac{\hbar^2}{2 M_i} \left[ \int\! \mathrm{d}x \, \chi_k^*(x; X) \frac{\partial^2}{\partial X_i^2} \chi_{k'}(x; X) \right] \varphi^\Psi_{k'}(X) \\ &\qquad - \sum_{k'} \sum_i \frac{\hbar^2}{M_i} \left[ \int\! \mathrm{d}x \, \chi_k^*(x; X) \frac{\partial}{\partial X_i} \chi_{k'}(x; X) \right] \left[ \frac{\partial}{\partial X_i} \varphi^\Psi_{k'}(X) \right] \\ &\qquad - \sum_i \frac{\hbar^2}{2 M_i} \left[ \frac{\partial^2}{\partial X_i^2} \varphi^\Psi_k(X) \right] \\ &\qquad + E_k(X) \varphi^\Psi_k(X). \end{aligned}$ We have used the product rule, which states that $\begin{aligned} & \frac{\partial}{\partial X_i} \chi_{k'}(x; X) \varphi^\Psi_{k'}(X) \\ &= \left[ \frac{\partial}{\partial X_i} \chi_{k'}(x; X) \right] \varphi^\Psi_{k'}(X) \\ &\qquad + \chi_{k'}(x; X) \left[ \frac{\partial}{\partial X_i} \varphi^\Psi_{k'}(X) \right] \end{aligned}$ and consequently $\begin{aligned} & \frac{\partial^2}{\partial X_i^2} \chi_{k'}(x; X) \varphi^\Psi_{k'}(X) \\ &= \left[ \frac{\partial^2}{\partial X_i^2} \chi_{k'}(x; X) \right] \varphi^\Psi_{k'}(X) \\ &\qquad + 2 \left[ \frac{\partial}{\partial X_i} \chi_{k'}(x; X) \right] \left[ \frac{\partial}{\partial X_i} \varphi^\Psi_{k'}(X) \right] \\ &\qquad + \chi_{k'}(x; X) \left[ \frac{\partial^2}{\partial X_i^2} \varphi^\Psi_{k'}(X) \right]. \end{aligned}$ More pertinently, we have used the continuously-varying parametric dependence of $| k ; X \rangle$ on $X$ to allow the kinetic energy operator to take its derivative remotely through the derivative of the delta function.

For convenience, we use the gradient vector $\mathbf{\nabla}$ with elements $\nabla_i = \frac{\partial}{\partial X_i},$ so that $\nabla^2 = \mathbf{\nabla} \cdot \mathbf{\nabla} = \sum_i \frac{\partial^2}{\partial^2 X_i},$ and we drop the unitful quantities to make the expressions below look clean. If this makes you feel dirty, don't hesitate to pencil them in where appropriate.

With this in mind, we can write $\langle X \, k | \hat{H} | \Psi \rangle = \sum_{k'} \left[ \left( -\nabla^2 + E_k(X) \right) \delta_{k k'} - 2 \tau^{(1)}_{k k'}(X) \cdot \mathbf{\nabla} - \tau^{(2)}_{k k'}(X) \right] \varphi^\Psi_{k'}(X)$ where $\tau^{(1)}_{k k'}(X) = \int\! \mathrm{d}x \, \chi_k^*(x; X) \mathbf{\nabla} \chi_{k'}(x; X)$ and $\tau^{(2)}_{k k'}(X) = \int\! \mathrm{d}x \, \chi_k^*(x; X) \nabla^2 \chi_{k'}(x; X)$ are the non-adiabatic couplings and the terms containing them are the non-adiabatic coupling terms (NACTs).

Because the derivative operator is antihermitian, we find that $\tau^{(1)}_{k k'}(X) = -\tau^{(1) *}_{k' k}(X)$, so $\tau^{(1)}_{k k',i}(X)$ is a skew-Hermitian matrix (in $k$ and $k'$). A consequence of this is that all its diagonal terms vanish: $\tau^{(1)}_{k k}(X) = 0$. Because the second derivative operator is Hermitian, we also find that $\tau^{(2)}_{k k'}(X) = \tau^{(2) *}_{k' k}(X)$, so $\tau^{(2)}_{k k'}(X)$ is a Hermitian matrix. Hence, all its diagonal terms $\tau^{(2)}_{k k}(X)$ are real.

Note how the terms in the big square brackets smell like a Hamiltonian for the slow degrees of freedom, parameterized by $k$ and $k'$, and expressed in the position representation. If we define the matrix $\hat{\mathbf{H}}$ with elements $\hat{H}_{k k'}$ that have the position representation $\langle X | \hat{H}_{k k'} = \left[ \left( -\nabla^2 + E_k(X) \right) \delta_{k k'} - 2 \tau^{(1)}_{k k'}(X) \cdot \mathbf{\nabla} - \tau^{(2)}_{k k'}(X) \right] \langle X |,$ then the overall Hamiltonian $\hat{H}$ looks like a matrix of Hamiltonians for the slow degrees of freedom, indexed by the surfaces. On the diagonal, we have simply $-\nabla^2 + E_k(X) - \tau^{(2)}_{k k}(X),$ where the last two terms are plain old potentials. On the off-diagonal, we instead have $-2 \tau^{(1)}_{k k'}(X) \cdot \mathbf{\nabla} - \tau^{(2)}_{k k'}(X),$ which is a bit strange, because instead of a second derivative, it has first derivatives.

Nevertheless, this has the effect of turning the single Schrödinger equation $\hat{H} | n \rangle = E_n | n \rangle$ into a collection of coupled differential equations, indexed by $k$: $\sum_{k'} D_{k k'} \varphi^n_{k'}(X) = E_n \varphi^n_k(X),$ where we have given a name to the differential operator $D_{k k'} = \left( -\nabla^2 + E_k(X) \right) \delta_{k k'} - 2 \tau^{(1)}_{k k'}(X) \cdot \mathbf{\nabla} - \tau^{(2)}_{k k'}(X)$ as a shorthand. In the matrix picture, this looks like $\begin{pmatrix} D_{1,1} & D_{1,2} & \cdots \\ D_{2,1} & D_{2,2} & \cdots \\ \vdots & \vdots & \ddots \end{pmatrix} \begin{pmatrix} \varphi^n_1(X) \\ \varphi^n_2(X) \\ \vdots \end{pmatrix} = E_n \begin{pmatrix} \varphi^n_1(X) \\ \varphi^n_2(X) \\ \vdots \end{pmatrix}.$ If one is able to find the functions $\varphi^n_k(X)$ that simultaneously satisfy these equations, one can then assemble the eigenfunction $\langle X \, x | n \rangle = \sum_k \varphi^n_k(X) \chi_k(x; X)$ of the full Hamiltonian $\hat{H}$.

Before we continue, a few brief words about the Hamiltonians $\hat{H}_{k k'}$. It is tempting to say that these are partial matrix elements of $\hat{H}$ in the $| k ; X \rangle$ basis, but that direction is full of potential pitfalls. For starters, which basis do we mean? After all, there is a different one for each $X$, and no matter which one we pick, it would be a mistake to claim that $\langle k ; X | \hat{H} | k' ; X \rangle$ is the object of interest, since its position representation $\langle X' | \langle k ; X | \hat{H} | k' ; X \rangle$ is not useful for us. We could also try $\begin{aligned} \langle X \, k | \hat{H} | k' ; X \rangle &= \int\! \mathrm{d}x \, \langle k ; X | x \rangle \langle X \, x | \hat{H} | k' ; X \rangle \\ &= \int\! \mathrm{d}x \, \langle k ; X | x \rangle \langle X \, x | \hat{K}^\mathrm{s} | k' ; X \rangle \\ &\qquad + \int\! \mathrm{d}x \, \langle k ; X | x \rangle \langle X \, x | \hat{H}^\mathrm{f}(X) | k' ; X \rangle \\ &= \int\! \mathrm{d}x \, \langle k ; X | x \rangle \langle x | k' ; X \rangle \langle X | \hat{K}^\mathrm{s} \\ &\qquad + \int\! \mathrm{d}x \, \langle k ; X | x \rangle \langle x | k' ; X \rangle E_{k'}(X) \langle X | \\ &= \langle X | \left( \hat{K}^\mathrm{s} + E_{k'}(X) \right) \delta_{k k'}, \end{aligned}$ which is definitely not what we wanted. No, this sort of thinking just will not do.

Born–Oppenheimer approximation

Now that we have the Hamiltonian in the adiabatic representation, all that remains is to assume that the non-adiabatic couplings are sufficiently small that neglecting the NACTs entirely is a good approximation. This leaves us with a Hamiltonian that is diagonal in surfaces: $\langle X | \hat{H}_{k k'} = \left( -\nabla^2 + E_k(X) \right) \delta_{k k'} \langle X |.$ In other words, each $\hat{H}_{k k}$ is the complete Hamiltonian for the slow degrees of freedom on surface $k$.

It is then clear that the resulting collection of differential equations is completely uncoupled: $\left( -\nabla^2 + E_k(X) \right) \varphi^n_k(X) = E_n \varphi^n_k(X),$ and they may be solved independently. In the matrix picture, that's $\begin{pmatrix} D_{1,1} & 0 & \cdots \\ 0 & D_{2,2} & \cdots \\ \vdots & \vdots & \ddots \end{pmatrix} \begin{pmatrix} \varphi^n_1(X) \\ \varphi^n_2(X) \\ \vdots \end{pmatrix} = E_n \begin{pmatrix} \varphi^n_1(X) \\ \varphi^n_2(X) \\ \vdots \end{pmatrix}.$

It's somewhat peculiar that all of these seemingly independent Hamiltonians share the same eigenspectrum! It's easy to show that $\tau^{(1)}_{k k',i}(X) = 0$ implies that $\nabla_i \chi_k(x; X) = 0$. The integral $\int\! \mathrm{d}x \, \chi_k^*(x; X) \nabla_i \chi_{k'}(x; X)$ is the overlap between the familiar state $\chi_k^*(x; X)$ and the weird object $\nabla_i \chi_{k'}(x; X)$. Because the $\chi_k(x; X)$ form a complete basis, having $\nabla_i \chi_{k'}(x; X)$ be orthogonal to all $\chi_k(x; X)$ implies that $\nabla_i \chi_{k'}(x; X)$ is the zero element. Hence, the $X_i$ derivative of $\chi_k(x; X)$ is zero; if this is true for all $i$, $\chi_k(x; X)$ doesn't depend on $X$, so we can write just $\chi_k(x)$.

Now we can quickly show in two related ways that the above conclusion isn't a figment of our imagination: $\begin{aligned} \langle X \, k | \hat{H} | n \rangle &= -\int\! \mathrm{d}x \, \chi_k^*(x) \sum_{k'} \sum_i \int\! \mathrm{d}X' \, \delta^{(2)}_i(X - X') \chi_{k'}(x) \varphi^n_{k'}(X') \\ &\qquad + E_k(X) \varphi^n_k(X) \\ &= \left( -\nabla^2 + E_k(X) \right) \varphi^n_k(X) \\ &= E_n \langle X \, k | n \rangle \end{aligned}$ and $\begin{aligned} \langle X \, x | \hat{H} | n \rangle &= \int\! \mathrm{d}X' \sum_k \langle X \, x | \hat{H} | X' \, k \rangle \langle X' \, k | n \rangle \\ &= \sum_k \int\! \mathrm{d}X' \, \langle X \, x | \hat{K}^\mathrm{s} | X' \, k \rangle \langle X' \, k | n \rangle \\ &\qquad + \sum_k \int\! \mathrm{d}X' \, \langle X \, x | \hat{H}^\mathrm{f}(X) | X' \, k \rangle \langle X' \, k | n \rangle \\ &= \sum_k \int\! \mathrm{d}X' \, \langle X | \hat{K}^\mathrm{s} | X' \rangle \langle x | k \rangle \langle X' \, k | n \rangle \\ &\qquad + \sum_k E_k(X) \int\! \mathrm{d}X' \, \langle X | X' \rangle \langle x | k \rangle \langle X' \, k | n \rangle \\ &= -\sum_k \sum_i \int\! \mathrm{d}X' \, \delta^{(2)}_i(X - X') \langle x | k \rangle \varphi^n_k(X') \\ &\qquad + \sum_k E_k(X) \langle x | k \rangle \varphi^n_k(X) \\ &= \sum_k \left( -\nabla^2 + E_k(X) \right) \varphi^n_k(X) \chi_k(x) \\ &= E_n \sum_k \varphi^n_k(X) \chi_k(x) \\ &= E_n \langle X \, x | n \rangle. \end{aligned}$ In fact, because all the $\varphi^n_k(X)$ are degenerate, any linear combination $\sum_k \beta_k \varphi^n_k(X) \chi_k(x)$ is also an eigenstate of $\hat{H}$, including those that only include a single term.

(I think that the above implies that $\mathbf{\nabla} \hat{H}^\mathrm{f}(X) = 0$, so $\mathbf{\nabla} E_k(X) = 0$. Proving this or giving a counterexample is left as an exercise for the reader.)

However, in practice the NACTs don't disappear by themselves. If that were the case, the word "approximation" wouldn't appear on this page. Instead, we create a new Hamiltonian $\hat{H}'$ which has the same diagonal elements as $\hat{H}$, but with the couplings artificially set to zero. In this more realistic scenario, it's not the case that all the surfaces are identical. Still, because they are not coupled in this approximation, they may be dealt with independently.

When the non-adiabatic couplings are sufficiently strong that they can't be neglected, more sophisticated methods are necessary to treat more than one surface at a time; for example, switching to a diabatic representation. This is commonly termed going "beyond the Born–Oppenheimer approximation", and is beyond the scope of the present work.