Mass–Radius Sequences

An equation of state does not predict one star but a family of equilibria. Varying the central density produces a curve in the mass–radius plane whose shape reveals the low-density envelope, the stiffness of the core, possible phase-transition branches, and the largest static mass. Constructing that curve is a computational problem with unusually transparent physics: every point must satisfy the same TOV boundary-value structure, yet numerical choices near the center, surface, and EOS joins can manufacture apparent branches if they are not controlled.

Required background. The Einstein Equation supplies relativistic field equations; Killing Fields and Isometry Groups supplies the static spherical reduction.

Helpful background. Flat-Background Hilbert Stress-Energy from Metric Variation fixes matter conventions; Relativistic Matter and Compact Objects supplies EOS context; Computational Relativity develops convergence and reproducibility practice.

One central density, one equilibrium

Assume a cold barotropic EOS in the form $p(\epsilon)$, with a stated crust prescription and a causal validity interval. For each central energy density $\epsilon_c$, evaluate $p_c=p(\epsilon_c)$ and integrate

$$\frac{dm}{dr}=4\pi r^2\epsilon, \qquad \frac{dp}{dr} =-\frac{G(\epsilon+p)(m+4\pi r^3p)} {r^2(1-2Gm/r)}.$$

Regular center series provide starting data at $r=r_0\gt0$. The first zero-pressure point defines the circumferential radius $R$, and $M=m(R)$. Repeating this calculation over an ordered grid in $\epsilon_c$ produces $M(\epsilon_c)$ and $R(\epsilon_c)$. A logarithmic central-density coordinate usually resolves both diffuse low-mass stars and compact high-density models more efficiently than a uniform grid.

The curve is parameterized, not merely drawn. Two configurations can have the same $M$ but different $R$, central density, composition, and stability. A sharp first-order phase transition may produce a disconnected stable branch, but an apparent back-bending can also come from inconsistent interpolation. The representative piecewise-polytropic construction and its density joins are specified in Read et al. 2009, §III, Eqs. (1)–(6); any imported parameters must be translated to explicit $G$ and a single energy-density convention.

A worked relativistic polytrope

Take a cold teaching EOS with baryon reference rest-mass density $\rho_0=m_b n$,

$$p=K\rho_0^2,\qquad \epsilon=\rho_0+p,\qquad K\gt0.$$

The internal-energy term is essential: replacing $\epsilon$ by $\rho_0$ in a relativistic calculation changes the model and violates $d\epsilon=(\epsilon+p)\,d\rho_0/\rho_0$. Here the sound speed is $c_s^2=2K\rho_0/(1+2K\rho_0)\lt1$.

Introduce $q=K\rho_0$ and dimensionless quantities

$$\bar r=r\sqrt{\frac{G}{K}},\qquad \bar m=\frac{G^{3/2}m}{\sqrt K},\qquad \bar p=Kp=q^2,\qquad \bar\epsilon=K\epsilon=q+q^2.$$

The barred TOV equations have the displayed form above with $G$ replaced by one. The relativistic enthalpy is particularly simple: $H=\ln(1+2q)$. Thus the numerical solution needs only the central value $q_c$; changing $K$ rescales every mass and radius by $\sqrt K$ at fixed $G$, while leaving $GM/R=\bar M/\bar R$ unchanged. No neutron-star density scale or composition has been fitted.

The downloadable Python calculation integrates each model in two ways: using radius and using enthalpy as the independent variable. It also integrates baryonic mass, $d\bar M_b/d\bar r=4\pi\bar r^2q/\sqrt{1-2\bar m/\bar r}$. Selected results are:

$q_c$	$\bar M$	$\bar R$	$\bar M_b$	$GM/R$
$0.01$	$0.023095$	$1.219832$	$0.023315$	$0.018933$
$0.05$	$0.086159$	$1.108344$	$0.089626$	$0.077737$
$0.10$	$0.126957$	$1.004735$	$0.135313$	$0.126359$
$0.20$	$0.157377$	$0.865793$	$0.171753$	$0.181772$
$0.32$	$0.163727$	$0.762308$	$0.179861$	$0.214778$
$0.50$	$0.158431$	$0.668271$	$0.172651$	$0.237077$
$1.00$	$0.136846$	$0.551113$	$0.142386$	$0.248308$

Mass first rises and then falls, while compactness continues to grow across these samples. Refining the extremum gives $q_c\simeq0.31824$, $\bar M_{\max}\simeq0.163728$, and $\bar R\simeq0.763515$. The largest compactness in the table therefore belongs to a configuration with less mass than the maximum. These are equilibrium results; the script does not calculate radial eigenfrequencies or assign stability labels.

To reproduce the table, save the linked file and run:

python3 verify-stellar-teaching-models.py > stellar-teaching-results.json

