Hybrid GRMHD and Force-Free Simulations of Black Hole Accretion

We present a new approach for stably evolving general relativistic magnetohydrodynamic (GRMHD) simulations in regions where the magnetization $\sigma=b^2/\rho c^2$ becomes large. GRMHD codes typically struggle to evolve plasma above $\sigma\approx100$ in simulations of black hole accretion. To ensure stability, GRMHD codes will inject mass density artificially to the simulation as necessary to keep the magnetization below a ceiling value $\sigma_{\rm max}$. We propose an alternative approach where the simulation transitions to solving the equations of general relativistic force-free electrodynamics (GRFFE) above a magnetization $\sigma_{\rm trans}$. We augment the GRFFE equations in the highly magnetized region with approximate equations to evolve the decoupled field-parallel velocity, plasma energy density, and plasma mass density. Our hybrid scheme is explicit and easily added to the framework of standard-volume GRMHD codes. We present a variety of tests of our method, implemented in the GRMHD code KORAL, and we show first results from a 3D hybrid GRMHD+GRFFE simulation of a magnetically arrested disc (MAD) around a spinning black hole. Our hybrid MAD simulation closely matches the average properties of a standard GRMHD MAD simulation with the same initial conditions in low magnetization regions, but it achieves a magnetization $\sigma\approx10^6$ in the evacuated jet funnel. We present simulated horizon-scale images of both simulations at 230 GHz with the black hole mass and accretion rate matched to M87*. Images from the hybrid simulation are less affected by the choice of magnetization cutoff $\sigma_{\rm cut}$ imposed in radiative transfer than images from the standard GRMHD simulation.


INTRODUCTION
General Relativistic Magnetohydrodynamic (GRMHD) codes are now a standard tool for investigating accretion flows, jets, and outflows on black hole horizon scales.Most supermassive black holes accrete slowly, with an accretion rate  many orders of magnitude below the Eddington limit  Edd .The accreting plasma around these black holes is hot ( > 10 10 ), dilute, and optically thin (Yuan & Narayan 2014).These flows can be strongly magnetized and frequently launch jets by extracting the spin energy of the black hole (Blandford & Znajek 1977).
Millimeter-wave synchrotron emission from hot accretion flows around supermassive black holes is produced on horizon scales at radii  < 10  g , where  g =  / 2 is the black hole's gravitational radius.Using Very Long Baseline Interferometry (VLBI) at 230 GHz, the Event Horizon Telescope (EHT) has resolved and imaged this horizon-scale emission around the supermassive black holes with the largest apparent sizes: M87* (Event Horizon Telescope Collaboration et al. 2019b) and Sgr A* (Event Horizon Telescope Collaboration et al. 2022a).The EHT images of both sources have been extensively compared to simulated images from GRMHD simulation models made using polarized, relativistic ray-tracing radiative transfer codes (e.g.Dexter 2016;Moscibrodzka & Gammie 2018).Connecting EHT images to simulated data from GRMHD simulations has indicated that the emitting plasma around M87* co-rotates ★ E-mail: achael@princeton.edu with the black hole (Event Horizon Telescope Collaboration et al. 2019c) and that it is likely in the magnetically arrested (MAD) state of black hole accretion (Event Horizon Telescope Collaboration et al. 2021, 2023).EHT observations and simulations of Sgr A* also indicate that it is likely magnetically arrested (Event Horizon Telescope Collaboration et al. 2022b, 2024a,b).
The strongly magnetized regime is particularly interesting for understanding the physics of relativistic jet launching, plasma heating, and flaring processes around supermassive black holes.Unfortunately, numerical GRMHD codes struggle when the magnetic field energy density greatly exceeds the plasma rest mass density, or when the magnetization parameter  ≫ 1: Here,  2 is the magnetic energy density in the fluid rest frame in Heaviside-Lorentz units, 1 and  is the fluid mass density.In highly magnetized regions, the fluid rest mass and thermal energy make up a small contribution to the overall energy budget.Small numerical errors in the evolution of the overall energy-momentum can cascade to large errors in the fluid quantities, and the simulation can evolve to an unphysical state and crash.
The typical solution to the instability of GRMHD codes at high  1 Heaviside-Lorentz units are related to Gaussian units by a factor of 4  in the magnetic energy density, such that  2 HL =  2 G /4 .
is to introduce numerical "floors" on the mass density .In practice, this is usually achieved by limiting  <  max ≈ 100 throughout the simulation by constantly increasing the density in regions that would otherwise prefer to climb to higher magnetization.The details of this flooring procedure, such as which frame mass is added to, can affect the results of the simulation (Ressler et al. 2017).Furthermore, when producing simulated synchrotron images from GRMHD simulations for comparison to observations, it is necessary to remove emission from "floored" regions with an artificially high plasma density, or  >  max .In practice, most studies of GRMHD images more conservatively remove all emission from regions with  >  cut ≈ 1, as the evolution of the plasma internal energy and temperature is considered to be untrustworthy even below the ceiling value  max .The resulting millimeter-wavelength image structure can be sensitive to the chosen value of  cut (e.g.Chael et al. 2019;Zhang et al. 2024).This paper proposes an alternative approach to density floors in GRMHD simulations.Instead of artificially increasing the density to force the simulation to stay below a maximum magnetization  max , we present a method to smoothly connect the evolution of the GRMHD equations to their force-free limit in regions of high .The equations of general relativistic force-free electrodynamics (GRFFE) represent the dynamics of an infinitely conductive, zerorest-mass magnetized plasma in curved spacetime and are the limit of the GRMHD equations as  → ∞.The GRFFE equations are typically solved by evolving the magnetic and electric fields in some coordinate system (e.g.Komissarov 2002;McKinney 2006;Etienne et al. 2017;Mahlmann et al. 2021).However, McKinney (2006) showed that the GRFFE equations can be solved in standard finitevolume GRMHD codes with only small modifications, essentially removing the back-reaction of the fluid variables on the evolution of the magnetic field and plasma velocity.
In our approach, we adopt the method of McKinney (2006) to solve the equations of GRFFE in highly magnetized regions above a transition magnetization  >  trans .Since the GRFFE equations do not constrain the evolution of the fluid density and temperature, nor the fluid velocity parallel to the magnetic field lines, we solve additional equations to approximately evolve these quantities in the high- or GRFFE region.We then smoothly connect the GRFFE evolution region in the magnetized jet region of a black hole accretion simulation to standard GRMHD evolution in most of the simulation volume.
Other works have considered different methods for connecting GRMHD and GRFFE evolution in studies of neutron star magnetospheres.Paschalidis & Shapiro (2013) switch between GRMHD and GRFFE evolution at a fixed boundary, the neutron star surface; as they are interested primarily in the evolution of the neutron star magnetosphere, they do not attempt to evolve fluid quantities in this region.Parfrey & Tchekhovskoy (2017) take a different approach, adapting a GRMHD code as in this work to operate in the high  limit by manually fixing the fluid density  and internal energy  gas to pre-determined values in the magnetosphere and solving for the evolution of the electromagnetic field by the GRMHD equations on this fixed background.In our approach, we solve directly for the magnetic field, velocity, and plasma density and internal energy in both regions, under different equations.Like Parfrey & Tchekhovskoy (2017), we allow the boundary between the GRMHD and GRFFE regions to evolve dynamically throughout the simulation.We find that our approach can be straightforwardly applied to GRMHD simulations of MAD accretion discs, and that in these simulations we can stably evolve fluid with magnetizations  > 10 6 in the jet region close to the black hole.Our approach affects near-horizon images obtained from GRMHD simulations less than the standard floor ap-proach, as the density in the uncertain high- regions naturally falls to zero.
The paper is organized as follows.First, we review the GRMHD and GRFFE equations and their properties in section 2. Then in section 3 we review the standard finite-volume method for evolving the GRMHD equations, its extension to force-free evolution, and our new approach for coupling the two in different regions of a simulation.In section 4 we present several tests of our method, implemented in the GRMHD code KORAL (Sądowski et al. 2013).In section 5 we present a comparison between two 3D MAD simulations performed with a standard GRMHD method including density floors and our new hybrid scheme.We discuss and conclude in section 6.

