jesterTOV.eos module

The equation of state (EOS) module provides classes and functions for modeling neutron star matter equations of state.

Mathematical Background

The metamodel equation of state is based on the formulation from Margueron et al. (2021), which provides a flexible framework for modeling nuclear matter at various densities.

The total energy density is given by:

\[\varepsilon(n, \delta) = \varepsilon_{\text{kinetic}}(n, \delta) + \varepsilon_{\text{potential}}(n, \delta)\]

where \(n\) is the baryon number density and \(\delta = 1 - 2Y_p\) is the isospin asymmetry parameter.

API Reference

class Interpolate_EOS_model

Bases: object

Base class for interpolating equation of state (EOS) data.

This class provides the fundamental interpolation framework for converting tabulated EOS data (density, pressure, energy) into the auxiliary quantities needed for neutron star structure calculations using the TOV equations.

interpolate_eos(n, p, e)

Convert physical EOS quantities to geometric units and compute auxiliary quantities.

This method transforms the input EOS data from nuclear physics units to geometric units used in general relativity calculations, and computes derived quantities needed for the TOV equations.

Parameters:

• n (Float[Array, n_points]) – Number densities [\(\mathrm{fm}^{-3}\)]

• p (Float[Array, n_points]) – Pressure values [\(\mathrm{MeV} \, \mathrm{fm}^{-3}\)]

• e (Float[Array, n_points]) – Energy densities [\(\mathrm{MeV} \, \mathrm{fm}^{-3}\)]

Returns:

tuple –

A tuple containing (all in geometric units):

• ns: Number densities

• ps: Pressures

• hs: Specific enthalpy \(h = \int \frac{dp}{\varepsilon + p}\)

• es: Energy densities

• dloge_dlogps: Logarithmic derivative \(\frac{d\ln\varepsilon}{d\ln p}\)

class MetaModel_EOS_model(

kappas=(0.0, 0.0, 0.0, 0.0, 0.0, 0.0),

v_nq=[0.0, 0.0, 0.0, 0.0, 0.0],

b_sat=17.0,

b_sym=25.0,

nsat=0.16,

nmin_MM_nsat=0.75,

nmax_nsat=12,

ndat=200,

crust_name='DH',

max_n_crust_nsat=0.5,

min_n_crust_nsat=2e-13,

ndat_spline=10,

proton_fraction=None,

)

Bases: Interpolate_EOS_model

Meta-model equation of state for nuclear matter.

This class implements the meta-modeling approach for nuclear equation of state as described in Margueron et al. (Phys. Rev. C 103, 045803, 2021). The EOS is constructed by combining a realistic crust model with a meta-model for core densities based on nuclear empirical parameters (NEPs).

The meta-model uses a kinetic + potential energy decomposition:

\[\varepsilon(n, \delta) = \varepsilon_{\mathrm{kin}}(n, \delta) + \varepsilon_{\mathrm{pot}}(n, \delta)\]

where \(\delta = (n_n - n_p)/n\) is the isospin asymmetry parameter.

The kinetic part is based on a Thomas-Fermi gas with relativistic corrections, while the potential part uses a polynomial expansion around saturation density \(n_0\).

compute_b(delta)

compute_cs2(n, p, e, x, delta, f_1, f_star, f_star2, f_star3, b, v)

compute_energy(x, f_1, f_star, f_star2, f_star3, b, v)

Return type:

Array

compute_f_1(delta)

compute_f_star(delta)

compute_f_star2(delta)

compute_f_star3(delta)

compute_pressure(x, f_1, f_star, f_star2, f_star3, b, v)

Return type:

Array

compute_proton_fraction(coefficient_sym, n)

Compute proton fraction from beta-equilibrium condition.

This method solves the beta-equilibrium condition:

\[\mu_e + \mu_p - \mu_n = 0\]

where the chemical potentials are related to the EOS through:

