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

Motivation

• Physically admissible solutions require positive density and pressure and bounded material speeds, but preserving these constraints is difficult for high-order schemes.Highly diffusive remedies require fine grids and small time steps, motivating high-order constraints-preserving methods.Source: §1 · p. 2
• General equations of state are needed because the commonly used ideal equation of state poorly approximates relativistic flows, especially semi-relativistic and two-component fluids.The ideal equation of state is simple but derived from non-relativistic thermodynamics.Source: §1 · p. 2

Method

• The method uses a high-order Lax-Wendroff flux reconstruction approach that requires one time step for arbitrary-order accuracy, unlike Runge-Kutta methods.The approach is presented as an alternative for reducing data-sharing and memory costs associated with common time integrators.Source: §1 · p. 2

Contributions

• A new conservative-to-primitive conversion method is proposed for the equation of state in [50], supporting simulations with high Lorentz factors.The method addresses cases where the existing conversion approach does not succeed.Source: §1 · p. 2
• The scheme blends the high-order method with a constraints-preserving low-order method on sub-cells to ensure admissibility and control Gibbs oscillations.The low-order scheme’s constraints-preserving property is rigorously proved.Source: §1 · p. 2

Validation

• One- and two-dimensional tests assess accuracy and effectiveness for strong discontinuities, large Lorentz factors, and low density or pressure.The paper organizes validation across multiple equations of state and challenging relativistic flow regimes.Source: §1 · pp. 2–3

Governing equations

• Relativistic hydrodynamics requires an equation of state to recover primitive variables from conservative variables and close the conservation-law system.The conservative variables include energy and momentum quantities, while the primitive variables include density, velocity, and pressure.Source: §2 · p. 3

Equation of state

• The ideal equation of state is a poor relativistic approximation because its constant adiabatic index cannot represent the required variation across flow regimes.The paper therefore considers TM-EOS, IP-EOS, and RC-EOS alongside ID-EOS; the exact relativistic gas law is more expensive because it uses modified Bessel functions.Source: §2.1 · p. 4

Admissible set

• Physical admissibility requires positive density, pressure, and internal energy together with subluminal velocity, and leaving this set can violate hyperbolicity and cause computational failure.The admissible set is defined in terms of primitive variables, making direct verification computationally expensive.Source: §2.2 · p. 5
• An equivalent admissible-set representation uses conservative-variable constraint functions D(u) and q(u), which are directly computable and concave for the equations of state considered.This representation supports efficient admissibility checking and limiting without first recovering all primitive variables.Source: §2.2 · p. 5

Primitive-variable recovery

• For RC-EOS, the existing iterative recovery of Γ can converge incorrectly, so the paper instead solves S(Π)=0 for Π using Newton-Raphson iteration.S(Π) has a unique root in (E, ∞), and iterations initialized at Π0=E increase toward that root.Source: §2.2 · pp. 5–6
• The proposed RC-EOS recovery typically converges in four to five iterations even for low density, low pressure, and speeds near unity, although the authors cannot algebraically prove every iterate remains above E.A damped Newton-Raphson method is suggested as a safeguard.Source: §2.2 · p. 6

3 Numerical scheme

• The Lax-Wendroff flux reconstruction scheme evolves Gauss-Legendre solution-point degrees of freedom using a polynomial approximation and time-averaged fluxes.Flux reconstruction uses correction functions to form a continuous approximation from locally approximated time-average fluxes.Source: §3 · p. 7
• Finite-difference formulas approximate temporal flux derivatives instead of replacing them with spatial derivatives through the chain rule.This avoids the large computational cost associated with the chain-rule approach.Source: §3 · p. 7

3.1 Time average flux

• Non-admissible intermediate flux arguments are scaled into the admissible region before primitive-variable conversion and flux evaluation.The scaling is repeated across admissibility constraints and derivative orders while preserving admissible cell averages through convexity and concavity arguments.Source: §3.1 · p. 9

3.2 Blending of the scheme