Units and Definitions
Throughout, we use units normalized such that  =  = 1.We work in a fixed background metric   with signature (−, +, +, +), where we take the zeroth coordinate to be timelike at spatial infinity and we take Latin indices to refer to spatial components.Throughout, we use the ADM form of the metric: where the lapse , shift vector   , and spatial metric    are We frequently work in the normal observer frame   for a given coordinate system.The normal observer frame is defined by the covariant timelike vector   = (−, 0, 0, 0), so that We define a projection tensor    into the normal observer frame: Given a four-velocity   , we can then find the projected four-velocity in the normal observer frame ũ =      .The projected normalobserver four-velocity has components where the Lorentz factor  defined in the normal observer frame is We can also compute the Lorentz factor from the normal observer frame three-velocity ṽ = ũ / by the familiar expression where ṽ2 =    ṽ ṽ  . 2 The projection of the four-velocity into the normal observer frame, Equation 6, is invertible, so that given ũ or ṽ we can solve for  using Equation 8, and then determine the four-velocity from   = ũ +   .
An electromagnetic field is defined by its antisymmetric field strength (Faraday) tensor   , and its dual (Maxwell) tensor ★   . 3e decompose the field tensors into coordinate-dependent electric and magnetic field four-vectors in the normal observer frame: By definition, the contravariant normal observer electric and magnetic field vectors4 have no time component, E 0 = B 0 = 0, and they are both orthogonal to   .

GRMHD Equations
Here, we review the equations of ideal GRMHD following Gammie et al. (2003) and McKinney (2006).Ideal magnetohydrodynamics describes a single perfect fluid coupled to a degenerate electromagnetic field.The perfect fluid stress energy tensor is where   is the timelike four-velocity, ℎ is the relativistic enthalpy density, and  is the pressure, both measured in the fluid rest frame.
We adopt an adiabatic equation of state with index Γ such that  = (Γ − 1) gas , where  gas is the fluid internal energy.The enthalpy density ℎ is then where  is again the fluid rest-mass density.The dimensionless temperature is To take the ideal MHD approximation, we assume that the fluid's conductivity is infinite and hence the electric field is zero in the fluid rest frame,   =     = 0.The magnetic field in the fluid rest frame is generally nonzero and is given by which is orthogonal to   by the antisymmetry of ★   .The fluidframe magnetic field vector   is related to the normal observer frame magnetic field B  by a projection:5 The Faraday tensor in GRMHD has a nontrivial kernel; that is, there exists at least one frame   where     = 0, so that the electric field in that frame vanishes.In general, when   has a nontrivial kernel, the electromagnetic field is called degenerate.When at least some of the frames   in the kernel are timelike, the field is magnetically dominated.
When the electromagnetic field is both degenerate and magnetically dominated, the Maxwell tensor can be expressed simply in terms of the velocity   and fluid frame magnetic field   as The stress-energy tensor for a degenerate, magnetically dominated electromagnetic field is In the Appendix, section A, we review the forms of   , ★   , and   EM in terms of the normal observer frame fields.Given these definitions, the GRMHD equations of motion can be expressed simply as the conservation of the total stress-energy along with the conservation of the mass density current and the homogeneous Maxwell equation The inhomogeneous Maxwell equation ∇    = −  is not required to evolve the GRMHD system of equations, and it is instead taken to define the electric current   .

Force-Free Equations
When the electromagnetic stress-energy dominates over the fluid stress-energy,   EM ≫   fluid , the GRMHD equations for the electromagnetic field become independent of the fluid: These are the equations of force-free electrodynamics (GRFFE).
In both GRFFE and ideal GRMHD, the electromagnetic field is degenerate and magnetically dominated.One can show (e.g.McKinney 2006) that if there is at least one timelike frame   with vanishing electric field,   =     = 0, the electric field also vanishes in an infinite family of timelike frames.These frames are connected by Lorentz boosts along magnetic field lines, and Equation 16and Equation 17 are valid descriptions of the Maxwell and stress-energy tensors in any of these frames. 6In GRMHD, the fluid velocity picks out a unique   with vanishing electric field, but in GRFFE the evolution equations alone do not pick out a unique   .
Among the infinite number of frames with vanishing electric field, the drift frame   ⊥ is the unique frame with zero electric field where the Lorentz factor relative to the normal observer   is minimized: The Lorentz factor of the drift frame is It is clear from the form of Equation 24 that the field must be magnetically dominated for the drift velocity in the normal observer frame to be slower than the speed of light.Since all other frames with vanishing electric field have a larger Lorentz factor than the drift frame, the field must be magnetically dominated for any frame with vanishing   to be timelike.
Given a magnetically dominated degenerate field and a fixed coordinate system, the drift frame is the unique frame where the electric field vanishes, where the Lorentz factor relative to the normal observer is minimized, and which is orthogonal to the normal-observer magnetic field,  ⊥ B  =    ⊥ ★   = 0.Because the drift frame depends on the choice of normal-observer frame   , it is not a "physical" frame.Unlike the fluid frame in GRMHD, the drift frame and the drift frame Lorentz factor change depending on the choice of coordinate system.
All other frames with vanishing electric field can be parameterized by their three-velocity ṽ along the magnetic field line in the normal observer frame, which is encoded in the zeroth component of the magnetic field   : A general frame where the electric field   vanishes is thus related to the drift frame   ⊥ by a Lorentz boost along the normal observer magnetic field: To ensure that Equation 26 remains timelike, the parallel threevelocity has a maximum value Alternatively, we can parameterize the boost in terms of a parallel Lorentz factor  ≥ 1 such that the total Lorentz factor is In the Appendix, section B, we present a more detailed discussion of the field-parallel and field-perpendicular velocities.
The drift frame and the parallel boost Lorentz factor are coordinate-dependent quantities defined in terms of the normal observer   .In different coordinate systems for the same metric (e.g.Boyer-Lindquist and Kerr-Schild coordinates for the Kerr metric of a spinning black hole), identical degenerate EM fields will have different drift velocities, depending on   .When solving the forcefree equations 21 and 22 numerically in terms of a velocity   and magnetic field   , the drift frame velocity is typically used, since the GRFFE equations do not determine the velocity parallel to the magnetic field lines (McKinney 2006).

Finite Volume GRMHD Evolution
Finite-volume GRMHD codes define the fluid and magnetic field in terms of a vector of "primitive" variables P defined in each cell of a discretized grid.The eight primitive quantities in GRMHD are GRMHD codes usually use the normal observer frame velocity components ũ as a primitive quantity instead of   , since they can take any value from −∞ to ∞ and still produce a timelike   (McKinney & Gammie 2004).
The GRMHD equations in Equation 18-Equation 20 can then all be expressed in finite volume form where U are the "conserved" quantities, F are the fluxes between spatial cells, and S are source terms.The conserved quantities in GRMHD are where D =  is the density in the normal observer frame and Q  = −     =  0  is the energy flux four-vector in the normal observer frame.The fluxes are In ideal GRMHD without additional physics (e.g.coupling to radiation), the source terms S arise purely from the geometry; they are nonzero only for the stress-energy equations and involve products of the Christoffel symbols (see Gammie et al. 2003, Eq 4).
The GRMHD equations are hyperbolic and for a given coordinate system they may be evolved forward in the time coordinate  =  0 explicitly.A timestep Δ begins by interpolating the primitive quantities P to cell walls.The conserved quantities and source terms in each cell and the flux terms on each cell wall are then computed.The conserved quantities are updated by summing up the fluxes entering and leaving the cell, along with the geometric source terms, multiplied times Δ.
After the conserved quantities U are explicitly evolved forward in time with Equation 30, a GRMHD code must solve for the primitives P(U) in each cell.In GRMHD, the map from the conserved quantities U to the primitives P is not analytically invertable, so the conservative-to-primitive inversion must be done numerically.Noble et al. (2006) showed that the primitives P can all be expressed analytically in terms of the conserved quantities and the relativistic enthalpy  ≡  2 ℎ.Most GRMHD codes numerically solve for this single variable, , from the conserved quantities using a Newton-Raphson approach.We review the Noble et al. (2006) inversion procedure in the Appendix, section C Unfortunately, GRMHD numerical inversion can fail.When the magnetization  is large, the magnetic parts of the conserved quantities Q and U ≡ −  Q  dominate the contributions from the fluid energy momentum (Equation C7).As a result of truncation error, a GRMHD code can end up in a situation where no consistent solution for  can be found given the numerically evolved conserved quantities.GRMHD codes handle these failures by imposing artificial ceilings on , or equivalently floors on the density  These floors are handled differently in different codes, but they all have the effect of limiting the magnetization below some ceiling  <  max in some frame.Except for some very well-behaved problems,  max ≈ 100 in most applications in most GRMHD codes.
The region of parameter space where GRMHD conserved-toprimitive inversion fails, as  becomes large, is also where the GRMHD equations approach their force-free limit.Thus, instead of imposing floors on , we may hope to handle the failure of GRMHD codes at large  by transitioning to solving the GRFFE equations in regions of large .However, since the GRFFE equations only uniquely determine the drift frame velocity   ⊥ and its associated magnetic field   ⊥ , we will have to add additional equations to determine the evolution of the density, fluid internal energy, and parallel velocity.

