Bridging Scales in Black Hole Accretion and Feedback: Magnetized Bondi Accretion in 3D GRMHD
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 orders of magnitude in spatial scale, to simulate accretion onto a non-spinning SMBH from an external medium with Bondi radius , where 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 rather than and the accretion rate is consequently suppressed by over 2 orders of magnitude below the Bondi rate . We find continuous energy feedback from the accretion flow to the external medium at a level of . 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.
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 pc, which translates to for the SMBH in M87; here is the gravitational radius, where is the SMBH mass. Meanwhile, probing smaller scales, general relativistic magnetohydrodynamic (GRMHD) simulations focus on the inner few tens or hundreds of , 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 , so our unit of length and time . Gas temperature is a dimensionless quantity which is normalized by , where is the mean molecular weight and is the Boltzmann constant; thus gas pressure with the density. We denote the time average of a quantity as and the density-weighted, time- and shell-averaged value of as (Appendix A). Finally, the Bondi radius is , where is the sound speed in the external medium far outside , is the temperature there, and is the gas adiabatic index. The free-fall time at the Bondi radius is . For the models presented here, , (), , 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 down to as well as from out to . 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 and into a set of overlapping annuli evenly spaced in , each of which we label as zone-. For , zone- extends from an inner radius to an outer radius . We start by initializing and running zone-, the largest annulus, for some fraction of its free-fall time. We then simulate zone-, inheriting the final state of zone- 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
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 annuli with an outer radius of and a total runtime of ( for M87) at the Bondi radius. We confirmed convergence by checking that the results at remain unchanged even when we run up to . The resolution is per annulus giving a total simulation resolution of . Details of the set-up are presented in Appendix B.1.
Figures 1(a,b) show the radial profiles of -averaged gas density and temperature (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 in Figure 1(c) matches the analytical Bondi accretion rate very well, , across all radii within , 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 we can define a time-averaged energy inflow rate towards the BH, (Equation (A5)), but this is not useful for our purposes since it includes the rest mass energy of the accreted gas. Instead we consider (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 , we then obtain a feedback efficiency (Equation (A7))
(1) |
Following Appendix A, and noting that there is no magnetic field and the system is spherically symmetric, . At large radii, and , hence
(2) |
The agreement of the simulation results in Figure 1(d) with this theoretical expectation is excellent.
4 Magnetized Bondi Accretion
Accretion flows in galactic nuclei involve magnetized plasma, and it has been long known that magnetic fields strongly perturb the energetics (Shvartsman, 1971; Meszaros, 1975) and dynamics (Igumenshchev & Narayan, 2002; Pen et al., 2003; Lalakos et al., 2022) of Bondi accretion. Apart from a toy model with purely radial magnetic field and a non-spinning BH, where no dynamical interaction occurs between the fluid and the magnetic field (Gammie et al., 2003), there is no known analytical solution to the magnetized accretion problem. It can be studied only via numerical simulations. Our computational technique is well-suited for this.
To simulate the magnetohydrodynamic (MHD) problem, we use the same external density and temperature as adopted for the HD case (Section 3), but now we include an initial magnetic field in the -direction with plasma- where is the magnetic pressure. Our choice of is motivated by the Milky Way ISM (e.g., Ferrière, 2020; Guerra et al., 2023), but we have verified that the results remain unchanged if we set the initial (Appendix C). We use zones, which allows us to probe out to (and reach steady state out to nearly ). Our effective resolution is 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 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 , 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 , we were unable to obtain a reliable measurement of at 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 (solid blue) and mass outflow rate (dotted cyan), defined in Appendix A, as a function of radius. With increasing distance from the BH, both rates become much larger than the net , which is the difference of the two very similar numbers. It is therefore challenging to calculate . Here, we restrict ourselves to smaller radii, where we believe our estimates are reliable. At 333It is customary in GRMHD simulations not to measure 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 is a safe choice for our work. we find , 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 in M87* (when comparing the BH accretion rate /yr in Event Horizon Telescope Collaboration et al. (2021) with the Bondi accretion rate /yr in Russell et al. (2015)). It is also consistent with the observed accretion rate for Sgr A* (Wang et al., 2013).
We plot the time-averaged energy efficiency (Equation 1) from the simulation in Figure 1(d). There are two striking features. First, the radial profile of is nearly constant (ideally it should be perfectly constant), which indicates that the simulation energetics have converged well.444The simulations indicate a factor increase in close to the BH. We believe this is an artifact, caused by the need to apply numerical floors near the BH horizon. Second, 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 , which is lower than the values adopted in cosmological simulations (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 , and the normalized enclosed magnetic flux (Tchekhovskoy et al., 2011),
(3) |
We average the inverse of for stability, which is a commonly used technique. The profile we obtain has a flat plateau with a nearly constant value between . Over this extended volume, the field is in equi-partition with the gas and the magnetic flux varies , indicating that the organized poloidal field varies as . The pressure in the organized field is also in equi-partition: . We find that for radii , declines rapidly to well below unity and saturates at . 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 down to the BH horizon at . One minor difference is that at the BH horizon in our simulation, whereas typical GRMHD simulations which consider orbiting gas (instead of spherical infall) find for a non-spinning BH (Tchekhovskoy et al., 2012; Narayan et al., 2022).
5 Feedback via Reconnection-driven Convection
As noted, our simulation permits us for the first time to study the coupled BH feeding and feedback problem. Following the time evolution of the BH mass accretion rate (in units of the Bondi rate ) and the magnetic flux parameter at the horizon (top panel of Figure 2) reveals large amplitude fluctuations in both quantities around a fairly stable mean value, and . 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 . The other curves show individual physical components: (blue, left top panel), the energy carried in the fluid (as opposed to the magnetic field), (green, middle), the advected part of the fluid energy, and (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 , we find that very closely matches , 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 and a spinning BH (), Ressler et al. (2021) found that the electromagnetic exceeded . For the non-spinning case here, we find the opposite.
When we divide into an advective part and a convective part , we find that convection dominates by a substantial margin between and ; 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 , carrying the dissipated energy outward. Note that our definition of 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, , , (defined in Appendix A), in the logarithmic radius - polar angle ( - ) plane. The advective flux carries energy toward the BH near the midplane but away from the BH near the poles. When integrated over , the two zones nearly cancel. On the other hand, convection transfers energy outward at all and overwhelms the advective energy at all so that the net fluid flux 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 , as expected since convection requires gravity and buoyancy, while gravity becomes sub-dominant for . 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, . For spherical accretion without magnetic fields, we find perfect agreement with the GR analytical Bondi solution, i.e., for , the density and temperature of the simulated gas scale as , respectively, as expected, and the measured accretion rate also matches the analytical Bondi accretion rate . The thermal energy of the external gas, which is (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 changes to , and the mass accretion rate into the BH drops dramatically to , consistent with observations of M87* and Sgr A*. Equi-partition of magnetic and gas pressure is achieved () over the radius range , 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 , i.e., with an efficiency . Expressed in terms of the Bondi accretion rate, the energy outflow rate is . This positive feedback occurs even in the absence of BH spin or net gas rotation. The energy outflow is not strongly collimated; it is quasi-isotropic. Advection carries energy inwards in the mid-plane and transports energy outwards near the poles, while convection transfers energy outwards at all . We find that the level of energy feedback is significantly lower than that typically invoked in cosmological simulations, so we view the feedback estimate in the present work as representing a lower limit to energy outflow. Additional ingredients, especially BH spin and gas angular momentum, are needed to study the problem more fully.
The dominant outward energy transport mechanism is convection driven by magnetic reconnection in the magnetosphere near the BH. Even though a strong magnetic field is essential for this mechanism to work, the energy is carried almost entirely in the form of thermal and kinetic energy of the gas, while the electromagnetic energy flux is negligible. Thus the feedback mechanism seen here is different from the electromagnetically-dominated 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 . 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 written as (see Anile 1989; Komissarov 1999; Gammie et al. 2003 for details)
(A1) |
where is the four-velocity of the gas, is its rest mass density, is its internal energy density, and is a four-vector which describes the magnetic field in the fluid frame (Anile, 1989; Komissarov, 1999) with and . The gas pressure , where is the adiabatic index. The temperature is defined in relativistic units such that , i.e., . Given the initial conditions of a system, the GRMHD code evolves the equation of mass conservation, , and the equations of energy-momentum conservation, , along with the ideal MHD induction equation. The output data describe the time evolution of the primitive quantities, , , , , on each grid cell. At a given instant of time and at a given location in the grid, the quantity measures the rate at which rest mass flows out radially within a solid angle , where . The instantaneous net mass accretion rate through a sphere of radius is then:
(A2) |
where the negative sign is because we associate accretion with inflow, which corresponds to a negative value of . Because the accreting gas has large fluctuations in MHD (see the top panel in Figure 2), we are rarely interested in the instantaneous . 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 with a bar. Furthermore, when we quote a unique accretion rate, we use the value at , which we write as . 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:
(A3) |
Per this definition, (inflow) and (outflow). Following the convention used in this paper, and 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 -weighted shell average of the quantity as:
(A4) |
The one exception is , which we compute as a straight shell-average of the density without a further -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 measures the rate at which energy flows in radially within a solid angle . Thus the instantaneous energy inflow rate through a sphere at radius is given by
(A5) |
When energy flows toward the BH, 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:
(A6) |
and the corresponding efficiency as:
(A7) |
To study the contribution of the fluid alone (ignoring the magnetic field), the relevant terms in the stress-energy tensor are:
(A8) |
Then the time-averaged net fluid energy outflow rate is:
(A9) |
where the relativistic Bernoulli parameter is (e.g., Penna et al., 2013, but not including the magnetic field)
(A10) |
Following Ressler et al. (2021), we further sub-divide the fluid energy outflow rate into an advective part,
(A11) |
and a convective part,
(A12) |
In the same manner as in (A7), we define the corresponding efficiencies, , , .
In addition to the time-averaged and shell-integrated radial profiles defined above, we also define local time- and azimuth-averaged energy fluxes as a function of and :
(A13) | |||||
(A14) | |||||
(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 () covering the range . 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 , where the radial coordinate is , 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 , and , 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 . 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 such that it is interior to and constant outside of . 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 relation in the presence of strong magnetic fields. However, the results are insensitive to the choice of , 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 such that the plasma- parameter is approximately constant over all radii. The resulting magnetic field is vertical outside , and becomes more radial inside . The field is normalized with the parameter such that the initial , motivated by observations of the Milky Way ISM (Ferrière, 2020; Guerra et al., 2023). A per-zone random perturbation of maximum 5% is applied to the internal energy when initializing fluid at all radii, to help break axisymmetry. Since the magnetic field picks out a rotation axis for the resulting flow parallel to the coordinate axis, we choose a coordinate system which compresses the coordinate away from coordinate axes and toward the midplane. We adopt a simplified version of “funky modified Kerr-Schild” (FMKS) coordinates (see e.g. in the Appendix in Wong et al., 2021), eliminating the cylindrification term to obtain:
(B1) |
where is a normalization factor, , , , and . 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 , internal energy , temperature , and magnetization , 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 . While running each active annulus, we keep , and fixed at the inner and outer boundaries, as in the hydro case. In addition, we also keep fixed and make sure that we preserve .
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 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
We have presented our fiducial MHD simulation results in Section 4 for a resolution of per annulus with the FMKS coordinate system and initialized with at . 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. 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:
(C1) |
Here we use . In the case of , it is equivalent to the uniform grid of the eKS coordinates in Section B.1.
Six extra simulations were run with parameters outlined in Table 1. Using MKS coordinates and other set-up details the same as in our fiducial run, we ran simulations at (i) resolution per annulus and (ii) resolution per annulus. Run (iii) with resolution per annulus used FMKS but with a different which is less extreme than in the fiducial run. Run (iv) used a weaker initial magnetic field with and used MKS coordinates with resolution per annulus. Runs (v) and (vi) were initialized with a different initial density profile compared to the other runs and used MKS with 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 for .
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 of all five runs lies in the range of , 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 , where the highest value is from the lowest resolution () run (i). Even though initialized with a weaker magnetic field, run (iv) still shows magnetic field saturation at at radii . Near the horizon, run (iv) accreted the least amount of magnetic field as can be inferred from its larger and smaller . However, the difference is not large. The -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 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 . This demonstrates that our results are not dependent on the initial conditions.
Label | Resolution | Coordinates | Initial | Initial |
---|---|---|---|---|
, fmks (fiducial) | FMKS (, ) | |||
(i) | MKS () | |||
(ii) | MKS () | |||
(iii) , fmks | FMKS (, ) | |||
(iv) , | MKS () | |||
(v) , const | MKS () | piecewise constant | ||
(vi) , bondi | MKS () |
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