FORGE’D IN FIRE II: THE FORMATION OF MAGNETICALLY-DOMINATED QUASAR ACCRETION DISKS FROM COSMOLOGICAL INITIAL CONDITIONS

In a companion paper, we reported the formation of quasar accretion disks with inflow rates ∼ 10 M ⊙ yr − 1 down to < 300 Schwarzschild radii from cosmological radiation-magneto-thermochemical-hydrodynamical galaxy and star formation simulations. We see the formation of a well-defined, steady-state accretion disk which is stable against star formation at sub-pc scales. The disks are optically thick, with radiative cooling balancing accretion, but with properties that are distinct from those assumed in most previous accretion disk models. The pressure is strongly dominated by (primarily toroidal) magnetic fields, with a plasma 𝛽 ∼ 10 − 4 even in the disk midplane. They are qualitatively distinct from magnetically elevated or arrested disks. The disks are strongly turbulent, with trans-Alfvénic and highly super-sonic turbulence, and balance this via a cooling time that is short compared to the disk dynamical time, and can sustain highly super-Eddington accretion rates. Their surface and 3D densities at ∼ 10 3 − 10 5 gravitational radii are much lower than in a Shakura-Sunyaev disk, with important implications for their thermo-chemistry and stability. We show how the magnetic field strengths and geometries arise from rapid advection of flux with the inflow from much weaker galaxy-scale fields in these “flux-frozen” disks, and how this stabilizes the disk and gives rise to efficient torques. Re-simulating without magnetic fields produces catastrophic fragmentation with a vastly smaller, lower-(cid:164) 𝑀 Shakura-Sunyaev-like disk.

1. INTRODUCTION Accretion disks are important in a wide variety of astrophysical contexts, ranging from supermassive black hole (BH) growth and evolution to star and planet and satellite formation to X-ray binaries and neutron-star mergers.Around supermassive BHs in particular, these disks are believed to be the engine that powers quasars, the most luminous sources in the Universe (Schmidt 1963;Salpeter 1964), as well as less-luminous active galactic nuclei (AGN).As such, they funnel mass at enormous rates even exceeding ≳ 10 M ⊙ yr −1 to the BH, and ultimately provide most of the mass in SMBHs today (Soltan 1982).The radiation, outflows, and jets launched from the inner regions of such disks (Laor et al. 1997;Crenshaw et al. 2000;Dunn et al. 2010;Sturm et al. 2011;Faucher-Giguère & Quataert 2012;Faucher-Giguère et al. 2012;Zakamska et al. 2016;Williams et al. 2017) -collectively "AGN feedback" -are also widely believed to explain (Silk & Rees 1998;King 2003;Di Matteo et al. 2005;Murray et al. 2005;Hopkins et al. 2005a,b;Torrey et al. 2020) the observed remarkable correlations between BH and host galaxy properties (Magorrian et al. 1998;Ferrarese & Merritt 2000;Gebhardt et al. 2000;Hopkins et al. 2007b;Aller & Richstone 2007;Kormendy et al. 2011) and to dramatically influence galaxy formation and evolution (Croton et al. 2006;Hopkins et al. 2006Hopkins et al. , 2008;;Wellons et al. 2022;Mercedes-Feliz et al. 2023;Cochrane et al. 2023).Understanding the nature, origins, and dynamics of quasar accretion disks, therefore, remains a crucial challenge in theoretical astrophysics.
Since the seminal work by Shakura & Sunyaev (1973, SS73) (and others like Novikov & Thorne 1973), much of the work on quasar accretion disks has assumed as a starting point some variation of the "SS73 -disk" model: this takes disks to be geometrically thin (height  ≪ ), optically thick (black-body like), sub-sonically turbulent (sonic Mach number M  < 1), slowly cooling ( cool ≫  dyn = 1/Ω), thermal-pressure-dominated (plasma  ≡  2  / 2  > 1), radiatively efficient, well-ionized, and parameterized by an effectively constant- viscosity, where  ∼ () 2 / 2  < 1 represents some Maxwell or Reynolds stresses and the kinematic viscosity scales as  ≡  2  /Ω.Numerous variations have been introduced, including e.g.radiatively inefficient and/or advection-dominated, optically thin disks believed to be relevant for very low accretion rates (e.g.Narayan & Yi 1995); radiatively inefficient but still optically-thick "slim" disks at super-Eddington accretion rates (Paczyńsky & Wiita 1980); magnetically "elevated" disks with upper atmospheres (at multiple scale-heights above the midplane) or coronae with  |z| ≫H ≲ 1 (Miller & Stone 2000); magnetically "arrested" disks where magnetic pressure halts accretion near the in-nermost stable circular orbit (ISCO) at low accretion rates (Tchekhovskoy et al. 2011); gravitoturbulent disks relevant at low values of Toomre  (Gammie 2001); and many others (for reviews, see e.g.Pringle 1981;Frank et al. 2002;Abramowicz & Fragile 2013;Jafari 2019).Yet for typical quasars accreting around the Eddington limit, some form of SS73-like -model (whether "thin" or "slim" disk in flavor) is still overwhelmingly the "default" model of reference.Likewise, it is usually assumed (for typical quasars) that the effective viscosity in the disk is dominated by a combination of Maxwell and Reynolds stresses produced by the weak-field ( ≫ 1) magneto-rotational instability (MRI; Balbus & Hawley 1998).
This leaves some crucial questions unresolved, however.For one, it has been known for decades that if one simply extrapolates an SS73 disk to large radii with quasar-level luminosities, then outside a few hundred gravitational radii (∼ 10 −4 − 10 −3 pc, much smaller than typical ISM or even obscuring "torus" scales), it would naively become gravitationally unstable and should rapidly fragment rather than fueling the BH (Shlosman & Begelman 1989;Shlosman et al. 1990;Goodman 2003).Moreover, the properties of the disks in the models above (including both analytic models and traditional "accretion disk" simulations which can only extend to some modest number of gravitational radii from the BH), and even "which type of disk" one actually has, depend fundamentally on the "outer boundary conditions" set by larger-scale inflows into the accretion disk region.Most notably, the accretion rate  itself is simply taken as some constant input, and this has a major effect on the qualitative properties of the disks in the models above.But even for a fixed , one can imagine different distributions of angular momentum of inflowing material, which can produce qualitatively distinct phenomena (including e.g.warps, precessions, flips, or dynamical instabilities, if not highly coherent and close-to-circular; see Scheuer & Feiler 1996;Nayakshin 2005;Hobbs & Nayakshin 2009).And both the magnetic flux and geometry of the magnetic fields (e.g.primarily toroidal or poloidal, tangled or coherent) generally must be assumed.Likewise, many other possible boundary condition effects are often ignored -for example, the effects of global gravitational modes (e.g.coherent eccentric/lopsided disk modes) sourced by external perturbations or collective effects of stars at larger radii outside the disk (Hopkins & Quataert 2010a,b;Anglés-Alcázar et al. 2021).As a result, historical simulations of "strongly magnetized disks" (for example Gaburov et al. 2012;Forgan et al. 2017;Ju et al. 2017;White et al. 2019;Mishra et al. 2020;Kudoh et al. 2020), while crucial for understanding the internal evolution, structure, dynamics, and variability of such disks, must adopt critical parameters like the magnetic flux ad-hoc in their initial conditions and so cannot answer the question of whether or not such disks should or even could arise in real quasars and AGN.Doing so would require a self-consistent predictive model that follows the gas flows and magnetic field dynamics, star formation and feedback on much larger (galactic and intergalactic) scales, all the way down to the BH accretion disk scales.
Motivated by this, in Hopkins et al. (2023a) (henceforth Paper I) we presented the first simulations to follow all of these physical and dynamical effects in a single simulation around a SMBH, from cosmological simulations (using a super-Lagrangian hyper-refinement technique) down to scales of ∼ 80 au (less than ∼ 300  schw ).We observed the formation of a true ( ≫ 1) "accretion disk."In Paper I, we focused on the hierarchy of processes driving angular momentum loss and gas inflow from scales as large as ≳ Mpc onto the galaxy, through the galactic nucleus, the BH radius of influence (BHROI), torus-like regions, all the way down to the accretiondisk scales.We also studied how turbulence, magnetic fields, and radiation-hydrodynamics produced star formation on large scales and lead to the suppression of star formation on small scales (≲ 0.1−1 pc) around the SMBH, allowing for the conditions to transition naturally from "galactic-type" or "ISM-like" conditions at ≳ pc scales to "accretion-disk-like" at ≪ pc.
Here, in Paper II of the series, we focus on the emergent properties of the accretion disks in these simulations, and the physics that give rise to their key behaviors.Specifically, we show that the simulations naturally produce disks that are strongly magnetically dominated ( ≪ 1, with values much smaller in the midplane than usually assumed in historical models), specifically dominated by a toroidal magnetic field (but with substantial "turbulent," radial, and poloidal field components), with vigorous trans-Alfvénic, highly supersonic turbulence, large coherent eccentricities and coherent global modes, as well as gravito-turbulence and spiral arm-like structures.All together these produce rapid radiatively efficient and potentially super-Eddington accretion.We show that the fields are amplified by simple flux-freezing -or similarly, that the toroidal field dynamo is "closed" by rapid advection of new magnetic flux -with completely "normal" ISM magnetic field strengths (themselves built up from extremely small trace cosmological fields).We therefore refer to them as "fluxfrozen disks" for simplicity.In a companion paper (Hopkins 2023a; hereafter Paper III), we further demonstrate that these ideas are supported by a simple analytic accretion disk model.As such, while this is just one simulation, we might expect these behaviors to be quite common.Moreover, we show that the local turbulence may not, in fact, be dominated by the traditional weak-field MRI, but perhaps by distinct instabilities or variants of the MRI that arise when the magnetic fields are extremely strong.
In § 2 we summarize the numerical methods, and Table 1 defines some useful variables we will refer to throughout.In § 3 we summarize the basic conditions and properties of the ISM predicted on sub-pc scales to set the context here, including the connection to large radii in § 3.1, basic gas properties in § 3.2, and (lack of) star formation in § 3.3.In § 4 we examine the magnetic structure of the disks in more detail, discussing both the strength and detailed structure ( § 4.1) and physical origins ( § 4.2) of the strong magnetic fields.In § 5 we consider the same for the velocity field structure in the midplane ( § 5.1) and out-of-plane ( § 5.2), its relation to global coherent eccentric/lopsided disk modes ( § 5.3) and the details of the turbulent structure ( § 5.4), its physical origins/driving ( § 5.5), and the (relatively weak) role of turbulent resistivity ( § 5.6).In § 6 we discuss the vertical structure and profiles of various thermo-chemical and magnetic disk properties ( § 6.1) and their (weak) stratification ( § 6.2) as well as the physics behind this ( § 6.3).In § 7, we explore the physical torques and angular momentum exchange processes in the disk ( § 7.1) and their relation to different stresses including traditional Maxwell and Reynolds stresses ( § 7.2).With this in mind, § 8 compares to a variety of previous literature models of magnetized accretion disks, including magnetically arrested ( § 8.1), magnetically elevated ( § 8.2), galactic/star-forming ( § 8.3), toroidal-field dominated ( § 8.4), and decaying ( § 8.5) disks.We briefly describe how the disks are likely to be mis-aligned with the pre-existing BH spin in § 9.In § 10, we contrast a simulation without magnetic fields and describe how this leads to runaway nuclear star formation ( § 10.1), ordersof-magnitude lower accretion rates ( § 10.2), and a razor-thin, much smaller and lower-mass gravitoturbulent  disk ( § 10.3).We summarize and conclude in § 11.

METHODS 2.1. Physics & Resolution
The simulations studied here are presented and extensively described in Paper I. Briefly, we begin from a ∼ (100 cMpc) 3 cosmological periodic box at redshift  ∼ 100 with a primordial trace magnetic field, and follow it as a cosmological galaxy formation simulation following the full combined Feedback In Realistic Environments (FIRE; Hopkins et al. 2018Hopkins et al. , 2023b) ) and STARFORGE (Grudić et al. 2021;Guszejnov et al. 2021) physics in the code GIZMO1 (Hopkins 2015).We evolve the simulation with a modest refinement (target mass resolution ∼ 4000 M ⊙ in the galaxy, or spatial resolution ∼ 10 pc in the galaxy nucleus) until a redshift  ∼ 4.4 when a period of violent merging and starburst activity induces large inflows into the central ∼ kpc of the galaxy.We then initiate an additional hyper-refinement layer (as in e.g.Anglés-Alcázar et al. (2021)) to go to higher and higher resolution following the gas inflows, to reach sufficient resolution to resolve individual (proto)star formation, accretion and evolution and protostellar disk structure in the central ≲ 10 − 100 pc of the galaxy.We continue to refine to a target resolution of Δ < 0.01 M ⊙ in the central ≲ 10 pc, to follow gas inflows and disk formation down to < 300 Schwarzschild radii around the super-massive black hole of mass ∼ 1.3 × 10 7 M ⊙ .Figs. 1-2 show some illustrative images of the circum-BH disk which forms in the fiducial simulation.
The simulations include a wide range of physics including magnetic fields (see Fig. 3), using the high-order constrainedgradient method from Hopkins & Raives 2016; Hopkins 2016, with kinetic (anisotropic Braginskii viscosity and conduction) and non-ideal (Ohmic, ambipolar, Hall) effects (Su et al. 2017;Hopkins 2017); cosmic ray transport and coupling to gas dynamics (Hopkins et al. 2022b;Hopkins 2022;Hopkins et al. 2022c,d); self-gravity with adaptive softenings scaling with the resolution and high-order Hermite integrators capable of accurately integrating ≳ 10 5 orbits in hard binaries (Grudić & Hopkins 2020;Grudić et al. 2021;Grudić 2021;Hopkins et al. 2022a); metal enrichment and dust destruction/sublimation (Ma et al. 2017;Gandhi et al. 2022;Choban et al. 2022); super-massive black hole seed formation and growth via gravitational capture of gas (Hopkins et al. 2016;Shi et al. 2022;Wellons et al. 2022); (proto)star formation and accretion and explicit feedback from stars in the form of protostellar jets, main-sequence stellar mass-loss, radiation, and supernovae (Grudić et al. 2022;Guszejnov et al. 2022b,c,a).They include explicit multi-band M1 radiation-hydrodynamics with adaptive-wavelength bands (Hopkins et al. 2020a;Hopkins & Grudić 2019;Grudić et al. 2021) coupled explicitly to all the thermo-chemical processes, radiative cooling and thermo-chemistry incorporating cosmic backgrounds, radiation from local stars, re-radiated cooling radiation, dust, molecular, atomic, metal-line, and ionized species opacities and processes, cosmic rays, and other processes allowing us to robustly model the thermochemistry and opacities in gas with densities from  ∼ 10 −8 − 10 16 cm −3 and temperatures ∼ 1 − 10 10 K in a range of radiation and cosmic ray environments.A pure inflow/accretion boundary is enforced at ≈ 80 au from the central SMBH -we do not model any flux from e.g.jets or radiation emerging from the inner region, as these should depend on the accretion disk properties themselves.Fig. 2 illustrates some of the complex phase structure that emerges even in just the nuclear regions.
We stress that the entire simulation uses the identical, full physics -there is no discontinuous change in the equations integrated in space nor time.Instead, as described in Paper I, we evolve the full self-gravitating radiation-magnetothermochemical-hydrodynamics for all gas cells, and simply allow the code to form two distinct types of star particles: (1) FIRE "stellar population" particles which form from starforming gas in the low-resolution cells (resolution ≫ 1 M ⊙ ), and therefore sample an assumed stellar initial mass function -Image of the gas surface density in our fiducial cosmological simulation.We show projected gas density on a logarithmic scale (increasing dark-to-bright, dynamic range rescaled in each panel from a median   ∼ 10 19 cm −2 at the largest scale to ∼ 10 8 times larger at the smallest scale).Multiple scales are shown to illustrate the dynamic range of the simulation, which zooms in down to ≈ 80 au scales around a ∼ 10 7  ⊙ SMBH in the center of a massive galaxy at redshift  ≈ 4.4 in a ∼ (100 Mpc) 3 cosmological box.The simulations include explicit multi-band radiation-magnetohydrodynamics, detailed thermochemistry/cooling, self-gravity with resolved individual (proto)star formation, accretion, evolution, and feedback, and many other physical processes ( § 2).Tidally captured gas streams from an encounter with a massive star-forming cloud complex triggered via gravitational torques in a galaxy-scale merger fall into the BH radius of influence (BHROI) at a few ∼ pc and circularize at ∼ 0.1 − 1 pc to form an accretion disk which we follow down to ∼ 300 BH Schwarzschild radii.
(IMF) and calculate IMF-integrated rates for stellar feedback; and (2) STARFORGE "individual star" particles which form in the high-resolution cells (≪ 1 M ⊙ ) and therefore evolve along individual (proto)stellar+main sequence evolutionary tracks.
Some of these physics are not important on the scales we will study here, although they may play a crucial role in determining the boundary conditions via their role on larger scales.On all scales we study in detail in this paper, the refinement has reached the target resolution of < 0.01 M ⊙ (we briefly re-ran with refinement a factor ∼ 8 higher, and see no difference in our results), and in the densest regions just outside our inner boundary condition we reach local spatial resolution ∼ 10 − 20 au and time resolution as small as ∼ days.The most relevant physics on these scales are gravity, (ideal) MHD, and explicit radiation-thermodynamics.

Analysis and Definitions
Table 1 defines a number of variables we use throughout.In this manuscript, we will often refer to cylindrical radial/azimuthal/vertical coordinates , , , defined with respect to the angular momentum vector of the inner accretion disk (e.g.gas at  < 0.01 pc) and centered on the SMBH, so ẑ points along the angular momentum vector and φ points in the rotation direction.We distinguish the spherical radius/distance  (with spherical radial/polar/azimuthal angles , , φ defined in the same way so that  =  cos ) from the cylindrical radius .Given our focus, we will generally use the terms toroidal and azimuthal interchangeably (and likewise for poloidal and vertical).
The instantaneous values of fluctuating quantities like B(x, ) ≡ B(x, ) − ⟨B(x, )⟩ are defined by respect to their appropriately weighted averages ⟨B⟩ ≡ ⟨B(x, )⟩ x ≡ ( ∫  B(x, )  3 x)/( ∫   3 x) = (  B    Vol  )/(    Vol  ) at a given time (with  the chosen weight, and the summation over all cells ).For example, unless otherwise specified we will define volume-weighted averages in radial annuli, i.e.  = 1 for  or  (whichever is plotted) within some narrow logarithmic radial annulus of width ∼ 0.1 dex, and  = 0 outside the annulus (though we have checked that the exact choice of bin widths makes no appreciable differences to any plot here).We will sometimes compare mass-weighted ( = ) or other explicitly weighted distributions, where stated.We also define the corresponding weighted 90% (5 − 95%) or ±1  (16 − 84%) inclusion intervals, as e.g. the values of B above/below which correspond to the appropriate fraction of the total weight ( ∫   3 x).We note below that the stress tensor  can be written in terms of a mean and fluctuating component, e.g.⟨ kin   ⟩ = ⟨     ⟩ + ⟨     ⟩, so will specify throughout whether we refer to the "total" or "fluctuating" components.
As described below, we consider a few different definitions of the (gas mass) scale height  and show they give nearly identical results, including measuring the mass-weighted me-dian || in annuli, the mass-weighted rms ⟨ 2 ⟩ 1/2 , or fitting a vertical Gaussian or sech 2 profile to the density.Projected quantities like Σ gas are defined as the sum within cylindrical annuli.For vector/tensor quantities like B, v, , we follow usual convention and define them in terms of cylindrical components.However, we have remade all salient plots of these quantities, instead (1) defining in terms of the spherical components (e.g.⟨  ,  , φ ⟩ instead of  , ,  ); (2) plotting versus spherical radius  instead of cylindrical ; (3) allowing the ẑ axis to vary in annuli (defining it with respect to the angular momentum axis in each annulus instead of a global, fixed coordinate system); and (4) re-defining them in eccentric annuli or similarly subtracting the mean  = 1 non-axisymmetric component (defined and plotted below) from the fluctuating components.Unless we explicitly state otherwise for specific quantities in the text below, these choices do not qualitatively change any of our conclusions or comparisons.
Because we are interested in various non-equilibrium properties and dynamics and their time evolution, we will focus on specific representative times chosen after the simulation reaches steady-state, to represent the range of instantaneous behaviors.However we have surveyed hundreds of snapshots of the simulations and confirm that quantities studied here such as the accretion rates and mass profiles are representative of the range over all times after the simulation reaches its highest refinement level.We explicitly show this for several properties (showing their time evolution) below.For radial profiles of positive-definite quantities like  or ⟨|B| 2 ⟩, we show in Paper I and Paper III that their time-averaged behaviors over the en- We measure ⟨ |B| 2 ⟩ 1/2 in spherical annuli from our inner boundary at ∼ 80 au to > Mpc scales, and show the corresponding ⟨ 2  ⟩ 1/2 for the radial, toroidal, and poloidal components of the field (with directions defined relative to the angular momentum vector of the BH accretion disk at  < 0.1 pc).Inside ≪ 10 pc, the resolution is uniformly < 0.01  ⊙ and star formation follows the STARFORGE individual-resolved-stars physics.Outside of this radius the resolution is lower and the FIRE prescriptions form stellar population particles representing multiple stars sampling an assumed IMF (see § 2).Intergalactic sub-nG magnetic fields are amplified by the turbulent and galactic dynamo to "typical" ISM magnetic field strengths of ∼ 1 − 15 G at ∼ 1 − 10 kpc (labeled), but the fields are amplified and approach kG at  ∼ 100 au.At most radii the fields are roughly isotropic (representing the dis-ordered ISM and CGM without a preferred direction), with a mild radial bias in the infall region, before becoming predominantly toroidal where the clear rotating thin BH accretion disk forms (see Fig. 1).
tire simulation lie well within the scatter we plot at a given time and radius.For many signed quantities, like the toroidal and radial magnetic fields, time-averaging the profiles would artificially obscure the interesting behaviors.

BASIC PROPERTIES & CONDITIONS ON SUB-PC SCALES 3.1. Morphology and Connection to Larger Scales
Figs. 1-2 illustrate the nuclear disk, which forms from cosmological initial conditions.We see the disk forms within a chaotic, clumpy, massive high-redshift ( ≈ 4.5) starburst (galaxy-integrated SFR > 100 M ⊙ yr −1 ) galaxy, where a massive star-forming cloud complex has a close passage to the ∼ 10 7 M ⊙ SMBH in the galaxy nucleus.This outer galaxy is studied in Paper I. Some of the material from the cloud is tidally stripped by the SMBH around its radius of influence (BHROI, a few pc, interior to which the BH dominates the potential) and falls in initially in a radial stream (a tidal tail) and we see it circularize at ∼ 0.1 − 1 pc.As expected given the large cloud size and inhomogeneous structure, not all the inflowing gas has an identical impact parameter, so some falls in at slightly different angles (giving rise to warps and eccentricity in the outer disk).
Fig. 3 illustrates the scaling of the magnetic fields with BHcentric radius, showing the full dynamic range of the simulation.We will study scales within the disk below but, because we will argue below that the magnetic flux carried into the disk is important, we wish to highlight that the fields grow ultimately from sub-nG intergalactic magnetic fields, with ISM magnetic fields on scales ∼ 1 − 10 kpc, which have completely "typical" values of ∼ 1 − 15 G (Beck 2015), similar to those observed in the local ISM of the Milky Way (despite this be-ing a massive, high-redshift galaxy).The magnetic fields (and other properties like gas densities) grow relatively smoothly down to disk scales, without a sharp discontinuity at some particular radius.We also identify the range of scales where our simulation reaches maximum target resolution ( ≲ 10 pc), to make it clear that the entire dynamic range we study here uniformly has mass resolution < 0.01 M ⊙ .
Fig. 4 visualizes several different properties in a projected 2D image over a dynamic range of a factor ≳ 1000 in radius, now focusing just on the innermost regions of our simulation from 10 −3 − 1 pc where the disk forms.These include: 3D gas density  gas , temperature , plasma  ≡  2  / 2  , ratio of toroidal to total field strength | tor |/|B|, ratio of azimuthal to total velocity |  |/|v|, ratio of magnetic pressure  B to total non-rotational kinetic/turbulent ram pressure  turb , and "effective" total local Toomre  eff parameter (including ther-mal+magnetic+turbulent support).Fig. 5 complements this, showing the 1D (averaged in concentric radial shells) radial profiles of different properties out to ∼ pc scales.