Finite Volume GRFFE Evolution
Most force-free codes determine the evolution of the electromagnetic field   under the GRFFE equations Equation 22 using the electric and magnetic field components {B  , E  } as primitive variables.McKinney (2006) showed that it is possible to evolve the force-free equations of motion in the same framework as a standard GRMHD code, where the primitive variables are the magnetic field components   and drift frame velocity ũ ⊥ .In this approach to GRFFE, the primitive variables are P FFE = ũ ⊥ , B  .The associated conserved quantities, fluxes, and source terms are the same as in GRMHD, except that only  Unlike in GRMHD, the conserved-to-primitive inversion in this formulation of GRFFE is analytic.Specifically, McKinney (2006) showed that from the conserved quantities B  and Q  =  0  obtained from finite volume evolution, the normal observer electric field is and then the drift-frame velocity   C2).Namely, given the energy flux projected in the normal observer frame Q , the drift three-velocity is and the normal observer velocity is then ũ where  ⊥ is calculated with Equation 8.
Force-free electrodynamics can be implemented in this adapted GRMHD finite volume framework as long as the field remains magnetically dominated, B 2 > E 2 .In practice, a set of GRFFE initial conditions are not guaranteed to remain magnetically dominated as they evolve.In particular, in current sheets the field can become electrically dominated as the assumptions of GRFFE break down.We thus need to implement a ceiling on  ⊥ <  max to keep the field magnetically dominated and the evolution stable throughout the simulation region.Such a ceiling on  is also standard in GRMHD codes.

Decoupled GRFFE & Fluid Evolution
Next, we wish to define a strategy for fully evolving a fluid and electromagnetic field in the force-free limit, where the fluid does not back-react on the evolution of the field.Given a GRFFE solution for B  and   ⊥ , we can couple a non-interacting fluid to the force-free field by adding three evolution equations for the fluid mass density , the fluid energy density  gas , and the field-parallel velocity ṽ , which along with   ⊥ specifies the fluid velocity through Equation 26.Assuming we know   , to solve for the fluid density  we use the advection equation, as in GRMHD: To solve for the fluid energy density, we next assume that the fluid evolves adiabatically in the force-free region: In Equation 36,  is the entropy per unit mass: The corresponding conserved quantity for the entropy is S = D.
The adiabatic approximation minimizes the temperature of the fluid in the force-free regions when compared to a solution to the full GRMHD equations.Evolving the fluid adiabatically will not produce a physical solution to the GRMHD equations in the high  regions, and this choice ignores physical sources of dissipation that may exist in the force-free jet.The standard GRMHD approach of imposing a density floor in the jet tends to produce a jet interior that is relatively high-density and high-temperature; by contrast, the adiabatic GRFFE approach in this region tends to produce an evacuated, lower-temperature funnel (see section 5).
Advecting the density and entropy density with Equation 35and Equation 36 requires us to know both the drift velocity   ⊥ and the parallel velocity ṽ in Equation 26.One option is to keep ṽ = 0, minimizing the Lorentz factor of the fluid rest frame with respect to the normal observer.However, even when solving for ṽ⊥ from pure force-free electrodynamics, we may wish to solve for a nonzero ṽ .For instance, when coupling GRMHD and GRFFE equations in different regions of the simulation domain, we may have a nonzero ṽ on the boundary of the GRFFE region, so setting ṽ = 0 in the GRFFE region could introduce discontinuities in the velocity.Furthermore, the decomposition of the velocity into field-parallel and field-aligned components in GRFFE is coordinate-dependent, so zeroing ṽ in the GRFFE region will give different physical velocities when a simulation is run in different coordinate systems.
An exact equation for the evolution of the parallel momentum in GRMHD can be obtained by taking the dot product of the full GRMHD energy-momentum equation ∇    = 0 with the fluidframe magnetic field 4-vector   (see Camenzind 1986 and the Appendix, section D for a derivation): Equation 38 indicates that the parallel momentum along a magnetic field line is only changed by gas pressure gradients along the field line.The parallel momentum equation Equation 38can in principle be used to evolve the parallel velocity ṽ , since  0 =  ṽ √ B 2 / by Equation 25.Unfortunately, the gradient of pressure on the right of Equation 38 makes using Equation 38 to solve for ṽ in an explicit GRMHD code difficult.We can make progress by taking one of two limits.
In the adiabatic limit, we assume the entropy current is conserved (by Equation 36), and additionally that the entropy per unit mass is constant along field lines B     = 0.Under these assumptions, Equation 38 simplifies to where  is the enthalpy-per-unit-mass,  ≡ ℎ/.If we use Equation 39 together with the GRFFE equations to evolve the parallel velocity, we can obtain  0 as a conserved quantity along with the normal observer frame mass density D and entropy density S. To invert the system we then need to numerically solve for ṽ from the conserved quantities  0 , and  = S/D.We can do this numerical inversion using a similar Newton-Raphson method as in standard GRMHD codes (Noble et al. 2006), either iterating on ṽ directly or on  =  2 ℎ.
Equation 38 simplifies further in the cold limit, when we assume In the cold limit, we neglect acceleration of the parallel velocity from pressure gradients along the field line.In this approximation, we can evolve  0 directly as a conserved quantity.Fixing   ⊥ from the force-free equations, we can then exactly solve for ṽ from the known parameter (41) The cold limit thus provides an explicit prescription for obtaining ṽ by evolving one additional equation for  0 in addition to the force-free solution for B  and   ⊥ .Because of its computational simplicity and simple analytic inversion, we adopt the cold limit (Equation 40) in this paper.In the case of GRMHD simulations of black hole accretion, using the cold limit of the coupled GRFFE + parallel momentum + adiabatic equations in the high  jet region will not exactly solve the gas dynamics, but it should give a better approximation to the evolution of the gas density, temperature, and momentum in this region than in the standard treatment when density floors are imposed.

Hybrid Evolution
Finally, we propose a method to couple GRMHD to GRFFE in different regions of a finite volume simulation.In particular, we are interested in switching from solving the full GRMHD equations to solving the approximate system of GRFFE plus the parallel velocity and entropy advection equations described in subsection 3.3 in regions of high magnetization  ≫ 1.This method allows us to evolve the gas density to very small values with  ≫ 100, a typical "ceiling" value in GRMHD.
Our approach is conceptually simple.We keep the same 'primitive' expressions in each cell as in GRMHD (with the addition of the entropy per unit mass ): In addition to the standard GRMHD conserved quantities (Equation 31), we add additional "conserved" quantities for force-free evolution, corresponding to the three field-perpendicular momenta Q EM, =  0 EM, , the conserved quantity for the parallel velocity  0 , and the conserved entropy density S: If we use the cold limit for the parallel velocity (Equation 40), we set  = 1 in Equation 43.We explicitly evolve U forward in time in every cell using either the GRMHD or force-free equations, as appropriate for the given conserved quantity; B  is evolved identically for both sets of equations from the homogeneous Maxwell equation (Equation 20).Then, we determine which quantities to use in the inversion U → P based on the magnetization  of each cell during the last timestep.
If the magnetization was less than some critical value  <  trans at the last timestep, we use the standard GRMHD inversion procedure and obtain P MHD .If the magnetization was higher than some critical value  >  trans , we use the force-free conserved quantities to obtain the updated primitives P FFE following the method in subsection 3.2 and subsection 3.3.
In practice in a standard GRMHD torus simulation, the vast majority of the simulation cells are in the "GRMHD" regime and GRMHD U → P inversion takes place normally.In some cells concentrated in the simulation jet region,  ≫ 1, and we use the force-free inversion equations.We turn off density floors and magnetization ceilings entirely in the GRMHD region; whenever these would be needed, we change the equations and conserved quantities in the inversion U → P instead of injecting extra density.
Instead of switching discretely from one inversion scheme to another depending on , we can allow for a transition between the GRMHD and GRFFE solutions for the inversion by adding a mixing fraction  (): We determine  in the range [0, 1] by using the local value of  computed at the last time step in a given cell of the simulation.When  ≈ 0 the cell is inverted with GRMHD and when  ≈ 1 it is inverted with the GRFFE conserved quantities.We use a hyperbolic tangent transition function in ln , with a functional form: In Equation 45, the parameter  trans is the value of  where the transition from GRMHD to GRFFE is centered, and  controls the width of the transition.The limit  → 0 corresponds to an sharp transition exactly at  trans .
Using finite transition width between the GRMHD and GRFFE solutions for P(U) in Equation 44has the disadvantage of not exactly conserving energy and momentum under either the GRMHD or GRFFE equations in the transition region.However, it may be useful to smooth out any sharp discontinuities in the simulation from the boundary between the two regions.In practice, in GRMHD torus simulations, we find that  climbs so rapidly in the transition into the jet region that only a few cells along the jet sheath are not in either the limit  ≈ 0 or  ≈ 1.
For computational efficiency, and to ensure the bulk of the simulation either purely GRMHD or purely GRFFE, we use only GRMHD inversion when  < 1/64 and only force-free inversion when  > 63/64.This amounts to mixing the results on the two inversion procedures only between  low =  trans /   and  high =  trans   , where   = 3  7 /2 .We plot a typical curve  () with  trans = 50 and  = 0.1 in Figure 1.

