Bridging Scales in Black Hole Accretion and Feedback: Magnetized Bondi Accretion in 3D GRMHD

Fueling and feedback couple supermassive black holes (SMBHs) to their host galaxies across many orders of magnitude in spatial and temporal scales, making this problem notoriously challenging to simulate. We use a multi-zone computational method based on the general relativistic magnetohydrodynamic (GRMHD) code KHARMA that allows us to span 7 orders of magnitude in spatial scale, to simulate accretion onto a non-spinning SMBH from an external medium with a Bondi radius of R B ≈ 2 × 105 GM •/c 2, where M • is the SMBH mass. For the classic idealized Bondi problem, spherical gas accretion without magnetic fields, our simulation results agree very well with the general relativistic analytic solution. Meanwhile, when the accreting gas is magnetized, the SMBH magnetosphere becomes saturated with a strong magnetic field. The density profile varies as ∼r −1 rather than r −3/2 and the accretion rate Ṁ is consequently suppressed by over 2 orders of magnitude below the Bondi rate ṀB . We find continuous energy feedback from the accretion flow to the external medium at a level of ∼10−2Ṁc2∼5×10−5ṀBc2 . Energy transport across these widely disparate scales occurs via turbulent convection triggered by magnetic field reconnection near the SMBH. Thus, strong magnetic fields that accumulate on horizon scales transform the flow dynamics far from the SMBH and naturally explain observed extremely low accretion rates compared to the Bondi rate, as well as at least part of the energy feedback.


INTRODUCTION
Most nearby galaxies are found to harbor central supermassive black holes (SMBHs) whose masses are correlated with properties of the stellar components of their hosts (e.g., Magorrian et al. 1998;Ferrarese & Merritt 2000;Gebhardt et al. 2000;Kormendy & Ho 2013).However, the details of how gas flows into the galactic nucleus from large cosmic distances and how the SMBH in turn imparts feedback into the galaxy remain coupled unresolved problems.
Due to limitations in computational power and resolution, typically there has been a partitioning of scales in tackling the SMBH feeding and the feedback problem.In galaxy-scale simulations, the gas flows from kpc scales and the co-evolution of SMBH and galaxy are explicitly modeled (e.g., Sijacki et al. 2015;Rosas-Guevara et al. 2016;Weinberger et al. 2018;Ricarte et al. 2019;Ni et al. 2022;Wellons et al. 2023).However, even in idealized simulations of an isolated galactic nucleus, prescriptions need to be adopted as "sub-grid models" for including accretion and feedback (e.g., Li & Bryan 2014;Fiacconi et al. 2018;Anglés-Alcázar et al. 2021;Talbot et al. 2021;Weinberger et al. 2023) usually at a scale ≲ 10 pc, which translates to ∼ 3 × 10 4 r g for the ∼ 6.5 × 10 9 M ⊙ SMBH in M87; here r g = GM • /c 2 is the gravitational radius, where M • is the SMBH mass.Meanwhile, probing smaller scales, general relativistic magnetohydrodynamic (GRMHD) simulations focus on the inner few tens or hundreds of r g , self-consistently tracing gas flows and feedback on these limited spatial scales (e.g., Gammie et al. 2003;Tchekhovskoy et al. 2011;Porth et al. 2019;Chatterjee et al. 2023).Typically, GRMHD simulations use idealized settings and disregard the larger-scale cosmological environment.Bridging these vastly different scales while self-consistently following gas flows and incorporating both the GRMHD effects and cosmological evolution remains challenging.
Recently, some attempts to connect these disparate scales have been made using nested simulations where each smaller-scale simulation is initialized from the next larger-scale simulation (e.g., Hopkins & Quataert 2010;Ressler et al. 2020;Guo et al. 2023); or by adopting Lagrangian hyper-refinement methods (e.g., Anglés-Alcázar et al. 2021;Hopkins et al. 2023); or by pushing out the simulation regime to larger scales (e.g., Lalakos et al. 2022;Kaaz et al. 2023).While successful in studying the process of accretion, feedback from the SMBH has not been followed out to galaxy scales in prior work because either the communication between scales is only directed inwards or because of the inability to include the entire galaxy and the black hole.
We employ here a multi-zone computational method which represents a first attempt to not only follow accretion flows but also trace feedback across 7 orders of magnitude from the event horizon to galactic scales.With this method, we study purely hydrodynamic Bondi accretion as well as the magnetized Bondi problem.
The outline of this Letter is as follows: In Section 2 we give an overview of the numerical scheme that captures a wide dynamic range simultaneously.In Section 3, we present a purely hydrodynamic simulation and test our numerical scheme by comparing against the general relativistic (GR) analytical solution.In Section 4 we include magnetic fields for a more realistic representation of the environment of SMBHs and study how Bondi accretion is modified.We discuss our findings of feedback via convection in Section 5 and summarize conclusions of our study in Section 6.Additional details are presented in a set of short appendices: various GRMHD quantities are defined in Appendix A, the details of the simulation set-up are outlined in Appendix B, and a resolution and initial condition study is presented in Appendix C.
Throughout, we adopt units where GM • = c = 1, so our unit of length r g = 1 and time t g ≡ r g /c = 1.Gas temperature T is a dimensionless quantity which is normalized by µc2 /k B , where µ is the mean molecular weight and k B is the Boltzmann constant; thus gas pressure p g = ρT with ρ the density.We denote the time average of a quantity X as X and the density-weighted, time-and shell-averaged value of X as ⟨X⟩ (Appendix A).Finally, the Bondi radius is is the sound speed in the external medium far outside R B , T ∞ is the temperature there, and γ ad is the gas adiabatic index.The freefall time at the Bondi radius is t B ≡ (R B /r g ) 3/2 t g .For the models presented here, γ ad = 5/3, T ∞ = 3.4 × 10 −6 (≈ 2 × 10 7 K), R B = 1.8 × 10 5 r g , and we consider a non-spinning SMBH (Schwarzschild metric).

