# How Low Can Q Go?

John C. Forbes Center for Computational Astrophysics Flatiron Institute 162
5th Avenue New York, NY, 10010, USA John C. Forbes<EMAIL_ADDRESS>
###### Abstract

Gravitational instability plays a substantial role in the evolution of
galaxies. Various schemes to include it in galaxy evolution models exist,
generally assuming that the Toomre $Q$ parameter is self-regulated to
$Q_{\mathrm{crit}}$, the critical $Q$ dividing stable from unstable conditions
in a linear stability analysis. This assumption is in tension with
observational estimates of $Q$ that find values far below any plausible value
of $Q_{\mathrm{crit}}$. While the observations are subject to some
uncertainty, this tension can more easily be relieved on the theoretical side
by relaxing the common assumption that $Q\geq Q_{\mathrm{crit}}$. Based on
observations of both $z\sim 2$ disks and local face-on galaxies, we estimate
the effect of gravitational instability necessary to balance out every other
physical process that affects $Q$. In particular we find that the disk’s
response to low $Q$ values can be described by simple functions that depend
only on $Q$. These response functions allow galaxies to maintain $Q$ values
below $Q_{\mathrm{crit}}$ in equilibrium over a wide range of parameters.
Extremely low values of $Q$ are predicted when the gas surface density is
$\gtrsim 10^{3}M_{\odot}\ \mathrm{pc}^{-2}$, the rotation curve provides
minimal shear, the orbital time becomes long, and/or when the gas is much more
unstable than the stellar component. We suggest that these response functions
should be used in place of the $Q\geq Q_{\mathrm{crit}}$ ansatz.

††software: numpy (van der Walt et al., 2011; Harris et al., 2020), matplotlib
(Hunter, 2007), scipy (Virtanen et al., 2020)

##  1 Introduction

A variety of theoretical and observational results over the past two-plus
decades have suggested that gravitational instability, fundamentally like the
linear instability in axisymmetric self-gravitating disks discovered by Toomre
(1964), may play a substantial role in the evolution of galaxies. $Q$
summarizes the stability of such disks, where $Q$ below some critical value
$Q_{\mathrm{crit}}$ of order 1 indicates that the disk will be subject to
clumping or the formation of spiral arms. In local galaxies where surface
densities and kinematics of both gaseous and stellar components of galaxies
may be measured, the Toomre $Q$ parameter is observed to be in the
neighborhood of 2 across many different galaxies (Leroy et al., 2008). At
$z\sim 2$, gas velocity dispersions (e.g. Kassin et al., 2012; Genzel et al.,
2014; Wisnioski et al., 2015; Förster Schreiber et al., 2018; Übler et al.,
2019; Price et al., 2020) far in excess of the level likely to be sustainable
by supernova feedback alone (Joung et al., 2009; Krumholz et al., 2018; Brucy
et al., 2020), along with clumpy UV morphologies (Elmegreen et al., 2007) have
also suggested that gravitational instability had a substantial role to play.
Indeed Genzel et al. (2014) has reported estimates of the Toomre $Q$ well
below $Q_{\mathrm{crit}}$. There are also a number of local examples of $Q$
falling below $Q_{\mathrm{crit}}$ (Westfall et al., 2014) from the DiskMass
Survey (Bershady et al., 2010).

In modeling these disks it has been convenient to assume that
$Q=Q_{\mathrm{crit}}$ (e.g. Thompson et al., 2005; Dekel et al., 2009a;
Krumholz & Burkert, 2010; Forbes et al., 2012; Faucher-Giguère et al., 2013;
Krumholz & Burkhart, 2016; Hayward & Hopkins, 2017; Krumholz et al., 2018), or
at least that $Q$ may not fall below $Q_{\mathrm{crit}}$ (Forbes et al., 2014;
Rathaus & Sternberg, 2016; Forbes et al., 2019). Physically the idea is that
when $Q$ falls below $Q_{\mathrm{crit}}$, the resulting fragmentation,
formation of clumps, star formation, and/or torques exerted by the resulting
inhomogeneous gravitational potential, will tend to increase the velocity
dispersion and drive $Q$ back to $Q_{\mathrm{crit}}$. This assumption has
enabled, under various additional assumptions, the calculation of star
formation rates, mass loading factors, mass transport rates, viscous evolution
timescales, and contributions to the turbulent energy of the disk. This
process does occur in simulations, where many groups have found that galactic
disks tend to self-regulate to $Q\sim 1$ (Bournaud & Elmegreen, 2009; Goldbaum
et al., 2015, 2016; Behrendt et al., 2016; Gurvich et al., 2020) in an
azimuthally-averaged sense.

The exact value of $Q_{\mathrm{crit}}$ is complicated by several factors.
First, galactic disks have multiple components (gas and stars), and indeed
more if one requires that each component be describable by a single
temperature or velocity dispersion. Each quasi-isothermal component is
described by its own $Q$ parameter,

$Q_{i}=\frac{\kappa\sigma_{i}}{\pi G\Sigma_{i}}$ (1)

where $\kappa=\sqrt{2(\beta+1)}\Omega$ is the epicyclic frequency,
$\Omega=v_{\mathrm{circ}}/r$ is the orbital frequency of an object orbiting on
a circular orbit at velocity $v_{\mathrm{circ}}$, and $\beta=d\ln
v_{\mathrm{circ}}/d\ln r$ is the logarithmic derivative of
$v_{\mathrm{circ}}$. The velocity dispersion and column density of the $i$th
component are $\sigma_{i}$ and $\Sigma_{i}$ respectively. For a single-
component disk, equation 1 is exactly $Q$. The gravitational stability of a
gas-star mixture is approximately $Q^{-1}\sim Q_{g}^{-1}+Q_{*}^{-1}$ (Wang &
Silk, 1994). Intuitively the most unstable component dominates the stability
of the mixture. More accurate approximations (Romeo & Wiegert, 2011) account
for the finite thickness of the disk and deviations from the particularly
simple Wang & Silk (1994) approximation. With more components, the most
unstable mode (which sets $Q$) can be found numerically as well (Rafikov,
2001). Closer examination of the dispersion relation from which the $Q$
stability thresholds are derived has yielded the realization that the presence
of fast cooling in galactic disks means that $Q$ needs to be substantially
higher than the formal linear instability threshold, i.e $Q_{\mathrm{crit}}$
should be closer to 2 or 3 than 1, to realistically avoid Toomre-like
instabilities (Elmegreen, 2011). Perhaps because finite thickness and multi-
component corrections are often difficult to estimate directly from data,
especially outside the local Universe, and perhaps because even the gas-only
measurements are challenging with a wide variety of systematic uncertainties
(Davies et al., 2011), little concern has been devoted to the discrepancy
between $Q_{\mathrm{crit}}\approx 2-3$, and the observed
$Q_{\mathrm{gas}}\approx 0.2$ in $z\sim 2$ surveys.

In this work we relax the $Q\geq Q_{\mathrm{crit}}$ assumption and derive
equilibrium $Q$ values. Importantly these equilibrium values of $Q$ are not
about what $Q_{\mathrm{crit}}$ itself is, but rather how far below
$Q_{\mathrm{crit}}$ $Q$ may go. In the next section we show the key governing
equations and discuss the uncertainties in this formulation, then in section 3
we show how these uncertainties may be at least partially resolved by
estimating how quickly $Q$ is changing from everything besides gravitational
instability in the observed galaxies. The results are interpreted in section 4
and the implications are discussed in section 5. We conclude in section 6.

##  2 Governing Equations

Our goal is to find $Q_{\mathrm{eq}}$, the value of $Q$ at which

$\frac{dQ}{dt}=0,$ (2)

in the regime where gravitational instability is active, i.e. $Q\leq
Q_{\mathrm{crit}}$. Equation 2 may be re-expressed using the chain rule, and
by dividing up the resulting terms into those arising from gravitational
instability, which we denote “GI,” and all other terms, which we call
“sources.”

$\frac{dQ}{dt}=0=\frac{dQ}{dt}\Big{|}_{\mathrm{GI}}+\frac{dQ}{dt}\Big{|}_{\mathrm{sources}}=\frac{\partial
Q}{\partial\Sigma}\frac{\partial\Sigma}{\partial t}+\frac{\partial
Q}{\partial\sigma}\frac{\partial\sigma}{\partial t}+\frac{\partial
Q}{\partial\Sigma_{*}}\frac{\partial\Sigma_{*}}{\partial t}+\frac{\partial
Q}{\partial\sigma_{*}}\frac{\partial\sigma_{*}}{\partial t}$ (3)

In principle more terms can be added to the chain rule expansion, e.g.
derivatives with respect to other components of $Q$ like $v_{\mathrm{circ}}$,
but we assume that these terms are small (Krumholz & Burkert, 2010; Forbes et
al., 2012). The main issue in evaluating the terms in Equation 3 is deciding
which formulation of $Q$ itself to use. It is plausible that only
$Q_{\mathrm{gas}}$ matters, especially in the gas-rich $z\sim 2$ disks on
which we are focused, in which case the terms with $\Sigma_{*}$ and
$\sigma_{*}$ would be approximately zero. Unless otherwise noted, we will
adopt the Romeo & Wiegert (2011) approximation to $Q$ since it accounts for
the two components we consider here and for finite thickness effects. This
approximation has the disadvantage of involving a conditional depending on
which component of the disk contributes more to $Q$, so we will carry through
the partial derivatives with respect to $Q$ without explicitly evaluating them
analytically.

###  2.1 Source Terms

The source terms come from changes in the state of the disk that are unrelated
to the direct effects of gravitational instability. We enumerate these terms
as follows, recalling that the quantities in parentheses following each
$\partial Q/\partial X$ are just $\partial X/\partial t$ arising from any
physical process besides gravitational instability (see Eq. 3).

