## Finite Spherical Square Well with spin-orbital potential

Last time, we studied the finite spherical square well and calculated the energy levels for difference angular momentum. In that case, the good quantum numbers (quantities that is commute with the Hamiltonian, i.e. conserved, invariance with time) are the principle quantum number $n$, which dictated the number of node and the angular momentum $l$, and the spin state $s$. A eigen state can be denoted as $|nl m_l s m_s\rangle$. The degeneracy for each state is $2(2l+1)$.

With the spin-orbital potential

$\displaystyle V_{so} = a L \cdot S$

The angular momentum $m_l$ and spin $m_s$ are not good quantum numbers. The sum, total angular momentum $j = l + s$ and it z-component $m_j$ are. The new eigen state is $|n j m_j l s\rangle$.

The Hamiltonian inside the well is

$\displaystyle H = H_0 + V_{so} = -\frac{\hbar^2}{2m} \frac{1}{r^2} \frac{d}{dr}r^2\frac{d}{dr} - |V_0| + \frac{\hbar^2}{2m} \frac{1}{r^2}L^2 + a L\cdot S$

The energy for $H_0$ dose not change, as the $L^2$ is still a good quantum number. The spin-orbital coupling is evaluated like

$\displaystyle L\cdot S = \frac{1}{2}(J^2 - L^2 - S^2)$

then,

$\displaystyle \langle n j m_j l s| L\cdot S |n j m_j l s\rangle = \frac{1}{2} ( j(j+1) - l(l+1) - s(s-1) )$

Since the spin has only 2 spin states, thus,

$\displaystyle \langle L\cdot S \rangle \\ = \frac{1}{2} l , j = l+\frac{1}{2} \\ = -\frac{1}{2}(l+1), j = l - \frac{1}{2}$

Note that, the mean energy of the L-states are unchanged. The degeneracy for $j_+ = l+1/2$ is $2(l+1)$, for $j_- = l-1/2$ is $2l$. Also, the s-states are not affected, because $l = 0$.

Thus, the wave vector inside the well $k$ becomes

$\displaystyle k = \sqrt{\frac{2m}{\hbar^2 }(|V_0|- \beta \langle L\cdot S\rangle) - |E|)}$

the wave vector outside the well $k'$

$\displaystyle k' = i \kappa = i \sqrt{\frac{2m}{\hbar^2 }(\beta \langle L\cdot S\rangle) + |E|)}$

And the rest is matching the boundary condition.

We can also add one more term, an additional $\gamma L^2$ term so that the centroid of the L-state is shifted. The energy for this potential is $\gamma l(l+1)$.

Here is the energy levels for $\beta = -1$ and $\gamma = -0.5$.

We can see, we reproduced the shell ordering 0s1/2, 0p3/2, 0p1/2, 0d5/2, 1s1/2, 0d3/2, 0f7/2, 1p3/2, 1p1/2, 0f5/2….

The magic number 2, 8, 20, and 40 are reproduced. May be, if I adjust the strength of the spin-orbit coupling, my get the magic number 28 and 50. Since the energy level of the s-states are fixed, I need to adjust the shift of the centroid of the others, and the splitting.

Here is the result of $\beta = -1.2, \gamma = -0.2$

It seems that the magic number of 28, which is from the isolation of 0f7/2, is recreated. And the isolation of 0g9/2, created the magic number of 50.

## 1-D Harmonic Oscillator

It is quite surprised that This is not here!

The Hamiltonian is

$\displaystyle H=\frac{P^2}{2m} + \frac{m \omega^2}{2} X^2$

Lets define creation operator $A^\dagger$ and destruction operator $A$, or the ladder operators, such that,

$\displaystyle A^\dagger =\sqrt{\frac{1}{2m}} P + i \sqrt{\frac{m \omega^2}{2}} X$

$\displaystyle A =\sqrt{\frac{1}{2m}} P - i \sqrt{\frac{m \omega^2}{2}} X$

Then we can see

$\displaystyle H = A^\dagger A + \frac{\hbar \omega}{2} = A A^\dagger - \frac{\hbar \omega}{2}$

and the commutative relations

$\displaystyle A^\dagger A - A A^\dagger = -\hbar \omega$

$\displaystyle HA^\dagger - A^\dagger H = \hbar \omega A^\dagger$

