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

The boundary conditions are continuity and differential continuity.

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

$\displaystyle \left(\frac{d}{dr}j_l(kr)\right)_{r=a} = A \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,

## Hartree-Fock for Helium excited state II

This times, we will show the energy level diagram for helium atom. We already calculated the ground state, the 1s2s singlet and triplet excited states. We will calculate higher excited states, for example, 1s3s, 1s2p or 1s3p, etc, and included in the diagram.

The most difficult part is the evaluation of the mutual interaction matrix element

$\displaystyle G_{ij}^{hk} = \langle b_i(x) b_h(y) | \frac{1}{r_{xy}} | b_j(x) b_k(y) \rangle$

The angular integral was evaluated in the last post. And the radial integral is done using Mathematica.

I will use hydrogen 1s, 2s, 2p, 3s, 3p, 4s, and 4p states for basis. Thus, the diagram will contain some of the excited states from some possible combination. The 1s4s and 1s4p state will not be calculated, because there is no room for mixing with higher excited states but only for lower excited states. This will make the eigen state be unbound.

During the calculation, one of the tricky point is the identify of the n-th s-state. It is because the mixing is always there, and the mixing can be very large. In order to determine the principle quantum number, we have to check the wave function can see how many peaks in there.

Here is the energy level diagram obtained by Hartree-Fock method using limited hydrogen wave functions as basis set.

## Evaluation of mutual interaction

We are going to evaluation the integral

$\displaystyle G_{ij}^{hk} = \langle b_i(x) b_h(y) | \frac{1}{r_{xy}} | b_j(x) b_k(y) \rangle$

Recalling the multi-pole expansion,

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

and the basis

$b_{nlm}(\vec{r}) = R_{nl}(r) Y_{lm}(\Omega)$

Set an averaged basis

$\displaystyle b_{nl}(\vec{r}) = R_{nl}(r) \frac{1}{\sqrt{2l+1}} \sum_{m=-1}^{l}Y_{lm}(\Omega)$

$\displaystyle \Gamma_{ij}^{hk}(l) = \frac{4 \pi}{2l+1} \sum_{m=-l}^{l} \frac{1}{\sqrt{(2l_i+1)(2l_j+1)(2l_h+1)(2l_k+1)}} \\ \sum_{m_i,m_j}\int Y_{l_i m_i}^{*}(\Omega_1) Y_{l m}^{*}(\Omega_1) Y_{l_j m_j}(\Omega_1) d \Omega_1 \\ \sum_{m_h,m_k} \int Y_{l_h m_k}^{*}(\Omega_2) Y_{l m}(\Omega_2) Y_{l_k m_k}(\Omega_2) d \Omega_2$

In the angular integrals, using Wigner 3-j symbol and the integral

$\displaystyle\int Y_{l_1 m_1}(\Omega) Y_{l m}(\Omega) Y_{l_2 m_2}(\Omega) d \Omega \\= \frac{\sqrt{(2l_1+1)(2l+1)(2l_2+1)}}{\sqrt{4\pi}} \begin{pmatrix} l_1 & l & l_2 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_1 & l & l_2 \\ m_1 & m & m_2 \end{pmatrix}$

$Y_{lm}^{*}(\Omega) = (-1)^m Y_{l(-m)}(\Omega)$

$\displaystyle \begin{pmatrix} l_1 & l & l_2 \\ m_1 & m & m_2 \end{pmatrix} = (-1)^{l_1+l+l_2} \begin{pmatrix} l_1 & l & l_2 \\ -m_1 & -m & -m_2 \end{pmatrix} \\ =(-1)^{l_1+l+l_2} \begin{pmatrix} l_1 & l_2 & l \\ m_1 & m_2 & m \end{pmatrix}$

we have

$\displaystyle\int Y_{l_i m_i}^{*}(\Omega) Y_{l m}^{*}(\Omega) Y_{l_j m_j}(\Omega) d \Omega \\= (-1)^{m_j} \frac{\sqrt{(2l_i+1)(2l+1)(2l_j+1)}}{\sqrt{4\pi}} \begin{pmatrix} l_i & l & l_j \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_i & l & l_j \\ m_i & m & -m_j \end{pmatrix}$

$\displaystyle\int Y_{l_h m_h}^{*}(\Omega) Y_{l m}(\Omega) Y_{l_k m_k}(\Omega) d \Omega \\= (-1)^{m_h} \frac{\sqrt{(2l_h+1)(2l+1)(2l_k+1)}}{\sqrt{4\pi}} \begin{pmatrix} l_h & l & l_k \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_h & l & l_k \\ -m_h & m & m_k \end{pmatrix}$

Thus,

$\displaystyle \Gamma_{ij}^{hk}(l) = \sum_{m=-l}^{l} \sum_{m_i,m_j } (-1)^{m_j} \begin{pmatrix} l_i & l_j & l \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_i & l_j & l \\ m_i & -m_j & m \end{pmatrix} \\ \sum_{m_h,m_k} (-1)^{m_h} \begin{pmatrix} l_h & l_k & l \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_h & l_k & l \\ -m_h & m_k & m \end{pmatrix}$

