Conservative Relativistic Hydrodynamics

Conservative relativistic hydrodynamics rewrites local stress-energy and baryon conservation so finite-volume updates handle shocks with the correct weak solution. Geometry supplies nonconservative source terms, and every time step requires recovery of physical fluid variables from conserved ones. The method is credible only when smooth convergence, discontinuous Riemann problems, integrated balances, and source-term equilibria are all tested. This page develops the discretization; microphysical closure and astrophysical interpretation remain with Relativistic Matter and Compact Objects.

Helpful background. Formulations of General Relativity supplies the spacetime split, Ordinary Differential Equations supplies time integration, Relativistic Matter and Compact Objects supplies fluid closure, Computational Relativity supplies evidence standards, and Fluid and Kinetic Theory supplies balance laws and characteristic waves.

Valencia conserved variables

For Eulerian three-velocity $v^i$, rest density $\rho$, specific internal energy $\epsilon$, and pressure $p$, define

$$v_i=\gamma_{ij}v^j, \qquad W=\left(1-\gamma_{ij}v^iv^j\right)^{-1/2}, \qquad h=1+\epsilon+\frac{p}{\rho}$$

in units with $c=1$. A standard undensitized conserved state $\mathbf U=(D,S_j,\tau)$ is

$$D=\rho W, \qquad S_i=\rho hW^2v_i, \qquad \tau=\rho hW^2-p-D.$$

With metric determinant $\gamma$, the balance law has the schematic form

$$\partial_t(\sqrt\gamma\,\mathbf U) +\partial_i(\sqrt\gamma\,\mathbf F^i)= \sqrt\gamma\,\mathbf S.$$

Writing the transport velocity as $v_{\rm tr}^i=\alpha v^i-\beta^i$, the flux already multiplied by the appropriate lapse factors is

$$\mathbf F^i=( Dv_{\rm tr}^i,\; S_jv_{\rm tr}^i+\alpha p\,\delta^i{}_j,\; \tau v_{\rm tr}^i+\alpha p\,v^i ).$$

The lapse $\alpha$, shift $\beta^i$, spatial metric, and extrinsic curvature enter these fluxes and the geometric sources. Record whether densitized or undensitized variables are stored and how $c$ and $G$ are set. Banyuls et al. 1997, pp. 223–225, Eqs. (5)–(25) derives the formulation and characteristic structure.

This numerical regime assumes an inviscid perfect fluid with a causal EOS and a valid continuum closure. Magnetic stresses, radiation exchange, physical viscosity, and kinetic scales are outside this page’s system; using the update where a fluid closure breaks down does not become valid merely because the conservative residual is small.

Why the conserved density is not the rest density

Write the fluid velocity as $u^\mu=W(n^\mu+v^\mu)$ with $n_\mu v^\mu=0$. Baryon conservation is $\nabla_\mu(\rho u^\mu)=0$, or $\partial_\mu(\sqrt{-g}\rho u^\mu)=0$. Because $\sqrt{-g}=\alpha\sqrt\gamma$, $u^t=W/\alpha$, and $u^i=W(v^i-\beta^i/\alpha)$, it becomes

$$\partial_t(\sqrt\gamma D) +\partial_i[\sqrt\gamma D(\alpha v^i-\beta^i)]=0.$$

Thus the cell integral of $\sqrt\gamma D$, not of $\rho$ alone, counts rest mass on a spatial slice. The Lorentz factor accounts for motion relative to the Eulerian observer, and $\sqrt\gamma$ converts coordinate volume to proper volume. This derivation also checks that the baryon source is zero even when the metric evolves.

In Minkowski spacetime, take an ideal gas with $\Gamma=5/3$, $\rho=1$, $p=0.2$, and $v^x=0.6$. Then $\epsilon=0.3$, $h=1.5$, $W=1.25$, and

$$(D,S_x,\tau)=(1.25,1.40625,0.89375), \qquad (F_D^x,F_{S_x}^x,F_\tau^x)=(0.75,1.04375,0.65625).$$

These numbers follow directly from the definitions and make a compact forward-map test. All densities are in a common arbitrary density unit; $\epsilon$, $h$, $W$, and velocities are dimensionless. In the rest limit, the same formulas give $D=\rho$, $S_i=0$, and $\tau=\rho\epsilon$, which checks the rest-mass subtraction.

Reconstruct, solve, and update

A high-resolution shock-capturing method reconstructs primitive or characteristic states to cell faces, solves an approximate Riemann problem, and updates cell averages. The reconstruction, limiter, Riemann solver, wave-speed estimate, time integrator, CFL factor, and source quadrature form one algorithm. Characteristic reconstruction can reduce oscillations but requires a correctly normalized eigenbasis; componentwise reconstruction is simpler but can mix waves.

For smooth solutions, high-order reconstruction and time stepping should reach their designed order. At shocks, monotonicity requirements reduce local order; global $L_1$ convergence is often first order. Positivity or atmosphere fixes are additional nonconservative operations and must be counted in mass, momentum, and energy budgets rather than hidden.

Geometric source terms require a well-balanced perspective. A stationary TOV star can drift even when fluxes are accurate because discrete pressure gradients and gravity fail to cancel. Track both terms and refine them together.