$\displaystyle HA - A H = -\hbar \omega A$

From the last two equations, we can see

$\displaystyle H |n\rangle = E_n |n\rangle \\ A H|n\rangle = (HA+\hbar \omega A) |n\rangle = E_n A|n\rangle \\ HA|n\rangle = (E_n-\hbar \omega) A|n\rangle$

This shows that, if $|n\rangle$ is a state with energy $E_n$, then $A|n\rangle$ is also a state with energy $E_n - \hbar \omega$. Similar for $A^\dagger |n\rangle$, the energy will be $E_n + \hbar \omega$.

Suppose the ground state is $|0\rangle$, which is is lowest energy state. If we apply the destruction operator, it will be gone, i.e.

$\displaystyle A|0 \rangle = 0 \\ A^\dagger A|0\rangle = 0 \\ \left(H- \frac{\hbar\omega}{2}\right)|0\rangle \\ H|0\rangle = \frac{\hbar\omega}{2}|0\rangle = E_0 |0\rangle$

Thus, the ground state energy is $\hbar \omega /2$!

We can donate

$\displaystyle A^\dagger |n\rangle = U_n |n+1\rangle \\ A |n\rangle = D_n |n-1\rangle$

Using $\displaystyle A^\dagger A - A A^\dagger = -\hbar \omega$, we have

$\displaystyle (A^\dagger A - A A^\dagger ) |n \rangle \\ = A^\dagger D_n |n-1 \rangle - A U_n|n+1 \rangle \\ = (U_{n-1}D_n- U_n D_{n+1} )|n\rangle = - \hbar \omega |n\rangle$

or

$\displaystyle U_n D_{n+1} = U_{n-1}D_n + \hbar \omega$

Consider above equation on $|0\rangle$, such that $A|0\rangle = 0$, thus,

$\displaystyle (A^\dagger A - A A^\dagger ) |0 \rangle = - U_0 D_{1} )|0\rangle = - \hbar \omega |0\rangle$

Thus, $U_0 D_1 = \hbar \omega$, then, in general, we have

$\displaystyle U_n D_{n+1} = (n+1) \hbar \omega$

Then

$\displaystyle A^\dagger A |n\rangle = U_{n-1} D_n |n\rangle \\ A A^\dagger |n\rangle = U_{n} D_{n+1} |n\rangle$

The simplest solution is

$U_n = \sqrt{\hbar \omega}\sqrt{n+1} \\ D_n =\sqrt{\hbar \omega}\sqrt{n}$

And, if we treat the state $|n\rangle$ contains $n$ oscillators, we have a number operator $N = A^\dagger A$, such that

$A^\dagger A |n \rangle = n \hbar \omega |n \rangle$

some people will normalized the ladder operators with $\sqrt{\hbar \omega}$. For me, I don’t want to make the starting point be so “artificial”. We can summarized the 1-D Harmonic oscillator as

## Finite Spherical Square Well

Hello, everyone, in order to calculate deuteron by Hartree-Fock method, I need a basis. The basis of infinite spherical square well is too “rigid”, that it has to “extension” to non-classical region. Beside of the basis of Wood-Saxon potential. The finite spherical square well is a good alternative. The radial equation is basically the same as finite spherical square well.

The potential is

$\displaystyle V(r) = -|V_0|, r\leq a$

$\displaystyle V(r) = 0, r > a$

Within the well, the wave vector is

$\displaystyle k = \sqrt{\frac{2m}{\hbar^2} (|V_0 |-|E|) }$

, outside the well, the wave vector is

$\displaystyle k' = i \kappa = i \sqrt{\frac{2m}{\hbar^2} |E| }$

The solution is spherical Bessel function. Since the Bessel function of the first and second kind are oscillating like sin or cosine function. To form a decay function when $r > a$, we need the Hankel function with complex argument.

$\displaystyle h_n( i \kappa r) = h_n(x) = - (i x)^n \left( \frac{1}{x} \frac{d}{dx}\right)^n \left(\frac{\exp{(-x)}}{x} \right)$

To make it real, we need a factor $i^n$.

The boundary conditions are continuity and differential continuity.

$\displaystyle j_l(ka) = A i^l h_l(i\kappa a)$