$\displaystyle\frac{dQ}{dt}\Big{|}_{\mathrm{sources}}=$
$\displaystyle\frac{\partial
Q}{\partial\Sigma}\left(\dot{\Sigma}_{\mathrm{accr}}-(f_{R}+\eta)\dot{\Sigma}_{\mathrm{SF}}\right)+$
(4) $\displaystyle\frac{\partial Q}{\partial\sigma}\left(\frac{\langle
p_{*}/m_{*}\rangle\dot{\Sigma}_{\mathrm{SF}}}{3\Sigma}-\frac{\sigma^{2}}{2H}\left(1-\frac{\sigma_{\mathrm{therm}}^{2}}{\sigma^{2}}\right)^{3/2}+\frac{\epsilon_{\mathrm{accr}}\dot{\Sigma}_{\mathrm{accr}}v_{\mathrm{circ}}^{2}}{3\Sigma\sigma}\right)+$
$\displaystyle\frac{\partial
Q}{\partial\Sigma_{*}}\left(f_{R}\dot{\Sigma}_{\mathrm{SF}}\right)+$
$\displaystyle\frac{\partial
Q}{\partial\sigma_{*}}\left(f_{R}\dot{\Sigma}_{\mathrm{SF}}\frac{\sigma^{2}-\sigma_{*}^{2}}{2\Sigma_{*}\sigma_{*}}\right)+$
$\displaystyle\frac{\partial Q}{\partial
Q_{*}}\left(\frac{Q_{*}-Q_{*,\mathrm{crit}}}{\gamma
t_{\mathrm{orb}}}\mathcal{I}_{Q_{*}}\right).$

The change in surface density per unit time at any point in the disk is the
sum of the effects of accretion, $\dot{\Sigma}_{\mathrm{accr}}$, and star
formation plus outflows, $(f_{R}+\eta)\dot{\Sigma}_{\mathrm{SF}}$. Accretion
is assumed to be set externally, e.g. by the physics of cosmological accretion
(e.g. Kereš et al., 2005; Dekel et al., 2009b), although in many galaxies it
may be more reasonable to tie it to past, even recent star formation via the
action of galactic fountains or star-formation induced condensation
(Oppenheimer et al., 2010; Hobbs et al., 2015). The effect of star formation
is reduced by a fraction $f_{R}\approx 0.5$ the fraction of material that
remains in long-lived stars or stellar remnants rather than being returned to
the ISM through stellar winds and supernovae (Tinsley, 1980; Leitner &
Kravtsov, 2011). Star formation is also assumed to remove gas from the galaxy
in galactic winds with mass loading factor $\eta$. In principle we could have
included a term proportional to $(1/2\pi r)(\partial\dot{M}/\partial r)$, the
change in surface density arising from a mass flow through the disk $\dot{M}$
when that flow rate varies with radius. We assume, however, that $\dot{M}$ is
primarily set by gravitational instability in this regime where $Q\leq
Q_{\mathrm{crit}}$, so this term is not included among the source terms
$dQ/dt|_{\mathrm{sources}}$.

The velocity dispersion of the gas, similarly, may be affected by
gravitational stability, but that effect is not included among the source
terms in Eq. 4. Instead, we include the effects from star formation feedback
(first term), turbulent dissipation (second term), and direct accretion (third
term). The star formation feedback term includes the factor $\langle
p_{*}/m_{*}\rangle$, the amount of momentum injected by supernovae per unit
mass formed in stars. Generally simulations find values of around $3000\
\mathrm{km}\ \mathrm{s}^{-1}$ (e.g. Kim & Ostriker, 2015), but for
sufficiently high gas fractions and long dynamical times, Orr et al. (2022)
find that this number can be reduced substantially by the breakout of
superbubbles from the disk before they finish their energy-driven phase. We
adopt the Orr et al. (2022) values of $\langle p_{*}/m_{*}\rangle$ unless
otherwise noted111In particular we evaluate their equation 8 for the velocity
of the superbubble at $t=t_{\mathrm{SNe}}$, the time that supernovae continue
to go off after a burst of star formation. If this velocity is greater than
the local velocity dispersion, we evaluate their Equation 19 for the momentum
per unit mass of stars formed. We do not adopt their simplifications like
assuming that $Q_{\mathrm{gas}}$ is 1.. The dissipation term assumes that
turbulent energy dissipates on a scaleheight crossing time (Stone et al.,
1998; Mac Low et al., 1998), and as in Forbes et al. (2014) this dissipation
is assumed to gradually approach zero as
$\sigma\rightarrow\sigma_{\mathrm{therm}}=8\ \mathrm{km}\ \mathrm{s}^{-1}$,
the sound speed of the warm neutral medium. We adopt the expression for $H$
from Ostriker & Kim (2022) (see below). For the contribution from direct
accretion, we adopt the formula suggested by Klessen & Hennebelle (2010) with
$\epsilon_{\mathrm{accr}}=10\%$ as measured in Illustris TNG50 by Forbes et
al. (2022).

The remaining three terms summarize the effect of changes to the stars, namely
changes to $\Sigma_{*}$ and $\sigma_{*}$, on $Q$. Just as star formation
removes gas, it adds to $\Sigma_{*}$. The stellar population as a whole also
loses specific kinetic energy by the addition of new stars that have the same
velocity dispersion as the (dynamically colder) gas. Stars may also heat up
via their interactions with spiral arms (Sellwood & Carlberg, 1984; Carlberg &
Sellwood, 1985) and other structures in the disk like molecular clouds (Lacey,
1984). We assume that these effects follow the Sellwood & Carlberg (1984) and
Carlberg & Sellwood (1985) prescription that $Q_{*}$ increases in proportion
to its difference from $Q_{*,\mathrm{crit}}$ per unit orbital time, so long as
$Q_{*}<Q_{*,\mathrm{crit}}$. This last requirement is encapsulated in the
indicator function $I_{Q_{*}}$, which is 1 when $Q_{*}<Q_{*,\mathrm{crit}}$
and zero otherwise. We adopt $\gamma=4$ and $Q_{*,\mathrm{crit}}=2$. Note that
this rate of change in $Q_{*}$ due to spiral heating could be decomposed into
its effects on $\Sigma_{*}$ and $\sigma_{*}$ (e.g. Forbes et al., 2014), but
for our purposes here all we care about is its effects on $Q$, so such a
decomposition is not necessary.

We also need to specify exactly how we estimate $H$, the star formation rate
$\dot{\Sigma}_{\mathrm{SF}}$, and the accretion rate
$\dot{\Sigma}_{\mathrm{accr}}$. Following Ostriker & Kim (2022), we adopt

$H=\frac{2\sigma^{2}}{\pi G\Sigma+((\pi
G\Sigma)^{2}+32\pi\zeta_{d}G\rho_{\mathrm{sd}}\sigma^{2})^{1/2}},$ (5)

where $\zeta_{d}\approx 0.33$ depending on the exact shape of the vertical
density profile, and $\rho_{\mathrm{sd}}$ is the density of stars plus dark
matter in the midplane. For convenience we adopt the simple assumptions
mentioned in Ostriker et al. (2010),

$\rho_{\mathrm{sd}}\approx\frac{v_{\mathrm{circ}}^{2}/r^{2}}{4\pi G}+\frac{\pi
G\Sigma_{*}^{2}}{2\sigma_{*}^{2}},$ (6)

where the first term comes from assuming that $v_{\mathrm{circ}}$ from dark
matter is constant with radius, and the second comes from taking the stars to
be locally isothermal with little modification in their vertical profile from
the self-gravity of the gas. Caution is warranted in adopting both of these in
high-z gas-rich disks, given that more explicit mass models find that the
disks are largely baryon-dominated (Genzel et al., 2017; Price et al., 2021),
and the gas fractions are high enough (Tacchella et al., 2015) that the gas
will likely have some effect on the stellar scaleheight. However, these same
effects that make Equation 6’s estimate of $\rho_{\mathrm{sd}}$ slightly worse
also make $\rho_{\mathrm{sd}}$ less relevant in calculating $H$ (see Equation
5).

Next we adopt the Krumholz (2012) star formation law, in which the star
formation rate surface density is proportional to the surface density of
molecular gas, $f_{\mathrm{H}_{2}}\Sigma$, times an efficiency factor
$\epsilon_{\mathrm{ff}}\approx 0.01$ (Krumholz & Tan, 2007; Krumholz, 2012),
divided by an estimate of the local effective freefall time.

$\dot{\Sigma}_{\mathrm{SF}}=\epsilon_{\mathrm{ff}}f_{\mathrm{H}_{2}}\Sigma/t_{\mathrm{SF}},$
(7)

The star formation timescale is the shorter of the following two timescales,

$t_{\mathrm{SF}}=\begin{cases}t_{\mathrm{GMC}}=\frac{\pi^{1/4}}{\sqrt{8}}\frac{\sigma}{G(\Sigma_{0}^{3}\Sigma)^{1/4}},&\mathrm{if}\
t_{\mathrm{GMC}}<t_{Q}\\\
t_{\mathrm{Q}}=\frac{1}{\Omega}\sqrt{\frac{3\pi^{2}Q^{2}}{32(\beta+1)\phi_{P}}},&\mathrm{if}\
t_{\mathrm{GMC}}>t_{Q}.\end{cases},$ (8)

where $\phi_{P}$ is a constant that we set to 3, and $\Sigma_{0}$ is $85\
M_{\odot}\ \mathrm{pc}^{-2}$. The molecular gas fraction $f_{\mathrm{H}_{2}}$
is set according to the Krumholz (2013) model, which is metallicity-dependent,
but in the regime of interest where $Q\leq Q_{\mathrm{crit}}$, usually
$f_{\mathrm{H}_{2}}\approx 1$. Finally, the accretion rate is set following
Bouché et al. (2010), who specify simple powerlaw accretion rates as a
function of halo mass and redshift. Halo masses are estimated by inverting the
median Moster et al. (2013) stellar mass-halo mass relation, and the radial
dependence of the accretion rate surface is taken to follow a simple
exponential with a scale length of 0.1 times the Virial radius (Forbes et al.,
2019).