IMPLEMENTATION AND CODE TESTS
Here we present details of our implementation of both GRFFE and hybrid GRMHD+GRFFE evolution in the code KORAL (Sądowski et al. 2013).We then present a number of 1D and 2D code tests of both pure GRFFE and GRMHD+GRFFE hybrid problems in both flat space and the Kerr geometry.

Implementation Details
We have implemented both the GRFFE solution method of McKinney (2006) and its extension to hybrid GRFFE+GRMHD evolution in the GRMHD code KORAL (Sądowski et al. 2013).KORAL was based on the HARM GRMHD code (Gammie et al. 2003;Noble et al. 2006), but is extensively modified for radiative GRMHD simulations of black hole accretion (Sądowski et al. 2014(Sądowski et al. , 2015) ) as well as two-temperature radiative simulations (Sądowski et al. 2017;Chael et al. 2017).The KORAL code has been extensively tested on standard GRMHD test problems (Sądowski et al. 2014) and compared to other GRMHD codes in the EHT code comparison project (Porth et al. 2019).KORAL works on a regular grid in arbitrary spacetime coordinates; all quantities (including the magnetic field) are cell-centered.In most applications, KORAL interpolates primitives P to cell walls with the second order piecewise parabolic method (PPM).Fluxes at cell walls are evaluated with the local Lax-Friedrichs (LLF) method, and evolution in time is performed with a second-order Runge-Kutta method.KORAL ensures the magnetic field remains divergence-free by using the FluxCT constrained transport algorithm (Tóth 2000).
In standard GRMHD applications, KORAL applies ceilings on the magnetization  in the normal observer frame, but an option to apply ceilings in the drift velocity frame (Ressler et al. 2017) is also implemented.A typical value of the ceiling magnetization value used in GRMHD evolution in KORAL is  max ≈ 100.In the GRFFE and hybrid GRMHD+GRFFE code tests presented here, we completely turn off the standard  max ceiling in KORAL and instead transition to force-free evolution at a given  trans , using either a sharp transition or applying a smoothing function  (), as in Equation 44.In hybrid GRMHD+GRFFE evolution, we update the local  () of a cell once per Runge-Kutta sub-timestep.As discussed above, when using a continuous mixing function  (), we apply upper and lower limits on ,  high and  low , above and below which quantities are evolved only according to the GRFFE or GRMHD equations, respectively.In all applications we must retain a ceiling on the maximum Lorentz factor; we typically set  max = 100 unless otherwise specified.

GRFFE Tests
We first test our implementation of pure GRFFE evolution in KORAL with the standard set of 1D test problems in flat spacetime introduced in Komissarov (2002).These tests consist of (1) a nonlinear fast wave propagating to the right at the speed of light,  = 1; (2) an nonlinear Alfvén wave propagating left with speed  = −0.5;(3) a nonlinear degenerate Alfvén wave propagating to the right with speed  = 0.5; (4) a superposition of a stationary Alfvén wave and two nonlinear fast waves propagating left and right at  = 1; (5) a "breakdown" test that evolves to a state where the field becomes electrically dominated.See Paschalidis & Shapiro (2013) for clear and detailed implementations of the initial conditions in all five problems.
The spatial domain is  ∈ [−2, 2] for the fast wave and Alfvén tests,  ∈ [−1, 1] for the three wave test, and  ∈ [−0.5, 0.5] for the breakdown test.We used outflow boundary conditions in all cases, and a fiducial resolution of   = 256 in all tests.We ran all tests to a maximum time of  max = 1.5.
Because these tests are pure FFE problems, we turned the parallel momentum solver off and only evolved the drift velocity   ⊥ and magnetic field B  .The initial density and internal energy are advected adiabatically, but do not back-react on the magnetic field and velocity.We used first-order spatial reconstruction with a monotonized central (MC) limiter instead of second-order PPM reconstruction, as first-order reconstruction is often more robust for propagating sharp discontinuities (Gammie et al. 2003).
We present our results for all five test problems in Figure 2; these results may be compared to results from other GRFFE codes, including Figure 1 of McKinney (2006) and Figure 1 of Paschalidis & Shapiro (2013).For the fast wave, Alfvén and three-wave tests, we plot a given component of the electric or magnetic field evolved by KORAL over the analytic solution at a given time; in all cases, the numerical solution lines up well with the analytic expectation.In the breakdown test (which does not have an analytic solution) we plot  2 −  2 , indicating the degree of magnetic domination of the electromagnetic field.As reported by Komissarov (2002), we see that for the given initial condition the electric field strength approaches the magnetic field strength at  ≈ 0.02, after which the FFE solver must begin applying a ceiling on the Lorentz factor .

Hybrid GRMHD+GRFFE tests
Next we turn to testing our implementation of matching the GRFFE solver in part of the spatial domain with standard GRMHD in the rest of the domain.For a GRMHD problem (like the Bondi problem in subsubsection 4.3.2), the GRFFE solution will only be approximate and should match the expected GRMHD solution only in regions of high ; similarly, for a GRFFE problem setup, we expect any GRMHD solution to diverge from the GRFFE solution in regions of low .Thus, how well any hybrid solution matches an analytic result will depend on the choice of  trans and desired error tolerance.In these tests, we primarily aim to show that our implementation of hybrid GRMHD+GRFFE does not introduce unexpected artifacts at the boundary between the GRFFE and GRMHD domains.

Linear Waves
We first consider an Alfvén wave in flat space in the linear regime, moving rightward from a region where we solve the MHD equations  Komissarov (2002).For each problem, dashed lines represent the initial solution, solid lines indicate the analytic solution at the given time, and symbols indicate the numerical solution for the reported magnetic or electric field component.In the breakdown test (bottom row), there is no analytic result, but our solution closely matches the result from Komissarov (2002) indicating that that at time  = 0.02 the electric field strength is nearly equal to the magnetic field strength.
in a region of strong magnetic field into a region where we solve the approximate set of FFE and decoupled parallel momentum, density and entropy equations.For our initial conditions, we fix a background B-field  0 = 1 along the −direction.We fix the magnetization  = 250, so  =  2 0 / = 0.4, and we fix the initial temperature Θ gas = 0.2 and adiabatic index Γ = 4/3.
Alfvén waves propagating on this MHD background have a velocity   = √︁   /(1 +   ) = 0.996, where   =  2 /ℎ = 138.8 is the magnetization defined relative to the enthalpy.Our initial pulse propagates to the right with this velocity in the MHD region and at the speed of light in the FFE region.We define our initial pulse in  48.
the domain −0.5 <  < 0.5: while for  > 0.5 and  < −0.5,   =  0 and   =   = 0. We fix  = 10 −3 .The drift velocity of the initial pulse is Rather than determine which solver we use based on the local , in this problem we transition from GRMHD to GRFFE based on the spatial coordinate .We thus replace  () in Equation 44with  () defined by where  trans = 0.5 and  = 0.1.Again, we use MHD inversion exclusively when  < 1/64 and FFE inversion exclusively when  > 63/64.Our spatial domain is  ∈ [−2, 2], and we use a fiducial resolution of   = 512.We run the simulation to a maximum time of  max = 1.5.In the FFE solution, we use the cold approximation to determine the parallel momentum.We use second-order PPM spatial reconstruction, and limit our maximum Lorentz factor to  <  max = 1000.We present our results in Figure 3, where we plot the − and − components of the electric and magnetic field (orthogonal to the background magnetic field in the −direction) after evolving the system to  max = 1.5.The numerical solution lines up well with the analytic expectation of a propagating linear pulse and does not show any artifacts from passing through the transition between the MHD inversion and FFE inversion at  = 0.5.

