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

# FORGE’d in FIRE III: The IMF in Quasar Accretion Disks from STARFORGE

Philip F. Hopkins1,∗ Michael Y. Grudić2,† Kyle Kremer1,† Stella S. R. Offner3
Dávid Guszejnov3,† Anna L. Rosen4 1TAPIR, Mailcode 350-17, California
Institute of Technology, Pasadena, CA 91125, USA  2Carnegie Observatories, 813
Santa Barbara St, Pasadena, CA 91101, USA  3Department of Astronomy,
University of Texas at Austin, TX 78712, USA  4Center for Astronomy & Space
Sciences, University of California, San Diego, La Jolla, CA 92093, USA

###### Abstract

Recently, we demonstrated self-consistent formation of strongly-magnetized
quasar accretion disks (QADs) from cosmological radiation-magnetohydrodynamic-
thermochemical galaxy-star formation simulations, including the full STARFORGE
physics shown previously to produce a reasonable IMF under typical ISM
conditions. Here we study star formation and the stellar IMF in QADs, on
scales from $100\,$au to $10\,$pc from the SMBH. We show it is critical to
include physics often previously neglected, including magnetic fields,
radiation, and (proto)stellar feedback. Closer to the SMBH, star formation is
suppressed, but the (rare) stars that do form exhibit top-heavy IMFs. Stars
can form only in special locations (e.g. magnetic field switches) in the outer
QAD. Protostars accrete their natal cores rapidly but then dynamically
decouple from the gas and “wander,” ceasing accretion on timescales $\sim
100\,$yr. Their jets control initial core accretion, but the ejecta are “swept
up” into the larger-scale QAD flow without much dynamical effect. The strong
tidal environment strongly suppresses common-core multiplicity. The IMF shape
depends sensitively on un-resolved dynamics of protostellar disks (PSDs), as
the global dynamical times can become incredibly short ($\ll$ yr) and tidal
fields are incredibly strong, so whether PSDs can efficiently transport
angular momentum or fragment catastrophically at $\lesssim 10\,$au scales
requires novel PSD simulations to properly address. Most analytic IMF models
and analogies with planet formation in PSDs fail qualitatively to explain the
simulation IMFs, though we discuss a couple of viable models.

###### keywords:

stars: formation — galaxies: star formation — galaxies: starburst — quasars:
general — accretion, accretion disks — ISM: general

##  1 Introduction

How stars form, and what shapes their masses and other characteristics, are
fundamental questions that inform almost every subfield of astrophysics. A key
parameterization is the “initial mass function” (IMF) of stars - the
distribution of stellar masses at the time of their formation.
Observationally, the IMF exhibits remarkable uniformity across different
observed environments (for reviews, see Bastian et al., 2010; Kroupa et al.,
2011; Offner et al., 2014; Hopkins, 2018). In the local Solar neighborhood and
similar “typical” galactic environments, a universal IMF is perhaps not
surprising, as the properties of progenitor giant molecular clouds (GMCs) and
the typical interstellar medium (ISM) fall along well-defined scaling
relations within a relatively narrow parameter space (Larson, 1981; Solomon et
al., 1987; Rosolowsky et al., 2003; Bolatto et al., 2008; Goodman et al.,
2009; Heyer et al., 2009; Pokhrel et al., 2020). But this makes it all the
more important to explore extreme environments which lie well outside of this
space. And there are some hints of (modest) IMF variation in special
environments: for example towards more bottom-heavy (with a lower-mass
turnover) in some massive galaxy bulges (van Dokkum & Conroy, 2011; Spiniello
et al., 2012), or towards more top-heavy (“standard” but with a slightly
shallower-than-Salpeter slope) in a couple of special environments very close
to the supermassive black hole (SMBH) Sgr A∗ in our Galactic center (Nayakshin
& Sunyaev, 2005; Paumard et al., 2006; Maness et al., 2007; Espinoza et al.,
2009; Bartko et al., 2010; Löckmann et al., 2010; Lu et al., 2013; Hosek et
al., 2019) and perhaps (much more tentatively) in similar environments around
luminous quasars (e.g. Toyouchi et al., 2022; Fan & Wu, 2023).

These Galactic center observations, in particular, have prompted a number of
theoretical studies (e.g. Nayakshin et al., 2007; Bonnell & Rice, 2008; Hobbs
& Nayakshin, 2009; Hocuk & Spaans, 2011; Hocuk et al., 2012; Mapelli et al.,
2012; Alig et al., 2013; Davies & Lin, 2020), which have largely argued that
the “excess” population of intermediate-mass stars in that neighborhood formed
in some kind of accretion disk at $\lesssim 0.2\,$pc from the SMBH during an
earlier accretion episode. Moreover, in recent years there has been a surge of
theoretical interest in star formation in quasar111In this manuscript we will
use the term “quasar” loosely to refer to any SMBH accreting at high
luminosities, as opposed to e.g. the traditional optical luminosity and/or
spectral criteria. accretion disks (QADs) and the circum-quasar medium
(CQM),222In this paper, for the sake of clarity we will refer to the medium on
scales $\ll 10\,$pc around the SMBH as the “circum-quasar medium” (CQM), and
within the CQM we refer to the global, nuclear disk of gas and stars orbiting
the SMBH in the galaxy center as the “quasar accretion disk” (QAD). The CQM we
are careful to distinguish from the more “typical” interstellar medium (ISM or
LISM) and giant molecular clouds (GMCs) observed in the Solar neighborhood. We
are also careful to distinguish individual circum-stellar or proto-planetary
or proto-stellar disks (themselves potentially embedded in the larger CQM or
QAD) – which we will refer to as “protostellar disks” (PSDs) throughout – from
the QAD itself. When we refer to a global “position” or “radius” $R$, we refer
always to a coordinate frame in which the SMBH is at the origin (e.g. a radius
$R$ from the SMBH), and we often refer to a cylindrical coordinate frame based
on the orientation of the inner QAD (so $\phi$ is the prograde azimuthal angle
rotating with the inner QAD and $z$ the vertical/angular momentum direction).
as a potential means to form hyper-massive and/or infinitely-long lived
accreting stars or compact binaries which can act as a source for
gravitational wave events (for LIGO/LISA) or tidal disruption events (TDEs) or
other types of astrophysical transients (Cuadra et al., 2009; Stone et al.,
2017; Tagawa et al., 2020; Baruteau et al., 2011; Bortolas, 2022)

But star formation – even in “normal” Solar-neighborhood environments – is a
highly non-linear, chaotic, multi-physics problem which necessitates detailed
numerical simulations (see discussion in McKee & Ostriker, 2007; Offner et
al., 2014). It has only recently become possible for “first principles”
simulations of star formation under Solar neighborhood/local ISM (LISM)
conditions to reproduce the full behavior of the IMF after “completion” of
star formation events in a GMC, and those simulations have argued that this
depends on physics including self-gravity; hydrodynamics and magnetic fields;
explicit radiation-hydrodynamics accounting for multiple wavelengths and self-
consistently able to handle the transition between optically thin and thick
cooling; detailed (potentially non-equilibrium) thermo-chemistry of dusty,
partially-ionized atomic and molecular gas; self-consistent driving of
turbulence, GMC surface densities, and other large-scale outer boundary or
initial conditions; and detailed feedback from proto-stars and (pre)-main-
sequence stars including (proto)stellar jets, stellar radiation, and stellar
mass-loss (Bate, 2012; Lee & Hennebelle, 2019; Hennebelle et al., 2020;
Guszejnov et al., 2021; Grudić et al., 2021; Lee et al., 2021; Grudić et al.,
2022; Guszejnov et al., 2022a, b, 2023). Meanwhile, simulating the outer
regions of QADs (where star formation may be possible) clearly requires all of
these physics as well, with the structure of the QAD sensitive to magnetic
field structure (both strengths and geometry/topology), and QADs clearly
necessitating treatments that can handle radiation-pressure-dominated limits;
dust sublimation and its effects on thermo-chemistry; and the ability to
follow global modes such as spiral arms and gravito-turbulence (see Davis &
Tchekhovskoy 2020 for a recent review).

This highlights two major areas for improvement on previous-generation
simulations of star formation in QAD/CQM environments. (1) All of the previous
QAD IMF studies cited above used ad-hoc initial/boundary conditions, as they
were modeling some portion of an arbitrarily analytically constructed QAD, as
opposed to predicting these conditions self-consistently from simulations of
larger scales actually producing inflows and forming the QAD. (2) None of the
QAD IMF studies cited above included a combination of physics and numerical
methods which has been actually shown to reproduce the Solar neighborhood IMF.
Indeed, the vast majority of them considered only hydrodynamics (i.e. ignored
magnetic fields), with idealized equilibrium thermochemistry (or even further
simplified ideal gas equations-of-state) neglecting explicit multi-band
radiation hydrodynamics, and did not include any kind of explicit stellar
evolution or feedback. It is well-established that this combination of physics
alone will, even in Solar-neighborhood environments, produce a systematically
incorrect IMF, incorrect star formation rates/efficiencies, and far too much
cloud-to-cloud or cluster-to-cluster IMF variation compared to observations
(Guszejnov et al., 2016; Guszejnov et al., 2017a, b, 2018b, 2019, 2020).

In this paper, we therefore present a first attempt to address both of these
issues. Regarding issue (1), in Hopkins et al. (2023a) (hereafter Paper I), we
presented the first fully-cosmological galaxy and star formation simulations
which employed an adaptive hyper-refinement technique to self-consistently
resolve the formation of a QAD (the properties of which were studied in detail
in Hopkins et al. 2023b, hereafter Paper II) around a SMBH of $\sim
10^{7}\,M_{\odot}$ during a bright quasar phase (resolving the interior
structure of the QAD down to scales of $<80$ au, or $\sim 300$ Schwarzchild
radii, from the SMBH). Regarding issue (2), these simulations utilize the
STARFORGE modules and physics, including all of the specific physics discussed
above (e.g. detailed non-equilibrium ion/atomic/molecular/dust thermo-
chemistry coupled to multi-band explicit radiation transport; magnetic fields;
proto-stellar and main-sequence stellar evolution with feedback from jets,
radiation, winds, and supernovae), which has been explicitly shown to
reproduce reasonable IMFs with minimal variation and higher-order properties
(e.g. clustering, metallicity, star formation efficiencies) plausibly
consistent with observational constraints for “typical” Solar-neighborhood-
like environments (Guszejnov et al., 2022a, b, 2023). Furthermore, in Paper I
and Paper II, we showed that we resolve the transition to a weakly star-
forming (i.e. not catastrophically-fragmenting) QAD, as required to sustain
high accretion rates for quasar activity, and that the QAD is strongly
magnetized with magnetic fields playing a fundamental role in the dynamics in
a manner qualitatively distinct from the assumptions in the historical
numerical and analytic studies cited above. This provides yet more reason to
further explore the IMF in these conditions. Here, we present a first
systematic study of the stellar IMF in both the QAD and CQM environments in
these simulations, and explore how it depends qualitatively on physics not
included in most previous studies.

In § 2, we review the numerical methods used in these simulations (§ 2.1),
with special attention to how un-resolved PSDs/accretion flows and their
accretion/feedback processes are treated around individual sink particles (§
2.2). In § 3, we present and review from Papers I & II some of the qualitative
properties of the QAD and CQM and how they differ from the LISM, including:
their surface density, density, and opacity scales (§ 3.1); magnetic field
strengths (§ 3.2); dust, thermal gas, and radiation temperature structures (§
3.3); turbulent properties (§ 3.4); and tidal fields (§ 3.5). In § 4 we
present and discuss the key results of the simulations. We first focus on the
sites of star formation (§ 4.1), in infalling CQM filaments (§ 4.1.1) and the
outer QAD (§ 4.1.2); we then discuss how these special regions are able to
overcome magnetic support and collapse to form stars (§ 4.2). We explore the
timescales over which stars form and accrete in § 4.3, and show how this
relates to the formation of a large population of “wandering stars” in § 4.3.1
and has crucial effects on stellar accretion and feedback in § 4.3.2. We
discuss the global influence of (proto)stellar feedback on the environment in
§ 4.4. We then consider the IMF directly (§ 4.5), showing its general
properties (§ 4.5.1), dependence on distance from the SMBH (§ 4.5.3) and
evolution in time (§ 4.5.4), exploring how it depends on physics such as the
un-resolved treatment of PSDs on $\sim$ au scales (§ 4.5.2, with further
discussion of their fragmentation and the system-versus-single star IMF in §
4.5.5) and magnetic fields (§ 4.5.6). We briefly discuss multiplicity in §
4.5.7. We then discuss the physics driving the IMF shape (§ 4.6), specifically
reviewing a large number of simple analytic scalings that have been proposed
to explain the dependence of the IMF on environment which are clearly not
describing our simulations (§ 4.6.1) as well as physical explanations which
appear to reproduce the simulations in both the weak feedback/slow PSD-
accretion limit (§ 4.6.2) and strong feedback/fast PSD-accretion limit (§
4.6.3). With this in mind we then review the extent to which star formation in
the QAD is (or is not) analogous to planet formation in a “typical” PSD (§
4.7). We discuss some caveats of resolution in § 5, and then conclude in § 6.

Figure 1: Images of the projected gas density in our simulation (§ 2) at one
moment in time (redshift $z\approx 4.4$) when we “zoom in,” projected face-on
to the QSO accretion disk (QAD). Color encodes surface density increasing
dark-to-light on a logarithmic scale (each panel rescaled owing to different
dynamic range) – a median pixel in the largest-scale panel has $N_{H}\sim
10^{19}\,{\rm cm^{-2}}$ (density $n_{H}\sim 10^{-5}\,{\rm cm^{-3}}$), while in
the smallest-scale panel $N_{H}\sim 10^{27}\,{\rm cm^{-2}}$ ($n_{H}\sim
10^{12}\,{\rm cm^{-3}}$). We label the different characteristic scales. From
the cosmic web, a galaxy forms on $\sim$ kpc scales at high redshifts, with
rapid star formation (galactic SFR $>200\,{\rm M_{\odot}\,yr^{-1}}$) in a
merging, clumpy, gas-rich galaxy. A large star-forming cloud complex of mass
$\gtrsim 10^{8}\,M_{\odot}$ has a close passage to a pre-existing $\sim
10^{7}\,M_{\odot}$ SMBH, and is partially tidally disrupted around the BH
radius of influence (BHROI) at a few pc. The bound mass falls in as a tidally
compressed stream/filamentary structure in the circum-quasar medium (CQM),
which circularizes at $\sim 0.1\,$pc to form the QAD. The QAD persists down to
our inner boundary condition at $r=80\,$au distance from the SMBH. In this
paper we study SF interior to the BHROI, where the simulation resolves
individual (proto)stars.  Figure 2: Illustration of the gas phase structure in
the same snapshot as Fig. 1, projected edge-on to the QAD. We show distance
from the SMBH along the QAD midplane ($\hat{R}$ direction; horizontal axis),
in a wedge of azimuthal angle with $|\sin{\phi}|<0.3$, versus the vertical
distance along the angular momentum direction ($\hat{z}$; vertical axis), with
both $r$ and $|z|$ stretched logarithmically. Colors encode gas temperature:
$T<10^{3}\,$K (green), $10^{3}<T<10^{4}$ K (yellow), $10^{4}<T<10^{5}\,$K
(magenta), $10^{5}<T<10^{6}\,$K (purple), $T>10^{6}\,$K (cyan). We see multi-
phase self-gravitating gas clumps/filaments appear on all scales from $\sim
10\,$kpc (“galactic” GMCs) to $\sim 0.1\,$pc, but the structure becomes
smoother and more homogeneous with a warmer $\gtrsim 10^{3}\,$K temperature
within the QAD itself (see § 3). We label the range of radii corresponding to
the QAD and CQM, and denote where the resolution is sufficiently high for our
STARFORGE resolved-star-formation model. At larger radii, individual stars are
un-resolved and the simulation uses the FIRE model for un-resolved
fragmentation to an assumed IMF, so this part of the simulation informs our
initial/boundary conditions but cannot be used to predict the IMF.