$\displaystyle \left(\frac{d}{dr}j_l(kr)\right)_{r=a} = A i^l \left(\frac{d}{dr}h_l(i\kappa a) \right)_{r=a}$

These two conditions solved for two parameters $A$ and $E$. However, I cannot find an analytical solution to the energy $E$.

If the potential depth is 60 MeV for proton. Radius is 1 fm for a light nuclei. Set $\hbar = 1, m = 1$.

The result is follow,

The 1st column is 1s, 1p, 1d, 1f, and 1g. The 2nd column is 2s, 2p, 2d, 2f, and 2g. The 3rd column is 3s, 3p, and 3d.

By compare with infinite square well,

The energies are lower in finite well, because the wave functions can spread-out to non-classical region, so that the wave length is longer and energy is lower.

In this example, the 3d and 2g orbits are bounded (of course, all orbits in infinite well are bound). This is not because of the depth of the well, but the boundary of the well. In other word, to bring down an orbit, the wave function has to spread out, that is connected with the neutron-halo.

## Wave function in momentum space

The wave function often calculated in spatial coordinate. However, in experimental point of view, the momentum distribution can be extracted directly from the experimental data.

The conversion between momentum space and position space is the Fourier transform

$\displaystyle \phi(\vec{k}) = \frac{1}{\sqrt{2\pi}^3} \int Exp\left(-i \vec{k}\cdot \vec{r} \right) \phi(\vec{r}) d\vec{r}$

Using the plane wave expansion

$\displaystyle Exp(i k\cdot r) = \sum_{l=0}^\infty (2l+1) i^l j_l(kr) P_l(\hat{k}\cdot\hat{r})$

or

$\displaystyle Exp(i k\cdot r) = 4\pi \sum_{l=0}^\infty \sum_{m=-l}^{l} i^l j_l(kr) Y_{lm}(\Omega_k) Y_{lm}^{*}(\Omega_r)$

Thus,

$\displaystyle \phi(\vec{k}) = \frac{4\pi}{\sqrt{2\pi}^3} \sum_{l=0}^\infty (-i)^l \sum_{m=-l}^{l} \int j_l(k r) Y_{lm}(\Omega_k) Y_{lm}^*(\Omega) \phi(\vec{r}) r^2 dr d\Omega$

where $j_l (x)$ is spherical Bessel function. Usually, due to conservation of angular momentum, the angular part can be separated from the spatial part. Let assume the wave function in position space is

$\phi(\vec{r}) = \psi(r) Y_{l_r m_r}(\Omega)$

$\phi(\vec{k}) = \psi(p) Y_{l_k m_k}(\Omega_k)$

Then we have

$\displaystyle \psi(k) = \frac{4\pi}{\sqrt{(2\pi)^3}} \int j_l(k r) \psi(r) r^2 dr$

$\displaystyle \phi(\vec{k}) = \psi(k) (-i)^l Y_{lm}(\Omega_k)$

where $l = l_r = l_k, m=m_r = m_k$, due to the orthogonality of spherical harmonics.

For s, p, d, f-state, the spherical Bessel function is

$\displaystyle j_0(x) = \frac{\sin(x)}{x}$

$\displaystyle j_1(x) = \frac{\sin(x)}{x^2} - \frac{\cos(x)}{x}$

$\displaystyle j_2(x) = \sin(x)\frac{2-x^2}{x^3} - \cos(x)\frac{3}{x^2}$

$\displaystyle j_3(x) = \sin(x)\frac{15-6x^2}{x^4} - \cos(x)\frac{15-x^2}{x^3}$

For Hydrogen-like wave function, the non-normalized momentum distribution is

$\displaystyle \psi_{10}(k) = \frac{4 Z^{5/2}}{(k^2 + Z^2)^2}$

$\displaystyle \psi_{20}(k) = \frac{16 \sqrt{2}Z^{5/2}(4k^2-Z^2)}{(4k^2 + Z^2)^3}$

$\displaystyle \psi_{21}(k) =\sqrt{\frac{2}{3}}\frac{64 k Z^{7/2}}{(4k^2 + Z^2)^3}$

$\displaystyle \psi_{30}(k) = \frac{36 \sqrt{3} Z^{5/2} (81k^3-30k^2Z^2+Z^4)}{(9k^2 + Z^2)^4}$