\[\mu_p - \mu_n = \frac{\partial \varepsilon}{\partial x_p} = -4 E_{\mathrm{sym}} (1 - 2x_p)\]

and the electron chemical potential is \(\mu_e = \hbar c (3\pi^2 x_p n)^{1/3}\).

Parameters:

• coefficient_sym (list) – Symmetry energy expansion coefficients.

• n (Float[Array, "n_points"]) – Number density [\(\mathrm{fm}^{-3}\)].

Return type:

Float[Array, 'n_points']

Returns:

Float[Array, “n_points”] – Proton fraction \(x_p = n_p/n\) as a function of density.

compute_v(v_sat, v_sym2, delta)

Return type:

Array

compute_x(n)

construct_eos(NEP_dict)

Construct the complete equation of state from nuclear empirical parameters.

This method builds the full EOS by combining the crust model with the meta-model core, ensuring thermodynamic consistency and causality.

Parameters:

NEP_dict (dict) –

Nuclear empirical parameters including:

• E_sat: Saturation energy per nucleon [\(\mathrm{MeV}\)] (default: -16.0)

• K_sat: Incompressibility at saturation [\(\mathrm{MeV}\)]

• Q_sat, Z_sat: Higher-order saturation parameters [\(\mathrm{MeV}\)]

• E_sym: Symmetry energy [\(\mathrm{MeV}\)]

• L_sym: Symmetry energy slope [\(\mathrm{MeV}\)]

• K_sym, Q_sym, Z_sym: Higher-order symmetry parameters [\(\mathrm{MeV}\)]

Return type:

tuple

Returns:

tuple –

Complete EOS data containing:

• ns: Number densities [geometric units]

• ps: Pressures [geometric units]

• hs: Specific enthalpies [geometric units]

• es: Energy densities [geometric units]

• dloge_dlogps: Logarithmic derivative \(\frac{d\ln\varepsilon}{d\ln p}\)

• mu: Chemical potential [geometric units]

• cs2: Speed of sound squared \(c_s^2 = \frac{dp}{d\varepsilon}\)

esym(coefficient_sym, x)

u(x, b, alpha)

class MetaModel_with_CSE_EOS_model(

nsat=0.16,

nmin_MM_nsat=0.75,

nmax_nsat=12,

ndat_metamodel=100,

ndat_CSE=100,

**metamodel_kwargs,

)

Bases: Interpolate_EOS_model

Meta-model EOS combined with piecewise speed-of-sound extensions (CSE).

This class extends the meta-model approach by allowing for piecewise-constant speed-of-sound extensions at high densities. This is useful for modeling phase transitions or exotic matter components in neutron star cores that may not be captured by the meta-model polynomial expansions.

The EOS is constructed in two regions:

1. Low-to-intermediate density: Meta-model approach (crust + core)

2. High density: Speed-of-sound extension scheme

construct_eos(NEP_dict, ngrids, cs2grids)

Construct the EOS

Parameters:

• NEP_dict (dict) – Dictionary with the NEP keys to be passed to the metamodel EOS class.

• ngrids (Float[Array, n_grid_point]) – Density grid points of densities for the CSE part of the EOS.

• cs2grids (Float[Array, n_grid_point]) – Speed-of-sound squared grid points of densities for the CSE part of the EOS.

Return type:

tuple

Returns:

tuple – EOS quantities (see Interpolate_EOS_model), as well as the chemical potential and speed of sound.

class MetaModel_with_peakCSE_EOS_model(

nsat=0.16,

nmin_MM_nsat=0.75,

nmax_nsat=12,

ndat_metamodel=100,

ndat_CSE=100,

**metamodel_kwargs,

)

Bases: Interpolate_EOS_model

Meta-model EOS with Gaussian peak Constant Speed-of-sound Extensions (peakCSE).

This class implements a sophisticated CSE parametrization based on the peakCSE model, which combines a Gaussian peak structure with logistic growth to model phase transitions while ensuring asymptotic consistency with perturbative QCD (pQCD) at the highest densities.