• The high-order update is blended with a low-order sub-cell scheme to control Gibbs oscillations and preserve admissibility near discontinuities.Each element is divided into N + 1 sub-elements, and the low-order update uses a first-order finite-volume method with Rusanov flux.Source: §3.2 · pp. 10–11
• The first-order Rusanov update is admissibility-preserving for ID-EOS, TM-EOS, IP-EOS, and RC-EOS under CFL-type restrictions, assuming the previous solution is admissible.At inter-element faces, extremal solution points require a separate treatment because the blended inter-element flux can violate the low-order update’s admissibility property.Source: §3.2 · p. 11

4 Numerical simulations

• Numerical simulations demonstrate the scheme’s robustness and capability to capture different wave structures.The simulations use RHDTenkai.jl, developed with Tenkai.jl as a library.Source: §4 · p. 12

4.1.1 Accuracy test

• The scheme converges at order O(∆x)^(N+1) for N = 3 and N = 4 across all tested equations of state.Source: §4.1.1 · p. 13

4.1.2 1-D Riemann problem 1

• It effectively captures shock waves, contact discontinuities, and rarefaction waves in the first Riemann problem across the tested equations of state.The simulations use 500 cells and polynomial degrees N = 3, 4, with reference solutions from exact or fine-mesh Rusanov calculations.Source: §4.1.2 · pp. 13, 15
• ID-EOS with γ = 4/3 more closely approximates TM-EOS, IP-EOS, and RC-EOS than γ = 5/3 in relativistic-temperature regions.The other three equations of state also produce closely aligned solutions, reflecting similar specific-enthalpy distributions.Source: §4.1.2 · pp. 15–16

4.1.3–4.1.4 1-D Riemann problems

• For the second and third Riemann problems, TM-EOS, IP-EOS, and RC-EOS remain close, while ID-EOS with γ = 4/3 better matches them than γ = 5/3.The γ = 5/3 case produces a more elongated rarefaction wave in the second problem because of its higher sound speed.Source: §4.1.3 · p. 16, §4.1.4 · p. 16

4.1.5 1-D density perturbation problem

• The density-perturbation test resolves the smooth sinusoidal profile together with rarefaction, contact, and shock waves for every equation of state.The equivalence of the TM-EOS, IP-EOS, and RC-EOS solution graphs is also observed.Source: §4.1.5 · p. 18

4.1.6 1-D blast wave problem

• The blast-wave test captures all structures for every equation of state, while TM-EOS, IP-EOS, and RC-EOS remain close but differ from ID-EOS in zoomed regions.Source: §4.1.6 · p. 18

• Fourth-order accuracy for N = 3 and fifth-order accuracy for N = 4 are observed across all tested equations of state.Source: §4.2.1 · p. 19

2-D Riemann problem 1

• The first two-dimensional Riemann problem forms a jet-like structure, a mushroom-like lower-left structure, and curved shocks that the scheme captures effectively.The ID-EOS case with γ = 5⁄3 produces faster curved shocks, while the other equations of state give similar solutions.Source: §4.2.2 · p. 22, Figure 9 · p. 22, Figure 10 · p. 23

2-D Riemann problem 4

• The fourth Riemann problem evolves into four rarefaction waves and two symmetric shocks that the scheme captures for all equations of state.Results are compared along a diagonal cut through the domain.Source: §4.2.5 · p. 25, Figure 13 · p. 26, Figure 14 · p. 27

2-D relativistic jet

• The high-speed jet test captures the Mach shock at the beam head and other waves for all equations of state.The test uses a near-light-speed jet, strong relativistic shocks, shear waves, interface instabilities, and an ultra-relativistic region.Source: §4.2.7 · p. 27, Figure 19 · p. 32

2-D bubble shock interaction

• The bubble-shock tests capture the deformed bubble and waves generated during and after interaction for both density cases and all tested equations of state.Reflected waves heat the bubble again after striking the reflective boundaries.Source: §4.2.8 · pp. 27–28, Figure 15 · p. 28, Figure 16 · p. 29

2-D double Mach reflection

