A star of free neutrons can weigh at most 0.71 suns

The reference, in posts/refs/the-free-neutron-star-limit.py, integrates the same physics a different way: in cgs units, with the radius as the independent variable and the Fermi momentum as the unknown, by scipy's DOP853 at a relative tolerance of $10^{-12}$, and finds the maximum with minimize_scalar. It gives $M_{\text{max}}=0.71018029M_{\odot }$ at $R=9.161465$ km. The literature value is from G. Narain, J. Schaffner-Bielich and I. N. Mishustin, "Compact stars made of fermionic dark matter", Phys. Rev. D 74, 063003 (2006), arXiv:astro-ph/0605724: from their Eqs. (31) and (32), a star of free neutrons has "the maximum mass $M_{\text{max}}=0.710M_{\odot }$ at a radius of $R_{\text{min}}=9.192$ km, which match very well the original results of Oppenheimer and Volkoff" (Phys. Rev. 55, 374, 1939).

The program takes 400 steps per star, and its maximum agrees with scipy's to $2\times 10^{-8}$ of a solar mass; the tests allow $10^{-6}$, relative.

The question

A neutron star is held up by the pressure of its neutrons. Neutrons are fermions: no two can share a state, so packing them tighter forces them into states of higher momentum, and that momentum pushes outward even at zero temperature. In 1939 Oppenheimer and Volkoff asked how heavy a star held up this way alone can be, ignoring the nuclear force between the neutrons. In general relativity pressure is itself a source of gravity, so a star that squeezes its center harder to push back also pulls harder, and past some central density the extra pull wins.

The equation of state

A cold gas of free neutrons is described by one number, the Fermi momentum $p_F$, written as $x=p_F/m_{n}c$. The energy density $\epsilon$ (including the neutrons' rest mass) and the pressure $P$ are

with ${\epsilon }_0=m_n^4{c}^5/8{\pi }^2{\hslash }^3$. The bracket in $P$ is a difference of nearly equal terms when $x$ is small, so below $x=0.01$ the program uses the first two terms of its series, $\frac{8}{15}{x}^5(1-\frac{5}{14}{x}^2)$, which agree with the full formula to $2\times 10^{-9}$ there.

The structure equations

In geometric units ($G=c=1$, lengths in km) the Tolman-Oppenheimer-Volkoff equations give the mass $m(r)$ inside radius $r$ and the pressure there. The program writes them with the log-enthalpy $h=\ln (\mu /m_n{c}^2)$ as the independent variable, where $\mu$ is the energy it costs to add one neutron. For free neutrons $\mu =m_n{c}^2\sqrt{1+x^2}$, so $x=\sqrt{e^{2h}-1}$ and the equation of state needs no inverting. The surface, where the pressure is zero, is $h=0$ exactly, so every star is integrated over a known interval, from $h_c$ at the center to 0; the equations are the two at the top of the page (Lindblom, ApJ 398, 569, 1992).

Near the center $r$ grows like $\sqrt{h_c-h}$, and near the surface the density falls like $h^{3/2}$; a fixed-step method follows neither well. So the program steps evenly in a variable $t$ that runs from 0 at the center to 1 at the surface, defined by $\sqrt{h_c-h}=\sqrt{h_c}t(2-t)$. Near the center that makes $r$ a straight line in $t$, $r\approx \sqrt{h_c-h}\sqrt{3/2\pi ({\epsilon }_c+3P_c)}$, and near the surface $h\approx 2h_c(1-t{)}^2$, so $x\propto 1-t$ and the density is a polynomial in $t$ there too. The integration starts at $t=0.01$ from those two leading terms, with $m=\frac{4\pi }{3}{\epsilon }_c{r}^3$, because a fourth-order step across the center itself, where $r$ and $m$ are both zero, is only second-order accurate. With 400 steps a star's radius and mass agree with the scipy twin to about $3\times 10^{-7}$, the digits the twin was recorded to; with 100 steps, to $3\times 10^{-5}$.

The sequence of stars

Each central density gives one star. A light star is large: at $x_c=0.5$ (central Fermi momentum half of $m_{n}c$) it has 0.594 solar masses in a radius of 13.5 km. Before running: squeeze the center harder, to $x_c=2$. Is the star heavier?

It is lighter, 0.476 solar masses in 5.08 km. Somewhere in between the mass passed through a maximum. The next cell draws the whole sequence and then finds the peak by golden-section search on $\ln x_c$.

The heaviest star has $M=0.7102$ solar masses and $R=9.16$ km (m_max, r_max), at a central Fermi momentum $x_c=0.834$ (x_max) and a central density of $4.2\times 10^{15}$ g/cm³ (rho_c). Narain, Schaffner-Bielich and Mishustin (2006) quote 9.192 km for the radius, 0.3% larger than both this program and the scipy twin; their number comes from a scaling formula (their Eq. (32)) whose coefficient, 3.367, they give to four digits.

Along the sequence, the stars to the right of the peak (larger, lighter, less dense) are stable: add a little mass and the star settles slightly smaller. The stars to its left are not, and a small squeeze makes them collapse. So 0.71 solar masses is the most a star of free neutrons can hold. Heavier neutron stars need the force between neutrons, which this model leaves out, to be repulsive at these densities.