Mathematical Framework: The speed of sound squared is parametrized as:

\[c^2_s &= c^2_{s,{\rm break}} + \frac{\frac{1}{3} - c^2_{s,{\rm break}}}{1 + e^{-l_{\rm sig}(n - n_{\rm sig})}} + c^2_{s,{\rm peak}}e^{-\frac{1}{2}\left(\frac{n - n_{\rm peak}}{\sigma_{\rm peak}}\right)^2}\]

Reference: Greif:2018njt, arXiv:1812.08188

Note

The peakCSE model provides greater physical realism than simple piecewise-constant CSE by incorporating smooth transitions and theoretically motivated high-density behavior.

construct_eos(NEP_dict, peakCSE_dict)

Construct the complete EOS using meta-model + peakCSE extensions.

This method builds the full EOS by combining the meta-model approach with peakCSE extensions that model phase transitions through Gaussian peaks and approach the pQCD conformal limit at high densities.

Parameters:

• NEP_dict (dict) – Nuclear empirical parameters for meta-model construction. Must include ‘nbreak’ key specifying the transition density between meta-model and peakCSE regions. See MetaModel_EOS_model.construct_eos for complete NEP parameter descriptions.

• peakCSE_dict (dict) –

  peakCSE model parameters defining the high-density behavior:

  • gaussian_peak (float): Amplitude \(A\) of the Gaussian peak

  • gaussian_mu (float): Peak location \(\mu\) [\(\mathrm{fm}^{-3}\)]

  • gaussian_sigma (float): Peak width \(\sigma\) [\(\mathrm{fm}^{-3}\)]

  • logit_growth_rate (float): Growth rate \(k\) for pQCD approach

  • logit_midpoint (float): Midpoint density \(n_{\mathrm{mid}}\) for logistic transition

Returns:

tuple –

Complete EOS data containing:

• ns: Number densities [geometric units]

• ps: Pressures [geometric units]

• hs: Specific enthalpies [geometric units]

• es: Energy densities [geometric units]

• dloge_dlogps: Logarithmic derivative \(\frac{d\ln\varepsilon}{d\ln p}\)

• mu: Chemical potential [geometric units]

• cs2: Speed of sound squared including peakCSE structure

Note

The peakCSE speed of sound follows: \(c^2_s &= c^2_{s,{\rm break}} + \frac{\frac{1}{3} - c^2_{s,{\rm break}}}{1 + e^{-l_{\rm sig}(n - n_{\rm sig})}} + c^2_{s,{\rm peak}}e^{-\frac{1}{2}\left(\frac{n - n_{\rm peak}}{\sigma_{\rm peak}}\right)^2}\)

This ensures smooth transitions, realistic phase transition modeling, and asymptotic consistency with the pQCD conformal limit \(c_s^2 = 1/3\).

offset_calc(nbreak, cs2_break, peakCSE_dict)

construct_family(eos, ndat=50, min_nsat=2)

Solve the TOV equations and generate mass-radius-tidal deformability relations.

This function constructs a neutron star family by solving the Tolman-Oppenheimer-Volkoff (TOV) equations for a range of central pressures. The TOV equations describe the hydrostatic equilibrium of a spherically symmetric, static star:

\[\begin{split}\frac{dm}{dr} &= 4\pi r^2 \varepsilon(r) \\ \frac{dp}{dr} &= -\frac{[\varepsilon(r) + p(r)][m(r) + 4\pi r^3 p(r)]}{r[r - 2m(r)]}\end{split}\]

Parameters:

• eos (tuple) – Tuple of the EOS data (ns, ps, hs, es, dloge_dlogps).

• ndat (int, optional) – Number of datapoints used when constructing the central pressure grid. Defaults to 50.

