## Scattering with Coulomb potential

I only stay the result in here. The result is very similar to the case without Coulomb potential.

The general solution for Coulomb potential is

$\displaystyle \Phi_c( r \rightarrow \infty )= e^{i \vec{k} \cdot \vec{r}} e^{-i \eta \ln(k(r-z))} + f_c(\theta) \frac{e^{ikr}}{r} e^{-i \eta \ln(k(r-z))}$

$\displaystyle f_c(\theta) = - \frac{\eta}{ 2k \sin^2(\theta/2)} e^{-i \eta \ln(\sin^2(\theta/2))} e^{2i \sigma_0 }$,

$\displaystyle \sigma_l = \arg( \Gamma(1+l+i\eta) )$

In partial wave,

$\displaystyle \Phi_c = A \frac{1}{2kr} \sum_{l=0}^\infty (2l+1) i^{l+1} P_l(\cos\theta) e^{i\sigma_l} \left( H_l^{-} - H_l^+\right)$

Scattering with Coulomb and short-range potential, the potential is

$\displaystyle V(r) = V_c(r) + V_N(r)$

The solution takes the form

$\Psi(r) = \Phi_c + \Psi_N$

The asymptotic form of $\Psi_N$ is

$\displaystyle \Psi_N \rightarrow A f_N(\theta) \frac{e^{ikr}}{r} e^{-i\eta \ln(2kr)}$

$\displaystyle f_N(\theta) = \frac{1}{2ik} \sum_{l=0}^\infty (2l+1) P_l(\cos\theta) e^{2i\sigma_l} (S_l -1 )$

$S_l = e^{2i \delta_l}$

$\displaystyle \tan \delta_l = - \frac{kR F_l'(kR) - F_l(kR) L^I}{kR G_l'(kR) -G_l(kR) L^I}$

The cross section is

$\displaystyle \frac{d\sigma_c}{d\Omega} = |f_C(\theta) + f_N(\theta)|^2$

## Coulomb wave function (II)

In the previous post, we tried to derived to Coulomb wave function, and the regular Coulomb wave function $F_l(x)$ is derived, except for the normalization constant, and the mysterious Coulomb phase shift.

From this arxiv article, the Coulomb “Hankel” function is

$\displaystyle H_L^{pm}(x) = D_L^{\pm} x^{L+1} e^{\pm i x} U(L+1 \pm i \eta, 2L+2, - \pm 2 i x)$

$\displaystyle D_L^{\pm} = (-\pm 2i)^{2L+1} \frac{\Gamma(L+1\pm i \eta)}{ C_L \Gamma(2L+2)}$

$\displaystyle C_L = z^L \frac{\sqrt{\Gamma(L+1+ i \eta) \Gamma(L+1- i \eta)}}{ e^{\eta \pi/2}\Gamma(2L+2)}$

where $U(a,b,z)$ is the confluent hypergeometric function of the second kind.

The regular Coulomb wave function is

$\displaystyle F_L(x) = C_L x^{L+1} e^{\pm i x} ~_1F_1(L+1 \pm i \eta, 2L+2, - \pm 2 i x) = \frac{1}{2i} \left( H_L^+ - H_L^- \right)$

I am fail to prove the last equality. The irregular Coulomb wave function is

$\displaystyle G_L(x) = \frac{1}{2} \left( H_L^+ + H_L^- \right)$

In Mathematica 11.2, the $U(a,b,z)$ is fail to evaluate when $\eta > 2$ for $l = 0$, and $\eta > 0.1$ for $l = 5$.

Update 20200420, settign WorkingPrecision to 45 can solve the problem. Here are some plots for the Coulomb wave functions.

When $\eta$ getting larger, the wave function pushed further.

## Coulomb wave function

The Coulomb wave function is the solution of a REPULSIVE Coulomb potential. The Schrödinger equation is

$\displaystyle \left( -\frac{\hbar^2}{2m} \nabla^2 + \frac{Q_1Q_2 e^2}{r} \right) \Phi = E \Phi$

where $e^2 = 1.44~\textrm{MeV.fm}$, separate the radial and angular part as usual,

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

setting

$\displaystyle k^2 = \frac{2mE}{\hbar^2}, \eta =\frac{k}{2E} Q_1Q_2 e^2$

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

using the usual substitution $R = u / r$

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

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

Setting $x = kr$

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

The short range behaviour is approximated as

