Born-Oppenheimer separation in Dirac bra-ket notation Most derivations I have seen of the Born-Oppenheimer approximation are made using wave-functions. To understand it better, I was trying to write a derivation using Dirac notation, but I am stuck. I am going to post what I have done so far so you guys can help me out.
The hamiltonian of the molecule can be written as the sum of two parts, $H_\text{mol} = H_\text{el} + H_\text{nuc}$, where $H_\text{nuc}$ is the hamiltonian of the nuclei by themselves, and $H_\text{el}$ is the hamiltonian of the electrons interacting with the nuclei:
$$H_\text{el} = T_\text{el} + V_\text{el-el} + V_\text{el-nuc}$$
$$H_\text{nuc} = T_\text{nuc} + V_\text{nuc-nuc}$$
We want to find the energy levels of the molecule. That is, we want to solve $H_{\text{mol}} |\mathcal{E}\rangle =\mathcal{E} |\mathcal{E}\rangle$, where $\mathcal{E}$ and $|\mathcal{E}\rangle$ are the eigenenergy and corresponding eigenket of the molecule.
The state space of the molecule can be separated into electronic and nuclear parts: $\mathcal{S}_\text{mol} = \mathcal{S}_\text{el} \otimes \mathcal{S}_\text{nuc}$. Let $|R\rangle \in \mathcal{S}_\text{nuc}$ be the position basis of the nuclei, where $R$ denotes the coordinates of all the nuclei. 
Define $H_{\text{el}}(R) = \langle R| H_{\text{el}} |R\rangle$, which is an operator in $\mathcal{S}_{\text{el}}$. Let $E_a(R)$ and $|E_a(R)\rangle$ be the eigenvalues and corresponding eigenkets of $H_{\text{el}}(R)$ in $\mathcal{S}_{\text{el}}$, so that $H_{\text{el}}(R) |E_a(R)\rangle = E_a(R) |E_a(R)\rangle$. For each $R$, the kets $|E_a(R)\rangle$ make a basis for $\mathcal{S}_{\text{el}}$.
The set of kets $|E_a(R)\rangle |R\rangle \in \mathcal{S}_{\text{mol}}$ for all $R$ and $a$ then make a basis for $\mathcal{S}_{\text{mol}}$. Using this basis, the state of the molecule can be written:
$$|\psi \rangle =\sum _a \int \chi_a(R) |E_a(R)\rangle |R\rangle dR$$
where $\chi_a(R)$ is an amplitude function.
Note that:
$$H_\text{el} |E_a(R)\rangle |R\rangle = [ H_\text{el}(R) |E_a(R)\rangle ] |R\rangle = E_a(R) |E_a(R)\rangle |R\rangle $$
Therefore, the molecular eigenproblem $H_{\text{mol}} |\psi \rangle =\mathcal{E} |\psi \rangle$ can be written:
$$\sum_a \int (E_a + T_\text{nuc} + V_\text{nuc-nuc} - \mathcal{E}) \chi_a(R) |E_a(R)\rangle |R\rangle dR = 0$$
Multiplying by $\langle E_a(R)|$ on the left:
$$\int (E_a + T_\text{nuc} + V_\text{nuc-nuc} - \mathcal{E}) \chi_a(R) |R\rangle dR = 0$$
At last, we define a ket $|\chi_a\rangle \in \mathcal{S}_\text{nuc}$ such that $\chi_a(R) = \langle R | \chi_a \rangle$:
$$|\chi_a\rangle := \int \chi_a(R) |R\rangle dR$$
Then we can write:
$$(T_\text{nuc} + V_\text{nuc-nuc} + E_a - \mathcal{E}) | \chi_a \rangle = 0$$
HERE ENDS MY DERIVATION SO FAR.
I must have done something wrong somewhere, because the final equation that I obtain is, as far as I can tell, the Born-Oppenheimer approximation, but I am obtaining it here as an exact equation. What did I do wrong?
Also, if anyone knows of some reference, textbook or paper, that deals with the Born-Oppenheimer approximation in Dirac notation, please post it.
 A: Your mistake is to multiply the equation
$$\sum_a \int dR \left[ (E_a + T_\text{nuc} + V_\text{nuc-nuc} - \mathcal{E}) \chi_a(R) |E_a(R)\rangle |R\rangle\right] = 0$$
on the left by $\langle E_a(R)|$. The left-hand side is a vector in $\mathcal{S}_\text{mol}$, and while you can project it down to $\mathcal{S}_\text{nuc}$ by multiplying it on the left with a bra from $\mathcal{S}_\text{el}$, the variables $R$ and $a$ are mute in the sense that they do not have a value outside the integral / summation.
Thus, you should multiply on the left by the (slightly but importantly different) bra $\langle E_{a'}(R')|$:
$$\langle E_{a'}(R')|\sum_a \int dR \left[ (E_a + T_\text{nuc} + V_\text{nuc-nuc} - \mathcal{E}) \chi_a(R) |E_a(R)\rangle |R\rangle\right] = 0.$$
Equivalently,
$$(E_a + T_\text{nuc} + V_\text{nuc-nuc} - \mathcal{E}) \sum_a \int dR \left[ \langle E_{a'}(R')|E_a(R)\rangle  \right] \chi_a(R)|R\rangle= 0.$$
(Note that here the dependence of $E_a$ on $R$ is eliminated by making it the operator $E_a=\int dR E_a(R)|R\rangle\langle R|$ on $\mathcal{S}_\text{nuc}$.)
Now, if the vectors $|E_{a}(R)\rangle$ and $|E_{a'}(R')\rangle$ are orthogonal, i.e. if 
$$\langle E_{a'}(R')|E_a(R)\rangle=\delta_{aa'},$$ 
then you can eliminate both the summation and the integral to arrive at the Born-Oppenheimer nuclear TISE. However, there isn't really any reason for this to happen, since the two are electronic eigenstates for different hamiltonians. It is these off-diagonal terms that must be neglected in the Born-Oppenheimer approximation; assuming they are zero is equivalent to postulating the tensor-product structure of the total eigenstates which is at the heart of the BOA. 
In practice, of course, these are close to zero since the molecule doesn't stretch or bend that much and the electronic eigenstates' structure is not altered too much. The point is then to remember that these are approximately zero.
