Shocks, Weak Solutions, and Entropy Conditions

Compressive waves in relativistic Euler flow can steepen until derivatives diverge even though density, pressure, and velocity remain finite. Beyond that time the equations must be interpreted weakly: integral conservation selects the allowed jump relations, and an entropy condition excludes unphysical branches. This page derives those conditions and reports an executed, reproducible special-relativistic shock-tube verification against an independently implemented exact Riemann solution.

Required background. Hilbert Stress–Energy from Metric Variation supplies stress–energy conservation; Events, Intervals, and Minkowski Geometry supplies relativistic velocity; and Differential Equations supplies weak solutions and characteristics.

Helpful background. The Einstein Equation gives gravitational coupling; Relativistic Matter and Compact Objects gives the matter context; and Computational Relativity provides numerical evolution methods.

Conservation across a discontinuity

Assume that the shock carries no separate surface stress or particle reservoir. The continuum inputs are baryon conservation, stress–energy conservation, and an equation of state:

$$\nabla_a(n u^a)=0, \qquad \nabla_aT^{ab}=0, \qquad T^{ab}=(\varepsilon+p)u^a u^b+p g^{ab}.$$

Let a smooth hypersurface $\Sigma$ separate left and right states and let $s_a$ be any nonzero normal covector. Integrating over a thin pillbox and taking its thickness to zero gives the Rankine–Hugoniot conditions

$$\bigl[n u^a\bigr]s_a=0, \qquad \bigl[T^{ab}\bigr]s_a=0,$$

where $[Q]=Q_R-Q_L$. Scaling $s_a$ changes neither condition. In one-dimensional Minkowski coordinates, a balance law

$$\partial_tU+\partial_xF(U)=0$$

with a jump moving at speed $V_s$ obeys

$$V_s[U]=[F].$$

For one-dimensional special-relativistic hydrodynamics, let $\rho$ be rest-mass density, $e$ specific internal energy, $h=1+e+p/\rho$ specific enthalpy, and $W=(1-v^2)^{-1/2}$. A rest-mass-subtracted conservative state and its flux are

$$U= \begin{pmatrix} D\\ S\\ \tau \end{pmatrix} = \begin{pmatrix} \rho W\\ \rho hW^2v\\ \rho hW^2-p-D \end{pmatrix}, \qquad F= \begin{pmatrix} Dv\\ Sv+p\\ S-Dv \end{pmatrix}.$$

The total energy is $E=\tau+D$ and has flux $S$; mixing the $E$ and $\tau$ conventions would corrupt the energy jump. Andersson and Comer develop the related Eulerian variables and conservative-to-primitive map in Andersson and Comer 2021, §§ 11.2–11.3, Eqs. (11.21)–(11.50). Conservation alone admits compressive shocks, unphysical expansion shocks, and contact discontinuities unless the characteristic family and entropy behavior are also examined. Martí and Müller describe the hyperbolic wave structure in Martí and Müller 2015, § 3.2.1 and give the exact relativistic Riemann construction in § 8.5, pp. 128–133, Eqs. (179)–(200).

Entropy admissibility and wave families

For a smooth adiabatic perfect fluid, entropy per baryon satisfies $u^a\nabla_as=0$. Across a physical shock, irreversible microscopic processes occur in a layer omitted by the ideal equations, and downstream entropy must exceed upstream entropy. Equivalently, the entropy current satisfies a distributional inequality

$$\nabla_a(n s u^a)\geq0.$$

In one dimension choose the normal covector $s_a=(-V_s,1)$. Baryon conservation makes $j=n u^as_a$ the same on both sides. The entropy inequality at the discontinuity is then $j[s]\geq0$. For $j\gt0$, matter flows from left to right and requires $s_R\geq s_L$; if the flow reverses, the upstream and downstream labels reverse too. A contact has $j=0$ and permits an entropy jump without entropy production.

For a genuinely nonlinear characteristic family, the Lax compressivity condition states that its characteristics enter a shock from both sides: $\lambda_k(U_R)\lt V_s\lt\lambda_k(U_L)$. A rarefaction has the opposite structure: it is a continuous self-similar fan in which the relevant characteristic speed spreads monotonically. A contact discontinuity moves with the fluid; pressure and normal velocity remain continuous, while density and entropy may jump.

The entropy condition is not an extra conserved equation. It is an inequality inherited from the dissipative regularization as viscosity and heat conduction tend to zero. Different regularizations can matter for non-convex equations of state, phase boundaries, or undercompressive shocks. For a convex gamma-law gas, however, the standard relativistic Riemann solution consists of two nonlinear waves separated by a contact. The distinction between convex and non-convex wave patterns is reviewed in Martí and Müller 2015, § 3.2.3, and the rarefaction, shock, and matching branches are explicit in §§ 8.5.1–8.5.3, Eqs. (186)–(200).

Exact shock-tube benchmark

A demanding reference problem uses flat spacetime, domain $0\leq x\leq1$, discontinuity at $x=0.5$, outflow boundaries, and gamma-law equation of state