NUMERICAL METHODS
Our numerical scheme utilizes the GRMHD code KHARMA1 , a performance-portable C++ implementation based on iharm3D (Prather et al. 2021); iharm3D is itself an extension of HARM, an efficient secondorder conservative finite-volume scheme for solving MHD equations on Eulerian meshes in stationary curved space-times (Gammie et al. 2003).KHARMA offers a more flexible, portable, and scalable implementation suitable for multiple uses, by leveraging the Parthenon Adaptive Mesh Refinement Framework and the Kokkos programming model (Grete et al. 2023;Trott et al. 2022).
The approach adopted thus far to handle a wide range of spatial and temporal scales has been to run a series of nested simulations sequentially, starting from the outside and proceeding down to the BH (e.g., Ressler et al. 2020;Guo et al. 2023).This set-up is appropriate for studying the feeding of the BH from the external galactic medium (or stellar winds), but it cannot tackle feedback from the BH to the galaxy.To handle both processes, we run nested simulations from radii r ≫ R B down to r g as well as from r g out to r ≫ R B .In addition, we cycle in and out hundreds of times to ensure self-consistent, two-way communication between the BH and the external medium.
We use a grid based on spherical coordinates and split the simulation volume between r = r g and r ≫ R B into a set of n overlapping annuli evenly spaced in log(r), each of which we label as zone-i.For i ∈ {0, 1, ..., (n − 1)}, zone-i extends from an inner radius r i,in = 8 i r g to an outer radius r i,out = 8 i+2 r g .We start by initializing and running zone-(n − 1), the largest annulus, for some fraction of its free-fall time.We then simulate zone-(n − 2), inheriting the final state of zone-(n−1) as the initial condition over half the new domain.This continues inward to zone-0 and back out in a sort of "V-cycle."Our method can be thought of as defining a global domain, but "pausing" everything outside the particular active region, shifting the active region up and down the range of scales a large number of times (in our fiducial MHD run, the V-cycle was traversed 254 times) until the full system converges to a steady state.(f) Figure 1.Time-averaged radial profiles of selected quantities for our HD (red) and MHD (black) Bondi accretion simulations with RB = 1.8 × 10 5 rg (gray vertical dashed line).In this and forthcoming plots, negative values are plotted with dotted rather than solid lines.While most naturally expressed in gravitational units (lower x-axis), we also present distances scaled to M87 (upper x-axis).Where available, analytic solutions of relativistic Bondi accretion are plotted with thick gray lines.Magnetic fields lead to a shallower density slope, a dramatically suppressed accretion rate, and positive energy feedback.
We first test the accuracy of our method by numerically simulating hydrodynamic (HD) spherical Bondi (1952) accretion on a Schwarzschild BH; a relativistic analytical solution is available for this problem (Michel 1972;Shapiro & Teukolsky 1983).We use a total of n = 7 annuli with an outer radius of 8 8 ∼ 2 × 10 7 r g and a total runtime of t run (R B ) ∼ 9 t B (∼ 0.7 Myr for M87) at the Bondi radius.We confirmed convergence by checking that the results at t run = 9 t B remain unchanged even when we run up to t run > 56 t B .The resolution is 128 3 per annulus giving a total simulation resolution of 512 × 128 × 128.Details of the set-up are presented in Appendix B.1.
Figures 1(a,b) show the radial profiles of t, θ, φaveraged gas density ρ and temperature T (red lines).The simulation results show excellent agreement with the GR Bondi analytic solution (thick gray line) over 7 orders of magnitude in radius.Further, the timeaveraged accretion rate Ṁ (r) in Figure 1(c) matches the analytical Bondi accretion rate ṀB very well, Ṁ / ṀB = 1, across all radii within R B , confirming a steady state.We note that Guo et al. (2023) have previously sim-ulated Bondi accretion under Newtonian gravity spanning 5 orders of magnitude.Here we have successfully reproduced the GR version over a larger range of scales.
Analogous to Ṁ (r) we can define a time-averaged energy inflow rate towards the BH, Ė(r) (Equation (A5)), but this is not useful for our purposes since it includes the rest mass energy of the accreted gas.Instead we consider Ėnet (r) ≡ ( Ṁ − Ė) (Equation (A6)), which removes the rest mass contribution and is flipped in sign such that it is the net rate of outflow of energy (i.e., feedback) from the BH to large radii.Normalizing by Ṁ 10 ≡ Ṁ (r = 10 r g ), we then obtain a feedback efficiency η tot (Equation (A7)) (1) Following Appendix A, and noting that there is no magnetic field and the system is spherically symmetric, (2) The agreement of the simulation results in Figure 1(d) with this theoretical expectation is excellent.