###  2.2 Gravitational Instability

In equilibrium the source terms discussed in the previous subsection will be
cancelled by the effects of gravitational instability, which simultaneously
induces radial transport and heats the disk. One possibility following stellar
N-body simulations (Sellwood & Carlberg, 1984; Carlberg & Sellwood, 1985) is
to take

$\frac{dQ}{dt}\Big{|}_{\mathrm{GI}}=\frac{\mathcal{G}(Q)}{t_{\mathrm{orb}}},$
(9)

where $\mathcal{G}$ is some dimensionless function, and $t_{\mathrm{orb}}=2\pi
r/v_{\mathrm{circ}}$ is the local orbital time. Sellwood & Carlberg (1984)
found $\mathcal{G}=(1/4)(Q_{\mathrm{crit}}-Q)$ when $Q<Q_{\mathrm{crit}}$ and
0 otherwise. This illustrates some of the properties we might expect of
$\mathcal{G}$, namely that $\mathcal{G}(Q=Q_{\mathrm{crit}})\approx 0$, then
monotonically increases as $Q$ decreases. Arguably $\mathcal{G}(Q)$ must also
depend only on $Q$ and not on any other parameters, since $Q$ should totally
encapsulate the degree of instability. In general this reasoning is not
totally certain because the physics of the non-linear evolution of the
stability is not just a function of $Q$, but also the physics of feedback
(Goldbaum et al., 2015, 2016; Behrendt et al., 2016).

An alternative starting point is to assume that the GI term can be modeled as
a viscosity (Balbus & Papaloizou, 1999; Gammie, 2001; Lodato & Rice, 2005;
Krumholz & Burkert, 2010), in which case we can write out the corresponding
rate of change in $Q$ explicitly, equivalently as a function of viscosity, a
Shakura & Sunyaev (1973) $\alpha$ parameter, or a vertically integrated torque
$\mathcal{T}$. These quantities are related via

$\nu=\alpha\sigma^{2}/\Omega$ (10)

which serves as a definition of $\alpha$, and

$\mathcal{T}=2\pi r\Sigma v_{\mathrm{circ}}(\beta-1)\nu$ (11)

(see Shu, 1992, for instance). The change in $Q$ due to a local torque is then

$\frac{dQ}{dt}\Big{|}_{\mathrm{GI}}\approx\frac{\partial
Q}{\partial\sigma}\frac{(\beta-1)v_{\mathrm{circ}}\mathcal{T}}{6\pi
r^{3}\Sigma\sigma}=\frac{\partial
Q}{\partial\sigma}\frac{1}{3}(\beta-1)^{2}\alpha\sigma\Omega,$ (12)

Equation 12 is an approximation because terms related to the advection of
energy through the disk that contribute to $\partial\sigma/\partial t$ have
been dropped, as has a term related to the accumulation (or loss) of mass
because of differential mass flow, proportional to
$\partial^{2}\mathcal{T}/\partial r^{2}$. The two energy terms, proportional
to $\partial\mathcal{T}/\partial r^{2}$ and
$\mathcal{T}\partial\sigma/\partial r$, as well as the mass flow term, are all
zero in the equilibrium solutions found by Krumholz & Burkert (2010) and
Krumholz et al. (2018), and regardless are non-local in the sense that they
depend on the structure of the disk instead of conditions at one radius. The
Krumholz & Burkert (2010) solution was obtained by solving the radial
differential equation for $\mathcal{T}$ under the assumption that
$Q=Q_{\mathrm{crit}}$ exactly (see also Rafikov, 2009; Rice et al., 2011;
Rafikov, 2015, in the context of protoplanetary disks). In other words, the
torque/viscosity/$\alpha$ is set based on global conditions in the disk in
order to maintain $Q=Q_{\mathrm{crit}}$.

One possibility at this stage is to specify a particular $\mathcal{G}(Q)$,
then solve for $\mathcal{T}$ everywhere in the disk simultaneously via a
similar ordinary differential equation in $\mathcal{T}(r)$. This approach was
taken by Forbes et al. (2014) to compute the radial motions and heating of
stars subject to spiral arms in a way that obeyed the Sellwood & Carlberg
(1984) result, and by Krumholz & Forbes (2015) to illustrate numerically how a
disk would evolve towards $Q=Q_{\mathrm{crit}}$. When the advection terms are
dropped, no differential equation is needed and one need only specify an
explicit form of $\mathcal{T}$ or $\alpha$ dependent on local conditions. In
this case we would arrive at

$\mathcal{G}=t_{\mathrm{orb}}\frac{dQ}{dt}\Big{|}_{\mathrm{GI}}\approx\frac{\partial
Q}{\partial\sigma}\frac{2\pi}{3}(\beta-1)^{2}\alpha\sigma\sim\frac{2\pi}{3}(\beta-1)^{2}\alpha
Q$ (13)

where the last approximation applies in the limit where $Q$ is set mostly by
the instability of the gas. This form makes it clear that $\mathcal{G}$
plausibly really is just a function of $Q$ (and the approximately fixed
rotation curve), provided that $\alpha$ itself is only a function of $Q$.

In the protoplanetary disk literature, a variety of $Q$-dependent values of
$\alpha$ have been proposed, though all are fairly ad-hoc, essentially aiming
to satisfy the expectations we had for $\mathcal{G}$, namely that it’s close
to zero at $Q=Q_{\mathrm{crit}}$ and rises monotonically as $Q$ decreases.
Examples include

$\alpha_{\mathrm{Zhu}}=\exp{\left(-Q^{4}\right)}$ (14)

(Zhu et al., 2010b, a) and

$\alpha_{LP}\propto\left(\left(\frac{Q_{\mathrm{crit}}}{Q}\right)^{2}-1\right)$
(15)

(Lin & Pringle, 1987, 1990). Note that these prescriptions for $\alpha(Q)$ are
separate from the prescriptions that arise from a calculation that assumes
$Q=Q_{\mathrm{crit}}$ exactly, which also produce estimates of $\alpha$ (e.g.
Pringle, 1981; Gammie, 2001; Rafikov, 2009; Krumholz & Burkert, 2010), since
these values only apply at the critical value of $Q$ (Kratter & Lodato, 2016).
A compilation of $\alpha$ from simulations (Kratter et al., 2008) provides a
smooth transition in $\alpha$ as $Q$ goes from 2 to 1, with some additional
dependence on the ratio of the disk mass to the stellar mass, but then assumes
a saturation below $Q=1$ since the disk begins to fragment. There is
presumably no recovering from this in a protoplanetary disk, but stellar
feedback and the ability of stars as a large population to self-regulate their
own value of $Q$ mean that galaxies can plausibly equilibrate in this regime.
However, given the variety of $\alpha$ prescriptions available and their
acknowledged ad-hoc nature, the next step is to compare to data.

##  3 Comparison to Data

Figure 1: Data (forcing, $\mathcal{S}$, points) versus theory (GI response, $\mathcal{G}$, lines). Every panel shows some version of $\mathcal{S}=-t_{\mathrm{orb}}dQ/dt|_{\mathrm{sources}}$ as a function of some version of $Q$. The top panel shows a schematic of what we expect each panel to look like: gravitational instability has no effect for $Q>Q_{\mathrm{crit}}$, then acts to increase $Q$ at some unknown rate as a function of $Q$ below $Q_{\mathrm{crit}}$. In equilibrium, the source terms (estimated from the data) will exactly cancel this effect. The points show radius-by-radius estimates of $\mathcal{S}$, i.e. the negative rate of change of $Q$ per dynamical time owing to star formation, accretion, and turbulent dissipation. The red squares are from the DiskMass sample at $z\sim 0$, and the blue circles are from the SINS $z\sim 2$ sample. The lines show prescriptions for $\mathcal{G}=t_{\mathrm{orb}}dQ/dt|_{\mathrm{GI}}$, the rate of change of $Q$ per dynamical time owing to gravitational instability. From highest to lowest at $Q=2$, we show lines corresponding to constant $\alpha=3$ (gray), $\alpha\propto Q^{-0.75}$ (purple), Lin & Pringle (1987) (maroon), $\alpha\propto(2/Q-1)$ (blue), and Zhu et al. (2010a) (red). Uncertainties shown by the colorbars are estimated by varying each point according to Table 1 and taking the standard deviation of the perturbed values of $\mathcal{S}$. The range of $Q_{\mathrm{crit}}$ suggested by Elmegreen (2011), from 2-3, is highlighted as a vertical gray bar in the plots, with higher values of $Q$ shown in a lighter gray. Typical reported statistical errors and estimated systematic errors for the SINS sample are shown as the errorbars in the top left panel. Ingredient for Estimating $\mathcal{S}$ | Perturbation for DiskMass Data | Perturbation for SINS Data
---|---|---
$\Sigma$ | 0.1 dex | 0.3 dex
$\sigma$ | 0.08 dex times log-uniform from 0.5 to 1 | 0.15 dex
$\Sigma_{*}$ | 0.18 dex | 0.3 dex
$\sigma_{*}$ | 0.1 dex times log-uniform from 1 to 2 | 0.48 dex
$v_{\mathrm{circ}}$ | 0.04 dex | 0.15 dex
$\beta=d\ln v_{\mathrm{circ}}/d\ln r$ | 0.1 | 0.1
$\dot{\Sigma}_{\mathrm{SF}}$ | 0.3 dex | 0.3 dex
$Z/Z_{\odot}$ | 0.3 dex | 0.3 dex
$\rho_{\mathrm{dm}}$ | 0.3 dex | 1 dex
$\langle p_{*}/m_{*}\rangle$ | 0.11 dex | 0.3 dex
$\eta$ | 0.3 dex | 1 dex
Table 1: Perturbations to the observational data to estimate the size of each point’s systematic uncertainty. Most perturbations are multiplicative with the exception of those for $\beta$. $\beta$ is additionally clipped to be between -0.5 and 1. $\alpha=3\mathcal{G}/(2\pi(\beta-1)^{2}Q)$ | Reference | Color in Figure 1
---|---|---
3 | - | Gray
$3/Q^{0.75}$ | This work | Purple
$2^{2}/Q^{2}-1$ |  Lin & Pringle (1987), Lin & Pringle (1990) | Maroon
2(2/Q - 1) | This work | Blue
$\exp{\left(-Q^{4}\right)}$ |  Zhu et al. (2010a), Zhu et al. (2010b) | Red
Table 2: Values of $\alpha$ viscosity and the color of their corresponding
lines in Figure 1