$\displaystyle \frac{d^2u(x)}{dx^2} = \frac{l(l+1)}{x^2}u(x) \Rightarrow u(x) \approx x^{l+1}~\textrm{or}~ x^{-l}$

and for long range behaviour, it should approach Riccati–Bessel functions $\hat{j}_l, \hat{n}_l$ with phase shift.

Set $u(x) = x^{l+1} \exp(i x) y(x)$, the equation of $u(x)$ becomes,

$\displaystyle x \frac{d^2y}{dx^2} + (2L+2 + 2ix) \frac{dy}{dx} + (2i(L+1) - 2\eta) y = 0$

Change of variable $z = - 2i x$,

$\displaystyle z \frac{d^2y}{dz^2} + (2L+2 - z ) \frac{dy}{dx} - (L+1 + i \eta) y = 0$

This is our friend Laguerre polynomial again!!! with $\alpha = 2L+1, n = -(L+1+i\eta)$. Since the $n$ is not an integer anymore, we go to a more general case, that is the Kummer’s equation,

$\displaystyle z \frac{d^2 w}{dz} + (b-z) \frac{dw}{dz} - a w = 0$

The solution of Kummer’s equation is the confluent hypergeometric function

$w(z) = _1F_1(a, b, z)$

Thus, the solution for the radial function is

$\displaystyle u_l(x) =A x^{l+1} e^{ix} _1F_1(L+1+i \eta, 2L+2, -2 i x )$

where $A$ is a normalization constant by compare the long range behaviour with Riccati-Bessel function. The full solution is,

$\displaystyle u_l(x) = \\ F_l(x) = \frac{2^l e^{-\pi \eta/2} |\Gamma(l+1+i\eta)| }{(2l+1)!} x^{l+1} e^{ix}~_1F_1(L+1+i \eta, 2L+2, -2 i x )$

At long range,

$\displaystyle F_l(x \rightarrow \infty) = \sin \left( x - l \frac{\pi}{2} - \eta \log(2x) + \sigma_l \right)$

where $\sigma_l = \arg( \Gamma(l+1+i\eta) )$ is the Coulomb phase shift.

Using Kummer’s transform

$\displaystyle _1F_1(a,b,z) = e^z ~_1F_1(b-a, b, -z)$

The solution can be written as

$\displaystyle u_l(x) = \\ F_l(x) = \frac{2^l e^{-\pi \eta/2} |\Gamma(l+1+i\eta)| }{(2l+1)!} x^{l+1} e^{-ix}~_1F_1(L+1-i \eta, 2L+2, 2 i x )$

Another solution should behave like $\cos$ and unbound at $x = 0$.

Set $u(x) = x^{-l} \exp(i x) y(x)$, the equation of $u(x)$ becomes,

$\displaystyle x \frac{d^2y}{dx^2} + (-2L + 2ix) \frac{dy}{dx} + (-2i L - 2\eta) y = 0$

Change of variable $z = -2ix$,

$\displaystyle z \frac{d^2y}{dz^2} + (-2L -z ) \frac{dy}{dx} - ( - L + i \eta) y = 0$

$\displaystyle u_l(x) =A x^{-l} e^{ix}~ _1F_1(-L+i \eta, -2L, -2 i x )$

However, $_1F_1(a, b, z)$ does not exist for non-position b.

We can also transform into Whittaker’s equation by change of variable $z = 2ix$

$\displaystyle -4 \frac{d^2u(x)}{dz^2} + \left( 1 - \frac{-4l(l+1)}{z^2} - \frac{4i\eta}{z} \right) u(x) =0$

$\displaystyle \frac{d^2u(x)}{dz^2} + \left( -\frac{1}{4} +\frac{i\eta}{z}- \frac{l(l+1)}{z^2} \right) u(x) =0$

using $\displaystyle l(l+1) = (l+1/2)^2 - 1/4$

$\displaystyle \frac{d^2u(x)}{dz^2} + \left( -\frac{1}{4} +\frac{i\eta}{z} + \frac{ 1/4 - (l+1/2)^2}{z^2} \right) u(x) =0$

This is the Whittaker’s equation with $\kappa = i \eta, \mu = l+1/2$

The solutions are

$u_l(x) = e^{-ix} (2ix)^{l+1} ~_1F_1(l+1-i \eta, 2l+2, 2ix)$

