The Piecewise Parabolic Method for Multidimensional Relativistic Fluid Dynamics

1. INTRODUCTION

• Relativistic astrophysical systems can involve flow velocities ≳0.99c, making special- or general-relativistic modeling necessary.Examples include jets, accretion flows, pulsar winds, gamma-ray bursts, and heavy-ion collisions.Source: §1 · p. 1
• Rapid spatial and temporal changes make relativistic problems difficult to study numerically and motivate more efficient algorithms.Source: §1 · p. 1
• High-resolution shock-capturing schemes use conservative fluid equations and Riemann problems to model steep gradients and discontinuities robustly.These methods provide at least second-order accuracy in smooth flow regions.Source: §1 · p. 2
• The paper extends the Piecewise Parabolic Method from one-dimensional relativistic flows to multidimensional relativistic hydrodynamics.The relativistic formulation couples normal and transverse velocities through the Lorentz factor and specific enthalpy.Source: §1 · p. 2
• The method combines third-order spatial interpolation, second-order characteristic time evolution, a two-shock iterative Riemann solver, and unsplit corner-transport upwinding.The formulation is presented for general orthogonal curvilinear coordinates.Source: §1 · p. 2

• The fluid equations enforce particle-number and energy-momentum conservation in divergence form in the laboratory frame.Source: §2 · p. 2
• The laboratory-frame variables are connected to rest-frame thermodynamic quantities through the Lorentz factor, including proper density ρ, specific enthalpy h, and pressure p.The specific enthalpy is related to the proper internal energy density e.Source: §2 · p. 3
• An equation of state closes the conservation system by expressing specific enthalpy as h = h(p, τ), with τ = 1/ρ, and defines the sound speed.The formulation assumes causal equations of state with cs < c, preserving hyperbolicity.Source: §2 · p. 3
• Relativistic states use conservative variables U = (D, m1, m2, m3, E) and primitive variables V = (ρ, v1, v2, v3, p).Momentum and velocity components are projections along the coordinate axes.Source: §2 · p. 3
• The inverse map V(U) is nonlinear because the Lorentz factor depends on pressure, so pressure is recovered from an implicit equation.The pressure dependence also enters τ and, generally, the specific enthalpy h(p, τ(p)).Source: §2 · pp. 3–4
• Newton-Raphson pressure recovery requires the derivative of the implicit equation, and the resulting U-to-V map is solved at every time step for the numerical method.The derivative is expressed using derivatives of the specific enthalpy.Source: §2 · p. 4

3. DESCRIPTION OF THE ALGORITHM

• The scheme integrates conservation laws over curvilinear control volumes, using generalized volume coordinates, face geometry, volumetric centroids, and geometrical source terms.The integral equations conserve mass, momentum, and energy; the conservation-law expressions are exact before numerical approximations are introduced.Source: §3 · pp. 5–6

3.1. Evolution in Multi-dimensions

• The multidimensional evolution uses Strang-split homogeneous and source operators, with a spatially unsplit corner-coupled update and second-order temporal accuracy.The homogeneous update proceeds through predictor states, transverse corrections, fully corner-coupled states, and Riemann-solver fluxes; geometrical sources use a second-order Runge-Kutta step and the timestep obeys a CFL condition.Source: §3.1 · pp. 7–9

3.2. The predictor step

• The predictor reconstructs fourth-order interface values from volume averages, defines piecewise parabolic cell profiles, and forms time-centered left and right states using characteristic information.Reconstruction is performed in the volume coordinate ξ; characteristic fields include two genuinely nonlinear and three linearly degenerate families, with states obtained by averaging profiles over characteristic domains of dependence.Source: §3.2 · pp. 9–11
• The predictor step constructs effective left and right states by retaining positive-speed characteristics for left states and negative-speed characteristics for right states.This selection is performed by expanding the matrix A(V) in eigenvalues and bi-orthonormal left and right eigenvectors.Source: §3.2 · pp. 11–12

3.3. Solution of the Riemann Problem

• The Riemann solver models a self-similar three-wave pattern and approximates rarefaction waves as shocks, while retaining contact discontinuity conditions and relativistic shock jumps.The contact discontinuity preserves normal velocity and pressure, whereas density and tangential velocities may jump.Source: §3.3 · p. 12
• Pressure matching across the contact discontinuity is solved with Newton-Raphson iteration, converging when relative pressure changes reach 10^-4-10^-6, typically within 6-7 iterations.The solution is parameterized by post-shock pressure, with the contact condition imposed through equality of left and right velocities.Source: §3.3 · p. 14

3.4. Equation of State

• The equation-of-state treatment contrasts the exact but computationally expensive relativistic perfect gas with analytic alternatives and addresses mass-flux singularities as pressure jumps vanish.For most presented equations of state, factoring [p] from [hτ] yields stable alternative expressions; the relativistic perfect-gas case requires a limiting value.Source: §3.4 · pp. 16–17
• The proposed TM equation of state satisfies Taub’s inequality at all temperatures, has correct limiting values, and is analytically solvable at nearly no computational cost.Its effective adiabatic index differs by less than 4% from the relativistic perfect-gas value.Source: §3.4 · pp. 16–17

4.1. Relativistic Shock Tubes

• The method reproduces relativistic shock-tube wave structures, with below-5% L−1 error in one case, but thin shells and tangential velocities remain challenging.At 400 zones, one test captures shocks and contacts over a few zones; another underestimates post-shock density by about 25%, improving to about 6% when resolution doubles.Source: §4.1 · pp. 18–19
• Large tangential velocity lowers post-shock pressure, normal velocity, shock strength, and mass flux, while right-state tangential velocity can produce denser, thicker shells.These effects arise from relativistic coupling between enthalpy and the Lorentz factor.Source: §4.1 · p. 19