##  2 Methods

###  2.1 Review and Key Physics

The simulations studied here are presented and extensively described in Paper
I. Briefly, we begin from a $\sim(100\,{\rm cMpc})^{3}$ cosmological periodic
box at redshift $z\sim 100$ with a primordial trace magnetic field, and follow
it as a cosmological galaxy formation simulation following the combined
Feedback In Realistic Environments (FIRE) project methods from Hopkins et al.
(2018, 2023d) and STARFORGE physics treatment from Grudić et al. (2021);
Guszejnov et al. (2021). At a redshift $z\sim 4.5$ when a period of violent
activity induces large inflows into the central $\sim\,$kpc of the galaxy, we
hyper-refine to go to higher and higher resolution, reaching sufficiently high
resolution to resolve individual (proto)star formation, accretion and
evolution and PSD structure in the central $\lesssim 10-100\,$pc of the
galaxy. We continue to integrate self-consistently for a mix of “individual
star” (STARFORGE) and “stellar population” (FIRE) particles given the
methodology in Paper I, and continue to refine to a target resolution of
$\Delta m<0.01\,{\rm M_{\odot}}$ in the central $\lesssim 1-10\,$pc, to follow
gas inflows and QAD formation down to $<300$ Schwarzschild radii around the
super-massive black hole of mass $\sim 1.3\times 10^{7}\,{\rm M_{\odot}}$.
Fig. 1 shows an image of the gas properties on some of the wide range of
scales simultaneously resolved in the simulation, illustrating the circum-BH
QADs which form on sub-pc scales and form the focus of our study here.

The simulations include a wide range of physics including magnetic fields
(using the high-order constrained-gradient method from Hopkins & Raives 2016;
Hopkins 2016), with kinetic (anisotropic Braginskii viscosity and conduction)
effects (Su et al., 2017; Hopkins, 2017); a variable local cosmic-ray
background (Hopkins et al., 2022a; Hopkins, 2023; Hopkins et al., 2022c;
Hopkins et al., 2022d) here modeled using the simple sub-grid method from
Hopkins et al. (2023e); self-gravity with adaptive self-consistent softenings
scaling with the resolution and high-order Hermite integrators capable of
accurately integrating $\gtrsim 10^{5}$ orbits in hard binaries (Grudić &
Hopkins, 2020; Grudić et al., 2021; Grudić, 2021; Hopkins et al., 2023f);
metal enrichment and dust and dust destruction/sublimation (Ma et al., 2017;
Gandhi et al., 2022; Choban et al., 2022); super-massive black hole seed
formation and growth via gravitational capture of gas (Hopkins et al., 2016;
Shi et al., 2022; Wellons et al., 2023); (proto)star formation and accretion
and explicit feedback from stars in the form of protostellar jets, main-
sequence stellar mass-loss, radiation, and supernovae (Grudić et al., 2022;
Guszejnov et al., 2022a, b, 2023). The simulations evolve explicit multi-band
radiation-hydrodynamics with adaptive-wavelength bands (Hopkins et al., 2020;
Hopkins & Grudić, 2019; Grudić et al., 2021) coupled explicitly to all the
thermo-chemical processes, together with radiative cooling and thermo-
chemistry incorporating cosmic backgrounds, radiation from local stars, re-
radiated cooling radiation, dust, molecular, atomic, metal-line, and ionized
species opacities and processes, cosmic rays, and other processes. This allows
us to self-consistently model the thermochemistry and opacities in gas with
densities from densities $n\sim 10^{-8}-10^{16}\,{\rm cm^{-3}}$ and
temperatures $\sim 1-10^{10}\,$K in a range of radiation and cosmic ray
environments (see Paper I). Fig. 2 illustrates some of the complex phase
structure which emerges even in just the nuclear regions. As we show below,
the most important physics on these scales are gravity, MHD, and radiation-
thermodynamics, coupled to accretion onto and feedback from (proto)-stars.

Protostellar sink particle formation and accretion follow the methods in
Grudić et al. (2021) for STARFORGE. Briefly, sink particles form with a mass
of their parent gas cell, if and only if (1) they are above some minimum
density; (2) they represent the density maximum among interacting gas
neighbors; (3) there is no other sink within the interaction radius; (4) the
cell density is increasing and velocity divergence is negative; (5) the cell
is fully bound/self-gravitating at the resolution scale, including its
turbulent, magnetic, and thermal support; (6) the free-fall time for collapse
is faster than the tidal timescale or accretion/orbital timescale onto any
other sink; and (7) the tidal tensor at the cell center-of-mass is fully
compressive (possesses three negative eigenvalues), which ensures it is the
local gravitational center/focus and can collapse even in extreme tidal
environments like those here. Once formed, sinks capture and remove from the
gas flow neighboring gas cells if and only if: (1) those cells fall within the
sink capture radius ($=80\,$au); (2) the gas cell is bound to the sink
(including its kinetic, thermal, and magnetic energy relative to the sink);
(3) the cell possesses less angular momentum than a circular orbit around the
sink at its separation; (4) its apocentric orbital two-body radius also falls
inside $R_{\rm sink}$; (5) the volume of the cell is smaller than the volume
enclosed by $R_{\rm sink}$; and (6) the gas cell is not eligible for capture
by any other sink with a shorter free-fall time onto that sink.

On all scales we study in detail in this paper, the refinement has reached the
target resolution of $<0.01\,{\rm M_{\odot}}$, and so the simulation is fully
in the STARFORGE limit (with no FIRE “stellar population” particles). These
and other resolution properties are shown explicitly in Fig. 3. This is
necessary for any meaningful predictions about the IMF, since the FIRE model
used at much lower resolution in the Galactic radii $\gg 10\,$pc assumes a
well-sampled Milky Way IMF, from which to model stellar populations. In the
most-dense regions around collapsing proto-stars, this mass resolution
corresponds to a local spatial resolution $\sim 10\,{\rm au}$ and time
resolution as small as $\sim$ days. We discuss some potential effects of
resolution on our conclusions in § 5.

###  2.2 Sub-Grid PSD/Reservoir Accretion Rates

In the simulations here, when gas is captured by a (proto)stellar sink
particle, it is first added to a “reservoir” of mass $M_{\rm PSD}$ which
represents the sub-grid (by definition un-resolved in the simulation) portion
of the PSD/accretion flow. This then accretes onto the (proto)star of mass
$M_{\ast}$ at some (necessarily) prescribed rate $\dot{M}_{\ast}=f_{\rm
acc}\,M_{\rm PSD}/t_{\rm dep}$, where $t_{\rm dep}$ is the assigned depletion
time, and $f_{\rm acc}$ is the fraction which actually is retained in the
(proto)star as opposed to ejected in outflows/protostellar jets. In previous
studies of GMCs with typical Solar-neighborhood-like conditions, Guszejnov et
al. (2018b, 2021, 2022a); Grudić et al. (2021) showed that predictions for
stellar properties (the IMF, multiplicity) as well as cloud-scale properties
(global star formation efficiencies, cloud dynamics) were relatively
insensitive to the details of the model for $t_{\rm dep}$, so long as $t_{\rm
dep}$ is relatively “short” compared to the global timescales of interest
(e.g. cloud lifetimes, timescales for massive stars to accrete all of their
mass from large scales around their initial core, etc., of order $\sim$ Myr).
The default STARFORGE model is therefore a relatively simple prescription
based on an isothermal Shu (1977) sphere in LISM conditions, which – inserting
typical values for our resolution and the assumed default simulation
parameters – gives a more or less constant $t_{\rm dep}$ between
$2000-5000\,{\rm yr}$. This easily satisfies the “sufficiently fast’ condition
in LISM conditions.

However, as we note below, on the scales we follow, the global
dynamical/orbital times around the SMBH can be as short as months or even
days. This also necessarily means any PSD/PPD/accretion reservoir at this
radius must have an even shorter internal dynamical time or else it would be
tidally disrupted. So assuming the same “default” $t_{\rm dep}$ leads to
accretion timescales from the PSD/reservoir to the star which are very large
compared to the global dynamical time, and thus mass can rapidly accumulate
and the “reservoir” mass $M_{\rm PSD}$ can greatly exceed $M_{\ast}$, which is
probably unphysical. We therefore denote this as the “slow” accretion regime.
Given that this is, by definition, unresolved behavior in our simulations, a
detailed study of the “correct” accretion rates or $t_{\rm dep}$ – and an
ultimate determination of whether these PSDs fragment efficiently (indeed
whether a PSD even forms at all, or whether the accretion flow on small scales
remains in some sort of quasi-spherical isothermal collapse) – necessarily
requires much higher-resolution idealized simulations which resolve the PSD
scales in detail given the conditions predicted by our simulations here.

Therefore we bracket the range of interesting behaviors by re-running our
fiducial simulation (restarting just before the first resolved STARFORGE sink
actually appears) with identical physics except assuming $t_{\rm
dep}\equiv{\rm MIN}[\,50\,{\rm yr}\ ,\ 100\,t_{\rm dyn,s}^{\rm max}]$ where
$t_{\rm dyn,s}^{\rm max}\equiv t_{\rm tidal}\equiv(\|{\bf T}\|^{2}/6)^{-1/4}$.
Here $\|{\bf T}\|$ is the Frobenius norm of the gravitational tidal tensor
${\bf T}\equiv\nabla{\bf g}$ and the normalization is defined such that in a
Keplerian potential (as we approximately have here with the SMBH dominant),
$(\|{\bf T}\|^{2}/6)^{-1/4}\approx 1/\Omega\approx(G\,M_{\rm
BH}/R^{3})^{-1/2}$. As discussed below $t_{\rm tidal}$ is roughly the maximum
possible dynamical time of the un-resolved PSD/accretion flow. We denote this
as the “fast” accretion regime.

We have also experimented (re-running for a short period of time) with models
using even shorter/faster $t_{\rm dep}$ (or even longer/slower than our “slow”
model), as well as a couple of intermediate cases. We find that on the scales
resolved in our simulations the key behaviors split into just two regimes,
essentially depending on whether we assume the un-resolved accretion occurs
rapidly (“fast”) compared to the global QAD/CQM dynamical time, or slowly
(“slow”), but are insensitive to the details of $t_{\rm dep}$ within each
regime. We therefore treat these two simulations as representative of these
two broad limiting cases.

Figure 3: Effective resolution (§ 2) of the simulation as a function of radial
distance from the SMBH ($r$). We show the minimum (dotted black), median
(solid), $50\%$ and $90\%$ inclusion intervals (shaded) of the mass resolution
(gas cell mass $\Delta m$; left), spatial resolution (equivalent cell size
$\Delta x\equiv(\Delta m/\rho)^{1/3}$; center), and time resolution (numerical
timestep $\Delta t$; right), at each $r$. For reference, we compare the total
enclosed gas mass ($M_{\rm gas}(<r)$), radius ($r$), dynamical time ($t_{\rm
dyn}\equiv\Omega^{-1}=r/V_{c}(r)$), and how long the simulation was run after
reaching maximum refinement at each radius (“duration”). We also show the
resolution test with additional refinement in the QAD from Paper II (dotted
magenta), run for much shorter duration. We denote where the simulation treats
star formation using the un-resolved FIRE “single stellar population particle”
(SSP) approach (sampling an assumed IMF), versus STARFORGE individual star
particles. The scales of interest here (interior to the BHROI at a few pc) are
uniformly resolved at the target $\Delta m\sim 0.001-0.01\,M_{\odot}$, with
spatial (time) resolution reaching $\sim 10$ au (hours). The duration at
highest resolution $\sim 10^{4}$ yr corresponds to $\sim 10^{5}$ orbital
dynamical times of the innermost QAD. Note the external/galactic tidal field
$\sim G\,M_{\rm enc}/r^{3}\sim\Omega^{2}$ can be inferred here: it is $\sim
10^{17}$ times stronger in the inner QAD compared to the LISM.  Figure 4:
Characteristic mass, velocity, and inflow/outflow rate scales in the
simulations (§ 3). Left: Circular velocity profile $V_{\rm
c}\equiv\sqrt{G\,M_{\rm enc}(<r)/r}$ versus distance from the SMBH $r$, with
contributions from different mass components. We focus on radii where the
STARFORGE physics applies. Note the SMBH excision radius truncates the gas
mass distribution at $\lesssim 0.001\,$pc, and the contribution from dark
matter (while included) is negligible at these radii. We see a clear
transition from a stellar-mass dominated galactic nucleus through the BHROI at
a few pc, with the stellar density rapidly shut-down at sub-pc scales as SF is
suppressed overall in the QAD (Paper I). Note the characteristic velocity
scales $\gg 100\,{\rm km\,s^{-1}}$ in the QAD/CQM. Right: Mass flow rate
$\dot{M}$, showing the total inflow rate $\dot{M}_{\rm in}$ and total outflow
rate $\dot{M}_{\rm out}$ through each annulus, and the cumulative SFR
(averaged over the previous $t_{\rm dyn}$) summed within each radius
$\dot{M}_{\ast}(<r)$. Lack of spherical symmetry and non-equilibrium dynamics
mean that inflow and outflow co-exist and can change relative sign, but
despite this a remarkably stable order-of-magnitude $\dot{M}_{\rm in}\sim
10-100\,{\rm M_{\odot}\,yr^{-1}}$ persists (see Paper II). We see the
suppression of SF relative to inflow explicitly on small scales.  Figure 5:
Properties of the gas in the QAD/CQM versus radius (§ 3). In each panel, we
compare the simulations to the analytic prediction for a standard Shakura &
Sunyaev (1973, SS73) $\alpha$-disk QAD ($\alpha=0.1$) extrapolated to the same
radii. Top Left: Projected gas and stellar surface density $\Sigma$ in
cylindrical shells. We see the transitions from Fig. 4. The surface densities
are uniformly $\gtrsim 10^{4}\,{\rm M_{\odot}\,pc^{-2}}$, much higher than
LISM GMCs, and the entire QAD is optically thick to its own cooling radiation.
Top Right: Gas density $\rho$ versus radius $r$. We show volume and mass-
weighted means (line) and $90\%$ inclusion intervals (shaded; mean can exceed
this if the distribution has large “tails”). Middle Left: Temperatures: we
plot volume and mass-weighted means and $90\%$ inclusion radii for gas
temperature, effective radiation temperature (integrating over all RHD bands
evolved), and dust temperature (dynamically evolved). Owing to optical-depth
and heating, all temperatures exceed the CMB temperature ($\approx 15\,$K at
this redshift) at radii here, with the temperatures in the QAD approaching
blackbody equilibrium. Differences in opacities and surface densities explains
the difference with the SS73 model. The dust sublimates in-code at $T_{\rm
dust}\gtrsim 1500\,$K. Middle Right: Magnetic field strengths: energy-weighted
field strength and radial/toroidal/poloidal components. In and outside the CQM
fields are crudely isotropic with a mild radial bias, but the mean toroidal
component comes to dominate in the QAD (Paper II). Bottom Left: Toomre $Q$
parameter accounting for thermal ($Q_{\rm therm}$), magnetic ($Q_{\rm mag}$),
turbulent ($Q_{\rm turb}$), or combined support of the gas. The outer CQM has
turbulent $Q\sim 1$ with thermal $Q<1$, but even the thermal-only $Q\gg 1$ in
the QAD. Bottom Right: Sonic ($\mathcal{M}_{s}\equiv\delta v_{\rm
turb}/c_{s}$) and Alfvénic ($\mathcal{M}_{A}\equiv\delta v_{\rm turb}/v_{A}$)
Mach numbers (mass-weighted), for each component of the random motions in the
midplane. The velocity dispersions are broadly isotropic and highly
supersonic, with a gradual trend from super-Alfvénic to sub-Alfvénic at
decreasing $r$. While $\Sigma$, $\rho$, and $|{\bf B}|$ are vastly larger than
in the LISM, the QAD is vastly lower-density compared to the SS73 model, owing
to strong magnetic support (Paper III; note even though $|{\bf B}|$ is lower
than in SS73, $v_{A}$ and $v_{A}/c_{s}$ are orders-of-magnitude larger, owing
to the different densities). This means $Q$ is larger here by a factor
$\gtrsim 10^{5}$ and $\mathcal{M}_{s}$ larger by a factor $\gtrsim 100$.

