General Relativity

Governing equations

Let g_μν be the components of a stationary metric with time coordinate x⁰, so ∂₀g_μν = 0. We can define the lapse α = (-g⁰⁰)^-1/2. Then the unit future-pointing timelike normal to surfaces of constant x⁰ has components n_μ = -αδ_0μ. Projection onto the spatial part of the tangent space can be accomplished with the components j_μν = g_μν + n_μn_ν.

When working with general relativity, the primitive variables in Athena++ are taken to be

• comoving rest mass density ρ,
• gas pressure p_gas,
• spatial components of the projected 4-velocity ũⁱ = jⁱ_μu^μ, and
• comoving passive scalar concentration χ (one for each species).

In MHD there will also be the magnetic field components Bⁱ = (★F)ⁱ⁰, where F is the electromagnetic field tensor. Note that these are not the projected variables ℬⁱ = -jⁱ_μn_ν(★F)^μν = αBⁱ.

The following quantities are useful to define:

• projected field components b^μ = u_ν(★F)^νμ, in particular
  • b⁰ = g_iμBⁱu^μ,
  • bⁱ = (Bⁱ + b⁰uⁱ)/u⁰;
• magnetic pressure p_mag = b_μb^μ/2;
• total pressure p_tot = p_gas + p_mag;
• adiabatic index Γ;
• gas enthalpy w_gas = ρ + Γ/(Γ−1) × p_gas;
• magnetic enthalpy w_mag = 2p_mag;
• total enthalpy w_tot = w_gas + w_mag;
• stress-energy components T^μν = w_totu^μu^ν + p_totg^μν − b^μb^ν;
• metric determinant g;
• connection coefficients Γ^σ_μν = 1/2 × g^σλ(∂_μg_νλ + ∂_νg_μλ − ∂_λg_μν).

The evolution equations are

• ∂₀ ((-g)^1/2ρu⁰) + ∂_j ((-g)^1/2ρu^j) = 0,
• ∂₀ ((-g)^1/2T⁰_μ) + ∂_j((-g)^1/2T^j_μ) = (-g)^1/2T^ν_σΓ^σ_μν,
• ∂₀ ((-g)^1/2(★F)ⁱ⁰) + ∂_j ((-g)^1/2(★F)^ij) = 0, and
• ∂₀ ((-g)^1/2χρu⁰) + ∂_j ((-g)^1/2χρu^j) = 0 (one equation for each species).

Athena++ considers the conserved variables to be

• conserved density D = ρu⁰,
• energy density E = T⁰₀,
• momentum density M_i = T⁰_i, and
• conserved species density S = χρu⁰ (one for each species),

as well as the magnetic fields Bⁱ. In particular, the factors of (-g)^1/2 are not included in the output data, nor should they be included when setting conserved variables in problem generators. Also, these are coordinate-frame quantities, not quantities in the normal frame often seen in 3+1 decompositions. This is the difference between components A^μ (coordinate frame) and -n_μA^μ and jⁱ_μA^μ (normal frame).

Considerations for GR in Athena++

To enable GR, add the -g flag when configuring. The -s flag should not be used, as this is reserved for SR without GR. There are several important choices one must make, detailed here.

Coordinates (--coord option)

Coordinate systems for GR are entirely distinct from those used in SR or Newtonian simulations. In particular, cartesian, cylindrical, and spherical_polar coordinates do not function in GR. Instead choices include:

• minkowski (t,x,y,z): ds² = -dt² + dx² + dy² + dz²;
• schwarzschild (t,r,θ,ϕ): ds² = -(1−2M/r)dt² + (1−2M/r)^-1dr² + r²dθ² + r²sin²(θ)dϕ²;
• kerr-schild (t,r,θ,ϕ): ds² = g_μνdx^μdx^ν, with:
  • Σ = r² + a²cos²(θ),
  • g₀₀ = -(1−2Mr/Σ),
  • g₀₁ = 2Mr/Σ,
  • g₀₃ = -(2Mar/Σ)sin²(θ),
  • g₁₁ = 1 + 2Mr/Σ,
  • g₁₃ = -(1 + 2Mr/Σ)a sin²(θ),
  • g₂₂ = Σ,
  • g₃₃ = (r² + a² + (2Ma²r/Σ)sin²(θ))sin²(θ),
  • g₀₂ = g₁₂ = g₂₃ = 0;
• gr_user: user-defined coordinates.

Kerr-Schild coordinates are a standard choice for black holes, since they are horizon penetrating.

Note that all primitive and conserved variables are in the coordinate basis. The natural basis vectors are not normalized, even in cases where it would be easy to do so such as Schwarzschild.

Equation of state (--eos option)

GR in Athena++ only works with an adiabatic (gamma-law) equation of state.

Riemann solver (--flux option)

In order of decreasing numerical diffusion, the Riemann solvers compatible with GR are llf (local Lax-Friedrichs), hlle, hllc, and hlld. HLLC only works for hydrodynamics, and HLLD only works for MHD. In addition, HLLC and HLLD can only work using frame transformations at interfaces to use SR-compatible algorithms (see the Athena++ GRMHD method paper), and thus they require the -t option. Frame transformations can be either on or off when using LLF or HLLE.

At the polar axis coordinate singularity in Schwarzschild and Kerr-Schild coordinates, frame transformations become singular. The Riemann solver will use the appropriate non-transforming method (LLF in that case, HLLE in all other cases) at those interfaces automatically.

Mesh refinement

Athena++'s mesh refinement, both static and adaptive, is designed to be fully compatible with GR. Note however that GR simulations often involve spherical coordinate systems, and polar boundaries require careful consideration.

Timestep size

The size of a timestep is set by the shortest light crossing time in the grid. It does not change even if the sound (hydrodynamics) or fast magnetosonic (MHD) speeds are considerably slower than unity (i.e. the speed of light). The x¹-width of cells is defined to be ∫(g₁₁)^1/2dx¹ along the line of constant x² and x³ from one face to the other passing through the cell center. The x²- and x³-widths are defined likewise. In coordinate systems where the integral has no simple antiderivative, lower bounds are used instead. Also, in complex coordinate systems the cell center may be defined as the arithmetic mean of the face coordinates rather than the halfway point in volume.

Floors and ceilings on variables

As with SR, the complexity of the equations of GR can prove troublesome for robustness. There are a number of ways in which a simulation can produce problematic values, often manifesting in recovering negative densities, negative pressures, or 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 = max(dfloor, rho_min × (x 1)rho_pow)
rho_min	dfloor
rho_pow	0
pfloor	≈10-35	pgas,floor = max(pfloor, pgas_min × (x 1)pgas_pow)
pgas_min	pfloor
pgas_pow	0
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

The power-law scalings of ρ_floor and p_gas,floor are meant to be used in spherical-like coordinates, where x¹ is the radial coordinate. These scaling components of the floors are ignored if the exponents are 0.

The Lorentz factor is in the normal frame: γ = -n_μu^μ = αu⁰ = (1 + g_ijũⁱũ^j)^1/2. If a value γ > γ_max is found, the velocity components ũⁱ are each multiplied by ((γ_max²−1)/(γ²−1))^1/2 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 ∂_i ((-g)^1/2Bⁱ) = 0.