recent simulations of even more typical Solar-neighborhood GMCs; see Mac Low &
Klessen, 2004; Vázquez-Semadeni et al., 2011; Chen & Ostriker, 2014; Wurster
et al., 2021; Sadanari et al., 2022).

The Hall effect is slightly more complicated. We do see some places where Hall
MHD could in principle be important (at least in the sense that the fastest
Hall timescales – those associated with fast whistler waves – become shorter
than other resolved timescales like Alfvén or turbulent crossing times at our
resolution). However this only occurs in a very small subset of the gas in the
simulation, which has extremely small ion fractions $x_{i}\sim n_{\rm
ion}/n_{\rm neutral}\sim 10^{-17}-10^{-15}$ (vastly lower than the diffuse
CQM). But there are four numerical challenges involved here: (1) gas only ever
appears to reach these conditions on extremely small scales of order one
resolution element around sink particles (at the very highest densities before
being accreted), so coherent Hall MHD effects are not well-resolved; (2) this
usually involves very thermally-cold cells dropping to around the temperature
“floor” set by the CMB and other physics, but these cells always have plasma
$\beta$ extremely small (often $\ll 10^{-4}$), so accurately integrating their
thermal evolution is extremely challenging (very small integration errors in
their kinetic/magnetic energies can easily dominate over the physical thermal
energy); (3) these fast Whistler waves impose a very short timestep
requirement, so in our default simulations we choose to ignore this effect as
it limits our simulation runtime; and (4) it is well known that the Hall
effect is very sensitive to the thermo-chemistry and especially dust
properties of the PSDs (Bai & Stone, 2017; Zhao et al., 2020; Lee et al.,
2021; Tsukamoto et al., 2017), which are very different here compared to
“typical” PSDs in the Solar neighborhood, so it is by no means clear if our
chemical network is totally adequate for these extremes. So while it is
conceivable that these effects could be important in the most-dense PSDs which
form here, following Hall effects believably clearly requires much more
highly-resolved simulations.

We conclude that our “fast” sub-grid accretion model is probably more
reasonable than the “slow” model, at least at lowest-order. But saying
anything more detailed requires future work. In particular, this would be much
better studied by idealized simulations of a single collapsing core, in the
more extreme environments here. Critical for this problem as compared to more
“typical” Solar-neighborhood environments is to include (1) analytically (as
in e.g. a shearing-box simulation) the extremely strong tidal/shear
background, (2) the very strong fluxes not just of mass (with $\dot{M}$
reaching $\gtrsim 1\,{\rm M_{\odot}\,yr^{-1}}$ in the rapid collapse phase)
but also magnetic flux from the ambient ionized medium and thermal energy
(accounting for the accreting gas being warmer), (3) the much stronger local
interstellar radiation field and cosmic ray background (usually taken as
constant with some LISM values, but these can be scaled to a given location in
the QAD/CQM using the values at a given $R$ given in Paper I), and (4) some
model for dust sublimation/destruction which allows one to evolve a high-
metallicity but dust-poor ambient medium.

Figure 23: IMF for the “fast” (top) and “slow” (bottom) PSD accretion model
for all stars interior to the BHROI as Fig. 18 (left) and versus radius as
Fig. 21 (right), but in simulations where we repeat our hyper-refinement
experiment keeping all other physics identical but removing magnetic fields
from the initial conditions per § 4.5.6. These are normalized by $N_{s}^{-1}$:
the absolute SFR and number of stars is $\gg 100$ times larger at $r\ll 1$ pc
compared to the runs with magnetic fields. Per Paper I and expectations from §
3, artificially removing magnetic fields qualitatively changes the properties
of the QAD+CQM, leading to runaway fragmentation on sub-pc scales, orders-of-
magnitude higher SFRs, and a much smaller, thinner, less-massive QAD with a
factor $\sim 500$ lower accretion rate onto the SMBH. But despite the
enhancement of global fragmentation and star formation and lower gas masses in
the QAD, the predicted IMF actually becomes much more top-heavy absent
magnetic fields (in both PSD accretion-rate scenarios). Moreover, the
predicted IMFs become monotonically more top-heavy at larger radii $r$
interior to the BHROI (perhaps simply limited by the small QAD mass as
$r\rightarrow 0$ in these simulations). This demonstrates the critical role of
magnetic fields suppressing not just the global fragmentation/SFR (discussed
in Paper I), but also the masses of individual cores and their ability to
accrete.

####  4.5.6 Dependence on Magnetic Fields

In Paper I & Paper II, we compared a re-run of this simulation without
magnetic fields, and showed that this produces radically different global
properties. In particular, the suppression of SF and gravitoturbulent
fragmentation of the QAD depends strongly on the magnetic fields. In the full-
physics simulation, magnetic pressure dominates the CQM and QAD pressure, and
so removing magnetic fields leads to a much thinner, less massive, smaller-
spatial-scale, rapidly gas-depleted QAD, with the SFR at CQM scales larger by
factors $\sim 10^{3}-10^{5}$, and the accretion rates onto the SMBH decreasing
by factors of $\sim 10^{2}-10^{3}$. So clearly magnetic fields play a crucial
role in the global dynamics on these scales.

In Fig. 23, we show how removing MHD influences the predicted IMF in these
simulations. Even on timescales as short as a few hundred years after the
highest refinement level is reached (shorter than the duration of our full-
physics run at its highest resolution), we see that absent magnetic fields,
the IMF “runs away” to much more top-heavy values. Even with the “fast” PSD
accretion model (which produces a mildly top-heavy IMF in the full-physics
case), the IMF absent magnetic fields shifts to resemble a second lognormal
peak at $\gtrsim 10\,M_{\odot}$. Also important, without MHD the IMF becomes
similarly top-heavy at all radii where we can predict it (where we employ full
STARFORGE physics). Recall, in the full-physics case, the IMF becomes
progressively less top-heavy at larger radii, resembling a “typical” IMF in
the more ambient ISM at $R\gg 1\,$pc. Absent MHD, even at $>\,$pc scales, the
IMF appears to be getting systematically more top-heavy at larger-$r$.

The fact that magnetic fields in our simulations help to suppress global
fragmentation and star formation rates overall is qualitatively consistent
(albeit more dramatic, owing to much stronger fields) with most previous
studies of star formation in Solar-neighborhood-like ISM conditions (see
Krumholz & Federrath, 2019, and references therein). However, studies of the
IMF in LISM conditions (including those with the same STARFORGE physics used
here) have generally concluded that magnetic fields have a relatively weak net
effect on the IMF of the stars that do form (see Myers et al., 2014; Lee &
Hennebelle, 2019; Krumholz & Federrath, 2019; Sharda et al., 2020), with the
sign of the net effect varying as magnetic fields can both suppress
fragmentation (producing larger “cores”) and also accretion (suppressing
growth after initial collapse). But our results are more similar to Hocuk et
al. (2012), who considered idealized simulations of a $\sim 800\,M_{\odot}$
cloud at a distance of $\sim 10$ pc from a $\sim 10^{7}\,M_{\odot}$ BH,777We
do caution that while the tidal environment considered is similar, the
simulations in Hocuk et al. (2012) differ fairly radically from those here (§
3) in that their simulations considered orders-of-magnitude smaller mass,
density, column density, temperature, and magnetic field values. and showed
adding strong magnetic fields shifted their predicted turnover mass from $\sim
1\,M_{\odot}$ to $\sim 0.1\,M_{\odot}$ (see also Li et al. 2010). Those
authors noted that under these conditions, the magnetic fields strongly
suppress accretion, and can even lower the thermal Jeans mass while also
suppressing fragmentation by restricting fragmentation to a limited set of
more extreme over-dense regions (with higher $\rho$ so lower $M_{\rm Jeans}$)
– both effects we also see and describe above. The more extreme difference we
see between our magnetized and non-magnetized runs likely owes to the more
extreme magnetic field strengths.

Figure 24: Resolved multiplicity statistics (§ 4.5.7). Top: For the “slow” PSD
accretion simulation, we plot the multiplicity fraction $f$ of stars versus
primary sink mass $M_{1}$, including all mass ratios, selecting stars with
formation distances from the SMBH $r$. We compare LISM measurements compiled
in Offner et al. (2023), similar to STARFORGE simulations of typical LISM GMCs
(Guszejnov et al., 2023). Bottom: Semi-major axis ($a$) distribution of the
same multiples (given the low $f$, these are almost all binaries). Arrows
indicate a typical 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})$
for a binary (of mass $M_{\rm s}$ equal to the median mass of the binaries in
each line) in the SMBH potential. Note that the “tails” of stars with larger
$a$ largely correspond to more-massive stars or stars near the outer edge of
the radii ranges, for which $r_{\rm t}$ is larger. We plot LISM observations
from Raghavan et al. (2010) after completeness correction in Moe & Di Stefano
(2017). Shaded range shows sink gravitational softening radii (harder binaries
are un-resolved). At radii $\gtrsim\,$pc where the tidal effects of the SMBH
become relatively weaker, multiplicity statistics begin to approach the LISM
results (though still suppressed), but at smaller radii, multiplicity is very
strongly suppressed. The suppression occurs both at formation (fewer “common
core” binaries form), and then via subsequent evolution/tidal disruption. We
see slightly older stars, even $>300\,$yr or $\sim 3\,\Omega^{-1}$ at
$r\lesssim 0.1\,$pc, have even lower multiplicity, and the semi-major axis
distribution is truncated at the tidal radii as expected. Trends are similar
for the “fast” accretion case.

####  4.5.7 Multiplicity

We see strongly suppressed multiplicity in the QAD environment. This is true
whether we attempt to visually identify binaries or use more sophisticated
linking algorithms such as those in Guszejnov et al. (2023). We emphasize this
is quite different from simulations with the same physics and code in
“typical” Solar-neighborhood ISM GMCs (Bate, 2012; Lee et al., 2019; Guszejnov
et al., 2023; Kuruwita & Haugbølle, 2023). But it is not surprising: even if
binaries can form, unequal-mass stellar binaries which have a separation
comparable to or larger than the 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})$ will be tidally sheared apart in these
extreme environments. And binary formation, at least “common core” binaries,
should be suppressed in these environments by multiple effects: (1) the warm
gas temperatures and strong magnetic fields raise the Jeans/magnetic critical
mass as discussed above and in § 4.6 below; and (2) the strong magnetic fields
and high ionization fraction lead to rapid accretion/magnetic braking and
inhibit formation of PSDs (§ 4.5.5).