Benchmark ladder

Use this one-dimensional special-relativistic shock tube as a fully specified first record. On $x\in[-0.5,0.5]$ in Minkowski spacetime, set a $\Gamma=5/3$ ideal-gas EOS and the left and right primitive states

$$(\rho,p,v^x)_L=(1,1,0), \qquad (\rho,p,v^x)_R=(0.125,0.1,0),$$

with the discontinuity at $x=0$, transverse velocity zero, and zero-gradient outflow boundaries. Use fifth-order WENO reconstruction, an HLL flux, SSPRK3, CFL $0.4$, and $N=200,400,800,1600$ uniform cells to $t=0.4$. At every $N$, repeat primitive-recovery tolerances $10^{-8},10^{-10},10^{-12}$. The reference is an independently evaluated exact self-similar special-relativistic Riemann solution, integrated to exact cell averages.

Compare density using

$$E_1(\rho)= \frac{\sum_j\Delta x\,|\rho_j-\bar\rho_{{\rm ex},j}|} {\sum_j\Delta x\,|\bar\rho_{{\rm ex},j}|},$$

and locate the rarefaction edges, contact, and shock from the exact solution rather than from visual alignment. Accept the $N=1600$, tolerance-$10^{-12}$ run only if $E_1(\rho)\leq1.5\times10^{-2}$, every wave location is within two finest cells of the reference, the boundary-flux-corrected baryon-budget residual is at most $10^{-10}$, and no density, pressure, subluminality, recovery, or floor failure occurs. Require the discontinuity-dominated $E_1$ sequence to approach first order and the result to change by at most $10^{-4}$ when tightening the recovery tolerance from $10^{-10}$ to $10^{-12}$.

Observed discrepancy: not available—specification only.

Complement the shock record with smooth linear sound waves in special relativity. At three or more resolutions, compare phase and amplitude with the analytic solution and demonstrate the designed order before boundaries matter. Additional relativistic shock tubes should span stronger shocks and rarefactions; report $L_1$ errors, overshoots, and physical-state failures.

Font 2008, §2.1.3 and §4 reviews conservative GR hydrodynamics and standard numerical tests. Add a stationary TOV star and spherical Bondi accretion on a fixed black-hole background to exercise geometric sources. Monitor central density, mass flux, baryon mass, and constraint-independent metric inputs.

Vary reconstruction, Riemann solver, and CFL factor one at a time. The exact shock tube should agree across robust solver choices after refinement, whereas excessive diffusion changes contact width in a predictable way. For Bondi flow, move inner and outer boundaries and compare the conserved accretion rate on multiple spheres. For the star, measure pulsation frequencies as well as secular drift; strong dissipation can suppress drift by damping the physical response and is not automatically an improvement.

For a fixed coordinate control volume $\Omega$, integrate the update using the outward physical unit normal $s_i$ and physical boundary area $dA$:

$$\Delta\int_\Omega \sqrt\gamma\,\mathbf U\,d^3x =-\int_{t_0}^{t_1}\oint_{\partial\Omega}\mathbf F^i s_i\,dA\,dt +\int_{t_0}^{t_1}\int_\Omega\sqrt\gamma\,\mathbf S\,d^3x\,dt.$$

Evaluate this budget with stored boundary fluxes and sources. Baryon number has no geometric creation term; energy and momentum in a general dynamical spacetime require careful interpretation, but the discrete balance law can still be checked.

Environment and failure criteria

Pin equations and variable definitions, EOS, reconstruction and limiter, Riemann solver, wave speeds, primitive recovery, floors, source discretization, metric data, grid, boundaries, integrator, CFL factor, precision, and code revision. A clean interface could be hydro-benchmark --config shock-tube-1.toml --output run/, emitting cell averages, primitives, face fluxes, sources, budgets, failure counters, refinement sweeps, and environment metadata. This specifies a rerun record rather than claiming an executable is present.

Reject a result with negative density or pressure outside a documented fallback, wrong shock speed, nonconvergent smooth-wave phase, unexplained baryon creation, or a stationary solution whose drift lacks a resolution trend. Failure modes include inconsistent densitization, using superluminal wave estimates incorrectly, limiter activation in smooth extrema, atmosphere domination, and geometric sources evaluated at a lower effective order.

Diagnostic

A shock tube matches the exact wave positions, but total baryon mass grows at each step even though boundary flux is zero. The growth decreases only when the atmosphere floor is lowered, not when resolution increases. Is the conservative scheme verified?

Solution

No. The finite-volume flux update may be conservative, but the atmosphere reset is injecting baryons and dominates the budget. Record reset mass explicitly, sweep the floor and threshold independently, and demonstrate a limit in which the physical solution and integrated balance are stable. Conservation claims must include every nonconservative repair.

References

• Banyuls, F., Font, J. A., Ibáñez, J. M., Martí, J. M., and Miralles, J. A. (1997). “Numerical 3+1 general relativistic hydrodynamics: A local characteristic approach.” Astrophysical Journal 476, 221–231. DOI.
• Font, J. A. (2008). “Numerical hydrodynamics and magnetohydrodynamics in general relativity.” Living Reviews in Relativity 11, 7. DOI.