$u_l(x) = e^{-ix} (2ix)^{l+1} U(l+1-i \eta, 2l+2, 2ix)$

I still cannot get the second solution, the $G_l$. According to Wolfram, https://mathworld.wolfram.com/CoulombWaveFunction.html

$\displaystyle G_l(x) = \frac{2}{\eta C_0^2(\eta)} F_l(x) \left( \log(2x) + \frac{q_l(\eta)}{p_l(\eta)} \right) + \frac{x^{-l}}{(2l+1) C_l(\eta)} \sum_{K=-l}^\infty a_k^l(\eta) x^{K+l}$,

where $q_l, p_l, a_k^l$ are defined inAbramowitz, M. and Stegun, I. A. (Eds.). “Coulomb Wave Functions.” Ch. 14 in Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 9th printing. New York: Dover, pp. 537-544, 1972.

## Coulomb Energy for two protons system

From this post and this post,  we already have the framework to calculate the energy. The 1-body intersection is zero, as we do not interested on the kinematics energy. And also, suppose we already have the solution for the wave function for a fixed mean-field. Although it is unrealistic, it is a good starting point.

In fact, we have to solve the Schrödinger equation (using Hartree-Fock or numerical solution) with a realistic NN-interaction with and without Coulomb potential, and subtract the difference to get the Coulomb energy.

We have to evaluate the integral

$\displaystyle W = \left<12||12\right> - \left<12||21\right>$

When the 2 protons are in same orbital, the total wavefunction is

$\displaystyle \Psi(1,2) = \frac{1}{\sqrt{2}}\left(\uparrow \downarrow - \downarrow \uparrow \right) \phi(r_1) \phi(r_2)$

The exchange term is zero as the Coulomb operator does not act on the spin part.

$\displaystyle \left<12||12\right> = \sum_{L}^{\infty} \sum_{M=-L}^{L} \frac{4\pi}{2L+1}\int_{0}^{\infty} \int_{0}^{\infty} \phi^2(x) \phi^2(y) \frac{r_<^L}{r_>^{L+1}} x^2 y^2 dx dy \\ \int Y_{lm}^*(1) Y_{LM}^*(1) Y_{lm}(1) d\Omega_1 \int Y_{lm}^*(2) Y_{LM}(2) Y_{lm}(2) d\Omega_2$

For the angular part, from this post or this post, we have

$\displaystyle \int Y_{l_1m_1}^* Y_{LM}^* Y_{l_2m_2} d\Omega = \sqrt{\frac{(2L+1)(2l_1+1)}{4\pi(2l_2+1)}} C_{L0l_10}^{l_20} C_{LMl_1m_1}^{l_2m_2}$

For $l_1 = l_2$, the Clebsh-Gordon coefficient dictated that $2l \geq L$, $L = even$, $M = 0$.

We use the potential is 3D spherical infinite well with radius $a = 2.2$ fm. the wave function is spherical Bessel function ,

$\phi_n(r ) = \begin{matrix} A j_n(k r) & r < a \\ 0 & r > a \end{matrix}, j_n(ka) = 0$,

where $A$ is the normalization factor.

The first non-zero root of spherical Bessel Function take the form,

$r_0(n) = 2.96519 + 0.701921\sqrt{n} + 0.993301 n, n < 500$

The energy for the $n$-th s-state is

$\displaystyle E_n = \frac{\hbar k^2}{2 m_p^2}, k a = r_0(n)$

In Mathematica, the radial calculation can be done as

a=2.2;
R[r_]:=Piecewise[{{SphericalBesselJ[n, r0[n] r], 0 < r < a}, {0, r > 0}}];
A=(NIntegrate[(R[r])^2 r^2, {r, 0, ∞}]])^(-½);
w = NIntegrate[(R[x]R[y])^2 y^L/x^(L+1) x^2 y^2 , {x, 0, ∞}, {y,0, x}] + NIntegrate[(R[x]R[y])^2 x^L/y^(L+1) x^2 y^2 , {x, 0, ∞}, {y,x, ∞}]
W=1.44 A^4 w

The result for first n s-orbital is

I also extract the radius for uniformly charged sphere. The trend is reasonable as higher principle number, the wave function is close to the boundary $a = 2.2 \textrm{fm}$.

Next time, we will work on Woods-Saxon potential.

## Coulomb energy for generic charge distribution

Classically, for an arbitrary charge distribution $\rho(r,\theta, \phi)$, the potential energy at position $\vec{r}_0$ is