The 3-j symbol restricted that

$l_i+l_j+l = even$

$l_h+l_k+l = even$

$m_i-m_j + m = 0, |m_i| \leq l_i , |m_j| \leq l_j$

$-m_h + m_k +m = 0, |m_h| \leq l_h , |m_k| \leq l_k$

I guess it is the most simplified formula for the angular part.

The total integral is

$\displaystyle G_{ij}^{hk} = \sum_{l=0}^\infty \Gamma_{ij}^{hk}(l) \langle R_i(x) R_h(y) | \frac{r_<^l}{r_>^{l+1}} |R_j(x) R_k(y) \rangle$

The angular integral imposes condition for $l$.

I am not sure this is a correct way to treat the problem.

First, the averaged basis is still an energy eigen state. It is not the eigen state for the angular part. So, this averaging could introduces an error and we should reminded that this is an approximation. But in the perturbation view point, this averaged basis is still valid.

Second thing is, the sum

$\displaystyle P_{(l_i m_i) (l_j m_j)}^{(l_h m_h)(l_k m_k)}(l) = \sum_{m=-l}^{l} \langle Y_{l_i m_i}|Y_{lm}^*|Y_{l_j m_j}\rangle \langle Y_{l_h m_h}|Y_{lm}|Y_{l_k m_k}\rangle$

is not symmetry for exchange of $i, j, h,k$ in general. For example,

$\displaystyle \frac{-1}{3}=P_{(1-1) (00)}^{(11)(00)}(1) {\neq} P_{(00)(1-1)}^{(11)(00)}(1) = 0$

This is a very uprising result that the mutual interaction dependent on the magnetic quantum number. Thus, in detail, we should use $n l m m_s$ as a basis.

Third, the sum $P_{ij}^{hk}(l)$ is depend on $l$. The mutual interaction require us to sum all possible $l$.

Fourth, the coupling between 1s2p triplet state, the total spin is $S = 1$, total L is $L = 1$, and the total angular momentum can be  $J = 0, 1, 2$. In our treatment, we did not coupled the angular momentum in the calculation explicitly. In fact, in the integral of the spherical harmonic, the coordinate are integrated separately, and the coupling seem to be calculated implicitly. I am not sure how to couple two spherical harmonics with two coordinates.

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

## Using generating function to calculate integrals of Laguerre polynomial

The generating function of Laguerre polynomial is

$\displaystyle \sum_n^\infty L_n^\alpha(x) t^n = \frac{1}{(1-t)^{\alpha+1}} \exp \left(-\frac{tx}{1-t}\right) =G_\alpha(t, x)$

We are going to evaluate the integral

$\displaystyle \int_{0}^{\infty} x^{\alpha +k} e^{-x} L_m^{\beta}(x) L_n^{\alpha}(x) dx$

Consider this integral

$\displaystyle \int_0^\infty x^{\alpha+k} e^{-x} G_\alpha(t,x) G_{\beta}(s,x) dx$

It is equal to

$\displaystyle \int_0^\infty x^{\alpha+k} e^{-x} \sum_{n,m}^{\infty} L_n^\alpha(x) L_m^\beta(x) t^n s^m dx \\ = \sum_{n,m}^{\infty} t^n s^m \int_0^\infty x^{\alpha+k} e^{-x} L_n^\alpha(x) L_m^\beta(x) dx \\ = \frac{1}{(1-t)^{\alpha+1}} \frac{1}{(1-s)^{\beta+1}} \int_0^\infty x^{\alpha+k} e^{-x} e^{-\frac{t x}{1-t} - \frac{s x}{1-s}} dx$

It is easy to calculate the last integral,

$\displaystyle \int_0^\infty x^{\alpha+k} e^{-x} e^{-\frac{t x}{1-t} - \frac{s x}{1-s}} dx = \int_0^\infty x^{\alpha+k} \exp \left(-\frac{1- t s}{(1-t)(1-s)}x \right) dx$

Replace variable,

$\displaystyle \int_0^\infty x^{\alpha+k} e^{-x} e^{-\frac{t x}{1-t} - \frac{s x}{1-s}} dx \\ = \frac{(1-t)^{(\alpha+1+k)}(1-s)^{(\alpha+1+k)}}{(1-t s)^{(\alpha+1+k)}}\int_0^\infty y^{\alpha+k} e^{-y} dy$

The integral becomes

$\displaystyle \sum_{n,m}^{\infty} t^n s^m \int_0^\infty x^{\alpha+k} e^{-x} L_n^\alpha(x) L_m^\beta(x) dx = \frac{(1-t)^{k}(1-s)^{(\alpha-\beta+k)}}{(1-t s)^{(\alpha+1+k)}} (\alpha+k)!$

Expand the right hand side,

$\displaystyle (1-t)^k = \sum_i^k (-1)^i \begin{pmatrix} k \\ i \end{pmatrix} t^i$

$\displaystyle (1-s)^{(\alpha-\beta+k)} = \sum_i^{\alpha-\beta+k} (-1)^i \begin{pmatrix} \alpha-\beta+k \\ i \end{pmatrix} s^i$