##  3 Summary of Key Differences Between the CQM/QAD and the “Normal” ISM and
$\alpha$-Disk QAD Models

Before going forward, it is useful to review and differentiate the extreme
conditions here – on sub-pc scales around an accreting quasar – from even the
ISM of galactic nuclei (typically scales $\lesssim$ kpc), let alone the Solar
neighborhood LISM (the focus of the vast majority of IMF literature). To aid
in this we refer to the “diffuse”/volume-filling gas phases in the high-
resolution region of our simulations interior to the BHROI at $\lesssim 5\,$pc
scales as the CQM, and the material circularized into an accretion disk around
the SMBH on $\lesssim 0.1$ pc scales as the QAD.

The properties discussed below are all studied in detail in Paper I and Paper
II, to which we refer interested readers for details. Our intent here is to
simply provide context for the results below. A visual illustration is
provided in Figs. 1 & 2. Some zeroth-order quantitative properties of the
medium are presented in Figs. 3, 4, & 5 (similar properties, but with more
focus on larger radii, are shown in Paper I).

Some useful “global” numbers including the enclosed gas mass and dynamical
times $t_{\rm dyn}\equiv r/V_{c}(r)$ as a function of radius $r$ are presented
in Fig. 3, alongside the relevant simulation resolution. The inflow rates,
outflow rates, and star formation rates averaged over the last dynamical time
of the simulation are shown in Fig. 4, along with the relevant contributions
to $V_{c}(r)$ or equivalently $\Omega\equiv V_{c}/r$ or the global tidal field
$\sim\Omega^{2}$. As is immediately obvious from that Figure, we are focused
in this paper exclusively on scales interior to the BH radius of influence
(BHROI; defined as $R_{\rm BHROI}\sim GM_{\rm BH}/\sigma_{\rm gal}^{2}$ in a
galactic potential with velocity dispersion $\sigma_{\rm gal}$) at several pc,
interior to which (by definition) the BH dominates the global gravitational
forces. As shown in Paper II, properties such as the inflow rates in Fig. 4
are time-steady over the entire duration of our simulation $\sim
10^{5}\,t_{\rm dyn}^{\rm inner}$.

For reference, in Fig. 5 we also compare the simulation properties to the
scalings predicted by the traditional Shakura & Sunyaev (1973) (henceforth
SS73) $\alpha$-disk QAD model, assuming the typical $\alpha=0.1$, with the
same $\dot{M}$ as the simulations from Fig. 4. This comparison, and the
fundamental differences in accretion disk properties between the magnetically-
dominated disks seen in the simulations and SS73-like models (which
fundamentally assume relatively weak magnetic support, with Alfvén speeds much
smaller than thermal sound speeds), are discussed at length in Paper II and
Hopkins et al. (2023c) (from which the scalings in Fig. 5 are taken,
henceforth Paper III). We refer to those papers for discussion of the physics,
but we show this because a number of recent analytic studies assuming
properties of stars forming in QADs (see references in § 1) have assumed the
QAD to have properties given by an SS73-like model, which Fig. 5 shows can be
orders-of-magnitude different from the conditions here, in ways that have
important consequences for star formation.

###  3.1 High (Surface and 3D) Densities and Optical Depths

Per Fig. 5, on sub-pc scales, the surface density of gas in the QAD/CQM
(averaged in cylindrical or spherical shells) scales approximately as
$\langle\Sigma_{\rm gas}\rangle\sim 2\times 10^{5}\,({\rm pc}/R)\,{\rm
M_{\odot}\,{\rm pc}^{-2}}$ as a function of BH-centric radius $R$, and the
mean 3D gas density in the star-forming “midplane” scales as $\sim 3\times
10^{6}\,{\rm cm^{-3}}\,({\rm pc}/R)^{2}$ (in a volume-weighted sense, or about
$\sim 30-100$ times larger in a mass-weighted sense). The related “total
surface density” $\Sigma_{\rm eff}\sim M_{\rm enc}(<R)/\pi\,R^{2}$ or
effective acceleration scale $a_{\rm eff}\sim G\,M_{\rm
enc}/R^{2}\sim\pi\,G\,\Sigma_{\rm eff}\sim V_{c}^{2}/R$, includes stars and
the SMBH so is even larger at $\sim 2\times 10^{6}\,{\rm
M_{\odot}\,pc^{-2}}\,(1+{\rm pc^{2}}/R^{2})$, orders-of-magnitude larger than
the critical value where simulations generically find that stellar feedback is
unable to suppress SF, drive strong winds, or have other appreciable influence
on the medium (see Grudić et al., 2018a; Grudić et al., 2019; Grudić et al.,
2020; Hopkins et al., 2022b, and references therein).

Thus the least-dense CQM structure that could be considered an “overdensity”
in any meaningful sense in our simulations is already $\sim 5000$ times higher
column density than a typical Solar-circle GMC, and $\gtrsim 1000$ times more
dense in 3D than a typical protostellar core. This, in turn, means that the
optical depth scales are also large. Even in the least-dense outer regions,
the “diffuse” CQM is opaque to optical/UV/EUV radiation and already marginally
optically thick to IR cooling radiation. As shown in Paper I and noted below
this leads to a different thermal structure of this QAD/CQM, and means that we
need to consider the possibility of optically-thick cooling physics
everywhere.

That said, the densities are also much lower than those of a canonical SS73
disk with the same $\dot{M}$. As discussed at length in Paper II and Paper
III, this owes to the combination of strong magnetic fields and turbulence
“puffing up” the disk (making it much thicker) but also providing strong
Maxwell and Reynolds stresses, so the same $\dot{M}$ corresponds to much lower
surface densities and optical depths. As shown in Paper II, this means that
non-LTE and even multi-phase structure can still be present in the outer QAD,
and the disk is much less strongly radiation-pressure dominated compared to
those models.

###  3.2 Extremely Strong Magnetic Field Scales

The typical magnetic field strengths in the CQM scale as $|{\bf B}|\sim
0.1\,{\rm G}\,({\rm pc}/R)^{4/3}$ (with a factor of a few scatter; Fig. 5). In
other words, we are following star formation at radii where the ambient gas
before collapse into a PSD can already have magnetic field strengths reaching
up to kiloGauss. These are extremely strong in an absolute sense, but also in
a relative sense in terms of the plasma $\beta\equiv c_{s}^{2}/v_{A}^{2}$ or
ratio of magnetic field strength to gravity. As shown in Paper I and discussed
at length in Paper II, $\beta$ is extremely small in these simulations, often
$\sim 10^{-4}$ in the diffuse CQM – much smaller than typical GMCs. This is
also orders-of-magnitude smaller than assumed in the SS73 model (which by
definition assumes $\beta\gg 1$).333As discussed in Paper III, the absolute
value of $|{\bf B}|$ is actually larger in an SS73 disk, because the densities
are so enormous by comparison to the simulations here that producing even a
modest Maxwell stress to explain $\alpha\sim 0.1$ requires extremely large
$|{\bf B}|$. But the Alfvén speeds $v_{A}$ and “magnetic support” of the disk
(e.g. $v_{A}/V_{c}$) are vastly smaller in SS73 compared to the simulations
here).

Moreover, the magnetic critical mass using the median $|{\bf B}|$ is actually
larger than the entire QAD/CQM mass on sub-pc scales (shown explicitly in
Paper I and below in § 4.6.1). Similarly, the magnetic Toomre $Q_{\rm mag}$
parameter is $\gtrsim 100$ everywhere we analyze in the QAD and $\gtrsim 1$
throughout the entire CQM outside the QAD (Fig. 5). This means (see Paper I)
that the field is strong enough to ensure the QAD is stable against gravito-
turbulence in the classical sense even where the thermal-only Toomre $Q_{\rm
thermal}$ is modest ($\sim 10$ at the largest QAD radii, though it drops
rapidly at larger $R$ in the CQM). As shown in Paper I, this strongly inhibits
star formation on these smaller (QAD) scales relative to simulations which
simply do not include magnetic fields. And we see that for a “typical”
thermal-pressure supported SS73 disk, the predicted $Q_{\rm SS73}$ is $\sim
10^{5}-10^{7}$ times smaller ($\ll 1$ at all radii simulated here). And
indeed, absent magnetic fields, the QAD fragments catastrophically at $\sim$
pc scales (Paper I & Paper II). So it is clearly important to include them in
our study here, though they do not completely suppress star formation.

As shown in Paper II and visually below, the fields within the QAD are
primarily toroidal/azimuthal, albeit with non-negligible radial and
poloidal/vertical components. This means radial collapse (e.g. Toomre-style
radial modes in the QAD) is especially inhibited. However there can be some
regions of locally lower $|{\bf B}|$, and some regions where the field
switches direction (owing to its being the flux-frozen relics of originally
more isotropically turbulent ISM fields amplified as they fall around the BH),
at which point collapse along field lines is less prohibitive, as we show
below. In the CQM, where the inflow has the form of an infalling filamentary,
tidally-disrupted molecular cloud complex, the magnetic fields are more
isotropically turbulent/disordered (see Paper II for details), with a mild
radial bias, more akin to LISM molecular clouds in geometry.

###  3.3 Warm Thermal Dust, Gas, & Radiation Temperatures

Turning to the thermal properties, the diffuse QAD and CQM do not at first
glance appear so unusual – most of the QAD mass is in a WNM-like phase of
atomic gas with temperatures of a few thousand Kelvin and free electron
fractions $\sim 1-10\%$, and metallicities close to ($\sim 1/2$) Solar (see
Paper I, and Figs. 2 & 5). As shown in Paper II and Paper III, the cooling
times are short compared to dynamical times ($t_{\rm cool}\ll t_{\rm dyn}$)
throughout the CQM and QAD, like in the LISM. However, given the extremely
high densities and optical depths noted above, the dust and radiation
temperatures and gas temperature of the midplane phases are tightly coupled,
and scale as approximately $T_{\rm rad}\sim T_{\rm dust}\sim 300\,{\rm
K}\,({\rm pc}/R)^{1/2}$ (Fig. 5). There is some gas cooler than this but not
much; and note that since this is at $z\approx 4.4$, the temperatures
essentially never get cooler than the CMB at $\sim 15\,$K. Thus there is
little “cold” phase in the classic ISM sense. Indeed, even much of the
molecular gas in the diffuse QAD (where present) is actually quite warm at
$\gtrsim 1000\,$K, and a factor of a few colder in the outer CQM.

The warm temperatures owe to a combination of effects: optically-thick cooling
as mentioned above, but also the fact that this is occurring around a quasar
in the center of a high-redshift massive starburst galaxy. The star formation
rate inside of the central couple hundred pc is $\gtrsim 100\,{\rm
M_{\odot}\,yr^{-1}}$, completely typical for a high-redshift extremely
luminous quasar but obviously vastly larger than the Milky Way. Moreover, the
accretion luminosity of the quasar itself is enormous and is balancing the
local accretion rate through each annulus, $\sigma_{\rm B}\,T_{\rm
eff}^{4}\sim\dot{M}\,\Omega^{2}$, naturally giving a much larger radiation
density. Indeed the interstellar radiation field or ISRF in the CQM here has
an energy density $\sim 10^{7}\,({\rm pc}/R)^{2}\,{\rm eV\,cm^{-3}}$. The
estimated cosmic ray energy density is also similarly elevated (with both
protons and electrons well into the calorimetric limit by any reasonable
estimation).

Moreover, because of its coupling to these large temperatures, the dust begins
to sublimate in the QAD interior to $R\lesssim 0.1\,$pc. So in the inner
regions we study, the medium has only gas-phase metals, and the cooling
becomes again quite different.

Note that the thermal structure is not so dramatically different from SS73:
the differences primarily owe to different optical depths and opacity
structure of the disk owing to its different surface densities and NLTE and
multi-phase structure. But in SS73, $t_{\rm cool}\gg t_{\rm dyn}$ by
assumption, even though the effective photospheric disk temperatures are not
so different, owing to the very different optical depths and densities/masses
of the disks.

###  3.4 Large Turbulent Velocities

The turbulence in the QAD/CQM is highly super-sonic (sonic Mach numbers
$\mathcal{M}_{s}\sim 30-300$, related physically to the rapid cooling noted
above) and trans-Alfvénic (Alfvén Mach numbers $\sim 0.3-3$), shown in Fig. 5.
In these dimensionless terms this is not so different from some cold ISM
phases. However recalling the large temperature and magnetic field scales,
this translates to turbulent velocities at the driving scale (of order the QAD
scale height $H\sim(0.1-0.3)\,R$) of $\gtrsim 200\,{\rm km\,s^{-1}}$ at all
radii we study, rising to as large as $\sim 1000\,{\rm km\,s^{-1}}$ in the
innermost regions at $\lesssim 0.001\,{\rm pc}$. So there can be extremely
large post-shock temperatures and compression ratios. Again this differs
dramatically from an SS73 disk, where (by assumption) the turbulence is always
subsonic.

###  3.5 Extremely Strong Tidal Forces

Of special importance in a circum-BH environment, the tidal field is extremely
strong (Fig. 4). Since the large-scale potential is dominated by the SMBH on
the scales we study, the tidal tensor scales as $|{\bf T}|\equiv|\nabla{\bf
g}|\equiv|\nabla\otimes\nabla\Phi|\sim G\,M_{\rm BH}/R^{3}\sim\Omega^{2}$ in
terms of the CQM/QAD dynamical frequency $\Omega\equiv V_{\rm c}/R\sim
1/t_{\rm dyn}$ (where $t_{\rm dyn}$ is the CQM/QAD dynamical time). Using the
simulation $M_{\rm BH}\approx 1.3\times 10^{7}\,{\rm M_{\odot}}$, we see that
the dynamical time $1/\Omega\sim 4000\,{\rm yr}\,(R/{\rm pc})^{3/2}$, so it
reaches as short as $\sim 1\,$month at $R\sim 0.001\,{\rm pc}$ – compare this
to $1/\Omega\sim 40\,{\rm Myr}$ in the Solar neighborhood. This in turn means
that the tidal field $\sim\Omega^{2}$ is stronger than that in the Solar
neighborhood by a factor of $\sim 10^{8}\,({\rm pc}/R)^{-3}$, i.e. $\sim
10^{8}-10^{17}$ times stronger at the radii we consider.

This tidal field has immediate consequences for the structure of PSDs/PPDs in
the CQM. Consider a PSD with system mass $M_{\rm s}\equiv M_{\ast}+M_{\rm
PSD}$, and assume its potential is quasi-Keplerian for simplicity. At a
distance $R$ from the SMBH, the PSD has a tidal radius $r_{\rm t}=R\,(M_{\rm
s}/4\,M_{\rm BH})^{1/3}\sim 500\,{\rm au}\,(M_{\rm s}/{\rm
M_{\odot}})^{1/3}\,(R/{\rm pc})$, or $G\,M_{\rm s}/r_{\rm t}^{3}\leq
4\,\Omega^{2}$. Outside of this radius it will be rapidly tidally disrupted.
Equivalently, the maximum PSD dynamical time $t_{\rm dyn,\,s}$ or minimum
$\Omega_{\rm s}$ is $\Omega_{\rm s}^{\rm min}=1/t_{\rm dyn,\,s}^{\rm
max}\equiv\sqrt{G\,M_{\rm s}/r_{\rm t}^{3}}=2\,\Omega$, so $t_{\rm
dyn,\,s}^{\rm max}\approx 2000\,{\rm yr}\,(R/{\rm pc})^{3/2}$ ($\sim 20\,$days
at the smallest $R$ here).

