# Why does exchange integral 1s2p not vanish for Helium excited state?

I am studying first-order perturbation theory applied to the first excited state of Helium. I am facing the following problem:

I understand that the Coulomb integral for all 2p functions will not be $$0$$, but for me, all the exchange integrals involving the 2p functions will be $$0$$. For example, look at this two integrals for the 1s2pz state:

1s function: $$1s=\sqrt{\frac{Z^3}{a_0^3 π}}e^{\frac{-Z}{a_0}r}$$ 2pz function: $$2p_z=\frac{1}4\sqrt{\frac{Z^5}{2a_0^5 π}}re^{\frac{-Z}{2a_0}r}cosθ$$

Coulomb integral: $$J_{p_z}=H_{(1s2pz;1s2pz)} =∫∫\frac{1s(1)2p_z(2)e^21s(1)2p_z (2)}{4πε_0r_{12}}dr_1 dr_2$$ $$J_{p_z}=\frac{Z^8 e^2}{128a_0^8 π^3ε_0} ∫_0^∞e^{{\frac{-2Z}{a_0}}r_1} r_1^2 dr_1 ∫_0^∞e^{{\frac{-2Z}{a_0}}r_2} r_2^4 dr_2 ∫_0^π\frac{senθ_1 dθ_1}{\sqrt{r_1^2+r_2^2-2r_1 r_2 cosθ_1}} ∫_0^{2π}dϕ_1 ∫_0^πcos^2 θ_2 senθ_2 dθ_2 ∫_0^{2π}dϕ_2$$ Exchange integral:

$$K_{p_z}=H_{(1s2pz;1s2pz)} =∫∫\frac{1s(1)2p_z(2)e^22p_z(1)1s(2)}{4πε_0r_{12}}dr_1 dr_2$$ $$K_{p_z}=\frac{Z^8 e^2}{128a_0^8 π^3ε_0} ∫_0^∞e^{{\frac{-3Z}{2a_0}}r_1} r_1^3 dr_1 ∫_0^∞e^{{\frac{-3Z}{2a_0}}r_2} r_2^3 dr_2 ∫_0^π\frac{cosθsenθ_1 dθ_1}{\sqrt{r_1^2+r_2^2-2r_1 r_2 cosθ_1}} ∫_0^{2π}dϕ_1 ∫_0^πcos θ_2 senθ_2 dθ_2 ∫_0^{2π}dϕ_2$$

In all cases, the electron 1 was placed along the z-axis for its $$θ$$ angle to be equal to the angle between the two electrons. So, the exchange integral (K) will be zero because the integral: $$∫_0^πcos θ_2 senθ_2 dθ_2=-\frac{(cos^2π-cos^20)}2=0$$ And the Coulomb Integral will not vanish because: $$∫_0^πcos^2 θ_2 senθ_2 dθ_2=-\frac{(cos^3π-cos^30)}3=\frac{2}3$$ But in all the materials I have read, both Coulomb and exchange integrals are not zero, but they do not show the calculations. I have searched a lot on the internet, and I still have not found where my mistake is.

Introduction

I start off this derivation more generally and analyize the configuration where $$(n_A,\ell_A) \neq (n_B,\ell_B)$$. By $$E_{ n_A\ell_Am_A,n_B\ell_Bm_B}^{Z_1,Z_2}$$. Later, I restrict my attention to the case where $$(n_A,\ell_A) = (1,0)$$.

Throughout this answer, the first two (of four) integrals in the equations tend to be the direct integrals, while the second two (of four) integrals tend to be the exchange integrals.

Results