$\displaystyle (1-ts)^{-(\alpha+1+k)} = \sum_i^\infty (-1)^i \begin{pmatrix} -(\alpha+k+1) \\ i \end{pmatrix} (t s)^i = \sum_i^\infty \begin{pmatrix} \alpha+k+i\\ i \end{pmatrix} (t s)^i$

We used

$\displaystyle \begin{pmatrix} m \\ n \end{pmatrix} = \prod_{i=0}^{n-1} \frac{m-i}{n-i}$

$\displaystyle (-1)^n\begin{pmatrix} -m \\ n \end{pmatrix} = \prod_{i=0}^{n-1} \frac{m+i}{n-i} = \frac{(m+n-1)!}{n!(m-1)!}$

At the last step. The product

$\displaystyle \frac{(1-t)^{k}(1-s)^{(\alpha-\beta+k)}}{(1-t s)^{(\alpha+1+k)}} \\ = \sum_a^k \sum_{b}^{\alpha-\beta+k} \sum_c^\infty (-1)^{(a+b)} \begin{pmatrix} k \\ a \end{pmatrix} \begin{pmatrix} \alpha-\beta+k \\ b \end{pmatrix} \begin{pmatrix}\alpha+k+c \\ c \end{pmatrix} t^{(a+c)} s^{(b+c)}$

By comparing the order of $t^n s^m$, the integral can be evaluated. With $a+c = n, b+c = m$

$\displaystyle \int_0^\infty x^{\alpha+k} e^{-x} L_n^\alpha(x) L_m^\beta(x) dx \\ = (\alpha+k)! \sum_a^k \sum_{b}^{\alpha-\beta+k} \sum_c^\infty (-1)^{(a+b)} \begin{pmatrix} k \\ a \end{pmatrix} \begin{pmatrix} \alpha-\beta+k \\ b \end{pmatrix} \begin{pmatrix} \alpha+k+c \\ c \end{pmatrix}$

When $\alpha = \beta$

$\displaystyle \int_0^\infty x^{\alpha+k} e^{-x} L_n^\alpha(x) L_m^\alpha(x) dx \\ = (\alpha+k)! \sum_a^k \sum_{b}^{k} \sum_c^\infty (-1)^{(a+b)} \begin{pmatrix} k \\ a \end{pmatrix} \begin{pmatrix} k \\ b \end{pmatrix} \begin{pmatrix} \alpha+k+c \\ c \end{pmatrix}$

For example, $k = 0$,

$\displaystyle \frac{(1-s)^{(\alpha-\beta)}}{(1-t s)^{(\alpha+1)}} = \sum_{b}^{\alpha-\beta} \sum_c^\infty (-1)^{b} \begin{pmatrix} \alpha-\beta \\ b \end{pmatrix} \begin{pmatrix} \alpha+c \\ c \end{pmatrix} t^{c} s^{(b+c)}$

When $\alpha \neq \beta, , n= m = 1$, only $b = 0, c = 1$

$\displaystyle \int_0^\infty x^{\alpha} e^{-x} L_1^\alpha(x) L_1^\beta(x) dx = (\alpha)! \begin{pmatrix} \alpha-\beta \\ 0 \end{pmatrix} \begin{pmatrix} \alpha+1 \\ 1 \end{pmatrix} = (\alpha+1)!$

When $\alpha = \beta, , n = m$, only $b = 0, c = n$

$\displaystyle \int_0^\infty x^{\alpha} e^{-x} (L_n^\alpha(x))^2 dx = (\alpha)! \begin{pmatrix} \alpha+n \\ n \end{pmatrix} = \frac{(\alpha+n)!}{n!}$

For $k = 1, \alpha = \beta, n= m$, there are two terms, $a = b = 0, c = n$  and $a = b = 1, c = n -1$

$\displaystyle \int_0^\infty x^{\alpha+1} e^{-x} (L_n^\alpha(x))^2 \\ = (\alpha+1)! \left(\begin{pmatrix} \alpha+n+1 \\ n \end{pmatrix} + \begin{pmatrix} \alpha+n \\ n-1 \end{pmatrix} \right) = \frac{(\alpha+n)!}{n!}(2n+\alpha+1)$

We recovered all formula. Using generating function, we don’t have to worry about the condition of $n, m, \alpha, \beta, k$. We only need to fulfill the sum with condition $a +c = n, b+c = m$ and $a, b, c \geq 0$.

## Qt Script Engine

In order to use “any” function for fitting algorithm, the Qt script engine is a very powerful tool. And the basic usage is very simple.

QScriptEngine engine;
QString function = "p1*x + pow(x,2)";
double p1 = 3.2;
double x = 5;
engine.globalObject().setProperty("p1", p1);
engine.globalObject().setProperty("x", x);
ans = engine.evaluate(function_str).toNumber();

The above code will evaluate the expression $p_1 x + x^2$. Although the running time will be a bit slower then hard code, the advantage is that, the expression can be read from a file.