Radial Trends and Basic Scalings in the Disk
On sub-pc scales, Figs.4-5 show: 1. Inflow is systematically larger than outflow, and star formation is largely negligible (note the three can coexist, as the medium is clearly not spherically symmetric, nor in strict long-term equilibrium), with an orderof-magnitude constant  in ∼ 20−30 M ⊙ yr −1 sustained into the central ∼ 80 au around the SMBH.
2. The gas densities increase towards the center, with a slightly more shallow profile than a singular isothermal sphere ( ∝  −2 ), with clumpiness evident in the visual projection and in the mass-weighted spherical profile.
3. The radiation and thermal energy densities are in rough equilibrium (with the radiation, dust, and gas kinetic temperatures coming into increasing equilibrium at  ≪ pc, and the dust mostly sublimated at ≲ pc as shown in Paper I), but with both of their energy densities well below the magnetic energy density in the disk, with  ranging from 10 −6 − 10 −2 depending on the local phase.The temperature is dominated by cold and warm atomic media (with some warm ionized gas), and rises weakly towards the center.
4. The kinetic energy of the gas is uniformly larger than magnetic, but most of that on scales ≪ pc is the disk rotation, i.e. |v| ∼   ∼   (), and the disk is rotation-dominated within ≲ 0.1 pc.The remaining non-rotational or "turbulent" kinetic energy density is more comparable to magnetic energy densities, with  B / turb ranging from ∼ 0.01 − 10 depending on the local phase sub-structure and density of the gas.In other words, the turbulence is broadly trans-Alfvénic.
5. The magnetic fields are stronger in the center, with |B| increasing slightly steeper than  −1 , and the toroidal component of the field dominating inside the radii where the disk is ordered and rotation-dominated.
6.A combination of turbulence and magnetic pressure support the measured disk scale height (i.e. the disk is in approximate vertical equilibrium) with / ∼ 0.1 Fig. 4.-Wedge plot showing various disk properties (colormaps labeled) including: 3D gas density  gas , gas temperature , plasma  ≡  2  / 2  , ratio of toroidal to total field strength |  tor |/|B|, ratio of azimuthal to total velocity |  |/|v|, ratio of magnetic pressure  B to total non-rotational kinetic/turbulent ram pressure  turb , and "effective" total (including thermal+magnetic+turbulent support) local Toomre  eff parameter.The plot shows cylindrical coordinates , , but with the cylindrical radius  stretched on a log scale from ∼ 80 au to ∼ 2 pc (labeled), to show the behavior over a large range of scales.All quantities are mass-weighted averages within each image pixel, in a slice through the disk midplane (|  |/ < 0.1).Densities increase inwards; temperatures vary but tend to be mostly cold-to-warm at these scales;  ≪ 1 and decreases in the inner disk  ≲ 0.01 pc; the toroidal field and azimuthal (rotational) velocities dominate inside ≲ 0.1 pc; the magnetic pressure is crudely of order turbulent pressure but varies locally from ∼ 0.01 − 10 times its value; and the effective stability parameter  eff ≫ 1 on all these scales.in the central regions (and quasi-spherical structure at ∼ pc, outside the rotation-dominated region).At all radii ≪ pc, the "effective' Toomre  parameter including thermal+magnetic+turbulent support is large.But even the pure-thermal Toomre  thermal ∼ 10 from  ∼ 10 −2 − 10 −1 pc inside the disk and rises rapidly to  thermal ≫ 10 at  < 0.01 pc.
Given that the accretion rates at ∼ 80 au (our inner boundary) correspond to super-Eddington accretion if they remained constant down to horizon scales with a fixed radiative efficiency of   = 0.1, predicting in detail the quasar luminosities requires radiation-GRMHD simulations which can extend the disk simulations here to those scales (e.g.Jiang et al. 2019).If we assume that the accretion becomes radiatively inefficient (or strong outflows suppress  on near-horizon scales) so that the luminosity remains limited to Eddington (Abramowicz et al. 1988), we would predict  bol ∼  Edd ∼ 10 45 erg s −1 , if on the other hand the radiative efficiency remains high (  ∼ 0.1) and  remains constant to the horizon (essentially an upper limit), we would predict  bol ∼ 10 47 erg s −1 .At the redshift  ≈ 4.5 here, this range brackets the "knee" of the observed bolometric quasar luminosity function (Hopkins et al. 2007a;Shen et al. 2020), so in terms of luminosities, this should (  > 0), and star formation rate (stellar mass formed/accreted within the last dynamical time, in each shell).We see consistent inflow at up to ∼ 20 M ⊙ yr −1 .Top right: Toomre  parameter for the gas alone, separating thermal support (  ), magnetic (  ), turbulent ( turb ) and total ( 2  , eff ≡  2  +  2  +  2 turb ).The system is stable on all these scales, but the relatively modest thermal-only  produces some of the large-scale gravitoturbulence and some (very slow/weak) star formation.Middle left: Gas density (lines show mean, shaded range ∼ 90% inclusion interval) weighted by mass (spikes correspond to denser arms/clumps) or volume.Middle right: Total contribution to the stress/pressure tensor (Frobenius norm of each tensor, lines show volume-weighted mean, shaded range ∼ 90% interval), for kinetic ( v v), magnetic (|B| 2 I/8  − BB/4 ), thermal (   I), and radiation ( ∫ d   D  ) energies.The thermal and radiation terms are in rough equilibrium as expected (the disk is optically thick) but are strongly dominated by magnetic and kinetic terms.Note the kinetic energy is primarily rotation of the disk, the turbulent kinetic energy is order-of-magnitude comparable to magnetic (turbulence is broadly trans-Alfvénic).Bottom left: Scale height / directly measured (median or rms , as labeled), and expected (≈ /  ) from turbulent/magnetic/thermal support.Given the above, we see turbulent+magnetic support clearly dominate the disks' vertical structure.Bottom right: Magnetic field strength (rms volume-weighted value in solid, shaded shows 90% range for total field strength), for total and by toroidal/poloidal/radial component.The fields rise to very large values at small , and are toroidal-dominated in the disk but with non-negligible poloidal+radial terms.In a companion paper (Paper III) we note that all of these scalings can be reasonably approximated with a simple analytic similarity model.correspond to a "typical" quasar at the redshifts simulated.
Because we will refer to it below, we note that Paper III compares these profiles to those predicted for a Shakura & Sunyaev (1973) or SS73 -disk with the same accretion rate , to show they differ by orders of magnitude.Briefly, SS73 and other "weakly-magnetized" models assume  ≳ 1, so the vertical support in the outer disk comes only from thermal pressure (/ ∼   /  ), with an effective viscosity ( ∝    ) provided by some Reynolds/Maxwell stresses with  ∼  2 turb / 2  ∼  2  / 2  ∼ 1/ ≪ 1 leading to inflow (  ∝  Σ gas ).So /,   /  , and  turb /  are predicted by SS73 to be smaller by factors of ∼ 100 − 1000 than the values here, while (for the same2 ) Σ gas would be larger in SS73 by factors ∼ 10 4 − 10 6 .Correspondingly, the midplane density ( ∝ Σ gas /) is larger and effective Toomre  (∝  eff Ω/ Σ gas ∝  Ω 2 /Σ gas ) is smaller by factors ∼ 10 6 −10 8 in SS73 compared to these simulations.

Stability and (Lack of) Star Formation
As noted above, the SFR interior to ≲ pc is much smaller than inflow rates.The physics of this suppression is discussed in Paper I, but as we will see below in our tests without MHD, magnetic fields play an important role (raising , preventing local collapse perpendicular to the mean toroidal field, and promoting faster torques and more rapid accretion).For our purposes here, we show therein and in Hopkins (2023b) (where the properties of the few stars that do form in the outer accretion disk are studied) that the density or total mass of stars on these scales is very small compared to the gas mass, and that the stars contribute negligibly to the dynamics or stresses (momentum/energy flux or turbulence driving) or heating in the disk either via their gravitational influence or via their stellar "feedback" effects (jets, winds, radiation).Indeed, if we re-start the simulations from a snapshot and simply delete all stars at < 1 pc (and disable new star formation at these radii) and run for ∼ 10 − 20 dynamical times, we see no appreciable difference in any of the properties we study on these scales.
So we are justified in neglecting stars and stellar feedback in our discussion of the accretion disk structure, on these scales.This is distinct from larger, more "ISM-like" scales, where stars dominate the mass and dynamics (see e.g.Hopkins & Quataert 2011;Anglés-Alcázar et al. 2021).

STRUCTURE AND ORIGINS OF THE MAGNETIC FIELDS 4.1. Overview of Properties & Transition to Toroidal
Now we turn to study of the magnetic fields in the simulations.Fig. 3 and its "zoomed-in" version, Fig. 5, show that the total magnetic field strength rises to values approaching ∼ kG at ≲ 10 −3 pc, somewhat steeper than |B| ∝  −1 , from ∼ pc scales.We see directly in Fig. 3 that something like ⟨|B|⟩ ∝  −1 also describes the field strength (very crudely) at > pc scales out to ≳ kpc.Defining the poloidal, toroidal, and radial components of the field,3 we see in Figs. 3, 4, & 5 a clear transition from an isotropic or slightly radially-biased field at ≳ 0.1 pc (which Fig. 3 shows is true at much larger radii as well) to a toroidal-dominated field at smaller radii, coinciding with the visually well-ordered disk in Fig. 1.
Fig. 6 plots the projected structure of the field lines in the disk plane (taking a slice through the midplane, so in cylindrical - coordinates considering a wedge of some width in ||/) and edge-on (in cylindrical - coordinates considering a wedge of some width in ), in one example time snapshot within the ordered disk ( ≲ 0.1 pc).Fig. 7 considers the face-on field structure at three different scales  ≲ (0.01, 0.1, 1) pc.Fig. 8 shows the edge-on field-line structure on the smaller scales ≲ 0.01 pc at two different times, to illustrate how they can change, and Fig. 9 shows the field lines overplotted on an edge-on density map to see how illustrate the relation to the density substructure within the disk.
In the face-on projections, we see the fields are fairly well ordered, with a clear transition from more radial field lines pointing along the direction of gas infall onto the nuclear disk, circularizing where the disk forms, to become toroidal, with increasing order in the toroidal field at smaller radii (by eye, it becomes closer to purely azimuthal).We also see clear repeating sign flips in the toroidal field joined by field reversals (occasionally breaking off into large-scale loop-type structures), with some intermediate turbulent zones.However, the inflow continues to be traced even as the toroidal field strengthens as we see a spiral-type structure (i.e.non-negligible coherent radial field components pointing inwards).We also plainly visually see an anti-correlation in the signs of   and   , as expected if the toroidal field is sourced by radial flux.
In the edge-on projections we see less large-scale coherence.In the midplane there is a (much weaker) coherent radial/vertical field component in the - projection, and there is some vaguely "jet like" vertical bipolar field in - (the vertical/conical fields at small  with coherent   and   but oppositely-signed   above/below the disk).Interestingly, looking at the sign of the toroidal field in the edge-on projection, we see that there can be sign flips of  tor at different vertical heights, as well as at different radial intervals; but, these are generally at || ≳  -i.e.outside the body of the disk (with coherence lengths ≳ ).It is worth noting that there is no sign flip of   across the midplane, as is often seen if   were sourced by a much stronger mean poloidal field (although there can be exceptions to this).We see clear evidence for some mode structure in   and   interior to the disk, with wavelength ∼ .
In Fig. 10 we examine the magnetic field profile on these scales somewhat more quantitatively.Here, we consider a 1D profile (averaged in spherical shells), plotting both the mass-weighted mean-field values of the cylindrical , ,  components, as well as their ∼ 1 range.4We plot these both in absolute units, as well as in units of the Alfvén speed  ,  ≡   / √︁ 4 , relative to the circular velocity For each, colors denote the sign of toroidal/azimuthal (  ; top), radial (  ; middle), or poloidal/vertical (  ; bottom) field.Our sign convention is that ẑ and φ point in the direction of the angular momentum vector and direction of co-rotation of the inner gas disk (so   > 0 is prograde), and R points away from the BH.We see large-scale ordered structure tracing the radial inflows from scales outside the BHROI being wrapped around the center as the disk circularizes at ∼ 0.1 pc, with sign changes on large scales determined by the large-scale inflow structure (and reversals or occasional loops forming at the sign change boundaries).There is also an obvious anti-correlation in the sign of   and   .Fig. 7.-Field lines as Fig. 6, in face-on projection of the midplane on three different spatial scales: ±1 pc (left), ±0.1 pc (middle), ±0.01 pc (right).We more clearly see the field lines transition from a somewhat less coherent and radially-biased flow with the filamentary inflow on the largest ∼ 1 pc scales to the more ordered and more clearly azimuthal field geometry on ≪ 0.01 pc scales.We also see the expected   −   anti-correlation on all scales, with a somewhat less correlated/more turbulent   .Fig. 8.-Field lines as Fig. 6, in edge-on projection on smaller ±0.01 pc scales, at two different times (earlier at left, later at right), separated by ∼ 2 × 10 4 dynamical times at the innermost boundary radius.The times are chosen to represent times before (left) and after (right) the sign-flip in   on small scales accretes through our inner boundary at < 80 au, leaving a somewhat larger coherent toroidal field at the late time (but still with successive sign-flips on larger scales).In both cases we see some   −   anti-correlation but at the later time it is less obvious (but would become more obvious if we subtract the coherent mean-field   and plot the fluctuations    ≡   − ⟨   () ⟩ and    ).Coherent modes with wavenumber   ∼   ∼ 1/ are obvious and especially prominent in   at the later time.We denote the  = 0 and  = 0 axes with the thin white dotted line for clarity: we see no evidence for a sign switch in   at  = 0 as would be predicted if the toroidal   were ultimately sourced by a dominant mean-field   .Fig. 9.-Field lines (as Fig. 8; grey) superimposed on an edge-on gas density projection (colors,  as labeled) in - coordinates, in a narrow wedge (sin  < 0.1).There is a large-scale asymmetry in the disk density, reflecting the coherent large-scale eccentric disk (Figs.1-2).Insets show higher-resolution examples of a couple of regions (top-left inset zooms in on the midplane gas just below in the image).There is considerable density substructure, consistent with highly super-sonic and trans-Alfvénic turbulence in the disk, and tangled internal disk fields.We see clear radial/vertical field structure and cells/loops with scale length ∼ .Approaching  → 0, the disk becomes relatively thick, while being thinner at intermediate radii.Eqs.1-4).These successfully reproduce the key behaviors in the simulations.We compare the mass-weighted simulation mean (thick line) value and 90% inclusion range (shaded) of |B| in spherical annuli to the value |B predicted | predicted by the models (labeled), versus radius .We show one early time (top) and one late time (bottom), matching the times in Figs. 8 & 10.These different scalings, which make slightly different assumptions about how the gas is compressed as it accretes, all predict |B| within factors ∼ 2 without a large systematic residual or offset.
c ≡ √︁   enc (< )/ at each radius.We see the same rise in the field values as Fig. 5, and the increasing preference for toroidal fields within the disk as above, but now more quantitatively see the transition from a primarily turbulent field (the rms/dispersion component larger than mean) to more coherent (mean similar to or even larger than dispersion, for   ).The mean ⟨  ⟩ dominates at the inner disk, but the "turbulent" or rms components are not vastly smaller.There appears to be a robust sort of "hierarchy" of the different field components in the inner disk, with the mean toroidal and vertical fields are the strongest and weakest, respectively, with the fluctuating toroidal secondstrongest, followed by the fluctuating radial and vertical fields.This is independent of whether we define these components in a global cylindrical coordinate system (shown), or a spherical coordinate system, or a spatially-variable coordinate system where we rotate each annulus independently to correspond to the angular momentum axis of gas just in that annulus, and/or whether we subtract the mean  = 1 component in each an-   11, but comparing alternative previously proposed analytic models for amplification of |B|, which are not based on advection/freezing of toroidal/radial magnetic flux (see § 4.2.2, scalings (2), (3), (5), and (6)).We see that these do not accurately describe the simulations, especially noting the much larger vertical axis range here compared to Fig. 11.The model used for most previous discussions of toroidal magnetically dominated disks is the predicted saturation value of the linear (local, unstratified) MRI at  Pessah & Psaltis (2005), but we see that the mean fields here are already well above this value at all nuclear radii, the scatter is very large, and there is a systematic trend in the median with radius (so no simple re-normalization leads to agreement).The same is true for constant- models or models with |B| ∝ Σ gas sometimes invoked in the literature.A model where   saturates at ∼  K fares somewhat better but still shows a systemic deviation with radius and a normalization offset from the simulations.nulus (to control for a coherent eccentric mode).We also see that the sign flips in the disk clearly evolve, as mass is accreted through the disk, but the magnitude of the total field stays broadly consistent over tens of thousands of disk dynamical times.
As discussed above ( § 3.2), the magnetic pressure dominates the vertical support of the disk but with O (1) contribution from trans-Alfvénic turbulence: this is reflected both in direct comparison of the turbulent velocity components to the various  ,  below; or from the kinetic energy densities, or relative contribution to the effective stability parameter , or comparison of the disk / to   /  in Fig. 5.In other words,   2 turb ∼ ⟨|B| 2 ⟩ broadly speaking, as we quantify in more detail below.Of course, both   and  turb are much larger in the simulations than in an SS73 disk which assumes  ≫ 1.As expected for a mean-field-dominated disk in steady-state with most of the mass in midplane inflow, this traces a similar power-law |B| ∝  2/3 .Bottom: Variation in |B| with  at a given radius (narrow radial annuli at the different  shown).To compare we normalize |B| and  to their mean within each annulus.Here variation in |B| in an annulus is not dominated by global compression of the mean field, so is not expected to obey the same powerlaw, and instead varies weakly (∝  0.2−0.3 ), as expected if fluctuations are generally modest relative to the mean field.The tail to high /⟨( ) ⟩ is dominated by rare sites of star-formation, so this reflects the physics of local collapse (along field lines).√︃  3 inner /  BH , with  inner ≡ 80 au (the time zero-point is arbitrary).We zoom in to show times a few thousand dynamical times before and after a sign flip event where the innermost gas with ⟨  ⟩ > 0, as evident in e.g.Fig. 6, is accreted, with gas from larger radii with ⟨  ⟩ < 0 moving in towards these smaller radii.We see that the magnitude of both ⟨ |B| 2 ⟩ 1/2 and | ⟨  ⟩ | recover quickly after sign-flip events and remain stable for ≳ 10 4 dynamical times.Bottom: Accretion rate into the central < 80 au, averaged in Δ ∼ 2 yr (∼ 50 Ω −1 inner ) increments.

