r"""
Post-TOV equation solver with beyond-GR corrections.
This module extends the standard TOV equations to include phenomenological
modifications that parameterize deviations from General Relativity. The
modifications are implemented through additional sigma terms in the pressure
gradient equation.
**Units:** All calculations are performed in geometric units where :math:`G = c = 1`.
**Reference:** Yagi & Yunes, Phys. Rev. D 88, 023009 (2013)
"""
import jax
import jax.numpy as jnp
from diffrax import diffeqsolve, ODETerm, Dopri5, SaveAt, PIDController
from jesterTOV import utils
from jesterTOV.tov.base import TOVSolverBase
from jesterTOV.tov.data_classes import EOSData, TOVSolution
def _sigma_func(p, e, m, r, lambda_BL, lambda_DY, lambda_HB, gamma, alpha, beta):
r"""
Compute the non-GR correction term sigma for post-TOV equations.
This function implements various phenomenological modifications to
General Relativity that appear as additional terms in the TOV equations.
The corrections are parameterized by several coupling constants.
The sigma function includes:
- **Bowers-and-Liang**: :math:`\sigma_{\mathrm{BL}} = -\frac{\lambda_{\mathrm{BL}}\epsilon^2 r^2}{3}\left(1 + \frac{3p}{\epsilon}\right)\frac{\left(1 + \frac{p}{\epsilon}\right)}{\left(1 - \frac{2M}{r}\right)}`
- **Horvat et al.**: :math:`\sigma_{\mathrm{DY}} = \lambda_{\mathrm{DY}} \frac{2m}{r} p`
- **Cosenza et al.**: :math:`\sigma_{\mathrm{HB}} = -(\frac{1}{\lambda_{\mathrm{HB}}} - 1) \frac{r}{2} \frac{dp}{dr}`
- **Post-Newtonian**: :math:`\sigma_{\mathrm{PN}} = \gamma \frac{2m}{r} p \tanh(\alpha(\frac{m}{r} - \beta))`
Parameters
----------
p : float
Pressure at current radius
e : float
Energy density at current radius
m : float
Enclosed mass at current radius
r : float
Current radius
lambda_BL : float
Bowers-Liang coupling parameter
lambda_DY : float
Doneva-Yazadjiev coupling parameter
lambda_HB : float
Herrera-Barreto coupling parameter
gamma : float
Post-Newtonian amplitude parameter
alpha : float
Post-Newtonian steepness parameter
beta : float
Post-Newtonian transition point parameter
Returns
-------
float
Total sigma correction term
"""
# Metric coefficient A = 1/(1-2m/r)
A = 1.0 / (1.0 - 2.0 * m / r)
dpdr = -(e + p) * (m + 4.0 * jnp.pi * r * r * r * p) / r / r * A
sigma = 0.0
# models reviewed in https://doi.org/10.1140/epjc/s10052-020-8361-4
# in Eq. 12, the power of epsilon should be 2
sigma += -lambda_BL * r * r / 3.0 * (e + 3.0 * p) * (e + p) * A
sigma += lambda_DY * 2.0 * m / r * p
sigma += -(1.0 / lambda_HB - 1.0) * r / 2.0 * dpdr
# post-Newtonian modification
sigma += gamma * 2.0 * m / r * p * jnp.tanh(alpha * (m / r - beta))
return sigma
def _tov_ode(h, y, eos):
r"""
Post-TOV ordinary differential equation system.
Includes beyond-GR corrections through the sigma function.
Parameters
----------
h : float
Enthalpy (independent variable)
y : tuple
State vector (r, m, H, β)
eos : dict
EOS interpolation data plus modification parameters
Returns
-------
tuple
Derivatives (dr/dh, dm/dh, dH/dh, dβ/dh)
"""
# fetch the eos arrays
ps = eos["p"]
hs = eos["h"]
es = eos["e"]
dloge_dlogps = eos["dloge_dlogp"]
# actual equations
r, m, H, b = y
e = utils.interp_in_logspace(h, hs, es)
p = utils.interp_in_logspace(h, hs, ps)
dedp = e / p * jnp.interp(h, hs, dloge_dlogps)
# evalute the sigma and dsigmadp
sigma = _sigma_func(
p,
e,
m,
r,
eos["lambda_BL"],
eos["lambda_DY"],
eos["lambda_HB"],
eos["gamma"],
eos["alpha"],
eos["beta"],
)
dsigmadp_fn = jax.grad(_sigma_func, argnums=0) # Gradient w.r.t. p
dsigmade_fn = jax.grad(_sigma_func, argnums=1) # Gradient w.r.t. e
dsigmadp = dsigmadp_fn(
p,
e,
m,
r,
eos["lambda_BL"],
eos["lambda_DY"],
eos["lambda_HB"],
eos["gamma"],
eos["alpha"],
eos["beta"],
)
dsigmadp += dedp * dsigmade_fn(
p,
e,
m,
r,
eos["lambda_BL"],
eos["lambda_DY"],
eos["lambda_HB"],
eos["gamma"],
eos["alpha"],
eos["beta"],
)
A = 1.0 / (1.0 - 2.0 * m / r)
# terms for radius and mass
dpdr = -(e + p) * (m + 4.0 * jnp.pi * r * r * r * p) / r / r * A
dpdr += -2.0 * sigma / r # adding non-GR contribution
# terms for tidal deformability
C1 = 2.0 / r + A * (2.0 * m / (r * r) + 4.0 * jnp.pi * r * (p - e))
C0 = A * (
-6 / (r * r)
+ 4.0 * jnp.pi * (e + p) * (1.0 + dedp) / (1.0 - dsigmadp)
+ 4.0 * jnp.pi * (4.0 * e + 8.0 * p)
+ 16.0 * jnp.pi * sigma
) - jnp.power(2.0 * (m + 4.0 * jnp.pi * r * r * r * p) / (r * (r - 2.0 * m)), 2.0)
drdh = (e + p) / dpdr
dmdh = 4.0 * jnp.pi * r * r * e * drdh
dHdh = b * drdh
dbdh = -(C0 * H + C1 * b) * drdh
dydt = drdh, dmdh, dHdh, dbdh
return dydt
def _calc_k2(R, M, H, b):
r"""
Calculate the second Love number k₂ for tidal deformability.
Parameters
----------
R : float
Neutron star radius [geometric units]
M : float
Neutron star mass [geometric units]
H : float
Auxiliary tidal variable at surface
b : float
Auxiliary tidal variable β at surface
Returns
-------
float
Second Love number k₂
"""
y = R * b / H
C = M / R
num = (
(8.0 / 5.0)
* jnp.power(1 - 2 * C, 2.0)
* jnp.power(C, 5.0)
* (2 * C * (y - 1) - y + 2)
)
den = (
2
* C
* (
4 * (y + 1) * jnp.power(C, 4)
+ (6 * y - 4) * jnp.power(C, 3)
+ (26 - 22 * y) * C * C
+ 3 * (5 * y - 8) * C
- 3 * y
+ 6
)
)
den -= (
3
* jnp.power(1 - 2 * C, 2)
* (2 * C * (y - 1) - y + 2)
* jnp.log(1.0 / (1 - 2 * C))
)
return num / den
[docs]
class AnisotropyTOVSolver(TOVSolverBase):
"""
Post-TOV solver with phenomenological beyond-GR corrections.
Solves post-TOV equations with correction term sigma:
.. math::
\\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}
The sigma function includes:
- Bowers-and-Liang corrections (lambda_BL)
- Horvat et al. corrections (lambda_DY)
- Cosenza et al. corrections (lambda_HB)
- Post-Newtonian corrections (gamma, alpha, beta)
"""
[docs]
def solve(
self, eos_data: EOSData, pc: float, tov_params: dict[str, float]
) -> TOVSolution:
r"""
Solve post-TOV equations for given central pressure.
This function integrates the post-TOV equations that include beyond-GR
corrections. The integration procedure is identical to the standard TOV case,
but the differential equations include additional sigma terms.
Args:
eos_data: EOS quantities in geometric units
pc: Central pressure [geometric units]
tov_params: Beyond-GR coupling parameters, as returned by
:meth:`~jesterTOV.tov.base.TOVSolverBase.fetch_params`.
Must contain all keys listed in :meth:`get_required_parameters`.
Returns:
TOVSolution: Mass, radius, and Love number in geometric units.
Returns NaN values on solver failure (JAX-compatible).
Notes:
The modifications affect the stellar structure but the same integration
method and boundary conditions as the standard TOV case are used.
"""
lambda_BL = tov_params["lambda_BL"]
lambda_DY = tov_params["lambda_DY"]
lambda_HB = tov_params["lambda_HB"]
gamma = tov_params["gamma"]
alpha = tov_params["alpha"]
beta = tov_params["beta"]
# Convert EOSData to dictionary for ODE solver
eos_dict = {
"p": eos_data.ps,
"h": eos_data.hs,
"e": eos_data.es,
"dloge_dlogp": eos_data.dloge_dlogps,
# Add modification parameters
"lambda_BL": lambda_BL,
"lambda_DY": lambda_DY,
"lambda_HB": lambda_HB,
"gamma": gamma,
"alpha": alpha,
"beta": beta,
}
# Extract EOS arrays
ps = eos_data.ps
hs = eos_data.hs
es = eos_data.es
dloge_dlogps = eos_data.dloge_dlogps
# Central values and initial conditions
hc = utils.interp_in_logspace(pc, ps, hs)
ec = utils.interp_in_logspace(hc, hs, es)
dedp_c = ec / pc * jnp.interp(hc, hs, dloge_dlogps)
dhdp_c = 1.0 / (ec + pc)
dedh_c = dedp_c / dhdp_c
# Initial values using series expansion near center
dh = -1e-3 * hc
h0 = hc + dh
r0 = jnp.sqrt(3.0 * (-dh) / 2.0 / jnp.pi / (ec + 3.0 * pc))
r0 *= 1.0 - 0.25 * (ec - 3.0 * pc - 0.6 * dedh_c) * (-dh) / (ec + 3.0 * pc)
m0 = 4.0 * jnp.pi * ec * jnp.power(r0, 3.0) / 3.0
m0 *= 1.0 - 0.6 * dedh_c * (-dh) / ec
H0 = r0 * r0
b0 = 2.0 * r0
y0 = (r0, m0, H0, b0)
sol = diffeqsolve(
ODETerm(_tov_ode),
Dopri5(scan_kind="bounded"),
t0=h0,
t1=0,
dt0=dh,
y0=y0,
args=eos_dict,
saveat=SaveAt(t1=True),
stepsize_controller=PIDController(rtol=1e-5, atol=1e-6),
throw=False,
)
# Extract solution values (throw=False guarantees ys is populated)
R = sol.ys[0][-1] # type: ignore[index]
M = sol.ys[1][-1] # type: ignore[index]
H = sol.ys[2][-1] # type: ignore[index]
b = sol.ys[3][-1] # type: ignore[index]
k2 = _calc_k2(R, M, H, b)
return TOVSolution(M=M, R=R, k2=k2)
[docs]
def get_required_parameters(self) -> list[str]:
"""
Post-TOV requires 6 additional theory parameters.
Returns:
list[str]: ["lambda_BL", "lambda_DY", "lambda_HB", "gamma", "alpha", "beta"]
"""
return ["lambda_BL", "lambda_DY", "lambda_HB", "gamma", "alpha", "beta"]