$\displaystyle \psi_{31}(k) = \frac{144 \sqrt{6} Z^{7/2} (9k^3-kZ^2)}{(9k^2 + Z^2)^4}$

Here is the plot for momentum distribution $(\psi_{nl}(k))^2 k^2$.

It is interesting that, the number of node decrease with higher angular momentum. But be-aware that it is only in atomic case, not a universal.

The higher the principle quantum number, the smaller of the spread of momentum. This is because, the spread of position wave function getting larger, and the uncertainty in momentum space will be smaller. This is a universal principle.

We also plot the Hydrogen radial function in here $\psi(r)^2 x^2$, for reference,

## Complete derivation from Schrodinger equation to Laguerre equation for Coulomb potential

The Hamiltonian is

$\displaystyle H = -\frac{\hbar^2}{2m}\nabla^2 - \frac{Ze^2}{4\pi\epsilon_0r}$

$\displaystyle \left( -\frac{\hbar^2}{2m}\frac{1}{r^2}\left(\frac{d}{dr} r^2 \frac{d}{dr} \right) - \frac{Ze^2}{4\pi\epsilon_0 r} + \frac{\hbar^2}{2m} \frac{l(l+1)}{r^2} \right) R(r) = E R(r)$

rearrange, using Atomic unit

$\displaystyle \left( -\frac{1}{2}\frac{1}{r^2}\left(\frac{d}{dr} r^2 \frac{d}{dr} \right) - E -\frac{Z}{ r} + \frac{1}{2} \frac{l(l+1)}{r^2} \right) R(r) = 0$

Since the normalization condition of $R(r)$ is $\int_0^\infty R(r)R(r) r^2 dr = 1$, it is natural to define $u(r) = r R(r)$

$\displaystyle R = \frac{u}{r}, \frac{d}{dr}\left(r^2\frac{dR}{dr}\right) = r \frac{d^2u}{dr^2}$

$\displaystyle -\frac{1}{2}\frac{d^2u}{dr^2} - \left( E + \frac{Z}{ r} - \frac{l(l+1)}{2r^2} \right) u(r)= 0$

Substitute the eigen energy $E_n = - \frac{Z^2}{2n^2}$

$\displaystyle \frac{d^2u}{dr^2} + \left( - \frac{Z^2}{n^2} + \frac{2Z}{ r} - \frac{l(l+1)}{r^2} \right) u(r)= 0$

Set

$\displaystyle x = \frac{2Z}{n} r, \frac{d}{dr} = \frac{d}{dx} \frac{2Z}{n}$

$\displaystyle \left(\frac{2Z}{n}\right)^2\frac{d^2u(x)}{dx^2} - \left( - \frac{Z^2}{n^2} + \frac{4Z^2}{n x} - \frac{4Z^2}{n^2}\frac{l(l+1)}{x^2} \right) u(x)= 0$

Take out the $4Z^2/n^2$

$\displaystyle \frac{d^2u(x)}{dx^2} - \left(\frac{n}{ x} - \frac{1}{4} - \frac{l(l+1)}{x^2} \right) u(x)= 0$

Now, we know the short and long-range behaviour of $u(x)$, assume it to be

$u(x) = \exp\left(-\frac{x}{2} \right) x^{l+1} y(x)$

Then the equation of $y(x)$ is

$\displaystyle x \frac{d^2y}{dx^2} + \left( 2(l+1) - x\right) \frac{dy}{dx} +\left( n - l- 1 \right) y(x) = 0$

This is the Laguerre equation

$\displaystyle x \frac{d^2y}{dx^2} + \left( \alpha+1 - x\right) \frac{dy}{dx} + m y(x) = 0$

where $\alpha = 2l+1$ and $m = n-l -1$. Therefore, the solution of the radial equation is,

$R(r) = \frac{1}{r} A \exp\left(-\frac{x}{2} \right) x^{l+1} L_{n-l-1}^{2l+1}(x)$

where $A$ is normalization factor.

Notice that the Laguerre polynomial is only defined for $m \geq 0$, thus, $n > l$.

## Hartree-Fock method for Helium excited state

There are two excited wavefunctions, singlet and triplet state. In the singlet state, the total spin is zero, while in the triplet state, the total spin is 1.

