Global GRMHD Simulations

• Global GRMHD simulations are first-principles computational models that solve covariant relativistic MHD equations in curved spacetimes.
• They employ advanced numerical algorithms, adaptive grids, and high-order solvers to resolve turbulence, MRI, and dynamic jet-launching processes.
• Applications span black hole accretion, jet-disk coupling, supernovae, and neutron star systems, providing insights into high-energy astrophysical phenomena.

Global general relativistic magnetohydrodynamics (GRMHD) simulations provide a first-principles framework for predicting the plasma dynamics, energetics, and emergent phenomenology of strongly magnetized flows in astrophysical systems where both relativistic gravity and magnetic fields are essential. Utilizing discretized solutions of the covariant relativistic MHD equations in curved spacetimes, these global simulations resolve the nonlinear coupling between MHD turbulence, accretion, jet launching, radiation, and dissipative microphysics. They have established themselves as the foundational method for investigating black hole accretion disks, jet-disk coupling, neutron star accretion, core-collapse supernovae, binary black hole mergers in circumbinary disks, and collapsar engines.

1. Fundamental Equations and Physical Scope

Global GRMHD simulations solve the equations of ideal relativistic MHD on a curved spacetime background, typically either a fixed (Kerr or Schwarzschild) or dynamical metric in the BSSN formalism. The system consists of:

• Covariant mass conservation: $\nabla_\mu(\rho u^\mu) = 0$
• Conservation of energy-momentum: $\nabla_\mu T^{\mu\nu} = 0$
• Homogeneous Maxwell induction: $\nabla_\mu (^*F^{\mu\nu}) = 0$ where $$T^{\mu\nu} = (\rho + u_g + p_g + b^2)u^\mu u^\nu + (p_g + \tfrac12 b^2)g^{\mu\nu} - b^\mu b^\nu$$ for the total stress-energy tensor in ideal MHD, and $b^\mu = (^*F^{\mu\nu}u_\nu)$ encodes the comoving magnetic field (Santhiya et al., 1 Dec 2025, Ball et al., 2017, Issa et al., 3 Oct 2024, Shiokawa et al., 2011).

For applications involving radiation or microphysical processes, additional evolution equations for the radiation stress-energy tensor $R^{\mu}_{\ \nu}$ using M1 closure (Mishra et al., 2016, Sadowski et al., 2014), neutrino transport (Issa et al., 3 Oct 2024), or other source terms are included.

2. Numerical Algorithms and Grid Structures

Key features of global GRMHD codes include:

• Horizon-penetrating curvilinear coordinates (e.g., Kerr–Schild, modified Boyer–Lindquist) to regularize the spacetime near compact objects (Ball et al., 2017, Liska et al., 2018, Li, 14 Aug 2025).
• Adaptive mesh refinement or static multiblock grids for efficient resolution near disks, jets, or collapsing cores (Santhiya et al., 1 Dec 2025, Issa et al., 3 Oct 2024, Nathanail et al., 2020).
• High-order finite-volume or (more recently) discontinuous spectral solvers with approximate Riemann solvers (HLL, HLLE, LLF) for balancing accuracy and shock-capturing (Li, 14 Aug 2025, Das et al., 25 Nov 2024, Nathanail et al., 2020).
• Constrained transport schemes to maintain $\nabla\cdot\mathbf{B} = 0$ to machine precision (Li, 14 Aug 2025, Santhiya et al., 1 Dec 2025, Nathanail et al., 2020).
• Advanced inversion algorithms for converting between conservative and primitive variables, often via multidimensional Newton–Raphson solvers (Santhiya et al., 1 Dec 2025, Das et al., 25 Nov 2024).
• GPU-acceleration for massively parallel 3D simulations, enabling problem sizes up to $\sim 10^{9}$ cells (Liska et al., 2018, Shankar et al., 15 Apr 2025, Issa et al., 3 Oct 2024).