$\displaystyle V(r_0) = \int\int\int{\frac{\rho(r,\theta,\phi)}{|\vec{r}-\vec{r}_0|} r^2 \sin\theta d\phi d\theta dr}$

Suppose the charge distribution is spherical symmetric, i.e. $\rho(r,\theta,\phi) = \rho(r )$. And the distance $|\vec{r}-\vec{r}_0|^2 = r^2 + r_0^2 - 2 rr_0 \cos\theta$.

$\displaystyle V(r_0) = 2\pi\int{\rho(r)}r^2\int{\frac{1}{\sqrt{r^2 + r_0^2 - 2 rr_0 \cos\theta}}\sin\theta d\theta} dr$

The angular part,

$\displaystyle \int_{0}^{\pi}{\frac{1}{\sqrt{r^2 + r_0^2 - 2 rr_0 \cos\theta}}\sin\theta d\theta } \\ = \int_{1}^{-1}\frac{1}{\sqrt{r^2+r_0^2 - 2 r r_0 x}} dx \\ = \frac{-1}{2 rr_0} \int_{1}^{-1} \frac{d(r^2+r_0^2 - 2rr_0 x)}{\sqrt{r^2+r_0^2 - 2 r r_0 x}} \\ = \frac{1}{2 rr_0} \int_{(r-r_0)^2}^{(r+r_0)^2} y^{-\frac{1}{2}} dy = \frac{r + r_0 - |r - r_0|}{rr_0}$

The last step is a bit tricky. So, the angular part is

$\displaystyle \int_{0}^{\pi}{\frac{1}{\sqrt{r^2 + r_0^2 - 2 rr_0 \cos\theta}}\sin\theta d\theta } = \frac{r + r_0 - |r - r_0|}{rr_0} = \frac{2}{r_>}$

Recall that

$\displaystyle \frac{1}{|\vec{r}-\vec{r}_0|} = \sum_{L}^{\infty}\sum_{M=-L}^{L} \frac{4\pi}{2L+1} \frac{r_<^L}{r_>^{L+1}} Y_{LM}^*(\Omega)Y_{LM}(\Omega_0)$

For $L = 0$, that reduced to $4\pi/r_>$.

Thus, for spherical arbitrary charge distribution, the Coulomb potential is,

$\displaystyle V(r_0) = 4\pi\int{\rho(r)} \frac{r^2}{r_>} dr$

Lets redo the Coulomb energy for the uniform charge sphere here.

The charge distribution is

$\displaystyle \rho_0(x) = Z \frac{3}{4\pi} \frac{1}{R^3} , x < R$

The potential at position $y$ is

$\displaystyle V(y) = 4\pi \int_0^y \rho_0 (x) \frac{x^2}{y}dx + 4\pi \int_y^\infty \rho_0(x) x dx$

For $y < R$

$\displaystyle V(y) = Z \frac{3}{R^3} \left( \int_0^y \frac{x^2}{y}dx + \int_y^R x dx \right) = Z \frac{3}{R^3} \left( \frac{y^2}{3} + \frac{R^2}{2} - \frac{y^2}{2} \right) \\ = \frac{Z}{R} \left( \frac{3}{2} - \frac{y^2}{2 R^2} \right)$

For $y > R$

$\displaystyle V(y) = 4\pi \int_0^R \rho_0 (x) \frac{x^2}{y}dx = \frac{Z}{y}$

The result is the same.

The Classical Coulomb energy for another charge distribution $\rho_0(\vec{r}_0)$ is

$\displaystyle W = \int V(r_0) \rho_0(\vec{r}_0) d\vec{r}_0$

Quantum Mechanically, the Coulomb energy between two wavefunction,

$\displaystyle W = \left<\psi_1(r_1) \psi_2(r_2)|V_{12}| \psi_1(r_1) \psi_2(r_2) \right> \\ = \int e \psi_1^2(r_1) \left(\int \psi_2^2(r_2) \frac{e}{|\vec{r}_1 - \vec{r}_2|} d\vec{r}_2 \right) d\vec{r}_1 \\ = \int \rho_1(r_1) V(r_1) d\vec{r}_1$

We can see the expression is consistent with the Classical case.

For fermion N-body, the total wave function (in the 3rd part of this post) is,