Bondi Accretion
In this test we simulate spherically symmetric Bondi accretion onto a Schwarzschild ( * = 0) black hole, which is a standard test problem in GRMHD codes (e.g. de Villiers et al. 2003;Porth et al. 2017).We set up an analytic Bondi flow in Boyer-Lindquist coordinates with a sonic point at  c = 8  g , an accretion rate  = −1 in code units, and adiabatic index Γ = 4/3.We then solve for the initial density, internal energy, and velocity for the Bondi solution as a function of radius in the equatorial plane numerically following the description in Rezzolla & Zanotti (2013), section 11.4.
A magnetic field parallel to the radial infall does not affect the Bondi solution (Porth et al. 2017).We set up an initial radial magnetic field, Where the velocity and density at the sonic point are From top to bottom, we show profiles of the plasma density , Boyer-Lindquist radial velocity   , dimensionless temperature Θ gas , and magnetization  for the spherically symmetric Bondi accretion setup described in this section.The solid black line indicates the initial condition; the Bondi problem should remain stationary in time.Blue circles show the numerical solution after 100  g for a standard GRMHD method, grey triangles show the results from a GRMHD simulation which evolves the fluid energy adiabatically, and orange xs show the numerical solution at the same time from our hybrid approach.In the hybrid approach, we transition between GRMHD and GRFFE solutions at  trans = 50, indicated by the horizontal dashed line in the bottom panel; this transition magnetization corresponds to a stable transition radius of  ≈ 6.3  g , indicated by the vertical dashed line in all panels.
We fix the magnetization at the sonic point to be  c = 25.The magnetization at the horizon climbs to  ≈ 1000.
The Bondi solution should remain stationary in time.We evolve the initial conditions numerically in 1D in modified Kerr-Schild (KS) coordinates in the domain / g ∈ [1.8, 100].We use a logarithmic grid in ; the resolution is   = 128 in our fiducial test.While we run the simulation in KS coordinates, we present our results in Boyer-Lindquist coordinates in Figure 4. We use outflow boundary conditions and second order PPM reconstruction, and we run the problem to  max = 100  g .
In Figure 4 we compare results from our new hybrid GRMHD+GRFFE method to a standard GRMHD run of the same problem using KORAL.When testing the hybrid method, we fixed the transition point at  trans = 50 in Equation 45, with width  = 0.05.This choice puts the transition around  ≈ 6.3  g in our problem setup, indicated by the vertical dashed line in Figure 4.For the given set-up, standard GRMHD is well behaved even up to  = 1000 at the horizon, so we we do not implement any ceiling on  in the GRMHD comparison run.This problem tests the parallel momentum solver in the force-  52) between the solution at  = 100  g and the analytic solution for the density , internal energy  gas , and Boyer-Lindquist radial velocity   as a function of the number of radial cells   .We plot the  1 error for the standard MHD run in blue circles and the error for the hybrid GRMHD+GRFFE run in orange triangles.Results from the adiabatic GRMHD simulation are shown in grey.
free region, as all of velocity is parallel to the radial magnetic field line and hence   ⊥ = 0.Because the dimensionless temperature is relatively high close to the horizon, Θ gas ≈ 0.1, we use the adiabatic approximation (Equation 39) instead of the cold approximation (Equation 40) for solving for the parallel velocity in the force-free region.
In Figure 4 we show the initial conditions and the GRMHD and hybrid GRMHD+GRFFE solutions after evolving the system for  = 100  g for the plasma density , radial four-velocity   , temperature Θ gas and magnetization  from our fiducial run with resolution   = 128.The standard MHD solution begins to deviate from the analytic stationary solution close to the horizon, as  > ∼ 100; by contrast, the hybrid solution more closely tracks the analytic expectation in this region as the code switches from solving the GRMHD to GRFFE equations, including the approximate equations for the parallel velocity and adiabatic internal energy evolution.In particular, the MHD solution under-predicts the infall velocity and over-predicts the plasma temperature Θ g in the high magnetization region close to the horizon, while the hybrid solution more closely tracks the ground truth.We also show results from an GRMHD run where the conserved entropy S is used instead of the energy Q 0 ; when adiabatic evolution is enforced, the standard GRMHD run recovers the gas temperature significantly more accurately (Porth et al. 2017), but it still performs worse than the hybrid run in recovering  and   .
In Figure 5 we perform a convergence test for both methods.We plot the normalized  1 error for the density , internal energy  int and velocity   of the numerical solution at  = 100  g compared to the stationary ground truth initial condition.The  1 error for a given quantity (, ) is defined, where we sum over all radial points   in the domain.Figure 5 indicates that the accuracy of both the standard GRMHD and hybrid GRMHD+GRFFE approaches converges at second order in the spatial resolution, but the absolute error of the hybrid set-up is lower for all three quantities than the MHD solution.As seen in Figure 4, the standard MHD approach introduces particularly large error in the internal energy  int when compared to the error in the density or velocity.By contrast, the hybrid approach features a similar magnitude of the relative error in all three quantities for a given resolution   .
The GRMHD simulation run with adiabatic evolution enforced has an absolute error in all three quantities that is between the standard GRMHD and hybrid simulations.The GRMHD solution is stable at the given magnetization  c = 25 and does not need to impose ceilings on the magnetization close to the horizon.However, if we increase the magnetization further beyond  c = 100, standard GRMHD begins to break down close to the horizon, while our new hybrid method remains well-behaved.

BZ Monopole
For our final test problem we simulate a 2D force-free monopole in the Kerr metric for a black hole with dimensionless spin  * = 0.5.We set up an initial radial monopolar field in the Kerr spacetime.As the field evolves under the GRFFE equations, it is quickly spun up by the black hole and develops a significant toroidal component.The wound-up magnetic field achieves an angular velocity Ω  = 0.5Ω  , where Ω  = /2 + is the angular frequency of the horizon, thus maximally extracting energy from the black hole spin by the Blandford & Znajek (1977) process.The BZ monopole is a standard test problem for GRFFE codes (e.g.Komissarov 2002;McKinney 2006;Paschalidis & Shapiro 2013;Mahlmann et al. 2021); it can also be evolved to relatively high  in GRMHD codes (e.g.Tchekhovskoy et al. 2009).
We work in Kerr-Schild coordinates (, , , ) in two dimensions.Our initial monopolar magnetic field is set by a toroidal vector potential   : The initial radial magnetic field in the normal observer frame is B  =     / √ −.We set up an initial density profile (): where  init = 1000.The atmosphere is initially isothermal with dimensionless temperature Θ gas = 10 −4 , and the adiabatic index is Γ = 4/3.The initial velocities in the normal observer frame are all zero.In the left plot, we show angular slices of the covariant toroidal magnetic field   =    ★  0 at radii  = 2  g (grey),  = 5  g (blue), and  = 10  g (orange) for the hybrid monopole simulation shown in Figure 6 at  = 50  g .In the right panel, we show profiles of the fieldline angular frequency Ω  at the same radii, normalized by the horizon angular frequency Ω  .The analytic solution from Blandford & Znajek (1977) for each radius is plotted in a dotted line.
We evolve the problem forward in time on a regular 2D grid in modified Kerr-Schild coordinates (,  1 ,  2 , ), defined in KORAL by For this test we use  0 = 0 and ℎ 0 = 0.8.We use   = 192 radial cells in the range  ∈ 0.7 + , 250  g and   = 128 cells in polar angle between  2 ∈ [0.001, 0.999].We run the simulation to a maximum time  max = 50  g .
For this problem, we compare a simulation run with pure GRFFE to a hybrid simulation that transitions from GRFFE to GRMHD evolution at   = 100, with a width  = 0.1 in Equation 45.In the GRMHD region of the hybrid simulation, we enforce entropy conservation (replacing Q 0 with S in the GRMHD conserved-to-primitive inversion, section C).In both simulations, in the FFE region we solve for the field-parallel momentum under the cold approximation, Equation 40.We found that using the parallel momentum solver is essential even in the inner GRFFE region of the hybrid simulation to correctly match the boundary conditions set by GRMHD evolution in the outer region.In both simulations we use outflow boundary conditions and first-order spatial reconstruction with monotonized central (MC) limiter.We set the maximum Lorentz factor  max = 2000.
As the monopolar field winds up and begins extracting energy from the black hole, it launches an outflow, driving plasma from the central region and expanding the transition radius where  =  trans = 100.In Figure 6 we show 2D poloidal slices of the magnetization , normal observer Lorentz factor  in KS coordinates, and field-line angular frequency Ω  for both the GRFFE and hybrid GRFFE+GRMHD simulations.For an axisymmetric degenerate electromagnetic field, the fieldline angular frequency can be computed from the plasma velocity and magnetic field components as The results shown in Figure 6 are not identical between the completely force-free simulation and hybrid GRMHD+GRFFE simulation.The differences are most pronounced at the transition region between GRMHD and GRFFE evolution at the  trans = 100 contour, which slightly lags the outer radius where the field lines are wound up by the black hole (the outer radius with Ω  ≈ 0.5Ω  ).The pure force-free simulation achieves a slightly higher Lorentz factor than the hybrid simulation, and the  = 100 contour is at a slightly smaller radius in the force-free simulation at  = 50 g .Compared to the force-free simulation, the hybrid simulation has a slightly ( < ∼ 5%) smaller radial velocity at and just outside the transition surface but nearly the same radial velocity deeper in the force-free region; as a result, some gas piles up outside the expanding wind front in the hybrid simulation, slightly lowering  compared to the pure force-free result.
Despite these differences, the overall structure of all three quantities in both simulations is similar, particularly interior to the transition radius in the force-free region at  ≈ 25  g .Thus, the boundary between GRMHD evolution and GRFFE evolution in our hybrid scheme does not adversely affect the force-free evolution inside the expanding region of wound-up field lines in this set-up.
In Figure 7 we show angular profiles of the poloidal covariant toroidal field   and field line angular speed Ω  for the hybrid simulation compared to the analytic Blandford & Znajek (1977) solution at different radii.The simulation results closely track the analytic solution; in particular, the fieldline angular speed achieves the expected optimal value for energy extraction, Ω  = 0.5Ω  .

