Equations of State and Microphysics Interfaces

An equation-of-state interface is part of the numerical algorithm, not a passive lookup. Relativistic simulations may query pressure, energy, entropy, chemical potentials, sound speed, and derivatives as functions of density, temperature, and composition. Interpolation must preserve thermodynamic consistency closely enough for primitive recovery and wave speeds, and every out-of-range query needs an explicit policy. A table can be smooth in values yet inconsistent in derivatives, producing entropy drift or superluminal characteristics.

Helpful background. Formulations of General Relativity supplies stress-energy coupling, Ordinary Differential Equations supplies interpolation and nonlinear inversion concepts, Relativistic Matter and Compact Objects supplies microphysical closure, Computational Relativity supplies evidence standards, and Fluid and Kinetic Theory supplies thermodynamic identities and characteristic speeds.

Define one canonical state interface

Choose independent variables, commonly rest-mass density $\rho$, temperature $T$, and electron fraction $Y_e$. Record units, logarithmic transforms, baryon-mass convention, rest-mass subtraction, and whether tabulated energy is specific, volumetric, or per baryon. A canonical query returns values and declared derivatives from one table version and interpolation scheme.

For specific entropy $s$ and specific internal energy $\epsilon$, define the composition conjugate $\mu_Y=(\partial\epsilon/\partial Y_e)_{s,\rho}$. The first law is then

$$d\epsilon=T\,ds-p\,d\left(\frac{1}{\rho}\right)+\mu_Y\,dY_e,$$

where $\epsilon$, $s$, and $\mu_Y$ are specific quantities per unit rest mass, while $p$ is pressure and $T$ is temperature. A table that instead stores energy, entropy, or chemical potentials per baryon must convert all terms with the same baryon-mass convention; $\mu_Y$ is not automatically the electron chemical potential supplied by a microphysical table. Thermodynamic fields are not independent; deriving all outputs from an interpolated potential can improve consistency, while independently interpolating each column may violate Maxwell relations. CompOSE Core Team 2022, §§3.6, 4.1.1–4.1.2, 4.2.2, 7 and Appendix A.1 documents table content, interpolation, and thermodynamic conventions.

Derivatives control dynamics

Primitive recovery and Riemann solvers need derivatives such as $\partial p/\partial\rho$ and $\partial p/\partial\epsilon$ under stated held-fixed variables. Compute the relativistic adiabatic sound speed from a thermodynamically consistent combination and require

$$0\leq c_s^2\leq 1$$

in units with $c=1$, except where the input model itself violates the physical restriction and is rejected. Finite-difference the returned pressure on an independent dense sample and compare with derivative queries. Repeat near table edges and phase transitions.

High-order unconstrained interpolation can overshoot even when all nodes are physical. Monotone methods reduce overshoot but may lower smoothness. Swesty 1996, §§II–IV develops thermodynamically consistent interpolation and shows why matching tabulated values alone is insufficient.

The sound-speed calculation can be derived from the same first law. At fixed entropy and composition, $d\epsilon=p\,d\rho/\rho^2$. Define

$$\chi=\left(\frac{\partial p}{\partial\rho}\right)_{\epsilon,Y_e}, \qquad \kappa=\left(\frac{\partial p}{\partial\epsilon}\right)_{\rho,Y_e}, \qquad h=1+\epsilon+\frac p\rho.$$

Since the total rest-frame energy density is $e=\rho(1+\epsilon)$ and $de=h\,d\rho$ along that adiabat,

$$c_s^2=\left(\frac{\partial p}{\partial e}\right)_{s,Y_e} =\frac{\chi+p\kappa/\rho^2}{h}.$$

For $p=(\Gamma-1)\rho\epsilon$, this reduces to $c_s^2=\Gamma p/(\rho h)$. It tends to the Newtonian $\Gamma p/\rho$ when $\epsilon,p/\rho\ll1$, and to $\Gamma-1$ when thermal energy dominates. This explains the benchmark’s $2/3$ upper limit for $\Gamma=5/3$. The held-fixed composition assumption describes sound waves faster than weak-reaction equilibration; an equilibrium-composition sound speed is a different derivative.

At a first-order phase boundary, derivatives may be discontinuous. The interface must state whether the table contains a Maxwell or Gibbs construction and how interpolation respects it. Do not claim spectral convergence of fluid variables through a genuine non-smooth phase surface.

Benchmarks and range policy

Use an analytic $\Gamma=5/3$ ideal gas in units $c=k_B=m_b=1$:

$$p=\rho T, \qquad \epsilon=\frac{T}{\Gamma-1}, \qquad s=\frac{\ln(T/T_*)}{\Gamma-1}-\ln(\rho/\rho_*), \qquad \mu_Y=0.$$

