General Relativistic Magnetohydrodynamics

• General Relativistic Magnetohydrodynamics (GRMHD) is a framework that unifies magnetized plasma dynamics with general relativity to model high-energy astrophysical systems.
• It utilizes a 3+1 flux-conservative formulation and covariant methods to accurately simulate the complex interactions in curved spacetime.
• Advanced numerical schemes like constrained transport and robust primitive recovery ensure precise enforcement of divergence-free constraints and handle non-ideal effects.

General Relativistic Magnetohydrodynamic (GRMHD) Framework refers to the mathematical and numerical formulation of magnetized plasma dynamics in a curved, Lorentzian spacetime, typically as governed by the Einstein-Maxwell-MHD equations. The GRMHD model unifies the conservation laws for fluid matter and electromagnetic fields within the context of general relativity and has become the standard for describing high-energy astrophysical systems, such as accretion flows onto black holes, magnetized neutron stars, and relativistic jets.

1. Covariant Fundamentals and Symmetry Structure

The foundational GRMHD equations are formulated on a four-dimensional Lorentzian manifold $(M,g)$ of signature $(-,+,+,+)$. The essential conservation laws are:

• Baryon number conservation: $\nabla_{\mu}(n\,u^{\mu})=0$,
• Energy-momentum conservation: $\nabla_{\mu}T^{\mu\nu}=0$, where $T^{\mu\nu}$ is the total (fluid + electromagnetic) stress-energy tensor,
• Maxwell’s equations: $\nabla_{\nu} F^{*\mu\nu}=0$ in the ideal MHD (comoving electric field vanishes: $F_{\mu\nu}u^{\nu}=0$).

For stationary and axisymmetric spacetimes with Killing fields $\xi$ (timelike, $\mathcal{L}_{\xi}g=0$) and $\chi$ (spacelike, $\mathcal{L}_{\chi}g=0$), the metric can be described with adapted coordinates $(t, x^1, x^2, \varphi)$ and the electromagnetic field tensor $F$ decomposed uniquely in terms of scalar potentials $\Phi$ (electric), $\Psi$ (magnetic), and $I$ (toroidal current) as

$$F = d\Phi \wedge \xi + d\Psi \wedge \chi + \frac{I}{\sigma} \epsilon(\xi, \chi, \cdot, \cdot),$$

where $\sigma = V X + W^2$ with $V = -\xi\cdot\xi$, $W = \xi\cdot\chi$, $X = \chi\cdot\chi$, and $\epsilon$ is the volume form (Gourgoulhon et al., 2011).

2. Conservative 3+1 Formulation and Variable Sets

In the $3+1$ split, the spacetime metric is written as

$$ds^2 = -\alpha^2 dt^2 + \gamma_{ij} (dx^i + \beta^i dt)(dx^j + \beta^j dt),$$

where $\alpha$ is the lapse, $\beta^i$ is the shift, and $\gamma_{ij}$ is the spatial metric. The equations are cast into flux-conservative Valencia form: $$\partial_t (\sqrt{\gamma} U) + \partial_i (\sqrt{\gamma} F^i) = \sqrt{\gamma} S,$$ with conserved variables $U = [D, S_j, \tau, B^j]$, where

• $D = \sqrt{\gamma} \rho W$,
• $$S_j = \sqrt{\gamma}[(\rho h + b^2)W^2 v_j - \alpha b^0 b_j]$$,
• $$\tau = \sqrt{\gamma}[(\rho h + b^2)W^2 - (p + \frac{1}{2}b^2) - (\alpha b^0)^2] - D$$,
• $B^j$ the Eulerian magnetic field (Moesta et al., 2013, Etienne et al., 2015, Cipolletta et al., 2019, Kalinani et al., 2021, White et al., 2015).

Fluxes and sources depend nonlinearly on the primitive variables $(\rho, \epsilon, v^i, B^i)$; thus, robust recovery/inversion schemes are essential for stable evolution at high Lorentz factors and magnetization.

3. Geometrical Equilibrium Equations and Master Potentials

