Thermodynamic Consistency of Equations of State

A thermodynamic equation of state is more than a collection of pressure and energy values. Its fields must be derivatives of one potential, so the first law, Maxwell identities, stability conditions, and composition constraints agree at every state. A table can interpolate each column smoothly yet violate those relations between grid points. The safest construction interpolates a generating potential and derives dependent quantities consistently, then tests analytic limits and derivative convergence.

Required background. Hilbert Stress–Energy from Metric Variation supplies the energy and pressure entering the material stress tensor.

Helpful background. The Einstein Equation supplies gravitational use of the equation of state; Relativistic Matter and Compact Objects gives the volume context; and Computational Relativity supplies table interpolation and evolution settings.

One potential generates the state

Let $\mathcal F(T,n_i)$ be Helmholtz free-energy density as a function of temperature and independent species densities. Free energy is the useful starting point when temperature and particle densities are the table axes: its slopes tell us the entropy and the energy cost of adding a particle. Its derivatives define

$$s=-\left(\frac{\partial\mathcal F}{\partial T}\right)_{n_i}, \qquad \mu_i=\left(\frac{\partial\mathcal F}{\partial n_i}\right)_{T,n_{j\ne i}},$$ $$p=\sum_i n_i\mu_i-\mathcal F, \qquad \varepsilon=\mathcal F+Ts.$$

Differentiation gives the first law and Gibbs–Duhem relation,

$$d\varepsilon=T ds+\sum_i\mu_i dn_i, \qquad dp=s dT+\sum_i n_i d\mu_i.$$

Mixed derivatives commute wherever $\mathcal F$ is twice differentiable. For example,

$$\left(\frac{\partial s}{\partial n_i}\right)_T =-\left(\frac{\partial\mu_i}{\partial T}\right)_{n_j}.$$

This Maxwell identity is an independent derivative check. Oertel and collaborators summarize the thermodynamic variables and equilibrium constraints in Oertel et al. 2017, §§II.A–II.B.2.

Many astrophysical tables use baryon density $n_B$, temperature $T$, and charge or electron fraction $Y_e$. Then composition derivatives must specify what is held fixed, and chemical potentials may need a constrained Legendre transformation. Beta equilibrium is not a universal table identity: it is imposed only when weak reactions act faster than the dynamics and the relevant neutrino chemical potential is known.

Units and zero points are part of the state definition. A table may report energy per baryon with or without a reference rest mass, mass density rather than baryon density, and entropy per baryon rather than entropy density. Converting one column without converting its conjugate derivative breaks the first law. The electron, photon, or neutrino contribution may be included, excluded, or supplied as a separate module. Before combining tables, reconstruct one test state from the documented conventions and verify pressure, enthalpy, and chemical-potential identities in a common unit system.

Stability and causal derivatives

Local thermal stability requires positive heat capacity at fixed volume and composition,

$$c_V=T\left(\frac{\partial s}{\partial T}\right)_{n_i}\gt0.$$

Mechanical and diffusive stability constrain the Hessian of the free energy with respect to density and composition. Along a one-component isotherm, positive compressibility requires $(\partial p/\partial n)_T\gt0$. A negative eigenvalue of the full susceptibility matrix indicates instability toward a density or composition fluctuation even when the one-dimensional pressure curve rises.

The hydrodynamic sound derivative is adiabatic and follows the dynamical composition constraint,

$$c_s^2=\left(\frac{\partial p}{\partial\varepsilon}\right)_{s/n_B,Y_i}.$$

It is not generally equal to an isothermal derivative or a finite difference taken along a table axis. Stable perfect-fluid matter requires $c_s^2\geq0$, and causal acoustic propagation requires $c_s^2\leq1$. These are not complete causality tests for magnetic, elastic, multifluid, or dissipative closures. A first-order coexistence region is non-smooth and must be treated with its phase construction rather than forcing a high-order derivative through the kink.

For a concrete derivative calculation, let $a=\mathcal F/n$ and hold composition fixed. Use subscripts for partial derivatives. Since entropy per baryon is $\sigma=-a_T$, moving along an isentrope gives

$$\left.\frac{dT}{dn}\right|_\sigma=-\frac{a_{nT}}{a_{TT}}, \qquad c_V=-nT a_{TT}.$$

The pressure derivative along that same path is therefore

$$\left.\frac{dp}{dn}\right|_\sigma =2na_n+n^2a_{nn}-n^2\frac{a_{nT}^2}{a_{TT}}, \qquad c_s^2=\frac{n}{\varepsilon+p} \left.\frac{dp}{dn}\right|_\sigma .$$