MAGNETIZED BONDI ACCRETION
Accretion flows in galactic nuclei involve magnetized plasma, and it has been long known that magnetic fields strongly perturb the energetics (Shvartsman 1971;Meszaros 1975) and dynamics (Igumenshchev & Narayan 2002;Pen et al. 2003;Lalakos et al. 2022) of Bondi accretion.Apart from a toy model with purely radial magnetic field and a non-spinning BH, where no dynamical interaction occurs between the fluid and the magnetic field (Gammie et al. 2003), there is no known analytical solution to the magnetized accretion problem.It can be studied only via numerical simulations.Our computational technique is well-suited for this.
To simulate the magnetohydrodynamic (MHD) problem, we use the same external density and temperature as adopted for the HD case (Section 3), but now we include an initial magnetic field in the z-direction with plasma-β ≡ p g /p b ≈ 1 where p b is the magnetic pressure.Our choice of β is motivated by the Milky Way ISM (e.g., Ferrière 2020;Guerra et al. 2023), but we have verified that the results remain unchanged if we set the initial β ≈ 10 (Appendix C).We use n = 8 zones, which allows us to probe out to r ∼ 10 8 r g (and reach steady state out to nearly 10 7 r g ).Our effective resolution is 576 × 128 × 128 over the whole domain.Additional setup details are given in Appendix B.2.
Unlike the laminar flow that we find in the HD simulation, the MHD simulation gives a dynamic and turbulent flow which is driven by magnetic stresses and reconnection.We need to time-average the results to obtain meaningful steady state profiles.The black lines in Figure 1 show time-averaged radial profiles of a number of quantities after the simulation was run for a physical duration of t run (R B ) ∼ 9t B at the Bondi radius.The temperature profile is fairly similar to the HD Bondi case (compare black and red lines), but the density profile is significantly modified due to the presence of the magnetic field.The radial scaling is now ρ ∝ r −1 , which is similar to the slope reported in several previous GRMHD simulations of rotating accretion flows (e.g., Ressler et al. 2020;Chatterjee & Narayan 2022) and even some hydrodynamical simulations (Guo et al. 2023).It is also consistent with observational constraints in M87* and Sgr A* from X-ray and EHT observations (e.g.Chatterjee & Narayan 2022).There is a curious dip in the density and a corresponding bump in the temperature near the Bondi radius, which we explain in Section 5.
Despite running the simulation for an effective time of ∼ 9t B , we were unable to obtain a reliable measurement of Ṁ at r ≳ 10 2 r g as seen in Figure 1(c), unlike the HD case.This is because the medium exhibits violent convective turbulence (Section 5), which produces large random velocities that overwhelm the mean accretion velocity.To illustrate the difficulty, we show the average mass inflow rate Ṁin (solid blue) and mass outflow rate Ṁout (dotted cyan), defined in Appendix A, as a function of radius.With increasing distance from the BH, both rates become much larger than the net Ṁ , which is the difference of the two very similar numbers.It is therefore challenging to calculate Ṁ .Here, we restrict ourselves to smaller radii, where we believe our Ṁ estimates are reliable.At r = 10 r g3 we find Ṁ 10 ≈ 0.005 ṀB , i.e., the accretion rate is suppressed by more than a factor of 100 relative to the HD Bondi case.This estimate is consistent with the observed accretion rate of ∼ 10 −3 − 10  We plot the time-averaged energy efficiency η tot (r) (Equation 1) from the simulation in Figure 1(d).There are two striking features.First, the radial profile of η tot is nearly constant (ideally it should be perfectly constant), which indicates that the simulation energetics have converged well.4Second, η tot is positive, i.e., energy flows out from the accretion flow to the external medium.In other words, the simulation exhibits bona fide feedback.It is worth emphasizing that we are able to study feedback only because our numerical technique is specifically designed to allow two-way communication between the BH and the external medium (Section 2).The amount of feedback energy is ∼ 10 −2 Ṁ 10 c 2 ∼ 5 × 10 −5 ṀB c 2 , which is lower than the values adopted in cosmological simulations ≳ 10 −3 ṀB c 2 (e.g., Booth & Schaye 2009;Weinberger et al. 2018;Davé et al. 2019;Ni et al. 2022).However, the feedback in our runs operates without BH spin and without net rotation in the accreting gas 5 .We speculate therefore, that it represents the minimum feedback that any hot quasi-spherical magnetized accretion flow will produce.The actual feedback will likely be larger when more favorable conditions are present (e.g., a spinning BH).
The radial profiles of two magnetic field related quantities are shown in the bottom panels of Figure 1: the plasma β, and the normalized enclosed magnetic flux ϕ b (Tchekhovskoy et al. 2011), (3) We average the inverse of β for stability, which is a commonly used technique.The profile we obtain has a flat plateau with a nearly constant value β ∼ 3 between r = 10 2 − 10 5 .Over this extended volume, the field is in equi-partition with the gas and the magnetic flux ϕ b varies ∝ r, indicating that the organized poloidal field varies as B pol ∝ r −1 .The pressure in the organized field is also in equi-partition: B 2 pol ∝ ρT ∝ r −2 .We find that for radii r < 10 2 , β declines rapidly to well below unity and ϕ b saturates at ∼ 30.Our results are very reminiscent of the magnetically arrested disk (MAD) model of MHD accretion (Tchekhovskoy et al. 2011;Igumenshchev et al. 2003;Narayan et al. 2003, see also Bisnovatyi-Kogan & Ruzmaikin 1974).In previous work, Ressler et al. (2020) found that the entire stellar wind-fed accretion flow in their model of Sgr A* approached the MAD state.In our case, MAD accretion extends from the Bondi radius R B ≈ 1.8 × 10 5 r g down to the BH horizon at 2 r g .One minor difference is that ϕ b ∼ 30 at the BH horizon in our simulation, whereas typical GRMHD simulations which consider orbiting gas (instead of spherical infall) find ϕ b ∼ 50 for a non-spinning BH (Tchekhovskoy et al. 2012;Narayan et al. 2022).

