Hybrid GRMHD+GRFFE: Jet Simulation Method

• The paper introduces a hybrid scheme that transitions between GRMHD and GRFFE based on local magnetization without using artificial density floors.
• The method applies a smooth mixing function and analytic inversion to accurately capture jet funnel dynamics and preserve physical fidelity.
• Implementation in the KORAL code shows enhanced stability and robust synthetic imaging of highly magnetized jets.

Hybrid GRMHD+GRFFE refers to a numerical methodology designed for black hole accretion simulations that stably unifies general relativistic magnetohydrodynamics (GRMHD) with general relativistic force-free electrodynamics (GRFFE). This scheme explicitly evolves standard GRMHD variables in low-magnetization ($\sigma$) regions and switches to a GRFFE-based evolution when $\sigma$ exceeds a chosen threshold, allowing accurate treatment of extremely magnetized (jet/funnel) zones without artificial density floors or unphysical magnetization ceilings. The hybrid approach is formulated to be modular, explicit, and compatible with finite-volume GRMHD code architectures, as demonstrated in the code KORAL (Chael, 2024).

1. Theoretical Framework

Standard GRMHD describes a magnetized relativistic fluid with eight primitive variables per cell, $P = [\rho, u_{\rm int}, \tilde u^i, B^i]$, comprising rest-mass density, internal energy, spatial momentum, and magnetic field. The evolved conserved variables are $U = \sqrt{-g}\,\alpha [D, Q_0, Q_i, B^i]$, where $D = \gamma \rho$, $\gamma$ is the Lorentz factor, and $Q_\mu = -\eta_\nu T^\nu{}_\mu$ with $$T^{\mu\nu}=T^{\mu\nu}_{\rm fluid} + T^{\mu\nu}_{\rm EM}$$.

In low-magnetization regions ($\sigma = b^2/\rho \ll 1$), the explicit finite-volume update of the conserved vector proceeds with standard GRMHD fluxes and sources, using the ideal, adiabatic GRMHD equations: \begin{align*} &\nabla_\mu(\rho u^\mu) = 0, \ &\nabla_\mu T^{\mu\nu} = 0, \ &\nabla_\mu {}^*F^{\mu\nu} = 0\,, \end{align*} wherein the stress-energy and Faraday tensor take their usual forms.

In regimes of high magnetization ($\sigma \gg 1$), the force-free limit ($\sigma$0, $\sigma$1) becomes more appropriate. GRFFE evolution uses primitive variables $\sigma$2, where $\sigma$3 is the drift velocity perpendicular to the magnetic field, and the stress-energy tensor depends solely on electromagnetic contributions. In this limit,

$\sigma$4

Inversion between conserved and primitive variables is analytic in GRFFE.

2. Transition and Mixing Strategy

The hybrid approach replaces the traditional GRMHD magnetization ceiling with a continuous transition from GRMHD to GRFFE governed by the cell-wise magnetization parameter,

$\sigma$5

where $\sigma$6. When $\sigma$7 exceeds a prescribed threshold $\sigma$8, the evolution transitions to GRFFE. No artificial density or internal energy floors are imposed, and the standard GRMHD upper limit $\sigma$9 is disabled.

The matching between regimes is achieved by the smooth mixing function

$P = [\rho, u_{\rm int}, \tilde u^i, B^i]$0

where $P = [\rho, u_{\rm int}, \tilde u^i, B^i]$1 sets the transition width. Pragmatically, pure GRMHD inversion is applied for $P = [\rho, u_{\rm int}, \tilde u^i, B^i]$2, and pure GRFFE for $P = [\rho, u_{\rm int}, \tilde u^i, B^i]$3. Intermediate values blend primitive variables as $P = [\rho, u_{\rm int}, \tilde u^i, B^i]$4.

In implementation, cell classification and mixing are recalculated at each substep of the time integration (Runge–Kutta). This avoids any discontinuous physical interface and minimizes numerical artifacts.

3. Augmented Fluid Evolution in Force-Free Regions

Although the electromagnetic field equations dominate in the high-$P = [\rho, u_{\rm int}, \tilde u^i, B^i]$5 zones, passive tracking of fluid quantities is preserved via advection: \begin{align*} &\nabla_\mu(\rho u^\mu) = 0, \ &\nabla_\mu(\rho s u^\mu) = 0, \end{align*} with entropy $P = [\rho, u_{\rm int}, \tilde u^i, B^i]$6 and $P = [\rho, u_{\rm int}, \tilde u^i, B^i]$7. The field-parallel velocity $P = [\rho, u_{\rm int}, \tilde u^i, B^i]$8 is evolved in the “cold” limit by enforcing

$P = [\rho, u_{\rm int}, \tilde u^i, B^i]$9