As noted above, these extreme (short) timescales force us to consider
different prescriptions for the depletion timescale $t_{\rm dep}$ from the un-
resolved PSDs onto (proto)stars. But they also have important implications for
PSD structure and stability and stellar evolution, which we will discuss in
more detail below. And the extreme tidal fields can, of course, be prohibitive
of star formation and accretion under many circumstances.

##  4 Results & Discussion

Having described the simulation methods (§ 2) and emergent properties of the
CQM (§ 3) in which we will study SF, we now present the results for CQM SF.

Figure 6: Image illustrating the locations of star formation in the CQM/QAD (§
4.1). We plot the projected density of gas (face-on to the inner QAD) akin to
Fig. 1, and superimpose locations of individual (proto)stars (points). Stars
are colored by “dynamical age”: time since formation $t_{\ast}<t_{\rm dyn}$ at
their current radius (dark blue), or $<10\,t_{\rm dyn}$ (light blue) or older
(white). Inset images highly two sub-regions. SF primarily occurs in two
regions here: (1) the free-falling, tidally compressed filament fueling the
QAD (tidally stripped from larger-scale massive cloud-complexes) in the CQM,
and (2) in the outer QAD. The filamentary regions are analogous to LISM-like
dense filamentary star-forming regions, though they exhibit structure on much
smaller absolute physical owing to the extreme densities and tidal fields
(individual “cores” and shock sub-structures have scales $<0.001$ pc, as
compared to $\sim 0.1$ pc).  Figure 7: As Fig. 6, but focusing on the QAD
(projected face-on and edge-on). SF occurs around where gas circularizes and
is stabilized by magnetic fields but is still appreciably gravito-turbulent
with a modest thermal-only Toomre $Q_{\rm thermal}\sim 1-10$ (Fig. 5). In the
inner QAD at $\ll 0.01\,$pc, where even the thermal-only $Q_{\rm thermal}$
rises to $\gg 10^{2}-10^{3}$, star formation becomes nearly impossible: all
stars at these radii are dynamically “old,” and formed at larger radii on
highly-eccentric orbits.  Figure 8: Images highlighting the multi-phase
structure of the gas over the same range of scales as Figs. 6-7. We plot
projected density edge-on, colored by gas temperature: $\ll 10^{4}\,$K
(magenta), $\sim 10^{4}\,$K (green), and $\gg 10^{4}$ K (red). Warm/hot, low
density gas is present in the CQM from shocks and supernovae explosions, and
vents from the inner QAD in a jet-like structure. The outer CQM also exhibits
colder “clouds” embedded in warm gas, but the tidally-disrupted GMC-like
complex is much more visible as the extended cold “stream,” with some material
free-falling onto the QAD, until it is heated via accretion and inefficient
radiative cooling to temperatures $\gtrsim 1000\,$K.

###  4.1 Where do stars form in the QAD & CQM?

Having shown the zeroth-order properties and large-scale structure of the
medium, Fig. 6 “zooms in” to show the gas and stars in the CQM, and Fig. 7
does the same in the QAD, at a time near the end of the evolved simulation
duration (well after it reaches a quasi-steady-state in the QAD). We highlight
some key features (some of which are also clear in the gas-only morphology
shown in Figs. 1 & 2). First, despite the large magnetic Toomre $Q$ parameter
(see Paper I), the gas exhibits significant inhomogeneity, although it becomes
“more smooth” on smaller scales. Fig. 8 shows that this is reflected as well
in the gas phase structure. Second, the positions of the (proto)stars do not
appear particularly correlated with the gas structure. We discuss this in more
detail below in § 4.3.1, but it fundamentally relates to the fact that the
dynamical times at these radii are vastly shorter than even proto-stellar
evolution times (§ 3), so stars can migrate far from their “birth” locations
in timescales as short as years, in the inner QAD.

If we wish to understand “where star formation occurs,” it is more useful to
look at stars which have just formed. Highlighting the youngest stars/sinks in
Figs. 6-7, we see that there are two qualitatively distinct “sites” of star
formation.

####  4.1.1 Star Formation in Infalling CQM Filaments

At radii between $\sim 0.1-10$ pc, star formation is largely occurring within
relatively large-scale filamentary structures composed of gas captured by the
BH, falling onto it in tidal streams, highlighted in Figs. 6 & 8. On these
scales there is no QAD to speak of: as discussed in Paper I, the galaxy is a
chaotic, highly turbulent, clumpy, massive, gas-rich system undergoing a major
merger at $z\sim 4.4$, and the quasar event here is related to a close passage
(induced by the merger) of one extremely massive molecular cloud complex (gas
mass $>10^{8}\,M_{\odot}$) by the galactic nucleus containing a SMBH with mass
$\sim 10^{7}\,M_{\odot}$. Some material is tidally stripped off of the cloud
and captured in this event, so it initially free-falls onto the BH from the BH
radius of influence (where the BH begins to dominate the potential) at a few
pc. This forms the “parent” filament.

This infalling filament is already part of a massive star-forming complex when
it is captured from the ISM, so it is not surprising that star formation
continues. As discussed in Paper II, as it falls in, it is tidally stretched
in the radial/infall direction while being compressed in both perpendicular
directions, so it becomes more visibly elongated/filamentary but the local 3D
density of the dense clumps within it is not suppressed and is actually net
tidally enhanced/compressed. More formally, it already has a sufficiently low
virial/Toomre parameter ($Q_{\rm eff}\sim 1$ even including turbulent+magnetic
support at $r\gtrsim 1\,$pc, with thermal-only $Q_{\rm thermal}\ll 1$; see
Fig. 5), and sufficiently rapid cooling ($t_{\rm cool}\ll t_{\rm dyn}$), and
relatively low optical depth to its own cooling radiation, such that we should
expect fragmentation. Likewise while the magnetic fields are strong, the
turbulence is still trans-to-super Alfvénic (Fig. 5), and the magnetic field
geometry is not strongly toroidal on these scales, but more quasi-isotropic
(meaning the fields will be less efficient at preventing collapse).

In this regime, therefore, the star formation is at least “expected,” and
reasonably analogous to star formation in dense filamentary structures in the
LISM (see e.g. Pineda et al., 2023, for a review). We stress that absolute
values of most parameters here – densities, magnetic field strengths,
turbulent velocities, etc. – are orders-of-magnitude different from the LISM
(§ 3), so the problem is rescaled and not exactly analogous. For example, (1)
the high surface density/acceleration scales mean we should be in the regime
where stellar feedback has a weak effect on the dynamics (Grudić et al.,
2018a) and we confirm this below; (2) the high optical depths and strong
radiation background lead to $T_{\rm dust}\sim T_{\rm rad}\sim T_{\rm gas}\sim
100-1000$ K; (3) the absolute size scales of structures are much smaller, with
for example the sonic scale of the turbulence (scale below which rms turbulent
velocity fluctuations become sub-sonic, and density fluctuations become small,
approximately the $\sim H/\mathcal{M}_{s}^{2}$ in a supersonically turbulent
disk) and resolved well-separated density inhomogeneities/structures extending
to $<0.001\,$pc.

####  4.1.2 Star Formation in the Rotating QAD

At radii $\ll 0.1\,$pc, the focus of Fig. 7, the infalling gas circularizes
and forms a QAD that persists in to the smallest radii we can resolve ($\ll
0.001\,$pc, with no obvious reason it could not continue down to scales of
order the horizon around the SMBH; see Paper III). As discussed at length in
Paper I and Paper II and briefly reviewed above, star formation becomes
strongly suppressed at small radii below these scales, for several reasons.
Strong magnetic fields produce rapid accretion through an annulus while
providing a “magnetic” $Q\gg 1$ or magnetic critical mass much larger than the
entire QAD mass, and even the thermal-only $Q$ rises to $\sim 10-30$ on scales
down to $R\sim 0.005\,$pc (before rising very rapidly to even much larger
values at still smaller $R$), the turbulence becomes sub-Alfvénic, the
temperatures of the dust and radiation and “cold/warm” phases continue to rise
to $\gtrsim 10^{3}-10^{4}\,$K and most of the dust sublimates inside of radii
$\lesssim 0.01\,$pc (see Paper I). As a result, the global SFR integrated
within the QAD+CQM drops precipitously and at all radial annuli $\ll$ pc it is
more than an order of magnitude smaller than inflow/accretion rates onto the
SMBH, and the CQM becomes a gas-dominated QAD in the usual sense (Paper I,
Paper II). However, this does not mean there is zero star formation, just that
in e.g. Paper II where we study the properties of the QAD, we can neglect star
formation and stellar feedback and stellar dynamics as significant
perturbations to the gas dynamics/inflow/accretion physics.

Indeed, we see here that there are still some dense star-forming structures in
this QAD zone. They tend to occur at larger radii within the QAD:
unsurprisingly, there is essentially no detectable star formation at the
smallest radii $\lesssim 0.01\,$pc, where accretion maintains the QAD at
temperatures $>10^{4}$ K with a thermal-only Toomre $Q$ parameter of $\gtrsim
1000$ (Fig. 5). Stars might still dynamically make their way into these radii,
however, as we discuss below. Star formation occurs primarily in the disk
midplane as expected, and in radius it occurs vaguely around spiral-arm like
structures which are clearly “shearing out” in the outer disk, but we show
below that in detail the picture is more complicated. It is clear that stars
form even within these structures in some preferred locations of the QAD (the
youngest stars are not uniformly distributed with radius or even azimuthal
angle). We discuss this further below.

Figure 9: Comparison of where stars form in the QAD to the magnetic field
structure. Left: Magnetic field lines (in-plane projection, so
$B_{\phi}-B_{R}$), overlaid on a colormap showing the sign of the toroidal
magnetic field ($B_{\phi}>0$ points in the direction of rotation), on
different scales (top-to-bottom). To emphasize where stars form we show only
dynamically young stars ($t_{\ast}<t_{\rm dyn}(r)$). Right: Same, but overlaid
on the gas density projection from Fig. 7. While there is some density sub-
structure throughout the QAD, we see a striking correspondence between where
stars form and “switches” in the mean toroidal field sign. Per Paper II, these
“switches” arise from flux-freezing of turbulent ISM fields which are
stretched into the dominant toroidal component and advect inwards with the
gas. The strong coherent toroidal field strongly suppresses gravitational
collapse in the vertical/radial directions, while azimuthal modes are sheared
out by the extreme differential rotation, so together this prevents star
formation throughout most of the QAD. However, at the “switches,” the field is
both locally weaker and oriented primarily in the radial direction, so radial
(Toomre-like) collapse becomes possible.  Figure 10: Comparison of where stars
form in the QAD to the velocity field structure; we show young stars and
density/field line maps as Fig. 9 but color code by radial velocity $v_{r}$.
If we plot the velocity stream lines directly, no small-scale structure is
visible, as the velocity is totally dominated by rotation at these scales, so
instead we plot the streamlines of the velocity after subtracting the expected
mean circular velocity: $\Delta{\bf v}\equiv{\bf v}-V_{c}(r)\,\hat{\phi}$. To
leading order the large-scale residual and $v_{r}<0$, $>0$ sign pattern is
dominated by the global eccentricity of the QAD (visible in the density map),
though in the inner QAD the velocity perturbations from individual spiral arms
becomes visible. In either case, there is no correlation between these large-
scale velocity modes and star formation (indeed SF appears in especially
smooth regions of the velocity flow), so SF is not primarily driven by spiral
modes in the QAD, but rather by the magnetic structure in Fig. 9.  Figure 11:
Magnetic (top) and (rotation-subtracted $\Delta{\bf v}$) velocity (bottom)
field line structure as Figs. 9-10 overlaid on gas density in the inner QAD
where there is no SF in the last $\sim t_{\rm dyn}$ here. The coherence of
$B_{\phi}$ is visible, and prevents the spiral overdensities (which clearly
correlate with the $\Delta{\bf v}$ field) from non-linearly collapsing.
Figure 12: Magnetic and velocity fields face-on in the QAD as Fig. 9-10, but
in an edge-on projection in cylindrical coordinates ($\hat{R}$-$\hat{z}$
plane, selecting a wedge with $|\sin{\phi}|<0.3$). The field lines are again
in the projected plane ($B_{R}-B_{z}$, $v_{R}-v_{z}$). The dynamically-young
stars clearly form in the denser midplane. QAD is thick and turbulent, with
largest eddies of size $\sim H\sim 0.1\,R$, which are trans-Alfvénic but
highly super-sonic so generate strong density enhancements in the shocks seen,
but these are rapidly sheared out into the diagonal “fans” seen here as they
form.  Figure 13: Edge-on projection of the QAD on larger scales compared to
Fig. 12, now in a Cartesian projection (not taking a wedge in cylindrical
coordinates), and showing the dynamically “older” stars (white stars at
left/top), as well as all gas cells which are primarily composed of
(proto)stellar ejects (winds and jets from material accreted by sink
particles; cyan points at right). This allows us to better see that while the
youngest stars form in the midplane, the older stars and their ejecta are out-
of-midplane with an asymmetric skew “above” the QAD. This is real and owes to
the continuous precession of both the QAD eccentricity axis and the QAD plane
as the angular momentum vector of the “accretion stream” (visible at top and
Figs. 6) wraps around the QAD at a continuously varying angle (as expected,
since the disrupted cloud/stream is still moving), on a timescale of order the
dynamical time at the outer QAD radius ($\sim 100\,$yr). The dynamically older
stars and their ejecta reflect the previous QAD orbits.

###  4.2 Overcoming Support to Form Stars: The Critical Role of Toroidal
Magnetic Fields and Tidal Forces

As discussed above (§ 4.1.1), in the “infall/filament” CQM region, just like
in a typical GMC, there is not really a dramatic “barrier” to fragmentation
that needs to be overcome. As such, the formation of bound/collapsing sub-
clumps can occur via much the same processes: trans/super-Alfvénic, highly
super-sonic turbulence creates a spectrum of dense sub-structures that can be
internally Jeans/Toomre unstable and collapse/fragment (see e.g. Hopkins,
2013b, and references therein for a theoretical review).

In the disky, rotating QAD region, per § 4.1.2, there is such a barrier, which
prevents most of the gas from fragmenting. While there are many contributing
factors, as reviewed above and shown in detail in Paper I, the most important
factor is the combination of strong toroidal magnetic fields in a strong tidal
environment. This is illustrated in Figs. 9 & 10 which plot the face-on
magnetic and velocity field lines in the QAD alongside the gas density field
and locations of the most-recently-formed stars, with Figs. 11 & 12 showing
the same for the inner QAD and in edge-on projection, respectively.

Without magnetic fields, much of the QAD outside of $\gg 0.001\,$pc would
easily meet the conditions for efficient gravito-turbulent fragmentation
(Paardekooper, 2012; Meru & Bate, 2012; Hopkins & Christiansen, 2013; Deng et
al., 2017): it features a cooling time short compared to the dynamical time
($t_{\rm cool}\sim\beta\,t_{\rm dyn}\sim(10^{-5}-10^{-2})\,t_{\rm dyn}$; see
Paper II), supersonic turbulence, and the thermal-only Toomre $Q\sim 10$ is
modest. And indeed we confirmed this in Paper I, where we compared simulations
with magnetic fields artificially removed, which experienced runaway
fragmentation and star formation on these scales. A strong mean toroidal field
directly prevents collapse in the vertical and radial directions (i.e.
suppresses axisymmetric/Toomre-like modes); collapse would in principle be
possible in the azimuthal direction but any modes/structures with finite
radial extent or radial velocities are rapidly sheared apart in this direction
owing to the strong tidal field. For the typical field strengths here, this
would be sufficient to totally suppress gravito-turbulent collapse if e.g. the
field were idealized as a coherent, cylindrical pure-azimuthal mean field
everywhere (Lizano et al., 2010; Lin, 2014; Riols & Latter, 2016; Forgan et
al., 2017). The strong field also enables strong torques which accelerate
accretion and prevent mass from “piling up” in the QAD (Paper II).