• The double Mach reflection produces no Kelvin–Helmholtz instability with ID-EOS γ = 1.4, whereas the other equations of state develop such instabilities.The simulations use 960 × 240 cells, N = 4, and are run to t = 4.Source: §4.2.9 · p. 30

Kelvin–Helmholtz instability

• The Kelvin–Helmholtz test captures small-scale instabilities for all equations of state despite their difficulty for diffusive schemes.The reported simulations use 640 × 320 cells with N = 4.Source: §4.2.10 · pp. 30, 32

• The method combines a high-order scheme with a constraints-preserving low-order subcell scheme to ensure admissibility and control Gibbs oscillations.The low-order scheme's constraints-preserving property is rigorously proved.Source: §5 · p. 32
• The work develops a high-order Lax–Wendroff flux reconstruction scheme for relativistic hydrodynamics with several general equations of state.The admissible-region characterization is changed to concave constraints computable directly from conservative variables.Source: §5 · p. 32
• An alternative pressure-based conservative-to-primitive conversion is proposed because the existing RC-EOS procedure does not meet expectations in some cases.Source: §5 · p. 32
• The accuracy tests and challenging simulations demonstrate robustness for high Lorentz factors, low density or pressure, rarefaction, strong shocks, and other discontinuities.Source: Figure 20 · p. 33, Figure 21 · p. 33, Figure 22 · p. 33

• Animations showing temporal evolutions for some simulations can be viewed separately.Source: Additional data · p. 37

Constraints preserving nature of first-order finite volume scheme

• The proof establishes admissibility for TM-EOS, IP-EOS, and RC-EOS, while the ID-EOS case was proved previously.The argument targets the admissibility constraints of the first-order finite-volume method.Source: A Constraints preserving nature of… · pp. 37–38
• The Rusanov numerical flux and eigenvalue bounds provide the ingredients for the first-order update analysis.The spectral radius satisfies Λ < 1 because the sound speed satisfies cs < 1.Source: A Constraints preserving nature of… · pp. 37–39
• For TM-EOS, IP-EOS, and RC-EOS, the required positivity inequalities follow from p > 0 and, where needed, ρ > 0.The EOS-specific calculations establish positive terms for each of the three general equations of state.Source: A Constraints preserving nature of… · pp. 38–39
• Positive quadrature weights, spectral-radius bounds, and a CFL-type restriction make the update coefficients positive.The resulting update is a convex combination in the convex admissible set U′ad.Source: A Constraints preserving nature of… · p. 40

Proof of S′(Π) > 0 for Π ≥E

• R(h) is decreasing on (1, ∞) and satisfies R(1) = 8/5.These properties imply that R(h) remains greater than one throughout the interval.Source: B Proof of S′(Π) > 0 for Π ≥E · p. 40
• R(h) > 1 on (1, ∞), giving S′(Π) > 0 for Π ≥ E.The positivity conclusion follows directly from the bound on R(h).Source: B Proof of S′(Π) > 0 for Π ≥E · p. 40

Julia code for conversion to primitive variables in the case of RC-EOS

• Algorithm 2 is the conservative-to-primitive conversion procedure for RC-EOS.Source: C Julia code for conversion to… · p. 41
• The algorithm is presented specifically for conversion in the RC-EOS case.Source: C Julia code for conversion to… · p. 41
• The section associates the RC-EOS conversion procedure with executable Julia code.Source: C Julia code for conversion to… · p. 41

Comparison of conservative to primitive conversion for ideal equation of state

• The compared methods solve different scalar formulations for conversion: a pressure equation from [17] versus a velocity polynomial from.Both formulations are solved using iterative root-finding methods.Source: D Comparison of conservative to… · p. 42
• In Example 2, the velocity-function values are approximately equal while the Φ-function values differ significantly, favoring the new method from.The passage reports this discrepancy as evidence that the new conversion is more accurate.Source: D Comparison of conservative to… · p. 43
• In Example 3, Φ(pnew) = −3.557432126655158e −133 and Φ(pold) = 2.3753539135640267e −9, while the Ξ values agree to machine precision.The reported Φ discrepancy again supports the greater accuracy of the new method from.Source: D Comparison of conservative to… · p. 43