In any case, the Slater determinant is

$\displaystyle \Psi(x,y) = \frac{1}{\sqrt{2}} \begin{vmatrix} \phi_1(x) & \phi_2(x) \\ \phi_1(y) & \phi_2(y) \end{vmatrix}$

For spin-singlet state, we set,

$\phi_1(q) = \uparrow \phi_{1s}(r), \phi_2(q) = \downarrow \phi_{2s}(r)$

For triplet state,

$\phi_1(q) = \uparrow \phi_{1s}(r), \phi_2(q) = \uparrow \phi_{2s}(r)$

Since the Hamitonian is spin independent, we only take care of the spatial part.

For the triplet state, if we put in the wavefucntion in the Hamiltonian. Using the basis set approximation, The Fock matrix is

$\displaystyle F_{ij} = H_{ij} + \sum_{hk}^M a_{\nu h}^*a_{\nu k} (G_{ij}^{hk} - \delta(\sigma_h, \sigma_k)G_{ik}^{hj})$

The factor $latex \delta(\sigma_h, \sigma_k)$ is due to the spin component.

From the beginning, the Fock matrix is

$\displaystyle \sum_{j}^M F_{ij} a_{\mu j} = \sum_{j}^M H_{ij} a_{\mu j} + \sum_{j}^M a_{\mu j} \sum_{hk}^M a_{\nu h}^*a_{\nu k} G_{ij}^{hk} - \sum_{j}^M a_{\nu j} \sum_{hk}^M a_{\nu h}^*a_{\mu k} G_{ij}^{hk}$

We can see the exchange term is difference. We can interchange $k \leftrightarrow j$

$\displaystyle \sum_{j}^M a_{\nu j} \sum_{hk}^M a_{\nu h}^*a_{\mu k} G_{ij}^{hk} = \sum_{k}^M a_{\nu k} \sum_{hj}^M a_{\nu h}^*a_{\mu j} G_{ik}^{hj}$

and pull the $a_{\mu j}$ out.

The calculation is straight forward. First, we calculate the 1s state due to the mean field created by the 2s orbital, then calculate the 2s state due to the mean field by the 1s orbital. And then repeat until converged.

Using the Hydrogen 1s, 2s, 3s, and 4s orbtial, we calculated

1s KE+PE = -2 Hartree = -54.42 eV
2s KE+PE = -0.4349 Hartree = -11.83 eV
The direct term = 0.2800 Hartree = 7.62 eV
The exchange term = 0.016 Hartree = 0.44 eV

The total energy = -2.171 Hartree = -59.07 eV.

Compare to the calculated ground state energy ( -77.11 eV), the excitation energy is 18.04 eV. And the ionization energy is 4.65 eV.

The experimental value for the triplet state excitation energy is 19.77 eV, and the ionization energy is 4.70 eV.

Consider the ground state energy is different from the experimental value by 2 eV. the calculation is fairly good.

For the singlet state, because the spin state of the two electron are opposite, the Fock matrix becomes Hartree matrix.

$\displaystyle F_{ij} = H_{ij} + \sum_{hk}^M a_{\nu h}^*a_{\nu k} G_{ij}^{hk}$

The calculated result is,

1s KE+PE = -2.000 Hartree = -54.42 eV
2s KE+PE = -0.400 Hartree = -10.89 eV
The direct term = 0.2531 Hartree = 6.89 eV
The exchange term = 0.040 Hartree = 1.09 eV

The total energy = -2.147 Hartree = -58.43 eV.

Compare to the calculated ground state energy ( -77.11 eV), the excitation energy is 18.68 eV. And the ionization energy is 4.00 eV.

The experimental value for the triplet state excitation energy is 20.55 eV, and the ionization energy is 3.92 eV.

It is interesting to compare with other formalism. Many other textbook use singlet wavefunction,

$\displaystyle \Psi_S(x,y) = \frac{1}{\sqrt{2}}( \phi_1(x) \phi_2(y) + \phi_2(x) \phi_1(y)) \frac{1}{\sqrt{2}}(\uparrow \downarrow - \downarrow \uparrow)$

The triplet state is the same,