However, although there is a strong toroidal preference, the magnetic fields
here are not perfectly cylindrical, uniform, azimuthal structures. In Paper II
and Fig. 12 we show for example that there is also a clear turbulent field
within the QAD midplane with fluctuating components not too much smaller than
the mean toroidal field (see Fig. 5), but this generally does not change its
qualitative ability to resist collapse (the mean field is still strong and
primarily toroidal). More relevant, at certain radial annuli in the QAD (which
move inwards with the accreting gas as it advects), the mean toroidal field
experiences coherent sign flips (Fig. 9, see Paper II for details). This can
occur semi-stochastically if the field is driven by instabilities like the MRI
(Kudoh et al., 2020), but we show in Paper II that here it owes to an even
simpler explanation: the mean field “dynamo” is powered primarily by advection
of magnetic flux as material accretes from larger radii. As such, initially
tangled fields from ISM scales (where the turbulence is super-Alfvénic) are
tidally stretched and shorn into a radial (at large BH-centric radii) then
toroidal (inside the QAD) configurations with the sign flips and “zones” of
weaker field reflecting the initial small-scale ISM turbulent magnetic field
conditions.

We can see directly in Fig. 9 that star formation in the QAD is clearly
strongly associated with these “flips” in the toroidal field. The young sinks
are not homogeneously distributed throughout the QAD, and they do not
correspond to obvious global structures in the velocity field in Fig. 10 or
Fig. 11. There we plot the residual velocity field lines, after subtracting
the mean rotational motion, to highlight deviations from circular orbits,
which are strongly dominated by the global coherent eccentric structure of the
disk and the spiral arms, but these structures do not actually appear to play
the key role in local SF in the QAD. Likewise, in Fig. 12, while we see the
stars form vaguely around the midplane/width of the disk (with $|z|\lesssim
H\sim v_{A}/\Omega$), there is no sharp association with some structure in
either the $v_{R}-v_{z}$ or $B_{R}-B_{z}$ fields, as there is for the
$B_{\phi}$ field.

In Fig. 9, we see that where the mean field is coherent, fragmentation is
almost totally suppressed, as we expect from the arguments above. But near the
sign flips, there must be regions where the toroidal field vanishes entirely
(the fields become primarily radial, in this case), and we see gas collapse
radially along the mean field lines to form structures, just like we would
expect for simple axisymmetric Toomre instabilities in a gravitoturbulent QAD.
In these regions, the field strength is also somewhat weaker, but this is a
smaller effect (the magnetic critical mass, using just $|{\bf B}|$, is often
still large compared to the condensing gas mass), so the more important effect
is the local field geometry, which when non-toroidal means that gas can
collapse along field lines in the strongly-magnetized QAD (and thus not feel
the strong magnetic pressure which would otherwise resist collapse). In
contrast, in the CQM at $r\gtrsim 1\,$pc, where the field is more
isotropically turbulent and the gas is not in an ordered disk, there are
essentially “switches” located throughout, so there are always directions for
collapse to proceed even if the local magnetic field value is strong/magnetic
critical mass is large.

Note that while we see essentially all SF in the QAD associated with switches
in the polarity of the toroidal field, the converse is not true: not every
polarity switch produces star formation. In the inner disk, even the thermal-
only $Q$ is much too large. The SF is most vigorous in the outer disk “flips,”
where $Q_{\rm thermal}$ is relatively low, and the most intense episodes
(evident in Fig. 9) tend to coincide with the overlap of a polarity switch and
spiral overdensity in the outermost QAD.

Figure 14: Growth of massive stars. We select the $\sim 50$ most massive stars
which form at radii $0.02<r<0.05\,$pc in the “slow” PSD accretion simulation,
and plot total mass $M_{\rm s}$ versus time since the sink particle was
created in the simulation ($\Delta t\equiv t-t_{\rm sink}$). Linewidths are
ordered by the time when the star reaches $1/2$ of its final mass. Since this
simulation has weak (proto)stellar accretion/feedback and more massive stars
compared to the “fast” PSD accretion run, this represents an “upper limit” to
the duration of star formation. We see that almost all of the growth occurs in
an extremely short period of time (§ 4.3), with the rapid rise in $M_{\rm s}$
occurring on of order the dynamical time at this radius $t_{\rm
dyn}=\Omega^{-1}(r)\sim 50\,$yr (median time to $1/2$ or $1/6$ of the final
mass is $\sim 3$ or $\sim 1\,t_{\rm dyn}$). Even those stars with “late”
growth often grow via a delayed secondary “burst” with duration $\lesssim
t_{\rm dyn}$. In this strong tidal environment, growth is limited by the
dynamical separation of the initial cores/sinks and surrounding gas and the
extremely short dynamical times relative to normal pre-main sequence stellar
evolution timescales.

###  4.3 Timescales of Star Formation and Growth

In Fig. 14, we examine how rapidly collapse and star formation actually occurs
in these environments. Specifically, we can examine how the total mass of a
typical sink particle grows in time at different initial radii around the
SMBH. Recall from § 2.1, sinks form with the mass of their parent gas cell,
$\sim(1-5)\times 10^{-3}\,{\rm M_{\odot}}$ here, then capture only gas which
is bound and has essentially un-resolved orbits around the sink. This is added
to a reservoir representing any un-resolved accretion flow or a PSD (or, if
the PSD fragments, un-resolved multiples and/or planets) which can also
accrete onto the proto-star or be ejected (representing protostellar jets,
magnetic outflows from the PSD, main-sequence mass-loss, etc.). Rather than
focus on the growth of just the “star” itself, which depends on these un-
resolved assumptions about transfer between “reservoir” and “star” labels
within the sink (see § 2.2), we first examine the total sink masses which
represent the total accreted (minus expelled) mass, and represent the resolved
dynamics.

We see the sink masses rise initially extremely rapidly, approaching a maximum
on a timescale as short as years, even at $\sim 1\,$pc from the BH where the
dynamical time $t_{\rm dyn}\sim 1/\Omega\sim 4000\,{\rm yr}\,(R/{\rm
pc})^{3/2}$. Then they flatten or asymptote to a final mass, or even begin to
decline (as e.g. some mass is lost from the systems owing to jets). Note that
this initial rise time being much faster than the local dynamical time should
not be taken too literally and does not actually mean cores form much faster
than $\sim 1/\Omega$ (though they can in principle, e.g. in strong shocks).
Rather, what we often see is that the dense gas condensations (which will form
sinks) form as described in § 4.2, on timescales broadly $\sim 1/\Omega$, but
it is (by construction) only after a clear collapsing core has formed that the
code will (numerically) convert the gas cell at the density maximum (roughly
the “center” of the condensation/core) into a sink, which can then very
rapidly numerically accrete the neighbor cells which were already part of its
bound, resolved collapsing core mass at the time of the initial sink particle
“promotion.”

In Fig. 14, we see that this qualitative behavior is independent of the
initial radius of star formation over the $\ll 10\,$pc scales we follow (i.e.
independent of whether the SF occurs in the “infall/filament” CQM or rotating
QAD region; § 4.1), and independent of the sub-grid prescription for accretion
from the “reservoir” to the protostar (§ 2.2).

The important thing is that after this “initial core” is rapidly accreted and
converted into a sink, there is very little on-going growth/accretion from the
medium, even for the most massive stars formed here. This is very much unlike
the situation for massive stars in STARFORGE simulations (with identical
physics and numerics) of “typical” Solar neighborhood GMCs, where Grudić et
al. (2022) showed the most massive stars in particular have the most extended
accretion histories often extending over an appreciable fraction of the GMC
lifetime, i.e. on the order of the large-scale dynamical time or multiple Myr
(also, notably longer than the protostellar lifetimes of such massive stars in
isolation). Physically, this relates directly to the dynamics of the stellar
orbits, and their accretion/mass ejection efficiencies after formation, which
we discuss below. It also means that the IMF in our QAD/CQM environment is
more directly “set at birth” by said QAD/CQM conditions, as compared to being
more strongly modified or self-regulated via subsequent accretion in more
typical GMC environments.

Figure 15: Distribution of velocity differences $\delta{\bf v}\equiv{\bf
v}_{\ast}-\langle v_{\rm gas}({\bf r}_{\ast})\rangle$ between individual stars
(${\bf v}_{\ast}$) and their surrounding gas ($\langle v_{\rm gas}({\bf
r}_{\ast})\rangle$, defined as the mass-weighted mean velocity of gas within a
resolution element $\Delta x$ surrounding each star). We plot both the signed
$\delta{\bf v}$ (left; divided into radial/azimuthal/vertical components) and
absolute value $|\delta{\bf v}|$, both for very young stars (top; time since
sink formation $<100\,$yr) and all stars (bottom). Even within $<100\,$yr
after sink formation, stars have acquired a typical velocity difference of
order $\sim 10\%$ of the circular velocity $V_{c}$ at their formation radius
from the gas around them. Per § 4.3.1 this must happen owing to the extremely
short dynamical times at these radii: stars and gas cannot be co-orbital owing
to pressure/magnetic/radiation forces on the gas, which generate a velocity
difference $|\delta{\bf v}|\sim(H/R)\,V_{c}$ on a timescale
$\lesssim\mathcal{O}(\Omega^{-1})$. Physically, this amounts to several
thousands (hundreds) of ${\rm km\,s^{-1}}$ in the QAD (CQM). This also
explains the very rapid suppression of accretion in Fig. 14. The older stars
with velocity offsets $\sim V_{c}$ are formed on highly eccentric orbits and
now “plunging” near pericenter through the inner QAD, whose precession (Fig.
13) produces velocity offsets $\mathcal{O}(V_{c})$.

####  4.3.1 The Formation of “Wandering Stars”

The suppression of stellar accretion just after cores form (§ 4.3) is closely
related to the phenomenon of “wandering stars” which we see are ubiquitous in
these simulations. Specifically, we see, following the positions of
progressively “older” stars (in e.g. Figs. 6, 7, 12, and more directly in Fig.
13), that very quickly after formation (on timescales as short as $\lesssim$
years in the QAD, or hundreds of years at $\sim 1\,$pc) stars “separate” from
the dense gas structures in which they form. On somewhat longer timescales
$\gtrsim\Omega^{-1}$, they end up in completely distinct orbits from the dense
gas. This is illustrated in Fig. 13 where we plot all the stars (including
more of the old stars), in the QAD (where the dynamical times are shorter so
this is even more evident), where we see that the ongoing precession of the
QAD has left behind a population of dynamically older stars whose orbits take
them further out of the QAD plane the older the star is (in units of
$\Omega^{-1}$), even though the youngest/newest sinks are clearly forming in
the midplane (Fig. 12). We stress that these are not stars scattered by
dynamical processes (e.g. interacting triples, sub-cluster merging) as often
occurs in LISM GMCs (Grudić et al., 2018b; Guszejnov et al., 2022a; Farias et
al., 2023): if it were so, we would see a much more isotropic and symmetric
stellar distribution, and it would occur over many dynamical times, not
$\mathcal{O}(1)$ dynamical time. Instead we see clear “one sided” examples
where old stars are all in the orbits that the QAD/CQM gas used to occupy but
the gas and stars have drifted systematically relative to one another.

We can also directly test this by measuring the relative motions between stars
and the gas surrounding them in Fig. 15, and validate it by checking in the
formal orbital parameters of the stars. There, we clearly see strong
systematic velocity offsets between stars and neighboring gas. For
dynamically-old stars, this is less surprising, but we see even for the
youngest stars (we obtain similar results in Fig. 15 whether we define this by
$t_{\ast}<\Omega^{-1}(r)$ or $t_{\ast}<100\,{\rm yr}$) a clear offset, in
their orbits. And we can further see this reflected in the “cometary trails”
made by stellar ejecta (jets and winds) as the stars traverse through the
gaseous medium as shown in Fig. 13 and discussed further below.

This is further enhanced by the gravitational structures giving rise to
overdense modes in the QAD. Broadly speaking, we see modes on two scales:
first, a very obvious global lopsided coherent eccentric ($m=1$) mode at all
radii (see details in Paper I). Such modes are generically expected around
quasi-Keplerian potentials like accreting SMBHs (Sellwood & Merritt, 1994;
Jacobs & Sellwood, 2001; Sambhus & Sridhar, 2002; Hopkins & Quataert, 2010a,
2011a, 2011b), and the same physics necessarily produces a phase offset
between the patterns for the collisional and collisionless components, which
contributes (coherently) to the velocity drift between stars and gas as soon
as stars form (see e.g. Noguchi, 1988; Wada, 1994; Barnes & Hernquist, 1996;
Hopkins & Quataert, 2011b; Berentzen et al., 2007). This global eccentricity
also naturally explains why the radial velocity offsets in Fig. 15 are
especially notable. On small scales, we see local collapsing dense gravito-
turbulent structures in the star-forming zones (§ 4.2). But these are at least
partially density waves and shocks (they are formally more related to “slow
modes” as in Tremaine 2001; Hopkins & Quataert 2010b; Hopkins 2010), so
material compressed in them “moves through them” at a large fraction of the
Keplerian velocity. As such, the timescale for separation is of order the size
of the structure $\sim\ell$, divided by some appreciably fraction of the
Keplerian speed $\sim\alpha\,V_{\rm c}$, or $\sim(\ell/\alpha\,r)\,t_{\rm
dyn}$. This is quite familiar from the well-studied problem (both
observationally and theoretically) of star formation in GMCs and especially
spiral arms in the galaxy (e.g. Tasker & Tan, 2009; Foyle et al., 2010; Leroy
et al., 2013; Colombo et al., 2014; Grasha et al., 2019; Meidt et al., 2022;
Lee et al., 2022). As has been studied for decades (Sellwood & Carlberg, 1984;
Young & Scoville, 1991), density modes lead to compression and shocks, which
produces rapid star formation, then stars drift through and separate from the
overdense gaseous arm region.

The physics causing this is actually quite straightforward. Recall, the gas
feels strong anisotropic stresses from pressure forces and, in particular,
magnetic fields on these scales (Paper II). As cores collapse into protostars,
they effectively de-couple from the background gas MHD forces, and become
“collisionless” (or at least “pressure-free”) in their gravitational dynamics.
This means that as soon as this force diminishes, they must begin to separate
from the gas positions, with typical relative velocities of order the ratio of
pressure forces to gravitational forces times the circular velocity – in
practice, up to tens of percent of the circular velocity at a given radius
(tens to hundreds of ${\rm km\,s^{-1}}$; see Fig. 5 and Paper II). This is
consistent with the behavior of the young stars in Fig. 15: we see a
characteristic velocity offset of order $\sim 0.1\,V_{c}$ imprinted at
formation. With time, a combination of eccentric orbits (e.g. the stars on
“plunging” orbits through the inner QAD), precession of the QAD, systemic
drift, and turbulence/new shocks creating new or shifting the locations of
filamentary structures in the outer CQM, lead to non-linear offsets in
velocities/orbits. This in turn means the local star-gas velocity difference
for stars more than a few dynamical times old is more like $\sim V_{c}$
(order-unity, in a relative sense).444Interestingly, unlike the simple case of
e.g. a planet in a nearly-circular PSD, in many cases the star-gas velocity
offset is such that the gas is actually locally moving faster in the mutual
direction of motion than the star, so the ejecta can be swept “forward.” This
can of course occur for many reasons, depending on e.g. where the stars are in
their non-linearly eccentric orbits and the precession of the QAD.