$$p=(\Gamma-1)\rho e, \qquad \Gamma=\frac53.$$

Take primitive initial data, rationalizing the source’s printed $p_L=13.33333$ as $40/3$,

$$(\rho,p,v)_L=\left(10,\frac{40}{3},0\right), \qquad (\rho,p,v)_R=(1,10^{-6},0).$$

The exact self-similar solution has a left-family rarefaction, a contact, and a right-going shock. The rarefaction’s head travels left, but its tail moves right for these data; the family label does not mean that every ray has negative coordinate speed. An exact relativistic Riemann solver determines the common star-region pressure by matching the velocity behind the left and right waves, then reconstructs density and wave speeds. It is the benchmark; a high-resolution numerical solution should not be compared only with another discretization.

The exact calculation gives

$$p_*=1.44794515599414, \qquad v_*=0.714020700932636,$$ $$\rho_{L*}=2.63929554954561, \qquad \rho_{R*}=5.07077596424752, \qquad V_s=0.828398034190528.$$

The normalized residuals of $V_s[U]-[F]$ are below $2\times10^{-12}$ componentwise. At $t=0.4$, the rarefaction head, rarefaction tail, contact, and shock are at $x=0.21355405$, $0.56689452$, $0.78560828$, and $0.83135921$. The shell speed, shell width divided by time, shock speed, and compression agree with the rounded published anchors $0.714020$, $0.114378$, $0.828398$, and $5.070776$. The source and benchmark are Martí and Müller 2015, §§ 6.3–6.3.1, pp. 76–81, Table 3.

The spacetime diagram turns these numbers into a wave pattern. Read horizontally at $t=0.4$: the expanding fan comes first, then the left star state, the contact, the compressed right star state, and the shock. The two star states share $p_*$ and $v_*$ but have different densities.

Quantitative spacetime wave fan for the stated gamma-law shock tube. Each boundary follows $x=0.5+\lambda t$; the dotted horizontal line samples the solution at $t=0.4$. Positions and slopes come from the retained exact Riemann solution. Rays inside the rarefaction indicate the continuous range of similarity speeds, rather than additional discontinuities.

The numerical contact and shock locations reported below are explicitly target-assisted diagnostics, not blind feature detections. Each is the unique midpoint-level crossing inside a window centered on the corresponding exact position, with a quoted uncertainty of $\pm\Delta x$. Requiring exactly one crossing makes an ambiguous feature inside that window fail rather than silently selecting the nearest candidate; oscillations outside the window or below the chosen midpoint level require separate diagnostics.

Reproducible numerical verification

The retained verification source is dependency-free JavaScript executed with Node.js 24.18.0 using IEEE binary64 arithmetic. It evolves finite-volume cell averages with primitive-variable piecewise-linear reconstruction, a monotonized-central limiter, an HLLE flux, and second-order strong-stability-preserving Runge–Kutta time stepping. Primitive recovery uses a safeguarded Newton–bisection pressure solve and requires both a normalized pressure residual and a conserved-state round-trip postcondition. No atmosphere or primitive-variable floor is applied: an inadmissible reconstruction falls back locally and is counted, while an inadmissible cell, exhausted solve, or failed postcondition terminates the run.

The shock convergence calculation uses $N=200,400,800$ uniform cells, constant outflow boundaries, and Courant factors $0.4$ and $0.2$. It compares numerical cell averages with exact cell averages, splitting quadrature at every wave. The conserved-density error and observed order are

$$E_1(N)=\frac1N\sum_{j=1}^N|D_j-D_j^{\mathrm{exact}}|, \qquad q_N=\log_2\!\frac{E_1(N)}{E_1(2N)}.$$

Courant factor	$N$	$E_1(D)$	$q_N$
0.4	200	0.0799949	0.940
0.4	400	0.0417104	0.779
0.4	800	0.0243073	—
0.2	200	0.0764177	0.939
0.2	400	0.0398695	0.784
0.2	800	0.0231612	—

The errors decrease under both refinements. Their pre-asymptotic orders lie below one because the diffused contact and shock dominate the global norm. A separate periodic entropy wave,

$$\rho(x,t)=1+0.2\sin[2\pi(x-0.3t)], \qquad p=1, \qquad v=0.3,$$

tests the smooth regime. For $N=50,100,200,400$, its $E_1(D)$ values are $7.47193\times10^{-4}$, $1.90671\times10^{-4}$, $4.82236\times10^{-5}$, and $1.19141\times10^{-5}$. The conserved fields $D,S,\tau$ and recovered density all give orders $1.970$, $1.983$, and $2.017$ to the shown precision. Pressure and velocity are exact invariants of this test; their largest deviations are respectively below $2.8\times10^{-15}$ and $7.3\times10^{-16}$. This full-state check separates the second-order smooth discretization from the lower global order of a discontinuous solution.

Halving the Courant factor is also compared directly at fixed spatial resolution. For each $Q\in\{D,S,\tau,\rho\}$, define

$$R_Q(N)=\frac{\lVert Q_{0.4}-Q_{0.2}\rVert_1}{E_{1,0.4}(Q)}.$$