$\displaystyle \Psi_T(x,y) = \frac{1}{\sqrt{2}}( \phi_1(x) \phi_2(y) - \phi_2(x) \phi_1(y)) \frac{1}{\sqrt{2}}(\uparrow \uparrow)$

Since the Hamiltonian is spin independent, we can see, when we calculate the total energy of the singlet state and triplet state are

$E_S = E_0 + J + K$

$E_T = E_0 + J - K$

Thus, the energy difference between singlet and triplet state is 2 exchange term.

In our calculation, the wave function of the singlet state is difference. We directly calculate from the Slater determinant. And the Fock matrix is lack of the exchange term due to the opposite spin. However, the triplet state is the same, the exchange term is 0.44 eV. If we do not do Hartree-Fock on the singlet state, but calculate the singlet state using the exchange term, the singlet state energy would be -58.19 eV. Comparing with the Hartree-Fock result (-58.43 eV), I would say they agree. The small difference is due to the difference of the mean fields.

## Hartree method for Helium ground state

After long preparation, I am ready to do this problem.

The two electron in the helium ground state occupy same spacial orbital but difference spin. Thus, the total wavefunction is

$\displaystyle \Psi(x,y) = \frac{1}{\sqrt{2}}(\uparrow \downarrow - \downarrow \uparrow) \psi(x) \psi(y)$

Since the Coulomb potential is spin-independent, the Hartree-Fock method reduce to Hartree method. The Hartree operator is

$F(x) = H(x) + \langle \psi(y)|G(x,y) |\psi(y) \rangle$

where the single-particle Hamiltonian and mutual interaction are

$\displaystyle H(x) = -\frac{\hbar^2}{2m} \nabla^2 - \frac{Ze^2}{4\pi\epsilon_0 x} = -\frac{1}{2}\nabla^2 - \frac{Z}{x}$

$\displaystyle G(x,y) = \frac{e^2}{4\pi\epsilon_0|x-y|} = \frac{1}{|x-y|}$

In the last step, we use atomic unit, such that $\hbar = 1, m=1, e^2 = 4\pi\epsilon_0$. And the energy is in unit of Hartree, $1 \textrm{H} = 27.2114 \textrm{eV}$.

We are going to use Hydrogen-like orbital as a basis set.

$\displaystyle b_i(r) = R_{nl}(r)Y_{lm}(\Omega) \\= \sqrt{\frac{(n-l-1)!Z}{n^2(n+l)!}}e^{-\frac{Z}{n}r} \left( \frac{2Z}{n}r \right)^{l+1} L_{n-l-1}^{2l+1}\left( \frac{2Z}{n} r \right) \frac{1}{r} Y_{lm}(\Omega)$

I like the left the $1/r$, because in the integration $\int b^2 r^2 dr$, the $r^2$ can be cancelled. Also, the $i = nlm$ is a compact index of the orbital.

Using basis set expansion, we need to calculate the matrix elements of

$\displaystyle H_{ij}=\langle b_i(x) |H(x)|b_j(x)\rangle = -\delta \frac{Z^2}{2n^2}$

$\displaystyle G_{ij}^{hk} = \langle b_i(x) b_h(y) |G(x,y) |b_j(x) b_k(y) \rangle$

Now, we will concentrate on evaluate the mutual interaction integral.

Using the well-known expansion,

$\displaystyle G(x,y) = \frac{1}{|x-y|}=\frac{1}{r_{12}} = \sum_{l=0}^{\infty} \sum_{m=-l}^{l} \frac{4\pi}{2l+1} \frac{r_<^l}{r_>^{l+1}} Y_{lm}^{*}(\Omega_1)Y_{lm}(\Omega_2)$

The angular integral

$\displaystyle \langle Y_i(x) Y_h(y)| Y_{lm}^{*}(x) Y_{lm}(y) | Y_j(x) Y_k(y) \rangle \\ = \big( \int Y_i^{*}(x) Y_{lm}^{*}(x) Y_j(x) dx \big) \big( \int Y_h^{*}(y) Y_{lm}(y) Y_k(y) dy \big)$

where the integral $\int dx = \int_{0}^{\pi} \int_{0}^{2\pi} \sin(\theta_x) d\theta_x d\phi_x$.

From this post, the triplet integral of spherical harmonic is easy to compute.