We now turn to two published datasets where we can estimate $Q$ and the
magnitude of each source term, and where we know there are cases where $Q$ is
observed to be below plausible values of $Q_{\mathrm{crit}}$. These are the
SINS sample of galaxies (Förster Schreiber et al., 2009), including fits to
the dynamics and mass profiles of 19 galaxies at $z\sim 2$ (Genzel et al.,
2014; Tacchella et al., 2015), and the DiskMass sample (Bershady et al., 2010;
Martinsson et al., 2013b; Westfall et al., 2014) of 30 local nearly face-on
galaxies. Our goal is to estimate the value of $dQ/dt|_{\mathrm{sources}}$ for
each measured annulus in these two galaxy samples, and check whether they can
be cancelled out by $dQ/dt|_{\mathrm{GI}}$ below $Q_{\mathrm{crit}}$.

The SINS data are largely based on H$\alpha$ line emission whose flux,
centroid, and width can be measured. From this a dynamical model is
constructed assuming a velocity dispersion that is constant in radius, and
assuming that the H$\alpha$ emission itself comes from star formation, from
which a surface density of gas can be inferred via inversion of the Kennicutt-
Schmidt relation. We enumerate the published values of $v_{\phi}(r)$,
$\Sigma(r)$, and $Q_{\mathrm{gas}}(r)$, the latter two of which also include
errorbars based on uncertainties in the H$\alpha$ flux. Each galaxy contains
$\sim 20$ radial bins along the main axis of the galaxy for $\Sigma$ and
$Q_{\mathrm{gas}}$. We take each of these points and plug them in to Equation
4 assuming $\eta=1$, $Z=0.5Z_{\odot}$, $\langle p_{*}/m_{*}\rangle$ from Orr
et al. (2022), $\epsilon_{\mathrm{accr}}=0.1$, and
$\dot{\Sigma}_{\mathrm{accr}}$ assuming the exponential scalelength of
accretion is 10% of the Virial radius (Forbes et al., 2019), where the Virial
radius and halo mass are estimated by inverting the median relation between
$M_{*}$ and $M_{h}$ from Moster et al. (2013). The normalization of
$\dot{\Sigma}_{\mathrm{accr}}$ uses the simple powerlaw scaling with halo mass
and $(1+z)$ from Bouché et al. (2010) with an efficiency of 10%. The stellar
surface densities are taken from the single- or two-component surface density
models of Tacchella et al. (2015) derived from H-band photometry. We adopt the
two-component model where available, and the one-component model otherwise. If
neither model is available for a given galaxy, it is excluded from the sample.
We also check that the adopted values of $Q_{\mathrm{gas}}$, $\Sigma$, and
$v_{\mathrm{circ}}$ yield a constant $\sigma$, and exclude any annuli where
this is not the case owing to digitization errors. Since the quoted $\Sigma$
in these observations is inferred from the H$\alpha$ luminosity by assuming a
particular redshift-dependent molecular gas depletion time, but we use the
Krumholz & Dekel (2012) star formation law (see Section 2.1), we use the
quoted $\Sigma$ values to find the original estimated star formation rate
surface density, then use the Krumholz & Dekel (2012) star formation law to
infer $\Sigma$. In other words, we keep the estimate of the star formation
rate surface density unchanged, but infer a different underlying gas column
density. The effect of this change is among the sources of systematic error
that we explore.

Finally, in order to estimate a two-component $Q$, which we will call
$Q_{\mathrm{2c}}$, we need to assume something about the velocity dispersion
of the stars, $\sigma_{*}$. We adopt a smooth interpolation between the case
that $\sigma_{*}\approx\sigma$ as the result of recent star formation, and
$Q_{*}=Q_{*,\mathrm{crit}}=2$, from the action of spiral arm heating. In
particular,

$\frac{\sigma_{*}}{\sigma}=\frac{\mathcal{I}_{\mathrm{spiral}}\sigma_{*,\mathrm{crit}}+\sqrt{\mathcal{I}_{\mathrm{spiral}}\sigma_{*,\mathrm{crit}}^{2}+(\phi^{2}+2\phi\mathcal{I}_{\mathrm{spiral}})\sigma^{2}}}{(\phi+2\mathcal{I}_{\mathrm{spiral}})\sigma}$
(16)

Here $\phi=f_{R}\gamma\ \mathrm{sSFR}\ t_{\mathrm{orb}}$,
$\sigma_{*,\mathrm{crit}}=Q_{*,\mathrm{crit}}\pi G\Sigma_{*}/\kappa$, and
$\mathcal{I}_{\mathrm{spiral}}=1$ when $\sigma<\sigma_{*,\mathrm{crit}}$ and
zero otherwise. This formula is derived by assuming that
$\partial\sigma_{*}/\partial t=0$. In future work we will validate and refine
this estimate, but for now we adopt it and keep this in mind as a source of
systematic uncertainty for our analysis of the SINS galaxies.

In the DiskMass sample, we include all radii that have estimated values of
$\sigma_{z}$, namely the stellar velocity dispersion in the direction
perpendicular to the plane of the disk presented in Martinsson et al. (2013a),
and line-of-sight velocity dispersions of [OIII]$\lambda 5007\rm{\AA}$ as
presented in Martinsson et al. (2013b). For each annulus, we also incorporate
the estimated HI and $\mathrm{H}_{2}$ surface densities from the Westerbork
survey and Spitzer 24 $\mu$m photometry respectively (Bershady et al., 2010).
We estimate the local star formation rate surface density as
$\Sigma_{\mathrm{H}_{2}}/(2\ \mathrm{Gyr})$ (Bigiel et al., 2011). Stellar
mass surface densities are estimated based on K-band photometry, and circular
velocities are adopted from the dynamical modelling (Case Ia) of Martinsson et
al. (2013a). The DiskMass survey corrects their observed line of sight stellar
velocity dispersions to estimate $\sigma_{z}$, but typically the radial
component of the stellar velocity dispersion is larger by up to a factor of
two depending on the physical mechanism by which the stars were heated (e.g.
Lacey, 1984), and it is the radial component of the velocity dispersion that
enters $Q_{*}$. This is another source of systematic uncertainty.

When estimating $Q_{\mathrm{2c}}$ within the DiskMass collaboration (Westfall
et al., 2014), the authors take the observed ionized gas velocity dispersions
and divide by two to estimate the cold gas velocity dispersion (Andersen et
al., 2006; Caldú-Primo et al., 2013). This difference in velocity dispersion
also arises in simulations with resolved multi-phase ISMs (Jeffreson et al.,
2022). In our treatment where we take the gas to be a single component with a
single velocity dispersion, it is not obvious which velocity dispersion to
adopt in general. For the DiskMass sample, usually the molecular gas is
subdominant in terms of mass, so the velocity dispersion of the disk may be
better-represented by the velocity dispersion of the warmer component. Girard
et al. (2021) argue that for $z\sim 2$ galaxy analogues observed at $z\sim 0$,
there is a factor of $>2$ offset between molecular and ionized gas velocity
dispersions. This is more serious since the $z\sim 2$ galaxies are dominated
by molecular gas. On the other hand Übler et al. (2021) found in an analysis
of TNG50 data that analyzing the simulation in the same way as the
observations yields an offset that varies from galaxy to galaxy, with the
observation-like pipeline yielding higher velocity dispersions, though often
just a few km s-1, but occasionally up to a factor of two. By default we
therefore do not adjust the ionized velocity dispersions, but acknowledge that
the velocity dispersions may plausibly be lower by a factor of 2. Keep in mind
that this plausible adjustment factor may be different between the two
samples.

Having now collected all of the assumptions necessary to estimate
$dQ/dt|_{\mathrm{sources}}$ for every annulus in these two samples, we now
plot the results. Figure 1 shows 12 variations of $\mathcal{S}$ vs. $Q$, and a
schematic representation of what we might expect this plot to look like a
priori. We define $\mathcal{S}\equiv-
t_{\mathrm{orb}}dQ/dt|_{\mathrm{sources}}$, which is the rate of change of $Q$
per local dynamical time. The minus sign allows us to directly compare
$\mathcal{S}$ to possibilities for $\mathcal{G}$, the rate of change of $Q$
per dynamical time from the effects of gravitational instability, since in
equilibrium $\mathcal{S}=\mathcal{G}$. Schematically, we expect that
$\mathcal{S}=\mathcal{G}=0$ when $Q>Q_{\mathrm{crit}}$, and for $\mathcal{S}$
and $\mathcal{G}$ to trace out some positive function for
$Q<Q_{\mathrm{crit}}$. Each column shows a different version combination of
$\mathcal{S}$ and $Q$ – the subscripts indicate which version of $Q$ is used
to evaluate either $Q$ itself or in the case of $\mathcal{S}$, the partial
derivatives of $Q$ in Equation 4. Note that to evaluate
$\mathcal{S}_{\mathrm{gas}}$, the last three terms of Equation 4 are all zero
because $Q_{\mathrm{gas}}$ has no explicit dependence on the stellar
component. Each row shows an alternative calculation of $\mathcal{S}$ based on
a plausible change in our assumptions. The second row shows the effect of
ignoring the Orr et al. (2022) model and simply using $\langle
p_{*}/m_{*}\rangle=3000\ \mathrm{km}\ \mathrm{s}^{-1}$, the third row shows a
systematic reduction in observed values of $\sigma$ by a factor of 2, and the
final row shows the adoption of the Genzel et al. (2014) values of $\Sigma$
based on the observed H$\alpha$ flux for those galaxies (the DiskMass sample
is unchanged between the first and fourth row).

