Hierarchical Equation-of-State Reconstruction

Hierarchical EOS reconstruction treats the cold equation of state as a latent function shared by neutron stars while keeping each star’s mass, spin, temperature approximation, geometry, and measurement nuisance parameters distinct. The challenge is to impose thermodynamic stability and causality without allowing the parameterization or its prior to decide the answer in advance—especially when sharp phase transitions are possible. This page compares the principal representations, derives the common hierarchical likelihood, and states a reproducible validation route. Dated inputs use an evidence cutoff of 2026-08-18.

Required background. Equivalence principles and gravitational measurement supplies the separation between invariant stellar parameters and messenger-specific observables.

Helpful background. Stress–energy and conservation supplies the EOS variables; the Einstein equation supplies the stellar map; relativistic matter and compact objects develops equilibrium and tidal response; gravitational radiation and waveforms supplies tidal likelihoods; and computational relativity supplies solver, emulator, and sampling validation.

Thermodynamic map from a latent EOS to a star

For cold catalyzed matter, a barotrope may be represented by pressure $p(\varepsilon)$ as a function of total energy density. For the corresponding single-fluid barotropic model, mechanical stability and causal acoustic propagation require, away from a discontinuity,

$$0\le c_s^2=\frac{dp}{d\varepsilon}\le1.$$

This derivative follows the equilibrium composition. Oscillations faster than chemical equilibration use a frozen-composition derivative instead; passing the barotropic bound does not verify every mode of a reacting or multifluid material. Thermodynamic consistency also fixes baryon density. At zero temperature,

$$d\varepsilon=\frac{\varepsilon+p}{n}\,dn, \qquad \frac{dn}{n}=\frac{d\varepsilon}{\varepsilon+p}.$$

Thus a pressure curve and one matching value of $n$ determine chemical potential $\mu=(\varepsilon+p)/n$ and baryon density along a phase. A sharp first-order transition at pressure $p_t$ has an energy-density jump $\Delta\varepsilon$. For a Maxwell construction, equality of chemical potentials gives

$$\frac{n_+}{n_-} =\frac{\varepsilon_++p_t}{\varepsilon_-+p_t}.$$

The two endpoint densities belong to different phases at the same pressure; they are not a smooth stellar layer with an arbitrarily small sound speed. Mixed-phase constructions require their own surface-tension and charge-neutrality assumptions. The phase-transition construction explains the distinction.

Given central pressure and a low-density crust match, the Tolman–Oppenheimer–Volkoff and tidal perturbation equations define a deterministic map

$$\mathcal F_{\rm star}:p(\varepsilon) \longmapsto \{M(R),\Lambda(M),M_{\max},\ldots\}.$$

The map is used only on radially stable branches unless the likelihood explicitly models metastable or disconnected families. Numerical interpolation must preserve monotonicity and the first law; a smooth pressure interpolant with noisy derivatives can create spurious sound-speed excursions.

What different EOS representations permit

Representation	Native variables	Strength	Characteristic failure
Piecewise polytrope	$p=K_i\rho^{\Gamma_i}$ between density knots	Fast and interpretable; broad stiffness changes with enough segments	Fixed knots and few segments smear a sharp transition; priors on $\Gamma_i$ induce nonuniform pressure priors
Spectral adiabatic index	Smooth coefficients for $\Gamma(p)=d\ln p/d\ln n$	Efficient approximation to smooth causal barotropes	Truncated smooth basis can ring or exclude discontinuities
Sound-speed segments	$c_s^2(\varepsilon)$ plus transition jumps	Stability and causality are local; peaks and soft regions are explicit	Knot and smoothing priors can manufacture preferred feature widths
Gaussian process	Correlated prior over transformed $p(\varepsilon)$ or $c_s^2$	Flexible, with tunable correlation structure and training information	Kernel, mean, training ensemble, and conditioning can still suppress out-of-family behavior
Explicit phase-transition model	$p_t$, $\Delta\varepsilon$, and high-density phase parameters	Directly represents latent heat and disconnected branches	Strong Occam and label effects; branch assignment may be multimodal