FEEDBACK VIA RECONNECTION-DRIVEN CONVECTION
As noted, our simulation permits us for the first time to study the coupled BH feeding and feedback problem.Following the time evolution of the BH mass accretion rate Ṁ (in units of the Bondi rate ṀB ) and the magnetic flux parameter ϕ b at the horizon (top panel of Figure 2) reveals large amplitude fluctuations in both quantities around a fairly stable mean value, Ṁ 10 ∼ 0.005 ṀB and ϕ b ∼ 30.The time scale of the largest fluctuations is a few Bondi times, as expected for feeding from the Bondi radius.The cause of the fluctuations is violent turbulence, evidence for which can be seen in the multi-scale snapshots in the middle and lower rows of Figure 2.
To investigate the mechanism driving the turbulence, we explore the physics of the outward flow of energy.The black curves in the upper panels of Figure 3 show the radial profile of the efficiency parameter η tot .The other curves show individual physical components: η fl (blue, left top panel), the energy carried in the fluid (as opposed to the magnetic field), η adv (green, middle), the advected part of the fluid energy, and η conv (red, right), the part of the fluid energy transport that can be attributed to convection.The definitions of these components are given in Appendix A and are identical to those used by Ressler et al. (2021).
At all radii r ≲ R B , we find that η fl very closely matches η tot , i.e., essentially all the energy flux is carried by the fluid, reminiscent of the convection-dominated accretion model of Narayan et al. 2000;Quataert & Gruzinov 2000), while the electromagnetic contribution is neg-ligible.The electromagnetic energy transport would likely become important if the BH had spin, since field lines near the horizon would be frame dragged and produce a relativistic jet.In a smaller scale simulation with R B = GM • /c 2 s,∞ = 100 and a spinning BH (a * = 0.9375), Ressler et al. (2021) found that the electromagnetic η exceeded η fl .For the non-spinning case here, we find the opposite.
When we divide η fl into an advective part η adv and a convective part η conv , we find that convection dominates by a substantial margin between r = r g and ∼ 0.1R B ; energy transport is principally via convection.From the work of Igumenshchev & Narayan (2002), the convection is initiated by magnetic reconnection near the BH, which transforms the dissipated magnetic energy into thermal energy (entropy) and kinetic energy (Ripperda et al. 2022).Convective turbulence develops and extends all the way to r ∼ 0.1R B , carrying the dissipated energy outward.Note that our definition of η conv closely follows Ressler et al. (2021).However, it is possible that some fraction of what we call convection corresponds to a violently varying wind (e.g., Yuan et al. 2015).Interestingly, although magnetic reconnection drives the convection, the field does not transport the energy because field lines remain largely fixed in radius in the MAD state.
The bottom three panels in Figure 3 present detailed information on the fluid, advective and convective energy fluxes, F fl E , F adv E , F conv E (defined in Appendix A), in the logarithmic radius -polar angle (log 10 r -θ) plane.The advective flux carries energy toward the BH near the midplane but away from the BH near the poles.When integrated over θ, the two zones nearly cancel.On the other hand, convection transfers energy outward at all θ and overwhelms the advective energy at all r ≲ 0.1R B so that the net fluid flux F fl E is outward everywhere.While Ressler et al. (2021) concluded that convection was inefficient in their spinning BH set-up, we find exactly the opposite for our non-spinning BH.Pen et al. (2003) found energy inflow rather than outflow in their non-relativistic simulation, which they termed "frustrated convection."The difference from our result might be because their magnetic field did not reach the saturated MAD limit and/or their spatial dynamic range was only a few tens.
Finally, we see in Figure 3 that the convective energy transport ceases at R B , as expected since convection requires gravity and buoyancy, while gravity becomes subdominant for r > R B .Thus, in our model, the feedback energy carried outward by convection is dumped near the Bondi radius and has nowhere to go.This explains the bump in the temperature profile of the MHD run near the Bondi radius in Figure 1, and the corresponding dip in density (required by pressure balance).