Each point in Figure 1 is shaded according to its estimated systematic error,
with larger errorbars shown as lighter colors, and uncertainties greater than
5 also using a much smaller symbol size. The red squares show the DiskMass
sample, and the blue circles show the SINS sample. The errorbars in the top
left panel show the median statistical and systematic uncertainties for the
SINS sample. For each point the systematic error is estimated by re-
calculating $\mathcal{S}$ 10 times with perturbed values of the default input
parameters and recording the standard deviation of the samples. A more
principled approach would be to fit the data with an explicit model and
representing the systematic offsets as parameters to be fit drawn from their
own distributions per annulus and per galaxy, but short of that we are able to
see which points are reliable in this space, and which are more uncertain. The
adopted random perturbations are shown in Table 1. Remarkably much of the
DiskMass sample, particularly the annuli with $Q_{\mathrm{gas}}\gtrsim 3$, is
subject to large uncertainties in $\mathcal{S}_{\mathrm{gas}}$. These arise
largely from the modest uncertainty in the star formation rate, since in this
regime we expect no reduction in $\langle p_{*}/m_{*}\rangle$, and the low
velocity dispersions of the gas make it quite sensitive to the effects of
feedback (see also Section 4). Meanwhile $\mathcal{S}_{\mathrm{2c}}$ for this
sample is insensitive to these uncertainties because $Q_{\mathrm{2c}}$ is
dominated by the stellar component.

For comparison to the data, we plot a number of curves for $\mathcal{G}$,
namely plausible but ad-hoc descriptions of how gravitational instability will
behave as a function of $Q$. The lines are identical in all 8 panels where
they are shown in Figure 1. Many of them arise from a particular formulation
of $\alpha$, but notably there is not a strict one-to-one translation between
$\alpha$ and the two-component version of $\mathcal{G}$, which we could call
$\mathcal{G}_{\mathrm{2c}}$. That is, once the stellar component is included,
an $\alpha$ viscosity that acts by heating and transporting only gas will have
a different effect on $Q_{\mathrm{2c}}$ depending on a wide range of other
parameters. Therefore when we show the same line in the middle panel as we
showed in the left panel, $\alpha$ itself is not the same. Generally the same
$\alpha$ would lead to lower values of $\mathcal{G}_{\mathrm{2c}}$ than
$\mathcal{G}_{\mathrm{gas}}$. Recall from Equation 13 that
$\mathcal{G}_{\mathrm{gas}}=2\pi(\beta-1)^{2}\alpha Q_{\mathrm{gas}}/3$, and
this is how we translate the values of $\alpha$ in Table 2 to the lines shown
in Figure 1.

##  4 Interpretation

We are searching for the function $\mathcal{G}(Q)$ that balances out every
other effect on $Q$ so that $dQ/dt\approx 0$, and which approaches zero for
$Q\gtrsim Q_{\mathrm{crit}}$, whose exact value is also not known. It is also
not obvious what version of $Q$ to use, i.e. the gas-only version, or the
multi-component version. This is why in Figure 1 we show both
$\mathcal{S}_{\mathrm{gas}}$ vs $Q_{\mathrm{gas}}$ and
$\mathcal{S}_{\mathrm{2c}}$ vs $Q_{\mathrm{2c}}$. This ambiguity arises
because $Q$ and its various versions only specify the stability of the disk in
the linear regime with no cooling, whereas per Elmegreen (2011), realistic
disks are never formally stable. We therefore examine the evidence for a
purely gas instability-driven effect, and alternatively an instability in
which the stars also participate.

In the gas-driven case, the relevant panels of Figure 1 are the first column,
showing $\mathcal{S}_{\mathrm{gas}}$ as a function of $Q_{\mathrm{gas}}$. The
data as a whole, largely dominated by the SINS data given the systematic
uncertainties for disks in the DiskMass regime, show a pattern similar to what
we might expect a priori. The data have $\mathcal{S}$ convincingly above zero
until $Q_{\mathrm{gas}}\gtrsim 3$, at which point the errorbars become too
large to reliably discern much of a trend. When $Q_{\mathrm{gas}}\lesssim 2$,
the data are close to a line corresponding to
$\alpha=3/Q_{\mathrm{gas}}^{0.75}$.

In the two-component case, the data quite clearly show that
$\mathcal{S}\approx 0$ for $Q_{\mathrm{2c}}\gtrsim 2$, a regime well-populated
in the DiskMass sample. The SINS data display a larger scatter than in the
gas-driven case, although this may be driven by a combination of greater
systematic uncertainties given the higher redshift and the fact that
$\sigma_{*}$ is not measured but guessed in our analysis. Roughly speaking the
data follow the blue line, corresponding to
$\mathcal{G}_{\mathrm{2c}}=(4\pi/3)(\beta-1)^{2}Q_{\mathrm{2c}}(2/Q_{\mathrm{2c}}-1)$.
In this version we would also expect that $\mathcal{S}_{\mathrm{gas}}$ should
be zero when the instability is inactive, i.e. when
$Q_{\mathrm{2c}}>Q_{\mathrm{crit}}$. Showing this relationship in the third
column of Figure 1, we see that on balance $\mathcal{S}_{\mathrm{gas}}$ is
above zero, though this hinges on only a handful of points with substantial
errorbars.

While these datasets do not immediately tell us exactly which version of $Q$
and $\mathcal{G}$ to use, each prescription is similar in the low-$Q$ regime,
since in these samples low $Q_{\mathrm{2c}}$ only occurs when
$Q_{\mathrm{2c}}$ is dominated by the gas component. Since this is the regime
that originally motivated our investigation into how low $Q$ could go, it is
worth examining $\mathcal{S}$ in this regime to better understand how exactly
the data end up showing $\mathcal{S}\sim 5$. In this limit, namely that $Q$ is
set almost entirely by the gas, $f_{\mathrm{H}_{2}}\approx 1$, and
$t_{\mathrm{SF}}\approx t_{Q}$,

$\displaystyle t_{\mathrm{orb}}\frac{dQ}{dt}\Big{|}_{\mathrm{sources}}\sim$
$\displaystyle
2\epsilon_{\mathrm{ff}}\sqrt{\frac{32(\beta+1)\phi_{P}}{3}}\left(f_{R}+\eta+\frac{1}{3}\frac{\tilde{Q}_{\mathrm{SN}}}{Q}\right)$
(17)
$\displaystyle+\frac{\dot{\Sigma}_{\mathrm{accr}}t_{\mathrm{orb}}}{\Sigma}Q\left(\epsilon_{\mathrm{accr}}\frac{\tilde{Q}_{\mathrm{orb}}^{2}}{Q^{2}}-1\right)$
$\displaystyle-\pi\sqrt{\frac{\beta+1}{2}}\left(1+\left(1+4\zeta_{d}Q^{2}\left(\frac{1}{(\beta+1)}+\frac{8\pi}{Q_{*}^{2}}\right)\right)^{1/2}\right)\left(1-\frac{\tilde{Q}_{\mathrm{therm}}^{2}}{Q^{2}}\right)^{3/2}$

To non-dimensionalize the equation, we have defined several $Q$-like
quantities here, where in each case the velocity dispersion has been replaced
with a different relevant velocity. In particular,
$\tilde{Q}_{\mathrm{orb}}\equiv v_{\mathrm{orb}}\kappa/(\pi G\Sigma)$ and
$\tilde{Q}_{\mathrm{SN}}\equiv\langle p_{*}/m_{*}\rangle\kappa/(\pi G\Sigma)$.
In both cases, the ratio of these $\tilde{Q}$’s to $Q$ is often $\gg 1$,
though given the large velocity dispersions in these galaxies
$\tilde{Q}_{\mathrm{orb}}/Q$ may be less than 1 in galaxy centers. As in
Forbes et al. (2014), we include a factor in the final term containing
$\tilde{Q}_{\mathrm{therm}}\equiv\sigma_{\mathrm{therm}}\kappa/(\pi G\Sigma)$
to reduce the dissipation term to zero as
$\sigma\rightarrow\sigma_{\mathrm{therm}}$, i.e. where the velocity dispersion
arises purely from the thermal soundspeed of the Warm Neutral Medium.

Figure 2: The dependence of star formation, turbulent dissipation, and
gravitational instability on $Q$. The dashed lines showing the effect of
gravitational instability as encapsulated by
$\mathcal{G}=t_{\mathrm{orb}}dQ/dt|_{\mathrm{GI}}$ are the same in each panel,
and vary quite dramatically as a function of $\beta$. We have adopted a
particular form of $\mathcal{G}$ referred to in the text as the two-component
case. In equilibrium, $\mathcal{G}=\mathcal{S}$, shown as the solid lines.
These equilibria for each value of $\beta$ are highlighted as three circles in
each panel. The total sources $\mathcal{S}$ are the sum of the effect of star
formation (dot-dashed lines), which always acts to increase $Q$, and
dissipation (dotted lines), which always act to decrease $Q$.

In this limit, the loss term ranges from about 2-7, with a dependence on $Q$
that is at most linear under the right circumstances. Meanwhile star
formation, both via the removal of gas and the injection of momentum, always
opposes the loss term in a $Q$-dependent way, so that star formation becomes
more effective at lower $Q$. In some models this is the end of the story –
$\epsilon_{\mathrm{ff}}$ changes such that these two terms exactly cancel at
$Q=Q_{\mathrm{crit}}$ (e.g. Faucher-Giguère et al., 2013), but if
$\epsilon_{\mathrm{ff}}$ is constant (Krumholz & Tan, 2007; Krumholz, 2012),
they will only cancel at one value of $Q$ depending on the gas fraction, the
mass loading factor, and the momentum injection from SN feedback. The
remaining term representing the effects of accretion may be positive or
negative, since energy, but also mass, is added to the system. While in
principle it appears that this term could dominate in the outer parts of disks
where $\tilde{Q}_{\mathrm{orb}}/Q\gg 1$, $t_{\mathrm{orb}}$ is long, and
$\Sigma$ is low, in that regime $\sigma$ is usually already at the thermal
sound speed when $Q=Q_{\mathrm{crit}}$, meaning that in practice this term
does not appear to matter very often.

