Neutron pairing in a neutron-star crust

Below the outer crust of a neutron star, nuclei are so neutron-rich that neutrons drip out of them and fill the space between as a gas. That gas is expected to be a superfluid: neutrons of opposite spin pair up, as electrons do in a superconductor, and the energy needed to break a pair is twice the pairing gap $\Delta$. Deeper down, the nuclei are thought to flatten into sheets (the "lasagna" phase). This post takes one such sheet, a slab of nuclear matter, in one dimension: a cell 30 fm wide with the slab in the middle and neutron gas on either side.

The method is the Hartree-Fock-Bogoliubov one, cut down to its pairing part. A neutron's energy in the cell is $h=-\frac{{\hslash }^2}{2m}\frac{d^2}{d{x}^2}+V(x)$, with $V$ a Woods-Saxon well 40 MeV deep and 10 fm wide standing for the slab. Pairing mixes a neutron state ($u$) with a hole state ($v$), which gives the Bogoliubov-de Gennes (BdG) equations above. The pairing field $\Delta (x)$ depends on the solutions through the anomalous density $\kappa (x)=-{\sum }_{E<0}u(x)v(x)$, so the equations are solved over and over until $\Delta$ stops changing. In one dimension a contact force gives a finite $\kappa$ without a cutoff, so no regularization is needed.

posts/refs/pairing-in-the-neutron-star-crust.py solves the same cell with numpy (the matrix in block form, not interleaved, and 100 times tighter convergence), and solves the uniform gas from the BCS gap equation without diagonalizing anything.

The slab cell

The kinetic energy is the three-point second difference, so each site couples to its neighbors with $\tau ={\hslash }^2/(2m\Delta x^2)$. The cell has hard walls at both ends. The iteration starts from $\Delta =2$ MeV everywhere and stops when no site's $\Delta$ moves by more than $10^{-6}$ MeV in one step. With the mixing, what is left is a few $10^{-6}$ MeV: the gas gap comes out 0.915113 against the numpy solver's 0.915110, which ran to $10^{-8}$. The tests allow $10^{-5}$.

The BCS gap equation in the uniform gas

To check the solver apart from the slab, the program sets $V=0$ and joins the ends of the cell into a ring. Every site is then the same, $\Delta$ is uniform, and the BdG equations separate into plane waves of momentum $k$, each with quasiparticle energy $E_k=\sqrt{{\xi }_k^2+{\Delta }^2}$ and ${\xi }_k=2\tau (1-\cos k\Delta x)-\mu$. The self-consistency condition becomes the BCS gap equation on this lattice,

which the reference script solves with a root finder: $\Delta =0.9610894$ MeV. The ring run, iterated to $10^{-9}$, lands on it to $10^{-8}$.

In the slab cell the gas-side gap, 0.9151 MeV, is 5% below this uniform value. The gas there is 10 fm wide on each side, bounded by a hard wall and by the slab, not an infinite gas. The gap there peaks at 1.012 MeV about 4 fm from the wall, falls to zero at the wall, and falls toward the slab. Inside the slab, where the neutrons are 40 MeV deeper in the well and far denser, it drops to 0.152 MeV. The figure shows the gap and the neutron density across the cell.