Source code for jesterTOV.eos.metamodel.base
r"""Meta-model equation of state for nuclear matter."""
import jax.numpy as jnp
from jax.scipy.special import factorial
from jaxtyping import Array, Float, Int
from jesterTOV import utils
from jesterTOV.eos.base import Interpolate_EOS_model
from jesterTOV.eos.crust import Crust
from jesterTOV.tov.data_classes import EOSData
from jesterTOV.logging_config import get_logger
logger = get_logger("jester")
[docs]
class MetaModel_EOS_model(Interpolate_EOS_model):
r"""
Meta-model equation of state for nuclear matter.
This class implements the meta-modeling approach for nuclear equation of state
as described in Somasundaram 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).
For a more detailed introduction into the metamodel approach, see Margueron et al. (Phys.Rev.C 97 (2018) 2, 025805).
The meta-model uses a kinetic + potential energy decomposition:
.. math::
\varepsilon(n, \delta) = \varepsilon_{\mathrm{kin}}(n, \delta) + \varepsilon_{\mathrm{pot}}(n, \delta)
where :math:`\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 :math:`n_{\rm{sat}}`,
taking into account non-quadratic contributions to the symmetry energy.
"""
[docs]
def __init__(
self,
kappas: tuple[Float, Float, Float, Float, Float, Float] = (
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
),
v_nq: list[float] = [0.0, 0.0, 0.0, 0.0, 0.0],
b_sat: Float = 17.0,
b_sym: Float = 25.0,
# density parameters
nsat: Float = 0.16,
nmin_MM_nsat: Float = 0.75,
nmax_nsat: Float = 12,
ndat: Int = 200,
# crust parameters
crust_name: str = "DH",
max_n_crust_nsat: Float = 0.5,
min_n_crust_nsat: Float = 2e-13,
ndat_spline: Int = 10,
# proton fraction
proton_fraction: bool | float | None = None,
):
r"""
Initialize the meta-model EOS with nuclear empirical parameters.
The meta-model approach parameterizes nuclear matter using empirical parameters
measured from finite nuclei and infinite nuclear matter calculations. This
implementation combines a realistic crust model with a meta-model description
of the core using nuclear empirical parameters (NEPs).
**Reference:** Somasundaram et al., Phys. Rev. C 103, 045803 (2021)
**Physical Framework:**
The meta-model decomposes the energy density as:
.. math::
\varepsilon(n, \delta) = \varepsilon_{\mathrm{kin}}(n, \delta) + \varepsilon_{\mathrm{pot}}(n, \delta)
where :math:`\delta = (n_n - n_p)/n` is the isospin asymmetry parameter.
Args:
kappas (tuple[Float, Float, Float, Float, Float, Float], optional):
Meta-model expansion coefficients :math:`(\kappa_{\mathrm{sat}}, \kappa_{\mathrm{sat2}}, \kappa_{\mathrm{sat3}}, \kappa_{\mathrm{NM}}, \kappa_{\mathrm{NM2}}, \kappa_{\mathrm{NM3}})`.
Controls the density dependence of kinetic energy corrections in the Thomas-Fermi gas approximation.
These parameters modify the kinetic energy beyond the non-relativistic limit.
Defaults to (0.0, 0.0, 0.0, 0.0, 0.0, 0.0).
v_nq (list[float], optional):
Quartic isospin coefficients :math:`v_{\mathrm{nq},\alpha}` for symmetry energy expansion up to :math:`\delta^4` terms.
These control the high-order isospin dependence of the potential energy contribution.
Defaults to [0.0, 0.0, 0.0, 0.0, 0.0].
b_sat (Float, optional):
Saturation parameter :math:`b_{\mathrm{sat}}` controlling potential energy cutoff function.
Higher values lead to sharper exponential cutoffs at high density. Defaults to 17.0.
b_sym (Float, optional):
Symmetry parameter :math:`b_{\mathrm{sym}}` for isospin-dependent cutoff corrections.
Modifies the cutoff parameter as :math:`b = b_{\mathrm{sat}} + b_{\mathrm{sym}} \delta^2`. Defaults to 25.0.
nsat (Float, optional):
Nuclear saturation density :math:`n_0` [:math:`\mathrm{fm}^{-3}`]. The density at which
symmetric nuclear matter reaches its minimum energy per nucleon. Defaults to 0.16.
nmin_MM_nsat (Float, optional):
Starting density for meta-model region as fraction of :math:`n_0`.
Must be above the crust-core transition density. Defaults to 0.75 (= 0.12/0.16).
nmax_nsat (Float, optional):
Maximum density for EOS construction in units of :math:`n_0`.
Determines the high-density reach of the neutron star model. Defaults to 12.
ndat (Int, optional):
Number of density points for meta-model region discretization.
Higher values provide smoother interpolation at computational cost. Defaults to 200.
crust_name (str, optional):
Crust model name (e.g., 'DH', 'BPS') or path to custom .npz file containing crust EOS data.
The crust provides low-density EOS data below nuclear saturation. Defaults to 'DH'.
max_n_crust_nsat (Float, optional):
Maximum crust density as fraction of :math:`n_0`. Defines the crust-core
transition region where spline matching occurs. Defaults to 0.5.
ndat_spline (Int, optional):
Number of points for smooth spline interpolation across crust-core transition.
Ensures thermodynamic consistency and causality preservation. Defaults to 10.
proton_fraction (bool | float | None, optional):
Proton fraction treatment strategy:
- None: Calculate :math:`\beta`-equilibrium (charge neutrality + weak equilibrium)
- float: Use fixed proton fraction value throughout the star
- bool: Use simplified uniform composition model
:math:`\beta`-equilibrium is the physical condition for neutron star matter. Defaults to None.
Note:
The meta-model uses a Thomas-Fermi kinetic energy approximation with relativistic
corrections controlled by the :math:`\kappa` parameters, combined with a potential
energy expansion around saturation density with an exponential cutoff at high densities.
This approach provides a flexible framework for exploring nuclear physics uncertainties
in neutron star structure calculations.
"""
# Save as attributes
self.nsat = nsat
self.v_nq = jnp.array(v_nq)
self.b_sat = b_sat
self.b_sym = b_sym
# NOTE: the expansion is truncated at N=4, consistent with Somasundaram et al. and Margueron et al.
self.N = 4
self.nmin_MM_nsat = nmin_MM_nsat
self.nmax_nsat = nmax_nsat
self.ndat = ndat
self.max_n_crust_nsat = max_n_crust_nsat
self.min_n_crust_nsat = min_n_crust_nsat
self.ndat_spline = ndat_spline
if isinstance(proton_fraction, float):
self.proton_fraction_val = proton_fraction
self.proton_fraction = lambda _v_sat, _v_sym2, _n: self.proton_fraction_val
else:
self.proton_fraction = (
lambda v_sat, v_sym2, n: self.compute_proton_fraction(v_sat, v_sym2, n)
)
# Constructions
assert (
len(kappas) == 6
), "kappas must be a tuple of 6 values: kappa_sat, kappa_sat2, kappa_sat3, kappa_NM, kappa_NM2, kappa_NM3"
(
self.kappa_sat,
self.kappa_sat2,
self.kappa_sat3,
self.kappa_NM,
self.kappa_NM2,
self.kappa_NM3,
) = kappas
self.kappa_sym = self.kappa_NM - self.kappa_sat
self.kappa_sym2 = self.kappa_NM2 - self.kappa_sat2
self.kappa_sym3 = self.kappa_NM3 - self.kappa_sat3
# t_sat or TFGsat is the kinetic energy per nucleons in SM and at saturation, see just after eq (13) in the Margueron paper
self.t_sat = (
3
* utils.hbar**2
/ (10 * utils.m)
* (3 * jnp.pi**2 * self.nsat / 2) ** (2 / 3)
)
# These coefficients are defined in Eq. (B3) of Somasundaram et al
self.v_sat_0_no_NEP = -self.t_sat * (
1 + self.kappa_sat + self.kappa_sat2 + self.kappa_sat3
)
self.v_sat_1_no_NEP = -self.t_sat * (
2 + 5 * self.kappa_sat + 8 * self.kappa_sat2 + 11 * self.kappa_sat3
)
self.v_sat_2_no_NEP = (
-2
* self.t_sat
* (-1 + 5 * self.kappa_sat + 20 * self.kappa_sat2 + 44 * self.kappa_sat3)
)
self.v_sat_3_no_NEP = (
-2
* self.t_sat
* (4 - 5 * self.kappa_sat + 40 * self.kappa_sat2 + 220 * self.kappa_sat3)
)
self.v_sat_4_no_NEP = (
-8
* self.t_sat
* (-7 + 5 * self.kappa_sat - 10 * self.kappa_sat2 + 110 * self.kappa_sat3)
)
# These are the difference between the PNM and SNM coefficients in Appendix B of Somasundaram et al,
# Eq. (B4) - Eq. (B3)
self.v_sym2_0_no_NEP = (
-self.t_sat
* (
2 ** (2 / 3) * (1 + self.kappa_NM + self.kappa_NM2 + self.kappa_NM3)
- (1 + self.kappa_sat + self.kappa_sat2 + self.kappa_sat3)
)
- self.v_nq[0]
)
self.v_sym2_1_no_NEP = (
-self.t_sat
* (
2 ** (2 / 3)
* (2 + 5 * self.kappa_NM + 8 * self.kappa_NM2 + 11 * self.kappa_NM3)
- (2 + 5 * self.kappa_sat + 8 * self.kappa_sat2 + 11 * self.kappa_sat3)
)
- self.v_nq[1]
)
# NOTE the factor at the front, so the power of two is correct here
self.v_sym2_2_no_NEP = (
-2
* self.t_sat
* (
2 ** (2 / 3)
* (-1 + 5 * self.kappa_NM + 20 * self.kappa_NM2 + 44 * self.kappa_NM3)
- (
-1
+ 5 * self.kappa_sat
+ 20 * self.kappa_sat2
+ 44 * self.kappa_sat3
)
)
- self.v_nq[2]
)
# NOTE the factor at the front, so the power of two is correct here
self.v_sym2_3_no_NEP = (
-2
* self.t_sat
* (
2 ** (2 / 3)
* (4 - 5 * self.kappa_NM + 40 * self.kappa_NM2 + 220 * self.kappa_NM3)
- (
4
- 5 * self.kappa_sat
+ 40 * self.kappa_sat2
+ 220 * self.kappa_sat3
)
)
- self.v_nq[3]
)
# NOTE the factor at the front, so the power of two is correct here
self.v_sym2_4_no_NEP = (
-8
* self.t_sat
* (
2 ** (2 / 3)
* (-7 + 5 * self.kappa_NM - 10 * self.kappa_NM2 + 110 * self.kappa_NM3)
- (
-7
+ 5 * self.kappa_sat
- 10 * self.kappa_sat2
+ 110 * self.kappa_sat3
)
)
- self.v_nq[4]
)
# Load and preprocess the crust
crust = Crust(
crust_name,
min_density=min_n_crust_nsat * nsat,
max_density=max_n_crust_nsat * nsat,
filter_zero_pressure=True,
)
self.ns_crust: Float[Array, "n_crust"] = crust.n
self.ps_crust: Float[Array, "n_crust"] = crust.p
self.es_crust: Float[Array, "n_crust"] = crust.e
self.mu_lowest: Float = crust.mu_lowest
self.cs2_crust: Float[Array, "n_crust"] = crust.cs2
# Make sure the metamodel starts above the crust
self.max_n_crust = self.ns_crust[-1]
# Create density arrays
self.nmax = nmax_nsat * self.nsat
self.ndat = ndat
self.nmin_MM = self.nmin_MM_nsat * self.nsat
self.n_metamodel = jnp.logspace(
jnp.log10(self.nmin_MM), jnp.log10(self.nmax), self.ndat, endpoint=False
)
self.ns_spline = jnp.append(self.ns_crust, self.n_metamodel)
self.n_connection = jnp.linspace(
self.max_n_crust + 1e-5, self.nmin_MM, self.ndat_spline, endpoint=False
)
[docs]
def construct_eos(self, params: dict[str, float]) -> EOSData:
r"""
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.
Args:
params: Nuclear empirical parameters including:
- **E_sat**: Saturation energy per nucleon [:math:`\mathrm{MeV}`] (default: -16.0)
- **K_sat**: Incompressibility at saturation [:math:`\mathrm{MeV}`]
- **Q_sat, Z_sat**: Higher-order saturation parameters [:math:`\mathrm{MeV}`]
- **E_sym**: Symmetry energy [:math:`\mathrm{MeV}`]
- **L_sym**: Symmetry energy slope [:math:`\mathrm{MeV}`]
- **K_sym, Q_sym, Z_sym**: Higher-order symmetry parameters [:math:`\mathrm{MeV}`]
Returns:
EOSData: Complete EOS with all required arrays
"""
# Use the old parameter name internally for compatibility
NEP_dict = params
E_sat = NEP_dict["E_sat"]
K_sat = NEP_dict["K_sat"]
Q_sat = NEP_dict["Q_sat"]
Z_sat = NEP_dict["Z_sat"]
E_sym = NEP_dict["E_sym"]
L_sym = NEP_dict["L_sym"]
K_sym = NEP_dict["K_sym"]
Q_sym = NEP_dict["Q_sym"]
Z_sym = NEP_dict["Z_sym"]
# Potential energy: see Appendix B of Somasundaram et al. for details
v_sat = jnp.array(
[
E_sat + self.v_sat_0_no_NEP,
0.0 + self.v_sat_1_no_NEP,
K_sat + self.v_sat_2_no_NEP,
Q_sat + self.v_sat_3_no_NEP,
Z_sat + self.v_sat_4_no_NEP,
]
)
v_sym2 = jnp.array(
[
E_sym + self.v_sym2_0_no_NEP,
L_sym + self.v_sym2_1_no_NEP,
K_sym + self.v_sym2_2_no_NEP,
Q_sym + self.v_sym2_3_no_NEP,
Z_sym + self.v_sym2_4_no_NEP,
]
)
# Auxiliaries first
x = self.compute_x(self.n_metamodel)
proton_fraction = self.proton_fraction(v_sat, v_sym2, self.n_metamodel)
delta = 1 - 2 * proton_fraction
f_1 = self.compute_f_1(delta)
f_star = self.compute_f_star(delta)
f_star2 = self.compute_f_star2(delta)
f_star3 = self.compute_f_star3(delta)
v = self.compute_v(v_sat, v_sym2, delta)
b = self.compute_b(delta)
u = self.compute_u(x, b)
# Other quantities
# NOTE: The following functions only compute pressure and energy density for the nucleonic part of the EOS
# The electron contribution is not added here.
# Instead, the electron contribution is added in the compute_cs2 function, which is then integrated
# to obtain the total energy density and pressure later on
p_metamodel = self.compute_pressure_nucleons(
x, f_1, f_star, f_star2, f_star3, b, v, u
)
e_metamodel = self.compute_energy_nucleons(
x, f_1, f_star, f_star2, f_star3, v, u
)
# Get cs2 for the metamodel
cs2_metamodel = self.compute_cs2(
self.n_metamodel,
p_metamodel,
e_metamodel,
x,
delta,
f_1,
f_star,
f_star2,
f_star3,
b,
v,
u,
)
# Spline for speed of sound for the connection region
cs2_spline = jnp.append(jnp.array(self.cs2_crust), cs2_metamodel)
cs2_connection = utils.cubic_spline(
self.n_connection, self.ns_spline, cs2_spline
)
cs2_connection = jnp.clip(cs2_connection, 1e-5, 1.0)
# Concatenate the arrays
n = jnp.concatenate([self.ns_crust, self.n_connection, self.n_metamodel])
cs2 = jnp.concatenate(
[jnp.array(self.cs2_crust), cs2_connection, cs2_metamodel]
)
# Compute pressure and energy from chemical potential and initialize the parent class with it
log_mu = utils.cumtrapz(cs2, jnp.log(n)) + jnp.log(self.mu_lowest)
mu = jnp.exp(log_mu)
p = utils.cumtrapz(cs2 * mu, n) + self.ps_crust[0]
e = mu * n - p
ns, ps, hs, es, dloge_dlogps = self.interpolate_eos(n, p, e)
# Count points where symmetry energy is negative (unphysical region)
esym_vals = self.esym(v_sat, v_sym2, self.n_metamodel)
n_esym_violations = jnp.sum(esym_vals < 0.0)
extra_constraints = {"n_esym_violations": n_esym_violations}
return EOSData(
ns=ns,
ps=ps,
hs=hs,
es=es,
dloge_dlogps=dloge_dlogps,
cs2=cs2,
mu=mu,
extra_constraints=extra_constraints,
)
[docs]
def get_required_parameters(self) -> list[str]:
"""
Return list of nuclear empirical parameters required by MetaModel.
Returns:
list[str]: ["E_sat", "K_sat", "Q_sat", "Z_sat", "E_sym", "L_sym", "K_sym", "Q_sym", "Z_sym"]
"""
return [
"E_sat",
"K_sat",
"Q_sat",
"Z_sat",
"E_sym",
"L_sym",
"K_sym",
"Q_sym",
"Z_sym",
]
#################
### AUXILIARY ###
#################
[docs]
def u(self, x: Array, b: Array, alpha: Int):
"Defined right after equation (40) in Margueron et al. Note that Somasundaram et al. has more than 1 b parameter, compared to Margueron et al., but the definition of u is similar. For a single u value"
return 1 - ((-3 * x) ** (self.N + 1 - alpha) * jnp.exp(-b * (1 + 3 * x)))
[docs]
def compute_x(self, n: Array):
"Defined right before equation (2) in Margueron et al."
return (n - self.nsat) / (3 * self.nsat)
[docs]
def compute_f_1(self, delta: Array | float):
"Equation (14) in Margueron et al."
return (1 + delta) ** (5 / 3) + (1 - delta) ** (5 / 3)
[docs]
def compute_f_star(self, delta: Array | float):
return (self.kappa_sat + self.kappa_sym * delta) * (1 + delta) ** (5 / 3) + (
self.kappa_sat - self.kappa_sym * delta
) * (1 - delta) ** (5 / 3)
[docs]
def compute_f_star2(self, delta: Array | float):
return (self.kappa_sat2 + self.kappa_sym2 * delta) * (1 + delta) ** (5 / 3) + (
self.kappa_sat2 - self.kappa_sym2 * delta
) * (1 - delta) ** (5 / 3)
[docs]
def compute_f_star3(self, delta: Array | float):
return (self.kappa_sat3 + self.kappa_sym3 * delta) * (1 + delta) ** (5 / 3) + (
self.kappa_sat3 - self.kappa_sym3 * delta
) * (1 - delta) ** (5 / 3)
[docs]
def compute_u(self, x: Array, b: Array):
"Computes an array of N+1 u values for all alphas"
return jnp.array([self.u(x, b, alpha) for alpha in range(self.N + 1)])
[docs]
def compute_v(self, v_sat: Array, v_sym2: Array, delta: Array | float) -> Array:
"""Compute an array of length N+1 containing the v coefficients for the potential energy"""
return jnp.array(
[
v_sat[alpha] + v_sym2[alpha] * delta**2 + self.v_nq[alpha] * delta**4
for alpha in range(self.N + 1)
]
)
[docs]
def compute_energy_nucleons(
self,
x: Array,
f_1: Array,
f_star: Array,
f_star2: Array,
f_star3: Array,
v: Array,
u: Array,
) -> Array:
r"""Compute the nucleon-only energy per nucleon :math:`e(n, \delta)` [MeV].
This returns the kinetic plus potential energy per nucleon for the nuclear
matter part only — electrons are not included. The result is an intermediate
quantity passed to :meth:`compute_cs2`, which adds the electron contribution
before computing the sound speed. The final energy density stored in
:class:`~jesterTOV.tov.data_classes.EOSData` is obtained by integrating the
total cs2 (nucleons + electrons) via the thermodynamic identity, so it
consistently includes electrons.
"""
# Kinetic energy
prefac = self.t_sat / 2 * (1 + 3 * x) ** (2 / 3)
linear = (1 + 3 * x) * f_star
quadratic = (1 + 3 * x) ** 2 * f_star2
cubic = (1 + 3 * x) ** 3 * f_star3
kinetic_energy = prefac * (f_1 + linear + quadratic + cubic)
# Potential energy
potential_energy = 0
for alpha in range(self.N + 1):
potential_energy += v[alpha] / (factorial(alpha)) * x**alpha * u[alpha]
return kinetic_energy + potential_energy
[docs]
def esym(self, v_sat: Array, v_sym2: Array, n: Array) -> Array:
r"""Symmetry energy :math:`e_{\rm sym}(n)` [MeV].
Computed as the difference in energy per nucleon between pure neutron matter
(:math:`\delta = 1`) and symmetric nuclear matter (:math:`\delta = 0`):
.. math::
e_{\rm sym}(n) = e(n, \delta=1) - e(n, \delta=0)
This self-consistent approach uses the full metamodel kinetic and potential
energy (including the u-function corrections) rather than the NEP polynomial
approximation, giving more accurate proton fractions especially below
saturation density.
Args:
v_sat: Potential energy coefficients for isospin-symmetric matter, shape (N+1,).
v_sym2: Potential energy coefficients for the isospin-asymmetric part, shape (N+1,).
n: Number density array [:math:`\mathrm{fm}^{-3}`].
Returns:
Symmetry energy at each density point [MeV].
"""
x = self.compute_x(n)
# Pure neutron matter (delta = 1)
e_pnm = self.compute_energy_nucleons(
x,
self.compute_f_1(1.0),
self.compute_f_star(1.0),
self.compute_f_star2(1.0),
self.compute_f_star3(1.0),
self.compute_v(v_sat, v_sym2, 1.0),
self.compute_u(x, self.compute_b(1.0)),
)
# Symmetric nuclear matter (delta = 0)
e_snm = self.compute_energy_nucleons(
x,
self.compute_f_1(0.0),
self.compute_f_star(0.0),
self.compute_f_star2(0.0),
self.compute_f_star3(0.0),
self.compute_v(v_sat, v_sym2, 0.0),
self.compute_u(x, self.compute_b(0.0)),
)
return e_pnm - e_snm
[docs]
def compute_pressure_nucleons(
self,
x: Array,
f_1: Array,
f_star: Array,
f_star2: Array,
f_star3: Array,
b: Array,
v: Array,
u: Array,
) -> Array:
r"""Compute the nucleon-only pressure [MeV fm\ :sup:`-3`].
Returns the pressure of nuclear matter (protons + neutrons) only:
.. math::
P_{\rm nuc}(n, \delta) = n^2 \frac{\partial e(n, \delta)}{\partial n}
Electrons are not included here. This is an intermediate quantity used
inside :meth:`compute_cs2` to evaluate the nuclear incompressibility
:math:`K_{\rm nuc}`. The total cs2 from :meth:`compute_cs2` adds the
electron incompressibility and enthalpy, and integrating that total cs2
via the thermodynamic identity yields the final pressure stored in
:class:`~jesterTOV.tov.data_classes.EOSData`, which does include electrons.
"""
# Contribution from kinetic energy
p_kin = (
1
/ 3
* self.nsat
* self.t_sat
* (1 + 3 * x) ** (5 / 3)
* (
f_1
+ 5 / 2 * (1 + 3 * x) * f_star
+ 4 * (1 + 3 * x) ** 2 * f_star2
+ 11 / 2 * (1 + 3 * x) ** 3 * f_star3
)
)
# Contribution from potential energy
p_pot = 0
for alpha in range(1, self.N + 1):
fac1 = alpha * u[alpha]
fac2 = (self.N + 1 - alpha - 3 * b * x) * (u[alpha] - 1)
p_pot += v[alpha] / (factorial(alpha)) * x ** (alpha - 1) * (fac1 + fac2)
p_pot = p_pot - v[0] * (-3) ** (self.N + 1) * x**self.N * (
self.N + 1 - 3 * b * x
) * jnp.exp(-b * (1 + 3 * x))
p_pot = p_pot * (1 / 3) * self.nsat * (1 + 3 * x) ** 2
return p_pot + p_kin
[docs]
def compute_cs2(
self,
n: Array,
p: Array,
e: Array,
x: Array,
delta: Array | float,
f_1: Array,
f_star: Array,
f_star2: Array,
f_star3: Array,
b: Array,
v: Array,
u: Array,
):
### Compute incompressibility
# Kinetic part
K_kin = (
self.t_sat
* (1 + 3 * x) ** (2 / 3)
* (
-f_1
+ 5 * (1 + 3 * x) * f_star
+ 20 * (1 + 3 * x) ** 2 * f_star2
+ 44 * (1 + 3 * x) ** 3 * f_star3
)
)
# Build the contributions from the potential part
K_pot = 0
# alpha = 0 term
K_pot += (
v[0]
* (-(self.N + 1) * (self.N) + 6 * b * x * (self.N + 1) - 9 * x**2 * b**2)
* ((-3) ** (self.N + 1) * x ** (self.N - 1) * jnp.exp(-b * (1 + 3 * x)))
)
# alpha = 1 term
K_pot += (
2
* v[1]
* (self.N - 3 * b * x)
* (-((-3) ** (self.N)) * x ** (self.N - 1) * jnp.exp(-b * (1 + 3 * x)))
)
K_pot += (
v[1]
* (-(self.N) * (self.N - 1) + 6 * b * x * (self.N) - 9 * x**2 * b**2)
* ((-3) ** (self.N) * x ** (self.N - 1) * jnp.exp(-b * (1 + 3 * x)))
)
# Summation over alpha >= 2 terms
for alpha in range(2, self.N + 1):
# x times u prime
x_up = (self.N + 1 - alpha - 3 * b * x) * (u[alpha] - 1)
# x squared times u double prime
x2_upp = (
-(self.N + 1 - alpha) * (self.N - alpha)
+ 6 * b * x * (self.N + 1 - alpha)
- 9 * x**2 * b**2
) * (1 - u[alpha])
K_pot = K_pot + v[alpha] / (factorial(alpha)) * x ** (alpha - 2) * (
alpha * (alpha - 1) * u[alpha] + 2 * alpha * x_up + x2_upp
)
# Multiply by the overall prefactor
K_pot *= (1 + 3 * x) ** 2
# Sum the kinetic and potential contributions, and add the contribution from the pressure
K = K_kin + K_pot + 18 / n * p
# Add electron contributions
K_Fe = (
(3.0 * jnp.pi**2 / 2.0 * n) ** (1.0 / 3.0)
* utils.hbarc
* (1.0 - delta) ** (1.0 / 3.0)
)
x_e = K_Fe / utils.m_e
C = utils.m_e**4 / (8.0 * jnp.pi**2) / utils.hbarc**3
f = x_e * (1 + 2 * x_e**2) * jnp.sqrt(1 + x_e**2) - jnp.arcsinh(x_e)
epsilon_electron = C * f
p_electron = -epsilon_electron + 8.0 / 3.0 * C * x_e**3 * jnp.sqrt(1 + x_e**2)
K_electron = 8 * C / n * x_e**3 * (3 + 4 * x_e**2) / (
jnp.sqrt(1 + x_e**2)
) - 9 / n * (epsilon_electron + p_electron)
# Sum together:
K_tot = K + K_electron
# Finally, get cs2:
chi = K_tot / 9.0
total_energy_density = (e + utils.m) * n + epsilon_electron
total_pressure = p + p_electron
h_tot = (total_energy_density + total_pressure) / n
cs2 = chi / h_tot
return cs2
[docs]
def compute_proton_fraction(
self, v_sat: Array, v_sym2: Array, n: Array
) -> Float[Array, "n_points"]:
r"""
Compute proton fraction from beta-equilibrium condition.
Solves the beta-equilibrium condition:
.. math::
\mu_n - \mu_p = \mu_e
using the quadratic approximation :math:`\mu_n - \mu_p \approx 4 e_{\rm sym}(n) \delta`
together with the self-consistent symmetry energy from :meth:`esym`. Protons and
neutrons are treated as having equal rest mass (:math:`m_n = m_p`).
While this does not perform an exact root-finding step, it provides a good approximation
that is much more efficient to evaluate during sampling than a full root-finder.
Args:
v_sat: Potential energy coefficients for isospin-symmetric matter, shape (N+1,).
v_sym2: Potential energy coefficients for the isospin-asymmetric part, shape (N+1,).
n: Number density [:math:`\mathrm{fm}^{-3}`].
Returns:
Float[Array, "n_points"]: Proton fraction :math:`x_p = n_p/n` as a function of density.
"""
esym = self.esym(v_sat, v_sym2, n)
a = 8.0 * esym
b = jnp.zeros(shape=n.shape)
c = utils.hbarc * jnp.power(3.0 * jnp.pi**2 * n, 1.0 / 3.0)
d = -4.0 * esym - (utils.m_n - utils.m_p)
coeffs = jnp.stack([a, b, c, d], axis=1)
ys = utils.cubic_root_for_proton_fraction(coeffs)
proton_fraction = ys**3 # type: ignore[operator]
return proton_fraction