The key difference between the situation in the Solar neighborhood and the
simulations here is simply that the timescales for these offsets to occur, and
the magnitude of the gravitational velocities, are wildly different. Stars
drift out of e.g. spiral arms in the LISM on timescales of tens to hundreds of
Myr ($\sim 1/\Omega$ at Solar circle Galactocentric radii) – very long
compared to the protostellar collapse/accretion timescales of massive stars.
But here, separation on a fraction of $\sim 1/\Omega$ means timescales of
years, much shorter than typical protostellar evolution timescales. And in the
QAD, star-gas velocity separation at even $\sim 10\%$ of $V_{c}$ does not mean
a relative velocity of $\sim 10\,{\rm km\,s^{-1}}$ as in the LISM, but more
like several hundred ${\rm km\,s^{-1}}$.

Figure 16: Extreme examples of “cometary trails” of (proto)stellar
outflows/jets. We show images of the total gas density in a couple of the
filamentary regions from Fig. 6, and highlight gas cells (cyan, as Fig. 13)
whose mass is dominated by the direct ejecta from a couple individual
protostars in each sub-region. Owing to the rapid motion of the stars relative
to the gas even for quite young stars (Fig. 15), the ejected jet material
immediately encounters an extremely high ram and magnetic pressure
($\sim\rho\,V_{c}^{2}$) which produces a bow shock/stream/tail morphology.

####  4.3.2 Consequences for Stellar Accretion and Feedback?

The separation between stellar positions and their formation sites in gas (§
4.3.1) has immediate important consequences for their ability to continue
accreting gas, which is reflected in the time histories in § 4.3. Consider the
Bondi-Hoyle-Lyttleton capture rate onto a core (moving through a statistically
homogeneous, isotropic ambient medium with no self-gravity and no other
forces), which scales roughly as $\dot{M}\sim 4\pi\,G^{2}\,M_{\rm
core}^{2}\,\rho_{\rm gas}/v_{\rm eff}^{3}$ where $\rho_{\rm gas}$ is the
ambient gas density and $v_{\rm eff}^{2}\sim
c_{s}^{2}+\langle(1-\cos^{2}{\theta})\,\rangle\,v_{A}^{2}+|{\bf v}_{\rm
core}-\langle{\bf v}_{\rm gas}\rangle|^{2}+|\delta{\bf v}_{\rm turb}|^{2}$
collects contributions from the thermal sound speed ($c_{s}$), magnetic
support (the Alfvén speed $v_{A}$, where
$\langle(1-\cos^{2}{\theta})\,\rangle$ depends on the geometry and is $\sim
2/3$ for tangled fields, see Lee et al. 2014), the bulk relative velocity
between core and gas $|{\bf v}_{\rm core}-\langle{\bf v}_{\rm gas}\rangle|$,
and small-scale turbulence $|\delta{\bf v}_{\rm turb}|$.

If we assumed a strictly laminar, non-magnetized, non-turbulent, thermal-
pressure supported QAD with zero global/eccentric structure and zero drift
velocity or asymmetric drift/pressure support of the gas and the stars
perfectly co-orbital with the gas on aligned co-rotating circular orbits, then
$v_{\rm eff}\sim c_{s}$. But in our simulations, the thermal sound speed is
small compared to $v_{A}$ by a typical factor of $\sim 100$ (i.e. $\beta\sim
10^{-6}-10^{-2}$ is very small). The turbulence is trans-Alfvénic, so
$|\delta{\bf v}_{\rm turb}|\sim v_{A}$. And the effects described in § 4.3.1
each naturally lead to systematic/drift velocity offsets of $|{\bf v}_{\rm
core}-\langle{\bf v}_{\rm gas}\rangle|\sim(H/R)\,V_{c}\sim
v_{A}\sim|\delta{\bf v}_{\rm turb}|$ (see Fig. 15) and add coherently. So we
expect $v_{\rm eff}\sim{\rm a\ few}\times v_{A}$. And since the accretion rate
is suppressed by $v_{\rm eff}^{-3}$, this means that the accretion rates will
typically be suppressed relative to the laminar, co-moving expectation by a
factor of $\sim(0.01-0.1)\times\beta^{3/2}\sim 10^{-8}-10^{-5}$. Thus almost
as soon as the cores begin to “de-couple” after their density becomes
sufficiently high as described in § 4.3.1, the accretion rates will drop by
orders of magnitude (and we have even neglected the additional linear
suppression factor that will arise from $\rho_{\rm gas}$ decreasing as the
stars move away from the local density maxima where they initially formed).

Worse yet, stars from outer radii on eccentric orbits intersect the “inner”
QAD (at their pericentric passage) with $\mathcal{O}(1)$ angle between their
eccentric velocity vector and the local mean gas orbital velocity vector,
owing to the fact that the QAD+CQM structure, eccentricity, and orbital plane
are not perfectly independent of radius $r$. So in these intersection cases
although the density rises, the velocity also rises to
$\sim\mathcal{O}(1)\times V_{c}(r)$, which net suppresses the accretion rate
by an additional factor of $\sim 10^{-4}-10^{-3}$ or so.

The same relative velocity argument also immediately explains why we see the
bow shocks and almost-cometary “trails” left behind by ejecta from
(proto)stars, which we show some typical examples of in Fig. 16. As soon as
material is ejected from a star, it re-couples to these ambient RMHD forces,
introducing a drift velocity relative to the star which ranges from $\sim
10\%$ to a couple times the circular velocity. Since the star always feels
some coherent gas velocity from its point of view (analogous to e.g. a
planetesimal or rocky body in a circum-stellar PSD), the ejecta are “swept up”
in the stellar frame. Of course some examples of this (albeit much less common
and with much lower absolute velocity scales) are known in protostars and YSOs
(as well as AGB stars), such as L1551 IRS 5 or PV Cep (Goodman & Arce, 2004),
but the phenomenon here is more pronounced and ubiquitous owing to the extreme
dynamical conditions.

Figure 17: Distribution of stellar ejecta (jets/outflows/winds) in the broader
QAD/CQM, as Fig. 16 but focusing on larger scales and showing more of the
ejecta mass. Left: Gas density map as Fig. 6, but with blue colorscale overlay
highlighting gas weighted by the mass fraction of ejecta cells (direct ejecta
cells “spawned” by protostellar particles to represent winds and jets)
recently merged into the “ambient” (standard) gas cells. These correspond to
regions of the QAD/CQM with particularly recent very massive star formation.
Right: Zooming in on the QAD from left, and overplotting all un-merged ejecta
cells as Fig. 16. In both cases the ambient gas interacting with or recently
affected by ejecta show no deviation in their kinematics or density structure
from the larger-scale “background” flow, clearly indicating that stellar
feedback has very little dynamical effect on the CQM/QAD. This is expected (§
4.4), owing to the extremely deep potential combined with suppressed/very low
overall SFR.

###  4.4 The (Weak) Effects of Stellar Feedback On the Environment

The simulations exhibit very weak global effects of stellar feedback on
environment. In Fig. 17, we see that stellar ejecta in both the CQM
filamentary zone and disky QAD zone are simply mixed/incorporated into the
surrounding gas with orbits typical of the ambient gas that has not been
directly inflenced by feedback. We have evolved the system for many local
dynamical times (up to $\gtrsim 10^{5}\,\Omega^{-1}$ at the innermost radii
simulated) so this remains true in a “local steady state” sense. And we show
in Paper I and Paper II that the contribution of feedback on sub-pc scales to
the turbulent or dynamical or binding energy budget of the QAD itself and/or
infalling gas into the QAD is negligible.

This is anticipated by many previous studies of stellar feedback processes in
circum-quasar environments (see Thompson et al. 2005; Wada et al. 2009;
Hopkins et al. 2016; Grudić et al. 2019; Kawakatu et al. 2020; Anglés-Alcázar
et al. 2021, or for a review, Hopkins et al. 2022b). So we only briefly review
the reasons here. To begin, note that inside of $r\sim(0.01,\,0.1,\,1)\,{\rm
pc}$, the star formation+stellar accretion rate is
$\dot{M}_{\ast}(<r)\sim(10^{-4},\,0.1,\,10)\,{\rm M_{\odot}\,yr^{-1}}$, while
the gas inflow rate is $\dot{M}_{\rm in}\sim(40,\,70,\,100)\,{\rm
M_{\odot}\,yr^{-1}}$ (see § 3 or Paper I), so clearly the SFR and stellar
mass-loss/jets (even if mass-loss returned the entire accreted stellar mass
budget) are small perturbations to the inflow mass budget.

In terms of momentum/force balance, recall that for our “slow” PSD accretion
models, accretion from the PSD to the star (which is what actually determines
the feedback rates) is very slow compared to the local dynamical time, so
feedback is strongly suppressed. We therefore instead consider the opposite
“fast” PSD accretion regime, where one can approximate accretion from PSD to
star as instantaneous so that for calculating IMF-integrated fluxes from e.g.
jets, we simply have $\dot{M}_{\rm jets}\sim f_{\rm jet}\,\dot{M}_{\ast}$.
With that in mind, we can directly compute the momentum flux in jets
$\dot{P}_{\rm jets}\sim\langle f_{\rm jet}\,\dot{M}_{\ast}\,v_{\rm
jet}\rangle$ (with $f_{\rm jet}\sim 0.3$ the mass fraction going into jets and
$\langle v_{\rm jet}\rangle\sim 0.3\,\sqrt{G\,M_{\odot}/R_{\odot}}$, roughly),
or from main sequence winds plus radiation pressure $\dot{P}_{\ast}\sim
10^{3}\,L_{\odot}\,M_{\ast}/(M_{\odot}\,c)$, and compare to the gravitational
force on the gas $\dot{P}_{\rm grav}\sim G\,M_{\rm BH}\,M_{\rm
gas}(<r)/r^{2}$, and obtain $\dot{P}_{\rm jets}/\dot{P}_{\rm
grav}\sim(10^{-7},\,10^{-4},\,0.01)$ and $\dot{P}_{\ast}/\dot{P}_{\rm
grav}\sim(4\times 10^{-8},\,10^{-5},\,5\times 10^{-4})$, at
$r\sim(0.01,\,0.1,\,1)\,{\rm pc}$. This is consistent with numerous studies
(see references above and Fall et al. 2010; Colín et al. 2013; Gavagnin et al.
2017; Grudić et al. 2018a; Kim et al. 2018) which have shown that simply
balancing the momentum injection rate for young stellar populations against
the gravitational force per unit area (here $\sim|{\bf a}_{\rm
grav}|\,\Sigma_{\rm gas}\sim G\,M_{\rm bh}\,\Sigma_{\rm gas}/r^{2}$) means
that feedback becomes sub-dominant when the acceleration scale exceeds $|{\bf
a}_{\rm grav}|\gg 10^{-8}\,{\rm cm\,s^{-2}}$ (equivalent to an “effective
surface density” $(M_{\rm bh}+M_{\rm gas}+M_{\rm other})/r^{2}\gg 1000\,{\rm
M_{\odot}\,pc^{-2}}$). But here $|{\bf a}_{\rm grav}|\sim G\,M_{\rm
bh}/r^{2}\sim 2\times 10^{-4}\,({\rm pc}/R)^{2}$ exceeds this critical
threshold by enormous ratios of $\gtrsim 10^{4}$. In other words, in the
vicinity of the SMBH, the gravitational forces competing against feedback are
equivalent to those of a cloud with effective surface density $\sim
10^{7}\,{\rm M_{\odot}\,pc^{-2}}\times({\rm pc}/R)^{2}$.

Figure 18: Stellar IMFs (§ 4.5) in the simulation CQM+QAD. Here we plot the
(normalized) mass function $N_{s}^{-1}dN_{s}/d\ln{M_{s}}$ (including the total
mass of the sink, both that assigned to “star” and accreted “reservoir”),
including all stars within the CQM and QAD (radii within the BHROI, $r\lesssim
5\,$pc), at the final time to which the simulation is run ($\sim 10^{5}$
dynamical times at the innermost radii). The normalization of the IMFs is
arbitrary. We include all stars formed during the simulation after it reaches
the highest refinement level (so we exclude contamination by low-resolution
from earlier pre-refinement stages, though this has only a small effect here).
We compare the IMF of just the dynamically-young stars with
$t_{\ast}<\Omega^{-1}$. We compare simulations which adopt the “fast” (left)
or “slow” (right) PSD-to-protostar accretion model (§ 2.2): as discussed in §
4.5.2, this has an obvious effect on how “top-heavy” the IMF is, where the
“slow” model, by suppressing accretion from PSD to star also suppresses
(proto)stellar feedback and leads to accretion of more initial bound clump
mass (hence a more top-heavy IMF). We compare to the canonical observed
“universal” single-star IMFs from Kroupa (2001) and Chabrier (2005), in their
standard form (with a Salpeter 1955-like high-mass slope $\alpha=2.3$) or with
a modified top-heavy high-mass slope $\alpha=1.8$ (§ 4.5.1). The lack of the
“bump” at $\sim 10-20\,M_{\odot}$ for young stars in the “slow” accretion run
owes to the lack of recent star formation in the outer QAD zone (where the IMF
is most top-heavy) at this particular moment.  Figure 19: IMFs as Fig. 18, but
comparing with the Kroupa (2001) and Chabrier (2005) LISM system IMFs as well
as single-star IMFs. As we show below, fragmentation and multiplicity are
strongly suppressed by the extreme tidal fields in the simulations. This may
explain the modest shift in the peak of the IMF for the “fast” accretion
model, but both the high-mass slope and IMF peak/turnover mass are shifted
higher than even the LISM system IMF.

In terms of energetics and heating, noting that the circular velocities at
these radii are $V_{c}\sim(2500,\,750,\,250)\,{\rm km\,s^{-1}}$, and that the
gravitational luminosity (energy released by accretion, $\dot{E}_{\rm
grav}\sim\dot{M}_{\rm in}\,V_{c}^{2}$), turbulent dissipation rate
($\dot{E}_{\rm cool}\sim M_{\rm gas}\,v_{\rm turb}^{2}\,\Omega$), and
radiative cooling losses ($\dot{E}_{\rm cool}$) are all similar (as expected
for an accretion-regulated system), we have $\dot{E}_{\rm
grav}\sim(10^{44},\,3\times 10^{43},\,5\times 10^{42})\,{\rm erg\,s^{-1}}$ at
$r\sim(0.01,\,0.1,\,1)\,{\rm pc}$. We can directly compare this to the kinetic
luminosity of jets (or O/B winds, which are comparable) $\dot{E}_{\rm
jets}\sim\langle f_{\rm get}\,\dot{M}_{\ast}v_{\rm
jet}^{2}\rangle\sim\sim(5\times 10^{35},\,5\times 10^{38},\,5\times
10^{40})\,{\rm erg\,s^{-1}}$, and see that the latter is negligible in the
global energy budget. Essentially, the escape velocities from the SMBH
potential are larger than typical jet velocities, so the SFR would have to be
much larger than QAD inflow rates (the opposite of what we see) for the jets
to impact the kinematics globally. For comparison, a typical GMC would have a
gravitational luminosity of $\sim 10^{37}\,{\rm erg\,s^{-1}}\,(M_{\rm
cloud}/10^{6}\,M_{\odot})^{5/4}\,(\Sigma_{\rm cloud}/100\,{\rm
M_{\odot}\,pc^{-2}})^{5/4}$, so again the extreme potential of the SMBH makes
this more akin to the gravitational dynamics of a “cloud” with effective
surface densities $\gtrsim 10^{7}\,{\rm M_{\odot}\,pc^{-2}}$.