Fig. 24 shows this quantitatively: we plot the multiplicity fraction, defined
using the algorithm to identify all bound multiples (including hierarchical
multiples) from Bate (2012) as applied in Guszejnov et al. (2023).888Note that
we do not apply the additional mass-ratio and age cuts considered in Guszejnov
et al. (2023). The age cuts of $\gtrsim$ Myr would not make sense given the
timescales of interest here; the mass-ratio cuts make little difference to the
result but will slightly decrease the multiplicity fractions if we did apply
them. First, we show the multiple fraction (defined here simply as the
fraction of stars in bound multiples),999Because the multiplicity is low, most
of the multiples are simple binaries, not higher-order multiples, so it makes
little difference here if we consider a binary fraction versus some higher-
order multiplicity fraction statistics as discussed in Guszejnov et al.
(2023). versus stellar mass, for stars at various BH-centric radii $R$ and
with various ages since the initial sink formation time. For reference, we
also plot the observational compilation in the LISM with which the simulations
in Guszejnov et al. (2023) agreed. Recall, these used the same STARFORGE
physics to simulate “typical” Solar neighborhood GMC/LISM environments, and
therein there, the multiplicity fraction rises from $\sim 0.5$ at $\sim
1\,M_{\odot}$ to $\sim 0.9$ at $>10\,M_{\odot}$, much larger than we see.

Clearly there is some radial trend where at larger radii, the simulations
begin to approach the “typical” ISM, though still with a notable suppression,
which again is expected from the simple tidal radius argument above. We also
see clear evidence for this in the semi-major axis ($a$) distribution: at e.g.
$R<0.1\,$pc, considering binaries $>300\,$yr old, the only surviving (small
number) of examples have semi-major axes $a<100\,$au, more or less exactly the
tidal radius predicted for this $R$ and typical stellar masses $\sim
10-30\,M_{\odot}$ in the simulation at those radii (this is the example with
“slow” PSD accretion). Even at $R\sim 1$ pc, the semi-major axis distribution
is sharply suppressed at $a$ larger than a few thousand au, again roughly the
tidal radius at these BH-centric radii. For comparison, the observed semi-
major axis distribution in the LISM (for Solar-type and massive primaries) is
flat up to at least $a\sim 10^{4}$ au.

So clearly some of the suppression of multiplicity owes to the strong tidal
forces, but this is probably not all of it. If we simply take the Solar-
neighborhood binary fractions and truncate them by eliminating any binaries
with $a\gtrsim r_{\rm t}$ (as well as removing any which would be numerically-
unresolved with $a\lesssim 20\,$au), we would still predict a “surviving”
binary fraction of $\sim 20-30\%$ in the massive $\gtrsim 10\,M_{\odot}$
stars, for example. But the simulations predict a binary fraction $\lesssim
5-10\%$. This additional suppression of stellar multiplicity could occur if a
substantial fraction of binaries at $\lesssim 20$au separation form through
fragmentation at larger separations and then migrate inwards through gas or
dynamical friction (Offner et al., 2016; Lee et al., 2019; Kuruwita &
Haugbølle, 2023). Here, such binaries would be tidally disrupted before
migration could occur. Alternatively, the remaining suppression owes to the
physics described above suppressing the incidence of fragmentation (both core
and disk).

Briefly, we do not see any obvious unusual behavior in the mass ratio
distribution, but given the combination of small statistics and resolution
limitations (and incompleteness for observational mass-ratio distributions),
this is not a particularly significant statement.

Of course, we cannot resolve compact binaries that could form in principle via
fragmentation in some dense PSD or collapsing cores at scales below $\sim
10\,$au in these simulations. Those are particularly interesting for e.g. TDE
and LIGO/LISA sources as noted in § 1. Predictions for such compact binaries
will fundamentally require dedicated, much higher-resolution simulations of
the formation of individual systems with the extreme boundary/initial
conditions provided by the simulations here, as argued for in § 4.5.5.

###  4.6 Physics: Why Is the IMF Top-Heavy?

What physics drives the IMF to be top-heavy in the simulations?

Figure 25: Different characteristic mass scales versus distance to the SMBH (§
4.6), defined as the mass-weighted median of said mass scale at the end of the
“fast” PSD accretion simulation, in radial annuli. Using the mean or other
weighting, or comparing different times, or the “slow” PSD simulation, makes
little difference. We compare the total enclosed gas mass $M_{\rm gas}(<r)$
and a number of other mass scales which might set the IMF mass scale. None of
these alone appears able to entirely explain the IMF trends seen in the
simulations § 4.6.1. However some combination may be able to do so for a
couple scenarios discussed in § 4.6.2-4.6.3.

####  4.6.1 Scalings Which Cannot Alone Explain the Simulation IMF Trends

First, it is worth discussing many candidate explanations that can be ruled
out as the sole explanation for the IMF variation we see. In Fig. 25, we show
the typical values of many different characteristic mass scales, which one
might imagine setting the IMF in principle. None of these exhibits the
behaviors needed to fully explain the variation in IMF behavior we see, but we
briefly discuss them in turn. To do so, it is helpful to divide them into a
few different categories.

(1) Jeans-Type Masses: First, one can consider some Jeans-type length
$\lambda_{J}\sim v_{0}/\sqrt{G\,\rho}$ and corresponding Jeans mass
$M_{J}\sim(\lambda_{J}/2)^{3}\,(4\pi/3)\,\rho$, which have been argued to be
proportional to the IMF mean or turnover mass (Bate & Bonnell, 2005). One can
then generalize the classic thermal Jeans $v_{0}\rightarrow c_{s}$ mass, to a
magnetic Jeans $v_{0}\rightarrow v_{A}$, or turbulent Jeans
$v_{0}\rightarrow\langle\delta v_{\rm turb}^{2}(\lambda)\rangle^{1/2}$ mass.
We see that the magnetic and turbulent Jeans masses are vastly too-large at
all radii to represent some mean mass in the IMF – indeed, they exceed the
total disk mass $M_{\rm gas}(<r)$, which means such collapse would not
generically occur. And the thermal Jeans mass is much too high at large $r$
and increases with $r$ (the opposite of the IMF characteristic mass). We can
also define the magnetic critical mass $M_{B}\sim 0.1\,\Phi_{B}/G^{1/2}$. If
we assume $\Phi_{B}\sim|{\bf B}|\,\pi\,R^{2}$,101010For the actual plot we
calculate a more appropriate effective $\int{\bf B}\cdot{\rm d}{\bf A}$ with
the projection area defined as the normal to the mean ${\bf B}$ field in
cylindrical annuli. then $M_{B}\sim M_{\rm gas}\,(M_{J,\,B}/M_{\rm
gas})^{1/3}$, so this is also much too large ($>M_{\rm gas}$) at all radii of
interest, and it (like $M_{J,\,T}$) increases with $r$. It is also well-known
that these hypotheses would predict vastly more-dramatic IMF variation in the
Solar neighborhood compared to either that observed or seen in simulations
with the full physics (e.g. stellar feedback) here (Guszejnov et al., 2017b,
2019).

(2) Opacity/Heating-Limited Masses: We can also consider the classic opacity
limit (Rees, 1976) and related arguments, where a critical mass scale is set
by some point where heating balances cooling and effective cooling times
become long compared to collapse times, preventing fragmentation. If we just
take the naive scalings for the LISM, assuming opacities are dominated by dust
only and that compression dominates heating, as in Jappsen et al. (2005); Bate
(2009); Lee & Hennebelle (2018), we obtain the opacity limit (LISM) mass. As
shown in Grudić & Hopkins (2023), one can more carefully generalize this to
cases with arbitrary gas-phase plus dust-phase opacities and for the limit
where external radiative heating dominates over compressional heating – using
these scalings (e.g. Eqs. 19 & 21 therein) with the exact compressional and
radiative heating rates and opacities taken directly from the simulations
gives the opacity limit (consistent) mass. One can then scale to a first
Larson core mass as shown therein (Eq. 23 therein) assuming quasi-adiabatic
collapse past the limit of inefficient cooling (Lee & Hennebelle, 2018).
Alternatively, but related in a physical sense, some models predict that the
IMF turnover mass is strictly self-regulated by protostellar feedback from
radiation (Krumholz, 2011; Guszejnov et al., 2016). Alone, these models
predict an extremely weak dependence of the turnover mass on local conditions
such as $T$ and $\rho$, as e.g. $M_{\rm PSFB}\propto(\rho\,T)^{-1/18}$. Most
of these masses scale far too weakly with radius $r$ to explain the IMFs in
our simulations, moreover they are much smaller than the mass scales where we
see interesting IMF variation. The one potentially interesting case is the
consistent opacity limit mass – this rises at small radii owing to the
rapidly-rising gas temperatures and gas-phase opacities (even though the dust
is sublimated, the gas opacities can be quite large here), combined with the
stronger sensitivity to conditions in the radiation-dominated regime as
derived in Grudić & Hopkins (2023). But given the (again very low) mass scales
of (e.g. $\sim 0.03\,{\rm M_{\odot}}$ at $r\sim 0.1\,$pc), this is likely at
most relevant for explaining the weak variation in the very low-mass end of
the IMF.

