FORGE'd in FIRE: Resolving the End of Star Formation and Structure of AGN Accretion Disks from Cosmological Initial Conditions

It has recently become possible to zoom-in from cosmological to sub-pc scales in galaxy simulations to follow accretion onto supermassive black holes (SMBHs). However, at some point the approximations used on ISM scales (e.g. optically-thin cooling and stellar-population-integrated star formation [SF] and feedback [FB]) break down. We therefore present the first cosmological radiation-magnetohydrodynamic (RMHD) simulation which self-consistently combines the FIRE physics (relevant on galactic/ISM scales where SF/FB are ensemble-averaged) and STARFORGE physics (relevant on small scales where we track individual (proto)stellar formation and evolution), together with explicit RMHD (including non-ideal MHD and multi-band M1-RHD) which self-consistently treats both optically-thick and thin regimes. This allows us to span scales from ~100 Mpc down to<100 au (~300 Schwarzschild radii) around a SMBH at a time where it accretes as a bright quasar, in a single simulation. We show that accretion rates up to $\sim 10-100\,{\rm M_{\odot}\,yr^{-1}}$ can be sustained into the accretion disk at $\ll 10^{3}\,R_{\rm schw}$, with gravitational torques between stars and gas dominating on sub-kpc scales until star formation is shut down on sub-pc scales by a combination of optical depth to cooling and strong magnetic fields. There is an intermediate-scale, flux-frozen disk which is gravitoturbulent and stabilized by magnetic pressure sustaining strong turbulence and inflow with persistent spiral modes. In this paper we focus on how gas gets into the small-scale disk, and how star formation is efficiently suppressed.

1. INTRODUCTION The origins and growth of super-massive black holes (SMBHs) represents one of the most important open problems in extragalactic astrophysics.Most sufficiently-massive galaxies host SMBHs whose masses correlate with various host galaxy bulge properties and reach masses as large as ∼ 10 10  ⊙ (Magorrian et al. 1998;Ferrarese & Merritt 2000;Gebhardt et al. 2000;Hopkins et al. 2007a,b;Aller & Richstone 2007;Kormendy et al. 2011; for a review see Kormendy & Ho 2013).Many constraints indicate that most of this BH mass is assembled via accretion of gas in a few bright quasar phases (Soltan 1982;Salucci et al. 1999;Yu & Tremaine 2002;Hopkins et al. 2006c), giving rise to a picture of "co-evolution" between galaxies and active galactic nuclei (AGN) or quasars (Merloni & Heinz 2008).Understanding this "co-evolution" has crucial consequences far beyond the BHs themselves, for example in the form of AGN "feedback" launching galactic winds (Silk & Rees 1998;King 2003;Di Matteo et al. 2005;Murray et al. 2005;Hopkins et al. 2005a,b;Debuhr et al. 2010;Faucher-Giguère & Quataert 2012;Torrey et al. 2020), regulating galaxy masses (Croton et al. 2006;Hopkins et al. 2006aHopkins et al. ,b, 2008a)), and changing the structure of the circum-galactic or

