Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
License: CC BY 4.0
arXiv:2310.19135v2 [astro-ph.HE] 01 Dec 2023

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

Hyerin Cho (조혜린) Center for Astrophysics |||| Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA Ben S. Prather CCS-2, Los Alamos National Laboratory, PO Box 1663, Los Alamos, NM 87545, USA Ramesh Narayan Center for Astrophysics |||| Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA Priyamvada Natarajan Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA Department of Astronomy, Yale University, Kline Tower, 266 Whitney Avenue, New Haven, CT 06511, USA Department of Physics, Yale University, P.O. Box 208121, New Haven, CT 06520, USA Kung-Yi Su Center for Astrophysics |||| Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA Angelo Ricarte Center for Astrophysics |||| Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA Koushik Chatterjee Center for Astrophysics |||| Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA
Abstract

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 magneto-hydrodynamic (GRMHD) code KHARMA that allows us to span 7777 orders of magnitude in spatial scale, to simulate accretion onto a non-spinning SMBH from an external medium with Bondi radius RB2×105GM/c2subscript𝑅𝐵2superscript105𝐺subscript𝑀superscript𝑐2R_{B}\approx 2\times 10^{5}\,GM_{\bullet}/c^{2}italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT ≈ 2 × 10 start_POSTSUPERSCRIPT 5 end_POSTSUPERSCRIPT italic_G italic_M start_POSTSUBSCRIPT ∙ end_POSTSUBSCRIPT / italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, where Msubscript𝑀M_{\bullet}italic_M start_POSTSUBSCRIPT ∙ end_POSTSUBSCRIPT 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 r1similar-toabsentsuperscript𝑟1\sim r^{-1}∼ italic_r start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT rather than r3/2superscript𝑟32r^{-3/2}italic_r start_POSTSUPERSCRIPT - 3 / 2 end_POSTSUPERSCRIPT and the accretion rate M˙˙𝑀\dot{M}over˙ start_ARG italic_M end_ARG is consequently suppressed by over 2 orders of magnitude below the Bondi rate M˙Bsubscript˙𝑀𝐵\dot{M}_{B}over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT. We find continuous energy feedback from the accretion flow to the external medium at a level of 102M˙c25×105M˙Bc2similar-toabsentsuperscript102˙𝑀superscript𝑐2similar-to5superscript105subscript˙𝑀𝐵superscript𝑐2\sim 10^{-2}\dot{M}c^{2}\sim 5\times 10^{-5}\dot{M}_{B}c^{2}∼ 10 start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT over˙ start_ARG italic_M end_ARG italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ∼ 5 × 10 start_POSTSUPERSCRIPT - 5 end_POSTSUPERSCRIPT over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. 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.

Accretion (14), Active galactic nuclei (16), Bondi accretion (174), Schwarzschild black holes (1433), Supermassive black holes (1663), Magnetohydrodynamical simulations (1966)