(3) Global Fragmentation Mass Scales: There are some characteristic global
mass scales in fragmentation cascades. One is the total mass $M_{\rm
gas}(<r)$. This and similarly any “tidally limited” mass, e.g. the classic
Roche or Hill mass (making the analogy of star formation in the QAD around the
SMBH to planet formation in a PSD around a single star), $M_{\rm
Hill}\sim\Sigma_{\rm gas}\,H^{2}$, all scale in the opposite manner desired,
increasing monotonically with radius $r$. Alternatively the Toomre mass scales
as the Hill/Roche/Tidal mass multiplied by $\sim Q^{2}$ – this makes it vastly
larger than $M_{\rm gas}(<r)$ and the characteristic IMF scales.
Unsurprisingly, on large scales these have more to do with the mass scales of
GMCs/star-forming cloud complexes than individual stars. In a turbulent
cascade, there are two additional commonly-cited mass scales. First is the
“sonic mass” or turbulent Bonner-Ebert sphere mass $M_{\rm sonic}\sim
c_{s}\,R_{\rm sonic}/G$ (where $R_{\rm sonic}$ is the scale where the
turbulence becomes sub-sonic; Padoan & Nordlund 2002; Hopkins 2012b, 2013a,
2013c; Guszejnov & Hopkins 2015) – this is fairly constant around $\sim
1\,M_{\odot}$ at large scales but again increases with $r$ from $r\sim
0.01-1\,$pc. Second is the “maximally-unstable turbulent scale” defined by the
first-crossing distribution (Hopkins, 2012a, 2013b) (Eq. 44, combined with 16
& 24 in Hopkins 2013b), analogous to the Hill/tidal scale but in an
inhomogeneous turbulent medium (and identical to these and the “turbulent
Toomre mass” in a supersonically-turbulent, $Q_{\rm turb}\sim 1$, weakly-
magnetized medium). Again this increases with $r$ (and at the largest radii is
much too large to represent the IMF).

(4) Accretion-Based Models: We have extensively discussed (§ 4.3, 4.3.1,
4.3.2, 4.5.6) how accretion is strongly suppressed in these stars, beyond
their initial collapse masses. So they are accreting less, not more, from
their environments than comparably massive stars in typical GMCs. Thus any
explanation for the IMF relying on significant accretion beyond the initial
collapsing core (e.g. competitive accretion) is immediately ruled out.

(5) Metallicity-Driven or Strictly Cooling-Rate-Driven Models: Various models
have proposed heuristically that the IMF could depend on some simple power-law
function of metallicity, or otherwise on the ratio of cooling to dynamical
time (Cappellari et al., 2012; Spiniello et al., 2012; Narayanan & Davé, 2013;
Martín-Navarro et al., 2015), although more detailed assessments such as
Offner et al. (2014) have noted that observations require any such dependence
be weak and/or highly indirect (as in e.g. the opacity limit, which can depend
on metallicity, but very weakly, and depends on other parameters as well).
However, as shown in Paper I, the entire galactic nuclear region for which we
model the IMF is quite chemically homogeneous here: there is no appreciable
difference in abundances, nor in the ratio $t_{\rm cool}/t_{\rm dyn}\ll 1$,
over the dynamic range where we see strong IMF variation.

####  4.6.2 A Potential Explanation for the IMF in the “Slow” Accretion Limit:
(Gravito)Turbulent Fragmentation in a Quasi-Stable, Quasi-Keplerian, Trans-
Alfvénic QAD

One possible explanation for the behavior of the IMF in our “slow” accretion
simulations is suggested by the remarkable similarity between the IMF shape
where it is top heavy, and the IMF of planets formed by gravito-turbulent
fragmentation in a quasi-Keplerian disk in Hopkins & Christiansen (2013)
(after, of course, scaling the absolute units of the disk to the QAD here).
That study considered turbulent fragmentation (generalizing models for star
formation in the ISM from e.g. Hennebelle & Chabrier 2009; Hopkins 2012b), in
a Keplerian disk, with a quasi-isothermal EOS ($t_{\rm cool}\ll t_{\rm dyn}$),
modest thermal Toomre $0.1\lesssim Q\lesssim 100$, and modestly compressible
turbulence (compressive Mach number of $\sim 1$ producing $\mathcal{O}(1)$
local density fluctuations in shocks). While subsequent observations have
argued that these assumptions probably do not apply to most LISM PSDs
(Flaherty et al., 2015, 2017; Andrews et al., 2018), they are an excellent
description of the QAD (see Paper I & Paper II). Moreover, unlike in LISM GMC
star formation, we can safely ignore the effects of stellar feedback both
locally and globally regulating SF in this particular simulation.

Essentially, what happens in this model (and can be more generally derived
from the models in Hopkins (2013b)), is as follows. As discussed above (§
4.6.1) there are two key turbulent fragmentation scales: the sonic mass
$M_{\rm sonic}$ and maximally-unstable turbulent mass $M_{\rm max}$, which
define the potential dynamic range of the turbulent fragmentation cascade.
When “effective” compressive Mach number $\mathcal{M}_{c,\,\rm eff}$ (defined
in practice as directly related to the magnitude of rms density fluctuations
generated by turbulence, $\delta\rho/\rho\sim\mathcal{M}_{c,\,\rm eff}$) is
large ($\mathcal{M}_{c,\,\rm eff}\gg 1$), then the entire fragmentation
cascade is realized: objects form initially at the first-crossing $M_{\rm
max}$ scale and hierarchically fragment rapidly into progressively smaller
sub-structures (e.g. clouds/clumps/cores, etc.) producing scale-free mass
functions and correlation functions all the way down to the last-
crossing/sonic scale $M_{\rm sonic}$ (see e.g. Guszejnov et al. 2018a), with
the dynamic range between the two scaling as $\sim\mathcal{M}_{c,\,\rm
eff}^{4}/Q$. In the diffuse ISM, this can qualitatively explain the MFs of
GMCs/clumps/cores and the stellar IMF down to the mean/turnover mass $\sim
M_{\rm sonic}$ (see references above). But when $\mathcal{M}_{c,\,\rm eff}$
becomes $\lesssim 1$ with $Q_{\rm thermal}+Q_{\rm magnetic}\gtrsim 1$, as
considered in Hopkins & Christiansen (2013) the IMF becomes strongly-peaked
around $M_{\rm max}$, because essentially all mass scales are formally
“stable” if the disk were laminar, but turbulent fluctuations mean there is a
small (but non-vanishing) probability of a sufficiently overdense fluctuation
which can collapse and that probability is (by definition) maximized at
$M_{\rm max}$ (it is generically more difficult to generate a fluctuation with
sufficient mass to collapse on progressively smaller scales). So in practice,
we might expect the mean sink mass to “jump” at this critical point from
$M_{\rm sonic}\sim 1-3\,{\rm M_{\odot}}$ at $r\gtrsim\,$pc to $M_{\rm max}\sim
30-100\,{\rm M_{\odot}}$ in Fig. 25 at $r\ll\,$pc, which agrees at least
crudely with the sharp change in behavior seen at the outer QAD (Fig. 21).

This transition occurs fairly abruptly because of the same physics discussed
in § 4.2. At radii outside the QAD in the CQM, the turbulence is super-
Alfvénic, highly super-sonic, and the magnetic fields are isotropically
tangled (see Paper I & Paper II). In this regime, as discussed at length in
Hopkins & Christiansen (2013) and Hopkins (2013b) and shown in idealized
numerical simulations in Kowal et al. (2007); Lemaster & Stone (2009); Kritsuk
et al. (2011); Molina et al. (2012), magnetic fields only modify the effective
compressive Mach number by a geometric pre-factor $\mathcal{M}_{c,\,\rm
eff}\rightarrow b\,\mathcal{M}_{\rm sonic}/\sqrt{3}$ (with $b\sim 1/2$) as
$\beta\rightarrow 0$, because the fields cannot resist compressive motions
along the field lines. But in the QAD, the turbulence rapidly drops to sub-
Alfvénic (Paper I) and the field becomes dominated by the mean toroidal field,
which strongly suppresses the compressibility of the gas (suppressing
$\mathcal{M}_{c,\,\rm eff}$ by a factor $\sim c_{s}/v_{A}\sim\beta^{1/2}$ in
the perpendicular radial/vertical directions; see Hopkins 2013b § 6 for a
theoretical discussion or Beattie et al. 2021 for simulation examples), while
of course strong shear in the QAD prevents azimuthal compressions (suppressing
$\mathcal{M}_{c,\,\rm eff}$ in that direction by a factor $\sim
c_{s}/\kappa\,H\sim c_{s}/v_{A}$ for the magnetically-supported disks here,
see Hopkins & Christiansen 2013).111111Another way of saying this is that the
sonic mass/length defined purely in terms of the sonic Mach number
$\mathcal{M}_{s}$ is no longer relevant in the QAD, because in turbulent
fragmentation models this is just a proxy for the scale below which large
compressions are rare/difficult, so the “effective” sonic mass jumps to a
value $\gtrsim M_{\rm max}$. So $\mathcal{M}_{c,\,\rm eff}$ effectively drops
to $\sim\mathcal{M}_{A}\lesssim 1$, and the IMF changes abruptly.

####  4.6.3 A Potential Explanation for the IMF in the “Fast” Accretion Limit:
Feedback-Regulated Stellar Masses

In the “fast” accretion model, the behavior of the IMF is somewhat more
subtle. As discussed in § 4.5.2, the key difference is that in this case,
stellar feedback – primarily in the form of radiation from
accretion+contraction luminosity and protostellar jets, which scale with the
sub-grid mass and accretion rate onto the “star” within the sink – is able to
act because mass can transfer from sink to “star” rapidly compared to the
local QAD/CQM dynamical time.