Cosmic Rays
Dynamically-evolved with LEBRON approximation, coupled to chemistry, sourced from fast shocks from SNe and stellar mass-loss.SSP Particles FIRE: Formation in self-gravitating, Jeans-unstable gas, sampling IMF, when cell resolution > 1  ⊙ .SSP Feedback FIRE: Main-sequence IMF-sampled tracks: radiation, stellar mass-loss (O/B & AGB), supernovae (Types I & II), cosmic rays.Star Particles STARFORGE: Formation in self-gravitating, isolated, resolved Larson cores inside own Hill sphere, accrete bound mass, resolution < 1  ⊙ .Star Feedback STARFORGE: Protostellar & main-sequence single-star tracks with accretion, radiation, jets, surface mass-loss, end-of-life explosions.Supermassive BH Live sink particles formed dynamically.Refinement centers on ∼ 1.3 × 10 7  ⊙ BH, accretion at < 300  schw .Kawakatu & Wada 2008;Nayakshin & King 2007)1 However, the last couple of decades have seen considerable progress on this front at scales ≳ 0.1 − 1 pc.Initial analytic arguments (Shlosman et al. 1989), followed by semi-idealized numerical simulations of different "levels" of the scale hierarchy (Escala et al. 2004;Escala 2007;Mayer et al. 2007;Wise et al. 2008;Levine et al. 2008;Hopkins & Quataert 2010b;Costa et al. 2022), and then simulations using "super-Lagrangian" or "hyper-refinement" techniques to probe small scales (Curtis & Sĳacki 2015;Prieto & Escala 2016;Prieto et al. 2017;Bourne & Sĳacki 2017;Su et al. 2019Su et al. , 2020Su et al. , 2021;;Franchini et al. 2022;Talbot et al. 2022;Sivasankaran et al. 2022) within larger boxes eventually reaching up to cosmological scales (Beckmann et al. 2019;Bourne et al. 2019;Bourne & Sĳacki 2021;Anglés-Alcázar et al. 2021) have led to a robust emergent picture wherein on large scales, gravitational torques between non-axisymmetric structures (including e.g.mergers, bars, large clumps, lopsided/warped disks; see Hopkins & Quataert 2010b), especially between collisionless and collisional components of the galaxy (e.g.torques from stars on gas not only driving angular momentum exchange but inducing shocks which can orders-of-magnitude enhance inflow rates over classical single-component disk models; see Hopkins & Quataert 2011b) can produce inflows of gas on timescales of order the dynamical time, ensuring some can reach sub-pc radii without turning into stars (references above and e.g.Levine et al. 2008Levine et al. , 2010;;Hopkins & Quataert 2011a;Hopkins et al. 2012aHopkins et al. , 2016)).While these represent an enormous progress, there are still many open questions and key issues unresolved by these simulations.In particular, it has not yet been possible to "bridge the gap" between these (≳ pc) scales and the traditional ( ≫ 1) accretion disk.This is not just a question of dynamic range, but of physics: the physics believed to drive accretion on small scales -physics like the magneto-rotational instability (MRI) -are qualitatively different from the physics of gravitational torques on larger scales.And it is by no means clear what physically occurs when the different physics most relevant on different scales intersect.On large scales ≳ 10 − 100 pc, simulations of high-redshift quasar fueling require cosmological dynamics following optically-thin cooling from dusty ionized, atomic, and molecular gas, self-gravity, and detailed models of star formation and stellar feedback which model the formation and collective effects of entire stellar populations (spanning the range of the entire stellar initial mass function [IMF]), including their radiation, acceleration of cosmic rays, massloss, and supernovae.On smaller scales ≲ 10 pc, simulations of star formation need to follow individual stars and protostars as they form, accrete, and grow, while injecting feedback in the form of jets, winds, radiation, and (eventually) supernovae, all in a dusty medium which spans both optically-thin and optically-thick cooling.At even smaller scales around a SMBH (≲ 10 4   , where   =   BH / 2 ) traditional "accretion disk" simulations must be able to accurately evolve radiation-magneto-hydrodynamics, with global simulations that can accurately follow the growth of the magneto-rotational instability even in warped or irregular disks, radiation-pressure dominated fluids with explicit radiation-dynamics (accounting for finite-speed-of-light effects), with opacities dominated by partially-ionized (largely dust-free) gas, and gravity integrators which must follow huge numbers of orbits accurately.
As a result, there have not been simulations that can span all three of these regimes simultaneously and self-consistently.Even today, very few codes include all of the physics listed for even just one of the three scale regimes described above, let alone two or all three.So simulations using super-Lagrangian hyper-refinement have generally either (a) had to stop at some radius or resolution where the physics prescriptions simply cease to make sense (e.g. at ∼ pc scales, for simulations with traditional "galaxy-scale" cooling, star formation and feedback prescriptions as in e.g.Anglés-Alcázar et al. 2021); or (b) consider only restricted special cases like accretion onto lowredshift SMBHs at extremely low accretion rates (≲ 10 −4 times the Eddington limit) in gas-poor ellipticals in "hot halos" (where star formation and many other physical processes above can be neglected relatively "safely"; as in Guo et al. 2022), or (c) simply neglect most of the physics above even on scales where it could be important.
In this paper, we present the first simulation to span all three of these regimes including all of the physics above.The key to this is to leverage a suite of physics that has been developed and extensively studied in the GIZMO code (Hopkins 2015(Hopkins , 2017a) ) over the last several years.On large scales, all of the physics above (and more) has been developed into a physics suite as part of the Feedback In Realistic Environments (FIRE) (Hopkins et al. 2014(Hopkins et al. , 2018b(Hopkins et al. , 2022a) ) project, designed fundamentally for simulations of galaxies on scales where stars can be treated as ensemble stellar populations, so star formation occurs in environments which should fragment to form stars and produces "stellar population particles" which represent many stars that can act on the environment in the form of radiation, cosmic rays, mass-loss, and supernovae.In parallel, we have also developed a suite of physics as part of the STARFORGE project (Grudić et al. 2021;Guszejnov et al. 2021), designed for simulations which resolve individual star formation, where sink particles form representing individual (proto)stars which then follow individual (proto)stellar evolution tracks as they grow, accrete, evolve, ending up on the main sequence, and eventually ending their main sequence lives as remnants or SNe, explicitly modeling jets, radiation, mass-loss, and supernovae.As a part of this, we have developed gravity and radiation-magnetohydrodynamics solvers which have been applied to high accuracy to evolving e.g. the MRI in global disk simulations (Gaburov et al. 2012;Hopkins & Raives 2016;Hopkins 2016Hopkins , 2017b;;Deng et al. 2019Deng et al. , 2020b;;Deng & Ogilvie 2022), dynamics of strongly radiationpressure dominated fluids (Hopkins & Grudić 2019;Hopkins et al. 2020aHopkins et al. , 2021a;;Williamson et al. 2020aWilliamson et al. , 2022;;Lupi et al. 2022;Braspenning et al. 2022;Soliman & Hopkins 2022) including radiation-pressure-dominated AGN accretion disks, and accurately evolving individual gravitational orbits (allowing for "hard" N-body dynamics) for up to millions of orbits (Grudić & Hopkins 2020;Grudić et al. 2021;Guszejnov et al. 2022a;Hopkins et al. 2022c).Crucially, the physics of all of the above are built in a modular fashion in the code, allowing for cross-compatibility -this allows us to evolve all of the relevant physics simultaneously for the first time.
These physics and numerical methods allow us to "zoom in" from truly cosmological initial conditions down to < 300  s around a super-massive BH during an extremely highaccretion-rate quasar episode, and to see the formation of the true accretion disk and cessation of star formation on sufficiently small scales in a self-consistent manner.In § 2, we summarize the numerical methods and physics included ( § 2.1) including both the FIRE ( § 2.2) and STARFORGE ( § 2.3) regimes, and initial conditions ( § 2.4) and architecture ( § 2.5) of the fiducial simulation studied here.In § 3 we study the results of the simulation (including some variants with different physics): we describe the qualitatively different behaviors over the vast hierarchy of scales ( § 3.1), our effective resolution ( § 3.2), the (gas/stellar/dark matter) mass density and accretion rate profiles ( § 3.3), the plasma and thermodynamic properties on these scales ( § 3.4), dynamics of fragmentation and star formation and its suppression at small radii ( § 3.5), and the torques driving inflow ( § 3.6).In § 4 we contrast a simulation that ignores magnetic fields entirely, and in § 5 we summarize the scales where different physics "ingredients" play a crucial role.In § 6 we compare to previous work in different regimes from galactic ( § 6.1) through accretion disk ( § 6.4) scales.We summarize our conclusions in § 7.

Overview and Common Physics
Fundamentally, the simulation suite presented here combines two well-tested numerical physics implementations: the Feedback In Realistic Environments (FIRE) physics (specifically the FIRE-3 version from Hopkins et al. 2022a), and STARFORGE physics (Grudić et al. 2021).Both of these physics modules have been extensively tested in the literature2 so we will only summarize what is included and refer to the relevant methods papers for each, in order to focus on what is novel here (how the two are integrated within our refinement scheme).An even more succinct high-level overview is provided in Table 1.
All of the relevant physics are implemented in the code GIZMO3 (Hopkins 2015).
The simulations evolve the radiation-magneto-hydrodynamics (RMHD) equations, using the meshless finite-mass MFM scheme (a mesh-free Lagrangian Godunov method).The ideal MHD equations are numerically integrated as described in Hopkins & Raives (2016); Hopkins (2016) using the constrained-gradient method from Hopkins (2016) for greater accuracy4, with the addition of non-ideal terms including fully-anisotropic Spitzer-Braginskii conduction and viscosity (implemented as in Hopkins 2017b; Su et al. 2017;Hopkins et al. 2020b, including all appropriate terms needed to self-consistently apply them at arbitrary values of temperature or plasma  and in both saturated and unsaturated limits), as well as ambipolar diffusion, the Hall effect, and Ohmic resistivity (Hopkins 2017b).The RHD equations are integrated using the M1 moments method (as tested and implemented in Lupi et al. 2018Lupi et al. , 2021Lupi et al. , 2022;;Hopkins et al. 2020a;Hopkins & Grudić 2019;Williamson et al. 2020aWilliamson et al. , 2022;;Bonnerot et al. 2021) for each of five bands (H ionizing, FUV/photo-electric, NUV, optical-NIR, and an adaptive-wavelength blackbody FIR band).5As described in Grudić et al. (2021), this includes the ability to evolve the effective wavelength or temperature of the IR radiation field so as to accurately handle wavelength/temperature-dependent opacities and emission from wavelengths of ∼ 0.1−1000 m.Also as described therein, we separately evolve the gas, dust, and radiation temperatures and different radiation bands, with the appropriate physical coupling/exchange terms between these, so that the code can self-consistently handle both limits where the various temperatures are arbitrarily de-coupled from one another and limits where they become closely-coupled (e.g.Bonnerot et al. 2021;Grudić et al. 2021).Note that compared to previous STARFORGE or FIRE RHD simulations, we greatly increase the reduced speed of light, with most runs here using c = 0.1 , though we have tested runs for a limited time with c = 0.01  and = 1  (i.e.no reduced speed of light at all) to validate that the radiation properties in the galaxy nucleus at ≪ 100 pc scales (the regime of greatest We see structure on all scales, with a chaotic, cold, disordered morphology on most scales until an ordered disk forms from capture of gas from a passage of a giant molecular cloud complex (itself triggered by an ongoing galaxy merger in the rapidly-accreting proto-galaxy), forming the accretion disk at ≲ 0.1 pc.interest) are converged.Gravity is solved with an adaptive Lagrangian force softening matching hydrodynamic and force resolution for gas cells, and fixed softenings specified below for the collisionless particles, using a fifth-order Hermite integrator designed to accurately integrate "hard" gravitational encounters (e.g.close binaries) for the entire duration of our simulation (Grudić et al. 2021).We explicitly follow the enrichment, chemistry, and dynamics of 11 abundances (H, He, Z, C, N, O, Ne, Mg, Si, S, Ca, Fe; Colbrook et al. 2017), allowing for micro-physical and turbulent/Reynolds diffusion (Escala et al. 2018), as well as a set of tracer species.
Thermo-chemistry is integrated with a detailed solver de-scribed in Grudić et al. (2021); Hopkins et al. (2022a): we follow all of the expected dominant processes at temperatures of ∼ 1 − 10 10 K including explicit non-equilibrium ionization/atomic/molecular chemistry as well as molecular, fine-structure, photo-electric, dust, ionization, cosmic ray, and other heating/cooling processes (including the effects of the meta-galactic background from Faucher-Giguère 2020, with self-shielding).Crucially, the explicit radiationhydrodynamics is coupled directly to the thermo-chemistry: radiative heating, ionization, and related processes scale directly from the local (evolved) radiation field, and cooling radiation is not simply "lost" (as assumed in many imple- mentations of optically-thin cooling), but is emitted back into the evolved RHD bands appropriately (see Grudić et al. 2021, for various tests demonstrating that this accurately captures the transition between optically thin and thick cooling regimes).We assume a dust-to-gas ratio which scales as  dg = 0.01 (/ ⊙ ) exp (− dust /1500 K), i.e. a standard dustto-metals ratio at low dust temperatures with dust destruction above a critical dust temperature.This allows us to capture the most important dust transition at small radii in these simulations, namely dust destruction within the QSO sublimation radius.We stress that the thermo-chemistry modules are designed to self-consistently include essentially all processes which dominate radiative cooling and opacities from densi-ties  ≪ 10 −10 cm −3 through  ≫ 10 15 cm −3 in proto-stellar disks (with or without dust).We separately account for the dust and gas opacities in each of the ionizing, photo-electric, NUV, optical-NIR, and gray-body IR bands, calculated as an appropriate function of the (distinct) dust and gas temperature and radiation temperature in each band including dust opacities from Semenov et al. (2003), bound-free/ionizing, freefree, Rayleigh and Thompson opacities for free  − , HI, HII, H − , H 2 , CO, and partially-ionized heavy elements evolved, with the abundances of each of these species calculated in the chemical network (see e.  2, but in stars.The chaotic merger morphology at a few kpc, and clumpy, highly asymmetric stellar morphology driving gravitational torques on the gas is evident on all scales.In most panels we show a continuous projection of stellar density, but in the last panel this breaks down (the inter-stellar separation is no longer much smaller than a pixel) so we show individual O-stars.We do not show images at ≪ 1 pc because there is a negligible stellar mass compared to the gas on these scales.2020a), and the (dust-dominated) protostellar disk and molecular cloud limits (Grudić et al. 2021), we have also validated our opacities against those tabulated in Lenzuni et al. (1991) for metal-free gas with densities  ∼ 10 12 − 10 16 cm −3 .
We emphasize that all the physics and numerical methods above, including gravity, radiation transport, MHD, and thermochemistry, apply always and everywhere in the simulation: there is no distinction between FIRE and STARFORGE treatments.6The only difference between the FIRE and STAR-FORGE limits in our simulations lies in how we treat "stars."Specifically, when a gas cell is eligible for "star formation", we must decide whether to convert it into a "single stellar population (SSP) particle" which represents a statistically-sampled ensemble of multiple stars (the FIRE limit, relevant at lower resolution/large cell masses) or to convert it into a "sink/single (proto)star particle" which represents a single (proto)star (the STARFORGE limit, relevant at higher resolution/small cell 6 Again note that some previous FIRE simulations using the "default" model in Hopkins et al. (2018b) employ a simpler approximate LEBRON radiationhydrodynamics solver, and FIRE-2 used a simpler thermochemistry module compared to the FIRE-3 version here.However these simplifications are not designed for handling extremely high densities or optically thin-to-thick transitions, so we adopt the more accurate M1 RHD and FIRE-3/STARFORGE thermochemistry detailed above.But we stress that these RMHD and thermochemistry modules have been used (and compared to those simpler modules) in multiple previous FIRE studies (e.g.Hopkins et al. 2022aHopkins et al. , 2020a;;Hopkins 2022;Hopkins et al. 2022f;Shi et al. 2022;Schauer et al. 2022, and references therein) as well as STARFORGE (Guszejnov et al. 2021;Grudić et al. 2021Grudić et al. , 2022b;;Guszejnov et al. 2022b,c,a), and many additional simulations using GIZMO, referenced above.masses).Those particles each then use their distinct appropriate (SSP or single-star) evolutionary tracks to calculate the mass/momentum/energy/cosmic ray/photon fluxes which are deposited back onto the grid, at which point that injected material is again evolved identically according to the algorithms above.Thus both operate simultaneously in the simulation.

FIRE Treatment of Stars (Relevant for Large Cell
Masses) FIRE was designed for galaxy-scale simulations, with resolution sufficient to resolve some phase structure in the ISM, but insufficient to resolve individual proto-star formation and stellar growth/evolution histories.As such, we apply the FIRE treatment of stars when the resolution is still sufficiently low (cell mass > 1  ⊙ ), using the FIRE-3 implementation in (Hopkins et al. 2022a) and summarized above.In this limit, gas is eligible for star formation if it is locally self-gravitating at the resolution scale, Jeans unstable, and in a converging flow, as in Hopkins et al. (2013b); Grudić et al. (2018).The intention here is not to resolve e.g.local "peaks" in the density field which will become individual stars, but rather to identify "patches" of the ISM where the fragmentation cascade becomes unresolved, so the gas should continue to fragment and ultimately form a population of stars.As such, for the cells which meet this criterion, we assume fragmentation on a dynamical time (per Hopkins et al. 2018b) and convert them into "single stellar population (SSP) particles" -i.e.collisionless particles which represent ensemble populations of stars which sample an assumed universal stellar IMF.
Once formed, these SSP particles evolve as detailed in Hopkins et al. (2022a) according to explicit stellar evolution models from the 2021 version of STARBURST99 (Leitherer et al. 2014), and return metals, mass, momentum, and energy to the ISM via resolved individual SNe (both Ia & core-collapse) and O/B and AGB mass-loss as in Hopkins et al. (2018a), with radiative heating and momentum fluxes determined from the stellar population spectra as in Hopkins et al. (2020a), appropriate for a Kroupa (2001) IMF.Cosmic rays are injected in fast stellar wind or SNe shocks as described in Hopkins et al. (2022a), using the approximate method from Hopkins et al. (2022b).Once injected onto the grid, gas/radiation/cosmic rays obey the physics in § 2.1.
To deal with intermediate resolution cases, we employ the stochastic IMF sampling scheme from Ma et al. (2015); Su et al. (2018b); Wheeler et al. (2019); Grudić & Hopkins (2019): when a SSP particle forms, we draw a quantized number of massive stars from an appropriate binomial distribution, from which the relevant feedback properties particular to (rare) massive stars (e.g.core-collapse SNe, ionizing radiation) scale.This allows us to at least approximately apply SSP particles down to gas mass resolution ∼ 1 − 10  ⊙ , where the resolution is still too poor to explicitly resolve individual (proto)star formation, but so high that the expected number of massive stars per star particle is ≪ 1 and as such discreteness effects could be important (see discussion in Ma et al. 2015Ma et al. , 2016bMa et al. , 2020b;;Grudić & Hopkins 2019).

STARFORGE Treatment of Stars (Relevant at Small
Cell Masses) STARFORGE, on the other hand, was designed for simulations which resolve individual (proto)star formation and evolution, e.g.simulations of individual molecular clouds, clumps, star clusters, or (proto)stellar disks.In this limit, each sink represents a single star, which obviously means it cannot meaningfully represent systems with resolution poorer than ≳ 1  ⊙ .As such, we apply the STARFORGE treatment of stars when the resolution is sufficiently high (cell mass < 1  ⊙ ).In this limit, gas is eligible for (proto)star formation if it meets a standard but more stringent set of seed criteria described in Grudić et al. (2021), including a strict virial/selfgravity, Jeans, converging flow, fully-compressive tides, and local density/potential maximum criteria as well as restricting to gas cells without a pre-existing neighboring sink and requiring their collapse time be much shorter than the infall time onto the nearest sink (whatever its distance).If all of these criteria are met, a cell is immediately converted into a sink or individual star particle.
Once formed, each sink accretes gas that is bound (accounting for thermal, kinetic, and magnetic energies) to it, and whose current and circularized radii fall within the sink radius (set comparable to the force softening), following a standard strict sink accretion model validated in a variety of idealized accretion problems (details in Grudić et al. 2021).The sinks evolve along combined proto and main-sequence stellar evolution tracks, explicitly following the stellar evolution physics versus time (e.g.contraction, heating, different burning stages) allowing for the dynamic accretion rate in every timestep (Offner et al. 2009).In the proto-stellar evolution stage, sinks radiate in all bands with the appropriate effective temperature, and launch collimated jets with a mass-loading proportional to the surface accretion rate onto the star and a launch velocity comparable to the escape velocity from the protostellar surface (details in Grudić et al. 2021;Guszejnov et al. 2021;Grudić et al. 2022b).Main-sequence stars continue to emit jets and accretion luminosity if accretion continues, and radiate in all bands following their stellar evolutiontrack calculated effective temperatures and full spectra, while also emitting continuous stellar surface winds (assumed to be isotropic in the frame of the star, with a continuous mainsequence mass-loss rate given by Grudić et al. 2022b Eq. 1), as a function of the instantaneous stellar luminosity.At the end of their main-sequence lifetime stars can, if sufficiently massive, explode as SNe with 10 51 erg of ejecta kinetic energy.Again we refer to Grudić et al. (2021) for details.
Importantly, we note that these physics have been validated by direct comparison to dense molecular gas properties, the observed stellar IMF and multiplicity statistics, and star cluster properties (Guszejnov et al. 2021(Guszejnov et al. , 2022b;;Grudić et al. 2022b;Lane et al. 2022), for typical Milky Way-like galaxy conditions.In Hopkins (2023a) (henceforth Paper III) we will study the predicted IMF from the simulations here around the SMBH; for now we note that although it becomes top-heavy very close to the SMBH (≪ 1 pc), this analysis appears to justify the choice of a more universal Kroupa (2001)-like IMF at the larger radii ≳ 10 pc assumed by our FIRE SSP particles.

Initial Conditions and Refinement Choices
Our initial condition (see Fig. 2) is a fully-cosmological "zoom-in" simulation, evolving a large box from redshifts  ≳ 100, with resolution concentrated in a ∼ 10 Mpc comoving volume centered on a "target" halo of interest (specifically, halo "A1" aka "m12z4" studied in Feldmann et al. 2016Feldmann et al. , 2017;;Oklopčić et al. 2017;Anglés-Alcázar et al. 2017c;Hopkins et al. 2020b;Ma et al. 2021;Wellons et al. 2022).While there are many smaller galaxies in that volume, for the sake of clarity we focus just on the properties of the "primary" (i.e.best-resolved) galaxy in the volume.The dynamic refinement scheme employed here is numerically identical to that used in many previous GIZMO studies,7 including examples which have refined to similar resolution around single or binary SMBHs, just without the use of the hybrid FIRE-STARFORGE physics described above (see Orr et al. 2018;Hopkins et al. 2018b;Su et al. 2019Su et al. , 2020Su et al. , 2021;;Benincasa et al. 2020;Anglés-Alcázar et al. 2021;Franchini et al. 2022).
The simulation is run from  ∼ 100 down to some redshift (here  < 4), using a dynamic refinement scheme that adaptively varies the mass resolution, as a function of the minimum of either the thermal Jeans mass (ensuring it is resolved by ∼ 100 cells) or a function of distance to the nearest SMBH particle (progressively moving from minimum refinement at > 100 kpc from the nearest SMBH to maximum refinement at < 10 kpc), between an imposed minimum refinement cell mass of ≈ 4000  ⊙ and maximum refinement mass of ≈ 10 6  ⊙ .To ensure even low-density, non-self-gravitating multi-phase structure is reasonably resolved, once refined a cell cannot be de-refined unless it escapes far from the galaxy (though 7 Briefly, gas cells are assigned a target resolution Δ target as a function of e.g.BH-centric distance: if they deviate by more than a tolerance factor ∼ 2 shown below, they are refined in pairs (or de-refined via merging with their nearest neighbor if both are eligible for de-refinement) maintaining all numerically-integrated conserved quantities (mass, momentum, energy, etc.) to machine precision, with the refinement into pairs of equal-mass Lagrangian cells separated along the axis of lowest pre-existing neighbor cell number density (akin to an approximate local Voronoi-type remeshing).For details see Hopkins (2015); Su et al. (2019); Anglés-Alcázar et al. (2021).The full algorithm is provided in the public GIZMO code. in practice, given the short simulation duration after reaching maximum resolution and relatively weak winds from the central sub-pc scales, very few cells of high refinement escape to interact with low-resolution cells).This allows us to resolve the galaxy at ∼ 4000  ⊙ resolution through formation and initial growth of the SMBH, through a redshift of  ≈ 4.
We then select a specific time  0 (at a redshift  0 ≈ 4.4) from this original simulation just before a period where it (at its relatively modest resolution) identified rapid quasar-level SMBH growth, with one clearly-dominant SMBH particle within the galaxy nucleus.We will show the galaxy properties at this time in great detail below, but to summarize at  ≈ 4.4 it has a dark matter halo mass of ∼ 3 × 10 12 M ⊙ inside  < 250 kpc, and a galaxy stellar mass of 2 × 10 10 M ⊙ (and very similar gas mass) inside < 10 kpc (with stellar half-mass radius of ∼ 1.5 kpc), and a nuclear SMBH mass of ∼ 10 7 M ⊙ .We then re-start the simulation from this time, with an additional refinement layer: on top of the refinement scheme above, we apply a multiplier  ( ≡ |x−x SMBH |, − 0 ) to the "target" mass resolution, which is a continuous function of  with 0 <  ≤ 1, where at small radii  ∝  3 and at radii ≳ 1 kpc,  = 1 exactly.8We reduce the minimum allowed/target cell mass to Δ ≲ 0.01  ⊙ .
To avoid pathological behaviors, this refinement layer does not "instantly" activate but appears as a smooth function of time  −  0 , designed such that each concentric radius  (beginning at ∼ 1 kpc interior to which the refinement begins) is evolved for a few local dynamical times before the next "layer" of refinement is applied interior to this.Specifically, the refinement is still continuous in both space and time.For each decade in  (between  max and  min =  max /10) until we reach our maximum target resolution, if the desired final refinement function  (,  → ∞) ∝    , we take  ∝  ( ) with () =   (1 − exp {( −   /Δ)})/(1 − exp {(  −   /Δ)}) after time   equal to   of the previous (decade-larger) interval, with Δ ∼ 1/Ω (at the midpoint of the annulus) and   −   ∼ 5 Δ, with the normalization of  set to the value from the previous annulus at  max to ensure continuity, and  = constant at  <  min for  <   .In § 3.2 below, we explicitly show and discuss the effective final "ramp" refinement functions in both space and time.We also enforce a requirement that no cell can refine or de-refine within less than 30 numerical timesteps of its last refinement/de-refinement (where the numerical timesteps ensure this is some multiple of the relevant signal-crossing timescale).This gradual "ramp" of the refinement in both space and time is designed to reduce initial noise and spurious features imprinted by e. Fig. 4.-Illustration of the time evolution of the main galaxy in our simulation before our hyper-refinement.We plot the galaxy-integrated SFR  * , and the sub-grid estimated BH accretion rate (BHAR) from the model -which scales approximately as  subgrid ∼   gas Ω at the low-resolution limit of ∼ 10 − 100 pc before hyper-refinement is turned on -as a function of time prior to the time of refinement (at  ∼ 4.4).The inset shows the resolved gas inflow rate into the central sink resolution around the SMBH of < 80 au, as a function of time in units of the dynamical time  dyn ≡ 1/Ω at this resolution scale (∼ 80 au), over the final ∼ 1500 yr of our simulation duration (well after it reaches the maximum refinement level everywhere).Though this is a very short relative timescale (compared to the order Hubble time evolution on large scales), we see that the inflow rate into the inner accretion disk is quite stable over tens of thousands of dynamical times in the center, for a given set of conditions at larger radii.We also see extremely high inflow rates, as expected based on the high nuclear gas masses and densities in the "parent" simulation which motivated the choice of this particular moment in time to "zoom in." olution by the radii and times desired for zooming in on the formation of the quasar accretion disk.In summary, this means that the total duration of the simulation after the beginning of the "hyper-refinement" period is a couple of Myr, but after the highest resolution level (resolving orbits at ≲ 100 au around the SMBH with an orbital dynamical time of ∼ 10−20 days) is reached at a redshift closer to  ∼ 4, we evolve for ∼ 10 4 years.

Types of Cells/Particles
In summary, our simulations include five types of cells or particles: 1. Gas/Radiation Cells: These define the effective mesh on which the equations of radiation-magneto-hydrodynamics are solved, including thermochemistry and all the physics described above.The mesh resolution (for gravity and all other forces) is adaptive with spatial resolution Δ = (/Δ) 1/3 , and Δ ranges smoothly from ≲ 0.01  ⊙ in the high-resolution region around the SMBH (after the "hyperrefinement" phase begins) to a median of ∼ 5000  ⊙ in the ∼ kpc-scale ISM of galaxies to ∼ 10 6  ⊙ in the diffuse IGM.The physics/equations integrated for the gas are independent of resolution -the choice of FIRE or STAR-FORGE physics only appears as a choice of whether cells should form SSP or sink/single-star particles.
2. Dark Matter Particles: Dark matter is represented in standard fashion by collisionless particles which interact only via gravity.In the low-resolution cosmological box (well outside of the high-resolution region) these particles have lower resolution in factor-of-two increments (with the poorest resolution ∼ 4 × 1010  ⊙ in the ∼ 100 Mpc region), but following standard practice we have confirmed that within ∼ 500 kpc of the "target" galaxy of interest, there are zero low-resolution dark matter particles.The high-resolution dark matter particles have Δ ∼ 10 6  ⊙ and a force softening equivalent to Δ ≈ 200 pc.While crucial for cosmological evolution and galactic-scale dynamics, the dark matter contributes negligibly to the nuclear dynamics (contributing just ∼ 6% of the total mass inside < 200 pc and a vanishingly small fraction of the mass inside < 20 pc),9 and the force softening is sufficiently large that the worstscale -body deflection from an encounter between a highresolution baryonic cell and DM particle would be no larger than the acceleration/deflection for a gas cell with density  ≲ 1 cm −3 (vastly lower density/acceleration scales than those of interest in the galaxy nucleus).103. SMBH Particles: SMBHs are represented by collisionless sink particles.In the "pre-simulation" phase (running from redshift  ∼ 100 to the time  0 ( ∼ 4.4) when we begin our hyper-refinement), the BHs are formed and evolve according to the default sub-grid FIRE-3 seeding, dynamics, accretion, and feedback models described in Hopkins et al. (2022a).But this is only relevant in that it gives us a plausible initial condition for our hyper-resolution run.Once the hyper-refinement phase begins, we disable all "sub-grid" models for BH growth and accretion: the BHs are represented by sinks that follow normal gravitational dynamics.Any cells/particles which fall inside of the SMBH sink radius set to ≈ 80 au ∼ 300   (where   = 2   BH / 2 ) are immediately captured (removed from the domain and added to the sink mass).11We do not include SMBH "feedback" from within the sink radius during this phase.For simplicity, we choose a time where the primary galaxy of interest contains just one SMBH: a sink with mass ≈ 1.3 × 10 7  ⊙ .
4. Single Stellar Population (SSP) Particles: Unresolved fragmentation to stars formed in the "FIRE limit" ( § 2.2) produce SSP particles which represent stellar populations.The mass of an SSP particle is fixed to the mass of the gas cell from which it forms, so varies smoothly with radius accordingly.The force softening is fixed for each SSP particle after it forms but depends on its mass with the scaling Δ ≈ 0.6 pc (Δ/100 M ⊙ ) 1/3 , equivalent to a preformation gas density ≈ 2 × 10 4 cm −3 typical of where cells of this resolution form stars.12These particles act on the surrounding medium via continuous stellar mass-loss, radiation, and SNe, sampling from a universal stellar IMF, as described in § 2.2.
5. Sink/Single (Proto)Star Particles: Resolved collapse to individual stars in the "STARFORGE limit" ( § 2.3) produces sink or individual (proto)star particles, each of which represents a single stellar object or system.The mass of a sink reflects the current mass of the star: at formation, the initial mass equals that of its progenitor gas cell but it grows via resolved accretion.The force softening/sink radius for these is set to a fixed value with a Plummerequivalent  ≈ 25 au (Δ ≈ 40 au).These stars follow individual evolution tracks and act on the surrounding medium via radiation, jets, mass-loss, and SNe as described in § 2.2.
3. RESULTS We summarize some of the results for our fiducial simulation in Figs.2-3, showing images of the simulation on scales from > Mpc to < 100 au.Specifically, Figs. 2 & 3 show images of the projected simulation gas and stellar mass densities, viewed from the same viewing angle (chosen to be face-on to the accretion disk in the center), on a range of scales.Fig. 4 illustrates the pre-history of the large-scale "parent" simulation, for reference.

Different Characteristic Scales/Regimes
We can clearly see in Figs.2-3 that the dynamic range spanned by the simulation is enormous -a factor of > 109 in black-hole centric radii, and more like ∼ 10 13 if we compare our smallest spatial resolution at radii ∼ 10 −3 − 10 −2 pc from the SMBH, to the size of our entire cosmological box.Fig. 5 shows an alternative illustration, plotting the gas and stars on a logarithmic scale and identifying the different scales with the labels below, in an attempt to visualize the qualitative structures and phases of gas on each scale.It is difficult to actually describe so many orders of magnitude in scale at once, so we break the scales from 10 −3 − 10 6 pc in SMBH-centric radius down into each order-of-magnitude and both assign a characteristic label for these scales and describe some of the key physics and processes occurring.From largest to smallest scales, we follow gas inflows as follows: • IGM → CGM: On scales ≫ 100 kpc, the IGM is "cool" (temperatures ∼ 10 4 K), diffuse ( ≪ 10 −2   cm −3 ), quasi-spherical (/ ∼ 1), dark-matter dominated, weakly-magnetized ( plasma ≡  thermal / magnetic =    /(|B| 2 /8) ≫ 100, with || ∼ 1 − 10 nG), with weak outflows and strong primarily-radial (so Π  ≡ ⟨     ⟩ dominates the kinetic stress tensor) super-sonic inflows of ∼ 300 M ⊙ yr −1 onto the halo.Essentially gas is in freefall collapsing with dark matter via the cosmic web.Since this has been well-studied and resolved in many previous simulations, it is not our goal to study this regime in detail here, but what we see is consistent with the most previous studies with the FIRE simulations (Hafen et al. 2019(Hafen et al. , 2020;;Butsky et al. 2020;Hopkins et al. 2021b;Ponnada et al. 2022;Ji et al. 2020Ji et al. , 2021;;Li et al. 2021;Stern et al. 2021;Esmerian et al. 2021;Kim et al. 2022;Butsky et al. 2022) et al. ( 2022) is several orders of magnitude smaller than the typical turbulent dissipation rates in the simulation.(defined so  = 0 corresponds to the midplane of the inner disk), in a wedge of azimuthal opening angle sin  < 0.3.For gas (left) colors denote different phases  < 10 3 K (green), 10 3 <  < 10 4 K (yellow), 10 4 <  < 10 5 K (magenta), 10 5 <  < 10 6 K (purple),  > 10 6 K (cyan).We see other galaxies on IGM scales, the virialized CGM with accretion in warm clumps/filaments, the highly clumpy/inhomogeneous/asymmetric and multi-phase structure in the galaxy and (thermally colder, primarily atomic+molecular) galaxy nucleus, settling into the more ordered (but still visibly thick and turbulent) non-star forming disk and BH accretion disk on sub-pc scales.
as well as results from other codes and semi-analytic models (Hummels et al. 2019;Pandya et al. 2022) and standard observational inferences (see Tumlinson et al. 2017;Chen et al. 2020;Lan & Prochaska 2020, reviews).
• CGM → Galactic ISM: On scales ∼ 10 − 100 kpc, the volume-filling gas in the CGM is shock-heated to virial temperatures ∼ 10 6 K, with  plasma ∼ 100 and trans-sonic or sub-sonic turbulence, mostly ionized, with thermal pressure comparable to the total pressure and gravity.But the gas is multi-phase, with accretion and outflows of comparable magnitude, with outflows prominent in the diffuse/volume filling phases and inflows dominated by accretion of "cool" (≲ 10 5 K) gas along filaments with densities ∼ 100 times larger than the median background hot gas, lower  plasma ∼ 10, and velocities of order the free-fall speed in a darkmatter dominated potential.This is essentially the classic "cold flows in hot halos" picture, again consistent with many previous theoretical studies (Kereš et al. 2005;Dekel & Birnboim 2006;Brooks et al. 2009;Kereš et al. 2009a,b;Faucher-Giguère et al. 2011;Sĳacki et al. 2012;Kereš et al. 2012;Vogelsberger et al. 2012;Stern et al. 2020) and more recent observations (Ribaudo et al. 2011;Kacprzak et al. 2012;Vayner et al. 2022).
• Galactic ISM → Galactic Core/Proto-Bulge: On scales 1 − 10 kpc in the galaxy, the gas is highly multi-phase with self-shielding of the UV radiation field (Σ gas ≳ 10 M ⊙ pc −2 ) allowing formation of "cold" neutral medium (CNM) and molecular medium with  ≪ 10 4 K, alongside hot gas with  ≳ 10 7 K from SNe, while gas densities range from ≲ 10 −2   cm −3 to ≫ 10   cm −3 in cold cloud complexes (and  plasma similarly ranges from ∼ 0.1 in cold phases to ∼ 1 − 10 in warm phases and ∼ 100 in the most diffuse volume-filling phases, and || ≳ a few G).These cold complexes maintain most of the SF, with a SFR inside < 10 kpc of ∼ 50 − 100 M ⊙ yr −1 (over the last ∼ 100 Myr).
The potential becomes dominated by stars inside a few kpc (the galaxy effective radius).The turbulence is mildly supersonic (sonic M  ∼ 1− a few) in a volume-averaged sense (with the volume-average dominated by warm ionized media [WIM] and warm neutral media [WNM] at ∼ 10 4 K), but highly super-sonic (M  ∼ 10 − 100) in the "cold" phases.Most of the gas is atomic or molecular.While turbulence maintains an effective volume-averaged  ∼ a few as at all larger radii, the thermal Toomre  parameter drops to ≪ 1 in the cold phases, in particular, meaning that fragmentation via self-gravity is rapidly promoted, with the characteristic "most unstable" fragment masses expected to contain most of the power in the fragment mass spectrum (e.g. the largest self-gravitating complexes) ranging from ∼ 10 7 to a few 10 9  ⊙ (larger than in low-redshift galaxies, owing to the massive gas content of this dense, high-redshift galaxy, similar to complexes observed at high redshift).Again this is broadly consistent with previous theoretical (Noguchi 1999;Bournaud et al. 2008;Agertz et al. 2009;Dekel et al. 2009;Ceverino et al. 2010;Hopkins et al. 2012b;Oklopčić et al. 2017) and observational (Elmegreen et al. 2004;Martínez-Sansigre et al. 2009;Kriek et al. 2009;Daddi et al. 2010;Förster Schreiber et al. 2011;Newman et al. 2012) studies of massive star-forming and quasar-host galaxies at redshifts  ≳ 2. The system is extremely inhomogeneous, with non-axisymmetric mode amplitudes | 1 | ∼ 0.1 − 1 and large clump and cloud complexes and star clusters visible.
At the time of this particular simulation, the torques from ∼ 1 − 10 kpc clearly involve large non-axisymmetries which are visually dominated by a large minor merger (with the companion at ∼ 10 kpc, having just passed pericenter).
• Galactic Core/Proto-Bulge → Galactic Nucleus: On scales ∼ 0.1 − 1 kpc, most of the "Galactic ISM" intuition applies, and a significant fraction of the SFR (∼ 20 M ⊙ yr −1 averaged over the last ∼ 100 Myr, but closer to ∼ 100 M ⊙ yr −1 in the last ∼ 1 Myr) comes from these radii.The gas is primarily molecular (by mass) and the turbulence is supersonic and super-Alfvénic (M  ∼ M  ∼ 3 − 10 ).The potential is deeper, and the density and surface density scales are higher, with the surface density approaching the scales ∼ 10 3 M ⊙ pc −2 at which stellar feedback becomes less efficient (Fall et al. 2010;Grudić et al. 2018;Grudić & Hopkins 2019;Grudić et al. 2020;Kim et al. 2018b;Hopkins et al. 2022e) so the outflows weaken again relative to inflow, and the most massive cloud complexes are more like ∼ 10 6 − 10 7 M ⊙ .For the first time the radial Π  component does not strongly dominate the stress, as various clumps and accreting gas are not on primarily radial orbits as they are when accreting at larger radii but on quasiisotropic (often tangential) orbits, having quasi-circularized though with different angular momenta (as there is no coherent disk).Accretion is strongly dominated by gravitational torques with | 1 | ∼ 1 -i.e.order-unity asymmetries in the potential dominated by the stellar structure (since this dominates the mass) leading to the gas structures shocking and losing angular momentum on a timescale comparable to the dynamical/orbital time (see Levine et al. 2008;Hopkins & Quataert 2010b, 2011b;Hopkins et al. 2016;Anglés-Alcázar et al. 2013, 2017a, 2021;Prieto & Escala 2016;Prieto et al. 2017).The system begins to be optically thick to cooling radiation in NIR/optical/NUV/UV bands, so the IR radiation energy density begins to rise.
• Galactic Nucleus → Black Hole Radius of Influence (BHROI): On scales ∼ 10 − 100 pc, the increasing density and surface density scale ( ∼ 10 2 − 10 6   cm −3 , Σ gas ∼ 10 3 − 10 4 M ⊙ pc −2 ) means that the gas cools rapidly and the "hot" and "warm ionized" phases vanish rapidly so by 10 pc the gas has an average temperature ≲ 1000 K (primarily in relatively "warm" molecular gas).This also means  plasma drops from ∼ 1 at the outer range of these radii to ≲ 0.01 at the inner radii.The system has moved above the critical Σ eff at which stellar feedback becomes highly inefficient and we see outflows diminish.But strong instability, fragmentation to more "GMC-like" mass scales, inhomogeneity with | 1 | ∼ 1, and highly super-sonic (and still super-Alfvénic) turbulence (M s ∼ 30 even in the volumefilling phases) persists.The region is still actively starforming but given the smaller area and mass contributes only ∼ a few M ⊙ yr −1 in "steady-state", but this is spiking to ∼ 10 − 100 M ⊙ yr −1 in the ∼ 1 − 10 Myr immediately preceding the snapshots analyzed owing to the rapid inflows.
At the time of our hyper-refinement, the inflows through these scales are dominated by one large gaseous complex (itself torqued by the ongoing mergers above) which has a close passage with the BHROI13 allowing the BH to tidally capture the gas -gravitational torques still clearly dominate in this regime.The system is beginning to become more optically thick at some wavelengths but still has cooling times much shorter than dynamical times and is not in a blackbody like state (the dust, radiation, and gas temperatures are all significantly different).
• BHROI → "Torus": On scales ∼ 1 − 10 pc, the BH begins to dominate the potential, though stars still strongly dominate over gas in the local fluctuations in the potential (since the density of stars is much higher than gas).
Because the system is now "fully" optically thick to its own cooling radiation, we begin to see a clear inversion of the density-temperature relation, with denser gas being warmer (in both its kinetic, gas, and radiation temperatures, even though these are not yet all in equilibrium with one another) in a quasi-adiabatic manner (as opposed to the usual case at larger radii where denser gas is colder), with  plasma ∼ 0.01.The densities in the midplane and dense gas phases begin to exceed ≳ 10 6   cm −3 , at which point the dust temperature starts to couple appreciably to the gas kinetic temperature so the two begin to approach one another, but the large inhomogeneity of the medium and much shorter dynamical times (compared to e.g. the conventional case in molecular clouds) mean this coupling is still relatively weak/gradual and incomplete.Despite inflow rates still as large as ∼ 100 M ⊙ yr −1 , the SFRs averaged over the last (100, 10, 1) Myr are ∼ (0.5, 5, 25) M ⊙ yr −1 , indicating its highly non-global-equilbrium nature.Turbulence remains highly super-sonic (in primarily warm molecular gas) but becomes only mildly super-Alfvénic with M  ∼ 1 − 3. Again, gravitational torques clearly dominate the visual structure (with large coherent gas asymmetries and | 1 | ∼ 1): the major change from larger radii is that instead of being incoherent/clumpy structures, the increasingly Keplerian nature of the potential means that coherent, non-linear  = 1-like perturbations related to torques between gas and stars dominate (Tremaine 2001;Zakamska & Tremaine 2004;Hopkins & Quataert 2010b, 2011b,a;Hopkins 2010).
• "Torus" → Non-Star-Forming Disk: As shown in § 3.4, on scales ∼ 0.1 − 1 pc, the temperatures rise to a few 10 3 K, and molecules begin to dissociate into atomic gas at these warmer temperatures (though a non-negligible molecular fraction remains).The gas is still dust-laden (we are outside the sublimation radius).We see an important morphological transition in Figs.2-3: the stream of gas tidally torn from the external gas complex and fueling the accretion disk begins to circularize and self-intersect, forming a more coherent disk.Compared to larger radii, we see ( § 3.5) an accompanying rapid rise in  turb and  mag .At  ≳ pc, we still have  thermal ≲ 1 (and a cooling time short compared to the dynamical time), and as shown in previous studies (see e.g.Hopkins & Christiansen 2013 and discussion below) a disk with these conditions is still unstable to gravitational fragmentation within "patches" even if it is statistically marginally stable with turbulent+magnetic support (the system still has  ≪ 1, with trans-Alfvénic turbulence), so some star formation continues.Indeed, the SFR averaged on (100, 10, 1, 0.1) Myr timescales inside ∼ 1 pc is ∼ (0.01, 0.1, 1, 10) M ⊙ yr −1 (to the extent that it can be defined in any meaningful way on these small timescales), as the individual protostars and main-sequence stars formed since the gas arrived at these radii are still accreting.However, as we approach ≲ 0.1 pc, the system dramatically changes and star formation effectively ceases.
• Non-Star-Forming Disk → "Accretion Disk": (∼ 0.01 − 0.1 pc) Just outside ∼ 0.1 pc, a crucial transition occurs as  thermal increases to ≳ 1 with  mag ≫ 1 dominated by increasingly-organized toroidal fields.Meanwhile, the characteristic maximal fragment mass ∼  Σ gas  2 starts to drop into the stellar mass range.As a result (discussed in detail below), star formation shuts down.The SFR inside < 0.1 pc averaged over the entire duration of our simulation is ≲ 10 −3 M ⊙ yr −1 , compared to inflow rates of ∼ 10 − 100 M ⊙ yr −1 .With this transition, the disk mass inside < 0.1 pc is now locally gas-dominated instead of stellardominated, so gravitational torques rapidly become inefficient (Hopkins & Quataert 2011b), but we still see  = 1 modes propagate into these radii and gravito-turbulent behavior, but now a combination of Maxwell and Reynolds stresses take over as the dominant provider of the torques maintaining a similar global bulk inflow rate.As we go to smaller scales still, the deepening potential means the scale height becomes somewhat smaller and the disk becomes increasingly well-ordered.The disk is strongly-magnetized, with  plasma ≲ 10 −3 , and the turbulence becomes modestly sub-Alfvénic at smaller radii.
• "Accretion Disk" → ISCO: On scales ≪ 0.01 pc (≪ 10 4 gravitational radii), the disk is essentially in the regime of a "traditional" -like accretion disk in many ways.It is optically thick, geometrically thin or "slim," radiating in increasingly black-body-like fashion, nearly-Keplerian and close to circular, gravitationally stable ( thermal ≳ 1 with  mag ≫ 1), with maximum fragmentation mass scale ≲ 1  ⊙ so it is not able to fragment efficiently at all.But there are many important differences between the disk here and what is usually assumed in accretion disk studies (to be studied in detail in Hopkins 2023b, henceforth Paper II).Even though the optical depth is large, the effective black-body cooling time is much shorter than the dynamical time (by a factor ∼ 10 −3 ), and the turbulence is supersonic, so it maintains a quasi-isothermal, relatively cool global structure.The disk is strongly magnetized with  plasma ∼ 10 −4 (|| ≳ 100 G primarily toroidal fields), sustained by flux-freezing from the flux it is fed from the ISM (hence a "flux-frozen" and/or "flux-fed" accretion disk) and modestly sub-Alfvénic (hence highly-supersonic) turbulence.

3.2.
Effective Resolution of the Simulation Fig. 6 shows the effective mass, spatial, and time resolution of the fiducial simulations as a function of BH-centric radius  at the times where we study it.This reflects the target resolution discussed in § 2 above.In brief: at galactic radii (≲ 10 kpc) the resolution is uniformly better than ∼ 10 4  ⊙ as given by our target refinement criterion before the "hyperrefinement" is activated; this then achieves the desired radial refinement with Δ smoothly decreasing from ∼ kpc to ∼ pc scales before we saturate at our target resolution inside ≲ 1 pc of Δ < 0.01  ⊙ .This also lets us clearly identify where the simulation lies in different "limits" with regard to star formation per § 2.2-2.3: at > 10 pc scales, the resolution is always in the "FIRE" limit (forming SSP particles), and at -Effective resolution of the simulation as a function of radial distance from the BH ().We show the median, 50% and 90% inclusion intervals (shaded) of the mass resolution (gas cell mass Δ; left), spatial resolution (equivalent cell size Δ ≡ (Δ/) 1/3 ), and time resolution (numerical timestep Δ), at each .The bulk of the galaxy is resolved with Δ < 10 4  ⊙ ; from 1 − 100 pc the resolution is rapidly refined as  decreases, with the target resolution of Δ < 0.01  ⊙ and Δ < 10 −3 pc reached for the particles at  ≲ pc scales.For Δ, Δ, and Δ we compare the global enclosed gas mass  gas ( ), radius , and dynamical time  dyn at each .In Δ, we also denote the region where the simulation will form STARFORGE single-star/sink particles, as opposed to FIRE SSP particles.In Δ, we also denote as "duration" approximately how long the simulation was run after reaching the maximum refinement level at each radius.
< 1 pc scales, the resolution is always in the "STARFORGE" limit (forming single-star particles).As shown in Grudić & Hopkins (2019), the "transition" between these limits, with mass resolution ∼ 1 − 10 M ⊙ , is awkward in the sense that individual stars and the IMF are clearly not well-resolved, but neither can one effectively sample massive stars from the IMF with a single particle (and resulting detailed stellar feedback efficiencies can depend on the statistical method used for IMF sampling at the factor ∼ 2 − 3 level),14 so we intentionally design our refinement scheme to interpolate through this intermediate resolution over a narrow range of radii ( § 2.4).We can also compare to the total enclosed gas mass in the simulation, to demonstrate that there are always  ≫ 1 gas cells in each radial annulus.
Likewise, we see the spatial resolution at small radii is uniformly ≪  (the scale height of the gas at that radius, defined below), reaching Δ ∼ 10 −5 pc, and the time resolution is always much shorter than the local dynamical time.Our timesteps reach extremely small values Δ < 1 day in the central regions at the maximum refinement level: even with hierarchical timestepping (obviously necessary for such large dynamic range) this limits how long the simulations can be run.

Mass and Accretion Rate Profiles
Figs. 7 & 8 more quantitatively examines the radial profiles of various quantities related to the mass and mass flows: the circular velocity (defined as  c ≡ √︁   enc (< )/) and its contribution from the SMBH, gas, stars, and dark matter; the radial profile of surface density Σ gas and mid-plane threedimensional density , and the inflow and outflow rates  through each annulus.15 Because  2  at some  is just proportional to the enclosed mass,16 we can clearly read off from Fig. 7 where different components dominate the potential and the local matter distribution.The BH dominates inside the BHROI at a few pc, and we see the local density is gas-dominated only interior to the radii where star formation shuts down (≪ 0.2 pc here), while stars dominate the local density from ∼ 0.2 pc to ∼ 2 kpc, and dark matter dominates the density at much larger scales ≫ kpc 14 We technically address this by simply adopting properties for the IMF sampling for whatever mass star "should" exist according to the draw from the IMF as per Su et al. (2018b); Wheeler et al. (2019), but limiting the returned mass to the total star particle mass to retain mass conservation.More sophisticated schemes exist, which will be explored in future work (Steinwandel et al., in prep.), but in practice this only occurs for a negligibly small fraction of star particles around ∼ 10 pc in the simulation here, and as a result we see no effects re-running briefly with a simpler model where we adopt IMF-averaged properties without sampling for these particles. 15 where Θ(  ) = 0 or 1 for  < 0,  > 1 identifies all the inflowing gas, and likewise for  out (with Θ(  )), in some sufficiently narrow annulus Δ (here ∼ 0.1 dex).This allows for there to be both inflow and outflow through a given radius, given the lack of spherical symmetry.We also define the SFR as the total mass locked into stars (via new (proto)star formation and accretion onto stars) in the last dynamical time in each annulus, as this gives a better sense of the "instantaneous" rate at small radii where the dynamical time is much shorter than usual fixed timescales like ∼ 10 − 100 Myr used observationally to define "recent" star formation.
16 Note at very large radii outside the halo virial radius,   defined here as √︁   enc / begins to rise as it should while the (mean) dark matter and gas densities approach constant values near cosmic mean, so  enc ∝  3 .10 −3 10 −2 10 −1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 Distance from BH r [pc] √︁   enc (<  )/ versus spherical distance from the SMBH , with contributions from different mass components.The BH excision radius truncates the gas mass distribution at ≲ 0.001 pc.We clearly see a transition from dark matter dominating outside the galaxy, stars dominating within the galaxy, until the BHROI at ∼ a few pc, and the cessation of star formation giving a gas-dominated disk at ≲ 0.1 pc.Top Right: Projected mean gas, stellar, dark matter, and "star formation rate" (defined by mass of stars formed in the last  dyn ≡ 1/Ω at each radius) surface density Σ, in cylindrical shells of different projected radii  (plotted down to radii where we have at least > 1000 particles of each type).We see the same transitions.While there are large deviations at any radius, Σ gas ∝  −1 approximates the scaling reasonably well from ∼ 10 −3 − 10 6 pc.In this and subsequent plots we independently determine the "midplane" in each annulus from the angular momentum vector of the gas in the annulus.The suppression of Σ SFR at small radii is discussed below ( § 3.5.3).Bottom Left: Three-dimensional gas density  versus spherical radius .We show both the volume-weighted density mean (line) in shells and mass-weighted mean (shaded range shows the volume-weighted 90% inclusion interval; the mean can exceed this upper limit if the distribution has large tails).The gas density profile is crudely isothermal (in a very order-of-magnitude sense),  ∝  −2 , from ∼ 10 −3 − 10 6 pc.Mass-weighted densities are systematically enhanced relative to volume-weighted but also vary more owing to clumping and phase structure.Bottom Right: Mass flow rate , showing the total inflow rate  in and total outflow rate  out through each annulus, and the cumulative SFR summed within each radius  * (<  ), versus spherical .Grey bar shows the 50% (solid; ∼ 18 − 30 M ⊙ yr −1 ) and 90% (dotted; ∼ 15 − 73 M ⊙ yr −1 ) range of inflow rates at < 100 au over the duration of our highest-resolution simulation.Lack of spherical symmetry and non-equilibrium dynamics mean that inflow and outflow co-exist and can change relative sign (e.g. the "outflow" at ≲ pc scales is mostly just coherent eccentric motion), but despite this a remarkably stable order-of-magnitude inflow rate to the SMBH of  in ∼ 10 − 100 M ⊙ yr −1 persists.
(while the steep profile of stars and gas mean dark matter is completely negligible at smaller radii, for any reasonable dark matter profile).While there are clearly very large local fluctuations in gas density, it (rather remarkably) appears to follow an approximately isothermal-sphere-like  gas ∝  −2 profile on average over nine decades in radius.17This leads to a gas surface density profile scaling approximately as Σ gas ∝  −1 (an actual power-law fit gives a very slightly shallower slope).
The density profiles and   profiles of stars, gas, and dark matter also clearly demonstrate that, despite the chaotic dynamics of the galaxy at this time, there is a reasonably well defined galaxy "center" (with the SMBH within it) on ≳ 10 pc scales (i.e. the stellar and gas and dark matter densities rise monotonically down to these radii), interior to which the SMBH itself dominates the potential and defines the "center."As discussed in Anglés-Alcázar et al. ( 2021); Hopkins et al. (2022e) and Byrne et al. (2023a), this is important for establishing the pre-conditions for a rapid accretion event, enabling gas flows to efficiently reach the SMBH radius of influence from much larger scales, so it should not be surprising to see .Left: "Depletion" (or accretion) timescales for inflow/accretion, outflow, and star formation, in different radial annuli, defined as shown (with Δ  ≡   / ln ), in units of the galactic dynamical time ( dyn ≡ Ω −1 ) at each radius.Star formation is efficient at Galactic radii from ∼ pc to ∼ kpc, as expected.Inflows are dynamical at essentially all radii, and produce net inflow at sub-kpc scales during this strong-inflow phase.Right: Gas mass-weighted mean radial (  ), azimuthal/rotational (  , with the φ axis defined in each annulus by the net angular momentum vector as Fig. 7), and polar (  ) velocities, and velocity dispersions   , in units of   .We see a mix of inflow/outflow motions (large   ) dominate at large radii (≫ kpc), partial rotation in a kinematically hot/thick galaxy configuration at Galactic radii (∼ 0.1 − 10 kpc), and near free-fall at the radii interior to the BHROI where the GMC-like gas complex fueling the SMBH is tidally disrupted and captured (∼ 1 − 10 pc), which then circularizes at sub-pc scales to form a kinematically cold, rotation-dominated disk at ≲ 0.1 pc scales.
here (given our selection of a massive accretion event to "zoom in" onto).
Recall, the duration of our simulation at its highest resolution is still short compared to the global dynamical/evolution timescales on ≳ 1 − 10 pc scales, so (as expected) these profiles are robust in time over the duration of the simulations.Even at the smallest radii, where we run for many dynamical times, after the initial refinement period they remain consistent within the scatter shown, as they are determined by the boundary conditions from larger radii.
In terms of accretion rates, we also see a surprisingly closeto-constant  in () from radii of ∼ Mpc down to ≲ 10 −3 pc.This is especially surprising given (a) the wildly different characteristic dynamical times on these scales, and (b) as noted from the morphologies above and some kinematic discussion below, that many radii are strongly out-of-equilibrium.The latter does produce some of the large "wiggles" in  in , but seems to produce much more dramatic variation in the outflow rates  out at different radii.That is consistent with e.g. the behavior seen in Anglés-Alcázar et al. ( 2021), especially when we consider the time variability shown in Fig. 7 over the dynamical time at each radius, but we caution that we focus on much smaller spatial scales for a much shorter overall period of time, compared to their study.Interestingly however, comparing the different simulations considered in Anglés-Alcázar et al. (2021), the (weak) variation in  in () we see is most similar to their "full-QSO" simulation (the simulation with the largest sustained inflow, most similar to the case here).That suggests the radial and time variability may be much larger at lower accretion rates (which is plausible, as e.g.star formation and outflows and other potential "bottlenecks" may play a much larger role limiting gas supply at low ).We do, on average, see some systematic decline in  in from the largest to smallest radii, as expected (material can "stall" and simply cease inflow without efficient angular momentum transport mechanisms, or be ejected in outflows, or go into star formation, at each radius), but this is weak, especially at the smallest radii ≪ pc where star formation has ceased (again notably weaker than simulations modeling systems with orders-of-magnitude lower mass inflow rates like M87, see e.g.Guo et al. 2022).And we see outflow rates are order-of-magnitude comparable to inflow rates at most radii; but even where  out >  in locally (which again clearly indicates out-of-equilibrium behavior, and is much more transient in these simulations) inflows are sustained over the entire duration of the simulation.It is also the case that even at large radii both the inflow/outflow rates are generally much larger than star formation rates within a given annulus (except right around ∼ 1 kpc), as shown in Fig. 8, owing to feedback self-consistently regulating star formation to be relatively slow (with an average efficiency ∼ 1% per free-fall time; see Hopkins et al. 2011;Orr et al. 2018), as expected from previous studies of gas-rich, star-forming galaxies.Finally we can also see where star formation ceases at small radii in both Figs.7  & 8.
Looking at different times in our simulations in Fig. 4, we see that while there are some significant (factor of a few to order of magnitude) variations in the accretion rates into the central < 80 au over the duration of the simulation, the accretion rates at these radii are quite slowly-evolving in a dynamical sense (these variations occur over tens of thousands of local dynamical times  dyn ∼ 1/Ω at the smallest radii).Thus it is reasonable to consider the inner regions to be in some kind of statistical quasi-steady-state in terms of accretion and dynamics, at given large-scale time in the galaxy.Of course, there are other peaks in the accretion rate into the central ∼ 10 − 100 pc in the simulation before hyper-refinement; we select the strongest inflow for refinement, but in future work hope to explore these other episodes.

Plasma and Thermo-Chemical Properties
In Fig. 9 we illustrate some of the multi-phase structure of the simulation more explicitly, alongside the (highly inhomogenous) stellar distribution.For comparison, Fig. 10 illustrates the average radial profiles (smoothing out these local variations) in various plasma and thermo-chemical properties of Thermal properties: we plot the volume-weighted and mass-weighted mean gas temperatures (and their range, shaded) in each annulus, together with the effective radiation temperature (integrating over all bands evolved, but excluding the CMB which dominates at  ≳ 200 pc) and dust temperature (also dynamically evolved), all versus distance from the BH in spherical annuli.Top Right: Chemical properties: average free electron (  ), ionized hydrogen ( HII ), molecular hydrogen mass (  H 2 ), atomic hydrogen ( HI ), and metal (, so solar ∼ 0.01) mass fractions.We see the highly multi-phase, optically-thin medium collapse into a cool atomic+molecular medium with warm dust in the galactic nucleus, then eventually see the system converge towards black-body like behavior with  rad ∼  dust ∼  gas in the center, with the dust sublimating.The scatter at a given  in e.g.  is large, and shown below.Bottom Left: Magnetic field strengths: we plot the energy-weighted rms field strength, and radial/toroidal/poloidal components.To crude approximation ⟨ |B| ⟩ 1/2 ∝  −1 from ∼ 10 −3 − 10 6 pc, with crudely isotropic (some mildly radial-dominant reflecting inflows) fields at most radii until the inner disk forms, and the field becomes primarily toroidal (see Paper II). Bottom Right: Radial profile (compensated by  2 to make comparison easier) of different components of the pressure/stress tensor (for anisotropic components, we plot ∥ ⟨⟩ ∥ for the stress ).At all radii, the kinetic energy density/ram pressure is important, and it is largely isotropic with mild radial bias at most radii until the disk forms where it becomes tangential.We see the transition from plasma  ≫ 1 at large radii to ≪ 1 at small radii.Radiation (again excluding the CMB term here) is generally sub-dominant at all radii, and other stresses (viscous, cosmic ray) are even smaller.the medium.
We see that the mean temperature jumps from typical ∼ 10 4 K IGM values to much warmer ≳ 10 6 K (comparable to the virial temperature) inside the virial radius of the dark matter halo as the gas shocks, although as shown in Fig. 9 much of the accretion onto the galaxy can still be in the form of warm/cold filaments and clumps.At these large radii the radiation and dust temperatures are largely determined by the CMB (if we specifically ignore the CMB, the radiation temperature of the residual radiation instead just reflects the meta-galactic UV background), because the medium is optically thin without significant sources on these scales.As expected, the gas is largely a mix of ionized and atomic phases.Inside the galaxy, we see an even more dramatic multi-phase structure (evident in e.g. the separation between mass and volume-weighted mean temperature), with large amounts of gas at ∼ 10 4 − 10 5 K (both ionized and warm atomic), but some cold neutral and molecular star-forming gas and a hot phase at ≳ 10 6 K.In the galaxy nucleus at scales, ≲ 100 pc, the mean temperature drops as the densities are so high that there is very little hot phase, and the medium becomes primarily molecular.As we go to smaller radii and the star formation rate density becomes higher and infrared optical depths become appreciable, we see the dust temperature rise from CMB values up to ∼ 100 K, very similar to the typical values observed in low-redshift starburst nuclei and circum-nuclear disks around AGN where the molecular gas and SFR densities are comparable to those predicted here (see e.g.Narayanan et al. 2005;Evans et al. 2006;Iono et al. 2007;Hopkins et al. 2008b,c;Casey et al. 2009;Wang et al. 2008;Izumi et al. 2016;Lelli et al. 2022, and references therein).We also see a significant "warm" molecular component at ∼ 1000 K begin to appear at ∼ pc.At  ≲ pc, most of the medium is at warm phase temperatures ∼ 10 3 − 10 4 K, and molecules begin to be dissociated again, and at ≲ 0.1 pc the optical depths to cooling radiation become large while the densities are sufficiently large (≫ 10 6 cm −3 ) that the dust and gas and radiation temperatures all begin to couple to one another (rapidly converging to broadly similar values by ∼ 0.01 pc).
Meanwhile, similar to what we saw with the density field, the mean magnetic fields follow a profile ⟨|B|⟩ ∝  −1 (becoming slightly steeper in the CGM/IGM, see Ponnada et al. 2022).18Note that, because of the extreme dynamic range here, the fact that |B| scales slightly steeper than  −1 , while  scales slightly shallower than  −2 means that the mean ideal-MHD Alfvén speed (  = (|B| 2 /4) 1/2 ) is not exactly constant but increases gradually from tens of km s −1 on "galactic" scales ∼ 1 − 10 kpc to hundreds of km s −1 at scales ≪ 0.01 pc, but this is consistent with a very weak trend ⟨  ⟩ ∝  −0.15 or so.At radii ≫ pc we see large variance in |B| reflecting the multi-phase structure of the gas.We also see that the energyweighted typical plasma  ≡  2  / 2  ≫ 1, as expected and observed in the ISM and CGM of typical galaxies (e.Where we see the mean temperature/thermal pressure of the medium drop sharply as the mass collapses into cold phases at ≲ 100 pc, we therefore see a sharp transition from a total-energy or volumeweighted  ≫ 1 to  ≪ 1 at smaller radii.The magnetic field evolution through this point is relatively smooth; it is the thermal phase structure which changes much more rapidly.We can also see that outside of the nuclear disk at ≳ 0.1 pc, the magnetic fields are order-of-magnitude isotropic (no single component strongly dominates |B|) and exhibit a large dispersion, reflecting the disordered and super-Alfvén turbulent morphology of the gas and multi-phase structure (there is at some radii a small bias towards radial B fields, reflecting inflows and outflows).Inside ≲ 0.1 pc where the ordered, thin, nuclear disk forms, we see this produces a much more ordered, predominantly toroidal field.This will be studied in more detail in Paper II, where we will examine the time-dependence of the field and its amplification mechanisms in detail.
At all radii, the kinetic energy density of gas is non-negligible (whether primarily ordered or disordered), as expected.The radiation energy density is always sub-dominant to kinetic, magnetic+thermal, and gravitational energy densities (∼   2  ) -this is expected at large radii where the galaxy is optically thin, but is surprising at the smallest radii, where again it will be discussed in greater detail in Paper II. However the radiation energy density we see in the simulation at small radii is expected (it is approximately that of a simple blackbody if we set the cooling luminosity equal to the accretion luminosity ∼  in  2  () at each ) -it is simply that the magnetic and kinetic energy densities are much larger.The relative composition of the radiation energy density is unsurprising: at large radii the broad NUV band dominates as expected for an optically thin young stellar population, whereas at small radii our adaptive "IR" (but really just any re-radiated light) band dominates when the medium becomes optically thick to NUV and optical emission.The cosmic ray energy density is also small (except perhaps at the very largest radii) compared to others in such a dense environment, as expected and discussed in more detail below.
Briefly, it is worth noting that for a reasonable estimate of the UV luminosity from un-resolved (≪ 100 au) scales around the SMBH given the accretion rate here, we might expect the broad line region (BLR) to reside at radii ∼ 20−150 light-days (Kaspi et al. 2005), or ∼ 0.01 − 0.1 pc.It is notable that the gravitational velocities at these radii are ≳ 1000 km s −1 , and (per Fig. 12) the disk covering fraction or / is relatively large ∼ 0.03 − 0.1 and increasing at these radii, where it is broadly comparable to the fraction of the UV/optical quasar continuum emitted in the broad lines (Vanden Berk et al. 2001;Richards et al. 2006).This is highly suggestive, but more quantitative comparisons (and conclusions related to the physical nature of the BLR "clouds" in these simulations) will require detailed post-processing radiative line transfer, which we hope to explore in future work.
As noted previously in § 3.3, the profiles here remain stable in time (well within their fairly large scatter) over the duration of the highest-resolution simulation, though they can, of course, evolve on larger radii (pre-refinement) on much longer timescales (of order many galaxy dynamical times).The time evolution of the innermost magnetic field structure, and its relation to amplification mechanisms, will be studied in Paper II.

Star Formation & Fragmentation Dynamics On
Different Scales We next turn to examining more dynamical properties of the simulation, in order to better understand what drives fragmentation and star formation (or the lack thereof) and inflows on various scales.

Definitions of "Disk" Dynamical Properties
Fig. 12 shows radial profiles (as Figs.7-10) for different dynamical properties of the simulation.We show the Toomre  parameter of the gas, defined as   ≈   /  Σ disk where   =   for the thermal  thermal , =  A for the magnetic  mag , =  turb 19 for the turbulent  turb , and  2 eff =  2  +  2 A +  2 turb for the total effective  eff .We also show how the different 19 More formally we follow Orr et al. (2019Orr et al. ( , 2021) ) and define  using the expressions for a multi-component disk from e.g.Romeo 1992, using the appropriate mass-weighted integrals over the distribution function (similar to defining ⟨   ⟩ −1 ∼ (Δ) −1 ∫ (/  )   ) to define the dispersion/sound speed in the gas (since the system is multi-phase).Note this thermal/magnetic/turbulent components contribute to the vertical support of the gas and the gas scale height, the sonic M  ≡  turb /  and Alfvénic M A ≡  turb / A Mach numbers, and the characteristic fragmentation scales of the disk determined by the characteristic maximum/dominant fragment mass ∼  Σ gas  2  (see Hopkins 2013b) and minimal Jeans mass ∼ (/6)  3   −3/2  −1/2 .Unless otherwise specified, these are mass-weighted averages in each radial annulus.

Fragmentation and Star Formation
In the CGM/IGM, we see the gas is thermally stable against self-gravity ( therm ≳ 1), the turbulence is trans-sonic (or sub-sonic in the hottest phases), and the gas is quasi-spherical ( ∼ ), all as expected.On galaxy scales from ∼ 1 pc to ∼ 10 kpc, the gas is not thermally stable, but has a  therm ≪ 1, so it fragments and forms stars, which in turn maintain supersonic turbulence (with turbulent dispersion dominating the effective scale-height, i.e.  ∼  turb, z /Ω) with an approximately constant (self-regulating) turbulent  turb ∼ 1, as observed in both nearby galaxies (Leroy et al. 2008) and highredshift starburst and quasar host systems (Forster Schreiber et al. 2009;Fisher et al. 2022;Reichardt Chu et al. 2022), and seen in previous simulations with similar physics (Hopkins et al. 2011(Hopkins et al. , 2023;;Orr et al. 2019Orr et al. , 2021)).The galaxy size (∼ kpc) and compactness are also reasonably similar to observed massive galaxies at these redshifts (compare Bezanson et al. 2009;Damjanov et al. 2009;van Dokkum et al. 2010;Hopkins et al. 2009cHopkins et al. , 2010)).The characteristic fragment mass at ∼ kpc scales is relatively large, ∼ 10 8  ⊙ , and this corresponds to the mass of the large star-forming cloud complexes or "clumps" seen in the gas morphology in Fig. 2 -these are more massive than Milky Way GMCs as expected because, for constant  turb , the clump mass scales as the gas fraction ∝  3 gas and this galaxy (with a gas fraction of ∼ 30% at ∼ 1 kpc) is a factor of ∼ 6 more gas-rich than the Milky Way (so we expect a maximal clump mass ∼ 200 times larger than the largest GMC complexes in the Milky Way), as studied in more detail for similar systems in Oklopčić et al. (2017).Given the large gas fractions and  turb , the system is still thick, with / ∼ 0.3 − 0.5 at these radii.All of these behaviors are consistent with many previous studies of the star-forming ISM in both idealized and cosmological galaxy simulations (Noguchi 1999;Bournaud et al. 2007;Ceverino et al. 2010;Hopkins et al. 2012bHopkins et al. , 2013a)), and again generically expected.As the ISM becomes thermally cold at ≪ 100 pc, we see the turbulence go from super-Alfvénic to trans or even mildly sub-Alfvénic, and highly super-sonic, and   contributes negligibly (even in a volume-weighted sense) to the vertical disk support.
It is important to stress that even though star-forming systems with thermal  therm ≪ 1 can and do self-regulate to turbulent (and magnetic)  turb ∼ 1, these are not locally stable against fragmentation and star formation (see references above).In fact, Hopkins & Christiansen (2013) show that such systems will always produce more fragmentation on small scales as  turb increases if  therm remains constant, owing to shocks and compressions generating locally overdense regions with  eff ≪ 1 (including the local turbulent, magnetic, gives significantly lower  (so is more conservative for our purposes hereessentially measuring  in the "most unstable" gas phases) than using e.g. a thermally-weighted average or a simple rms value of , which can be strongly biased by outflow motions or small amounts of gas in hot phases populating the "tails."Thermal Jeans Mass π c 3 s /6 G 3/2 ρ 1/2 Magnetic Jeans Mass π v 3 A /6 G 3/2 ρ 1/2 Total Mgas(< r) Fig. 12.-Top Left: Toomre  parameter of the gas (accounting for a multi-component potential) in annuli , accounting for the thermal ( therm ), magnetic ( mag ), turbulent ( turb ), or combined support of the gas.The ISM/CGM are semi-stable as expected, the galactic ISM is thermally unstable and with marginal turbulent stability, until ≲ 0.1 pc when stability leads to cessation of star formation.Top Right: Scale height / of the gas (mass-weighted), directly measured in each annulus as the median or rms |  | after rotating to the angular momentum axis in that annulus, compared to different velocities (thermal sound speed   , Alfvén   , or vertical turbulent   ) relative to   .At most scales, kinetic/turbulent support dominates, with trans-sonic thermal support in the CGM and magnetic support taking over at ≲ 0.01 pc.Bottom Left: Sonic (M  ≡  turb /  ) and Alfvénic (M  ≡  turb /  ) Mach numbers in each annulus (mass-weighted), for each component of the random motions.The velocity dispersions are broadly isotropic at most radii but inflow/outflow leads to a mild radial bias at ≳ kpc scales.The CGM/IGM are trans-sonic (sub-sonic in the diffuse gas, but the average here is dominated by dense substructure), galactic and smaller scales highly super-sonic; we see a clear trend of  decreasing at small  so the accretion disk is modestly sub/trans-Alfvénic.Bottom Right: Enclosed gas mass inside , versus characteristic maximum gravitational fragmentation (Hill/Toomre) mass ∼ Σ gas  2 , maximal turbulent fragmentation mass from Hopkins (2013b), and mass-weighted (so biased to much lower values versus volume or thermal-energy weighted) thermal or magnetic Jeans masses.At galactic radii, the gas-rich galaxy produces massive clump complexes with masses ≳ 10 9  ⊙ .Between ∼ 0.1 − 1 pc, the thermal Jeans mass approaches Σ gas  2 (equivalent to  thermal ≳ 1), and the magnetic Jeans mass exceeds  gas (<  ) (equivalent to all possible scales being magnetically sub-critical).and thermal energy densities).This is of course implicitly necessary for star formation to explain their turbulent selfregulation.

The Cessation of Star Formation at Small Radii
At smaller radii inside the BHROI, we see (1)  begins to rise for all components, owing to the steep rise in Ω, and the in particular the thermal  therm ≳ 1; (2) the disk begins to become thinner (as for relatively slowly-varying  A ,  c begins to rise); (3) correspondingly the turbulence becomes somewhat weaker (more sub-Alfvénic); (4) as the optical depth increases and gas becomes more thermally homogeneous at warm temperatures (Fig. 10) the minimum Jeans mass sta-bilizes20 while the characteristic upper fragmentation mass21 decreases into the stellar-mass range and becomes comparable to the minimum thermal Jeans mass (implying all scales are thermally stable); (5) the magnetic field becomes more ordered and dominated by a coherent toroidal component as the disk becomes more organized; (6) the "magnetic Jeans mass" becomes larger than the enclosed gas mass, so all scales are magnetically sub-critical.The combination of these effects simple mean or volume-weighted median is much larger), so we can see the behavior of the coldest, most fragmentation-prone phases at each radius.Note also that in a single-component, homogeneous slab disk with only thermal support, the ratio of thermal Jeans mass to Σ gas  2 scales as ∼  3/2 , but we see that this precise correspondence is broken for the more complex systems here.
21 This is defined as ∼ Σ gas  2 , dimensionally the same as the mostunstable mass or Toomre mass or maximum Hill mass in a disk, but with the O (1) pre-factor taken from Hopkins & Christiansen (2013) as the value above which the mass function of fragments becomes exponentially suppressed.leads to the cessation of star formation.
When  therm ≳ 1, the system is nominally locally "stable" in a formal sense.Still, since the cooling time is short compared to the free-fall time at these radii, we see  therm remains modest (not ≫ 10, for example) at all but the smallest radii ≲ 0.01 pc.As a result, one might expect an intermediate gravitoturbulent regime (e.g.Gammie 2001), and indeed the gas morphology in Fig. 2 appears consistent with this.While fragmentation in the gravitoturbulent regime is not "catastrophic" in the same sense as in the galactic regime (where  thermal ≪ 1 and gas locally fragments on its free-fall time), it could in principle still produce efficient fragmentation if we neglected some of the other effects above ( -Angular momentum direction of the gas (cos  ≡ ĵ ≡ ĵ • ĵinner ), relative to the mean angular momentum direction of the inner disk j inner (averaged within < 0.01 pc), as a function of BH-centric radius .We clearly see the aligned, dynamically cold disk (cos  ≈ 1, with relatively little variation) in the central ≲ 0.1 pc, with an un-aligned, and much more spherical/kinematically hot (broad distribution of cos ) at larger radii.
lence with these fields can create a dynamo or locally mix/reorder the field lines in a manner which further suppresses local collapse (Deng et al. 2020b;Riols et al. 2021).
A closely related and more formal statement of magnetic stability comes from our comparison of the magnetic Jeans mass23  B J ∝  3  / 3/2  1/2 to the total enclosed mass  gas (< ).While the magnetic Jeans mass in and of itself is not a strong determinant of star formation in the same way as the thermal Jeans mass, and the "magnetic  parameter"  mag ≫ 1 likewise does not alone formally ensure local stability (Lynden-Bell 1966), when  B J exceeds  gas (< ) (which occurs here at ≪ 1 pc), it is equivalent to the statement (for a homogeneous disk or spheroid) that any perturbation of any wavelength/size ≤  is magnetically sub-critical (has a massto-flux ratio sufficiently low that it cannot collapse), i.e. that fragmentation is strongly suppressed (Armitage 2015).
A second, less important but still non-negligible barrier to fragmentation at these radii is the strong torques which we see are producing angular momentum loss and inspiral on a timescale of order the orbital time (discussed in detail below, but this can be read directly off from Fig. 8, or inferred from Fig. 7 by simply noting  in ∼ 0.1 Σ gas  2 ).Akin to the problem of giant planet formation which was the focus of many historical gravito-turbulence experiments above (see e.g.Armitage 2015; Kratter & Lodato 2016), even if we neglect the magnetic fields, vigorous gravitoturbulence would lead to an initial collapse of a perturbation by a factor of a few in density, at which point (since the cooling time, while shorter than dynamical, is not completely negligible, and the clumps of order the most unstable wavelength in size are optically thick to their cooling radiation) its cooling would proceed more resisted by magnetic fields; however any such mode with a finite radial width/wavenumber -i.e. containing finite mass -will be strongly resisted by differential rotation/shear and rapidly sheared out.
23 We plot magnetic Jeans mass  B J rather than magnetic critical mass  Φ ≡  Φ Φ  / 1/2 =  Φ |B| / 1/2 because the latter can only be defined by reference to a specific area A or size, whereas  B J can be defined locally.But note that the dimensionless ratio  B J / for some volume enclosing mass  is simply proportional to (  Φ / ) 3 , so for our purposes the two can be treated the same as we only focus on their relative scaling.
slowly and it would contract quasi-adiabatically on a cooling time, but this must occur before inspiral.For planet formation one might have an inspiral time of millions of orbits in the gas-rich disk; here, one has only a few orbits.As a result, we see that with magnetic fields present so fragmentation is already suppressed, most of the mildly-overdense clumps that do form (e.g. one evident in Fig. 2) spiral inwards and are tidally sheared out upon reaching smaller radii (or simply accrete into our sink particle SMBH) before they can reach even order-of-magnitude overdensities (let alone become anywhere near dense enough to approach star formation).
It is worth noting that none of the above effects completely eliminate all fragmentation and star formation.We only see that occur on even smaller scales,  ≪ 0.01 pc, where the thermal Toomre  thermal rises extremely rapidly to values ≳ 1000 by  ≲ 0.001 pc (much more strongly suppressing gravitoturbulence).However, it is sufficient to ensure the star formation rate and total gas accretion rates onto stars are negligible compared to the gas inflow rates, and therefore that star formation (as well as stellar feedback, at least for the duration of this simulation at highest resolution) plays an essentially negligible role in the global dynamics of the system on ≪ pc scales.More detailed properties of the star formation at these innermost radii (including the IMF), exploring how the rare stars that do form are influenced by their environment (and how their feedback does or does not influence that environment locally) will be studied in Paper III.
We stress that this suppression of star formation is only possible because of the unique circum-SMBH environment, and is not a generic property of magnetically-dominated ( ≪ 1) media.The suppression of fragmentation by strong shear (Ω and  rising rapidly), and the increasing order of the toroidaldominated mean-field (which efficiently suppresses collapse in the radial and vertical directions) are directly related to the increasingly-strong differential rotation of the disk in an external (BH-dominated) potential, and as discussed in Paper III do not arise generically in environments like the circumquasar gas outside the BH radius of influence at radii ≫ pc, let alone more typical GMCs.

Different Contributions to the Torques
In Fig. 15, we now turn to understanding the torques driving gas inflows in more detail.We first simply plot the actual torques in the simulation.For every gas cell, we calculate the specific torque vector  ≡ r × a, where a is the acceleration from various sources (recording the value directly computed in-code), and we consider the component along the existing specific angular momentum direction j ≡ r × v (where r is defined as the vector distance to the SMBH), and plot this in units of  2  = |r|   () Ω(), so that a value of  • ĵ =   2  corresponds to the torque removing all of the angular momentum from an initially circular orbit in a time Δ ≈ ( Ω) −1 .We separately quantify this for the acceleration from MHD forces (the Riemann problem in the code), cosmic ray forces (using the full expressions from Hopkins et al. 2022d,f, which allow for both the tight-coupling and free-streaming limits), radiation forces (likewise allowing for both limits and anisotropic radiation tensors following Hopkins et al. 2020a), and gravitational forces.
Closely related to this, we quantify different components of the stress tensor in the code.The momentum equation solved in the simulation can be written:  (v)/ +∇ •  * = S, where  ) in shells as a function of distance .We restrict to cool gas with  < 10 4.5 K but the trends are qualitatively similar regardless.Shaded region shows the ∼ 90% inclusion interval.We plot the torques directly from the simulation from gravitational forces, radiation pressure forces, and MHD (magnetic, thermal, turbulent) forces.At all radii torques are efficient, implying angular momentum loss in a few orbital times.At large radii gravitational torques from stars on gas, plus MHD torques driven by stellar feedback, dominate.Where there are few stars, magnetic and turbulent/Reynolds torques take over.Top Right: Fractional (Fourier) mode amplitudes of asymmetric modes in the face-on projected gas surface density (Σ ( , ) ≡ Σ 0 (1 +    cos (  +  0,  ) )) in cylindrical annuli.At radii ≫ 0.1 pc there are order-unity asymmetries, often dominated by the lowest- modes (global asymmetries rather than just small-scale clumping).At small, increasingly Keplerian radii |   | decreases but only modestly, to ∼ 0.1 at ≲ 10 −3 pc.Bottom Left: Profile of different (volume-weighted) mean components of the magnetic stress tensor  mag ≡ (1/4  ) ( |B| 2 I − BB/2), in spherical coordinates, versus radius.We normalize to the mean value of the total stress tensor  at each radius.Solid (dotted) line correspond to  > 0 ( < 0).Bottom Right: Same, for the (total) kinetic stress  kin ≡  v v (note this is distinct from e.g. a Reynolds stress).
in the source term S includes e.g.non-hyperbolic terms from radiation and cosmic rays (see references above) and other terms relevant in the weak-coupling limit, and  * is the stress tensor which can be decomposed into: (1)  kin +  mag +  therm +  visc +  cr +  rad +  grav representing the sum of kinetic (including turbulent)  kin , magnetic  mag , thermal  therm , viscous  visc , cosmic ray  cr , radiation  rad stress tensors constituting the usual "total stress tensor"  =  internal , plus gravitational  grav forces.
In order to better understand the origin of the gravi-tational and kinetic stresses in the disk plane in particular, it is also helpful to quantify the degree of nonaxisymmetry of the system.Noting that the total surface mass density Σ tot (, ) within each cylindrical annulus  can be Fourier decomposed into Σ tot (, ) = ⟨Σ tot ()⟩ 1 + ∞ =1   () cos ( [ −  0,  ()]) , we extract the coefficients |  ()|, and plot the first few coefficients   as a measure of the global asymmetry.The behavior is similar for higher- modes, but  = 1 is most relevant for linear global gravitational instabilities interior to the BHROI; see Hopkins & Quataert 2010b, 2011b;Hopkins et al. 2009a,e.The first thing to note is that the torques are large, in a dimensionless sense, | • ĵ| ∼ 0.1  2  , i.e. the timescale for angular momentum loss of an initially-circular orbit is just a couple of orbital times ( orbit = 2/Ω).This is expected from the very large  in in Fig. 7 and shown in Fig. 8: we anticipate  in ∼  gas () | • ĵ|/(   ) ∼ (| • ĵ|/ 2  )  Σ gas    ∼ 10 − 100  ⊙ yr −1 at these radii (inserting typical values from Fig. 7 & Fig. 15 in the final evaluation).This means that accretion is fundamentally dynamical here, occurring on of order the dynamical time, as opposed to a slow, secular, viscoustype process as often assumed for much lower accretion-rate systems.
From a cursory examination of the components of  * extracted directly from the simulation, or from the torques in Fig. 15 (where various torques fall below the plotted range), or our discussion above of relevant physics and scalings, it is easy to confirm that the physical viscosity ( visc ), cosmic ray ( cr ), radiation ( rad ), and pure-thermal (isotropic by definition  therm ) terms in the stress tensor contribute negligibly to the torques at essentially all radii modeled here.There are really only three contributions of broad importance: the "gravitational torque" (r × g), and the MHD torques arising from a combination of magnetic ( mag ) and kinetic or Reynolds-like ( kin ) stresses.

The Gravitational & "Stellar Feedback" Torques
On scales ≳ pc, we see that gravitational torques are important (even if not always dominant) for the dynamics of angular momentum exchange, with large-amplitude |  | ∼ O (1) asymmetries (obvious in the visual morphology of gas and stars) producing strong torques.This is studied on these scales in much greater detail in Anglés-Alcázar et al. (2021), who show that such torques act most efficiently with gas forced into shocks and dissipation (allowing gas orbits to decay rapidly) via asymmetries in the stellar distribution which dominates the local mass density (hence the actual masses exerting the torques in Fig. 15).This does mean that even when gravitational torques are dominant on these scales, the MHD torque is generally order-of-magnitude comparable (the torques induce shocks which have comparable amplitude and [usually] opposite sign, as we see).As shown in Anglés-Alcázar et al. (2021), phenomena like the sign flip we see in Fig. 15 at ∼ 0.5 − 10 kpc are often transient and can flip back-and-forth, with the time-averaged effect of these torques on these scales being to reduce gas angular momentum.This also agrees with the results of previous "nuclear zoom-in" simulations (Levine et al. 2008;Prieto & Escala 2016;Prieto et al. 2017;Anglés-Alcázar et al. 2021) and both idealized simulations of small scales in gas+stellar nuclear disks (Hopkins & Quataert 2010b;Hopkins et al. 2016;Williamson et al. 2022), observations of nearby nuclear disks (Lauer et al. 2002;Hopkins & Quataert 2010a;Querejeta et al. 2016) as well as galaxy-scale simulations of galaxy mergers, strong bars, and large clump-type perturbations, which were the first to describe this gravitational torques process as uniquely efficient in mixed (collisional+collisionless) systems (Barnes & Hernquist 1991, 1996;Hopkins et al. 2009d,b).
On galactic scales ≳ kpc, we also see comparable and sometimes dominant MHD torques to gravitational torques, which relate to strong shocks sometimes (as noted above) driven by gravitational motions (e.g.infall/accretion, bar or clumpinduced shocks) but sometimes also due to e.g.strong shocks owing to stellar feedback events (motions that can be traced directly back to e.g.superbubbles and outflows).The latter, in particular, is directly related to the large scatter in the magnitude of the MHD torques at  ≫ pc -cells with || ≫  2  for example almost entirely owe to material being strongly accelerated by supernovae shocks and superbubbles/winds.24This is again consistent with previous studies and very similar to the results in e.g.Prieto & Escala (2016); Prieto et al. (2017).If feedback is "self-regulating" on these scales (e.g.on-average balances gravitational collapse and maintains a super-sonic turbulent  turb ∼ 1), then it is (by definition) true that the non-circular (or non-hydrostatic) motions generated by such feedback should be comparable to those generated by gravity (and generally have the opposite sign).This is closely related to the question of "what powers the turbulence" in the ISM (gravity or stellar feedback), where self-regulating models necessarily predict that if inflow is balanced by outflow and star formation (as it is here on super-kpc scales; see Figs. 7 &  8) the two should be comparable (Orr et al. 2020).
Note that this physical connection, as well as the large scatter, make it difficult to un-ambiguously determine precisely "how much inflow" is driven by each mechanism, as transient events like shocks could have large integrated effects on gas orbits.In future work, one can envision following the Lagrangian evolution of individual gas parcels which ultimately accrete onto the SMBH over cosmic time, and asking the question (distinct from what we plot here) of how they lost angular momentum (and whether they preferentially sample lower angular-momentum material in the first place) at different times and locations on their way to being accreted (akin to the studies on larger scales in e.g.Anglés-Alcázar et al. 2017b;Hafen et al. 2019Hafen et al. , 2020)).
Regardless of this, at smaller radii ≪ pc, with star formation efficiently shut down, we see a transition around ∼ 0.3 − 0.6 pc (Fig. 7) from the local mass density and gravitational field being dominated by stars at larger radii, to entirely gas-dominated at smaller radii.This dramatically reduces the efficacy of gravitational torques and (of course) any torques owing directly or indirectly to stellar feedback.As noted above and shown more formally in Hopkins & Quataert (2011b), the leading-order gravitational torque arises when one has a two-component system with a collisional, dissipative gas component being acted upon by a dominant collisionless component (e.g.stars).When the collisionless component becomes small, so the gas disk is effectively "one-component," the strength of the torques (averaged over the orbit of some gas parcel) drops dramatically, until it eventually is reduced to the small and higher-order resonant-only contributions given by Kalnajs (1971).The analytic prediction from Hopkins & Quataert (2011b) is that the torque drops ∝  * ≡ Σ * /(Σ * + Σ gas ) (in terms of the mean gas and stellar mass densities in an annulus of some radius  from the BH) when Σ * becomes smaller, which agrees well with the trend we see in the "transition region" (comparing Fig. 15 & Fig. 7).
Note that this does not mean that the  = 1 modes themselves cease; as shown in Fig. 15 they propagate towards small , albeit with decreasing amplitude.Indeed Hopkins & Quataert (2010b,a); Hopkins (2010) showed that this lopsided disk mode, if excited at large radius where the disk is marginally self-gravitating, can excite a response at all radii  → 0. However, the point in Hopkins & Quataert (2011b) is that for a given asymmetric mode amplitude | =1 |, the effective net torque on gas is weaker if  * is smaller.Moreover, as shown in Hopkins (2010), if the mass profile of the collisionless component ceases to rise sufficiently steeply towards  → 0 (e.g. for Σ * ∝  −  with  ≲ 0.5), then there is a refraction barrier (akin to an inner Lindblad resonance in many ways) across which the torque switches sign, exactly as we see in Fig. 15.25 If gravitational torques were the only mechanism for angular momentum transfer, then as speculated in Hopkins & Quataert (2011b), this barrier would create a trap and "pileup" of gas between ∼ 0.1 − 1 pc, until it became so dense it would necessarily fragment and form stars, until sufficient stars formed (assuming a large gas supply continued to flow in) to steepen the profile and reverse the barrier (making the stars-on-gas gravitational torque strong again), moving the barrier gradually inwards and building a steep stellar cusp.And this is what appears to happen in the simulations in Anglés-Alcázar et al. ( 2021), as discussed below.
However, here we see something at first apparently rather remarkable (though perhaps not on further reflection): at the radii where the gravitational torque becomes highly inefficient, the MHD torques "take over," with the same sign and broadly similar magnitude.