Read et al. 2009, §§ III–V and Appendix C, Table III constructs the widely used piecewise-polytropic map and quantifies its stellar approximation error. Lindblom 2010, §§ II–IV develops a spectral representation of the adiabatic index. Landry and Essick 2019, §§ II–III provides a nonparametric Gaussian-process route. These are coordinate systems on function space, not neutral alternatives: “uniform coefficients” means a different prior in pressure, radius, maximum mass, and phase structure for each.

Every analysis should therefore show prior-predictive bands for $p(\varepsilon)$, $c_s^2$, $R(M)$, $\Lambda(M)$, $M_{\max}$, and the frequency and size of phase transitions. Legred et al. 2021, §§ II–IV provides a concrete nonparametric example of propagating an added pulsar radius likelihood through high-density quantities. If the posterior resembles a narrow prior-predictive band in the density interval of interest, the data have not independently established that feature.

The hierarchical likelihood

Let $\phi$ describe the common EOS and $\xi_i$ the latent parameters of object or event $i$. If $d_i$ are conditionally independent datasets after shared calibration and theory nuisances $\zeta$ have been exposed, the posterior is schematically

$$p(\phi,\zeta,\{\xi_i\}\mid\{d_i\}) \propto \pi(\phi,\zeta) \prod_i \mathcal L_i\!\left[d_i\mid \mathcal F_i(\phi,\xi_i,\zeta)\right] \pi_i(\xi_i\mid\phi,\zeta).$$

For a population selected by detection and conditioned on its observed size $N$, its factor additionally contains the appropriate selection normalization $\alpha(\phi,\lambda)^{-N}$. An analysis that also infers the event rate uses the corresponding Poisson count term rather than treating this conditional factor as the whole likelihood. Targeted NICER observations, radio timing samples, and gravitational-wave searches have different selection processes; one denominator cannot be copied across them. A radio mass likelihood used inside a NICER analysis must not be multiplied a second time in the joint EOS fit.

Recovering likelihood information from posterior samples

Suppose an event analysis used an interim prior $\pi_0(\theta)$ and returned samples from $q(\theta)=\mathcal L(d\mid\theta)\pi_0(\theta)/Z_0$. To evaluate an alternative, absolutely continuous population density $\pi(\theta\mid\phi)$, integrate

$$\begin{aligned} \mathcal L_{\rm pop}(\phi) &=\int \mathcal L(d\mid\theta)\pi(\theta\mid\phi)\,d\theta\\ &=Z_0\,\mathbb E_q\!\left[ \frac{\pi(\theta\mid\phi)}{\pi_0(\theta)}\right]. \end{aligned}$$

The evidence $Z_0$ is constant in $\phi$, so posterior samples can estimate the average weight. They cannot restore regions where $\pi_0=0$, and a few dominant weights indicate poor overlap. The diagnostic $N_{\rm eff}=(\sum_s w_s)^2/\sum_s w_s^2$ measures this weight concentration, before accounting for sample autocorrelation.

A deterministic EOS curve in a continuous $M$–$R$ space needs extra care: it has zero area under an ordinary two-dimensional posterior. One evaluates a reconstructed likelihood along the curve and integrates over a declared central-pressure or mass measure, including branch weights. Testing whether samples land exactly on the curve is not a valid reweighting algorithm. A finite model-discrepancy width is an alternative only when that width is physically justified.

The latent map can be multimodal. The same measured mass may lie on different stable branches when a strong transition produces twins. Sampling only one central-density interval or ordering stars by central density can silently discard a branch. Branch labels should be marginalized, and the prior probability of forming or observing each branch must not be inferred from equilibrium existence alone.

Model uncertainty includes crust matching, phase construction, gravity theory, stellar rotation, tidal waveform, X-ray atmosphere and spot geometry, nuclear-theory covariance, and population distributions. Some nuisances are correlated across objects: an X-ray effective-area calibration may be target-specific or shared by epoch, while one chiral interaction ensemble affects the entire low-density EOS. The conditional-independence graph should determine which nuisance is copied and which is shared.

