General-Relativistic Hydrodynamical Simulations

• General-relativistic hydrodynamical simulations are computational techniques that model relativistic fluids in curved spacetime by solving conservation laws coupled with Einstein's field equations.
• Advanced numerical methods, including high-resolution shock-capturing schemes and flux-conservative formulations, accurately resolve shocks, turbulence, and strong spacetime curvature effects near compact objects.
• These simulations provide insights into accretion dynamics, torus formation, and turbulence, aiding predictions of electromagnetic and gravitational wave signals in extreme astrophysical environments.

General-relativistic hydrodynamical simulations are computational studies of self-gravitating, relativistic fluids—typically inviscid or with prescribed microphysical transport—on backgrounds with strong spacetime curvature. These simulations solve the covariant conservation laws of mass and energy-momentum for a perfect (or sometimes magnetized) fluid, coupled where necessary to the evolution of the spacetime metric via the Einstein field equations. Applications include accretion flows onto compact objects, stellar collapse, neutron star post-merger dynamics, and turbulence in strongly curved spacetime regions. Modern codes employ high-resolution shock-capturing schemes, flux-conservative formulations, and specialized coordinate systems and boundary conditions to handle horizons and singular geometries.

1. Covariant Framework and Equations of Motion

The foundation of general-relativistic hydrodynamics (GRHD) is the set of conservation equations:

• Mass conservation: $\nabla_\mu(\rho\,u^\mu) = 0$,
• Energy–momentum conservation: $\nabla_\mu T^{\mu\nu} = 0$,

where $\rho$ is the rest-mass density, $u^\mu$ the fluid four-velocity, and $T^{\mu\nu} = \rho h\,u^\mu u^\nu + p\,g^{\mu\nu}$ with specific enthalpy $h = 1 + \epsilon + p/\rho$, specific internal energy $\epsilon$, and fluid pressure $p$.

An equation of state (EOS), typically ideal $\Gamma$-law: $p = (\Gamma-1)\rho\,\epsilon$, specifies thermodynamics. Relativistic sound speed is $c_s^2 \equiv \Gamma p / (\rho h)$.

The flux-conservative formulation is most widely adopted (Valencia formalism). In a $3+1$ split:

• Conserved variables: $D = \rho W$, $S_j = \rho h W^2 v_j$, $\tau = \rho h W^2 - p - D$, where $W = \alpha u^0 = (1-v^i v_i)^{-1/2}$ is the Lorentz factor and $v^i = (u^i + \beta^i u^0)/(\alpha u^0)$ is the 3-velocity (see (Mach et al., 2018, Montero et al., 2013)).

2. Coordinate Systems and Metric Decomposition

To resolve the physics near horizons and avoid coordinate singularities, simulations employ specialized coordinate choices:

• Horizon-penetrating Eddington-Finkelstein coordinates:

$$ds^2 = -(1-2M/r)\,dt^2 + (4M/r)\,dt\,dr + (1+2M/r)\,dr^2 + r^2\,(d\theta^2 + \sin^2\theta\,d\phi^2)$$

allow the computational grid's inner boundary $r_{\rm in}$ to be set inside the event horizon ($r_{\rm in} < 2M$), avoiding artificial reflections or nonphysical constraints at the horizon (Mach et al., 2018).

• Spherical polar, curvilinear, or reference-metric approaches treat singular $1/r$ and $\cot\theta$ structures analytically in metric Christoffel symbols, eliminating the need for special regularization (Montero et al., 2013).

Typical $3+1$ split includes lapse $\alpha$, shift $\beta^i$, spatial 3-metric $\gamma_{ij}$, enabling flux/balance-law formulation.

3. Numerical Methods: Shock Capturing and Stability

GRHD codes employ high-resolution shock-capturing (HRSC) schemes:

• Riemann solvers: HLLE (robust, diffusive; used in strong gravity), HLLC (resolves contacts; preferred where sharp features are present).
• Reconstruction: Second-order (minmod, MC), third-order PPM, or WENO for primitive variables and velocities, often reconstructing on $$[\rho, \epsilon, \sqrt{\gamma^{ii} v_i v^i}\,{\rm sgn}(v_i)]$$ rather than directly on $v_i$ (Mach et al., 2018).
• Time integration: Method-of-lines (MoL), typically with TVD Runge-Kutta (second or third order), ensures stability across shocks (Mach et al., 2018).
• CFL criterion: $\Delta t$ chosen so that $\max(|\lambda_{\max}|)\,\Delta t/\Delta x \leq 0.5$–$0.8$, where $\lambda_{\max}$ are the local characteristic speeds.