SUMMARY AND CONCLUSIONS
We have simulated general relativistic accretion on a non-spinning SMBH using a numerical technique which allows us to maintain two-way communication across the entire range of spatial scales from the BH horizon to beyond the Bondi radius, R B ≈ 1.8 × 10 5 r g .For spherical accretion without magnetic fields, we find perfect agreement with the GR analytical Bondi solution, i.e., for r < R B , the density and temperature of the simulated gas scale as ρ ∝ r −3/2 , T ∝ r −1 respectively, as expected, and the measured accretion rate also matches the analytical Bondi accretion rate ṀB .The thermal energy of the external gas, which is 1.5/R B ≈ 10 −5 (Equation 2) of rest-mass energy, is advected into the BH and there is no feedback in the opposite direction.
Adding a magnetic field of realistic strength to the external gas changes the results drastically.The density slope inside R B changes to ρ ∝ r −1 , and the mass accretion rate Ṁ into the BH drops dramatically to ≈ 5×10 −3 ṀB , consistent with observations of M87* and Sgr A*.Equi-partition of magnetic and gas pressure is achieved (β ∼ 3) over the radius range 10 2 r g ≲ r ≲ R B , and the magnetic pressure becomes much more dominant close to the BH.Except for the lack of net rotation, the solution closely resembles the magnetically arrested disk (MAD) model.Significantly, the magnetized simulation shows net energy feedback from the accretion flow to the external medium at a rate ≈ 10 −2 Ṁ c 2 , i.e., with an efficiency η tot ∼ 0.01.Expressed in terms of the Bondi accretion rate, the energy outflow rate is ∼ 5 × 10 −5 ṀB c 2 .This positive feedback occurs even in the absence of BH spin or net gas rotation.The energy outflow is not strongly collimated; it is quasi-isotropic.Advection carries energy inwards in the mid-plane and transports energy outwards near the poles, while convection transfers energy outwards at all θ.We find that the level of energy feedback is significantly lower than that typically invoked in cosmological simulations, so we view the feedback estimate in the present work as representing a lower limit to energy outflow.Additional ingredients, especially BH spin and gas angular momentum, are needed to study the problem more fully.The dominant outward energy transport mechanism is convection driven by magnetic reconnection in the magnetosphere near the BH.Even though a strong magnetic field is essential for this mechanism to work, the energy is carried almost entirely in the form of thermal and kinetic energy of the gas, while the electromagnetic energy flux is negligible.Thus the feedback mechanism seen here is different from the electromagnetically-dominated Cho et al.
outflows found in relativistic jets and rotation-driven disk winds.
These new insights are possible because our numerical technique permits two-way communication between the BH and the external medium over many decades of spatial and temporal scales.In future work, we plan to extend our investigation to the case of a spinning BH that also includes the angular momentum of the gas.and the corresponding efficiency η tot as: To study the contribution of the fluid alone (ignoring the magnetic field), the relevant terms in the stress-energy tensor are: Then the time-averaged net fluid energy outflow rate is: where the relativistic Bernoulli parameter is (e.g., Penna et al. 2013, but not including the magnetic field) In the same manner as in (A7), we define the corresponding efficiencies, η fl (r), η adv (r), η conv (r).
In addition to the time-averaged and shell-integrated radial profiles Ė(r) defined above, we also define local timeand azimuth-averaged energy fluxes as a function of r and θ: F adv E (r, θ) ≡ (1/2π) Be ρu r √ −g dφ, (A14) In the hydrodynamic simulation described in Section 3, we use a domain with 7 zones (i = 0 − 6) covering the range [1, 1.7 × 10 7 ] r g .Zone-6 is initialized with a static fluid of constant density and temperature, and no perturbations, so as to recover the analytic Bondi solution without rotation.We use spherical Kerr-Schild coordinates with coordinate four-vector (x t , x r , x θ , x φ ), where the radial coordinate is x r ≡ log r, sometimes referred to as "eKS" coordinates.Rather than allow vacuum to arise in the inner regions before the initial material from zone-6 accretes, new areas are filled with material before they are first simulated, with the same density and temperature as the innermost nonzero density zone present when they are initialized.Initialization in this way prevents disruption of the simulation, but bears no resemblance to the Bondi solution, validating that the scheme can converge to the correct result regardless of its initialization.After the first pass, no further initialization is needed.During the rest of the simulation, each active annulus receives half its initial data from the previous active annulus and half from the last output of the current active annulus.In the process, the inner and outer radial boundaries of the active annulus end up with new values of ρ, u and u µ , which are treated as fixed boundary conditions while the annulus is evolved.Once the system has reached its final steady state, conservation of mass and energy between annuli are satisfied in a time-averaged sense.