Origins of the Mean Field in Flux-Freezing/Advection of Flux with Accreting Gas
We find that the qualitative behavior of the dominant mean toroidal field ⟨  ⟩ can be understood primarily from simple flux-freezing considerations.Below, we test and validate this in more detail, but first let us describe the qualitative scenario and key behaviors in the simulations.To begin, recall that the gas forming the disk is tidally captured from a close passage by a molecular cloud complex to the BHROI (Figs. 1-2), so behaves as an initially "cold" (weakly pressurized) tidal filament/stream, akin to satellite galaxy encounters on large scales (Hernquist & Mihos 1995;Bullock & Johnston 2005;Younger et al. 2008;Moster et al. 2010) and similarly analogous to some of the behaviors seen in simulations of magnetic fields in stellar tidal disruption events (TDEs) (Bonnerot et al. -Mean toroidal ⟨  ⟩ in azimuthal annuli (as Fig. 10), versus the expected leading-order term from the induction equation for amplification of the toroidal field via advection/stretching of radial magnetic flux   (⟨  ⟩/ ≈ (    )/ + ...), multiplied by the dynamical time Ω −1 at each radius.We show two times (top and bottom), before and after the sign flip in Fig. 14.We see a more quantitative version of the   −   anti-correlation and that this leading-order average roughly predicts the correct mean toroidal field.Together with the previous figures, this suggests that the strong toroidal field arises primarily from flux conservation (advection of radial flux "closes the dynamo loop").

2017; Guillochon & McCourt 2017).
As shown in Fig. 3, the gas at large radii (e.g. in the sub-kpc galactic nucleus) has broadly isotropic turbulent fields, with magnetic energy density a few percent of the kinetic energy density (see also Fig. 9 in Paper I), as expected for the supersonic dynamo in both idealized (Federrath et al. 2014;Rieder & Teyssier 2017) and multi-phase galaxy formation simulations (Martin-Alvarez et al. 2022;Guszejnov et al. 2020;Seta & Federrath 2022).When some portion of this is captured and falls in (with tangential velocity below  c ), it is tidally stretched into a radial stream of length ℓ and width  (as we see occurring in Figs.6-7).For   and   , the perpendicular areas -(Δ) ( Δ) ∼ ℓ  and (Δ) (Δ) ∼ ℓ , respectivelyare increased by the stretching5 so these components are weakened; in contrast,   is amplified, owing to the perpendicular compression (perpendicular area ( Δ) (Δ) ∼  2 ).Equivalently we can think of the field lines as being "stretched" in the radial direction as they are dragged.The expected amplification of   ranges between  −0 (for e.g. an infalling clump with  ≪ , or no radial motion) to  −2 (maximal case where  ∼  and pure-radial inflow), depending on the efficiency of the perpendicular compression.This agrees with the behavior seen in Figs. 5 & 10.
The infalling gas has non-zero impact parameter  and circularizes at ∼ 0.1 pc (plainly visible in Figs. 2 & 4).The "initially" radial field therefore follows the gas flow and wraps/winds up to become toroidal (Fig. 10).Equivalently the compression ratios rotate as shear now means that the elongation direction is azimuthal, so the mean azimuthal field ⟨  ⟩ is amplified strongly while the mean radial and vertical fields grow less rapidly (remaining sub-dominant to their fluctuating field components).Since the field was initially tangled and isotropic (with |  | ≳ |⟨  ⟩| for all components of B) at radii outside of the BHROI, there are sign flips in the "initial" field which is now stretched into some mean ⟨  ⟩ as it was radially accreted -these become the successive sign flips in the radial direction.Indeed, following the fluid over time in Fig. 8, we note below that the sign flips in   simply reflect these frozen-in trends and are advected inwards with the fluid as it accretes.In steady-state, even if there is some damping of the coherent toroidal field owing to turbulent resistivity or buoyant escape, ⟨  ⟩ is constantly replenished by the steady supply of radial magnetic flux into and through the disk.
For flux-freezing, the compression in the  direction suggests  ∼   ∝ 1/ (Machida et al. 2006).Meanwhile direct analysis of the simulations or simple analytic considerations give weak compression in the  direction.6As shown below, a yet simpler isotropic-flux-freezing expectation  ∝  2/3 works equally well in explaining the evolution of the field strengths both in time (following a Lagrangian parcel) or space (fitting the radial profile).Checking the normalization of |B|, if we fit  ∝  2/3 to the simulations using the midplane values of ⟨|B| 2 ⟩ 1/2 and ⟨ midplane ⟩ (Figs. 3 & 5), we obtain |B| ∼ (2 − 8) G (/cm −3 ) 2/3 (depending on how we weight it; or |B| ∼ (2 − 6) G (/10 kpc) −1 ) -i.e. at Galactic radii this extrapolates to typical mundane values of |B|, consistent with our direct estimates from the simulations in Fig. 3.Moreover, as discussed below and in Paper III, the absolute field strengths here (in Gauss) are actually smaller than in an SS73-like  ∼ 0.1 disk with the same  (even though  ≪ 1 is much smaller).Thus no extreme fields at large radii are needed to sustain these strong fields in the disk.This is consistent with the global field geometry, and at least qualitatively explains the sign flips in the toroidal field in successive radial annuli (as this reflects field reversals from the turbulent fields at much larger radii, before amplification).The idea also explains the lack of systematic sign flips in ⟨  ⟩ at  = 0 (i.e.reversals of ⟨  ⟩ as one vertically crosses the midplane), as well as the broad anti-correlation between   and   .The relative amplification versus suppression above also explains why we see a dominant mean ⟨  ⟩ component with much weaker mean ⟨  ⟩ and ⟨  ⟩ components (Fig. 10; note also that ⟨  ⟩ seems to grow somewhat more slowly  and  2  defined in Pessah & Psaltis (2005) which separate the different regimes of instability labeled (for an analytic, linear, laminar, unstratified system).The bottom-left quadrant corresponds to the classical (linear, unstratified, weak-field) MRI, top-left to a stable regime, while the top-right and bottom-right to the "Type II" and "Type III" (bouyancy-related) instabilities defined therein (or the suprathermal slow mode instability [SSMI] and suprathermal hybrid mode instability [SHMI], respectively, in Das et al. 2018).We see that the volume-filling midplane phases in the simulations reside well into the Type II/III (SSMI/SHMI) regime, with the largest wavelengths (∼ ) broadly similar to the dividing line between those modes, i.e. the toroidal field is stronger than the maximum amplitude at which the linear MRI grows under the usual laminar, unstratified disk assumptions ( max, MRI ,  ∼ √    K , as shown in Fig. 12 above).Colder-phase gas (not shown here) shifts diagonally to the upper right, further into the Type II/SSMI region.
than ⟨  ⟩).Instead, the radial and poloidal fields are more dominated by their turbulent/fluctuating components, whose amplitudes are of order | ,  | ∼ |( ,  /  ) ⟨  ⟩| -i.e.consistent with field lines being stretched, distorted, and perhaps modestly amplified by turbulence within the disk (as we appear to see occurring in Figs.8-9).We stress that, as discussed in the numerical tests in Paper I, microphysical resistivity is not expected to play a significant role here owing to the fact that the gas is still relatively diffuse and, in the nuclear regions of importance, still highly ionized with ionization fraction ≳ 0.01 (so ideal MHD remains a good approximation).This owes to a combination of high temperatures (≳ 1000 K), dust destruction, high cosmic-ray energy densities and ionization rates (even in the Milky Way center, these exceed typical Solar circle values by several orders of magnitude; see Indriolo et al. 2015), and high interstellar radiation field densities due to the concentrated intense star formation (≳ 100 M ⊙ yr −1 in the central ∼ 100 pc; see Paper I).This is unlike accretion of magnetic fields onto a protostar via a protostellar disk, where ionized fractions are expected to be in the range ≲ 10 −17 − 10 −15 in much of the disk.7In § 5.6, we discuss the role of turbulent resistivity in detail and show that it is also sub-dominant: damping of the dominant radial and toroidal fields via turbulent resistivity is generically slower than their growth via flux-freezing and advection, though such damping may be important for the less-coherent poloidal field.Thus our discussion in this section does not assume weak turbulence and indeed we have assumed trans-Alfvénic or even modestly super-Alfvénic (and highly super-sonic) turbulence could be present throughout ( § 5.6), so long as the turbulence is sub-virial ( turb ≲   ), so that the turbulent coherence length is smaller than the characteristic radial distance ≳  over which the field lines are being stretched as part of the tidal inflow/stream.

Flux-Freezing and Advected Flux
We now consider various quantitative tests of the picture described above in § 4.2.1, to validate more rigorously that it is a reasonable description of the simulations.Moreover, one might imagine several alternative scenarios/models that attempt to explain or predict strong magnetic fields in an accretion disk, but we find most of these do not reproduce the behaviors seen in our simulation.Consider the following possibilities: 1. Flux-freezing of the radial/azimuthal magnetic flux, the scenario from § 4.2.1.Here ⟨  ⟩ is sourced by flux-frozen accreted fields, and the mean toroidal field originates from advection of radial flux with the radially-infalling gas from outside the BHROI.Its evolution following the description in § 4.2.1 (Eqs.1-4 below).
One could instead assume that the magnetic field was dominated by a mean poloidal/vertical/dipolar field ⟨  ⟩, which was flux-frozen as mass advected inwards in a laminar disk (e.g.torqued by magnetic braking).Then the fixed fluxto-mass ratio would predict |B pred | ∼ |⟨  ⟩| ∝ Σ gas in the disk.

Traditional MRI-Driven Fields (following Begelman & Pringle 2007
).If the fields (including the mean field) were primarily driven by the traditional weak-field MRI (sourced as usually assumed by some initial poloidal field), then this would predict that the maximum saturation value of the magnetic field strength should be given by the value above which the MRI ceases to grow efficiently.For the analytic models of magnetically-dominated disks in Begelman & Pringle (2007), this is taken to be K (the value derived for linear MRI growth in a laminar, unstratified, local analysis in Pessah & Psaltis 2005).8 4. The Small-Scale Super-Sonic, Rapidly-Cooling Turbulent Dynamo.In the standard super-sonic turbulent dynamo, the fields saturate with magnetic energy a few percent the turbulent kinetic energy (i.e.Alfvén Mach numbers ∼ 3 − 10), with isotropic, tangled fields, and fluctuating components much larger than mean (|B pred | ≫ |B|).
5. The Small-Scale Sub-Sonic, Gravitational/Protostellar Dynamo.Various models for the local dynamo in more slowly Begelman & Armitage 2023).This correction generally makes a small (tens of percent) difference, but even in the annulus where the correction is maximized (near ∼ 0.002 pc where  ln   / ln  is close to −1) the correction never exceeds a factor of ∼ 2.5, so does not change any of our conclusions.cooling dense gas in molecular clouds/clumps/cores collapsing to proto-stellar disks have argued for saturation at fixed  ∼ 1, giving  pred  ∼   (Mocz et al. 2017), again with isotropically tangled fields (|B pred | ≫ |B|).
6. "Arrested" Fields.If the field simply saturated at a value where it would dynamically arrest further inflow (as in e.g."magnetically arrested disks", discussed further below), we might expect  pred  ∼   First, we note that the measurements of the different mean and fluctuating components (e.g. which components are dominant where, and the ratio of mean-to-fluctuating component amplitudes), as well as the presence/absence of different sign flips, shown in Figs.6-10, are all consistent with our favored flux-freezing picture (i), as described in § 4.2.1.The field geometry we see is immediately inconsistent with model (ii) (which assumes the field is dominated by a mean poloidal/vertical field) and models (iii), (iv) and (v) (which all predict that the mean toroidal field should be much smaller than the fluctuating field).Moreover, phenomena such as the sign flips are either not predicted or predicted to have qualitatively different behaviors in models (ii)- (vi).
Second, in Figs.11-12, we specifically compare the measured |B|() in the simulation to the value predicted by the different simple assumptions/models above (at two different times).To begin, in Fig. 11 we compare a group of models for |B| ∼ |⟨  ⟩| motivated by the simple flux-freezing considerations of model (i) as described in § 4.2.1.We see that these reasonably reproduce the absolute magnitude and radial trend of |B| with relatively little scatter.If we compare or (both noted in § 4.2.1), the  pred ∝  −1 scaling clearly exhibits even smaller scatter and more accurate prediction of the mean |B| compared to  pred ∝  −1 .We obtain almost as good a fit with the simpler expression (3) (shown explicitly in Fig. 13).9We can also assume ⟨|  |⟩ ∝ is the scale height set by the magnetic pressure (specifically focusing on the toroidal component providing the vertical support).This gives which also provides a good fit.
To compare, Fig. 12 plots the predicted value of |B pred | from models (ii), (iii), (v), and (vi).We see that all of these models fail to correctly predict the magnitude of the typical fields and their dependence on radius within the disk (not just the field geometries and more detailed structure).Models (ii), (iii), and (v) fare especially poorly, predicting huge variations in the local value of |B| at a given annulus that we do not see in the simulations, predicting the incorrect radial trend (there is a systematic trend in the mean offset from the actual simulation |B|), and predicting the incorrect normalization of |B|.Model (vi) fares somewhat better but is still notably offset from the actual field strengths (and the fact that the disk is actually accreting appears to immediately contradict model (vi)).Model (iv) is not shown here because we have not explicitly separated the turbulent velocity fields/kinetic energies, but we show below ( § 5.4) that the typical Alfvén Mach numbers in the disk at ≲ 0.1 pc are M  ∼ 0.3 − 1 -i.e. the saturation magnetic field strength in the simulation is an order-of-magnitude larger than model (iv) would predict (and, again, the field geometry is completely different).Nonetheless, it is worth noting model (iv) appears perfectly reasonable as a description of both the field geometry and strength at much larger radii  ≫ pc, far outside the disk (in the ISM).Thus we see that the model (i) variants clearly provide a much better fit to the simulated values of |B|, compared to other hypotheses (ii)-(vi) above.
Third, in Fig. 13 we have followed the Lagrangian time evolution of individual fluid elements (since this is a Lagrangian code, this is numerically trivial), and verified that they obey the approximate model-(i) scalings and behaviors in time, as well as in space.This is expected if the disk is in steady state with steady inwards accretion, and ⟨⟩ or  are monotonic functions of , but it is important to validate.It also allows us to confirm that the evolution occurs continuously as the material advects, and not, for example, only after it reaches some critical radius.We can also immediately confirm that the mean toroidal field is not amplified from some trace/turbulent seed field as models (iii), (iv), and (vi) predict.
Fourth, by examining different snapshots, we have verified that the sign flips in   move with Lagrangian fluid elements over time exactly as expected in model (i), as opposed to oscillating as they would if they arose from e.g.instabilities like the MRI (model (iii)) or the turbulent dynamo (models (iv) and (v)).We see this directly in Figs. 8 and Fig. 10, as well as via the fact that when following Lagrangian parcels we do not see sign flips.We illustrate this physics more explicitly in Fig. 14, where we follow the gas at the inner radii of the simulation just before it is accreted.This Figure illustrates a number of important properties: (1) that the magnetic field strength, toroidal field prominence, and accretion rate into < 80 au are stable (to within a factor of a couple) over tens of thousands of dynamical times at our inner radii; (2) that sign flips occur, on a timescale comparable to the accretion timescale (∼  in / gas (< )), as new gas moves from larger radii into the annulus of interest; (3) that the field strength at a given radius is robust to these flips and restores quickly "through" the flip.We also see this in Fig. 15, where at two different times more closely separated than the times in Fig. 10 we more plainly see the sign reversals in ⟨  ⟩ systematically propagating inwards with the gas.As noted above, the form of the sign flips is also qualitatively inconsistent with model (ii).
Fifth, Fig. 15 compares the mean toroidal field ⟨  ⟩ with its expected growth if it was ultimately sourced primarily from advection of radial flux -as predicted by model (i) -in a close-to-Keplerian disk.In particular, we approximate the induction equation for ⟨  ⟩,   ⟨  ⟩ = ⟨  (    −     )⟩ + ⟨  (    −     )⟩, to include only the term that accounts for the stretching of the radial field by the toroidal flow   ⟨  ⟩ ≈   ⟨    ⟩.Note the pure radial   transport term   ⟨    ⟩ is usually smaller, but not always negligible, while the vertical flux divergence     vanishes and the vertical inflow term     is small.In any case we see this approximation works well at describing quantitatively the sign flips and trends in ⟨  ⟩.This behavior is generally distinct from the predictions of the alternative models (ii)-(vi) above.Moreover, if we multiply   ⟨    ⟩ by the characteristic dynamical time  dyn = Ω −1 , this appears to provide a remarkably good order-of-magnitude estimate of the saturation ⟨  ⟩this is expected if either (i) the accretion is dynamical (so the amplification time is limited to some multiple of Ω −1 ) or (ii) if the disk is trans-Alfvénically turbulent (so  ∼  turb /Ω by definition in vertical equilibrium, and the turbulent magnetic dissipation time is ∼ Ω −1 ) or (iii) if the midplane toroidal flux is lost via buoyancy on the vertical buoyancy timescale (a few to tens of times Ω −1 ).In any of these cases, provided the disk support is dominated by a mean toroidal field and Maxwell stresses with trans-Alfvénic turbulence (as we have here), the dimensional expectation for the rate at which flux is "lost" through a fixed Eulerian annulus is similar, and we can think of it as being "replenished" by advection of new ra-dial+toroidal flux with the inflow from larger radii.So we can say -effectively equivalently to our description above -that the dynamo is "closed" by advection of mean radial+toroidal magnetic flux from the inflowing gas.
Together, all of these comparisons quantitatively (and qualitatively) support our argument from § 4.2.1 that the mean field is ultimately sourced via flux-freezing (model (i) here), rather than via some other scenario (e.g.models (ii)-(vi)).

Discriminating Between Models for the "Turbulent" Field
Now we consider the sub-dominant, but still not negligible, fluctuating or "turbulent" field components B.Given the observed strength of the turbulent velocity fields, which are trans-Alfvénic and crudely isotropic, the typical magnitude of the fluctuating field components   in Fig. 10 are consistent with the usual trans/sub-Alfvénic relation ||/|B| ∼ ||/  , as expected.This implies that the origins of the turbulent/fluctuating magnetic-field components are likely related to the origins of the turbulence in the disk, which we will examine in detail in the next section.
However, one important question is whether we are in the regime of the "traditional" MRI.As shown in Fig. 12, the typical magnetic-field strength in the disk (both |⟨  ⟩| and |B|) is an order-of-magnitude or more larger than the characteristic strength   ∼ √    K at which the linear growth rate of the MRI is usually assumed to vanish following the analytic analysis in e.g.Pessah et al. (2006).And Figs. 6-9 do not show obvious magnetic morphological signatures of the MRI (e.g.channel modes).We further show below that the ratio of Maxwell to Reynolds stress differs significantly from commonly-quoted saturated weak-seed-field MRI simulation results through much of the disk (Brandenburg et al. 1995).
But it is well known that for fields stronger than the commonly-quoted MRI limit (  ≳ √    K ), magnetized disks are still unstable -the instabilities can simply change in character (Pitts & Tayler 1985;Terquem & Papaloizou 1996;Kim & Ostriker 2000;Pessah & Psaltis 2005;Hirabayashi & Hoshino 2016;Das et al. 2018).We therefore follow Pessah & Psaltis (2005) and plot the simulations on the "mode diagram" shown in e.g.Fig. 3 therein, in our Fig.16.Specifically, those authors show that the parameter space of different characteristic unstable modes -at least in a simplified ana-Fig.17.-Velocity streamlines in a cylindrical wedge inside < 0.1 pc, as in Fig. 6, face-on in the midplane (left, showing   -  ) and edge-on (right, showing   -  ), colored by the sign of the (spherical)   ≡ v • r (where the rest frame of the system is defined as the velocity of the SMBH).Compared to the magnetic field, the velocity fields are much more clearly dominated by coherent rotation in the disk (there are no sign flips), and vertical infall onto the disk.Gas pileups at spiral shocks forming the arms in Figs.1-2 are obvious, with increasing order inside < 0.01 pc.The next-order inflow/outflow motion in the midplane is clearly dominated by a coherent, large-scale eccentric disk ( = 1) "slow mode" which propagates inwards from large  where the disk is self-gravitating, but the shocks and turbulence break the exact cancellation of this motion and lead to net angular momentum loss.In the edge-on midplane, we more clearly see the turbulent internal disk dynamics.The eccentric disk slowly precesses around the SMBH lagging the path of the infalling cloud from which it is forming, as expected.Fig. 18.-Velocity streamlines (as Fig. 17), but now with the lines colorcoded by the sign of the radial velocity (inflow/outflow) and super-imposed on a face-on projected density map of the gas in the disk midplane (as Figs.1-2) and on somewhat larger scales.This highlights the flows that form the disk from captured gas initially beyond the BHROI in a massive star-forming complex.As generically expected in tidal capture/disruption events, a largescale stream falls in primarily radially, shocking on self-intersection, and circularizes to form the disk, with an order-unity fraction of the complex mass not bound and instead in a large unbound tidal stream or outflow "fan."The inflow, being tidally compressed, forms the denser filamentary structure while the outflow becomes lower-density gas.Within the inflowing filament in this "free fall" zone before it circularizes, the flow is not highly turbulent: there is both outflow and inflow (so the variance in   can be large, at a given ), but these clearly reflect large-scale, geometrically separate structures rather than turbulent local motions.lytic linear stability analysis of a laminar, unstratified diskcan be, to leading order, represented in a two-dimensional plot of dimensionless vertical wavenumber k ≡     ,  /  (with  ,  ≡ |  |/ √︁ 4 ) versus a measure of the dimensionless toroidal magnetic field strength  ,  /  .The parameter space is then divided into four characteristic regimes based on the intersection of two critical wavenumbers: k1 Modes with k > MAX( k1  , k2  ) are stable; those with k2  < k < k1  are unstable to the traditional MRI (with the near-vertical part of k2 defining the upper limit  ,  < √    K above which, for the conditions considered in Pessah & Psaltis 2005, the MRI growth rate vanishes);10 those with k1  < k < k2  are unstable to a second or "Type II" instability; and those with k < MIN( k1  , k2  ) are unstable to a third or "Type III" instability.11Note that the "Type II" and "Type III" instabilities in Pessah & Psaltis (2005) are also called axisymmetric toroidal buoyancy (ATB) mode[s] in Kim & Ostriker (2000), or superthermal slow mode instability (SSMI) and suprathermal hybrid mode instability (SHMI) in Das et al. (2018) -while there are subtle but important differences in these analyses, the order-of-magnitude dividing criteria between the different instability regimes and 10 If we retain the radial stratification term B ≡  ln   / ln , then k2 is modified to 16, the effect is small here: the k2  boundary between MRI and "Type III" shifts upwards by a factor ∼ 1.1 − 1.4 in each panel, which has no effect on our conclusions.Although B is broadly similar to −1 on average in Fig. 3, evaluating this correction term in factor ∼ 2 radial intervals from ∼ 80 au to ∼ 10 pc, we find that |1 + B | −1/2 never exceeds a factor of ∼ 3.
11 Technically, in Pessah & Psaltis (2005) there is a small stable "strip" just interior to the boundaries of the "Type III" parameter space.However on the (large) dynamic range plotted in Fig. 16, this occupies a negligibly small fraction of the parameter space and is unimportant for our comparison.Moreover, Das et al. (2018) argue this strip vanishes (the full Type-III regime is unstable) in a global mode analysis.key behaviors are, for our purposes, identical.We plot the simulation gas in the disk, at each of several radii, on this diagram, assuming that the characteristic minimum wavenumber of interest (and wavenumber containing most of the power) is   ∼ 2/.We see that the simulations lie solidly in the "Type II/III range," even if we focus only on the warmer (higher-) volume-filling phases of the gas in the disk where it is multi-phase (in the colder, denser gas,  is even smaller and the simulations lie even further from the traditional MRI regime).This confirms our intuition and quantitative statement that  ,  ≫ √    K above.These specific modes are fundamentally related to radial magnetic buoyancy (see references above), operating near the midplane, but as noted in Pessah & Psaltis (2005) they generically involve comparable in-plane and vertical displacements.Moreover, given that the disk is vertically stratified at some level (see § 5.2), additional buoyancy modes in a manner potentially similar to that discussed in idealized simulation studies such as Johansen & Levin (2008).While the analytic models of such modes are somewhat less clear in the regime here (non-constant  ≪ 1, non-isothermal d/d|| > 0, strongly differentially rotating), dimensional considerations and simpler versions of said instabilities suggest that the fastest-possible growth timescales (a few times the vertical Alfvén crossing time, ∼ a few Ω −1 ) at the characteristic scales ∼  should be order-of-magnitude similar to the Type II/III modes above (see e.g.Foglizzo & Tagger 1994;Vishniac 1995;Rodrigues et al. 2016;Salvesen et al. 2016a, and references therein).Together this would naturally explain why we see broadly similar turbulent vertical and radial components |  | ∼ |  | (Fig. 10) with similar-scale structures (Fig. 6).The Type III/SHMI instability in particular is also robust to the different vertical and radial density/pressure/magnetic stratification terms and range of mode propagation angles considered in Pessah & Psaltis (2005); Das et al. (2018).It depends on differential rotation in a similar manner to the traditional MRI and its linear eigenvectors feature a broadly similar structure: notably the linear Type III/SHMI instability, like the linear MRI, always produces both Maxwell and Reynolds stresses which transport angular momentum outwards (dominated by the  component), although the linear mode Reynolds-to-Maxwell ratio can be higher or lower than the traditional MRI over the parameter space spanned by the simulations.And its growth rate peaks at relatively long wavelengths.In fact, if we insert values of the simulation parameters (including the radial gradients in  and   ) from Fig. 16 into the equations from Pessah & Psaltis (2005), we find the simulations can often be in a parameter space which produces an even faster-growing, longer-wavelength version of their Type III instability (Maxwell-stress dominated, with peak growth rate of ∼ (| ,  |/ K ) 1/3 Ω at wavenumber  ∼ 1/ ∼ Ω/ ,  ).Altogether, this suggests these modes could play an important role in driving turbulence and angular momentum transport, but clearly further non-linear simulation studies are needed.12Briefly, the fact that the simulation modes with wavelength ∼  reside order-of-magnitude around the Type II/III (or SHMI/SSMI) dividing line ( k1  ) is not actually surprising: when  10), but  ,  is reduced by a factor of ∼ 3 − 10.This places modes with wavelength ∼  more firmly in the lower-wavelength or "Type III/SHMI" regime (these modes more strongly depend on, and interact with, the differential rotation).But if we considered a somewhat larger wavenumber |  | > 1/, then the vertical position of the simulations would instead shift upwards by a corresponding factor towards the more local-mode "Type II/SSMI" regime.
Briefly, some other instabilities may be less likely to drive the magnetic fluctuations we see.As discussed in a number of previous studies Machida et al. (2006); Begelman & Pringle (2007); Oda et al. (2009); Sądowski (2016); Habibi & Abbassi (2019), magnetically-dominated disks such as these are generically stable against the usual (linear) viscous and thermal instabilities.The disks here are also (linearly) stable against the short-wavelength (k ⊥ B) Parker-like magnetic convective/Rayleigh-Taylor/interchange modes, as these require  ln (||/)/ < 0 (for  > 0 and  ≪ 1; Tayler 1973;Terquem & Papaloizou 1996;Kim et al. 2002), i.e. that the magnetic scale-height is smaller than the density scaleheight, which we plainly show in § 5.2 is not satisfied.And although the long-wavelength Parker instability requires only a vertically-decreasing |B| ( ln  2 / < 0), the characteristic wavelengths  ≳  crit ∼ (5 − 10) ×  (Parker 1966;Kim et al. 1997;Lee & Hong 2007) are extremely large (> , given the relatively large  here), and the trans-Alfvénic magnetic fluctuations are much larger than commonly-quoted thresholds above which the instability may be strongly suppressed (see Kim & Ryu 2001, and references therein), so it is not clear if these modes can actually exist (at a minimum, a global analytic treatment is required).
Of course, all of these analytic "dividing lines" are predicated on analytic linear stability analysis, with a number of simplifying assumptions (e.g. that the disk is azimuthally symmetric, laminar, and adiabatic, and that various vertical and radial stratification terms can be neglected).It is not obvious, therefore, how much can be applied to simulations like ours with complicated stratification, fully-developed strong turbulence, cooling, and highly non-linear modes.And other instabilities or variants of those discussed above may be present as well.For example, as discussed in Begelman & Armitage (2023), the "traditional" MRI can persist in a supra-thermal form if  ln   / ln  is very close to −1 (interestingly similar to the "average" slope we see in Fig. 3), although from the linear analysis in Pessah & Psaltis (2005) this would require |1 +  ln   / ln | <  ( K / ,  ) 2 ∼  (/) −2 ≪ 0.01 − 0.1 for the disk parameters here, so if it is occurring here it may be a transient phenomenon in space and/or time.It is difficult to speculate further -our intent here is to motivate more exploration in both analytic studies and non-linear but idealized numerical simulations of magnetic instabilities in the strong toroidal field parameter space of interest here, and to highlight that the field strengths here are not in fact restricted to the sometimes-quoted value of   ≲ √    K (as assumed in e.g.Begelman & Pringle 2007).

