The Surprisingly Small Impact of Magnetic Fields On The Inner Accretion Flow of Sagittarius A* Fueled By Stellar Winds
MMon. Not. R. Astron. Soc. , 000–000 (0000) Printed 15 January 2020 (MN L A TEX style file v2.2)
The Surprisingly Small Impact of Magnetic Fields On The InnerAccretion Flow of Sagittarius A* Fueled By Stellar Winds
S. M. Ressler , , E. Quataert , J. M. Stone Departments of Astronomy & Physics, Theoretical Astrophysics Center, University of California, Berkeley, CA 94720 Kavli Institute for Theoretical Physics, University of California Santa Barbara, Santa Barbara, CA 93107 Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544
15 January 2020
ABSTRACT
We study the flow structure in 3D magnetohydrodynamic (MHD) simulations of accretiononto Sagittarius A* via the magnetized winds of the orbiting Wolf-Rayet stars. These simula-tions cover over 3 orders of magnitude in radius to reach ≈
300 gravitational radii, with onlyone poorly constrained parameter (the magnetic field in the stellar winds). Even for windswith relatively weak magnetic fields (e.g., plasma β ∼ ), flux freezing / compression in theinflowing gas amplifies the field to β ∼ few well before it reaches the event horizon. Overall,the dynamics, accretion rate, and spherically averaged flow profiles (e.g., density, velocity)in our MHD simulations are remarkably similar to analogous hydrodynamic simulations. Weattribute this to the broad distribution of angular momentum provided by the stellar winds,which sources accretion even absent much angular momentum transport. We find that themagneto-rotational instability is not important because of i) strong magnetic fields that areamplified by flux freezing / compression, and ii) the rapid inflow / outflow times of the gas andine ffi cient radiative cooling preclude circularization. The primary e ff ect of magnetic fields isthat they drive a polar outflow that is absent in hydrodynamics. The dynamical state of theaccretion flow found in our simulations is unlike the rotationally supported tori used as initialconditions in horizon scale simulations, which could have implications for models being usedto interpret Event Horizon Telescope and GRAVITY observations of Sgr A*. Key words:
Galaxy: centre – accretion, accretion discs –hydrodynamics – stars: Wolf-Rayet– X-rays: ISM – black hole physics
The accretion system immediately surrounding Sagittarius A* (SgrA*), the supermassive black hole in the centre of The Milky Way,o ff ers an unparalleled view of the diverse physical processes at playin galactic nuclei. Compared to other active galactic nuclei, the lu-minosity of the black hole is strikingly small, only ∼ − times theEddington limit, this places it firmly into the regime of the well-studied Radiatively Ine ffi cient Accretion Flow (RIAF) models (seeYuan & Narayan 2014 for an extensive review). The proximity ofthe Galactic Centre allows for the environment immediately sur-rounding the black hole to be spatially resolved, including (cid:38) ff et al. 2003), and the ionized mini-spirals streaming inwards sur-rounded by the cold, molecular circumnuclear disc. Direct con-straints on the near horizon environment are now possible with thedetection of several localized infrared flares orbiting the black holewithin ∼
10 gravitational radii ( r g ≡ GM / c , where M is the mass ofthe black hole, G is the gravitational constant, and c is the speed oflight) by GRAVITY (Gravity Collaboration et al. 2018) and the first resolved mm images by the Event Horizon Telescope (Doelemanet al. 2009; Event Horizon Telescope Collaboration et al. 2019a,b)soon to come. With such a wealth of observational data, Sgr A*can be used as a test-bed of accretion models in a way that no othersystem can.It is generally believed that the black hole’s gas supply isprimarily set by the stellar winds of the ∼
30 Wolf-Rayet (WR)stars orbiting at distances of ∼ ∼ keV temperatures, producingX-rays around the Bondi radius that are well resolved by Chan-dra (Bagano ff et al. 2003). However, a spherical Bondi estimatevastly over-predicts the observed Faraday rotation of the linearlypolarized radio emission (Agol 2000; Quataert & Gruzinov 2000a;Bower et al. 2003; Marrone et al. 2007). Instead, only a small frac-tion (cid:46) − of this gas reaches the horizon. It is this material thatproduces the X-ray and infrared flares as well as the 230 GHz emis-sion targeted by EHT.What exactly prevents most of the material at the Bondi ra-dius from accreting is still an open debate. Several viable modelshave been proposed, including those that appeal to strong outflows c (cid:13) a r X i v : . [ a s t r o - ph . H E ] J a n S. M. Ressler, E. Quataert, J. M. Stone (Blandford & Begelman 1999) and those that appeal to convec-tive instabilities that trap gas in circulating eddies (Stone, Pringle& Begelman 1999; Narayan, Igumenshchev & Abramowicz 2000;Quataert & Gruzinov 2000b; Igumenshchev & Narayan 2002; Pen,Matzner & Wong 2003). The range of models corresponds to a de-pendence of density on radius between the two extremes of r − / and r − / , with the combination of multiple observational estimatesat ∼ ff erent radii supporting r − in the inner regions of the flow(Gillessen et al. 2019) with a potential break near the Bondi ra-dius (Wang et al. 2013). Another key consideration is the angularmomentum of the gas being fed at large radii. In the absence ofmagnetic fields or other processes, gas in axisymmetric flows canonly accrete if it has a specific angular momentum (roughly) lessthan the Keplerian value at the event horizon. On the other handthe magnetorotational instability (MRI; Balbus & Hawley 1991)can amplify an initially weak field causing gas to accrete whilealso driving strong magnetically dominated outflows in the polarregions.Most simulations of accretion onto low luminosity AGN op-erate either explicitly or implicitly on the assumption that the MRIis the primary driver of accretion. For instance, General Relativis-tic Magnetohydrodynamic (GRMHD) simulations used to modelthe horizon-scale accretion flow in the Galactic Centre (e.g., DeVilliers & Hawley 2003; Gammie, McKinney & T´oth 2003; McK-inney & Gammie 2004; Mo´scibrodzka et al. 2009; Narayan et al.2012; Sa¸dowski et al. 2013; Mo´scibrodzka et al. 2014; Chan et al.2015; Event Horizon Telescope Collaboration et al. 2019c; see alsoPorth et al. 2019 for a recent GRMHD code comparison) almostuniformly start from equilibrium tori (e.g. Fishbone & Moncrief1976; Penna, Kulkarni & Narayan 2013) seeded with weak mag-netic fields that are unstable to the MRI. No low angular momentumgas is initially present. In this picture, understanding the physics ofthe MRI and how it depends on physical parameters like the netvertical flux in the disc or numerical parameters like resolution isessential for understanding accretion physics.Sgr A* is unique among AGN in that we can plausibly ex-pect to directly model the accretion of gas from large radii whereit is originally sourced by the winds of the Wolf-Rayet stars. Sincethe hydrodynamic properties of these winds (Martins et al. 2007;Yusef-Zadeh et al. 2015) as well as the orbits of the star themselves(Paumard et al. 2006; Lu et al. 2009) can be reasonably estimatedfrom observations, the freedom in our modeling is limited mainlyto the magnetic properties of the winds, which are less well known.In principle, a simulation covering a large enough dynamical rangein radius could self consistently track the gas from the stellar windsas it falls into the black hole, determining the dominant physicalprocesses responsible for accretion and directly connecting the ac-cretion rate, density profile, and outflow properties of the system tothe observations at parsec scales.With this motivation, Cuadra et al. (2005, 2006); Cuadra,Nayakshin & Martins (2008) studied wind-fed accretion in theGalactic Centre with a realistic treatment of stellar winds andCuadra, Nayakshin & Wang (2015); Russell, Wang & Cuadra(2017) added a “subgrid” model to study how feedback from theblack hole a ff ects the X-ray emission. In Ressler, Quataert & Stone(2018) (RQS18) we built on this key earlier work by treating thewinds of the WR stars as source terms of mass, momentum, andenergy in hydrodynamic simulations encompassing the radial rangespanning from ∼ ∼ × − pc ( ∼ r g ). One key resultof RQS18 was that even in hydrodynamic simulations the accretionrate onto the black hole is significant and comparable to previousobservational estimates (e.g., Marrone et al. 2007; Shcherbakov & Bagano ff ∼ ff ect magneticfields would have. Would the rotating gas be unstable to the MRI?Would the MRI growth time be short enough compared to the in-flow / outflow time in order to significantly e ff ect the flow structure?If so, how is the net accretion rate altered? How significant are largescale magnetic torques in transporting angular momentum? Theseand more are the questions we address in this work.Ressler, Quataert & Stone (2019) (RQS19) presented amethodology for modeling the accretion of magnetized stellarwinds by introducing additional source terms to account for theazimuthal field in each wind. In that work, we showed that a singlesimulation of fueling Sgr A* with magnetized winds can satisfy anumber of observational constraints, providing a convincing argu-ment that our model is a reasonable representation of the accretionflow in the Galactic Centre. First, our simulations reproduce the to-tal X-ray luminosity observed by Chandra (Bagano ff et al. 2003),meaning that we capture at least a majority of the hot, di ff use gas atlarge radii. Second, our simulations reproduce the r − density scal-ing inferred from observations that were taken over a large radialrange (Gillessen et al. 2019), implying that we are capturing a ma-jority of the gas at all radii and that our inflow / outflow rates havethe right radial dependence. Third, our simulations can reproducethe magnitude of the RM of both the magnetar (produced at r (cid:38) . r (cid:46) − pc, Mar-rone et al. 2007), demonstrating that our calculated magnetic fieldstrengths are reasonable at both small and large scales. Fourth, oursimulations can plausibly explain the time variability of the RMof Sgr A* (Bower et al. 2018), the time variability of the magne-tar’s RM, as well as the time variable part of its dispersion measure(Desvignes et al. 2018). In this work, we study the dynamics ofthis model in more detail, with the primary focus of determiningthe degree to which magnetic fields alter the flow structure seen inpurely hydrodynamic simulations (e.g., Cuadra, Nayakshin & Mar-tins 2008, RQS18).One key open question regarding the horizon scale accretionflow onto Sgr A* is whether or not it is magnetically arrested.This state can occur when coherent magnetic flux is consistentlyaccreted onto the black hole and amplified by (e.g.) flux freez-ing (Shvartsman 1971) to the point that the magnetic pressure be-comes large enough to halt the inflow of matter. This configurationis generally referred to as a “Magnetically Arrested Disc” (MAD;Narayan, Igumenshchev & Abramowicz 2003). Simulations showthat accretion in the MAD state are much more time variable thantheir Standard and Normal Evolution (SANE) counterparts, havemuch stronger jets, and the bulk of accretion occurs along thintransient streams that are able to penetrate to the horizon (Igumen-shchev, Narayan & Abramowicz 2003; Tchekhovskoy, Narayan& McKinney 2011). The periodicity of the polarization vectorsof the localized infrared flares detected by GRAVITY favors thepresence of a strong, coherent, vertical magnetic field at horizonscales (Gravity Collaboration et al. 2018), large enough to poten-tially be in a MAD state. MAD models are also favored over SANEmodels of the emission from the supermassive black hole in M87because they more naturally account for the energetics of the jet c (cid:13) , 000–000 eeding Sagittarius A* (Event Horizon Telescope Collaboration et al. 2019c). Simulationsof magnetized wind accretion in the Galactic Centre are uniquelyequipped to address the question of whether or not the winds of theWR stars can provide enough coherent magnetic flux for Sgr A* tobecome MAD.This paper is organized as follows. § § § § § § Our simulations use the conservative, grid-based code
Athena++ coupled with the model for magnetized winds outlined in RQS19.This model is an extension of the purely hydrodynamic wind modelpresented in RQS18 and treats the winds of the WR stars as sourcesof mass, momentum, energy, and magnetic field that move on fixedKeplerian orbits. The hydrodynamic properties of the winds are pa-rameterised by their mass loss rates, ˙ M w and their wind speeds, v w .The magnetic fields of the winds are purely toroidal as defined withrespect to the spin axes of the stars and have magnitudes set by theparameter β w , defined by the ratio between the ram pressure of thewind and its magnetic pressure at the equator (a ratio that is inde-pendent of radius in an ideal stellar wind).In brief, the equations solved in our simulations are ∂ρ∂ t + ∇ · ( ρ v ) = f ˙ ρ w ∂ ( ρ v ) ∂ t + ∇ · (cid:32) P tot I + ρ vv − BB π (cid:33) = − ρ GM BH r ˆ r + f ˙ ρ w (cid:104) v w , net (cid:105) ∂ ( E ) ∂ t + ∇ · [( E + P tot ) v − v · BB ] = − ρ GM BH r v · ˆ r + (cid:104) ˙ E B (cid:105) + f ˙ ρ w (cid:104)| v w , net | (cid:105) − Q − ∂ B ∂ t − ∇ × ( v × B ) = ∇ × (cid:16) ˜ E w (cid:17) , (1)where ρ is the mass density, v is the velocity vector, B is the mag-netic field vector, E = / ρ v + P / ( γ − + B / (8 π ) is the totalenergy, γ = / P tot = P + B / (8 π ) is the total pressure including both thermal andmagnetic contributions, Q − is the cooling rate per unit volume dueto radiative losses caused by optically thin bremsstrahlung and linecooling (using Z = Z (cid:12) and X = f is the fraction of the cell byvolume contained in the wind, ˙ ρ w = ˙ M w / V w , V w = π/ r w , v w , net is the wind speed in the fixed frame of the grid, (cid:104)(cid:105) denotes a vol-ume average over the cell, ˙ E B is the magnetic energy source termprovided by the winds, and ˜ E w is the average of the wind sourceelectric field, E w , over the appropriate cell edge (see Equations 22-24 of Stone et al. 2008). Each ‘wind’ has a radius of r w ≈ √ ∆ x ,where ∆ x is the edge length of the cell containing the centre of thestar. Athena++ is rewrite of the widely used
Athena code (Stone et al.2008) in the c ++ language. For the latest version of Athena++ , seehttps: // princetonuniversity.github.io / athena / . To test that our implementation of the source terms drives a magne-tized wind with the desired properties, we place a stationary windin the centre of a 3D, 1 pc grid and run for 4 wind crossing times.The mass-loss rate of the wind is ˙ M w = − M (cid:12) / yr, the wind speedis v w = / s, and radiative cooling is disabled. We choose β w =
100 to ensure that the magnetic field is non-negligible butrelatively weak.The left panel of Figure 1 shows that the angle-averaged ϕ component of the magnetic field matches the analytic expecta-tion, scaling as r − as determined by flux conservation. The othercomponents of the field are negligible. The right panel of Figure1 shows the dependence of B ϕ on polar angle θ at a distance of10 times the wind radius. Since the ϕ source term in the induc-tion equation is ∝ sin( θ ), a dynamically unimportant magnetic fieldwould also be ∝ sin( θ ). For β w = P m ∝ sin( θ ) pushes the gas away from the midplane andtowards the poles. This leads to the field being slightly lower thanprescribed in the midplane and slightly higher near the poles, by afactor of (cid:46) θ ) dependence of the sourceterm in the induction equation was chosen simply because it van-ishes at 0 and π and not based on detailed modeling of the angularstructure of MHD winds.The angular structure of the wind seen in Figure 2 becomeseven more pronounced for β w =
10, where the magnetic pressureis now 10% of the ram pressure. Here the wind becomes highlycollimated, as shown by the left panel of Figure 2, where the densityis now concentrated at the poles and the magnetic field is roughlyindependent of polar angle. At the same time, the total mass outflowrate and angle averaged wind speeds are still in good agreementwith the intended v w and ˙ M w as shown in the right panel of Figure2. As noted in RQS19, when β w is further decreased to (cid:46) β w (cid:38)
10 and focus primarily on β w (cid:62) . We show later that our results are insensitive to β w for β w ∼ − . / Initial Conditions
We use a base grid in Cartesian coordinates of 128 covering a boxsize of (2 pc) with an additional 9 levels of nested static mesh re-finement. This doubles the e ff ective resolution every factor of ∼ ≈ × − pc. No additional mesh refinement is used(e.g.) near the stellar wind source terms. The inner boundary ofour simulation is approximately spherical, with a radius equal totwice the length of an edge of the smallest cubic cell, correspond-ing to r in ≈ × − pc ≈ . × − (cid:48)(cid:48) ≈ r g . All cells with Note that in RQS18 and RQS19 it was stated that the box size of oursimulations was (1 pc) . This was an error. The box size in both works wasactually (2 pc) , as it is here.c (cid:13)000
We use a base grid in Cartesian coordinates of 128 covering a boxsize of (2 pc) with an additional 9 levels of nested static mesh re-finement. This doubles the e ff ective resolution every factor of ∼ ≈ × − pc. No additional mesh refinement is used(e.g.) near the stellar wind source terms. The inner boundary ofour simulation is approximately spherical, with a radius equal totwice the length of an edge of the smallest cubic cell, correspond-ing to r in ≈ × − pc ≈ . × − (cid:48)(cid:48) ≈ r g . All cells with Note that in RQS18 and RQS19 it was stated that the box size of oursimulations was (1 pc) . This was an error. The box size in both works wasactually (2 pc) , as it is here.c (cid:13)000 , 000–000 S. M. Ressler, E. Quataert, J. M. Stone r / r w r B / B , an B /max B B ( r = 10 r w ) B ( r = 10 r w ) sin( ) Figure 1.
Left: Angle averaged azimuthal magnetic field, (cid:104) B ϕ (cid:105) , normalized to the r (cid:38) r w analytic expectation for β w = (solid red) and to its peak value(dashed blue). The agreement with the analytic solution is excellent. Right: θ dependence of B ϕ at 10 wind radii (solid green) compared to sin( θ ) (dashedpurple). The magnetic field is slightly more spread out in θ than sin( θ ) because the imbalanced magnetic pressure tends to push the gas towards the poles. Thise ff ect is more extreme for β w =
10 (Figure 2). B ( r = 10 r w ) B ( r = 10 r w ) sin( ) log ( / hydro ) r / r w B / B , an M / M w v r / v w Figure 2.
Left: θ dependence of the azimuthal magnetic field, B ϕ (solid green), and the logarithm of the mass density divided by the purely hydrodynamicsolution, log ( ρ/ρ hydro ) (dotted orange), both evaluated at 10 wind radii for β w =
10 and compared with sin( θ ) (dashed purple). Right: angle averaged (cid:104) B ϕ (cid:105) (solid red), accretion rate normalized to the expected value, (cid:104) ˙ M (cid:105) / ˙ M w (dashed blue), and radial velocity normalized to the expected value, (cid:104) v r (cid:105) / v w (dottedblack). Compared to the β w =
100 case in Figure 1, the magnetic field is now strong enough to collimate the wind, enhancing the density by almost a factor of100 at the poles. Despite this, the net accretion rate and wind speed are still consistent with the input parameters. Nonetheless, we focus on β w (cid:62)
100 for oursimulations to avoid the collimating e ff ect of magnetic fields on the stellar winds. centre points within this radius are set to have zero velocity andfloored density / pressure, yet the magnetic field is allowed to freelyevolve. Even though we generally expect the solution just outsidethis radius in our simulations to have large inwards radial speeds,we chose to set the velocity to zero inside r in for simplicity andhave found that it does not strongly a ff ect our results. First, we havetested that simulations using this inner boundary condition are ableto successfully reproduce a spherical Bondi flow even when thesonic radius is smaller than r in (that is, when gas within the innerboundary radius is in causal contact with gas at larger radii), andsecond, the simulations presented in RQS18 using the same innerboundary showed that − v r just outside the inner boundary was stillable to reach the local sound speed (see Figure 11 in that work).The outer boundary of each of our simulations is set at the faces ofthe computational box using “outflow” conditions, where primitivevariables are simply copied from the nearest grid cell into the ghostzones. Our simulations use the Harten-Lax-van Leer + Einfeldt(HLLE; Einfeldt 1988) Riemann solver and 2nd order piece-wiselinear reconstruction on the primitive variables.For the WR stars, we use the orbits, mass loss rates, and windspeeds exactly as described in RQS18, drawing primarily fromMartins et al. (2007), Cuadra, Nayakshin & Martins (2008), Pau-mard et al. (2006), and Gillessen et al. (2017). These values di ff erslightly from those used in RQS19, where we modified the massloss rates and wind speeds of four stars (within reasonable system-atic observational uncertainties) to show that our simulations couldreproduce the observed RM of the Galactic Centre magnetar. Eachwind is given a randomly chosen direction for its spin axis that de-termines the azimuthal direction for the magnetic field; this randomselection is made only once for each star so that each simulation werun has the same set of spin axes. Note that RQS19 used a di ff er-ent set of spin axes, but we have found our results insensitive tothis choice. Here we also use a value of r in that is a factor of 2 c (cid:13) , 000–000 eeding Sagittarius A* Figure 3.
2D slice in the plane of the sky of electron number density overplotted with projected magnetic field lines for the inner ∼ β w = β w = (top right), β w = (bottom left), and hydrodynamic (bottom right) simulations. Each ‘star’ in our simulation provides a purely toroidalmagnetic field with direction determined by the random, independently chosen spin axes of the stars. No significant di ff erence is seen in the simulations at thisscale because the magnetic fields are relatively weak compared to the ram and thermal pressures of the gas. smaller than RQS19. We ran a total of 5 simulations; four in MHDwith β w = , , and 10 , and one in hydrodynamics (i.e., β w → ∞ ).In § β w =
10 in our modelbecome highly collimated. We have found that this collimationhas nontrivial e ff ects on the resulting dynamics of the inner accre-tion flow (in particular, altering the angular momentum directionat small radii) in a way that makes separating the e ff ects of largemagnetic fields from this extra hydrodynamic consideration di ffi -cult. Furthermore, the precise nature of this collimation dependson our choice of angular dependence of ˜ E w in Equation (1), whichwas arbitrary. Thus we do not find it instructive to include β w = β w = − simulations are consistent with those de-rived from the β w =
10 simulations that we have run.We initialize each simulation with floored density and pres-sure, zero velocity, and no magnetic field, starting at 1.1 kyr in thepast. Here, for consistency with RQS18 January 1, 2017 is definedas the present day, i.e., t =
0. The simulations are run for 1.3 kyr toa final time of t f =
200 yr. In Appendix B we argue that our resultsare independent of the arbitrary choice of starting our simulations1.1 kyr in the past by comparing with simulations that start 9 kyr inthe past. Finally, we use floors on the density and pressure (seeRQS18), and a ceiling on the Alfven speed (which is e ff ectivelyan additional floor on density; see RQS19).To ensure that our results are well converged, we ran an ad-ditional simulation that used a factor of 4 finer resolution within ∼ Figure 3 shows a 1 pc
2D slice in the plane of the sky (centredon the black hole) of the electron number density over-plotted withmagnetic field lines for β w = , , and 10 compared to thehydrodynamic case. Magnetic fields do not significantly alter thedynamics at this scale because even for β w =
100 the magneticpressure in the winds is insignificant compared to their ram pres-sure. Thus, all panels are nearly identical. Slices of the temperaturesshow similarly small di ff erences from the right panel of Figure 7 inRQS18 and are thus not included here. Figure 3 shows that, as de-sired, the magnetic fields lines wrap around the “stars,” which showup as dense circles typically surrounded by bow shocks. Again, c (cid:13)000
100 the magneticpressure in the winds is insignificant compared to their ram pres-sure. Thus, all panels are nearly identical. Slices of the temperaturesshow similarly small di ff erences from the right panel of Figure 7 inRQS18 and are thus not included here. Figure 3 shows that, as de-sired, the magnetic fields lines wrap around the “stars,” which showup as dense circles typically surrounded by bow shocks. Again, c (cid:13)000 , 000–000 S. M. Ressler, E. Quataert, J. M. Stone
Figure 4.
2D slice in the plane of the sky of electron number density overplotted with projected magnetic field lines (top four panels) and temperature (bottomfour panels) for the inner ∼ β w = and β w = simulationsare more dynamically important and thus clear di ff erences are seen in the density distribution compared to the hydrodynamic case. In addition, the larger thefield strength in the winds, the larger the spatial scale over which the field lines are coherent. c (cid:13) , 000–000 eeding Sagittarius A* since the field is not dynamically important at this scale the fieldlines for di ff erent β w ’s all have essentially the same geometry.The top four panels of Figure 4 again show 2D slices of theelectron number density overplotted with field lines but on a scaleof ∼ β w = run still looks similar to the hydrodynamic simulation,the β w = and (particularly) β w =
100 runs show significantdi ff erences. This is because, as we shall show, the field in the lattercases starts to become dynamically important at this scale. The fieldlines become increasingly ordered with decreasing β w , going frommostly tangled for β w = (where the field is easily dragged alongwith the flow) to mostly coherent for β w =
100 (where the field canresist the gas motion). This will have important implications for thefield geometry at small radii in § β w = simulation alone, a prominent large scale,hot, collimated outflow can be seen at particular times (at t = ∼ T (cid:62) × K and spanning t = −
540 yr to t = −
480 yr in 20 yrincrements. In the initial frame, no clear outflow structure is seen,only strong shocks between winds. As time progresses, however, athin (cid:62) K outflow appears coming out from the right side of theblack hole, ∼ parallel to the angular momentum direction at thistime. This outflow is magnetically driven and originates at smallradii, as we shall show in the next section.Figure 6 shows the angle-averaged root-mean-squared (rms)magnetic field strength and plasma β ≡ P / P m , where P m is themagnetic pressure, as a function of radius for di ff erent values of β w . As expected, at large radii ( (cid:38) . β scalesimply as (cid:112) /β w and β w , respectively. At small radii ( (cid:46) − pc),however, there is a much weaker dependence of the rms field and β on β w . In fact, both β w = and β w = reach β ∼ ∼
1G field strengths by 10 − pc. Even in the β w = case, the fieldstrength ( β ) at small radii is only a factor of ∼ β w =
10 simulation. This is why RQS19 found thatthe rotation measure of Sgr A* and the net vertical flux threadingthe inner boundary of the simulation were roughly independent of β w , since both quantities are set by the field at the innermost radii.Though this result might seem like a clear signature of a field regu-lated by the magnetorotational instability (MRI), we argue in § ff erences at the ∼ / temperature (Figure 4) and the fact that the flowcan reach β of ∼ a few over orders of magnitude in radius (Figure6), the radially averaged gas properties in the MHD simulations re-main strikingly similar to the hydrodynamic results at all radii evenfor β w = β w = net accretion rate including both inflow andoutflow components, i.e., (cid:104) ˙ M (cid:105) = (cid:104) πρ v r r (cid:105) . Though there can be aslarge as a factor of three di ff erence in accretion rate (correspondingto a di ff erence in density at small radii) at specific times, on av-erage, the accretion rate through the inner boundary is unchangedby the presence of the magnetic fields, falling between ∼ × − M (cid:12) / yr . The di ff erences in the average sound speed andradial velocity are negligible. We have tested that this result alsoholds for di ff erent values of the inner boundary radius.The net accretion rates shown in Figure 7 (and in all our sim-ulations) are negative and roughly constant in radius from the innerboundary out to r ∼ − pc. Then it rises in magnitude between r ∼ − pc and r ∼ − pc with a sign that fluctuates with time.Finally for r (cid:38) − pc, it is positive and increasing with radius. Anet accretion rate that is constant in radius is expected for a flow insteady-state in the absence of source terms. Our simulations, how-ever, have a time-variable source of mass describing the contribu-tions of stellar winds, depending on the stellar wind properties andstellar locations, the latter of which change as the stars proceedalong their orbits. In the limit of a large number of stars, this timedependence can be small if at each radial distance from the blackhole there are always a similar number of stellar winds (or at least asimilar total mass-loss rate). This is roughly the case for the stellarwinds in our simulations for 10 − pc (cid:46) r (cid:46) . − pc (cid:46) r (cid:46) . t = × − pc (cid:46) r (cid:46) − pcso the source term in mass is time variable in this region. Becauseof this, the regions between 10 − pc (cid:46) r (cid:46) − pc never reach asteady state but instead depend on the time-dependent location ofthe stars and the properties of their winds, both of which are ob-servationally constrained. For reference, the mass-weighted inflowtime, r / (cid:104) v r ( v r < (cid:105) ρ , is shorter than the simulation run time for allradii (cid:46) . (cid:46) − pc (cid:46) r (cid:46) − pc would have reached inflow equi-librium. The time-dependence of the location of the stellar winds,however, results in the magnitude of the accretion rate in this re-gion increasing with radius though temporally fluctuating in sign.For smaller radii, however, with r (cid:46) − pc, there are no signifi-cant source terms and the dynamical time is short compared to thetime-scale for the temporal variability of the stellar winds sourcingthe flow so that the flow reaches a negative accretion rate that isconstant with radius. Finally, the angle-averaged flow properties at r (cid:38) . ff erences in the angle-averaged radial profilesof fluid quantities (Figure 13), Figure 8 shows the various time andangle-averaged components of the outwards radial force balancinggravity for our β w = simulation. This includes the thermal pres-sure force, the Lorentz force, the centrifugal force, and the radialram pressure force. The magnetic field accounts for only (cid:46) ∼
40% each and the radial ram pressure force ac-counting for ∼ β , which takes into account onlythermal pressure, is ∼ ff ect of the magnetic field is reduced because of the largecentrifugal and ram pressure contributions. Due to the chaotic nature of our simulations, the instantaneous value ofthe accretion rate at t = (cid:13)000
40% each and the radial ram pressure force ac-counting for ∼ β , which takes into account onlythermal pressure, is ∼ ff ect of the magnetic field is reduced because of the largecentrifugal and ram pressure contributions. Due to the chaotic nature of our simulations, the instantaneous value ofthe accretion rate at t = (cid:13)000 , 000–000 S. M. Ressler, E. Quataert, J. M. Stone
Figure 5.
Time series of jet formation over the course of 150 years in the β w = simulation. Plotted are 2D gas temperature slices in the plane of the sky forthe inner ∼ δ t is the time elapsed since the first snapshot.Time proceeds clockwise starting from the upper left panel. The color scale di ff ers from that used in Figure 4 and was chosen to particularly highlight thehighest temperatures. As time progresses, a collimated, high temperature outflow emerges asymmetrically from small radii until it reaches r ∼ β w = simulation and not at all in the higher β w simulations. To facilitate analysis of the accretion flow at small radii it is usefulto define time intervals over which the angular momentum vectorof the gas is relatively constant in time. Due to the stochastic natureof the simulations, this occurs at di ff erent times for each run, oftennot centred at t =
0. The purpose of this analysis, however, is to un-derstand the general properties of the accretion disc, outflow, andmagnetic field structure, not to make overly specific predictions forthe present day. We expect that the intervals we choose are repre-sentative of the general accretion flow dynamics and structure.Figure 9 shows the three components of the angle and radius-averaged (over the interval r = × − pc and r = × − pc) angular momentum direction vector as a function of time forour four simulations. We use this information to choose our par-ticular choice of time intervals for averaging the flow structure:[100 , , , − ,
0] yrs for β w = , , and the hydrodynamic simulation, respectively.All of these intervals have angular momentum directions that areapproximately constant in time and nearly aligned with the stellardisc containing about half of the WR stars. The angular momen-tum of the accretion flow is aligned with that of the stellar disc most of the time, though for the β w = simulation it has morefrequent and larger deviations from the stellar disc than in the hy-drodynamic simulation. The most significant of these is seen near t =
0, where the angular momentum of the gas in the β w = sim-ulation is nearly anti-aligned with the stellar disk for a brief ∼ magnitude of the angle and time-averagedangular momentum is similar for all simulations, being ∼ l kep for r (cid:46) . l in MHD di ff ers at most by 20% from that inhydrodynamics as shown in Figure 7).Defining a new ‘ z (cid:48) direction as the direction of the time av-eraged angular momentum vector, Figure 10 and Figure 11 showtime series of the midplane ( θ = π/
2) mass density on ∼ β w = simulations, respectively. A time series for the β w = ( β w = ) simulation is not shown but it looks qualitatively verysimilar to the β w = (hydrodynamic) case. Figure 12 shows theBernoulli parameter, a measure of how bound the gas is to the blackhole, in the same frame for all four simulations at a representativetime. These figures show that the majority of the unbound, highangular momentum gas in the midplane at large radii is providedby the closest one or two stellar winds (namely, those of E23 / IRS c (cid:13) , 000–000 eeding Sagittarius A* r (pc)10 B ( G) w = 10 w = 10 w = 10 Figure 6.
Comparison between the root-mean-squared magnetic fieldstrength (red), (cid:112) (cid:104) B (cid:105) , and plasma β (blue), (cid:104) P (cid:105) / (cid:104) P m (cid:105) , for di ff erent valuesof β w , which quantifies the magnetization of the WR stellar winds. Eventhough the field strength varies by 2 orders of magnitude at large radii (cor-responding to a 4 orders of magnitude di ff erence in β ), the field strengths atsmall radii are all within a factor of (cid:46) β ’s within a factor of (cid:46) β reaches ∼ a few. / IRS 16C) as they orbit the black hole in all simula-tions. As this material streams inwards, however a clear di ff erenceis seen in the behavior at smaller radii in the di ff erent runs. In thehydrodynamic case, each fluid element largely conserves its an-gular momentum and energy, thus remaining unbound. Gravity isonly strong enough to bend the inflowing streams of gas around theblack hole until they emerge on the other side as a spray of outflowthat sends the gas out to larger radii without much accretion. Thehigh angular momentum gas does not spend enough time at smallradii to circularize or form a disc; instead the supply of matter atsmall radii is continually being lost and replenished. The same istrue of weakly magnetized simulations (i.e., β w = and higher).In MHD with strong magnetic fields (i.e., β w = and β w = ), however, this picture is di ff erent. Now the strong fields( β ∼ a few at small radii, Figure 6) are able to e ffi ciently removesome angular momentum and energy from the gas via large-scaletorques. This results in the originally unbound material becomingbound as it falls in so that its trajectory alters to form an inwardspiral that ultimately accretes instead of spraying out the other sideto large radii. The main di ff erence, however, as we shall argue in § ϕ and time-averaged accretion rates for our four simula-tions while Figure 14 shows the same for the Bernoulli parameter.The hydrodynamic midplane structure described above results in anet outflow of high angular momentum, modestly unbound mate-rial in the midplane, while low angular momentum material freelyfalls in along the poles. The polar inflow also contains some higher r (pc)10 MHDHydro c s | v r | M (>0) M (<0) l t (yr)0.000.250.500.751.001.251.501.752.00 | M | ( M / y r ) Figure 7.
Comparison between the β w = MHD simulation (solid) and apurely hydrodynamic simulation (dashed). Top: present day angle-averagedmass density, ρ (M (cid:12) / pc ), sound speed, c s (pc / kyr), radial velocity, | v r | (pc / kyr), specific angular momentum, l (pc / kyr), and net mass accretionrate, ˙ M (M (cid:12) / kyr), as a function of distance from the black hole. Positive netaccretion rates are green, while negative net accretion rates are orange. Bot-tom: Mass accretion rate as a function of time measured at ≈ ≈ r g . Despite the relatively large magnetization of the stellar winds, the mag-netic field has an almost negligible e ff ect on the radial profiles. The smalldi ff erence in density (and hence, accretion rate) is caused by the slightlydi ff erent time dependence of the accretion rate leading to a di ff erent real-ization of the flow at t = β w . angular momentum, unbound (Be / | Φ | (cid:46) − but (cid:62)
0) materialthat eventually hits a centrifugal barrier and turns aside and addsto the midplane outflow. For β w = , where the field is relativelyweak, the same structure is seen. However, for β w = and 10 ,the hydrodynamic accretion structure is completely reversed. Forthese more magnetized flows, not only is there net inflow of bound,Be / | Φ | <
0, material in the midplane, but the energy released fromthe gas as it loses angular momentum due to magnetic torques is de-posited into an unbound, Be / | Φ | ∼
10 polar outflow. As evidencedby the fact that the net accretion rates are comparable in both cases,this outflow is similar to the one present in the hydrodynamic sim-ulation but redirected from the midplane to the poles.An additional consequence of the di ff erent poloidal dynamicsis that the β w = and 10 simulations display a stronger densitycontrast between the midplane and polar regions compared to the β w = and hydrodynamic simulations. This is seen in Figure 15,which plots the time and ϕ -averaged density folded over the mid- c (cid:13)000
10 polar outflow. As evidencedby the fact that the net accretion rates are comparable in both cases,this outflow is similar to the one present in the hydrodynamic sim-ulation but redirected from the midplane to the poles.An additional consequence of the di ff erent poloidal dynamicsis that the β w = and 10 simulations display a stronger densitycontrast between the midplane and polar regions compared to the β w = and hydrodynamic simulations. This is seen in Figure 15,which plots the time and ϕ -averaged density folded over the mid- c (cid:13)000 , 000–000 S. M. Ressler, E. Quataert, J. M. Stone r (pc)0.00.20.40.60.81.0 thermal pressureLorentzram pressure( v r )centrifugal Figure 8.
Various components of the radial force exerted on a parcel ofgas relative to the gravitational force, ρ GM / r , as a function of radius inour β w = simulation. Plotted are the angle and time-averaged radialcomponent of the pressure gradient (solid), − ∂ P /∂ r , Lorentz force (dashed),ˆ r · [( ∇ × B ) / (4 π )] × B , the v r portion of the advection derivative (dotted), − ρ v r ∂ v r /∂ r , and centrifugal force (dot-dashed), ρ v ϕ / r . The Lorentz forceis ≈
10% of the gravitational force for most radii, comparable to the “rampressure force” of the v r component of the fluid velocity. Thermal pressureand rotation each balance about 40% of gravity, providing a majority of theradial support.
600 400 200 0 2001.00.50.00.51.0 w = 10
600 400 200 0 2001.00.50.00.51.0 w = 10 L x / L L y / L L z / L
600 400 200 0 2001.00.50.00.51.0 w = 10
600 400 200 0 2001.00.50.00.51.0 hydro t (yr) Figure 9.
Angular momentum direction as a function of time for the gas inour β w = , , and hydrodynamic simulations, averaged in radiusand angle over 5 × − pc to 3 × − pc. The blue shaded regions repre-sent the time intervals that we choose to analyze the inner accretion flow,over which the angular momentum vector is relatively stable. Dashed linesrepresent the three components of the angular momentum direction vectorof the clockwise stellar disc (Paumard et al. 2006; Lu et al. 2009). plane at r = ff erent simulations. Though the “disc”of gas is still quite thick, the equatorial to polar density contrastin the most magnetized case is now a factor of ∼ ∼ In order to quantify the relative contribution of the magnetic field toangular momentum transport, we calculate the Shakura & Sunyaev(1973) α -viscosity from our simulations. We do this in the sameframes defined by the angular momentum direction during the timeintervals shown in Figure 9 as described in the preceding subsec-tion. We follow Jiang, Stone & Davis (2017) by defining the timeand angle averaged Reynold’s stress S h ≡ (cid:104) ρ v r v ϕ sin( θ ) (cid:105) − (cid:104) ρ v r (cid:105)(cid:104) v ϕ sin( θ ) (cid:105) (2)and Maxwell stress S m ≡ (cid:104) B r B ϕ sin( θ ) (cid:105) . (3)Then the Shakura & Sunyaev (1973) α viscosities are simply α h = S h / P and α m = S m / P , where we have chosen to use the thermalpressure instead of the total (thermal plus magnetic) pressure in thedenominator for fair comparison between hydrodynamic and MHDsimulations.The resulting α ’s for each simulation are plotted in Figure16. In the hydrodynamic simulation, the total stress is by defini-tion equal to the Reynolds stress. This nonzero stress even withoutmagnetic fields or other sources of viscosity can be understood byconsidering the inflow / outflow structure seen in Figure 13. Accre-tion occurs via low angular momentum (i.e., low v ϕ ) material inthe polar regions where the ϕ -averaged v r is large (i.e., close tofree-fall) and negative while the midplane consists of high angularmomentum (i.e., large v ϕ ) material with smaller in magnitude andpositive ϕ -averaged v r . Thus, (cid:104) ρ v r v ϕ sin( θ ) (cid:105) is significantly di ff er-ent than (cid:104) ρ v r (cid:105)(cid:104) v ϕ sin( θ ) (cid:105) ∝ ˙ Ml , leading to a large α h . Thus, α h is notpredominantly a turbulent viscosity but a measure of the fact thatthere is a superposition of two types of flows: low angular momen-tum accretion and high angular momentum outflow. For β w = ,the Maxwell stress provided by the magnetic field is comparableto Reynolds stress and both work together to transport angular mo-mentum inwards. This picture is altered for β w = and β w = ,where the total stress is a competition between a large Maxwellstress and a non-negligible, negative Reynolds stress (where a neg-ative stress implies transport of angular momentum inwards). Forthese more magnetized flows, the magnetic field is strong enoughto resist being wound up in the ϕ direction, providing significanttorque to rotating gas as it falls in.In all cases, the total stress, α tot = α h + α m , is similar for r (cid:46) − pc, varying between ∼ . − .
2. This simply reflectsthe overall dynamical similarity of the flows independent of β w . Ina steady state accretion flow, the total stress can be written as (fromequations 2 and 3) α tot = F J − (cid:104) ˙ M (cid:105)(cid:104) l (cid:105) π r (cid:104) P (cid:105) , (4)where F J = (cid:104) ρ v r v ϕ sin( θ ) (cid:105) − (cid:104) π r B r B ϕ sin( θ ) (cid:105) = is the constantflux of angular momentum and l is the specific angular momentum.The constant F J is set by the accretion rate and angular momentumat the inner boundary and is generally small. Thus, since ˙ M and l are relatively unchanged in an angle averaged sense going fromhydrodynamics to MHD, the total stress is unchanged. We have shown (Figure 6) that the magnitude of the magnetic fieldat small radii is only weakly dependent on β w , the parameter gov-erning the strength of the magnetic field in the stellar winds. A c (cid:13) , 000–000 eeding Sagittarius A* Figure 10.
Sequence of midplane slices of the mass density weighted by radius separated by 25 year intervals for the hydrodynamic simulation. Here δ t is thetime elapsed since the first snapshot and time proceeds clockwise starting from the upper left. Material provided by the nearest two stellar winds (E20 / IRS16C in the upper left quadrant and E23 / IRS 16SW in the upper right quadrant) streams inward but mostly has too much angular momentum to accrete withoutany redistribution of angular momentum. Instead, the streams of material ultimately hit a centrifugal boundary and then “spray” outwards on the opposite sideof the black hole from which they approached. The bulk of the gas does not circularise nor form a steady disc. natural mechanism to explain this is the magnetorotational instabil-ity, which can amplify an arbitrarily small magnetic field to reach β (cid:46)
10 in di ff erentially rotating flows, such as we have here. How-ever, we have also shown that the gas in our simulations never cir-cularizes and therefore does not spend many orbits at small radii.Consequently, there is not su ffi cient time in a Lagrangian sense forthe MRI to grow.We can further evaluate the role of the MRI by using an esti-mate of the fastest growing wavelength for perturbations given by λ MRI ,θ ≈ π | B θ | (cid:112) πρ Ω , (5)where Ω ≡ v ϕ / ( r sin( θ )) is the rotational frequency. At least twocriteria need to be met in order for the MRI to operate in numericalsimulations: 1) λ MRI ,θ needs to be resolved, that is, the cell length ∆ x needs to be (cid:28) λ MRI ,θ , and 2) λ MRI ,θ needs to be smaller than thescale height of the disc, otherwise the perturbations have no roomto grow.Figure 17 shows λ MRI ,θ in the midplane of the disc comparedto the scale height of the disc, defined as H ≡ r (cid:104) ρ | θ − π/ |(cid:105) / (cid:104) ρ (cid:105) ,and the resolution of our grid. We find that λ MRI ,θ is su ffi cientlyresolved at all radii but that it is larger than the scale height for all of our MHD simulations. This implies that even if the gas inour simulations were to circularize, which we reiterate does not infact occur, the MRI would have no room to operate. Therefore, weconclude that the MRI is not an important source of magnetic fieldamplification or angular momentum transport in our simulations.This finding, along with the fact that β ∼ a few at small radii inour lower β w runs, is also observed in MAD simulations (Igumen-shchev, Narayan & Abramowicz 2003; White, Stone & Quataert2019). Instead of the MRI, we explain the saturation of the mag-netic field at small radii displayed in Figure 6 with simple com-pression / flux freezing. An initially weak field at large radii will becompressed as it is pulled inwards by the bulk motion of the gas.It will continue to do so until β ∼ a few, when the field becomesdynamically important and starts to resist the fluid motion. At thispoint, the field maintains β ∼ a few as it continues to accrete. Forsmall β w , this happens at large radii, while as β w increases the fieldreaches β ∼ a few at progressively smaller radii. If we were able toreach even smaller radii with our simulations, we predict that eventhe β w = run will ultimately reach β of ∼ a few and the fieldwould become dynamically important (see also § c (cid:13)000
10 in di ff erentially rotating flows, such as we have here. How-ever, we have also shown that the gas in our simulations never cir-cularizes and therefore does not spend many orbits at small radii.Consequently, there is not su ffi cient time in a Lagrangian sense forthe MRI to grow.We can further evaluate the role of the MRI by using an esti-mate of the fastest growing wavelength for perturbations given by λ MRI ,θ ≈ π | B θ | (cid:112) πρ Ω , (5)where Ω ≡ v ϕ / ( r sin( θ )) is the rotational frequency. At least twocriteria need to be met in order for the MRI to operate in numericalsimulations: 1) λ MRI ,θ needs to be resolved, that is, the cell length ∆ x needs to be (cid:28) λ MRI ,θ , and 2) λ MRI ,θ needs to be smaller than thescale height of the disc, otherwise the perturbations have no roomto grow.Figure 17 shows λ MRI ,θ in the midplane of the disc comparedto the scale height of the disc, defined as H ≡ r (cid:104) ρ | θ − π/ |(cid:105) / (cid:104) ρ (cid:105) ,and the resolution of our grid. We find that λ MRI ,θ is su ffi cientlyresolved at all radii but that it is larger than the scale height for all of our MHD simulations. This implies that even if the gas inour simulations were to circularize, which we reiterate does not infact occur, the MRI would have no room to operate. Therefore, weconclude that the MRI is not an important source of magnetic fieldamplification or angular momentum transport in our simulations.This finding, along with the fact that β ∼ a few at small radii inour lower β w runs, is also observed in MAD simulations (Igumen-shchev, Narayan & Abramowicz 2003; White, Stone & Quataert2019). Instead of the MRI, we explain the saturation of the mag-netic field at small radii displayed in Figure 6 with simple com-pression / flux freezing. An initially weak field at large radii will becompressed as it is pulled inwards by the bulk motion of the gas.It will continue to do so until β ∼ a few, when the field becomesdynamically important and starts to resist the fluid motion. At thispoint, the field maintains β ∼ a few as it continues to accrete. Forsmall β w , this happens at large radii, while as β w increases the fieldreaches β ∼ a few at progressively smaller radii. If we were able toreach even smaller radii with our simulations, we predict that eventhe β w = run will ultimately reach β of ∼ a few and the fieldwould become dynamically important (see also § c (cid:13)000 , 000–000 S. M. Ressler, E. Quataert, J. M. Stone
Figure 11.
Sequence of midplane slices of the mass density weighted by radius separated by 25 year intervals for the β w = simulation. Time proceedsclockwise starting from the upper left panel. Note that δ t = ff erent absolute time relative to Figure 10. Instead of simply streaming in and “spraying”outwards as seen in the hydrodynamic case (Figure 10), strong magnetic fields are able to redirect the outflowing, high angular momentum gas towards thepolar regions so that the midplane slice pictured here is mostly comprised of spiraling inflow. As in the hydrodynamic simulation, gas does not truly circulariseinto a disc but either accretes or outflows after only (cid:46) a few orbits around the black hole. While magnetic fields do provide a non-negligible torque that canremove angular momentum from the gas, this torque has limited time to operate and does not significantly modify the accretion rate, which is very similarin the hydrodynamic and MHD simulations (see Figure 7). The two stellar winds providing most of the material in this plot are E20 / IRS 16C in bottom leftquadrant and E23 / IRS 16SW in the upper left quadrant.
We now turn our attention to the structure of the magnetic fieldat small radii. RQS19 predicted that the amount of magnetic fluxultimately threading the inner radii of the domain, φ in , is roughlyinsensitive to β w and roughly constant in time, falling between ≈ − ≈
50 (Tchekhovskoy, Narayan & McKinney 2011).To demonstrate this result more clearly, Figure 18 plots φ in , definedas φ in ≡ / (cid:82) | B r | r d Ω r (cid:113) | ˙ M | v kep (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) r = r in , (6)for our three MHD simulations. Across four orders of magnitude in β w , φ in varies by only a few, and for each run it is roughly constantin time. Averaged over the interval (-100 yr, 100 yr), the values are4.4, 3.5, and 1.1 for β w of 10 , 10 , and 10 , respectively. These dif-ferences in φ in are even smaller when extrapolated to smaller radii, which we do in § β w , φ in at the horizon will be around the β w = value shown in Fig-ure 18, independent of β w . The result that φ in becomes quasi-steadydespite the fact that matter is continually being accreted (bottompanel of Figure 7) is noteworthy. We hypothesize that this is a con-sequence of the magnetic field being accreted changing directionwith time so that the incoming field reconnects with the field in theboundary in a way that regulates the value of φ in . If instead the in-coming field had the same orientation at all times, φ in would show acontinual rise until the field threading the boundary became strongenough to arrest accretion. Alternative possibilities include that theoutflow preferentially removes magnetic fields, or that a balance ofadvection and di ff usion regulates the value of φ in (as seen in simu-lations of magnetically “elevated” discs, Zhu & Stone 2018; Mishraet al. 2019).It is important to note, however, that the amount of net mag-netic flux threading the event horizon required for a simulation toreach the MAD state in GR ( φ in ≈
50) is not necessarily the same forthe Newtonian simulations we have here. GR e ff ects cause gravity c (cid:13) , 000–000 eeding Sagittarius A* Figure 12.
Midplane contours of the Bernoulli parameter, Be ≡ | v | / + γ/ ( γ − P /ρ − GM / r , divided by the gravitational potential, | Φ | = GM / r , in the fourdi ff erent simulations at representative times. Orange denotes bound material while purple denotes unbound material. Absent magnetic fields, the relativelyhigh angular momentum gas provided by the nearby stellar winds is mostly unbound with too much angular momentum to accrete (see also Figure 10). Strongmagnetic fields (as present in the β w = and 10 simulations), however, can torque the gas enough that it loses some angular momentum and becomesmoderately bound to the black hole. The energy released by this process drives polar outflow. near the event horizon to be e ff ectively stronger, requiring moremagnetic flux (i.e. more magnetic pressure) to arrest the accre-tion flow. In Newtonian simulations the threshold value for φ in islikely lower. To e ff ectively arrest the flow, the magnetic field mustbe strong enough to exert an outward radial force that is at leastas large the radial ram pressure, ρ v r , and perhaps as large as thegravitational force. Conservatively, then, if we assume that a MADstate is reached when the magnetic pressure at the inner boundaryroughly balances gravity, then by equating the gradient of the mag-netic pressure with ρ GM / r one finds a rough threshold value of φ in ∼ π (cid:112) v kep / v r (cid:12)(cid:12)(cid:12) r = r in , on the order of 6 −
10 if v r is a little less thanfree fall at r = r in as it is here. Note that in deriving this thresholdvalue on φ in we have assumed that B scales roughly as r − (Fig-ure 6) and neglected the contribution of magnetic tension. Relaxingthese assumptions, the dashed line in Figure 8 plots the outwardLorentz force (including both magnetic pressure and magnetic ten-sion forces, ∝ [ {∇× B }× B ] · ˆ r ) relative to the gravitational force. Wefind that the Lorentz force provided by the magnetic field is a factorof 10 smaller than the gravitational force at the inner boundary eventhough φ in is near the previous simple estimate of the MAD thresh-old. On the other hand, Figure 8 also shows that the Lorentz forceis as large as or larger than the v r ram pressure. So while our sim- ulations do not appear to be fully magnetically arrested based onthe fact that the accretion rate is comparable in both hydrodynamicand MHD simulations (Figure 7), the amount of magnetic flux atthe inner boundary must be only modestly less than the MAD limit.One possible concern is that just outside of the inner boundaryour simulations are unresolved, where r / ∆ x (cid:38) ∆ x is the sizeof an edge of a cubic cell. This could potentially lead to numeri-cal di ff usion of magnetic fields and prevent a larger amount of fluxfrom accumulating. Two things suggest that this locally limited res-olution is not e ff ecting our calculation of φ in . First, φ ( r ), i.e., Equa-tion (6) evaluated at a radius r instead of r in , is roughly independentof r in the β w = simulation, equal to φ in . This includes largerradii that are much better resolved where the typical grid spacing is r / ∆ x ∼ β w = simulation with an inner boundary radius that was a factorof 8 times the size of the smallest cell edge instead of our fiducialfactor of 2, meaning that the region just outside the boundary was4 times as well resolved. This additional simulation displayed ap-proximately the same value of φ in as a simulation with the same r in but only 2 cells per r in . Ultimately, much better resolved simula-tions will be required to definitively assess the impact of numericaldi ff usion on the values of φ in determined here. c (cid:13)000
10 if v r is a little less thanfree fall at r = r in as it is here. Note that in deriving this thresholdvalue on φ in we have assumed that B scales roughly as r − (Fig-ure 6) and neglected the contribution of magnetic tension. Relaxingthese assumptions, the dashed line in Figure 8 plots the outwardLorentz force (including both magnetic pressure and magnetic ten-sion forces, ∝ [ {∇× B }× B ] · ˆ r ) relative to the gravitational force. Wefind that the Lorentz force provided by the magnetic field is a factorof 10 smaller than the gravitational force at the inner boundary eventhough φ in is near the previous simple estimate of the MAD thresh-old. On the other hand, Figure 8 also shows that the Lorentz forceis as large as or larger than the v r ram pressure. So while our sim- ulations do not appear to be fully magnetically arrested based onthe fact that the accretion rate is comparable in both hydrodynamicand MHD simulations (Figure 7), the amount of magnetic flux atthe inner boundary must be only modestly less than the MAD limit.One possible concern is that just outside of the inner boundaryour simulations are unresolved, where r / ∆ x (cid:38) ∆ x is the sizeof an edge of a cubic cell. This could potentially lead to numeri-cal di ff usion of magnetic fields and prevent a larger amount of fluxfrom accumulating. Two things suggest that this locally limited res-olution is not e ff ecting our calculation of φ in . First, φ ( r ), i.e., Equa-tion (6) evaluated at a radius r instead of r in , is roughly independentof r in the β w = simulation, equal to φ in . This includes largerradii that are much better resolved where the typical grid spacing is r / ∆ x ∼ β w = simulation with an inner boundary radius that was a factorof 8 times the size of the smallest cell edge instead of our fiducialfactor of 2, meaning that the region just outside the boundary was4 times as well resolved. This additional simulation displayed ap-proximately the same value of φ in as a simulation with the same r in but only 2 cells per r in . Ultimately, much better resolved simula-tions will be required to definitively assess the impact of numericaldi ff usion on the values of φ in determined here. c (cid:13)000 , 000–000 S. M. Ressler, E. Quataert, J. M. Stone
Figure 13.
Time and ϕ -averaged accretion rate on mpc scales for each simulation, where the z direction is defined as the net angular momentum of the gas(Figure 9). Red represents inflow, blue represents outflow, and the accretion rate has been folded over the midplane and normalized such that the absolutevalue at 5 mpc is unity. For su ffi ciently large magnetic fields, the inflow / outflow structure seen in the hydrodynamic case is reversed, because the field is strongenough to redirect the outflow and confine it to the polar regions. We quantify the relative strength of the vertical magneticfield by computing the ratio between the magnitude of the averagemagnetic field vector, | B | , to the root-mean-squared field strength, (cid:112) (cid:104) B (cid:105) . For a completely vertical field this quantity would be 1,while for a completely toroidal or random field it would be 0. Fig-ure 19 plots (cid:104)| B |(cid:105) / (cid:112) (cid:104) B (cid:105) averaged over angle and the inner 5 × − pc to 3 × − pc in radius for β w = , , and 10 , where we findthat the relative strength of the ordered field increases with decreas-ing β w . This same trend is seen in the poloidal field lines (Figure 20where they are plotted on top of mass density), where the directionof the field goes from mostly random at β w = , to nearly verticalat β w = . The weaker the magnetic field, the more it is able tobe twisted by the motion of the gas and lose its original structure.The quantities plotted in Figures 19 and 20 do not e ff ectivelyprobe the ϕ component of the field, which in principle could besignificant. To quantify this, we compare the relative strength ofthe mean B ϕ to the mean B r and B θ field components. We define an‘antisymmetric’ average of B r as (cid:104) ˜ B r (cid:105) = t (cid:90) t π (cid:90) π/ (cid:90) B r d θ d ϕ dt − t (cid:90) t π (cid:90) π (cid:90) π/ B r d θ d ϕ dt , (7)where t and t are the endpoints of the time interval for aver- aging. The minus sign in Equation (7) prevents the radial fieldfrom averaging to zero over all angles. For β w = , the toroidalfield dominates with (cid:104) B ϕ (cid:105) / (cid:16) (cid:104) ˜ B r (cid:105) + (cid:104) B θ (cid:105) + (cid:104) B ϕ (cid:105) (cid:17) ≈ . − r (cid:46) × − pc because the field is weak enough to be com-pletely stretched out by the orbital motion of gas. For β w = ,on the other hand, the field is able to resist the orbital motion (seenalso in the torque that it exerts; Figure 16) and retain a predomi-nantly poloidal structure, with (cid:104) B ϕ (cid:105) / (cid:16) (cid:104) ˜ B r (cid:105) + (cid:104) B θ (cid:105) + (cid:104) B ϕ (cid:105) (cid:17) (cid:46) . r (cid:46) × − pc. This is even more true for β w = , which has (cid:104) B ϕ (cid:105) / (cid:16) (cid:104) ˜ B r (cid:105) + (cid:104) B θ (cid:105) + (cid:104) B ϕ (cid:105) (cid:17) (cid:46) . r (cid:46) × − pc. Thus far we have presented seemingly paradoxical results. On onehand, for su ffi ciently magnetized stellar winds (e.g., β w = , ),the magnetic field at small radii reaches near equi-partition with theplasma, achieving β ∼ a few, reversing the polar inflow seen in hy-drodynamic simulations, and driving accretion in the midplane. Onthe other hand, the net accretion rate through the inner boundaryand the radially averaged fluid quantities are largely una ff ected bythe presence of magnetic fields. How can this be? In the conven-tional picture of MRI driven accretion from a rotationally supported c (cid:13) , 000–000 eeding Sagittarius A* Figure 14.
Time and ϕ -averaged Bernoulli parameter, Be, normalized to the gravitational potential, | Φ | = GM / r on mpc scales for each simulation, where the z direction is defined as the net angular momentum of the gas (Figure 9). Orange represents bound, purple represents unbound, and the Bernoulli parameter hasbeen folded over the equator. Without magnetic fields, the material is slightly unbound throughout the domain except for some slightly bound material nearthe polar axis. Magnetic fields provide torque, releasing energy from the high angular momentum gas in the midplane and depositing it in the polar outflow(Figure 13). ( n o r m a li z e d ) w = 10 w = 10 w = 10 hydro Figure 15.
Time and ϕ averaged mass density at r = θ = π/ torus, it would require an improbable cooincidence, where the mid-plane accretion driven by the MRI exactly equals the original hy-drodynamic polar accretion despite the fact that they are governedby di ff erent physical considerations. As we have shown, however,our simulations do not fit this conventional picture. The gas withsignificant angular momentum clearly does not circularize into aconfiguration where the velocity is primarily in the azimuthal di-rection (e.g., Figures 10 and 11), but instead retains significant ra-dial velocity of order free fall throughout the domain. Simply put,gas accreting from large radii is quick to either flow through the in-ner boundary or flow right back out. Magnetic fields are not strongenough to modify these flows by more than order unity even at β ∼
1. Moreover, even in the hydrodynamic simulation, inflow is not oc-curring only in the poles as Figure 13 would imply but at all polarangles. It is only in an azimuthally-averaged sense that v r is positiveand small in the midplane because there is also significant outflowpresent (at di ff erent ϕ ). The primary role of magnetic fields, then, isnot to drive accretion but to redirect the outflow from the midplaneto the pole. This means that 1) the same physical processes governaccretion in the hydrodynamic and MHD simulations and 2) the netaccretion rate is essentially determined by hydrodynamic consider- c (cid:13)000
1. Moreover, even in the hydrodynamic simulation, inflow is not oc-curring only in the poles as Figure 13 would imply but at all polarangles. It is only in an azimuthally-averaged sense that v r is positiveand small in the midplane because there is also significant outflowpresent (at di ff erent ϕ ). The primary role of magnetic fields, then, isnot to drive accretion but to redirect the outflow from the midplaneto the pole. This means that 1) the same physical processes governaccretion in the hydrodynamic and MHD simulations and 2) the netaccretion rate is essentially determined by hydrodynamic consider- c (cid:13)000 , 000–000 S. M. Ressler, E. Quataert, J. M. Stone w =10 mr tot w =10 mr tot w =10 mr tot hydro rtot r (pc) Figure 16.
Comparison between the time averaged Maxwell (blue solid),Reynolds (dashed yellow), and total (dotted black) α viscosities as de-fined in § § α m ≈ .
2, with the Reynolds stress becoming negativeto compensate. r (pc)10 MRI , , w = 10 MRI , , w = 10 MRI , , w = 10 xH Figure 17.
Wavelength of the most unstable mode for the MRI computed inthe midplane for our three MHD simulations as compared to the scale heightof the disc, H , and the resolution of our grid, ∆ x . λ MRI ,θ is well resolved butis larger than the scale height of the disc. The MRI is suppressed by thestrong β ∼ few magnetic field (Figure 6) produced by compression as thegas flows in. ations, namely, the distribution of angular momentum at large radii,a quantity set by the winds of the WR stars.The lack of circularization in our simulations, the crucial fac-tor in determining this accretion structure, is at least in part due toradiative cooling being ine ffi cient at removing dissipated energy inthe gas streamers seen in Figures 10 and 11. As the gas comes inalong nearly parabolic orbits it heats up and (because it can’t cool)expands outward, making it more di ffi cult for it to circularize. Thisis analogous to the di ffi culty that simulations of tidal disruptionevents have in forming a circular accretion disc (for recent discus-sion, see, e.g., Stone et al. 2019; Lu & Bonnerot 2019).
400 300 200 100 0 100 200 t (yr) i n w = 10 w = 10 w = 10 Figure 18.
Magnetic flux threading the inner boundary as a function of timefor our three MHD simulations, φ in , in units such that the MAD value is ∼
50 in GR (see Equation 6). Orders of magnitude di ff erence β w correspondsto only a factor of ∼ few di ff erence in φ in because the field strength at smallradii is only weakly sensitive to β w (Figure 6). In all cases φ in is (cid:46)
10% ofthe GR MAD limit, but the β w = simulation is near or has reached theexpected Newtonian MAD limit of (cid:46)
10 appropriate for these simulations. t (yr)0.00.10.20.30.40.5 | B |/ B w = 10 w = 10 w = 10 Figure 19.
Magnitude of the time and angle-averaged magnetic field vectordivided by the rms magnitude of the field for β w = (solid blue), β w = (dashed orange), and β w = (dotted green). This measures the degreeto which the magnetic field is ordered, and increases with decreasing β w because stronger fields are able to resist fluid motion and more e ff ectivelyretain a coherent structure. The simulations we have performed, while modeling a radial rangeof just over 3 orders of magnitude, are not able to penetrate all theway to the event horizon of Sgr A* but have inner boundary radiistill a few hundred times farther out. Thus, it is important for usto understand how the artificially large inner boundary of our sim-ulation (which acts as the black hole) a ff ects the results. By vary-ing the inner boundary, RQS18 showed that the predicted accretionrate through the inner boundary in our hydrodynamic simulationsis ˙ M ≈ . × − ( r in / r G ) / M (cid:12) / yr, where the dependence on r in isset by the distribution of accretion rate with angular momentum atlarge radii; for a smaller inner boundary radius, less material has an-gular momentum low enough to ultimately accrete. This predictedaccretion rate is consistent with both observational constraints andemission models (see RQS18).In MHD, we have shown that even for strong magnetic fieldsthe radially averaged fluid variables are mostly unchanged goingfrom hydrodynamics to MHD (Figure 7), including the accretionrate. Thus the above relation between ˙ M and r in still holds. As wehave argued in the preceding section, this counter-intuitive result is c (cid:13) , 000–000 eeding Sagittarius A* Figure 20.
Time and ϕ -averaged mass density weighted by radius in a co-ordinate system such that the z -direction is aligned with the angular mo-mentum of the gas (Figure 9), normalized, and overplotted with magneticfield lines for β w = (top), 10 (middle), and 10 (bottom). The strongerthe field, the more it is able to resist being dragged along by the randommotions of the flow and retain a coherent structure a consequence of the fact that the supply of infalling gas at smallradii is still mostly set by the distribution of accretion rate with an-gular momentum at large radii. The gas provided by nearby stellarwinds has a typical distribution of d ˙ M / dl ≈ const. which results in˙ M in ∝ √ r (see Appendix A in RQS18). Note that since the windsemit at all angles, this is the distribution of infalling gas for boththe poles and the midplane. Figure 21 confirms this expectation,showing that ˙ M ∝ √ r for the inflow in both hydrodynamics (in thepolar region) and β w = MHD (in the midplane).Figure 6 shows that β tends to decrease with decreasing radiusuntil it reaches ∼ a few, at which point it becomes independent ofradius. For β w = β is ≈ β w = it decreases from β ≈
200 at large radii to β ≈ r ≈ × − pc and remains constant for r (cid:46) × − pc,while for β w = it decreases from β ≈ × to β ≈ r (pc)10 M ( M / k y r ) w = 10 , = 90hydro, = 170 r r (pc)10 M ( M / k y r ) w = 10 , = 170hydro, = 90 r Figure 21.
Top: Comparison between the time and ϕ -averaged accretionrate in the midplane for the β w = simulation (solid blue) vs. the polefor the hydrodynamic simulation (dashed orange). Both of these regions aredominated by inflow and show ˙ M ∝ √ r as expected for the d ˙ M / dl ≈ const.distribution provided by the stellar winds. Bottom: Comparison between thetime and ϕ -averaged accretion rate in the pole for the β w = simulation(solid blue) vs. the midplane for the hydrodynamic simulation (dashed or-ange). Both of these regions have a net outflow and show ˙ M ∝ r , implyinga roughly constant velocity outflow since ρ ∝ r − (Figure 7). inner boundary. It is natural to suppose that if the inner boundaryradius of the simulation was reduced then the β w = run wouldultimately also reach β ≈
2. Regrettably this is not something wecan test with our current computational resources; however, we can increase the size of the inner boundary and infer how β dependson r in in the same way that we used to extrapolate ˙ M in RQS18.Doing so we estimate that all models with β w (cid:54) will reach β of ∼ a few by 2 r g (the event horizon radius of a non-rotating blackhole). Thus β w ∼ is a critical value that determines whetheror not the horizon scale accretion flow will more closely resemblethe hydrodynamic simulations ( β w (cid:38) ) or the more magnetizedwind simulations ( β w (cid:46) ).Similar behavior is seen with the magnetic flux threading theinner boundary. Figure 22 shows the time-averaged φ in for β w = and β w = and four values of the inner boundary radius. As wasthe case for β , φ in is independent of r in for β w = . This is againbecause the β w = simulation has already reached β ∼ few at r (cid:29) r in . Since ˙ M ∝ r / in v kep ( r in ) ∝ r − / in , Equation (6) gives φ in ≈ const. For β w = on the other hand, φ in increases with decreasing r in . Empirically, we find in Figure 22 that for β w = , φ in ∝ r . in , c (cid:13)000
2. Regrettably this is not something wecan test with our current computational resources; however, we can increase the size of the inner boundary and infer how β dependson r in in the same way that we used to extrapolate ˙ M in RQS18.Doing so we estimate that all models with β w (cid:54) will reach β of ∼ a few by 2 r g (the event horizon radius of a non-rotating blackhole). Thus β w ∼ is a critical value that determines whetheror not the horizon scale accretion flow will more closely resemblethe hydrodynamic simulations ( β w (cid:38) ) or the more magnetizedwind simulations ( β w (cid:46) ).Similar behavior is seen with the magnetic flux threading theinner boundary. Figure 22 shows the time-averaged φ in for β w = and β w = and four values of the inner boundary radius. As wasthe case for β , φ in is independent of r in for β w = . This is againbecause the β w = simulation has already reached β ∼ few at r (cid:29) r in . Since ˙ M ∝ r / in v kep ( r in ) ∝ r − / in , Equation (6) gives φ in ≈ const. For β w = on the other hand, φ in increases with decreasing r in . Empirically, we find in Figure 22 that for β w = , φ in ∝ r . in , c (cid:13)000 , 000–000 S. M. Ressler, E. Quataert, J. M. Stone r in (pc)10 i n r w = 10 w = 10 Figure 22.
Magnetic flux threading the inner boundary of our simulation, φ in (Equation 6), averaged over the time interval (-100 yr,100 yr) and plottedas a function of inner boundary radius for β w = (blue solid) and β w = (orange dashed). φ in is independent of r in for β w = where β hasreached ∼ few (Figure 6), while it increases with decreasing r in for β w = because β is decreasing with decreasing r in . For β w = we find φ in ˜ ∝ r . and thus expect it to reach ∼ β w = value) for an inner boundaryat the event horizon. predicting that it will reach ≈ r in = × − pc ≈ r g . At thatpoint, we expect φ in to stop increasing in the same way that φ in isindependent of r in for β w = . In analyzing our simulations, we have found it instructive to com-pare and contrast our results with previous simulations in the lit-erature that considered the problem of accretion onto Sgr A* andrelated systems via large scale feeding. In this section we do so fortwo key works.
Proga & Begelman (2003a,b) (hereafter PB03A and PB03B) pre-sented the results of 2D inviscid hydrodynamic (PB03A) and MHD(PB03B) simulations of accretion onto supermassive black holes asfed by gas with a θ -dependent distribution of angular momentum atlarge radii. This approach di ff ers from the standard method of ini-tializing simulations with equilibrium tori without any feeding atlarge radii and is perhaps a better approximation of the feeding ofgas via stellar winds in the Galactic Centre. In fact, in many waysthe results of PB03A and PB03B are strikingly similar to the resultsof RQS18 and those presented here. We both find that accretion inhydrodynamic simulations occurs via low angular momentum gasfalling in mostly along the polar regions while the higher angularmomentum midplane is (on average) outflowing. We both also findthat this structure is reversed in MHD for su ffi ciently large mag-netic fields, with the low angular momentum polar inflow gettingquenched by magnetically driven polar outflow while gas in themidplane accretes. PB03B, however, found that the accretion ratein the MHD case was significantly reduced compared to the hy-drodynamic case because the induced midplane accretion was notenough to compensate for the loss of polar inflow. In this work, onthe other hand, the midplane accretion in MHD seems to roughly equal the original hydrodynamic polar inflow so that the net accre-tion rate is relatively the same in MHD and hydrodynamics.The key di ff erence lies in the structure of the high angularmomentum gas in the midplane. PB03A found that this gas wasable to circularise and build up into a nearly constant angular mo-mentum torus that blocked the inward flow of gas for polar anglesclose to the equator. We find that the high angular momentum gasin our hydrodynamic simulation never circularises but mostly flowsright back out after falling in to small radii. This is more easily ac-complished in 3D where flow streams can avoid intersecting; in 2Daxisymmetry (used in PB03A and PB03B), collisions between theinfalling and outflowing high angular momentum gas are unavoid-able and can dissipate radial kinetic energy and e ffi ciently circu-larise the material. Because of this circularisation in PB03A, byadding even a weak magnetic field, the MRI is able to grow as thegas in the torus orbits and becomes the dominant driver of accre-tion. Thus, the accretion rate in PB03B is mostly set by completelydi ff erent physical considerations (the properties of the MRI) thanin PB03A (the availability of low / zero angular momentum gas). Inour simulations, however, even in MHD the dominant source of ac-cretion is still the supply of low angular momentum gas with anorder unity correction for global torques provided by strong mag-netic fields that have been compressed to β ∼ few at small radii.This means that the local supply of mass available to accrete is setmostly by hydrodynamic considerations (i.e., the distribution of an-gular momentum vs. accretion rate provided by the nearest stellarwinds).One of the main conclusions of PB03B was that the MRIdriven accretion seen in their simulations was roughly independentof the angular momentum distribution of material sourced at largeradii, ultimately resembling simulations that are initialized with anequilibrium torus seeded by a weak magnetic field. This served aspartial motivation and justification for future work to mostly ignorelarge radii and instead focus on horizon scale ( (cid:46) − r g ) sim-ulations starting from equilibrium tori. Our results, however, sug-gest that when a more complicated treatment of accretion sourcedby stellar winds in full 3D is considered, the properties of the accre-tion flow at small radii are very strongly influenced. We suspect thatthe major source of this di ff erence is the non-axisymmetric natureof how the gas is fed by the winds, which inhibits circularization(see Janiuk, Proga & Kurosawa 2008). Pang et al. (2011) (P11) presented 3D MHD simulations in whichgas is fed through the outer boundary with uniform magnetic field,spherically symmetric density and pressure, and purely rotationalvelocity such that the specific angular momentum l , varies as sin( θ ).This set up is quite similar to PB03B except for the field geome-try (uniform vs. radial in PB03B), the field strength, and the ad-dition of a third dimension. Unlike PB03B, however, P11 foundthat global magnetic torques and not the MRI were the governingphysical mechanisms driving accretion in their simulations. Thisdi ff erence relative to PB03B is probably a consequence of the ini-tial magnetic field in P11 being much stronger, with the initial β being ∼ in P11 compared to β ∼ − in PB03B. Thiscauses the magnetic field to become dynamically important beforethe gas can circularize (if it ever would have) and also suppressesthe MRI. Dissipation of the field also leads to an unstable entropyprofile, driving convection. A steady state is reached in which thegas is in near hydrostatic equilibrium, slowly falling inwards withmagnetic pressure resisting the upward buoyancy force. Several as- c (cid:13) , 000–000 eeding Sagittarius A* pects of the P11 simulations are similar to what we have found inours. Both show a lack of circularization, both have the MRI sup-pressed by strong magnetic fields, and both find a density profileof ∼ r − with a corresponding ˙ M ∝ √ r in relationship. At the sametime, the accretion flow structures are very di ff erent in the two setsof calculations. Unlike P11, the gas in our simulations is not hy-drostatic, because the ram pressure, ρ v is comparable to or largerthan the magnetic and thermal pressures throughout the domain.We also find significant and coherent outflow, something absent inP11. The root cause of these di ff erences are related to the morecomplicated, asymmetric way that the winds of the WR stars supplygas (and magnetic field) to the black hole. While both sets of simu-lations can contain relatively large and coherent magnetic fields atlarge radii, the steady state of P11 is one in which the gas is beingsourced in an approximately spherically symmetric way with rota-tion playing only a minor role. This is because after the initial tran-sient in which the sourced gas first free falls and then transitions toa PB03A-like configuration, the build up of gas at small radii pro-vides radial pressure support for the gas at large radii, significantlyincreasing the time it takes to accrete. At this point the magnetictorques have enough time to remove a large amount of the angularmomentum at large radii, ultimately resulting in a quasi-sphericalsteady-state. In contrast, because the feeding in our simulations oc-curs in more of a stream-like manner (Figures 10 and 11), we haveno build-up of gas to provide radial pressure support. Instead, ra-dial velocities remain large and thus the e ff ect of even significant( β ∼ few) magnetic torques are limited by the short inflow / outflowtimes. This means that rotation of gas is important throughout oursimulations, with the distribution of angular momentum being theprimary determinant of the accretion rate. The main properties of nearly all GRMHD simulations used tomodel the Galactic Centre are governed by the evolution of theMRI. The supply of gas is determined by an initial rotating toruswhile low angular momentum material is absent. Our results sug-gest that for the Galactic Centre it may be critical to consider a moredetailed model for how the gas is fed into the domain, particularlywith respect to the distribution of angular momentum coming infrom larger radii.One large remaining uncertainty is how strong the outflowsare from near the horizon and whether they significantly modifythe dynamics at ∼ r g found here. Figure 5 shows that, at times,we do see strong outflows that can modify the gas out to ∼ . M ∝ √ r in , the energy released in outflowsshould scale with the inner boundary as ∝ ˙ M v kep ∝ r − / in , meaningthat the strength of this outflow would be a factor of (cid:38) √ ≈ ∼ r g hasbeen interpreted as the results of an orbitting “hot spot” embeddedin a face-on rotating flow threaded by a magnetic field primarily inthe vertical direction. Qualitatively, the geometry of the magneticfield at small radii in our β w = simulation agrees with thispicture (Figures 20), with the poloidal field being larger than B ϕ .On the other hand, the angular momentum direction of the inner accretion flow in Figure 9 is rarely as face-on as that of the best-fitorbit of the three flares ( L z / L ≈ . ± . / black hole (which they assume to be aligned),and find that an inclination angle roughly aligned with the clock-wise stellar disc is preferred. This is consistent with measurementsof the position angle of the 86 GHz and X-ray emission performedby the VLBA on a scale of ∼
10s of r g (Bower et al. 2014) and by Chandra on a much larger scale of ∼
1” (Wang et al. 2013), respec-tively. Our results are in good agreement with these observations,as the angular momentum of our innermost accretion flow is typ-ically aligned with the stellar disc (Figure 9). Forthcoming highersensitivity EHT measurements will be important for resolving thediscrepancy with the leading interpretation of the GRAVITY data.
We have presented the results of 3D simulations of accretion ontothe supermassive black hole in the Galactic Centre fueled by mag-netized stellar winds. Our simulations span a large radial range,having an outer boundary of 1 pc and an inner boundary of ∼ × − pc ( ∼ r g ), with approximately logarithmic resolution inbetween. The mass loss rates, wind speeds, and orbits of the stellarwind source terms that represent the ∼
30 WR stars are largely con-strained by observations while the relative strength of the magneticfield in each wind is parameterized by a single parameter β w , de-fined as the ratio between the ram pressure and midplane magneticpressure of the wind. In previous work, we have shown that oursimulations naturally reproduce many of the observational proper-ties of Sgr A* such as an accretion rate that is much less than theBondi estimate, a density profile ˜ ∝ r − , a total X-ray luminosityconsistent with Chandra measurements, and the rotation measureof Sgr A*. In the present paper we have focussed on the dynamicsof accretion onto Sgr A* from magnetized stellar winds.Our most significant and a priori surprising result is that theaccretion rate onto the black hole, as well as the radial profiles ofmass density, temperature, and velocity are set mostly by hydrody-namic considerations (Figure 7). This is true even when plasma β is as low as ≈ l ≈
0. Thisbroad range of angular momentum is a consequence of the fact thatthe WR stellar wind speeds ( ∼ / s) are comparable to theirorbital speeds. As a result, the stellar winds provide enough low an-gular momentum material to result in an extrapolated accretion ratethat is in good agreement with previous estimates for Sgr A*. Withmagnetic fields, global torques provide only order unity correctionsto this picture, with the accretion rate still mostly being determinedby the supply of low angular momentum gas. This is a consequenceof the fact that the high angular momentum material in our simu-lations does not circularize but mostly flows in and out with largeenough radial velocity that the inflow / outflow times are short com-pared to the time scale for magnetic stresses to redistribute angularmomentum.Simulations with strong magnetic fields at small radii dohowever di ff er from hydrodynamic simulations in one important c (cid:13)000
0. Thisbroad range of angular momentum is a consequence of the fact thatthe WR stellar wind speeds ( ∼ / s) are comparable to theirorbital speeds. As a result, the stellar winds provide enough low an-gular momentum material to result in an extrapolated accretion ratethat is in good agreement with previous estimates for Sgr A*. Withmagnetic fields, global torques provide only order unity correctionsto this picture, with the accretion rate still mostly being determinedby the supply of low angular momentum gas. This is a consequenceof the fact that the high angular momentum material in our simu-lations does not circularize but mostly flows in and out with largeenough radial velocity that the inflow / outflow times are short com-pared to the time scale for magnetic stresses to redistribute angularmomentum.Simulations with strong magnetic fields at small radii dohowever di ff er from hydrodynamic simulations in one important c (cid:13)000 , 000–000 S. M. Ressler, E. Quataert, J. M. Stone way. Hydrodynamic simulations are dominated by inflow along thepoles, while the midplane is on average outflowing but composedof both inflow and outflow components at di ff erent θ and ϕ . By con-trast, MHD simulations are dominated by inflow in the midplane,while the polar regions are on average outflowing but composed ofboth inflow and outflow components at di ff erent θ and ϕ . This isa consequence of the β ∼ few magnetic fields redirecting the highangular momentum outflow away from the midplane.We find that the magnetic field increases rapidly with radiusso that β tends to eventually saturate at small radii to a value oforder unity independent of β w (Figure 6). This growth of the fieldis caused by advection / compression as gas falls inwards and not bythe MRI. There is neither su ffi cient time for the MRI to grow be-fore gas is accreted or advected to larger radii, nor is there su ffi cientspace for the instability to grow because flux freezing builds up afield for which the most unstable MRI wavelength is comparableto or larger than the disc scale height (Figure 17). Thus the con-ventional MRI-driven torus simulations that dominate the literaturedo not appear to have reasonable initial conditions for studying ac-cretion in the Galactic Centre, at least on the scales that we cansimulate here.Elaborating on the result first presented RQS19, we haveshown that our model predicts that the magnetic flux ultimatelythreading the event horizon, φ in , will be on the order of 5, indepen-dent of β w for β w (cid:46) (Figure 22). This prediction relies on ex-trapolation to smaller radii, ignores the e ff ects of GR, and assumesthat the scaling between φ in and the inner boundary radius that wefound (Figure 22) holds at smaller radii than our simulations probe.Not accounting for GR e ff ects, this amount of flux threading thehorizon is potentially large enough to induce a MAD state near thehorizon, with the outward Lorentz force reaching (cid:46)
10% of the in-ward gravitational force in our simulations (Figure 8), which is (cid:38) the ram pressure due to v r . It is worth noting that our lower β w simu-lations do in fact display many similar properties to GRMHD MADaccretion flows, for instance, the MRI is suppressed by strong ver-tical fields, the poloidal component of the field dominates over thetoroidal component, and the angular velocity of the gas is roughlyhalf the Keplerian value (Narayan et al. 2012), though we findthe latter to be true in both MHD and hydrodynamic simulations.We also find that φ in is relatively independent of time (Figure 18).There are, however, a number of ways in which the simulations inthe present work do not appear to be fully magnetically arrested.For example, the hydrodynamic and MHD simulations show sim-ilar accretion rates and radial profiles (Figures 7 and 13), whichexplicitly demonstrates that the magnetic field is not dynamicallycritical for establishing the flow properties. In addition, our simu-lations find very di ff erent angular distributions of density (Figure15) and plasma β from traditional MAD simulations (and, in fact,most GRMHD simulations). Because of the significant amount oflow angular momentum material provided by the stellar winds, thepolar regions are only a factor or 3–5 less dense and have only aslightly lower β than the midplane. In contrast, a typical GRMHDMAD simulation would show near-vacuum, highly magnetized po-lar regions with orders of magnitude less mass than in the mid-plane, where most of the mass is condensed into a relatively thindisc compressed by the magnetic field. The evacuation of the fun-nel is caused by both the strong jets present in MADs (Igumen-shchev, Narayan & Abramowicz 2003; Tchekhovskoy, Narayan &McKinney 2011) as well as the fact that the rotating hydrodynamicconfigurations used as initial conditions in GRMHD simulationsgenerally do not allow material close to the poles (Abramowicz,Jaroszynski & Sikora 1978; Kozlowski, Jaroszynski & Abramow- icz 1978). We do not at all, however, rule out that by the time thegas reaches the event horizon of the black hole that a MAD statecould be reached. This may also depend on whether or not the blackhole is rotating, since a rotating hole provides an additional sourceof energy for outflows that could strongly impact the dynamics inthe polar region.For su ffi ciently magnetized winds (i.e., β w = here), themagnetically driven, polar outflow can, at times, reach scales aslarge as ∼ >
10 times stronger in asimulation that reached the event horizon. This is even without con-sidering the rotation of the black hole itself, which can also be ane ffi cient mechanism for driving magnetized jets (Hawley & Krolik2006). Though there is no clear signature of a jet in the Galac-tic Centre, strong outflows from Sgr A* have been invoked as onepossible explanation for the recent ALMA observations that showhighly blue-shifted emission from unbound gas in a narrow cone(Royster et al. 2019).The magnetic field structure at small radii depends on the pa-rameter β w (Figure 20). For smaller β w (10 and to a lesser extent,10 ) the field is strong enough to resist being wound up in the ϕ direction and remains mostly poloidal at small radii. For larger β w ( (cid:38) ), the field is easily dragged along with the motion of thegas so that it becomes predominantly toroidal by the time β reachesorder unity. The leading interpretation of the GRAVITY observa-tions of astrometric motion of the IR emission during Sgr A* flares(Gravity Collaboration et al. 2018) requires that the horizon scalemagnetic field be mostly perpendicular to the angular momentumof the gas. We find qualitative agreement with this result in oursimulations that have more magnetized winds. A more quantita-tive comparison to the observations using full radiative transfer inGRMHD simulations using such a field as initial conditions willrequire additional work.Cuadra, Nayakshin & Martins (2008) found that the windsof only 3 WR stars (E20 / IRS 16C, E23 / IRS 16SW, and E39 / IRS16NE) dominated the t = ’ orbital configuration. This is both because ofthe proximity and relatively slow wind speeds ( ∼
600 km / s) of thesewinds. Since we adopted the ‘ ’ configuration from Cuadra,Nayakshin & Martins (2008) with only slight changes, it is not sur-prising that the same three stellar winds seem to be the most impor-tant for determining the properties of the innermost accretion flowin our simulations (e.g., Figures 10 and 11). Future observationsthat place stronger constraints on the mass-loss rates, wind speeds,and especially the magnetic field strengths of these stars would thusgo a long way towards reducing the uncertainty in our calculation.Several observations suggest that gas surrounding Sgr A* isaligned with the clockwise stellar disc both near the horizon andjust inside the Bondi radius (Wang et al. 2013; Bower et al. 2014;Psaltis et al. 2015; though see also Gravity Collaboration et al.2018). Our simulations are consistent with this result for β w (cid:62) (Figure 9) but not for β w =
10 due to wind collimation altering thedistribution of angular momentum in the winds (not shown). If alarge fraction of the accreting gas (and associated magnetic field) Unlike the particle based calculation of Cuadra, Nayakshin & Martins(2008), we do not have a rigorous way to track the gas from each individ-ual wind in our current implementation. We can only infer which stellarwinds dominate the accretion budget from, e.g., the poloidal and toroidalanimations. c (cid:13) , 000–000 eeding Sagittarius A* were sourced from material outside the region where the major-ity of the WR stars reside, then it would also be unlikely for itsangular momentum to coincide with the stellar disc. In all of oursimulations, the direction of the angular momentum of the inneraccretion flow is not strictly constant in time over the ∼ must be tilted with respect to the spin of the black hole at least mod-erately often since the time scale for the spin of the black hole tochange is much longer than 1000 yr. Simulations of tilted accre-tion discs (like those of Fragile & Anninos 2005; Liska et al. 2017;Hawley & Krolik 2019) are thus likely necessary for horizon scalemodeling of Sgr A*.Our results could have a significant impact on current state ofthe art models of horizon scale accretion onto Sgr A*. GRMHDsimulations to date almost universally rely on the MRI as themechanism to drive accretion. It is not clear how much the re-sults of these simulations and their observational consequencesmight change using the dynamically di ff erent flow structure foundhere. For instance, if the disc is less turbulent without the MRI,how does this e ff ect the time-variability properties of the emission?Would nearly empty, magnetically dominated jets still be robustlypresent in GRMHD and does this depend on black hole spin andhorizon-scale flux in the same way as in current simulations (e.g.Tchekhovskoy, Narayan & McKinney 2011)? Such questions andmore will be important to answer in order to further our understand-ing of the emission from Sgr A* and other low luminosity AGN. ACKNOWLEDGMENTS
We thank the referee, Ramesh Narayan for a careful and thoroughreading of this manuscript. We thank F. Yusef-Zadeh for usefuldiscussions, as well as all the members of the horizon collabo-ration, http: // horizon.astro.illinois.edu. SR was supported by theGordon and Betty Moore Foundation through Grant GBMF7392for part of the duration of this work. EQ thanks the Princeton As-trophysical Sciences department and the theoretical astrophysicsgroup and Moore Distinguished Scholar program at Caltech fortheir hospitality and support. This work was supported in part byNSF grants AST 13-33612, AST 1715054, AST-1715277, Chan-dra theory grant TM7-18006X from the Smithsonian Institution, aSimons Investigator award from the Simons Foundation, and bythe NSF through an XSEDE computational time allocation TG-AST090038 on SDSC Comet. This work was made possible bycomputing time granted by UCB on the Savio cluster.
REFERENCES
Abramowicz M., Jaroszynski M., Sikora M., 1978, A & A, 63,221Agol E., 2000, ApJL, 538, L121Bagano ff F. K. et al., 2003, ApJ, 591, 891Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214Blandford R. D., Begelman M. C., 1999, MNRAS, 303, L1Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433Bower G. C. et al., 2014, ApJ, 790, 1Bower G. C., Wright M. C. H., Falcke H., Backer D. C., 2003,ApJ, 588, 331Bower G. C., et al., 2018 Calder´on D., Cuadra J., Schartmann M., Burkert A., Russell C.M. P., 2019, arXiv e-prints, arXiv:1910.06976Chan C.-K., Psaltis D., ¨Ozel F., Narayan R., Sad¸owski A., 2015,ApJ, 799, 1Cuadra 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, 358Cuadra J., Nayakshin S., Wang Q. D., 2015, MNRAS, 450, 277De Villiers J.-P., Hawley J. F., 2003, ApJ, 589, 458Desvignes G. et al., 2018, ApJL, 852, L12Doeleman S. et al., 2009, in Astronomy, Vol. 2010, astro2010:The Astronomy and Astrophysics Decadal Survey, p. 68Eatough R. P. et al., 2013, Nature, 501, 391Einfeldt B., 1988, SIAM Journal on Numerical Analysis, 25, 294Event Horizon Telescope Collaboration et al., 2019a, ApJL, 875,L1Event Horizon Telescope Collaboration et al., 2019b, ApJL, 875,L2Event Horizon Telescope Collaboration et al., 2019c, ApJL, 875,L5Fishbone L. G., Moncrief V., 1976, ApJ, 207, 962Fragile P. C., Anninos P., 2005, ApJ, 623, 347Gammie C. F., McKinney J. C., T´oth G., 2003, ApJ, 589, 444Gillessen S. et al., 2017, ApJ, 837, 30Gillessen S. et al., 2019, ApJ, 871, 126Gravity Collaboration et al., 2018, A&A, 618, L10Hawley J. F., Krolik J. H., 2006, ApJ, 641, 103Hawley J. F., Krolik J. H., 2019, ApJ, 878, 149Igumenshchev I. V., Narayan R., 2002, ApJ, 566, 137Igumenshchev I. V., Narayan R., Abramowicz M. A., 2003, ApJ,592, 1042Janiuk A., Proga D., Kurosawa R., 2008, ApJ, 681, 58Jiang Y.-F., Stone J., Davis S. W., 2017, arXiv e-prints,arXiv:1709.02845Kozlowski M., Jaroszynski M., Abramowicz M. A., 1978, A& A,63, 209Lemaster M. N., Stone J. M., Gardiner T. A., 2007, ApJ, 662, 582Liska M., Hesp C., Tchekhovskoy A., Ingram A., van der Klis M.,Marko ff S., 2017, ArXiv e-printsLu J. R., Ghez A. M., Hornstein S. D., Morris M. R., BecklinE. E., Matthews K., 2009, ApJ, 690, 1463Lu W., Bonnerot C., 2019, arXiv e-prints, arXiv:1904.12018Marrone D. P., Moran J. M., Zhao J.-H., Rao R., 2007, ApJL, 654,L57Martins F., Genzel R., Hillier D. J., Eisenhauer F., Paumard T.,Gillessen S., Ott T., Trippe S., 2007, A & A, 468, 233McKinney J. C., Gammie C. F., 2004, ApJ, 611, 977Mishra B., Begelman M. C., Armitage P. J., Simon J. B., 2019,arXiv e-prints, arXiv:1907.08995Mo´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., S ¨A dowski A., Penna R. F., Kulkarni A. K., 2012,MNRAS, 426, 3241 c (cid:13)000
Abramowicz M., Jaroszynski M., Sikora M., 1978, A & A, 63,221Agol E., 2000, ApJL, 538, L121Bagano ff F. K. et al., 2003, ApJ, 591, 891Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214Blandford R. D., Begelman M. C., 1999, MNRAS, 303, L1Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433Bower G. C. et al., 2014, ApJ, 790, 1Bower G. C., Wright M. C. H., Falcke H., Backer D. C., 2003,ApJ, 588, 331Bower G. C., et al., 2018 Calder´on D., Cuadra J., Schartmann M., Burkert A., Russell C.M. P., 2019, arXiv e-prints, arXiv:1910.06976Chan C.-K., Psaltis D., ¨Ozel F., Narayan R., Sad¸owski A., 2015,ApJ, 799, 1Cuadra 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, 358Cuadra J., Nayakshin S., Wang Q. D., 2015, MNRAS, 450, 277De Villiers J.-P., Hawley J. F., 2003, ApJ, 589, 458Desvignes G. et al., 2018, ApJL, 852, L12Doeleman S. et al., 2009, in Astronomy, Vol. 2010, astro2010:The Astronomy and Astrophysics Decadal Survey, p. 68Eatough R. P. et al., 2013, Nature, 501, 391Einfeldt B., 1988, SIAM Journal on Numerical Analysis, 25, 294Event Horizon Telescope Collaboration et al., 2019a, ApJL, 875,L1Event Horizon Telescope Collaboration et al., 2019b, ApJL, 875,L2Event Horizon Telescope Collaboration et al., 2019c, ApJL, 875,L5Fishbone L. G., Moncrief V., 1976, ApJ, 207, 962Fragile P. C., Anninos P., 2005, ApJ, 623, 347Gammie C. F., McKinney J. C., T´oth G., 2003, ApJ, 589, 444Gillessen S. et al., 2017, ApJ, 837, 30Gillessen S. et al., 2019, ApJ, 871, 126Gravity Collaboration et al., 2018, A&A, 618, L10Hawley J. F., Krolik J. H., 2006, ApJ, 641, 103Hawley J. F., Krolik J. H., 2019, ApJ, 878, 149Igumenshchev I. V., Narayan R., 2002, ApJ, 566, 137Igumenshchev I. V., Narayan R., Abramowicz M. A., 2003, ApJ,592, 1042Janiuk A., Proga D., Kurosawa R., 2008, ApJ, 681, 58Jiang Y.-F., Stone J., Davis S. W., 2017, arXiv e-prints,arXiv:1709.02845Kozlowski M., Jaroszynski M., Abramowicz M. A., 1978, A& A,63, 209Lemaster M. N., Stone J. M., Gardiner T. A., 2007, ApJ, 662, 582Liska M., Hesp C., Tchekhovskoy A., Ingram A., van der Klis M.,Marko ff S., 2017, ArXiv e-printsLu J. R., Ghez A. M., Hornstein S. D., Morris M. R., BecklinE. E., Matthews K., 2009, ApJ, 690, 1463Lu W., Bonnerot C., 2019, arXiv e-prints, arXiv:1904.12018Marrone D. P., Moran J. M., Zhao J.-H., Rao R., 2007, ApJL, 654,L57Martins F., Genzel R., Hillier D. J., Eisenhauer F., Paumard T.,Gillessen S., Ott T., Trippe S., 2007, A & A, 468, 233McKinney J. C., Gammie C. F., 2004, ApJ, 611, 977Mishra B., Begelman M. C., Armitage P. J., Simon J. B., 2019,arXiv e-prints, arXiv:1907.08995Mo´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., S ¨A dowski A., Penna R. F., Kulkarni A. K., 2012,MNRAS, 426, 3241 c (cid:13)000 , 000–000 S. M. Ressler, E. Quataert, J. M. Stone
Pang B., Pen U.-L., Matzner C. D., Green S. R., Liebend¨orfer M.,2011, MNRAS, 415, 1228Parker E. N., 1965, Space Sci. Rev., 4, 666Paumard T. et al., 2006, ApJ, 643, 1011Pen 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, arXiv e-prints, arXiv:1904.04923Proga D., Begelman M. C., 2003a, ApJ, 582, 69Proga D., Begelman M. C., 2003b, ApJ, 592, 767Psaltis D., Narayan R., Fish V. L., Broderick A. E., Loeb A.,Doeleman S. S., 2015, ApJ, 798, 15Quataert E., Gruzinov A., 2000a, ApJ, 545, 842Quataert E., Gruzinov A., 2000b, ApJ, 539, 809Ressler S. M., Quataert E., Stone J. M., 2018, MNRASRessler S. M., Quataert E., Stone J. M., 2019, MNRAS, 482, L123Ressler S. M., Tchekhovskoy A., Quataert E., Gammie C. F.,2017, MNRAS, 467, 3604Royster M. J., Yusef-Zadeh F., Wardle M., Kunneriath D., CottonW., Roberts D. A., 2019, ApJ, 872, 2Russell C. M. P., Wang Q. D., Cuadra J., 2017, MNRAS, 464,4958Sakurai T., 1985, A&A, 152, 121Sa¸dowski A., Narayan R., Penna R., Zhu Y., 2013, MNRAS, 436,3856Shakura N. I., Sunyaev R. A., 1973, A & A, 24, 337Shcherbakov R. V., Bagano ff F. K., 2010, ApJ, 716, 504Shvartsman V. F., 1971, Soviet Atron., 15, 377Stone J. M., Gardiner T. A., Teuben P., Hawley J. F., Simon J. B.,2008, ApJS, 178, 137Stone J. M., Pringle J. E., Begelman M. C., 1999, MNRAS, 310,1002Stone N. C., Kesden M., Cheng R. M., van Velzen S., 2019, Gen-eral Relativity and Gravitation, 51, 30Tchekhovskoy A., Narayan R., McKinney J. C., 2011, MNRAS,418, L79Teyssier R., 2002, A & A, 385, 337Wang Q. D. et al., 2013, Science, 341, 981White C. J., Stone J. M., Quataert E., 2019, The AstrophysicalJournal, 874, 168Yuan F., Narayan R., 2014, AR&A, 52, 529Yusef-Zadeh F., Bushouse H., Sch¨odel R., Wardle M., Cotton W.,Roberts D. A., Nogueras-Lara F., Gallego-Cano E., 2015, ApJ,809, 10Zhu Z., Stone J. M., 2018, ApJ, 857, 34
APPENDIX A: RESOLUTION STUDY
We have argued in the main text in § ff ecting our results, we performedtwo additional simulation with β w =
10 and an inner boundary ra-dius of r in ≈ . × − pc. The first simulation was run for 1.025kyr until t = − .
075 kyr with our usual base resolution of 128 cellswith 8 levels of mesh refinement that increase the resolution by afactor of 2 each time the radius decreases by a factor of 2. The sec-ond simulation increased the resolution by a factor of four within ∼ .
007 pc and thus spans many orbital times for the small radiiof interest. Both calculations were done without radiative cooling r (pc)10
4x res.fiducial c s | v r | M (>0) M (<0) Figure A1.
Demonstration of convergence for β w =
10 simulations. Dashedlines are from a simulation at the fiducial resolution while solid lines arefrom a simulation with a factor of 4 higher resolution for r (cid:46) ρ (black), sound speed, c s (blue), radial ve-locity, v r (red), and net accretion rate, ˙ M (green for positive, orange fornegative), are essentially identical after being run for 25 yr (one orbital timeat ≈ .
007 pc) and thus converged at the fiducial resolution used throughoutthis work. in order to reduce the computational cost. Figure A1 shows thatthe resulting radial profiles of angle-averaged fluid quantities in thetwo simulations are nearly identical. Thus we are confident that thegeneral properties of our simulations are well converged and notlimited by resolution.
APPENDIX B: SIMULATIONS WITH LONGER RUNTIMES
In this Appendix we show that our results are not sensitive to thestart time ( t ) of our simulations, that is, the the length of time werun our simulations before the present day at t =
0. In principle, t should be chosen to represent the typical duration of the WR phase( ∼
100 kyr), however, such a simulation would be extremely ex-pensive and perhaps unnecessarily so if the resulting dynamics ofthe accretion flow are mostly sensitive to the t = ∼ ffi cient to study accretion at t =
0. As this manuscriptwas in press, however, Calder´on et al. (2019) (C19) released newresults of conservative, grid-based, hydrodynamic simulations ofaccretion via stellar winds onto Sgr A* that suggest otherwise. Us-ing the
RAMSES code (Teyssier 2002) with the same wind speeds,orbits, and mass-loss rates of the WR stars used in the present work,they studied the e ff ect of starting their simulations further back intime. The results of their control run with t = − . t = − . ∼ t =
0. Presumably the longerthe simulation is run the more massive the disc becomes. Based onthis finding, they conclude that simulations using t = − . c (cid:13) , 000–000 eeding Sagittarius A* In light of this result, we revisited the question of how sen-sitive our results are to the choice of t by running an additional t = − t f = . ∼ ∼ t = − . r in ≈ × − pc. Figure B1 shows the accre-tion rates through r ≈ × − pc for the t = − t = − . t = t = − . t = − t = t = − t further back in time than -1.1kyr. Though we cannot run a simulation going as far back as the ∼ -100 kyr timescale appropriate for the lifetime of a typical WRstar, we would not expect to find any significant di ff erences fromthe t = − t = − . ff erence in our implementation seems to bethe specific treatment of the stellar winds: C19 use the techniquedescribed in Lemaster, Stone & Gardiner (2007) that overwritescertain cells with the analytic solution of a spherical wind whilewe treat the winds as source terms spread over a certain numberof cells (see §
2) that ultimately drive a spherical wind. Overall,our resolution is a factor of ∼ t (kyr) | M | ( M / y r ) t = 1.1 kyr t = 9 kyr Figure B1.
Accretion rate vs. time at r ≈ × − pc for hydrodynamicsimulations with t = − . t = − ff ect on the resulting ac-cretion rate around t =
0. This demonstrates that t = -1.1 kyr is su ffi cientlyfar back in time before the present day to start our simulations.c (cid:13)000