B.2. Magnetohydrodynamic Simulation Setup
The magnetized simulation in Section 4 is run over a larger domain with 8 zones [1, 1.3 × 10 8 ] r g .Zone-7 is initialized with a similar density and temperature as in the hydrodynamical model, but the two differ in several ways.The density is now initialized in the same way in all zones, to follow ρ init (r) = ρ 0 (r + R B )/r such that it is ∝ r −1 interior to R B and constant outside of R B .This is motivated by simulations (e.g.Ressler et al. 2020;Chatterjee & Narayan 2022;Guo et al. 2023) and observations (Chatterjee & Narayan 2022) that find a ρ ∝ r −1 relation in the presence of strong magnetic fields.However, the results are insensitive to the choice of ρ init , as discussed in Appendix C. The temperature is initialized to the Bondi analytical solution, and the velocities are initialized to zero.The magnetic fields are initialized with a vector potential A φ (r, θ) = b z r(r + R B ) sin 2 θ/2 such that the plasma-β parameter is approximately constant over all radii.The resulting magnetic field is vertical outside R B , and becomes more radial inside R B .The field is normalized with the parameter b z such that the initial β init (r) ∼ 1, motivated by observations of the Milky Way ISM (Ferrière 2020;Guerra et al. 2023).A per-zone random perturbation of maximum 5% is applied to the internal energy when initializing fluid at all radii, to help break axisymmetry.Since the magnetic field picks out a rotation axis for the resulting flow parallel to the coordinate axis, we choose a coordinate system which compresses the θ coordinate away from coordinate axes and toward the midplane.We adopt a simplified version of "funky modified Kerr-Schild" (FMKS) coordinates (see e.g.θ j in the Appendix in Wong et al. 2021), eliminating the cylindrification term to obtain: where N ≡ π/2(1 + χ −α t /(1 + α)) −1 is a normalization factor, y = 2x θ − 1, χ t = 0.8, α = 16, and x θ ∈ [0, 1].The simplified system retains the primary benefit of FMKS coordinates: zones near the coordinate pole are widened, relaxing the Courant condition and allowing a greater simulation timestep, by a factor of 2 or more, dramatically reducing simulation cost.As we are interested primarily in steady state behavior rather than resolving particular structures, e.g., a stable jet, the cylindrification is unnecessary for our case.
In addition to changing the coordinates, we lower the Courant factor for safety.We also introduce dynamical floors: in order to maintain ρ > 10 −6 r −3/2 , internal energy u > 10 −8 r −5/2 , temperature u/ρ < 100, and magnetization b 2 /ρ < 100, we introduce material and energy in the coordinate frame.We also reduce the fluid velocity when the Lorentz factor measured by the Eulerian observer is larger than γ max = 10.While running each active annulus, we keep ρ, u and u µ fixed at the inner and outer boundaries, as in the hydro case.In addition, we also keep b µ fixed and make sure that we preserve ∇ • B = 0.
We have validated the 'multi-zone' method described in Section 2 in the presence of strong magnetic fields, by comparing a smaller 4-zone simulation with R B ≈ 460 against a single-zone simulation of the same problem.We obtained substantially similar averaged radial profiles.These comparisons, along with detailed descriptions of the method, boundary conditions, etc., will be presented in a forthcoming longer paper.