Comparison to Field Strengths in "Traditional" 𝛼 Disks
As shown in Paper III, if we take the scalings for the outer accretion disk from Shakura & Sunyaev (1973) for an SS73-like "weakly-magnetized" ( ≪ 1) -disk, then building up a sufficient Maxwell stress to produce the canonical  ∼  2  / 2  ∼ 0.1 in such a disk would actually require magnetic fields whose absolute strength (in Gauss) is approximately a factor of ∼ 10 larger than those seen in our simulations (e.g. from SS73 Eq. 2.19 therein,  SS73 ∼ 5000 G (/0.001 pc) −1.3 for  ∼ 20 − 30 M ⊙ yr −1 and  BH ∼ 10 7 M ⊙ ).Essentially, "weakly-magnetized" models such as SS73, as well as the vast majority of historical accretion disk simulations, make the implicit assumption that the disk "initially" formed with negligible vertical magnetic support (i.e.vanishingly small or strictly vertical magnetic fields).This implies the disk would collapse to much smaller scale-heights / ∼   /  , with midplane densities  a factor of ∼ 10 6 − 10 8 larger than those seen here ( § 3.2).This in turn would require some process like the MRI to very efficiently amplify |B| up to quite large absolute values to produce even a modest Alfvén speed   ∼ |B|/ √︁ 4 , as needed to produce a Maxwell stress that is any appreciable fraction of   2  .This of course has important consequences for the observational properties of disks like those simulated here: even though   /  and   /  are much larger here than in a traditional SS73-like disk, |B| can actually be significantly smaller.It also re-emphasizes that the absolute field strengths in our simulations are not particularly extreme or implausible ( § 4.1).And it further suggests that if one "initially" forms the disk from gas with more realistic ISM magnetic field strengths (with non-trivial toroidal and/or radial fields), it will likely reach a magnetically-dominated state akin to the simulations here, well before it could actually collapse to extreme densities like those assumed for "weakly-magnetized" SS73-like disks.

STRUCTURE OF THE VELOCITY FIELDS
5.1.Overview In Figs.17-19 we plot the structure of the velocity fields (face-on and edge-on).We clearly see radial infall at large , circularizing and forming the disk (Figs.17-18).Face-on, the disk appears highly ordered, but with an obvious coherent eccentric disk mode present in Fig. 17.Edge-on, we see vertical inflow onto the disk, with a thick turbulent midplane layer evident in Fig. 19.Note that edge-on, since we show below this is essentially the same as a plot of the velocity fluctuations.This shows, as we saw above in the morphology and B, that the midplane is not a uniform, perfectly rigid layer but features a complex internal density structure with many warps and even streams with somewhat different orientations at large radii.Fig. 20 plots the radial profile of the different velocity components (in absolute units and relative to  c ()), showing both the mean and fluctuating velocities.
As discussed in Paper I, at radii ≫ pc, outside the BHROI, the velocity structure in the nucleus is primarily turbulent with quasi-isotropic velocity fields -we see this already at ≳ 0.3 pc in Figs.17-20, where the incoherent or "turbulent" velocity components dominate over the mean velocities, and are comparable both to one-another and to the circular velocity 1 pc, we see the disk clearly in kinematic space, with ⟨  ⟩ ≈  c much larger than other components, which have a dispersion ∼ 0.1   .The kinematics of the disk are much more coherent compared to B above: there is a global smooth rotation-dominated flow with no sign flips (i.e.all the material is rotating the same direction, as expected).

Vertical Structure
Vertically we see infall from above/below onto the disk, at most radii; but we stress this does not dominate the total inflow rate, most of which occurs through the thicker and denser midplane.This may change if we modeled emergent "feedback" (e.g.jets or high-velocity outflows) from the unresolved portions of the accretion disk at < 80 au.
We discuss stratification in greater detail below, but we do not see any evidence for coherent stratification of the velocity structure within the disk (at − ≲  ≲ ).

The Eccentric Disk
Following the system in time, we see the initial inflows as the gas captured first falls into the center and circularizes to form an eccentric, lopsided structure with a clear shock/compression region as orbits self-intersect (Fig. 18).This "settles" over some tens of dynamical times into a smoother structure with more coherent velocity structure visible in Fig. 17.The eccentricity does not disappear even over many thousands of dynamical times.Rather, the system settles into a coherent  = 1 eccentric disk mode, where we clearly see the eccentricities of individual gas orbits align, and the entire eccentric structure precesses with a slow pattern speed (pattern speed Ω  ∼ Ω( max ) where  max ∼ pc, so at smaller radii Ω  ≪ Ω( ≪  max )).The amplitude of the eccentricity declines weakly as  → 0. Fig. 21 plots the amplitude of the eccentric ( = 1) and higher- modes in the gas surface density in the face-on surface density projection of the disk, as a function of radius, which demonstrates this explicitly.
These are exactly the predicted structures for "slow"  = 1 modes in nearly-Keplerian potentials, which are well-studied and (at least in linear theory) can persist indefinitely (Tremaine 2001;Bacon et al. 2001;Hopkins 2010).It is important to note that these are unique among global/large-scale modes in the orbit structure: other modes, e.g. = 2 bars or higher  > 2, are damped as  → 0, but for slow  = 1 modes there is essentially zero energetic cost of the mode propagating inwards as it just involves coherent alignment of already-eccentric orbits.The pattern speed is set by the "driving" of the eccentric mode: namely, the motion of the parent gas complex which is being tidally stripped by the SMBH to fuel the accretion event and form the disk in the first place.That complex both torques the disk directly (as its mass is larger than that of the SMBH and most of the complex lies outside the BHROI), and provides the newly-infalling gas which follows the trajectory of its parent cloud (itself on an effectively hyperbolic orbit) as it passes through some impact parameter or pericenter  complex .The cloud complex therefore drives a characteristic (lagging) pattern frequency Ω eff ∼ ||/ complex ∼  galaxy / complex ∼ 200 km s −1 /a few pc ∼ 10 −12 s −1 , i.e. precession on a ∼ 10 4 yr timescale.
This means that in an instantaneous sense, for a given gas parcel in the outer disk, its radial (inflow/outflow) velocity is dominated by where it is instantaneously in its eccentric orbit.This is clear in the face-on projections of the velocity streamlines in Figs.17-18.But to leading order, these orbits are of course closed and their radial flow cancels over the course of the full orbit, so we stress that this should not be conflated with the systematic or net inflows feeding accretion or outflows ejecting material from the nucleus.Instead, as expected, the combination of some shocks/dissipation (which break the exact symmetries of the eccentric orbits), plus non-zero local turbulent/Reynolds/Maxwell stresses, means that there is a non-zero torque on the gas which causes a systematic net inflow/accretion of gas.But it is worth noting that the leadingorder description of the outermost disk is not a perturbed circular disk, but a perturbed eccentric disk.However, the fact that there is also clearly less-coherent/global more "turbulence"like cells or eddies in the inner disk ( ≲ 0.01 pc) with coher- We can attempt to separate the coherent eccentric motion in our definition of v (see § 2.2), by subtracting the best fit  = 1 component from   (, ).We have done so (fitting this independently in each radial annulus, which we caution may over-estimate the coherent component), and find that it has a modest effect most notably on   in the outer disk (0.01 pc ≲  ≲ 0.1 pc), reducing it by a factor of ∼ 2 (and a much smaller effect at smaller or larger radii, as expected).Interestingly after doing so, the residual velocity fluctuations in the outer disk become closer to isotropic, suggesting the more traditionally "turbulent" component is indeed close-to-isotropic.We show Fig. 20.-Mean velocity with its dispersion (top), and separation of the mean ⟨v  ⟩ and fluctuating v  velocity components in both units of   (middle) and km s −1 (bottom), in the same style and showing the same two times (left and right) as Fig. 10.We see quantitatively that the velocity fields are more ordered than the magnetic fields.Inside the disk (≲ 0.1 pc) the velocity is clearly dominated by rotational motion, with the turbulent/fluctuating components broadly all similar to each other in magnitude at v  ∼ 0.1   .Almost always we see ⟨  ⟩ < 0, indicating inflow, as expected, with modest amplitude.For   we show ⟨  • /|  | ⟩, so that a value of > 0 or < 0 indicates vertical outflow or inflow respectively (regardless of whether the gas parcel happens to be "above" or "below" the midplane), showing weak vertical inflows as well.below ( § 10) that this is quite different from the case where we re-run without magnetic fields, where the disk is much thinner and exhibits much more extreme anisotropic structure and higher eccentricities, with much weaker Reynolds stresses and lower inflow rates.Of course, in either case, the eccentric motions have no measureable effect on the vertical   .

Turbulent/Velocity Fluctuation Structure
From Fig. 20 we see the turbulence or more general velocity fluctuations, while sub-dominant to rotation, are still vigorous, with ⟨ 2 ⟩ 1/2 ∼ 0.1  c ∼   .Comparison of Figs. 10 & 20 immediately shows turbulent ram pressure is comparable to magnetic pressure and, from Figs. 4-5, we see both are much larger than thermal or radiation pressure within the disk (i.e. the turbulence is broadly trans-Alfvénic, but highly supersonic).Fig. 22 shows this more explicitly, plotting the typical sonic and Alfvénic Mach numbers of the velocity fluctuations in different radial annuli.
Figs. 20 & 22 also show that the typical velocity fluctuation  is generally larger than the mean velocities in  or  directions, and not strongly anisotropic (|  | ∼ |  | ∼ |  |, to within a factor of ∼ 2 − 3 or so, though note that the apparent transient dominance of e.g.  at small  in Fig. 20 owes as much to the presence of coherent warps/bends in the disk as to actual "local" small-scale vertical turbulence).This promotes strong mixing and turbulent structure within − <  < , contributing to the complex edge-on structure (as compared to well ordered face-on structure), as well as the Reynolds stresses (analyzed below).
Of course, the velocity fluctuations here are much stronger than in a typical SS73-like -disk, where (by assumption) the sonic Mach number M  ∼ ||/  ∼  1/2 < 1 (i.e. the non-circular motions are always sub-sonic).5.5.What Drives the Turbulence (Or sets its Amplitude)?5.5.1.Does Anything Actually Need to Drive The Turbulence Interior to the Disk?Given these vigorous super-sonic and quasi-isotropic velocity dispersions, it is natural to ask what their "driving mechanism" may be.However, it is not necessarily obvious that a driving mechanism for the velocity fluctuations -in the usual specific sense of some local instability constantly amplifying turbulent modes on some driving scale ≲  -is strictly necessary.If the turbulence is Alfvénic, simple arguments (Völk & Aplers 1973) show that the effective "adiabatic index" of a group of passively-advected13 Alfvén wave packets is 3/2 (i.e.we expect |B| ∝  3/4 moving with a Lagrangian parcel, or equivalently d  ⟨B 2 ⟩ ∼ −(3/2) ⟨B 2 ⟩ (∇ • ⟨v⟩)).
The challenge here is the turbulent dissipation time: if the turbulence dissipates on a timescale ∼  turb, diss , then the expressions above should be modified to d  ⟨ v 2 ⟩ ∼ − ⟨ v 2 ⟩ (∇ • ⟨v⟩) − ⟨ v 2 ⟩/ turb, diss (where  can be in the range from 3/2 to 5/3 depending on the regime, as described above).Noting that turb, diss , so it is clear that the scalings outlined in the previous paragraph require that the accretion/advection timescale is comparable to, or faster than, the turbulent damping timescale at each radius.For hydrodynamic, supersonic turbulence with  cool ≪  dyn , the expected dissipation time is of order a couple of turnover times at the driving scale  turb, diss ∼  turnover ∼  drive / drive ∼ /|v turb | ∼ Ω −1 ∼  dyn (Thompson et al. 2005;Pan & Scannapieco 2010;Hopkins 2013); the expectation is similar for magnetized trans-Alfvénic but subsonic/incompressible turbulence (Goldreich & Sridhar 1995).But the inflow time ∼  gas /  is ∼ ( K /|⟨  ⟩|) Ω −1 ∼ (  / K ) −2  dyn .While the accretion here is very fast compared to more traditional -disk models (as we discuss below), inserting typical values from e.g.Fig. 20 into the above still gives an inflow time ∼ 30 − 100  dyn or ∼ 5 − 15  orbit .So we would naively expect the gas to sit too long at a given , meaning the fluctuations will damp faster than they are amplified via compression.
However, we stress that we use the term "turbulence" rather loosely here, encompassing any non-zero rms fluctuating motions/fields v ≡ v − ⟨v⟩.So some of the power in these fluctuations could be in structures which are not strictly "turbulent" in the classical sense (i.e.not part of some cascade) and therefore potentially dissipate on slower timescales.For example, in addition to the eccentric motions discussed above, recent idealized studies (Skalidis et al. 2021(Skalidis et al. , 2022;;Beattie et al. 2020Beattie et al. , 2022) ) of randomly-forced super-sonic but sub/trans-Alfvénic motions (akin to what we see here; Fig. 22) have argued that most of the kinetic energy at the driving scale ends up neither going into shocks (as it would with weaker fields) nor a classical MHD cascade (as in the sub-sonic case) but into structures such as compressive transverse non-linear magnetosonic modes (e.g.ram pressure of flows transverse to B with |B • B| ∼   2 ), which do not dissipate rapidly.So the damping time for these types of motions could be much longer and potentially remain in balance with amplification via inflow, obviating the need for an explicit "driving" mechanism.A more detailed analysis of the dissipation and power spectrum/structure functions of the velocity field is therefore an interesting subject for future work.

Viable Driving Mechanisms
That said, there are multiple viable driving mechanisms operating here in principle.As expected given the lack of e.g.strong star formation and therefore stellar feedback from within the disk, the ultimate source of energy for the turbulence on these scales is the gravitational energy of the gas.We see this reflected in the fact that the inflows obey the usual Reynolds/Maxwell energy balance condition Σ gas  2  Ω ∼  Ω 2 (at least at the order-of-magnitude level; this is discussed in detail below).But this does not immediately provide a local mechanism.
One obvious possibility is the strong global eccentric mode discussed above (quantified in Fig. 21).Per Noguchi (1988); Barnes & Hernquist (1996); Hopkins & Quataert (2011), we would expect this to excite non-circular motions of the order   / c ∼ | =1 | ∼ 0.1, which is not far from what we see.Given that the mode is a slow mode (pattern speed small compared to  c ), gas will intersect the compressive orbit crossing/pileup locations in Fig. 17 once per orbit, so the turbulence can be rapidly "rejuvenated" on a timescale by definition comparable to its decay/damping time, even in the fast-decay ( decay ∼ 1/Ω) limit.
Another possibility is local gravito-turbulence.This is distinct from the global-eccentric-mode driving discussed above in that it depends on the local self-gravity of the gas collapsing on itself and driving high- ( ≳ 1/) modes in Fig. 21, while the global mode is driven by an external perturbation (in this case, the self-sustaining asymmetric distribution at large-).This would naturally give rise to turbulence of the form Σ gas  2  Ω ∼  Ω 2 (Gammie 2001).Both of these are likely playing some role, especially at larger radii (≫ 0.01 pc), where Reynolds stresses dominate over Maxwell stresses and the turbulence is modestly super-Alfvénic.This is further reinforced by our comparisons in § 10 below to simulations without magnetic fields.However, at smaller radii in the disk (≲ 0.03 pc), it seems unlikely that these dominate the turbulence at saturation, for several reasons.The growth/driving rates of these modes are independent of B and neither therefore provides an obvious reason the turbulence should saturate at trans-Alfvénic values in Fig. 22. Moreover we show below that the Reynolds stresses at these radii are much weaker (and often opposite in sign) if we re-run without magnetic fields, whereas the gravitoturbulence and disk asymmetry are stronger.For the eccentric mode, examining the spatial pattern of the turbulent torques, we see a hint at the largest radii that these may trace the eccentric shock but comparing Fig. 17  It is possible the vertical turbulence is instead driven by the weak vertical inflows but again this would not explain why it is roughly isotropic with the radial/azimuthal turbulence, and moreover it is challenging to explain how the turbulent velocity would be much larger than the mean inflow velocity as seen at most radii in Fig. 20.For gravitoturbulence, as discussed in Paper I, at small radii  thermal is too large, and the magnetic field strength much too large, to actually sustain vigorous gravitoturbulence (see Kim & Ostriker 2001;Lizano et al. 2010;Riols & Latter 2016;Forgan et al. 2017).Further, such vigorous gravitoturbulence would also be inconsistent with the low star formation rates we actually see at these radii.
As discussed in detail in § 4.2.3, the "traditional" weak-field MRI -at least in the simplest linear form usually invokeddoes not appear viable as the primary turbulent driving mechanism.Specifically, as shown above (there and § 4.2.2), the toroidal   is sufficiently large that the system should be stable against the standard MRI in the usual idealized context where it has been studied (i.e. for the mean-field parameters here, the standard weak seed-field MRI modes should generally be quenched and linearly stable, at least for a laminar, azimuthally-symmetric, unstratified system; see Terquem & Papaloizou 1996;Kim & Ostriker 2000;Pessah & Psaltis 2005;Lin 2014;Hirabayashi & Hoshino 2016;Das et al. 2018).Also when we analyze the structure of the torques in more detail below, neither the mode morphology nor ratio of Reynolds to Maxwell stresses agrees quantitatively with the predictions of some previous saturated MRI simulations in similarly simplified/idealized shearing boxes (compare e.g.Brandenburg et al. 1995;Pessah et al. 2006).But as reviewed in § 4.2.3, while this particular form of the MRI is stabilized under such conditions at sufficiently high   , this just means that different linear instabilities appear.Specifically, even under those conditions, at   above the threshold where the MRI growth rate vanishes, in a  ≪ 1 disk, the Type II/III (Pessah & Psaltis 2005) or ATB (Kim & Ostriker 2000) or SSMI/SHMI (Das et al. 2018) modes appear (see Fig. 16).As shown in those studies, these modes are essentially radial buoyancy modes driven by the competition between dominant magnetic pressure support of the gas and the Keplerian gravity/differential rotation around the BH.For the conditions here, these have rapid growth rates (∼ (  / K ) Ω for the largest-wavelength unstable modes with  ∼ 1/), are at least plausibly consistent with the turbulent morphologies and Reynolds stresses below, and -given their dependence on the mean-field direction/structure -would plausibly saturate at trans-Alfvénic turbulence.Likewise vertical buoyancy modes, if unstable, should have broadly similar growth rates for modes with  ∼ 1/, potentially explaining (at least at the order-of-magnitude level) the broadly isotropic turbulent velocity and magnetic field structure.And modified (supra-thermal) versions of the MRI could exist where  ln   / ln  is very close to −1 (Begelman & Armitage 2023).Of course, it is also plausible that other, distinct instabilities appear, or even that some alternative modified versions of the MRI persist, under the more complicated conditions present in our simulations (e.g.allowing for eccentricity, stratification, different equations-of-state, turbulence, non-linear perturbations, etc.).It therefore seems likely that some unstable modes are at least present, and could play an important role in the accretion dynamics.
Unfortunately, identifying a "unique" driver in fully nonlinear, saturated, multi-physics and multi-scale simulations like these (as opposed to more idealized controlled numerical experiments) is challenging, and will require additional work with different simulations motivated by the discussion above to isolate the distinct mechanisms discussed above in more con-trolled circumstances.Importantly, however, this has allowed us to identify a relatively small list of likely processes which are occurring, in order to guide said simulations.We also wish to stress that there is no smooth or laminar "initial" phase from which the turbulence and/or magnetic field strengths in the disk amplify -in idealized simulations, such a phase is an artifact of artificially smooth initial conditions.Instead, the system here begins from a highly out-of-equilibrium gas configuration as the disk forms from infall, with strong turbulence and magnetic flux already well above the commonly-quoted MRI saturation threshold and inhomogeneity in the velocity, density, temperature, and magnetic fields, before relaxing into a steady-state strongly-turbulent disk.
In potential future, more idealized/controlled numerical experiments with initial and boundary conditions motivated by these simulations, our discussion above also demonstrates that global simulations may be necessary to properly capture the magnetic curvature effects and low- "Type III" as well as eccentric modes, as well as the inflow of gas and magnetic flux.But our discussion above also shows that such studies could be made more simplified in various ways -e.g. by replacing the explicit radiation-hydrodynamics and thermochemistry with more idealized models, or neglecting star formation and stellar dynamics -and still be used to explore e.g. the effects of self-gravity and different instabilities.It would be particularly interesting to study these in the context of angular momentum transfer in the disk.While this has been studied more extensively (under more idealized conditions) for the weak-field MRI, gravito-turbulence, and global gravity (eccentric/spiral) modes (all of which are known to produce efficient angular momentum exchange), there has been much less study of the non-linear buoyancy modes we discuss above (let alone other instabilities that may be present here).
At least in linear theory, the longer-wavelength "Type III/SHMI" instabilities (active at wavenumbers |  | ≲  1  ) appear to depend on differential rotation in a similar manner to the MRI (Pessah & Psaltis 2005), and the linear mode eigenvectors of the Type III/SHMI instabilities are similar to the (compressible) MRI both in the relative magnitude of the density/velocity/magnetic field perturbation components and in the fact that both the resulting Maxwell and Reynolds stresses always transport angular momentum outwards (although the Maxwell-to-Reynolds ratio can vary with respect to the MRI).This suggests that their non-linear behavior and ability to produce efficient angular momentum transport may also be similar to the MRI.On the other hand, the short-wavelength "Type II/SSMI" modes are more local and persist without differential rotation (and their linear eigenvectors involve weak Maxwell/Reynolds stresses with angular momentum transport of opposing signs), so may be more akin to e.g.convective instabilities in their non-linear behavior.But further study will be needed for more rigorous conclusions.