Dated datasets, calibration, and reproduction

At the cutoff, a representative joint dataset names Australia Telescope National Facility 2026, ATNF v2.8.1 release-history row and the primary timing likelihoods for heavy pulsars; target-specific NICER/HEASARC ObsID lists and exact CALDB/response, screening and background products for an explicitly declared subset of the five peer-reviewed targets J0030+0451, J0740+6620, J0437−4715, J1231−1411, and J0614−3329; the GWOSC GW170817 event release and current GWTC-5.0 catalog context with event-level strain, calibration and selection metadata; and the final PREX-II/CREX publications plus a versioned chiral-EFT calculation. The cutoff-date J2124−3358 preprint is a separate optional input, not part of that peer-reviewed set; its atmosphere and prior sensitivity are recorded in the cutoff-date NICER target inventory. “NICER radius,” “GW tides,” or “nuclear prior” is not a dataset identifier.

Selections include radio survey and timing follow-up, X-ray good-time and detector cuts, gravitational-wave search threshold and data-quality segments, and nuclear kinematic cuts. Likelihood versions, interim priors, population assumptions, waveform/atmosphere alternatives, response calibration, and correlated systematic uncertainty are archived alongside the input products.

A bounded reproduction pins the TOV/Love solver, EOS interpolation, crust table, parameterization code, likelihood adapters, sampler, dependency environment, random seeds, and all data-product checksums. It verifies a reference polytrope, mass–radius and Love convergence under tighter tolerances, thermodynamic identities, stable-branch detection, and identical likelihood normalization on an analytic Gaussian toy problem. Then it performs synthetic-data injection and recovery for both a smooth EOS and one with a density jump, repeats the result under at least two representation families, checks effective sample size after reweighting, and reports prior-to-posterior information by density.

This page documents the method and does not claim a new joint inference. A numerical constraint is accepted only when it survives parameterization, prior, calibration, messenger-removal, and solver-resolution alternatives within a declared domain.

Check your understanding

On one EOS segment let $p=p_0+s(\varepsilon-\varepsilon_0)$ with constant $s$. Which values of $s$ satisfy local stability and causality? Why does $s=0$ require separate bookkeeping in an inverse map $\varepsilon(p)$?

Solution

Here $dp/d\varepsilon=s$, so local mechanical stability and relativistic causality give

$$0\le s\le1.$$

For $s>0$, the segment is invertible:

$$\varepsilon(p)=\varepsilon_0+\frac{p-p_0}{s}.$$

At $s=0$, pressure is constant while energy density may jump. Division by $s$ is invalid; the transition must instead record the two endpoint densities or latent heat $\Delta\varepsilon$ at $p_t$. Replacing it by an arbitrarily small positive $s$ changes the stellar model and can hide a disconnected stable branch.

References

• Australia Telescope National Facility. “ATNF Pulsar Catalogue.” PSRCAT version 2.8.1, released 18 June 2026. Official version history.
• Landry, Philippe, and Reed Essick. “Nonparametric Inference of the Neutron Star Equation of State from Gravitational Wave Observations.” Physical Review D 99 (2019): 084049. DOI.
• Legred, Isaac, Katerina Chatziioannou, Reed Essick, Sophia Han, and Philippe Landry. “Impact of the PSR J0740+6620 Radius Constraint on the Properties of High-Density Matter.” Physical Review D 104 (2021): 063003. DOI.
• Lindblom, Lee. “Spectral Representations of Neutron-Star Equations of State.” Physical Review D 82 (2010): 103011. DOI.
• Read, Jocelyn S., Benjamin D. Lackey, Benjamin J. Owen, and John L. Friedman. “Constraints on a Phenomenologically Parametrized Neutron-Star Equation of State.” Physical Review D 79 (2009): 124032. DOI.