For stationary, axisymmetric configurations, all equilibrium equations can be consolidated into a master nonlinear second-order partial differential equation (the relativistic Soloviev transfield equation): $$A \Delta \Upsilon + \frac{n}{h} \gamma \nabla\left(\frac{h}{\gamma}\right) \cdot \nabla\Upsilon + (\gamma \gamma' - \mu_0^{-1} h^{-1} \beta^2)\, \nabla\Upsilon \cdot \nabla\Upsilon + \cdots\ d\Upsilon = 0,$$ where $\Upsilon$ is the "master potential", $A$ involves combinations of the spacetime metric scalars and prescribed functions $\alpha(\Upsilon)$, $\beta(\Upsilon)$, $\gamma(\Upsilon)$, representing the freedom in equilibrium "surface functions" (e.g., field rotation, current) (Gourgoulhon et al., 2011).

Special subcases of $\Upsilon$ yield various equilibrium constructions:

• The relativistic Grad-Shafranov equation for poloidal fields,
• Purely toroidal field equilibria,
• The relativistic Stokes equation for unmagnetized flows.

4. Numerical Schemes and Divergence-Free Enforcement

Modern GRMHD codes employ a combination of high-order finite-volume shock-capturing methods for robust evolution:

• Reconstruction: TVD, PPM, WENO, ENO, MP5, etc., for primitive state interpolation at cell faces.
• Riemann solvers: HLLE, HLLC/HLLD (with local Lorentz-frame transformation), Roe-type, for approximate solution of the local MHD wave structure.
• Time integration: Method-of-lines with strong-stability-preserving Runge-Kutta schemes.
• Constrained transport (CT): Staggered-mesh discrete curl update for $B^i$, ensuring divergence-free constraint $\partial_i(\sqrt{\gamma} B^i) = 0$ to machine precision (Moesta et al., 2013, White et al., 2015, Etienne et al., 2015, Prather et al., 2021, Cipolletta et al., 2019).
• Alternative divergence control: Hyperbolic divergence cleaning (Dedner) and 8-wave (Powell) schemes, including generalized versions suitable for mesh-free methods (Fedrigo et al., 18 Jun 2025, Moesta et al., 2013).

Development in mesh-free (MFV, MFM) discretizations extends GRMHD to particle-based codes, coupling conservative update equations to generalized Riemann solvers and divergence cleaning (Fedrigo et al., 18 Jun 2025).

5. Non-Ideal Effects and Generalized Ohm’s Law

The ideal GRMHD limit is valid when resistive, Hall, and electron-inertial terms remain negligible compared to the leading-order MHD terms. For two-fluid (non-ideal) generalizations, the generalized Ohm law reads

$$\frac{h J^{\mu}}{n^2 e} + \frac{2\Delta h}{n^2 e} U^{\mu} = F^{\mu\nu} U_{\nu} - \eta \big[ J^{\mu} - q_e(1+\Theta)U^{\mu} \big] + \frac{1}{4 n e} \nabla_{\nu}[(h_+ - h_-) (U^\mu J^\nu + J^\mu U^\nu)],$$

with $\Theta$ the energy-exchange term determined by closure (Koide, 2020, Gorard et al., 29 Oct 2025). At sufficiently large scales and small resistivities, the single-fluid ideal MHD equations are recovered, with electrons’ infinite mobility enforcing $F_{\mu\nu} u_e^\nu \rightarrow 0$ (Gorard et al., 29 Oct 2025).

For applications near black holes and in neutron star interiors, the plasma skin depth, reconnection layer thickness, and electron inertia set the limits for the validity of ideal GRMHD (Koide, 2020, Andersson et al., 2021). Global flows in, e.g., M87* accretion disks are classically in the ideal regime, though local reconnection may require explicit non-ideal treatment.

6. Code Architectures, Primitive Recovery, and Exascale Implementation

GRMHD simulations are enabled by large-scale, highly parallel code architectures:

• Structured-grid codes: e.g., GRHydro and IllinoisGRMHD (block-structured AMR, BSSN spacetime, vector-potential CT) (Moesta et al., 2013, Etienne et al., 2015, Cipolletta et al., 2019, Etienne et al., 2010).
• GPU and exascale codes: GRaM-X (GPU AMR with AMReX, hybrid MPI/OpenMP/CUDA, Z4c spacetime, WENO5) (Shankar et al., 2022); GR-Athena++ (oct-tree AMR, task DAG, Z4c) (Daszuta et al., 7 Jun 2024).
• Mesh-free: GIZMO GRMHD (kernel-weighted MFV/MFM, divergence cleaning) (Fedrigo et al., 18 Jun 2025).
• Discontinuous spectral: mapped Chebyshev-Fourier grids for exponential spatial convergence (Li, 14 Aug 2025).
• Radiation-MHD: Athena++ extensions with finite-solid-angle discretization, implicit matter–radiation coupling (White et al., 2023).

A major bottleneck is conservative-to-primitive variable inversion, especially in regions of high magnetization and Lorentz factor. Newer robust schemes (e.g., RePrimAnd) guarantee unique, physically admissible inversion, outperforming standard multidimensional Newton–Raphson methods in extreme regimes (Kalinani et al., 2021). Primitive inversion in multifluid models is purely hydrodynamic and decoupled from magnetization (Gorard et al., 29 Oct 2025).

7. Extensions and Hybrid Regimes

Recent works extend the framework:

• Hybrid GRMHD+GRFFE: Explicitly blend ideal GRMHD with force-free electrodynamics in highly magnetized ($\sigma \gg 1$) regions, eliminating the need for unphysical density floors, crucial for jet funnel modeling (Chael, 1 Apr 2024).
• Multifluid and Non-Ideal MHD: Generalized systems can treat electron inertia, Hall and battery effects, and support strictly hyperbolic well-posedness—even in regimes inaccessible to single-fluid GRMHD (e.g., $\sigma \gtrsim 10^8$, $W \gtrsim 10^4$) (Gorard et al., 29 Oct 2025, Koide, 2020, Andersson et al., 2021).
• Reference-metric and Non-Cartesian Coordinates: SphericalNR implements Valencia+CT GRMHD in a reference-metric BSSN/CCZ4 setting for full dynamical evolution in spherical coordinates, with orthonormal variable representations (Mewes et al., 2020).

A summary of implementation approaches and features:

Codebase / Framework	Features	Notes
GRHydro, IllinoisGRMHD	BSSN dynamical spacetime, CT, AMR	Open-source, widely validated (Moesta et al., 2013, Etienne et al., 2015)
Athena++, GR-Athena++	Advanced Riemann, AMR, CT	Frame-transforms for HLLC/HLLD, exascale scaling (White et al., 2015, Daszuta et al., 7 Jun 2024)
iharm3D	Arbitrary stationary spacetime	HARM algorithm, Flux-CT, soA layout (Prather et al., 2021)
GIZMO-GRMHD	Mesh-free MFV/MFM, divergence cleaning	GR extension of MFV/MFM methods (Fedrigo et al., 18 Jun 2025)
SphericalNR	Reference-metric, spherical coords	Orthonormal storage, AMR, vector potential (Mewes et al., 2020)
GRaM-X	GPU, AMReX, Z4c, WENO, tabulated EoS	GPU exascale AMR (Shankar et al., 2022)
RePrimAnd	Robust primitive recovery	Unique bracketing, error handling in Spritz (Kalinani et al., 2021)
Hybrid GRMHD+GRFFE	Smooth $\sigma$-driven switch	Jet funnel simulations, removes floor-induced artifacts (Chael, 1 Apr 2024)

• "Magnetohydrodynamics in stationary and axisymmetric spacetimes: a fully covariant approach" (Gourgoulhon et al., 2011)
• "The physics of non-ideal general relativistic magnetohydrodynamics" (Andersson et al., 2021)
• "Implementing a new recovery scheme for primitive variables in the general relativistic magnetohydrodynamic code Spritz" (Kalinani et al., 2021)
• "GRHydro: A new open source general-relativistic magnetohydrodynamics code for the Einstein Toolkit" (Moesta et al., 2013)
• "IllinoisGRMHD: An Open-Source, User-Friendly GRMHD Code for Dynamical Spacetimes" (Etienne et al., 2015)
• "An Extension of the Athena++ Code Framework for GRMHD Based on Advanced Riemann Solvers and Staggered-Mesh Constrained Transport" (White et al., 2015)
• "A general relativistic magnetohydrodynamics extension to mesh-less schemes in the code GIZMO" (Fedrigo et al., 18 Jun 2025)
• "Hybrid GRMHD and Force-Free Simulations of Black Hole Accretion" (Chael, 1 Apr 2024)
• "Beyond GRMHD: A Robust Numerical Scheme for Extended, Non-Ideal General Relativistic Multifluid Simulations" (Gorard et al., 29 Oct 2025)
• "Generalized general-relativistic magnetohydrodynamic equations for plasmas of active galactic nuclei in the era of the Event Horizon Telescope" (Koide, 2020)
• "SphericalNR: A new dynamical spacetime and general relativistic MHD evolution framework for the Einstein Toolkit" (Mewes et al., 2020)
• "GR-Athena++: magnetohydrodynamical evolution with dynamical space-time" (Daszuta et al., 7 Jun 2024)