The gas-dominated limit shown in Equation 17 is useful in understanding Figure
2 where we show how the source terms tend to vary as a function of $Q$ for a
few different parameter combinations, along with the other side of the
equation, $\mathcal{G}=t_{\mathrm{orb}}(dQ/dt)|_{\mathrm{GI}}$. For the latter
we have set $\mathcal{G}=(4\pi/3)(\beta-1)^{2}Q(2/Q-1)$, which is a reasonable
description of the data. The different colors correspond to different values
of the rotation curve slope $\beta$: black, red, and blue are respectively
$-0.45$, $0.25$, and $0.95$. These values have a huge effect on $\mathcal{G}$,
but only modest effects on the source terms. This is because the viscosity is
much less effective when there is little inherent shear in the disk (as
$\beta\rightarrow 1$ and the disk approaches solid body rotation).

Generally the source terms increase with increasing $Q$ as the influence of
star formation wanes and the loss term becomes large. The source term on its
own can be zero when star formation is important, and the loss term can be
zero when $Q$ is so low that, at fixed values of $v_{\mathrm{circ}}$, $r$,
$\beta$, and $\Sigma$, the velocity dispersion reaches the thermal soundspeed
of the disk (taken here to be 8 km s-1). For low values of $\beta$ where the
disk experiences substantial shear, $\mathcal{G}$ is large enough to intersect
the source terms at $Q$ only modestly less than $Q_{\mathrm{crit}}$, whereas
at high values of $\beta$, it takes much lower values of $Q$ before
$\mathcal{G}$ can counteract the source terms for this choice of viscosity.
Note also that the intersections of the source and GI terms are almost always
stable, since for example to the left of the equilibrium point,
$\mathcal{S}<\mathcal{G}$, so $dQ/dt>0$, which would move the disk back
towards the equilibrium point.

With this understanding of the gas-dominated regime, it is clear that the SINS
data simply represent a set of galaxies where the turbulent dissipation rate
is large enough that it cannot be balanced by the effects of star formation,
even when stellar feedback is able to add a momentum per unit mass of 3000 km
s-1 (see second row of Figure 1) or when the velocity dispersions are
substantially over-estimated (see the third row of Figure 1). This is very
similar to previous results that have shown that high velocity dispersions
cannot be maintained by supernovae alone (e.g. Joung et al., 2009; Brucy et
al., 2020), but here we show that the effect of gravitational instability is
likely not so strong that $Q\rightarrow Q_{\mathrm{crit}}$. Instead, $Q$
remains in equilibrium at values far less than 1.

##  5 Implications

Figure 3: Equilibrium values of $Q_{\mathrm{2c}}$. Each panel shows a slice
through parameter space. If a parameter is not being varied, it is set as
follows: $\beta=0$, $v_{\mathrm{circ}}=200\ \mathrm{km}\ \mathrm{s}^{-1}$,
$r=5\ \mathrm{kpc}$, $\Sigma=\Sigma_{*}=1000M_{\odot}\ \mathrm{pc}^{-2}$. The
red region indicates that $Q$ cannot fall below $Q_{\mathrm{crit}}$ because
$\sigma$ is too low. Substantial parts of parameter space admit equilibrium
values of $Q$ well below $Q_{\mathrm{crit}}$, particularly for high values of
$\Sigma$ and $\beta$.

By adopting a value of $\mathcal{G}(Q)$ that balances out the observed source
terms, encapsulated as $\mathcal{S}$, we can check how varying the conditions
in the disk will change $Q_{\mathrm{eq}}$, the value of $Q$ where
$\mathcal{S}=\mathcal{G}$. We compute this solution numerically for
$Q_{\mathrm{2c}}$, adopting Equation 4 for the sources and
$\mathcal{G}=(4\pi/3)(\beta-1)^{2}Q(2/Q-1)$ for the effect of gravitational
instability. We then fix all but two values of $r$, $\sigma_{*}/\sigma$,
$\Sigma$, $\Sigma_{*}$, $v_{\mathrm{circ}}$, and $\beta$, and plot the
resulting equilibrium $Q_{\mathrm{2c}}$ in Figure 3. The fixed values are
$r=5\ \mathrm{kpc}$, $\sigma_{*}/\sigma=1.1$, $\Sigma_{*}=10^{3}M_{\odot}\
\mathrm{pc}^{-2}$, $v_{\mathrm{circ}}=200\ \mathrm{km}\ \mathrm{s}^{-1}$, and
$\beta=0.1$. Large radii, small $v_{\mathrm{circ}}$, large $\Sigma_{*}$, or
large $\beta$, will all, so long as $\Sigma\gtrsim 10M_{\odot}\
\mathrm{pc}^{-2}$, produce $Q_{\mathrm{eq}}$ noticeably lower than the assumed
$Q_{\mathrm{crit}}=2$. The hatched region in red indicates that in that part
of parameter space, $\sigma$ has reached the thermal sound speed at
$Q=Q_{\mathrm{crit}}$, so $Q$ cannot fall below $Q_{\mathrm{crit}}$ there (and
indeed these regions are surrounded by regions of $Q_{\mathrm{eq}}\approx
Q_{\mathrm{crit}}$).

The patterns in Figure 3 can largely be understood by referring back to Figure
2 and Equation 17. The simplest to understand is $\beta$, where increasing
values decrease the shear in the disk, leading to weaker effects on $dQ/dt$
for a given value of $\alpha$. Interestingly this behavior with $\beta$
suggests that the phenomenon we are studying here, namely equilibrium values
of $Q$ well below $Q_{\mathrm{crit}}$, is much harder to generate in
protoplanetary disks than galactic disks, since in the former shear is always
strong given the Keplerian rotation curves of $\beta=-1/2$. Variation in $r$
is a stand-in for variation in $\Omega$ and hence $\kappa$, which scale every
version of $Q$ or $\tilde{Q}$ by the same factor, essentially stretching or
narrowing many of the terms shown in Figure 2 along the $x$-axis. Higher
values of $\sigma_{*}/\sigma$ also tend to lower $Q_{\mathrm{eq}}$, sometimes
quite sharply, as $Q_{*}$ becomes less important in determining
$Q_{\mathrm{2c}}$. Essentially once $Q_{*}$ is important in setting
$Q_{\mathrm{2c}}$, the partial derivatives of $Q$ with respect to $\Sigma$ and
$\sigma$ become much smaller, and so the dramatic effects of star formation,
feedback, and turbulent dissipation, are much weaker (see Equation 4). A
similar effect is visible in the lower right panel showing variation in
$\Sigma_{*}$.

While it is straightforward to compute $Q_{\mathrm{eq}}$ numerically once we
adopt a particular function to represent $\mathcal{G}$ and have chosen a
particular value of $Q_{\mathrm{crit}}$, it bears repeating that the data we
have compiled here do not make it obvious precisely what $Q_{\mathrm{crit}}$
is, nor the value of $\mathcal{G}$ near $Q_{\mathrm{crit}}$, nor even which
version of $Q$ to use. This dataset may be augmented by other surveys and by
simulations, and the analysis may be improved by casting it as a hierarchical
Bayesian inference problem.

The ambiguities discussed here lead to uncertainty in how to implement our
results here into simplified disk models that rely on $Q_{\mathrm{crit}}$ in
some way (e.g. Krumholz & Burkert, 2010; Forbes et al., 2012; Genel et al.,
2012; Cacciato et al., 2012; Forbes et al., 2014; Porter et al., 2014;
Krumholz & Forbes, 2015; Stevens et al., 2016; Krumholz & Burkhart, 2016;
Krumholz et al., 2018; Forbes et al., 2019). Clearly imposing that $Q\geq
Q_{\mathrm{crit}}$ everywhere at all times is not adequate. However, the two
different models for gravitational instability-induced changes in $Q$, namely
the gas-only version and the multi-component version, may imply somewhat
different prescriptions in simplified models. The gas-only version is likely
the simplest, since there there is a direct one-to-one relationship between
$\mathcal{G}$ and $\alpha$, so one could simply adopt
$\alpha_{\mathrm{GI}}\approx 3/Q_{\mathrm{gas}}^{0.75}$ when
$Q_{\mathrm{gas}}\lesssim 2.5$ and $\approx 0$ otherwise. This $\alpha$ would
have a local effect on $d\sigma/dt$ and $\dot{M}$, the rate of mass flow
through the disk. The multi-component version where
$\mathcal{G}_{\mathrm{2c}}=(4\pi/3)(\beta-1)^{2}Q_{\mathrm{2c}}(2/Q_{\mathrm{2c}}-1)$
is more ambiguous because a single value of $\alpha$ leads to a plethora of
values of $\mathcal{G}_{\mathrm{2c}}$ depending on local conditions, and in
general there is no guarantee that any value of $\alpha$ will be large enough
to generate the required $\mathcal{G}_{\mathrm{2c}}$ when the stellar
component is important in setting $Q_{\mathrm{2c}}$. Essentially because stars
may participate in the instability in this formulation, they too may be
altered by $\mathcal{G}$. The flatness of the data in the third column of
Figure 1 suggests that one way to implement this may be to set
$\alpha_{\mathrm{GI}}\propto 1/Q_{\mathrm{2c}}$ when
$Q_{\mathrm{2c}}<Q_{\mathrm{crit}}$ and 0 otherwise, then calculating the
$\mathcal{G}_{*}\equiv t_{\mathrm{orb}}dQ_{*}/dt|_{\mathrm{GI}}$ necessary to
reach the given $\mathcal{G}_{\mathrm{2c}}$. Note that in general
$\mathcal{G}_{\mathrm{2c}}\neq\mathcal{G}_{\mathrm{gas}}+\mathcal{G}_{*}$
because