$\displaystyle \int Y_h^{*}(y) Y_{lm}(y) Y_k(y) dy = \sqrt{\frac{(2l+1)(2l_k+1)}{4\pi (2l_h+1)}} C_{l0l_k0}^{l_h0} C_{lm l_km_k}^{l_hm_h}$

The Clebsch-Gordon coefficient imposed a restriction on $l,m$.

$\displaystyle \langle R_i(x) R_h(y)| \frac{r_<^l}{r_>^{l+1}} | R_j(x) R_k(y) \rangle \\ = \int_0^{\infty} \int_{0}^{\infty} R_i(x) R_h(y) \frac{r_<^l}{r_>^{l+1}} R_j(x) R_k(y) y^2 x^2 dy dx \\ = \int_0^{\infty} R_i(x) R_j(x) \\ \left( \int_{0}^{x} R_h(y) R_k(y) \frac{y^l}{x^{l+1}} y^2dy + \int_{x}^{\infty} R_h(x)R_k(x) \frac{x^l}{y^{l+1}} y^2 dy \right) x^2 dx$

The algebraic calculation of the integral is complicated, but after the restriction of $l$ from the Clebsch-Gordon coefficient, only few terms need to be calculated.

The general consideration is done. now, we use the first 2 even states as a basis set.

$\displaystyle b_{1s}(r) = R_{10}(r)Y_{00}(\Omega) = 2Z^{3/2}e^{-Zr}Y_{00}(\Omega)$

$\displaystyle b_{2s}(r) = R_{20}(r)Y_{00}(\Omega) = \frac{1}{\sqrt{8}}Z^{3/2}(2-Zr)e^{-Zr/2}Y_{00}(\Omega)$

These are both s-state orbital. Thus, the Clebsch-Gordon coefficient

$\displaystyle C_{lm l_k m_k}^{l_h m_h} = C_{lm00}^{00}$

The radial sum only has 1 term. And the mutual interaction becomes

$\displaystyle G(x,y) = \frac{1}{|x-y|}=\frac{1}{r_{12}} = 4\pi \frac{1}{r_>} Y_{00}^{*}(\Omega_1)Y_{00}(\Omega_2)$

The angular part

$\displaystyle \langle Y_i(x) Y_h(y)| Y_{lm}^{*}(x) Y_{lm}(y) | Y_j(x) Y_k(y) \rangle = \frac{1}{4\pi}$

Thus, the mutual interaction energy is

$G_{ij}^{hk} = \displaystyle \langle b_i(x) b_h(y) |G(x,y) |b_j(x) b_k(y) \rangle = \langle R_i(x) R_h(y)| \frac{1}{r_>} | R_j(x) R_k(y) \rangle$

$G_{ij}^{hk} = \displaystyle \langle R_i(x) R_h(y)| \frac{1}{r_>} | R_j(x) R_k(y) \rangle \\ \begin{pmatrix} G_{11}^{hk} & G_{12}^{hk} \\ G_{21}^{hk} & G_{22}^{hk} \end{pmatrix} = \begin{pmatrix} \begin{pmatrix} G_{11}^{11} & G_{11}^{12} \\ G_{11}^{21} & G_{11}^{22} \end{pmatrix} & \begin{pmatrix} G_{12}^{11} & G_{12}^{12} \\ G_{12}^{21} & G_{12}^{22} \end{pmatrix} \\ \begin{pmatrix} G_{21}^{11} & G_{21}^{12} \\ G_{21}^{21} & G_{21}^{22} \end{pmatrix} & \begin{pmatrix} G_{22}^{11} & G_{22}^{12} \\ G_{22}^{21} & G_{22}^{22} \end{pmatrix} \end{pmatrix} \\= \begin{pmatrix} \begin{pmatrix} 1.25 & 0.17871 \\ 0.17871 & 0.419753 \end{pmatrix} & \begin{pmatrix} 0.17871 & 0.0438957 \\ 0.0439857 & 0.0171633 \end{pmatrix} \\ \begin{pmatrix} 0.17871 & 0.0438957 \\ 0.0438957 & 0.0171633 \end{pmatrix} & \begin{pmatrix} 0.419753 & 0.0171633 \\ 0.0171633 & 0.300781 \end{pmatrix} \end{pmatrix}$