At large radii in the CQM, the physics is probably similar to the LISM. The
sonic mass may play a role in setting the “turnover” mass (Guszejnov et al.,
2020), but this is modified by feedback: massive stars in particular self-
regulate via feedback making the high-mass slope less steep (Grudić et al.,
2022), and as many studies have shown protostellar jets shift the IMF peak
mass a factor of a few121212We stress that this is not “put in by hand”: the
jet mass-loading is just $\sim 30\%$ of the accreted mass (Grudić et al.,
2021). But as shown in these studies the jets efficiently entrain material
which would otherwise accrete onto the protostar. Here, while we showed the
jets have little effect on the global properties of the QAD/CQM, they by
definition can unbind a significant mass from accreting onto the protostar
since they are assigned a launch velocity comparable to the escape velocity
from the surface of the star, which of course must be larger than the escape
velocity from the local potential perturbation caused by the star itself on
larger (resolved) scales at $\sim 10-1000\,$au distances (even if it is less
than the escape velocity from the global potential of the SMBH, so the
outflows are not unbound from the global QAD/CQM. lower (Guszejnov et al.,
2021, 2022b; Lebreuilly et al., 2023), naturally explaining the difference in
the IMF shapes and peak masses between the “fast” and “slow” accretion models
(Figs. 18-21). Any modest residual difference between the IMF shape at
$\gtrsim$ pc scales and the “universal” IMF can be largely explained by the
suppression of fragmentation on smaller scales (§ 4.5.7), so the predicted IMF
is more like the system IMF (Fig. 20) – but we caution that we need higher
resolution to make that a more rigorous and quantitative statement (§ 5). We
do not see much evidence at the lowest masses for a role of a radially-
dependent opacity limit as in Fig. 25 (if anything, the “cutoff” mass of the
IMF may be slightly larger at the largest $r$, opposite this expectation), but
if the “first Larson core mass” is more relevant than the “opacity limit mass”
the weak evolution in Fig. 25 is mostly cancelled out and it would explain the
weak evolution in the lower-limit or “cutoff” mass at the lowest resolved
masses.

But clearly, something still happens around the outer QAD (see Fig. 21). At
that specific radius, the IMF is comparably top-heavy and shallow in slope to
the “slow accretion” simulation, though with perhaps the peak mass shifted
down by a factor of a few. The shift between “fast” and “slow” models again
can be attributed to feedback. But the stark “jump” to a more top-heavy shape
of the IMF in the outer QAD likely has a similar origin to that discussed in §
4.6.2 – this is a special environment where the effective compressibility of
the gas is strongly suppressed. This overall suppresses star formation, but
means that where stars do form, they are most likely to be at the largest
(most unstable) scales. This is further supported by the fact that we see the
sites of star formation in the outer QAD are the same for the the “fast” and
“slow” cases (§ 4.2).

The one behavior in the “fast” run which is somewhat more difficult to explain
with the arguments above is the apparent reversion to a “normal” IMF at
slightly smaller radii $0.02<r<0.05\,$pc, just before SF is completely
suppressed at smaller radii in the QAD (see Fig. 21). But as discussed in §
4.5.3, the statistics for this region are quite sparse (there are just a few
stars forming here) so it is difficult to draw a robust conclusion. But this
appears to be worth future investigation.

###  4.7 Is Star Formation in the QAD Analogous to Planet formation in a
“Typical” PSD?

It is common to draw the analogy between star formation in a QAD and planet
formation in “typical” PSDs (at least planet formation via gravitational
instability, as opposed to solid core formation), and indeed we did this to
some extent in § 4.6. As argued there, some parts of that analogy are
reasonable – after all, gravity is strictly scale-free, and the mass ratios
are conceivably similar in the QAD-star and PSD-planet cases. And indeed some
key aspects of gravito-turbulent fragmentation appear to cross-apply, at least
with a careful accounting for the correct “effective” compressive Mach numbers
(after accounting for the effects of magnetic fields and radiation
explicitly), as discussed in § 4.6.2.

However, most of our discussion above should make it clear that the analogy
should not be pushed too far. Much of the physics which plays a crucial role
(e.g. the thermo-chemistry) is not scale-free. And much of the physics that is
scale-free in principle (e.g. ideal MHD) occupies a qualitatively different
dimensionless parameter space (see § 3). Most obviously, the plasma $\beta$ in
the QADs here is $\beta_{\rm QAD}\sim 10^{-4}$, whereas in PSDs it is
generally assumed to be $\gg 1$ and in standard models for MRI-active PSDs,
the effective viscosity parameter $\alpha\sim\beta_{\rm PSD}^{-1}$ implying
$\beta_{\rm PSD}\gtrsim 10^{4}-10^{6}$ in observations (see Hughes et al.
2011; Flaherty et al. 2015, 2017; Andrews et al. 2018). Moreover in the QAD,
the magnetic fields are extremely well-coupled to the gas and well in the
ideal limit, with ionization fractions (given the warm temperatures and
extreme radiation environments) rising from $f_{\rm ion}\sim 0.01$ in the
outskirts to $\sim 1$ in the inner QAD. In contrast typical PSDs are expected
to have $f_{\rm ion}\ll 10^{-15}$ throughout much of the PSD, making magnetic
fields as a whole far less dynamically important and making Hall MHD and Ohmic
resistivity prominent (Wardle, 2007; Keith & Wardle, 2014; Tsukamoto et al.,
2015; Bai, 2017; Zhao et al., 2020; Lee et al., 2021).

Similarly, the QADs here are vastly more strongly turbulent even in a
dimensionless sense compared to observed PSDs. Typical sonic Mach numbers in
the QAD are $\mathcal{M}_{s}\sim 30-100$, and even the (much lower) Alfvén
Mach numbers are $\mathcal{M}_{A}\sim 1$, while for PSDs inferred Mach numbers
are $\mathcal{M}_{s}\sim 0.001-0.01$ (Hughes et al., 2011; Flaherty et al.,
2015, 2017; Andrews et al., 2018).

Related to this, the cooling physics is qualitatively different, not just in
the dominant cooling channels and opacities (though these are fundamentally
different, with dust sublimated in the QAD and orders-of-magnitude stronger
radiation fields and cosmic ray densities), but in the dimensionless ratio of
cooling time $t_{\rm cool}$ to dynamical time $t_{\rm dyn}$. For a PSD, except
perhaps for a very short transient period when they first form, one
generically expects (and again, observations appear to support) $t_{\rm
cool}\gg t_{\rm dyn}$, while for the QADs here cooling is efficient with
$t_{\rm cool}\ll t_{\rm dyn}$. This is one of the many causes of strong
gravito-turbulence and hence higher Mach numbers in the QADs, but also means
that the QADs can be multi-phase and feature radiative (quasi-isothermal)
shocks, and qualitatively different behaviors upon contraction.

Yet despite the much faster cooling in this sense in the QADs, the absolute
temperatures of the QAD are significantly warmer owing to the strong accretion
heating and extreme environment. Those higher QAD temperatures and accretion
rates also mean that the QADs feature orders-of-magnitude higher ratios of
radiation pressure to thermal gas pressure: in typical PSDs the radiation-to-
gas-thermal pressure ratios are $\sim 10^{-7}$ (e.g. Chiang & Youdin, 2010),
while in the QADs here they exceed unity throughout much of the QAD (and even
in the outermost QAD never drop much below $\sim 1$).

There is also the effect of (proto)-stellar feedback, e.g. in the form of
protostellar jets on these timescales, which even here (where it cannot
totally disrupt the background QAD/CQM) we showed has a qualitatively
important effect regulating the masses of the forming objects (§ 4.5.2). These
jets are much stronger and carry much more mass (fractionally speaking) than
what is generally assumed for planets (although of course the jet efficiencies
could vary with environment and magnetic field strengths).

And obviously, even if dust were not sublimated in the QADs, the Stokes
numbers are orders-of-magnitude different and efficient coagulation will not
occur as it does in PSDs, so there is also no direct analogy for planetesimal
formation, pebble accretion, the streaming instability (though other resonant
drag instabilities do exist in QAD environments, they behave in qualitatively
different ways; see Squire & Hopkins 2018b, a; Soliman & Hopkins 2022), or
core accretion.

These differences will also strongly modify the dynamics of the extant stars
and black holes in the QAD, via their coupling to the gaseous disk (for one
example, see § 4.3.1). But even for a pure $N$-body disk, i.e. after the gas
disk is fully consumed or expelled, there are some subtle but important
distinctions. In a planetary disk, the deviations from Keplerian orbits are,
to leading order, dominated by an $\mathcal{O}(1)$ number $N$ of large bodies
(e.g. Jupiter), while for a QAD (assuming the gas has been removed) they are
dominated by $N\gg 1$ broadly similar-mass stars and/or compact objects.
Moreover for the planetary disk, the total mass of orbiting bodies is
generally much less than the mass of the host star, while it is basically
guaranteed that in the outer CQM (approaching the BHROI), the total enclosed
mass of stars becomes comparable to the SMBH mass. These differences will
strongly influence and often qualitatively change the outcomes of phenomena
such as migration, resonant relaxation, propagation of $m=1$ modes and other
perturbations, and resonance and trojan trapping.

##  5 Potential Effects of Resolution

We discussed the typical simulation numerical resolution in Fig. 3 and § 2. Of
course, ideally one would always like to push to even higher-resolution.
However given the extremely-small timesteps of $\sim$ days involved in these
simulations (owing to the short dynamical and radiation timescales), it is not
possible to do so while running for even a modest fraction of the (already-
limited) run-time at highest resolution ($\sim 10^{4}\,{\rm yr}$; see Fig. 3).
It is also worth stressing that our “default” resolution is already
significantly finer than any of the historical, idealized circum-SMBH IMF
studies referenced in § 1.

More extensive STARFORGE resolution tests for LISM conditions (where such
tests are feasible) are presented in Guszejnov et al. (2021); Grudić et al.
(2021); Guszejnov et al. (2022b). These studies argued that the predicted IMFs
(the subject of our discussion in § 4.5) were reasonably well-converged for
${\rm M_{\odot}}\gtrsim 100\,\Delta m$ ($100\times$ the mass resolution), or
$\gtrsim 0.5\,{\rm M_{\odot}}$ here at the median resolution in the QAD/CQM.
And we stress that the global, integrated quantities in these simulations
which we study (e.g. the total mass of sinks, star formation efficiency,
locations within GMCs where most stars formed, bounded-ness of the stars,
etc., and the subjects of our discussion in § 3-4.4) are converged in these
tests even at much poorer resolution than our fiducial simulation, as low as
$\sim 10^{-2}-10^{-1}\,M_{\odot}$, suggesting the simulations should be robust
in these respects.

For the IMF specifically, our primarily focus on the “top-heavy” behavior in §
4.5 was at much larger masses going up to $\sim 100\,M_{\odot}$, where the key
behaviors (and e.g. Salpeter slope) are reasonably well-converged in the LISM-
like resolution tests in Guszejnov et al. (2021); Grudić et al. (2021);
Guszejnov et al. (2022b). Moreover, the offsets we study in the IMF (even
around the $\sim 0.1-1\,M_{\odot}$ mass range) are generally much larger than
the effects from incomplete convergence at this ($\sim$ few $\times
10^{-3}\,{\rm M_{\odot}}$) resolution in Grudić et al. (2021). For example,
the systematic differences we predict between the IMFs with the “fast” versus
“slow” accretion models (§ 4.5.2), or at different SMBH-centric radii (§
4.5.3) or in simulations with or without magnetic fields § 4.5.6, are all much
larger than the difference between the IMF at resolutions from $\sim
10^{-5}-10^{-2}\,{\rm M}_{\odot}$ in Grudić et al. (2021). Furthermore, there
are good physical reasons to think the convergence may be better here than in
the tests in small LISM clouds presented in those papers: for example, most of
the resolution dependence of the IMF in those tests owed to the ability to
resolve more fragmentation of individual cores into multiple (lower-mass)
stars at higher resolution. But owing the extreme tidal fields and magnetic
fields, we argued that such fragmentation is and must be very strongly
suppressed here (see § 4.5.7 & 4.6). Thus our key conclusions should be
robust.

That said, we have still attempted to check for additional resolution
dependence. In Paper II, we describe a set of resolution tests for both our
default simulation and the “no magnetic fields” simulations, where we apply an
additional refinement layer applied interior to $<0.1\,$pc, with mass
resolution reaching as small as $\sim 10^{-7}\,{\rm M}_{\odot}$ at $\lesssim
80\,$au, shown in Fig. 3. By necessity owing to the increased computational
expense this could only be run for $\sim 10^{3}$ dynamical times at the
innermost radii, or $\sim 100\,$yr, and the increased resolution was focused
on the smallest radii where the enclosed gas mass was small (reaching such
resolution inside the entire BHROI would require $\sim 10^{14}$ gas cells, an
infeasible number). The primary goal of that study was (as discussed in Paper
II) not to study the IMF, but to confirm that the global QAD properties in the
inner regions were converged. But this still allows us to confirm a few
important conclusions of ours: (1) the global properties of the QAD at
$\lesssim 0.1\,$pc are robust to resolution; (2) the strong suppression of SF
at small radii is also robust – we see zero sinks form for the duration of
this higher-resolution run interior to $\lesssim 0.05\,$pc; (3) the weak
effects of stellar ejecta which reach this radii, and the dynamics of the
(older) massive stars on eccentric orbits going through the inner disk, are
essentially identical. However, because of the refinement scheme adopted, this
only applies a modest additional refinement in the outermost QAD where star
formation occurs (factor of $\sim 4$ in mass resolution at $\sim 0.05\,$pc),
and none in the CQM (see Fig. 3 for details), and the duration of the
additional hyper-refinement ($\sim 100\,$yr) is just $\sim 1$ dynamical time
in the outer QAD at $\sim 0.1\,$pc. So this is not so useful for studying the
effects of resolution on the IMF and multiplicity statistics. We do see a few
sinks ($6$, specifically) form in the outer QAD ($\sim 0.05-0.1\,$pc) after
the additional refinement layer is applied, and can confirm (1) they form
exclusively at the toroidal magnetic field switch locations as we saw above (§
4.1); (2) include zero binaries, consistent with the low binary fraction seen
above (§ 4.5.7); and (3) are statistically consistent with being drawn from
the young ($<100\,$yr-old) IMF of stars formed at the same radii in our
“default” run (§ 4.5.3) – but we stress that the latter two statements in
particular are not highly statistically significant, given the tiny sample
size.

Unfortunately, refining even one-level higher in the CQM is extremely
computationally expensive, as the total gas mass enclosed in the BHROI ($\sim
10^{7}\,{\rm M_{\odot}}$) means this requires $\sim 10^{10}$ gas cells. We
attempted such a refinement but given the extreme cost could only run this
test for a small fraction of one dynamical time at these larger radii, so the
system could not relax and we do not consider this a particularly useful test.
So we suspect the sub-Solar IMF at these radii $\gtrsim 1\,$pc, where tidal
effects from the BH may not suppress fragmentation as severely and the IMF
begins to approach a more LISM-like form, to perhaps be the most uncertain of
the predictions we present here. In future work, dedicated hyper-refinement
centered on a smaller total gas mass (as opposed to refining on everything
within some radius of the SMBH) could be used to study better the behavior of
the IMF at these intermediate radii, as well as the IMF in the ISM of the
starburst galaxy (at $\sim$ kpc scales, where our current resolution and FIRE
treatment prevent any meaningful IMF predictions).

##  6 Conclusions

We have studied the stellar IMF in the outer regions of quasar accretion disks
and the circum-quasar medium (QADs and CQM), where the SMBH dominates the
potential and star formation is suppressed but still can occur at some level.
For the first time, we consider this with a combination of detailed physics,
which have been shown (in previous STARFORGE papers) to reproduce a plausible
observed IMF in LISM/“typical” Galactic GMC conditions, including: non-
equilibrium ion/atomic/molecular/dust thermo-chemistry coupled to explicit
multi-band radiation transport capable of handling both optically thick and
thin regimes, magnetic fields, self-gravity, and star formation and accretion
with detailed proto-stellar and main-sequence stellar evolutionary models
enabling feedback in the form of (proto)stellar jets, winds/mass-loss,
radiation, and supernovae. We also for the first time apply this to
simulations which self-consistently form a QAD around a SMBH from cosmological
initial conditions, so we do not have to make an ad-hoc assumption for
critical properties such as boundary/initial conditions, magnetic field
strengths/geometries, etc. Our important conclusions include:

  1. 1.

Star formation is possible in the CQM & outer QAD, though it is strongly
suppressed on sub-pc scales around the quasar by magnetic fields and optical
depth effects (Paper I), and almost completely vanishes on $<0.01\,$pc scales.
With the physics here the SFR is much less than the inflow rates to the SMBH
so stars and SF do not appreciably influence the structure or global dynamics
of the QAD or quasar accretion. It is possible that much later, very low-level
accretion could be driven by stellar mass-loss after the QAD is completely
consumed or blown away, but this would necessarily involve orders-of-magnitude
lower accretion rates.

  2. 2.

The sites of star formation are “special.” There are two “zones” of star
formation with different detailed structure. In the CQM inside the BHROI but
outside the circularized QAD, star formation occurs in fragmenting filamentary
inflows enhanced by strong super-sonic transverse tidal compression/shocks,
producing a population of stars on highly radial/eccentric orbits. Inside
where the QAD circularizes (the outer QAD) star formation occurs via gravito-
turbulence in special regions where toroidal magnetic field lines stretched
from the CQM reverse sign, hence where magnetic pressure is locally unable to
fully-stabilize the QAD against local collapse.

  3. 3.

The IMF is likely top-heavy at “intermediate” radii $\sim 0.1\,$pc between the
innermost CQM and outermost QAD. At larger radii approaching the BHROI, the
system IMF approaches the “typical” or “universal” IMF form (despite the fact
that the physical conditions at the BHROI are still orders-of-magnitude
different from typical Solar-neighborhood GMCs). At very small radii $\ll
0.01\,$pc essentially no star formation occurs. The variation of IMF versus
radius within the QAD (but outside this radius) depends on the unresolved sub-
grid assumptions regarding efficiency of accretion from PSDs onto stars in
these extreme conditions, but the IMF can at some point become more bottom-
heavy at small radii simply because the total enclosed QAD gas mass (and
corresponding Hill mass) at some $R$ become very small ($\ll 10\,M_{\odot}$)
and so massive star formation is not possible.

  4. 4.

How “top-heavy” the IMF actually is, and perhaps whether close binaries and/or
planets can form, depends sensitively on unresolved accretion from sink radii
to (proto)stars. Even in our extremely high-resolution simulations, we cannot
resolve the interior structure of accretion flows/PSDs around newly-forming
protostars interior to the numerical sink radii. Typical assumptions motivated
by PSDs in the Solar-neighborhood ISM naively predict timescales for accretion
over this range of scales which are many thousands of global orbital times in
the QAD – in other words, essentially no accretion onto the (proto)star. But
this is almost certainly unphysical, as the accretion “reservoir,” even if it
could form a PSD (which seems unlikely) must be far more compact, dense, warm,
and strongly magnetized if they form in the QAD owing to the extreme
environment. We show that assuming the accretion is “fast” (occurring on a few
sink dynamical times) leads to a far less extremely top-heavy system IMF. This
clearly motivates dedicated simulations of individual core and (proto)star
formation which can resolve scales $\ll$ au, with initial/boundary conditions
motivated by the simulations here.

  5. 5.

Accretion beyond the initial core is negligible. Beyond the initial tightly-
bound core collapse, stars accrete a negligible amount of gas from the
CQM/QAD, even integrating their trajectories over $\sim 10^{5}$ QAD dynamical
times during which the global QAD properties (e.g. the hypothetical mass
supply available to the star) are constant. If anything, the IMF becomes
slightly less top-heavy with time. All the initial core accretion occurs on
$\lesssim 1$ initial clump dynamical time. The accretion is strongly
suppressed, even in the absence of any (proto)stellar feedback, by the strong
magnetic fields, turbulence, radiation, warm gas pressures, and non-trivial
orbit structure (e.g. eccentric and spiral modes in the QAD as well as
stellar-gas orbit separation as detailed below). This is qualitatively
different from simulations with the same code and physics in Solar-
neighborhood GMCs, where the most massive stars become so by sustaining
accretion over $\gtrsim$ Myr timescales (Grudić et al., 2022; Guszejnov et
al., 2022a).

  6. 6.

Stars almost immediately dynamically decouple from gas after formation, and
follow orbits which are fundamentally distinct from the surrounding gas, with
typical systematic relative velocities between the stars and the immediately
surrounding ambient gas between $\sim 0.1-1$ times the Keplerian speed, i.e.
$\gtrsim 100-1000\,{\rm km\,s^{-1}}$ at the radii of interest. This is
expected and indeed inevitable, for fundamental dynamical reasons in any disk
which has non-negligible eccentricity, warps, precession,
pressure/magnetic/radiation/turbulent support, migration or sinking of clumps
and/or stars, etc. As a consequence, stars and gas are not co-orbital, which
has important effects on their post-initial collapse evolution, helping to
suppress accretion by orders of magnitude, and leaving ejecta/jet material
behind on cometary trails rather than allowing it to coherently locally
regulate the structure of the gas.

  7. 7.

Multiplicity on $\gg 10\,$au scales is strongly suppressed, owing to a
combination of inefficient fragmentation (the same physics that pushes the IMF
to be top-heavy) and the strong tidal environment disrupting any binaries that
are not very compact/hard. Whether compact binaries form is un-resolved here
and probably depends on the same uncertain physics of the un-resolved PSD
accretion/structure described above.

  8. 8.

Magnetic fields play a crucial role. The magnetic fields in the CQM and QAD
are extremely strong (plasma $\beta\sim 10^{-4}-10^{-2}$, field strengths
$\gg$ Gauss) and the ionization fractions extremely high ($\gtrsim 0.01$,
given the warm gas temperatures and extreme radiation environments) – hence
the magnetic fields much more strongly-coupled – compared to typical Galactic
GMCs and PSDs. We showed in Paper I and Paper II that the properties of the
CQM and QAD are fundamentally different without magnetic fields, including
(notably) that the system experiences runaway catastrophic fragmentation on
these scales if we (unphysically) remove said fields. So obviously they shape
the initial conditions in a key manner at all radii we consider (not just in
the QAD). We show here that the ensuing IMF absent magnetic fields is vastly
more top-heavy than the (already top-heavy) IMF with magnetic fields, and
becomes even more top-heavy with increasing distance from the SMBH, so
magnetic fields are also playing a key role regulating the masses of these
systems.

  9. 9.

Many proposed explanations/scalings for the IMF mass are ruled out at least as
an explanation for the simulations here. For example, simply assuming the IMF
turnover mass should scale with e.g. a thermal Jeans/Toomre/magnetic
critical/Hill/Roche/tidal/first Larson core/opacity limit mass, gives the
qualitatively incorrect scalings compared to the behavior we see. Likewise
explanations predicated on competitive accretion are not viable given the lack
of accretion noted above. We argue that the IMF shape in simulations with slow
PSD accretion (hence weak stellar feedback) are potentially explained by
analytic models for gravito-turbulent fragmentation in a disk with trans-
Alfvénic turbulence and efficient cooling from Hopkins & Christiansen (2013);
while in simulations with fast PSD accretion (hence more efficient stellar
feedback) are better explained by a feedback-regulated model as in Guszejnov
et al. (2022b).

There are many obvious interesting questions to explore in future work,
including more detailed exploration of properties such as multiplicity and
stellar correlation functions, and pushing further out in distance from the
SMBH at sufficiently high resolution to model the IMF in galactic nuclear
starburst environments. It would also of course be interesting to explore
similar simulations zooming in on different times and in different galaxies
accretion phases, to explore different conditions, such as star formation
around SMBHs more akin to Sgr A∗ (at much later cosmic times with vastly lower
gas supply to the nucleus even during the presumptive star-forming episode).
Indeed, although it is intriguing to note that the predicted IMFs here are top
heavy at a similar spatial scale around the SMBH to where the Galactic nucleus
is also observed to have a top-heavy IMF (and the IMF shapes are potentially
similar to those observational claims), we strongly caution that we are
simulating a vastly different environment from the Galactic nucleus. The
simulations here or others like them could also be used to explore predictions
and consequences for LIGO sources and TDEs (evolving the stars and their
relics dynamically in more idealized simulations or analytic calculations,
using their formation properties here as initial conditions). And of course,
as we have noted above, it would be especially interesting, and important for
informing future generations of circum-SMBH star formation simulations, to
consider much higher-resolution/smaller-scale simulations of individual
cores/accretion flows/PSDs and formation of individual massive stars with
initial/boundary conditions taken from the simulations here, to inform the
crucial uncertain sub-grid physics herein.

###### Acknowledgements.

Support for PFH was provided by NSF Research Grants 1911233, 20009234,
2108318, NSF CAREER grant 1455342, NASA grants 80NSSC18K0562, HST-AR-15800.
Numerical calculations were run on the Caltech compute cluster “Wheeler,”
allocations AST21010 and AST20016 supported by the NSF and TACC, and NASA HEC
SMD-16-7592. Support for MYG was provided by NASA through the NASA Hubble
Fellowship grant #HST-HF2-51479 awarded by the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Inc., for NASA, under contract NAS5-26555.

## References

  * Alig et al. (2013) Alig C., Schartmann M., Burkert A., Dolag K., 2013, ApJ, 771, 119
  * Andrews et al. (2018) Andrews S. M., et al., 2018, ApJ, 869, L41
  * Anglés-Alcázar et al. (2021) Anglés-Alcázar D., et al., 2021, ApJ, 917, 53
  * Bai (2017) Bai X.-N., 2017, ApJ, 845, 75
  * Bai & Stone (2017) Bai X.-N., Stone J. M., 2017, ApJ, 836, 46
  * Barnes & Hernquist (1996) Barnes J. E., Hernquist L., 1996, ApJ, 471, 115
  * Bartko et al. (2010) Bartko H., et al., 2010, ApJ, 708, 834
  * Baruteau et al. (2011) Baruteau C., Cuadra J., Lin D. N. C., 2011, ApJ, 726, 28
  * Bastian et al. (2010) Bastian N., Covey K. R., Meyer M. R., 2010, ARA&A, 48, 339
  * Bate (2009) Bate M. R., 2009, MNRAS, 392, 1363
  * Bate (2012) Bate M. R., 2012, MNRAS, 419, 3115
  * Bate & Bonnell (2005) Bate M. R., Bonnell I. A., 2005, MNRAS, 356, 1201
  * Beattie et al. (2021) Beattie J. R., Mocz P., Federrath C., Klessen R. S., 2021, arXiv e-prints, p. arXiv:2109.10470
  * Berentzen et al. (2007) Berentzen I., Shlosman I., Martinez-Valpuesta I., Heller C. H., 2007, ApJ, 666, 189
  * Bolatto et al. (2008) Bolatto A. D., Leroy A. K., Rosolowsky E., Walter F., Blitz L., 2008, ApJ, 686, 948
  * Bonnell & Rice (2008) Bonnell I. A., Rice W. K. M., 2008, Science, 321, 1060
  * Bortolas (2022) Bortolas E., 2022, MNRAS, 511, 2885
  * Cappellari et al. (2012) Cappellari M., et al., 2012, Nature, 484, 485
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Chabrier (2005) Chabrier G., 2005, in E. Corbelli, F. Palla, & H. Zinnecker ed., Astrophysics and Space Science Library Vol. 327, The Initial Mass Function 50 Years Later (Published by Springer, Dordrecht). pp 41–+ (arXiv:astro-ph/0409465)
  * Chen & Ostriker (2014) Chen C.-Y., Ostriker E. C., 2014, ApJ, 785, 69
  * Chiang & Youdin (2010) Chiang E., Youdin A. N., 2010, Annual Review of Earth and Planetary Sciences, 38, 493
  * Choban et al. (2022) Choban C. R., Kereš D., Hopkins P. F., Sandstrom K. M., Hayward C. C., Faucher-Giguère C.-A., 2022, MNRAS, 514, 4506
  * Colín et al. (2013) Colín P., Vázquez-Semadeni E., Gómez G. C., 2013, MNRAS, 435, 1701
  * Colombo et al. (2014) Colombo D., et al., 2014, ApJ, 784, 3
  * Cuadra et al. (2009) Cuadra J., Armitage P. J., Alexander R. D., Begelman M. C., 2009, MNRAS, 393, 1423
  * Davies & Lin (2020) Davies M. B., Lin D. N. C., 2020, MNRAS, 498, 3452
  * Davis & Tchekhovskoy (2020) Davis S. W., Tchekhovskoy A., 2020, ARA&A, 58, 407
  * Deng et al. (2017) Deng H., Mayer L., Meru F., 2017, ApJ, 847, 43
  * Deng et al. (2020) Deng H., Mayer L., Latter H., 2020, ApJ, 891, 154
  * Espinoza et al. (2009) Espinoza P., Selman F. J., Melnick J., 2009, A&A, 501, 563
  * Fall et al. (2010) Fall S. M., Krumholz M. R., Matzner C. D., 2010, ApJ, 710, L142
  * Fan & Wu (2023) Fan X., Wu Q., 2023, ApJ, 944, 159
  * Farias et al. (2023) Farias J. P., Offner S. S. R., Grudić M. Y., Guszejnov D., Rosen A. L., 2023, arXiv e-prints, p. arXiv:2309.11415
  * Flaherty et al. (2015) Flaherty K. M., Hughes A. M., Rosenfeld K. A., Andrews S. M., Chiang E., Simon J. B., Kerzner S., Wilner D. J., 2015, ApJ, 813, 99
  * Flaherty et al. (2017) Flaherty K. M., et al., 2017, ApJ, 843, 150
  * Forgan et al. (2017) Forgan D., Price D. J., Bonnell I., 2017, MNRAS, 466, 3406
  * Foyle et al. (2010) Foyle K., Rix H.-W., Walter F., Leroy A. K., 2010, ApJ, 725, 534
  * Gammie (2001) Gammie C. F., 2001, ApJ, 553, 174
  * Gandhi et al. (2022) Gandhi P. J., Wetzel A., Hopkins P. F., Shappee B. J., Wheeler C., Faucher-Giguère C.-A., 2022, MNRAS, 516, 1941
  * Gavagnin et al. (2017) Gavagnin E., Bleuler A., Rosdahl J., Teyssier R., 2017, MNRAS, 472, 4155
  * Goodman & Arce (2004) Goodman A. A., Arce H. G., 2004, ApJ, 608, 831
  * Goodman et al. (2009) Goodman A. A., Pineda J. E., Schnee S. L., 2009, ApJ, 692, 91
  * Grasha et al. (2019) Grasha K., et al., 2019, MNRAS, 483, 4707
  * Grudić (2021) Grudić M. Y., 2021, MNRAS, 507, 1064
  * Grudić & Hopkins (2020) Grudić M. Y., Hopkins P. F., 2020, MNRAS, 495, 4306
  * Grudić & Hopkins (2023) Grudić M. Y., Hopkins P. F., 2023, arXiv e-prints, p. arXiv:2308.16268
  * Grudić et al. (2018a) Grudić M. Y., Hopkins P. F., Faucher-Giguère C.-A., Quataert E., Murray N., Kereš D., 2018a, MNRAS, 475, 3511
  * Grudić et al. (2018b) Grudić M. Y., Guszejnov D., Hopkins P. F., Lamberts A., Boylan-Kolchin M., Murray N., Schmitz D., 2018b, MNRAS, 481, 688
  * Grudić et al. (2019) Grudić M. Y., Hopkins P. F., Quataert E., Murray N., 2019, MNRAS, 483, 5548
  * Grudić et al. (2020) Grudić M. Y., Boylan-Kolchin M., Faucher-Giguère C.-A., Hopkins P. F., 2020, MNRAS, 496, L127
  * Grudić et al. (2021) Grudić M. Y., Guszejnov D., Hopkins P. F., Offner S. S. R., Faucher-Giguère C.-A., 2021, MNRAS, 506, 2199
  * Grudić et al. (2022) Grudić M. Y., Guszejnov D., Offner S. S. R., Rosen A. L., Raju A. N., Faucher-Giguère C.-A., Hopkins P. F., 2022, MNRAS, 512, 216
  * Guszejnov & Hopkins (2015) Guszejnov D., Hopkins P. F., 2015, MNRAS, 450, 4137
  * Guszejnov et al. (2016) Guszejnov D., Krumholz M. R., Hopkins P. F., 2016, MNRAS, 458, 673
  * Guszejnov et al. (2017a) Guszejnov D., Hopkins P. F., Krumholz M. R., 2017a, MNRAS, 468, 4093
  * Guszejnov et al. (2017b) Guszejnov D., Hopkins P. F., Ma X., 2017b, MNRAS, 472, 2107
  * Guszejnov et al. (2018a) Guszejnov D., Hopkins P. F., Grudić M. Y., 2018a, MNRAS, 477, 5139
  * Guszejnov et al. (2018b) Guszejnov D., Hopkins P. F., Grudić M. Y., Krumholz M. R., Federrath C., 2018b, MNRAS, 480, 182
  * Guszejnov et al. (2019) Guszejnov D., Hopkins P. F., Graus A. S., 2019, MNRAS, 485, 4852
  * Guszejnov et al. (2020) Guszejnov D., Grudić M. Y., Hopkins P. F., Offner S. S. R., Faucher-Giguère C.-A., 2020, MNRAS, 496, 5072
  * Guszejnov et al. (2021) Guszejnov D., Grudić M. Y., Hopkins P. F., Offner S. S. R., Faucher-Giguère C.-A., 2021, MNRAS, 502, 3646
  * Guszejnov et al. (2022a) Guszejnov D., Markey C., Offner S. S. R., Grudić M. Y., Faucher-Giguère C.-A., Rosen A. L., Hopkins P. F., 2022a, MNRAS, 515, 167
  * Guszejnov et al. (2022b) Guszejnov D., Grudić M. Y., Offner S. S. R., Faucher-Giguère C.-A., Hopkins P. F., Rosen A. L., 2022b, MNRAS, 515, 4929
  * Guszejnov et al. (2023) Guszejnov D., Raju A. N., Offner S. S. R., Grudić M. Y., Faucher-Giguère C.-A., Hopkins P. F., Rosen A. L., 2023, MNRAS, 518, 4693
  * Hennebelle & Chabrier (2009) Hennebelle P., Chabrier G., 2009, ApJ, 702, 1428
  * Hennebelle et al. (2020) Hennebelle P., Commerçon B., Lee Y.-N., Chabrier G., 2020, ApJ, 904, 194
  * Heyer et al. (2009) Heyer M., Krawczyk C., Duval J., Jackson J. M., 2009, ApJ, 699, 1092
  * Hobbs & Nayakshin (2009) Hobbs A., Nayakshin S., 2009, MNRAS, 394, 191
  * Hocuk & Spaans (2011) Hocuk S., Spaans M., 2011, A&A, 536, A41
  * Hocuk et al. (2012) Hocuk S., Schleicher D. R. G., Spaans M., Cazaux S., 2012, A&A, 545, A46
  * Hopkins (2010) Hopkins P. F., 2010, arXiv e-prints, arXiv:1009.4702 [astro-ph],
  * Hopkins (2012a) Hopkins P. F., 2012a, MNRAS, 423, 2016
  * Hopkins (2012b) Hopkins P. F., 2012b, MNRAS, 423, 2037
  * Hopkins (2013a) Hopkins P. F., 2013a, MNRAS, 428, 1950
  * Hopkins (2013b) Hopkins P. F., 2013b, MNRAS, 430, 1653
  * Hopkins (2013c) Hopkins P. F., 2013c, MNRAS, 433, 170
  * Hopkins (2016) Hopkins P. F., 2016, MNRAS, 462, 576
  * Hopkins (2017) Hopkins P. F., 2017, MNRAS, 466, 3387
  * Hopkins (2018) Hopkins A. M., 2018, PASA, 35, e039
  * Hopkins (2023) Hopkins P. F., 2023, MNRAS, 518, 5882
  * Hopkins & Christiansen (2013) Hopkins P. F., Christiansen J. L., 2013, ApJ, 776, 48
  * Hopkins & Grudić (2019) Hopkins P. F., Grudić M. Y., 2019, MNRAS, 483, 4187
  * Hopkins & Quataert (2010a) Hopkins P. F., Quataert E., 2010a, MNRAS, 405, L41
  * Hopkins & Quataert (2010b) Hopkins P. F., Quataert E., 2010b, MNRAS, 407, 1529
  * Hopkins & Quataert (2011a) Hopkins P. F., Quataert E., 2011a, MNRAS, 411, L61
  * Hopkins & Quataert (2011b) Hopkins P. F., Quataert E., 2011b, MNRAS, 415, 1027
  * Hopkins & Raives (2016) Hopkins P. F., Raives M. J., 2016, MNRAS, 455, 51
  * Hopkins et al. (2016) Hopkins P. F., Torrey P., Faucher-Giguère C.-A., Quataert E., Murray N., 2016, MNRAS, 458, 816
  * Hopkins et al. (2018) Hopkins P. F., et al., 2018, MNRAS, 480, 800
  * Hopkins et al. (2020) Hopkins P. F., Grudić M. Y., Wetzel A., Kereš D., Faucher-Giguère C.-A., Ma X., Murray N., Butcher N., 2020, MNRAS, 491, 3702
  * Hopkins et al. (2022a) Hopkins P. F., Squire J., Butsky I. S., 2022a, MNRAS, 509, 3779
  * Hopkins et al. (2022b) Hopkins P. F., Wellons S., Anglés-Alcázar D., Faucher-Giguère C.-A., Grudić M. Y., 2022b, MNRAS, 510, 630
  * Hopkins et al. (2022c) Hopkins P. F., Butsky I. S., Panopoulou G. V., Ji S., Quataert E., Faucher-Giguère C.-A., Kereš D., 2022c, MNRAS, 516, 3470
  * Hopkins et al. (2022d) Hopkins P. F., Squire J., Butsky I. S., Ji S., 2022d, MNRAS, 517, 5413
  * Hopkins et al. (2023a) Hopkins P. F., et al., 2023a, arXiv e-prints, p. arXiv:2309.13115
  * Hopkins et al. (2023b) Hopkins P. F., et al., 2023b, arXiv e-prints, p. arXiv:2310.04506
  * Hopkins et al. (2023c) Hopkins P. F., et al., 2023c, arXiv e-prints, p. arXiv:2310.04507
  * Hopkins et al. (2023d) Hopkins P. F., et al., 2023d, MNRAS, 519, 3154
  * Hopkins et al. (2023e) Hopkins P. F., Butsky I. S., Ji S., Kereš D., 2023e, MNRAS, 522, 2936
  * Hopkins et al. (2023f) Hopkins P. F., Nadler E. O., Grudić M. Y., Shen X., Sands I., Jiang F., 2023f, MNRAS, 525, 5951
  * Hosek et al. (2019) Hosek Matthew W. J., Lu J. R., Anderson J., Najarro F., Ghez A. M., Morris M. R., Clarkson W. I., Albers S. M., 2019, ApJ, 870, 44
  * Hughes et al. (2011) Hughes A. M., Wilner D. J., Andrews S. M., Qi C., Hogerheijde M. R., 2011, ApJ, 727, 85
  * Jacobs & Sellwood (2001) Jacobs V., Sellwood J. A., 2001, ApJ, 555, L25
  * Jappsen et al. (2005) Jappsen A.-K., Klessen R. S., Larson R. B., Li Y., Mac Low M.-M., 2005, A&A, 435, 611
  * Kawakatu et al. (2020) Kawakatu N., Wada K., Ichikawa K., 2020, ApJ, 889, 84
  * Keith & Wardle (2014) Keith S. L., Wardle M., 2014, MNRAS, 440, 89
  * Kim et al. (2018) Kim J.-G., Kim W.-T., Ostriker E. C., 2018, ApJ, 859, 68
  * Kowal et al. (2007) Kowal G., Lazarian A., Beresnyak A., 2007, ApJ, 658, 423
  * Kratter & Lodato (2016) Kratter K. M., Lodato G., 2016, ARA&A, in press, arXiv:1603.01280,
  * Kritsuk et al. (2011) Kritsuk A. G., et al., 2011, ApJ, 737, 13
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Kroupa et al. (2011) Kroupa P., Weidner C., Pflamm-Altenburg J., Thies I., Dabringhausen J., Marks M., Maschberger T., 2011, Stellar Systems and Galactic Structure, V
  * Krumholz (2011) Krumholz M. R., 2011, ApJ, 743, 110
  * Krumholz & Federrath (2019) Krumholz M. R., Federrath C., 2019, Frontiers in Astronomy and Space Sciences, 6, 7
  * Kudoh et al. (2020) Kudoh Y., Wada K., Norman C., 2020, ApJ, 904, 9
  * Kuruwita & Haugbølle (2023) Kuruwita R. L., Haugbølle T., 2023, A&A, 674, A196
  * Larson (1981) Larson R. B., 1981, MNRAS, 194, 809
  * Lebreuilly et al. (2023) Lebreuilly U., Hennebelle P., Maury A., González M., Traficante A., Klessen R., Testi L., Molinari S., 2023, arXiv e-prints, p. arXiv:2309.05397
  * Lee & Hennebelle (2018) Lee Y.-N., Hennebelle P., 2018, A&A, 611, A89
  * Lee & Hennebelle (2019) Lee Y.-N., Hennebelle P., 2019, A&A, 622, A125
  * Lee et al. (2014) Lee A. T., Cunningham A. J., McKee C. F., Klein R. I., 2014, ApJ, 783, 50
  * Lee et al. (2019) Lee A. T., Offner S. S. R., Kratter K. M., Smullen R. A., Li P. S., 2019, ApJ, 887, 232
  * Lee et al. (2021) Lee Y.-N., Marchand P., Liu Y.-H., Hennebelle P., 2021, ApJ, 922, 36
  * Lee et al. (2022) Lee J. C., et al., 2022, arXiv e-prints, p. arXiv:2212.02667
  * Lemaster & Stone (2009) Lemaster M. N., Stone J. M., 2009, in Revista Mexicana de Astronomia y Astrofisica Conference Series. pp 243–+
  * Leroy et al. (2013) Leroy A. K., et al., 2013, AJ, 146, 19
  * Li et al. (2010) Li Z.-Y., Wang P., Abel T., Nakamura F., 2010, ApJ, 720, L26
  * Lin (2014) Lin M.-K., 2014, ApJ, 790, 13
  * Lizano et al. (2010) Lizano S., Galli D., Cai M. J., Adams F. C., 2010, ApJ, 724, 1561
  * Löckmann et al. (2010) Löckmann U., Baumgardt H., Kroupa P., 2010, MNRAS, 402, 519
  * Lu et al. (2013) Lu J. R., Do T., Ghez A. M., Morris M. R., Yelda S., Matthews K., 2013, ApJ, 764, 155
  * Ma et al. (2017) Ma X., Hopkins P. F., Wetzel A. R., Kirby E. N., Anglés-Alcázar D., Faucher-Giguère C.-A., Kereš D., Quataert E., 2017, MNRAS, 467, 2430
  * Mac Low & Klessen (2004) Mac Low M.-M., Klessen R. S., 2004, Reviews of Modern Physics, 76, 125
  * Maness et al. (2007) Maness H., et al., 2007, ApJ, 669, 1024
  * Mapelli et al. (2012) Mapelli M., Hayfield T., Mayer L., Wadsley J., 2012, ApJ, 749, 168
  * Martín-Navarro et al. (2015) Martín-Navarro I., et al., 2015, ApJ, 798, L4
  * McKee & Ostriker (2007) McKee C. F., Ostriker E. C., 2007, ARA&A, 45, 565
  * Meidt et al. (2022) Meidt S. E., et al., 2022, arXiv e-prints, p. arXiv:2212.06434
  * Meru & Bate (2012) Meru F., Bate M. R., 2012, MNRAS, in press, arXiv:1209.1107,
  * Moe & Di Stefano (2017) Moe M., Di Stefano R., 2017, ApJS, 230, 15
  * Molina et al. (2012) Molina F. Z., Glover S. C. O., Federrath C., Klessen R. S., 2012, MNRAS, 423, 2680
  * Myers et al. (2014) Myers A. T., Klein R. I., Krumholz M. R., McKee C. F., 2014, MNRAS, 439, 3420
  * Narayanan & Davé (2013) Narayanan D., Davé R., 2013, MNRAS, 436, 2892
  * Nayakshin & Sunyaev (2005) Nayakshin S., Sunyaev R., 2005, MNRAS, 364, L23
  * Nayakshin et al. (2007) Nayakshin S., Cuadra J., Springel V., 2007, MNRAS, 379, 21
  * Noguchi (1988) Noguchi M., 1988, A&A, 203, 259
  * Offner et al. (2014) Offner S. S. R., Clark P. C., Hennebelle P., Bastian N., Bate M. R., Hopkins P. F., Moraux E., Whitworth A. P., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. pp 53–75 (arXiv:1312.5326), doi:10.2458/azu_uapress_9780816531240-ch003
  * Offner et al. (2016) Offner S. S. R., Dunham M. M., Lee K. I., Arce H. G., Fielding D. B., 2016, ApJ, 827, L11
  * Offner et al. (2023) Offner S. S. R., Moe M., Kratter K. M., Sadavoy S. I., Jensen E. L. N., Tobin J. J., 2023, in Inutsuka S., Aikawa Y., Muto T., Tomida K., Tamura M., eds, Astronomical Society of the Pacific Conference Series Vol. 534, Protostars and Planets VII. p. 275 (arXiv:2203.10066), doi:10.48550/arXiv.2203.10066
  * Paardekooper (2012) Paardekooper S.-J., 2012, MNRAS, 421, 3286
  * Padoan & Nordlund (2002) Padoan P., Nordlund Å., 2002, ApJ, 576, 870
  * Paumard et al. (2006) Paumard T., et al., 2006, ApJ, 643, 1011
  * Pineda et al. (2023) Pineda J. E., et al., 2023, in Inutsuka S., Aikawa Y., Muto T., Tomida K., Tamura M., eds, Astronomical Society of the Pacific Conference Series Vol. 534, Protostars and Planets VII. p. 233 (arXiv:2205.03935), doi:10.48550/arXiv.2205.03935
  * Pokhrel et al. (2020) Pokhrel R., et al., 2020, ApJ, 896, 60
  * Raghavan et al. (2010) Raghavan D., et al., 2010, ApJS, 190, 1
  * Rees (1976) Rees M. J., 1976, MNRAS, 176, 483
  * Riols & Latter (2016) Riols A., Latter H., 2016, MNRAS, 460, 2223
  * Riols & Latter (2018) Riols A., Latter H., 2018, MNRAS, 474, 2212
  * Rosolowsky et al. (2003) Rosolowsky E., Engargiola G., Plambeck R., Blitz L., 2003, ApJ, 599, 258
  * Sadanari et al. (2022) Sadanari K. E., Omukai K., Sugimura K., Matsumoto T., Tomida K., 2022, arXiv e-prints, p. arXiv:2212.08075
  * Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
  * Sambhus & Sridhar (2002) Sambhus N., Sridhar S., 2002, A&A, 388, 766
  * Sellwood & Carlberg (1984) Sellwood J. A., Carlberg R. G., 1984, ApJ, 282, 61
  * Sellwood & Merritt (1994) Sellwood J. A., Merritt D., 1994, ApJ, 425, 530
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Sharda et al. (2020) Sharda P., Federrath C., Krumholz M. R., 2020, MNRAS, 497, 336
  * Shi et al. (2022) Shi Y., Kremer K., Grudić M. Y., Gerling-Dunsmore H. J., Hopkins P. F., 2022, MNRAS, submitted, arXiv e-prints, p. arXiv:2208.05025
  * Shu (1977) Shu F. H., 1977, ApJ, 214, 488
  * Soliman & Hopkins (2022) Soliman N. H., Hopkins P. F., 2022, arXiv e-prints, p. arXiv:2210.13517
  * Solomon et al. (1987) Solomon P. M., Rivolo A. R., Barrett J., Yahil A., 1987, ApJ, 319, 730
  * Spiniello et al. (2012) Spiniello C., Trager S. C., Koopmans L. V. E., Chen Y. P., 2012, ApJ, 753, L32
  * Squire & Hopkins (2018a) Squire J., Hopkins P. F., 2018a, MNRAS, 477, 5011
  * Squire & Hopkins (2018b) Squire J., Hopkins P. F., 2018b, ApJ, 856, L15
  * Stone et al. (2017) Stone N. C., Metzger B. D., Haiman Z., 2017, MNRAS, 464, 946
  * Su et al. (2017) Su K.-Y., Hopkins P. F., Hayward C. C., Faucher-Giguère C.-A., Kereš D., Ma X., Robles V. H., 2017, MNRAS, 471, 144
  * Tagawa et al. (2020) Tagawa H., Haiman Z., Kocsis B., 2020, ApJ, 898, 25
  * Tasker & Tan (2009) Tasker E. J., Tan J. C., 2009, ApJ, 700, 358
  * Thompson et al. (2005) Thompson T. A., Quataert E., Murray N., 2005, ApJ, 630, 167
  * Toyouchi et al. (2022) Toyouchi D., Inayoshi K., Ishigaki M. N., Tominaga N., 2022, MNRAS, 512, 2573
  * Tremaine (2001) Tremaine S., 2001, AJ, 121, 1776
  * Tsukamoto et al. (2015) Tsukamoto Y., Iwasaki K., Okuzumi S., Machida M. N., Inutsuka S., 2015, MNRAS, 452, 278
  * Tsukamoto et al. (2017) Tsukamoto Y., Okuzumi S., Iwasaki K., Machida M. N., Inutsuka S.-i., 2017, PASJ, 69, 95
  * Vázquez-Semadeni et al. (2011) Vázquez-Semadeni E., Banerjee R., Gómez G. C., Hennebelle P., Duffin D., Klessen R. S., 2011, MNRAS, 414, 2511
  * Wada (1994) Wada K., 1994, PASJ: Publications of the Astronomical Society of Japan (ISSN 0004-6264), 46, 165
  * Wada et al. (2009) Wada K., Papadopoulos P. P., Spaans M., 2009, ApJ, 702, 63
  * Wardle (2007) Wardle M., 2007, Astrophysics and Space Science, 311, 35
  * Wellons et al. (2023) Wellons S., et al., 2023, MNRAS, 520, 5394
  * Wurster et al. (2021) Wurster J., Bate M. R., Bonnell I. A., 2021, MNRAS, 507, 2354
  * Young & Scoville (1991) Young J. S., Scoville N. Z., 1991, ARA&A, 29, 581
  * Zhao et al. (2020) Zhao B., Caselli P., Li Z.-Y., Krasnopolsky R., Shang H., Lam K. H., 2020, MNRAS, 492, 3375
  * van Dokkum & Conroy (2011) van Dokkum P. G., Conroy C., 2011, ApJ, 735, L13+