$\mathcal{G}_{\mathrm{2c}}=t_{\mathrm{orb}}\left(\frac{\partial
Q_{\mathrm{2c}}}{\partial Q_{\mathrm{gas}}}\frac{\partial
Q_{\mathrm{gas}}}{\partial t}+\frac{\partial Q_{\mathrm{2c}}}{\partial
Q_{*}}\frac{\partial Q_{*}}{\partial t}+\frac{\partial
Q_{\mathrm{2c}}}{\partial s}\frac{\partial s}{\partial
t}\right)\Big{|}_{\mathrm{GI}},$ (18)

where $s$ is the ratio $\sigma/\sigma_{*}$. While backing out
$\mathcal{G}_{*}$ from this equation is slightly more arduous than just
subtracting off $\mathcal{G}_{\mathrm{gas}}$, it is straightforward to
compute. Moreover this multi-component formulation of the effect of
gravitational instability on $Q$ is more appealing theoretically simply
because the stability of a multi-component disk definitely depends on both
components in the linear regime.

##  6 Conclusion

We attempt to reconcile the common theoretical assumption that $Q\geq
Q_{\mathrm{crit}}$ with observations of $z\sim 2$ galaxies that show $Q\ll 1$
almost ubiquitously. This can be accomplished by relaxing this rigid
theoretical assumption slightly and setting the viscosity according to a
$Q$-dependent $\alpha$. Despite substantial observational uncertainties, a
simple analytic form $\alpha\approx 3/Q_{\mathrm{gas}}^{0.75}$ for a gas-only
version of the instability, or
$dQ/dt|_{\mathrm{GI}}=(4\pi/3)(\beta-1)^{2}(2-Q_{\mathrm{2c}})t_{\mathrm{orb}}^{-1}$
for the multi-component version, would cancel out the rate of change of
$Q_{\mathrm{2c}}$ due to other physical effects. Of these physical effects,
the most important is the high rate of turbulent dissipation, tempered by the
effects of star formation, implied by the high surface densities and velocity
dispersions of the $z\sim 2$ galaxies. This provides a simple explanation for
the low observed values of $Q$, and a path forward for improving the
implementation of the effects of gravitational instability in simplified
models.

We thank Rachel Somerville, Shy Genel, Mark Krumholz, Doug Lin, Chris Hayward,
and Julianne Dalcanton, for helpful conversations. JCF is supported by a
Flatiron Research Fellowship through the Flatiron Institute, a division of the
Simons Foundation.

