Born-Oppenheimer approximation and perturbation theory

In the book Molecular Physics by Demtroder there is an explanation of the Born-Oppenheimer approximation and the adiabatic approximation in terms of a perturbative series. The Hamiltonian is $H_0 + T_text{nuc} = H_0 + lambda W$, where $T_text{nuc}$ is the nuclear kinetic energy term. The solutions of $H_0$ are $phi_m(r, R)$ and a general wavefunction is written as a linear combination $psi(r, R) = sum_m chi_m(R) phi_m(r,R)$. The Schroedinger equation $(H_0 + lambda W)psi = Epsi$ is then solved by writing $E$ and $chi_m$ as a series expansion of $lambda$. According to the book, the energy to second order is
$$ E_n = E_n^{(0)} + W_{nn} + sum_{kneq n} frac{W_{nk}W_{kn}}{E_n^{(0)} – E_k^{(0)}}$$
where
$$ W_{nk} = int phi_n^{(0)*} T_text{nuc} phi_k^{(0)} dr.$$

Now I am lost on how this was obtained. They then relate the second term $W_{nn}$ to the adiabatic correction.

Since the energies are labelled, I would assume that then $psi$ would also have to be labelled by $n$.

I would appreciate any suggestions or explanations for this.

Edit: My attempt.

So if I use the $chi$ and $E$ expansion as well as their form of $psi$ then I obtain
$$
(H_0 + lambda W) psi = sum_m chi_m(R) E_m(R) phi_m(r,R) + lambda sum_m phi_m(r,R) W chi_m(R) + chi_m(R) W phi(r,R) = sum_m E chi_m(R)
phi_m(r,R) $$

Now I multiply $phi_k^*(r,R)$ and integrate over $r$ to obtain
$$
E_k(R) chi_k(R) + lambda W chi_k(R) + lambda sum_m chi_m(R) W_{mk} = E chi_k(R).
$$
I should have made clear that in the book the energy was expanded as
$$
E_n = E_n^{(0)} + lambda E_n^{(1)}…$$
Here the energy is labelled so the $psi$ need to be labelled but I will ignore this. If $psi$ is to be normalised then $chi_m(R)^*chi_m(R)$ over $m$ must sum to 1. The zeroth order is then
$$
E^{(0)} = E_k(R)
$$
which is not possible as it depends on $R$. I must be reading something wrong here.