What Is the "Turbulent Resistivity" and Why Does It
Not Destroy Flux-Freezing?Since we have shown that the magnetic fields in the accretion disk are largely "sourced" by simple flux-freezing/advection ( § 4) and that microphysical resistivity is negligible ( § 4.2.1), it is worth briefly discussing the role of the effective turbulent resistivity  turb , and whether this could destroy flux-freezing in the disk (see e.g.Lubow et al. 1994, for a detailed discussion of this for pure-poloidal fields, although this is the opposite of the regime of interest here).
Consider first the radial field component   , as we ar- We take the velocity dispersion   estimated as 1/2 the 16 − 84% range (to avoid being pulled by outliers) of   , weighted by gas mass, in radial annuli restricting to gas within |  |/ < 0.3 of the midplane, for each component of  in cylindrical coordinates.We then define the sonic Mach number M ,  ≡   /⟨  ⟩ and Alfvén M ,  ≡   /⟨  ⟩.The turbulence is broadly trans-Alfvénic and highly super-sonic (reflecting  ≪ 1 in the disk).In the inner disk, the turbulence becomes mildly sub-Alfvénic (  , ,  ∼ 0.3   , so  total ∼ 0.5   ) -this closely maps to the radii where the mean toroidal field ⟨  ⟩ becomes larger than the fluctuating field in Fig. 10, as expected.The radius of this transition can fluctuate by a factor of a few as new magnetic flux advects into and through the disk .gue this sources   , following a Lagrangian parcel of gas.The turbulent resistivity should damp the field as d turb    ∼ − turb  2     , where    gives some inverse gradient scale length of   and  turb ∼  turb / turb .Making the standard assumption that genuinely turbulent eddies with crossing time larger than ∼ Ω −1 will be sheared out (which appears to be valid in the simulations) we have14  turb  turb ∼ Ω so  turb ∼  2 turb /Ω.Meanwhile, flux-freezing will grow   on the inflow/compression timescale ∼ /⟨  ⟩, so15 d adv    ∼ −(⟨  ⟩/)   .For turbulence to completely 14 We will ignore O (1) prefactors here but note that  turb could be much smaller by a factor ∼ 10 − 100 than this estimate (see e.g.Salvesen et al. 2016a, who find a prefactor ∼ 0.3/Pr ∼ 0.01 − 0.1 in terms of the turbulent Prandtl number Pr).So we are effectively considering the maximum possible effect of turbulent resistivity.
15 Again we ignore order-unity coefficients but these can be derived from the induction equation for specific geometries, e.g.we show in § 4.2 that for a tidally captured radially infalling stream, the compression in the ẑ and φ directions ( −1   ⟨    ⟩ and   ⟨    ⟩) give a coefficient = 2. So, we are again being conservative and assuming a somewhat weaker growth via destroy flux-freezing, clearly we require With this in mind, is it helpful to divide the simulation into three regimes or "zones" based on distance to the SMBH, which exhibit three qualitatively different behaviors: 1.The outer/ISM/turbulent "zone" well outside the BHROI and accretion disk ( ≫ pc).Here the physics are "ISMlike" (per above): there is no disk or any coherent BH "accretion flow" to speak of (the BH by definition does not dominate the potential).Instead we see self-gravitating, collapsing, rapidly star-forming GMC-like gas complexes "stirred" by stellar feedback, with super-Alfvénic internal turbulence (these can globally lose angular momentum and get closer to the galactic nucleus, but are not systematically "accreting").In this regime (see e.g.Fig. 3 and Paper I) the magnetic fields are tangled and quasi-isotropic, mean fields are much smaller than fluctuating fields, and "global" compression (with coherent infall towards  → 0) is small compared to compression of gas internal to the cloud via shocks and self-gravity.So inserting values we naturally find |d turb    | ≫ |d adv    |: as discussed in Paper I, the clouds internally reach the expected saturation for the small-scale dynamo with turbulent amplification balancing turbulent resistivity (magnetic energy a few percent of turbulent energy).16 2. The intermediate/capture/free-fall zone inside the BHROI ( ∼ 0.1−a few pc).Here we see a complex is tidally disrupted on close passage to the SMBH, leading to a tidally elongated, radially free-falling stream of bound gas towards the SMBH (see § 3).We stress that while the rms ⟨  ⟩ 2 might appear large in e.g.Fig. 20 at these radii, a cursory examination of the visual morphology of the gas stream (Fig. 1) or velocity fields (Figs.17-18) or magnetic field lines (Fig. 7), or a more rigorous quantitative diagnostic of the local, genuinely turbulent components of v, immediately makes it obvious that the flow is locally coherent and strongly dominated by radial motion stretching the radial field lines, as opposed to turbulence ( § 4.1).The large rms ⟨  ⟩ 2 does not reflect local turbulence, but the fact that just like in any tidal disruption event, of order half the mass is unbound, so there is a free-falling stream inwards and a large unbound stream/fan being ejected.Within the accreting filament (which dominates the gas supply so is what actually matters here), the gas is effectively free-falling radially onto the SMBH, so   / ∼  freefall / ≳ Ω reaches its maximum possible value, while  turb is relatively weak (≪  freefall ), hence the net effect is stretching    → 1/.So in this regime, |d turb    | ≪ |d adv    | is easily ensured even if the gas supply at larger radii is turbulent.
3. The accretion disk zone at radii ≲ 0.1 pc, where the gas circularizes and an ordered accretion disk forms.This is the regime of interest and unlike the previous zones, the behavior is not trivial.Here, the relations noted above should hold for |d turb    | and |d adv    | so flux-freezing is destroyed flux-freezing by neglecting such order-unity factors.
16 It is worth noting that the bulk velocity of individual clouds/complexes is larger than their internal turbulence, as expected for any clouds that are smaller than the galaxy in which they are embedded.This is important insofar as it explains why on a close passage to the SMBH, complexes will be tidally disrupted and an O (1) fraction of the mass will be strongly bound and begin to radially free-fall with relatively small impact parameter.
if  2 turb ≫ |⟨  ⟩|  −2   Ω/.In our simulation  we clearly have    ∼ 1/ in this regime, which should be expected analytically from the boundary conditions forming the disk (from zone (ii) above), and -if flux-freezing is valid -will be maintained by the disk17 as  → 0. This means that the resistivity condition becomes  2 turb ≫ |⟨  ⟩|  K .But for a system dominated by Maxwell and Reynolds stresses we expect and measure in the simulation |⟨  ⟩| ∼ (Π Reynolds + Π Maxwell )/(  Ω) ∼ ( 2 turb +  2  )/ K (see § 7 and Figures therein below, as well as Paper III), so this is impossible to satisfy for  2  > 0, and hence flux-freezing/advection should not be strongly modified by turbulent resistivity (though the resistivity is likely to be a non-negligible correction).
Note that for the accretion disk zone (iii), if we had made the more common accretion disk assumptionopposite the behavior we actually see in the simulations -that the disk "begins" from negligible mean   (e.g. a strictly poloidal field with all of   sourced from the MRI), then    would necessarily be ∼  turb ∼ Ω/ turb .In this case, the condition Ω/ would simply reduce to |⟨  ⟩| ≪  K , which is always true in the disk.Thus there are two qualitatively distinct but each internally selfconsistent regimes: (i) the "traditional" weakly-magnetized assumption where   is sourced primarily by a poloidal field and turbulent resistivity invalidates flux-freezing, and (ii) the flux-frozen and therefore hyper-magnetized regime where there is some coherent   being "fed" to the disk, ensuring that flux-freezing remains a valid assumption throughout.18The regime of interest clearly depends on the disk outer boundary condition: thus even though the dynamic range of the "free-fall" zone in radius  is quite small compared to the "ISM zone" or "accretion disk zone," it plays a crucial role in establishing these conditions.
We can repeat these arguments for the toroidal   and reach similar conclusions in the disk zone where   is already dominant, but it is more interesting to consider the "sourcing" of   from   which we demonstrated in § 4.2.1.Motivated by Fig. 15, if we compare turbulent resistivity d    ∼ − turb  2     to the source term d    ∼   ⟨    ⟩ ∼ −Ω   , then (given that it is even more obvious that    ∼ 1/ is coherent in the disk zone) turbulent resistivity is a small effect so long as |⟨  ⟩| ≲ ( K / turb ) 2 |⟨  ⟩| ∼ 100 |⟨  ⟩| (inserting typical  turb from Fig. 20).This is easily satisfied in the disk (Fig. 10); indeed, |⟨  ⟩| never really reaches this upper limit because other terms dominate (e.g. the 17 If the disk begins from an initial condition with    ∼ 1/ with fluxfreezing valid, then we can approximate the evolution of    through the disk analytically via the stretching between two gas elements in the radial direction.The radial distance between two elements Δ evolves in a Lagrangian manner as they migrate radially (in the midplane of a steady azimuthally-symmetric accretion flow) as d  Δ =   ( + Δ) −   () ≈ Δ    , and noting in this limit we can replace d  with     , this becomes   (ln Δ) ≈   (ln [ −  ] ).So if we begin from    ∼ 1/,    ≲ 1/ will be preserved so long as |  | ∝   with  ≤ 1.This is easily satisfied in the simulation in a time-averaged sense (ignoring obvious transient features; see Fig. 20), and in the analytic model presented in Paper III this is automatically ensured not just for our default trans-sonic flow assumption ( ∼ 1/3 in the variables defined therein) but for any  >= −1/4 (where any physical solution requires  > 0, so this condition is always met).gas and field line advection/accretion time is short enough that this is balancing the induction term, rather than resistivity).
Distinct from the usual assumption for magnetized accretion disks, the relevant interactions here are between   and   : the poloidal field plays a sub-dominant dynamical role in either (though it does provide some additional stability and source for   and   , see Salvesen et al. 2016b).This is consistent with the relatively weak, more turbulence-dominated poloidal fields we see in the midplane (Fig. 10).However we caution that this does not mean the behaviors of interest here could be captured in a strictly 2D simulation: phenomena such as turbulence and vertical compression of the disk ( → 0 as  → 0, contributing to the flux-freezing amplification of both   and   ) still depend on 3D structure.

VERTICAL STRUCTURE & THERMO-CHEMICAL
PROPERTIES OF THE DISK Briefly, we examine the vertical and thermo-chemical structure of the disk.The details of the thermo-chemistry as a function of scale are discussed in much greater detail in Paper I to which we refer interested readers, because it is of great potential importance to the cessation of star formation in the disk.But, because  ≪ 1, the thermal properties of the disk actually play relatively little role in the accretion dynamics.We thus only briefly discuss them here insofar as they are useful to inform our understanding of the disk structural properties and to help distinguish the behaviors we observe from some other possible models (as with e.g. the amplitude of |B| as discussed above).

General Structure & Vertical Density Profiles
Recall, Fig. 4 shows various thermo-chemical properties such as the temperature and plasma  in face-on projection.Edge-on, Fig. 23 shows more detailed quantitative vertical profiles at a variety of BH-centric radii  (we use cylindrical , ,  coordinates, and azimuthally average over  in cylindrical rings).We see a broad radial trend in Fig. 4 where the gas becomes less dramatically multi-phase at smaller radii, and starts to become somewhat warmer in denser regions towards small , as expected in an optically thick disk (though much of the disk at the radii resolved here is still atomic with modest ionization fraction ∼ 1 − 10%; see Paper I).In Fig. 23, the obvious take-away is that the thermal and magnetic properties of the disk are quite weakly stratified.
In detail, in Fig. 23, we see that the vertical density profile follows a very typical "disk+halo/corona" profile.Specifically we can fit (, ) fairly well with the sum of a "disk" and "halo" component.The disk is Gaussian or sech 2 ,  disk ≈  mid sech 2 (/ H ), as is standard in the resolved/galactic disk literature (Van der Kruit 1988; the constant  H ∼ 1 depends on whether we define  as an rms, median, half-mass height, etc.).The halo follows a power-law  halo ∝  −  halo with lower normalization,  halo between 2 − 4, and  2 =  2 +  2 .This is quite generic to disks with extended gaseous halos, coronae, or outflows, and very similar to the profiles typically used to fit e.g.edge-on gas densities in galaxies or protostellar disks.The disk is clearly evident here, as it should be given its morphological presence in e.g.Fig. 2.There is significant skewness/asymmetry especially at large radii owing to the asymmetric inflow, but following Lagrangian parcels in time we typically see the vertical profile of "new" material settle into this profile over just a couple of dynamical times once it reaches its circularization radius (inside e.g.≲ 0.1 pc).

Weak or Inverse Stratification of |B|, 𝑇, and 𝛽
Comparing the magnetic field strength versus height  in Fig. 23, we see notably weaker vertical stratification compared to the gas density.At various radii we see that a ∼ 3 dex decrease in  from midplane to large heights || ≳ 2  corresponds to a < 1 dex change in |B|.This means that, as noted in § 4.2.2, although we can approximate the midplane scaling of ⟨|B|⟩ midplane ∝ ⟨⟩ 2/3 midplane , we cannot assume |B| ∝  2/3 everywhere outside the disk.Interestingly, we can actually model the vertical |B| profiles in Fig. 23 fairly well if we assume |B|() scales with the "halo" or "coronal" component of  as defined in § 6.1, as |B| ∝  2/3 halo -i.e. if the shape of |B|() only follows the diffuse/slowly varying "corona" component, with a normalization scaling as ⟨|B|⟩ ∝ ⟨ mid ⟩ 2/3 , such that the "disk" component and stratification effectively disappears when we consider |B|.As noted below, this is naturally connected to the strong turbulence in the disk, itself potentially related to strong magnetic buoyancy, convective, and other instabilities producing effective vertical motions.
Turning to the gas temperature , we again do not see any strong stratification within the disk in Fig. 23 -i.e. at ||/ ≲ 0.3 or so where we see clear structure in ().And the (weaker) stratification between the entire disk "zone" and extended halo at || >  is inverted -i.e.we see higher temperatures at large ||/.The absolute values of the temperatures we obtain depends on how we weight the average in each annulus (a generic expectation if there is any multi-phase structure, since most of the thermal energy often resides in phases that do not dominate the total gas mass), and the trend is more clear if we weight by thermal energy as compared to e.g.gas mass.19 The inverse stratification we see in Fig. 23 is generic to almost all diffuse coronae/halo type-systems, and is seen in e.g.AGN accretion disk coronae, proto-planetary disks, stellar atmospheres and their coronae, and galactic disk-halo (or CGM) interface regions.Indeed it has also been seen in other idealized simulations of magnetized disks not unlike those here (Kudoh et al. 2020).Like in all these other systems, in the more diffuse gas at || ≳ , the absorption optical depth is low and the gas is far from LTE.The diffuse gas cools relatively inefficiently and is heated by a wide range of non-equilibrium and/or non-local processes, including shocks and turbulent dissipation, magnetic reconnection, radiation from the disk, external radiation from starlight in the galaxy, the galactic cosmic ray background, and stellar feedback (e.g.jets, winds, supernovae) from the stars further out in the disk.
Combining these trends, we see that  is both ≪ 1 everywhere, and inversely or weakly stratified (with  lower near the midplane).Unsurprisingly the absolute value of the mean  and precise degree of inverse stratification depend systematically on the averaging method.For example, weighting by thermal energy will give the highest ⟨⟩, magnetic energy is biased towards the lowest ⟨⟩, gas mass gives results in between, or one could also use ⟨⟩ ≡ ⟨ 2  ⟩/⟨ 2  ⟩.But the qualitative results are the same for all weightings.For each, we mass or energy-average in cylindrical annuli at radius  (line colors, labeled) in factor ∼ 3 intervals from ∼ 100 au to ∼ 1 pc.The density is stratified as expected for / ∼ 0.1 − 0.3 at this range of radii; magnetic field strengths peak in the midplane but are more weakly stratified.The temperatures are coolest in the midplane and rise in the diffuse, more optically-thin corona, leading to the plasma  also rising (weakly) with scale height.Energy-weighted thermal averages appear more "noisy" owing to physical multi-phase structure (the average can be dominated by e.g.hot gas in shocks or SNe at large radii).In the density plot, we illustrate for one radius (∼ 5000  schw ) a decomposition into a two-component "disk" (Gaussian/exponential or  ∝ sech 2 (/ H  ) here; thin dotted) plus power-law "halo/coronal" ( ∝ ( 2 +  2 ) −  halo with  halo ∼ 3.5 here; thick dotted) profile.In this decomposition, the disk component with / ∼ 0.1 contains most of the mass, but there is clearly an extended vertical gas distribution beyond the extrapolation of the vertical profile from smaller .The vertical thermal structure is expected for a more tenuous corona/atmosphere above the disk, but the weak stratification within the disk of  and |B| owes to a combination of optical depth effects and strong super-sonic turbulence driving efficient vertical mixing ( § 6.3).

Physics of the Weak Disk Stratification
The lack of strong vertical stratification in |B| and  within the disk is a prediction distinct from e.g. a classic weaklymagnetized, weakly-turbulent Shakura & Sunyaev (1973) disk, but should be expected for the disks here, given the vigorous turbulence present.As noted above (Fig. 22), the disk is exhibits super-sonic and trans-Alfvénic, quasi-isotropic turbulence.This means, essentially by definition, that an order-unity fraction of the turbulent power will be in vertical flows with coherence length ∼ , crossing/turnover time ∼ / t ∼ Ω −1 ∼  dyn , and ram pressure comparable to or greater than the total (thermal+radiation+magnetic) pressure in the disk.
As discussed in more detail in Paper I and above, we also should not expect LTE to apply in the diffuse gas above the disk nor in the outermost regions of the disk.In the inner disk midplane, the behavior does begin to resemble more blackbody-like LTE behavior.As we showed in Paper I, the radiation temperature and gas temperature begin to converge in the dense gas with optical depth ≫ 1, and the gas radiation flux roughly balances the change in gas energy from accretion -but this really only applies to the most dense, nuclear midplane gas.And even that gas features rapid cooling ( cool ≪  dyn ), and is primarily atomic with abundances of important species such as H − and free electrons far from the naive Saha equation estimate, owing to the large deviations from LTE and large contribution to ionization and reactions from terms like external irradiation, shocks, cosmic rays, and other processes described above.

STRUCTURE OF THE TORQUES AND STRESSES
WITHIN THE DISK We now explore the origins and physics of the torques/stresses/angular momentum transfer processes in these disks in more detail.

