Skyrme Hartree-Fock of oxygen-16

A nucleus is up to a few hundred protons and neutrons (sixteen here) held together by the strong force, which has no simple closed form. Nuclear density functional theory replaces it with an energy that depends only on a few local densities: the number density ${\rho }_q$ of neutrons ($q=n$) and protons ($q=p$), their kinetic-energy densities ${\tau }_q$, and their spin-orbit densities ${\mathbf{J}}_q$. Skyrme's functional is the oldest of these, with about ten parameters fitted to nuclear data. Minimizing the energy gives each nucleon a Schrödinger equation (the second equation above) with a position-dependent effective mass $m_q^{\ast }$, a central potential $U_q$ and a spin-orbit field ${\mathbf{W}}_q$, all computed from the densities. That makes the problem self-consistent, and the program solves it by iteration: densities, then fields, then orbitals, then new densities.

${}^{16}$O has 8 protons and 8 neutrons, which fill the $1s_{1/2}$, $1p_{3/2}$ and $1p_{1/2}$ shells of each, so it is spherical and every orbital is the lowest of its $(l,j)$ block. The force is SLy5 (Chabanat et al., Nucl. Phys. A 643, 441 (1998)), with the parameters and conventions of skyrme_rpa (G. Colò, L. Cao, N. Van Giai and L. Capelli, Comput. Phys. Commun. 184, 142 (2013)), whose distribution includes a Hartree-Fock run of ${}^{16}$O with SLy5 as its test case. That output file is the published reference here. posts/refs/skyrme-hartree-fock-of-oxygen-16.py copies its numbers, describes rerunning the program, and holds a numpy twin of this solver that diagonalizes each block densely.

The discretization

Each orbital is $\varphi =\frac{u(r)}{r}[Y_l\otimes \chi {]}_j$, and on a mesh $r_i=ih$, $i=1\dots N$, the radial equation for $u$ is a symmetric tridiagonal matrix:

The first term uses $B$ at the half points, which keeps the matrix symmetric; the $B^{\prime }/r$ term comes from writing the equation for $u=r\varphi$ rather than $\varphi$. Each block's lowest eigenvector comes from inverse iteration: three solves of $(H-\sigma )w=u$ by the tridiagonal (Thomas) algorithm per pass, starting from the previous pass's vector, with the shift $\sigma$ 0.1 MeV below the previous orbital energy. On the first pass the eigenvalue itself is found first, by bisection on the Sturm count (the number of negative pivots of $H-x$ is the number of eigenvalues below $x$), so the iteration starts on the right state. Derivatives of the densities are centered differences with second-order one-sided ones at the ends, and integrals are trapezoidal. The protons also feel the Coulomb field of their own density, with exchange in the Slater approximation.

Every piece is second order in $h$, so the program solves the problem at $h=0.1$ and $h=0.05$ fm in a 15 fm box and extrapolates, $x\approx (4x_{h/2}-x_h)/3$.

Two meshes and the extrapolation

The iteration starts from a Woods-Saxon density and mixes half of each new density into the old one until the energy changes by less than $10^{-9}$ MeV, which takes 38 passes on each mesh. On both meshes the energy per nucleon, the three radii and the six orbital energies match the numpy twin to $7\times 10^{-9}$; the tests allow $10^{-7}$. Going from $h=0.1$ to $h=0.05$ fm moves the energy per nucleon by 0.016 MeV and the charge radius by 0.002 fm.

Against the published run

The extrapolated binding energy is $8.0250$ MeV per nucleon, against skyrme_rpa's $8.0261$: the two differ by $1.1\times 10^{-3}$ MeV per nucleon, 18 keV out of 128 MeV. skyrme_rpa's own answer moves by $1.2\times 10^{-4}$ MeV per nucleon between its 0.1 fm and 0.05 fm meshes (the reference script records the rerun), so most of the remaining difference is between the two programs' formulas, not their meshes. The radii agree to $2\times 10^{-4}$ fm and the six orbital energies to 3 keV. On the 0.1 fm mesh alone, this program's second-order scheme is off by 0.02 MeV per nucleon, which is why it extrapolates.

The figure shows the converged densities on the fine mesh.