Here $T_*$ and $\rho_*$ are fixed reference units, chosen to have numerical value one in this benchmark; the logarithms have dimensionless arguments and an additive entropy constant is immaterial. Tabulate this reference on numerical $\rho,T\in[10^{-6},10^2]$ with $Y_e\in\{0.1,0.3,0.5\}$, using $17,33,65,129$ equally spaced nodes in each of $\log_{10}(\rho/\rho_*)$ and $\log_{10}(T/T_*)$. Use the production interpolation and derivative path without analytic shortcuts. There is no spatial boundary; the closed table rectangle is the domain, and a query even one floating-point step outside any face must return a flagged range error rather than clamp or extrapolate. At each refinement, query the $3\times3$ Gauss–Legendre tensor points inside every table cell and all three composition slices. The analytic formulas above and their derivatives are the independent reference.

For each nonzero returned field or first derivative define a maximum relative error using the analytic magnitude as denominator. For entropy, which can cross zero, use $\max(1,|s_{\rm ex}|)$; for analytically zero composition derivatives and $\mu_Y$, normalize the absolute error by $1+\epsilon+p/\rho$. With $q=(\ln\rho,\ln T,Y_e)$, evaluate the differential first-law residual componentwise as

$$r_k=\partial_{q_k}\epsilon -T\partial_{q_k}s +p\partial_{q_k}\!\left(\frac{1}{\rho}\right) -\mu_Y\partial_{q_k}Y_e,$$

and normalize $|r_k|$ by the sum of the absolute magnitudes of its four displayed terms, using one only if that sum vanishes. Accept the $129^2\times3$ table only if the maximum value and derivative errors are at most $10^{-6}$, the normalized first-law residual is at most $10^{-8}$, the sound-speed relative error is at most $10^{-6}$, and every sample satisfies $0\leq c_s^2\lt 2/3$. Successive table refinements must reduce the maximum value error until these thresholds are reached, and all deliberately out-of-range queries must be flagged.

Observed discrepancy: not available—specification only.

Also test a Maxwell relation and round trips between $(\rho,T,Y_e)$ and alternative variables used by primitive recovery. Those checks should use derivatives assembled independently from the table’s returned derivative columns.

Then test the versioned production table using checksums and the provider’s metadata. Sample the full allowed range plus targeted neighborhoods of boundaries and phase transitions. Compare an independent reader or interpolation implementation. Couple the interface to a smooth sound wave, shock tube, and TOV star; track convergence, entropy, and stationarity.

Every out-of-range query is classified. Options include terminate, apply a documented analytic extension, clamp for a narrowly defined atmosphere state, or invoke a flagged fallback. Silent extrapolation is never acceptable. Count events, their locations, distance beyond the table, and changes to conserved quantities. The validity domain is exactly the table and declared extension region.

Archive representative raw table rows and the canonical unit-conversion tests with the results. A checksum proves file identity but not that the reader interpreted column order, logarithm base, missing-value markers, or composition convention correctly; hand-calculated queries test that semantic layer.

Environment and failure criteria

Pin table name, release, checksum, license record, units, independent variables, interpolation and derivative method, cache or accelerator representation, range policy, precision, analytic benchmark grid, and code revision. A clean interface could be eos-check --config compose-table.yaml --output run/, emitting query samples, value and derivative residuals, identity tests, sound speeds, range events, refinement sweeps, and environment metadata. This specifies reproducibility without claiming an implementation ran here.

Reject an interface with inconsistent units, failed first-law or Maxwell residuals beyond interpolation error, negative or superluminal sound speed, branch-crossing interpolation, or unreported clamping. Failure modes include subtracting rest mass twice, holding the wrong variable fixed in a derivative, cubic overshoot, stale cached table data, and a primitive solver requesting a thermodynamic region the table does not cover.

Diagnostic

Pressure values from two interpolators agree to $10^{-8}$ at all sampled nodes, but one gives negative $c_s^2$ between nodes near a phase transition and causes primitive-recovery failures. Are the interpolators equivalent?

Solution

No. Node agreement does not test derivatives or between-node thermodynamic consistency. Evaluate independent derivative and first-law residuals on a dense neighborhood, use an interpolation compatible with the phase construction, and propagate the difference through sound-wave and primitive-recovery benchmarks. Reject any region with unphysical sound speed unless the input EOS itself explicitly supplies and justifies it.

References

• CompOSE Core Team (S. Typel, M. Oertel, T. Klähn, et al.) (2022). CompOSE Reference Manual, Version 3.01. arXiv:2203.03209. arXiv.
• Swesty, F. D. (1996). “Thermodynamically consistent interpolation for equation of state tables.” Journal of Computational Physics 127, 118–127. DOI.