Magnetically Modified Spherical Accretion in GRMHD: Reconnection-Driven Convection and Jet Propagation
Sean M. Ressler, Eliot Quataert, Christopher J. White, Omer Blaes
MMon. Not. R. Astron. Soc. , 000–000 (0000) Printed 4 February 2021 (MN L A TEX style file v2.2)
Magnetically Modified Spherical Accretion in GRMHD:Reconnection-Driven Convection and Jet Propagation
S. M. Ressler , E. Quataert , , C. J. White , O. Blaes Kavli Institute for Theoretical Physics, University of California Santa Barbara, Santa Barbara, CA 93107 Departments of Astronomy & Physics, Theoretical Astrophysics Center, University of California, Berkeley, CA 94720 Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544 Department of Physics ,University of California, Santa Barbara, CA 93106, USA
ABSTRACT
We present 3D general relativistic magnetohydrodynamic(GRMHD) simulations of zero an-gular momentum accretion around a rapidly rotating black hole, modified by the presenceof initially uniform magnetic fields. We consider serveral angles between the magnetic fielddirection and the black hole spin. In the resulting flows, the midplane dynamics are governedby magnetic reconnection-driven turbulence in a magnetically arrested (or a nearly arrested)state. Electromagnetic jets with outflow e ffi ciencies ∼ ∼ ◦ with respect to blackhole spin, but this tilt does not depend strongly on the tilt of the initial magnetic field. A jetforms even when there is no initial net vertical magnetic flux since turbulent, horizon-scalefluctuations can generate a net vertical field locally. Peak jet power is obtained for an initialmagnetic field tilted by 40–80 ◦ with respect to the black hole spin because this maximizes theamount of magnetic flux that can reach the black hole. These simulations may be a reasonablemodel for low luminosity black hole accretion flows such as Sgr A* or M87. Key words: accretion, accretion discs – (magnetohydrodynamics) MHD – black hole physics– convection – Galaxy: centre – magnetic reconnection
Well-resolved, three dimensional general relativistic magnetohy-drodynamic (GRMHD) simulations of accreting black holes arenow commonplace (De Villiers & Hawley 2003; Gammie, McK-inney & T´oth 2003; White, Stone & Gammie 2016; Porth et al.2017; Liska et al. 2019a) and are frequently used to model low lu-minosity accretion flows such as those surrounding the supermas-sive black holes Sagittarius A* (Sgr A*) and the one at the cen-ter of M87 (Mo´scibrodzka et al. 2009; Mo´scibrodzka, Falcke &Shiokawa 2016; Chan et al. 2015; Event Horizon Telescope Col-laboration et al. 2019a; Chael, Narayan & Johnson 2019; Dexteret al. 2020a). The most common strategy is to initialize the sim-ulations with a torus in hydrodynamic equilibrium (e.g., Fishbone& Moncrief 1976; Penna, Kulkarni & Narayan 2013) and add a dy-namically unimportant magnetic field to seed the magnetorotationalinstability (MRI, Balbus & Hawley 1991) and drive accretion. Byvarying the size of the torus, the initial magnetic field geometry,and the magnitude of the black hole spin, a relatively diverse set ofsolutions emerge that can, in principle, be used to constrain prop-erties of observed systems. In order to do so precisely, one mustunderstand not only the dependence of the results on these particu-lar parameters, but also on the assumptions that go into each model.For instance, much recent work has been devoted to exploring the e ff ects of non-ideal physics expected to be important for lowluminosity accretion flows. This includes two-temperature e ff ects(Ressler et al. 2015; Sa¸dowski et al. 2017; Chael et al. 2018), non-thermal electron acceleration (Mao, Dexter & Quataert 2017; Ballet al. 2016; Chael, Narayan & Sad¸owski 2017), radiation (Ryan,Dolence & Gammie 2015; Sa¸dowski et al. 2017; Ryan et al. 2017;Chael, Narayan & Johnson 2019), magnetic resistivity (Ripperdaet al. 2019), and anisotropic conduction / viscosity (Chandra, Fou-cart & Gammie 2015; Foucart et al. 2017). These works build onearlier 1D semi-analytic models that more easily incorporated extraphysics (e.g., Narayan & Yi 1995; ¨Ozel, Psaltis & Narayan 2000;Yuan, Quataert & Narayan 2003; Johnson & Quataert 2007). An-other active research frontier is that of tilted accretion disks, wherethe spin of the black hole is misaligned with the angular momen-tum of the flow (Fragile et al. 2007; Liska et al. 2018, 2019b;White, Quataert & Blaes 2019; Liska et al. 2020). Tilted disksnot only have di ff erent dynamics than non-tilted disks (e.g., pre-cession, standing shocks), they can also produce qualitatively dif-ferent images / spectra (e.g. ellipsoidal morphologies, larger imagesize, White et al. 2020; Chatterjee et al. 2020).Less explored are the e ff ects of the torus initial conditions.This was the motivation for a series of multi-scale simulationsof Sgr A* (Ressler, Quataert & Stone 2018, 2020; Ressler et al. © a r X i v : . [ a s t r o - ph . H E ] F e b S. M. Ressler, E. Quataert, C. J. White, O. Blaes ff ect onthe dynamics and the flow became magnetically arrested (Narayan,Igumenshchev & Abramowicz 2003; Igumenshchev, Narayan &Abramowicz 2003) on horizon scales (at least for a =
0, where a is the dimensionless spin of the black hole). Although we do nottypically have such intricate knowledge of the stellar / gas popula-tion surrounding supermassive black holes in other galaxies, it isnot unreasonable to assume that some fraction of these objects arefed in a similar way.In this work we seek to investigate an alternative model oflow angular momentum, magnetized accretion flows in GRMHD.In particular, we consider the limiting case of a non-rotating flowat large radii with an initially uniform magnetic field, the MHDanalog of the classic force-free Wald (1974) and Bicak & Janis(1985) solutions. In non-relativistic simulations (Igumenshchev &Narayan 2002; Pen, Matzner & Wong 2003; Pang et al. 2011),such flows have been shown to have distinct characteristics fromthose with significant angular momentum, e.g., turbulence drivenby magnetic reconnection and nearly hydrostatic force balance.For a non-spinning black hole, the qualitative nature of these so-lutions likely carries over to full GRMHD, though with order unityquantitative corrections near the horizon. Conversely, for a rapidlyspinning black hole, the large supply of coherent magnetic fluxwill presumably drive a Blandford & Znajek (1977) jet that wouldback-react on the quasi-spherical inflow and could significantly al-ter the resulting dynamics. Additionally, frame dragging and mag-netic torque could lead to non-negligible azimuthal velocities nearthe horizon even if there is initially no rotation.One potentially important parameter in these simulations isthe angle that the initial magnetic field makes with the black holespin axis. Previous works that study tilt have used equilibrium toriseeded with various magnetic field geometries contained withinthose tori (e.g., Fragile et al. 2007; Liska et al. 2018, 2019b; White,Quataert & Blaes 2019; Liska et al. 2020). This means that theorientation of the magnetic field was tilted along with the angularmomentum of the gas, making it hard to disentangle the e ff ects ofmagnetic tilt vs. angular momentum tilt. Here, starting with no an-gular momentum, we can investigate the magnetic tilt in isolationin order to determine how much, if any, it e ff ects the orientation ofboth the jet and the accretion flow. In particular, we focus on fourmagnetic tilt angles, 0 ◦ , ◦ , ◦ , and 90 ◦ , with an initial plasma β of 100, a Bondi radius of 200 gravitational radii, and a dimension-less black hole spin a = . Since non-radiative GRMHD simulations are scale-free, we willgenerally use length and time units that scale with the mass of the black hole and set the speed of light, c , and the gravitational con-stant, G , equal to unity. Thus, in these units, the mass of the blackhole, M , is our standard unit for both length and time, though weoften use the gravitational radii r g = M as an equivalent notation.Our simulations are carried out in Athena++ (White, Stone &Gammie 2016; Stone et al. 2020 in press), a 3 dimensional grid-based code that solves the equations of conservative GRMHD inCartesian Kerr–Schild coordinates (CKS, Kerr 1963). These relateto the spherical Kerr–Schild r , θ, ϕ coordinates via x = r sin( θ ) cos( ϕ ) + a sin( θ ) sin( ϕ ) y = r sin( θ ) sin( ϕ ) − a sin( θ ) cos( ϕ ) z = r cos( θ ) , (1)where a is the dimensionless spin of the black hole. Our computa-tional domain is a rectangular box of size 3200 r g × r g × r g with base resolution 24 × ×
192 cells in the x , y , and z directions.The box is larger in the z -direction in order to ensure that any rela-tivistic jet that might form does not interact with the outer bound-ary. We use an additional 12 levels of nested static mesh refinement(SMR) to achieve approximately logarithmic spacing in r , usingmeshblocks of size 12 × ×
16 cells. The outer boundary of the firstlevel is located at | z | = r g and | x | , | y | = r g , while the outerboundary of the second level is located at | z | = r g and | x | , | y | = r g . The outer boundaries for the remaining n (cid:62) | x | , | y | , | z | = / n − r g . The finest grid spacing is achievedfor | x | , | y | , | z | < . r g , where ∆ x , ∆ y , ∆ z ≈ . r g , ensuringthat the event horizon is well resolved since ( ∆ x , ∆ y , ∆ z ) / r H (cid:28) r H = + √ − a is the event horizon radius. All of oursimulations use a = . r H ≈ . r g . Note that for | x | , | y | , | z | (cid:54) r g , the grid is equivalent to a (1600 r g ) box with a192 base resolution with 9 levels of extra mesh refinement. Ourresolution is comparable to or better than the highest resolution(192 ) spherical KS torus simulations presented in the EHT codecomparison project that were found to be converged (Porth et al.2019). Specifically, at ( r = r g , θ = π/ ϕ direction, while having a comparableresolution in the θ direction to the most midplane-refined of theirsimulations.The initial conditions for these simulations are inspired by pre-vious non-relativistic models of magnetically modified sphericalaccretion in the literature (e.g., Igumenshchev & Narayan 2002;Pen, Matzner & Wong 2003; Pang et al. 2011). Outside of r = r g ,we start with a uniform density, ρ =
1, and uniform pressure, P , such that the Bondi radius is r B ≡ M / c s , = M , where c s , = (cid:112) γ P /ρ is the non-relativistic, adiabatic sound speed.These non-relativistic expressions are appropriate for large Bondiradii where GR e ff ects are small. For simplicity, the gas is station-ary, with the spatial components of the four-velocity set to zero.Inside r = r g , the density and pressure are set to the numericalfloors and again the spatial components of the four-velocity are setto zero. We set the initial magnetic field by taking the curl of avector potential A µ , B i = √− g (cid:15) ijk ∂∂ x j A k , (2)where g = − (cid:15) ijk isthe Levi–Cevita symbol. All spatial components of A µ are zero ex-cept for the y -component (note that the t -component is irrelevant),which we set to be A y ∝ e − / r ) (cid:2) − z sin( ψ ) + x cos( ψ ) (cid:3) , with thenormalization set such that β ≡ P / P B =
100 at large radii, where © , 000–000 P B = b / b = b µ b µ , and b µ is the mag-netic field four vector in Lorentz–Heaviside units. This vector po-tential corresponds to a uniform field that is inclined by an angle of ψ with respect to the black hole spin axis for r > r g . For r (cid:46) r g the field exponentially decays to zero, meaning that the field lineswrap around r ≈ r < r H /
2, we set the density / pressure to the numeri-cal floors and the four-velocity to free fall. Within this region theinduction equation is solved without alteration. Although CKS co-ordinates are horizon-penetrating, there is still a coordinate singu-larity at z = (cid:112) x + y = | a | (i.e., r = θ = π/ | z | < − r g we set z = − r g in the calculation of the metric.Since this is well inside the event horizon, it has no e ff ect on ourresults. The outer boundary conditions at the box faces are fixed tothe initial conditions.The density floor is 10 − ( r / r g ) − / and the pressure floor is3 . × − ( r / r g ) − / , with σ ≡ b /ρ (cid:54) = σ max and β (cid:62) . β is the ratio between the thermal and magnetic pressures.Additionally, the velocity of the gas is limited such that the maxi-mum Lorentz factor is 50 in the CKS normal frame.We run simulations with ψ = ◦ , ◦ , ◦ , and 90 ◦ for 20,000 M , about 10 free-fall times at the Bondi radius. We set the adia-batic index to γ = / x , y , and z to spher-ical KS coordinates r , θ , and ϕ using the inverse of Equations (1).Four-vectors are converted using the Jacobian: A µ (KS) = ∂ x µ (KS) ∂ x ν (CKS) A ν (CKS) . (3)We will also make use of an orthonomal tetrad to measure localvelocities in the zero angular momentum (ZAMO) frame, which wedefine in Boyer–Lindquist coordinates (Boyer & Lindquist 1967) asan observer with four-velocity η BL µ = − (cid:112) − g tt BL , , , . (4)The corresponding tetrad is e µ t , BL = η µ BL e µ r , BL = , (cid:112) g BL rr , , e µθ, BL = , , (cid:113) g BL θθ , e µφ, BL = , , , (cid:113) g BL ϕϕ , (5)so that the physical velocities in the ith direction are V i = u µ KS e i , KS µ u µ KS e t , KS µ , (6)where e µν, KS is e µν, BL converted to KS coordinates. The ψ = ◦ simulation uses a slightly di ff erent grid than the other threesimulations in that the | z | > r g region is excluded. Within | x | , | y | , | z | < r g the grids are identical. Figures 1–4 show 2D contour slices of the mass density over-plotted with magnetic field lines at four di ff erent times for our ψ = ◦ , 30 ◦ , 60 ◦ , and 90 ◦ simulations, respectively. In all cases,the initially weak magnetic field is dragged inwards with the flowuntil it becomes dynamically important. Absent black hole spin,this would result in a ‘pinched’ geometry of the field where matteris pushed away from the magnetic poles and towards the magneticmidplane. In non-relativistic MHD this ultimately leads to a highlydisordered state where heating associated with the reconnectionof the pinched field drives convection (Igumenshchev & Narayan2002; Pen, Matzner & Wong 2003; Pang et al. 2011). Here, how-ever, the e ff ects of the rapidly spinning black hole significantly al-ter this scenario. For ψ = ◦ , ◦ and 60 ◦ , the accumulation of fluxthreading the horizon (particularly flux parallel to the black holespin) caused by the initially spherical accretion leads to the forma-tion of a magnetically dominated jet by ∼ M via the Blandford& Znajek (1977) mechanism in which the magnetic field is tightlywound up in the ϕ -direction and propels matter / energy outwards.We find that the jet is always aligned with the black hole spin for r (cid:46) r g , though at larger radii it can be significantly tilted (we dis-cuss jet orientation in more detail in §3.3). The strong outflow evac-uates the polar regions of matter but leaves a broad distribution ofgas in the midplane. Turbulence driven by magnetic reconnectionis visible at small ( r (cid:46) r g ) radii in Figures 1–4.The ψ = ◦ simulation shows similar behavior between5000–15000 M but much di ff erent behavior at earlier and latertimes. Initially, the dynamics are governed by reconnection driventurbulence with no preferred direction, the GRMHD analog to thenon-relativistic simulations mentioned above. This state persistsuntil turbulent fluctuations spontaneously generate a net locallyvertical magnetic field that accretes and powers a jet. However,since this vertical field is not being continuously supplied by gasat large radii, the jet dies o ff and the convection dominated stateresumes. We expect that if the simulation were to run for a longertime then another turbulently-generated, locally vertical field wouldbe created and the cycle would repeat.The simulations also have interesting azimuthal structure. Fig-ure 5 shows plasma β in the midplane of the ψ = ◦ simulation asan example at four di ff erent times (the other three simulations dis-play similar contours and are thus not shown). The flow structure ishighly nonaxisymmetric, with clear spiral structure caused by therotating spacetime of the central Kerr black hole. Such structuresare largest and most coherent in the ψ = ◦ simulation since thetwisting of the initially horizontal field by the black hole providesthe most spatially extended torque on the gas. In all cases, mattergenerally accretes in thin streams with higher β ( (cid:38)
1) regions sur-rounded by lower β ( (cid:46)
1) magnetically dominated regions.Figure 6 shows time and angle averages of several importantquantities for the four simulations during the interval 10000–15000 M . These include the mass density, ρ ; the midplane angular velocityrelative to the prograde Keplerian value measured in the ZAMOframe, V ϕ / V kep , where V kep is computed by transforming ˙ ϕ BLKep = / ( r / + a ) to the ZAMO frame; the radial velocity divided bythe sound speed measured in the ZAMO frame, V r / c s , where c s = (cid:112) γ P / w and w = ρ + γ/ ( γ − P ; temperature, T g = P /ρ (assumingionized hydrogen); entropy, s = ( γ − − log( P /ρ γ ); and plasma β . © , 000–000 S. M. Ressler, E. Quataert, C. J. White, O. Blaes
Figure 1.
2D slice through the x – z plane of mass density over-plotted with magnetic field lines at four di ff erent times in the ψ = ◦ simulation. Starting fromthe top left panel and proceeding clockwise, the panels represent 0, 4000, 10000, and 18000 M . The initially uniform field accretes through the event horizonand the rapidly spinning black hole drives a powerful jet. The accreting field also reconnects, heating the gas and driving turbulence. An animated version ofthis figure is available at https: // youtu.be / N7wz-vPLwLo.
Here we compute angle averages in KS coordinates by evaluating (cid:104) A (cid:105) X = (cid:33) AXd Ω (cid:33) Xd Ω , (7)where d Ω = √− g KS d ϕ d θ and g KS = − sin( θ ) [ r + a cos( θ ) ] . Thequantities plotted in Figure 6 are typically mass-weighted in orderto study the properties of the accreting midplane. We will discussthe properties of the jet later in §3.3. Qualitatively, the four simu-lations display the same basic trends in the mass-weighted radialprofiles. Within r (cid:46) r g , the magnetic field has been compressedto the point that the flow is magnetically dominated, with β (cid:46) ∼ r g instead of the hydro- dynamic Bondi value of ∼ r g . Note, however, that this Lorentzforce is still (cid:46) the thermal pressure force because the magnetic pres-sure and tension terms tend to act in opposite directions, reducingthe overall net magnetic force on the gas. The strong field is recon-necting and heating the gas, as evidenced by the entropy increasingwith decreasing r and the associated temperature increase beyondwhat one would predict from adiabatic compression alone. Inside The hydrodynamic Bondi radius for a = a = . It is somewhat surprising to see the entropy continue to increase afterthe fluid passes through the sonic point ( r (cid:46) r g ) and plunges through theevent horizon at r ≈ r g because n this narrow radial range one wouldexpect the gas to essentially be experiencing adiabatic free-fall. This may © , 000–000 Figure 2.
Same as Figure 1 but for ψ = ◦ . the Bondi radius (200 r g ) the mass density settles into an ρ ˜ ∝ r − profile. Finally, despite the stationary, spherical initial conditions,the frame-dragging e ff ects of the a = . r (cid:46) r g , reaching up to ∼ ff ect in the sense that u ϕ is nonzero and u ϕ is larger than what onewould expect from frame dragging of the gas alone (in which case V ϕ would be zero). Compared to the gas in the ψ = ◦ and ψ = ◦ simulations, the gas in the ψ = ◦ and ψ = ◦ simulations hasa shallower rotation curve for r (cid:38) r g , meaning that significantangular velocity is present out to slightly larger radii. This is likely be because 1) averaging e ff ects result in an average radial entropy profilethat is not necessarily representative of the entropy profile along individualstreamlines or 2) the thermodynamics in this region may be less reliablethan at larger radii because σ (cid:38) because the larger amount of field perpendicular to the spin axisprovides a longer ‘lever arm’ to torque the flow.The biggest outlier in Figure 6 is the ψ = ◦ simulation,which has a ∼ ∼ (cid:46) ∼ P /ρ γ for r (cid:46) r g . As we will show in the next subsection, this is because theamount of flux able to reach the black hole for ψ = ◦ is largerthan the other three simulations. We speculate as to why in §4.In terms of energy flux, we find that convection is ine ffi cient attransporting energy outwards. Figure 7 demonstrates this by plot-ting the average total flux of energy contained in matter, F E , Ma = − π (cid:16) r + a (cid:17) (cid:42)(cid:32) + γγ − P ρ (cid:33) ρ u r u t − ρ u r (cid:43) , (8)compared to the average flux of matter energy through laminar ad-vection, F adv , Ma = − π (cid:16) r + a (cid:17) (cid:42)(cid:32) + γγ − P ρ (cid:33) u t − (cid:43) (cid:104) ρ u r (cid:105) , (9) © , 000–000 S. M. Ressler, E. Quataert, C. J. White, O. Blaes
Figure 3.
Same as Figure 1 but for ψ = ◦ . and the average convective flux, F conv , Ma = F E , Ma − F adv , Ma , for allfour simulations. In averaging these terms, we exclude the low den-sity jet by focussing only on | θ − π/ | < θ H , where θ H = (cid:104)| θ − π/ |(cid:105) ρ is the half opening angle of the mass distribution. Energy flows in-wards predominantly by laminar advection, while convection pro-vides only a small outwards contribution with | F conv , Ma | a factor of3–10 times smaller than | F adv , Ma | except for r (cid:46) a few r g where thetwo fluxes can be comparable. The lack of a strong outwards con-vective flux explains why the density profile in Figure 6 is steeperthan the r − / predicted by models where convection is the domi-nant energy transport mechanism (e.g., Narayan, Igumenshchev &Abramowicz 2000; Quataert & Gruzinov 2000; Igumenshchev &Narayan 2002). We now turn our attention to the temporal variability of the foursimulations. A key quantity in magnetized black hole accretion is the dimensionless measure of the magnetic flux threading the hori-zon, φ BH ≡ √ π (cid:33) | B r | d Ω (cid:112) | ˙ M | , (10)where B r is the radial component of the magnetic three-vector and˙ M is the accretion rate. In GRMHD torus simulations this quan-tity seems to have a maximum possible value at φ BH ≈ φ BH for our four simulations. The ψ = ◦ curve dis-plays quintessential magnetically arrested disk (MAD) character-istics, with φ BH varying between 30–60 in periodic cycles of ac-cumulation and dissipation. The accretion rate in the second panel © , 000–000 Figure 4.
Same as Figure 1 but for ψ = ◦ . At 4000 M (top right panel), the gas is in a state of disordered convection that persists until a net vertical magneticfield is randomly generated locally and a jet is formed. An animated version of this figure is available at https: // youtu.be / L3apzf2cTP4. of Figure 8 varies accordingly (between 0.05–0.35 times the Bondirate) with the peaks in φ BH associated with valleys in ˙ M , and viceversa. In contrast, the ψ = ◦ simulation does not seem to reachthe MAD state, with a saturated φ BH of ∼
30 that is much less vari-able than the ψ = ◦ case. The accretion rate is also more constantin time than ψ = ◦ , though the time-averaged values in both sim-ulations are comparable, ∼ ψ = ◦ simulation has an even smaller amount of magnetic flux reachingthe black hole, φ BH ≈
20 and roughly constant in time. This makessense since the jet seen in Figure 4 pushes perpendicular to the ini-tial magnetic field and e ff ectively limits φ BH . Because of this, thetime-averaged accretion rate is the largest of the four simulations, ∼ ψ = ◦ , which shows the largest amount of magnetic flux reach-ing the black hole, φ BH ∼ M ∼ ψ = ◦ curves. Naively, we wouldhave expected the results for ψ = ◦ simulation to fall somewhere in between the ψ = ◦ and the ψ = ◦ results with a relativelysmall amount of magnetic flux reaching the black hole. Instead, itseems that supplying a field that is significantly tilted from verticalultimately results in more net flux (as measured in both the absolutesense and normalized to the square root of the accretion rate, φ BH )than a field completely aligned with the spin of the black hole. Wespeculate on why in §4.The third panel of Figure 8 plots the outflow e ffi ciency, η = ( ˙ E − ˙ M ) / | ˙ M | . As expected, the simulation with the largest φ BH showsthe largest e ffi ciency ( ∼
200 % for ψ = ◦ ), while the simulationwith the smallest φ BH shows the smallest e ffi ciency (either ∼ ∼
10 % for ψ = ◦ ). An e ffi ciency of 200 % can only occur if en-ergy is being extracted from the spin of the black hole , somethingobserved in previous highly magnetized GRMHD simulations (e.g., For η < © , 000–000 S. M. Ressler, E. Quataert, C. J. White, O. Blaes
Figure 5.
Midplane slice of plasma β at four di ff erent times in our ψ = ◦ simulation overplotted with magnetic field lines. Starting from the top left paneland proceeding clockwise, the panels represent 2000, 8000, 14000, and 20000 M . Accretion proceeds along thin, high density streams surrounded by highlymagnetized, low β regions. Contours of β in the other three simulations look similar. Tchekhovskoy, Narayan & McKinney 2011). The temporal vari-ability of the e ffi ciency in the ψ = ◦ simulation di ff ers from theother three simulations in that there are both quiescent phases andan active phase. The active phase (with η ∼ M (cid:46) t (cid:46) M , while the quiescent phases (with η ∼ t (cid:46) M and t (cid:38) M . The quiescent phases corre-spond to the highly convective state discussed in §3.1 where a lackof significant vertical flux prevents the formation of a jet, insteadfavoring disordered reconnection-driven turbulence (see top rightand bottom left panel of Figure 4). The active phase arises afterenough locally net vertical magnetic field is created via turbulenceto form a jet. In this section we focus specifically on the magnetically dominatedjet. Figures 9–12 show 2D slices at four di ff erent times of the mag- netization parameter b /ρ zoomed out to a 200 r g × r g box inthe x – z plane in the four simulations. At early times (i.e., 2000 M in the left panel of Figures 9–12), the jets are relatively symmetricabout both the spin axis of the black hole and the midplane. As timegoes on, however, this symmetry is broken by external kink insta-bilities that create more complicated structures (see, e.g. Bromberg& Tchekhovskoy 2016; Tchekhovskoy & Bromberg 2016; BarniolDuran, Tchekhovskoy & Giannios 2017). At some times the jet issignificantly more extended in either the upper or lower directions(e.g., the two rightmost panels of Figure 9 and Figure 10 as well asthe rightmost panel of Figure 11) and can favor the left or right sideof the domain (e.g., the two rightmost panels fo Figure 9 and thethree rightmost panels of Figure 10). The latter phenomenon is par-ticularly striking in the second panel of Figure 10 ( ψ = ◦ ), wherethe upper and lower jet make an approximately 90 ◦ angle with eachother. While the upper jet is aligned with the direction of the initialmagnetic field (inclined by 30 ◦ from the spin axis), the lower jet is © , 000–000 r [ r g ] V / V k e p ( = / ) V / V kep r = 0= 30 = 60= 9010 r [ r g ] r r T T = 0= 30 = 60= 9010 r [ r g ] s s = 0= 30 = 60= 90 Figure 6.
Angle and time averaged fluid quantities plotted versus radiusin each of our four simulations. Top: midplane angular velocity relative toKeplerian in the ZAMO frame, (cid:104) V ϕ ( θ = π/ (cid:105) ρ / V kep , and mass density, (cid:104) ρ (cid:105) .Middle: radial velocity divided by the sound speed in the ZAMO frame, (cid:104)M r (cid:105) ρ , and the temperature, (cid:104) T (cid:105) ρ . Bottom: entropy, (cid:104) ( γ − − log( P /ρ γ ) (cid:105) ρ ,and plasma β , (cid:104) β − (cid:105) − ρ . Note that log here represents the natural logarithm.The accretion flow in these simulations is transonic, but the sonic radii at ∼ r g is much closer to the event horizon than for the a = ∼ r g . This is primarily due tothe dynamically important, β (cid:46) r (cid:46) r g .Torque exerted by these fields being dragged by the rotating Kerr spacetimealso leads to significant ( ∼ . V kep ) rotational velocities near the horizon r (cid:46) r g . perpendicular to that direction and thus the alignment in the upperjet may just be a coincidence. That conclusion is supported by thefact that the ψ = ◦ jet also shows a non-neglible tilt away from thespin axis in the two rightmost panels of Figure 9 despite the fieldhaving been initially aligned with the spin axis. It is likely that thedegree to which the jet is tilted is determined by stochastic sym-metry breaking as it propagates through a relatively dense mediumand wobbles about via kink modes.The jets are also intermittent in the sense that there are spatialgaps between di ff erent sections (e.g., third panel of Figure 11 andFigure 12) and their spatial extent periodically grows and dimin-ishes. The intermittency is caused by the periodic fluctuations ofthe fundamental power source of the jet, the magnetic flux, seen inthe top panel of Figure 8. Without this power source, the cavitiesinitially excised by the jet are quickly overcome by the surroundinggas as it falls inwards.Another interesting aspect of Figures 9–12 are the bubbles thatpinch o ff from the main body of the jet (e.g., the middle two pan-els of Figure 12 and rightmost panel of Figure 9). These are lowdensity, magnetic pressure supported regions of tangled magneticfield that form as the jets propagate outwards through the ambientmedium. They buoyantly rise through the gas until they dissipate atlarger radii.In order to further quantify the properties of the jet, we de-fine it precisely as all regions with b / > . ρ (as in Liska et al.2019b). This allows us to assign radii to the upper and lower jets, r + / − jet , as the maximum radii for which b / > . ρ and θ < π/ + ) or θ > π/ r + jet is plotted vs. time for the foursimulations in the top panel Figure 13 (the qualitative behavior of r − jet is similar). Initially, the jets expand outwards at a fraction ofthe speed of light but eventually stall around 300–800 r g . At theseradii the surrounding gas is roughly uniform and stationary sincethe Bondi radius is 200 r g . After stalling, the jet radius oscillates intime, correlated directly with fluctuations in φ BH and η in Figure8. The di ff erences between the time averaged jet radius in the foursimulations are as expected based on the di ff erences in outflow ef-ficiency (Figure 8). The ψ = ◦ jet stalls at the smallest radius, ∼ r g , while the ψ = ◦ jet stalls at the largest radius, ∼ r g . The ψ = ◦ and ψ = ◦ jets stall somewhere in between, ∼ r g .We find that velocities inside the jet are only mildly rela-tivistic, reaching maximum Lorentz factors γ = u t ( − g tt ) − / (cid:46) u r / u t (cid:46) . c ) at a few hundred r g and then rapidly decelerating inradius near the edge of the jet.We can also define a half opening angle of the jet, α jet as thesolution to the equation A jet ( r )2 π − α jet (cid:90) √− gd θ = , (11)where A jet ( r ) ≡ / (cid:33) H ( b / − . ρ ) √− gd θ d ϕ is the e ff ectivecross-sectional area of one jet (hence the factor of 1 /
2) and H isthe Heaviside step function. The time-averaged α jet is plotted vs.radius in the bottom panel of Figure 13 for our simulations. Nearthe event horizon, the jets are quite wide, α jet ∼ ◦ as shownpreviously in Figures 1–4. As material in the jet flows outwards itbecomes confined to a narrower and narrower angle, decreasing ina way that is consistent with a parabolic jet profile ( A jet ˜ ∝ r ). Forexample, at r = r g α jet ∼ ◦ while at r = r g α jet ∼ a fewdegrees. © , 000–000 S. M. Ressler, E. Quataert, C. J. White, O. Blaes = 0 tot (neg.)conv (pos.) adv (neg.) = 30 = 60 = 90 r ( r g ) r F E r ( r g ) Figure 7.
Flux of energy in our four di ff erent simulations, averaged over time and angle broken down into the total (solid), convective (dashed), and advective(dotted) components. Here we restrict the averages to | θ − π/ | < θ H , where θ H is the half-angle associated with the scale height of the disc. Convectiontransports outwards only a fraction ( (cid:46) /
3) of the energy being advected inwards except at the innermost radii ( r (cid:46) r g ). The convective flux is particularlysmall in the ψ = ◦ simulation, (cid:46) Figure 14 shows spacetime diagrams of the location of thejets (at r = r g ) projected onto the x – y plane. More precisely, wedefine x + jet = (cid:33) xH ( b / − . ρ ) H ( θ < π/ √− gd θ d ϕ (cid:33) H ( b / − . ρ ) H ( θ < π/ √− gd θ d ϕ (12)and similarly for y + jet . The jets wobble around with time in spiral pat-terns characteristic of the external kink instability, which is knownto occur when a highly magnetized jet attempts to penetrate a densemedium (e.g., Bromberg & Tchekhovskoy 2016; Tchekhovskoy &Bromberg 2016; Barniol Duran, Tchekhovskoy & Giannios 2017)and is generally caused by toroidally dominant fields in collimatedjets (Begelman 1998; Lyubarskii 1999). As the jets in our simula-tions push against the surrounding gas, they decelerate, convertingsome of their poloidal magnetic field to toroidal field like a springbeing compressed, leading to the unstable configuration near the ra-dial edge of the jet. As a result, its electromagnetic energy is largelyconverted into internal and kinetic energy by heating up and accel-erating the surrounding gas. Figure 15 demonstrates this for the ψ = ◦ simulation by plotting the time averaged ˙ E vs. radius, bro-ken up into electromagnetic, kinetic, and thermal components. Thethermal component starts o ff as negligible but steadily increasesuntil it overtakes the electromagnetic component at ≈ r g (typi- cally ∼ a few 100 r g for the other three simulations with smaller jetpower). At this point the jet has blown out a low density, high tem-perature cavity of much larger spatial extent than its cross-section.This is shown in Figure 16, which plots the time averaged radialvelocity, temperature, β , and pressure in the x – z plane for the same ψ = ◦ simulation on an r ∼ r g scale. Compared to the sur-rounding medium, the cavity is a factor of ∼
10 times hotter and ∼
10 times less dense so that the pressures are comparable. Most ofthe gas in the cavity is unbound, with velocities exceeding free-fall.Also plotted in the top left panel of Figure 16 are velocitystreamlines. Material on the edge of the jet / cavity with ˙ r (cid:46) free-falltends to circulate back towards the midplane and feed the inflow.In practice, while this does not significantly alter the hydrodynamicproperties of the accretion flow near the Bondi radius, it can notice-ably change the field geometry there (see, e.g., the rightmost panelof Figure 11 where the field lines in the ψ = ◦ simulation bearalmost no resemblance to the initial conditions and look more likethe ψ = ◦ simulation). Because of this e ff ect we are less concernedthat the idealized, purely laminar magnetic fields lines in the initialconditions are artificially reducing turbulence in the midplane. © , 000–000 t [ M ] B H = 0= 30 = 60= 900 2500 5000 7500 10000 12500 15000 17500 20000 t [ M ] | M / M B o n d i | = 0= 30 = 60= 900 2500 5000 7500 10000 12500 15000 17500 20000 t [ M ] [ % ] = 0= 30= 60= 90 Figure 8.
Angle integrated quantities plotted vs. time for the four simula-tions. Top: Dimensionless measure of the amount of magnetic flux thread-ing the event horizon, φ BH . Middle: Net mass accretion rate through thehorizon normalized to the Bondi rate, | ˙ M / ˙ M Bondi | . Bottom: Outflow e ffi -ciency measured at r = r g , η ≡ | ˙ E − ˙ M | / | ˙ M | . The ψ = ◦ and ψ = ◦ showmagnetically arrested behavior, with flux dissipation events (i.e., rapid de-creases in φ BH ) followed by large spikes in accretion rates. The other twosimulations are less time variable in both magnetic flux and accretion rate,saturating at moderately smaller values of φ BH . η is dominated by the elec-tromagnetic energy of the jet, and reaches as high as 200–300% (meaningthat energy is extracted from the rotation of the black hole) for ψ = ◦ . Theother three simulations have more modest e ffi ciencies, ∼ ψ = ◦ simulation with no net initial verticalflux, though only for 5000 M (cid:46) M when a jet is present (see Figure4). At other times η is < ψ = ◦ OUTLIER
Naively, one would assume that the amount of magnetic fluxthreading the black hole (and thus jet power) would either decreasewith increasing ψ (for 0 ◦ (cid:54) ψ (cid:54) ◦ ) because of the decreasedsupply of initial vertical magnetic field or be relatively independentof ψ because φ BH would tend to saturate at some unique value de-termined by the balance between the outwards Lorentz force andinflow. Instead, we find that the ψ = ◦ simulation has the highestvalue of φ BH . To determine whether this result is robust or simplycaused by stochastic variation between the four simulations, we runan additional six simulations with ψ = ◦ , ◦ , . ◦ , ◦ , . ◦ ,and 75 ◦ with exactly half the resolution in each direction as themain simulations of this paper. The resulting φ BH as function of ψ is plotted in Figure 17 (including all 10 simulations), where φ BH isaveraged over intervals containing several limit cycles once the fluxhas roughly saturated. With the range in ψ now finely sampled, itis clear that simulations with 40 ◦ (cid:46) ψ (cid:46) ◦ have systematicallyhigher values of φ BH than the rest ( ∼ ∼ φ BH on ψ shown in Figure17, first consider the case of a non-spinning black hole. In the limitof large Bondi radius, the initially uniform, weak magnetic field be-come almost entirely radial by the time it reaches the event horizon,flipping sign across a current sheet in the plane perpendicular tothe initial field and containing the origin, as illustrated in Figure 18which plots b and magnetic field lines at an early time in our sim-ulations (at this early time the e ff ects of a (cid:44) | θ − π/ | < α mid that will only contain theentire angular extent of the initial current sheet if ψ < ( π/ − α mid )(assuming ψ < ◦ ). If the jet forms before the onset of turbulencethen the inflow in that case would plausibly be characterized bymore reconnection (and hence, more turbulence and more dissipa-tion of the magnetic field) than the case when the initial currentsheet is only partially contained in the inflow region (the rest beingblown away by the jet). Less dissipation in the more laminar flowswould allow for larger amount of magnetic flux to accumulate anda higher saturation value for φ BH . This argument breaks down ifthe field reconnects and starts driving turbulence before the jet hasa chance to form, which happens when ψ nearly approaches 90 ◦ since it takes longer for vertical field to accumulate via advection.Thus there is a ‘sweet spot’ around π/ − α mid < ψ < ◦ where wewould expect φ BH to be largest. Since we find ( π/ − α mid ) (cid:38) ◦ near the horizon (bottom panel of Figure 13), this would explainwhy the ψ = ◦ simulation has the highest value of φ BH . For most simulations we use the interval 5000 M –20000 M , but for someof the lower resolution simulations the flux takes longer to saturate and sowe use intervals with later start times. For example, for ψ = ◦ we use16000 M –25000 M © , 000–000 S. M. Ressler, E. Quataert, C. J. White, O. Blaes
Figure 9.
Azimuthal slice of the magnetization parameter b /ρ at four di ff erent times in our ψ = ◦ simulation. Starting from the top left panel and proceedingto the right, the panels represent 2000, 8000, 14000, and 20000 M . A highly magnetized, well-collimated jet is visible at all times, though its size andorientation are significantly time-variable. It is also highly asymmetric about both the midplane and the polar axes, with the upper ( z >
0) and lower ( z < ff erent behavior. Each can be tilted by as much as 30 ◦ from the black hole spin axis at certain times. An animatedversion of this figure is available at https: // youtu.be / c6tVZ5qyumA. In this section we speculate how the results of our simulations de-pend on certain input parameters, namely, the Bondi radius, the ini-tial β , and the maximum magnetization, σ max . Although it is beyondthe scope of this work to do a comprehensive parameter survey, wecan make reasonable extrapolations based on physical argumentsand the results of past work in the literature. Even so, this is nosubstitute for systematically varying these parameters in our simu-lations and thus the conclusions we draw in this section should betreated as conjecture. Our choice of r B = r g is much less than appropriate for re-alistic low luminosity AGN such as Sgr A* and M87. In practicethis limits the region of reconnection-driven turbulence in our sim-ulations to a relatively small radial range. Since an initially uni-form, weak magnetic field will grow in strength due to the com-pression caused by spherical inflow as ( r / r B ) − (dominated by theradial component), the magnetic pressure will grow as ( r / r B ) − while the thermal pressure will grow as ∼ r − / (assuming a Bonditype flow). This means that β ∼ r conv ≈ r B β / where the field will start reconnecting. This ex-pression gives r conv ≈ r g for our parameters, roughly agreeingwith the simulations (Figure 6). For much larger r B , r con v will alsobe much larger and so we expect that convection will be presentacross several orders of magnitude in radius. In this subsection, forthe purposes of extrapolation we assume that the basic power lawdependences of the flow properties that we observe in these sim-ulations where convection / turbulence occurs in a relatively narrowradial range will hold for simulations where convection / turbulenceoccurs in a more extended radial range. This assumption requiresthat a significant amount of net magnetic flux will survive transportthrough the larger convective regions and still reach the black hole.Because the competition between advection and di ff usion of mag-netic fields in such a regime is not well understood, it is not obvioushow well this assumption is justified. The crux of the issue is thatturbulent di ff usion could, in principle, decouple the motion of themagnetic flux from the mass flux so that infall does not necessarilyguarantee an accumulation of φ BH onto the black hole. © , 000–000 Figure 10.
Same as Figure 9 but for ψ = ◦ . φ BH vs. ψ Dependence
Also uncertain is whether or not the φ BH vs. ψ dependence shownin Figure 17 will hold for larger convective regions. Our hypothesisfor why the ψ = ◦ simulation has the largest φ BH (§4) is that alarge fraction of the initial current sheet is close enough to the polarregions to be partially blown away by the jet. If the onset of con-vection / turbulence begins at much larger radii this subtlety of themagnetic field geometry may be washed out by the time it reachesthe event horizon even if the overall net direction is preserved. An-other possibility is that the jet could stall at smaller radii than r conv and have less of an influence on flux transport / reconnection. Theseconcerns prevent us from determining whether the dependence of φ BH on ψ seen in our models is a robust feature of spherical infall ofinitially coherent magnetic fields or unique to relatively small r B .Future work varying r B could shed more light on this issue. Since we expect that the radial velocity near the event horizon toalways be ∝ free fall ∼ c , the ρ ˜ ∝ ρ ( r B / r ) observed in our sim-ulations implies ˙ M ∝ ρ r B . Since ˙ M Bondi ∝ ρ r / , we predictthat ˙ M / ˙ M Bondi ∝ r − / . For Sgr A*, where r B ≈ × r g and˙ M Bondi ≈ − ˙ M (cid:12) / yr (Bagano ff et al. 2003; Wang et al. 2013), thisextrapolation would imply that ˙ M ∼ . × − –6 × − ˙ M Bondi ∼ × − ˙ M (cid:12) / yr, a reasonable accretion rate given observa- tional constraints (Marrone et al. 2006) and previous simulation-based estimates (Shcherbakov & Bagano ff r B ≈ × r g and ˙ M Bondi ≈ M (cid:12) / yr (Di Matteo et al.2003), we extrapolate that ˙ M ∼ × − –4 × − ˙ M Bondi ∼ × − ˙ M (cid:12) / yr, a number in good agreement with the mean accre-tion rates required for GRMHD MAD models to reproduce the 230GHz flux Event Horizon Telescope Collaboration et al. (2019a) andconsistent with Faraday rotation estimates (Kuo et al. 2014, thoughsee Ricarte et al. 2020 for a discussion of why the rotation measuremight not be a reliable tracer for the accretion rate).Physically, the reason that this extrapolation predicts muchsmaller accretion rates than the Bondi value is not primarily dueto mass outflows (though these are present near the jet boundary)but the reduced radial velocity caused by the presence of magneticfields. The Bondi rate is appropriate for a spherically symmetricdistribution of gas under the influence of gravity alone and can thusbe thought of as an e ff ective upper limit on the accretion rate fora system that includes other forces (e.g., magnetic forces), irre-versible dissipation (e.g., shocks or turbulence), or non-ideal e ff ects(e.g., conduction). © , 000–000 S. M. Ressler, E. Quataert, C. J. White, O. Blaes
Figure 11.
Same as Figure 9 but for ψ = ◦ . The radial extent of the jet is regulated by the kink instability, whichcauses the outer parts of the jet to spiral around the axis (Figure14) and generally become disrupted (Figures 9–12). This motionserves to dissipate the magnetic energy of the jet, heating the sur-rounding gas and propelling it into a wider, uncollimated outflow.The transition happens when the timescale for kink mods to grow, t kink , is su ffi ciently small compared to the time it takes for a fluidelement to traverse the remaining height of the jet, t dyn . Using ananalytic model of a cylindrical, relativistic jet propagating throughan external medium motivated by their simulations, Bromberg &Tchekhovskoy (2016, see also Bromberg et al. 2011) estimate that t kink t dyn ˜ ∝ L jet ρ a r γ / , (13)where L jet is the luminosity of the jet, ρ a ( r ) is the density of theambient medium, and γ jet is the Lorentz factor of the jet. First con-sider r < r B , where ρ a ˜ ∝ r B / r . Assuming L jet ∝ ˙ E ∝ ˙ Mc ∝ r B ,Equation (13) predicts that the stability properties of the jet will beindependent of r B in in this region (as long as γ jet is roughly in-dependent of r B ) and t kink / t dyn ∝ ( r / r g ) − / . On the other hand, for r > r B , where the jets in our simulations stall, ρ a ≈ ρ and thus t kink / t dyn ∝ [ r / ( r B r g )] − / .In order to estimate where a jet would ultimately stall if r B was (cid:29) r g , we can thus solve t kink t dyn (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ρ a = ρ ( r B / r ) (cid:16) r = r jet , r B (cid:29) r g (cid:17) = t kink t dyn (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ρ a = ρ (cid:16) r = r g , r B = r g (cid:17) (14)for r jet , where we have used 400 r g as a characteristic jet stallingradius for our simulations (Figure 13) and approximated the densityprofile as a piecewise power law. Substituting Equation (13) for t kink / t dyn into Equation (14) we find r jet ≈ r g independent of r B ,which is valid as long as r B (cid:38) r g . For r B (cid:46) r g , we replace ρ a = ρ ( r B / r ) with ρ a = ρ in the right hand side of Equation(14) and obtain r jet ≈ r g (cid:112) r B / r g . In summary, we (roughly)predict r jet ≈ r g (cid:113) r B r g r B (cid:46) r g r g r B (cid:38) r g . (15)In Sgr A* and M87, r B (cid:29) r g , so we expect that the jetwould stall somewhere around ∼ r g , ± a few hundred r g . Thisassumes, of course, that the linear stability criterion (Equation 13)correctly predicts the scaling of the non-linear disruption of the jetby the kink instability as r B is varied. © , 000–000 Figure 12.
Same as Figure 9 but for ψ = ◦ . β All simulations that we have presented use an initial β = ffi ciently large ( (cid:38)
10 or so) so that the initial magneticfield is not dynamically important. This is because no matter howsmall the field is at t =
0, the roughly spherical accretion of frozen-in field lines will lead to a configuration of nearly radial magneticfield lines near the horizon with a dependence on radius of b ∝ r − and a magnitude that is growing with time (Bisnovatyi–Kogan &Ruzmaikin 1974). The field strength will continue to grow untilit becomes dynamically important, at which point the reconnectionof oppositely directed field lines heats the gas and drives turbulence(e.g., Shvartsman 1971; Meszaros 1975). In non-relativistic simu-lations, Pang et al. (2011) found that their results were relatively in-sensitive to β over the range 10–1000, with the density power lawindex p having only a weak dependence on this parameter, β − . . σ max As described in §2, we use a density floor to impose σ (cid:54) σ max = ff ect on the mainbody of the accretion flow (where σ is generally (cid:28) ff ectivelysets the value of σ at the base of the jets in our simulations. In re-ality σ in these regions would be larger, though its precise valueis unknown and may be set by pair-production processes (Bland- ford & Znajek 1977; Chen, Yuan & Yang 2018). A larger σ in thebase of the jet (i.e., a lower mass density) would lead to a largerLorentz factor throughout the main body of the jet, which coulda ff ect the kink stability criterion and, consequently, r jet . Equation(13) predicts, however, that the dependence of t kink / t dyn on γ jet isrelatively weak, ∝ γ − / , and the dependence on σ is probably evenweaker. For example, Chatterjee et al. (2019) found that adoptingvalue for σ max of 3, 10, 50, and 100 resulted in maximum Lorentzfactors of ≈
2, 5, 8, and 10, while McKinney (2006) found thata value of σ max = resulted in a maximum Lorentz factor of ≈
10. While both of these works assumed axisymmetry (thus notable to capture kink instabilities) and empty polar regions, theseassumptions would, if anything, result in a higher conversion e ffi -ciency from magnetic to kinetic energy and a stronger dependenceof γ jet on σ max than in our 3D, kink unstable jets that have to con-tinually fight against the ambient medium. This is encouraging butby no means conclusive evidence that r jet and the evolution of thekink instability in our simulations are not strongly sensitive to thearbitrary choice of σ max . © , 000–000 S. M. Ressler, E. Quataert, C. J. White, O. Blaes t [ M ] r + j e t [ r g ] = 0= 30= 60= 9010 r [ r g ] j e t [ d e g r ee s ] = 0= 30= 60= 90 Figure 13.
Upper jet radius, r + jet , as a function of time (top) and jet half-opening angle, α jet , as a function of radius (bottom) in our four simulations.Jets initially travel outwards at a constant speed that is a fraction of the speedof light but ultimately stall at ∼ r g , at which point they periodicallyfall back and then re-expand. The opening angles of the jets decrease withradius in a manner consistent with a parabolic shape.
150 100 50 0 50 100 150 x jet ( r = 200 r g ) [ r g ] y j e t ( r = r g ) [ r g ] = 0= 30 = 60= 90 Figure 14.
Space-time diagrams of the location of the jet at r = r g pro-jected onto the x − y plane in the four simulations. The kink instability causesthe jets to move around in spiral-like patterns as a function of time. Withthe largest jet power (Figure 8), the ψ = ◦ simulation also has smallestamplitude variations in jet location. r [ r g ] E Tot.EM KETh.
Figure 15.
Radial energy flux, ˙ E , broken up into total (solid), electromag-netic (dashed), kinetic (dotted), and thermal (dot-dashed) components forthe ψ = ◦ simulation. At small radii, the energy outflow is dominated bythe Poynting dominated jet. As the jet propagates outwards, its magneticenergy dissipates into kinetic and thermal energy. The latter dominates for r (cid:38) r g . The other three simulations show qualitatively similar behaviorand are thus not shown. The reconnection-driven convection present in our simulations is afeature of both convection-dominated Bondi flow (CDBF) models(Igumenshchev & Narayan 2002) and magnetically frustrated mod-els (Pen, Matzner & Wong 2003; Pang et al. 2011). The fundamen-tal di ff erence between these two classes of solutions is the amountof energy transported by convection, F conv . In CDBF models F conv is large and positive, dominating the total energy transport budgetand thus requiring ρ ∝ r − / in steady state. In the magneticallyfrustrated model, on the other hand, the magnitude of F conv is rel-atively small and can even be negative (that is, carrying energy in-wards). This allows for the power law index, p , of the radial densityprofile to be anywhere between − / − / F conv is suppressed because the magnetic field dampsthe convective motions, “frustrating” its ability to transport energy.In a survey over initial magnetic field strength and Bondi radius,Pang et al. (2011) found p ≈ − p ≈ − p = − / ψ = ◦ simula-tion looks quite similar to those in the Pen, Matzner & Wong (2003,their Figure 1) and Pang et al. (2011, their Figure 1) simulations,with turbulent motions creating fairly random structure with no pre-ferred direction. At other times in the same simulation (e.g., bottomright panel of Figure 4) and at all times (after the initial transientperiod) in the other simulations, however, the jets powered by thespinning black hole evacuate the polar regions of matter and local-ize the turbulence to the midplane. These black hole spin-poweredjets have no analog in non-relativistic MHD.While we do see a reduced inwards radial velocity comparedto the hydrodynamic case (as evidenced by the critical sonic pointmoving inwards from ∼ r g to ∼ r g (Figure 6), we do not find © , 000–000 Figure 16.
Time-averaged, y = ψ = ◦ simulation at large radii, where the jet breaks up due to the kink instability. Top left:magnitude of the radial velocity normalized to the free-fall speed, | ˙ r | / v ff , overplotted with velocity streamlines. Top right: gas temperature, T g . Bottom left:plasma β . Bottom right: pressure normalized to the initial pressure, P / P . The “wobbling” of the jet (Figure 14) caused by the kink instability carves out acavity of hot, magnetized gas with large radial velocities. This cavity is in rough pressure equilibrium with its surroundings. Gas near the edge of the jet / cavitycirculates back into the inflow, altering the supply of magnetic flux in a nonlinear way. The other three simulations show qualitatively similar behavior and arethus not shown. the ∼ zero radial velocities characteristic of the nearly hydrostaticCDBF and magnetically frustrated solutions. This is a strictly GRe ff ect, as all flows must pass through a sonic point before reach-ing the event horizon. The stronger outwards force provided by acombination of magnetic fields and pressure can only delay the in-evitable by moving the sonic point towards smaller radii. If we wereto extend our simulations to larger r B , the convective / turbulent zonewould also extend to larger radii and could plausibly reach a nearlyhydrostatic state in the region where GR e ff ects are negligible. / Jet Simulations
Two of our simulations ( ψ = ◦ and ψ = ◦ ) appear to be mag-netically arrested based on the traditional build-up and dissipationcycle seen in φ BH and the associated fluctuations in accretion ratecaused by the brief repulsion of matter from near the black hole.The other two simulations ( ψ = ◦ and ψ = ◦ ) do not appear tobe fully arrested since φ BH and ˙ M are much more stable with time,yet they are likely just under the threshold in magnetic flux to reachthe arrested state. Given that conclusion, it is instructive to compare our results to previous GRMHD MAD simulations (e.g., Igumen-shchev, Narayan & Abramowicz 2003; Tchekhovskoy, Narayan &McKinney 2011; Narayan et al. 2012; White, Stone & Quataert2019; Chael, Narayan & Johnson 2019). To first approximation,Figures 1–4, which show 2D slices of mass density and magneticfield lines, look remarkably similar to GRMHD MAD simulationsat most times. Both have evacuated, strongly magnetized polar re-gions with a large amount of coherent / vertical magnetic field lines.Both see accretion via thin, non-axisymmetric streams due to themagnetic field strongly compressing the inflow close to the horizon(e.g., Igumenshchev, Narayan & Abramowicz 2003 ; see also, ourFigure 5 which shows midplane slices of β for ψ = ◦ ). Thesesimilarities between torus-based MADs and our initially zero an-gular momentum simulations are, perhaps, not that surprising be-cause MAD flows are generally more sub-Keplerian than SANEflows (e.g., Narayan et al. 2012) and the magnetic fields draggedby the a = . © , 000–000 S. M. Ressler, E. Quataert, C. J. White, O. Blaes [degrees] B H high res.low res. Figure 17.
Dimensionless magnetic flux threading the black hole averagedover time (see main text), φ BH , vs. initial magnetic tilt angle ψ . Circlesrepresent simulations that employ the fiducial resolution described in §2,squares represent simulations that employ a resolution reduced by 2 in allthree directions, while error-bars represent the standard deviation of φ BH over the chosen time interval. There is a clear peak at ψ ∼ ◦ –80 ◦ that wediscuss in §4. gions are essentially vacuum. This allows for them to reach muchlarger distances and Lorentz factors than we see in our simulations(e.g., McKinney 2006; Chatterjee et al. 2019). In this respect, oursimulations are more similar to those performed by Bromberg &Tchekhovskoy (2016); Tchekhovskoy & Bromberg (2016); BarniolDuran, Tchekhovskoy & Giannios (2017), where the jets were em-bedded in spherical interstellar mediums of various density pro-files ρ − s . Bromberg & Tchekhovskoy (2016) provided an analyticalmodel for the jet as a function of height and found that for s < s ≈ r (cid:46) r Bondi and s ≈ r (cid:38) r Bondi ) do both of those things (Figure 13 and 9–12).Note that the value of s in this work for r (cid:46) r B is not a parameterbut is determined by the dynamics of the accretion flow.Simulations of tilted disks find that the jets at large radii tendto align perpendicular to the large-scale disc (e.g., White, Quataert& Blaes 2019; Liska et al. 2020). This is primarily caused by the jetpropagating through the pre-evacuated, low-density funnel carvedout by the initial torus. This region is the path of least resistancefor the Poynting flux to travel outwards (as opposed to penetrat-ing through the outer disc). Since our simulations start with a zeroangular momentum, spherical distribution of gas, the jets have noother option but to penetrate through the ambient medium. Thisprevents them from getting systematically channelled away fromthe black hole spin axis. We now compare our results to the simulations that inspired themin the first place, namely, the GRMHD wind-fed Sgr A* accre-tion simulations of Ressler et al. (2020), hereafter R20. Both setsof simulations display magnetically confined midplanes, evacuatedpolar regions, magnetically arrested behavior, and relatively low β ,coherent magnetic fields. Both also find ρ ˜ ∝ r − , though it is notclear that it is for the same reasons. In Ressler et al. (2020) the density profile was caused by the underlying distribution of angu-lar momentum with accretion rate in material being fed from largeradii, while here it is a result of ‘magnetically frustrated’ convec-tion in a pressure supported configuration. Why exactly the lat-ter results in an r − density profile is an open question. Gruzinov(2013) argued that if there was some physical process enforcinga constant momentum flux (i.e., ˙ P ∝ ρ v r r = const. in the non-relativistic regime) then ρ ∝ r − follows as a natural consequence.˙ P = const. can be achieved in a steady state only if gravity, pres-sure, and magnetic forces are perfectly balanced. While this doesseem to be the case in the non-relativistic simulations (e.g., Figure2 in Pen, Matzner & Wong 2003), in our simulations ˙ P increaseswith decreasing radius. The true reason for the r − density scalingmay just be the less satisfying observation that the simulations fallsomewhere between Bondi-type flows and Convection DominatedAccretion Flow (CDAF)-type flows so the density power law indexmust fall somewhere in between − / − /
2. It is also not clearif our simulations have su ffi ciently large dynamic range in radiusto determine the truly self-similar scaling of the density profile.The most obvious di ff erence between the accretion flow in thiswork and that found in R20 is the nature of the polar regions. SinceR20 used a =
0, the poles were evacuated due to a high concentra-tion of vertical magnetic fields but lacked the toroidal field strengthnecessary to drive powerful jets. As a result, the flow became quasi-spherical past a few 10s of r g . Furthermore, these magnetic poleschanged direction with time (corresponding to variations in the netfield being fed from larger radii), also tilting the gas away fromthe previous midplane. As we have shown, because the black holein our simulations is rapidly spinning, it strongly winds up thepoloidal magnetic field in the toroidal direction, producing jets thatreach to ∼ hundreds of r g . While the jets at large radii can be signif-icantly tilted with respect to the black hole spin due to propagatione ff ects and the kink instability, near the horizon the they are alwaysaligned with the spin axis regardless of initial magnetic field orien-tation. We expect that if the R20 simulations were run with a largeblack hole spin, they would undoubtably feature a jet that woulddisrupt the nature of the flow out to a relatively large radius andwould most likely be less prone to oscillations in the location ofthe midplane. As this work was in the latter stages of preparation, Kelly et al.(2020, hereafter K20) appeared on arXiv presenting the results of aparallel study of low angular momentum GRMHD accretion ontorotating black holes in the context of black hole merger remnants,focussing mainly on the orientation and power of the resulting jet.The initial conditions used in their simulations were essentially thesame as those used here: uniform pressure, density, and magneticfield of varying inclination angles with respect to the black holespin. The main di ff erences were that they 1) assumed a black holespin of a = .
69, lower than our a = . r B ≈ r g compared toour r B = r g , and 3) an adiabatic index appropriate for radiationpressure-supported gas, γ = /
3. Their fiducial model had r B ≈ r g In contrast to our results, K20 found that the jets in their sim-ulations only aligned with the spin axis of the black hole for r (cid:46) r g , while at larger radii the jets aligned with the initial magneticfield direction. Moreover, they found jet power to monotonicallydecrease with ψ and found no evidence of the ψ ∼ ◦ peak in η as we do (see bottom panel of Figure 8, Figure 17, and §4) seen © , 000–000 Figure 18.
2D slice in the x - z plane of b overplotted with magnetic field lines at an early time ( t = M ) in our four simulations. Spherical advection of themagnetic field causes a current sheet of b ≈ ψ = ◦ and ψ = ◦ ), it overlaps with the polar regions that are ultimately evacuated by jets. This meansthat if the jet forms before the first reconnection event (as it does for ψ = ◦ but not ψ = ◦ ), there can be a reduced amount of turbulence and an increase inthe amount of magnetic flux that reaches the event horizon. in our simulations. We believe that this is predominantly becausethe K20 simulations are in an entirely di ff erent parameter regimethan ours, where the dynamical range in radius is very small andany extrapolation in Bondi radius breaks down. This can be seenfrom their plots of various quantities vs. ‘ κ ’, their proxy for the en-tropy / temperature of the initial gas which relates directly to r B . For r B (cid:46) r g , their simulations are roughly independent of r B , whilefor r B (cid:38) r g they show an increase of | ˙ M | , ˙ E , and η with increasing r B . Moreover the simulations with smaller r B are much less variablein time than those with larger r B . Thus, there is a clear transition at r B ∼ r g ; the regime of r B (cid:46) r g is appropriate for some mergerscenarios (and perhaps other stellar mass sized black holes), whilethe regime of r B (cid:38) r g is appropriate for active galactic nuclei. Since K20 focused on the former while we focus on the latter, ourworks are complementary. We have presented GRMHD simulations of spherical accretiononto a rotating black hole modified by the presence of uniform,relatively weak magnetic fields of various tilt angles with respectto the black hole spin. As the field is advected inwards it forms cur-rent sheets and reconnects, dissipating magnetic energy and drivingturbulence (Figures 1–4). Similar to previous non-relativistic ‘mag-netically frustrated’ models, we find that convection does not e ffi -ciently transport energy outwards (Figure 7) and the mass densityscales as ˜ ∝ r − (Figure 6). Our simulations are either magnetically © , 000–000 S. M. Ressler, E. Quataert, C. J. White, O. Blaes arrested or very close to magnetically arrested, with the dimension-less magnetic flux, φ BH saturating at 20–80, often showing limitcycles of slow build-up followed by rapid dissipation (Figure 8).Torque caused by frame dragging of magnetic field lines threadingthe rapidly spinning black hole is enough to result in significant(up to ∼ ψ (the initial magnetic field direction), butrather results from symmetry breaking as the turbulently fed jetplows through the surrounding gas (Figures 9–12). At larger radiithe jet dynamics become dominated by the kink instability (Figure14), ultimately becoming disrupted at a few hundred r g (Figure 13)at which point they transition into an uncollimated outflow of hotgas reaching to ∼ r g (Figure 16). The dissipation mechanismof the kink modes is likely magnetic reconnection near the jet wallboundary (Bromberg et al. 2019), which could in turn accelerateparticles to high energy (Davelaar et al. 2020) and contribute tononthermal emission.The nature of our simulations is highly stochastic. Not onlyare the dynamics of the accretion flow governed by reconnection-driven turbulence, the dynamics and shape of the jet are governedby the kink instability and interaction with the turbulent ambientmedium. This is evident by the various contours of density, β , and b /ρ plotted in Figures 1–4, 5, and 9–12 which vary dramaticallyas a function of time and frequently display large global asymme-tries. Small changes in the initial conditions such as the tilt of themagnetic field ( ψ ) could therefore lead to completely di ff erent real-izations of the flow, accounting for much of the variation betweenthe simulations we have presented. All things considered, the dy-namics of the flow do not depend strongly on ψ except for a coupleof special cases.The first case is ψ = ◦ where there is exactly zero net ini-tial vertical flux (and presumably for ψ close to 90 ◦ ): the jet takeslonger to form, is less e ffi cient than the other simulations, and is in-termittent even at horizon scales; all because vertical flux can onlybe generated through turbulent motions. The second case is ψ ∼ ◦ –80 ◦ where there is a peak in the magnetic flux threading thehorizon, φ BH , (Figure 17) that is a factor of (cid:46) ψ (cid:46) ◦ and ψ (cid:38) ◦ . We suspect that this range inangles represents a “sweet spot” where there is slightly less turbu-lence inhibiting the accretion of magnetic flux because the large-scale current sheet is partially blown away by the polar outflow andyet there still enough initial vertical flux to form a jet before theonset of reconnection (see §4).The fact that even the ψ = ◦ simulation forms (at times) arelatively high e ffi ciency jet with η (cid:46)
20% (third panel of Figure 8)gives further evidence that a large-scale net vertical flux is not nec-essarily a prerequisite for Blandford & Znajek (1977) jet formation.This was previously shown for flows with larger angular momen-tum where an MRI-driven dynamo can convert toroidal magneticfield into poloidal magnetic field (Parfrey, Giannios & Beloborodov2015; Christie et al. 2019; Liska, Tchekhovskoy & Quataert 2020).Here we find that turbulent motions alone can produce enough lo-cal net vertical magnetic flux to power a jet from an already mostlypoloidal but zero net vertical flux field. These results suggest thatjets may be ubiquitous in systems with rapidly rotating black holes regardless of the magnetic field geometry in the accretion flow’ssource.The simulations presented in this work may be a reasonablemodel for Sgr A*. Magnetically arrested-type flows have been fa-vored in several recent analyses comparing models to observa-tions due to 1) the relatively high levels of linear polarization ob-served in the mm emission (models with less magnetic flux tendto require higher accretion rates that can depolarize the emission,Jim´enez-Rosales & Dexter 2018; Dexter et al. 2020a) and 2) therelatively large amount of vertical magnetic field required to re-produce the NIR flare polarization periodicity (Gravity Collabora-tion et al. 2018; Dexter et al. 2020b; Porth et al. 2020; GRAVITYCollaboration et al. 2020). Moreover, wind-fed GRMHD accretionsimulations that take into account the WR stellar winds primarilysourcing the accretion flow naturally evolve into an arrested state(Ressler et al. 2020). On the other hand, magnetically arrested mod-els are known to produce powerful jets (at least for rapidly spinningblack holes, e.g., Tchekhovskoy, Narayan & McKinney 2011) andyet there is no unambiguous detection of a collimated jet in SgrA*. Our simulations suggest, however, that this fact does not nec-essarily rule out the arrested model. Feeding by stellar winds canresult in a quasi-spherical distribution of gas outside of the innerhorizon-scale that pressure-confines the jet prompting dissipationthrough the kink instability. If this happens well inside the Bondiradius then larger scale observations may only be able to probe thehot, matter-dominated outflow blown out by the unstable jet ‘head’as it dissipates (e..g, Figure 16). In fact, recent work using ALMAmay have detected just that, finding evidence for a high-velocity,bipolar outflow with an opening angle of ∼ ◦ extending from ∼ × r g to ∼ r g (Royster et al. 2019; Yusef-Zadeh et al. 2020).Alternatively, Sgr A* could be slowly spinning or the jet may bepointed directly at us. Finally, the r − dependence of the mass den-sity that we find in our simulations (Figure 6) is consistent withobservations that use di ff erent methods to probe multiple di ff erentscales of the gas surrounding Sgr A* (Xu et al. 2006; Gillessenet al. 2019).It is less clear whether our simulations represent a viablemodel for M87 or other similar galactic nuclei. Resolved X-rayemission at the Bondi radius of M87 can be fit with an ρ ∝ r − distribution (Russell et al. 2015), as in our simulations, though itis unknown if this dependence persists to smaller radii. In order toreproduce the unambiguous, well-collimated jet observed out to ∼ –10 r g we would require the jets in our simulations to be sta-ble to kink instabilities out to radii (cid:38) times what we see here(Figure 13). Simply increasing the Bondi radius to a more realisticvalue may accomplish this, but based on an (admittedly uncertain)extrapolation we have argued that the radius at which the jet dis-rupts may be independent of Bondi radius (see §5.1.1 and the nextparagraph for a discussion). This could indicate that M87 is fed ina di ff erent manner than the Galactic Centre, with less low angularmomentum gas and more disk-like accretion. Indeed, even if theM87 black hole is predominantly fed by stellar winds (which isnot clear given the larger hot gas reservoir in the inner intraclustermedium of Virgo), the nature of this feeding must be di ff erent insome respects since the much larger (absolute) accretion rate re-quires a much larger number of winds to sustain. This is especiallytrue given that in Sgr A* only a handful of winds end up contribut-ing to the accretion rate (Loeb 2004; Cuadra, Nayakshin & Martins2008; Ressler, Quataert & Stone 2018).The self-consistent formation of the rising, buoyant bubbles(or cavities) in our simulations due to a combination of the kinkinstability and time variability of the jet (Figures 9–12) may be © , 000–000 broadly analogous to what occurs in galaxy clusters (e.g., Fabianet al. 2000; Blanton et al. 2010), though the spatial and tempo-ral scales we simulate are orders of magnitude smaller than thoseobserved in clusters. If the estimated scalings of our results with r B (§5.1.1) are correct, then even for much larger r B these cavitieswould never reach the sizes and distances required to match obser-vations. On the other hand, if the flow had more angular momentumso that the poles were pre-evacuated, the jet could potentially reachthe required distances before becoming kink unstable enough toproduce bubbles.For computational reasons we have focused on one particularBondi radius ( r B = r G ) that is orders of magnitude less thanthe actual Bondi radii for real systems, including both M87 andSgr A*. Using the dependence on radius of ρ and the kink stabilitycriteria, we extrapolated to larger r B (cid:29) r g (§5.1.1) to predictthat ˙ M / ˙ M Bondi ∝ r − / and r jet ∝ const, where r jet is the radius atwhich the jet disrupts by the kink instability. However, since wehave not rigorously tested these dependencies by running additionalsimulations with larger r B , they should be taken with a grain of salt.Future work can test these predictions using a parameter survey in r B or by using the actual Bondi radius in a multi-scale, multiplesimulation technique to extend the dynamic range in the manner ofR20.We have also focused on one particular value of σ max , the max-imum σ ≡ b /ρ allowed in the simulations, as well as one particularvalue of the initial β . It is unclear the degree to which our resultsare a ff ected by these choices. In §5.2 and §5.3 we argue that the de-pendence should be relatively weak based on previous simulationsin the literature, though none of those works had an entirely similarset-up as we have here. For β , the hope is that for any initial β (cid:29) β just taking longerto reach that state. For σ max , realistic astrophysical jets likely have σ max orders of magnitude larger than the σ max =
100 imposed here,which if the conversion of magnetic to kinetic energy in the jetswas 100% e ffi cient, this would lead to orders of magnitude largerLorentz factors. This is an important consideration since jets withlarger Lorentz factors could be even more unstable to kink modes(Equation 13). However, GRMHD jets tend to be rather ine ffi cientat converting magnetic to kinetic energy even in 2D where they arekink stable, which makes γ jet not as strongly dependent on σ max .We expect this to be even more true in the 3D, kink unstable jetswe have presented here and thus we can reasonably hope that themain conclusions drawn in this work are not strongly sensitive to σ max . ACKNOWLEDGMENTS
We thank S. Philippov for useful discussions, J. Stone and C. F.Gammie for comments on the manuscript, as well as all the mem-bers of the Horizon Collaboration, http: // horizon.astro.illinois.edu.SMR was supported by the Gordon and Betty Moore Founda-tion through Grant GBMF7392. SMR also thanks R. and D.Ressler for their generous hospitality during the writing of thismanuscript. This work was supported in part by NSF grants NSFPHY–1748958, AST–1715054, AST–1715277, a Simons Investi-gator award from the Simons Foundation, and by the NSF throughXSEDE computational time allocations TG–AST170012 and TG–AST200005 on Stampede2. This work was made possible by com-puting time granted by UCB on the Savio cluster. DATA AVAILABILITY
The data underlying this paper will be shared on reasonable requestto the corresponding author.
REFERENCES
Bagano ff F. K. et al., 2003, ApJ, 591, 891Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214Ball D., ¨Ozel F., Psaltis D., Chan C.-k., 2016, ApJ, 826, 77Barniol Duran R., Tchekhovskoy A., Giannios D., 2017, MNRAS,469, 4957Begelman M. C., 1998, ApJ, 493, 291Bicak J., Janis V., 1985, MNRAS, 212, 899Bisnovatyi–Kogan G. S., Ruzmaikin A. A., 1974, Ap&SS, 28, 45Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433Blanton E. L., Clarke T. E., Sarazin C. L., Randall S. W., Mc-Namara B. R., 2010, Proceedings of the National Academy ofScience, 107, 7174Boyer R. H., Lindquist R. W., 1967, Journal of MathematicalPhysics, 8, 265Bromberg O., Nakar E., Piran T., Sari R., 2011, ApJ, 740, 100Bromberg O., Singh C. B., Davelaar J., Philippov A. A., 2019,ApJ, 884, 39Bromberg O., Tchekhovskoy A., 2016, MNRAS, 456, 1739Calder´on D., Cuadra J., Schartmann M., Burkert A., Russell C.M. P., 2020, ApJL, 888, L2Chael A., Narayan R., Johnson M. D., 2019, MNRAS, 486, 2873Chael A., Rowan M., Narayan R., Johnson M., Sironi L., 2018,MNRAS, 478, 5209Chael A. A., Narayan R., Sad¸owski A., 2017, MNRAS, 470, 2367Chan C.-K., Psaltis D., ¨Ozel F., Narayan R., Sad¸owski A., 2015,ApJ, 799, 1Chandra M., Foucart F., Gammie C. F., 2015, ApJ, in prep.Chatterjee K., Liska M., Tchekhovskoy A., Marko ff S. B., 2019,MNRAS, 490, 2200Chatterjee K. et al., 2020, arXiv e-prints, arXiv:2002.08386Chen A. Y., Yuan Y., Yang H., 2018, ApJL, 863, L31Christie I. M., Lalakos A., Tchekhovskoy A., Fern´andez R., Fou-cart F., Quataert E., Kasen D., 2019, MNRAS, 490, 4811Cuadra J., Nayakshin S., Martins F., 2008, MNRAS, 383, 458Cuadra J., Nayakshin S., Springel V., Di Matteo T., 2005, MN-RAS, 360, L55Cuadra J., Nayakshin S., Springel V., Di Matteo T., 2006, MN-RAS, 366, 358Davelaar J., Philippov A. A., Bromberg O., Singh C. B., 2020,ApJL, 896, L31De Villiers J.-P., Hawley J. F., 2003, ApJ, 589, 458Dexter J. et al., 2020a, MNRASDexter J. et al., 2020b, MNRAS, 497, 4999Di Matteo T., Allen S. W., Fabian A. C., Wilson A. S., YoungA. J., 2003, ApJ, 582, 133Event Horizon Telescope Collaboration et al., 2019b, ApJL, 875,L1Event Horizon Telescope Collaboration et al., 2019a, ApJL, 875,L5Fabian A. C. et al., 2000, MNRAS, 318, L65Fishbone L. G., Moncrief V., 1976, ApJ, 207, 962Foucart F., Chandra M., Gammie C. F., Quataert E.,Tchekhovskoy A., 2017, MNRAS, 470, 2240Fragile P. C., Blaes O. M., Anninos P., Salmonson J. D., 2007,ApJ, 668, 417 © , 000–000 S. M. Ressler, E. Quataert, C. J. White, O. Blaes
Gammie C. F., McKinney J. C., T´oth G., 2003, ApJ, 589, 444Gillessen S. et al., 2019, ApJ, 871, 126Gravity Collaboration et al., 2018, A&A, 618, L10GRAVITY Collaboration et al., 2020, arXiv e-prints,arXiv:2009.01859Gruzinov A., 2001, arXiv e-prints, astroGruzinov A., 2013, arXiv e-prints, arXiv:1311.5813Hawley J. F., Smarr L. L., Wilson J. R., 1984, ApJ, 277, 296Igumenshchev I. V., Narayan R., 2002, ApJ, 566, 137Igumenshchev I. V., Narayan R., Abramowicz M. A., 2003, ApJ,592, 1042Jim´enez-Rosales A., Dexter J., 2018, MNRAS, 478, 1875Johnson B. M., Quataert E., 2007, ApJ, 660, 1273Kelly B. J., Eteinne Z. B., Golomb J., Schnittman J. D.,Baker J. G., Noble S. C., Ryan G., 2020, arXiv e-prints,arXiv:2010.11259Kerr R. P., 1963, Phys. Rev. Lett., 11, 237Kuo C. Y. et al., 2014, ApJL, 783, L33Liska M. et al., 2019a, arXiv e-prints, arXiv:1912.10192Liska M., Hesp C., Tchekhovskoy A., Ingram A., van der Klis M.,Marko ff S., 2018, MNRAS, 474, L81Liska M., Hesp C., Tchekhovskoy A., Ingram A., van der Klis M.,Marko ff S. B., Van Moer M., 2020, MNRASLiska M., Tchekhovskoy A., Ingram A., van der Klis M., 2019b,MNRAS, 487, 550Liska M., Tchekhovskoy A., Quataert E., 2020, MNRAS, 494,3656Loeb A., 2004, MNRAS, 350, 725Lyubarskii Y. E., 1999, MNRAS, 308, 1006Mao S. A., Dexter J., Quataert E., 2017, MNRAS, 466, 4307Marrone D. P., Moran J. M., Zhao J.-H., Rao R., 2006, ApJ, 640,308McKinney J. C., 2006, MNRAS, 368, 1561Meszaros P., 1975, A& A, 44, 59Michel F. C., 1972, Ap. Space Sci., 15, 153Mo´scibrodzka M., Falcke H., Shiokawa H., 2016, A&A, 586, A38Mo´scibrodzka M., Falcke H., Shiokawa H., Gammie C. F., 2014,A&A, 570, A7Mo´scibrodzka M., Gammie C. F., Dolence J. C., Shiokawa H.,Leung P. K., 2009, ApJ, 706, 497Narayan R., Igumenshchev I. V., Abramowicz M. A., 2000, ApJ,539, 798Narayan R., Igumenshchev I. V., Abramowicz M. A., 2003, PASJ,55, L69Narayan R., Sa¸dowski A., Penna R. F., Kulkarni A. K., 2012, MN-RAS, 426, 3241Narayan R., Yi I., 1995, ApJ, 444, 231¨Ozel F., Psaltis D., Narayan R., 2000, ApJ, 541, 234Pang B., Pen U.-L., Matzner C. D., Green S. R., Liebend¨orfer M.,2011, MNRAS, 415, 1228Parfrey K., Giannios D., Beloborodov A. M., 2015, MNRAS, 446,L61Pen U.-L., Matzner C. D., Wong S., 2003, ApJ, 596, L207Penna R. F., Kulkarni A., Narayan R., 2013, A& A, 559, A116Porth O. et al., 2019, ApJS, 243, 26Porth O., Mizuno Y., Younsi Z., Fromm C. M., 2020, arXiv e-prints, arXiv:2006.03658Porth O., Olivares H., Mizuno Y., Younsi Z., Rezzolla L., Mosci-brodzka M., Falcke H., Kramer M., 2017, Computational Astro-physics and Cosmology, 4, 1Quataert E., Gruzinov A., 2000, ApJ, 539, 809Ressler S. M., Quataert E., Stone J. M., 2018, MNRAS, 478, 3544 Ressler S. M., Quataert E., Stone J. M., 2020, MNRAS, 492, 3272Ressler S. M., Tchekhovskoy A., Quataert E., Chandra M., Gam-mie C. F., 2015, MNRAS, 454, 1848Ressler S. M., White C. J., Quataert E., Stone J. M., 2020, ApJL,896, L6Ricarte A., Prather B. S., Wong G. N., Narayan R., Gammie C.,Johnson M. D., 2020, MNRAS, 498, 5468Ripperda B. et al., 2019, ApJS, 244, 10Royster M. J., Yusef-Zadeh F., Wardle M., Kunneriath D., CottonW., Roberts D. A., 2019, ApJ, 872, 2Russell H. R., Fabian A. C., McNamara B. R., Broderick A. E.,2015, MNRAS, 451, 588Ryan B. R., Dolence J. C., Gammie C. F., 2015, ApJ, 807, 31Ryan B. R., Ressler S. M., Dolence J. C., Tchekhovskoy A., Gam-mie C., Quataert E., 2017, ApJL, 844, L24Sa¸dowski A., Wielgus M., Narayan R., Abarca D., McKinneyJ. C., Chael A., 2017, MNRAS, 466, 705Shcherbakov R. V., Bagano ff F. K., 2010, ApJ, 716, 504Shvartsman V. F., 1971, Soviet Atron., 15, 377Stone J. M., Tomida K., White C. J., Felker K. G., 2020 in press,ApJS, arXiv:2005.06651Tchekhovskoy A., Bromberg O., 2016, MNRAS, 461, L46Tchekhovskoy A., Narayan R., McKinney J. C., 2011, MNRAS,418, L79Wald R. M., 1974, Phys. Rev. D, 10, 1680Wang Q. D. et al., 2013, Science, 341, 981White C. J., Dexter J., Blaes O., Quataert E., 2020, ApJ, 894, 14White C. J., Quataert E., Blaes O., 2019, ApJ, 878, 51White C. J., Stone J. M., Gammie C. F., 2016, ApJS, 225, 22White C. J., Stone J. M., Quataert E., 2019, The AstrophysicalJournal, 874, 168Xu Y.-D., Narayan R., Quataert E., Yuan F., Bagano ff F. K., 2006,ApJ, 640, 319Yuan F., Quataert E., Narayan R., 2003, ApJ, 598, 301Yusef-Zadeh F., Royster M., Wardle M., Cotton W., KunneriathD., Heywood I., Michail J., 2020, MNRAS ©000