C. RESOLUTION AND INITIAL CONDITION STUDY
We have presented our fiducial MHD simulation results in Section 4 for a resolution of 128 3 per annulus with the FMKS coordinate system and initialized with ρ init ∝ r −1 at r < R B .Here we show that the results remain the same even when we choose different resolutions, coordinate systems, or initial conditions.First, we introduce a new "Modified Kerr-Schild" (MKS) coordinate system (see e.g.θ g in the Appendix of Wong et al. 2021) where the resolution is more focused in the midplane but the zones near the pole are not widened: Here we use h = 0.3.In the case of h = 1, it is equivalent to the uniform θ grid of the eKS coordinates in Section B.1.Six extra simulations were run with parameters outlined in Table 1.Using MKS coordinates and other set-up details the same as in our fiducial run, we ran simulations at (i) 48 3 resolution per annulus and (ii) 64 3 resolution per annulus.Run (iii) with 96 3 resolution per annulus used FMKS but with a different α = 14 which is less extreme than α = 16 in the fiducial run.Run (iv) used a weaker initial magnetic field with β ∼ 10 and used MKS coordinates with 64 3 resolution per annulus.Runs (v) and (vi) were initialized with a different initial density profile compared to the other runs and used MKS with 64 3 resolution.Run (v) was initialized with a piecewise constant density profile as in the HD run.Run (vi) was initialized with the analytical HD Bondi solution where the density scales as ρ init ∝ r −3/2 for r < R B .
Radial profiles of all the quantities of interest for the six extra simulations are compared with the fiducial run in Figure 4. Overall, the profiles are very similar between the various runs implying that our conclusions from the fiducial run hold true regardless of the choice of coordinate system, resolution, initial magnetic field strength, or initial density profile.The accretion rate at 10 r g of all five runs lies in the range of Ṁ ≈ (1 − 5) × 10 −3 ṀB , which is consistent with 2 −2 ṀB in M87* (when comparing the BH accretion rate Ṁ ∼ 10 −3 M ⊙ /yr in Event Horizon Telescope Collaboration et al. (2021) with the Bondi accretion rate ṀB ∼ 0.2 M ⊙ /yr in Russell et al. (2015)).It is also consistent with the observed accretion rate for Sgr A* 10 −4 − 10 −2 ṀB (Wang et al. 2013).