The last step uses the first law $d\varepsilon=(\varepsilon+p)\,dn/n$ at fixed $\sigma$ and composition. For a stable thermal state $a_{TT}\lt0$, so the additional thermal term in the adiabatic derivative is nonnegative. It is precisely the temperature rise during compression that an isothermal derivative misses. The formula requires nonzero heat capacity; the $T=0$ slice is treated by its cold EOS.

Analytic reference states

Two limits catch different mistakes. For a dilute classical nonrelativistic gas with free energy per particle

$$a=m+T\left[\ln\!\left(\frac{n}{n_Q(T)}\right)-1\right], \qquad n_Q(T)=g\left(\frac{mT}{2\pi}\right)^{3/2},$$

the derivative $p=n^2(\partial a/\partial n)_T$ gives $p=nT$. Directly, $\sigma=5/2-\ln(n/n_Q)$, so an isentrope has $T\propto n^{2/3}$. Thus $(dp/dn)_\sigma=5T/3$, whereas $(dp/dn)_T=T$. Including rest energy gives $\varepsilon=mn+3nT/2$ and $c_s^2=(5T/3)/(m+5T/2)\simeq5T/(3m)$ in the required $T\ll m$ regime. This benchmark tests thermal and density derivatives together; it must not be extrapolated to relativistic temperatures.

For a cold nonrelativistic ideal Fermi gas,

$$E_F=\frac{(3\pi^2n)^{2/3}}{2m}, \qquad \varepsilon=mn+\frac35nE_F, \qquad p=\frac25nE_F.$$

It checks the $n^{5/3}$ scaling, the zero-temperature limit, and consistency of rest-mass bookkeeping. An ultrarelativistic radiation component supplies a third check, $p=\varepsilon/3$ and $s=(\varepsilon+p)/T$.

Table workflow and convergence

The continuum input is a versioned potential or microphysical model on declared axes $(n_B,T,Y_e)$ with units, particle content, reference masses, and valid domain. The CompOSE format specifies thermodynamic quantities and consistency conventions in CompOSE Core Team 2022, Reference Manual v3.01, §§3.5–3.6, §4.1.2, and Appendix A.1 (PDF). The linked file is labeled version 3.01 on its title page despite the older version number in its filename.

For a numerical table study, record the actual interpreter and library versions, precision, source revision, table release, checksum, axis arrays, unit conversions, interpolation order, and derivative method. Interpolating one potential ensures derivative identities only when the derivatives are actually taken from that interpolant; it does not automatically preserve convexity or causal sound speed. Even a monotone interpolant of values may have unphysical second derivatives. The protocol below is a proposed table test, not a report of an executed production-table calculation.

Evaluate the ideal-gas and degenerate-gas reference tables on nested $32^3$, $64^3$, and $128^3$ grids over the same physical domain. Compare derived $p$, $s$, $\varepsilon$, $c_V$, and $c_s^2$ against analytic values in $L_1$ and $L_\infty$ norms. Report first-law and Maxwell residuals normalized to local thermodynamic scales. Separately refine each axis to identify whether density, temperature, or composition interpolation dominates. Halve finite-difference steps until truncation and roundoff exhibit a stable window.

One failure mode is independent interpolation of $p$ and $\varepsilon$ that violates the generating relation. Others include derivative oscillations from high-order splines, negative entropy or heat capacity, superluminal sound speed created between nodes, inconsistent rest-mass offsets, extrapolation outside the table domain, and differentiation across a phase boundary. A small value error does not guarantee a small derivative error; both must converge.

Exercises

Starting from a one-component free-energy density $\mathcal F(T,n)$, derive $p=n\mu-\mathcal F$ and the Maxwell relation between $s$ and $\mu$.

Solution

Extensivity gives $F=V\mathcal F(T,N/V)$. At fixed $T,N$,

$$p=-\left(\frac{\partial F}{\partial V}\right)_{T,N} =n\frac{\partial\mathcal F}{\partial n}-\mathcal F=n\mu-\mathcal F.$$

Since $s=-\partial_T\mathcal F$ and $\mu=\partial_n\mathcal F$, commuting mixed derivatives gives $(\partial s/\partial n)_T=-(\partial\mu/\partial T)_n$.

References

• CompOSE Core Team. CompOSE Reference Manual, Version 3.01, 1 April 2022. Official PDF.
• Oertel, Micaela, Matthias Hempel, Thomas Klähn, and Stefan Typel. “Equations of State for Supernovae and Compact Stars.” Reviews of Modern Physics 89, 015007 (2017). DOI.