General Relativity TOV Solver

General Relativity TOV Solver#

This page describes the standard Tolman-Oppenheimer-Volkoff (TOV) equations in general relativity and their numerical implementation in jester.

Background#

The structure of a spherically symmetric, non-rotating neutron star in hydrostatic equilibrium is governed by the TOV equations, first derived by Tolman [1] and Oppenheimer & Volkoff [2]. In geometric units (\(G = c = 1\)), the equations read

(1)#\[\frac{dp}{dr} = -\frac{(\varepsilon + p)(m + 4\pi r^3 p)}{r(r - 2m)}\]

(2)#\[\frac{dm}{dr} = 4\pi r^2 \varepsilon\]

where \(r\) is the circumferential radius, \(p(r)\) is the pressure, \(\varepsilon(r)\) is the energy density, and \(m(r)\) is the gravitational mass enclosed within radius \(r\). All calculations in jester are carried out in geometric units, and the results are converted to physical units (solar masses and kilometres) only at the final output stage.

Reformulation in terms of the specific enthalpy#

In practice it is numerically advantageous to integrate the TOV equations using the specific enthalpy

(3)#\[h = \int_0^p \frac{dp'}{\varepsilon(p') + p'}\]

as the independent variable rather than the radius \(r\). The enthalpy decreases monotonically from its central value \(h_c\) to zero at the stellar surface, which gives a natural stopping condition for the integration. Inverting the relation \(dh/dr = -(1/(\varepsilon + p))(dp/dr)\) and substituting Eq. (1) yields

(4)#\[\frac{dr}{dh} = -\frac{r(r - 2m)}{m + 4\pi r^3 p}\]

and consequently

(5)#\[\frac{dm}{dh} = 4\pi r^2 \varepsilon \frac{dr}{dh}\]

These are the two equations integrated by the GRTOVSolver in jester. The pressure \(p\) and energy density \(\varepsilon\) at a given enthalpy value are recovered by interpolating the pre-computed EOS table (stored in EOSData).

Tidal deformability#

In addition to the stellar structure, jester computes the tidal deformability \(\Lambda\), which is directly measurable from gravitational wave observations and is defined as

(6)#\[\Lambda = \frac{2}{3} k_2 C^{-5}\]

where \(C = M/R\) is the compactness and \(k_2\) is the second tidal Love number [3, 4]. To compute \(k_2\), a linear perturbation \(H(r)\) to the metric is evolved alongside the stellar structure. Expressed in terms of the enthalpy, the two auxiliary ODEs are

(7)#\[\frac{dH}{dh} = \beta \frac{dr}{dh}\]

(8)#\[\frac{d\beta}{dh} = -\left(C_0 H + C_1 \beta\right) \frac{dr}{dh}\]

where \(\beta \equiv dH/dr\) and the coefficients \(C_0\), \(C_1\) depend on the local stellar structure via

(9)#\[C_1 = \frac{2}{r} + A\!\left(\frac{2m}{r^2} + 4\pi r (p - \varepsilon)\right)\]

(10)#\[C_0 = A\!\left[-\frac{6}{r^2} + 4\pi\!\left(\varepsilon + p\right)\frac{d\varepsilon}{dp} + 4\pi\!\left(5\varepsilon + 9p\right)\right] - \left(\frac{2(m + 4\pi r^3 p)}{r(r - 2m)}\right)^2\]

with \(A = (1 - 2m/r)^{-1}\) being the metric function. At the stellar surface, the Love number \(k_2\) is extracted from the ratio \(y = R\,\beta(R)/H(R)\) through the closed-form expression (see [4])

(11)#\[k_2 = \frac{\tfrac{8}{5}C^5(1-2C)^2\bigl[2C(y-1) - y + 2\bigr]} {D(C, y)}\]

where the denominator \(D(C, y)\) is

(12)#\[\begin{split}\begin{aligned} D(C, y) &= 2C\!\left[4(y+1)C^4 + (6y-4)C^3 + (26-22y)C^2 + 3(5y-8)C - 3y + 6\right] \\ &\quad - 3(1-2C)^2\bigl[2C(y-1) - y + 2\bigr]\ln\!\frac{1}{1-2C} \end{aligned}\end{split}\]

Numerical implementation#

The four coupled ODEs — Eqs. (4), (5), (7), (8) — are integrated from \(h_c\) down to zero using the Dormand-Prince 5th-order Runge-Kutta method (Dopri5) provided by the Diffrax library, with an adaptive PID step-size controller (relative tolerance \(10^{-5}\), absolute tolerance \(10^{-6}\)).

Initial conditions. The state vector \((r, m, H, \beta)\) is initialised just off-centre at \(h_0 = h_c - 10^{-3} h_c\) using a series expansion valid near \(r = 0\):

(13)#\[r_0 \approx \sqrt{\frac{3(-\Delta h)}{2\pi(\varepsilon_c + 3p_c)}} \left[1 - \frac{(\varepsilon_c - 3p_c - 0.6\,d\varepsilon/dh|_c)\,(-\Delta h)} {4(\varepsilon_c + 3p_c)}\right]\]

(14)#\[m_0 \approx \frac{4\pi \varepsilon_c r_0^3}{3} \left[1 - \frac{0.6\,(d\varepsilon/dh|_c)\,(-\Delta h)}{\varepsilon_c}\right]\]

(15)#\[H_0 = r_0^2, \qquad \beta_0 = 2r_0\]

where \(\Delta h = h_0 - h_c = -10^{-3} h_c\).

Enthalpy clamping. Diffrax’s adaptive step controller can evaluate the ODE at trial points where \(h \leq 0\), causing undefined logarithmic operations in the EOS interpolation routines. jester guards against this by clamping the enthalpy to the minimum tabulated value before any EOS lookup, which freezes the ODE derivatives at their surface value without affecting the converged solution.

Family construction. To build a full mass-radius-tidal family, the TOV equations are solved for a logarithmically spaced grid of central pressures \(p_c\) running from a minimum density (set by min_nsat, in units of saturation density \(n_\mathrm{sat} = 0.16\,\mathrm{fm}^{-3}\)) up to the maximum pressure at which the EOS remains causal (\(c_s^2 < 1\)). JAX’s vmap is used to evaluate all central pressures in parallel, making the family construction efficient on GPU hardware. The stable branch of the \(M(R)\) curve (up to the maximum mass \(M_\mathrm{TOV}\)) is retained, and the family is re-sampled onto a uniform mass grid for downstream inference.

Further resources#

API reference: jesterTOV.tov.gr.GRTOVSolver

References

[1]

Richard C. Tolman. Static solutions of Einstein's field equations for spheres of fluid. Phys. Rev., 55:364–373, 1939. doi:10.1103/PhysRev.55.364.

[2]

J. R. Oppenheimer and G. M. Volkoff. On massive neutron cores. Phys. Rev., 55:374–381, 1939. doi:10.1103/PhysRev.55.374.

[3]

Eanna E. Flanagan and Tanja Hinderer. Constraining neutron star tidal Love numbers with gravitational wave detectors. Phys. Rev. D, 77:021502, 2008. arXiv:0709.1915, doi:10.1103/PhysRevD.77.021502.

[4] (1,2)

Tanja Hinderer. Tidal Love numbers of neutron stars. Astrophys. J., 677:1216–1220, 2008. [Erratum: Astrophys.J. 697, 964 (2009)]. arXiv:0711.2420, doi:10.1086/533487.