FIRST 3D MAD SIMULATIONS IN GRMHD AND HYBRID GRMHD+GRFFE
To test our new method for hybrid GRMHD+GRFFE simulations on a real black hole accretion problem, we performed two 3D simulations of magnetically arrested (MAD; Bisnovatyi-Kogan & Ruzmaikin 1976;Narayan et al. 2003) accretion discs around black holes with dimensionless spin  * = 0.5.One simulation ("MHD") was performed with the standard GRMHD approach in KORAL, with a ceiling on the magnetization  max = 50 imposed in the normal observer frame (McKinney et al. 2012).The other simulation ("Hybrid") was performed with our new hybrid approach with a transition between GRMHD and GRFFE evolution at  trans = 50 and with transition width  = 0.1.The setup of the other numerical floors and simulation parameters were otherwise identical in both simulations.In Appendix E, we examine the effect of choosing higher and lower values of  trans in hybrid simulations with the same parameters.60) in the GRMHD and Hybrid MAD simulations from  = 0 to  = 10 4  g .

Simulation setup
The parameters of the both simulations are reported in Table 1.Both simulations were conducted on a uniform grid in modified Kerr-Schild coordinates similar to those used in the monopole test (Equation 55).The inner edge of the simulation grid is at  = 1.5  g inside the event horizon  + ≈ 1.87  g ; the outer edge of the grid is at  = 1000 g .The resolution of the simulations were 160 × 128 × 98 in radius, polar angle, and azimuth.Both simulations assume a gas adiabatic index Γ = 13/9.
We initialized both simulations with a thick torus of plasma in hydrodynamic equilibrium following Fishbone & Moncrief (1976).The torus inner edge is at  in = 20  g and the maximum pressure is at  ,max = 42  g .The maximum density is set to  max = 1 in code units; as the GRMHD equations are scale free, this normalization has no effect on the evolution.The initial torus is threaded with a single loop of poloidal magnetic field from the vector potential (Narayan et al. 2022): where  mag = 400  g .The magnetic field in the initial loop is scaled such that the ratio of the maximum thermal pressure to maximum magnetic pressure in the disc is  max = 100.To initiate accretion, we seed the initial torus with Gaussian perturbations in the pressure with a fractional standard deviation of 2 percent relative to the equilibrium Fishbone & Moncrief (1976) value.
Both simulations used second order PPM spatial reconstruction and LLF fluxes.Outflowing boundary conditions were used at the inner and outer radii, and reflecting boundary conditions are imposed at the polar axes.To control numerical instability from material reflecting off of the polar axes, we smoothly interpolate the poloidal velocity   to zero across the two cells closest to the axis.

Simulation results
We ran both simulations to a final time  final = 10000  g and saved snapshot files every 10  g .In Figure 8, we show the evolution of the mass accretion rate and magnetic flux on the horizon from the beginning to the end of the simulation.At every time  we we calculate the mass accretion rate  and magnetic flux through the horizon Φ  by integrating over the horizon at  =  + : Note that we use   = ★   = B  / in calculating the magnetic flux.Given these quantities, the dimensionless magnetic flux or "MAD parameter" is We use the factor √ 4 to convert the magnetic field back to Gaussian units in computing .Under this convention, MAD accretion discs tend to saturate at a normalized magnetic flux  ≈ 50 (Tchekhovskoy et al. 2011).
Figure 8 shows that both the standard GRMHD and hybrid simulations take approximately 5000  g from the start of the simulation to reach a steady-state in both  and .In Table 1, we report the averaged values of  and  for both simulations over the second half of both simulations 5000 ≤ / g ≤ 10000.Both simulations saturate with an average value of  ≈ 55 and differ by ≈ 7%.
Notably, when the magnetic field is first building up on the black hole horizon around  = 5000  g , the two simulations show somewhat different behavior.The hybrid simulation climbs to higher values of the magnetization  up to  ≈ 100 before its first flux eruption event, and then settles down to steady state for the remainder of the simulation.In a future work we will investigate higher resolution hybrid simulations and run them to longer times to see if larger values  in the initial phase of MAD simulations is a generic feature of our hybrid method.
In Figure 9, we show snapshot data in the poloidal plane from the two simulations of three quantities: the density  (in dimensionless code units), the magnetization , and the ratio of the thermal to magnetic pressure , defined as In Figure 10  In all images, the cyan contour indicates the critical curve or "black hole shadow" (Bardeen 1973); the magenta curve indicates the image of the equatorial event horizon, or "inner shadow" feature (Chael et al. 2021).All images are displayed in a log color scale.
averaged in azimuth and over the time interval 7500 ≤ / g ≤ 10000.
In plotting the averaged profiles of  and , we average the magnetic pressure, density, and thermal pressure independently and then take the appropriate ratios in Equation 1 and Equation 61.
Both Figure 9 and Figure 10 show that the global structure of the two simulations within  < ∼ 50 g is substantially similar, particularly for low  regions in the bulk of the accretion disc and corona region.Inside the magnetized jet, the density in the hybrid simulation is allowed to fall to much lower values and the magnetization correspondingly rises to a maximum  ≈ 10 6 close to the black hole; in the standard GRMHD simulation, the magnetization  < ∼ 50 everywhere because our ceiling limits the magnetization below this value. 7n addition to lowering the overall value of  in the jet region, the standard GRMHD simulation also reaches much lower levels of  in this region and develops a signature halo of high density  around the black hole purely from floor material (upper left panel of Figure 10).In Figure 11, we directly compare  and  on different contours of constant radius  in the averaged data from the two simulations.The average values of  and  in the simulations match well at  = 100  g and around the equatorial plane at smaller radii  = 2  g and  = 10 g ; at the smaller radii close to the poles, the GRMHD simulation values of  and  level off at around 50 and 0.01, respectively.By contrast, in the hybrid simulation, the average  can rise to 10 6 near the black hole, with  falling correspondingly to ≈ 10 −6 .
In Figure 12 we show time and azimuthal average profiles of the dimensionless plasma temperature Θ gas from both simulations.GRMHD simulations typically show very high temperatures Θ gas ≈ 1 in the jet region, which we also see in our comparison GRMHD run.The core of the jet in our hybrid simulation is kept cooler by our imposition of adiabatic fluid evolution in the force-free region of the simulation; Θ gas ≈ 0.1.The exact profile of temperature in the transition region near  = 1 between the jet and disc can  (Ressler et al. 2017).As in standard GRMHD, our hybrid method shows a significant increase in the plasma temperature as  increases, with temperatures near the  = 1 surface slightly exceeding the values seen in standard GRMHD.The jet core is slightly cooler in the hybrid simulation than in the standard GRMHD approach.