$\displaystyle \Psi = \frac{1}{\sqrt{N!}}\sum_{n}^{N!} (-1)^{p_n} P_n \phi_a(1)\phi_b(2)...\phi_z(N) = (N!)^{-1/2}\sum_{n}^{N!} (-1)^{p_n} P_n \phi_H$

where $P_n$ is the permutation operator with parity $p_n$.

The Column energy is

$\displaystyle W = \left<\Psi| V_{12} + V_{23} +.... V_{N-1,N} |\Psi \right> = \frac{N(N-1)}{2} \left<\Psi|V_{12}|\Psi \right>$

$\displaystyle \left<\Psi|V_{12}|\Psi \right> = (N!)^{-1} \sum_{i}^{N!} \sum_{j}^{N!} (-1)^{p_i+p_j+2} \left< P_i(\phi_H^*) |V_{12}| P_j(\phi_H) \right>$

Since the Coulomb operator $V_{12}$ only act on particle 1 and 2, thus, the wave function for particle 3, 4,… N has to be the same. Thus, any none-zero product is exchange between particle 1 and 2. therefore, for any given permutation $P_i$ the other permutation must be $P_j = P_i - P_{12} P_i$

$\displaystyle W = (2(N-2)!)^{-1} \sum_{i}^{N!} \left< P_i(\phi_H^*) |V_{12}| P_i(\phi_H) \right> - \left< P_i(\phi_H^*) |V_{12}|P_{12}P_i(\phi_H) \right>$

In the $N!$ permutation, for a fixed orbital of particle 1 and 2, there are $(N-2)!$ way to permute the $N-2$ particle. Thus

$\displaystyle W = \frac{1}{2} \sum_{a}^{N}\sum_{b\neq a}^{N-1} \left< ab|| ab \right> - \left$

or

$\displaystyle W = \sum_{a}^{N}\sum_{b>a}^{N} \left< ab|| ab \right> - \left$

Here we use

$\displaystyle \left<\phi_a^*(1)\phi_b^*(2)|V_{12}|\phi_h(1)\phi_k(2) \right> = \left$

for $N=2$

$\displaystyle W = \sum_{a=1}^{2}\sum_{b=2}^{2} \left< ab|| ab \right> - \left = \left< 12|| 12 \right> - \left<12||21\right>$

which is the result we expect.

## Coulomb Displacement Energy

I went back to my home town for visa renewal and read few books on history ( american revolution, early islam, modern chinese). After I got the visa a back to US, a lot of work has to be done and waiting for me….

The Coulomb displacement energy is the energy difference between isobaric analog states with same isospin. For example, the ground state energy different between 15O and 15N. The only difference between 15O and 15N is the proton in 0p1/2 shell replaced by a 0p1/2 neutron. Both ground state is in isospin $T=1/2$. The Coulomb displacement energy is defined as,

$\displaystyle \Delta = BE(A, Z-1) - BE(A,Z) \\ ~~~~ = M(A,Z) - M(A,Z-1) + (m_n - m_p)$

The binding energies are

$BE(^{15}\textrm{O}) = 111.95535 ~\textrm{MeV}$
$BE(^{15}\textrm{N}) = 115.4919 ~\textrm{MeV}$

The difference of the binding energy is 3.54 MeV.

Let’s calculate the Coulomb energy based on proton charge distributed uniformly in a sphere with radiu $R=1.25 A^{1/3}$. For $A = 15, R = 3.083 \textrm{fm}$. The Coulomb energy for uniform sphere is

$\displaystyle E_c = \frac{3}{5}\frac{e^2}{R} Z(Z-1), e^2 = 1.44 ~\textrm{MeV} \cdot \textrm{fm}$

Thus the Coulomb energies of 15O and 15N are

$E_c(^{15}\textrm{O}) = 15.69 \textrm{MeV}$
$E_c(^{15}\textrm{N}) = 9.81 \textrm{MeV}$

The difference is 4.11 MeV, about 0.5 MeV difference than the binding energy difference!

## Simple model for 4He and NN-interaction

Starting from deuteron, the binding energy, or the p-n interaction is 2.2 MeV.

From triton, 3H, the total binding energy is 8.5 MeV, in which, there are only 3 interactions, two p-n and one n-n. Assume the p-n interaction does not change, the n-n interaction is 4.1 MeV.

The total binding energy of 3He is 7.7 MeV. The p-p interaction is 3.3 MeV.