Figure 2 .
Figure 2. Top: Accretion rate (black) and magnetic flux parameter (blue) as a function of time for our MHD simulation.The mean accretion rate Ṁ 10 ∼ 5 × 10 −3 ṀB and the mean magnetic flux at the horizon ϕ b ∼ 30 are plotted in gray and blue horizontal lines respectively.We average values over the second half of the simulations (pink background), when the simulation has reached steady state.Bottom: Slices of β (left) and ρ (right) in our simulation, spanning 8 orders of magnitude in spatial scale.Lines of constant magnetic flux representing φ-averaged magnetic fields are overlaid in black.

Figure 3 .
Figure 3. Top: Dissecting contributions to the total energy outflow (η tot ), including the fluid contribution (η fl ), the advective outflow efficiency (η adv ), and the convective outflow efficiency (η conv ).Since η fl ≈ η tot , the electromagnetic field does not contribute much directly to the energy outflow.Instead, the energy is transported primarily via convection until near the Bondi radius (gray vertical dashed line).Positive (negative) values are shown in solid (dotted) lines.Bottom: Azimuthally-averaged energy fluxes in the r, θ plane, normalized by the accretion rate at 10 rg.Outflowing energy is shown in red and inflowing energy in blue.

Table 1 .
Simulation set-up for different runs.