Constraints Preserving Lax-Wendroff Flux Reconstruction for Relativistic Hydrodynamics with General Equations of State

In the realm of relativistic astrophysics, the ideal equation of state with a constant adiabatic index provides a poor approximation due to its inconsistency with relativistic kinetic theory. However, it is a common practice to use it for relativistic fluid flow equations due to its simplicity. Here we develop a high-order Lax-Wendroff flux reconstruction method on Cartesian grids for solving relativistic hydrodynamics equations with several general equations of state available in the literature. We also study the conversion from conservative to primitive variables, which depends on the equation of state in use, and provide an alternative method of conversion when the existing approach does not succeed. For the admissibility of the solution, we blend the high-order method with a low-order method on sub-cells and prove its physical admissible property in the case of all the equations of state used here. Lastly, we validate the scheme by several test cases having strong discontinuities, large Lorentz factor, and low density or pressure in one and two dimensions.

💡 Research Summary

The paper addresses the challenge of developing a high‑order, constraint‑preserving numerical scheme for the relativistic hydrodynamics (RHD) equations when general equations of state (EOS) are employed. While the ideal EOS with a constant adiabatic index is widely used for its simplicity, it is known to be inconsistent with relativistic kinetic theory and therefore unsuitable for many astrophysical flows. The authors therefore consider three more realistic EOSs—TM‑EOS, IP‑EOS, and RC‑EOS—together with the ideal EOS (ID‑EOS) and aim to construct a method that works robustly for all of them.

A central difficulty lies in the conversion from conserved variables (density D, momentum m, total energy E) to primitive variables (rest‑mass density ρ, velocity v, pressure p). For the RC‑EOS the standard iterative approach fails to converge to the physically admissible root. The authors introduce a novel formulation based on the auxiliary variable Π = E + p. By expressing the EOS as p = ρ T(h) and using the relations Π = ρ h Γ² and h = p Π² − |m|² /D, they derive a single scalar nonlinear equation S(Π)=0. They prove that S(Π) is strictly increasing for Π ≥ E, that S(E) < 0 and S(Π) → +∞ as Π → +∞, guaranteeing a unique root Π* ∈ (E,∞). A Newton‑Raphson iteration with Π₀ = E is shown to produce monotone increasing iterates that converge to Π* in typically 4–5 steps, even for extreme low‑density, low‑pressure, or near‑light‑speed states. This provides a reliable, EOS‑independent primitive‑variable recovery algorithm.

For time integration the authors avoid multi‑stage Runge‑Kutta or multi‑step schemes, which require substantial inter‑processor communication and storage of several past solution stages. Instead they adopt a single‑step Lax‑Wendroff (LW) procedure, which yields arbitrary‑order temporal accuracy in one evaluation. They employ the Jacobian‑free Lax‑Wendroff Flux Reconstruction (LW‑FR) framework, originally introduced for non‑relativistic flows, and adapt it to the relativistic setting. The flux reconstruction builds a high‑order polynomial representation of the flux inside each cell using correction functions and numerical fluxes at cell interfaces, making the method well‑suited for modern vector processors and GPU architectures.

High‑order schemes alone cannot guarantee that the numerical solution stays within the physically admissible set
U_ad = { D > 0, ρ > 0, p > 0, |v| < 1 }. To enforce this, the authors construct a low‑order, constraint‑preserving scheme on sub‑cells. By reformulating the admissibility conditions as D > 0 and q := E − p + |m|² > 0, both of which are concave functions of the conserved variables, the low‑order method can be proved to keep the solution inside U_ad for any EOS considered. The final algorithm blends the high‑order LW‑FR solution with the low‑order admissible solution using a convex combination determined by a local scaling factor that detects non‑admissible flux arguments. This blending simultaneously suppresses Gibbs oscillations near discontinuities and restores physical admissibility.

The authors validate the method on a suite of one‑ and two‑dimensional test problems: relativistic shock tubes, blast waves with Lorentz factors up to ~100, interactions involving very low rest‑mass density or pressure, and multi‑dimensional jet‑like configurations. For each test, all four EOSs are exercised. The results demonstrate the expected order of convergence in smooth regions, sharp resolution of discontinuities without spurious oscillations, and robust primitive‑variable recovery even for the RC‑EOS where previous methods failed. The blending parameter remains modest, confirming that the high‑order component dominates except in the immediate vicinity of strong shocks.

In summary, the paper makes three major contributions: (1) a unified, provably convergent primitive‑variable recovery algorithm for a class of realistic relativistic EOSs, (2) a single‑step, high‑order Lax‑Wendroff flux‑reconstruction scheme that is computationally efficient and well‑suited to modern parallel hardware, and (3) a mathematically rigorous blending strategy that guarantees the solution stays within the admissible physical set for all considered EOSs. This work advances the state of the art in relativistic hydrodynamics simulations, offering a practical tool for high‑performance astrophysical modeling where both accuracy and physical fidelity are essential.