The Torques and Net Angular Momentum Change of the
Gas First, consider the torques themselves.We directly compute and record the "true" in-code instantaneous specific torque  ≡ r × a on every Lagrangian gas element in the simulations (where r is the vector distance from the center and a is the -Quantitative properties of the torques from Fig. 24.We plot the instantaneous value of the torques ( ≡ r × a) in cylindrical annuli (restricting to gas with |  | < 0.3 ), divided into two components (separating the acceleration a into said components as they are operator-split in-code): the "MHD" torque (all terms solved in the Riemann problem including magnetic+kinetic+thermal stresses), and the "gravitational" torque (from a grav = g).Other terms ( § 7.1) such as the torque from radiation pressure, cosmic rays, or microphysical viscosity fall below the plotted range and are negligible here.As in Fig. 24 we take the component of the torque in the angular momentum [AM] direction  • ĵ (so values < 0 [solid] indicate AM loss, while values > 0 [dotted] indicate AM gain) , and normalize to  2  ( ) so that a value ∼ 1 indicates an order-unity change in AM in a dynamical time.We show the mean (line) in each annulus and 90% range of individual cell values (shaded).While Paper I showed gravitational torques dominate at ≳ pc scales, MHD torques clearly dominate here on smaller scales within the disk, and generate net angular momentum loss within the disk.There is large cell-to-cell scatter in the instantaneous torque, consistent with the fluctuating picture in Fig. 24.We compare (green dashed) the instantaneous net inflow rate /3    Σ gas measured (Fig. 5): this should be ∼ |  |/ 2  for a thin, homogeneous, steady-state, axisymmetric disk on slowly-decaying circular orbits with constant .Given the large deviations from these assumptions the two agree fairly well.In a local, instantaneous sense, there is a broadly comparable volume with  • ĵ > 0 and  • ĵ < 0, though of course the mean torque causes angular momentum loss.This must be the case given the steady inflow and inwards radial flow seen above, and is consistent with the mean Reynolds and Maxwell-type stresses as we show below.Fig. 25 plots the azimuthally-and-vertically-averaged torques as a function of cylindrical radius , now separating by components.We again in-code separate different contributions to the torques, so this is decomposing exactly the contributions to the torque seen by the gas cells.Specifically, we divide the acceleration (and therefore resulting torque) into several terms: the gravitational term (from all gravitational forces), the "MHD" term (from all MHD forces), the radiative term (from radiation pressure/photon momentum), the cosmic ray term (from cosmic ray pressure forces/scattering) and the viscous term (from physical molecular, atomic, and Spitzer-Braginskii viscosity).Since the computation of these terms is operator-split in GIZMO, their decomposition is straightfor-ward here.21The dominant component on-average at  ≪ pc is clearly the "MHD torque," which includes all the terms from the Riemann problem solved in-code (essentially, the sum of torques from magnetic stresses, thermal pressure gradients, kinetic/Reynolds stresses or winds/outflows, etc).This shows a very large variance and (per Fig. 24) often local sign changes, but the net torque, azimuthally averaged, is negative (i.e.causing a net loss of angular momentum).The variation of the instantaneous torques in time (whether we consider a given location x or Lagrangian parcel) is comparable to their variation in space shown in Figs.24-25, and the variation of the spatially-averaged torque at a given annulus  is a factor of a few, as the various "bumps" and "dips" seen in Fig. 25 appear and disappear.
Given these broad fluctuations, the range of torques are consistent with the measured instantaneous inflow rates at a given time and radius, to within a factor of a few.But since, as noted above, neither is in exact steady state and both show inhomogeneity in space and time, we do not expect perfect agreement.The agreement improves of course, as it must, if we average over time as well as space.
The "gravitational torque" (torques arising from the gravitational forces themselves directly, in a non-axisymmetric potential) are much weaker but the next-strongest on average.At large radii ≳ pc, these torques actually dominate, and they are extensively discussed in Paper I as a result (see also Anglés-Alcázar et al. 2021).However, a dominant gravitational torque relies on (a) the non-BH contributions to the potential being non-negligible (i.e.terms that could conceivably be non-axisymmetric), and (b) the presence of a dominant collisionless (in this case stellar) disk, with mass much larger than the gas disk at a given radius, which can efficiently torque the gas disk (thus giving rise to torques that are stronger than would appear with a gas disk alone).The cessation of efficient star formation at scales ≲ pc therefore produces the transition to MHD torques dominating.In addition to being relatively weak, the gravitational torque is quite smooth (being dominated by the external  = 1 mode of the stellar disk at ≳ pc scales) and flips sign (so actually spins up the disk) at sub-pc scales (as predicted and discussed in Paper I).So it cannot be the dominant source of inflow.
The radiative torques from radiation pressure forces are much smaller at all radii plotted.Also not shown because they fall entirely off the plot, but computed nonetheless in our simulations, are the torques from the micro-physical Spitzer-Braginskii and atomic/molecular viscosities (many orders-ofmagnitude smaller than those here) and anisotropic cosmic ray forces.
It is also easy to verify -and indeed, is expected due to the very low values of  ≪ 1 in the disk -that the thermal pressure gradient contribution to the MHD torque, is negligible.So the net torque inside the disk is dominated by some combination of the magnetic and kinetic stresses.Given that the torques are dominated by a combination of magnetic and kinetic stresses, we now turn to examine those stresses directly.Recall, the gas momentum equation can be written  (v)/ + ∇ •  internal = S +  g, with S being a source term representing e.g.non-hyperbolic terms from CR and radiation partial-coupling (which are small in the disk region of interest) and g = −∇Φ the gravitational acceleration. internal ≡  kin +  mag +  therm +  visc +  cr +  rad represents the usual pressure tensor decomposed into kinetic ( kin ≡ vv), magnetic ( mag ≡ (|B| 2 I/2 − BB)/4), thermal/"hydrodynamic" ( therm =   B  I), viscous ( visc ≡ ( visc /3) (3 B B − I) (3 B B − I) : (∇v)), cosmic-ray ( cr ≡ ∫ p cr v cr (p cr )  cr (p cr )  3 p cr ) and radiation ( rad ≡ ∫  rad,  3 D  d) components respectively.As shown in Figs. 4, 5, 25, and discussed further in § 7.1, the thermal/hydrodynamic, viscous, cosmic-ray, and radiation stresses, as well as the source terms S, are small compared to the leading-order magnetic  mag and kinetic  kin terms (and of course gravity).Fig. 26 therefore plots the average value of each component of  kin ≡ vv and  mag ≡ (|B| 2 I/2−BB)/4 in radial annuli around the BH, where again we directly extract the stress tensor from the simulation at each time.
As expected (per Fig. 20), within the disk (≲ 0.1 pc), the rotational/centrifugal/angular momentum term ⟨ kin   ⟩ ≡ ⟨  2  ⟩ dominates, and provides the dominant support versus the radial gravitational force from the SMBH.The next-most-prominent term, at least in the inner disk, is the azimuthal/toroidal magnetic term ⟨ mag   ⟩ ≡ −⟨ 2  ⟩/8, which, as we showed above, dominates the internal disk pressure and provides most of the vertical support (Figs. 5 & 10).Then there are a group of broadly order-of-magnitude comparable terms (as anticipated from Figs. 10 & 20) including the , , , and  terms.These reflect the not-extremely-anisotropic fluctuating velocity and magnetic field terms seen in Figs. 10 & 20, and of course the  and  terms are just dominated by the normal radial and vertical dispersions.The  term is often significantly weaker in the disk.
Of course, these are the total stresses, so for example ⟨ kin   ⟩ = ⟨     ⟩ = ⟨⟩⟨  ⟩⟨  ⟩ + ⟨    ⟩ includes both the "mean field" term ⟨⟩⟨  ⟩⟨  ⟩ and the "fluctuating" or traditional Reynolds stress term ⟨ kin   ⟩ ≡ ⟨    ⟩ (where   ≡   − ⟨  ⟩).⟨ kin   ⟩ is negative therefore, as expected in any system with net inflow, because it is dominated by the mean term with ⟨  ⟩ > 0 (the rotational motion) and ⟨  ⟩ < 0 (inflow).We therefore also plot the separation of each component of  into fluctuating components, which allows us to more clearly see (1) the crudely isotropic turbulent fluctuations,22 (2) the fact that the Reynolds stress ⟨    ⟩ and both the total and fluctuating Maxwell stresses are almost always positive (indicating angular momentum loss, given our sign convention), and (3) that in the rotationally-dominated disk at ≪ 0.1 pc, the dominant torque should come from the traditional Maxwell stress ⟨ mag   ⟩.As discussed below ( § 7.2.2), in the inner disk the Maxwell stress is itself dominated by the mean component, but with non-negligible contribution from the fluctuating (⟨−    ⟩/4) component.The magnitude of the fluctuating components in both cases is 22 As discussed in § 5.3, if we subtract the best fit  = 1 component from each   (, ) to attempt to remove the effects of coherent eccentric motion, this has at most a modest (order-unity) effect reducing the kinetic ⟨  kin  ⟩ and ⟨  kin   ⟩, and almost no effect on the most relevant Reynolds stress ⟨  kin   ⟩.
consistent with the properties of the turbulence discussed in § 5.4 (e.g.⟨    ⟩ ∼ ⟨⟩ ⟨ 2 turb ⟩).Briefly, it is worth noting that the sign and efficiency of the Maxwell stresses are expected here, given the strong anticorrelation between   and   , as demonstrated in § 4. Recall that this results from the simple fact that the toroidal field is supplied by radial fields in the disk plane (given the induction equation; see Fig. 15 and  § 4.2), with simple fluxfreezing/advection and trans-Alfvénic turbulence explaining the mean and fluctuating field strengths.We stress that this is distinct from some historical models wherein   and   are sourced from some strong mean poloidal field ⟨  ⟩, in which case such an anti-correlation (ensuring inflow) is non-trivial.For the Reynolds stress in the outermost disk where it dominates, it is less obvious what exactly regulates the detailed quantitative properties of the stress, related to our discussion in § 5.4 regarding the uncertainty in what exactly drives the turbulence.If instabilities related to self-gravity like the global  = 1 modes and/or gravito-turbulence drive the turbulence, these would naturally produce the kind of in-plane motion with the correct sign (on average) of ⟨    ⟩, but other instabilities that could be present, such as the low- extensions of the MRI (see § 4.2.3), have not yet been well-studied in the non-linear regime.
For completeness, we note if we define the comoving MHD stress tensor T ≡  mag +  kin , there are three terms in Fig. 26 that can in principle give rise to torques in the disk plane.23First, the usual Maxwell/Reynolds stress (T   ), which we discussed above and will discuss further below.Second, the non-axisymmetric azimuthal term which can produce a local torque ∝ T   /.And third, the wind or convective term ∝ T   /.The azimuthal T   / term is usually neglected (even when |T   | itself is large) because by definition the net torque integrated over  must vanish at leading order for orbits that are approximately closed and/or circular (assuming small deviations from axisymmetry in the potential) in a disk which is not evolving on a timescale fast compared to the orbital time (Kalnajs 1971).However, Hopkins & Quataert (2011) showed that this term can be leading order in the net torque when the "gravitational torques" discussed above are important (or when non-axisymmetric terms in the potential become non-linearly large), because one can break the periodicity in  that causes the integral of T   / to vanish if the external potential induces orbit crossings and shocks or pileups in the gas.So we note it briefly here because at radii ≳ 0.1 pc, where the gravitational torques can be important (and where we see in e.g.Fig. 17 that the strong eccentric pattern produces a clear discontinuity in the velocity streamlines), this can actually play a leading-order role.However, in the primarily rotationallysupported inner disk at ≪ 0.1 pc, it becomes less important.
The vertical/convective term in T   / is usually considered when there is a strong outflow/wind removing angular momentum from the disk.But recall here (Figs.17-20) that the vertical velocity is primarily (weak) inflow.This means that in order to significantly lower the specific angular momentum of disk material and promote accretion, the inflowing gas would have to (1) carry a considerable fraction of the total accretion rate (which it does not, because we showed above the density structure means most of the inflow is through the gas in the midplane with || ≲ ), and (2) be significantly sub-Keplerian (which it is usually not, since by definition it tends to join the disk at its circularization radius).
So it seems clear that, on average, the most important MHD torques in the inner disk indeed arise predominantly from the usual Maxwell+Reynolds stress.

Relative Importance of the Maxwell vs. Reynolds Stresses, and Fluctuating vs. Mean-Field Components
In Fig. 27, we therefore plot the azimuthally-and-vertically averaged Reynolds stress    ≡ ⟨ kin   ⟩ ≡ ⟨ (  − ⟨  ⟩) (  − ⟨  ⟩)⟩ ≡ ⟨     ⟩ and Maxwell stress −   ≡ ⟨ mag   ⟩ ≡ −⟨    ⟩/4, as a function of radius from the BH, with their range and sign.Because the Maxwell term ∝ ⟨    ⟩ includes both a mean-field ⟨  ⟩⟨  ⟩ and fluctuating field ⟨    ⟩ component, we also plot both the total Maxwell and fluctuating component alone.We can immediately compare the normalization of these terms to the torques  in Fig. 25 noting that, if dominant, the sum of these terms in the units given (  2  ) should approximately equal the torque in units of  2  (up to an order one constant that depends on the averaging weighting and slope of the different terms).We see reasonably good agreement, with a factor of < 2 difference owing to how the different weighting of the means.And of course the magnitude of the stresses are much larger than in a SS73-like  disk where (by assumption) the Alfvén and turbulent velocities are much smaller than the thermal sound speed.
Next, we compare and see that at large radii ≫ 0.01 pc, the Reynolds stress is dominant by a factor of a couple to a few, while at smaller radii the two are comparable or the Maxwell stress is dominant by a similar factor.To understand this better we recall Fig. 22, which plots the Mach number M of the turbulence versus radius.We see that it is quasi-isotropic, trans-Alfvénic, and highly super-sonic (consistent with all our previous analysis of  and velocity fields directly).We note here though that although the trend of Alfvén Mach number M  with radius is weak, we see the turbulence transition from mildly sub-Alfvénic at the smallest radii, to modestly super-Alfvénic at large radii.This is broadly consistent with the trend in the ratio of Reynolds to Maxwell stresses., is what appears in the angular momentum (AM) equation, while ⟨    ⟩ includes additional pure advection terms that do not influence the AM (i.e.only ⟨     ⟩ gives rise to AM exchange).Each quantity is averaged over several snapshots around the time of interest and in annuli within the disk, excluding any dense star-forming clumps.The shaded range shows the ±1 range of | −     /4  |, etc., in all cells in each annulus.Reynolds stresses tend to be larger at radii ≳ 0.01 − 0.1 pc and especially outside the ordered disk, where the turbulence is weakly super-Alfvénic (Fig. 22).Within the inner disk, Maxwell stresses dominate by a factor of ∼ 10.The total/mean Maxwell stress is larger than the fluctuating component by factors of several, consistent with the strong mean fields in e.g.Fig. 15 but distinct from predictions of e.g. the weak-seed-field MRI.We also (top) plot the instantaneous inflow rate /3    Σ gas from Fig. 25.This agrees well with the sum of total Maxwell+Reynolds stresses, as expected if they dominate the angular momentum transport.
We also see that the although mean field Maxwell stress is usually dominant, the fluctuating field stress is non-negligible, especially at some times at very small radii.It will occasionally occur that an excess of toroidal magnetic flux will build up in the center, sometimes briefly suppressing accretion but then producing a large mean-field Maxwell stress and leading to a rapid accretion event.This leads to the fields being accreted inwards, changing the ratio of mean-to-fluctuating stresses briefly.But in general the behavior plotted holds, showing that the fluctuations   and   are sufficiently strongly anti-correlated that they can be comparable to the mean-field stress, even when ⟨|  |⟩ ≲ |⟨  ⟩|.
The "origins" of the angular momentum transport and stresses, therefore, are directly tied to (a) the origins of the strong magnetic fields, and (b) the origin of the turbulence within the disk, each of which was discussed above.
8. COMPARISON TO PREVIOUSLY-STUDIED "STRONGLY-MAGNETIZED" DISKS While the majority of the literature on quasar accretion disks has focused on disks with "weak" magnetic fields (e.g.magnetic pressure less than radiation or thermal pressure), there has been some discussion of disks in the strong-field limit.These can generally be divided into a couple of different categories, some of which are indeed closely related to the behaviors we report here in our "flux-frozen" (and flux-fed) disks, and some of which are not.We therefore find it helpful to briefly review these past models and distinguish the results from their predictions.Note that a more quantitative comparison to SS73-like "weakly-magnetized" ( ≫ 1) disks is given in Paper III (see also § 3.2 & 4.2.4).

Magnetically Arrested (MAD) Disks
Magnetically arrested (MAD) disks, in which accretion is halted by strong magnetic fields near the BH, are usually characterized by extremely strong poloidal magnetic flux (Bisnovatyi-Kogan & Ruzmaikin 1976;Narayan et al. 2003).We clearly see that in many (though not all) respects, the behaviors here are opposite those predicted in the MAD regime and seen in idealized simulations of MAD disks (e.g.Tchekhovskoy et al. 2010Tchekhovskoy et al. , 2011;;White et al. 2019;Xie & Zdziarski 2019).Most importantly: (1) accretion proceeds rapidly here, (2) the accretion is in fact aided by magnetic fields, (3) the fields are primarily toroidal, which has qualitatively different consequences at these scales, and (4) the disks are vigorously turbulent and cool efficiently (which also helps promote strong inflows, compared to the more well-studied MAD limit).This is not surprising, as it is trivial to verify that, everywhere we resolve, the magnetic field strengths are below the MAD limit of ⟨  ⟩ 2 mad /8 ∼   BH Σ gas /4  2 (this limit is effectively equivalent to  mad ,  ≳  K , for the poloidal field, but we see relatively weak poloidal fields and even including toroidal components we see   <  K everywhere).This is true both in the simulations here, and the analytic models in presented in Paper III which can be extrapolated to ISCO scales.As argued therein, if one starts from an outer disk boundary like that in our simulations, then even if we ignore the important geometric constraint that one cannot generate a dominant mean vertical field via flux freezing from field configurations like those in our simulations, it would require much stronger amplification of the mean field compared to what we see (  ∝ ⟨|B| 2 ⟩ ∝   with  > 2, as opposed to  ∼ 4/3 which we measure in the simulations here).In addition, the turbulence would have to be strongly suppressed to become become highly sub-Alfvénic ( turb ≪   ), in order to satisfy a MAD-like criterion  2 /8 ≳   BH Σ gas /4  2 as  → 0.
However, the argument that a disk should eventually enter the MAD limit somewhere outside the ISCO or event horizon is often phrased in terms of an incoming magnetic flux from much larger scales, as e.g.Φ mag ∼   2  ≳ Φ crit ∼ G pc 2 ( bh /10 7 M ⊙ ) 3/2 (  BH /  Edd ) 1/2 using the scalings from Xie & Zdziarski (2019).At a glance (comparing Fig. 10) it would naively seem that our simulations (and most observed AGN) exceed this limit at ISM scales.But we caution that this flux-based extrapolation makes several key assumptions which are not valid in the simulations here.Most importantly, it assumes (1) that the "seed" field is dominated by a coherent poloidal/vertical/dipole field (uniform B ≈ ⟨  ⟩ ẑ with coherence length ∼ ), which then (2) is amplified much more rapidly than other field components, assuming strictly homologous, laminar accretion with B = ⟨  ⟩ ẑ, so that ⟨  ⟩ ∝ 1/ 2 , with (3) a thermal-pressure dominated Shakura & Sunyaev (1973)  disk assumed to estimate Σ gas and the disk thickness / and pressure support level "needed" to arrest the disk, together with (4) assumed radial/toroidal field components |  |/|  | ∼ |  |/|  | ∼ / ≪ 1 that are sourced only via turbulence from the dominant poloidal field.But none of these conditions are satisfied in our simulations.As shown above, at large radii the "seed" flux is isotropic and turbulent (with |  | ≫ ⟨  ⟩; Figs. 3, 5, 6 & 10).Moreover in the "free-fall" region where gas is tidally captured by the SMBH and then circularizes to form the disk (see § 4.1), the structure of the tidal compression/expansion means that any mean poloidal ⟨  ⟩ is preferentially suppressed relative to the dominant ⟨  ⟩ or ⟨  ⟩ components, ensuring that |B| ≫ |  | ≫ ⟨  ⟩ in the disk (e.g.Fig. 10).Per § 5.6 and discussion in Lubow et al. (1994), this means that the vertical field   will be "locked" in a regime where turbulent flux transport and resistivity suppress the amplification of ⟨  ⟩ (qualitatively unlike   and   , which grow via flux-freezing/advection rather than being sourced from ⟨  ⟩), explaining why the mean ⟨  ⟩ grows much more slowly than ⟨  ⟩ and ⟨  ⟩ as  → 0 (opposite the MAD assumption).And we have shown that assuming a thermal-pressure-dominated ( ≫ 1) Shakura & Sunyaev (1973) disk would predict orders-of-magnitude different disk properties from those in our simulations (with e.g./ and Σ gas incorrect by factors of ∼ 300 and ∼ 10 4 , respectively; see Fig. 5).
Thus it is obvious that our disks are not magnetically arrested in practice, nor should they be given the physical conditions; however, it is certainly possible that if the densities of the inflows (hence accretion rates) dropped sufficiently at some later time, the system might transition to a MAD-like state (depending on whether the magnetic fields also declined, and whether or not the turbulence became less vigorous as the densities and accretion rates declined).It is also conceivable that very close to the BH (around the ISCO) the behavior of these disks could become more "MAD-like" (or otherwise truncated or inefficient at low accretion rates, see Hogg & Reynolds 2018;Datta et al. 2022).But again, the simple analytic models we develop in Paper III suggest that this would require some qualitative change in the fundamental scalings for the turbulence and magnetic field strengths, compared to the resolved behaviors in our simulations.

Magnetically Elevated or Levitated Disks
Magnetically levitated and/or elevated disks are disks in which magnetic fields are relatively weak in the midplane (with  ≫ 1), but become fractionally more important with  ≲ 1 in the tenuous gas at a few scale heights above the midplane (|| ≫ ).This can produce a variety of interesting behaviors with e.g.inflow/outflow along current sheets or angular momentum transport via MHD winds, but the bulk of the mass is still effectively a "classical" thermal-pressure-dominated  disk (Johansen & Levin 2008;Gaburov et al. 2012;Sądowski 2016;Mishra et al. 2020).Again our simulations are in an opposite regime: (1) most importantly we see  ≪ 1 in the midplane and in a gas-mass-and gas-density-weighted mean sense; (2) the dependence of  on both gas density ( is lower in denser gas) and scale height ( is weakly stratified and often increases in the midplane) is the opposite of that in a magnetically elevated disk; (3) the thermal structure of the disk is also opposite the elevated disk model (it is warmer above the midplane); (4) the velocity structure above the disk is also distinct.

Vertically Magnetized Star-Forming/Galactic Disks
Recently, Begelman & Silk (2023) suggested that stronglymagnetized galactic disks could drive BH accretion in protogalaxies on scales much larger than the traditional accretion disk.They focus primarily on galactic radii ≳ pc, where they assume a magnetic field that is dominated by a mean global vertical field (with smaller, turbulence-resistivity-dominated radial and toroidal components sourced from ⟨  ⟩), turbulent  ∼ 1 maintained by stellar feedback, and magnetic torques that dominate on these scales with turbulence that is sub-Alfvénic and/or sub-sonic relative to the dominant mean vertical field (specifically they require . While these larger scales are primarily discussed in Paper I, we note here that we do not see these conditions in our simulations.At star-forming galactic radii ≫ pc, Paper I showed that (1) magnetic forces are subdominant to gravity and turbulence/bulk motions; removing magnetic fields has almost no effect on the torques/inflow rates at  ≫ pc, which are dominated on these scales by a combination of gravitational torques and stellar-feedback-induced shocks; (2) there is no dominant coherent/uniform mean vertical field (consistent with observed galactic magnetic fields; Mao et al. 2010;Beck 2015;Jaffe 2019;Ma et al. 2020;Krause et al. 2020); and (3) the turbulence at  ≫ pc is super-Alfvénic.In the accretion disk studied here (radii ≪ pc), we also see a very different situation to that assumed in Begelman & Silk (2023), with fields that are primarily toroidal and set by flux freezing ( § 5.6),   sourced by   ( § 4.2), weak and incoherent vertical fields in the disk ( § 4), ⟨ ,  ⟩ ≪  turb (by factors ∼ 10 − 100; Figs. 10 & 20), thermal+magnetic  ≫ 1 and minimal star formation ( § 3), and qualitatively different scalings with radius.

Previous Models of Toroidal-Field-Dominated Disks
In previous literature, the models that come closest to capturing the properties we observe in our simulations are those of Begelman & Pringle (2007); Oda et al. (2009).These have been studied numerically in idealized simulations, which consider a relatively small section of the disk with simplified physics and fixed initial and boundary conditions, in e.g.Salvesen et al. (2016a); Kudoh et al. (2020) (see also Johansen &Levin 2008 andGaburov et al. 2012, though their assumptions, with  ∼ 1, may be more similar to magnetically elevated disks).The main similarity is that these models posit  ≪ 1 from a primarily toroidal magnetic field, which dominates over the midplane gas+radiation pressure.As we discuss in more detail in Paper III with a simple analytic model motivated by our simulations, it turns out that this similarity alone is sufficient to capture most of the crucial properties of the simulations here.
Nonetheless, there do still exist qualitative differences between the behaviors seen in our simulations and those models.Most importantly, those models implicitly assume the field is amplified from initially small values via the MRI and produces a toroidal   which (a) is dominated by its fluctuating components, and so (b) rapidly changes sign, even following a given Lagrangian parcel, (c) is dominant over the poloidal   (which sources it in the first place) but by a relatively small factor, and (d) saturates at a value of   ∼ √    K , above which the linear growth of the MRI is assumed to be suppressed following the analytic analysis in Pessah & Psaltis (2005).The final condition leads (e) to relatively modest  ∼ 0.1 at the radii of interest here (much higher than the values we see).This also leads (f) to the prediction that the maximum accretion rate that can be maintained is only just about the Eddington limit around supermassive black holes, so these studies focused on much lower-density, lower-accretion rate regimes and did not consider highly super-Eddington accretion.And (g) this means that the disks simulated here are orders-ofmagnitude more gravitationally stable (e.g.retain  ≫ 1 out to orders-of-magnitude larger radii from the BH) compared to the predictions in these studies.
There are other differences as well that could be important: the models in Begelman & Pringle (2007) (as well as Oda et al. 2009) made very different assumptions for the temperature and opacity structure and predict a stratified disk with a hotter midplane; however the actual simulations in Kudoh et al. (2020), which include dynamical cooling and heating (albeit with a simplified prescription compared to the detailed network here) predict an opacity and thermal structure much closer to what we see (inversely stratified, with a mostly-atomic cool midplane at these radii).Moreover these previous studies neglect the fact that these models, almost by necessity, predict highly super-sonic (M  ≫ 1) accretiondisk turbulence at high , which in turn relates to very efficient/rapid cooling ( cool / dyn ∼ M −2  ≪ 1).This in turn leads to some star-formation but avoids catastrophic gravitoturbulent fragmentation via magnetic support.Indeed, Begelman & Pringle (2007) essentially make the Shakura & Sunyaev (1973)-style assumption of a laminar disk with turbulent velocity  t ∼  1/2   ≪   ≪  c .And these analytic/idealized models (including ours, in Paper III) all assume quasi-circular disks, neglecting the potentially important role of the coherent eccentric, large-amplitude  = 1 modes we see here.
Fundamentally, the key physical difference is that our simulations do not need to amplify B from some weak/trace "seed" field via the MRI in order to achieve their "magnetically dominated" state.Rather, they begin from this state, as the accreted gas carries in sufficiently large magnetic flux, with the consequence that the gas initially in the disk already has a magnetic field well above the specific (analytic) saturation threshold for the MRI assumed in Begelman & Pringle (2007); Oda et al. (2009), namely  , flux−freezing ≫ (   K ) 1/2 .This, in turn, produces a variety of other consequences as described above, as well as different instabilities operating in the disk.
That said, some key conclusions are robust and confirmed here: the fact that  ≪ 1 can be supported down to the ISCO, in principle; that the fields do not decay; that these can produce accretion rates far larger than a thermal Shakura & Sunyaev (1973) disk; that they can sustain super-Eddington accretion; that the disks are much more stable at large radii than an equivalent Shakura & Sunyaev (1973) disk; and that the outer disk is primarily atomic and "cool" and thermally weaklystratified or even inversely stratified.These are all at least qualitatively robust conclusions comparing to the simulations in Kudoh et al. (2020) and analytic arguments in Begelman & Pringle (2007), even if the details and origins of the magnetic fields differ in some important respects.

Why Does the Strong Field Not Decay (or Buoyantly
Escape)? Related to the discussion above, there have also been some previous claims in the literature (Salvesen et al. 2016a;Fragile & Sądowski 2017) that a strong toroidal-field dominated disk can, under the right conditions, rapidly (in tens of orbits) evolve to  > 1, due to a combination of adiabatic expansion with outflows and buoyant escape.Clearly we do not see this "decay," even having run our simulations for ≫ 10 4 orbits at their innermost radii (see e.g.Fig. 14).It is also worth noting that the simulations of Kudoh et al. (2020) discussed above also did not see any such decay/escape.Moreover as discussed extensively above, we have checked not just that the field is maintained at a given Eulerian position, but following a Lagrangian parcel over time, we see it amplified (not decaying) according to simple theoretical expectations.
There are a several important and straightforward reasons why the arguments for the specific situation considered in e.g.Salvesen et al. (2016a); Fragile & Sądowski (2017) should not hold here.Importantly, those authors considered a qualitatively different parameter space, initial conditions, and boundary conditions from those here.In those papers, the initial "disk" is strictly hydrostatic, with no net accretion or torques, and rapidly expands/puffs up, producing outflows almost everywhere (even in the disk midplane) and weakening the magnetic fields primarily via adiabatic expansion.This is almost the exact opposite of the behavior here, where we see |B| increase because gas flows in becoming more dense .Moreover, the density and temperature scales are ordersof-magnitude different, and those idealized simulation models were strictly adiabatic -i.e. had no cooling -while we see very efficient cooling  cool ≪  dyn .This can both maintain low- and dissipate the thermal energy, which in the idealized models could cause the disk to puff up or drive strong buoyancy instabilities (see also Paper III).And they also only considered quite modest initial  > 0.1, so already close to  ∼ 1, while our disks essentially "begin" at orders-of-magnitude lower  (Fig. 4-5) via the transport of magnetic flux from a weakly magnetized ISM .This also qualitatively changes which instabilities can operate in the disk (as shown in Fig. 16): for example as discussed in § 4.2.3, the specific Parker-like vertical buoyancy modes discussed in Johansen & Levin (2008); Salvesen et al. (2016a) may operate on such large wavelengths (∼ 10  ≳ , given the large scale heights in Fig. 5) that they cannot fit within the disk, or they may be suppressed by the combination of strong trans-Alfvénic turbulence, nonnegligible radial and vertical fields, and low- (Horiuchi et al. 1988;Kim & Ryu 2001), but other radial buoyancy modes may appear which have qualitative different effects on the field ( § 4.2.3).
Moreover, both Salvesen et al. (2016a) and Fragile & Sądowski (2017) considered exclusively azimuthal fields, and noted that adding some poloidal or radial component would prevent the field decay, as was subsequently shown explicitly in Salvesen et al. (2016b).Here, we do see dominant toroidal fields, but with substantial (order-unity) projections into both the poloidal and radial directions (Fig. 10).Indeed, this is inevitable in our simulations, given the origin of the "seed" fields for the accretion disk from a quasi-isotropically-tangled field in the ISM.Thus our simulations actually reside "safely" in the regime where Salvesen et al. (2016b) argued they should not experience rapid decay/loss/escape.
It is also crucial to note that there is no source of magnetic flux in the idealized simulations of Salvesen et al. (2016a); Fragile & Sądowski (2017).Absent any source of new flux in the midplane, it seems plausible that escape of toroidal field from the midplane and/or turbulent resistivity could eventually weaken   , perhaps producing something more like the magnetically elevated disks in Johansen & Levin (2008) with a minimum in  at the midplane.And we are not arguing that magnetic buoyancy is totally negligible here -in contrast, in § 4.2.3, 5.5, & 6 we argued that buoyancy instabilities unique to low- disks might play an important role driving the turbulence in the disk and explaining the weak stratification we observe.However, we showed above in e.g.Fig. 15 and  § 4.1 that even if vertical buoyancy removed midplane toroidal field on the fastest possible timescale it can operate (a few Ω −1 ), this would be balanced by the growth of toroidal field just from new radial flux carried in with the midplane accretion flow (see also Shibata et al. 1990, who make a similar argument from both analytic considerations and idealized MHD simulations).
In brief, the midplane mean toroidal field in the simulations is constantly sourced by advection of new radial and toroidal flux, "closing the dynamo loop" and replenishing/maintaining the magnetic-field strength.We stress that this possibility was indeed anticipated by Salvesen et al. (2016a); Fragile & Sądowski (2017), as well as others such as Shibata et al. (1990); Johansen & Levin (2008); Kudoh et al. (2020), all of whom emphasized the critical importance of physically-motivated boundary conditions for the accretion disk from larger (ISM) radii and therefore the source of flux (the motivation for our simulations in this paper).Salvesen et al. (2016a) specifically stated that, given "favorable conditions" for the source of magnetic flux along with the accretion flow (which they noted could arise from the "external" field from the ISM/galactic scales), "a strongly magnetized disc would necessarily follow."Though they envisioned primarily poloidal external flux, while we find a mix of radial and poloidal flux (see Fig. 10), the overall effect will be similar, maintaining strong toroidal fields in the disk (Fragile & Sądowski 2017).Therefore, taken together, the stable behavior of the strong toroidal fields here is expected, and does not contradict the results from the more idealized test-problem simulations considered in Salvesen et al. (2016a); Fragile & Sądowski (2017).

ANGULAR MOMENTUM OF THE DISK RELATIVE
TO THE PRE-EXISTING BH SPIN Our simulations do not have sufficient resolution to follow the spin of the SMBH in detail.However, even in the lowresolution "progenitor" simulation, we do dynamically follow the total angular momentum (AM) accreted by the SMBH j BH (from its formation as a seed, with j BH updated whenever some gas is accreted following the detailed numerical description in Hopkins et al. 2023b).While imperfect, ĵBH should serve as at least a plausible guess for the spin direction of the SMBH.Given the scenario occurring in our simulation ( § 3) -where a turbulent massive star-forming cloud complex in a highly chaotic, clumpy high-redshift massive galaxy is partially tidally disrupted on close passage to a pre-existing SMBH in the nucleus -there is no reason to think that the accretion disk AM vector j gas should be preferentially aligned with this pre-existing BH AM vector (or the SMBH spin direction).
We can immediately verify this in the simulations.Compar-ing the vector directions ĵgas and ĵBH , the accretion disk at the times we follow is approximately ∼ 140 • misaligned from ĵBH interior to ≲ 0.1 pc (where the disk AM direction ĵgas is quite stable within different annuli, to within a few degrees).This means it is likely retrograde and out-of-plane with the preexisting SMBH spin.As the disk evolves, the AM direction will evolve too, tracing material further away from the SMBH at this time: at ∼ 1 pc (or ∼ 10 pc) the mean AM vector of the accreting gas is ∼ 100 • (or ∼ 130 • ) misaligned, so can vary by ∼ 40 • .The corresponding enclosed gas masses and accretion timescales at the current accretion rate into < 80 au are: ∼ (3 × 10 4 , 6 × 10 5 , 10 7 )  ⊙ and ∼ (10 3 , 2 × 10 4 , 3 × 10 5 ) yr for gas within ∼ (0.1, 1, 10) pc.
There is also no correlation between the direction ĵBH and the secondary preferred direction of the accretion disk: namely the direction of the semi-major axis in the plane of the disk (since it is coherently eccentric).But this is expected even if ĵBH and ĵgas were aligned, since the semi-major axis direction precesses on timescales of order the orbital time at the outer disk radius.
Our inner resolution scale (∼ 80 au or ∼ 300  schw ) is much larger than any scale where BH spin directly influences the dynamics (e.g. via Lense-Thirring precession), so this is not directly important for the scales we resolve.However, simulations of smaller scales using these results as outer boundary conditions should include such effects.Moreover, this is consistent with the well-established observational result that the "spin axis" inferred from AGN jet directions does not appear to correlate with gas AM on any resolved macroscopic scales in galaxies (see e.g.Schmitt et al. 1997;Kinney et al. 2000;Hopkins et al. 2012a,c;Davies et al. 2014;Reynolds 2021, and references therein).Indeed, the misalignment angles that we find, and their variation with scale and time on sub-pc scales, are consistent with those seen in previous simulations by Anglés-Alcázar et al. (2021), who also emphasize the lack of correlation between the angular momentum vector at ≲ 1 pc (at any given instant) and on much larger (≳ kpc) scales.
Together, this further suggests that extremely misaligned accretion -which can itself even further help in promoting very rapid accretion on scales smaller than those we resolve here (see e.g.Kaaz et al. 2022) -might be quite common in quasars whose large-scale fueling episodes resemble that here.
10. WHAT HAPPENS WITHOUT MAGNETIC FIELDS?Paper I discusses several numerical experiments in detail, including an example where we begin our "hyper-refinement" stage of the simulation without magnetic fields (but from otherwise identical initial conditions).There we showed that this produces relatively weak effects on ≫ pc scales (in the galactic ISM), consistent with the vast majority of previous simulation studies of magnetic fields in cosmological galaxy formation simulations (Su et al. 2017(Su et al. , 2018(Su et al. , 2019;;Hopkins et al. 2020b;Ji et al. 2020;Steinwandel et al. 2019Steinwandel et al. , 2022;;Martin-Alvarez et al. 2021;Ponnada et al. 2022;Whitworth et al. 2022).But as expected from our arguments above, the effects on the disk at scales ≪ pc can be dramatic.Here we explore the disks that form without MHD in more detail.
Briefly, we note the resolution and run-time of this simulation.We initially run the simulation with our "default" resolution of the full-physics simulation with MHD (mass resolution ∼ (1 − 5) × 10 −3  ⊙ in the high-resolution region inside ≲ 10 pc), but for a more limited time (equivalent to ∼ 3000 Ω −1 inner or ∼ 100 yr at the very highest resolution level, and ∼ 10 5 yr at intermediate resolution following the infall during refinement; see Paper I for details of the refinement scheme).This reduced time is chosen both because it is a counterfactual numerical experiment and because it is notably more numerically expensive compared to our default simulation, owing to the much more rapid fragmentation into extremely dense sub-clumps, the formation of extremely massive stars from said fragmentation (see Paper IV), and the razor-thin disk.This razor-thin disk, discussed below, means that even at our extremely high resolution, it is challenging to resolve the vertical scale-height of the inner disk at  ≪ 0.01 pc.We therefore re-start this simulation with an even further layer of refinement continuing to a minimum mass resolution of Δ min ≈ 0.003 M ⊙ (/0.1 pc) 2 at  < 0.1 pc, reaching a highest resolution of Δ < 10 −7 M ⊙ at our innermost radii < 100 au.This is extremely expensive computationally, so we evolve it for a much shorter time at its highest resolution (∼ 100 Ω −1 inner ), just to ensure that the properties at the innermost radii can be reliably modeled.

Runaway Star Formation
The major focus of our comparison in Paper I was to show that, absent magnetic fields, on sub-pc scales (where MHD torques take over as dominant from gravitational torques in our default simulation per Fig. 25) fragmentation and subsequent star formation run away catastrophically.Whereas in our default simulations (including magnetic fields) the fields stabilize the disk and lead to a sharp suppression of the SFR per unit area or volume at radii ≲ 0.1 − 1 pc, the volumetric SFR in a simulation without MHD continues to rise steeply as  → 0. Per Fig. 30, the SFR interior to < 1 pc (or < 0.1 pc), for example, rises from ∼ 10 M ⊙ yr −1 (∼ 0.1 M ⊙ yr −1 ) with MHD to ∼ 250 M ⊙ yr −1 (∼ 5 M ⊙ yr −1 ) without MHD.As discussed below, the SFR at < 1 pc significantly exceeds the total mass inflow rate into this annulus, strongly suppressing inflows to the BH accretion disk.
The physical reasons for this -again the main focus of Paper I -are straightforward.Most importantly, the lack of magnetic fields means the disk is no longer stabilized at radii ∼ 0.01 − 1 pc against catastrophic gravito-turbulent fragmentation.On top of this, the lack of strong magnetic stresses/torques means that gas cannot inflow efficiently so "piles up" at large radii, further accelerating fragmentation.This directly enhances the ratio of SFR to inflow rate at ∼ 0.1 − 1 pc, therefore further suppressing inflow to even smaller radii.Thus, without magnetic fields, we indeed see the classic hydrodynamic problem -reviewed in e.g.Shlosman & Begelman (1989); Shlosman et al. (1990); Goodman (2003) -of runaway rapid fragmentation in the outer disk region when inflow rates from the ISM are large.

Vastly Lower Accretion Rates
Partly as a result of the runaway star formation, but also because of inefficient MHD torques (discussed below), we see that the surviving inflow rates of gas into ≪ 1 pc are reduced dramatically without MHD.The time-averaged net gas inflow rate through our inner boundary at < 80 au is reduced from ∼ 20 − 30 M ⊙ yr −1 in our default simulation, to ≲ 0.1 M ⊙ yr −1 in the simulation without MHD (a factor of ∼ 200 − 500 reduction).But even this accretion rate appears to be falling towards the end of the relatively short duration of the rerun without MHD, because the torques in the inner disk seem to be spinning up the disk and create a growing hole in the center Fig. 28.-Face-on images showing gas surface density in a re-simulation of our fiducial simulation without magnetic fields (see § 10), as in Figs.1-2, at the latest time to which we are able to run the simulation.We see much more vigorous clumping/fragmentation at radii ∼ 0.01 − 1 pc, a more compact disk that only emerges interior to ≲ 0.01 pc, a series of tightly-wound  = 1 modes with large amplitudes, giving a "concentric ring" appearance, and a central hole in the disk that is expanding (larger than our inner accretion boundary by a factor of several).Fig. 29.-Edge-on images of gas surface density in our fiducial/"fullphysics" simulation (top) and re-simulation without MHD (bottom; as in Fig. 28).Cylindrical - coordinates are used to better see the disk structure.We clearly see that without magnetic fields, the disk is much more compact and razor thin, featuring almost no extended vertical atmosphere/halo/coronal structure.
of the disk.This can be seen by-eye in Figs.28-29.Indeed, from Fig. 30, comparing either the total inflow+outflow rate within a given annulus (  in +  out = 4 ⟨    2 ⟩), or from the total  kinetic stress (⟨ kin   ⟩ = ⟨     ⟩ ∼  K ⟨   ⟩), we clearly see that the mean radial velocity in the disk midplane is outward -i.e.we actually have a decretion disk at this time.Some gas still "leaks" through owing to non-equilibrium motions and vertical inflows joining onto the inner disk.This means the inflow rate is highly intermittent/bursty, dominated by occasional clumps that make it through the inner region.The vertical inflows contribute little in a time-averaged sense and would likely be ejected if we included some form of jets and/or harder radiation emitted by the unresolved disk interior to < 80 au, since they are primarily polar.Notably, if we evolved the disk longer, this net decretion, coupled to the larger SFR compared to gas inflow rate, would lead to even further depletion of the inner gas and therefore further suppression of the time-averaged gas accretion rate.Indeed, if we monitor the gas mass inside of ≲ 100 au, we see that after it initially rises in an inflow event, over the time period of the last ∼ 600 Ω −1 inner at  ∼ 80 au (just ∼ 20 yr) it drops precipitously by a factor of ∼ 3 (then it increases slightly as the ring of gas that builds up goes unstable and a clump falls back in towards the central < 100 au).In contrast, we showed these properties were stable in our full-physics simulations (with MHD) for timescales at least ≳ 10 5 Ω −1 inner (∼ 10 4 yr, the duration of our default simulation at highest resolution).So again without magnetic fields, we see orders-of-magnitude lower accretion rates, and it appears difficult to sustain near-Eddington (let alone significantly super-Eddington) inflows into and through the disk.
Of course it is possible that on much longer timescales (≳ Myr) the runaway star formation seen here could produce non-linearly different conditions that eventually lead to efficient accretion.But this does not change our generic conclusion that for a given set of initial/boundary conditions, magnetic fields play a critical role on these scales.That said, there is still clearly in Figs.28-29 some disk that forms.But it is visually obvious in Fig. 29 that this disk is radically different from that in our full-physics simulations.First, consider the disk size.Per Fig. 30, and by-eye in Fig. 29, we see that the surface density of the disk falls of rapidly, and the ratios / and  turb /  all increase rapidly, outside of  ≳ 0.01 pc -so the "outer" disk radius here is at least an order of magnitude smaller than in the full-physics case.This corresponds to the radii where the Toomre  ∼   /(  Σ) ∼ 140 (/0.01pc) −3/2 (Σ/10 6 M ⊙ pc −2 ) −1 (  /8 km s −1 ) falls to ≲ 10, where more catastrophic turbulent fragmentation should set in (as opposed to "gravito-turbulent" fragmentation which can at least maintain some semblance of disk structure; see e.g.Rice et al. 2005;Hopkins & Christiansen 2013 for examples) given the fact that the disk is quasi-isothermal over the limited dynamic range of radii resolved here at ∼ 8000 K.So it is not surprising that the disk can exist at these radii, given the (much lower) accretion rates seen here and its reasonably warm temperatures.
Second, we see that the disk is extremely thin, with / ≲ 0.01 through most of its resolved extent.We see subtle warps, and a roughly concentric-ring morphology owing to tightlywound modes discussed below.And not only is the disk scaleheight quite small, but whereas in our full-physics simulations we saw a fairly extended, more slowly-falling (power-lawlike) vertical atmosphere/corona above the disk (|| ≫ ; see Fig. 23), without MHD we see a much more stark (exponential or super-exponential) transition with extremely low densities outside the very narrow midplane (Fig. 30).This, of course, is not surprising, given the lack of the dominant pressure support from magnetic fields seen in our full-physics simulations, but further accelerates fragmentation and gravito-turbulence in the disk.
It is also worth noting that absent MHD, the disk is radiationpressure dominated.This is expected without magnetic pressure, given the accretion rates and other properties (see Fig. 5, or comparison to SS73-like disks in Paper III).
The  = 1 modes reach high (order-unity) amplitudes | 1 | ∼ 0.1 − 1, and are very tightly-wound with almost all the power at short characteristic wavelengths  ∼ 1/.This produces a morphology that resembles a series of narrow concentric rings of alternating high and low density.This is all expected for strong gravito-turbulence in a very thin ( ≪ ) differentially rotating disk (Gammie 2001;Paardekooper 2012;Meru & Bate 2012).Consistent with this, the velocity fluctuations are highly anisotropic: while the in-plane motions driven by these modes have relatively large |  | ∼ |  | ∼ | 1 |   ∼ 0.1   , the vertical motion is an order-of-magnitude smaller, |  | ∼ 0.01   .Again, this is consistent with previous idealized simulations in the regime without efficient magnetic fields and/or stellar feedback to isotropize the turbulence and maintain a thick disk (Hopkins et al. 2012b): "pure gravitoturbulent" driving is unable to maintain appreciable vertical dispersion.
Examining the kinetic stress tensor in Fig. 30 in more detail, it is worth noting that the majority of the  and  dispersion (⟨    ⟩ and ⟨    ⟩) can be directly attributed to the  = 1 modes: subtracting the best-fit  = 1 component from e.g.  (, ) at each  leads to an order-of-magnitude reduction in  kin  .The latter is the true "turbulent" component.But even this is order-of-magnitude larger than the much smaller vertical dispersion, as noted above.Turning to the ⟨ kin   ⟩ component of the Reynolds stress relevant for angular momentum transport, we see that not only is it order-ofmagnitude smaller than the Reynolds stress in our full-physics simulations (Fig. 26), but it often has the opposite sign, i.e.

⟨𝛿𝚷 kin
⟩ = ⟨    ⟩ < 0, driving outflow rather than inflow.Both the large eccentric fluctuations in   , and the sign variations in the Reynolds stress, are directly reflected in the inflow/outflow rates in Fig. 30.
In short, while we see the expected strongly non-linear  = 1 modes, in-plane supersonic turbulence, and efficient fragmentation, as expected for strong gravito-turbulence in a hydrodynamic disk with  cool ≪  dyn , the effective Shakura & Sunyaev (1973)-equivalent  parameter is often smaller than ∼ M 2  , and the Reynolds stresses can even produce net outflow.Taken together, this demonstrates a simple but important point: we cannot simply assume the properties of the disk without magnetic fields would be the same as that with magnetic fields, but with the fields simply "removed" (i.e. that we would have the same Reynolds stress) -the disk is non-linearly different in crucial aspects.
11. CONCLUSIONS Paper I presented the first numerical simulations to follow gas dynamics in a fully-cosmological setup from scales of ≳ Mpc down to < 80 au or < 300 Schwarzschild radii from a SMBH accreting in a quasar phase (Fig. 1).The simulations capture a diverse range of key physical effects relevant to evolution on these scales including multi-band radiation transport; coupling to non-LTE and non-equilibrium atomic/molecular/ionized gas thermo-chemistry and radiative cooling; resolved individual (proto-and main-sequence) star formation and stellar evolution, with associated "feedback" in the form of jets, radiation, stellar mass-loss, and supernovae; and magneto-hydrodynamics with kinetic and non-ideal effects, and amplification only from trace cosmological fields in the inter-galactic medium.In Paper I we studied the large-scale properties of these simulations and how star formation is suppressed close to the BH and a true quasar accretion disk forms on sub-pc scales (Figs. 2,4).In this paper, we present a detailed study of the accretion physics, structure, and origins of the magnetic fields in these simulations.Our key conclusions include: 1.The accretion disk is magnetically-dominated, with plasma  ∼ 10 −6 − 10 −2 even in the midplane.In the inner disk, the field is dominated by the mean toroidal ⟨  ⟩, but with non-negligible mean radial field ⟨  ⟩ and fluctuating toroidal/radial/vertical components | , ,  | (Figs. 5, 10).
2. The magnetic fields arise from simple flux-freezing, viz. the dynamo is closed by advection of radial flux.The field strengths are amplified smoothly from sub-nanoGauss intergalactic fields and typical few-microGauss interstellar fields at ∼ kpc scales, without some sharp change at smaller radii (Fig. 3).The accretion disk initially forms from capture of gas with tangled, quasi-isotropic fields in the galactic nucleus; these are stretched into radial fields as the gas initially falls through the BHROI, and then into toroidal fields as it circularizes (Figs. 11,12), so it is a "flux-fed" and "flux-frozen" disk in a general sense.This leaves characteristic imprints such as sign flips in the toroidal field that are advected with the gas as it accretes (Figs. 6,7,14,15).Given the rapid inflow rates and corresponding advection of flux, we expect and confirm that the toroidal fields do not damp, even over timescales ≳ 10 5 times the dynamical time at our innermost resolved radii (Figs. 14,15).
3. The toroidal field is stronger than the often-quoted limit for linear growth of the "traditional" MRI (  > √    K ; Figs. 12, 16).Instead, we see that the disk lies nominally in the regime of the related Type II/III or SSMI/SHMI magnetic instabilities from e.g.Kim & Ostriker (2000); Pessah & Psaltis (2005); Das et al. (2018), when analyzed according to those previous analytic linear studies.The fluctuating field shows coherent structure with most of the power on wavelengths of roughly the scale height  (Figs. 8, 9).4. The disk is strongly turbulent, with trans-Alfvénic and highly super-sonic, broadly isotropic velocity fluctuations (Figs. 20,22).The velocity and magnetic fluctuations are consistent with one another and the stresses.We see some coherent vertical infall into the disk (Figs. 17,19) but this is negligible compared to inflow through the disk.
5. The disk is coherently eccentric at large radii with a coherent dimensionless eccentricity of | 1 | ≳ 0.1 throughout much of the disk (Fig. 21).This is driven by the infall and non-Keplerian potential at larger radii outside the BHROI.
To first order this dominates the deviations from perfectly circular orbits (Fig. 17) and may produce shocks that drive some of the turbulence especially at the largest radii.But these effects do not dominate the net torque/angular momentum transport of the gas, especially at smaller radii (Figs. 24,26,27;§ 7.2 & 10).
6.The disk is weakly stratified.This is due to a combination of the vigorous highly super-sonic turbulent transport and non-LTE effects on the thermo-chemistry in the more tenuous disk atmosphere and outer disk.As a result, the magnetic-field strength declines only weakly away from the midplane, the temperature is inversely stratified (as expected for a corona/tenuous atmosphere), and the plasma  is weakly stratified but lowest in the disk midplane (Fig. 23).
7. Accretion is driven by a combination of Maxwell and Reynolds stresses.Within the accretion disk, the torques on the gas are dominated by a quasi-turbulent (Fig. 24) MHD torque (as opposed to e.g.torques from gravitational or radiation pressure forces; Fig. 25).Within this, we see that the usual Maxwell+Reynolds stresses are the dominant sources of angular momentum transport (Fig. 26).In the outer disk, where the turbulence is mildly super-Alfvénic, Reynolds stresses dominate, while in the inner disk, the mean-field Maxwell stress dominates, followed by the fluctuating Maxwell and then Reynolds stress, but the three terms are always within an order-of-magnitude of one another (Fig. 27).
8. The disk is not rapidly fragmenting.This is the main subject of Paper I so we refer to that study for details.But we confirm here (e.g.Figs. 4, 5) that the accretion disk is stable against catastrophic turbulent or gravito-turbulent fragmentation on all scales ≪ pc, and that the inflow rates are much larger than the star formation rates.As discussed in more detail in Paper I and Paper IV, on all scales of interest here, the mass loss, injection, and turbulent driving contributions from star formation and/or stellar feedback are completely negligible compared to the other terms we study.
The torques (Fig. 25), Maxwell & Reynolds stresses (Fig. 27), mean radial-flow velocities (Fig. 20), time-steady structure of the disk (Fig. 5), and directly measured inflow/accretion rates (Fig. 5) are all consistent with a sustained gas inflow rate that is remarkably constant in both space and time:  in ∼ 20 − 30 M ⊙ yr −1 is sustained in steady state from scales ≲ 80 au to ∼ 1 pc for the duration of our simulation (≳ 10 5 inner dynamical times, or ∼ 10 4 yr at the highest refinement level).Given the SMBH mass of ≃ 1.3 × 10 7 M ⊙ , this is up to ∼ 100 times the canonical optically thin electron-scattering Eddington mass accretion rate for a nominal radiative efficiency of   = 0.1.Whether this can be sustained to horizon scales (and whether the accretion is radiatively efficient on those scales), and the effects of whatever radiation emerges from the inner disk externally illuminating the outer disk, are important subjects for future study.
10.The accretion disk is likely mis-aligned with the BH spin.The quasar episode here is triggered by tidal interactions with passing, highly-turbulent giant molecular cloud complexes in a clumpy, turbulent, merging galaxy.As such, it is unsurprising that we find the inner accretion disk angular momentum has essentially no correlation with the angular momentum vector of previous generations of BH accretion at earlier cosmic times (which we use as a proxy for the BH spin direction).The disk here is mis-aligned by ∼ 140 • ( § 9), so is both retrograde and mis-aligned.
We further validate these results by comparing an equivalent simulation without any magnetic fields, which we show produces completely different results : it undergoes catastrophic gravito-turbulent fragmentation, with orders-of-magnitude higher star formation rates and ordersof-magnitude lower gas inflow rates; the gas inflow rate drops rapidly towards the center; the disk is razor-thin, with extremely anisotropic turbulence (strong in-plane gravitoturbulent modes but very weak vertical stirring/mixing); the Reynolds stresses are an order-of-magnitude weaker and often have the opposite sign (pushing the disk outwards); the  = 1 modes reach much stronger amplitudes and are tightly wound up to short radial wavelengths of order the pressure scalelength, leading to a concentric-ring-like morphology; and the disk mass and outer extent are reduced by more than an orderof-magnitude (with no real disk outside of ≳ 0.01 pc and the disk mass inside ≲ 0.001 pc that can be supported with  ≫ 1 reduced by orders-of-magnitude).So it is clear that magnetic fields play an absolutely fundamental role on these scales.
In a companion paper (Paper III), we also construct a simple self-similar analytic model for the flux-frozen stronglymagnetized, super-sonically turbulent disks that form consistently in the simulations.We show there that this can, at least qualitatively, reproduce the most important features of the simulations described above, although it is certainly an over-simplification that cannot capture all of the subtleties observed in the simulations (including e.g.behaviors which are clearly not strictly scale-free/power-law-like; see e.g.Fig. 5).In addition to providing some additional consistency checks and aid in interpreting the simulations, these models at least suggest that there is no obvious barrier to extrapolating the behavior seen on our resolved scales down to even smaller scales where dedicated GRMHD simulations are required.
There are many obvious ways in which to extend and im-prove on the simulations here in future work.In principle, one could imagine refining even further, to smaller radii.However, doing so while retaining all of the physics of star formation, molecular cooling, cosmological expansion (let alone the mass of all the gas, stars, and dark matter to ≳ Mpc scales) is both unnecessary and imposes a huge computational overhead.Moreover, at some scales additional physics (e.g.general relativistic effects) will become important.Our goal is therefore not to extend these simulations directly to the ISCO, but to provide new motivation for exploration of disks like that we see here in dedicated GRMHD accretion-disk simulations.Such simulations can reach from scales within the ISCO out to hundreds of gravitational radii -directly overlapping the innermost resolved radii here.So it would be possible to either directly take our inner boundary conditions to set up such smaller-scale simulations, or to use a form like the analytic models in Paper III to set up initial conditions for idealized disk simulations.These simulations can then be used to not only survey different parameter space within the broad category of magneticallydominated disks, but also to make first-principles predictions for the radiation and jets/outflows that should emerge from the accretion disk on scales ≪ 80 au, but which we cannot resolve here.Because this is a first experiment and these kinds of flux-frozen accretion disks have not been explored on such radii, it remains deeply ambiguous whether the disk here should, for example, be radiatively efficient or not, let alone what properties (and orientation) a jet should have.For these reasons and because our goal was to predict the accretion rates and disk outer boundary conditions in the first place, we again caution that these simulations take a simple accretion inner boundary, even though it seems likely that given the combination of strong fields and high accretion rates, significant jets or outflows and radiation must emerge from that boundary.Ultimately, those "feedback" properties could be re-introduced to simulations like those here -running, for example, with some injection/boundary conditions motivated by those smaller-scale simulations -to follow "back up" to larger scales.
Other important extensions of the work here include exploring the accretion disks that form in different galaxies and at different times.While all of our experiments here suggest that there is nothing "special" or "pathological" about the time chosen here for hyper-refinement (relative to any other time that features large inflow rates into the BHROI and therefore potential quasar-level activity), and as noted above the accretion rates and implied quasar luminosities correspond to quasars around the "knee" of the observed luminosity function at these redshifts ( ∼ 4 − 5), it is important to validate this directly.More important still would be to explore how the accretion disks behave in qualitatively different regimes: for example, at lower accretion rates, spanning the vast range between the most luminous sources (like the simulation here) through to "intermediate" luminosity quasars, then Seyferts, low-luminosity AGN, and ultimately, extremely low-accretion-rate systems like those in M87 or our own Galaxy.
Further, there are many more properties to study in the simulations here, which could enable a unique exploration and prediction space.This includes the structure of the obscuring torus, the nature of the broad-line region and narrow-line region transition, the IMF of stars forming in the circumquasar medium and quasar accretion disk, the consequences for transient and gravitational-wave sources from stars and stellar-mass black holes in the disk, predictions for observable signatures of the strongly magnetized accretion disk, and more.We hope to explore these areas and other phenomenology in future work.

Fig
Fig. 1.-Image of the gas surface density in our fiducial cosmological simulation.We show projected gas density on a logarithmic scale (increasing dark-to-bright, dynamic range rescaled in each panel from a median   ∼ 10 19 cm −2 at the largest scale to ∼ 10 8 times larger at the smallest scale).Multiple scales are shown to illustrate the dynamic range of the simulation, which zooms in down to ≈ 80 au scales around a ∼ 10 7  ⊙ SMBH in the center of a massive galaxy at redshift  ≈ 4.4 in a ∼ (100 Mpc) 3 cosmological box.The simulations include explicit multi-band radiation-magnetohydrodynamics, detailed thermochemistry/cooling, self-gravity with resolved individual (proto)star formation, accretion, evolution, and feedback, and many other physical processes ( § 2).Tidally captured gas streams from an encounter with a massive star-forming cloud complex triggered via gravitational torques in a galaxy-scale merger fall into the BH radius of influence (BHROI) at a few ∼ pc and circularize at ∼ 0.1 − 1 pc to form an accretion disk which we follow down to ∼ 300 BH Schwarzschild radii.

Fig. 2 .
Fig.2.-As Fig.1, showing face-on (top) and edge-on (bottom) projections (relative to the nuclear disk) with order-of-magnitude different spatial scales (dynamic range a factor of ∼ 100 in each panel, rescaled as Fig.1).We focus on sub-pc scales; for discussion of the dynamics on larger scales driving these flows see Paper I. The dashed circle at right denotes the inner accretion boundary at  < 80 au.We see the disk circularization radius from the inflowing filament.The disk is thin but has complex structure with spiral arms and multiple warps and some large-scale arms at different angles tracing new inflow with slightly different impact parameter.

Fig. 3 .
Fig. 3.-Magnetic field strengths as a function of BH-centric radius .We measure ⟨ |B| 2 ⟩ 1/2 in spherical annuli from our inner boundary at ∼ 80 au to > Mpc scales, and show the corresponding ⟨ 2  ⟩ 1/2 for the radial, toroidal, and poloidal components of the field (with directions defined relative to the angular momentum vector of the BH accretion disk at  < 0.1 pc).Inside ≪ 10 pc, the resolution is uniformly < 0.01  ⊙ and star formation Fig. 5.-Radial profiles of properties versus spherical BH-centric distance , from < 10 −3 pc to ∼ 1 pc (properties are mass-averaged in concentric shells unless otherwise stated; for discussion of properties on larger (extra)galactic scales, see Paper I). Top left: Inflow rate (summing all gas with   < 0), outflow rate

Fig. 6 .
Fig. 6.-The nuclear disk in a ±0.1 pc box, showing the magnetic field lines (grey) in a face-on slice through the midplane (left; lines show the in-plane   −   field, in gas with | / | < 0.1) or edge-on slice through the disk in cylindrical - coordinates (right; lines show the in-plane   −   field, in gas with | sin  | < 0.1).For each, colors denote the sign of toroidal/azimuthal (  ; top), radial (  ; middle), or poloidal/vertical (  ; bottom) field.Our sign convention is that ẑ and φ point in the direction of the angular momentum vector and direction of co-rotation of the inner gas disk (so   > 0 is prograde), and R points away from the BH.We see large-scale ordered structure tracing the radial inflows from scales outside the BHROI being wrapped around the center as the disk circularizes at ∼ 0.1 pc, with sign changes on large scales determined by the large-scale inflow structure (and reversals or occasional loops forming at the sign change boundaries).There is also an obvious anti-correlation in the sign of   and   .

Fig. 10 .
Fig. 10.-Magnetic field strength and Alfvén speed (defined by  ,  ≡   / √︁ 4  ) relative to circular velocity   ≡ √︁   enc (<  )/ in spherical annuli at BH-centric radii  out to ∼ 1 pc.We compare two different well-separated times (left and right) corresponding to the two times shown in Fig. 8. Top:  ,  /  , showing the (mass-weighted) mean ⟨  ⟩ (line) and ±1 (16 − 84%) range (shaded) in each annulus.Dotted green line shows  = 0 for reference.We clearly see that the toroidal field is dominated by the mean component in the ordered disk region ( ≲ 0.1 pc), with coherent sign flips, while the radial and poloidal fields are dominated by the fluctuating component (that vary in time).The toroidal field also becomes non-negligible compared to gravity at the smallest radii.Middle: | ,  |/  on a logarithmic scale, showing the mean | ⟨   ⟩ | (dotted) and fluctuating rms component    ≡  , 84% −  , 16% ∼ ⟨ |   − ⟨   ⟩ | 2 ⟩ 1/2 (solid; we use a percentile-based definition to avoid biasing the component from a few dense star-forming cells with much-larger |B|).We see only the toroidal field is dominated by the mean component, and even then always has a non-negligible fluctuating component; the fluctuating radial and poloidal fields are comparable to each other but smaller than    by a factor of a few in the inner disk (outside the disk region, |  | ≫ | ⟨ ⟩ | and all components are comparable, indicating a turbulence-dominated system).  has mean comparable to    in the inner regions, while | ⟨   ⟩ | ≪ |    |.Bottom: As middle, but plotting the components   in Gauss.The field strength rises crudely as ∝  −1 .

Fig. 11 .
Fig. 11.-Comparison of the analytic predictions for the amplification of B versus radius for the set of models motivated by simple flux-freezing considerations for the toroidal and radial flux, as described in the text ( § 4.2.1;Eqs.1-4).These successfully reproduce the key behaviors in the simulations.We compare the mass-weighted simulation mean (thick line) value and 90% inclusion range (shaded) of |B| in spherical annuli to the value |B predicted | predicted by the models (labeled), versus radius .We show one early time (top) and one late time (bottom), matching the times in Figs. 8 & 10.These different scalings, which make slightly different assumptions about how the gas is compressed as it accretes, all predict |B| within factors ∼ 2 without a large systematic residual or offset.

Fig. 12 .
Fig. 12.-As in Fig. 11, but comparing alternative previously proposed analytic models for amplification of |B|, which are not based on advection/freezing of toroidal/radial magnetic flux (see § 4.2.2, scalings (2), (3), (5), and (6)).We see that these do not accurately describe the simulations, especially noting the much larger vertical axis range here compared to Fig.11.The model used for most previous discussions of toroidal magnetically dominated disks is the predicted saturation value of the linear (local, unstratified) MRI at Fig.13.-Magnetic field strength |B| versus gas density .Top: All gas cells within  < 100 pc at a fixed time (Eulerian), showing mass-weighted median (solid line) and 90% interval (shaded).We compare an analytic |B| ∝  2/3 scaling motivated by flux-freezing of the mean midplane toroidal/radial fields as they accrete ( § 4.2.1-4.2.2).Middle: Evolution of |B| and  over time for fixed Lagrangian gas elements (points; the line shows running median).We select random cells that are close to being accreted (within  < 100 au) in the final simulation time  1 , and trace each back to an earlier time  0 (∼ 5000 Ω −1 inner dynamical times earlier), to plot |B( 1 ) |/|B( 0 ) | versus ( 1 )/( 0 ).As expected for a mean-field-dominated disk in steady-state with most of the mass in midplane inflow, this traces a similar power-law |B| ∝  2/3 .Bottom: Variation in |B| with  at a given radius (narrow radial annuli at the different  shown).To compare we normalize |B| and  to their mean within each annulus.Here variation in |B| in an annulus is not dominated by global compression of the mean field, so is not expected to obey the same powerlaw, and instead varies weakly (∝  0.2−0.3 ), as expected if fluctuations are generally modest relative to the mean field.The tail to high /⟨( ) ⟩ is dominated by rare sites of star-formation, so this reflects the physics of local collapse (along field lines).

Fig. 16 .
Fig.16.-Instability parameter space map, following e.g.Fig.3ofPessah & Psaltis (2005).We select gas in the warm, volume-filling phases (excluding gas with  < 1000 K) in a radial annulus (factor of ∼ 2 range of  centered on  = (0.001, 0.01, 0.1, 1) pc, corresponding to the different panels, as labeled), within the disk midplane (|  | <  at each ), and plot the volume-weighted histogram (with a linear scale, increasing from 0 in black to 1 in white) of the gas in the space of the toroidal  ,  /  ≡   /(  √︁ 4  ) versus vertical "effective minimum wavenumber"  eff    ,  /  ∼ (2  / ) (  /(  √︁ 4  ) ).Dashed and dotted white lines show the critical wavenumbers  1  and  2  defined inPessah & Psaltis (2005) which separate the different regimes of instability labeled (for an analytic, linear, laminar, unstratified system).The bottom-left quadrant corresponds to the classical (linear, unstratified, weak-field) MRI, top-left to a stable regime, while the top-right and bottom-right to the "Type II" and "Type III" (bouyancy-related) instabilities defined therein (or the suprathermal slow mode instability [SSMI] and suprathermal hybrid mode instability [SHMI], respectively, inDas et al. 2018).We see that the volume-filling midplane phases in the simulations reside well into the Type II/III (SSMI/SHMI) regime, with the largest wavelengths (∼ ) broadly similar to the dividing line between those modes, i.e. the toroidal field is stronger than the maximum amplitude at which the linear MRI grows under the usual laminar, unstratified disk assumptions( max, MRI

Fig. 19 .
Fig.19.-Velocity field line structure on top of a gas density projection, as Fig.9.We see clear radial motion through the midplane of the disk carrying most of the mass (more obviously showing mean velocity in this direction, as compared to the magnetic field line structure in Fig.9), and vertical inflow onto the disk from the tenuous atmosphere.Interior to the disk we see turbulent cells of size ∼ .
ence length ∼  in the midplane and |v  | ∼ |v  | is evident in the edge-on projections in Figs. 17 & 19.
Fig. 21.-Amplitude of the asymmetric modes in the gas density distribution: fitting the (face-on) projected surface gas density to Σ gas () = Σ 0 () [1 + ∞ =1   ( ) cos (  +  0,  [] ) ] and plotting |   |.Here |  1−8 | sums the first 8 modes in quadrature, showing this is dominated by the global  = 1 mode at all the radii shown.The relatively large values of   clearly indicate that the disk has strong spiral and coherent eccentric disk structure that can drive shocks and inflow, as clearly evident in the velocity field lines in Fig. 17.
Fig.22.-TurbulentMach numbers of the gas in radial annuli, at two times as Figs. 10 & 20.We take the velocity dispersion   estimated as 1/2 the 16 − 84% range (to avoid being pulled by outliers) of   , weighted by gas mass, in radial annuli restricting to gas within |  |/ < 0.3 of the midplane, for each component of  in cylindrical coordinates.We then define the sonic Mach number M ,  ≡   /⟨  ⟩ and Alfvén M ,  ≡   /⟨  ⟩.The turbulence is broadly trans-Alfvénic and highly super-sonic (reflecting  ≪ 1 in the disk).In the inner disk, the turbulence becomes mildly sub-Alfvénic (  , ,  ∼ 0.3   , so  total ∼ 0.5   ) -this closely maps to the radii where the mean toroidal field ⟨  ⟩ becomes larger than the fluctuating field in Fig.10, as expected.The radius of this transition can fluctuate by a factor of a few as new magnetic flux advects into and through the disk.
Fig. 23.-Vertical profile of disk properties ( § 6): gas density  (top left), magnetic field strength |B| (top right), temperature  (bottom left), and plasma  (bottom right).For each, we mass or energy-average in cylindrical annuli at radius  (line colors, labeled) in factor ∼ 3 intervals from ∼ 100 au to ∼ 1 pc.The density is stratified as expected for / ∼ 0.1 − 0.3 at this range of radii; magnetic field strengths peak in the midplane but are more weakly stratified.The temperatures are coolest in the midplane and rise in the diffuse, more optically-thin corona, leading to the plasma  also rising (weakly) with scale height.Energy-weighted thermal averages appear more "noisy" owing to physical multi-phase structure (the average can be dominated by e.g.hot gas in shocks or SNe at large radii).In the density plot, we illustrate for one radius (∼ 5000  schw ) a decomposition into a two-component "disk" (Gaussian/exponential or  ∝ sech 2 (/ H  ) here; thin dotted) plus power-law "halo/coronal" ( ∝ ( 2 +  2 ) −  halo with  halo ∼ 3.5 here; thick dotted) profile.In this decomposition, the disk component with / ∼ 0.1 contains most of the mass, but there is clearly an extended vertical gas distribution beyond the extrapolation of the vertical profile from smaller .The vertical thermal structure is expected for a more tenuous corona/atmosphere above the disk, but the weak stratification within the disk of  and |B| owes to a combination of optical depth effects and strong super-sonic turbulence driving efficient vertical mixing ( § 6.3).

Fig. 24
Fig. 24.-Map of the net torque  ≡ r × a (taken directly from the simulation) on each gas element  (inside  < 0.1 pc left, or < 0.01 pc right), in the direction of the local angular momentum vector j ≡ r × v: |  • ĵ|/ 2  (from −1 in blue to +1 in red).The image is a face-on projection of the disk midplane, as in Fig. 6-7.The torques are clearly dominated by a fluctuating component, involving non-axisymmetric/tightly-wound modes with wavenumbers   ∼ a few.This is broadly expected if a combination of Reynolds & Maxwell stresses dominate the torques, potentially originating from the instabilities considered in Fig. 16.
acceleration of the parcel).20We compare this to the in-code specific angular momentum j ≡ r × v, to calculate the change in the scalar specific angular momentum  • ĵ.It is convenient to express this in units of   c Ω =  2 c , with  • ĵ/ 2 c representing the fractional angular momentum loss of a quasicircular orbit per dynamical time.These are shown in a twodimensional projection in the disk in Fig. 24.We clearly see that the torques follow a complicated structure in space, with turbulent/fluctuating and non-axisymmetric modes dominating (with wavenumbers   ∼ a few and   ≳   ,   ≳   ).

7. 2 .
Contributions to the Magnetic and Kinetic Stresses 7.2.1.Which Components of the Kinetic/Magnetic Stress Regulate Angular Momentum Transfer?

c
Fig. 27.-Quantitative properties of the Maxwell & Reynolds stresses from Fig. 26 in more detail.We plot two times as in Figs. 10 & 20 (panels), showing the   component of the Maxwell stress (both the total −⟨     ⟩ and fluctuating −⟨       ⟩ component, as labeled) and the Reynolds stress (⟨     ⟩).Note that the total Maxwell stress, −⟨     /4  ⟩ (i.e.not just the fluctuating component), is what appears in the angular momentum (AM) equation, while ⟨    ⟩ includes additional pure advection terms that do not influence the AM (i.e.only ⟨     ⟩ gives rise to AM exchange).Each quantity is averaged over several snapshots around the time of interest and in annuli within the disk, excluding any dense star-forming clumps.The shaded range shows the ±1 range of | −     /4  |, etc., in all cells in each annulus.Reynolds stresses tend to be larger at radii ≳ 0.01 − 0.1 pc and especially outside the ordered disk, where the turbulence is weakly super-Alfvénic (Fig.22).Within the inner disk, Maxwell stresses dominate by a factor of ∼ 10.The total/mean Maxwell stress is larger than the fluctuating component by factors of several, consistent with the strong mean fields in e.g.Fig.15but distinct from predictions of e.g. the weak-seed-field MRI.We also (top) plot the instantaneous inflow rate /3    Σ gas from Fig.25.This agrees well with the sum of total Maxwell+Reynolds stresses, as expected if they dominate the angular momentum transport.