4.2. Two-dimensional Riemann Problem

• The two-dimensional Riemann test preserves symmetry to about 10−9 in density and captures stationary contact discontinuities without typical diffusive smearing.The calculation develops colliding shocks and a double-shocked, drop-like region along the main diagonal.Source: §4.2 · pp. 19–20

4.3. Relativistic Spherical Shock Reflection

• The spherical shock reflection test reports relative errors of 10.7% in density, 0.7% in velocity, and 16.4% in pressure, with symmetries retained to 10−7.The solver nevertheless exhibits the wall-heating phenomenon typical of such schemes.Source: §4.3 · p. 20

4.4.1. Axisymmetric Jet

• The axisymmetric jet develops shocks, vortices, mixing, and a complex cocoon, while its propagation efficiency is approximately 0.92 and its head velocity broadly follows theory.The head decelerates at 40≲t≲60 before becoming comparable with the theoretical prediction again; the model favorably compares with Mart´ı et al. (1997).Source: §4.4.1 · p. 21, §4.4 · p. 21

4.4.2. Gamma Ray Burst

• The gamma-ray-burst jet remains narrowly collimated and reaches a similar final distance to Zhang et al. (2003), despite being about 50% faster before 0.7 s.The simulation develops a low-density polar funnel and thin cocoon, but numerical diffusion limits small-scale resolution, especially at large radii.Source: §4.4.2 · pp. 22–23

5. CONCLUSIONS

• The scheme evolves conserved quantities with a finite-volume Godunov formulation, solves Riemann problems iteratively using a two-shock approximation, and incorporates general equations of state.The two-shock approximation reduces computational cost and is reasonable for weak rarefaction waves when explicit time steps are sufficiently small.Source: §5 · p. 23
• The new relativistic perfect-gas equation of state matches the exact Synge (1957) solution to better than 4% while reducing Riemann-solver cost.It also has the correct low- and high-temperature asymptotic behavior and is consistent with relativistic kinetic theory.Source: §5 · p. 23
• The unsplit corner-transport integrator preserves flow symmetries without the special procedures often needed by directionally split advection.Second-order temporal accuracy uses characteristic projection operators, while multidimensional flux calculation requires 4 Riemann problems in two dimensions and 12 in three dimensions per cell per step.Source: §5 · p. 23
• The implementation is verified for Lorentz factors exceeding 2·10^3 and applied to multidimensional relativistic jets in cylindrical and spherical geometries.The applications include a pressure-matched light jet with Mach number ≈11 and Γ≈10, and a stellar-atmosphere jet with Γ>50.Source: §5 · p. 23

A. ORTHOGONAL CURVILINEAR COORDINATES

• Geometrical source terms and finite-volume geometric quantities are provided for arbitrary orthogonal coordinates, including Cartesian, cylindrical, and spherical systems.The appendix gives explicit scale factors, face areas, and zone volumes needed by the discretization.Source: A. ORTHOGONAL CURVILINEAR COORDINATES · p. 24
• Cartesian coordinates have unit scale factors and an identically zero geometric source term, with cell volume Δx_iΔy_jΔz_k.Cylindrical and spherical systems use coordinate-dependent scale factors and corresponding geometric measures.Source: A. ORTHOGONAL CURVILINEAR COORDINATES · pp. 24–25
• In spherical coordinates, the scale factors are h_r=1, h_θ=r, and h_φ=r sin θ.These factors determine the spherical geometric source and finite-volume measures.Source: A. ORTHOGONAL CURVILINEAR COORDINATES · p. 25

B. RECONSTRUCTION PROCEDURE

• The reconstruction follows CW84 using volume-coordinate interpolation, with quartic interpolation supplying single-valued zone-interface estimates.The procedure uses neighboring cell averages to construct the parabolic representation.Source: B. RECONSTRUCTION PROCEDURE · p. 25
• Primitive-variable reconstruction is preferred because it preserves pressure positivity and helps enforce relativistic kinematic and thermodynamic constraints.The paper notes that formal accuracy would instead favor conservative-variable interpolation.Source: B. RECONSTRUCTION PROCEDURE · p. 25
• Interface values are constrained by extrema over neighboring cells, involving 9 zones in two dimensions and 27 zones in three dimensions.This limiting step bounds reconstructed interface values before the parabolic profile is finalized.Source: B. RECONSTRUCTION PROCEDURE · p. 26
• Monotonicity violations trigger first-order interpolation or interface adjustment while preserving the correct zone average.The limiter handles local extrema and internal extrema of the parabolic profile.Source: B. RECONSTRUCTION PROCEDURE · p. 26
• Velocity interpolation reverts to first order when the reconstructed total velocity would exceed neighboring-cell bounds, preventing violation of v_1^2+v_2^2+v_3^2<1.This fallback enforces the relativistic velocity constraint during reconstruction.Source: B. RECONSTRUCTION PROCEDURE · p. 26

• Strong shocks trigger profile flattening to prevent unphysical oscillations in the reconstructed interface values.Source: §B.1 · p. 27
• The multidimensional flattening parameter χ measures flow convergence at cell corners and ranges from 0 to 1.Setting χ = 0 reduces the method to first-order accuracy.Source: §B.1 · pp. 27–28
• Flattening is restricted to compressive regions and controlled by shock width and strength, with βmin = 0.75 and βmax = 0.85.The adopted limits are Zmax = 0.75 and Zmin = 0.25.Source: §B.1 · p. 27
• An explicit diffusive flux is added to the conservative update to further reduce oscillations behind strong shocks.Source: §B.1 · p. 27