Let's calculate the energy of where $$(n_A,\ell_A) \neq (n_B,\ell_B)$$. By $$E_{ n_A\ell_Am_A,n_B\ell_Bm_B}^{Z_1,Z_2}$$, I denote the energy of the configuration.
\begin{align*} &E_{ n_A\ell_Am_A,n_B\ell_Bm_B}^{Z_1,Z_2} \\ &= \frac{e^2}{ 4\,\pi\,\epsilon_o } \int_{\mathbb{R}^{3\times 2}} \Psi^*(\mathbf{r}_1, \mathbf{r}_2) \, \frac{1}{ r_{12}} \, \Psi (\mathbf{r}_1, \mathbf{r}_2) \, d^3\mathbf{r}_1\,d^3\mathbf{r}_2 \\ &= \frac{e^2}{ 4\,\pi\,\epsilon_o } \int_{\mathbb{R}^{3\times 2}} u^*_{n_A\ell_Am_A}(\mathbf{r}_1) \, u^*_{n_B\ell_Am_B}(\mathbf{r}_2) \, u_{n_A\ell_Am_A}(\mathbf{r}_1) \, u_{n_B\ell_Am_B}(\mathbf{r}_2) \, \frac{1}{ r_{12}}\,d^3\mathbf{r}_1\,d^3\mathbf{r}_2 \\ & \pm \frac{e^2}{ 4\,\pi\,\epsilon_o } \int_{\mathbb{R}^{3\times 2}} u^*_{n_A\ell_Am_A}(\mathbf{r}_1) \, u^*_{n_B\ell_Am_B}(\mathbf{r}_2) \, u_{n_B\ell_Am_B}(\mathbf{r}_1) \, u_{n_A\ell_Am_A}(\mathbf{r}_2) \, \frac{1}{ r_{12}}\,d^3\mathbf{r}_1\,d^3\mathbf{r}_2 \\ &= \frac{e^2}{ 4\,\pi\,\epsilon_o } \int_{\mathbb{R}^{3\times 2}} u^*_{n_A\ell_Am_A}(\mathbf{r}_1) \, u_{n_A\ell_Am_A}(\mathbf{r}_1) \, u^*_{n_B\ell_Bm_B}(\mathbf{r}_2) \, u_{n_B\ell_Bm_B}(\mathbf{r}_2) \, \frac{1}{ r_{12}}\,d^3\mathbf{r}_1\,d^3\mathbf{r}_2 \\ & \pm \frac{e^2}{ 4\,\pi\,\epsilon_o } \int_{\mathbb{R}^{3\times 2}} u^*_{n_A\ell_Am_A}(\mathbf{r}_1) \, u_{n_B\ell_Bm_B}(\mathbf{r}_1) \, u^*_{n_B\ell_Bm_B}(\mathbf{r}_2) \, u_{n_A\ell_Am_A}(\mathbf{r}_2) \, \frac{1}{ r_{12}}\,d^3\mathbf{r}_1\,d^3\mathbf{r}_2 \end{align*} Putting in the separable wave functions. \begin{align*} &E_{ n_A\ell_Am_A,n_B\ell_Bm_B}^{Z_1,Z_2} \\ &= \frac{e^2}{ 4\,\pi\,\epsilon_o } \int_{\mathbb{R}^{3\times 2}} R^*_{n_A\ell_A}(\mathbf{r}_1) Y^*_{\ell_Am_A}(\mathbf{r}_1) \, R_{n_A\ell_A}(\mathbf{r}_1) Y_{\ell_Am_A}(\mathbf{r}_1) \, R^*_{n_B\ell_B}(\mathbf{r}_2) Y^*_{\ell_Am_B}(\mathbf{r}_2) \, R_{n_B\ell_B}(\mathbf{r}_2) Y_{\ell_Bm_B}(\mathbf{r}_2) \, \frac{1}{ r_{12}}\,d^3\mathbf{r}_1\,d^3\mathbf{r}_2 \\ & \pm \frac{e^2}{ 4\,\pi\,\epsilon_o } \int_{\mathbb{R}^{3\times 2}} R^*_{n_A\ell_A }(\mathbf{r}_1) Y^*_{ \ell_Am_A}(\mathbf{r}_1) \, R_{n_B\ell_B }(\mathbf{r}_1) Y_{ \ell_Bm_B}(\mathbf{r}_1) \, R^*_{n_B\ell_B }(\mathbf{r}_2) Y^*_{ \ell_Bm_B}(\mathbf{r}_2) \, R_{n_A\ell_A }(\mathbf{r}_2) Y_{ \ell_Am_A}(\mathbf{r}_2) \, \frac{1}{ r_{12}}\,d^3\mathbf{r}_1\,d^3\mathbf{r}_2 \end{align*} The computation of this integral is well served by the expansion of $$1/r_{12}$$ in terms of spherical harmonics. $$\frac{1}{ r_{12}} = \begin{cases} \frac{1}{r_2}\sum\limits_{k=0}^\infty \frac{4\pi}{2k+1} \left(\frac{r_{1 } }{r_{ 2} } \right)^{k } \sum\limits_{q=-k}^{k} \,Y^{*}_{k,q}(\theta_1, \phi_1) \,Y _{k,q}(\theta_2, \phi_2), & \text{for}~r_2>r_1;~\text{and} \\ \frac{1}{r_1}\sum\limits_{k=0}^\infty \frac{4\pi}{2k+1} \left(\frac{r_{2 } }{r_{ 1} } \right)^{k }\sum\limits_{q=-k}^{k} Y^{*}_{k,q}(\theta_2, \phi_2) \,Y _{k,q}(\theta_1, \phi_1), & \text{for}~r_1>r_2. \end{cases}$$ \begin{align*} &E_{ n_A\ell_Am_A,n_B\ell_Bm_B}^{Z_1,Z_2} \\ &= \frac{e^2}{ 4\,\pi\,\epsilon_o } \sum\limits_{k=0}^\infty \frac{4\pi}{2k+1} \sum\limits_{q=-k}^{k} \, \int\limits_{0}^\infty \int\limits_{r_1=0}^{r_2} R^*_{n_A\ell_A}(\mathbf{r}_1) \, R_{n_A\ell_A}(r_1) \, R^*_{n_B\ell_B}(r_2) \, R_{n_B\ell_B}(r_2) \, r_{1 }^{k+2}\,r_{ 2}^{1-k } \, dr_1\,dr_2 \\ &\quad\times \int\limits_{0}^\pi \int\limits_{0}^{2\pi} \int\limits_{0}^\pi \int\limits_{0}^{2\pi} Y^*_{\ell_Am_A}(\theta_1,\phi_1) Y_{\ell_Am_A}(\theta_1,\phi_1) Y^*_{\ell_Am_B}(\theta_2,\phi_2) Y_{\ell_Bm_B}(\theta_2,\phi_2) \, Y^{*}_{k,q}(\theta_1, \phi_1) \, Y _{k,q}(\theta_2, \phi_2), \, \sin\theta_1\,d\theta_1\,d\phi_1 \, \sin\theta_2\,d\theta_2\,d\phi_2 \\ &+ \frac{e^2}{ 4\,\pi\,\epsilon_o } \sum\limits_{k=0}^\infty \frac{4\pi}{2k+1} \sum\limits_{q=-k}^{k} \int\limits_{0}^\infty \int\limits_{r_1=r_2}^{\infty} R^*_{n_A\ell_A}(r_1) \, R_{n_A\ell_A}(r_1) \, R^*_{n_B\ell_B}(r_2) \, R_{n_B\ell_B}(r_2) \, r_{2}^{k+2}\,r_{1}^{1-k} \, dr_1\,dr_2 \\ &\quad\times \int\limits_{0}^\pi \int\limits_{0}^{2\pi} \int\limits_{0}^\pi \int\limits_{0}^{2\pi} Y^*_{\ell_Am_A}(\theta_1,\phi_1) Y_{\ell_Am_A}(\theta_1,\phi_1) Y^*_{\ell_Am_B}(\theta_2,\phi_2) Y_{\ell_Bm_B}(\theta_2,\phi_2) \, Y^{*}_{k,q}(\theta_2, \phi_2) \,Y _{k,q}(\theta_1, \phi_1) \, \sin\theta_1\,d\theta_1\,d\phi_1 \, \sin\theta_2\,d\theta_2\,d\phi_2 \\ & \pm \frac{e^2}{ 4\,\pi\,\epsilon_o } \sum\limits_{k=0}^\infty \frac{4\pi}{2k+1} \sum\limits_{q=-k}^{k} \, \int\limits_{0}^\infty \int\limits_{r_1=0}^{r_2} R^*_{n_A\ell_A }(r_1) \, R_{n_B\ell_B }(r_1) \, R^*_{n_B\ell_B }(r_2) \, R_{n_A\ell_A }(r_2) \, r_{1}^{k+2}\,r_{2}^{1-k} \, dr_1\,dr_2 \\ &\quad\times \int\limits_{0}^\pi \int\limits_{0}^{2\pi} \int\limits_{0}^\pi \int\limits_{0}^{2\pi} Y^*_{ \ell_Am_A}(\theta_1,\phi_1) Y_{ \ell_Bm_B}(\theta_1,\phi_1) Y^*_{ \ell_Bm_B}(\theta_2,\phi_2) Y_{ \ell_Am_A}(\theta_2,\phi_2) \, Y^{*}_{k,q}(\theta_1, \phi_1) \,Y _{k,q}(\theta_2, \phi_2), \, \sin\theta_1\,d\theta_1\,d\phi_1 \, \sin\theta_2\,d\theta_2\,d\phi_2 \\ & \pm \frac{e^2}{ 4\,\pi\,\epsilon_o } \sum\limits_{k=0}^\infty \frac{4\pi}{2k+1} \sum\limits_{q=-k}^{k} \int\limits_{0}^\infty \int\limits_{r_1=r_2}^{\infty} R^*_{n_A\ell_A }(r_1) \, R_{n_B\ell_B }(r_1) \, R^*_{n_B\ell_B }(r_2) \, R_{n_A\ell_A }(r_2) \, r_{2}^{k+2}\,r_{1}^{1-k} \, dr_1\,dr_2 \\ &\quad\times \int\limits_{0}^\pi \int\limits_{0}^{2\pi} \int\limits_{0}^\pi \int\limits_{0}^{2\pi} Y^*_{ \ell_Am_A}(\theta_1,\phi_1) Y_{ \ell_Bm_B}(\theta_1,\phi_1) Y^*_{ \ell_Bm_B}(\theta_2,\phi_2) Y_{ \ell_Am_A}(\theta_2,\phi_2) \, Y^{*}_{k,q}(\theta_2, \phi_2) \, Y _{k,q}(\theta_1, \phi_1) \, \sin\theta_1\,d\theta_1\,d\phi_1 \, \sin\theta_2\,d\theta_2\,d\phi_2 \end{align*} Now, since \begin{align} & \int Y_{l_1 m_1}(\theta, \phi) Y_{l_2 m_2}(\theta, \phi) Y_{l_3 m_3}(\theta, \phi)\,\sin\theta\,\mathrm{d}\theta\,\mathrm{d}\phi \\ &\quad = \sqrt{\frac{(2l_1 + 1)(2l_2 + 1)(2l_3 + 1)}{4\pi}} \begin{pmatrix} l_1 & l_2 & l_3 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_1 & l_2 & l_3\\ m_1 & m_2 & m_3 \end{pmatrix} \end{align} we rewrite as \begin{align*} &E_{ n_A\ell_Am_A,n_B\ell_Bm_B}^{Z_1,Z_2} \\ &= \frac{e^2}{ 4\,\pi\,\epsilon_o } \sum\limits_{k=0}^\infty \frac{4\pi}{2k+1} \sum\limits_{q=-k}^{k} \, (-1)^{m_A+m_B+q} \int\limits_{0}^\infty \int\limits_{r_1=0}^{r_2} R^*_{n_A\ell_A}(\mathbf{r}_1) \, R_{n_A\ell_A}(r_1) \, R^*_{n_B\ell_B}(r_2) \, R_{n_B\ell_B}(r_2) \, r_{1 }^{k+2}\,r_{ 2}^{1-k } \, dr_1\,dr_2 \\ &\quad\times \int\limits_{0}^\pi \int\limits_{0}^{2\pi} \int\limits_{0}^\pi \int\limits_{0}^{2\pi} Y_{\ell_A,-m_A}(\theta_1,\phi_1) Y_{k,-q}(\theta_1, \phi_1) Y_{\ell_A,m_A}(\theta_1,\phi_1) Y_{\ell_A,-m_B}(\theta_2,\phi_2) Y_{k,q}(\theta_2, \phi_2), Y_{\ell_B,m_B}(\theta_2,\phi_2) \, \sin\theta_1\,d\theta_1\,d\phi_1 \, \sin\theta_2\,d\theta_2\,d\phi_2 \\ &+ \frac{e^2}{ 4\,\pi\,\epsilon_o } \sum\limits_{k=0}^\infty \frac{4\pi}{2k+1} \sum\limits_{q=-k}^{k} (-1)^{m_A+m_B+q} \int\limits_{0}^\infty \int\limits_{r_1=r_2}^{\infty} R^*_{n_A\ell_A}(r_1) \, R_{n_A\ell_A}(r_1) \, R^*_{n_B\ell_B}(r_2) \, R_{n_B\ell_B}(r_2) \, r_{2}^{k+2}\,r_{1}^{1-k} \, dr_1\,dr_2 \\ &\quad\times \int\limits_{0}^\pi \int\limits_{0}^{2\pi} \int\limits_{0}^\pi \int\limits_{0}^{2\pi} Y_{\ell_A,-m_A}(\theta_1,\phi_1) Y _{k,q}(\theta_1, \phi_1) Y_{\ell_Am_A}(\theta_1,\phi_1) Y_{\ell_A,-m_B}(\theta_2,\phi_2) Y_{k,-q}(\theta_2, \phi_2) Y_{\ell_Bm_B}(\theta_2,\phi_2) \, \sin\theta_1\,d\theta_1\,d\phi_1 \, \sin\theta_2\,d\theta_2\,d\phi_2 \\ & \pm \frac{e^2}{ 4\,\pi\,\epsilon_o } \sum\limits_{k=0}^\infty \frac{4\pi}{2k+1} \sum\limits_{q=-k}^{k} \, (-1)^{m_A+m_B+q} \int\limits_{0}^\infty \int\limits_{r_1=0}^{r_2} R^*_{n_A\ell_A }(r_1) \, R_{n_B\ell_B }(r_1) \, R^*_{n_B\ell_B }(r_2) \, R_{n_A\ell_A }(r_2) \, r_{1}^{k+2}\,r_{2}^{1-k} \, dr_1\,dr_2 \\ &\quad\times \int\limits_{0}^\pi \int\limits_{0}^{2\pi} \int\limits_{0}^\pi \int\limits_{0}^{2\pi} Y_{ \ell_A,-m_A}(\theta_1,\phi_1) Y_{k,-q}(\theta_1, \phi_1) Y_{ \ell_Bm_B}(\theta_1,\phi_1) Y_{ \ell_B,-m_B}(\theta_2,\phi_2) Y_{k,q}(\theta_2, \phi_2), Y_{ \ell_Am_A}(\theta_2,\phi_2) \, \sin\theta_1\,d\theta_1\,d\phi_1 \, \sin\theta_2\,d\theta_2\,d\phi_2 \\ & \pm \frac{e^2}{ 4\,\pi\,\epsilon_o } \sum\limits_{k=0}^\infty \frac{4\pi}{2k+1} \sum\limits_{q=-k}^{k} (-1)^{m_A+m_B+q} \int\limits_{0}^\infty \int\limits_{r_1=r_2}^{\infty} R^*_{n_A\ell_A }(r_1) \, R_{n_B\ell_B }(r_1) \, R^*_{n_B\ell_B }(r_2) \, R_{n_A\ell_A }(r_2) \, r_{2}^{k+2}\,r_{1}^{1-k} \, dr_1\,dr_2 \\ &\quad\times \int\limits_{0}^\pi \int\limits_{0}^{2\pi} \int\limits_{0}^\pi \int\limits_{0}^{2\pi} Y_{ \ell_A,-m_A}(\theta_1,\phi_1) Y_{k,q}(\theta_1, \phi_1) Y_{ \ell_B,m_B}(\theta_1,\phi_1) Y_{ \ell_B,-m_B}(\theta_2,\phi_2) Y_{k,-q}(\theta_2, \phi_2) Y_{ \ell_A,m_A}(\theta_2,\phi_2) \, \sin\theta_1\,d\theta_1\,d\phi_1 \, \sin\theta_2\,d\theta_2\,d\phi_2. \end{align*} Upon substitution, \begin{align*} &E_{ n_A\ell_Am_A,n_B\ell_Bm_B}^{Z_1,Z_2} \\ &= \frac{e^2}{ 4\,\pi\,\epsilon_o } \sum\limits_{k=0}^\infty \frac{4\pi}{2k+1} \sum\limits_{q=-k}^{k} \, (-1)^{m_A+m_B+q} \int\limits_{0}^\infty \int\limits_{r_1=0}^{r_2} R^*_{n_A\ell_A}(\mathbf{r}_1) \, R_{n_A\ell_A}(r_1) \, R^*_{n_B\ell_B}(r_2) \, R_{n_B\ell_B}(r_2) \, r_{1 }^{k+2}\,r_{ 2}^{1-k } \, dr_1\,dr_2 \\ &\quad\times \frac{(2\ell_A + 1) (2k + 1) (2\ell_B + 1)}{4\pi} \begin{pmatrix} \ell_A & k & \ell_A \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell_B & k & \ell_B \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell_A & k & \ell_A \\ -m_A & -q & m_A \end{pmatrix} \begin{pmatrix} \ell_B & k & \ell_B \\ -m_B & q & m_B \end{pmatrix} \\ &+ \frac{e^2}{ 4\,\pi\,\epsilon_o } \sum\limits_{k=0}^\infty \frac{4\pi}{2k+1} \sum\limits_{q=-k}^{k} (-1)^{m_A+m_B+q} \int\limits_{0}^\infty \int\limits_{r_1=r_2}^{\infty} R^*_{n_A\ell_A}(r_1) \, R_{n_A\ell_A}(r_1) \, R^*_{n_B\ell_B}(r_2) \, R_{n_B\ell_B}(r_2) \, r_{2}^{k+2}\,r_{1}^{1-k} \, dr_1\,dr_2 \\ &\quad\times \frac{(2\ell_A + 1)(2k + 1)(2\ell_B + 1)}{4\pi} \begin{pmatrix} \ell_A & k & \ell_A \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell_B & k & \ell_B \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell_A & k & \ell_A \\ -m_A & q & m_A \end{pmatrix} \begin{pmatrix} \ell_B & k & \ell_B \\ -m_B & -q & m_B \end{pmatrix} \\ & \pm \frac{e^2}{ 4\,\pi\,\epsilon_o } \sum\limits_{k=0}^\infty \frac{4\pi}{2k+1} \sum\limits_{q=-k}^{k} \, (-1)^{m_A+m_B+q} \int\limits_{0}^\infty \int\limits_{r_1=0}^{r_2} R^*_{n_A\ell_A }(r_1) \, R_{n_B\ell_B }(r_1) \, R^*_{n_B\ell_B }(r_2) \, R_{n_A\ell_A }(r_2) \, r_{1}^{k+2}\,r_{2}^{1-k} \, dr_1\,dr_2 \\ &\quad\times \frac{(2\ell_A + 1)(2k + 1)(2\ell_B + 1)}{4\pi} \begin{pmatrix} \ell_A & k & \ell_B \\ 0 & 0 & 0 \end{pmatrix}^2 \begin{pmatrix} \ell_A & k & \ell_B \\ -m_A & -q & m_B \end{pmatrix} \begin{pmatrix} \ell_A & k & \ell_B \\ -m_B & q & m_A \end{pmatrix} \\ & \pm \frac{e^2}{ 4\,\pi\,\epsilon_o } \sum\limits_{k=0}^\infty \frac{4\pi}{2k+1} \sum\limits_{q=-k}^{k} (-1)^{m_A+m_B+q} \int\limits_{0}^\infty \int\limits_{r_1=r_2}^{\infty} R^*_{n_A\ell_A }(r_1) \, R_{n_B\ell_B }(r_1) \, R^*_{n_B\ell_B }(r_2) \, R_{n_A\ell_A }(r_2) \, r_{2}^{k+2}\,r_{1}^{1-k} \, dr_1\,dr_2 \\ &\quad\times \frac{(2\ell_A + 1)(2k + 1)(2\ell_B + 1)}{4\pi} \begin{pmatrix} \ell_A & k & \ell_B \\ 0 & 0 & 0 \end{pmatrix}^2 \begin{pmatrix} \ell_A & k & \ell_B \\ -m_A & q & m_B \end{pmatrix} \begin{pmatrix} \ell_A & k & \ell_B \\ m_A & -q & -m_B \end{pmatrix} . \end{align*} From \url{https://en.wikipedia.org/wiki/3-j_symbol}, the Wigner $$3-j$$ symbol is zero unless (i) $$m + q -m = 0$$, and (ii) $$|k-\ell| \leq \ell \leq k+\ell$$. Consequently, the only values of $$k$$ and $$q$$ that result in a non-trivial definite intergal are $$q=0$$ and $$k=0,1,\ldots,\ell$$. I now employ the fact\footnote{From \url{https://mathworld.wolfram.com/Wigner3j-Symbol.html} } that for $$J = 2\ell+k$$ \begin{align*} &E_{ n_A\ell_Am_A,n_B\ell_Bm_B}^{Z_1,Z_2} \\ &= \frac{e^2}{ 4\,\pi\,\epsilon_o } (2\ell_A + 1) (2\ell_B + 1) \sum\limits_{k=0}^\infty (-1)^{m_A+m_B} \int\limits_{0}^\infty \int\limits_{r_1=0}^{r_2} R^*_{n_A\ell_A}(\mathbf{r}_1) \, R_{n_A\ell_A}(r_1) \, R^*_{n_B\ell_B}(r_2) \, R_{n_B\ell_B}(r_2) \, r_{1 }^{k+2}\,r_{ 2}^{1-k } \, dr_1\,dr_2 \\ &\quad\times \begin{pmatrix} \ell_A & k & \ell_A \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell_B & k & \ell_B \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell_A & k & \ell_A \\ -m_A & 0 & m_A \end{pmatrix} \begin{pmatrix} \ell_B & k & \ell_B \\ -m_B & & m_B \end{pmatrix} \\ &+ \frac{e^2}{ 4\,\pi\,\epsilon_o } (2\ell_A + 1) (2\ell_B + 1) \sum\limits_{k=0}^\infty (-1)^{m_A+m_B} \int\limits_{0}^\infty \int\limits_{r_1=r_2}^{\infty} R^*_{n_A\ell_A}(r_1) \, R_{n_A\ell_A}(r_1) \, R^*_{n_B\ell_B}(r_2) \, R_{n_B\ell_B}(r_2) \, r_{2}^{k+2}\,r_{1}^{1-k} \, dr_1\,dr_2 \\ &\quad\times \begin{pmatrix} \ell_A & k & \ell_A \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell_B & k & \ell_B \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell_A & k & \ell_A \\ -m_A & 0 & m_A \end{pmatrix} \begin{pmatrix} \ell_B & k & \ell_B \\ -m_B & 0 & m_B \end{pmatrix} \\ & \pm \frac{e^2}{ 4\,\pi\,\epsilon_o } (-1)^{\ell_A+k+q} \, (2\ell_A + 1) (2\ell_B + 1) \sum\limits_{k=0}^\infty \int\limits_{0}^\infty \int\limits_{r_1=0}^{r_2} R^*_{n_A\ell_A }(r_1) \, R_{n_B\ell_B }(r_1) \, R^*_{n_B\ell_B }(r_2) \, R_{n_A\ell_A }(r_2) \, r_{1}^{k+2}\,r_{2}^{1-k} \, dr_1\,dr_2 \\ &\quad\times \begin{pmatrix} \ell_A & k & \ell_B \\ 0 & 0 & 0 \end{pmatrix}^2 \begin{pmatrix} \ell_A & k & \ell_B \\ -m_A & m_A-m_B & m_B \end{pmatrix}^2 \\ & \pm \frac{e^2}{ 4\,\pi\,\epsilon_o } (-1)^{ \ell_A + k + \ell_B } \, (2\ell_A + 1) (2\ell_B + 1) \sum\limits_{k=0}^\infty \int\limits_{0}^\infty \int\limits_{r_1=r_2}^{\infty} R^*_{n_A\ell_A }(r_1) \, R_{n_B\ell_B }(r_1) \, R^*_{n_B\ell_B }(r_2) \, R_{n_A\ell_A }(r_2) \, r_{2}^{k+2}\,r_{1}^{1-k} \, dr_1\,dr_2 \\ &\quad\times \begin{pmatrix} \ell_A & k & \ell_B \\ 0 & 0 & 0 \end{pmatrix}^2 \begin{pmatrix} \ell_A & k & \ell_B \\ -m_A & m_A-m_B & m_B \end{pmatrix}^2 . \end{align*} And \begin{align*} &E_{ n_A\ell_Am_A,n_B\ell_Bm_B}^{Z_1,Z_2} \\ &= \frac{e^2}{ 4\,\pi\,\epsilon_o } (2\ell_A + 1) (2\ell_B + 1) \sum\limits_{k=0}^\infty (-1)^{m_A+m_B} \begin{pmatrix} \ell_A & k & \ell_A \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell_B & k & \ell_B \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell_A & k & \ell_A \\ -m_A & 0 & m_A \end{pmatrix} \begin{pmatrix} \ell_B & k & \ell_B \\ -m_B & & m_B \end{pmatrix} \\ &\quad \times \begin{bmatrix} +\int\limits_{0}^\infty \int\limits_{r_1=0}^{r_2} R^*_{n_A\ell_A}(\mathbf{r}_1) \, R_{n_A\ell_A}(r_1) \, R^*_{n_B\ell_B}(r_2) \, R_{n_B\ell_B}(r_2) \, r_{1 }^{k+2}\,r_{ 2}^{1-k } \, dr_1\,dr_2 \\ + \int\limits_{0}^\infty \int\limits_{r_1=r_2}^{\infty} R^*_{n_A\ell_A}(r_1) \, R_{n_A\ell_A}(r_1) \, R^*_{n_B\ell_B}(r_2) \, R_{n_B\ell_B}(r_2) \, r_{2}^{k+2}\,r_{1}^{1-k} \, dr_1\,dr_2 \end{bmatrix} \\ & \pm \frac{e^2}{ 4\,\pi\,\epsilon_o } (2\ell_A + 1) (2\ell_B + 1) \sum\limits_{k=0}^\infty (-1)^{ \ell_A + k + \ell_B } \begin{pmatrix} \ell_A & k & \ell_B \\ 0 & 0 & 0 \end{pmatrix}^2 \begin{pmatrix} \ell_A & k & \ell_B \\ -m_A & m_A-m_B & m_B \end{pmatrix}^2 \\ &\quad\times \begin{bmatrix} + \int\limits_{0}^\infty \int\limits_{r_1=0}^{r_2} R^*_{n_A\ell_A }(r_1) \, R_{n_B\ell_B }(r_1) \, R^*_{n_B\ell_B }(r_2) \, R_{n_A\ell_A }(r_2) \, r_{1}^{k+2}\,r_{2}^{1-k} \, dr_1\,dr_2 \\ + \int\limits_{0}^\infty \int\limits_{r_1=r_2}^{\infty} R^*_{n_A\ell_A }(r_1) \, R_{n_B\ell_B }(r_1) \, r_{1}^{1-k} \, dr_1 \, R^*_{n_B\ell_B }(r_2) \, R_{n_A\ell_A }(r_2) \, r_{2}^{k+2} \,dr_2 \end{bmatrix} . \end{align*} Now for the case of interest here which includes one electron in the 1s state, I find that \begin{align*} &E_{ 1,0,0;n_B,\ell_B,m_B}^{Z_1=2,Z_2=2} \\ &= \frac{e^2}{ 4\,\pi\,\epsilon_o } (2\ell_B + 1) \sum\limits_{k=0}^\infty (-1)^{ m_B} \begin{pmatrix} 0 & k & 0 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell_B & k & \ell_B \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} 0 & k & 0 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell_B & k & \ell_B \\ -m_B & 0 & m_B \end{pmatrix} \\ &\quad \times \begin{bmatrix} +\int\limits_{0}^\infty \int\limits_{r_1=0}^{r_2} R^*_{1,0}(\mathbf{r}_1) \, R_{1,0}(r_1) \, R^*_{n_B\ell_B}(r_2) \, R_{n_B\ell_B}(r_2) \, r_{1 }^{k+2}\,r_{ 2}^{1-k } \, dr_1\,dr_2 \\ + \int\limits_{0}^\infty \int\limits_{r_1=r_2}^{\infty} R^*_{1,0}(r_1) \, R_{1,0}(r_1) \, R^*_{n_B\ell_B}(r_2) \, R_{n_B\ell_B}(r_2) \, r_{2}^{k+2}\,r_{1}^{1-k} \, dr_1\,dr_2 \end{bmatrix} \\ & \pm \frac{e^2}{ 4\,\pi\,\epsilon_o } (2\ell_B + 1) \sum\limits_{k=0}^\infty (-1)^{ k + \ell_B } \begin{pmatrix} 0 & k & \ell_B \\ 0 & 0 & 0 \end{pmatrix}^2 \begin{pmatrix} 0 & k & \ell_B \\ 0 & -m_B & m_B \end{pmatrix}^2 \\ &\quad\times \begin{bmatrix} + \int\limits_{0}^\infty \int\limits_{r_1=0}^{r_2} R^*_{1,0}(r_1) \, R_{n_B\ell_B }(r_1) \, R^*_{n_B\ell_B }(r_2) \, R_{1,0 }(r_2) \, r_{1}^{k+2}\,r_{2}^{1-k} \, dr_1\,dr_2 \\ + \int\limits_{0}^\infty \int\limits_{r_1=r_2}^{\infty} R^*_{1,0 }(r_1) \, R_{n_B\ell_B }(r_1) \, r_{1}^{1-k} \, dr_1 \, R^*_{n_B\ell_B }(r_2) \, R_{1,0 }(r_2) \, r_{2}^{k+2} \,dr_2 \end{bmatrix} . \end{align*} According to the rules for the Wigner-$$3j$$ numbers, the only values of $$k$$ that result in non-trivial values here are $$k=0$$ for the direct integral, and $$k=\ell_B$$ for the exchange integral. Thus,
\begin{align*} &E_{ 1,0,0;n_B,\ell_B,m_B}^{Z_1,Z_2} \\ &= \frac{e^2}{ 4\,\pi\,\epsilon_o } (2\ell_B + 1) (-1)^{ m_B} \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell_B & 0 & \ell_B \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell_B & 0 & \ell_B \\ -m_B & 0 & m_B \end{pmatrix} \\ &\quad \times \begin{bmatrix} +\int\limits_{0}^\infty \int\limits_{r_1=0}^{r_2} R^*_{1,0}(\mathbf{r}_1) \, R_{1,0}(r_1) \, R^*_{n_B\ell_B}(r_2) \, R_{n_B\ell_B}(r_2) \, r_{1 }^{ 2}\,r_{ 2}^{1 } \, dr_1\,dr_2 \\ + \int\limits_{0}^\infty \int\limits_{r_1=r_2}^{\infty} R^*_{1,0}(r_1) \, R_{1,0}(r_1) \, R^*_{n_B\ell_B}(r_2) \, R_{n_B\ell_B}(r_2) \, r_{2}^{ 2}\,r_{1}^{1 } \, dr_1\,dr_2 \end{bmatrix} \\ & \pm \frac{e^2}{ 4\,\pi\,\epsilon_o } (2\ell_B + 1) \begin{pmatrix} 0 & \ell_B & \ell_B \\ 0 & 0 & 0 \end{pmatrix}^2 \begin{pmatrix} 0 & \ell_B & \ell_B \\ 0 & -m_B & m_B \end{pmatrix}^2 \\ &\quad\times \begin{bmatrix} + \int\limits_{0}^\infty \int\limits_{r_1=0}^{r_2} R^*_{1,0}(r_1) \, R_{n_B\ell_B }(r_1) \, R^*_{n_B\ell_B }(r_2) \, R_{1,0 }(r_2) \, r_{1}^{\ell_B+2}\,r_{2}^{1-\ell_B} \, dr_1\,dr_2 \\ + \int\limits_{0}^\infty \int\limits_{r_1=r_2}^{\infty} R^*_{1,0 }(r_1) \, R_{n_B\ell_B }(r_1) \, r_{1}^{1-\ell_B} \, dr_1 \, R^*_{n_B\ell_B }(r_2) \, R_{1,0 }(r_2) \, r_{2}^{\ell_B+2} \,dr_2 \end{bmatrix} . \end{align*} Reducing further \begin{align} \label{eq:helium-energy-levels-l2-equals-0} \boxed{ E_{ 1,0,0;n_B,\ell_B,m_B}^{Z_1,Z_2} = \begin{matrix} +\frac{e^2}{ 4\,\pi\,\epsilon_o } \int\limits_{0}^\infty \int\limits_{r_1=0}^{r_2} R^{Z_1=2}_{1,0}(\mathbf{r}_1) \, R^{Z_1=2}_{1,0}(r_1) \, R^{Z_2=2}_{n_B\ell_B}(r_2) \, R^{Z_2=2}_{n_B\ell_B}(r_2) \, r_{1 }^{ 2}\,r_{ 2}^{1 } \, dr_1\,dr_2 \\[2em] + \frac{e^2}{ 4\,\pi\,\epsilon_o } \int\limits_{0}^\infty \int\limits_{r_1=r_2}^{\infty} R^{Z_1}_{1,0}(r_1) \, R^{Z_1}_{1,0}(r_1) \, R^{Z_2}_{n_B\ell_B}(r_2) \, R^{Z_2}_{n_B\ell_B}(r_2) \, r_{2}^{ 2}\,r_{1}^{1 } \, dr_1\,dr_2 \\[2em] \pm \frac{e^2}{ 4\,\pi\,\epsilon_o } \dfrac{\int\limits_{0}^\infty \int\limits_{r_1=0}^{r_2} R^{Z_1}_{1,0}(r_1) \, R^{Z_2}_{n_B\ell_B }(r_1) \, R^{Z_2}_{n_B\ell_B }(r_2) \, R^{Z_1}_{1,0 }(r_2) \, r_{1}^{\ell_B+2}\,r_{2}^{1-\ell_B} \, dr_1\,dr_2 }{ (2\ell_B + 1) } \\[2em] \pm \frac{e^2}{ 4\,\pi\,\epsilon_o } \dfrac{\int\limits_{0}^\infty \int\limits_{r_1=r_2}^{\infty} R^{Z_1}_{1,0 }(r_1) \, R^{Z_2}_{n_B\ell_B }(r_1) \, r_{1}^{1-\ell_B} \, dr_1 \, R^{Z_2}_{n_B\ell_B }(r_2) \, R^{Z_2}{1,0 }(r_2) \, r_{2}^{\ell_B+2} \,dr_2 }{ (2\ell_B + 1) } \end{matrix} . } \end{align} In the boxed solution above, the first two (of four) integrals are the direct integrals, while the second two (of four) integrals are the exchange integrals.

Discussion

If you look at Appendix B of Atomic Physics by Mark Fox (Cambridge University Press, 2018), it appears that $$Z_1 = Z_2 = 2$$. However, in Atomic Physics by Christopher Foot, Foot argues that $$Z_1 = Z_2 = 2$$ is correct when calculating the 1s2s configuration. However, $$Z_1 = 2$$ and $$Z_2 = 1$$ is correct when, for example, calculating the 1s2p configuration. Further, Foot states (Oxford University Pres, page 55, 2005) that for the 1s2p configuration, upon plugging in the hydrogenic wave functions $$Z_1 = 2$$ and $$Z_2 = 1$$, that the two direct integrals yield $$-\frac{e^2}{ 8\,\pi\,\epsilon_o\,a_o }\frac{{Z_{2}}^5\,\left(108\,{Z_{1}}^2+16\,Z_{1}\,Z_{2}+{Z_{2}}^2\right)}{2\,{\left(2\,Z_{1}+Z_{2}\right)}^6} = -0.028~\text{eV}.$$ Using the boxed equation above, I obtain that same value. Foot also states (Oxford University Pres, page 56, 2005) that for the 1s2p configuration, upon plugging in the hydrogen wave functions ($$Z_1 = 2$$, $$Z_2 = 1$$), that the two exchange integrals yield $$+0.10$$ eV (page 55, 2005 edition). However, using the boxed equation above, I obtain a different value; namely $$\frac{e^2}{8\pi\epsilon_o\,a_o} \frac{2^{10}\,{Z_{1}}^2\,{Z_{2}}^2\,\left({Z_{1}}^2+5\,Z_{1}\,Z_{2}+{Z_{2}}^2\right)}{3^8\,{\left(Z_{1}+Z_{2}\right)}^5}' = 0.52~\text{eV}.$$ In both of my computations, I use the Bohr radius $$a_0 = 52.9~\text{pm}$$. I do not understand why I get a separation 5 times greater than what I see in the book and on energy-level diagrams for helium.