• min_nsat (int, optional) – Starting density for central pressure in numbers of \(n_0\) (assumed to be 0.16 \(\mathrm{fm}^{-3}\)). Defaults to 2.

Return type:

tuple[Float[Array, 'ndat'], Float[Array, 'ndat'], Float[Array, 'ndat'], Float[Array, 'ndat']]

Returns:

tuple –

A tuple containing:

• \(\log(p_c)\): Logarithm of central pressures [geometric units]

• \(M\): Gravitational masses [\(M_{\odot}\)]

• \(R\): Circumferential radii [\(\mathrm{km}\)]

• \(\Lambda\): Dimensionless tidal deformabilities

construct_family_ST(eos, ndat=50, min_nsat=2)

# TODO: (updated later)

Return type:

tuple[Float[Array, 'ndat'], Float[Array, 'ndat'], Float[Array, 'ndat'], Float[Array, 'ndat']]

construct_family_ST_sol(eos, ndat=1, min_nsat=2)

# TODO: complete the description Also output stellar structure solution via sol_iter (interior) and solext (exterior)

Return type:

tuple[Float[Array, 'ndat'], Float[Array, 'ndat'], Float[Array, 'ndat'], Float[Array, 'ndat'], Float[Array, 'ndat'], Float[Array, 'ndat']]

construct_family_nonGR(eos, ndat=50, min_nsat=2)

Solve modified TOV equations with beyond-GR corrections.

This function extends the standard TOV equations to include phenomenological modifications that parameterize deviations from General Relativity. The modified pressure gradient equation becomes:

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

where \(\sigma(r)\) contains the non-GR corrections parameterized by \(\lambda_{\mathrm{BL}}\), \(\lambda_{\mathrm{DY}}\), \(\lambda_{\mathrm{HB}}\), and post-Newtonian parameters \(\alpha\), \(\beta\), \(\gamma\).

Parameters:

• eos (tuple) – Extended EOS data including GR modification parameters.

• ndat (int, optional) – Number of datapoints for central pressure grid. Defaults to 50.

• min_nsat (int, optional) – Starting density in units of \(n_0\). Defaults to 2.

Return type:

tuple[Float[Array, 'ndat'], Float[Array, 'ndat'], Float[Array, 'ndat'], Float[Array, 'ndat']]

Returns:

tuple –

A tuple containing:

• \(\log(p_c)\): Logarithm of central pressures [geometric units]

• \(M\): Gravitational masses [\(M_{\odot}\)]

• \(R\): Circumferential radii [\(\mathrm{km}\)]

• \(\Lambda\): Dimensionless tidal deformabilities

load_crust(name)

Load neutron star crust equation of state data from the default directory.

This function loads pre-computed crust EOS data from tabulated files, supporting both built-in crust models (BPS, DH) and custom files.

Parameters:

name (str) – Name of the crust model to load (e.g., ‘BPS’, ‘DH’) or a filename with .npz extension for custom files.

Return type:

tuple[Array, Array, Array]

Returns:

tuple[Array, Array, Array] –

A tuple containing:

• Number densities \(n\) [\(\mathrm{fm}^{-3}\)]

• Pressures \(p\) [\(\mathrm{MeV} \, \mathrm{fm}^{-3}\)]

• Energy densities \(\varepsilon\) [\(\mathrm{MeV} \, \mathrm{fm}^{-3}\)]

Raises:

ValueError – If the specified crust model is not found in the default directory.

Note

Available built-in crust models can be found in the crust/ subdirectory.

locate_lowest_non_causal_point(cs2)

Find the first point where the equation of state becomes non-causal.

The speed of sound squared \(c_s^2 = dp/d\varepsilon\) must satisfy \(c_s^2 \leq 1\) (in units where \(c = 1\)) for causality. This function locates the first density where this condition is violated.

Parameters:

cs2 (Array) – Speed of sound squared values.

Returns:

int – Index of first non-causal point, or -1 if EOS is everywhere causal.