Notices that we neglected the 3-body force in 3H and 3He. And it is strange that the n-n and p-p interaction is stronger then p-n interaction.

In 4He, the total binding energy becomes 28.3 MeV. I try to decompose the energy in term of 2-body, 3-body, and 4-body interaction.

If we only assume 2-body interaction, the interaction strength from n-p, n-n, and p-p are insufficient. One way to look is the 1-particle separation energy.

The neutron separation energy is 20.6 MeV = 2(p-n) + (n-n).
The proton separation energy is 19.8 MeV = 2(p-n) + (p-p).
The total energy is 28.3 MeV = 4(p-n) + (n-n) + (p-p).

There is no solution for above 3 equations. Thus, only consider 2-body interaction is not enough.

The neutron separation energy is 20.6 MeV = 2(p-n) + (n-n) + 2(n-n-p) + (n-p-p)
The proton separation energy is 19.8 MeV = 2(p-n) +(p-p) + 2(n-p-p) + (n-n-p).
The total energy is 28.3 MeV = 4(p-n) + (n-n) + (p-p) + 2(n-n-p) + 2(n-p-p) + (n-n-p-p).

Assuming the 2-body terms are the same in 2H, 3H, and 3He, the (n-p-p) and (n-n-p) is 4.03 MeV, which is strange again, as the Coulomb repulsion should make the (n-p-p) interaction smaller then the (n-n-p) interaction. The (n-n-p-p) interaction is -4.03 MeV.

Lets also add 3-body force in 3H and 3He.

The neutron separation energy of 3H is 6.3 MeV = (p-n) + (n-n) + (n-n-p)
The toal energy of 3H is 8.5 MeV = 2(p-n) + (n-n) + (n-n-p)
The toal energy of 3He is 7.7 MeV = 2(p-n) + (p-p) + (n-p-p)
The neutron separation energy of 4He is 20.6 MeV = 2(p-n) + (n-n) + 2(n-n-p) + (n-p-p)
The proton separation energy is 4He 19.8 MeV = 2(p-n) +(p-p) + 2(n-p-p) + (n-n-p).
The total energy of 4He is 28.3 MeV = 4(p-n) + (n-n) + (p-p) + 2(n-n-p) + 2(n-p-p) + (n-n-p-p).

We have 6 equations, with 6 unknown [ (p-n) , (n-n) , (p-p) , (n-n-p), (n-p-p), and (n-n-p-p)]. Notice that the equation from proton separation energy of 3He is automatically satisfied. The solution is

(p-n) = 2.2 MeV
(p-p) = – 4.7 – (n-n) MeV
(n-n-p)  = 4.1 – (n-n) MeV
(n-p-p) = 8 + (n-n) MeV
(n-n-p-p) = 0 MeV

It is interesting that there is redundant equation. But still, the (p-n) interaction is 2.2 MeV, and 4-body (n-n-p-p) becomes 0 MeV. Also the (n-p-p) is more bound than (n-n-p) by 3.9 + 2(n-n) MeV. If (n-p-p) should be more unbound, than (n-n) must be negative and smaller than -1.95 MeV.

Since the interaction strength has to be on the s-orbit (mainly), by considering 2H, 3H, 3He, and 4He exhausted all possible equations (I think). We need other way to anchor either (n-n), (p-p), (n-p-p), and (n-n-p) interactions.

Use the Coulomb interaction, the Coulomb interaction should add -1.44 MeV on the NN pair (assuming the separation is a 1 fm). Lets assume the (n-n) – (p-p) = 1.44 MeV

(n-n) = -1.63 MeV
(p-p) = – 3.07 MeV
(n-n-p)  = 5.73 MeV
(n-p-p) = 6.37  MeV

The (p-p) is more unbound than (n-n) as expected, but the (n-p-p) is more bound than (n-n-p) by 0.64 MeV. This is surprising! We can also see that, the 3-body interaction play an important role in nuclear interaction.

According to this analysis, the main contribution of the binding energies of 3H and 3He are the 3-body force.

In 3H:  (n-n) + 2(n-p) + (n-n-p) = -1.6 + 4.4 + 5.7 = 8.5 MeV
In 3He: (p-p) + 2(n-p) + (n-p-p) = -3.1 + 4.4 + 6.4 = 7.7 MeV

