††thanks: ∗E-mail<EMAIL_ADDRESS>
, ††thanks: $\dagger$NASA Hubble Fellow

# FORGE’d in FIRE II: The Formation of Magnetically-Dominated Quasar Accretion
Disks from Cosmological Initial Conditions

Philip F. Hopkins1,∗ Jonathan Squire2 Kung-Yi Su3 Ulrich P. Steinwandel4 Kyle
Kremer1,† Yanlong Shi1 Michael Y. Grudić5,† Sarah Wellons6 Claude-André
Faucher-Giguère7 Daniel Anglés-Alcázar8,6 Norman Murray9 Eliot Quataert10
1TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA
91125, USA  2Physics Department, University of Otago, 730 Cumberland St.,
Dunedin 9016, New Zealand  3Black Hole Initiative, Harvard University, 20
Garden St, Cambridge, MA 02138, USA  4Center for Computational Astrophysics,
Flatiron Institute, 162 5th Ave., New York, NY 10010 USA  5Carnegie
Observatories, 813 Santa Barbara St, Pasadena, CA 91101, USA  6Department of
Astronomy, Van Vleck Observatory, Wesleyan University, 96 Foss Hill Drive,
Middletown, CT 06459, USA  7CIERA and Department of Physics and Astronomy,
Northwestern University, 1800 Sherman Ave, Evanston, IL 60201, USA
8Department of Physics, University of Connecticut, 196 Auditorium Road,
U-3046, Storrs, CT 06269, USA  9Canadian Institute for Theoretical
Astrophysics, University of Toronto, Toronto, ON M5S 3H8, Canada  10Department
of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

###### Abstract

In a companion paper, we reported the self-consistent formation of quasar
accretion disks with inflow rates $\sim 10\,{\rm M_{\odot}\,yr^{-1}}$ down to
$<300$ Schwarzschild radii from self-consistent 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 $\beta\sim 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 $\sim 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-$\dot{M}$
Shakura-Sunyaev-like disk.

###### keywords:

quasars: general — accretion, accretion disks — quasars: supermassive black
holes — galaxies: active — galaxies: evolution — galaxies: formation