Boundary conditions are generally outflow in radial domains and periodic or reflecting in angular directions. Density and pressure floors are imposed in low-density regions to avoid numerical failures.

3. Initial and Boundary Conditions for Astrophysical Systems

Initial conditions are set by constructing equilibrium tori (Fishbone–Moncrief), translationally mapped stellar models (for supernovae/collapsars), or vertically-integrated disks with prescribed magnetic topology:

• SANE (Standard and Normal Evolution): multiple poloidal loops with alternating polarity, yielding low net flux through the horizon.
• MAD (Magnetically Arrested Disk): a single large-scale poloidal loop produces a dynamical magnetization at the event horizon, characterized by a horizon magnetic flux $\Phi_{\rm BH}/\sqrt{\dot M} \gtrsim 50$ (Ball et al., 2017, Liska et al., 2018, Santhiya et al., 1 Dec 2025).
• "Sub-SANE"/multi-loop setups: multiple small-scale, same-polarity loops designed to erase memory and allow MRI-driven dynamo emergence (Santhiya et al., 1 Dec 2025).
• Purely toroidal or mixed configurations to paper in situ dynamo processes (Liska et al., 2018).
• For neutron stars, the magnetosphere and light-cylinder structure are imposed via analytic dipole solutions, and the disk is embedded as an equilibrium torus (Das et al., 25 Nov 2024).

For full dynamical-spacetime simulations, initial metrics are prescribed as either binary configurations (e.g., puncture or conformal-thin-sandwich) or dynamically evolved via BSSN/Z4c while evolving MHD fields (Li, 14 Aug 2025, Shankar et al., 15 Apr 2025, Farris et al., 2012). Density and pressure floors, as well as resistive switches (for hybrid GRMHD–force-free), are applied as needed (Chael, 1 Apr 2024).

4. Diagnostics, Turbulence, and MRI Resolution

The effective resolution of global GRMHD simulations is assessed with:

• Magnetorotational instability (MRI) quality factors: $$Q_{\theta} \approx 2\pi / (\Omega \Delta x^{\theta}) |B^{\theta}| / \sqrt{\rho}$$, $Q_{\phi}$ analogously.
• Convergence in volume-averaged diagnostics: plasma $\beta$, magnetization $\sigma$, field correlation lengths, and synthetic spectra (Shiokawa et al., 2011).
• Outflow diagnostics: time- and azimuthally averaged accretion rates, magnetic flux, jet/outflow power (Blandford–Znajek efficiency $\eta_{\rm jet}$), mass ejection criteria (e.g., $-h u_t > 1$ for unbound material) (Santhiya et al., 1 Dec 2025, Issa et al., 3 Oct 2024, Sadowski et al., 2014).
• Turbulence analysis: decomposition of the induction equation into advection, compression, and turbulent dynamo terms to quantify the sustenance of large-scale field cycles in the nonlinear regime (Santhiya et al., 1 Dec 2025).
• For reconnection studies, sheet-finding and plasmoid-tracking algorithms are applied to isolate current sheets and measure statistical distributions of $\beta$, $\sigma$, and guide field magnitudes in SANE vs. MAD topologies (Ball et al., 2017, Nathanail et al., 2020).

Global convergence studies confirm robust thermal structure and spectra for $N \gtrsim 200^3$, but full convergence of magnetic correlation lengths requires further increases in linear resolution (Shiokawa et al., 2011).

5. Key Physical Insights: Jet Launching, Dynamo, and Instabilities

Global GRMHD simulations have established several physically robust results:

• MRI-driven turbulence mediates angular momentum transport, allows for the in situ generation of large-scale poloidal flux via mean-field $\alpha$–$\Omega$ dynamo cycles, even from a purely toroidal seed field (Liska et al., 2018, Santhiya et al., 1 Dec 2025).
• Rapid emergence of large-scale poloidal loops results in inward advection of flux (on the accretion timescale), leading to eventual pile-up at the black hole, triggering the MAD state ($\Phi_{\rm BH}/\sqrt{\dot M} \sim 50$). The accretion flow transitions from SANE to MAD, with discernible changes in the magnetization profile, jet efficiency, and field geometry (Ball et al., 2017, Liska et al., 2018, Chael, 1 Apr 2024).
• Jet launching via the Blandford–Znajek mechanism: measured jet efficiency $\eta_{\rm jet}$ can exceed unity (i.e., $P_{\rm jet} > \dot M c^2$) in extreme MAD states, with jet Lorentz factors $\gamma\sim10$, parabolic collimation, and stability against large-scale instabilities over $>10^3 r_g$ (Liska et al., 2018).
• Sustained jet power and longevity require not only a sufficient supply of large-scale poloidal flux at the horizon but, crucially, high magnetic-field coherence, quantified via the signed-to-unsigned flux ratio $\mathcal{C}_{\rm BH}$; jets shut off for $\mathcal{C}_{\rm BH} \lesssim 0.6$ (Santhiya et al., 1 Dec 2025).
• Dissipation in current sheets and reconnection layers is ubiquitous: SANE flows produce sheets with $0.1\lesssim\beta\lesssim10^3$, $10^{-4}\lesssim\sigma\lesssim1$ and weak guide fields, while MAD flows exhibit stronger guide fields and higher magnetizations in sheets ($10^{-3}\lesssim\sigma\lesssim10$) (Ball et al., 2017). Plasmoid formation in these sheets drives episodic flaring activity and is a natural source of rapid electromagnetic variability (Nathanail et al., 2020).

6. Microphysics, Radiation, and Multiphysics Extensions

Several extensions of global GRMHD are critical for modeling observed phenomena:

• Radiation transport is commonly included via M1 moment schemes for photons (Mishra et al., 2016, Sadowski et al., 2014), or with two-moment gray schemes for neutrinos (Issa et al., 3 Oct 2024), both crucial for modeling thin, thermally unstable disks and collapsar disks.
• Subgrid dynamo models emulate unresolved 3D turbulence in axisymmetric simulations, allowing for sustained MRI turbulence over many thermal times (Sadowski et al., 2014).
• Hybrid GRMHD–force-free methods are deployed to treat highly magnetized jet funnels ($\sigma \gtrsim 10^6$), avoiding floor artifacts and enabling accurate postprocessing of horizon-scale images (Chael, 1 Apr 2024).
• Applications to neutron stars demand high-resolution treatment of disk–magnetosphere coupling, spot formation, and radiative transfer for pulse-profile modeling (Das et al., 25 Nov 2024).

7. Applications and Astrophysical Implications

Global GRMHD simulations have been deployed across a wide range of scenarios:

• Sgr A* and M87* disk and jet structure, direct EHT-relevant modeling, and comparison to field strengths inferred from VLBI polarimetry (Sen et al., 4 Oct 2025).
• Jet formation from collapsar disks and neutron-rich outflows relevant for $r$-process nucleosynthesis, with full neutrino transport and MAD transitions (Issa et al., 3 Oct 2024).
• Core-collapse supernovae, with parameter-space surveys in the rotation–field plane, showing thresholds for jet-dominated, magnetorotational explosions (Shankar et al., 15 Apr 2025).
• Jet modulation and variability from magnetic reconnection, plasmoid and flux-rope ejection, explaining impulsive multiwavelength flaring in AGN (Nathanail et al., 2020, Ball et al., 2017).
• Accretion in tilted, precessing, and binary black hole systems, enabling predictions for QPOs and multi-messenger counterparts (Teixeira et al., 2014, Farris et al., 2012).
• High-order discontinuous spectral methods now mature for dynamical-spacetime coupled GRMHD with exascale capabilities (Li, 14 Aug 2025).

Taken together, global GRMHD simulations provide a comprehensive, predictive platform for connecting plasma microphysics, MHD turbulence, field generation, jet launching, radiative signatures, and large-scale astrophysical observables across all compact-object contexts.