Worked on the algebra, when ever the difference  (n-n) – (p-p)  > 0.8 MeV, the (n-p-p) will be more bound that (n-n-p). Thus, the average protons separation should be more than 1.8 fm. I plot the interactions energies with the change of Coulomb energy below.

The (n-n) and (p-p) are isoscalar pair, where tensor force is zero. While the (n-p) quasi-deuteron is isovector pair. Thus, the difference between (n-n) and (n-p) reflect the tensor force in s-orbit, which is 3.8 MeV. In s-orbit, there is no spin-orbital interaction, therefore, we can regard the tensor force is 3.8 MeV for (n-p) isovector pair.

Following this method, may be, we can explore the NN interaction in more complex system, say the p-shell nuclei. need an automatic method. I wonder the above analysis agreed present interaction theory or not. If not, why?

## 3D Harmonic oscillator

Set $x = r/\alpha$The Schrodinger equation is

$\displaystyle \left(-\frac{\hbar^2}{2m} \nabla^2 + \frac{1}{2} m \omega^2 r^2 \right)\Psi = E \Psi$

in Cartesian coordinate, it is,

$\displaystyle -\frac{\hbar^2}{2m}\left( \frac{d^2}{dx^2}+\frac{d^2}{dy^2}+\frac{d^2}{dz^2} \right) \Psi + \frac{1}{2} m \omega^2 (x^2+y^2+z^2) \Psi = E \Psi$

We can set the wave function to be $\Psi(r, \Omega) = X(x) Y(y) Z(z)$

$\displaystyle \left( -\frac{\hbar^2}{2m}\frac{d^2X}{dx^2} + \frac{1}{2} m \omega^2 x^2 X \right) YZ + \\ \left( -\frac{\hbar^2}{2m}\frac{d^2Y}{dy^2} + \frac{1}{2} m \omega^2 y^2 Y \right) XZ + \\ \left( -\frac{\hbar^2}{2m}\frac{d^2Z}{dz^2} + \frac{1}{2} m \omega^2 z^2 Z \right) XY = E XYZ$

we can see, there are three repeated terms, we can set

$\displaystyle -\frac{\hbar^2}{2m}\frac{d^2X}{dx^2} + \frac{1}{2} m \omega^2 x^2 X = E_x X$

We decoupled the X, Y, Z. Each equation is a quadratic equation with energy

$\displaystyle E_x, E_y, E_z = \left(\frac{1}{2} + n \right) \hbar \omega$

and

$\displaystyle E_x + E_y + E_z = \left(\frac{3}{2} + n_x + n_y + n_z \right) \hbar \omega = \left(\frac{3}{2} + n \right) \hbar \omega = E$

The number of states for each energy level is

$\displaystyle C^{n_x+n_y+n_z+2}_2 = C^{n+2}_2 = \frac{(n+2)!}{n!2!}$

The first few numbers of states are 1, 3, 6, 10, 15, 21, 28, … The accumulated numbers of states are 1, 4, 10, 20, 35, 56, 84, … Due to the spin-state, the accumulated numbers of particles are 2, 8, 20, 40, 70, 112, 168, … The few magic numbers are reproduced.

The wave function is the product of the Hermite functions $H_n(x)$ and exponential function

$\Phi(x,y,z) = N H_{n_x} (x) H_{n_y}(y) H_{n_z}(z) \exp(-r^2/2)$

If we simply replace $(x,y,z) \rightarrow r( \cos(\phi) \sin(\theta), \sin(\phi) \sin(\theta) , cos(\theta) )$, we can see the ground state consists of s-orbit, the 1st excited state consists of p-orbit, and the 2nd excited state consists of d-orbit.

To see more clearly, we can project the function onto spherical harmonic, which is fixed angular momentum, i.e.

$\langle \Phi(x,y,z) | Y_{lm}(\theta,\phi) \rangle = C_{lm} R_{nl}(r)$.

where $C_{lm}$ is coefficients and $R_{nl}(r)$ is radial function.

To have a better understanding, the radial function has to be solved. The procedure is very similar to solving Coulomb potential.

$\displaystyle \left(-\frac{\hbar^2}{2m}\left(\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d}{dr}\right)\right) + \frac{\hbar^2}{2m} \frac{L^2}{r^2} + \frac{1}{2} m \omega^2 r^2 \right)\Psi = E \Psi$

separate the radial part and angular part.