along fluid worldlines, with $U = \sqrt{-g}\,\alpha [D, Q_0, Q_i, B^i]$0.

These augmentations yield analytic expressions for primitive recovery in force-free cells and preserve the physical fluid entropy and parallel velocity, should $U = \sqrt{-g}\,\alpha [D, Q_0, Q_i, B^i]$1 drop below the transition threshold.

4. Numerical Realization and Coupling

The scheme modifies only the primitive and conserved arrays, augmenting them to include $U = \sqrt{-g}\,\alpha [D, Q_0, Q_i, B^i]$2 (with $U = \sqrt{-g}\,\alpha [D, Q_0, Q_i, B^i]$3). The evolution algorithm for each time update is as follows:

1. Compute $U = \sqrt{-g}\,\alpha [D, Q_0, Q_i, B^i]$4 and $U = \sqrt{-g}\,\alpha [D, Q_0, Q_i, B^i]$5 for all cells.
2. Evolve all conserved variables with the appropriate finite-volume update (GRMHD or FFE fluxes/sources).
3. Invert updated conserved variables to primitives using the mixing prescription and regime-appropriate inversion:
   • GRMHD: Newton–Raphson on $U = \sqrt{-g}\,\alpha [D, Q_0, Q_i, B^i]$6 as in Noble et al. (2006).
   • GRFFE: Analytic inversion via $U = \sqrt{-g}\,\alpha [D, Q_0, Q_i, B^i]$7, $U = \sqrt{-g}\,\alpha [D, Q_0, Q_i, B^i]$8.
4. Magnetic field is updated via constrained transport (FluxCT) to ensure $U = \sqrt{-g}\,\alpha [D, Q_0, Q_i, B^i]$9.

The explicitness, use of standard finite-volume updates, and analytic inversions enable direct integration into existing divergence-controlled GRMHD infrastructures with minimal code disruption.

5. Implementation in KORAL and Practical Considerations

In the KORAL code, the primitive set is extended to $D = \gamma \rho$0, and the conserved array to include the FFE-specific variables. KORAL uses cell-centered PPM reconstruction, local Lax–Friedrichs (LLF) fluxes, second-order Runge–Kutta, and FluxCT for divergence control. The explicit hybrid scheme imposes no additional global CFL constraint, and the runtime overhead is minor, since only the highly evacuated jet funnel contains cells exceeding $D = \gamma \rho$1. For numerical stability, a maximum Lorentz factor $D = \gamma \rho$2 is enforced.

6. 3D Simulation Outcomes and Synthetic Imaging

Comparative 3D simulations of magnetically arrested disks (MADs) around spinning black holes ($D = \gamma \rho$3) using both standard GRMHD and the hybrid scheme (grid: $D = \gamma \rho$4) demonstrate:

• Identical disc/corona structure for $D = \gamma \rho$5, with averaged dimensionless flux $D = \gamma \rho$6–$D = \gamma \rho$7.
• The standard GRMHD run caps jet funnel magnetization at $D = \gamma \rho$8 with artificially high $D = \gamma \rho$9 due to floors, while the hybrid run achieves $\gamma$0 and $\gamma$1 becomes minuscule in the funnel (plasma $\gamma$2).
• Synthetic 230 GHz imaging (M87* parameters) via ipole shows that standard GRMHD images depend sensitively on the ad hoc radiative transfer cutoff $\gamma$3, with variations above $\gamma$4 due to “floor haze.” The hybrid simulation removes this sensitivity for $\gamma$5, producing robust funnel images.

The physical evacuation of plasma from the funnel in the hybrid approach yields synthetic observables that more faithfully reflect high-$\gamma$6 jet conditions than those produced under traditional floor/ceiling prescriptions.

7. Performance, Stability, and Integration Guidelines

Key advantages of the hybrid GRMHD+GRFFE method include:

• Stable evolution at $\gamma$7, eliminating the need for unphysical density floors and magnetization ceilings.
• No extra global Courant–Friedrichs–Lewy restriction.
• Increased physical fidelity in funnel plasma evolution, critical for robust high-frequency synthetic imaging and jet modeling.

For porting to other divergence-controlled GRMHD codes:

1. Extend primitive and conserved arrays as in KORAL.
2. Implement analytic GRFFE inversions (as per McKinney 2006).
3. Add passive advection of $\gamma$8, $\gamma$9, and $Q_\mu = -\eta_\nu T^\nu{}_\mu$0.
4. Evaluate $Q_\mu = -\eta_\nu T^\nu{}_\mu$1 and consistently apply the mixing of primitives.
5. Retain existing constrained transport for $Q_\mu = -\eta_\nu T^\nu{}_\mu$2.

This procedure enables broad adoption in codes requiring accurate treatment of highly magnetized relativistic jets with minimal changes to existing infrastructure (Chael, 2024).