We can easy to see that $G_{ij}^{hk} = G_{ji}^{hk} = G_{ij}^{kh} = G_{hk}^{ij} = G_{ji}^{kh}$. Thus, if we flatten the matrix of matrix, it is Hermitian, or symmetric.

Now, we can start doing the Hartree method.

The general solution of the wave function is

$\psi(x) = a_1 b_{1s}(x) + a_2 b_{2s}(x)$

The Hartree matrix is

$F_{ij} = H_{ij} + \sum_{h,k} a_h a_k G_{ij}^{hk}$

The first trial wave function are the Hydrogen-like orbital,

$\psi^{(0)}(x) = b_{1s}(r)$

$F_{ij}^{(0)} = \begin{pmatrix} -2 & 0 \\ 0 & -0.5 \end{pmatrix} + \begin{pmatrix} 1.25 & 0.17871 \\ 0.17817 & 0.419753 \end{pmatrix}$

Solve for eigen system, we have the energy after 1st trial,

$\epsilon^{(1)} = -0.794702 , (a_1^{(1)}, a_2^{(1)}) = (-0.970112, 0.242659)$

After 13th trial,

$\epsilon^{(13)} = -0.880049 , (a_1^{(13)}, a_2^{(13)}) = (-0.981015, 0.193931)$

$F_{ij}^{(13)} = \begin{pmatrix} -2 & 0 \\ 0 & -0.5 \end{pmatrix} + \begin{pmatrix} 1.15078 & 0.155932 \\ 0.155932 & 0.408748 \end{pmatrix}$

Thus, the mixing of the 2s state is only 3.7%.

Since the eigen energy contains the 1-body energy and 2-body energy. So, the total energy for 2 electrons is

$E_2 = 2 * \epsilon^{(13)} - G = -2.82364 \textrm{H} = -76.835 \textrm{eV}$

In which ,

$G = \langle \psi(x) \psi(y) |G(x,y) |\psi(x) \psi(y) \rangle = 1.06354 \textrm{H} = 28.9403 \textrm{eV}$

So the energies for

From He to He++.  $E_2 = -2.82364 \textrm{H} = -76.835 \textrm{eV}$
From He+ to He++, $E_1^+ = -Z^2/2 = 2 \textrm{H} = -54.422 \textrm{eV}$.
From He to He+, is $E_1 = E_2 - E_1^+ = -0.823635 \textrm{H} = -22.4123 \textrm{eV}$

The experimental 1 electron ionization energy for Helium atom is

$E_1(exp) = -0.90357 \textrm{H} = -24.587 \textrm{eV}$
$E_1^+(exp) = -1.99982 \textrm{H} = -54.418 \textrm{eV}$
$E_2(exp) = -2.90339 \textrm{H} = -79.005 \textrm{eV}$

The difference with experimental value is 2.175 eV. The following plot shows the Coulomb potential, the screening due to the existence of the other electron, the resultant mean field, the energy, and $r \psi(x)$

Usually, the Hartree method will under estimate the energy, because it neglected the correlation, for example, pairing and spin dependence. In our calculation, the $E_2$ energy is under estimated.

From the $(a_1^{(13)}, a_2^{(13)}) = (-0.981015, 0.193931)$, we can see, the mutual interaction between 1s and 2s state is attractive. While the interaction between 1s-1s and 2s-2s states are repulsive. The repulsive can be easily understood. But I am not sure how to explain the attractive between 1s-2s state.

Since the mass correction and the fine structure correction is in order of $10^{-3} \textrm{eV}$, so the missing 0.2 eV should be due to something else, for example, the incomplete basis set.

If the basis set only contain the 1s orbit, the mutual interaction is 1.25 Hartree = 34.014 eV. Thus, the mixing reduce the interaction by 5.07 eV, just for 3.7% mixing

I included the 3s state,

$\epsilon^{(13)} = -0.888475 , (a_1^{(13)}, a_2^{(13)}, a_3^{(13)}) = (0.981096, -0.181995, -0.06579)$

The mutual energy is further reduced to 1.05415 Hartree = 28.6848 eV. The $E_2 = -77.038 \textrm{eV}$. If 4s orbital included, the $E_2 = -77.1058 \textrm{eV}$. We can expect, if more orbital in included, the $E_2$ will approach to $E_2(exp)$.