But what about radiative heating? Summing the accretion/contraction/main
sequence luminosities of the stars, we obtain a total radiative luminosity of
stars $L_{\ast}\sim(10^{38},\,10^{41},\,10^{43})\,{\rm erg\,s^{-1}}$ at these
radii. Recall that at $\sim 1$ pc the QAD/CQM has an IR optical depth $\gg 1$,
so this will mostly be re-processed locally (hence considering the local
$L_{\ast}$ end $\dot{E}_{\rm grav}$ at each $R$, instead of external
illumination by either the $\sim$ kpc-scale galactic starburst or inner QAD,
whose luminosity here can approach $\sim 10^{47}\,{\rm erg\,s^{-1}}$). Thus at
$R\ll$ pc, the heating by stars is negligible in the total QAD/CQM thermal
budget, but at $R\gtrsim\,$pc, stellar heating becomes comparable to or larger
than accretion luminosities, contributing importantly to the thermal balance
of the medium. Indeed, it was noted in Paper I & Paper II that the typical CQM
temperatures at these larger radii (approaching the BHROI) were clearly
elevated compared to those predicted from models with only accretion
luminosity (e.g. traditional QAD models, extrapolated to $\gtrsim$ pc scales,
using the surface densities/optical depths actually appearing in the
simulations555In Fig. 5, this is evident in the somewhat hotter mass-weighted
(closer to midplane) temperatures compared to the SS73 model at $\gtrsim$ pc
scales, but it is worth noting that SS73 predicts a higher surface density and
optical depth as noted therein, so the extra heating from e.g. ambient stars
on gas is even more significant than this might naively imply at $r\gg
1\,$pc.).

Finally, outside of a few pc (the BHROI), stars begin to dominate the
potential (Paper I) and the SFRs continue to grow, with the medium resembling
a typical starburst environment, rather than QAD/CQM. There, stellar feedback
may play a much more significant role.

###  4.5 The IMF

####  4.5.1 Overview & Shape

Fig. 18 shows the resulting initial mass function of stars in the simulations,
at the latest times we evolve. Below, we break down the IMF by region, time,
and physics model in more detail, but to begin, we simply show it at the
latest possible time. Given the curves of growth in Fig. 14 above, it makes
very little difference if we plot the “present-day mass function” (PDMF) at
this time, or the “peak” IMF where we define the “mass” of each sink by the
maximum mass it achieves. This is because, as discussed above, the dynamical
timescales of interest and duration of the simulations are quite short
compared to stellar evolution timescales. Here we plot the “total sink mass”
defined as the total bound mass associated with the sink particle – both what
the code would actually call the “(proto)star” and the “circum-stellar PSD
mass” or “accretion reservoir.” As we discuss below, this is more robust,
given our uncertainty over the rate at which material accretes from reservoir
to “star.” Moreover, if there is unresolved fragmentation into e.g. close
binaries which cannot be captured explicitly in our simulations, this mass
will be robust to that effect.

Figure 20: IMF for the “fast” (left) and “slow” (right) PSD accretion models
as Fig. 18, comparing the “total sink mass” $M_{s}$ (explicitly-evolved
“stellar” mass $M_{\ast}$ plus the “reservoir” accreted by the sign but not
yet assigned to the star $M_{\rm PSD}$) to the mass assigned just to the
“star” $M_{\ast}$ (excluding the reservoir). We show this for all stars. As
expected, with the “fast” PSD accretion model (§ 4.5.2), the PSD mass rapidly
accretes on to the “star” so $M_{s}\approx M_{\ast}$, though at the smallest
radii even the very fast assumed accretion can lead to some differences. In
contrast, with the “slow” PSD accretion model nearly all ($>90\%$) of the
accreted sink mass stays in the PSD “reservoir” and does not reach the “star”
over the entire duration of the simulation ($\sim 10^{4}\,{\rm yr}$, or $\sim
10^{5}\,\Omega^{-1}$ at the innermost resolved QAD radii) meaning that there
is much weaker local stellar feedback/ejection/suppression of accretion of the
original core (§ 4.5.2), explaining the more top-heavy IMF.

Integrated over the nuclear region, we see star formation with (proto)stellar
system masses ranging from $<0.01\,{\rm M_{\odot}}$ to $\gtrsim 200\,{\rm
M_{\odot}}$ is possible. It is also a robust result that the IMF is “top-
heavy” (in the sense of having a shallower high-mass slope) when compared to
either the “mean” IMF or observed variation in IMFs in nearby clouds in
various compilations (Kroupa, 2001; Chabrier, 2003, 2005; Bastian et al.,
2010; Offner et al., 2014; Hopkins, 2018). It is less robust whether the
lower-mass end of the IMF is also “top heavy” (i.e. shifted towards higher
masses) – this appears to be more sensitive to our assumptions and where/when
we measure it (as well as whether we compare to the individual-star or system
IMFs in e.g. Chabrier 2003, 2005), as we discuss in more detail below.

It is interesting to note that this is at least qualitatively similar to the
inferred IMF in the sub-pc region around Sgr A∗ in the Galactic center
(Paumard et al., 2006; Bartko et al., 2010; Lu et al., 2013; Hosek et al.,
2019). In fact, fitting a slope to the “fast sub-grid accretion” simulation
gives a high-mass slope $dN/dM\propto M^{-\alpha}$ with $\alpha\sim 1.6$
(where Salpeter is $\sim 2.3$), quite similar to that inferred for the central
$<0.5\,$pc around Sgr A∗ in Lu et al. (2013) (and also broadly similar to the
slightly-less-top-heavy claim for Arches in Hosek et al. 2019). The “slow sub-
grid accretion” model predicts an even more top-heavy IMF, though still not
quite as top-heavy as the most top-heavy estimates from e.g. Bartko et al.
(2010) ($\alpha\sim 0.5$) for the Galactic center (though see Löckmann et al.,
2010). This qualitative agreement – that both are top-heavy – is encouraging,
though we caution against taking any quantitative comparison too literally.
Not only is the detailed IMF shape in the Galactic center still controversial,
and clearly physics-dependent here, but we very much are not attempting to
model anything like the conditions in the Milky Way center. While the SMBH
mass here is only a factor of $\sim 3$ larger than Sgr A∗, our default
simulation represents an extremely bright quasar in a massive starburst galaxy
at redshift $z\sim 4.4$, quite distinct from the event which formed the
observed stars around Sgr A∗. And even our default LISM-cloud STARFORGE
simulations often produce a slightly flat ($\alpha\sim 2$) IMF (Guszejnov et
al., 2022b). Still, at least some of the physics arguments below should apply
to both regimes.

####  4.5.2 System versus Single-Star IMFs and Dependence on the Sub-Grid
Model for Accretion from Reservoir to (Proto)Star

In Fig. 18, we clearly see systematic differences between the simulations
where we adopt a “fast” versus “slow” sub-grid model for accretion from the
captured PSD “reservoir” onto the (proto)star within the sink particles. Fig.
19 compares the same IMFs from the simulations to both the individual-star and
system IMFs from the Milky Way/LISM in Kroupa (2001) and Chabrier (2005). For
the “slow” accretion case, the IMF is top-heavy both in terms of the high-mass
slope and in terms of the location of the peak/turnover mass (here at $\sim
0.25-1\,{\rm M_{\odot}}$), even compared to the system IMF. For the “fast”
accretion case the results are a bit more ambiguous. The high-mass slope is
shallower than Salpeter, but as noted above not nearly as shallow as in the
“slow” accretion case. The IMF peak in the “fast” accretion case at $\sim
0.2\,{\rm M_{\odot}}$ is clearly somewhat higher than the LISM individual-star
IMF, but notably similar to the system IMF peak. This may suggest that the
major difference here is the suppression of fragmentation/binarity (an effect
we examine explicitly below). However, resolution may also be playing a role
for these more modest offsets at sink masses which are only $\sim 40-60$ times
our median mass resolution in the QAD (see § 5).

Turning to why this may occur, in Fig. 20 we can see that this choice of
accretion model has a huge effect on the in-simulation sub-grid model
breakdown between “star” versus total sink mass (including the “reservoir”),
as speculated in § 2. For the “fast” accretion model, where the timescale for
accretion from the reservoir to star is, by construction, always faster than
the local dynamical/tidal timescale, the “star” masses closely track the
“system” or “total” masses. For the “slow” accretion model, which uses a
constant accretion time $\sim 2000-5000\,$yr from PSD onto the star, often
longer than the local dynamical time in the CQM, we see a much larger offset,
where the “star” mass often is $\sim 30$ or so times smaller than its
“reservoir” mass. This trend is expected: the key difference is when $t_{\rm
acc}\ll t_{\rm dyn}\sim 1/\Omega$, the mass is rapidly transferred to the sink
(so sink and total masses are similar) while when $t_{\rm acc}\gg t_{\rm dyn}$
the mass simply “piles up” in the reservoir/PSD. This in itself is neither
surprising nor predictive: we are simply getting out what we put in, and any
value of $t_{\rm acc}$ is explicitly sub-grid and so the ratio of “sink” to
“reservoir” mass (or other potential sub-grid choices like whether or not we
assume that, below the resolution, a sink should or could fragment into close
binaries) would simply reflect the analytic prescription chosen rather than
any self-consistent prediction of the simulation.

But what is at least somewhat interesting and predictive is that, even if we
consider just the total sink masses, the “slow” accretion model is
systematically much more top-heavy in Figs. 18-20, whether we consider just
recently-formed sinks or all sinks formed since the highest resolution was
reached. If everything were perfectly identical in these simulations on all
scales, these total masses would be the same (the only difference would be in
the sub-grid mass breakdown). However we stress that in every measureable
global QAD/CQM property we considered in Paper I (on the global
thermochemistry and dynamics and star formation rate) or Paper II (on the
structure of the magnetized accretion flow to the SMBH), this choice of models
make no difference: the primary difference appears to be restricted to the
stars themselves, and must be manifesting on scales of order the individual
collapsing cores.

Figure 21: IMF for the “fast” (left) and “slow” (right) PSD accretion models
as Fig. 18, but dividing stars into different distances $r$ (in pc) from the
SMBH (§ 4.5.3). We compare all stars at the end of the simulation (top), and
young stars with time since formation $t_{\ast}<t_{\rm dyn}=\Omega^{-1}(r)$
(middle) or $t_{\ast}<100\,{\rm yr}$ (bottom). the relative normalization of
different curves reflects the different numbers of stars (we do not normalize
by $N_{s}^{-1}$), so we can compare how many stars form at different times and
places. At radii approaching the outer CQM/BHROI, $r\gtrsim$ pc, where tidal
forces from the SMBH become small compared to self-gravity, the IMF approaches
the “universal” IMF form shown, especially in the Salpeter slope (even though
ISM conditions still differ; § 3). At intermediate radii, particularly the
outer QAD around $\sim 0.1-0.5\,$pc, the IMF of young stars becomes notably
more top-heavy. At the smallest radii ($r\ll 0.1$ pc) whether it becomes more
top-heavy or bottom-heavy depends on the PSD accretion model assumptions. At
radii $r\lesssim 0.01\,$pc there is essentially no measurable star formation
in the $\sim 10^{5}$ dynamical times simulated here. The IMF of stars at a
given present-day position, including all stars, can be influenced by stars
not formed at those radii – e.g. the massive stars at small radii in the
“slow” run primarily formed in the top-heavy outer-QAD ($\sim 0.1-0.5\,$pc)
region, but on highly eccentric orbits and are now plunging through the QAD.

This suggests that, while it plays a relatively weak role on large scales,
feedback from (proto)stars still plays an important role in regulating the
total sink masses. In the “slow” accretion limit, $t_{\rm acc}\gg t_{\rm
dyn}$, when mass is accreted into the sink it simply sits in the “reservoir,”
doing nothing in our simulations, as all of the feedback processes we model
(protostellar jets, main-sequence mass-loss, radiation, supernovae) depend on
some combination of the mass of and mass accretion rate onto the actual sub-
grid star, not the reservoir mass. In the “fast” accretion limit, the usual
assumptions of efficient feedback (e.g. that an order-unity fraction of the
swallowed mass is almost-instantaneously returned to the ISM in the form of
energetic protostellar jets) hold.

So given the strong effects of stellar feedback well-studied under more
“typical” conditions, especially their known ability to regulate the massive
end of the IMF (Guszejnov et al., 2021, 2022b; Grudić et al., 2022) it is not
surprising that we see a less top-heavy IMF (suppression of the most massive
system masses) in the “fast” accretion model. The key is that since the
primary formation/accretion period for the stars is extremely short-lived or
rapid in terms of physical years (§ 4.3.2), the accretion onto the (proto)star
must be correspondingly fast to enable any effect of (proto)stellar feedback.

####  4.5.3 IMF versus Distance from the SMBH

Fig. 21 plots the IMF versus “birth” location (distance from the SMBH, i.e.
initial $R_{\rm initial}$ in the QAD/CQM). Note here we have to be careful to
separate the birth location from the instantaneous location of the star, for
the reasons in § 4.3.1 (stars move significantly in their orbits over the
simulation duration on these small scales). So for both the “fast” and “slow”
accretion models, we compare the IMFs at different radii for: (1) all sinks
formed since the simulations reach their maximum resolution, (2) “dynamically
young” sinks which formed less than one dynamical time ago
($t_{\ast}<\Omega^{-1}$, with $\Omega$ measured at the present radius $r$ of
the sink), and (3) young sinks with absolute formation age $t_{\ast}<100\,$yr.

At $R_{\rm initial}\gtrsim 0.5-1\,$pc, the IMF increasingly becomes “normal”
or only slightly top-heavy. In particular, for both the “fast” and “slow”
accretion cases, the high-mass slope appears to be rapidly converging to
Salpeter-like with increasing distance from the SMBH. Interestingly, for both,
at $0.5<r<1$ pc, the IMF peak mass agrees very well with the observed
individual-star IMF, while at somewhat larger radii ($\sim 1-5\,$pc) it
becomes mildly top-heavy (shifting by a factor of $\sim 2-3$ to higher
masses). It is possible this reflects real trends with e.g. the opacity limit
or the amount of self-shielding to the intense galactic radiation field, but
as discussed in § 5, this is precisely the region of the disk where resolution
concerns are maximized, and (owing to computational expense), we cannot run at
our highest resolution for very long (the dynamical time at $5\,$pc reaches
$\Omega^{-1}\sim 5\times 10^{4}\,$yr), so it is also not clear if the IMF has
actually reached steady-state at these very largest radii. As a result, we
urge against over-interpreting the any apparent trend going from $\sim
0.5-1\,$pc to $\sim 1-5\,$pc.

In any case, with these caveats in mind, it is clear that on scales $\gtrsim
0.5\,$pc well outside the QAD (but still in the CQM), the IMFs approach
something more like the “universal” IMF. This is noteworthy because the
conditions at these radii, as reviewed in § 3, are still radically different
from the Solar neighborhood ISM, with much higher gas densities and column
densities and optical depths, warmer molecular/dust/CMB temperatures, stronger
turbulence and magnetic fields, and still significant tidal effects from being
within the BHROI. Yet it appears that there are not at least radical
deviations in the IMF predicted. We discuss the implications for IMF physics
in § 4.6 below.

From $R_{\rm initial}\sim 0.05-0.5\,$pc (the location of the outer QAD) on the
other hand, the “top-heavy” behavior of the IMFs is quite dramatic. Even in
the “fast” accretion case, we see that (1) the high-mass slope is quite
shallow, between $\sim 1.6$ for the dynamically-young stars and $\sim 0$ for
the very youngest ($<100\,$yr-old) stars (though we caution the statistics are
poor for that latter category so this is driven by just a few sinks). And the
turnover mass is also shifted by a factor of a few to higher masses, despite
the mass resolution being finer and duration of the simulation much longer. It
is clear that especially for the “fast” accretion model, the “overall” top-
heavy behavior of the total IMFs seen in Figs. 18-20 is driven in large part
by the stars forming in the outer QAD (mostly around $\sim 0.1-0.3\,$pc).