##  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
$\gtrsim 10\,{\rm M_{\odot}\,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., 2007; Aller
& Richstone, 2007; Kormendy et al., 2011) and to dramatically influence galaxy
formation and evolution (Croton et al., 2006; Hopkins et al., 2006, 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 $\alpha$-disk” model: this
takes disks to be geometrically thin (height $H\ll R$), optically thick
(black-body like), sub-sonically turbulent (sonic Mach number
$\mathcal{M}_{s}<1$), slowly cooling ($t_{\rm cool}\gg t_{\rm dyn}=1/\Omega$),
thermal-pressure-dominated (plasma $\beta\equiv c_{s}^{2}/v_{A}^{2}>1$),
radiatively efficient, well-ionized, and parameterized by an effectively
constant-$\alpha$ viscosity, where $\alpha\sim(\delta v)^{2}/c_{s}^{2}<1$
represents some Maxwell or Reynolds stresses and the kinematic viscosity
scales as $\nu\equiv\alpha c_{s}^{2}/\Omega$. 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 $\beta_{\rm|z|\gg H}\lesssim 1$
(Miller & Stone, 2000); magnetically “arrested” disks where magnetic pressure
halts accretion near the innermost stable circular orbit (ISCO) at low
accretion rates (Tchekhovskoy et al., 2011); gravitoturbulent disks relevant
at low values of Toomre $Q$ (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 $\alpha$-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 ($\beta\gg 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
($\sim 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 $\dot{M}$
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 $\dot{M}$, 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,111In this manuscript, we
will often refer to cylindrical radial/azimuthal/vertical coordinates $R$,
$\phi$, $z$, defined with respect to the angular momentum vector of the inner
accretion disk and centered on the SMBH (so $\hat{z}$ points along the angular
momentum vector and $\hat{\phi}$ points in the rotation direction). Given our
focus, we will generally use the terms toroidal and azimuthal interchangeably
(and likewise for poloidal and vertical). We distinguish the spherical
radius/distance $r$ from the cylindrical radius $R$. 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.

Table 1: Some common variables used throughout this manuscript (others are defined throughout where relevant). $R$, $\phi$, $z$ | Cylindrical radial, azimuthal, and vertical coordinates (centered on the SMBH, aligned with the inner disk)
---|---
$r$ |  Spherical (3D) distance $r\equiv|{\bf r}|$ from the SMBH
${\bf B}$, $B_{i}$ |  Magnetic field ${\bf B}$ and components $B_{i}$ (e.g. radial, toroidal, poloidal components $B_{R}$, $B_{\phi}$, $B_{z}$)
${\bf v}$, $v_{i}$ |  Gas velocity ${\bf v}$ and components $v_{i}$ (e.g. radial, azimuthal, vertical components $v_{R}$, $v_{\phi}$, $v_{z}$)
$\delta{\bf B}$, $\delta{\bf v}$ |  Fluctuating part of ${\bf B}$ or ${\bf v}$ (value relative to mean in some annulus/region, e.g. $\delta{\bf B}\equiv{\bf B}-\langle{\bf B}\rangle$)
$\rho$, $n_{\rm gas}$, $\Sigma_{\rm gas}$ |  Gas 3D density $\rho$ or number density $n_{\rm gas}$ (in ${\rm particles\,cm^{-3}}$), and projected surface density $\Sigma_{\rm gas}$
$T$, $c_{s}$ |  Gas temperature $T$ ($T_{\rm rad}$, $T_{\rm dust}$ denote radiation and dust temperatures) and thermal sound speed $c_{s}\equiv\sqrt{k_{B}\,T/\mu\,m_{p}}$
$v_{A}$, $v_{A,\,i}$ |  Alfvén speed $v_{A}\equiv|{\bf B}|/\sqrt{4\pi\rho}$, and component-wise $v_{A,\,i}\equiv B_{i}/\sqrt{4\pi\rho}$
$\beta$, $\mathcal{M}_{s}$, $\mathcal{M}_{A}$ |  Plasma $\beta\equiv c_{s}^{2}/v_{A}^{2}$ parameter, sonic $\mathcal{M}_{s}\equiv|\delta{\bf v}|/c_{s}$ and Alfvén $\mathcal{M}_{A}\equiv|\delta{\bf v}|/v_{A}$ Mach numbers
$H$ |  Gas disk vertical scale-height $H$ (defined within a given annulus $R$)
${\bf j}$, $\boldsymbol{\tau}$ |  Specific angular momentum vector ${\bf j}\equiv{\bf r}\times{\bf v}$ and specific torque vector $\boldsymbol{\tau}\equiv{\bf r}\times{\bf a}$
$V_{c}$, $v_{\rm K}$ |  Total circular velocity $V_{c}^{2}\equiv G\,M_{\rm enclosed}(<r)/r$ (including all mass), Keplerian speed $v_{\rm K}\equiv G\,M_{\rm BH}/r$ (so $V_{c}\rightarrow v_{\rm K}$ as $r\rightarrow 0$)
$M_{\rm BH}$, $\dot{M}_{i}$ |  SMBH mass $M_{\rm BH}$ and inflow/outflow/SF rates $\dot{M}_{i}=$ $\dot{M}_{\rm in}$, $\dot{M}_{\rm out}$, $\dot{M}_{\ast}$, respectively
$\boldsymbol{\Pi}_{ij}^{a}$ |  Stress tensor $\boldsymbol{\Pi}$: tensor component $ij$, $a$ denotes different physical contributions (e.g. magnetic, kinetic, thermal $\boldsymbol{\Pi}^{\rm mag}$, $\boldsymbol{\Pi}^{\rm kin}$, $\boldsymbol{\Pi}^{\rm therm}$)
$Q_{i}$ |  Effective Toomre $Q$ parameter $Q_{i}\equiv\sigma_{i}\,\kappa_{\rm ep}/\pi\,G\,\Sigma_{\rm gas}$ where $\sigma_{i}=c_{s},\,v_{A},\,|\delta v_{R}|,\,(c_{s}^{2}+v_{A}^{2}+\delta v_{R}^{2})^{1/2}$ appear in
|  the thermal $Q_{\rm therm}$, magnetic $Q_{\rm mag}$, turbulent $Q_{\rm
turb}$, and “total” $Q_{\rm eff}$ parameters, respectively (and $\kappa_{\rm
ep}^{2}\equiv(2\Omega/R)\,\partial_{R}(R^{2}\Omega)$)
$\Omega$, $t_{\rm dyn}$, $t_{\rm cool}$ |  Orbital frequency $\Omega\equiv V_{c}/r$, dynamical time $t_{\rm dyn}\equiv\Omega^{-1}$, and gas cooling time $t_{\rm cool}$
$|a_{m}|$ |  Mode amplitude for non-axisymmetric modes with azimuthal integer wavenumber $m$ (e.g. eccentricity $|a_{1}|$)
$R_{g}$, $R_{\rm schw}$, $R_{\rm BHROI}$ |  BH gravitational radius $R_{g}\equiv G\,M_{\rm BH}/c^{2}$, Schwarzschild radius $R_{\rm schw}\equiv 2\,R_{g}$, and BH radius of influence
|  $R_{\rm BHROI}\equiv G\,M_{\rm BH}/\sigma_{\rm gal}^{2}$ (where
$\sigma_{\rm gal}$ represents the “parent” galactic velocity dispersion)

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
self-consistently around a SMBH, from cosmological simulations (using a super-
Lagrangian hyper-refinement technique) down to scales of $\sim 80\,$au (less
than $\sim 300\,R_{\rm schw}$). We observed the self-consistent formation of a
true ($Q\gg 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
$\gtrsim$ Mpc onto the galaxy, through the galactic nucleus, the BH radius of
influence (BHROI), torus-like regions, all the way down to the accretion-disk
scales. We also studied how turbulence, magnetic fields, and radiation-
hydrodynamics produced star formation on large scales and lead, self-
consistently, to the suppression of star formation on small scales ($\lesssim
0.1-1\,$pc) around the SMBH, allowing for the conditions to transition self-
consistently from “galactic-type” or “ISM-like” conditions at $\gtrsim$ pc
scales to “accretion-disk-like” at $\ll$ 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 ($\beta\ll 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 super-sonic 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 “flux-frozen 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 from Paper I, and Table 1 defines
some useful variables we will refer to throughout. In § 3 we review the basic
conditions and properties of the ISM predicted on sub-pc scales, many of which
are detailed in Paper I but helpful to set the context here (including
morphology and basic structure in § 3.1, thermo-chemical and stability
properties in § 3.2, and radial trends in § 3.3). Then, in § 4 we review the
physics of star formation in the disk vicinity, again from Paper I, including
the role of magnetic fields suppressing star formation in § 4.1, the weak
effects of stellar accretion and feedback on inflows in § 4.2, and its
relation to efficient angular momentum transfer in § 4.3. In § 5 we examine
the magnetized accretion disks formed in the simulations in more detail,
discussing both the strength and detailed structure (§ 5.1) and physical
origins (§ 5.2) of the strong magnetic fields. In § 6 we consider the same for
the velocity field structure in the midplane (§ 6.1) and out-of-plane (§ 6.2),
its relation to global coherent eccentric/lopsided disk modes (§ 6.3) and the
details of the turbulent structure (§ 6.4), its physical origins/driving (§
6.5), and the (relatively weak) role of turbulent resistivity (§ 6.6). In § 7
we discuss the vertical structure and profiles of various thermo-chemical and
magnetic disk properties (§ 7.1) and their (weak) stratification (§ 7.2) as
well as the physics behind this (§ 7.3). In § 8, we explore the physical
torques and angular momentum exchange processes in the disk (§ 8.1) and their
relation to different stresses including traditional Maxwell and Reynolds
stresses (§ 8.2). With this in mind, § 9 compares to a variety of previous
literature models of magnetized accretion disks, including magnetically
arrested (§ 9.1), magnetically elevated (§ 9.2), galactic/star-forming (§
9.3), toroidal-field dominated (§ 9.4), and decaying (§ 9.5) disks. We briefly
describe how the disks are likely to be mis-aligned with the pre-existing BH
spin in § 10. In § 11, we contrast a simulation without magnetic fields and
describe how this leads to runaway nuclear star formation (§ 11.1), orders-of-
magnitude lower accretion rates (§ 11.2), and a razor-thin, much smaller and
lower-mass gravitoturbulent $\alpha$ disk (§ 11.3). We summarize and conclude
in § 12.

Figure 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). Multiple scales are shown to illustrate the dynamic range of
the simulation, which zooms in down to $\approx 80\,$au scales around a $\sim
10^{7}\,M_{\odot}$ SMBH in the center of a massive galaxy at redshift
$z\approx 4.4$ in a $\sim(100\,{\rm 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
$\sim$ pc and circularize at $\sim 0.1-1\,$pc to form an accretion disk which
we follow down to $\sim 300$ BH Schwarzschild radii.  Figure 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. 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 $r<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.

##  2 Methods

The simulations studied here are presented and extensively described in Paper
I. Briefly, we begin from a $\sim(100\,{\rm cMpc})^{3}$ cosmological periodic
box at redshift $z\sim 100$ with a primordial trace magnetic field, and follow
it as a cosmological galaxy formation simulation following the Feedback In
Realistic Environments (FIRE) project methods from Hopkins et al. (2018,
2023b) in the code GIZMO222A public version of GIZMO is available at
http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html (Hopkins, 2015), until
a redshift $z\sim 4.4$ when a period of violent activity induces large inflows
into the central $\sim\,$kpc of the galaxy. We then hyper-refine to go to
higher and higher resolution as in Anglés-Alcázar et al. (2021), but self-
consistently switch at sufficiently high resolution to resolve individual
(proto)star formation, accretion and evolution and protostellar disk structure
in the central $\lesssim 10-100\,$pc of the galaxy using the STARFORGE project
methods from Grudić et al. (2021); Guszejnov et al. (2021). These can be
integrated self-consistently for a mix of “individual star” (STARFORGE) and
“stellar population” (FIRE) particles given the methodology in Paper I. We
continue to refine to a target resolution of $\Delta m<0.01\,{\rm M_{\odot}}$
in the central $\lesssim 10\,$pc, to follow gas inflows and disk formation
down to $<300\,$ Schwarzschild radii around the super-massive black hole of
mass $\sim 1.3\times 10^{7}\,{\rm M_{\odot}}$. 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 constrained-gradient 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; Hopkins et al., 2022d);
self-gravity with adaptive self-consistent softenings scaling with the
resolution and high-order Hermite integrators capable of accurately
integrating $\gtrsim 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 self-consistently model the thermochemistry and
opacities in gas with densities from $n\sim 10^{-8}-10^{16}\,{\rm cm^{-3}}$
and temperatures $\sim 1-10^{10}\,$K in a range of radiation and cosmic ray
environments (see Paper I). Fig. 2 illustrates some of the complex phase
structure that emerges even in just the nuclear regions.

As shown in Paper I, some of these physics are not important on the scales we
will study here, although they may play a crucial role in self-consistently
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\,{\rm M_{\odot}}$ (we briefly re-ran with refinement a
factor $\sim 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 $\sim 10-20\,{\rm au}$ and time resolution as small as
$\sim$ days. The most relevant physics on these scales are gravity, (ideal)
MHD, and explicit radiation-thermodynamics.

Figure 3: Magnetic field strengths as a function of BH-centric radius $r$. We
measure $\langle|{\bf B}|^{2}\rangle^{1/2}$ in spherical annuli from our inner
boundary at $\sim 80\,$au to $>$ Mpc scales, and show the corresponding
$\langle B_{i}^{2}\rangle^{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 $r<0.1\,$pc). Inside $\ll 10\,$pc,
the resolution is uniformly $<0.01\,M_{\odot}$ 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 $\sim
1-15\,{\rm\mu G}$ at $\sim 1-10\,$kpc (labeled), but the fields are amplified
and approach kG at $r\sim 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).  Figure 4: Wedge plot showing various disk properties (colormaps
labeled) including: 3D gas density $n_{\rm gas}$, gas temperature $T$, plasma
$\beta\equiv c_{s}^{2}/v_{A}^{2}$, ratio of toroidal to total field strength
$|B_{\rm tor}|/|{\bf B}|$, ratio of azimuthal to total velocity
$|v_{\phi}|/|{\bf v}|$, ratio of magnetic pressure $P_{\rm B}$ to total non-
rotational kinetic/turbulent ram pressure $P_{\rm turb}$, and “effective”
total (including thermal+magnetic+turbulent support) local Toomre $Q_{\rm
eff}$ parameter. The plot shows cylindrical coordinates $R,\,\phi$, but with
the cylindrical radius $R$ stretched on a log scale from $\sim 80\,$au to
$\sim 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 ($|z|/R<0.1$). Densities increase inwards;
temperatures vary but tend to be mostly cold-to-warm at these scales;
$\beta\ll 1$ and decreases in the inner disk $R\lesssim 0.01\,$pc; the
toroidal field and azimuthal (rotational) velocities dominate inside $\lesssim
0.1\,$pc; the magnetic pressure is crudely of order turbulent pressure but
varies locally from $\sim 0.01-10$ times its value; and the effective
stability parameter $Q_{\rm eff}\gg 1$ on all these scales.  Figure 5: Radial
profiles of properties versus spherical BH-centric distance $r$, from
$<10^{-3}$ pc to $\sim 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 $v_{r}<0$), outflow rate ($v_{r}>0)$), and star formation rate (stellar
mass formed/accreted within the last dynamical time, in each shell). We see
consistent inflow at up to $\sim 20\,{\rm M_{\odot}\,yr^{-1}}$. Top right:
Toomre $Q$ parameter for the gas alone, separating thermal support ($c_{s}$),
magnetic ($v_{A}$), turbulent ($v_{\rm turb}$) and total ($\delta v_{r,\,{\rm
eff}}^{2}\equiv c_{s}^{2}+v_{A}^{2}+v_{\rm turb}^{2}$). The system is stable
on all these scales, but the relatively modest thermal-only $Q$ produces some
of the large-scale gravitoturbulence and some (very slow/weak) star formation.
Middle left: Gas density (lines show mean, shaded range $\sim 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
$\sim 90\%$ interval), for kinetic ($\rho\,{\bf v}\,{\bf v}$), magnetic
($|{\bf B}|^{2}\,\mathbb{I}/8\pi-{\bf B}{\bf B}/4\pi$), thermal
($n\,k\,T\,\mathbb{I}$), and radiation ($\int{\rm
d}\nu\,e_{\nu}\,\mathbb{D}_{\nu}$) 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 $H/R$ directly measured (median or rms $z$, as
labeled), and expected ($\approx\sigma/V_{c}$) 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 $r$, 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.

##  3 Basic Properties & Conditions on Sub-pc Scales

###  3.1 Morphology and Physical Structure

Figs. 1-2 illustrate the nuclear disk, which forms self-consistently from
cosmological initial conditions. It is visually obvious that there is a well-
ordered but still eccentric and highly sub-structured rotating gas disk. Fig.
1 illustrates that the galactic nucleus transitions from the highly multi-
phase larger-scale ISM (where star formation and stellar feedback produce a
mix of cold/molecular star-forming regions, warm neutral and ionized phases,
hot bubbles, and super-bubbles of SNe remnants, etc.; see Paper I for more
explicit details of the phase structure) to a more uniformly cold/cool,
mostly-atomic medium. It also shows that while the disk is visually well-
ordered when viewed face-on, it has considerable internal turbulent sub-
structure, with obvious warps and breaks alongside large eddies that carry
material to a couple of scale-heights. The phase structure on scales $\gtrsim$
pc and this transition is studied in detail in Paper I.

Fig. 2 expands on Fig. 1 showing face-on and edge-on projections at a few of
the nuclear scales resolved, illustrating the dynamics and “external” feeding
of the disk, which again are discussed in detail in Paper I. To explain
briefly, because it will be important here, we see that a large star-forming
cloud complex has had a close passage to the $\sim 10^{7}\,{\rm M_{\odot}}$
SMBH in the galaxy nucleus. That passage itself is the result of gravitational
dynamical torques from an ongoing multiple merger in a highly turbulent,
chaotic, clumpy massive high-redshift galaxy with a galaxy-integrated SFR in
excess of $>100\,{\rm M_{\odot}\,yr^{-1}}$. Some of the material from the
cloud is tidally stripped by the SMBH, and falls in initially in a radial
stream (essentially a tidal tail), which falls into the BH radius of influence
(BHROI) interior to which the SMBH dominates the potential (a few pc), then
circularizes at $\sim 0.1-1\,$pc. We can see that, as expected given the large
cloud mass, size, and internal turbulent velocity dispersion, not all the
inflowing gas has an identical impact parameter, so some falls in at slightly
different angles and circularizes at different radii, giving rise to warps and
initial eccentricity in the outer disk.

###  3.2 Additional Disk Properties

Fig. 3 illustrates the scaling of the magnetic fields with BH-centric radius
in the simulations (discussed in detail below), showing the full dynamic range
of the simulation. We highlight (1) the sub-nG intergalactic magnetic fields
from which all other fields are self-consistently amplified; (2) the ISM
magnetic fields on scales $\sim 1-10\,$kpc, which have completely “typical”
values of $\sim 1-15\,{\rm\mu G}$, similar to those observed in the local ISM
of the Milky Way (despite this being a massive, high-redshift galaxy), and
form the ultimate source of the nuclear fields; and (3) how the magnetic field
in the galaxy very crudely scales as $|{\bf B}|\propto r^{-1}$ over a wide
range of scales (discussed below). We also highlight (see Paper I for details)
the range of scales where we reach our maximum target resolution (distances
$\lesssim 10\,$pc uniformly have mass resolution $<0.01\,{\rm M}_{\odot}$) and
exclusively follow star formation with the resolved individual-star STARFORGE
stellar evolution/formation models.

Fig. 4 visualizes several different properties in a projected 2D image over a
dynamic range of a factor $\gtrsim 1000$ in radius, now focusing just on the
innermost regions of our simulation from $10^{-3}-1$ pc (recall the simulation
box is $\sim 100\,$cMpc on a side). These include: 3D gas density $n_{\rm
gas}$, temperature $T$, plasma $\beta\equiv c_{s}^{2}/v_{A}^{2}$, ratio of
toroidal to total field strength $|B_{\rm tor}|/|{\bf B}|$, ratio of azimuthal
to total velocity $|v_{\phi}|/|{\bf v}|$, ratio of magnetic pressure $P_{\rm
B}$ to total non-rotational kinetic/turbulent ram pressure $P_{\rm turb}$, and
“effective” total local Toomre $Q_{\rm eff}$ parameter (including
thermal+magnetic+turbulent support). We discuss these properties in more
detail below. Fig. 5 complements this, showing the 1D (averaged in concentric
radial shells) radial profiles of different properties out to $\sim$ pc
scales, allowing us to see the properties in the innermost radii (of greatest
interest here) in more detail (as compared to e.g. Fig. 3, which shows the
full dynamic range of the simulation but necessarily makes it more difficult
to discern the details on smaller spatial scales).

Again, Paper I focuses on these and other properties on “galactic” scales from
$\sim$ pc to $\sim$ Mpc, but for our purposes, the nuclear region is evolving
so rapidly (dynamical times reaching below a year, as we show below) that we
can treat properties on much larger scales as essentially fixed on the
timescales of interest and over which we evolve the nuclear-disk dynamics
(although large-scale properties are important for determining the
initial/boundary conditions in the first place). As a result, we will focus
almost exclusively on sub-pc scales in this paper, but remind the reader that
the simulations self-consistently extend out to $\sim$ Mpc scales.

###  3.3 Radial Trends and Basic Scalings in the Disk

On sub-pc scales, Figs. 4-5 and our analysis in Paper I show:

  1. 1.

Inflow is systematically larger than outflow and star formation (note the
three can co-exist, as the medium is clearly not spherically symmetric, nor in
strict long-term equilibrium), with an order-of-magnitude constant
$\dot{M}_{\rm in}\sim 20-30\,{\rm M_{\odot}\,yr^{-1}}$ sustained into the
central $\sim 80$ au around the SMBH.

  2. 2.

The gas densities increase towards the center, with a slightly more shallow
profile than a singular isothermal sphere ($\rho\propto r^{-2}$), with
clumpiness evident in the visual projection and in the mass-weighted spherical
profile.

  3. 3.

The temperature is more uniformly dominated by cold and warm atomic media, and
rises weakly with higher densities towards the center (reversing the trend at
larger radii where the denser gas is cooler, as the disk is optically thick to
its own cooling radiation at these radii, see Paper I).

  4. 4.

The radiation and thermal energy densities are in rough equilibrium (with the
radiation, dust, and gas kinetic temperatures coming into increasing blackbody
equilibrium at $r\ll\,$pc, and the dust mostly sublimated at $\lesssim$ pc,
see Paper I), but with both of their energy densities well below the magnetic
energy density in the disk, with $\beta$ ranging from $10^{-6}-10^{-2}$
depending on the local phase (and weighting, as we discuss below).

  5. 5.

The kinetic energy of the gas is uniformly larger than magnetic, but most of
that on scales $\ll\,$pc is the disk rotation, i.e. $|{\bf v}|\sim
v_{\phi}\sim V_{c}(r)$, and the disk is rotation-dominated within $\lesssim
0.1\,$pc.

  6. 6.

The remaining non-rotational or “turbulent” kinetic energy density is more
comparable to magnetic energy densities, with $P_{\rm B}/P_{\rm turb}$ ranging
from $\sim 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, as we will
show more explicitly below.

  7. 7.

The magnetic fields are stronger in the center, with $|{\bf B}|$ increasing
slightly steeper than $r^{-1}$, and the toroidal component of the field
dominating inside the radii where the disk is ordered and rotation-dominated.

  8. 8.

A combination of turbulence and magnetic pressure support the observed disk
scale height (i.e. the disk is in approximate vertical equilibrium) with
$H/R\sim 0.1$ in the central regions (and quasi-spherical structure at $\sim$
pc, outside the rotation-dominated region).

  9. 9.

At all radii $\ll$ pc, the “effective’ Toomre $Q$ parameter including
thermal+magnetic+turbulent support is large (this is defined rigorously in
Paper I), dominated by magnetic+turbulent support as expected given the
arguments above. But even the pure-thermal Toomre $Q_{\rm thermal}\sim 10$
from $r\sim 10^{-2}-10^{-1}$ pc inside the disk and rises rapidly to $Q_{\rm
thermal}\gg 10$ at $r<0.01\,$pc.

  10. 10.

Unsurprisingly given that the SFR is $\ll\dot{M}_{\rm in}$, and shown
explicitly in Paper I, interior to $\ll$ pc, the (surface or 3D) density of
stars and dark matter is very small compared to gas. Thus for the local
dynamics, we can on these radii consider the system to be just that of
gas+SMBH (although the stars at larger radii can exert non-negligible
perturbations on the gas at these smaller radii, as we show below).

Paper III compares these profiles and others to those predicted for a standard
Shakura & Sunyaev (1973) or SS73 $\alpha$-disk with the same accretion rate
$\dot{M}$, to show they differ by orders of magnitude. Briefly, SS73 and other
“weakly-magnetized” models assume $\beta\gtrsim 1$, so the vertical support in
the outer disk comes only from thermal pressure ($H/R\sim c_{s}/V_{c}$), with
an effective viscosity ($\nu\propto\alpha\,c_{s}\,H$) provided by some
Reynolds/Maxwell stresses with $\alpha\sim\delta v_{\rm
turb}^{2}/c_{s}^{2}\sim v_{A}^{2}/c_{s}^{2}\sim 1/\beta\ll 1$ leading to
inflow ($\dot{M}\propto\nu\,\Sigma_{\rm gas}$). So $H/R$, $v_{A}/V_{c}$, and
$\delta v_{\rm turb}/V_{c}$ are predicted by SS73 to be smaller by factors of
$\sim 100-1000$ than the values here, while (for the same $\dot{M}$)
$\Sigma_{\rm gas}$ would be larger in SS73 by factors $\sim 10^{4}-10^{6}$.
Correspondingly, the midplane density ($\rho\propto\Sigma_{\rm gas}/H$) is
larger and effective Toomre $Q$ ($\propto\delta v_{\rm
eff}\,\Omega/G\,\Sigma_{\rm gas}\propto H\,\Omega^{2}/\Sigma_{\rm gas}$) is
smaller by factors $\sim 10^{6}-10^{8}$ in SS73 compared to the simulations.

##  4 Stability and (Lack of) Star Formation

As noted above, the SFR interior to $\ll$ pc is much smaller than inflow
rates, so star formation can be (to first order) neglected when considering
the dynamics of the gas locally, accretion onto the BH, and driving of
turbulence. The physics of how star formation “shuts down” on these scales is
largely the focus of Paper I, so we will not examine it in detail here. But
briefly, we wish to review how this occurs because it is important for
determining the physics of greatest interest here.

###  4.1 Magnetic Support and the Slowdown of Star Formation

Paper I demonstrates that two major factors contribute to the SFR slowing down
relative to the dynamical time on sub-pc scales in these simulations
(specifically at radii between $\sim 0.1-1\,$pc). First, the disk becomes
optically thick to its own cooling radiation, which leads to an inversion of
the density-temperature relation (denser gas becomes warmer, as opposed to
cooler in the more diffuse ISM phases at larger radii), quasi-adiabatic
behavior, and overall warmer temperatures of the “cold” phases (e.g. warm
molecular gas at $\gtrsim 1000\,$K). These effects increase the thermal Toomre
$Q$ parameter and Jeans mass (alongside the shear becoming stronger as well),
suppressing local fragmentation (see Fig. 5). But as shown in Paper I, this
alone is insufficient to strongly suppress SF at sub-pc scales and would, if
we remove magnetic fields, lead to a strongly gravitoturbulent disk that could
still strongly fragment down to more like $\sim 10^{-3}\,$pc.

The second factor is therefore magnetic support in the disk, which plays a
dual role. This gives a “magnetic Toomre parameter” $Q_{\rm mag}\sim
v_{A}\,\kappa/\pi\,G\,\Sigma\gg 1$ at these radii, and (related but not
equivalent) a “magnetic Jeans mass” and magnetic critical mass. These become
larger than the entire enclosed gas mass of the disk $M_{\Phi}>M_{\rm
gas}(<r)\sim 2\pi\,\Sigma\,r^{2}$ (again see Fig. 5 and Paper I). Particularly
because the fields are toroidal, meaning there is no way for the gas to
collapse without being strongly opposed by shear or magnetic pressure, this
strongly suppresses gravitoturbulent fragmentation (Kim & Ostriker, 2001;
Lizano et al., 2010; Riols & Latter, 2016; Forgan et al., 2017) as well as
thermal instabilities (Mishra et al., 2022).

Moreover, for a given surface density $\Sigma$, the stronger magnetic fields
produce orders-of-magnitude more efficient torques hence effective “viscosity”
and inflow rates $\dot{M}_{\rm in}$ (as shown below and in Paper I explicitly
by comparison to a simulation without magnetic fields; see also Oda et al.
2009). This thus increases the ratio of $\dot{M}_{\rm in}/\dot{M}_{\ast}$,
making star formation less significant for the dynamics, even if it were
present at the same level. Equivalently, this lowers the steady-state gas
surface density $\Sigma_{\rm gas}$ needed to sustain some inflow rate
$\dot{M}_{\rm in}$, so not only does the SFR become smaller relative to
$\dot{M}_{\rm in}$, but the Toomre $Q$ parameter is further increased.
Together, as shown in Paper III, this means $Q$ in the outer disk is larger in
the simulations by a factor of $\sim 10^{6}-10^{8}$ compared to what would be
predicted from the SS73 model with the same $\dot{M}$.

###  4.2 Weak Effects of Stars on Small Radii

We stress that there is still some star formation at these radii – indeed, the
physics of star formation on scales $\ll$ pc and resulting properties of said
stars is the primary focus of Hopkins (2023b) (henceforth Paper IV). However,
up to the time we have run our simulations (as much as $\sim 10^{5}$ dynamical
times at the smallest refinement radii), the SFR (defined here by the total
mass going into new protostars and/or accreted onto existing stars and/or
proto-stellar disks in some past number of dynamical times, divided by said
time) represents only a small fraction of the inflow rates through each radius
at radii $\ll\,$pc. Specifically, it is at least an order-of-magnitude smaller
even if we include all mass processed through the protostellar disks and stars
(including that eventually returned to the ISM). It is also worth noting that
given the small dynamical times in the inner disk, the system is evolved
nowhere near long enough for supernovae from the stars forming in said disk to
explode, so stellar feedback effects in the innermost disk are even weaker
than one would expect from “instantaneous yield” type models (see discussion
in Torrey et al., 2017).

###  4.3 Relation to Angular Momentum Transfer

We note all of this for two reasons. First, to justify our neglect of star
formation as a crucial dynamical process, driver of torques, etc., in the
disk. Indeed, if we re-start the simulations from a snapshot and entirely
disable either star formation/accretion or stellar feedback, and run for $\sim
10-20$ dynamical times of the inner disk, we see no difference in the bulk
disk dynamics or inflow rates (as we would expect from the simple arguments
above).

Second, because this leads to a clear division in the dynamics and behaviors
of the disk. At radii $\gtrsim 0.3-0.6\,$pc where star formation is not
strongly suppressed, the consequence is that the stellar mass density
$\Sigma_{\ast}$ is much larger than the gas density $\Sigma_{\rm gas}$.
Indeed, at all radii $\gtrsim 0.5\,$pc, the collisionless (stars and/or dark
matter) components dominate the density, while at all radii $\lesssim
0.5\,$pc, gas dominates (see Paper I). This means not only that the relevant
physics is fundamentally more “ISM-like” (with more multi-phase, actively
fragmenting and star-forming, feedback-dynamics-dominated gas) on larger
scales, but also that the dominant source of the torques and angular momentum
exchange is fundamentally different. On larger scales, since
$\Sigma_{\ast}\gg\Sigma_{\rm gas}$, the dominant angular momentum loss comes
from “gravitational torques” where non-axisymmetry in the collisionless
component produces a leading-order torque by forcing gas into dissipative
shocks (Hopkins & Quataert, 2011; Anglés-Alcázar et al., 2021). But this
ceases to be efficient when $\Sigma_{\ast}(r)\lesssim\Sigma_{\rm gas}(r)$, and
in practice we see a transition to “MHD torques” (Maxwell+Reynolds torques)
dominating precisely at that point, which we will show and discuss in detail
in this paper.

So there is a clear division in the dominant physics and angular-momentum
exchange processes, itself related directly to the magnetic-field strengths,
further motivating our attempt to understand their values below.

Figure 6: The nuclear disk in a $\pm 0.1\,$pc box, showing the magnetic field
lines (grey) in a face-on slice through the midplane (left; lines show the in-
plane $B_{R}-B_{\phi}$ field, in gas with $|z/R|<0.1$) or edge-on slice
through the disk in cylindrical $R$-$z$ coordinates (right; lines show the in-
plane $B_{R}-B_{z}$ field, in gas with $|\sin{\phi}|<0.1$). For each, colors
denote the sign of toroidal/azimuthal ($B_{\phi}$; top), radial ($B_{R}$;
middle), or poloidal/vertical ($B_{z}$; bottom) field. Our sign convention is
that $\hat{z}$ and $\hat{\phi}$ point in the direction of the angular momentum
vector and direction of co-rotation of the inner gas disk (so $B_{\phi}>0$ is
prograde), and $\hat{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 $\sim 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 $B_{\phi}$ and
$B_{R}$.  Figure 7: Field lines as Fig. 6, in face-on projection of the
midplane on three different spatial scales: $\pm 1\,$pc (left), $\pm 0.1\,$pc
(middle), $\pm 0.01\,$pc (right). We can more clearly see the field lines
transition from a somewhat less coherent and radially-biased flow with the
filamentary inflow on the largest $\sim 1\,$pc scales to the more ordered and
more clearly azimuthal field geometry on $\ll 0.01\,$pc scales. We also see
the expected $B_{\phi}-B_{R}$ anti-correlation on all scales, with a somewhat
less correlated/more turbulent $B_{z}$.  Figure 8: Field lines as Fig. 6, in
edge-on projection on smaller $\pm 0.01\,$pc scales, at two different times
(earlier at left, later at right), separated by $\sim 2\times 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 $B_{\phi}$ 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
$B_{\phi}-B_{R}$ anti-correlation but at the later time it is less obvious
(but would become more obvious if we subtract the coherent mean-field
$B_{\phi}$ and plot the fluctuations $\delta B_{\phi}\equiv B_{\phi}-\langle
B_{\phi}(R)\rangle$ and $\delta B_{R}$). Coherent modes with wavenumber
$k_{z}\sim k_{R}\sim 1/H$ are obvious and especially prominent in $B_{z}$ at
the later time. We denote the $z=0$ and $R=0$ axes with the thin white dotted
line for clarity: we see no evidence for a sign switch in $B_{\phi}$ at $z=0$
as would be predicted if the toroidal $B_{\phi}$ were ultimately sourced by a
dominant mean-field $B_{z}$.  Figure 9: Field lines (as Fig. 8; grey)
superimposed on an edge-on gas density projection (colors, $\rho$ as labeled)
in $R$-$z$ coordinates, in a narrow wedge ($\sin{\phi}<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 can see clear radial/vertical field structure and cells/loops
with scale length $\sim H$. Approaching $R\rightarrow 0$, the disk becomes
relatively thick, while being thinner at intermediate radii.  Figure 10:
Magnetic field strength and Alfvén speed (defined by $v_{A,\,i}\equiv
B_{i}/\sqrt{4\pi\,\rho}$) relative to circular velocity
$V_{c}\equiv\sqrt{G\,M_{\rm enc}(<r)/r}$ in spherical annuli at BH-centric
radii $r$ out to $\sim 1\,$pc. We compare two different well-separated times
(left and right) corresponding to the two times shown in Fig. 8. Top:
$v_{A,\,i}/V_{c}$, showing the (mass-weighted) mean $\langle B_{i}\rangle$
(line) and $\pm 1\sigma$ ($16-84\%$) range (shaded) in each annulus. Dotted
green line shows $B=0$ for reference. We can clearly see that the toroidal
field is dominated by the mean component in the ordered disk region
($r\lesssim 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: $|v_{A,\,i}|/V_{c}$ on a logarithmic scale, showing the mean
$|\langle B_{i}\rangle|$ (dotted) and fluctuating rms component $\delta
B_{i}\equiv B_{i,\,84\%}-B_{i,\,16\%}\sim\langle|B_{i}-\langle
B_{i}\rangle|^{2}\rangle^{1/2}$ (solid; we use a percentile-based definition
to avoid biasing the component from a few dense star-forming cells with much-
larger $|{\bf 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 $\delta B_{\phi}$ by a factor of a few in the inner disk (outside
the disk region, $|\delta B|\gg|\langle B\rangle|$ and all components are
comparable, indicating a turbulence-dominated system). $B_{R}$ has mean
comparable to $\delta B_{R}$ in the inner regions, while $|\langle
B_{z}\rangle|\ll|\delta B_{z}|$. Bottom: As middle, but plotting the
components $B_{i}$ in Gauss. The field strength rises crudely as $\propto
r^{-1}$.

Figure 11: Comparison of the analytic predictions for the amplification of
${\bf 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 (§ 5.2.1). 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 $|{\bf B}|$ in spherical annuli to the
value $|{\bf B}_{\rm predicted}|$ predicted by the models (labeled), versus
radius $r$. 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 $|{\bf B}|$ within factors $\sim 2$ without a large systematic
residual or offset.

Figure 12: As in Fig. 11, but comparing alternative previously proposed
analytic models for amplification of $|{\bf B}|$, which are not based on
advection/freezing of toroidal/radial magnetic flux (see § 5.2.2). 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 $v_{A}^{\rm
pred}\sim\sqrt{c_{s}\,v_{\rm K}}$ from 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-$\beta$ models or models with $|{\bf
B}|\propto\Sigma_{\rm gas}$ sometimes invoked in the literature. A model where
$v_{A}$ saturates at $\sim v_{\rm K}$ fares somewhat better but still shows a
systemic deviation with radius and a normalization offset from the
simulations.

Figure 13: Magnetic field strength $|{\bf B}|$ versus gas density $\rho$. Top:
All gas cells within $r<100\,$pc at a fixed time (Eulerian), showing mass-
weighted median (solid line) and $90\%$ interval (shaded). We compare an
analytic $|{\bf B}|\propto\rho^{2/3}$ scaling motivated by flux-freezing of
the mean midplane toroidal/radial fields as they accrete (§ 5.2.1-5.2.2).
Middle: Evolution of $|{\bf B}|$ and $\rho$ over time for fixed Lagrangian gas
elements (points; the line shows running median). We select random cells that
are close to being accreted (within $r<100\,$au) in the final simulation time
$t_{1}$, and trace each back to an earlier time $t_{0}$ ($\sim
5000\,\Omega_{\rm inner}^{-1}$ dynamical times earlier), to plot $|{\bf
B}(t_{1})|/|{\bf B}(t_{0})|$ versus $\rho(t_{1})/\rho(t_{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 $|{\bf B}|\propto\rho^{2/3}$. Bottom:
Variation in $|{\bf B}|$ with $\rho$ at a given radius (narrow radial annuli
at the different $r$ shown). To compare we normalize $|{\bf B}|$ and $\rho$ to
their mean within each annulus. Here variation in $|{\bf B}|$ in an annulus is
not dominated by global compression of the mean field, so is not expected to
obey the same power-law, and instead varies weakly ($\propto\rho^{0.2-0.3}$),
as expected if fluctuations are generally modest relative to the mean field.
The tail to high $\rho/\langle\rho(r)\rangle$ is dominated by rare sites of
star-formation, so this reflects the physics of local collapse (along field
lines).

Figure 14: Top: Gas mass-weighted magnetic field strength $\langle|{\bf
B}|^{2}\rangle^{1/2}$ and mean toroidal field $\langle B_{\phi}\rangle$ in the
innermost region of the simulation (radii $80\,{\rm au}<r<120\,{\rm au}$), as
a function of time. The time is in units of the dynamical time $1/\Omega_{\rm
inner}\equiv\sqrt{r_{\rm inner}^{3}/G\,M_{\rm enc}(<r_{\rm
inner})}\approx\sqrt{r_{\rm inner}^{3}/G\,M_{\rm BH}}$, with $r_{\rm
inner}\equiv 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 $\langle B_{\phi}\rangle>0$, as evident in e.g. Fig. 6,
is accreted, with gas from larger radii with $\langle B_{\phi}\rangle<0$
moving in towards these smaller radii. We see that the magnitude of both
$\langle|{\bf B}|^{2}\rangle^{1/2}$ and $|\langle B_{\phi}\rangle|$ recover
quickly after sign-flip events and remain stable for $\gtrsim 10^{4}$
dynamical times. Bottom: Accretion rate into the central $<80\,{\rm au}$,
averaged in $\Delta t\sim 2\,$yr ($\sim 50\,\Omega_{\rm inner}^{-1}$)
increments.

Figure 15: Mean toroidal $\langle B_{\phi}\rangle$ 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 $B_{R}$ ($\partial\langle B_{\phi}\rangle/\partial
t\approx\partial(v_{\phi}\,B_{R})/\partial R+...$), multiplied by the
dynamical time $\Omega^{-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 $B_{\phi}-B_{R}$ 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”).

##  5 Structure and Origins of the Magnetic Fields

###  5.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 $\sim\,$kG at $\lesssim 10^{-3}\,$pc, somewhat steeper
than $|{\bf B}|\propto r^{-1}$, from $\sim\,$pc scales. We see directly in
Fig. 3 that something like $\langle|{\bf B}|\rangle\propto r^{-1}$ also
describes the field strength (very crudely) at $>$ pc scales out to $\gtrsim$
kpc. Defining the poloidal, toroidal, and radial components of the
field,333Throughout this paper, we consider a fixed angular momentum axis at
any given time defined by the net angular momentum of gas interior to
$r<0.01\,$pc, though our results are not particularly sensitive to exactly
where we define this cutoff radius so long as it is within the visually-
obvious disk. The poloidal and toroidal and radial field components and
azimuthal, vertical, radial velocity components are defined with respect to
this. Our sign convention is such that for a toroidal/azimuthal field, a
positive value indicates prograde fields (aligned with the gas rotation). we
see in Figs. 3, 4, & 5 a clear transition from an isotropic or slightly
radially-biased field at $\gtrsim 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 $R$-$\phi$ coordinates
considering a wedge of some width in $|z|/R$) and edge-on (in cylindrical
$R$-$z$ coordinates considering a wedge of some width in $\phi$), in one
example time snapshot within the ordered disk ($R\lesssim 0.1\,$pc). Fig. 7
considers the face-on field structure at three different scales
$R\lesssim\,(0.01,\,0.1,\,1)$ pc. Fig. 8 shows the edge-on field-line
structure on the smaller scales $\lesssim 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 can also
plainly visually see an anti-correlation in the signs of $B_{\phi}$ and
$B_{R}$, 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
$R$-$z$ projection, and there is some vaguely “jet like” vertical bipolar
field in $R$-$z$. Interestingly, looking at the sign of the toroidal field in
the edge-on projection, we see that there can be sign flips of $B_{\rm tor}$
at different vertical heights, as well as at different radial intervals; but,
these are generally at $|z|\gtrsim H$ – i.e. outside the body of the disk
(with coherence lengths $\gtrsim H$). It is worth noting that there is no sign
flip of $B_{\phi}$ across the midplane, as is often seen if $B_{\phi}$ 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 $B_{R}$
and $B_{z}$ interior to the disk, with wavelength $\sim H$.

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 $R$,
$\phi$, $z$ components, as well as their $\sim 1\sigma$ range.444To be robust
against outliers from e.g. small-scale structure in small sub-volumes
containing nascent protostellar disks, for example, we define the $\sim
1\sigma$ value of the value of each magnetic component as $1/2$ of the
$16-84\%$ mass-weighted inclusion interval. If we simply consider the usual
RMS $\langle B_{i}^{2}\rangle^{1/2}$ we see similar trends, but slightly
systematically larger values, indicating that the “tails” of the distribution
are somewhat fatter than Gaussian, as well as slightly more noise owing to
some sub-structure. The exact quantitative values also vary somewhat if we
mass or volume or magnetic-energy weight the results, but this does not change
any of our relative comparisons or conclusions. We plot these both in absolute
units, as well as in units of the Alfvén speed $v_{A,\,i}\equiv
B_{i}/\sqrt{4\pi\,\rho}$, relative to the circular velocity $V_{\rm
c}\equiv\sqrt{G\,M_{\rm enc}(<r)/r}$ 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
$B_{\phi}$). The mean $\langle B_{\phi}\rangle$ 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 $|\langle B_{\phi}\rangle|\gtrsim|\delta B_{\phi}|\gtrsim|\delta
B_{R}|\sim|\delta B_{z}|\gtrsim|\langle B_{R}\rangle|\gtrsim|\langle
B_{z}\rangle|$, i.e. the mean toroidal and vertical fields are the strongest
and weakest, respectively, with the fluctuating toroidal second-strongest,
followed by the fluctuating radial and vertical fields. We can 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.3), the magnetic pressure dominates the vertical
support of the disk but with $\mathcal{O}(1)$ contribution from trans-Alfvénic
turbulence: this is reflected both in direct comparison of the turbulent
velocity components to the various $v_{A,\,i}$ below; or from the kinetic
energy densities, or relative contribution to the effective stability
parameter $Q$, or comparison of the disk $H/R$ to $v_{A}/V_{c}$ in Fig. 5. In
other words, $\rho\,v_{\rm turb}^{2}\sim\langle|{\bf B}|^{2}\rangle$ broadly
speaking, as we quantify in more detail below. Of course, both $v_{A}$ and
$v_{\rm turb}$ are much larger in the simulations than in an SS73 disk which
assumes $\beta\gg 1$.

Figure 16: Instability parameter space map, following e.g. Fig. 3 of Pessah &
Psaltis (2005). We select gas in the warm, volume-filling phases (excluding
gas with $T<1000\,$K) in a radial annulus (factor of $\sim 2$ range of $R$
centered on $R=(0.001,\,0.01,\,0.1,\,1)$ pc, corresponding to the different
panels, as labeled), within the disk midplane ($|z|<H$ at each $R$), 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
$v_{A,\,\phi}/c_{s}\equiv B_{\phi}/(c_{s}\,\sqrt{4\pi\,\rho})$ versus vertical
“effective minimum wavenumber” $k^{\rm
eff}_{z}\,r\,v_{A,\,z}/c_{s}\sim(2\pi\,r/H)\,(B_{z}/(c_{s}\,\sqrt{4\pi\,\rho}))$.
Dashed and dotted white lines show the critical wavenumbers $k_{z}^{c1}$ and
$k_{z}^{c2}$ 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 ($\sim H$) 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 ($v_{A,\,\phi}^{\rm max,\,MRI}\sim\sqrt{c_{s}\,v_{\rm 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.

###  5.2 Physical Origins of the Field Strength & Structure

####  5.2.1 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
$\langle B_{\phi}\rangle$ 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., 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
$V_{\rm c}$), it is tidally stretched into a radial stream of length $\ell$
and width $w$ (as we see occurring in Figs. 6-7). For $B_{z}$ and $B_{\phi}$,
the perpendicular areas – $(\Delta r)\,(r\,\Delta\phi)\sim\ell\,w$ and
$(\Delta r)\,(\Delta z)\sim\ell\,w$, respectively – are increased by the
stretching555For pure radial infall in a Keplerian potential, the
perpendicular extent $w$ is tidally compressed with a transverse acceleration,
$\Delta a_{\bot}=w\,\partial a_{\rm grav,\,\bot}/\partial w=-\Omega^{2}\,w$,
while the radial extent $\ell$ is stretched with a parallel acceleration,
$\Delta a_{\|}=\ell\,\partial a_{\rm
grav,\,\|}/\partial\ell=2\,\Omega^{2}\,\ell$. If one begins from small
coherent shear/expansion velocities, this means the perpendicular areas for
$B_{z}$ and $B_{\phi}$ ($\sim\ell\,w$) will increase while the perpendicular
area for $B_{R}$ ($\sim w^{2}$) decreases. so these components are weakened;
in contrast, $B_{R}$ is amplified, owing to the perpendicular compression
(perpendicular area $(r\,\Delta\phi)\,(\Delta z)\sim w^{2}$). Equivalently we
can think of the field lines as being “stretched” in the radial direction as
they are dragged. The expected amplification of $B_{R}$ ranges between
$r^{-0}$ (for e.g. an infalling clump with $w\ll r$, or no radial motion) to
$r^{-2}$ (maximal case where $w\sim r$ 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 $b$ and circularizes at $\sim
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 $\langle
B_{\phi}\rangle$ 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
$|\delta B_{i}|\gtrsim|\langle B_{i}\rangle|$ for all components of ${\bf B}$)
at radii outside of the BHROI, there are sign flips in the “initial” field
which is now stretched into some mean $\langle B_{\phi}\rangle$ 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 $B_{\phi}$ 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, $\langle B_{\phi}\rangle$ is constantly replenished by the
steady supply of radial magnetic flux into and through the disk.

For flux-freezing, the compression in the $z$ direction suggests $B\sim
B_{\phi}\propto 1/H$ (Machida et al., 2006). Meanwhile direct analysis of the
simulations or simple analytic considerations give weak compression in the $R$
direction.666This is expected if the solution is wind-like (constant $v_{r}$),
free-fall-like, spherical ($H\sim r$, akin to the Solar wind), or simply
circularizing while conserving specific angular momentum from a broad/flat
initial distribution (for some uniform initial bulk velocity dispersion of the
captured cloud, $dr\propto db$). As shown below, a yet simpler isotropic-flux-
freezing expectation $B\propto\rho^{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 $|{\bf
B}|$, if we fit $B\propto\rho^{2/3}$ to the simulations using the midplane
values of $\langle|{\bf B}|^{2}\rangle^{1/2}$ and $\langle\rho_{\rm
midplane}\rangle$ (Figs. 3 & 5), we obtain $|{\bf B}|\sim(2-8)\,{\rm\mu
G}\,(n/{\rm cm^{-3}})^{2/3}$ (depending on how we weight it; or $|{\bf
B}|\sim(2-6)\,{\rm\mu G}\,(r/10\,{\rm kpc})^{-1}$) – i.e. at Galactic radii
this extrapolates to typical mundane values of $|{\bf 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 $\alpha\sim 0.1$ disk with the same
$\dot{M}$ (even though $\beta\ll 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 $\langle B_{\phi}\rangle$ at $z=0$ (i.e. reversals of $\langle
B_{\phi}\rangle$ as one vertically crosses the midplane), as well as the broad
anti-correlation between $B_{\phi}$ and $B_{R}$. The relative amplification
versus suppression above also explains why we see a dominant mean $\langle
B_{\phi}\rangle$ component with much weaker mean $\langle B_{R}\rangle$ and
$\langle B_{z}\rangle$ components (Fig. 10; note also that $\langle
B_{z}\rangle$ seems to grow somewhat more slowly than $\langle B_{R}\rangle$).
Instead, the radial and poloidal fields are more dominated by their
turbulent/fluctuating components, whose amplitudes are of order $|\delta
B_{R,\,z}|\sim|\delta(v_{R,\,z}/v_{A})\,\langle B_{\phi}\rangle|$ – 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 $\gtrsim 0.01$ (so
ideal MHD remains a good approximation). This owes to a combination of high
temperatures ($\gtrsim 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 ($\gtrsim 100\,{\rm M_{\odot}\,yr^{-1}}$
in the central $\sim 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 $\lesssim 10^{-17}-10^{-15}$ in much
of the disk.777Calculating the generalized Elsasser number $N_{\rm E}\equiv
B^{2}/(\rho\,\eta\,\Omega)\sim v_{A}^{2}\,t_{\rm orbit}/\eta$, where $\eta$ is
the largest of e.g. ambipolar, Hall, or Ohmic resistivities, and plugging in
numbers for typical values in the simulations at $\lesssim$ pc scales gives
$N_{\rm E}\gtrsim 10^{13}$ through most of the inner disk, while $N_{\rm
E}\lesssim 1$ is required for resistivity to have a large global effect (or
for e.g. ambipolar diffusion to generate substantial “drift” between ions and
neutrals on the timescales of interest). That would generally require an ion
fraction $x_{i}\sim n_{i}/n_{\rm neutral}\ll 10^{-15}$ in the accretion disk.
In § 6.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.

####  5.2.2 Validating that the Mean Field is Indeed Driven by Flux-Freezing
and Advected Flux

We now consider various quantitative tests of the picture described above in §
5.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. 1.

Flux-freezing of the radial/azimuthal magnetic flux, the scenario from §
5.2.1. Here $\langle B_{\phi}\rangle$ 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 § 5.2.1.

  2. 2.

Flux-freezing of a dominant mean dipolar/poloidal field. One could instead
assume that the magnetic field was dominated by a mean
poloidal/vertical/dipolar field $\langle B_{z}\rangle$, which was flux-frozen
as mass advected inwards in a laminar disk (e.g. torqued by magnetic braking).
Then the fixed flux-to-mass ratio would predict $|{\bf B}^{\rm
pred}|\sim|\langle B_{z}\rangle|\propto\Sigma_{\rm gas}$ in the disk.

  3. 3.

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 $v_{A}^{\rm
pred}\approx\sqrt{c_{s}\,v_{\rm K}}$, or $B^{\rm
pred}\sim\sqrt{4\pi\,\rho\,c_{s}\,v_{\rm K}}$ (the value derived for linear
MRI growth in a laminar, unstratified, local analysis in Pessah & Psaltis
2005).888Technically, following the more general dispersion relation in Pessah
& Psaltis (2005), this should be modified to $v_{A}^{\rm
pred}\approx\sqrt{c_{s}\,v_{\rm K}/|1+\partial\ln{B}_{\phi}/\partial\ln{R}|}$
in the presence of radial magnetic stratification (see discussion in Begelman
& Armitage, 2023). This correction generally makes a small (tens of percent)
difference, but even in the annulus where the correction is maximized (near
$\sim 0.002$ pc where $\partial\ln{B}_{\phi}/\partial\ln{R}$ is close to $-1$)
the correction never exceeds a factor of $\sim 2.5$, so does not change any of
our conclusions.

  4. 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 $\sim 3-10$),
with isotropic, tangled fields, and fluctuating components much larger than
mean ($|\delta{\bf B}^{\rm pred}|\gg|{\bf B}|$).

  5. 5.

The Small-Scale Sub-Sonic, Gravitational/Protostellar Dynamo. Various models
for the local dynamo in more slowly cooling dense gas in molecular
clouds/clumps/cores collapsing to proto-stellar disks have argued for
saturation at fixed $\beta\sim 1$, giving $v_{A}^{\rm pred}\sim c_{s}$ (Mocz
et al., 2017), again with isotropically tangled fields ($|\delta{\bf B}^{\rm
pred}|\gg|{\bf B}|$).

  6. 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 $v_{A}^{\rm pred}\sim V_{c}$

First, of course, 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 §
5.2.1. Obviously, 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 $|{\bf B}|(r)$ 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 $|{\bf B}|\sim|\langle B_{\phi}\rangle|$
motivated by the simple flux-freezing considerations of model (i) as described
in § 5.2.1. We see that these reasonably reproduce the absolute magnitude and
radial trend of $|{\bf B}|$ with relatively little scatter. If we compare
$B\propto r^{-1}$ or $B\propto H^{-1}$ (both noted in § 5.2.1), the $B\propto
H^{-1}$ scaling clearly exhibits even smaller scatter and more accurate
prediction of the mean $|{\bf B}|$ compared to $B\propto r^{-1}$. We obtain
almost as good a fit with the simpler expression $B\propto\langle\rho_{\rm
mid}\rangle^{2/3}$ (shown explicitly in Fig. 13).999Note this does more
accurately predict $|{\bf B}|$ than assuming $|{\bf B}|\propto\rho^{2/3}$ for
all $\rho$ within a given disk annulus, as some of the dense, cold and/or hot,
diffuse gas phases actually have more similar $|{\bf B}|$ if they are at the
same radial annulus near the midplane, owing to the fact that their $|{\bf
B}|$ is dominated by the mean $\langle B_{\phi}\rangle$ component and
collapse/expansion can occur along these field lines. So locally on small
scales within the disk isotropic flux-freezing is not always a good
approximation, even if it is not a bad approximation for the global behavior
of the mass in the disk on large scales. We can also assume
$\langle|B_{\phi}|\rangle\propto H_{A}^{-1}$ where $H_{A}\approx\langle
v_{A,\,\phi}\rangle/\Omega\approx\langle|B_{\phi}|/(\Omega\,\sqrt{4\pi\,\rho})$
is the scale height set by the magnetic pressure (specifically focusing on the
toroidal component providing the vertical support). This gives
$|B_{\phi}|\propto\rho^{1/4}\,r^{-3/4}$, which also provides a good fit.

To compare, Fig. 12 plots the predicted value of $|{\bf B}^{\rm 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 $|{\bf 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 $|{\bf B}|$),
and predicting the incorrect normalization of $|{\bf B}|$. Model (vi) fares
somewhat better but is still notably offset from the actual field strengths
(and, of course, 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 (§ 6.4) that the typical Alfvén Mach numbers in the disk at
$\lesssim 0.1\,$pc are $\mathcal{M}_{A}\sim 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 $R\gg$ 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 $|{\bf 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 of course expected if the disk is in
steady state with steady inwards accretion, and $\langle\rho\rangle$ or $H$
are monotonic functions of $R$, 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 $B_{\phi}$ 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 ($\sim\dot{M}_{\rm in}/M_{\rm gas}(<r)$), 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 obviously see the
sign reversals in $\langle B_{\phi}\rangle$ 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 $\langle B_{\phi}\rangle$ 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 $\langle
B_{\phi}\rangle$, $\partial_{t}\langle
B_{\phi}\rangle=\langle\partial_{z}(B_{z}v_{\phi}-B_{\phi}v_{z})\rangle+\langle\partial_{R}(B_{R}v_{\phi}-B_{\phi}v_{R})\rangle$,
to include only the term that accounts for the stretching of the radial field
by the toroidal flow $\partial_{t}\langle
B_{\phi}\rangle\approx\partial_{R}\langle v_{\phi}\,B_{R}\rangle$. Note the
pure radial $B_{\phi}$ transport term $\partial_{R}\langle
B_{\phi}\,v_{R}\rangle$ is usually smaller, but not always negligible, while
the vertical flux divergence $B_{z}\,v_{\phi}$ vanishes and the vertical
inflow term $B_{\phi}\,v_{z}$ is small. In any case we see this approximation
works well at describing quantitatively the sign flips and trends in $\langle
B_{\phi}\rangle$. This behavior is generally distinct from the predictions of
the alternative models (ii)-(vi) above. Moreover, if we multiply
$\partial_{R}\langle B_{\phi}\,v_{R}\rangle$ by the characteristic dynamical
time $t_{\rm dyn}=\Omega^{-1}$, this appears to provide a remarkably good
order-of-magnitude estimate of the saturation $\langle B_{\phi}\rangle$ – this
is expected if either (i) the accretion is dynamical (so the amplification
time is limited to some multiple of $\Omega^{-1}$) or (ii) if the disk is
trans-Alfvénically turbulent (so $H\sim v_{\rm turb}/\Omega$ by definition in
vertical equilibrium, and the turbulent magnetic dissipation time is
$\sim\Omega^{-1}$) or (iii) if the midplane toroidal flux is lost via buoyancy
on the vertical buoyancy timescale (a few to tens of times $\Omega^{-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 radial+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 § 5.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)).

Figure 17: Velocity streamlines in a cylindrical wedge inside $<0.1\,$pc, as
in Fig. 6, face-on in the midplane (left, showing $v_{R}$-$v_{\phi}$) and
edge-on (right, showing $v_{R}$-$v_{z}$), colored by the sign of the
(spherical) $v_{r}\equiv{\bf v}\cdot\hat{\bf 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
($m=1$) “slow mode” which propagates inwards from large $r$ 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.  Figure 18: Velocity streamlines (as Fig.
17), but now with the lines color-coded 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 large-scale 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 $v_{R}$ can be large, at a given $R$), but these
clearly reflect large-scale, geometrically separate structures rather than
turbulent local motions.

####  5.2.3 Discriminating Between Models for the “Turbulent” Field

Now we consider the sub-dominant, but still not negligible, fluctuating or
“turbulent” field components $\delta{\bf 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 $\delta B_{i}$ in
Fig. 10 are consistent with the usual trans/sub-Alfvénic relation $|\delta
B|/|{\bf B}|\sim|\delta v|/v_{A}$, as expected. The origins of the
turbulent/fluctuating magnetic-field components are therefore directly 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 $|\langle B_{\phi}\rangle|$ and $|{\bf B}|$) is an order-of-
magnitude or more larger than the characteristic strength
$v_{A}\sim\sqrt{c_{s}\,v_{\rm 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 ($v_{A}\gtrsim\sqrt{c_{s}\,v_{\rm 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 analytic linear stability analysis of a laminar, unstratified disk
– can be, to leading order, represented in a two-dimensional plot of
dimensionless vertical wavenumber $\tilde{k}_{z}\equiv
k_{z}\,r\,v_{A,\,z}/c_{s}$ (with $v_{A,\,i}\equiv|B_{i}|/\sqrt{4\pi\,\rho}$)
versus a measure of the dimensionless toroidal magnetic field strength
$v_{A,\,\phi}/c_{s}$. The parameter space is then divided into four
characteristic regimes based on the intersection of two critical wavenumbers:
$\tilde{k}_{z}^{c1}\equiv
k_{z}^{c1}\,r\,v_{A,\,z}/c_{s}\approx(V_{c}/c_{s})\,\sqrt{-2\,{\rm
d}\ln{\Omega}/{\rm d}\ln{r}}\approx\sqrt{3}\,v_{\rm K}/c_{s}$, and
$\tilde{k}_{z}^{c2}\equiv
k_{z}^{c2}\,r\,v_{A,\,z}/c_{s}\approx(V_{c}/c_{s})\,\left[v_{A,\,\phi}^{4}/c_{s}^{2}\,V_{c}^{2}-\kappa^{2}/\Omega^{2}\right]^{1/2}\approx\sqrt{(v_{A,\,\phi}/c_{s})^{4}-(v_{\rm
K}/c_{s})^{2}}$. Modes with $\tilde{k}_{z}>{\rm
MAX}(\tilde{k}_{z}^{c1}\,,\,\tilde{k}_{z}^{c2})$ are stable; those with
$\tilde{k}_{z}^{c2}<\tilde{k}_{z}<\tilde{k}_{z}^{c1}$ are unstable to the
traditional MRI (with the near-vertical part of $\tilde{k}_{z}^{c2}$ defining
the upper limit $v_{A,\,\phi}<\sqrt{c_{s}\,v_{\rm K}}$ above which, for the
conditions considered in Pessah & Psaltis 2005, the MRI growth rate
vanishes);101010If we retain the radial stratification term
$\hat{B}_{\phi}\equiv\partial\ln{B_{\phi}}/\partial\ln{R}$, then
$\tilde{k}_{z}^{c2}$ is modified to
$\sqrt{(1+\hat{B}_{\phi})^{2}\,(v_{A,\,\phi}/c_{s})^{4}-(v_{\rm
K}/c_{s})^{2}}$, so the boundary for MRI-like behavior shifts to
$v_{A,\,\phi}<\sqrt{c_{s}\,v_{\rm K}/|1+\hat{B}_{\phi}|}$. If we include this
$|1+\hat{B}_{\phi}|^{-1/2}$ correction in Fig. 16, the effect is small here:
the $\tilde{k}_{z}^{c2}$ boundary between MRI and “Type III” shifts upwards by
a factor $\sim 1.1-1.4$ in each panel, which has no effect on our conclusions.
Although $\hat{B}_{\phi}$ is broadly similar to $-1$ on average in Fig. 3,
evaluating this correction term in factor $\sim 2$ radial intervals from $\sim
80$ au to $\sim 10\,$pc, we find that $|1+\hat{B}_{\phi}|^{-1/2}$ never
exceeds a factor of $\sim 3$. those with
$\tilde{k}_{z}^{c1}<\tilde{k}_{z}<\tilde{k}_{z}^{c2}$ are unstable to a second
or “Type II” instability; and those with $\tilde{k}_{z}<{\rm
MIN}(\tilde{k}_{z}^{c1}\,,\,\tilde{k}_{z}^{c2})$ are unstable to a third or
“Type III” instability.111111Technically, 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. Note 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 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
$k_{z}\sim 2\pi/H$. We see that the simulations lie solidly in the “Type
II/III range,” even if we focus only on the warmer (higher-$\beta$) volume-
filling phases of the gas in the disk where it is multi-phase (in the colder,
denser gas, $\beta$ is even smaller and the simulations lie even further from
the traditional MRI regime). This confirms our intuition and quantitative
statement that $v_{A,\,\phi}\gg\sqrt{c_{s}\,v_{\rm 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 § 6.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 $\beta\ll 1$, non-isothermal ${\rm d}T/{\rm
d}|z|>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, $\sim{\rm a\
few}\,\Omega^{-1}$) at the characteristic scales $\sim H$ 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 $|\delta B_{z}|\sim|\delta
B_{R}|$ (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 $R\phi$ 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 $\rho$ and $B_{\phi}$) 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 $\sim(|v_{A,\,z}|/v_{\rm K})^{1/3}\,\Omega$ at wavenumber
$k\sim 1/H\sim\Omega/v_{A,\,\phi}$). 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.121212In contrast, the Type II/SSMI instabilities operate on much
shorter wavelengths, do not depend on differential rotation, and their linear
eigenmodes feature weak Maxwell and Reynolds stresses with anti-
aligned/opposing angular momentum transport, so in this sense they are more
akin to local convective instabilities.

Briefly, the fact that the simulation modes with wavelength $\sim H$ reside
order-of-magnitude around the Type II/III (or SHMI/SSMI) dividing line
($\tilde{k}_{z}^{c1}$) is not actually surprising: when
$(2\pi\,r/H)\,(v_{A,\,z}/c_{s})\sim\tilde{k}_{z}^{c1}\sim\sqrt{3}\,v_{\rm
K}/c_{s}$ in a disk dominated by toroidal magnetic pressure with $H\sim
v_{A,\,\phi}/\Omega$, the Type II/III dividing line is equivalent to
$v_{A,\,\phi}\sim 4\,v_{A,\,z}$. So residing broadly “near” this line is
simply a statement that the typical $|B_{\phi}|$ is not orders-of-magnitude
larger than the typical $|B_{z}|$ in a cell. Note that we use the value of
$|B_{\phi}|$ and $|B_{z}|$ in each cell for this histogram (to define
$v_{A,\,\phi}$ and $v_{A,\,z}$) – if we instead replaced these values with the
mean field $|\langle B_{\phi}\rangle|$, $|\langle B_{z}\rangle|$ (closer to
the assumption in Pessah & Psaltis 2005) then $v_{A,\phi}$ and the horizontal
position of the simulation changes very little (since this is mean-field
dominated per Fig. 10), but $v_{A,\,z}$ is reduced by a factor of $\sim 3-10$.
This places modes with wavelength $\sim H$ 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 $|k_{z}|>1/H$, 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 (${\bf k}\perp{\bf B}$) Parker-like magnetic convective/Rayleigh-
Taylor/interchange modes, as these require $d\ln{(|B|/\rho)}/dz<0$ (for $z>0$
and $\beta\ll 1$; Tayler 1973; Terquem & Papaloizou 1996; Kim et al. 2002),
i.e. that the magnetic scale-height is smaller than the density scale-height,
which we plainly show in § 6.2 is not satisfied. And although the long-
wavelength Parker instability requires only a vertically-decreasing $|{\bf
B}|$ ($d\ln{B^{2}}/dz<0$), the characteristic wavelengths
$\lambda\gtrsim\lambda_{\rm crit}\sim(5-10)\times H$ (Parker, 1966; Kim et
al., 1997; Lee & Hong, 2007) are extremely large ($>R$, given the relatively
large $H$ 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
$\partial\ln{B_{\phi}}/\partial\ln{R}$ 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+\partial\ln{B_{\phi}}/\partial\ln{R}|<\beta\,(v_{\rm
K}/v_{A,\,\phi})^{2}\sim\beta\,(H/R)^{-2}\ll 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 $v_{A}\lesssim\sqrt{c_{s}\,v_{\rm K}}$ (as assumed in e.g. Begelman &
Pringle 2007).

####  5.2.4 Comparison to Field Strengths in “Traditional” $\alpha$ 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” ($\beta\ll
1$) $\alpha$-disk, then building up a sufficient Maxwell stress to produce the
canonical $\alpha\sim v_{A}^{2}/c_{s}^{2}\sim 0.1$ in such a disk would
actually require magnetic fields whose absolute strength (in Gauss) is
approximately a factor of $\sim 10$ larger than those seen in our simulations
(e.g. from SS73 Eq. 2.19 therein, $B_{\rm SS73}\sim 5000\,{\rm
G}\,(R/0.001\,{\rm pc})^{-1.3}$ for $\dot{M}\sim 20-30\,{\rm
M_{\odot}\,yr^{-1}}$ and $M_{\rm BH}\sim 10^{7}\,{\rm M_{\odot}}$).
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 $H/R\sim
c_{s}/V_{c}$, with midplane densities $\rho$ a factor of $\sim 10^{6}-10^{8}$
larger than those seen here (§ 3.3). This in turn would require some process
like the MRI to very efficiently amplify $|{\bf B}|$ up to quite large
absolute values to produce even a modest Alfvén speed $v_{A}\sim|{\bf
B}|/\sqrt{4\pi\,\rho}$, as needed to produce a Maxwell stress that is any
appreciable fraction of $\rho\,V_{c}^{2}$.

This of course has important consequences for the observational properties of
disks like those simulated here: even though $v_{A}/c_{s}$ and $v_{A}/V_{c}$
are much larger here than in a traditional SS73-like disk, $|{\bf B}|$ can
actually be significantly smaller. It also re-emphasizes that the absolute
field strengths in our simulations are not particularly extreme or implausible
(§ 5.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.

Figure 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 $\sim H$.

##  6 Structure of the Velocity Fields

###  6.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 $r$, 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. This shows, as we saw above in the morphology and ${\bf 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 $V_{\rm
c}(r)$), showing both the mean and fluctuating velocities.

As discussed in Paper I, at radii $\gg$ pc, outside the BHROI, the velocity
structure in the nucleus is primarily turbulent with quasi-isotropic velocity
fields – we see this already at $\gtrsim 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 $\delta v_{i}=\langle|{\bf v}_{i}-\langle{\bf
v}_{i}\rangle|^{2}\rangle^{1/2}\sim V_{\rm c}\gg\langle{\bf v}_{i}\rangle$. At
$\lesssim 0.1\,$pc, we see the disk clearly in kinematic space, with $\langle
v_{\phi}\rangle\approx V_{\rm c}$ much larger than other components, which
have a dispersion $\sim 0.1\,v_{\phi}$. The kinematics of the disk are much
more coherent compared to ${\bf 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).

Figure 20: Mean velocity with its dispersion (top), and separation of the mean
$\langle{\bf v}_{i}\rangle$ and fluctuating $\delta{\bf v}_{i}$ velocity
components in both units of $V_{c}$ (middle) and ${\rm 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 ($\lesssim 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 $\delta{\bf v}_{i}\sim
0.1\,V_{c}$. Almost always we see $\langle v_{R}\rangle<0$, indicating inflow,
as expected, with modest amplitude. For $v_{z}$ we show $\langle v_{z}\cdot
z/|z|\rangle$, 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.

###  6.2 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 un-resolved
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 $-H\lesssim z\lesssim H$).

###  6.3 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 $m=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
$\Omega_{p}\sim\Omega(r_{\rm max})$ where $r_{\rm max}\sim$ pc, so at smaller
radii $\Omega_{p}\ll\Omega(r\ll r_{\rm max})$). The amplitude of the
eccentricity declines very weakly as $r\rightarrow 0$. Fig. 21 plots the
amplitude of the eccentric ($m=1$) and higher-$m$ 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” $m=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. $m=2$ bars or higher $m>2$, are
damped as $r\rightarrow 0$, but for slow $m=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 $b_{\rm complex}$. The
cloud complex therefore drives a characteristic (lagging) pattern frequency
$\Omega_{\rm eff}\sim|\delta v|/b_{\rm complex}\sim\sigma_{\rm galaxy}/b_{\rm
complex}\sim 200\,{\rm km\,s^{-1}}/{\rm a\ few\ pc}\sim 10^{-12}\,{\rm
s^{-1}}$, i.e. precession on a $\sim 10^{4}\,$yr timescale.

This means that in an instantaneous sense, for any given gas parcel, its
radial (inflow/outflow) velocity is dominated by where it is instantaneously
in its eccentric orbit. This is clear in Fig. 17. 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 leading-order description of the disk is not a
perturbed circular disk, but a perturbed eccentric disk.

###  6.4 Turbulent Structure

From Fig. 20 we see the turbulence, while sub-dominant to rotation, is still
vigorous, with $\langle\delta v^{2}\rangle^{1/2}\sim 0.1\,V_{\rm c}\sim
v_{A}$. 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 super-sonic). Fig. 22 shows this more
explicitly, plotting the typical sonic and Alfvénic Mach numbers in different
radial annuli.

Figs. 20 & 22 also show that the typical turbulent velocity $\delta v_{\rm
turb}$ is generally larger than the mean velocities in $z$ or $R$ directions,
and not strongly anisotropic ($\delta v_{\phi}\sim\delta v_{R}\sim\delta
v_{z}$, to within a factor of $\sim 2-3$ or so, though note that the apparent
transient dominance of e.g. $\delta v_{z}$ at small $r$ 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 $-H<z<H$, 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 turbulence here is much stronger than in a typical SS73-like
$\alpha$-disk, where (by assumption) the sonic Mach number
$\mathcal{M}_{s}\sim\delta v/c_{s}\sim\alpha^{1/2}<1$ (i.e. the turbulence is
always sub-sonic).

Figure 21: Amplitude of the asymmetric modes in the gas density distribution:
fitting the (face-on) projected surface gas density to $\Sigma_{\rm
gas}(R)=\Sigma_{0}(R)\,[1+\sum_{m=1}^{\infty}\,a_{m}(r)\,\cos{(m\,\phi+\phi_{0,\,m}[R])}]$
and plotting $|a_{m}|$. Here $|a_{1-8}|$) sums the first 8 modes in
quadrature, showing this is dominated by the global $m=1$ mode at all the
radii shown. The relatively large values of $a_{m}$ 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.

###  6.5 What Drives the Turbulence (Or sets its Amplitude)?

####  6.5.1 Does Anything Actually Need to Drive The Turbulence Interior to
the Disk?

Given this vigorous super-sonic turbulence, it is natural to ask what the
“driving mechanism” may be. However, it is not necessarily obvious that a
driving mechanism 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-advected131313In the configuration we consider, wave
packets propagating radially inwards faster than the bulk accretion flow will
have $|\delta{\bf B}|$ decay relative to $|{\bf B}|$, but those moving along
the mean field or otherwise transverse, or with the accretion flow, will be
amplified. Alfvén wave packets is $3/2$ (i.e. we expect $|\delta{\bf
B}|\propto\rho^{3/4}$ moving with a Lagrangian parcel, or equivalently ${\rm
d}_{t}\langle\delta{\bf B}^{2}\rangle\sim-(3/2)\,\langle\delta{\bf
B}^{2}\rangle\,(\nabla\cdot\langle{\bf v}\rangle)$). Similarly, if we instead
consider purely hydrodynamic turbulence, many authors have shown the Euler
equations in globally compressed turbulence (if dissipation is slow) give rise
to a Lagrangian quasi-adiabatic behavior with $|\delta{\bf
v}|\propto\rho^{1/3}$ (i.e. ${\rm d}_{t}\langle\rho\,\delta{\bf
v}^{2}\rangle\sim-(5/3)\,\langle\rho\,\delta{\bf
v}^{2}\rangle\,(\nabla\cdot\langle{\bf v}\rangle)$; see e.g. Vazquez-Semadeni
& Gazol, 1995; Robertson & Goldreich, 2012). So given that we saw
approximately $|{\bf B}|\propto\rho^{2/3}$ for the mean-field from flux-
freezing above (${\rm d}_{t}\langle{\bf B}\rangle^{2}\sim-(4/3)\,\langle{\bf
B}\rangle^{2}\,(\nabla\cdot\langle{\bf v}\rangle)$), we would have
$|\delta{\bf B}|/|{\bf B}|\propto|\delta{\bf v}|/v_{A}\propto\rho^{1/12-1/6}$
increasing with a (weak) positive power of $\rho$ as gas accretes inwards and
$\rho$ increases. This could sustain $\mathcal{M}_{A}\sim 1$ in principle,
even with some damping/dissipation.

The challenge here is the turbulent dissipation time: if the turbulence
dissipates on a timescale $\sim t_{\rm turb,\,diss}$, then the expressions
above should be modified to ${\rm d}_{t}\langle\rho\,\delta{\bf
v}^{2}\rangle\sim-C\,\langle\rho\,\delta{\bf
v}^{2}\rangle\,(\nabla\cdot\langle{\bf v}\rangle)-\langle\rho\,\delta{\bf
v}^{2}\rangle/t_{\rm turb,\,diss}$ (where $C$ can be in the range from $3/2$
to $5/3$ depending on the regime, as described above). Noting that
$|\nabla\cdot\langle{\bf v}\rangle|\sim r/|\langle v_{r}\rangle|\sim(v_{\rm
K}/|\langle v_{r}\rangle|)\,\Omega$ is the inverse of the accretion timescale,
this becomes ${\rm d}_{t}\ln{\langle\rho\,\delta{\bf v}^{2}\rangle}\sim
C\,(v_{\rm K}/|\langle v_{r}\rangle|)\,\Omega-t_{\rm turb,\,diss}^{-1}$, 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
$t_{\rm cool}\ll t_{\rm dyn}$, the expected dissipation time is of order a
couple of turnover times at the driving scale $t_{\rm turb,\,diss}\sim t_{\rm
turnover}\sim L_{\rm drive}/v_{\rm drive}\sim H/|\delta{\bf v}_{\rm
turb}|\sim\Omega^{-1}\sim t_{\rm 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 $\sim M_{\rm gas}/\dot{M}$ is $\sim(v_{\rm
K}/|\langle v_{r}\rangle|)\,\Omega^{-1}\sim(v_{t}/v_{\rm K})^{-2}\,t_{\rm
dyn}$. While the accretion here is very fast compared to more traditional
$\alpha$-disk models (as we discuss below), inserting typical values from e.g.
Fig. 20 into the above still gives an inflow time $\sim 30-100\,t_{\rm dyn}$
or $\sim 5-15\,t_{\rm orbit}$. So we would naively expect the gas to sit too
long at a given $r$, 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 $\delta{\bf
v}\equiv{\bf v}-\langle{\bf v}\rangle$. 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, recent idealized studies
(Skalidis et al., 2021, 2022; Beattie et al., 2020, 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 ${\bf B}$ with $|{\bf B}\cdot\delta{\bf B}|\sim\rho\,\delta
v^{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.

####  6.5.2 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 $\Sigma_{\rm
gas}\,v_{t}^{2}\,\Omega\sim\dot{M}\,\Omega^{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 $v_{t}/V_{\rm c}\sim|a_{m=1}|\sim 0.1$, which is not far from
what we see. Given that the mode is a slow mode (pattern speed small compared
to $V_{\rm 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 ($t_{\rm decay}\sim 1/\Omega$) 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-$m$ ($k\gtrsim
1/H$) modes in Fig. 21, while the global mode is driven by an external
perturbation (in this case, the self-sustaining asymmetric distribution at
large-$r$). This would naturally give rise to turbulence of the form
$\Sigma_{\rm gas}\,v_{t}^{2}\,\Omega\sim\dot{M}\,\Omega^{2}$ (Gammie, 2001).

Both of these are likely playing some role, especially at larger radii ($\gg
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 § 11 below to simulations without magnetic fields. However, at
smaller radii in the disk ($\lesssim 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 ${\bf B}$ and neither therefore
provides an obvious reason the turbulence should saturate at trans-Alfvénic
values in Fig. 22. Moreover we show in Paper I and 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 clear that the torques in the
inner disk are dominated by less coherent (non-global) structures. Moreover,
such a mechanism would not efficiently drive out-of-plane turbulence (leading
to a large anisotropy, unlike what we see; see e.g. Hopkins et al. 2012b and
discussion below). 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 $Q_{\rm 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 § 5.2.3, the “traditional” weak-field MRI – at least
in the simplest linear form usually invoked – does not appear viable as the
primary turbulent driving mechanism. Specifically, as shown above (there and §
5.2.2), the toroidal $B_{\phi}$ 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 § 5.2.3, while this particular form
of the MRI is stabilized under such conditions at sufficiently high
$B_{\phi}$, this just means that different linear instabilities appear.
Specifically, even under those conditions, at $B_{\phi}$ above the threshold
where the MRI growth rate vanishes, in a $\beta\ll 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 ($\sim(v_{A}/v_{\rm K})\,\Omega$ for the largest-
wavelength unstable modes with $k\sim 1/H$), 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 $k\sim 1/H$,
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
$\partial\ln{B}_{\phi}/\partial\ln{R}$ is very close to $-1$ (Begelman &
Armitage, 2023). Of course, it is also plausible that other, distinct