230 GHz images and lightcurves
We next consider the implications of our proposed method for hybrid GRMHD+GRFFE simulations for simulated observables of black hole accretion flows.We produce simulated 230 GHZ images of M87* from both of our simulations using the GR radiative transfer code ipole (Moscibrodzka & Gammie 2018) over the time period 7500  g <  < 10000  g .In producing these images set the temperature of the emitting electrons using the standard Mościbrodzka et al.
(2016)  high - low prescription.The ion-to-electron temperature ratio  is taken as a function of the local : We fix  low = 1,  high = 20 in the images shown here.Given , the number density   and temperature   of the emitting electrons are calculated from the simulation mass density  and internal energy  gas as We assume the adiabatic index of electrons is Γ  = 4/3 and the adiabatic index of the ions is Γ  = 5/3.Given the electron number density and temperature, as well as the magnetic field strength and orientation, ipole solves the radiative transfer equation for synchrotron emission along the curved geodesic trajectories of light around our  * = 0.5 Kerr black hole to produce a simulated image.
Because GRMHD simulations are scale-free, we scale the black hole mass and accretion rate in our simulations to values appropriate for M87* (Event Horizon Telescope Collaboration et al. 2019a).We take  = 6.5 × 10 9  ⊙ and set the distance to M87* as  = 16.8Mpc (Event Horizon Telescope Collaboration et al. 2019d).We fix the inclination angle of the spin axis to   = 163 deg (Walker et al. 2018).We then scale the mass accretion rate independently in each simulation so that the median flux density at 230 GHz over all of the images we raytrace is 0.5 Jy.The field of view of our images is 160as (≈ 43  g ), and we use a pixel size of 0.5as.Synchrotron images from GRMHD simulations typically do not include contributions to the emission or absorption from plasma with a magnetization  >  cut , where typically  cut = 1 (e.g.Prather et al. 2023).The value of  cut is typically set lower than the simulation ceiling value  max ≈ 50−100 for two reasons.First, we wish to avoid contributions to the simulated image from regions where the plasma density is dominated by floor material.Second, it is frequently a concern that the plasma properties from GRMHD simulations may not be reliable even in regions of intermediate magnetization 1 <  <  max where the simulation does not require floors for numerical stability.In this work we follow Chael et al. (2019) in setting a fiducial  cut = 25 instead of  cut = 1.We find that in both simulations there is significant emission from intermediate magnetization regions 1 <  < 25 along the jet sheath, but there is no significant emission from the jet core at  > 25 in the hybrid simulation.
For both the GRMHD and hybrid simulation, we produce images with  cut = 1,  cut = 25, and  cut = ∞, where in the latter case we include synchrotron emission and absorption from all material in the simulation domain.In Figure 13 we show a comparison between image snapshots from both simulations at all three values of  cut .It is evident that the hybrid simulation shows no change in image structure between  cut = 25 and  cut = ∞; by contrast, in the MHD simulation the  cut shows a significant contribution (∼ 20% of the total flux density) from the floor material in the forward jet, forming a haze in front of the black hole shadow and inner shadow (Chael et al. 2021) features.In Figure 14 we show simulation lightcurves over the range 7500  g <  < 10000  g from both simulations at both values of  cut .For the hybrid simulation, the lightcurves for both  cut = 25 and  cut = ∞ are nearly identical, while the overall flux density of the GRMHD lightcurve is increased with the addition of floor material when increasing  cut from 25 to ∞.In addition to allowing GRMHD simulations to stably evolve Poynting dominated jets to high values of , our proposed hybrid GRMHD+GRFFE simulation method may remove a potential source of systematic error in producing simulated images from GRMHD simulation by removing the requirement to set a  cut value when performing radiative transfer.

DISCUSSION AND CONCLUSIONS
In this paper we have presented a proposed approach for stably evolving GRMHD simulations of black hole accretion to high values of the magnetization  by switching to an augmented set of force-free equations in this region.Our method evolves the force free equation in the finite volume formulation of McKinney (2006); it solves for the field-parallel velocity, plasma density, and energy density in the magnetically dominated region with an approximate set of decoupled equations that do not back-react on the electromagnetic field evolution.We propose a method for joining force-free evolution to standard GRMHD evolution at a specified transition  trans , either via a sharp transition or a smooth average of the recovered GRMHD primitives across a certain range of  values.
We tested our new hybrid GRMHD+GRFFE code on a set of test problems and showed that it can accurately evolve a highly magnetized plasma across the interface  trans ; in the Bondi test problem in particular (subsubsection 4.3.2),we see a smaller error in the plasma density, energy density, and velocity evolved under our hybrid scheme than in a standard GRMHD approach.We then compared our hybrid method with a standard GRMHD set-up in a magnetically arrested simulation of an accreting black hole with the same initial conditions and resolution.Our hybrid method produced a MAD disc with average properties similar to standard GRMHD in low magnetization regions, but the absence of hard ceilings on  allows the jet funnel to evacuate and the magnetization to reach  ≈ 10 6 .Unlike the GRMHD simulations which are contaminated by floor material in high-magnetization regions, the hybrid simulations produce simulated submillimeter synchrotron images that are insensitive to the  cut parameter above  cut ≈ 25.In other words, our hybrid method naturally produces an evacuated force-free jet funnel in GRMHD simulations that does not contribute to the synchrotron emission simulated for comparison to EHT images.
Our proposed hybrid method is a relatively straightforward addition to the framework of finite-volume GRMHD codes.Because the conserved-to-primitive inversion P FFE (U) for the force-free quantities is analytic in the cold approximation that we use, and because we only perform the additional GRFFE conserved to primitive inversion step in a small part of the simulation domain, our hybrid scheme does not significantly increase the computational expense of standard GRMHD at the resolution we considered for our 3D MAD simulations.The primary additional computational cost comes from computing fluxes for the additional P FFE quantities at each cell walls.In our current implementation we compute these F(P FFE ) everywhere, but in the future we plan to increase the simulation efficiency further by limiting the domain where we compute F(P FFE ) only to regions where they are required in a given time step.
Our method is not the first to consider matching GRMHD to GRFFE in high magnetization regions in a single simulation volume.Notably, Paschalidis & Shapiro (2013) considered hybrid GRMHD and GRFFE simulations of neutron star magnetospheres.They use a similar finite-volume approach motivated by McKinney (2006) to the method presented in this paper to evolve the force-free region.However, they do not evolve the field-parallel velocity or fluid quantities in the force free region, and they fix the transition between force-free and GRMHD evolution to the surface of the neutron star.By contrast, in our proposed method the boundary between the GRMHD and GRFFE simulation regions evolves with the simulation and propagating jet from the black hole; plasma is allowed to flow in between the two regions by solving the parallel velocity and continuity equations.
Later, Parfrey & Tchekhovskoy (2017, 2023) conducted 2D and 3D hybrid neutron star magnetosphere simulations using a GRMHD code adapted for force-free evolution in the high magnetization region.In the simulations of Parfrey & Tchekhovskoy (2017, 2023), the GRMHD equations are evolved everywhere, but the code damps the field-parallel velocity in the magnetosphere region.The region where the code transitions to force-free behavior is determined by a combination of an advected passive scalar distinguishing the initial accretion flow and magnetosphere and a function of the coordinates that goes to zero outside the light cylinder.Recently, Phillips & Komissarov (2023) presented a novel operator splitting method in special relativistic MHD that solves for the MHD evolution of the velocity and magnetic field over a timestep as a correction to the value obtained from force-free evolution.This method can stably evolve SRMHD problems to high magnetizations without requiring a pre-defined transition  trans between MHD and FFE regions, but it has not yet been applied to 3D GRMHD simulations.
By coupling GRMHD to GRFFE, the method proposed in this paper can stably evolve jets in GRMHD simulations of black hole accretion to a magnetization four orders of magnitude larger than the largest  achieved in standard GRMHD.Other techniques or improvements to the method presented here may further increase the reliability of GRMHD in the jet region and open up a range of new questions to investigation using GRMHD simulations, including; do evacuated jets show more pronounced limb brightening (e.g.Lu et al. 2023) in simulated VLBI images than standard GRMHD images (Chael et al. 2019;Fromm et al. 2022)?Do GRMHD jets evolved without density floors reach higher Lorentz factors or have significantly different shapes when evolved for long times, and how do they compare to observation of M87* (Nakamura et al. 2018;Park et al. 2019)?Can we add self-consistent pair-production (e.g.Broderick & Tchekhovskoy 2015;Wong et al. 2021) to simulations of GRMHD jets to fill in the core and investigate the importance of pairs to EHT images?As GRMHD simulations of black hole accretion become more sophisticated, increasing in resolution and adding additional radiative and thermodynamic physics, it is essential to revisit the standard approach to evolving the MHD equations in magnetically dominated regions.
Then, plugging Equation C5 into Equation C2 and Equation C4, we can show that Q =  ṽ + ṽ ⊥ + B 2 ṽ ⊥ , (C6) Again, we see that the electromagnetic parts of the conserved quantities Q  and U only depend on the perpendicular parts of the velocity.

