Special Relativity

Governing equations

Athena++ supports special-relativistic fluid dynamics. In this mode the primitive variables are

• comoving rest mass density ρ,
• gas pressure p_gas,
• spatial components uⁱ of the 4-velocity, and
• passive scalar dimensionless concentration χ (one for each species).

In MHD we also have a magnetic field B with components Bⁱ.

The conserved variables are

• coordinate-frame density D = γρ,
• total energy density E = γ²w_tot − p_tot − b^tb^t,
• momentum density M with components Mⁱ = γ²w_totvⁱ − b^tbⁱ, and
• coordinate-frame species density S = χD.

Here we define the following useful quantities:

• Lorentz factor γ = (1 + u¹u¹ + u²u² + u³u³)^1/2,
• 3-velocity components vⁱ = uⁱ/γ,
• projected 4-field with components b^t = γ v⋅B and bⁱ = Bⁱ/γ + b^tvⁱ,
• magnetic pressure p_mag = η_μνb^μb^ν/2,
• total pressure p_tot = p_gas + p_mag,
• adiabatic index Γ,
• comoving gas enthalpy w_gas = ρ + Γ/(Γ−1) × p_gas,
• comoving magnetic enthalpy w_mag = 2p_mag,
• comoving total enthalpy w_tot = w_gas + w_mag.

The evolution equations are then

• ∂_tD + ∂_j (γρv^j) = 0,
• ∂_tE + ∂_jM^j = 0,
• ∂_tMⁱ + ∂_j (γ²w_totvⁱv^j + p_totδ^ij − bⁱb^j) = 0,
• ∂_tBⁱ + ∂_j (γv^jbⁱ − γvⁱb^j) = 0, and
• ∂_tS + ∂_j (χγρv^j) = 0.

The momentum equation might have source terms on the right-hand side in non-Cartesian coordinates.

General code considerations

To enable special relativity, add the -s flag when configuring. This option will automatically select an appropriate Riemann solver and equation of state. The Riemann solver can be local Lax-Friedrichs, HLLE, HLLC (hydrodynamics only), or HLLD (MHD only). Note the Roe solver is unavailable for relativity. The equation of state must be adiabatic (i.e. gamma-law).

SR is compatible with the same set of coordinates as Newtonian physics: Cartesian, cylindrical, and spherical-polar. Note the coordinates labeled "Minkowski" are for use with general relativity only. For those with GR experience, it is worth noting that cylindrical and spherical-polar coordinates are implemented in a locally orthonormal way for SR, rather than in the natural coordinate way. For example, a purely radial velocity field that is everywhere half the speed of light will have constant components (0.5, 0, 0) in spherical coordinates, rather than spatially varying components and understood basis vectors that also vary in magnitude spatially. This is reflected in the output data.

When calculating timesteps in SR, Athena++ assumes a signal speed of the speed of light (i.e. unity). Therefore a fixed grid will have fixed timesteps given by the CFL number multiplied by the smallest cell width.

SR is fully compatible with static and adaptive mesh refinement.

Floors and ceilings on variables

Caution is advised when simulating extreme conditions (very high Lorentz factors, high magnetizations, etc.). Due to the extra nonlinear coupling in SR compared to Newtonian physics, there are more ways in which a simulation can go wrong, often manifesting in recovering superluminal velocities in the conserved-to-primitive variable inversion.

In order to help robustness, various limits are available in the <hydro> section of the input file, described below:

Name	Default Value	Description
dfloor	≈10-35	floor on ρ
pfloor	≈10-35	floor on pgas
sfloor	≈10-35	floor on χ
gamma_max	1000	ceiling on γ
sigma_max	0	ceiling on 2pmag/ρ; only used if positive
beta_min	0	floor on pgas/pmag; only used if positive

If Lorentz factors above γ_max are found, the primitive velocities are all rescaled by the same factor in order to make γ = γ_max. For the limits on plasma σ and β, ρ and p_gas are adjusted and the magnetic field is left alone in order to preserve ∇⋅B = 0.