Interestingly, if we restrict to young stars forming at the radii discussed so
far ($\gtrsim 0.05\,$pc), the differences between the “slow” and “fast”
accretion models are not nearly as dramatic as they appeared in Fig. 18. It is
still clearly the case that the “slow” accretion runs produce a more-top-heavy
IMF at each radius, but it suggests that at least some of the most extreme
time-and-location-integrated difference seen in Fig. 18 owes not to the local
physics of the IMF at a given point, but to the question of how many stars
(relatively speaking) formed at a given location (since the IMF is quite
sensitive to location).

At slightly smaller radii, $0.02\lesssim R_{\rm initial}\lesssim 0.05$ pc, we
see a more qualitative divergence in the IMFs. In the “fast” model, the IMF
becomes less top-heavy again, with very few massive stars. But in the “slow”
model, it remains comparably top-heavy to the IMF at slightly larger radii
$0.05-0.5$ pc. While the “more top heavy” qualitative character of the “slow”
model is expected, the more dramatic apparent difference here may have to do
as much with accidents of exactly where flips in the polarity of ${\bf B}$
occur (see § 4.2). It is worth noting that almost all the sinks at this radius
form towards the very outer edge of the “bin” (at more like $\sim
0.04-0.05\,$pc). But it may have interesting implications for the physics
driving the IMF in each case.

At $R_{\rm initial}\lesssim 0.03$ pc, there is essentially no sink formation.
Sinks can exist at these radii, if we consider “all” stars, but we see zero
very young sinks here, because sink formation is strongly suppressed by the
strong toroidal magnetic fields and rapidly-rising thermal Toomre $Q$ at radii
interior to this (see Fig. 5 & § 4.1.2-4.2, as well as Paper I for more
details). We therefore cannot define a true “initial” mass function here.
Rather, stars migrate into these radii from larger formation radii. In the
“slow” run, this leads to the “apparent” IMF (simply taking all stars,
remembering that they all have ages $<10^{4}\,$yr, which can be extremely
large compared to the dynamical time at $\lesssim 0.03\,$pc) being extremely
top-heavy. But as we discuss below, this must reflect more which stars make it
into these radii (and we similarly see the all-stars mass function at
$0.02<r<0.05\,$pc being more top-heavy than the young-star IMF at the same
radii when stars actually form).

It is important to note, in thinking about Fig. 21 in particular, that the
IMFs here are normalized by number fraction, so they can all be compared on
the same vertical axis. Of course, as is evident in Figs. 4-5 or Figs. 6-7
that vastly different numbers of stars are appear in each bin. Even in the
“all stars” MFs, interior to $\lesssim 0.05$ pc there are only a few stars, so
just a couple massive stars making their way to these smaller radii has a huge
effect on the IMF. Most of the stars in the QAD are clearly forming in the
$0.05<r<0.5$ radial range (see Fig. 7), while overall the total number of
stars/sinks in the high-resolution simulation region is dominated by stars in
the CQM at $\gtrsim 1\,$pc (Fig. 6).

Figure 22: IMF for the “slow” PSD accretion model versus distance from the
SMBH (as Fig. 21), comparing all stars formed after reaching the highest
refinement level (left) and dynamically young stars with $t_{\ast}<t_{\rm
dyn}(r)=\Omega^{-1}(r)$ (right), at four different times (early-to-late from
top to bottom). If we only consider dynamically-young stars – i.e. stars close
to their “birth sites” – the IMF at each $r$ is quite stable in time. Per §
4.5.4, the global evolution of the “all stars” IMF owes to two effects: (1)
stars moving from their birth radii (the dominant effect at $r\lesssim
0.05\,{\rm pc}$ and reason for the existence of any stars at $r\lesssim
0.02\,{\rm pc}$); and (2) change in the relative number of stars forming at
each radius with times (owing to e.g. the global mass flow, advection of the
sign flips in $B_{\phi}$, etc.). Results are similar for the “fast” case.

####  4.5.4 Evolution over Time

Fig. 22 plots the IMFs versus time in the simulations. We focus on the “slow”
model, since in Fig. 18 this appears to show a notable difference between the
IMF integrated over all time since the simulation began versus that of just
young sinks (though of course there has been no time to have objects migrate
off the main sequence). All our conclusions qualitatively cross-apply to the
“fast” model case, but the differences there are much smaller.

Specifically, in Fig. 22 we plot both the MF (versus radius of formation as
Fig. 21) of all sinks formed since the simulation reaches maximum refinement,
and just dynamically young sinks with $t_{\ast}<\Omega^{-1}$ at their radii,
at four equally-spaced times leading to the final time (what we plotted in
Fig. 21). Note that at much earlier times than the earliest already plotted,
there are simply not enough sinks to meaningfully sample the MFs.

If we first just focus on the dynamically-young stars, the real “IMF” here,
then we see that there is not actually much significant evolution in the MF of
sinks just formed at a given radius. There are some subtle effects, but most
of these are within the sampling/Poisson errors (binning by both time and
radius and age leads to only a few sinks in several of the bins). For example,
the IMF at $0.5<r<1$ pc may become slightly less top-heavy after the earlier
times, as the CQM on these radii relaxes somewhat and comes more into steady-
state after refinement is completed, but the effect is modest. Overall, this
is reasonable: as shown in Paper I and Paper II, the global QAD and CQM
properties are in reasonable steady-state for the duration of our simulations,
since these are short compared to the global dynamical time of the gas flows
from $\gtrsim 10\,$pc fueling the inflows within the BHROI.

Interestingly, however, if we examine the MF of all sinks versus time, there
appears to be more dramatic evolution of the MFs at $r<0.05\,$pc, in
particular. These become much more top-heavy as time goes by. Since especially
at $r<0.02\,$pc this is just a few sinks, we have checked that this evolution
is not driven by some sort of late accretion episodes onto these stars –
indeed, consistent with what we showed in § 4.3 and Fig. 14, these stars
barely accrete after their dynamically-young period. Rather, this trend comes
entirely from “new” stars (stars which had not previously been at
$r<0.02\,$pc) appearing at these radii despite forming at larger radii. Note
that at just larger radii, $0.05<r<0.5\,$pc, the IMF is quite top-heavy, so
almost any migration of stars from that radius will make the MF at smaller
radii more top-heavy. It is possible that the most-massive stars migrate
especially quickly owing to effects like dynamical friction, but as far as we
can tell, these stars are not “sinking” via something like dynamical friction
(in the sense of systematically migrating in on slowly-decaying nearly-
circular midplane orbits) but rather “passing through” the disk on highly-
eccentric, often inclined orbits (see Figs. 13 & 15).

####  4.5.5 What Properties do the “Reservoirs” Have, and Should They Accrete
Rapidly, or Fragment?

Given the difference in § 4.5.2 between the “fast” vs. “slow” sub-grid
accretion models for our sink particles, it is natural to ask what “should”
occur, and our simulations can provide some insight (although fundamentally,
the goal here it to motivate more detailed future work, as described below).
If accretion is fast ($t_{\rm acc}\lesssim t_{\rm dyn}$), then the mass will
be accreted from the sink radius (where it is gravitationally captured by the
sink in a resolved sense) onto the protostar rapidly, and there will not be
much reservoir or “PSD” to speak of, but if accretion is “slow” ($t_{\rm
acc}\gg t_{\rm dyn}$) then mass piles up in the reservoir by definition. One
might imagine this occurring if, for example, gas accreted into quasi-circular
orbits which form some kind of $\alpha$-disk around the true protostar,
accreting slowly onto said protostar with a very low $\alpha\ll 1$. As shown
in Fig. 20, if this occurred at the rates assumed in our “slow” accretion
scenario, it would imply a buildup of a reservoir or “PSD” mass $M_{\rm
PSD}\gg M_{\ast}$. So if this really were in a PSD and accretion were “slow,”
it seems likely such a system would fragment catastrophically into close
multiples.

To gain some insight into what is reasonable, recall (per § 3.5) that the
dynamical time in the CQM is $t_{\rm dyn}^{\rm
CQM}=1/\Omega=\sqrt{R^{3}/G\,M_{\rm BH}}\sim 5000\,{\rm yr}\,R_{\rm
pc}^{3/2}\,M_{\rm BH,\,7}^{-1/2}$ ($M_{\rm BH,\,7}\equiv M_{\rm
BH}/10^{7}\,{\rm M_{\odot}}$, $R_{\rm pc}\equiv R/{\rm pc}$), and the tidal
radius of a PSD+star system with system mass $M_{\rm s}\equiv M_{\ast}+M_{\rm
PSD}=M_{\rm s,\,\odot}\,{\rm M_{\odot}}$ is $r_{\rm t}=R\,(M_{\rm s}/4\,M_{\rm
BH})^{1/3}\sim 1000\,{\rm au}\,M_{\rm s,\,\odot}^{1/3}\,M_{7}^{-1/3}\,R_{\rm
pc}$, which gives a maximum PSD dynamical time of $t_{\rm dyn}^{\rm
s,\,max}\equiv\sqrt{r_{\rm t}^{3}/G\,M_{\rm s}}=1/2\,\Omega\sim 2500\,{\rm
yr}\,R_{\rm pc}^{3/2}\,M_{\rm BH,\,7}^{-1/2}$. Given the accretion criteria
adopted in our simulations (see § 2.1), it is implicitly required that any
accreted gas onto a a protostellar sink would have a circularization radius
less than $r_{\rm t}$ (since it must be bound to the sink, have apocentric
radius inside the sink radius, etc.). So the key question is whether or not
accretion in the PSDs is “dynamical” ($t_{\rm acc}\sim t_{\rm dyn}^{\rm
s,\,max}$) or “secular” ($t_{\rm acc}\gg t_{\rm dyn}^{\rm s,\,max}$). In
gravitationally-unstable, rapidly cooling, or strongly-magnetized PSDs, the
former (dynamical accretion) is usually true in idealized simulations (Gammie,
2001; Deng et al., 2017, 2020; Riols & Latter, 2016; Kratter & Lodato, 2016);
but in a traditional weakly-turbulent Shakura & Sunyaev (1973)-like
$\alpha$-disk PSD with $\alpha\ll 1$, the accretion is secular with $t_{\rm
acc}/t_{\rm dyn}\sim(c_{s}^{\rm PSD}/V_{c}^{\rm PSD})^{-2}\,\alpha^{-1}\gg 1$.

It is notable that our numerical sink accretion radii $r_{\rm sink}\sim
80-100\,$au are within $r_{\rm sink}<r_{\rm t}$ so long as $R_{\rm pc}\gtrsim
0.1\,M_{\rm s,\,\odot}^{-1/3}$ (and we show this explicitly below when we
discuss multiplicity), and our mass resolution is such that the more massive
PSDs would have as many as $\sim 10^{4}$ cell). So in principle, we should be
able to marginally resolve at least the formation of these PSDs or “stalled”
accretion reservoirs if the accretion were really “slow,” at least at larger
radii $\sim 1\,$pc. But we find that at any given time, it is very rare to
find any systems with extended resolved PSDs or quasi-hydrostatic accreting
gas “reservoirs,” in either the “fast” or “slow” $t_{\rm acc}$ simulations. We
see none in the few stars that form in our higher-resolution test in § 5. And
we show below that multiplicity is very strongly suppressed in the
environments here. This, coupled to the sink growth histories in Fig. 14,
suggests that at least on semi-resolved scales from $\sim 100-1000$ au, the
accretion is indeed “dynamical” given the physics in our simulations,666This
does not necessarily mean the accretion would be dynamical throughout the
entire life of the star and/or reservoir/PSD, as we are really only interested
in and able to follow the duration of the initial phase of rapid stellar
growth. A more depleted compact PSD or reservoir, with mass much less than the
star could of course persist much longer with much lower accretion rates, more
akin to typically observed PSDs in the Solar neighborhood with $\gtrsim
10^{5}-10^{6}\,$yr ages that have depleted to $M_{\rm PSD}\ll M_{\ast}$. and
that massive PSDs never really build up in the first place. In contrast, if
accretion were “slow” on resolved scales we would expect to see massive PSDs
in both simulations, but with an over-depleted “hole” in the center in our
“fast accretion” simulations.

That should not, perhaps, be surprising, for at least three reasons. First, as
noted in § 3 and discussed extensively in Paper II, the CQM is extremely
strongly-magnetized and well-ionized, with the colder gas accreting onto
protostars having typical $\beta$ as low as $\sim 10^{-4}$. Magnetic braking
should therefore produce efficient torques and prevent PSD formation, far
moreso than in “typical” (far less-strongly magnetized) GMC conditions.
Second, the early stages of even pure-hydrodynamic collapse follow more
closely a Shu (1977)-type quasi-spherical, quasi-isothermal collapse with
$\dot{M}_{\ast}\sim c_{s}^{3}/G$. Here, even if we neglect magnetic fields
$c_{s}$ is much larger through most of the CQM, and even in the coldest
phases, owing to the dense gas trapping radiation, with high radiation and gas
temperature floors more akin to the WNM, and dust being sublimated (see § 3),
producing larger accretion rates by a factor of $\sim 10^{4}-10^{6}$ compared
to a $\sim 10\,$K standard cold molecular medium in the Solar neighborhood.
Third, even if one tried to form a PSD, the massive PSDs of greatest interest
here would likely be highly unstable: for a PSD with most of $M_{\rm s}$ in
$M_{\rm PSD}$, with temperature in the PSD $T_{\rm PSD,\,gas}$ and size equal
to the tidal radius $r_{t}$, we would have $Q\sim 0.4\,(T_{\rm
PSD,\,gas}/1000\,{\rm K})^{1/2}\,R_{\rm pc}^{1/2}\,(M_{\rm s}/100\,{\rm
M_{\odot}})^{-1/3}\,M_{\rm bh,\,7}^{-1/6}$, with a cooling time at $r_{\rm t}$
of $t_{\rm cool}\ll\Omega^{-1}\sim t_{\rm dyn,\,PSD}$ (since the density of a
PSD at $r_{\rm t}$ is only a factor of $\sim Q$ larger than the midplane of
the CQM; see Paper II). So the PSDs would be violently gravitoturbulently
unstable without magnetic fields (Gammie, 2001; Forgan et al., 2017; Riols &
Latter, 2018).

It is important to note that all of these properties strongly differentiate
the circum-stellar PSDs around stars embedded in the QAD from circum-planetary
(or satellitesimal) disks around planets embedded in a “typical” circum-
stellar disk, a point we discuss further in § 4.7.

Regarding the strong magnetization, one might worry that non-ideal MHD effects
should lead to less efficient angular momentum transport (as commonly argued
for Solar-neighborhood PSDs). However, we note that (unlike STARFORGE
simulations of more typical Solar-neighborhood GMCs), we also do not see
prominent PSDs on these scales even in our simulations without magnetic
fields, discussed below. So it cannot be a purely magnetic effect. And in the
“diffuse” CQM (outside of PSDs) it is easy to verify that non-ideal MHD
effects are irrelevant: as shown in detail in Paper I, and discussed in § 3,
given the high temperatures ($\gtrsim 1000\,$K) of even molecular gas, extreme
irradiation and cosmic ray background in the vicinity of a bright quasar and
nuclear galactic starburst, the ionization fractions more closely resemble the
WNM and WIM ($\sim 1-100\%$ ionization fraction) so ideal MHD is actually an
excellent approximation. Moreover, as discussed in Paper I, we have
experimented with turning on and off different non-ideal MHD terms. Ohmic
resistivity and ambipolar diffusion have negligible effects here: the former
(Ohmic) because we never reach low enough neutral fractions and weak enough
fields for it to dominate among the non-ideal effects, the latter (ambipolar)
because even where it is the “most important” non-ideal MHD effect, the
relevant ambipolar diffusion timescales are vastly longer than the effective
turbulent diffusivity timescales (given the trans-Alfvénic and highly super-
sonic turbulence) on resolved scales here (consistent with a number of other