$\displaystyle \left( \frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d}{dr}\right) - \frac{l(l+1)}{r^2} - \frac{m^2 \omega^2}{\hbar^2} r^2\right) R = - (2n+3) \frac{m \omega}{\hbar} R$

Set $\alpha^2 = \frac{\hbar}{m \omega}$

$\displaystyle \left( \frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d}{dr}\right) + \frac{2n+3}{\alpha^2} - \frac{l(l+1)}{r^2} - \frac{r^2}{\alpha^4}\right) R = 0$

Set $x = r/\alpha$

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

Set $u(x) = x R(x)$

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

as usual, the short range behaviour is $r^{l+1}$, long range behaviour is $\exp(-x^2/2)$, as stated in the Cartesian coordinate. Thus, we set

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

$\displaystyle x\frac{d^2f}{dx^2} + (2(l+1) -2x^2) \frac{df}{dx} + 2x(n-l) f(x) = 0$

with change of variable $y = x^2$, the equation becomes

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

This is our friend, the Laguerre polynomial! In the Laguerre polynomial, $(n-l)/2 \geq 0$ must be non-negative integer. Now we set $k = (n-l)/2$, than the energy is

$\displaystyle E = \hbar \omega \left( 2k + l + \frac{3}{2} \right)$

In order to have $n, l, k$ are integer, when

$n = 0 \rightarrow l = 0 \\ n = 1 \rightarrow l = 1 \\ n= 2 \rightarrow l = 0, 2 \\ n = 3 \rightarrow l = 1, 3 \\ n = 4 \rightarrow l = 0,2,4$

The overall solution without a normalization factor is

$\displaystyle \Psi_{nlm}(r, \theta, \phi) = r^l \exp\left(-\frac{r^2}{2\alpha^2}\right) L_{k}^{l+\frac{1}{2}}\left( \frac{r^2}{\alpha^2} \right) Y_{lm}(\theta, \phi)$

The normalization constant can be calculate easily using the integral formula of Laguerre polynomial.

$\displaystyle \int R^2(r) r^2 dr = 1$

$\displaystyle N^2 \int r^{2l} \exp\left(-\frac{r^2}{2\alpha^2}\right) \left( L_{k}^{l+\frac{1}{2}}\left( \frac{r^2}{\alpha^2} \right)\right)^2 r^2 dr = 1$

change of variable $y = r^2/\alpha^2 \rightarrow r=\alpha \sqrt{y}$

$\displaystyle dr = \frac{\alpha}{2\sqrt{y}}dy$

$\displaystyle N^2 \frac{\alpha^{2l+3}}{2} \int y^{l+\frac{1}{2}} \exp(-y) \left( L_k^{l+\frac{1}{2}} \right)^2 dy = 1$

The integration is

$\displaystyle N^2 \frac{\alpha^{2l+3}}{2} \frac{(k+l+\frac{1}{2})!}{k!} = 1$

$\displaystyle N^2 = \frac{2}{\alpha^{2l+3}} \frac{k!}{(k+l+\frac{1}{2})!}$

we can use

$\displaystyle \left(n+\frac{1}{2}\right) = \frac{(2n+1)!}{2^{2n+1}n!}\sqrt{\pi}$

Thus,

$\displaystyle N^2 = \frac{1}{\sqrt{\pi}\alpha^{2l+3}} \frac{k! (k+l)! 2^{2k+2l+3}}{(2k+2l+1)!}$

replace $k = (n-l)/2$. The total wave function is

$\displaystyle \Psi_{nlm}(r, \theta, \phi) \\ =\sqrt{ \frac{1}{\sqrt{\pi}\alpha^{2l+3}} \frac{(\frac{n-l}{2})! (\frac{n+l}{2})! 2^{n+l+2}}{(n+l+1)!}} r^l \exp\left(-\frac{r^2}{2\alpha^2}\right) L_{k}^{l+\frac{1}{2}}\left( \frac{r^2}{\alpha^2} \right) Y_{lm}(\theta, \phi)$

Here are some drawing of the square of the wave functions. From the below is $n = 0, 1, 2, 3$, from left to right, are s-orbit, p-orbit, d-orbit, f-orbit.

With the LS coupling, the spatial function does not affected, unless the coupling has spatial dependence. With the LS coupling, the good quantum numbers are $n, l, j, m_j$

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

This is for atomic case, that the eigen energy is $E_n = - \frac{Z^2}{2n^2}$, and the Coulomb force is attractive.

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 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)$.