## References

  * Andersen et al. (2006) Andersen, D. R., Bershady, M. A., Sparke, L. S., et al. 2006, The Astrophysical Journal Supplement Series, 166, 505, doi: 10.1086/506609
  * Balbus & Papaloizou (1999) Balbus, S. A., & Papaloizou, J. C. B. 1999, The Astrophysical Journal, 521, 650, doi: 10.1086/307594
  * Behrendt et al. (2016) Behrendt, M., Burkert, A., & Schartmann, M. 2016, The Astrophysical Journal, 819, L2, doi: 10.3847/2041-8205/819/1/L2
  * Bershady et al. (2010) Bershady, M. A., Verheijen, M. A. W., Swaters, R. A., et al. 2010, The Astrophysical Journal, 716, 198, doi: 10.1088/0004-637X/716/1/198
  * Bigiel et al. (2011) Bigiel, F., Leroy, A. K., Walter, F., et al. 2011, The Astrophysical Journal, 730, L13, doi: 10.1088/2041-8205/730/2/L13
  * Bouché et al. (2010) Bouché, N., Dekel, A., Genzel, R., et al. 2010, The Astrophysical Journal, 718, 1001, doi: 10.1088/0004-637X/718/2/1001
  * Bournaud & Elmegreen (2009) Bournaud, F., & Elmegreen, B. G. 2009, The Astrophysical Journal, 694, L158, doi: 10.1088/0004-637X/694/2/L158
  * Brucy et al. (2020) Brucy, N., Hennebelle, P., Bournaud, F., & Colling, C. 2020, The Astrophysical Journal, 896, L34, doi: 10.3847/2041-8213/ab9830
  * Cacciato et al. (2012) Cacciato, M., Dekel, A., & Genel, S. 2012, Monthly Notices of the Royal Astronomical Society, 421, 818, doi: 10.1111/j.1365-2966.2011.20359.x
  * Caldú-Primo et al. (2013) Caldú-Primo, A., Schruba, A., Walter, F., et al. 2013, The Astronomical Journal, 146, 150, doi: 10.1088/0004-6256/146/6/150
  * Carlberg & Sellwood (1985) Carlberg, R. G., & Sellwood, J. A. 1985, The Astrophysical Journal, 292, 79, doi: 10.1086/163134
  * Davies et al. (2011) Davies, R., Förster Schreiber, N. M., Cresci, G., et al. 2011, The Astrophysical Journal, 741, 69, doi: 10.1088/0004-637X/741/2/69
  * Dekel et al. (2009a) Dekel, A., Sari, R., & Ceverino, D. 2009a, The Astrophysical Journal, 703, 785, doi: 10.1088/0004-637X/703/1/785
  * Dekel et al. (2009b) Dekel, A., Birnboim, Y., Engel, G., et al. 2009b, Nature, 457, 451, doi: 10.1038/nature07648
  * Elmegreen (2011) Elmegreen, B. G. 2011, The Astrophysical Journal, 737, 10, doi: 10.1088/0004-637X/737/1/10
  * Elmegreen et al. (2007) Elmegreen, D. M., Elmegreen, B. G., Ravindranath, S., & Coe, D. A. 2007, The Astrophysical Journal, 658, 763, doi: 10.1086/511667
  * Faucher-Giguère et al. (2013) Faucher-Giguère, C.-A., Quataert, E., & Hopkins, P. F. 2013, Monthly Notices of the Royal Astronomical Society, 433, 1970, doi: 10.1093/mnras/stt866
  * Forbes et al. (2012) Forbes, J., Krumholz, M., & Burkert, A. 2012, The Astrophysical Journal, 754, 48, doi: 10.1088/0004-637X/754/1/48
  * Forbes et al. (2014) Forbes, J. C., Krumholz, M. R., Burkert, A., & Dekel, A. 2014, Monthly Notices of the Royal Astronomical Society, 438, 1552, doi: 10.1093/mnras/stt2294
  * Forbes et al. (2019) Forbes, J. C., Krumholz, M. R., & Speagle, J. S. 2019, Monthly Notices of the Royal Astronomical Society, 487, 3581, doi: 10.1093/mnras/stz1473
  * Forbes et al. (2022) Forbes, J. C., Emami, R., Somerville, R. S., et al. 2022, Gas Accretion Can Drive Turbulence in Galaxies 2022arXiv220405344F
  * Förster Schreiber et al. (2009) Förster Schreiber, N. M., Genzel, R., Bouché, N., et al. 2009, The Astrophysical Journal, 706, 1364, doi: 10.1088/0004-637X/706/2/1364
  * Förster Schreiber et al. (2018) Förster Schreiber, N. M., Renzini, A., Mancini, C., et al. 2018, The Astrophysical Journal Supplement Series, 238, 21, doi: 10.3847/1538-4365/aadd49
  * Gammie (2001) Gammie, C. F. 2001, The Astrophysical Journal, 553, 174, doi: 10.1086/320631
  * Genel et al. (2012) Genel, S., Dekel, A., & Cacciato, M. 2012, Monthly Notices of the Royal Astronomical Society, 425, 788, doi: 10.1111/j.1365-2966.2012.21652.x
  * Genzel et al. (2014) Genzel, R., Förster Schreiber, N. M., Lang, P., et al. 2014, The Astrophysical Journal, 785, 75, doi: 10.1088/0004-637X/785/1/75
  * Genzel et al. (2017) Genzel, R., Förster Schreiber, N. M., Übler, H., et al. 2017, Nature, 543, 397, doi: 10.1038/nature21685
  * Girard et al. (2021) Girard, M., Fisher, D. B., Bolatto, A. D., et al. 2021, The Astrophysical Journal, 909, 12, doi: 10.3847/1538-4357/abd5b9
  * Goldbaum et al. (2015) Goldbaum, N. J., Krumholz, M. R., & Forbes, J. C. 2015, The Astrophysical Journal, 814, 131, doi: 10.1088/0004-637X/814/2/131
  * Goldbaum et al. (2016) —. 2016, The Astrophysical Journal, 827, 28, doi: 10.3847/0004-637X/827/1/28
  * Gurvich et al. (2020) Gurvich, A. B., Faucher-Giguère, C.-A., Richings, A. J., et al. 2020, Monthly Notices of the Royal Astronomical Society, 498, 3664, doi: 10.1093/mnras/staa2578
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hayward & Hopkins (2017) Hayward, C. C., & Hopkins, P. F. 2017, Monthly Notices of the Royal Astronomical Society, 465, 1682, doi: 10.1093/mnras/stw2888
  * Hobbs et al. (2015) Hobbs, A., Read, J., & Nicola, A. 2015, Monthly Notices of the Royal Astronomical Society, 452, 3593, doi: 10.1093/mnras/stv1469
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jeffreson et al. (2022) Jeffreson, S. M. R., Sun, J., & Wilson, C. D. 2022, Monthly Notices of the Royal Astronomical Society, 515, 1663, doi: 10.1093/mnras/stac1874
  * Joung et al. (2009) Joung, M. R., Mac Low, M.-M., & Bryan, G. L. 2009, The Astrophysical Journal, 704, 137, doi: 10.1088/0004-637X/704/1/137
  * Kassin et al. (2012) Kassin, S. A., Weiner, B. J., Faber, S. M., et al. 2012, The Astrophysical Journal, 758, 106, doi: 10.1088/0004-637X/758/2/106
  * Kereš et al. (2005) Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, Monthly Notices of the Royal Astronomical Society, 363, 2, doi: 10.1111/j.1365-2966.2005.09451.x
  * Kim & Ostriker (2015) Kim, C.-G., & Ostriker, E. C. 2015, The Astrophysical Journal, 802, 99, doi: 10.1088/0004-637X/802/2/99
  * Klessen & Hennebelle (2010) Klessen, R. S., & Hennebelle, P. 2010, Astronomy and Astrophysics, Volume 520, id.A17, 18 pp., 520, A17, doi: 10.1051/0004-6361/200913780
  * Kratter & Lodato (2016) Kratter, K., & Lodato, G. 2016, Annual Review of Astronomy and Astrophysics, 54, 271, doi: 10.1146/annurev-astro-081915-023307
  * Kratter et al. (2008) Kratter, K. M., Matzner, C. D., & Krumholz, M. R. 2008, The Astrophysical Journal, 681, 375, doi: 10.1086/587543
  * Krumholz & Burkert (2010) Krumholz, M., & Burkert, A. 2010, The Astrophysical Journal, 724, 895, doi: 10.1088/0004-637X/724/2/895
  * Krumholz (2012) Krumholz, M. R. 2012, The Astrophysical Journal, 759, 9, doi: 10.1088/0004-637X/759/1/9
  * Krumholz (2013) —. 2013, Monthly Notices of the Royal Astronomical Society, 436, 2747, doi: 10.1093/mnras/stt1780
  * Krumholz & Burkhart (2016) Krumholz, M. R., & Burkhart, B. 2016, Monthly Notices of the Royal Astronomical Society, 458, 1671, doi: 10.1093/mnras/stw434
  * Krumholz et al. (2018) Krumholz, M. R., Burkhart, B., Forbes, J. C., & Crocker, R. M. 2018, Monthly Notices of the Royal Astronomical Society, 477, 2716, doi: 10.1093/mnras/sty852
  * Krumholz & Dekel (2012) Krumholz, M. R., & Dekel, A. 2012, The Astrophysical Journal, 753, 16, doi: 10.1088/0004-637X/753/1/16
  * Krumholz & Forbes (2015) Krumholz, M. R., & Forbes, J. C. 2015, Astronomy and Computing, Volume 11, p. 1-17., 11, 1, doi: 10.1016/j.ascom.2015.02.005
  * Krumholz & Tan (2007) Krumholz, M. R., & Tan, J. C. 2007, The Astrophysical Journal, 654, 304, doi: 10.1086/509101
  * Lacey (1984) Lacey, C. G. 1984, Monthly Notices of the Royal Astronomical Society, 208, 687, doi: 10.1093/mnras/208.4.687
  * Leitner & Kravtsov (2011) Leitner, S. N., & Kravtsov, A. V. 2011, The Astrophysical Journal, 734, 48, doi: 10.1088/0004-637X/734/1/48
  * Leroy et al. (2008) Leroy, A. K., Walter, F., Brinks, E., et al. 2008, The Astronomical Journal, 136, 2782, doi: 10.1088/0004-6256/136/6/2782
  * Lin & Pringle (1987) Lin, D. N. C., & Pringle, J. E. 1987, The Astrophysical Journal Letters, 320, L87, doi: 10.1086/184981
  * Lin & Pringle (1990) —. 1990, The Astrophysical Journal, 358, 515, doi: 10.1086/169004
  * Lodato & Rice (2005) Lodato, G., & Rice, W. K. M. 2005, Monthly Notices of the Royal Astronomical Society, 358, 1489, doi: 10.1111/j.1365-2966.2005.08875.x
  * Mac Low et al. (1998) Mac Low, M.-M., Klessen, R. S., Burkert, A., & Smith, M. D. 1998, Physical Review Letters, 80, 2754, doi: 10.1103/PhysRevLett.80.2754
  * Martinsson et al. (2013a) Martinsson, T. P. K., Verheijen, M. A. W., Westfall, K. B., et al. 2013a, Astronomy & Astrophysics, Volume 557, id.A131, 52 pp., 557, A131, doi: 10.1051/0004-6361/201321390
  * Martinsson et al. (2013b) —. 2013b, Astronomy and Astrophysics, 557, A130, doi: 10.1051/0004-6361/201220515
  * Moster et al. (2013) Moster, B. P., Naab, T., & White, S. D. M. 2013, Monthly Notices of the Royal Astronomical Society, 428, 3121, doi: 10.1093/mnras/sts261
  * Oppenheimer et al. (2010) Oppenheimer, B. D., Davé, R., Kereš, D., et al. 2010, Monthly Notices of the Royal Astronomical Society, 406, 2325, doi: 10.1111/j.1365-2966.2010.16872.x
  * Orr et al. (2022) Orr, M. E., Fielding, D. B., Hayward, C. C., & Burkhart, B. 2022, The Astrophysical Journal, 932, 88, doi: 10.3847/1538-4357/ac6c26
  * Ostriker & Kim (2022) Ostriker, E. C., & Kim, C.-G. 2022, The Astrophysical Journal, 936, 137, doi: 10.3847/1538-4357/ac7de2
  * Ostriker et al. (2010) Ostriker, E. C., McKee, C. F., & Leroy, A. K. 2010, The Astrophysical Journal, 721, 975, doi: 10.1088/0004-637X/721/2/975;
  * Porter et al. (2014) Porter, L. A., Somerville, R. S., Primack, J. R., & Johansson, P. H. 2014, Monthly Notices of the Royal Astronomical Society, 444, 942, doi: 10.1093/mnras/stu1434
  * Price et al. (2020) Price, S. H., Kriek, M., Barro, G., et al. 2020, The Astrophysical Journal, 894, 91, doi: 10.3847/1538-4357/ab7990
  * Price et al. (2021) Price, S. H., Shimizu, T. T., Genzel, R., et al. 2021, The Astrophysical Journal, 922, 143, doi: 10.3847/1538-4357/ac22ad
  * Pringle (1981) Pringle, J. E. 1981, Annual Review of Astronomy and Astrophysics, 19, 137, doi: 10.1146/annurev.aa.19.090181.001033
  * Rafikov (2001) Rafikov, R. R. 2001, Monthly Notices of the Royal Astronomical Society, 323, 445, doi: 10.1046/j.1365-8711.2001.04201.x
  * Rafikov (2009) —. 2009, The Astrophysical Journal, 704, 281, doi: 10.1088/0004-637X/704/1/281
  * Rafikov (2015) —. 2015, The Astrophysical Journal, 804, 62, doi: 10.1088/0004-637X/804/1/62
  * Rathaus & Sternberg (2016) Rathaus, B., & Sternberg, A. 2016, Monthly Notices of the Royal Astronomical Society, 458, 3168, doi: 10.1093/mnras/stw416
  * Rice et al. (2011) Rice, W. K. M., Armitage, P. J., Mamatsashvili, G. R., Lodato, G., & Clarke, C. J. 2011, Monthly Notices of the Royal Astronomical Society, 418, 1356, doi: 10.1111/j.1365-2966.2011.19586.x
  * Romeo & Wiegert (2011) Romeo, A. B., & Wiegert, J. 2011, Monthly Notices of the Royal Astronomical Society, 416, 1191, doi: 10.1111/j.1365-2966.2011.19120.x
  * Sellwood & Carlberg (1984) Sellwood, J. A., & Carlberg, R. G. 1984, The Astrophysical Journal, 282, 61, doi: 10.1086/162176
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, Astronomy and Astrophysics, Vol. 500, p. 33-51 (2009), 500, 33
  * Shu (1992) Shu, F. H. 1992, The Physics of Astrophysics. Volume II: Gas Dynamics.
  * Stevens et al. (2016) Stevens, A. R. H., Croton, D. J., & Mutch, S. J. 2016, Monthly Notices of the Royal Astronomical Society, 461, 859, doi: 10.1093/mnras/stw1332
  * Stone et al. (1998) Stone, J. M., Ostriker, E. C., & Gammie, C. F. 1998, The Astrophysical Journal Letters, 508, L99, doi: 10.1086/311718
  * Tacchella et al. (2015) Tacchella, S., Lang, P., Carollo, C. M., et al. 2015, The Astrophysical Journal, 802, 101, doi: 10.1088/0004-637X/802/2/101
  * Thompson et al. (2005) Thompson, T. A., Quataert, E., & Murray, N. 2005, The Astrophysical Journal, 630, 167, doi: 10.1086/431923
  * Tinsley (1980) Tinsley, B. M. 1980, Fundamentals of Cosmic Physics, 5, 287
  * Toomre (1964) Toomre, A. 1964, The Astrophysical Journal, 139, 1217, doi: 10.1086/147861
  * Übler et al. (2019) Übler, H., Genzel, R., Wisnioski, E., et al. 2019, The Astrophysical Journal, 880, 48, doi: 10.3847/1538-4357/ab27cc
  * Übler et al. (2021) Übler, H., Genel, S., Sternberg, A., et al. 2021, Monthly Notices of the Royal Astronomical Society, 500, 4597, doi: 10.1093/mnras/staa3464
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wang & Silk (1994) Wang, B., & Silk, J. 1994, The Astrophysical Journal, 427, 759, doi: 10.1086/174182
  * Westfall et al. (2014) Westfall, K. B., Andersen, D. R., Bershady, M. A., et al. 2014, The Astrophysical Journal, 785, 43, doi: 10.1088/0004-637X/785/1/43
  * Wisnioski et al. (2015) Wisnioski, E., Förster Schreiber, N. M., Wuyts, S., et al. 2015, The Astrophysical Journal, 799, 209, doi: 10.1088/0004-637X/799/2/209
  * Zhu et al. (2010a) Zhu, Z., Hartmann, L., & Gammie, C. 2010a, The Astrophysical Journal, 713, 1143, doi: 10.1088/0004-637X/713/2/1143
  * Zhu et al. (2010b) Zhu, Z., Hartmann, L., Gammie, C. F., et al. 2010b, The Astrophysical Journal, 713, 1134, doi: 10.1088/0004-637X/713/2/1134