The MHD Torques
In Paper II, we will study the torques in the inner accretion disk in much greater detail, in order to understand the origins of the strong toroidal field, relative role of the Maxwell versus Reynolds torques, their dominant components and fluctuations in time and space, physical origin and ultimate energy sources, and how they are dynamically maintaining accretion.Here, we simply wish to identify and summarize some basic properties of the "MHD torques" across a broad range of radii.
As noted above, we see in Fig. 15 that at radii ∼ 1 − 1000 pc, the gravitational torques play an important role, and the MHD torques are heavily influenced by stellar feedback, consistent with previous work ( § 6 below).There are occasional, usually transient, exceptions, when e.g.strong shocks are induced in mergers and the shocks instantaneously dominate the torque -though as discussed above in such a situation the angular momentum exchange in the shock can be itself ultimately driven/determined by the gravitational forces (or stellar superbubbles and winds), so the "labeling" can be somewhat ambiguous (Hopkins & Quataert 2011b).
But to better understand the structure of these torques in either case, we plot the different components (in spherical coordinates centered on the BH) of the magnetic  mag and kinetic  kin stess tensors versus radius in Fig. 15.We normalize the components to the magnitude (Frobenius norm) of the sum internal stress tensor ∥ internal ∥ ≡ ∥ kin +  mag +  therm +  visc +  cr +  rad ∥ (i.e. the total stress ignoring the source terms so assuming the tightly-coupling limit for cosmic rays and radiation, and excluding gravitational forces, so representing the internal forces from the gas).At all radii, we see the kinetic components sum close to unity, i.e. represent a dominant term in the total.26Other than a small range of radii in the CGM, where the thermal pressure contribution to the stress is comparable to the kinetic (as we showed in Fig. 12, where the turbulence is trans-sonic), the other (non-kinetic, non-magnetic) terms in the stress are generally fractionally small.At large radii ≳ pc, the kinetic terms are quasi-isotropic (with a mild radial bias from inflow/outflow motion), and dominated by random motions (e.g. , with the mixed/off-diagonal terms ( kin   ,  kin   ,  kin   ) having essentially random (rapidly alternating) signs as they average out to smaller values than the diagonal terms.All of this is consistent with incoherent (e.g.turbulent and/or feedbackdominated) motions at a non-negligible fraction of the circular velocity -i.e.roughly as expected from the virial theorem -at similar Mach numbers for the components shown in Fig. 12.At ≪ 1 pc we clearly see the ordered disk form: the azimuthal (rotational) component    dominates the total stress, this component is itself strongly dominated by its mean/coherent component (   ≈ ⟨⟩⟨  ⟩ 2 ), and the radial and azimuthal terms become sub-dominant by a factor of ∼ 100 (implying turbulent/incoherent velocities more like ∼ 0.1   ).
We will analyze these terms to study e.g. the Reynolds stresses within the accretion disk in Paper II, but note that what is plotted here, for the sake of comparing the entire stress tensor and radial range, is not the Reynolds stress.Specifically, components like  kin   ≡ ⟨     ⟩ are defined as the average of the total values of the relevant velocity components like   , whereas the Reynolds stress is defined in terms of the incoherent components   =   − ⟨  ⟩.So the fact that, for example,  kin   < 0 here at all radii ≲ 1 pc simply means that there is, at this snapshot in time, net inflow through all radii ≲ pc in the disk, because it is dominated by its coherent components  kin   ∼ ⟨⟩⟨  ⟩⟨  ⟩ (and ⟨  ⟩ > 0 by definition of our coordinate convention for the inner, rotating disk, while ⟨  ⟩ < 0 denotes inflow).
For the magnetic stresses  mag , we see the corresponding expected behavior: overall ∥ mag ∥ is a fractionally small contribution to ∥ internal ∥ at large radii ≫ pc where magnetic field effects on the dynamics are small and  is increasingly large, and at radii ≳ pc the magnetic fields are quasi-isotropic/tangled (with again a mild radial bias), and magnitudes relative to velocity consistent with the Alfvén Mach numbers in Fig. 12.At sub-pc scales, we see the toroidal magnetic field becomes dominant and at ≲ 0.01 pc begins to contribute at up to an orderunity level to the total stress, though it is still sub-dominant to rotational support of the disk.This component is dominated by the coherent field, while the others remain dominated by largely incoherent fields (more detailed analysis in Paper II).
Here the   component is a more traditional Maxwell stress, though it is dominated by a mix of coherent and incoherent components depending on exactly which radius we analyze; we will study this in detail in Paper II but for our purposes here we can note (a) the sign is positive, which in the convention here means it is transporting angular momentum outwards and therefore promoting inflow; (b) the fractional magnitude (assuming the disk is orbiting at ∼   ) is comparable to the values needed to explain the fractional MHD torques and inflow rates in Figs. 15 & 7;and (c) at the smallest radii, the magnitude | mag   | is comparable to its kinetic counterpart even with the kinetic term defined in terms of the bulk/coherent components, meaning that the Maxwell stress must be at least comparable to the Reynolds stress (if not larger) at these radii.
Why do these torques appear to smoothly "take over" from the gravitational torques with broadly similar magnitude?This might at first to appear to require some sort of "conspiracy," but closer examination of Fig. 15 suggests a more mundane explanation, namely that this is required for continuity/steadystate.First, upon examination of Fig. 15 we see that it is not the case that there is some sort of "exact" boundary-condition matching occurring here.Both the gravitational torque and (especially) MHD torque have huge fluctuations in magnitude, so the "transition" between one and the other dominating is actually spread out, in a local sense (e.g.considering different narrow annuli in solid angle from the SMBH), over at least an order of magnitude in radius (again see Paper II for more details of the distribution of torques and stresses at these radii, in particular).The transition only appears "narrow" because we follow and plot so many orders of magnitude in radius.Second, this large scatter also means there are large fluctuations where one or the other dominates outside/inside this radius.Third, we see that even the instantaneous mass-weighed mean specific torque is very clearly not precisely constant as a function of radius, but fluctuates by an order of magnitude: the boundary between the mean torque being MHD-dominated and gravity-dominated is one such example (there is a factor ∼ 10 fluctuation in their sum between ∼ 0.1 − 10 pc).It is true that the nearest "peaks" in the mean specific torque on either side of ∼ 1 pc happen to have remarkably similar amplitude in this particular snapshot, but comparing other snapshots even these peaks are only comparable in an order-of-magnitude sense.With this in mind, it is much easier to understand.Continuity means that if there were a sudden change in the torque efficiency around ∼ 1 pc, mass would either "pile up" (or be evacuated), which would, for most reasonable models for the origin of the MHD torques (e.g.those discussed in Paper II in detail) and for the gravitational torque models lead to a corresponding increase (decrease) in the torques, until  in ∼ constant with radius was in approximate steady-state.While it is true in principle that a "sharp" discontinuity in the torques could be balanced (for the same  in ) by a similar discontinuity in the gas surface density, the fact that the transition is "smeared out" in both space and time by large local fluctuations and turbulence means that this cannot reasonably be self-sustaining (so Σ gas must be smooth, hence the specific torques being smooth).In summary, the "transition" being continuous is a statistical, order-of-magnitude statement over a fairly wide range of radii.

SIMULATION WITHOUT MAGNETIC FIELDS
In this section, we compare a simulation run without magnetic fields.We begin from the same initial condition/snapshot used for "zooming in" as our fiducial simulation, at the same time when we would normally begin our hyper-refinement pro-cess.But now we remove magnetic fields.This technically means the "initial condition" for the zoom-in is slightly out-ofequilibrium, but recall as shown above (1) it is already a highly non-equilibrium system, (2) on large scales outside of where the iterative hyper-refinement procedure begins, the magnetic fields are not as dynamically important ( is larger), and (3) we evolve each hierarchical level of the hyper-refinement several dynamical times before allowing a subsequent level of refinement so that each level can re-equilibrate, so this is given time to occur in these runs.
Altogether, this experiment appears to support all of the statements above regarding the role of magnetic fields.
Figs. 16 & 17 repeat some of our earlier comparisons in e.g.Fig. 2-3 and Figs.7-15, respectively, but for this simulation without magnetic fields.Unsurprisingly, on large scales (again, outside of where the iterative hyper-refinement procedure begins), the change is small.Indeed at all  ≫ 1 pc, the morphologies, star formation rates, gas and stellar densities,  = 1 mode amplitudes, scale-heights, thermochemical properties of the gas (phase distribution, temperatures, ionization states), and strength of gravitational torques are basically the same as in our "default" simulation with MHD.However we caution that this simulation is only run a short time, so this could simply arise at the largest scales because the system has no time to come to a new equilibrium -for galactic scales ≫ 100 pc, explicit studies of cosmological simulations with and without magnetic fields are more informative.But at sub-pc scales, we can immediately see some major qualitative differences appear.
Visually, we can see much stronger fragmentation setting in on scales ≲ 0.1 − 0.5 pc.This is also immediately evident in the surface density of star formation, which in the runs with MHD "cuts off" and is strongly suppressed at ≪ 1 pc, while without MHD it continues to rise monotonically to small radii.The integrated SFR (averaged over a few dynamical times) inside of < (0.1, 1, 10) pc rises from ∼ (0.1, 10, 30) M ⊙ yr −1 with MHD to ∼ (5, 120, 300) M ⊙ yr −1 without.The surface density of stars begins to rapidly rise on small scales as gas is depleted: whereas with MHD we saw gas dominate the local density over stars within  ≲ 0.1 pc, we now see stellar densities increase rapidly and beginning to dominate after a few tens of dynamical times at all radii  → 0.
Explaining this rapid enhancement of fragmentation, we see that without magnetic support, the disk / becomes smaller/thinner (without MHD) inside of  ≪ pc, by a factor of ∼ 10 − 30.There is still highly super-sonic turbulence supporting it, but it is well-established that absent magnetic fields, a disk supported by stronger super-sonic turbulence will actually have more rapid fragmentation (Hopkins 2012a).We see this manifest in much larger   especially for  ≫ 1 without MHD, as gravo-turbulent fragmentation runs away (boosting the negative gravitational torque at small radii).We also see that without magnetic fields, the hydrodynamic torques on sub-pc scales are much weaker than the magnetized "MHD torques" -moreover the sign of the hydrodynamic torques without MHD is actually opposite (they are net moving material outward).Thus while the stellar densities are still relatively low at  ≪ pc, this creates a "bottleneck" or "pileup" of material at these radii, which further assists the runaway fragmentation.We also show in Paper III that this leads to an even more top-heavy stellar IMF at these radii.The inefficient hydrodynamic torques without magnetic fields, coupled to efficient fragmentation, mean that the inflow rate to the BH Fig. 16.-Images of a re-run of our fiducial simulation without magnetic fields ( § 4).We show gas face-on (top) and edge-on (bottom) as Fig. 2-13.We overlay the star particles which form (blue points) to emphasize that the extreme "clumpiness" in the gas is a real effect: these are collapsing dense gas clouds which rapidly form stars and produce runaway star formation on sub-pc scales.A much smaller (spatially and in mass/surface density) inner non-star-forming disk remains, but it is truncated at the radii where the thermal-only Toomre  parameter falls below  ≪ 10 (∼ 0.01 pc; see Fig. 12).
is strongly suppressed at ≪ pc radii (especially at the smallest radii we follow, ∼ 10 −3 pc).The net inflow rate into our central resolution element -i.e. total mass growth rate of the sink interior to < 80 au is lower by a factor of ∼ 200 − 300 on average over the duration of the more limited no-MHD run.27So while still non-zero owing to some transient and non-spherically-homogeneous structure (and indeed still fairly large in an absolute sense), the accretion rates without MHD are dramatically suppressed relative to those with magnetic fields present, at least until a much larger stellar density is able to build up (beyond the duration of our simulation to explore).We will study the more detailed consequences for the disk at ≪ pc scales in this simulation in Paper II.

SCALES WHERE DIFFERENT SIMULATION PHYSICS "INGREDIENTS" BECOME IMPORTANT
From the above, we can refer back to Table 1 and review the major physical "ingredients" in these simulations, in order to discuss where each plays a crucial role in determining the dynamics of the system.We emphasize this because of course, certain physics will always be important by definition if one is interested in them for their own sake, or their influence on certain observations (e.g. one could have dynamically negligible magnetic fields but they would still be "important" to predict Zeeman observations).This is summarized in Table 2 and Fig. 18.

Gravity & Collisional vs. Collisionless Dynamics
The importance of gravitational dynamics is self-evident.At radii ≳ 0.01 pc, self-gravity is essential to follow the formation of the galaxy, inflows, feedback, and (especially crucial even in idealized simulations of a "patch" of this medium) fragmentation to form multi-phase ISM structure and stars.Without this, no meaningful predictions for inflow rates to the SMBH can be made, since this is the primary "competitor" with inflow to determine whether or not gas can actually reach the BH (not to mention how it qualitatively changes much of the dynamics).The presence of stars (and at larger radii, dark matter) also means one must be able to integrate collisional+collisionless systems simultaneously.
At smaller radii, where star formation has ceased and the potential is dominated by the SMBH, accurate gravitational orbit integration is obviously necessary: certain numerical methods for example cannot accurately integrate warped or precessing disks, or nearly-Keplerian cold disks, for many orbits before spurious numerical torques or "grid alignment effects" (in e.g.fixed-mesh codes or many smoothed-particle hydrodynamics methods) will destroy or artificially grid-align the disks (for extensive discussion of this and validation of the methods here in test problems, see Gaburov & Nitadori 2011;Hopkins 2015;Zhu & Li 2016;Deng et al. 2017Deng et al. , 2020aDeng et al. , 2021;;Deng & Ogilvie 2022;Hubber et al. 2018;Fletcher et al. 2019;Bonetti et al. 2020;Yamamoto et al. 2021;Franchini et al. 2022;Bortolas et al. 2022).Here our high-order Hermite integrator provides the ability to, for example, reasonably integrate a hard stellar binary for ≫ 10 5 orbital times in a strong tidal field (Grudić et al. 2021), much longer than necessary given the duration that we actually run our simulations to at their highest refinement level.But more importantly, we see that self-gravity is not negligible even at radii ∼ 1000  g .The spiral arms and  = 1 modes seen plainly in Fig. 2 and discussed above can play an important role in the dynamics even for a gaseous disk-to-BH mass  disk (< )/ BH ≪ 1.

Magnetic Fields
As expected based on most previous studies on galactic scales, at ≳ 100 pc magnetic fields play a relatively minimal role in the dynamics or gas thermodynamics (Kim & Ostriker 2015;Su et al. 2017Su et al. , 2018aSu et al. , 2019;;Hopkins et al. 2020b;Ji et al. 2020;Steinwandel et al. 2019Steinwandel et al. , 2022;;Martin-Alvarez et al. 2021;Ponnada et al. 2022;Whitworth et al. 2022).Even on scales from ∼ 1 pc to ∼ 100 pc, we see no evidence that the magnetic fields play a major role in the overall gas dynamics, and re-starting our simulation and re-running for ∼ 100 dynamical times on these scales without magnetic fields (without refining down to ≪ 1 pc scales) produces no major qualitative differences in our predictions for these scales, despite a plasma  ≪ 1.This is also expected, based on many previous studies of magnetic fields in the cold, neutral ISM on similar scales, where despite  ≪ 1 (because the gas is thermally cold), the magnetic pressure is still sub-dominant to other forms of pressure such as the "turbulent pressure" or (in the outer ISM) cosmic ray pressure (Federrath et al. 2014;Martin-Alvarez et al. 2018, 2022b;Guszejnov et al. 2020aGuszejnov et al. , 2022c;;Hopkins et al. 2020bHopkins et al. , 2022a;;Grudić et al. 2022b;Seta & Federrath 2022).Equivalently, we see above that the turbulence is still super-Alfvénic on these scales.In such situations the magnetic fields are largely "passively" tracing the local gas dynamics, rather than controlling it.There can of course still be indirect effects via smaller-scale dynamics (e.g.magnetic fields modifying the IMF which in turn modifies feedback; see Guszejnov et al. 2020b, 2022c andreferences therein).On smaller scales, however, we clearly see a reversal in this situation: the magnetic field strengths continue to grow, magnetic pressure dominates the vertical disk support and torques (discussed in greater detail in Paper II), the turbulence becomes trans or sub-Alfvénic, and so magnetic fields become essential to the dynamics.
The role of non-ideal effects is more minor.Though formally included, atomic & molecular viscosities and conductivities are everywhere negligible compared to numerical diffusion (and other physical processes), as expected.Anisotropic Braginskii viscosity and conductivity, given their strong temperature dependence, are only expected to be important in the most diffuse, hot phases of the CGM and ISM, and even there have relatively small effects (Su et al. 2017), so while included here we do not expect it to change any of our conclusions if they were excluded (and we see the viscous stress tensor is almost always relatively small).
In highly-neutral gas, we have re-run our simulation briefly turning on and off Ohmic resistivity, the Hall effect, and am-bipolar diffusion in turn28 to examine their relative importance.A conservative estimate of this comes from comparing the relevant timescale or timestep ∝  2 /  for process  on scale  to the other code timescales/timesteps (for e.g.other diffusion processes, sound-crossing, Alfvén-wave crossing, etc.).We see that Ohmic resistivity is never dominant given the density and ionization fractions we resolve.Ambipolar diffusion can be the most important effect of these three non-ideal terms in the least-dense but still overwhelmingly-neutral phases of gas (e.g.ISM-like molecular phases), but we see that including or excluding it has almost no effect on the global dynamics in the simulation, because even in the regions where it dominates over Hall and Ohmic terms, the ambipolar diffusion time is almost always much longer (often by orders of magnitude) than other transport process timescales such as the turbulent dissipation/reconnection timescale (O (/ turb )), as seen in most modern idealized simulations of protostellar core collapse (Chen & Ostriker 2014;Wurster et al. 2021) and studies of individual GMCs with realistic star formation and feedback (Mac Low & Klessen 2004;Vázquez-Semadeni et al. 2011;Sadanari et al. 2022).The Hall term is most potentially interesting: we do see a regime where the Hall term is dominant among non-ideal MHD effects and where the relevant timescale is shorter than other resolved timescales.This specifically occurs in extremely-dense gas forming protostellar disks at our highest resolution level (∼ 10 au distance from sink particles) in the star-forming disk (e.g.mostly at ≳ 1 pc) where the local densities are ≫ 10 12 cm −3 (about a million times higher than the average at those radii; Fig. 7) and the ion fractions can (locally) become extremely small (typically   ≡  ion / neutral ∼ 10 −17 ).This is not surprising: indeed, studies have shown that Hall effects can be important for the dynamics of proto-stellar disks and planetary disks on these spatial and ionization fraction scales (Bai & Stone 2017;Zhao et al. 2020;Lee et al. 2021;Tsukamoto et al. 2017).However, this is clearly not directly important for the global, quasar accretion disk-scale dynamics (the fraction of the total gas mass or volume in such protostellar disks is, as we showed above, negligibly small at these radii).Though the Hall effect could in principle indirectly alter e.g. the IMF of stars if it helps regulate accretion through individual protostellar disks onto the protostars themselves via Hall MRI, it would require much higher resolution (compared to our simulation) within the individual protostellar disks to resolve this (see Paper III).Crucially, within the global quasar accretion disk on sub-pc scales, even though the densities are very large, the warmer temperatures mean that the ionization fractions are vastly larger, ∼ 0.01 as shown in Fig. 10.This means that the characteristic timescales for Hall MHD effects within the quasar accretion disk and ISM as a whole are typically ∼ 11 − 15 orders of magnitude longer than the disk dynamical time at the radii we model here, so can be safely neglected.29Fig. 18.-Cartoon illustrating the hierarchy of scales as Fig. 5, with a heuristic description of the process driving fastest angular momentum loss on each scale.We show the image from our simulation, size scale, and descriptor as Fig. 5, together with a list of characteristic processes that drive angular momentum loss on these scales in a "slow" or "secular" fashion (timescales much longer than the dynamical time) or in a "fast" or "dynamical" fashion (timescales of order dynamical times).See discussion in § 5-6.Illustrations of numerical simulations for each "fast" scale/process are shown, taken from simulations first presented in Hopkins & Quataert (2010b); Hopkins et al. (2014); Torrey et al. (2017) for ≳ pc scales and from the simulations here on sub-pc scales.

Cosmic Rays
We see fairly minor effects turning on/off explicit cosmic ray transport, or switching between the simpler sub-grid model from Hopkins et al. (2022b) and the more detailed and physically-derived explicit cosmic ray dynamics models developed in Hopkins et al. (2022d); Hopkins (2022); Hopkins et al. (2022f,g), or even simply assuming a uniform cosmic ray background for purposes of ionization rate calculations.This is expected given the detailed studies of cosmic ray dynamics in e.g.Su et al. (2019Su et al. ( , 2020Su et al. ( , 2021)); Chan et al. (2019Chan et al. ( , 2022)); Hopkins et al. (2020bHopkins et al. ( , 2021b,c,d);,c,d);Ji et al. (2020Ji et al. ( , 2021)); Buck et al. (2020); Peschken et al. (2022);Martin-Alvarez et al. (2022a) as well as observational constraints from starburst galaxies (Lacki et al. 2011;Griffin et al. 2016;Zhang et al. 2019;Heesen 2021), which show that for starburst systems and massive high-redshift galaxies the cosmic ray energy is lost via catastrophic and Coulomb+ionization interactions on a timescale short compared to other timescales of interest (i.e. the galaxies are approximate proton calorimeters).It is precisely the opposite regime: tenuous CGM/IGM gas around review), we therefore neglect the Hall term in our default simulation and only include it in these tests run for a shorter time period.low-redshift dwarf and ∼  * galaxies, where CRs are seen in the studies above to have the largest effects (where they are observed to escape from galaxies into the CGM efficiently; see Lacki et al. 2011;Rojas-Bravo & Araya 2016;Lopez et al. 2018;Persic & Rephaeli 2022;Butsky et al. 2022).This is not surprising, based both on the more detailed studies above, but also simple analytic considerations.If we assume CRs are injected in the midplane and treat the gas as a uniform slab with tangled magnetic fields, then assuming a typical scattering rate similar to the constraints from the Solar system applies (see Hopkins et al. 2022f, for details), then given the density profiles in Fig. 7 at initial injection radii  inj ≲ 10 kpc catastrophic losses would remove all of the proton energy before propagation to a distance ≲ 0.3  inj .If we further assume the injection is proportional to the SNe rate assuming a steady-state SNe rate proportional to the star formation rate, itself scaling as some efficiency  SF per freefall time, this would produce a steady-state CR energy density in the ISM (again, taking the profiles from Fig. 7) of ∼ 30 eV cm −3 ( SF /0.01) (kpc/ inj ) -reasonably similar to what we measure in the simulation and much less than any of the dominant energy densities that we plot at any radii  ≲ 100 kpc in Fig. 10.And this is essentially an upper limit to the CR energy density, as it ignores other loss terms from trapping, denser sub-structure, advection or streaming.

Radiation Transport & Thermo-Chemistry
Clearly some cooling physics is important on all scales we simulate, but again the nature of that cooling and the role of radiation changes with scale.On large scales ≳ 10 pc, the cooling can be well-approximated as optically-thin (the usual approximation in galaxy formation simulations), although the actual chemistry can be enormously complex (as e.g. the models here account for a huge range of atomic and molecular and ionization and dust and other processes, with a multi-band radiation background, interactions with cosmicrays, non-equilibrium photo-chemistry, etc.).On these scales radiation is important for determining self-consistent ionization and photo-heating and radiation pressure dynamics from stars within the galaxy but its global effects can be captured reasonably well by simple approximations such as the LE-BRON method (see Hopkins et al. 2020a), and the gas cooling radiation itself can be largely neglected in the dynamics.
On the smallest scales we resolve, cooling is still important -in fact the disk has a cooling time short compared to its dynamical time even at the smallest scales we resolve (discussed in more detail in Paper II), so one cannot simply approximate the disk as strictly adiabatic.But the chemistry becomes substantially less complex as dust is sublimated, molecules dissociated, and in general the system becomes more and more locally black-body-like (and eventually at sufficiently small radii in the accretion disk the medium will be largely ionized, with chemistry relevant for second-order [though still potentially important] effects like metal line absorption, see e.g.Proga 2007;Jiang et al. 2016 and references therein).In this regime radiation is dominated by the cooling radiation itself, although it can increasingly be approximated via simple blackbody or gray-body approximations and the photon mean-free paths become short, so methods like flux-limited diffusion or other even simpler analytic radiation treatments may be valid (as in e.g.Thompson et al. 2005;Rafikov 2007;Derdzinski & Mayer 2022), and one could even approximate the disk as e.g.locally isothermal (or following some effective adiabatic index, with a mean temperature that depends on the distance from the BH, as is common in e.g.shearing-box simulations).
The complexity is maximized "in between," here from radii ∼ 0.01 − 10 pc (or more generally given how this should depend on the opacities, from surface densities Σ gas ∼ 10 4 − 10 7 M ⊙ pc −2 ).Here the system is optically thick to its cooling radiation, but not so optically thick that one can treat the radiation/dust/gas temperatures as in strict local thermodynamic equilibrium (LTE); more complex species such as dust, atoms, and molecules are all present in varying abundances (also not in LTE); and processes such as H − Kramers opacity -often neglected in both galaxy-scale simulations and accretion disk simulations -can dominate the opacities (this is dominant over a significant fraction of the dynamic range at ≪ 1 pc, given the high atomic abundance and relatively high free electron fraction producing significant H − ).Notably, the opacities would be orders-of-magnitude incorrect if we simply ignored dust destruction or neglected H − or detailed gas-phase opacities, and the temperature and cooling rates would also be order-of-magnitude incorrect if we simply assumed LTE and that all (radiation/dust/gas) temperatures were in equilibrium.30Above, we discuss how this produces 30 Notably, Derdzinski & Mayer (2022) do point out the importance of these some substantial differences (most notably, shutting down star formation), compared to previous studies which treated the transition region more simply.

Star Formation (Sink Particle and Unresolved) &
Stellar Feedback Star formation is clearly important on scales ∼ 0.1 − 10 4 pc.On much larger CGM/IGM scales we do not expect star formation to occur given the low densities; on much smaller scales we self-consistently see it suppressed so it can be again neglected.On scales ≫ pc, we see that the characteristic "fragmentation mass" ∼ Σ gas  2 gas ∼ Σ gas (/Ω) 2 expected in any turbulent fragmentation cascade (Hopkins 2012a) (akin to the "Toomre mass" for a marginally-stable disk) is ≫ 10 4 M ⊙ , so the stellar IMF should be well-sampled by the typical fragmenting clouds containing most of the mass and most of the star formation (Evans 1999;Vázquez-Semadeni et al. 2003;Hopkins 2012bHopkins , 2013a;;Grudić & Hopkins 2019).This means on these scales, approximating the star formation via "galactictype" models wherein one simply seeks to identify fragmenting sub-regions and, from them, statistically samples some IMF, is reasonable.Of course, one could always imagine "zooming in" to some sub-region (e.g. an individual GMC) on these large scales and studying it with resolved star formation models as in STARFORGE, but this would be more akin to standard studies of star formation in isolated clouds in the typical ISM (Guszejnov et al. 2021;Grudić et al. 2021;Guszejnov et al. 2022b,c,a) and is not strictly necessary for recovery of the large-scale dynamics (though it could of course be indirectly important via calibration of the "correct" IMF to use and other related properties of the stars themselves).Our preliminary analysis in Paper III of the resolved IMF at smaller radii also supports the use of a statically-sampled universal IMF on these larger scales.
"Resolved" star formation physics is therefore strictly necessary over a relatively narrow range of intermediate radii, here primarily from ∼ 0.1 − 1 pc.Still it plays a crucial role in the simulation here of actually allowing us to validate that SF should indeed cease at ≪ 0.1 pc.The galactic-type models cannot truly self-consistently predict this, since the SFR for a small "patch" of the ISM with certain properties is assumed, not self-consistently resolved.On these scales the characteristic upper limit of the fragmentation mass is still relatively large, ≳ 100  ⊙ , but not so massive that a well-sampled IMF can be assumed.But the thermal conditions of the gas and its increasing warmth and magnetic support mean that the sites of individual star formation are highly constrained, as we discuss above and in more detail in Paper III.

COMPARISON TO PREVIOUS RESULTS
6.1.Galactic Scales (≳ 100 pc) As noted in § 5, on relatively large scales ≳ 100 pc, which one could reasonably call "galactic," our results are opacity terms in quasar accretion disks at broadly similar radii, though the parameter space they explore is rather distinct from that here, and they adopt a simpler analytic disk model with opacity fitting functions calibrated for proto-stellar disks in Lin & Papaloizou (1985); Bell & Lin (1994).However direct comparison with e.g.their Table 1 or Figures 1-3 shows that the explicit chemical network and non-equilibrium dynamics evolved in our simulationimportant for capturing conditions in the sub-pc AGN disk that are not analogous to protoplanetary disks at a similar density and temperature (e.g.much shorter dynamical times, stronger radiation and cosmic ray fields, sublimation of dust grains, higher accretion rates) -can produce order-of-magnitude quantitative differences in the detailed opacities.broadly consistent with previous FIRE studies of massive, high-redshift galaxies.But it is worth reiterating some basic conclusions (some of which are summarized in Fig. 18): (1) the galaxies are very much not in steady-state or equilibrium, with large clumps, mergers, and feedback-driven perturbations to the potential of order unity (Ma et al. 2015(Ma et al. , 2018a;;Oklopčić et al. 2017;Price et al. 2017); (2) they feature prominent "cold flows" in the halo contributing substantial "cold mode" accretion onto the galaxy (Feldmann et al. 2016;Faucher-Giguère et al. 2015, 2016;Sravan et al. 2016;Chen et al. 2020); (3) "gravitational torques" play a key role in the dynamics of angular momentum exchange, with the gas predominantly being forced into shocks and dissipation by asymmetries in the stars (Anglés-Alcázar et al. 2017c;Ma et al. 2021;Trapp et al. 2022); (4) the galactic ISM is highly multi-phase and unstable, with relatively short-lived structures (Orr et al. 2017;Kim et al. 2018a;Ma et al. 2020b;Smith et al. 2019); (5) the plasma  ≫ 1 except in the cold-phase ISM, where magnetic pressure is larger than thermal but still order-of-magnitude sub-dominant to turbulent energy densities (so magnetic fields are not dynamically dominant; Su et al. 2017Su et al. , 2018aSu et al. , 2019;;Guszejnov et al. 2017Guszejnov et al. , 2019Guszejnov et al. , 2020a;;Hopkins et al. 2020b); (6) stellar feedback rapidly becomes less efficient above a critical acceleration scale ∼ ⟨ * / * ⟩ ∼ 10 −8 cm 2 s −1 (corresponding to a total-enclosed-mass effective surface density  enc /  2 ≳ 10 3 M ⊙ pc −2 ; Grudić et al. 2018Grudić et al. , 2019Grudić et al. , 2022a;;Ma et al. 2020a;Shi et al. 2021;Byrne et al. 2023b); (7) star formation is rapid but inflows are dynamical, with  ∼  gas Ω, and co-exist with outflows (as there is no spherical symmetry) so can still out-pace star formation (Sparre et al. 2017;Orr et al. 2021;Flores Velázquez et al. 2021;Ma et al. 2018b).
There are novel advantages of the simulation here.For one, it includes a variety of physics not included in all previous FIRE studies: e.g.magnetic fields with non-ideal MHD, detailed thermochemical treatments of non-equilibrium chemistry and opacities for the highly optically-thick and/or dustfree regimes, explicit M1 radiation hydrodynamics (as compared to simpler RHD treatments).For another, it reaches significantly higher resolution than some previous FIRE studies cited above, with mass resolution ∼ 10 3 − 10 4 M ⊙ throughout the galaxy.This allows us to confirm that, at least in the simulation run up to the time of "hyper-refinement" in the nucleus, these additional physics and numerics improvements do not appear to have a major qualitative effect on any of the conclusions of those previous papers regarding global galaxy properties.The weak effects of these physics on gross properties at large scales had been noted before (Hopkins et al. 2018b(Hopkins et al. , 2020a(Hopkins et al. ,b, 2022a;;Su et al. 2017Su et al. , 2019;;Wheeler et al. 2019), but those studies focused on lower-mass systems, so we extend them here.
However, the simulation here also has a serious and obvious disadvantage compared to previous studies: we only simulate one case study, and once we turn on hyper-refinement, only simulate it for a very short cosmic time.So studies of galaxy-scale properties are generally better-served by dedicated simulations without hyper-refinement (using e.g.subgrid models for BH accretion and feedback) that can evolve for longer times.The primary purpose of our including these scales in our simulation here is to generate self-consistent initial and boundary conditions for smaller scales, and of course to inform and refine such sub-grid models for future studies.

Galactic Nuclei Scales (∼ 1 − 100 pc)
On scales ∼ 1 − 100 pc, our conclusions are largely similar to those in the dedicated hyper-refinement experiments in Anglés-Alcázar et al. (2021) (themselves in key respects similar to previous studies such as Levine et al. 2008Levine et al. , 2010;;Wada et al. 2009;Hopkins & Quataert 2010b;Hopkins et al. 2016;Torrey et al. 2017;Prieto & Escala 2016;Prieto et al. 2017 or other more idealized nuclear simulations in e.g.Emsellem et al. 2015;Beckmann et al. 2019;Sivasankaran et al. 2022;see Anglés-Alcázar et al. 2021 § 5.1 for a summary).For example (see also Fig. 18): (1) gravitational torques between gas and stars (largely stars at similar radii to the gas) again dominate the accretion physics (even more strongly than on galactic scales); (2) angular momentum support is the key "barrier" to inflows and accretion; the accretion is qualitatively distinct from a radial or Bondi or turbulent accretion problem, and application of Bondi-Hoyle-Lyttleton-type accretion-rate estimators based on the gas properties (as opposed to those which account for effects like gravitational torques, star formation, and stellar feedback; see Wada et al. 2009;Hopkins & Quataert 2011b;Hopkins et al. 2022e) at these scales gives an accretion rate which is typically incorrect by ∼ 4 − 8 orders of magnitude; (3) the ISM is highly multi-phase and unstable and rapidly star-forming, with most of the gas mass at these radii in cold/warm neutral phases, but (4) accretion is dynamical owing to said gravitational torques.
As noted in § 3.3, stellar-feedback driven wind/outflow rates decline interior to ≲ 100 pc, We note that this is evident even before our refinement begins (so is a largely resolutionindependent statement) and appears both in our simulations with and without magnetic fields.The same effect is seen in previous zoom-in simulations (Levine et al. 2008;Anglés-Alcázar et al. 2021) as well as more idealized simulations of galactic nuclei (Wada et al. 2009;Hopkins et al. 2016 and dense molecular clouds (Geen et al. 2017;Kim et al. 2018b;Grudić et al. 2018).As predicted analytically Fall et al. (2010); Grudić et al. (2019Grudić et al. ( , 2020)), this occurs when the gravitational acceleration scale |a| =  2  / exceeds the IMF-averaged momentum flux per unit mass of a young stellar population (∼ 10 −7 cm s −2 ), so stellar feedback becomes inefficient at wind-launching.
The much more detailed physics (MHD both ideal and nonideal, multi-band explicit RHD, expanded opacities, individual star formation/evolution) and higher resolution (∼ 3 ordersof-magnitude improved) here do lead to some differences of potential importance for observational diagnostics on these scales, but do not change the key qualitative physics of accretion, outflows, and star formation above.Given the cold ISM prominence, we see  ≪ 1, but magnetic fields do not dominate the torques and are sub-dominant to turbulent and "gravitational" pressure and do not strongly alter the dynamics on these scales (much like in GMCs and HI filaments in the "normal" ISM of  ∼ 0 galaxies; see Su et al. 2017;Martin-Alvarez et al. 2018, 2021, 2022b;Guszejnov et al. 2020a;Benincasa et al. 2020).The inner galaxy begins to become optically thick to its own cooling radiation as we reach Σ gas ≳ 10 3 − 10 4 M ⊙ inside of ∼ 10 − 100 pc, and this does modify the phase structure: there is a large amount of warm molecular gas at ∼ 10 3 K, the dust temperature rises to ∼ 100 K (well above the CMB temperature at this redshift), and the ISRF becomes strongly dominated by the re-radiated/IR radiation.These are notably similar to observed properties of gas at similar densities in the nuclei of local starburst galaxies such as Arp 220, NGC 6240 and others (Tacconi et al. 1999;Lonsdale et al. 2003;Evans et al. 2006;Iono et al. 2007;Mason et al. 2006;Greve et al. 2009;Ott et al. 2011;Scoville et al. 2017) as well as inferences in high-redshift quasar hosts (Casey et al. 2009;Riechers et al. 2009;Younger et al. 2009;Wang et al. 2008;Izumi et al. 2016;Lelli et al. 2022).In future work, we will investigate their consequences for observables as well as whether they have any impact on the stellar IMF, but for now, they do not appear to qualitatively change the most important conclusions from Anglés-Alcázar et al. ( 2021) regarding accretion.
Given this, the primary purpose of our extended physics and resolution on these scales is to (1) test and validate the conclusions of Anglés-Alcázar et al. ( 2021) with more detailed simulations accounting for a range of physics neglected therein; (2) make more accurate predictions for observables and future sub-grid models on these scales; (3) enable exploration of more detailed quantities like the IMF; and (4) to provide self-consistent initial and boundary conditions for even smaller scales.
6.3.Approaching The Accretion Disk (∼ 0.01 − 1 pc) On scales ≪ 1 pc, however, we see significant deviations from the behavior seen in Anglés-Alcázar et al. ( 2021) and other similar studies described above (and see also Kawakatu & Wada 2008;Hopkins & Quataert 2010a,b;Hopkins et al. 2016;Wada et al. 2009;Schartmann et al. 2010;Hobbs et al. 2011;Izumi et al. 2016;Williamson et al. 2020b;Kawakatu et al. 2020).There are several closely-related key qualitative differences.Perhaps most importantly, we see star formation shut down, as the magnetic  mag ≫ 1 and thermal  thermal ≳ 1 as the system becomes more optically-thick and magneticallydominated.This obviously involves both the MHD and RHD physics (coupled explicitly to the thermo-chemistry) here, as well as a resolved individual star model which can detect and mass-resolve individual stellar-mass patches that might be collapsing (or not).In contrast, in simulations like Anglés-Alcázar et al. (2021) and almost all the other simulation examples in § 6.2 above, a simple "galaxy-scale" sub-grid star formation prescription was adopted at all scales and magnetic fields were neglected, which meant star formation could not "cease" in this manner, and as a result, continued to be efficient (even growing in efficiency at smaller and smaller scales) at all resolved scales therein.
This immediately produces important consequences.With star formation suppressed, we see a transition at ≲ 0.5 pc where the local mass density becomes gas-dominated, and gravitational torques (whose efficiency depends strongly on there being a dominant collisionless component of the local mass density to drive shocks in the gas) become less efficient (and even reverse sign).However, Maxwell and Reynolds torques from the strongly-magnetized, gravito-turbulent disk take over and continue efficient inflow (Fig. 18).The detailed structure of the disk, its turbulence and magnetic fields, their origins, and how they drive accretion, will be the subject of detailed study in Paper II, but depend directly on magnetic fields.Strong  = 1 modes persist and a lopsided disk with clear spiral structure forms, and some star formation does occur, which will be studied in Paper III, but the SFR inside these radii is small compared to inflow rates.
As noted above, there is a body of work with overlapping physics and results here in historical idealized simulations of nuclear "torus" scales around AGN, such as those in Kawakatu & Wada (2008) Kawakatu et al. (2020).These studies generally were more akin to Anglés-Alcázar et al. (2021) in that they included a more limited range of physics (often, but not always, neglecting MHD and RHD, and in all cases using a much simpler thermo-chemical network compared to that here).But most crucially, these were like Anglés-Alcázar et al. (2021) in using idealized star formation prescriptions with some statisticallyaveraged star formation rate per free-fall time above some density put in "by hand" as opposed to explicitly resolving individual stars and star formation.As such, their conclusions, like Anglés-Alcázar et al. (2021) as summarized above, have a great deal in common with ours on larger scales but diverge from ours when star formation shuts down.
There has however been some work specifically using simulations designed to resolve individual star formation to predict e.g. the IMF of stars forming in circum-nuclear disks (Nayakshin & Sunyaev 2005;Nayakshin et al. 2007;Klessen et al. 2007;Hopkins 2013c;Bonnell & Rice 2008;Alexander et al. 2008;Hobbs & Nayakshin 2009;Frazer & Heitsch 2019), to which we will compare in more detail in Paper III.These studies have found some similar conclusions to those here, e.g. that coherent lopsided modes in disks are ubiquitous, as expected from analytic considerations as discussed above (and other conclusions specific to the IMF, like its being somewhat top-heavy, see Paper III).But again these usually neglected physics such as magnetic fields and self-consistent radiation hydrodynamics tied to the thermo-chemistry of molecular/atomic/neutral phases -all crucial, as we argued above, to follow the self-consistent suppression of star formation and transition in structure on these scales.Even more importantly, these simulations in the past have been much more limited in the dynamic range of scales around the SMBH which they could probe, so needed to adopt somewhat ad-hoc initial and outer boundary conditions at ∼ pc, and therefore could not self-consistently predict the transition between star-forming and accretion disk.Of equal importance, all of those IMF studies referenced above explored parameter space orders-ofmagnitude distinct from that here, with much lower gas masses and densities (in most cases because they were designed to specifically understand the sub-pc stellar disk around Sgr A * , rather than the most luminous quasar environments like our study here).
Alternatively, some recent studies have extended accretion disk models "outwards" to these scales (Namekata & Umemura 2016;Chen et al. 2023) to explore fragmentation.But again, these simulations necessarily focused on a small dynamic range with specific initial/boundary conditions and physics, so were not attempting to link different scales in the same way as we do here.
As such, we stress that these different types of intermediatescale simulations are highly complementary to studies like those here.Our hope is that a study like this provides additional motivation and improved understanding of the necessary initial and boundary conditions and choice of "physics included" in these sorts of idealized, more-restricted-in-scale nuclear simulations, in the future.6.4.Within the Accretion Disk (≲ 0.01 pc) By the time we get to the more traditional "accretion disk" scales at ≲ 0.01 pc, star formation is inefficient, so the dominant physical ingredients are broadly similar to those traditionally invoked in AGN accretion disk simulations: ideal MHD and radiation-hydrodynamics in a nearly-Keplerian potential.
Still, we see a several key qualitative differences between our predictions and the assumptions of the vast majority of existing quasar accretion disk simulation literature (although some recent work shows striking similarity, see e.g.Kudoh et al. 2020), which will be studied in detail in Paper II so we only briefly review them here.Most of these have to do with the initial/boundary conditions.A strong  = 1 coherent eccentric gas disk mode persists, induced (and propagating inwards) by the asymmetry from large radii (Hopkins 2010).Self-gravity is not completely negligible and some non-zero star formation persists (Paper III), along with some gravitoturbulence (Paper II).The disk is still predominantly neutral even at these outer radii (though this will change when the disk temperature rises to ≫ 10 4 K at smaller radii) and can still efficiently cool, with a cooling time still less than its dynamical time, so the gas cannot be treated as adiabatic, and the opacities include important contributions from mostly-neutral gas contributors like H − , usually ignored in accretion disk simulations (which generally, if they follow explicit RHD, assume just some combination of free-free/Compton and metal line opacities).As a result, the turbulence is vigorous: trans-Alfvénic and highly super-sonic.And perhaps most importantly, the disk is strongly magnetized as a result of flux-freezing from the magnetic flux being fed to it from the ISM (as detailed in Paper II), sustaining a "flux-frozen" or "flux-fed" disk with plasma  ≪ 1, even in the midplane.Again, the hope is that the simulations here will provide additional motivation for new generations of accretion-disk simulations exploring these rather distinct portions of parameter space.
7. CONCLUSIONS We present novel simulations which utilize the galaxy-scale cosmological physics of the FIRE simulations to inform the small-scale physics of individual star formation and stellar evolution of STARFORGE.This allows us to run a cosmological simulation employing a super-Lagrangian refinement technique to reach ∼ 10 −4 pc resolution in a ∼ (100 cMpc) 3 box, i.e. the equivalent of a ≳ (10 12 ) 3 uniform-resolution simulation.More importantly, we incorporate a range of physics including (non-ideal and kinetic) magneto-hydrodynamics; self-gravity, star formation, stellar evolution, and (proto)stellar feedback (including jets, main-sequence mass-loss, multiband radiation, core-collapse and Ia supernovae); explicit multi-band radiation-MHD (with separately evolved dust, gas, and radiation field temperatures/bands); and detailed thermochemistry accounting for a huge range of processes and opacities including dust-gas coupling, sublimation, non-equilibrium atomic and molecular chemistry, metal lines, H − and others directly coupled to the RMHD solver.This allows us, for the first time, to self-consistently treat both the limits of "traditional" accretion disk simulations and traditional "ISM-scale" simulations and the transition in between, all in the same simulation.
Our most notable conclusions from this particular study include: • Magnetic fields play a key role: We show that magnetic fields are critical for a wide range of effects on sub-pc scales within the accretion disk, ranging from maintaining efficient torques and high inflow rates, explaining the scale heights and vertical profiles of the disk structure, the outer size/boundary of the accretion disk, and perhaps most importantly the suppression of star formation at sub-pc scales.Without magnetic fields, a disk still forms, but it is an order of magnitude or more smaller in spatial scale and mass, and produces factor of ≳ 100 times lower accretion rates into the SMBH, with runaway fragmentation and ordersof-magnitude larger nuclear SFRs on ∼ 0.1 − 10 pc scales.
The accretion disk that forms also has qualitatively different structures as a result of magnetic flux-freezing and flux-feeding from the ISM, to be studied in detail in Paper II.
• Quasar-Level Inflow Rates are Plausible and Can be Maintained: Extending previous studies like those in Anglés-Alcázar et al. ( 2021) with both much higher resolution and more detailed micro-physics relevant on small scales, we confirm that strong torques on sub-kpc scales can maintain inflow rates as large as ≳ 10 M ⊙ yr −1 into a QSO accretion disk at < 80 au, for extended periods of time (hundreds of thousands of dynamical times at the smallest radii simulated here).On scales ∼ pc-kpc these are dominated by "gravitational torques" in a multi-component (gas+stellar) disk, inducing strong shocks and inflow.On sub-pc scales where star formation becomes inefficient these become weak but strong MHD torques (Maxwell+Reynolds stress) in a turbulent, strongly-magnetized outer flux-fed accretion disk take over and are able to sustain such large inflow rates down to the smallest resolved radii here, well within the "traditional accretion disk" range of scales.
• Suppression of star formation: On sub-pc scales, star formation is strongly suppressed.Models which simply assume some fixed star formation efficiency per free-fall time in sufficiently dense and/or self-gravitating gas (the standard on galactic scales and in our FIRE simulations) would not necessarily fully capture or be able to predict this effect a priori, but on these scales the simulations here explicitly resolve individual proto-stellar cores (at < 0.01  ⊙ mass resolution) with models for resolved single-star formation from STARFORGE that would, if star formation were "missed" incorrectly, allow the gas to collapse to infinitely high densities.Just as important, for the first time this prediction is made using physics and numerical methods which have been explicitly shown to reproduce reasonable observed star formation efficiencies, stellar masses/IMFs, and stellar multiplicity distributions under typical Solar-neighborhood ISM/GMC conditions.With these physics, we show that a combination of increasing optical depths producing warmer gas in the galactic nucleus, plus (crucially) strong toroidal magnetic fields raising the magnetic critical mass to be larger than the disk mass and strongly suppressing gravitoturbulent fragmentation, leads to a dramatic, almost complete suppression of star formation at distances ≪ pc from the SMBH.
There are many properties which could and should be studied in more detail in the simulations here.In Paper II, we explore the structure of the strongly-magnetized and flux-frozen/fed accretion disk, nature of the MHD (Maxwell/Reynolds) torques, origin and dynamics of these strong fields on ≪ pc scales, and their consequences for accretion disk theory.In Paper III we explore the detailed predictions for the dynamics and process of star formation and the resulting IMF of stars in the hyper-resolved inner region around the QSO accretion disk (the "circum-quasar medium").There are obvious extensions of the work here and therein modeling many different observables, contrasting how, for example, the much more detailed radiation-magneto-thermochemistry models produce different predictions for dust and atomic and molecular gas properties (and hence observables) from galactic nuclei and the quasar "torus" region, compared to previousgeneration simulations with more simplified physics.
In future work, there are also multiple ways one might extend the actual simulation work here.The most obvious is to consider other initial conditions of different galaxies at different times, to explore other regimes.Perhaps the biggest caveat here is that (owing to the computational expense of these simulations) we have studied just one case, so it is not obvious how much of our conclusions can be generalized to very different conditions with e.g.much lower accretion rates, let alone extremely low-accretion rate systems like M87 or Sgr A * .Another obvious limitation of the work here is that we do not include any "outward" fluxes from the un-resolved accretion disk at < 80 au, e.g.radiation or jets from the inner disk.Our hope is that, given the unexpected properties of the disks which form here, this work will first motivate smallerscale simulations of accretion disks (going down to the ISCO) with outer boundary conditions broadly similar to our inner boundary conditions, and those can provide some motivation for including such "inner disk feedback" prescriptions in a subsequent generation of simulations.
In principle, one could also attempt to improve the simulations here even further in resolution or run-time.However the simulations here already push the boundaries of what is possible, and it is not obvious such a strategy would be most efficient.Instead, we argued that conclusions on large (galactic) scales (where one would ideally like to run the simulations for much longer) are better studied in simulations without "hyperrefinement" (but using simulations like those here to inform the sub-grid models for BH accretion and feedback).Mean-while one can much more efficiently explore the physics and parameter space of accretion disk physics with classical, dedicated accretion disk simulations.But in those simulations, the initial and boundary conditions are largely arbitrary, so our goal here is to provide predictive values for those, motivating qualitatively new parameter space to be studied in future work.
Fig. 1.-Series of images of the projected gas density in our simulation ( § 2.4) at one moment in time at redshift  ≈ 4 typical of when we analyze it.Color encodes surface density increasing black-to-white on a logarithmic scale (each panel rescaled owing to the different dynamic range) -a median pixel in the largest-scale panel (top-left) has column   ∼ 10 19 cm −2 (density   ∼ 10 −5 cm −3 ), while in the smallest-scale panel   ∼ 10 27 cm −2 (  ∼ 10 12 cm −3 ).We see structure on all scales, with a chaotic, cold, disordered morphology on most scales until an ordered disk forms from capture of gas from a passage of a giant molecular cloud complex (itself triggered by an ongoing galaxy merger in the rapidly-accreting proto-galaxy), forming the accretion disk at ≲ 0.1 pc.
Fig.2.-As Fig.1, but tiling the images so more structure can be seen and identifying each with the heuristic label appropriate to the range of scales shown, per § 3.1.In order, each image zooms in by a factor of 10 around the previous image, with side-length  = (1000, 100, 10) kpc (top), (1000, 100, 10) pc (middle), (1, 0.1, 0.01) pc (bottom).The projection here is chosen to be face-on to the innermost central disk.Note the "hole" in the latter inside  ≲ 80 au is caused by our inner accretion boundary (dashed circle).
Fig.3.-AsFig.2, but in stars.The chaotic merger morphology at a few kpc, and clumpy, highly asymmetric stellar morphology driving gravitational torques on the gas is evident on all scales.In most panels we show a continuous projection of stellar density, but in the last panel this breaks down (the inter-stellar separation is no longer much smaller than a pixel) so we show individual O-stars.We do not show images at ≪ 1 pc because there is a negligible stellar mass compared to the gas on these scales.

Fig. 5 .
Fig. 5.-Images of the gas (left) and stars (right) with different scales and their approximate naming label conventions from § 3.1 shown.The images show BH-centric radius  increasing from bottom-to-top (the vertical axis) on a logarithmic scale as labeled.The horizontal axis shows cos  ≡ / from −1 to +1(defined so  = 0 corresponds to the midplane of the inner disk), in a wedge of azimuthal opening angle sin  < 0.3.For gas (left) colors denote different phases  < 10 3 K (green), 10 3 <  < 10 4 K (yellow), 10 4 <  < 10 5 K (magenta), 10 5 <  < 10 6 K (purple),  > 10 6 K (cyan).We see other galaxies on IGM scales, the virialized CGM with accretion in warm clumps/filaments, the highly clumpy/inhomogeneous/asymmetric and multi-phase structure in the galaxy and (thermally colder, primarily atomic+molecular) galaxy nucleus, settling into the more ordered (but still visibly thick and turbulent) non-star forming disk and BH accretion disk on sub-pc scales.
Fig.6.-Effectiveresolution of the simulation as a function of radial distance from the BH ().We show the median, 50% and 90% inclusion intervals (shaded) of the mass resolution (gas cell mass Δ; left), spatial resolution (equivalent cell size Δ ≡ (Δ/) 1/3 ), and time resolution (numerical timestep Δ), at each .The bulk of the galaxy is resolved with Δ < 10 4  ⊙ ; from 1 − 100 pc the resolution is rapidly refined as  decreases, with the target resolution of Δ < 0.01  ⊙ and Δ < 10 −3 pc reached for the particles at  ≲ pc scales.For Δ, Δ, and Δ we compare the global enclosed gas mass  gas ( ), radius , and dynamical time  dyn at each .In Δ, we also denote the region where the simulation will form STARFORGE single-star/sink particles, as opposed to FIRE SSP particles.In Δ, we also denote as "duration" approximately how long the simulation was run after reaching the maximum refinement level at each radius.
Fig. 7.-Top Left: Circular velocity profile  c ≡√︁   enc (<  )/ versus spherical distance from the SMBH , with contributions from different mass components.The BH excision radius truncates the gas mass distribution at ≲ 0.001 pc.We clearly see a transition from dark matter dominating outside the galaxy, stars dominating within the galaxy, until the BHROI at ∼ a few pc, and the cessation of star formation giving a gas-dominated disk at ≲ 0.1 pc.Top Right: Projected mean gas, stellar, dark matter, and "star formation rate" (defined by mass of stars formed in the last  dyn ≡ 1/Ω at each radius) surface density Σ, in cylindrical shells of different projected radii  (plotted down to radii where we have at least > 1000 particles of each type).We see the same transitions.While there are large deviations at any radius, Σ gas ∝  −1 approximates the scaling reasonably well from ∼ 10 −3 − 10 6 pc.In this and subsequent plots we independently determine the "midplane" in each annulus from the angular momentum vector of the gas in the annulus.The suppression of Σ SFR at small radii is discussed below ( § 3.5.3).Bottom Left: Three-dimensional gas density  versus spherical radius .We show both the volume-weighted density mean (line) in shells and mass-weighted mean (shaded range shows the volume-weighted 90% inclusion interval; the mean can exceed this upper limit if the distribution has large tails).The gas density profile is crudely isothermal (in a very order-of-magnitude sense),  ∝  −2 , from ∼ 10 −3 − 10 6 pc.Mass-weighted densities are systematically enhanced relative to volume-weighted but also vary more owing to clumping and phase structure.Bottom Right: Mass flow rate , showing the total inflow rate  in and total outflow rate  out through each annulus, and the cumulative SFR summed within each radius  * (<  ), versus spherical .Grey bar shows the 50% (solid; ∼ 18 − 30 M ⊙ yr −1 ) and 90% (dotted; ∼ 15 − 73 M ⊙ yr −1 ) range of inflow rates at < 100 au over the duration of our highest-resolution simulation.Lack of spherical symmetry and non-equilibrium dynamics mean that inflow and outflow co-exist and can change relative sign (e.g. the "outflow" at ≲ pc scales is mostly just coherent eccentric motion), but despite this a remarkably stable order-of-magnitude inflow rate to the SMBH of  in ∼ 10 − 100 M ⊙ yr −1 persists.
Fig.8.-Radial profiles of related timescales and bulk flows as Fig.7.Left: "Depletion" (or accretion) timescales for inflow/accretion, outflow, and star formation, in different radial annuli, defined as shown (with Δ  ≡   / ln ), in units of the galactic dynamical time ( dyn ≡ Ω −1 ) at each radius.Star formation is efficient at Galactic radii from ∼ pc to ∼ kpc, as expected.Inflows are dynamical at essentially all radii, and produce net inflow at sub-kpc scales during this strong-inflow phase.Right: Gas mass-weighted mean radial (  ), azimuthal/rotational (  , with the φ axis defined in each annulus by the net angular momentum vector as Fig.7), and polar (  ) velocities, and velocity dispersions   , in units of   .We see a mix of inflow/outflow motions (large   ) dominate at large radii (≫ kpc), partial rotation in a kinematically hot/thick galaxy configuration at Galactic radii (∼ 0.1 − 10 kpc), and near free-fall at the radii interior to the BHROI where the GMC-like gas complex fueling the SMBH is tidally disrupted and captured (∼ 1 − 10 pc), which then circularizes at sub-pc scales to form a kinematically cold, rotation-dominated disk at ≲ 0.1 pc scales.

Fig. 9 .
Fig.9.-Images of the gas phase structure as Fig.5, but showing distance from the BH along the accretion-disk midplane r direction (horizontal axis, in a wedge with azimuthal angle | sin  | < 0.3) versus ẑ (vertical axis), with both  and |  | stretched logarithmically (so both expand at the same rate with distance from the origin).Here we can see the highly-globally-asymmetric (and multi-phase) structures driving inflow on a wide range of scales.

Fig. 11 .
Fig. 11.-Two-dimensional projections of the mass-weighted mean thermochemical properties (temperature  in K; density  compensated by  1.7 for the sake of visualization, in units of cm −3 kpc 1.7 ; free electron fraction   ; plasma ) from Fig. 10.Horizontal axis shows distance  from the SMBH, log-scaled, along a wedge with opening angle | sin  | < 0.3, as Fig. 9, while the vertical axis shows the height / as Fig. 5.The wedge is oriented along the axis of the inner accretion disk.The radial trends from Fig. 10 are evident as is highly multi-phase structure at most radii.
Fig.13.-AsFig.2, but showing the edge-on projection to the inner disk at the same time.We can more clearly see the formation of the thin disk at sub-pc scales from capture via filamentary, misaligned inflow from a highly chaotic ISM/CGM at larger radii.
Fig.14.-Angular momentum direction of the gas (cos  ≡ ĵ ≡ ĵ • ĵinner ), relative to the mean angular momentum direction of the inner disk j inner (averaged within < 0.01 pc), as a function of BH-centric radius .We clearly see the aligned, dynamically cold disk (cos  ≈ 1, with relatively little variation) in the central ≲ 0.1 pc, with an un-aligned, and much more spherical/kinematically hot (broad distribution of cos ) at larger radii.
Fig.15.-TopLeft: Instantaneous (gas mass-weighted) mean torques  in the direction of the mean angular momentum vector ĵ within each radial annulus (normalized to    Ω =  2  ) in shells as a function of distance .We restrict to cool gas with  < 10 4.5 K but the trends are qualitatively similar regardless.Shaded region shows the ∼ 90% inclusion interval.We plot the torques directly from the simulation from gravitational forces, radiation pressure forces, and MHD (magnetic, thermal, turbulent) forces.At all radii torques are efficient, implying angular momentum loss in a few orbital times.At large radii gravitational torques from stars on gas, plus MHD torques driven by stellar feedback, dominate.Where there are few stars, magnetic and turbulent/Reynolds torques take over.Top Right: Fractional (Fourier) mode amplitudes of asymmetric modes in the face-on projected gas surface density (Σ ( , ) ≡ Σ 0 (1 +    cos (  +  0,  ) )) in cylindrical annuli.At radii ≫ 0.1 pc there are order-unity asymmetries, often dominated by the lowest- modes (global asymmetries rather than just small-scale clumping).At small, increasingly Keplerian radii |   | decreases but only modestly, to ∼ 0.1 at ≲ 10 −3 pc.Bottom Left: Profile of different (volume-weighted) mean components of the magnetic stress tensor  mag ≡ (1/4  ) ( |B| 2 I − BB/2), in spherical coordinates, versus radius.We normalize to the mean value of the total stress tensor  at each radius.Solid (dotted) line correspond to  > 0 ( < 0).Bottom Right: Same, for the (total) kinetic stress  kin ≡  v v (note this is distinct from e.g. a Reynolds stress).

TABLE 1
Summary of physics included in our default simulation.