The retained same-cell differences and the largest ratio over the four fields are

$N$	$\lVert\Delta D\rVert_1$	$\lVert\Delta S\rVert_1$	$\lVert\Delta\tau\rVert_1$	$\lVert\Delta\rho\rVert_1$	$\max_Q R_Q$
200	0.00449562	0.00625759	0.00487721	0.00368755	0.0694
400	0.00219680	0.00351726	0.00284457	0.00180231	0.0830
800	0.00131479	0.00179267	0.00134369	0.00105808	0.0646

All ratios remain below the declared $0.10$ gate. This is a timestep-sensitivity test, not a measurement of temporal order for a discontinuous solution.

Across every shock run, density and pressure remain positive, $|v|<0.721$, no reconstruction fallback occurs, and boundary-flux-corrected conservation residuals remain below $7\times10^{-14}$. The largest normalized primitive-recovery pressure residual is $1.84\times10^{-13}$ and the largest conserved-state round-trip residual is $2.85\times10^{-13}$, within gates of $2\times10^{-13}$ and $10^{-12}$. The exact star-pressure solve converges with a pre-alignment velocity mismatch below $2.3\times10^{-16}$. The exact rarefaction preserves $\ln(p/\rho^\Gamma)$ to $4.5\times10^{-16}$, while it increases by $11.48$ across the right shock. The bundle deliberately rejects negative pressure, superluminal primitive velocity, a non-dominant conserved state, and an expansion state passed to the shock branch. All 143 declared checks pass in the retained machine-readable results.

The public rerun is independent of the site’s private source repository. Start in an empty directory and download the complete single-file program together with the retained reference result:

mkdir srhd-shock-tube-rerun
cd srhd-shock-tube-rerun
curl --fail --location --output verify-srhd-shock-tube.mjs \
  https://blackhole.org/computations/compact-objects/verify-srhd-shock-tube.mjs
curl --fail --location --output reference-srhd-shock-tube-results.json \
  https://blackhole.org/computations/compact-objects/srhd-shock-tube-results.json
shasum -a 256 -c <<'SHA256'
4c4c0e42201118db57eeff7ab384b8cedfe29fd5eca35582ee76c1af2fc210fc  verify-srhd-shock-tube.mjs
e046d224941c675e88e2268d4aace3ad7be3961c88babd1e7e2d6cffd7d43f82  reference-srhd-shock-tube-results.json
SHA256

Use Node.js 24.18.0 to match the retained environment. The program needs only the Node.js standard library. Its portable no-argument mode recomputes the exact solution, all numerical refinements, every failure injection, and all 143 gates, then writes a fresh report to standard output:

node --version
node verify-srhd-shock-tube.mjs > rerun-srhd-shock-tube-results.json
node -e 'const fs=require("node:fs");const r=JSON.parse(fs.readFileSync("rerun-srhd-shock-tube-results.json","utf8"));if(r.summary.status!=="pass"||r.summary.failed!==0||r.summary.passed!==143||r.summary.total!==143)process.exit(1);console.log(`PASS ${r.summary.passed}/${r.summary.total} scientific checks`);'

The fresh report records its actual runtime and download location, so it is not expected to be byte-identical to the retained reference JSON. Its source hash, 143 named checks, physical outputs, tolerances, convergence results, and failure diagnostics are the comparison objects. No other script, data file, package, or environment file is required.

Failure modes are diagnostic. A negative reconstructed pressure signals an inadmissible state or failed primitive recovery; excessive contact smearing points to flux or limiter diffusion; oscillations behind the shock indicate inadequate limiting; a wrong wave order usually indicates a sign or relativistic-velocity error; apparent superconvergence often comes from comparing point values or shifting the exact solution to match the computed shock. Near vacuum, conservative-to-primitive inversion becomes ill-conditioned and requires an explicit admissibility or vacuum treatment. If floors are used, their mass and energy contributions must be measured; floors are not a mathematical requirement of every method.

Exercises

For a scalar conservation law with left and right constants, derive $V_s[U]=[F]$ by integrating over a thin spacetime control volume crossed by the discontinuity.

Solution

Integrate $\partial_tU+\partial_xF=0$ over a rectangle whose moving internal boundary is $x=V_st$. Flux through the two sides of that boundary must cancel. Per unit time, the boundary sweeps volume $V_s$, contributing $V_s(U_R-U_L)$, while the net physical flux is $F_R-F_L$. Equality gives

$$V_s(U_R-U_L)=F_R-F_L.$$

Applied componentwise to baryon number and stress–energy, this is the coordinate form of the covariant pillbox conditions.

References

• Andersson, Nils, and Gregory L. Comer. “Relativistic Fluid Dynamics: Physics for Many Different Scales.” Living Reviews in Relativity 24, 3 (2021). DOI.
• Martí, José María, and Ewald Müller. “Grid-Based Methods in Relativistic Hydrodynamics and Magnetohydrodynamics.” Living Reviews in Computational Astrophysics 1, 3 (2015). DOI.