1 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 10less-than-or-similar-toabsent10\lesssim 10≲ 10 pc, which translates to 3×104rgsimilar-toabsent3superscript104subscript𝑟𝑔\sim 3\times 10^{4}r_{g}∼ 3 × 10 start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT for the 6.5×109Msimilar-toabsent6.5superscript109subscript𝑀direct-product\sim 6.5\times 10^{9}\,M_{\odot}∼ 6.5 × 10 start_POSTSUPERSCRIPT 9 end_POSTSUPERSCRIPT italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT SMBH in M87; here rg=GM/c2subscript𝑟𝑔𝐺subscript𝑀superscript𝑐2r_{g}=GM_{\bullet}/c^{2}italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT = italic_G italic_M start_POSTSUBSCRIPT ∙ end_POSTSUBSCRIPT / italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT is the gravitational radius, where Msubscript𝑀M_{\bullet}italic_M start_POSTSUBSCRIPT ∙ end_POSTSUBSCRIPT is the SMBH mass. Meanwhile, probing smaller scales, general relativistic magnetohydrodynamic (GRMHD) simulations focus on the inner few tens or hundreds of rgsubscript𝑟𝑔r_{g}italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT, 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𝐺subscript𝑀𝑐1GM_{\bullet}=c=1italic_G italic_M start_POSTSUBSCRIPT ∙ end_POSTSUBSCRIPT = italic_c = 1, so our unit of length rg=1subscript𝑟𝑔1r_{g}=1italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT = 1 and time tgrg/c=1subscript𝑡𝑔subscript𝑟𝑔𝑐1t_{g}\equiv r_{g}/c=1italic_t start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT ≡ italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT / italic_c = 1. Gas temperature T𝑇Titalic_T is a dimensionless quantity which is normalized by μc2/kB𝜇superscript𝑐2subscript𝑘𝐵\mu c^{2}/k_{B}italic_μ italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT / italic_k start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT, where μ𝜇\muitalic_μ is the mean molecular weight and kBsubscript𝑘𝐵k_{B}italic_k start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT is the Boltzmann constant; thus gas pressure pg=ρTsubscript𝑝𝑔𝜌𝑇p_{g}=\rho Titalic_p start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT = italic_ρ italic_T with ρ𝜌\rhoitalic_ρ the density. We denote the time average of a quantity X𝑋Xitalic_X as X¯¯𝑋\overline{X}over¯ start_ARG italic_X end_ARG and the density-weighted, time- and shell-averaged value of X𝑋Xitalic_X as Xdelimited-⟨⟩𝑋\langle X\rangle⟨ italic_X ⟩ (Appendix A). Finally, the Bondi radius is RBGM/cs,2subscript𝑅𝐵𝐺subscript𝑀superscriptsubscript𝑐𝑠2R_{B}\equiv GM_{\bullet}/c_{s,\infty}^{2}italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT ≡ italic_G italic_M start_POSTSUBSCRIPT ∙ end_POSTSUBSCRIPT / italic_c start_POSTSUBSCRIPT italic_s , ∞ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, where cs,=γadTsubscript𝑐𝑠subscript𝛾adsubscript𝑇c_{s,\infty}=\sqrt{\gamma_{\rm ad}T_{\infty}}italic_c start_POSTSUBSCRIPT italic_s , ∞ end_POSTSUBSCRIPT = square-root start_ARG italic_γ start_POSTSUBSCRIPT roman_ad end_POSTSUBSCRIPT italic_T start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT end_ARG is the sound speed in the external medium far outside RBsubscript𝑅𝐵R_{B}italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT, Tsubscript𝑇T_{\infty}italic_T start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT is the temperature there, and γadsubscript𝛾ad\gamma_{\rm ad}italic_γ start_POSTSUBSCRIPT roman_ad end_POSTSUBSCRIPT is the gas adiabatic index. The free-fall time at the Bondi radius is tB(RB/rg)3/2tgsubscript𝑡𝐵superscriptsubscript𝑅𝐵subscript𝑟𝑔32subscript𝑡𝑔t_{B}\equiv(R_{B}/r_{g})^{3/2}t_{g}italic_t start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT ≡ ( italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT / italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT 3 / 2 end_POSTSUPERSCRIPT italic_t start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT. For the models presented here, γad=5/3subscript𝛾ad53\gamma_{\rm ad}=5/3italic_γ start_POSTSUBSCRIPT roman_ad end_POSTSUBSCRIPT = 5 / 3, T=3.4×106subscript𝑇3.4superscript106T_{\infty}=3.4\times 10^{-6}italic_T start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT = 3.4 × 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT (2×107Kabsent2superscript107K\approx 2\times 10^{7}\,{\rm K}≈ 2 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT roman_K), RB=1.8×105rgsubscript𝑅𝐵1.8superscript105subscript𝑟𝑔R_{B}={1.8\times 10^{5}}\,r_{g}italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT = 1.8 × 10 start_POSTSUPERSCRIPT 5 end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT, and we consider a non-spinning SMBH (Schwarzschild metric).

2 Numerical Methods

Our numerical scheme utilizes the GRMHD code KHARMA111https://github.com/AFD-Illinois/kharma, a performance-portable C++ implementation based on iharm3D (Prather et al., 2021); iharm3D is itself an extension of HARM, an efficient second-order 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 rRBmuch-greater-than𝑟subscript𝑅𝐵r\gg R_{B}italic_r ≫ italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT down to rgsubscript𝑟𝑔r_{g}italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT as well as from rgsubscript𝑟𝑔r_{g}italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT out to rRBmuch-greater-than𝑟subscript𝑅𝐵r\gg R_{B}italic_r ≫ italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT. 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=rg𝑟subscript𝑟𝑔r=r_{g}italic_r = italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT and rRBmuch-greater-than𝑟subscript𝑅𝐵r\gg R_{B}italic_r ≫ italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT into a set of n𝑛nitalic_n overlapping annuli evenly spaced in log(r)𝑟\log(r)roman_log ( italic_r ), each of which we label as zone-i𝑖iitalic_i. For i{0,1,,(n1)}𝑖01𝑛1i\in\{0,1,...,(n-1)\}italic_i ∈ { 0 , 1 , … , ( italic_n - 1 ) }, zone-i𝑖iitalic_i extends from an inner radius ri,in=8irgsubscript𝑟𝑖insuperscript8𝑖subscript𝑟𝑔r_{i,{\rm in}}=8^{i}\,r_{g}italic_r start_POSTSUBSCRIPT italic_i , roman_in end_POSTSUBSCRIPT = 8 start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT to an outer radius ri,out=8i+2rgsubscript𝑟𝑖outsuperscript8𝑖2subscript𝑟𝑔r_{i,{\rm out}}=8^{i+2}\,r_{g}italic_r start_POSTSUBSCRIPT italic_i , roman_out end_POSTSUBSCRIPT = 8 start_POSTSUPERSCRIPT italic_i + 2 end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT. We start by initializing and running zone-(n1)𝑛1(n-1)( italic_n - 1 ), the largest annulus, for some fraction of its free-fall time. We then simulate zone-(n2)𝑛2(n-2)( italic_n - 2 ), inheriting the final state of zone-(n1)𝑛1(n-1)( italic_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.222A few details of the numerical method are given in Section B. A complete description is deferred to a future paper.

3 Hydrodynamic Bondi accretion

Figure 1: Time-averaged radial profiles of selected quantities for our HD (red) and MHD (black) Bondi accretion simulations with RB=1.8×105rgsubscript𝑅𝐵1.8superscript105subscript𝑟𝑔R_{B}={1.8\times 10^{5}}\,r_{g}italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT = 1.8 × 10 start_POSTSUPERSCRIPT 5 end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT (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𝑛7n=7italic_n = 7 annuli with an outer radius of 882×107rgsimilar-tosuperscript882superscript107subscript𝑟𝑔8^{8}\sim 2\times 10^{7}\,r_{g}8 start_POSTSUPERSCRIPT 8 end_POSTSUPERSCRIPT ∼ 2 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT and a total runtime of trun(RB)9tBsimilar-tosubscript𝑡runsubscript𝑅𝐵9subscript𝑡𝐵t_{\rm run}(R_{B})\sim 9\,t_{B}italic_t start_POSTSUBSCRIPT roman_run end_POSTSUBSCRIPT ( italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT ) ∼ 9 italic_t start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT (0.7Myrsimilar-toabsent0.7Myr\sim 0.7\,{\rm Myr}∼ 0.7 roman_Myr for M87) at the Bondi radius. We confirmed convergence by checking that the results at trun=9tBsubscript𝑡run9subscript𝑡𝐵t_{\rm run}=9\,t_{B}italic_t start_POSTSUBSCRIPT roman_run end_POSTSUBSCRIPT = 9 italic_t start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT remain unchanged even when we run up to trun>56tBsubscript𝑡run56subscript𝑡𝐵t_{\rm run}>56\,t_{B}italic_t start_POSTSUBSCRIPT roman_run end_POSTSUBSCRIPT > 56 italic_t start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT. The resolution is 1283superscript1283128^{3}128 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT per annulus giving a total simulation resolution of 512×128×128512128128512\times 128\times 128512 × 128 × 128. Details of the set-up are presented in Appendix B.1.

Figures 1(a,b) show the radial profiles of t,θ,φ𝑡𝜃𝜑t,\theta,\varphiitalic_t , italic_θ , italic_φ-averaged gas density ρ𝜌\rhoitalic_ρ and temperature T𝑇Titalic_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 time-averaged accretion rate M˙(r)¯¯˙𝑀𝑟\overline{\dot{M}(r)}over¯ start_ARG over˙ start_ARG italic_M end_ARG ( italic_r ) end_ARG in Figure 1(c) matches the analytical Bondi accretion rate M˙Bsubscript˙𝑀𝐵\dot{M}_{B}over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT very well, M˙¯/M˙B=1¯˙𝑀subscript˙𝑀𝐵1\overline{\dot{M}}/\dot{M}_{B}=1over¯ start_ARG over˙ start_ARG italic_M end_ARG end_ARG / over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT = 1, across all radii within RBsubscript𝑅𝐵R_{B}italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT, confirming a steady state. We note that Guo et al. (2023) have previously simulated 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 M˙(r)¯¯˙𝑀𝑟\overline{\dot{M}(r)}over¯ start_ARG over˙ start_ARG italic_M end_ARG ( italic_r ) end_ARG we can define a time-averaged energy inflow rate towards the BH, E˙(r)¯¯˙𝐸𝑟\overline{\dot{E}(r)}over¯ start_ARG over˙ start_ARG italic_E end_ARG ( italic_r ) end_ARG (Equation (A5)), but this is not useful for our purposes since it includes the rest mass energy of the accreted gas. Instead we consider E˙net(r)(M˙E˙)¯subscript˙𝐸net𝑟¯˙𝑀˙𝐸\dot{E}_{\rm net}(r)\equiv\overline{(\dot{M}-\dot{E})}over˙ start_ARG italic_E end_ARG start_POSTSUBSCRIPT roman_net end_POSTSUBSCRIPT ( italic_r ) ≡ over¯ start_ARG ( over˙ start_ARG italic_M end_ARG - over˙ start_ARG italic_E end_ARG ) end_ARG (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 M˙¯10M˙¯(r=10rg)subscript¯˙𝑀10¯˙𝑀𝑟10subscript𝑟𝑔\overline{\dot{M}}_{10}\equiv\overline{\dot{M}}(r=10\,r_{g})over¯ start_ARG over˙ start_ARG italic_M end_ARG end_ARG start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT ≡ over¯ start_ARG over˙ start_ARG italic_M end_ARG end_ARG ( italic_r = 10 italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT ), we then obtain a feedback efficiency ηtotsuperscript𝜂tot\eta^{\rm tot}italic_η start_POSTSUPERSCRIPT roman_tot end_POSTSUPERSCRIPT (Equation (A7))

ηtot(M˙E˙¯)/M˙¯10.superscript𝜂tot¯˙𝑀˙𝐸subscript¯˙𝑀10\eta^{\rm tot}\equiv(\overline{\dot{M}-\dot{E}})/\overline{\dot{M}}_{10}\,.italic_η start_POSTSUPERSCRIPT roman_tot end_POSTSUPERSCRIPT ≡ ( over¯ start_ARG over˙ start_ARG italic_M end_ARG - over˙ start_ARG italic_E end_ARG end_ARG ) / over¯ start_ARG over˙ start_ARG italic_M end_ARG end_ARG start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT . (1)

Following Appendix A, and noting that there is no magnetic field and the system is spherically symmetric, ηtot1+(1+γadu/ρ)utsuperscript𝜂tot11subscript𝛾ad𝑢𝜌subscript𝑢𝑡\eta^{\rm tot}\approx 1+(1+\gamma_{\rm ad}u/\rho)u_{t}italic_η start_POSTSUPERSCRIPT roman_tot end_POSTSUPERSCRIPT ≈ 1 + ( 1 + italic_γ start_POSTSUBSCRIPT roman_ad end_POSTSUBSCRIPT italic_u / italic_ρ ) italic_u start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT. At large radii, uρT/(γad1)𝑢𝜌subscript𝑇subscript𝛾ad1u\to\rho T_{\infty}/(\gamma_{\rm ad}-1)italic_u → italic_ρ italic_T start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT / ( italic_γ start_POSTSUBSCRIPT roman_ad end_POSTSUBSCRIPT - 1 ) and ut1subscript𝑢𝑡1u_{t}\to-1italic_u start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT → - 1, hence

ηtotγadT/(γad1)=1.5/RB8.4×106.superscript𝜂totsubscript𝛾adsubscript𝑇subscript𝛾ad11.5subscript𝑅𝐵8.4superscript106\eta^{\rm tot}\rightarrow-\gamma_{\rm ad}T_{\infty}/(\gamma_{\rm ad}-1)=-1.5/R% _{B}\approx-8.4\times 10^{-6}.italic_η start_POSTSUPERSCRIPT roman_tot end_POSTSUPERSCRIPT → - italic_γ start_POSTSUBSCRIPT roman_ad end_POSTSUBSCRIPT italic_T start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT / ( italic_γ start_POSTSUBSCRIPT roman_ad end_POSTSUBSCRIPT - 1 ) = - 1.5 / italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT ≈ - 8.4 × 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT . (2)

The agreement of the simulation results in Figure 1(d) with this theoretical expectation is excellent.

4 Magnetized Bondi Accretion

Figure 2: Top: Accretion rate (black) and magnetic flux parameter (blue) as a function of time for our MHD simulation. The mean accretion rate M˙¯105×103M˙Bsimilar-tosubscript¯˙𝑀105superscript103subscript˙𝑀𝐵\overline{\dot{M}}_{10}\sim 5\times 10^{-3}\dot{M}_{B}over¯ start_ARG over˙ start_ARG italic_M end_ARG end_ARG start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT ∼ 5 × 10 start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT and the mean magnetic flux at the horizon ϕb¯30similar-to¯subscriptitalic-ϕ𝑏30\overline{\phi_{b}}\sim 30over¯ start_ARG italic_ϕ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT end_ARG ∼ 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 β𝛽\betaitalic_β (left) and ρ𝜌\rhoitalic_ρ (right) in our simulation, spanning 8 orders of magnitude in spatial scale. Lines of constant magnetic flux representing φ𝜑\varphiitalic_φ-averaged magnetic fields are overlaid in black.

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𝑧zitalic_z-direction with plasma-βpg/pb1𝛽subscript𝑝𝑔subscript𝑝𝑏1\beta\equiv p_{g}/p_{b}\approx 1italic_β ≡ italic_p start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT / italic_p start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT ≈ 1 where pbsubscript𝑝𝑏p_{b}italic_p start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT is the magnetic pressure. Our choice of β𝛽\betaitalic_β 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𝛽10\beta\approx 10italic_β ≈ 10 (Appendix C). We use n=8𝑛8n=8italic_n = 8 zones, which allows us to probe out to r108rgsimilar-to𝑟superscript108subscript𝑟𝑔r\sim 10^{8}r_{g}italic_r ∼ 10 start_POSTSUPERSCRIPT 8 end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT (and reach steady state out to nearly 107rgsuperscript107subscript𝑟𝑔10^{7}r_{g}10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT). Our effective resolution is 576×128×128576128128576\times 128\times 128576 × 128 × 128 over the whole domain. Additional set-up 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 trun(RB)9tBsimilar-tosubscript𝑡runsubscript𝑅𝐵9subscript𝑡𝐵t_{\rm run}(R_{B})\sim 9t_{B}italic_t start_POSTSUBSCRIPT roman_run end_POSTSUBSCRIPT ( italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT ) ∼ 9 italic_t start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT 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 ρr1proportional-to𝜌superscript𝑟1\rho\propto r^{-1}italic_ρ ∝ italic_r start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT, 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 9tBsimilar-toabsent9subscript𝑡𝐵\sim 9t_{B}∼ 9 italic_t start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT, we were unable to obtain a reliable measurement of M˙¯¯˙𝑀\overline{\dot{M}}over¯ start_ARG over˙ start_ARG italic_M end_ARG end_ARG at r102rggreater-than-or-equivalent-to𝑟superscript102subscript𝑟𝑔r\gtrsim 10^{2}\ r_{g}italic_r ≳ 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT 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 M˙insubscript˙𝑀in\dot{M}_{\rm in}over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT (solid blue) and mass outflow rate M˙outsubscript˙𝑀out\dot{M}_{\rm out}over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT (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 M˙¯¯˙𝑀\overline{\dot{M}}over¯ start_ARG over˙ start_ARG italic_M end_ARG end_ARG, which is the difference of the two very similar numbers. It is therefore challenging to calculate M˙¯¯˙𝑀\overline{\dot{M}}over¯ start_ARG over˙ start_ARG italic_M end_ARG end_ARG. Here, we restrict ourselves to smaller radii, where we believe our M˙¯¯˙𝑀\overline{\dot{M}}over¯ start_ARG over˙ start_ARG italic_M end_ARG end_ARG estimates are reliable. At r=10rg𝑟10subscript𝑟𝑔r=10\,r_{g}italic_r = 10 italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT 333It is customary in GRMHD simulations not to measure M˙˙𝑀\dot{M}over˙ start_ARG italic_M end_ARG at the horizon, where numerical floors might affect the fidelity of the computed density, but to measure it at a larger radius. We determined that r=10rg𝑟10subscript𝑟𝑔r=10r_{g}italic_r = 10 italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT is a safe choice for our work. we find M˙¯100.005M˙Bsubscript¯˙𝑀100.005subscript˙𝑀𝐵\overline{\dot{M}}_{10}\approx 0.005\dot{M}_{B}over¯ start_ARG over˙ start_ARG italic_M end_ARG end_ARG start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT ≈ 0.005 over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT, 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 103102M˙Bsimilar-toabsentsuperscript103superscript102subscript˙𝑀𝐵\sim 10^{-3}-10^{-2}\dot{M}_{B}∼ 10 start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT - 10 start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT in M87* (when comparing the BH accretion rate M˙103Msimilar-to˙𝑀superscript103subscript𝑀direct-product\dot{M}\sim 10^{-3}\,M_{\odot}over˙ start_ARG italic_M end_ARG ∼ 10 start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT/yr in Event Horizon Telescope Collaboration et al. (2021) with the Bondi accretion rate M˙B0.2Msimilar-tosubscript˙𝑀𝐵0.2subscript𝑀direct-product\dot{M}_{B}\sim 0.2\,M_{\odot}over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT ∼ 0.2 italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT/yr in Russell et al. (2015)). It is also consistent with the observed accretion rate for Sgr A* 104102M˙Bsuperscript104superscript102subscript˙𝑀𝐵10^{-4}-10^{-2}\,\dot{M}_{B}10 start_POSTSUPERSCRIPT - 4 end_POSTSUPERSCRIPT - 10 start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT (Wang et al., 2013).

We plot the time-averaged energy efficiency ηtot(r)superscript𝜂tot𝑟\eta^{\rm tot}(r)italic_η start_POSTSUPERSCRIPT roman_tot end_POSTSUPERSCRIPT ( italic_r ) (Equation 1) from the simulation in Figure 1(d). There are two striking features. First, the radial profile of ηtotsuperscript𝜂tot\eta^{\rm tot}italic_η start_POSTSUPERSCRIPT roman_tot end_POSTSUPERSCRIPT is nearly constant (ideally it should be perfectly constant), which indicates that the simulation energetics have converged well.444The simulations indicate a factor 2similar-toabsent2\sim 2∼ 2 increase in ηtotsuperscript𝜂tot\eta^{\rm tot}italic_η start_POSTSUPERSCRIPT roman_tot end_POSTSUPERSCRIPT close to the BH. We believe this is an artifact, caused by the need to apply numerical floors near the BH horizon. Second, ηtotsuperscript𝜂tot\eta^{\rm tot}italic_η start_POSTSUPERSCRIPT roman_tot end_POSTSUPERSCRIPT 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 102M˙¯10c25×105M˙Bc2similar-toabsentsuperscript102subscript¯˙𝑀10superscript𝑐2similar-to5superscript105subscript˙𝑀𝐵superscript𝑐2\sim 10^{-2}\overline{\dot{M}}_{10}c^{2}\sim 5\times 10^{-5}\,\dot{M}_{B}c^{2}∼ 10 start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT over¯ start_ARG over˙ start_ARG italic_M end_ARG end_ARG start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ∼ 5 × 10 start_POSTSUPERSCRIPT - 5 end_POSTSUPERSCRIPT over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, which is lower than the values adopted in cosmological simulations 103M˙Bc2greater-than-or-equivalent-toabsentsuperscript103subscript˙𝑀𝐵superscript𝑐2\gtrsim 10^{-3}\,\dot{M}_{B}c^{2}≳ 10 start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT (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 gas555Because of turbulent fluctuations, the gas at each radius tends to have a random amount of rotation, at the level of a few percent of Kepler. But the rotation is too weak to have any dynamical effect.. 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 β𝛽\betaitalic_β, and the normalized enclosed magnetic flux ϕbsubscriptitalic-ϕ𝑏\phi_{b}italic_ϕ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT (Tchekhovskoy et al., 2011),

ϕb(r)π/M˙¯10|Br|g𝑑θ𝑑φ.subscriptitalic-ϕ𝑏𝑟𝜋subscript¯˙𝑀10double-integralsuperscript𝐵𝑟𝑔differential-d𝜃differential-d𝜑\phi_{b}(r)\equiv\sqrt{\pi\big{/}\overline{\dot{M}}_{10}}\iint|B^{r}|\sqrt{-g}% \,d\theta\,d\varphi.italic_ϕ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT ( italic_r ) ≡ square-root start_ARG italic_π / over¯ start_ARG over˙ start_ARG italic_M end_ARG end_ARG start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT end_ARG ∬ | italic_B start_POSTSUPERSCRIPT italic_r end_POSTSUPERSCRIPT | square-root start_ARG - italic_g end_ARG italic_d italic_θ italic_d italic_φ . (3)

We average the inverse of β𝛽\betaitalic_β for stability, which is a commonly used technique. The profile we obtain has a flat plateau with a nearly constant value β3similar-to𝛽3\beta\sim 3italic_β ∼ 3 between r=102105𝑟superscript102superscript105r=10^{2}-10^{5}italic_r = 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - 10 start_POSTSUPERSCRIPT 5 end_POSTSUPERSCRIPT. Over this extended volume, the field is in equi-partition with the gas and the magnetic flux ϕbsubscriptitalic-ϕ𝑏\phi_{b}italic_ϕ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT varies rproportional-toabsent𝑟\propto r∝ italic_r, indicating that the organized poloidal field varies as Bpolr1proportional-tosubscript𝐵polsuperscript𝑟1B_{\rm pol}\propto r^{-1}italic_B start_POSTSUBSCRIPT roman_pol end_POSTSUBSCRIPT ∝ italic_r start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT. The pressure in the organized field is also in equi-partition: Bpol2ρTr2proportional-tosuperscriptsubscript𝐵pol2𝜌𝑇proportional-tosuperscript𝑟2B_{\rm pol}^{2}\propto\rho T\propto r^{-2}italic_B start_POSTSUBSCRIPT roman_pol end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ∝ italic_ρ italic_T ∝ italic_r start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT. We find that for radii r<102𝑟superscript102r<10^{2}italic_r < 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, β𝛽\betaitalic_β declines rapidly to well below unity and ϕbsubscriptitalic-ϕ𝑏\phi_{b}italic_ϕ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT saturates at 30similar-toabsent30\sim 30∼ 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 RB1.8×105rgsubscript𝑅𝐵1.8superscript105subscript𝑟𝑔R_{B}\approx{1.8\times 10^{5}}\,r_{g}italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT ≈ 1.8 × 10 start_POSTSUPERSCRIPT 5 end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT down to the BH horizon at 2rg2subscript𝑟𝑔2\,r_{g}2 italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT. One minor difference is that ϕb30similar-tosubscriptitalic-ϕ𝑏30\phi_{b}\sim 30italic_ϕ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT ∼ 30 at the BH horizon in our simulation, whereas typical GRMHD simulations which consider orbiting gas (instead of spherical infall) find ϕb50similar-tosubscriptitalic-ϕ𝑏50\phi_{b}\sim 50italic_ϕ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT ∼ 50 for a non-spinning BH (Tchekhovskoy et al., 2012; Narayan et al., 2022).

5 Feedback via Reconnection-driven Convection

Figure 3:
Refer to caption
Figure 4: Top: Dissecting contributions to the total energy outflow (ηtotsuperscript𝜂tot\eta^{\rm tot}italic_η start_POSTSUPERSCRIPT roman_tot end_POSTSUPERSCRIPT), including the fluid contribution (ηflsuperscript𝜂fl\eta^{\rm fl}italic_η start_POSTSUPERSCRIPT roman_fl end_POSTSUPERSCRIPT), the advective outflow efficiency (ηadvsuperscript𝜂adv\eta^{\rm adv}italic_η start_POSTSUPERSCRIPT roman_adv end_POSTSUPERSCRIPT), and the convective outflow efficiency (ηconvsuperscript𝜂conv\eta^{\rm conv}italic_η start_POSTSUPERSCRIPT roman_conv end_POSTSUPERSCRIPT). Since ηflηtotsuperscript𝜂flsuperscript𝜂tot\eta^{\rm fl}\approx\eta^{\rm tot}italic_η start_POSTSUPERSCRIPT roman_fl end_POSTSUPERSCRIPT ≈ italic_η start_POSTSUPERSCRIPT roman_tot end_POSTSUPERSCRIPT, 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,θ𝑟𝜃r,\thetaitalic_r , italic_θ plane, normalized by the accretion rate at 10rg10subscript𝑟𝑔10\;r_{g}10 italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT. Outflowing energy is shown in red and inflowing energy in blue.

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 M˙˙𝑀\dot{M}over˙ start_ARG italic_M end_ARG (in units of the Bondi rate M˙Bsubscript˙𝑀𝐵\dot{M}_{B}over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT) and the magnetic flux parameter ϕbsubscriptitalic-ϕ𝑏\phi_{b}italic_ϕ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT at the horizon (top panel of Figure 2) reveals large amplitude fluctuations in both quantities around a fairly stable mean value, M˙¯100.005M˙Bsimilar-tosubscript¯˙𝑀100.005subscript˙𝑀𝐵\overline{\dot{M}}_{\rm 10}\sim 0.005\dot{M}_{B}over¯ start_ARG over˙ start_ARG italic_M end_ARG end_ARG start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT ∼ 0.005 over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT and ϕ¯b30similar-tosubscript¯italic-ϕ𝑏30\overline{\phi}_{b}\sim 30over¯ start_ARG italic_ϕ end_ARG start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT ∼ 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 4 show the radial profile of the efficiency parameter ηtotsuperscript𝜂tot\eta^{\rm tot}italic_η start_POSTSUPERSCRIPT roman_tot end_POSTSUPERSCRIPT. The other curves show individual physical components: ηflsuperscript𝜂fl\eta^{\rm fl}italic_η start_POSTSUPERSCRIPT roman_fl end_POSTSUPERSCRIPT (blue, left top panel), the energy carried in the fluid (as opposed to the magnetic field), ηadvsuperscript𝜂adv\eta^{\rm adv}italic_η start_POSTSUPERSCRIPT roman_adv end_POSTSUPERSCRIPT (green, middle), the advected part of the fluid energy, and ηconvsuperscript𝜂conv\eta^{\rm conv}italic_η start_POSTSUPERSCRIPT roman_conv end_POSTSUPERSCRIPT (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 rRBless-than-or-similar-to𝑟subscript𝑅𝐵r\lesssim R_{B}italic_r ≲ italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT, we find that ηflsuperscript𝜂fl\eta^{\rm fl}italic_η start_POSTSUPERSCRIPT roman_fl end_POSTSUPERSCRIPT very closely matches ηtotsuperscript𝜂tot\eta^{\rm tot}italic_η start_POSTSUPERSCRIPT roman_tot end_POSTSUPERSCRIPT, 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 negligible. 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 RB=GM/cs,2=100subscript𝑅𝐵𝐺subscript𝑀superscriptsubscript𝑐𝑠2100R_{B}=GM_{\bullet}/c_{s,\infty}^{2}=100italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT = italic_G italic_M start_POSTSUBSCRIPT ∙ end_POSTSUBSCRIPT / italic_c start_POSTSUBSCRIPT italic_s , ∞ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT = 100 and a spinning BH (a*=0.9375subscript𝑎0.9375a_{*}=0.9375italic_a start_POSTSUBSCRIPT * end_POSTSUBSCRIPT = 0.9375), Ressler et al. (2021) found that the electromagnetic η𝜂\etaitalic_η exceeded ηflsuperscript𝜂fl\eta^{\rm fl}italic_η start_POSTSUPERSCRIPT roman_fl end_POSTSUPERSCRIPT. For the non-spinning case here, we find the opposite.

When we divide ηflsuperscript𝜂fl\eta^{\rm fl}italic_η start_POSTSUPERSCRIPT roman_fl end_POSTSUPERSCRIPT into an advective part ηadvsuperscript𝜂adv\eta^{\rm adv}italic_η start_POSTSUPERSCRIPT roman_adv end_POSTSUPERSCRIPT and a convective part ηconvsuperscript𝜂conv\eta^{\rm conv}italic_η start_POSTSUPERSCRIPT roman_conv end_POSTSUPERSCRIPT, we find that convection dominates by a substantial margin between r=rg𝑟subscript𝑟𝑔r=r_{g}italic_r = italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT and 0.1RBsimilar-toabsent0.1subscript𝑅𝐵\sim 0.1R_{B}∼ 0.1 italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT; 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 r0.1RBsimilar-to𝑟0.1subscript𝑅𝐵r\sim 0.1R_{B}italic_r ∼ 0.1 italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT, carrying the dissipated energy outward. Note that our definition of ηconvsuperscript𝜂conv\eta^{\rm conv}italic_η start_POSTSUPERSCRIPT roman_conv end_POSTSUPERSCRIPT 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 4 present detailed information on the fluid, advective and convective energy fluxes, FEflsuperscriptsubscript𝐹𝐸flF_{E}^{\rm fl}italic_F start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_fl end_POSTSUPERSCRIPT, FEadvsuperscriptsubscript𝐹𝐸advF_{E}^{\rm adv}italic_F start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_adv end_POSTSUPERSCRIPT, FEconvsuperscriptsubscript𝐹𝐸convF_{E}^{\rm conv}italic_F start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_conv end_POSTSUPERSCRIPT (defined in Appendix A), in the logarithmic radius - polar angle (log10rsubscript10𝑟\log_{10}{r}roman_log start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT italic_r - θ𝜃\thetaitalic_θ) plane. The advective flux carries energy toward the BH near the midplane but away from the BH near the poles. When integrated over θ𝜃\thetaitalic_θ, the two zones nearly cancel. On the other hand, convection transfers energy outward at all θ𝜃\thetaitalic_θ and overwhelms the advective energy at all r0.1RBless-than-or-similar-to𝑟0.1subscript𝑅𝐵r\lesssim 0.1R_{B}italic_r ≲ 0.1 italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT so that the net fluid flux FEflsuperscriptsubscript𝐹𝐸flF_{E}^{\rm fl}italic_F start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_fl end_POSTSUPERSCRIPT 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 4 that the convective energy transport ceases at RBsubscript𝑅𝐵R_{B}italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT, as expected since convection requires gravity and buoyancy, while gravity becomes sub-dominant for r>RB𝑟subscript𝑅𝐵r>R_{B}italic_r > italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT. 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).

6 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, RB1.8×105rgsubscript𝑅𝐵1.8superscript105subscript𝑟𝑔R_{B}\approx{1.8\times 10^{5}}r_{g}italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT ≈ 1.8 × 10 start_POSTSUPERSCRIPT 5 end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT. For spherical accretion without magnetic fields, we find perfect agreement with the GR analytical Bondi solution, i.e., for r<RB𝑟subscript𝑅𝐵r<R_{B}italic_r < italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT, the density and temperature of the simulated gas scale as ρr3/2proportional-to𝜌superscript𝑟32\rho\propto r^{-3/2}italic_ρ ∝ italic_r start_POSTSUPERSCRIPT - 3 / 2 end_POSTSUPERSCRIPT, Tr1proportional-to𝑇superscript𝑟1T\propto r^{-1}italic_T ∝ italic_r start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT respectively, as expected, and the measured accretion rate also matches the analytical Bondi accretion rate M˙Bsubscript˙𝑀𝐵\dot{M}_{B}over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT. The thermal energy of the external gas, which is 1.5/RB1051.5subscript𝑅𝐵superscript1051.5/R_{B}\approx 10^{-5}1.5 / italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT ≈ 10 start_POSTSUPERSCRIPT - 5 end_POSTSUPERSCRIPT (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 RBsubscript𝑅𝐵R_{B}italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT changes to ρr1proportional-to𝜌superscript𝑟1\rho\propto r^{-1}italic_ρ ∝ italic_r start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT, and the mass accretion rate M˙˙𝑀\dot{M}over˙ start_ARG italic_M end_ARG into the BH drops dramatically to 5×103M˙Babsent5superscript103subscript˙𝑀𝐵\approx 5\times 10^{-3}\dot{M}_{B}≈ 5 × 10 start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT, consistent with observations of M87* and Sgr A*. Equi-partition of magnetic and gas pressure is achieved (β3similar-to𝛽3\beta\sim 3italic_β ∼ 3) over the radius range 102rgrRBless-than-or-similar-tosuperscript102subscript𝑟𝑔𝑟less-than-or-similar-tosubscript𝑅𝐵10^{2}r_{g}\lesssim r\lesssim R_{B}10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT ≲ italic_r ≲ italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT, 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 102M˙c2absentsuperscript102˙𝑀superscript𝑐2\approx 10^{-2}\dot{M}c^{2}≈ 10 start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT over˙ start_ARG italic_M end_ARG italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, i.e., with an efficiency ηtot0.01similar-tosuperscript𝜂tot0.01\eta^{\rm tot}\sim 0.01italic_η start_POSTSUPERSCRIPT roman_tot end_POSTSUPERSCRIPT ∼ 0.01. Expressed in terms of the Bondi accretion rate, the energy outflow rate is 5×105M˙Bc2similar-toabsent5superscript105subscript˙𝑀𝐵superscript𝑐2\sim 5\times 10^{-5}\,\dot{M}_{B}c^{2}∼ 5 × 10 start_POSTSUPERSCRIPT - 5 end_POSTSUPERSCRIPT over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. 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 θ𝜃\thetaitalic_θ. 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 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.

Acknowledgements

The authors HC, RN, PN, KS, AR, and KC acknowledge support from the Gordon and Betty Moore Foundation and the John Templeton Foundation via grants to the Black Hole Initiative at Harvard University. This work used Delta at the University of Illinois at Urbana-Champaign through allocation PHY230079 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296. This research used resources provided to BP by the Los Alamos National Laboratory Institutional Computing Program, which is supported by the U.S. Department of Energy National Nuclear Security Administration under Contract No. 89233218CNA000001. We thank the referee for a thorough reading of the manuscript and constructive comments that have helped improve it.

Appendix A GRMHD Primer and Definitions

Grid-based GRMHD codes such as KHARMA used in this work evolve a magnetized fluid in a stationary spacetime (often a BH) described by a spacetime metric gμνsubscript𝑔𝜇𝜈g_{\mu\nu}italic_g start_POSTSUBSCRIPT italic_μ italic_ν end_POSTSUBSCRIPT. We focus on a non-spinning BH, and hence use the Schwarzschild metric. The magnetized fluid is treated within the ideal MHD approximation with the stress-energy tensor Tνμsubscriptsuperscript𝑇𝜇𝜈T^{\mu}_{\nu}italic_T start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT written as (see Anile 1989; Komissarov 1999; Gammie et al. 2003 for details)

Tνμ=(ρ+γadu+b2)uμuν+[pg+(b2/2)]δνμbμbν,subscriptsuperscript𝑇𝜇𝜈𝜌subscript𝛾ad𝑢superscript𝑏2superscript𝑢𝜇subscript𝑢𝜈delimited-[]subscript𝑝𝑔superscript𝑏22subscriptsuperscript𝛿𝜇𝜈superscript𝑏𝜇subscript𝑏𝜈T^{\mu}_{\nu}=\left(\rho+\gamma_{\rm ad}u+b^{2}\right)u^{\mu}u_{\nu}+\left[p_{% g}+(b^{2}/2)\right]\delta^{\mu}_{\nu}-b^{\mu}b_{\nu},italic_T start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT = ( italic_ρ + italic_γ start_POSTSUBSCRIPT roman_ad end_POSTSUBSCRIPT italic_u + italic_b start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) italic_u start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT italic_u start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT + [ italic_p start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT + ( italic_b start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT / 2 ) ] italic_δ start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT - italic_b start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT italic_b start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT , (A1)

where uμsuperscript𝑢𝜇u^{\mu}italic_u start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT is the four-velocity of the gas, ρ𝜌\rhoitalic_ρ is its rest mass density, u𝑢uitalic_u is its internal energy density, and bμsuperscript𝑏𝜇b^{\mu}italic_b start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT is a four-vector which describes the magnetic field in the fluid frame (Anile, 1989; Komissarov, 1999) with bμuμ=0superscript𝑏𝜇subscript𝑢𝜇0b^{\mu}u_{\mu}=0italic_b start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT italic_u start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT = 0 and bμbμ=b2superscript𝑏𝜇subscript𝑏𝜇superscript𝑏2b^{\mu}b_{\mu}=b^{2}italic_b start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT italic_b start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT = italic_b start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. The gas pressure pg=(γad1)usubscript𝑝𝑔subscript𝛾ad1𝑢p_{g}=(\gamma_{\rm ad}-1)uitalic_p start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT = ( italic_γ start_POSTSUBSCRIPT roman_ad end_POSTSUBSCRIPT - 1 ) italic_u, where γad=5/3subscript𝛾ad53\gamma_{\rm ad}=5/3italic_γ start_POSTSUBSCRIPT roman_ad end_POSTSUBSCRIPT = 5 / 3 is the adiabatic index. The temperature T𝑇Titalic_T is defined in relativistic units such that pg=ρTsubscript𝑝𝑔𝜌𝑇p_{g}=\rho Titalic_p start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT = italic_ρ italic_T, i.e., T=(γad1)u/ρ𝑇subscript𝛾ad1𝑢𝜌T=(\gamma_{\rm ad}-1)u/\rhoitalic_T = ( italic_γ start_POSTSUBSCRIPT roman_ad end_POSTSUBSCRIPT - 1 ) italic_u / italic_ρ. Given the initial conditions of a system, the GRMHD code evolves the equation of mass conservation, (ρuμ);μ=0(\rho u^{\mu})_{;\mu}=0( italic_ρ italic_u start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT ) start_POSTSUBSCRIPT ; italic_μ end_POSTSUBSCRIPT = 0, and the equations of energy-momentum conservation, Tν;μμ=0superscriptsubscript𝑇𝜈𝜇𝜇0T_{\ \nu;\mu}^{\mu}=0italic_T start_POSTSUBSCRIPT italic_ν ; italic_μ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT = 0, along with the ideal MHD induction equation. The output data describe the time evolution of the primitive quantities, ρ𝜌\rhoitalic_ρ, u𝑢uitalic_u, uμsuperscript𝑢𝜇u^{\mu}italic_u start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT, bμsuperscript𝑏𝜇b^{\mu}italic_b start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT, on each grid cell. At a given instant of time and at a given location in the grid, the quantity ρurgdθdφ𝜌superscript𝑢𝑟𝑔𝑑𝜃𝑑𝜑\rho u^{r}\sqrt{-g}\,d\theta d\varphiitalic_ρ italic_u start_POSTSUPERSCRIPT italic_r end_POSTSUPERSCRIPT square-root start_ARG - italic_g end_ARG italic_d italic_θ italic_d italic_φ measures the rate at which rest mass flows out radially within a solid angle dθdφ𝑑𝜃𝑑𝜑d\theta d\varphiitalic_d italic_θ italic_d italic_φ, where g=det(gμν)𝑔detsubscript𝑔𝜇𝜈g={\rm det}(g_{\mu\nu})italic_g = roman_det ( italic_g start_POSTSUBSCRIPT italic_μ italic_ν end_POSTSUBSCRIPT ). The instantaneous net mass accretion rate through a sphere of radius r𝑟ritalic_r is then:

M˙(r)ρurg𝑑θ𝑑φ,˙𝑀𝑟double-integral𝜌superscript𝑢𝑟𝑔differential-d𝜃differential-d𝜑\dot{M}(r)\equiv-\iint\rho u^{r}\sqrt{-g}\,d\theta\,d\varphi,over˙ start_ARG italic_M end_ARG ( italic_r ) ≡ - ∬ italic_ρ italic_u start_POSTSUPERSCRIPT italic_r end_POSTSUPERSCRIPT square-root start_ARG - italic_g end_ARG italic_d italic_θ italic_d italic_φ , (A2)

where the negative sign is because we associate accretion with inflow, which corresponds to a negative value of ursuperscript𝑢𝑟u^{r}italic_u start_POSTSUPERSCRIPT italic_r end_POSTSUPERSCRIPT. Because the accreting gas has large fluctuations in MHD (see the top panel in Figure 2), we are rarely interested in the instantaneous M˙˙𝑀\dot{M}over˙ start_ARG italic_M end_ARG. Rather, we average over an extended period of time (e.g., the latter half of the entire simulation), and consider the time-averaged mass accretion rate, which we write as M˙¯(r)¯˙𝑀𝑟\overline{\dot{M}}(r)over¯ start_ARG over˙ start_ARG italic_M end_ARG end_ARG ( italic_r ) with a bar. Furthermore, when we quote a unique accretion rate, we use the value at r=10rg𝑟10subscript𝑟𝑔r=10r_{g}italic_r = 10 italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT, which we write as M˙¯10subscript¯˙𝑀10\overline{\dot{M}}_{10}over¯ start_ARG over˙ start_ARG italic_M end_ARG end_ARG start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT. Sometimes it is useful to distinguish between mass that is flowing in at a given radius versus that flowing out (e.g., Fig. 1). In this case, we separately define:

M˙in(r)ρMin(ur,0)g𝑑θ𝑑φ¯,M˙out(r)ρMax(ur,0)g𝑑θ𝑑φ¯.formulae-sequencesubscript˙𝑀in𝑟¯double-integral𝜌Minsuperscript𝑢𝑟0𝑔differential-d𝜃differential-d𝜑subscript˙𝑀out𝑟¯double-integral𝜌Maxsuperscript𝑢𝑟0𝑔differential-d𝜃differential-d𝜑\dot{M}_{\rm in}(r)\equiv\overline{-\iint\rho\ {\rm Min}(u^{r},0)\sqrt{-g}\,d% \theta\,d\varphi},\qquad\dot{M}_{\rm out}(r)\equiv\overline{-\iint\rho\ {\rm Max% }(u^{r},0)\sqrt{-g}\,d\theta\,d\varphi}.over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT ( italic_r ) ≡ over¯ start_ARG - ∬ italic_ρ roman_Min ( italic_u start_POSTSUPERSCRIPT italic_r end_POSTSUPERSCRIPT , 0 ) square-root start_ARG - italic_g end_ARG italic_d italic_θ italic_d italic_φ end_ARG , over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT ( italic_r ) ≡ over¯ start_ARG - ∬ italic_ρ roman_Max ( italic_u start_POSTSUPERSCRIPT italic_r end_POSTSUPERSCRIPT , 0 ) square-root start_ARG - italic_g end_ARG italic_d italic_θ italic_d italic_φ end_ARG . (A3)

Per this definition, M˙in0subscript˙𝑀in0\dot{M}_{\rm in}\geq 0over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT ≥ 0 (inflow) and M˙out0subscript˙𝑀out0\dot{M}_{\rm out}\leq 0over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT ≤ 0 (outflow). Following the convention used in this paper, M˙insubscript˙𝑀in\dot{M}_{\rm in}over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT roman_in end_POSTSUBSCRIPT and M˙outsubscript˙𝑀out\dot{M}_{\rm out}over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT roman_out end_POSTSUBSCRIPT are shown in Figure 1 as solid lines (positive values) and dotted lines (negative values), respectively.

For the analyses in the paper, we define time-averaged and ρ𝜌\rhoitalic_ρ-weighted shell average of the quantity X𝑋Xitalic_X as:

X=(Xρg𝑑θ𝑑φρg𝑑θ𝑑φ)¯.delimited-⟨⟩𝑋¯double-integral𝑋𝜌𝑔differential-d𝜃differential-d𝜑double-integral𝜌𝑔differential-d𝜃differential-d𝜑\langle X\rangle=\overline{\left(\frac{\iint X\rho\sqrt{-g}~{}d\theta\,d% \varphi}{\iint\rho\sqrt{-g}~{}d\theta\,d\varphi}\right)}.⟨ italic_X ⟩ = over¯ start_ARG ( divide start_ARG ∬ italic_X italic_ρ square-root start_ARG - italic_g end_ARG italic_d italic_θ italic_d italic_φ end_ARG start_ARG ∬ italic_ρ square-root start_ARG - italic_g end_ARG italic_d italic_θ italic_d italic_φ end_ARG ) end_ARG . (A4)

The one exception is ρdelimited-⟨⟩𝜌\langle\rho\rangle⟨ italic_ρ ⟩, which we compute as a straight shell-average of the density without a further ρ𝜌\rhoitalic_ρ-weighting. In the shell averages, we leave out one closest cell from each of the two poles because these are affected by the polar boundary conditions.

For computing the flow of energy, we note that Ttrgdθdφsuperscriptsubscript𝑇𝑡𝑟𝑔𝑑𝜃𝑑𝜑T_{t}^{r}\sqrt{-g}\,d\theta d\varphiitalic_T start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_r end_POSTSUPERSCRIPT square-root start_ARG - italic_g end_ARG italic_d italic_θ italic_d italic_φ measures the rate at which energy flows in radially within a solid angle dθdφ𝑑𝜃𝑑𝜑d\theta d\varphiitalic_d italic_θ italic_d italic_φ. Thus the instantaneous energy inflow rate through a sphere at radius r𝑟ritalic_r is given by

E˙(r)Ttrg𝑑θ𝑑φ.˙𝐸𝑟double-integralsubscriptsuperscript𝑇𝑟𝑡𝑔differential-d𝜃differential-d𝜑\dot{E}(r)\equiv\iint T^{r}_{t}\sqrt{-g}\,d\theta\,d\varphi\,.over˙ start_ARG italic_E end_ARG ( italic_r ) ≡ ∬ italic_T start_POSTSUPERSCRIPT italic_r end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT square-root start_ARG - italic_g end_ARG italic_d italic_θ italic_d italic_φ . (A5)

When energy flows toward the BH, E˙˙𝐸\dot{E}over˙ start_ARG italic_E end_ARG is positive. For studying feedback we are interested in the flow of energy away from the BH. We are also not interested in the flow of rest mass energy, which has no effect on the thermodynamics. We thus define the time-averaged net rate of outflow of “useful” energy as:

E˙net(r)(M˙E˙)¯=(Ttr+ρur)¯g𝑑θ𝑑φ,subscript˙𝐸net𝑟¯˙𝑀˙𝐸double-integral¯subscriptsuperscript𝑇𝑟𝑡𝜌superscript𝑢𝑟𝑔differential-d𝜃differential-d𝜑\dot{E}_{\rm net}(r)\equiv\overline{(\dot{M}-\dot{E})}=-\iint\overline{(T^{r}_% {t}+\rho u^{r})}\sqrt{-g}\,d\theta\,d\varphi\,,over˙ start_ARG italic_E end_ARG start_POSTSUBSCRIPT roman_net end_POSTSUBSCRIPT ( italic_r ) ≡ over¯ start_ARG ( over˙ start_ARG italic_M end_ARG - over˙ start_ARG italic_E end_ARG ) end_ARG = - ∬ over¯ start_ARG ( italic_T start_POSTSUPERSCRIPT italic_r end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT + italic_ρ italic_u start_POSTSUPERSCRIPT italic_r end_POSTSUPERSCRIPT ) end_ARG square-root start_ARG - italic_g end_ARG italic_d italic_θ italic_d italic_φ , (A6)

and the corresponding efficiency ηtotsuperscript𝜂tot\eta^{\rm tot}italic_η start_POSTSUPERSCRIPT roman_tot end_POSTSUPERSCRIPT as:

ηtot(r)=E˙net(r)/M˙¯10=(M˙E˙)¯/M˙¯10.superscript𝜂tot𝑟subscript˙𝐸net𝑟subscript¯˙𝑀10¯˙𝑀˙𝐸subscript¯˙𝑀10\eta^{\rm tot}(r)=\dot{E}_{\rm net}(r)/\overline{\dot{M}}_{10}=\overline{(\dot% {M}-\dot{E})}/\overline{\dot{M}}_{10}\,.italic_η start_POSTSUPERSCRIPT roman_tot end_POSTSUPERSCRIPT ( italic_r ) = over˙ start_ARG italic_E end_ARG start_POSTSUBSCRIPT roman_net end_POSTSUBSCRIPT ( italic_r ) / over¯ start_ARG over˙ start_ARG italic_M end_ARG end_ARG start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT = over¯ start_ARG ( over˙ start_ARG italic_M end_ARG - over˙ start_ARG italic_E end_ARG ) end_ARG / over¯ start_ARG over˙ start_ARG italic_M end_ARG end_ARG start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT . (A7)

To study the contribution of the fluid alone (ignoring the magnetic field), the relevant terms in the stress-energy tensor are:

(Tfl)tr(ρ+γadu)urut.subscriptsuperscriptsuperscript𝑇fl𝑟𝑡𝜌subscript𝛾ad𝑢superscript𝑢𝑟subscript𝑢𝑡\left(T^{\rm fl}\right)^{r}_{t}\equiv(\rho+\gamma_{\rm ad}u)u^{r}u_{t}.( italic_T start_POSTSUPERSCRIPT roman_fl end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT italic_r end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ≡ ( italic_ρ + italic_γ start_POSTSUBSCRIPT roman_ad end_POSTSUBSCRIPT italic_u ) italic_u start_POSTSUPERSCRIPT italic_r end_POSTSUPERSCRIPT italic_u start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT . (A8)

Then the time-averaged net fluid energy outflow rate is:

E˙netfl(r)(Tfl+ρur)tr¯g𝑑θ𝑑φ=Beρur¯g𝑑θ𝑑φ,subscriptsuperscript˙𝐸flnet𝑟double-integral¯subscriptsuperscriptsuperscript𝑇fl𝜌superscript𝑢𝑟𝑟𝑡𝑔differential-d𝜃differential-d𝜑double-integral¯𝐵𝑒𝜌superscript𝑢𝑟𝑔differential-d𝜃differential-d𝜑\dot{E}^{\rm fl}_{\rm net}(r)\equiv-\iint\overline{\left(T^{\rm fl}+\rho u^{r}% \right)^{r}_{t}}\sqrt{-g}\,d\theta\,d\varphi=\iint\overline{Be~{}\rho u^{r}}% \sqrt{-g}\,d\theta\,d\varphi\,,over˙ start_ARG italic_E end_ARG start_POSTSUPERSCRIPT roman_fl end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_net end_POSTSUBSCRIPT ( italic_r ) ≡ - ∬ over¯ start_ARG ( italic_T start_POSTSUPERSCRIPT roman_fl end_POSTSUPERSCRIPT + italic_ρ italic_u start_POSTSUPERSCRIPT italic_r end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT italic_r end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT end_ARG square-root start_ARG - italic_g end_ARG italic_d italic_θ italic_d italic_φ = ∬ over¯ start_ARG italic_B italic_e italic_ρ italic_u start_POSTSUPERSCRIPT italic_r end_POSTSUPERSCRIPT end_ARG square-root start_ARG - italic_g end_ARG italic_d italic_θ italic_d italic_φ , (A9)

where the relativistic Bernoulli parameter is (e.g., Penna et al., 2013, but not including the magnetic field)

Be=(1+γadu/ρ)ut1.𝐵𝑒1subscript𝛾ad𝑢𝜌subscript𝑢𝑡1Be=-\left(1+\gamma_{\rm ad}u/\rho\right)u_{t}-1.italic_B italic_e = - ( 1 + italic_γ start_POSTSUBSCRIPT roman_ad end_POSTSUBSCRIPT italic_u / italic_ρ ) italic_u start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT - 1 . (A10)

Following Ressler et al. (2021), we further sub-divide the fluid energy outflow rate into an advective part,

E˙netadv(r)Be¯ρur¯g𝑑θ𝑑φ,subscriptsuperscript˙𝐸advnet𝑟double-integral¯𝐵𝑒¯𝜌superscript𝑢𝑟𝑔differential-d𝜃differential-d𝜑\dot{E}^{\rm adv}_{\rm net}(r)\equiv\iint\overline{Be}~{}\overline{\rho u^{r}}% \sqrt{-g}\,d\theta\,d\varphi\,,over˙ start_ARG italic_E end_ARG start_POSTSUPERSCRIPT roman_adv end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_net end_POSTSUBSCRIPT ( italic_r ) ≡ ∬ over¯ start_ARG italic_B italic_e end_ARG over¯ start_ARG italic_ρ italic_u start_POSTSUPERSCRIPT italic_r end_POSTSUPERSCRIPT end_ARG square-root start_ARG - italic_g end_ARG italic_d italic_θ italic_d italic_φ , (A11)

and a convective part,

E˙netconv(r)E˙netfl(r)E˙netadv(r).subscriptsuperscript˙𝐸convnet𝑟subscriptsuperscript˙𝐸flnet𝑟subscriptsuperscript˙𝐸advnet𝑟\dot{E}^{\rm conv}_{\rm net}(r)\equiv\dot{E}^{\rm fl}_{\rm net}(r)-\dot{E}^{% \rm adv}_{\rm net}(r)\,.over˙ start_ARG italic_E end_ARG start_POSTSUPERSCRIPT roman_conv end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_net end_POSTSUBSCRIPT ( italic_r ) ≡ over˙ start_ARG italic_E end_ARG start_POSTSUPERSCRIPT roman_fl end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_net end_POSTSUBSCRIPT ( italic_r ) - over˙ start_ARG italic_E end_ARG start_POSTSUPERSCRIPT roman_adv end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_net end_POSTSUBSCRIPT ( italic_r ) . (A12)

In the same manner as in (A7), we define the corresponding efficiencies, ηfl(r)superscript𝜂fl𝑟\eta^{\rm fl}(r)italic_η start_POSTSUPERSCRIPT roman_fl end_POSTSUPERSCRIPT ( italic_r ), ηadv(r)superscript𝜂adv𝑟\eta^{\rm adv}(r)italic_η start_POSTSUPERSCRIPT roman_adv end_POSTSUPERSCRIPT ( italic_r ), ηconv(r)superscript𝜂conv𝑟\eta^{\rm conv}(r)italic_η start_POSTSUPERSCRIPT roman_conv end_POSTSUPERSCRIPT ( italic_r ).

In addition to the time-averaged and shell-integrated radial profiles E˙(r)˙𝐸𝑟\dot{E}(r)over˙ start_ARG italic_E end_ARG ( italic_r ) defined above, we also define local time- and azimuth-averaged energy fluxes as a function of r𝑟ritalic_r and θ𝜃\thetaitalic_θ:

FEfl(r,θ)superscriptsubscript𝐹𝐸fl𝑟𝜃\displaystyle F_{E}^{\rm fl}(r,\theta)italic_F start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_fl end_POSTSUPERSCRIPT ( italic_r , italic_θ ) \displaystyle\equiv (1/2π)Beρur¯g𝑑φ,12𝜋¯𝐵𝑒𝜌superscript𝑢𝑟𝑔differential-d𝜑\displaystyle(1/2\pi)\int\overline{Be~{}\rho u^{r}}\sqrt{-g}\,d\varphi,( 1 / 2 italic_π ) ∫ over¯ start_ARG italic_B italic_e italic_ρ italic_u start_POSTSUPERSCRIPT italic_r end_POSTSUPERSCRIPT end_ARG square-root start_ARG - italic_g end_ARG italic_d italic_φ , (A13)
FEadv(r,θ)superscriptsubscript𝐹𝐸adv𝑟𝜃\displaystyle F_{E}^{\rm adv}(r,\theta)italic_F start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_adv end_POSTSUPERSCRIPT ( italic_r , italic_θ ) \displaystyle\equiv (1/2π)Be¯ρur¯g𝑑φ,12𝜋¯𝐵𝑒¯𝜌superscript𝑢𝑟𝑔differential-d𝜑\displaystyle(1/2\pi)\int\overline{Be}~{}\overline{\rho u^{r}}\sqrt{-g}\,d\varphi,( 1 / 2 italic_π ) ∫ over¯ start_ARG italic_B italic_e end_ARG over¯ start_ARG italic_ρ italic_u start_POSTSUPERSCRIPT italic_r end_POSTSUPERSCRIPT end_ARG square-root start_ARG - italic_g end_ARG italic_d italic_φ , (A14)
FEconv(r,θ)superscriptsubscript𝐹𝐸conv𝑟𝜃\displaystyle F_{E}^{\rm conv}(r,\theta)italic_F start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_conv end_POSTSUPERSCRIPT ( italic_r , italic_θ ) \displaystyle\equiv FEfl(r,θ)FEadv(r,θ).superscriptsubscript𝐹𝐸fl𝑟𝜃superscriptsubscript𝐹𝐸adv𝑟𝜃\displaystyle F_{E}^{\rm fl}(r,\theta)-F_{E}^{\rm adv}(r,\theta).italic_F start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_fl end_POSTSUPERSCRIPT ( italic_r , italic_θ ) - italic_F start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_adv end_POSTSUPERSCRIPT ( italic_r , italic_θ ) . (A15)

Appendix B Numerical Set-up

B.1 Hydrodynamic Simulation Setup

In the hydrodynamic simulation described in Section 3, we use a domain with 7 zones (i=06𝑖06i=0-6italic_i = 0 - 6) covering the range [1,1.7×107]rg11.7superscript107subscript𝑟𝑔[1,1.7\times 10^{7}]\;r_{g}[ 1 , 1.7 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT ] italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT. 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 (xt,xr,xθ,xφ)superscript𝑥𝑡superscript𝑥𝑟superscript𝑥𝜃superscript𝑥𝜑(x^{t},x^{r},x^{\theta},x^{\varphi})( italic_x start_POSTSUPERSCRIPT italic_t end_POSTSUPERSCRIPT , italic_x start_POSTSUPERSCRIPT italic_r end_POSTSUPERSCRIPT , italic_x start_POSTSUPERSCRIPT italic_θ end_POSTSUPERSCRIPT , italic_x start_POSTSUPERSCRIPT italic_φ end_POSTSUPERSCRIPT ), where the radial coordinate is xrlogrsuperscript𝑥𝑟𝑟x^{r}\equiv\log{r}italic_x start_POSTSUPERSCRIPT italic_r end_POSTSUPERSCRIPT ≡ roman_log italic_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 ρ𝜌\rhoitalic_ρ, u𝑢uitalic_u and uμsuperscript𝑢𝜇u^{\mu}italic_u start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT, 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×108]rg11.3superscript108subscript𝑟𝑔[1,1.3\times 10^{8}]\;r_{g}[ 1 , 1.3 × 10 start_POSTSUPERSCRIPT 8 end_POSTSUPERSCRIPT ] italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT. 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+RB)/rsubscript𝜌init𝑟subscript𝜌0𝑟subscript𝑅𝐵𝑟\rho_{\rm init}(r)=\rho_{0}(r+R_{B})/ritalic_ρ start_POSTSUBSCRIPT roman_init end_POSTSUBSCRIPT ( italic_r ) = italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ( italic_r + italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT ) / italic_r such that it is r1proportional-toabsentsuperscript𝑟1\propto r^{-1}∝ italic_r start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT interior to RBsubscript𝑅𝐵R_{B}italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT and constant outside of RBsubscript𝑅𝐵R_{B}italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT. 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 ρr1proportional-to𝜌superscript𝑟1\rho\propto r^{-1}italic_ρ ∝ italic_r start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT relation in the presence of strong magnetic fields. However, the results are insensitive to the choice of ρinitsubscript𝜌init\rho_{\rm init}italic_ρ start_POSTSUBSCRIPT roman_init end_POSTSUBSCRIPT, 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,θ)=bzr(r+RB)sin2θ/2subscript𝐴𝜑𝑟𝜃subscript𝑏𝑧𝑟𝑟subscript𝑅𝐵superscript2𝜃2A_{\varphi}(r,\theta)=b_{z}r(r+R_{B})\sin^{2}{\theta}/2italic_A start_POSTSUBSCRIPT italic_φ end_POSTSUBSCRIPT ( italic_r , italic_θ ) = italic_b start_POSTSUBSCRIPT italic_z end_POSTSUBSCRIPT italic_r ( italic_r + italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT ) roman_sin start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_θ / 2 such that the plasma-β𝛽\betaitalic_β parameter is approximately constant over all radii. The resulting magnetic field is vertical outside RBsubscript𝑅𝐵R_{B}italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT, and becomes more radial inside RBsubscript𝑅𝐵R_{B}italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT. The field is normalized with the parameter bzsubscript𝑏𝑧b_{z}italic_b start_POSTSUBSCRIPT italic_z end_POSTSUBSCRIPT such that the initial βinit(r)1similar-tosubscript𝛽init𝑟1\beta_{\rm init}(r)\sim 1italic_β start_POSTSUBSCRIPT roman_init end_POSTSUBSCRIPT ( italic_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 θ𝜃\thetaitalic_θ coordinate away from coordinate axes and toward the midplane. We adopt a simplified version of “funky modified Kerr-Schild” (FMKS) coordinates (see e.g. θjsubscript𝜃𝑗\theta_{j}italic_θ start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT in the Appendix in Wong et al., 2021), eliminating the cylindrification term to obtain:

θ=Ny(1+(y/χt)αα+1)+π2,𝜃𝑁𝑦1superscript𝑦subscript𝜒𝑡𝛼𝛼1𝜋2\theta=Ny\left(1+\frac{(y/\chi_{t})^{\alpha}}{\alpha+1}\right)+\frac{\pi}{2},italic_θ = italic_N italic_y ( 1 + divide start_ARG ( italic_y / italic_χ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT end_ARG start_ARG italic_α + 1 end_ARG ) + divide start_ARG italic_π end_ARG start_ARG 2 end_ARG , (B1)

where Nπ/2(1+χtα/(1+α))1𝑁𝜋2superscript1superscriptsubscript𝜒𝑡𝛼1𝛼1N\equiv\pi/2(1+\chi_{t}^{-\alpha}/(1+\alpha))^{-1}italic_N ≡ italic_π / 2 ( 1 + italic_χ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT start_POSTSUPERSCRIPT - italic_α end_POSTSUPERSCRIPT / ( 1 + italic_α ) ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT is a normalization factor, y=2xθ1𝑦2superscript𝑥𝜃1y=2x^{\theta}-1italic_y = 2 italic_x start_POSTSUPERSCRIPT italic_θ end_POSTSUPERSCRIPT - 1, χt=0.8subscript𝜒𝑡0.8\chi_{t}=0.8italic_χ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT = 0.8, α=16𝛼16\alpha=16italic_α = 16, and xθ[0,1]superscript𝑥𝜃01x^{\theta}\in[0,1]italic_x start_POSTSUPERSCRIPT italic_θ end_POSTSUPERSCRIPT ∈ [ 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 ρ>106r3/2𝜌superscript106superscript𝑟32\rho>10^{-6}\,r^{-3/2}italic_ρ > 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT italic_r start_POSTSUPERSCRIPT - 3 / 2 end_POSTSUPERSCRIPT, internal energy u>108r5/2𝑢superscript108superscript𝑟52u>10^{-8}\,r^{-5/2}italic_u > 10 start_POSTSUPERSCRIPT - 8 end_POSTSUPERSCRIPT italic_r start_POSTSUPERSCRIPT - 5 / 2 end_POSTSUPERSCRIPT, temperature u/ρ<100𝑢𝜌100u/\rho<100italic_u / italic_ρ < 100, and magnetization b2/ρ<100superscript𝑏2𝜌100b^{2}/\rho<100italic_b start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT / italic_ρ < 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=10subscript𝛾max10\gamma_{\rm max}=10italic_γ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT = 10. While running each active annulus, we keep ρ𝜌\rhoitalic_ρ, u𝑢uitalic_u and uμsuperscript𝑢𝜇u^{\mu}italic_u start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT fixed at the inner and outer boundaries, as in the hydro case. In addition, we also keep bμsuperscript𝑏𝜇b^{\mu}italic_b start_POSTSUPERSCRIPT italic_μ end_POSTSUPERSCRIPT fixed and make sure that we preserve B=0𝐵0\nabla\cdot B=0∇ ⋅ italic_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 RB460subscript𝑅𝐵460R_{B}\approx 460italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT ≈ 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.

Appendix C Resolution and initial condition study

Figure 5: Comparison between simulations with varying resolutions, coordinate systems, and initial magnetic strength, described in Table 1. Each profile is averaged between trun(RB)=4.5tB9tBsubscript𝑡runsubscript𝑅𝐵4.5subscript𝑡𝐵9subscript𝑡𝐵t_{\rm run}(R_{B})=4.5\,t_{B}-9\,t_{B}italic_t start_POSTSUBSCRIPT roman_run end_POSTSUBSCRIPT ( italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT ) = 4.5 italic_t start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT - 9 italic_t start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT measured at the Bondi radius (gray vertical dashed line). Positive (negative) values are shown in solid (dotted) lines. Overall, we find good agreement between these simulations.

We have presented our fiducial MHD simulation results in Section 4 for a resolution of 1283superscript1283128^{3}128 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT per annulus with the FMKS coordinate system and initialized with ρinitr1proportional-tosubscript𝜌initsuperscript𝑟1\rho_{\rm init}\propto r^{-1}italic_ρ start_POSTSUBSCRIPT roman_init end_POSTSUBSCRIPT ∝ italic_r start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT at r<RB𝑟subscript𝑅𝐵r<R_{B}italic_r < italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT. 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. θgsubscript𝜃𝑔\theta_{g}italic_θ start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT 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:

θ=πxθ+12(1h)sin(2πxθ)𝜃𝜋superscript𝑥𝜃1212𝜋superscript𝑥𝜃\theta=\pi x^{\theta}+\frac{1}{2}(1-h)\sin{(2\pi x^{\theta})}italic_θ = italic_π italic_x start_POSTSUPERSCRIPT italic_θ end_POSTSUPERSCRIPT + divide start_ARG 1 end_ARG start_ARG 2 end_ARG ( 1 - italic_h ) roman_sin ( 2 italic_π italic_x start_POSTSUPERSCRIPT italic_θ end_POSTSUPERSCRIPT ) (C1)

Here we use h=0.30.3h=0.3italic_h = 0.3. In the case of h=11h=1italic_h = 1, it is equivalent to the uniform θ𝜃\thetaitalic_θ 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) 483superscript48348^{3}48 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT resolution per annulus and (ii) 643superscript64364^{3}64 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT resolution per annulus. Run (iii) with 963superscript96396^{3}96 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT resolution per annulus used FMKS but with a different α=14𝛼14\alpha=14italic_α = 14 which is less extreme than α=16𝛼16\alpha=16italic_α = 16 in the fiducial run. Run (iv) used a weaker initial magnetic field with β10similar-to𝛽10\beta\sim 10italic_β ∼ 10 and used MKS coordinates with 643superscript64364^{3}64 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT resolution per annulus. Runs (v) and (vi) were initialized with a different initial density profile compared to the other runs and used MKS with 643superscript64364^{3}64 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT 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 ρinitr3/2proportional-tosubscript𝜌initsuperscript𝑟32\rho_{\rm init}\propto r^{-3/2}italic_ρ start_POSTSUBSCRIPT roman_init end_POSTSUBSCRIPT ∝ italic_r start_POSTSUPERSCRIPT - 3 / 2 end_POSTSUPERSCRIPT for r<RB𝑟subscript𝑅𝐵r<R_{B}italic_r < italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT.

Radial profiles of all the quantities of interest for the six extra simulations are compared with the fiducial run in Figure 5. 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 10rg10subscript𝑟𝑔10\,r_{g}10 italic_r start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT of all five runs lies in the range of M˙¯(15)×103M˙B¯˙𝑀15superscript103subscript˙𝑀𝐵\overline{\dot{M}}\approx(1-5)\times 10^{-3}\,\dot{M}_{B}over¯ start_ARG over˙ start_ARG italic_M end_ARG end_ARG ≈ ( 1 - 5 ) × 10 start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT, which is consistent with observations of M87* (Event Horizon Telescope Collaboration et al., 2021; Russell et al., 2015) and Sgr A* (Wang et al., 2013). The outflow efficiency is around ηtot13%superscript𝜂tot1percent3\eta^{\rm tot}\approx 1-3\,\%italic_η start_POSTSUPERSCRIPT roman_tot end_POSTSUPERSCRIPT ≈ 1 - 3 %, where the highest value is from the lowest resolution (483superscript48348^{3}48 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT) run (i). Even though initialized with a weaker magnetic field, run (iv) still shows magnetic field saturation at β3similar-to𝛽3\beta\sim 3italic_β ∼ 3 at radii r=102104𝑟superscript102superscript104r=10^{2}-10^{4}italic_r = 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - 10 start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT. Near the horizon, run (iv) accreted the least amount of magnetic field as can be inferred from its larger β𝛽\betaitalic_β and smaller ϕbsubscriptitalic-ϕ𝑏\phi_{b}italic_ϕ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT. However, the difference is not large. The β𝛽\betaitalic_β-parameter decreases more steeply towards the BH for run (iii) and the fiducial run, both of which use FMKS coordinates. We suspect this is because the FMKS coordinates are characterized by larger cells near the poles where β1superscript𝛽1\beta^{-1}italic_β start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT tends to be larger, and this might have affected the shell average. Runs (v) and (vi), despite their very different initial density profiles, converge to a similar final density profile with ρr1proportional-to𝜌superscript𝑟1\rho\propto r^{-1}italic_ρ ∝ italic_r start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT. This demonstrates that our results are not dependent on the initial conditions.

Label Resolution Coordinates Initial β𝛽\betaitalic_β Initial ρ𝜌\rhoitalic_ρ (r<RB)𝑟subscript𝑅𝐵(r<R_{B})( italic_r < italic_R start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT )
1283superscript1283128^{3}128 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT, fmks (fiducial) 1283superscript1283128^{3}128 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT FMKS (χt=0.8subscript𝜒𝑡0.8\chi_{t}=0.8italic_χ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT = 0.8, α=16𝛼16\alpha=16italic_α = 16) 1similar-toabsent1\sim 1∼ 1 r1superscript𝑟1r^{-1}italic_r start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
(i) 483superscript48348^{3}48 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT 483superscript48348^{3}48 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT MKS (h=0.30.3h=0.3italic_h = 0.3) 1similar-toabsent1\sim 1∼ 1 r1superscript𝑟1r^{-1}italic_r start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
(ii) 643superscript64364^{3}64 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT 643superscript64364^{3}64 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT MKS (h=0.30.3h=0.3italic_h = 0.3) 1similar-toabsent1\sim 1∼ 1 r1superscript𝑟1r^{-1}italic_r start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
(iii) 963superscript96396^{3}96 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT, fmks 963superscript96396^{3}96 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT FMKS (χt=0.8subscript𝜒𝑡0.8\chi_{t}=0.8italic_χ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT = 0.8, α=14𝛼14\alpha=14italic_α = 14) 1similar-toabsent1\sim 1∼ 1 r1superscript𝑟1r^{-1}italic_r start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
(iv) 643superscript64364^{3}64 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT, β10𝛽10\beta 10italic_β 10 643superscript64364^{3}64 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT MKS (h=0.30.3h=0.3italic_h = 0.3) 10similar-toabsent10\sim 10∼ 10 r1superscript𝑟1r^{-1}italic_r start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
(v) 643superscript64364^{3}64 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT, const 643superscript64364^{3}64 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT MKS (h=0.30.3h=0.3italic_h = 0.3) 1similar-toabsent1\sim 1∼ 1 piecewise constant
(vi) 643superscript64364^{3}64 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT, bondi 643superscript64364^{3}64 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT MKS (h=0.30.3h=0.3italic_h = 0.3) 1similar-toabsent1\sim 1∼ 1 r3/2superscript𝑟32r^{-3/2}italic_r start_POSTSUPERSCRIPT - 3 / 2 end_POSTSUPERSCRIPT
Table 1: Simulation set-up for different runs.

References

  • Anglés-Alcázar et al. (2021) Anglés-Alcázar, D., Quataert, E., Hopkins, P. F., et al. 2021, ApJ, 917, 53, doi: 10.3847/1538-4357/ac09e8
  • Anile (1989) Anile, A. M. 1989, Relativistic fluids and magneto-fluids: With applications in astrophysics and plasma physics
  • Bisnovatyi-Kogan & Ruzmaikin (1974) Bisnovatyi-Kogan, G. S., & Ruzmaikin, A. A. 1974, Ap&SS, 28, 45, doi: 10.1007/BF00642237
  • Bondi (1952) Bondi, H. 1952, MNRAS, 112, 195, doi: 10.1093/mnras/112.2.195
  • Booth & Schaye (2009) Booth, C. M., & Schaye, J. 2009, MNRAS, 398, 53, doi: 10.1111/j.1365-2966.2009.15043.x
  • Chatterjee & Narayan (2022) Chatterjee, K., & Narayan, R. 2022, ApJ, 941, 30, doi: 10.3847/1538-4357/ac9d97
  • Chatterjee et al. (2023) Chatterjee, K., Chael, A., Tiede, P., et al. 2023, Galaxies, 11, 38, doi: 10.3390/galaxies11020038
  • Davé et al. (2019) Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827, doi: 10.1093/mnras/stz937
  • Event Horizon Telescope Collaboration et al. (2021) Event Horizon Telescope Collaboration, Akiyama, K., Algaba, J. C., et al. 2021, ApJ, 910, L13, doi: 10.3847/2041-8213/abe4de
  • Ferrarese & Merritt (2000) Ferrarese, L., & Merritt, D. 2000, ApJ, 539, L9, doi: 10.1086/312838
  • Ferrière (2020) Ferrière, K. 2020, Plasma Physics and Controlled Fusion, 62, 014014, doi: 10.1088/1361-6587/ab49eb
  • Fiacconi et al. (2018) Fiacconi, D., Sijacki, D., & Pringle, J. E. 2018, MNRAS, 477, 3807, doi: 10.1093/mnras/sty893
  • Gammie et al. (2003) Gammie, C. F., McKinney, J. C., & Tóth, G. 2003, ApJ, 589, 444, doi: 10.1086/374594
  • Gebhardt et al. (2000) Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJ, 539, L13, doi: 10.1086/312840
  • Grete et al. (2023) Grete, P., Dolence, J. C., Miller, J. M., et al. 2023, 37, 465, doi: 10.1177/10943420221143775
  • Guerra et al. (2023) Guerra, J. A., Lopez-Rodriguez, E., Chuss, D. T., Butterfield, N. O., & Schmelz, J. T. 2023, AJ, 166, 37, doi: 10.3847/1538-3881/acdacd
  • Guo et al. (2023) Guo, M., Stone, J. M., Kim, C.-G., & Quataert, E. 2023, ApJ, 946, 26, doi: 10.3847/1538-4357/acb81e
  • Hopkins & Quataert (2010) Hopkins, P. F., & Quataert, E. 2010, MNRAS, 407, 1529, doi: 10.1111/j.1365-2966.2010.17064.x
  • Hopkins et al. (2023) Hopkins, P. F., Grudic, M. Y., Su, K.-Y., et al. 2023, arXiv e-prints, arXiv:2309.13115. https://arxiv.org/abs/2309.13115
  • Igumenshchev & Narayan (2002) Igumenshchev, I. V., & Narayan, R. 2002, ApJ, 566, 137, doi: 10.1086/338077
  • Igumenshchev et al. (2003) Igumenshchev, I. V., Narayan, R., & Abramowicz, M. A. 2003, ApJ, 592, 1042, doi: 10.1086/375769
  • Kaaz et al. (2023) Kaaz, N., Murguia-Berthier, A., Chatterjee, K., Liska, M. T. P., & Tchekhovskoy, A. 2023, ApJ, 950, 31, doi: 10.3847/1538-4357/acc7a1
  • Komissarov (1999) Komissarov, S. S. 1999, MNRAS, 303, 343, doi: 10.1046/j.1365-8711.1999.02244.x
  • Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511, doi: 10.1146/annurev-astro-082708-101811
  • Lalakos et al. (2022) Lalakos, A., Gottlieb, O., Kaaz, N., et al. 2022, ApJ, 936, L5, doi: 10.3847/2041-8213/ac7bed
  • Li & Bryan (2014) Li, Y., & Bryan, G. L. 2014, ApJ, 789, 153, doi: 10.1088/0004-637X/789/2/153
  • Magorrian et al. (1998) Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285, doi: 10.1086/300353
  • Meszaros (1975) Meszaros, P. 1975, A&A, 44, 59
  • Michel (1972) Michel, F. C. 1972, Ap&SS, 15, 153, doi: 10.1007/BF00649949
  • Narayan et al. (2022) Narayan, R., Chael, A., Chatterjee, K., Ricarte, A., & Curd, B. 2022, MNRAS, 511, 3795, doi: 10.1093/mnras/stac285
  • Narayan et al. (2000) Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A. 2000, ApJ, 539, 798, doi: 10.1086/309268
  • Narayan et al. (2003) —. 2003, PASJ, 55, L69, doi: 10.1093/pasj/55.6.L69
  • Ni et al. (2022) Ni, Y., Di Matteo, T., Bird, S., et al. 2022, MNRAS, 513, 670, doi: 10.1093/mnras/stac351
  • Pen et al. (2003) Pen, U.-L., Matzner, C. D., & Wong, S. 2003, ApJ, 596, L207, doi: 10.1086/379339
  • Penna et al. (2013) Penna, R. F., Kulkarni, A., & Narayan, R. 2013, A&A, 559, A116, doi: 10.1051/0004-6361/201219666
  • Porth et al. (2019) Porth, O., Chatterjee, K., Narayan, R., et al. 2019, ApJS, 243, 26, doi: 10.3847/1538-4365/ab29fd
  • Prather et al. (2021) Prather, B. S., Wong, G. N., Dhruv, V., et al. 2021, 6, 3336, doi: 10.21105/joss.03336
  • Quataert & Gruzinov (2000) Quataert, E., & Gruzinov, A. 2000, ApJ, 539, 809, doi: 10.1086/309267
  • Ressler et al. (2021) Ressler, S. M., Quataert, E., White, C. J., & Blaes, O. 2021, MNRAS, 504, 6076, doi: 10.1093/mnras/stab311
  • Ressler et al. (2020) Ressler, S. M., White, C. J., Quataert, E., & Stone, J. M. 2020, ApJ, 896, L6, doi: 10.3847/2041-8213/ab9532
  • Ricarte et al. (2019) Ricarte, A., Tremmel, M., Natarajan, P., & Quinn, T. 2019, MNRAS, 489, 802, doi: 10.1093/mnras/stz2161
  • Ripperda et al. (2022) Ripperda, B., Liska, M., Chatterjee, K., et al. 2022, ApJ, 924, L32, doi: 10.3847/2041-8213/ac46a1
  • Rosas-Guevara et al. (2016) Rosas-Guevara, Y., Bower, R. G., Schaye, J., et al. 2016, MNRAS, 462, 190, doi: 10.1093/mnras/stw1679
  • Russell et al. (2015) Russell, H. R., Fabian, A. C., McNamara, B. R., & Broderick, A. E. 2015, MNRAS, 451, 588, doi: 10.1093/mnras/stv954
  • Shapiro & Teukolsky (1983) Shapiro, S. L., & Teukolsky, S. A. 1983, Black holes, white dwarfs, and neutron stars : the physics of compact objects
  • Shvartsman (1971) Shvartsman, V. F. 1971, Soviet Ast., 15, 377
  • Sijacki et al. (2015) Sijacki, D., Vogelsberger, M., Genel, S., et al. 2015, MNRAS, 452, 575, doi: 10.1093/mnras/stv1340
  • Talbot et al. (2021) Talbot, R. Y., Bourne, M. A., & Sijacki, D. 2021, MNRAS, 504, 3619, doi: 10.1093/mnras/stab804
  • Tchekhovskoy et al. (2012) Tchekhovskoy, A., McKinney, J. C., & Narayan, R. 2012, in Journal of Physics Conference Series, Vol. 372, Journal of Physics Conference Series, 012040, doi: 10.1088/1742-6596/372/1/012040
  • Tchekhovskoy et al. (2011) Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2011, MNRAS, 418, L79, doi: 10.1111/j.1745-3933.2011.01147.x
  • Trott et al. (2022) Trott, C. R., Lebrun-Grandié, D., Arndt, D., et al. 2022, 33, 805, doi: 10.1109/TPDS.2021.3097283
  • Wang et al. (2013) Wang, Q. D., Nowak, M. A., Markoff, S. B., et al. 2013, Science, 341, 981, doi: 10.1126/science.1240755
  • Weinberger et al. (2018) Weinberger, R., Springel, V., Pakmor, R., et al. 2018, MNRAS, 479, 4056, doi: 10.1093/mnras/sty1733
  • Weinberger et al. (2023) Weinberger, R., Su, K.-Y., Ehlert, K., et al. 2023, MNRAS, 523, 1104, doi: 10.1093/mnras/stad1396
  • Wellons et al. (2023) Wellons, S., Faucher-Giguère, C.-A., Hopkins, P. F., et al. 2023, MNRAS, 520, 5394, doi: 10.1093/mnras/stad511
  • Wong et al. (2021) Wong, G. N., Du, Y., Prather, B. S., & Gammie, C. F. 2021, ApJ, 914, 55, doi: 10.3847/1538-4357/abf8b8
  • Yuan et al. (2015) Yuan, F., Gan, Z., Narayan, R., et al. 2015, ApJ, 804, 101, doi: 10.1088/0004-637X/804/2/101