It uses Python 3.9 or later and only the standard library. Its JSON records the actual environment, center and surface prescriptions, results at three resolutions, and pass/fail checks. Refinement from 2048 to 4096 steps and agreement between the two coordinate methods give relative differences below $2\times10^{-6}$ for $M$, $R$, and $M_b$ at all listed central densities. This is a convergence check on this smooth model, not an error estimate for a tabulated nuclear EOS.

Two physical checks add information that solver agreement alone would miss. In the dilute Newtonian limit the EOS is an index-one polytrope: $\bar R\to\sqrt{\pi/2}$ and $\bar M/q_c\to\sqrt{2\pi}$. At $q_c=10^{-5}$, the computed relative differences from these limiting values are $2.8\times10^{-5}$ and $8.4\times10^{-5}$, respectively; they decrease as $q_c$ decreases. The stellar first law also gives $dM/dM_b=\sqrt{1-2GM/R}$ for this EOS, because its zero-pressure chemical potential is $m_b$. At $q_c=0.1$, a centered sequence derivative agrees to $4.4\times10^{-7}$ relative accuracy after refinement.

Extending the calculation to a tabulated EOS

A production sequence must specify the EOS release, units, independent thermodynamic variables, crust join, central-density interval, surface condition, and observables. Its reproducible implementation must also fix the solver, error tolerances, interpolation algorithm, event-location rule, and floating-point precision.

A robust implementation evolves pressure or relativistic enthalpy. Enthalpy often handles the surface more cleanly because

$$H(p)=\int_0^p\frac{dp'}{\epsilon(p')+p'}$$

vanishes at $R$. Thermodynamic table interpolation must preserve monotonicity and should reproduce derivatives used for sound speed. Interpolating $p$ and $\epsilon$ independently can violate the first law even when the plotted curve looks smooth.

The central-density grid is adaptively refined near extrema and rapid phase changes. For each model, record the surface event residual, final ODE residual, mass and radius, baryon number if available, maximum compactness, and whether the EOS remained inside its tabulated domain. The sequence terminates rather than extrapolates silently beyond that domain.

Branches, extrema, and what they imply

For ordinary one-parameter nonrotating barotropes, the first maximum of $M(\epsilon_c)$ commonly coincides with a zero-frequency radial mode. That is a stability result, not a plotting convention; its hypotheses are developed in Radial Stability and the Turning-Point Criterion. Additional extrema can arise from phase transitions and require an explicit mode count before being labeled stable or unstable.

The qualitative dependence on the EOS is reviewed in Lattimer and Prakash 2016, §§2.1–2.2 and §4.2: greater pressure over the densities sampled by a star generally increases its radius or maximum supported mass, but no single pressure value controls the full curve. Low-density crust physics affects radii more strongly for low-mass configurations, whereas the high-mass endpoint probes higher-density pressure. The original neutron-core sequence in Oppenheimer and Volkoff 1939, §II, Eqs. (7)–(10), and Fig. 1 is historically important but uses an idealized free-neutron EOS, not a modern mass bound.

Validation, uncertainty, and failure modes

The primary benchmark is an independent integration of a smooth relativistic polytrope, compared at fixed central density. A constant-density solution is a useful sign and surface check even though it is acausal. Convergence requires simultaneously decreasing ODE tolerances, refining the central-density grid, and refining or changing EOS interpolation. Report norms or fractional changes in $M$, $R$, and the location of each extremum, not only visual overlap.

Important failure modes include starting too far from the center, stepping over the surface, interpolating across a discontinuity as though it were smooth, confusing rest-mass and energy density, mixing cgs and geometrized units, and declaring an endpoint because the table ended. Uncertainty from EOS parameters is distinct from numerical truncation error. A credible sequence propagates both and states which one dominates each quoted feature. A particularly deceptive failure mode is a smooth-looking but thermodynamically inconsistent interpolant, because ordinary ODE residuals will not expose it.

Check a reported maximum

A calculation reports $M_{\max}$ to four significant digits. Halving the ODE tolerance changes it by $0.01\%$, but doubling the number of central-density samples shifts it by $1.2\%$. Which error controls the reported maximum, and what should be done?

Solution

The central-density sampling error controls the extremum: the ODE solution at each sampled point is already much more accurate than the discrete location of the maximum. Refine adaptively around the sign change of $dM/d\epsilon_c$ or fit a local convergent interpolant, then repeat with a denser grid until the maximum’s mass and central density stabilize. Four significant digits are not supported while the dominant change is at the percent level.

References

• Lattimer, James M., and Madappa Prakash. “The Equation of State of Hot, Dense Matter and Neutron Stars.” Physics Reports 621 (2016): 127–164. Open preprint.
• Oppenheimer, J. Robert, and George M. Volkoff. “On Massive Neutron Cores.” Physical Review 55 (1939): 374–381. DOI.
• Read, Jocelyn S., Benjamin D. Lackey, Benjamin J. Owen, and John L. Friedman. “Constraints on a Phenomenologically Parametrized Neutron-Star Equation of State.” Physical Review D 79 (2009): 124032. DOI.