APPENDIX D: GRMHD PARALLEL MOMENTUM EQUATION
In this section we derive the general equation for the conservation of stress-energy parallel to the magnetic field in GRMHD, Equation 38.This derivation is based on a similar derivation in the Appendix of Camenzind (1986).First, we contract the fluid velocity   with the homogeneous Maxwell equation for GRMHD (Equation 20): where we used the fact that     = 0 and     = −1.Therefore, for degenerate magnetically dominated fields where the acceleration along the velocity direction is   ≡   ∇    .Next, we contract   with the conservation of stress-energy equation (Equation 18).We find first that contracting   with the divergence of   EM gives zero in all cases: In the second line we again used the fact that     = 0 and in the third line we used the relation for the divergence of   from Maxwell's equations, Equation D2.
Next we contract   with the divergence of  which we can easily verify from the definition of the entropy, Equation 37.If we assume the evolution of the gas is adiabatic and the entropy-per-particle is constant along field lines then   ∇   = 0. We can then use the relation Equation D6 to replace the ∇   term in Equation 38, and we find that which is Equation 39.
To take the cold limit of Equation D5, we simply assume |  ∇  | << |ℎ ∇    , |.In this case, which is the form of the equation for the parallel velocity we use in the force-free region for the MAD simulations reported in this paper.
APPENDIX E: SIMULATIONS WITH DIFFERENT  TRANS .
In this Appendix, we show results from an additional two MAD simulations run with our new hybrid method, varying the transition magnetization  trans .The simulation parameters and initial conditions are the same as the hybrid simulation in section 5, including the initial random pressure perturbations.In addition to the fiducial hybrid simulation with  trans = 50 presented in the main text, we ran additional simulations with  trans = 25 and  trans = 100.Our simulation set-up is stable at all three values of  trans .In Figure D1, we show the accretion rate  and MAD parameter  from  = 7500  g to  = 10000  g for the three simulations.The median  and  for the  trans = 25 and  trans = 100 simulations both differ from the values from the fiducial  trans = 50 simulation by less than 10%.
In Figure D2, we show radial slices of the averaged magnetization and plasma- in the three hybrid simulations.All three hybrid simulations achieve very high  (low ) in the jet region close to the black hole.However, the highest value of  (lowest value of ) achieved in the jet depends on the choice of  trans , with the  trans = 100 simulation reaching the highest overall  (lowest overall ) in the jet region.All three simulations have closely matched  and  profiles with each other and with the fiducial MHD simulation outside the  trans threshold.This paper has been typeset from a T E X/L A T E X file prepared by the author.

𝜇⊥
can be obtained from Equation 23.In fact, the McKinney (2006) solution Equation B1 is equal to what is obtained by setting  = 0 in the usual Noble et al. (2006) GRMHD inversion procedure (Equation

Figure 1 .
Figure 1.Plot of the mixing function  ( ) for the GRMHD simulations in this paper.We center the transition at  trans = 50, and use width  = 0.1 in Equation 45.Below  ( ) = 1/64 and above  ( ) = 63/64, we switch to entirely GRMHD or GRFFE evolution, respectively.

Figure 2 .
Figure 2. Special relativistic FFE tests fromKomissarov (2002).For each problem, dashed lines represent the initial solution, solid lines indicate the analytic solution at the given time, and symbols indicate the numerical solution for the reported magnetic or electric field component.In the breakdown test (bottom row), there is no analytic result, but our solution closely matches the result fromKomissarov (2002) indicating that that at time  = 0.02 the electric field strength is nearly equal to the magnetic field strength.

Figure 3 .
Figure 3. Linear Alfvén pulse test.In the top row, we show the − components of the magnetic and electric field in blue and red, respectively; in the bottom row, we show the − components in the same colors.Dashed lines indicate the initial condition; the solid lines indicate the analytic expectation for the solution at  = 1.5 and the open markers show the numerical solution.The numerical solution transitions from standard MHD inversion to the left of the dashed vertical line at  = 0.5 to FFE inversion to the right of  = 0.5, with the transition happening over a finite width using the mixing prescription in Equation 44, replacing  ( ) with  (  ) from Equation48.

Figure 4 .
Figure 4. Schwarzschild Bondi Accretion Test Problem.From top to bottom, we show profiles of the plasma density , Boyer-Lindquist radial velocity   , dimensionless temperature Θ gas , and magnetization  for the spherically symmetric Bondi accretion setup described in this section.The solid black line indicates the initial condition; the Bondi problem should remain stationary in time.Blue circles show the numerical solution after 100  g for a standard GRMHD method, grey triangles show the results from a GRMHD simulation which evolves the fluid energy adiabatically, and orange xs show the numerical solution at the same time from our hybrid approach.In the hybrid approach, we transition between GRMHD and GRFFE solutions at  trans = 50, indicated by the horizontal dashed line in the bottom panel; this transition magnetization corresponds to a stable transition radius of  ≈ 6.3  g , indicated by the vertical dashed line in all panels.

Figure 5 .
Figure 5. Convergence test for the Bondi problem.From left to right, we plot the  1 error (Equation52) between the solution at  = 100  g and the analytic solution for the density , internal energy  gas , and Boyer-Lindquist radial velocity   as a function of the number of radial cells   .We plot the  1 error for the standard MHD run in blue circles and the error for the hybrid GRMHD+GRFFE run in orange triangles.Results from the adiabatic GRMHD simulation are shown in grey.

Figure 6 .Figure 7 .
Figure 6.2D BZ Monopole Test.From top to bottom, we show poloidal profiles of the magnetization , Lorentz factor  in KS coordinates, and fieldline angular speed Ω  /Ω  for the monopole simulations presented in this section with dimensionless black hole spin  * = 0.5.We show results from the end of the simulation at  = 50 g .In the left column, we plot results from a purely force-free simulation, and in the right column we show results from a hybrid GRMHD+GRFFE simulation with the same initial conditions.The magenta contour in all plots shows the  = 100 surface, which is the center of the transition function  trans = 100 in the hybrid simulation.

Figure 8 .
Figure 8. Evolution of the mass accretion rate  (top) and dimensionless magnetic flux through the horizon  (bottom; Equation60) in the GRMHD and Hybrid MAD simulations from  = 0 to  = 10 4  g .

Figure 9 .Figure 10 .Figure 11 .Figure 12 .
Figure 9. Snapshot poloidal slices from the GRMHD (left column) and Hybrid GRMHD+GRFFE (right column) MAD simulations.From top to bottom, we plot the plasma mass density  in code units; the magnetization ; and the ratio of thermal to magnetic pressure .Both simulation snapshots were taken at  = 8550  g .In all plots we present snapshot data from a slice of constant azimuthal angle  in Kerr-Schild coordinates.Contours of constant   , the azimuthal component of the vector potential, are shown in white.The black disc indicates the black hole horizon, and the magenta contour in all plots indicates the  = 1 surface.

Figure 13 .
Figure 13.230 GHz images from both simulations with parameters (black hole mass, total flux density) scaled to match EHT observations of M87* (Event Horizon Telescope Collaboration et al. 2019a).The top column shows images from the standard GRMHD simulation taken at  = 8880  g ; the left column shows the image raytraced zeroing out emissivities wherever  >  cut = 25, while in the right column no regions have their emissivities set to zero in the radiative transfer.The bottom row shows corresponding images from a snapshot of the Hybrid GRMHD+GRFFE simulation at  = 8550  g .In all images, the cyan contour indicates the critical curve or "black hole shadow"(Bardeen 1973); the magenta curve indicates the image of the equatorial event horizon, or "inner shadow" feature(Chael et al. 2021).All images are displayed in a log color scale.
∇    fluid =   ∇  ℎ    +   ∇   = ℎ    +   ∇  , (D4)where we again used the fact that   is orthogonal to   , and the definition of the acceleration   .So, finally, from the conservation of total stress-energy and the relation Equation D2:0 =   ∇    EM +   fluid =   ∇   + ℎ ∇    ,(D5)which gives the GRMHD relation for the field-parallel velocity used in the main text, Equation 38.To obtain the adiabatic limit of Equation D5, we use the thermodynamic relation:

Table 1 .
,max / g   ×   ×    min / g ,  max / g Summary of 3D MAD simulations presented in this work.The mass accretion rate through the horizon  and the normalized magnetic flux  were averaged over the range [5000, 10000] / g .
Chael et al. 2019;simulation lightcurves over the range 7500  g <  < 10000  g with black hole mass mass and accretion rate scaled to match 2017 EHT observations of M87*.The blue lines represent lightcurves from images raytraced from the standard GRMHD simulation; the orange lines indicate data from the Hybrid GRMHD+GRFFE simulation.Dotted lines indicate that the lightcurves from images generated with  cut = 1; solid lines indicate the lightcurves from images generated with  cut = 25; dashed lines indicate lightcurves from images generated with no  cut applied.beimportant in determining observable image features and spectra from synchrotron emission (e.g.Chael et al. 2019; Event Horizon  Telescope Collaboration et al. 2019c).It is still not entirely clear if the rapid increase in temperature with  in this region is entirely physical, or if the rapid decrease in plasma density causes issues in the numerical treatment of dissipation