Boundary conditions are crucial. Horizon boundary uses outflow/extrapolation (minimal influence on external flow). External boundary often imposes analytic inflow/outflow conditions matching Michel or Fishbone-Moncrief torus solutions, with latitude-dependent angular momentum enforced via $v_\phi = \ell(\alpha - \beta^r v^r)$ and $v_\theta$ either fixed or free (Mach et al., 2018).

4. Initial Conditions and Model Parameters

Initial data typically consist of:

• Exact transonic Michel or Fishbone-Moncrief background flows for radial inflow.
• Angular momentum profiles added via $l(r,\theta) = l_0 [1 - |\cos\theta|]$ outside transition radius $r_*$. $l_0$ is chosen for a desired circularization radius $r_C$ using $l_0^2 = M r_C / (1-2M/r_C)^2$ (see Table 1 in (Mach et al., 2018)).
• Asymptotic sound speed $c_{s,\infty}$ sets Bondi radius $r_B = M / c_{s,\infty}^2$.
• Models are run with varying $r_S'/r_B$ and $c_{s,\infty}$ to explore the transition from Bondi-like inflow to thick torus formation.

Perturbed models include Michel solution for $\rho,\epsilon,v_r$, $v_\theta=0$, $v_\phi$ set by angular momentum prescription.

5. Physical Results: Accretion Flows and Turbulence

Simulations reveal key physical effects of relativistic, low-angular-momentum accretion onto black holes:

• Even modest latitude-dependent angular momentum dramatically suppresses accretion rate to $\lesssim 30\%$ of the radial (Michel) value (Mach et al., 2018).
• For sufficiently low $c_{s,\infty}$ ($\lesssim 0.022$), an equatorial, geometrically thick torus forms outside the event horizon, stalling equatorial accretion; mass influx predominantly occurs along the poles via low-density funnels (time-averaged polar fraction >80–90% of accretion; see Figs. 20_dens–35_dens_free).
• The system evolves from axisymmetry to fully turbulent flow patterns, with angular momentum patches and accretion rate variability (Figs. 30_ang_free, 35_ang_free). Turbulence develops after $10^5$–$10^6\ M$ of evolution, depending on $c_{s,\infty}$.
• Shock-capturing codes resolve the attenuation/oscillation of mass and angular momentum fluxes ($\dot{m}(t), \dot{L}(t)$; see Figs. M20–M35, L20–L35). When $v_\theta$ is free at the outer boundary, accretion can transiently overshoot the Michel rate before decaying to a small fraction.
• Equatorial matter encounters a centrifugal barrier, preventing direct accretion; instead, buildup of high-entropy, turbulent torus supports non-radial flows.

6. Code Validation: Convergence and Benchmarks

Code accuracy is established by:

• Comparison with exact steady-state solutions: Michel radial accretion and stationary Fishbone-Moncrief tori.
• Second-order convergence demonstrated across the transonic point: rescaled pointwise errors in $\rho$ decrease as $N_r^{-2}$ with resolution ((Mach et al., 2018), Fig. 0a).
• $L^1$ errors in $\rho(t)$ measured at fixed intervals; global error norms scale as $N^{-2}$ (Fig. 0c).

These validation strategies confirm both the mathematical fidelity and the physical reliability of the simulation results.

7. Outlook: Extensions and Future Applications

Results from these GRHD simulations elucidate accretion morphologies distinct from Newtonian predictions, reinforce the suppression of mass accretion by angular momentum injection, and capture the turbulent dynamics of thick tori and polar inflow (Mach et al., 2018). The use of horizon-penetrating coordinates ensures that boundaries do not artificially affect the external flow.

To achieve realistic predictions for jet launching, angular-momentum transport, and energetic feedback, future extensions must incorporate:

• Viscosity, e.g., via GR Israel–Stewart formalism, to model dissipation and angular momentum exchange.
• Magnetic fields, requiring relativistic magnetohydrodynamics (GRMHD) to self-consistently include MRI-driven turbulence and jet formation.
• Microphysical processes and radiative cooling to match observed accretion flows.

These developments are essential for interpreting electromagnetic and gravitational wave signals from black hole environments and for bridging simulations to observable phenomena across the relativistic astrophysics domain.