Mon. Not. R. Astron. Soc. 000, 000–000 (0000) Printed 31 May 2010 (MN LATEX style file v2.2)
How Do Massive Black Holes Get Their Gas?
Philip F. Hopkins1 & Eliot Quataert1
1Department of Astronomy and Theoretical Astrophysics Center, University of California Berkeley, Berkeley, CA 94720
Submitted to MNRAS, Dec. 14, 2009
ABSTRACT
We use multi-scale smoothed particle hydrodynamic simulations to study the inflow of gas
from galactic scales ( 10kpc) down to . 0:1pc, at which point the gas begins to resemble
a traditional, Keplerian accretion disk. The key ingredients of the simulations are gas, stars,
black holes (BHs), self-gravity, star formation, and stellar feedback (via a subgrid model);
BH feedback is not included. We use  100 simulations to survey a large parameter space of
galaxy properties and subgrid models for the interstellar medium physics. We generate initial
conditions for our simulations of galactic nuclei (. 300 pc) using galaxy scale simulations,
including both major galaxy mergers and isolated bar-(un)stable disk galaxies. For sufficiently
gas-rich, disk-dominated systems, we find that a series of gravitational instabilities generates
large accretion rates of up to  1 10M yr
1 onto the BH (i.e., at . 0:1 pc); this is com-
parable to what is needed to fuel the most luminous quasars. The BH accretion rate is highly
time variable for a given set of conditions in the galaxy at  kpc. At radii & 10 pc, our sim-
ulations resemble the “bars within bars” model of Shlosman et al, but we show that the gas
can have a diverse array of morphologies, including spirals, rings, clumps, and bars; the duty
cycle of these features is modest, complicating attempts to correlate BH accretion with the
morphology of gas in galactic nuclei. At  1 10 pc, the gravitational potential becomes
dominated by the BH and bar-like modes are no longer present. However, we show that the
gas can become unstable to a standing, eccentric disk or a single-armed spiral mode (m = 1),
in which the stars and gas precess at different rates, driving the gas to sub-pc scales (again
for sufficiently gas-rich, disk-dominated systems). A proper treatment of this mode requires
including star formation and the self-gravity of both the stars and gas (which has not been
the case in many previous calculations). Our simulations predict a correlation between the
BH accretion rate and the star formation rate at different galactic radii. We find that nuclear
star formation is more tightly coupled to AGN activity than the global star formation rate of a
galaxy, but a reasonable correlation remains even for the latter.
Key words: galaxies: active — quasars: general — galaxies: evolution — cosmology: theory
1 INTRODUCTION
The inflow of gas into the central parts of galaxies plays a crit-
ical role in galaxy formation, ultimately generating phenomena
as diverse as bulges and spheroidal galaxies, starbursts and ultra-
luminous infrared galaxies (ULIRGs), nuclear stellar clusters, and
accretion onto super-massive black holes (BHs). The discovery, in
the past decade, of tight correlations between black hole mass and
host spheroid properties including mass (Kormendy & Richstone
1995; Magorrian et al. 1998), velocity dispersion (Ferrarese & Mer-
ritt 2000; Gebhardt et al. 2000), and binding energy or potential
well depth (Hopkins et al. 2007; Aller & Richstone 2007) implies
that these phenomena are tightly coupled.
It has long been realized that bright, high-Eddington ratio ac-
 E-mail:phopkins@astro.berkeley.edu
cretion (i.e., a quasar) dominates the accumulation of mass in the
supermassive BH population (Soltan 1982; Salucci et al. 1999;
Shankar et al. 2004). In order to explain the existence of black holes
with masses  109M, the amount of gas required is comparable
to that contained in entire large galaxies. Given the short lifetime of
the quasar phase . 108 yr (Martini 2004), the processes of interest
must deliver a galaxy’s worth of gas to the inner regions of a galaxy
on a relatively short timescale.
There is also compelling evidence that quasar activity is pre-
ceded and/or accompanied by a period of intense star formation in
galactic nuclei (Sanders et al. 1988a,b; Dasyra et al. 2007; Kauff-
mann et al. 2003). The observed properties of bulges at z 0 inde-
pendently require that dissipative processes (gas inflow) dominate
the formation and structure of the innerkpc (Ostriker 1980; Carl-
berg 1986; Gunn 1987; Kormendy 1989; Hernquist et al. 1993).
Hopkins et al. (2009a,d, 2008a) showed that this inner dissipational
c
0000 RAS
ar
X
iv
:0
91
2.
32
57
v2
[
as
tro
-p
h.C
O]
2
8 M
ay
20
10

2 Hopkins and Quataert
component can constitute a large fraction 530% of the galaxy’s
mass, with stellar (and at some point probably gas) surface densi-
ties reaching  101112 M kpc
2.
On large (galactic) scales, several viable processes for initi-
ating such inflows are well-known. Major galaxy-galaxy mergers
produce strong non-axisymmetric disturbances of the constituent
galaxies; such disturbances may also be produced in some mi-
nor mergers and/or globally self-gravitating isolated galactic disks.
Observationally, major mergers are associated with enhancements
in star formation in ULIRGs, sub-millimeter galaxies, and pairs
more generally (e.g. Sanders & Mirabel 1996; Schweizer 1998; Jo-
gee 2006; Dasyra et al. 2006; Woods et al. 2006; Veilleux et al.
2009). Numerical simulations of mergers have shown that when
such events occur in gas-rich galaxies, resonant tidal torques lead to
rapid inflow of gas into the central kpc (Hernquist 1989; Barnes
& Hernquist 1991, 1996). The resulting high gas densities trigger
starbursts (Mihos & Hernquist 1994, 1996), and are presumed to
feed rapid black hole growth. Feedback from the starburst and a
central active galactic nucleus (AGN) may also be important, both
for regulating the BH’s growth (Di Matteo et al. 2005; Hopkins
et al. 2005; DeBuhr et al. 2009; Johansson et al. 2009b) and for
shutting down future star formation (Springel et al. 2005a; Johans-
son et al. 2009a; see, however, DeBuhr et al. 2009).
However, the physics of how gas is transported from  1 kpc
to much smaller scales remains uncertain (e.g., Goodman 2003).
Typically, once gas reaches sub-kpc scales, the large-scale torques
produced by a merger and/or large-scale bar/spiral become less ef-
ficient. In the case of stellar bars or spiral waves, there can even be
a “hard” barrier to further inflow in the form of an inner Linblad
resonance, if the system has a non-trivial bulge. In mergers, the co-
alescence of the two systems generates perturbations on all scales,
and so allows gas to move through the resonances, but the pertur-
bations relax rapidly on small scales, often before gas can inflow.
Local viscous stresses – which are believed to dominate angu-
lar momentum transport near the central BH (e.g., Balbus & Haw-
ley 1998) – are inefficient at radii & 0:01 0:1 pc (e.g., Shlosman
& Begelman 1989; Goodman 2003; Thompson et al. 2005). It is
in principle possible that some molecular clouds could be scattered
onto very low angular momentum orbits, but even the optimistic
fueling rates from this process are generally insufficient to produce
luminous quasars (see e.g. Hopkins & Hernquist 2006; Kawakatu
& Wada 2008; Nayakshin & King 2007). As a consequence, many
models invoke some form of gravitational torques (“bars within
bars”; Shlosman et al. 1989) to continue transport to smaller radii.
As gas is driven into the central kpc by large-scale torques, it will
cool rapidly into a disky structure; if this gas reservoir is massive
enough, the gas will be self-gravitating and thus again vulnerable
to global instabilities (e.g., the well-known bar and/or spiral wave
instabilities) that can drive some of the gas to yet smaller radii.
To date, numerical simulations have seen the formation of
such secondary bars in some circumstances, such as in adaptive
mesh refinement (AMR) simulations of galaxy formation (Wise
et al. 2007; Levine et al. 2008; Escala 2007) or particle-splitting
smoothed particle hydrodynamics (SPH) simulations of some ide-
alized systems (Escala et al. 2004; Mayer et al. 2007). These stud-
ies have served as a critical “proof of concept.” However, these ex-
amples have generally been limited by computational expense to
studying a single system at one instant in its evolution, and thus it
is difficult to assess how the sub-pc dynamics depends on the large
parameter space of possible inflow conditions from large radii and
galaxy structural parameters.
Alternatively, some simulations simply take an assumed
small-scale structure and/or fixed inflow rate as an initial/boundary
condition, and study the resulting gas dynamics at small radii (e.g.
Schartmann et al. 2009; Dotti et al. 2009; Wada & Norman 2002;
Wada et al. 2009). These studies have greatly informed our under-
standing of nuclear obscuration on small scales (the “torus”), and
the role of stellar feedback in determining the structure and dy-
namics of the gas at these radii; it is, however, unclear how to relate
this small-scale dynamics to the larger-scale properties of the host
galaxy. This is critical for understanding black hole growth and nu-
clear star formation in the broader context of galaxy formation.
Observationally, a long standing puzzle has been that many
systems, especially those with weaker inflows on large scales (e.g.
bar or spiral wave-unstable disks with some bulge, as opposed to
major mergers), exhibit no secondary instabilities at  0:11 kpc
– in several cases, torques clearly reverse sign inside these radii
(Block et al. 2001; García-Burillo et al. 2005). Whether this is
generic, or the consequence of a low duty cycle, or the result of
the large-scale inflows simply being too weak in these cases, is not
clear. Moreover, even among systems that do show nuclear asym-
metries, and that clearly exhibit enhanced star formation and lu-
minous AGN, the observed features at smaller radii are very often
not traditional bars. Rather, they exhibit a diverse morphology, with
spirals quite common, along with nuclear rings, barred rings, occa-
sional one or three-armed modes, and some clumpy/irregular struc-
tures (Martini & Pogge 1999; Peletier et al. 1999; Knapen et al.
2000; Laine et al. 2002; Knapen et al. 2002; Greene et al. 2009).
Even if secondary bars or spirals are present at intermediate
radii  10 100 pc, it has long been recognized that they will
cease to be important at yet smaller scales, when the potential be-
comes quasi-Keplerian and the global self-gravity of the gas less
important; this occurs as one approaches the BH radius of influ-
ence, which is 10pc in typical L galaxies (Athanassoula et al.
1983, 2005; Shlosman et al. 1989; Heller et al. 2001; Begelman &
Shlosman 2009). Indeed, in previous simulations and most analytic
calculations, the “bars-within-bars” model appears to break down
at these scales (see e.g. Jogee 2006, and references therein). How-
ever, local angular momentum transport is still very inefficient at
 10 pc, and the gas is still locally self-gravitating, and so should be
able to form stars rapidly (e.g., Thompson et al. 2005). Understand-
ing the physics of inflow through these last few pc, especially in a
consistent model that connects to gas on galactic scales ( 10kpc),
remains one of the key open questions in our understanding of mas-
sive BH growth.
In this paper, we present a suite of multi-scale hydrodynamic
simulations that follow gravitational torques and gas inflow from
the kpc scales of galaxy-wide events through to < 0:1pc where the
material begins to form a standard thin accretion disk. These sim-
ulations include gas cooling, star formation, and self-gravity; feed-
back from supernovae and stellar winds is crudely accounted for via
a subgrid model. In order to isolate the physics of angular momen-
tum transport, we do not include BH feedback in our calculations.
We systematically survey a large range of galaxy properties (e.g.,
gas fraction and bulge to disk ratio) and gas thermodynamics, in
order to understand how these influence the dynamics, inflow rates,
and observational properties of gas on small scales in galactic nu-
clei ( 0:1 100 pc). Our focus in this paper is on the results of
most observational interest: what absolute inflow rates, star forma-
tion rates, and gas/stellar surface density profiles result from sec-
ondary gravitational instabilities? What is their effective duty cy-
cle? And what range of observational morphologies are predicted?
In a future paper (Paper II) we will present a more detailed compar-
ison between our numerical results and analytic models of inflow
c
0000 RAS, MNRAS 000, 000–000

How Do Massive Black Holes Get Their Gas? 3
and angular momentum transport induced by non-axisymmetric in-
stabilities in galactic nuclei.
The remainder of this paper is organized as follows. In § 2 we
describe our simulation methodology, which consists of two levels
of “re-simulations” using initial conditions motivated by galactic-
scale simulations. In § 3 we present an overview of our results and
show how a series of gravitational processes leads to gas transport
from galactic scales to sub-pc scales. In §4 we quantify the resulting
inflow rates and gas properties as a function of time and radius in
the simulations. §5 summarizes the conditions required for global
gravitational instability and significant gas inflow. In § 6 we show
how the physics of accretion induced by gravitational instabilities
leads to a correlation (with significant scatter) between star for-
mation at different radii and BH accretion; we also compare these
results to observations. In § 7 we summarize our results and discuss
a number of their implications and several additional observational
tests. Further numerical details and tests of our methodology are
discussed in § A. In § B we show how the subgrid model of the
ISM we use influences our results.
2 METHODOLOGY
We use a suite of hydrodynamic simulations to study the physics
of gas inflow from  10kpc to  0:1pc in galactic nuclei. In or-
der to probe the very large range in spatial and mass scales, we
carry out a series of “re-simulations.” First, we simulate the dynam-
ics on galaxy scales. Specifically, we use representative examples
of gas-rich galaxy-galaxy merger simulations and isolated, moder-
ately bar-unstable disk simulations. These are well-resolved down
to  100 500pc. We use the conditions at these radii (at several
times) as the initial conditions for intermediate-scale re-simulations
of the sub-kpc dynamics. In these re-simulations, the smaller vol-
ume is simulated at higher resolution, allowing us to resolve the
subsequent dynamics down to 10pc scales – these re-simulations
approximate the nearly instantaneous behavior of the gas on sub-
kpc scales in response to the conditions atkpc set by galaxy-scale
dynamics. We then repeat our re-simulation method to follow the
dynamics down to sub-pc scales where the gas begins to form a
standard accretion disk.
Our re-simulations are not intended to provide an exact re-
alization of the small-scale dynamics of the larger-scale simula-
tion that motivated the initial conditions of each re-simulation (in
the manner of particle-splitting or adaptive-mesh refinement tech-
niques). Rather, our goal is to identify the dominant mechanism(s)
of angular momentum transport in galactic nuclei and what param-
eters they depend on. This approach clearly has limitations, espe-
cially at the outer boundaries of the simulations; however, it also
has a major advantage. By not requiring the conditions at small
radii to be uniquely set by a larger-scale “parent” simulation, we
can run a series of simulations with otherwise identical conditions
(on that scale) but systematically vary one parameter (e.g., gas frac-
tion or ISM model) over a large dynamic range. This allows us to
identify the physics and galaxy properties that have the biggest ef-
fect on gas inflow in galactic nuclei. As we will show, the diversity
of behaviors seen in the simulations, and desire to marginalize over
the uncertain ISM physics, makes such a parameter survey critical.
This methodology is discussed in more detail below. First, we
describe the physics in our simulations, in particular our treatment
of gas cooling, star formation, and feedback from supernovae and
young stars (§ 2.1). We then summarize the galaxy-scale simula-
tions that are used to motivate the initial conditions for subsequent
re-simulations (§ 2.2). The intermediate-scale re-simulations, and
the methodology used to construct their initial conditions, are dis-
cussed in § 2.3. Finally, we discuss the nuclear-scale resimulations,
which are themselves motivated by the intermediate-scale resimu-
lations (§ 2.4).
2.1 Gas Physics, Star Formation, and Stellar Feedback
The simulations were performed with the parallel TreeSPH code
GADGET-3 (Springel 2005), based on a fully conservative formula-
tion of smoothed particle hydrodynamics (SPH), which conserves
energy and entropy simultaneously even when smoothing lengths
evolve adaptively (see e.g., Springel & Hernquist 2002; Hernquist
1993; O’Shea et al. 2005). The detailed numerical methodology
is described in Springel (2005), Springel & Hernquist (2003), and
Springel et al. (2005b).
The simulations include supermassive black holes (BHs) as
additional collisionless particles at the centers of all progenitor
galaxies. In our calculations the BH’s only dynamical role is via
its gravitational influence on the smallest scales  1 10 pc. In
particular, to cleanly isolate the physics of gas inflow, we do not in-
clude the subgrid models for BH accretion and feedback that have
been used in previous works (e.g., Springel et al. 2005b). During
a galaxy merger, the BHs in each galactic nucleus are assumed to
coalesce and form a single BH at the center of mass of the system
once they are within a single SPH smoothing length of one another
and are moving at a relative speed lower than both the local gas
sound speed and relative escape velocities.
In our models, stars form from the gas using a prescription
motivated by the observed Kennicutt (1998) relation. Specifically,
we use a star formation rate per unit volume _ / 3=2 with the
normalization chosen so that a Milky-way like galaxy has a total
star formation rate of about 1M yr
1.
The precise slope, normalization, and scatter of the Schmidt-
Kennicutt relation, and even whether or not such a relation is gen-
erally applicable, are somewhat uncertain on the smallest spatial
scales we model here. This is especially true when the dynamical
times become short relative to the main-sequence stellar lifetime
(tdyn  105 106 yr in the smallest regions simulated). Nonethe-
less, there is some observational and physical motivation for the
“standard” parameters we have adopted, even at high surface den-
sities. For the densest star forming galaxies, observational studies
favor a logarithmic slope' 1:7 for the relation between _ and g
(Bouché et al. 2007), not that different from what our model imple-
ments. In addition, Tan et al. (2006) and Krumholz & Tan (2007)
show that local observations imply a constant star formation effi-
ciency in units of the dynamical time (i.e. _ / 1:5) at all densities
observed, n 1016 cm3 – the highest gas densities in these stud-
ies are comparable to the highest gas densities in our simulations
( 108 M of gas inside  10 pc). Finally, Davies et al. (2007) &
Hicks et al. (2009) estimate the star formation rate (SFR) and gas
surface densities in AGN on exactly the small scales of interest here
( 110pc); they find a SFR-density relation continuous with that
implied at “normal” galaxy densities.
To understand the possible impact of uncertainties in the
Schmidt-Kennicutt relation on our conclusions, we have adjusted
the slope d ln _=dln adopted in our simulations between 1:02:0
in a small set of test runs, fixing the star formation rate at MW-
like surface densities of ' 109 M kpc
2 where the observational
constraints are tight. This amounts to varying the absolute star for-
mation efficiency on the smallest resolved scales by a factor of
& 100; qualitatively, this could presumably mimic a wide variety
of different physics associated with stellar feedback and star for-
c
0000 RAS, MNRAS 000, 000–000

4 Hopkins and Quataert
10-2 100 102 104 106 108 1010
Bulk Average Density n [cm-3]
0
50
100
150
200
250
Ef
fe
cti
ve
S
ou
nd
S
pe
ed
c
s
[km
s-1
]
qeos=1qeos=0.25
qeos=0.125
qeos=0
Adiabatic
(No Cooling)
Figure 1. Effective equation of state of the ISM in our simulations; this accounts for the effects of stellar feedback that are not resolved in our calculations. We
plot the effective sound speed cs (i.e., the turbulent speed generated by feedback) versus average ISM density n. The qeos = 0 case is an isothermal “floor” at
cs = 10kms1. The qeos = 1 line is the “maximal feedback” model in Springel & Hernquist (2003), in which the ISM is multi-phase above a minimum n and
all supernova energy goes into pressurizing the “diffuse” medium. Intermediate qeos interpolate between the two (eq. 1). For each qeos, we show an analytic
curve for the equilibrium cs(n), and simulation results (colored points). We compile measurements of the turbulent or non-thermal velocities in observed
systems (black points): local ULIRGs (Downes & Solomon 1998, circles; large/small for the inner/outer observed radii in each), the nuclei of merging LIRGs
(Bryant & Scoville 1999, star), the core of M82 (Westmoquette et al. 2007, ), the central 400 pc of NGC 6240 (Tacconi et al. 1999, triangle), nuclear
clumps in NGC 6240 (Iono et al. 2007, pentagon), and the maser disk of NGC 3079 (Kondratko et al. 2005, small squares). The canonical spiral disk velocity
dispersion is also shown (large square). At high-redshifts, we show recent observations of “normal” star-forming disks from Förster Schreiber et al. (2006,
z 2; diamond) and Lemoine-Busserolle et al. (2009, z = 34; asterisk), and luminous sub-millimeter galaxies from Tacconi et al. (2006, inverted triangle).
The observations favor a median qeos  0:1250:3, which we adopt as our “fiducial” models.
mation. This variation can, unsurprisingly, have a dramatic affect
on the quasi-equilibrium gas densities at small radii, which are set
by gas inflow balancing star formation. However, even over this
large range of star formation efficiencies, the qualitative behavior
of the angular momentum transport and gas inflow does not change
dramatically; the gas dynamics in a low-star formation efficiency
run is similar to that in a run with much higher initial gas content
but also higher star formation efficiency. As a result, although the
absolute star formation efficiency is clearly somewhat uncertain,
we do not believe that this qualitatively affects our conclusions.
The largest uncertainties in our modeling stem from the treat-
ment of the interstellar medium (ISM) gas physics and the im-
pact of stellar feedback on the ISM. Our simulations are rela-
tively coarse and average over many star-forming clumps, HII re-
gions, supernova remnants, etc. As a result, the simulations use a
sub-resolution model of a multiphase ISM in which the gas has a
sound speed much larger than its true thermal velocity (Springel
& Hernquist 2003). Our assumption is that the large-scale grav-
itational torques produced by bars, spiral waves, and other non-
axisymmetric features, will not depend critically on the small-scale
structure of the ISM; although we believe that this is qualitatively
correct, more detailed calculations will be required to ultimately
assess this assumption. The key role of stellar feedback in this
model is to suppress the runaway fragmentation and clumping of
gas on small scales. In reality, this likely occurs via turbulence gen-
erated by stellar feedback and via the disruption of star clusters and
molecular clouds (e.g., Murray et al. 2009). In our model, all of
this physics is “accounted for” by the large effective sound speed,
which increases the Jeans and Toomre masses, thus suppressing the
formation of small-scale structure.
Figure 1 shows the range of effective sound speeds cs in our
calculations as a function of the ISM density n, compared to a num-
ber of observational constraints (large symbols). The solid lines in
Figure 1 represent analytic approximations to the equation of state,
while the small colored points are representative results from sim-
ulations that also include adiabatic cooling/heating and shock heat-
ing. We adopt a parameterization of the sound speed as a function
of density – i.e. the effective equation of state for the ISM gas –
c
0000 RAS, MNRAS 000, 000–000

How Do Massive Black Holes Get Their Gas? 5
following Springel et al. (2005b); Springel & Hernquist (2005);
Robertson et al. (2006a,b). With this model, we can interpolate
freely between two extremes using a parameter qeos. At one ex-
treme, the gas has an effective sound speed of 10kms1, motivated
by, e.g., the observed turbulent velocity in atomic gas in nearby
spirals or the sound speed of low density photo-ionized gas; this
is the “no-feedback” case with qeos = 0. The opposite extreme,
qeos = 1, represents the “maximal feedback” sub-resolution model
of Springel et al. (2005b), based on the multiphase ISM model
of McKee & Ostriker (1977); in this case, 100% of the energy
from supernovae is assumed to stir up the ISM. This equation of
state is substantially stiffer, with effective sound speeds as high as
 200kms1. Note that at the highest densities, cs begins to de-
cline in all of the models (albeit slowly), as the efficiency of star
formation asymptotes but cooling rates continue to increase.
Varying qeos between these two extremes amounts to varying
the effective sound speed of the ISM, with the interpolation
cs =
q
qeos c2s [q = 1] + (1qeos)c2s [q = 0] : (1)
The resulting sound speeds for qeos = 0:125 and 0:25 are shown in
Figure 1; these correspond to more moderate values of cs  30
100kms1 for the densities of interest.
Figure 1 compares these models to observations of the turbu-
lent (non-thermal) velocities in atomic and molecular gas in a num-
ber of systems (large symbols). At low mean densities, n . 0:3
1cm3, the turbulent velocity in nearby spirals is  10kms1.
Downes & Solomon (1998) present a detailed study of a number
of luminous local starbursts that have significantly higher mean
densities; they decompose the observed molecular line profiles into
bulk (e.g., rotation) and turbulent motions. We plot their determi-
nation of the mean density and turbulent velocities in each system
at several radii. We also show the results of similar observations
of the core of M82 (Westmoquette et al. 2007), additional nearby
luminous infrared galaxies (Bryant & Scoville 1999), NGC 6240
(Tacconi et al. 1999; Iono et al. 2007), and luminous starbursts at
high redshift, z 23 (Förster Schreiber et al. 2006; Tacconi et al.
2006; Lemoine-Busserolle et al. 2009); at the highest densities, we
also show the random velocities observed in the nuclear maser disk
in the nearby Seyfert 2 galaxy NGC3079 (Kondratko et al. 2005).
The observational results in Figure 1 favor models with qeos 
0:1 0:3, albeit with significant scatter. We thus take these values
of qeos as our “standard” choices, although we have carried out nu-
merical experiments over the entire range qeos = 01. Note that the
observations clearly do not support a simple no-feedback (qeos = 0)
model. Within the range qeos  0:10:3, our results on AGN fuel-
ing are not particularly sensitive to the precise value of qeos. More-
over, the functional form cs(n) is also not crucial: simulations using
a constant cs ' 50kms1 yield similar results. However, our simu-
lations with qeos = 0 and qeos = 1 predict results that are inconsistent
with observations of galactic nuclei – thus, our results themselves
favor qeos  0:10:3 (see Appendix B).
For the gas densities of interest in this paper, the precise form
of the cooling law does not significantly affect our conclusions.
This is because the cooling time is almost always much shorter
than the local dynamical time (typical tcool  106104 tdyn). As
a result the "sound speed" of the gas is nearly always pinned to
the subgrid ‘turbulent’ value discussed above (this is why the nu-
merical points in Fig. 1 are so close to the analytic models). This
is true even when the gas is optically thick to the infrared radia-
tion produced by dust, as can readily occur in the central  100
pc: the cooling time (diffusion time) is still much less than the dy-
namical time in the optically thick limit for the radii that we re-
solve (e.g., Thompson et al. 2005). We have, in fact, experimented
with alternative cooling rate prescriptions: including or excluding
metal-line cooling, uniformly increasing or decreasing the cooling
rate by a factor of  3, and in the most extreme case, assuming
instantaneous gas cooling (any gas parcel above the cooling floor
is assumed to immediately radiate the excess energy in a single
timestep). We do not see any significant changes in our results with
these variations, simply because the gas always cools rapidly in
our calculations (in contrast, in regimes such as the -disk where
the cooling time is comparable to the dynamical time, the details
of the cooling function can have a significant effect; see Gammie
2001; Nayakshin et al. 2007; Cossins et al. 2009). If the effective
minimum cs comes from turbulent velocities, then the ÒeffectiveÓ
cooling time around this SSoor should be given by the turbulent de-
cay time, which can be comparable to the dynamical time (Begel-
man & Shlosman 2009); this is not included in our calculations. In
the presence of such an effective cooling time, local gravitational
instability may lead to tightly wound spirals as opposed to frag-
mentation into star-forming clumps. These could be important for
angular momentum transport at some radii.
To conclude our discussion of the ISM physics in our sim-
ulations, it is important to reiterate that the key role of the sub-
resolution sound speed cs is that it determines the local Jeans and
Toomre criteria, and thus the physical scale on which gravitational
physics dominates. The Jeans mass for a disk of surface density
 and sound speed cs is given by MJ = ( c4s )=(4G
2 ). For the
outer regions of a galactic disk with  108109 M kpc
2 and
cs  10kms1, MJ  106 M, comparable to that of a molecular
cloud; the corresponding Jeans length is tens of pc, comparable to
that of massive molecular cloud complexes in galaxies. Thus our
sub-resolution model is effectively averaging over discrete molec-
ular clouds and star clusters in galaxies. Large-scale inflows can
increase the surface density in the central regions of galaxies, but
cs also rises. In our models with qeos ' 0:10:3, the Jeans mass re-
mains roughly similar down to pc scales, but as a result the Jeans
length is significantly smaller in galactic nuclei where the ambient
density is much higher. These physical mass and size-scales mo-
tivate the numerical resolution in our simulations; in all cases, we
ensure that the resolution is sufficient to formally resolve the Jeans
mass and length. Higher resolution simulations may be numerically
achievable, but can provide only minimal gains in the “reality” of
the simulation without a corresponding increase in the sophistica-
tion of the ISM model.
2.2 Large Scale Galaxy Mergers and Bars: 100 kpc to 100 pc
Our galaxy-scale simulations motivate the initial conditions chosen
for the smaller-scale re-simulation calculations described in §2.3 &
2.4. The galaxy-scale simulations include isolated disks (both glob-
ally stable and bar unstable) and galaxy-galaxy mergers. We will
ultimately focus on a few representative examples, but we chose
those having surveyed a large parameter space. These simulations
and the methodology used for building the initial galaxies are de-
scribed in more detail in a series of papers (see e.g. Di Matteo et al.
2005; Robertson et al. 2006b; Cox et al. 2006; Younger et al. 2008).
We briefly review the key points here.
For each simulation, we generate one or two stable, isolated
disk galaxies, each with an extended dark matter halo with a Hern-
quist (1990) profile, an exponential disk of gas and stars, and an
optional stellar bulge. The initial systems are chosen to be consis-
tent with the observed baryonic Tully-Fisher relation and estimated
halo-galaxy mass scaling laws (Bell & de Jong 2001; Kormendy &
c
0000 RAS, MNRAS 000, 000–000

6 Hopkins and Quataert
Table 1. Galaxy-Scale Simulations
Simulation  qeos fgas h=R B=T (< R) d(0) b(0) Mtot(< R) [M]
Name [pc] (500 pc) (500 pc) (300 pc) (1 kpc) [M kpc
2] (300 pc) (1 kpc)
(1) (2) (3) (4) (5) (6) (7) (8) (9)
Merger & Isolated Galaxy Simulations: Parameter Studies
— 20,50,100,150 0.0-1.0 0.05-1.0 0.05-0.35 0.1-0.9 0.01-0.5 1.0e9-1.0e11 1.0e9-1.0e12 1.0e8-1.0e11 1.0e8-1.0e11
Typical Gas-Rich  L Merger: Initial Conditions
b3ex(ic) 10a, 50 0.25 0.80 0.25 0.8 0.4 2.5e8 2.5e8 1.4e8 1.9e9
Typical Gas-Rich  L Merger: At Coalescence
b3ex(co) 10a, 50 0.25 0.45 0.33 0.15 0.25 2.0e10 2.6e9 4.4e9 2.5e10
Typical Gas-Rich  L Merger: 108 yr Post-Coalescence
b3ex(po) 10a, 50 0.25 0.08 0.37 0.10 0.15 2.3e11 9.1e9 2.8e10 4.6e10
Bar-Unstable  L Disk: Initial Conditions
barex(ic) 10a, 50 0.25 0.40 0.15 0.55 0.35 1.7e8 2.4e8 1.5e8 1.9e9
Bar-Unstable  L Disk: At Peak of Inflow
barex(pk) 10a, 50 0.25 0.48 0.20 0.13 0.23 1.1e9 2.5e8 5.2e8 3.8e9
Bar-Unstable  L Disk: After Bar Relaxation
barex(re) 10a, 50 0.25 0.04 0.11 0.16 0.18 1.2e10 5.5e9 5.6e9 1.5e10
Parameters describing representative examples of our galaxy-scale simulations of galaxy-galaxy mergers and unstable isolated disks. The top row gives the
range spanned in each parameter across our suite of simulations. Subsequent rows pertain to specific examples at chosen times in their evolution. (1)
Simulation name/ID. (2) Minimum smoothing length (in pc). (3) Equation of state parameter (Figure 1). (4) Gas fraction Mgas=(Mgas + M; disk) of the
disky/cold component (inside the given radius). (5) Scale height of the disky/cold component within the given radius (median jzj=R; the plane of the disk is
defined by the total angular momentum vector). (6) Bulge-to-total mass ratio inside a given radius (bulge here includes all spherical components: stellar
bulge, dark matter halo, and black hole). (7) Maximum surface density of the disky/cold component (gas plus stars) – for an exponential profile, the surface
density is nearly constant at radii below the disk scale length. Otherwise, averaged over a couple times our minimum smoothing length. (8) Average surface
density inside 300 pc for the bulge component. (9) Total mass enclosed inside a given radius. All simulations include black holes, but these are dynamically
unimportant on these scales.
a For each of the two simulations here, we have also run three ultra high-resolution simulations which also act as a moderate-resolution intermediate-scale
simulations ( 107 particles). They are identical in initial conditions to the standard merger and isolated disk run here, with initial gas fraction equal to,
one-half, and one-quarter that shown (six simulations in total).
Freeman 2004; Mandelbaum et al. 2006, and references therein).
The galaxies have total masses Mvir = V 3vir=(10GH[z]) for an initial
redshift z, with the baryonic disk having a mass fraction md (typi-
cally md ' 0:041) relative to the total mass. The system has an ini-
tial bulge-to-total baryonic mass ratio B=T , and the disk has initial
gas fraction fgas. The dark matter halos are assigned a concentration
parameter scaled as in Robertson et al. (2006b) for the galaxy mass
and redshift following Bullock et al. (2001). Disk scale lengths are
set in accordance with the above scaling laws.
In previous papers (referenced above), a large suite of these
simulations have been presented, with several hundred simula-
tions of varying equation of state, numerical resolution, merger
orbital parameters, structural properties (e.g. profile shapes, initial
bulge-to-disk ratios, and scale lengths), initial gas fractions, and
halo concentrations. In this suite, galaxies have baryonic masses
 1081013 M and gas fractions fgas = 01; mergers spanning
mass ratios from 1:1 to 1:20, and isolated disks have Toomre Q
parameters from 0:110.
In this work, we focus on galaxies with baryonic masses
 1011 M. Based on the survey above, we select a representa-
tive simulation of a gas-rich major merger and one of an isolated
disk, to provide the basis for our subsequent re-simulations. Some
of the salient parameters of these simulations are given in Table 1.
The merger is equal-mass, with 1011 M galaxies that have gas
fractions of  40% at the time of the merger/coalescence. The or-
bit is a moderately tilted prograde, parabolic case (orbit e in Cox
et al. 2006). Together this makes for a fairly violent, gas rich major
merger, representative of many of our other gas-rich major merger
simulations at both low and high redshifts. The isolated system is a
1011 M disk with fgas = 0:4, B=T = 0:3, scale length h = 3:2kpc,
and md = 0:041; it has a Toomre Q of order one and develops a
moderate bar (amplitude  15%), but the gas encounters an inner
Linblad resonance at 1-2 kpc.
Small variations in the orbits or the structural properties of the
galaxies will change the details of the tidal and bar features on large
scales. However, the precise details of these large-scale simulations
are not important for our study of the dynamics on small scales (see
Appendix A). Rather, the small-scale dynamics depends on global
parameters such as the gas mass channeled into the central region,
relative to the pre-existing bulge, disk, and black hole mass. In these
respects, we have chosen the simulations summarized in Table 1 to
be representative of a broad class of gas-rich systems.
Our galaxy-scale simulations have spatial resolution – grav-
itational softening length and minimum adaptive SPH smoothing
length – of 50pc. In the suite described above, the resolution scales
with galaxy mass and is 50100pc for M  1011 M systems,
but in a subset of higher-resolution cases is as small as 20pc. In
Hopkins et al. (2008b) and Hopkins et al. (2009a) we have demon-
strated that this resolution is sufficient to properly resolve not only
the mass fractions but also the spatial extent of the “starburst”
formed from gas which loses angular momentum in a merger or
via a strong bar instability. However, to assess how much of this
gas can ultimately fuel a central BH requires that we determine the
dynamics on even smaller spatial scales.
c
0000 RAS, MNRAS 000, 000–000

How Do Massive Black Holes Get Their Gas? 7
2.3 Intermediate Scales: Re-Simulating from 1 kpc to 10 pc
In order to follow the behavior of gas inflow on smaller scales, we
re-simulate the central regions of interest at higher resolution, in a
series of progressively smaller-scale runs. We begin by selecting a
number of representative outputs from the galaxy-scale simulations
described above, near the peak of activity. We select several snap-
shots in the gas-rich merger at key epochs: early close encounters
of the two galaxies, just at nuclear coalescence (which is the peak
of star formation in the nuclear region), and at the “end” of the
merger (roughly  108 yr after the final coalescence). We also se-
lect snapshots typical of isolated, moderately bar-unstable systems,
at times where a bar and some inflow has developed; for compar-
ison, we also consider a fully stable (pre-bar) galaxy disk. In each
case, we focus on the central 0:11 kpc region, which includes the
majority of the gas that has been driven in from larger scales. Some
of the representative properties of these snapshots, at these scales
and times, are outlined in Table 1.
Our approach to re-simulating the nuclear region is to use
the larger-scale simulations to motivate the initial conditions of a
smaller scale calculation (a “zoom-in” or “re-simulation”). We do
so by de-composing the potential, density, and velocity distribu-
tions of the gas, stars, and dark matter at a given time in the larger-
scale simulation using the basis expansion proposed in Hernquist &
Ostriker (1992). This allows us to not only re-construct a smoothed
density profile, but also to include the asymmetric structures from
the larger-scale simulation (if desired) and to define where the po-
tential is noise-dominated.1 From these stellar, gas, and dark matter
distributions, we re-populate the gas and stars in the central regions
(the scale we wish to re-simulate; generally out to an outer radius of
 1 2kpc) and use this as the initial condition for a new simula-
tion that we run for several local dynamical times. To be conserva-
tive, we typically initialize only a small amount of gas in the inner
parts of the re-simulation2, since the larger-scale simulation from
which the initial condition is drawn has little information about the
gas properties on small scales; in Appendix A we show that the
subsequent dynamics does not depend significantly on these details
of the initial conditions.
We have carried out a total of  50 simulations at these inter-
mediate scales, which together span a wide range in the key param-
eters: the equation of state of the gas and the relative mass fraction
in a pre-existing bulge, gas disk, and stellar disk. Table 2 summa-
rizes the key properties of physical importance in several of these
simulations (some numerical studies and surveys of initial condi-
tion, which turn out to have little effect on the key results, are not
listed). Because our general approach is to systematically survey
the initial conditions, we do not identify every simulation in Table 2
with an exact snapshot from Table 1; rather, they should be consid-
ered a systematic parameter survey of possible intermediate-scale
conditions, motivated by the typical range of sub-kpc conditions
seen in our galaxy-scale simulations. Dark matter is present, but is
1 As one continues the expansion to include arbitrarily small-scale modes,
the best-fit mode amplitude will eventually yield an amplitude consistent
with the shot noise in the simulation, roughly at the scale of the median
inter-particle spacing; we discard higher mode numbers as they are particle-
noise dominated.
2 To avoid numerical effects from a step-function cutoff in the mass profile
at small radii, we typically truncate the gas mass profile with a gas /
(R=R0)2 power law inside of a radius R0  35 in the parent simulation
( is the minimum smoothing length). Gas within this radius (< 1% of the
re-simulated mass) is initialized with circular orbits.
dynamically irrelevant at these scales. We have also varied the mass
of the central black hole, but at these scales it is still dynamically
unimportant. Although the initial conditions for our calculations
are drawn from galaxy-scale simulations, the dynamics on small-
scales depends primarily on a few key properties of the simulation
( fgas and B=T ), and is thus insensitive to many of the details of the
galaxy on larger scales.
Our intermediate scale simulations typically involve '
106 particles, with a force resolution of a few pc and a particle mass
of  104 M. The duration of the simulation is  10
7 108 yr –
this is many dynamical times at small radii, but small compared to
the dynamical time at larger radii.3 These re-simulations can thus
be thought of as a probe of the instantaneous behavior of the gas at
small radii given the inflow conditions set at larger radii.
We discuss a number of numerical tests and variations about
this basic methodology in Appendix A. Specifically, we show that
our results are not sensitive to including properties of the larger-
scale galaxy in which the simulation should be embedded, such
as the  kpc-scale tidal potential. They are also not sensitive to
whether we initialize an axisymmetric gas/potential distribution, or
whether the initial condition includes the non-axisymmetric modes
present in the larger-scale simulation. The reason is that instabili-
ties due to self-gravity can grow exponentially from shot-noise in
the simulation even given an initially axisymmetric structure. Thus
the presence of initial asymmetries on  100 pc scales does not
have a significant effect on the resulting transport of gas to smaller
radii; the transport is determined by the presence (or absence) of
internal instabilities in the gas on small scales. Finally, because the
“initial” densities on small scales are intentionally initialized to be
low relative to their values later in the simulation, our results are
not particularly sensitive to how we initialize gas at small radii in
the re-simulation; this is important because there is no reliable in-
formation about these small-scales in our larger-scale simulation.
To check that our re-simulation approach has not introduced
any artificial behavior, we have run a small number of ultra-high
resolution galaxy scale simulations, the inner properties of which
can be compared to the intermediate scale re-simulations summa-
rized here. For these very high resolution calculations there is con-
tinuous inflow from large scales, so they can be self-consistently
evolved for many dynamical times. The expense of these cal-
culations, however, limits the survey of initial conditions possi-
ble. We have 6 such simulations: three mergers, identical in mass
and geometry to our canonical case in Table 1, but with initial
fgas = 0:2; 0:4; 0:8, and  107 particles, which gives SPH smooth-
ing lengths of  10pc. While not quite as high-resolution as our
re-simulation runs, these provide an important check on the results
of the latter and are run self-consistently for 4 109 yr. We will
show results from these ultra-high resolution simulations at several
points; we find that they are quite similar to our re-simulation runs,
thus supporting the methodology used for most of our calculations.
The very high resolution merger calculations also allow us to
follow the binary BH pair to much smaller separation (our assump-
tions lead to rapid merger below  10pc in these calculations). We
confirm, in these cases, that the BH-BH merger precedes most of
the gas inflows at & 10 pc, so that the assumption that the binary
has merged is probably reasonable for our re-simulation calcula-
3 To ensure there are no later-time phenomena of interest, and to study the
relaxed structure of the stellar remnant produced by each re-simulation, we
evolve most for & 2108 years, by which time all the gas is exhausted. We
find that there is no qualitatively new behavior at these later times.
c
0000 RAS, MNRAS 000, 000–000

8 Hopkins and Quataert
Table 2. Intermediate-Scale Resimulations ( 101000 pc)
Simulation  qeos fgas h=R B=T (< R) d(0) b(0) Mtot(< R) [M]
Name [pc] (500 pc) (500 pc) 100 pc 300 pc [M kpc
2] 100 pc 300 pc
(1) (2) (3) (4) (5) (6) (7) (8) (9)
If9b5a 1.0 0b,0.175,0.25 0.90 0.30 0.5 0.15 1.0e10 1.1e10 5.4e8 2.9e9
If9b5thin 1.0 0.125,0.25 0.90 0.08,0.16 0.4 0.2 1.0e10 1.1e10 5.4e8 2.9e9
If9b5res 0.3,1,3,10 0.125 0.90 0.30 0.5 0.15 1.0e10 1.1e10 5.4e8 2.9e9
If9b5q 1.0 0,0c,0.125,0.25,0.5,1 0.90 0.30 0.5 0.15 1.0e10 1.1e10 5.4e8 2.9e9
Ilowresq 3.0 0,0.25,1 0.95 0.27 0.0 0.0 6.0e10 0.0 1.2e9 4.9e9
If1b1late 1.0 0.125 0.096 0.25 0.06 0.03 1.0e11 1.1e10 3.4e9 2.2e10
If1b0late 1.0 0.25 0.091 0.28 0.002 0.003 6.0e10 5.0e8 1.3e9 5.3e9
If1b0lateLd 2.0 0.25 0.091 0.28 0.005 0.01 6.0e10 5.0e8 7.4e9 8.5e9
If3b3mid 1.0 0.125 0.34 0.25 0.3 0.15 3.1e10 2.0e10 1.3e9 7.2e9
If3b3midRge 1.0 0.125 0.45,0.20,0.05 0.3,0.5 0.26 0.12 3.6e10 2.0e10 1.5e9 9.2e9
If1b3Lmid 1.0 0.25 0.10 0.30 0.07 0.10 6.4e10 1.6e9 1.4e9 5.7e9
If1b3LmidLd 2.0 0.25 0.10 0.30 0.15 0.25 6.4e10 1.6e9 8.4e9 1.1e10
If5b4mbul 1.0 0.25 0.50 0.24 0.40 0.55 7.4e8 1.9e8 0.3e8 1.3e8
If5b8mbul 1.0 0.25 0.50 0.24 0.80 0.90 7.4e8 3.1e9 1.2e8 5.3e8
If9b1lowm 1.0 0.25 0.90 0.15 0.03 0.06 6.4e9 1.2e7 1.3e8 5.4e8
If3b9dsk 1.0 0.25 0.32 0.22 0.95 0.80 1.6e9 1.1e11 2.1e9 6.2e9
If3b9dskLd 2.0 0.25 0.32 0.22 0.71 0.62 1.6e9 1.1e11 8.0e9 9.3e9
IfXb2gas 1.0 0.25 0.16,0.32,0.50 0.26 0.16 0.08 3.6e10 1.0e10 1.3e9 7.7e9
Inf28b2 1.0 0.20 0.20,0.80 0.25 0.51 0.28 6.3e9 3.0e8 3.6e8 1.7e9
Inf28b4 1.0 0.20 0.20,0.80 0.25 0.67 0.43 3.3e9 3.0e8 2.7e8 1.1e9
Inf28b6 1.0 0.20 0.20,0.80 0.25 0.78 0.56 3.3e9 6.0e8 4.0e8 1.4e9
Inf28b8 1.0 0.20 0.20.0.80 0.25 0.92 0.81 1.0e9 6.0e8 3.4e8 9.5e8
Inf2b9 1.0 0.20 0.20 0.25 0.96 0.89 5.2e8 6.0e8 3.2e8 8.6e8
Inf28b2hf 0.3 0.125 0.20,0.80 0.25 0.38 0.21 5.5e9 5.0e7 2.3e8 1.1e9
Inf8b2hrf 0.3 0.20 0.80 0.25 0.17 0.11 5.5e9 2.5e7 2.2e8 1.1e9
Inf28b9hf 0.3 0.125 0.20,0.80 0.25 0.81 0.75 5.5e9 5.0e10 2.1e9 1.1e10
Parameters describing our re-simulations of the 0:011 kpc regions from galaxy scale simulations. Parameters separated by commas denote
simulations with otherwise identical initial conditions, re-run with the specified parameter varied. (1) Simulation name/ID. (2) Minimum
smoothing length (in pc). (3) Equation of state parameter (Figure 1). (4) Initial gas fraction of the disky/cold component (inside the given radius).
(5) Initial scale height of the disk component (inside the given radius). (6) Initial bulge-to-total mass ratio inside a given radius (again, bulge
refers to all spherical components). (7) Initial maximum surface density of the disky/cold component (gas plus stars). (8) Initial average surface
density inside 300 pc for the bulge component (9) Initial total mass enclosed inside a given radius. All simulations include BHs and dark matter,
but these are dynamically unimportant on these scales.
a A series of 7 runs testing different means of constructing initial conditions, described in Appendix A.
b Isothermal equation of state, but with a large cs = 50kms1 cooling “floor.”
c Cooling allowed down to 100K, i.e. cs = 1kms1.
d Somewhat larger-scale simulation (between “galaxy scale” and standard “intermediate scale”). Instead of B=T (< R) and Mtot(< R) being
evaluated at 100 pc and 300pc, they are here evaluated at 500pc and 1kpc, respectively.
e Series where the gas disk profile is allowed to vary independent of the stellar disk profile. The gas has exponential, power-law, and truncated
power-law profiles, with varying concentrations with respect to the disk (for example including an extended gas “reservoir” at a distance  2
times the regular nuclear stellar disk length, with surface density profile / R1).
f Very high-resolution simulations which also act as a moderate-resolution nuclear-scale simulations (2107 particles; gas particle mass
 500M). A series of 6 galaxy-scale runs with very high ( 10pc) resolution, used as moderate-resolution intermediate-scale simulations, are
also described in the text.
tions. Moreover, the gas mass at  10 pc is large ( MBH) in the
merger simulations. Thus if gas-rich reservoirs indeed drive rapid
BH-BH coalescence, the rapid merger of the two BHs should be
a reasonable assumption on all of the scales that we simulate (see
e.g. Escala et al. 2004; Perets et al. 2007; Mayer et al. 2007; Perets
& Alexander 2008; Dotti et al. 2009; Cuadra et al. 2009).
2.4 Nuclear Scales: From 10 pc to 0.1 pc
The characteristic initial scale-lengths of the nuclear disks in our
intermediate scale calculations are  0:2 0:5kpc. As we discuss
in §3, if the gas fraction is sufficiently large, instabilities quickly
develop that transport material down to  110pc, near the reso-
lution limit of our intermediate scale calculations. Material begins
to pile up at these radii because the BH mass dominates the po-
tential and the efficiency of large-scale modes decreases at small
radii. In order to understand the dynamics on yet smaller scales, we
therefore repeat our “re-simulation” methodology once more. The
approach is identical to that described above, but this time using
the intermediate-scale simulations with resolution of. 10pc as our
“parent simulation” from which to motivate the initial conditions.
We again carried out  50 such simulations, typically with
 106 particles, and a force/spatial resolution of  0:1 pc (particle
mass  100M). The properties of these simulations are summa-
rized in Table 3. The simulations are evolved for  105 106 yr;
this is large compared to the dynamical time at the smallest radii
 0:1 pc, but very small relative to the dynamical time of the larger-
c
0000 RAS, MNRAS 000, 000–000

How Do Massive Black Holes Get Their Gas? 9
Table 3. Nuclear-Scale Resimulations ( 0:1100 pc)
Simulation  qeos fgas h=R MBH Mb=cl(< R) d(0) Md(< R) [M]
Name [pc] (100 pc) (100 pc) [M] 10 pc 50 pc [M kpc
2] 1 pc 10 pc 50 pc
(1) (2) (3) (4) (5) (6) (7) (8) (9)
Nf8h1c0 0.3 0.0,0.25 0.75 0.16,0.26 2.9e7 1.8e5 3.0e5 7.5e10 1.0e5 1.5e6 1.0e8
Nf8h1c0hola 0.3 0.25 0.75 0.24 2.9e7 1.8e5 3.0e5 7.4e10 1.0e3 1.2e6 1.0e8
Nf5h2c1 0.1 0.25 0.50 0.24 3.0e8 1.2e7 8.0e7 1.2e11 3.6e5 2.4e7 1.6e8
Nf28h1c1 0.1 0.25 0.19,0.75 0.27 3.0e7 1.3e7 2.8e7 7.5e10 3.6e5 2.4e7 1.6e8
Nf7h12c0dsk 0.1 0.25 0.70 0.23 2.9e7,3.0e8 0.0 0.0 4.5e9 1.0e4 0.9e6 6.0e6
Nf5h1c0 0.1 0.25 0.48 0.25 2.9e7 0.0 0.0 1.2e11 3.0e5 2.4e7 1.6e8
Nf5h1c2 0.1 0.25 0.48 0.25 2.9e7 3.0e7 9.0e7 1.2e11 3.0e5 2.4e7 1.6e8
Nf8h1c0thin 0.1 0.125 0.75 0.16,0.27 3.0e7 1.8e5 0.3e7 7.6e10 2.0e5 1.5e7 1.0e8
Nf5h1c1thin2 0.1 0.125,0.25 0.50 0.07,0.14 3.0e7 3.0e5 1.2e7 1.2e11 2.0e5 2.4e7 1.6e8
Nf8h1c1qs 0.1 0,0.125,0.25,1 0.75 0.28 3.0e7 3.0e6 1.4e7 7.3e10 2.0e5 1.7e7 1.2e8
Nf8h2c1 0.1 0.25 0.75 0.2,0.33 3.0e8 3.6e6 1.4e7 7.6e10 2.0e5 1.9e7 1.2e8
Nf1h1c1low 0.1 0.25 0.08 0.25 3.0e7 3.5e6 1.4e7 1.7e11 4.7e5 3.7e7 2.5e8
Nf3h1c1mid 0.1 0.20 0.26 0.23 3.0e7 3.5e6 1.4e7 2.1e11 6.0e5 4.6e7 3.0e8
Nf6h12c2dsk 0.1 0.25 0.57 0.22 2.9e7,3.0e8 7.2e6 2.6e7 4.4e9 1.1e5 8.1e6 3.4e7
Nf8h1c3dskM 0.1 0.25 0.75 0.30 2.9e7 1.6e7 7.1e7 7.4e10 3.5e5 3.0e7 1.7e8
Nf8h1c1dens 0.1 0.25 0.75 0.25 3.0e7 3.6e6 1.4e7 3.8e11 1.1e6 8.1e7 5.3e8
Nf8h1c1ICsb 0.1 0.25 0.75 0.28 3.0e7 3.1e6 1.4e7 7.3e10 2.0e5 1.7e7 1.2e8
Nf8h1c1thin 0.1 0.125,0.18 0.75 0.08,0.17 3.0e7 3.1e6 1.4e7 7.3e10 2.0e5 1.7e7 1.2e8
Nf2h2b2 0.1 0.20 0.20 0.25 3.0e7 6.5e6 2.0e7 3.0e11 8.1e4 7.0e7 4.2e8
Nf8h2b2 0.1 0.20 0.80 0.25 3.0e7 1.3e7 4.0e7 1.5e11 4.7e4 3.5e7 2.1e8
Nf2h2b4 0.1 0.20 0.20 0.25 3.0e7 6.5e6 2.1e7 1.5e11 4.3e4 3.5e7 2.1e8
Nf8h2b4 0.1 0.20 0.80 0.25 3.0e7 1.3e7 4.0e7 7.7e10 2.4e4 1.7e7 1.1e8
Nf2h2b5 0.1 0.20 0.20 0.25 3.0e7 6.5e6 2.1e7 7.7e10 2.1e4 1.8e7 1.1e8
Nf28h2b6 0.1 0.20 0.20,0.80 0.25 3.0e7 1.1e7 3.7e7 3.8e10 1.1e4 8.8e6 5.3e7
Nf8h2b8 0.1 0.20 0.80 0.25 3.0e7 1.3e7 4.1e7 1.5e10 4.7e3 3.5e6 2.1e7
Nf28h2b9 0.1 0.20 0.20,0.80 0.25 3.0e7 1.1e7 3.8e7 9.6e9 2.7e3 2.2e6 1.3e7
Nf8h2b1h 0.015 0.20 0.80 0.25 3.0e7 6.4e6 2.0e7 1.5e11 5.2e4 3.5e7 2.1e8
Nf8h2b3L 0.1 0.20 0.80 0.25 3.0e7 2.3e7 1.1e8 1.5e11 9.3e3 3.9e7 4.9e8
Nf8h2b4q 0.1 0,0.02,0.06,0.12, 0.80 0.25 3.0e7 1.3e7 4.0e7 7.7e10 2.4e4 1.7e7 1.1e8
0.25,0.35,0.5,0.7,1
Parameters describing our nuclear-scale re-simulations of the sub-100 pc regions from intermediate-scale simulations. (1) Simulation name/ID. (2)
Minimum smoothing length (in pc). (3) Equation of state parameter (Figure 1). (4) Initial gas fraction of the disky/cold component. (5) Initial scale
height of the disky component. (6) Black hole mass (M). (7) Initial bulge or nuclear stellar cluster mass, inside the given radius. (8) Initial
maximum surface density of the disky/cold component (gas plus stars). (9) Initial mass of the disky component (gas plus stars) inside a given radius
(does not include the BH mass or, if significant, nuclear star cluster/bulge mass). Dark matter is insignificant on these scales.
a Central “hole” is extended in disk out to 5pc.
b Simulations with no central deficit of matter; the initial density from the larger-scale simulation is extrapolated in to r! 0. Also expanded into series
of initial conditions, as described in Appendix A.
scale simulations from which the initial conditions are drawn. The
characteristic spatial scale of the re-simulated material is initially
 1030pc. As described in Appendix A, we carried out a num-
ber of numerical tests of the robustness of these simulations.
At radii  pc, the parameters that determine the dynamics are
largely the equation of state of the gas, the mass of the BH, the
mass of the nuclear disk formed by the inflow from larger scales,
and the gas fraction of that nuclear disk. Since the BH dominates
the spherical component of the potential at these radii, the “bulge”
mass at these radii is only of secondary importance; we include it
but find that it makes little difference.
As in §2.3, we have checked the results of these “re-
simulations” by carrying out a small subset of ultra-high resolu-
tion runs. These extend from  0:3 1000 pc and follow inflow
from larger scales deep into the potential of the BH; because they
resolve larger spatial scales than our typical "nuclear scale" simula-
tion, these can be run self-consistently for 2 108 yr. Specifically,
we have five such high-resolution intermediate scale simulations
(see Table 2), three with initial fgas = 0:8 (a low, intermediate, and
high B=T case), and two with fgas = 0:2 (low and high B=T ). They
have 107 particles and gravitational softening lengths of 0:3pc.
We show the results from these runs explicitly at several points;
we find that they are completely consistent with our survey of re-
simulations, which cover a larger parameter space of galaxy/BH
properties, but are more limited in dynamic range.
It is important to note up-front that our simplified treatment of
the ISM physics becomes particularly suspect on nuclear scales 
pc. At these radii, our assumption that we can average over the dy-
namics of stellar winds, supernovae, HII regions, etc. and define an
effective ISM equation of state may break down. Nonetheless, we
believe that the efficient angular momentum transport found here is
likely generic, so long as some of the gas is prevented from forming
stars and the gas fraction is sufficiently high that instabilities gen-
erated by self-gravity are initiated. The fact that the main sequence
lifetime of a massive star is longer than the local dynamical time
on small scales probably increases the efficacy of stellar feedback
c
0000 RAS, MNRAS 000, 000–000

10 Hopkins and Quataert
and decreases the fraction of the gas turned into stars per dynamical
time (Murray et al. 2009).
At scales 0:1pc, the potential is fully Keplerian and viscous
heating is sufficient to stabilize the disk against its own self-gravity
(i.e., Q & 1) (Goodman 2003). At these radii, the system begins to
approach a traditional accretion disk. Given the cessation of star
formation and the deep potential well of the BH, we assume that
the inflow rate at  0:1 pc is a reasonable proxy for the true accre-
tion rate onto the BH. Because our simulations are not well-suited
to describe the physics of the disk on scales . 0:1 pc, we do not
perform a further “zoom in.”
3 OVERVIEW: FROM KPC TO SUB-PC SCALES
Using the numerical simulations described in §2, we now describe
how gas is transported fromkpc scales topc scales. Initially,
our discussion is somewhat qualitative; we focus on emphasizing
the key physics at play and our key results. In § 5 we discuss the
relevant stability criteria more quantitatively, and outline some spe-
cific criteria necessary for “interesting” gas inflow.
Figure 2 shows an illustrative example of the results of our
re-simulations on various scales. We plot gas surface density maps,
with color encoding the gas effective sound speed, from scales of
 100kpc to < 1pc. The simulation in this case is a fairly gas-rich
major merger ( fgas  3040% at the time of the final coalescence)
of two 5 1010 M baryonic mass galaxies. The smaller scale re-
simulations were carried out just after the coalescence of the two
nuclei, which is near the peak of star formation activity, but when
the system is still quite gas rich. The initial systems had pre-existing
bulges of  1=3 the disk mass and BHs initialized on the MBH
relation ( 107 M). Each panel is rotated so that the view is close
to “face-on” with respect to the total angular momentum vector
of the gas plotted in the image. Viewed edge-on, much of the gas
forms a modestly thin (H=R . 0:3) disky distribution at all radii.
3.1 Large Scales: Mergers and Bars
From  100kpc to & 100pc, the gas flows are well-resolved by
our galaxy-scale simulation. In mergers, the final collision of the
two galaxies yields strong torques that efficiently cause most of the
gas to flow to the center on a timescale approximately equal to a
few dynamical times. This process has been described in detail in
e.g. Hopkins et al. (2009b), and references therein, but for com-
pleteness we briefly summarize the important physics. The sec-
ondary/merging galaxy does not directly torque the gas. Rather,
the torques on the gas are dominated by local torques from stars
originally in the same disk as the gas. The merger induces triaxial
“sloshes” and bar-like structures in the stars, i.e., non-axisymmetric
modes, supported by radial and/or random orbits. These are also
induced in the gas. However, because the gas is dissipational, the
gaseous modes slightly lead those in the stars. The stellar distur-
bance, being physically close to, trailing, and in near-resonance
with the gas, produces a strong torque that removes angular mo-
mentum from the gas and easily dominates the total torque (torques
directly from the secondary galaxy, from the primary halo, or hy-
drodynamic torques from shocks or internal clump collisions, are
all . 10% effects; see e.g. Barnes & Hernquist 1996; Barnes 1998;
Hopkins et al. 2009b).
After the two galactic nuclei coalesce, the disturbances in the
stellar component of the galaxy relax away in a number of crossing
times. Until they relax, gas inflows continue. Moreover, the coales-
cence of the nuclei completes much more rapidly, in a timescale
close to a single crossing time, at least in a major merger. Thus, a
significant fraction of the gas inflows can occur in the background
of a rapidly relaxing stellar potential, in the wake of the nuclear
coalescence. This is the stage illustrated in Figure 2.
The gas that loses angular momentum flows in to radii 
0:51kpc (for an  L system), where it participates in a nuclear
starburst and builds a dense central stellar mass concentration, criti-
cal for establishing the structural properties and size of the remnant
spheroid. At these scales, the system is often gas-dominated for a
short period of time owing to these inflows (provided the merger
is sufficiently gas-rich). However, as the gas forms stars, the cen-
tral region will quickly become more stellar-dominated; because
these stars form out of the gaseous disk, in the relaxing potential
— they are not themselves violently relaxed. This is important for
the subsequent evolution of the system because of the presence of
disk instabilities that would be suppressed by a larger dispersion-
supported (spherical) component in the very central region.
The general scenario summarized here can be applied not just
to major mergers, but also to minor mergers, fly-by encounters, and
even sufficiently bar-unstable stellar disks. The details will be dif-
ferent, but the qualitative steps above, and the exchange of angular
momentum between gas and stars, is robust, ultimately leading to
inflow to sub-kpc scales. The subsequent evolution depends largely
how much material is efficiently channeled to small radii (relative
to the bulge and BH mass), not on how that material gets there.
3.2 Morphology and Gas Transport From 1 kpc to 10 pc
3.2.1 General Behavior
The gas infalling from large radii begins to “pile up” at radii
 0:1 1 kpc, rather than continuously flowing in to yet smaller
radii. This is because the torques from the disturbances at large
radii become less efficient at small radii. This happens for three
reasons: (1) In the merger context, the stellar perturbation at small
radii relaxes after coalescence, decreasing the efficiency of gas in-
flow. (2) The rapid gas inflow implies that the system becomes in-
creasingly gas-dominated at radii  100 pc, even if the initial disk
gas fraction is low, 0:1. Because the primary angular momentum
sink of the gas is the local stars, when the system becomes locally
gas-dominated, angular momentum transfer is actually less efficient
(see Hopkins et al. 2009b,f). (3) The gas can encounter the equiv-
alent of an inner Linblad resonance. This is especially important
for unstable gas bars, minor mergers, and disturbances induced by
early passages. For the case of coalescence following major merg-
ers, the disturbance is not a single mode, but a series of modes at all
scales. As such, there is often no formal inner Linblad resonance or
angular momentum barrier (each mode may have such a barrier, but
these are spread over a wide range of scales; there is thus a means
to overcome the barrier associated with any single mode).
Figure 2 shows the outcome of gas pile-up in the central
kpc using an intermediate-scale re-simulation (middle row). In this
case, the intermediate scale simulation is a high-resolution re-
simulation of the larger-scale gas distribution at a given epoch in
a gas-rich major merger. The gas density reached from the larger-
scale inflows is quite large –  1010 M worth of gas has formed
a disky component with a scale length of 0:3kpc and an average
surface density of  1010 M kpc
2. This is a large fraction of the
galaxy mass – larger than the pre-existing bulge within these radii.
The small-scale gas disk is therefore strongly self-gravitating. In-
deed, we see from Figure 2 that it quickly develops unstable, non-
axisymmetric modes.4 This is essentially the “bars within bars”
4 These modes develop almost identically even if we initialize the re-
simulation to be perfectly smooth and remove all external tidal forces (see
c
0000 RAS, MNRAS 000, 000–000

How Do Massive Black Holes Get Their Gas? 11








30 kpc
Galaxy-Scale Simulation:








7.2 kpc










3.6 kpc








600 pc
Intermediate-Scale Re-Simulation:










180 pc








60 pc








15 pc
Nuclear-Scale Re-Simulation:








6 pc










2 pc
Effective cs [km s-1]
10 200
Figure 2. Example of our multi-scale simulations used to follow gas flows from  100kpc to  0:1pc. Each row is a separate simulation, with the initial
conditions for the intermediate and nuclear-scale simulations taken from the output of the larger-scale runs in the row above it. Each panel shows the projected
gas density (intensity) and effective sound speed (color; blue is gas with an effective cs  10kms1, through yellow at  100 200kms1). Each image
is rotated to project the gas density “face on” relative to its angular momentum vector. From top left to bottom right, panels zoom in to the nuclear region
around the BH, with resolution spanning a factor 106 in radius. Top: Large-scale gas-rich galaxy-galaxy major merger simulation, just after the coalescence
of the two nuclei (run b3ex(co) in Table 1). The apparent second nucleus is actually a clump formed from gravitational instability. Middle: A higher-resolution
re-simulation of the conditions in the central kpc (run If3b3midRg in Table 2). Despite the fact that the background potential is largely relaxed on these scales,
the very large gas inflows lead to a strongly self-gravitating disk on  0:5kpc scales that develops a strong spiral instability, leading to efficient angular
momentum transport to  10pc. Again, some clumping appears (there is only one nucleus). Bottom: High resolution re-simulation of the central  30 pc of
the intermediate-scale simulation, with a resolution  0:1pc (run Nf5h1c2 in Table 3). The potential is quasi-Keplerian, suppressing traditional bar and spiral
instabilities, but the large inflows lead to a self-gravitating system that develops a standing eccentric disk mode (single-armed m = 1). The stellar and gaseous
eccentric disks precess relative to one another on  110pc scales and drive efficient inflows of  10M yr
1 into the central 0:1pc.
c
0000 RAS, MNRAS 000, 000–000

12 Hopkins and Quataert
Figure 3. Images of the morphologies of the secondary instabilities that develop on intermediate scales ( 100 pc), as a consequence of galactic-scale torques
bringing gas into the central kpc. Each panel shows projected gas density (intensity) and local star formation rate (color, from blue through yellow) in the
central kpc. We sort the simulations by the parameters that have the largest impact on the dynamics: the bulge-to-disk ratio (disk being both stellar and
gaseous) within 300pc and gas fraction Mgas=(Mgas + M; disk): from top – bulge-dominated systems with B=T & 0:8 to bottom – disk-dominated systems
with B=T . 0:1; and from left – systems with fgas . 0:1 to right – systems with fgas & 0:8. Global instabilities become more prominent with decreasing B=T ;
local instabilities (clumping/fragmentation and/or more tightly wound spirals) become more prominent with increasing fgas. The simulations at fixed B/T
and/or fgas are similar, but include variation in e.g. the initial mass profile shapes; this contributes to the dramatic diversity of morphologies on nuclear scales.
Almost all systems that are not entirely bulge-dominated develop strong large-scale instabilities that efficiently transport gas to 10pc. Observationally, these
would be categorized as a mix of bars, spirals, nuclear rings, crossed rings, fragmentation/clump instabilities, and single or triple-armed spirals.
c
0000 RAS, MNRAS 000, 000–000

How Do Massive Black Holes Get Their Gas? 13
scenario predicted by Shlosman et al. (1989), although the mor-
phology of the system is clearly not a simple bar; this is an impor-
tant point to which we return below. The strength of the modes that
develop in this re-simulation depend on the fact that there is star
formation in the gas – as the gas turns into stars in situ, those stars
develop non-axisymmetric modes, and the two precess relative to
one another. As in the galactic-scale torques discussed above, this
produces particularly efficient angular momentum transport. These
processes ultimately lead to a significant amount of gas flowing
down to  10pc.
3.2.2 Diversity in Morphologies and Inflow Strengths
For our fiducial parameterization of the ISM equation of state, the
disk-to-bulge ratio and gas fraction on ' 0:10:3 kpc scales have
the largest influence on the dynamics and angular momentum trans-
port in our intermediate-scale simulations (for discussion of the role
of qEOS and sub-resolution physics, see Appendix B). Figure 3 il-
lustrates this by showing the 100 pc scale morphology of a repre-
sentative subset of our intermediate scale re-simulations, each after
a couple of local dynamical times of evolution ( 107 108 yr).
These are sorted by gas fraction fg and B=T .5
As Figure 3 shows, systems with very large B=T & 0:9 (top
row) are globally stable, as expected analytically. In the extremely
gas-rich, large B=T , case (top right), local Toomre-scale instabili-
ties develop – if such small clumps were infinitely long-lived, their
orbits would decay via dynamical friction and allow for some trans-
port of gas to the center. However, because the clumps are dense,
they quickly turn into stars – such a mechanism is largely a means
to move stars to the galactic nucleus, not gas (and in fact leads to
nuclear stellar clusters much larger than those observed; see Ap-
pendix B for a more detailed discussion). Previous claims that such
clump-sinking could efficiently fuel BH growth (see e.g. the dis-
cussion in Wada et al. 2009; Kawakatu & Wada 2008, and refer-
ences therein) have neglected star formation and have thus dra-
matically over-estimated the inflow of gas via this process. Lo-
cal clumping instabilities like these do, of course, occur, forming
molecular clouds and star clusters. However, there are both theo-
retical and observational arguments that suggest that they quickly
dissolve on a few local dynamical times, probably via feedback
from some combination of stellar winds, HII regions, and radiation
pressure (Larson 1981; Blitz 1993; Krumholz et al. 2006; Allen
et al. 2007; Bonnell et al. 2006; Elmegreen 2007). We also note
that, under certain conditions, such local instabilities could instead
produce small scale, tightly wound spiral waves instead of leading
to fragmentation (e.g. Lodato & Rice 2004; Rice et al. 2005; Boley
et al. 2006); however, this generally occurs when the cooling time
is comparable to or larger than the dynamical time, which is not the
case in these simulations (although it will also depend on the tur-
bulent decay time which, if sufficiently large, can allow turbulence
to suppress runaway clumping and star formation).
Figure 3 demonstrates that the strength of global non-
axisymmetric modes increases dramatically with decreasing B=T .
In fact, as soon as the bulge and disk are comparable (second row
from top), prominent modes appear. In the disk-dominated cases
(bottom row), very large angular momentum transport occurs even
Appendix A). It is thus not sensitive to the larger-scale environment; rather,
the system is simply strongly globally unstable.
5 The other parameters of the simulations are not all identical, representing
the properties of the simulations from which they are selected, but the key
qualitative behavior depends primarily on these two parameters.
in cases without much total gas. Of course, at each B=T , increasing
the gas content makes the system more vulnerable to local instabil-
ities as well – usually making the overall inflow more clumpy and
time-variable. In several cases, the large-scale mode (e.g. a spiral
arm) becomes self-gravitating and globally fragments (not neces-
sarily into smaller sub-units, but as a whole), leading to a major co-
alescence – what is almost (dynamically speaking) a scaled-down
merger in the central regions!
Figure 3 also demonstrates an important point seen in Fig-
ure 2. Although the instabilities seen in our simulations qualita-
tively resemble the “bars within bars” idea of secondary instabili-
ties once the gas density is sufficiently high, the morphologies vary
widely, and are not restricted to traditional bars (although these cer-
tainly do appear). At similar gas fraction and B=T , we find that the
strength of angular momentum transfer is generally similar, but we
also find that the visual morphology (and the precise modes im-
portant for transport) can vary widely, depending on time and on
the details of the gas, stellar disk, bulge, and halo profiles, and the
precise equation of state. Thus, global quantities such as the mass
profile and accretion rate are comparatively robust, but the obser-
vational classification of these systems would vary widely.6 Indeed,
Figure 3 shows traditional bars and spirals, nuclear rings, crossed
or barred rings, single or three-armed systems, flocculent disks, and
clumpy, irregular morphologies. These all appear, with no obvious
preference for one or another as a whole, in our simulations.
This feature of our simulations may account for a number of
observational results in the literature. For example, surveys of AGN
have often found that although nuclear bars on these scales only ap-
pear in some fraction of sources (not necessarily much larger than
the fraction of non-AGN in which they appear), there are ubiquitous
asymmetric gas structures of some sort, similar to those modeled
here (see e.g. Martini & Pogge 1999; García-Burillo et al. 2005;
Haan et al. 2009; Krips et al. 2007; Laine et al. 2002; Peletier et al.
1999; Sakamoto et al. 1999). We discuss this further in § 7.
3.3 Towards the Accretion Disk: Eccentric Disks at  Parsec
From 500 to 10pc, our intermediate-scale simulations success-
fully demonstrate efficient angular momentum transport via a wide
range of gravitationally unstable modes. Near the smallest radii in
these simulations, however, the systems encounter yet another an-
gular momentum barrier. At that point, gas has reached the BH ra-
dius of influence, i.e., the BH begins to contribute non-trivially to
the potential, which becomes quasi-Keplerian. This halts further
inflow because the disk is no longer strongly self-gravitating and
so is less susceptible to global modes. In addition, the gas gen-
erally encounters an inner Linblad resonance associated with the
intermediate-scale bar.
Indeed, it is widely appreciated that both “bars within bars”
and the direct or induced torques due to perturbations from merg-
ers, close passages, and large-scale bars do not produce efficient an-
gular momentum transfer at radii 0:110pc (see e.g. the discus-
sion in Athanassoula et al. 1983, 2005; Shlosman et al. 1989; Heller
et al. 2001; Begelman & Shlosman 2009, and references therein).
6 In detail, the angular momentum transport does depend on the structure of
the unstable mode/perturbation. The precise dependence will be discussed
in Paper II. At the dimensional level, the torques and inflow scale in the
same manner independent of the detailed mode morphology, so long as the
stellar torques are sufficiently strong to cause orbit crossing and shocks in
the gas. The details of the specific modes driving such shocks amounts to
numerical factors of a few in the torques.
c
0000 RAS, MNRAS 000, 000–000

14 Hopkins and Quataert
Figure 4. Images of the instabilities that develop in our small-scale ( 10 pc) nuclear re-simulations with 0:1pc resolution. Each panel shows projected gas
density (intensity) and local star formation rate (color, from blue through yellow) in the central kpc. Each of these simulations can be thought of as a simulation
of the corresponding nuclear scales from Figure 3. The simulations extend into the BH radius of influence. The primary parameters of importance are the ratio
of the gas mass to the BH mass (or BH plus bulge/star cluster mass, when the latter is present) and the gas fraction in the disky component: top to bottom
is decreasing BH/stellar mass while the disk gas fraction increases from left to right. A strong m = 1 mode is generic for reasonable BH/stellar mass and
gas fractions – this corresponds to an eccentric, globally precessing (non-winding) disk (or single-armed spiral), a mode that is special to the quasi-Keplerian
potential. The resulting torques drive inflows of up to 10M yr
1 at < 0:1pc scales (Figure 5), sufficient to fuel a luminous quasar.
As a consequence this is often considered the most difficult-to-
explain regime of gas inflow and angular momentum transport.
We find, however, that efficient angular momentum transfer
continues in our smaller-scale re-simulations for sufficiently gas-
rich, disk-dominated systems. Figure 4 shows this with images
spanning a range of fg and B=T ; we will discuss this physics in
more detail in §4 & 5. These results show that, so long as the gas
density is sufficiently large (relative to the BH mass), the system
develops a precessing eccentric disk (an m = 1 mode) that drives
gas down to sub-pc scales  0:1 pc. As before, stars rapidly form
out of the disk, leading to a similar mode in both the stars and
gas; these modes precess about the BH relative to one another with
c
0000 RAS, MNRAS 000, 000–000

How Do Massive Black Holes Get Their Gas? 15
0.0 0.5 1.0 1.5 2.0 2.5 3.0
t [Gyr]
0
100
200
300
G
as
In
flo
w
Ra
te
d
M
/d
t
[M
O •
yr
-
1 ]
R < 300 pc
0 2 4 6 8 10 12 14
t [Myr]
-50
0
50
100
150
200
250
300

R < 10 pc
0.0 0.2 0.4 0.6 0.8 1.0 1.2
t [Myr]
-2
0
2
4
6
8
10
12

R < 0.1 pc
Figure 5. Time-dependent inflow rate through various radii, for the example simulation shown in Figure 2. Main: Gas inflow rate into the central 300pc in
the large-scale galaxy merger simulation, as a function of time since the beginning of the merger. The large spike follows the coalescence of the two nuclei.
Dotted line shows the instantaneous inflow rates; solid line shows the inflow rate averaged over five local dynamical times. Top Right: Inflow rates into the
central  10 pc in the intermediate scale re-simulation for a small time interval just after coalescence. Top Left: Inflow rates into 0:1 pc in our smallest-scale
resimulation, initialized with conditions from the previous intermediate-scale simulation after it has evolved for several local dynamical times. The accretion
rate reaches 110M yr1, sufficient to fuel a luminous quasar. Note the episodic, highly variable nature of accretion/inflow at each scale.
slightly different pattern speeds, leading to crossing orbits, dissi-
pation of energy and angular momentum in the gas, and thus net
inflow. This is, once again, an instability that depends primarily on
the presence of sufficient gas at small radii in the first place. We
find that this condition is met at some point(s) in time in all of our
simulations with significant gas mass and instability at  100pc,
i.e., in those simulations that meet the “bars within bars” criteria in
§3.2 and Figure 3.
4 INFLOW RATES AND GAS PROPERTIES
Figure 5 shows the time-dependent inflow rates through several an-
nuli, for the same set of multi-scale re-simulations shown in Fig-
ure 2. These demonstrate and quantify the general scenario sum-
marized above: on large scales, coalescence and the final stages of
the merger drive a large quantity of gas into the central few hun-
dred pc, with inflow rates  100 300M yr
1. The total dura-
tion of this phase is a few times 108 yr. During this time, the gas
accumulates at small radii  100 pc; simulating the dynamics on
this scale for  107 yr, we find that secondary gravitational insta-
bilities develop that drive further inflows into the central  10pc,
with inflow rates  10 200M yr
1. Zooming in yet again dur-
ing one epoch of significant inflow to . 10 pc, our smallest-scale
simulations resolve the rapid formation of an eccentric nuclear disk
around the central BH; the accretion rate into the central 0:1 pc,
which is likely a reasonable proxy for the accretion rate onto the
BH, reaches  1 10M yr
1. This is sufficient to fuel a lumi-
nous AGN at the Eddington rate.
Figure 5 also demonstrates that the small-scale accretion rate
can be highly time-variable. This is in part a consequence of the
accretion of individual clumps/clouds (e.g., Fig. 4), but is also a
consequence of the fact that gravitationally unstable perturbations
rapidly grow, dissipate, and generate other structures; depending
on, e.g., the state of precession of the stellar versus gaseous disk,
the system can transition between inflow and outflow at a given
radius. Even on the largest scales, the inflow is still highly vari-
able, although is coherent over a time much longer than the dy-
namical time because it is driven by the global torques involved in
the merger. Because of the variability in _M on different scales, we
do not expect every merger (or isolated galaxy with a large-scale
bar), at every time, to exhibit significant inflow from large scales
all the way down to the BH. This is important – after all, a large
fraction of observed mergers are not bright quasars. In addition, the
large variation in the physical conditions at small radii for a rela-
tively fixed set of conditions at larger radii demonstrates that great
care must be taken when trying to correlate the galactic structure at
c
0000 RAS, MNRAS 000, 000–000

How Do Massive Black Holes Get Their Gas? 17
0.1 1 10 100R [pc]
-400
-200
0
200
〈τ(R)
〉 [km
2 s-
2 ]
Gravitational:
from stars
from gas
large-scale tidesHydrodynamic
Figure 7. Radial profile of the specific torque in a representative nuclear
scale simulation (Nf8h1c1qs in Table 3) during a stage of rapid BH ac-
cretion. The torque per unit mass acting on the gas is averaged over narrow
radial annuli at one time; negative torques decrease the gas angular momen-
tum. The total gravitational torque (black line) is divided into the torque
from stars in the same disk as the gas (purple), gas (red), and the contri-
bution from the large-scale tidal field (orange). The hydrodynamic torques
from pressure forces and artificial viscosity (green) are comparatively small.
Inflow at all radii is dominated by gravitational torques, largely from nearby
stars in the same m = 1 pattern as the gas.
 , where the z axis is defined to be the angular momentum vector
of the disk (torques in the other directions are negligible, so this is
almost identical to plotting j j). The torques from pressure forces
are similarly defined by
P; j 
1
Mgas[
R]
Z R0+
R
R0
gas(r) [raP]z dr dcosd
=
1
Mgas[
R]
X
i in
R
mi [ri
1
(ri)
rP(ri)]z (3)
where the pressures and densities are determined in standard fash-
ion from the SPH quantities. Note that the BH can and does move
(with typical amplitude R  0:1R[Md(< R)  0:1MBH]) in re-
sponse to the m = 1 modes on nuclear scales; but since we are in-
terested in the inflow onto the BH itself, we evaluate these torques
about its instantaneous position (rather than, say, the center-of-mass
of the system). The results are qualitatively similar in Figure 7 if we
fix the central position (albeit quantitatively altered within 10pc
at a factor  2 level), but are less directly relevant to interpreting
inflow onto the BH.
The qualitative behavior shown in Figure 7 is representa-
tive of all of our simulations. At essentially all radii, gravitational
torques dominate hydrodynamic torques. Moreover, the gravita-
tional torques themselves are dominated by torques from stars, not
the torques of the gas on itself; specifically, the stars that are impor-
tant are in the same asymmetric perturbation as the gas. This is also
the case on larger scales (& 100pc), for both mergers and barred
systems (see Barnes & Hernquist 1996; Hopkins et al. 2009b).
Torques from the spherical component of the system (e.g. the halo
and/or bulge stars) are negligible at all radii. Unsurprisingly, the
torques from the large-scale tidal field, defined here as the torques
from the 100pc scales of the parent simulation from which the
initial conditions of this re-simulation were drawn, become signifi-
cant only at the outer boundary of our re-simulation (see Appendix
A2). Figure 7 shows that there are several sign changes in the torque
profile, reflecting the specific state of the system at this time; the
torque is time variable, but we find there are sign changes as well
in a time-averaged sense. Overall, the net rate of change of the an-
gular momentum of the gas very closely tracks the time-averaged
gravitational torque from the stars. Hydrodynamic torques never
induce very strong torques (greater than those shown), whereas the
stellar gravitational torques can be a factor of  10 100 larger
than in Figure 7 at some times.
The fact that the torques on the gas are dominated by stars
is robust, and occurs for two reasons. First, the stars contribute
significantly to the mass on the scales & 0:1 pc that we focus on
here (where star formation can occur). For typical star formation
efficiencies, it is difficult to have the gas mass  50% of the to-
tal mass for a reasonable fraction of the lifetime of the system. On
large scales, galaxies are known to not be so gas rich (although they
certainly can reach 50% gas fraction). On small scales, if the star
formation efficiency is a few percent per dynamical time, then even
a pure gas inflow from larger scales is likely to only remain gas-
dominated for  10 local dynamical times (on the smallest scales,
 106 yr) – thus for the majority of the time during which inflow is
continued, the system will contain a significant stellar mass. This is
consistent with direct observations of sub-kpc regions of starburst
galaxies (Downes & Solomon 1998; Bryant & Scoville 1999) and
the  1 10 pc nuclear scales around AGN (Hicks et al. 2009).
On sufficiently small scales, . 0:01 0:1 pc, star formation will
become inefficient, but at precisely those radii, we have (by defini-
tion) essentially reached the -disk.
Even if the gas mass is comparable to or greater than the stellar
mass, the self-torque of the gas on itself is much weaker than the
torques from the stars on the gas. This has been demonstrated in
detail in the case of large-scale perturbations from galaxy mergers,
bars, and spiral waves (Noguchi 1988; Barnes & Hernquist 1996;
Barnes 1998; Berentzen et al. 2007; Hopkins et al. 2009b). We
show that it is also true for smaller radii in Figure 7. We discuss
this in detail in a subsequent paper (Paper II; in preparation), but
briefly outline the key physics here.
It is well-known that the magnitude of the self-torque in a pure
gaseous or pure stellar disk is second-order in the mode amplitude
(/ jaj2, where jaj is the fractional magnitude of the asymmetry;
see e.g. Lynden-Bell & Kalnajs 1972). Moreover, Kalnajs (1971)
showed (with the exact solution for a logarithmic spiral) that the
torques are further suppressed by powers of  kR and m, where
k (m) is the radial (azimuthal) wavenumber of the mode. This is
despite the fact that the instantaneous torque on gas/stars moving
through a perturbation has no such suppression and is linear in jaj.
For small mode amplitudes, the gas orbits periodically (in the test-
particle limit) in response to the perturbation and the positive and
negative linear terms exactly cancel because of the symmetry on
either “side” of the mode. In a mixed stellar+gas system, however,
the two perturbations are not exactly in phase because of dissipation
in the gas. Perhaps most importantly, the gas can be driven into
shocks by a sufficiently strong perturbation in the stars, leading to
strong dissipation. As a result, the symmetric response of the gas
is broken, and the net torque from stars on the gas is linear in jaj
and largest in the global limit. At the order of magnitude level, the
inflow rate is thus given by
_Mgrav  jajgas R
2
: (4)
Note that expressed in terms of a viscosity  such that _M '
3gas, equation (4) corresponds to  ' (jaj=3)VcR, where Vc =
R
is the local circular velocity. Equivalently, the inflow time of
the gas is jaj1
1. Typical values of jaj are 0:010:3 in our
simulations when significant inflow is present (see §5). We stress
c
0000 RAS, MNRAS 000, 000–000

18 Hopkins and Quataert
0.1 1


8
9
10
11
12
log
(Σ)
[M O
•
kpc
-
2 ]
Nuclear-ScaleRe-Simulation
10 100








Intermediate-ScaleRe-Simulation
GasStars
1000








Galaxy-ScaleSimulation
0.1 1

1
10
100
SF
R(<
R)
[M O •

yr-1
]
10 100





1000





0.1 1

1
10
100
dM
/dt
[M O •

yr-1
] |a|ΣgasR2Ω
10 100





1000





0.1 1
0.01
0.10
1.00
z 0

/ R
10 100




1000




0.1 1R [pc]
0.1
1
10
Q ga
s(R)
10 100R [pc]




1000R [pc]




Figure 8. Properties of representative galactic (right panel), intermediate
(middle), and nuclear-scale (left) simulations as a function of radius at a
single time. The simulations are the same as in Figures 2 & 5 and the time
chosen corresponds to near the peak of activity, with significant gas inflow
to small radii. Note that the initial/boundary conditions for the smaller-scale
simulations were taken from the larger-scale simulations shown, so the re-
sults are reasonably continuous at the radial boundaries (despite the fact
that these do not represent a self-consistent single simulation). Top: Gas
and stellar surface density. Second-from-Top: Enclosed star formation rate
in a given radius. Middle: Instantaneous inflow/outflow rate through the an-
nulus (black line); note that there are several sign changes as a function of
radius. Also shown is the order of magnitude estimate of the accretion rate
produced by gravitational torques from a non-axisymmetric mode of am-
plitude jaj, _M  jajgas R2
(dashed red line; § 4). Second-from-Bottom:
Scale height z0=R of the gas. Bottom: Toomre Q parameter of the gas.
that, as described above, equation 4 differs fundamentally from the
dimensional expectation for self-torques, which is second-order in
jaj and suppressed by terms in jkRj and m.
Figure 8 shows a number of properties of our fiducial sim-
ulations from Figures 2 & 5, for each of the three spatial scales
we consider: galactic ( kpc), intermediate ( 10 100 pc), and
0.1 1 10 100 1000R [pc]

0.1
1
10
100
[dM/
dt] /
[ |a| Σ
gas
R2 Ω

]
Figure 9. Comparison of the true instantaneous inflow rate in the simula-
tions (dM=dt) to the simple dimensional scaling jajgas R2
(R) expected
for pure gravitational torques, where jaj is the magnitude of the strongest
non-axisymmetric mode at a given radius R (directly measured in the simu-
lations; with typical values 0:010:3). Each solid line denotes a different
simulation, with dotted lines showing the radii near the boundaries of our
re-simulations, where the exact profile shape is suspect; colors correspond
to initial gas fractions as in Figure 13. Thick dot-dashed lines correspond to
ultra-high resolution simulations. The simulations are shown during the ac-
tive/strong inflow phases, as in Figure 5. The inflow rate in the simulations
is consistent with being produced by gravitational torques over 5 orders of
magnitude in radius, with relatively little scatter (0:30:5dex at all radii).
nuclear ( 0:1 10 pc). We show the stellar and gaseous surface
densities, the local star formation rate, the gas inflow rate in com-
parison to equation 4, and the vertical scale height and Toomre Q
of the gas disk. We stress that Figure 8 shows three independent
simulations, but that the smaller-scale simulations were initialized
with properties drawn from the larger-scale simulations; this ex-
plains the approximate, but not exact, continuity between the dif-
ferent scales. Although Figure 8 is quite instructive, it can also be
misleading to focus on the results of an individual simulation. The
reason is that there is large variation in time and potentially large
scatter introduced by modest differences in galaxy properties (due
to, e.g., large-scale fragmentation of a spiral arm biasing the re-
sults of a given simulation). For this reason we believe that it is
extremely important to consider the statistical properties of a large
number of simulations that vary the precise initial conditions, gas
fraction, bulge to disk ratio, etc. In what follows we now discuss the
statistical properties of the simulations summarized in Tables 1-3.
In doing so, we will show a number of Figures that contain the re-
sults of many of our simulations. In these Figures, the critical point
to focus on is less the results of any given simulation (which can
be hard to pick out), but rather the ensemble properties of all of
the simulations (e.g., median, scatter, ...). Recall that the parameter
space spanned by our simulations includes systematic surveys in
the gas fraction, bulge (or BH)-to-disk ratio, and sub-grid equation
of state of the gas.
Figure 9 compares the order of magnitude estimate in equa-
tion 4 with the true inflow rate in the simulations, as a function
of radius; as in Figure 8, we show the results of simulations on
all three of the spatial scales we have simulated, with dotted lines
showing results near the radial boundary of a given calculation. In
calculating _Mgrav in equation 4, we approximate jaj using the mag-
nitude of the larger of the m = 1 and m = 2 Fourier component of
the surface density distribution within R (see § 5). The results in
c
0000 RAS, MNRAS 000, 000–000

How Do Massive Black Holes Get Their Gas? 19
Figure 9 are at the peak of activity for each simulation with sig-
nificant inflow, i.e., when _M is near-maximum. Modulo a constant
numerical prefactor, the simple scaling in equation (4) works well
over the entire dynamic range of our simulations, with a scatter of
 0:3 0:5dex at any radius. The same equation also works for
our ultra-high resolution simulations, which provide a smooth in-
terpolation over the boundaries in our typical “re-simulation” cal-
culations. The agreement in Figure 9 demonstrates clearly that the
inflow rates do not follow the second-order expectations from self-
torques in the absence of dissipation. In Paper II, we show that more
detailed modeling of the gas response (accounting for e.g. the de-
tails of the perturbation structure) can both predict the overall am-
plitude of _M in Figure 9 and reduce the scatter to < 0:3dex.
It is also important to note that the inflow rate in equation 4
is independent of the sound speed cs, i.e., of the thermodynamics
of the gas. Rather, the inflow rate is determined by the magni-
tude of the non-axisymmetric gravitational potential. This is one
of the reasons why we believe that our results are likely to be rel-
atively robust, even given the significant uncertainties in the ISM
physics. By contrast, if the transport were dominated by local vis-
cous stresses, the accretion rate would be given by _Mvisc 3gas,
where  ' c2s=
. For large-scale spiral waves, the angular mo-
mentum flux associated with the wave corresponds to an effective
viscosity that is at most   csr (Goodman 2003). Comparing this
with the inflow rate induced by orbit crossing (eq. [4]) highlights
that the dominant effect of the non-axisymmetric potential is that
the stellar potential induces crossing orbits in the gas; the direct
transport by spiral waves in the gas is smaller by at least a factor
of  cs=Vc. We find that in simulations that are initially 100% gas
(even those initialized with already-strong mode amplitudes), the
torques tend to remain fairly weak for a few dynamical times un-
til the stellar component is non-negligible, at which point strong
shocks appear, the stellar and gaseous modes develop a phase-shift,
and the inflow rates rise sharply. The pure gas limit is, however,
more sensitive to the thermodynamics of the gas and should be
studied in more detail in future work. In our calculations, we find
that it is critical to include stars, gas, and star formation to prop-
erly capture inflow in galactic nuclei: it is the interplay between
the stars and gas, and the competition between star formation and
inflow, that sets the net accretion rate at small radii onto the BH.
Our interpretation and analysis of the inflow rates assumes that
the gas resides in a self-gravitating, at least modestly thin, disk.
To demonstrate this explicitly, Figure 10 shows the vertical scale
height and Toomre Q parameter of the gas as a function of radius.
For each simulation, at a given radius R we determine the scale
height z0 of the gas by fitting the gas density distribution to a Gaus-
sian with dispersion z0 (after projecting to the plane perpendicular
to the net angular momentum vector). At large radii, there is scatter
owing to the difference between mergers (with sizeable z0=R) and
unstable disks. From 1100pc, z0=R ranges from 0:050:25;
this is determined by our assumed ISM equation of state, as well as
kinematic features such as twists and misalignments. At the small-
est radii z0=R begins to decrease rapidly because the potential is
increasingly dominated by the BH (the circular velocity increases
but the maximum sound speeds are finite). Note that, at all radii, the
minimum SPH smoothing length hsml . z0, so that the disk thick-
ness is at least modestly resolved. Our simulated disks are always
reasonably thin because the cooling time is always short compared
to the dynamical time, so that the sound speed in the simulation
reliably tracks that of our subgrid model.
Figure 10 also shows the Toomre Q parameter of the gas, de-
fined as Q cs=Gd , where  is the local epicyclic frequency

0.01
0.10
1.00
z 0

/ R
0.1 1 10 100 1000
R [pc]
0.1
1
10
Q ga
s(R)
Figure 10. Properties of the gas disks in our simulations. Top: Scale height
z0=R of the gas as a function of radius, for simulations with significant in-
flows. Bottom: Toomre Q parameter of the gas. The inflows are robustly
self-gravitating (Q 1) and disky, with z0=R. 0:2. Each solid line denotes
a different simulation, with dotted lines showing the radii near the bound-
aries of our re-simulations, where the exact profile shape is suspect; colors
correspond to initial gas fractions as in Figure 13. Thick lines correspond
to ultra-high resolution simulations. The range in behavior emphasizes the
importance of a large, statistically representative suite of simulations.
and d is the total surface density of the disk (gas plus stars). With
some scatter, systems largely lie near Q  1. This is broadly ex-
pected for self-regulating disks; in the simulations, Q 1 is a con-
sequence of star formation, the assumed stellar feedback, and grav-
itational dynamics.
4.1 Gas Density Profiles
Figure 11 shows the gas surface density as a function of time, at
different radii, in the high-inflow simulation from Figures 2 & 5.8
Comparison of Figures 5 & 11 clearly indicates that there is a rea-
sonable correlation between high inflow rates and high gas surface
densities at the same radii. However, the correlation is not necessar-
ily one-to-one because of the complex time-dependent dynamics on
small scales. In some cases, the inflow leads the high surface den-
sity (e.g., 0:4Myr at 0:1pc or  24Myr at 10pc), in which case
the large _M causes the high , rather than the other way around.
Occasionally, _M and  can even be anti-correlated (e.g., at 10Myr
on 10pc scales).
Figure 12 shows the gas density profiles as a function of ra-
dius for our ensemble of simulations (including the individual sim-
ulations shown in Figures 5 & 6); for the gas-rich systems in the
left panel, these surface density profiles are at the peak of activ-
ity, i.e., when _M is relatively high. During these active phases, the
8 Note that for the re-simulations, the density at small r is intentionally
initialized to be extremely small, but rapidly rises at very early times.
c
0000 RAS, MNRAS 000, 000–000

20 Hopkins and Quataert
0.0 0.5 1.0 1.5 2.0 2.5 3.0
t [Gyr]
6
7
8
9
10

1 kpc300 pc
0 2 4 6 8 10 12 14
t [Myr]
9.0
9.5
10.0
10.5
11.0
log
( Σ g
as[R
] ) [
M O •
kpc
-
2 ]
100 pc30 pc10 pc
0.0 0.2 0.4 0.6 0.8 1.0 1.2
t [Myr]
10.5
11.0
11.5
12.0

2 pc0.5 pc0.1 pc
Figure 11. Gas surface density gas at different radii as a function of time,
for the simulations shown in Figure 5. Top: Large scale galaxy merger
simulation. The box shows the time interval re-simulated in the panel be-
low. Middle: Intermediate-scale re-simulation. The box again shows the
time interval for the smaller scale re-simulation. Bottom: Nuclear-scale re-
simulation. Comparing with Figure 5, high inflow rates are generally corre-
lated with higher local surface densities, but the relation is not one-to-one.
gas reaches a quasi-equilibrium density distribution. For a given
annulus, the net accretion rate to smaller radii is typically a small
fraction of the total rate at which gas is supplied to that annulus.
Thus to first approximation, the surface density can be estimated
by considering the competition between star formation and inflow
in a given annulus. If the star formation rate surface density scales
/ 3=2gas , as in Kennicutt (1998), then setting the star formation rate
inside some radius equal to the time-averaged inflow rate _Min gives
an expected hgasi / _M
2=3
in R
4=3. Given _Min / R1=2, which is rea-
sonably consistent with Figure 5, this implies gas / R1, similar
to the results for the gas rich systems in Figures 11-12 (although
slightly steeper than the numerical results). There is of course sig-
nificant variation: we show the results from our large suite of simu-
lations to emphasize the variation introduced by initial and bound-
ary conditions, treatment of star formation and gas physics, etc., but
also to highlight the robust average behavior and the dependence on
global quantities such as gas fraction and inflow rate at larger radii
(see Fig. 13 for the color scheme).
Note that the gas surface densities achieved on small scales
can be extremely large,  1011 1012 M kpc
2 in the central
 0:110pc; this is comparable to, or even somewhat larger than,
the highest-surface density stellar systems known (Hopkins et al.
2009e; Hopkins & Hernquist 2009). And it is a factor of  104
larger than the initial gas surface density at small radii (shown
in Fig. 12 for comparison)! Assuming that a significant fraction
of the gas eventually turns into stars, the relic stellar density pro-
files expected from these gas-rich simulations are similar to the ob-
served profiles of nearby “cusp” ellipticals (e.g., Kormendy et al.
2009), which are indeed believed to be direct descendants of gas-
rich mergers. However, estimating the “final” stellar profile at small
radii, where our simulations are run for only a modest number of
local dynamical times, requires a careful correction for the effects
of duty cycle, so we defer a more detailed comparison to future
work.
Figure 12 (right panel) also shows the gas surface density pro-
files for cases in which the initial large-scale gas inflows are not
strong (e.g., Fig. 6). Most of these simulations represent either ex-
tremely gas-poor mergers or (more commonly) systems that are
weakly bar unstable or bar unstable but with significant bulge com-
ponents (such that bars form efficiently, but inflow stalls at a large
inner Linblad resonance). In such cases, the large-scale dynamics
cannot increase the gas density at sub-kpc scales very far above the
initial value of  108109 M kpc
2. Without this enhancement,
the system is stable against secondary instabilities, and so there is
little gas inflow at small radii (e.g., Fig. 6).
In both the strong and weak inflow cases in Figure 12, our
small subset of ultra-high resolution simulations (dot-dashed thick
lines; 6 galaxy-scale runs that resolve down to  10pc, and 5
intermediate-scale runs that extend to < 1 pc) are fully consistent
with the larger suite of lower resolution re-simulation runs. These
high resolution simulations have continuous inflow from large to
small radii, and can therefore be run self-consistently for a much
larger number of dynamical times than our typical resimulation.
In spite of this advantage, we find that our resimulation technique
yields very similar results.
5 CONDITIONS FOR (IN)STABILITY
In the preceding sections, we have presented results for cases in
which gas inflow from large scales both is, and is not, sufficient
to trigger secondary instabilities at small radii, leading to efficient
gas inflow down to  0:1pc. Here we provide a more quantitative
assessment of the conditions under which there is significant inflow
of gas to sub-pc scales.
Figure 13 shows how several different measures of the effi-
ciency of angular momentum transport and gas inflow depend on
the initial gas fraction fgas and the “bulge” to total ratio B=T – the
same parameters that strongly influence the morphology of the gas
(Fig. 3 & 4). We show results for a few different equations of state
qeos (see Fig.1). The quantities we use to define the efficiency of gas
inflow are the gas mass within 10 and 0:1pc Mgas(< R) (top panel),
the inflow rates within the same radii (middle panel) and the for-
mal amplitude of the non-axisymmetric gravitational perturbations
at  1 and  100 pc (bottom panel). For the inflow rate (middle
c
0000 RAS, MNRAS 000, 000–000

How Do Massive Black Holes Get Their Gas? 21
0.1 1 10 100 1000
R [pc]
107
108
109
1010
1011
1012
1013
Σ ga
s
[M O •

kpc
-
2 ]
Gas-Rich MergersInitial Gas Profile
0.1 1 10 100 1000
R [pc]
107
108
109
1010
1011
1012
1013

Weakly Barred Disks
Figure 12. Left: Gas surface density profiles during the strong inflow phases, i.e., when the accretion rate at small radii peaks. Each solid line denotes a
different simulation, with dotted lines showing the radii near the boundaries of our re-simulations, where the exact profile shape is suspect; colors correspond
to initial gas fractions as in Figure 13. Thick lines are the ultra-high resolution simulations. Each simulation consists of re-simulations of the nuclear dynamics
during gas-rich galaxy-galaxy mergers, as in Figures 5 & 11. Dashed red line shows the initial gas density profile of the large-scale simulations (run b3ex(ic)
in Table 1). The gas density profile is quasi-steady over the entire active phase, with star formation offsetting inflow; the time variation in gas for one example
is shown in Figure 11. Right: Gas surface density profiles for simulations of isolated barred-disk galaxies and the corresponding re-simulations. In cases with
a very small pre-existing bulge, the central gas density can increase by an order of magnitude, but not the 24 orders of magnitude seen in the left panel.
panel), we show both the time average and “peak” values to convey
the variability as a function of time in the simulation.
The mode amplitude is measured in the standard fashion from
the surface density distribution at a given radius, using
jam(R; t)j=
jj
R 2
0 (R;) exp(im)djj
R 2
0 (R;)d
: (5)
We measure the amplitude from the stellar disk since it is slightly
more robust to local clumping, but using the gas surface density
gives similar results. For simplicity, we show the results for the
largest-amplitude mode in each case, but this is almost always m =
1 in the nuclear-scale simulations (right panels) and an m = 2 bar
or m = 1 spiral in the intermediate-scale simulations (left panels).
We measure the relevant amplitude at radii slightly larger than the
radii where the inflow is measured, since it is this non-axisymmetry
that drives the inflow (it is also somewhat less noisy at larger radii).
The “bulge-to-total” ratio in Figure 13 is defined as the mass
fraction in all spherical components within the specified radius,
since these all serve to suppress instability in the same manner (e.g.
any black hole, bulge, halo, and/or nuclear star cluster would qual-
ify). Gas fraction is defined as the gas fraction of the disky compo-
nent within the same radius.
Figure 13 shows clearly that there is a very strong trend of
inflow strength with B=T , in the expected sense. Strongly disk-
dominated systems have characteristic inflow rates 3 orders of mag-
nitude higher than strongly “bulge-dominated” systems; likewise
mode amplitudes and gas masses inside small radii are higher by
12 orders of magnitude.
There are several interesting additional results in Figure 13.
First, within the range of qeos that we consider (motivated by the
observations in Figure 1), the inflow properties do not depend sig-
nificantly on the ISM physics. This is consistent with our previous
arguments; so long as the sub-resolution model is sufficient to pre-
vent catastrophic fragmentation and runaway star formation, gravi-
tational torques that are reasonably independent of the ISM micro-
physics dominate the angular momentum redistribution. Note that
there is also no distinguishable difference between the results of
our ultra-high resolution runs (shown with large symbols) and our
standard resimulations. Second, the inflow properties depend much
more strongly on B/T than on the initial gas fraction; this is because
the stellar disk generates the same instabilities as the gas disk. In
fact, the collisional nature of gas means it cannot support certain
strong modes (Toomre 1977), so the transport can be less efficient
in the fully gas-dominated limit (the gas-rich systems modeled here
often do not develop strong modes until a non-trivial fraction of the
gas has turned into stars).
Third, there is considerable scatter at all B=T . This is in no
small part due to the time-variable nature of the small scale dynam-
ics (§3). In addition, however, a quantity such as B=T is not a global
invariant; it depends on both position and time in these simulations.
As such, some simulations with a given B=T at the specified radius
can have a larger or smaller B=T at some larger radius, leading to
more or less efficient inflow from those scales that in turn leads to
more or less efficient inflow to smaller scales. Identifying B=T at
a given radius and time is thus a crude measure of stability. More-
over, the torques and inflow rates depend on the precise structure of
the modes involved; a significant portion of the scatter in Figure 13
reflects the difference between e.g. bars, spiral arms, and loosely or
tightly-wrapped waves. This large scatter highlights the importance
of using a large suite of simulations to conduct parameter surveys
and build statistical samples – any individual system, studied for a
small amount of time, can be highly non-representative.
5.1 Intermediate (sub-kpc) Scales
For the intermediate scale results shown in the left panel of Fig-
ure 13, we can gain some analytic understanding of the behavior
by following Shlosman et al. (1989), who discuss the criteria for
secondary instabilities (“bars within bars”). It is well established
c
0000 RAS, MNRAS 000, 000–000

22 Hopkins and Quataert
0.0 0.2 0.4 0.6 0.8 1.0

5
6
7
8
log
[ M
ga
s(<R
) ]
[M O •
]
R < 10 pc
fgas(Initial): 0.1 0.2 0.3 0.5 0.8 0.9
ISM Equation of State:
qeos=0.20-0.35qeos=0.125-0.20qeos=0.05-0.125
0.0 0.2 0.4 0.6 0.8 1.0

0.1
1
10
100
Inf
low
dM
/dt
[M
O •
yr-
1 ]
R < 10 pc
0.0 0.2 0.4 0.6 0.8 1.0
B/T (<300 pc)
0.01
0.1
1
〈Mo
de
Am
plit
ud
e〉
R < 100 pc
0.0 0.2 0.4 0.6 0.8 1.0

2
3
4
5

R < 0.1 pc
0.0 0.2 0.4 0.6 0.8 1.0

0.01
0.1
1
10

R < 0.1 pc
Time-AveragedMaximum Sampled
0.0 0.2 0.4 0.6 0.8 1.0
B/T (<10 pc)
0.01
0.1

R < 1 pc
Figure 13. Inflow strengths as a function of bulge to total ratio (B/T) for several initial gas fractions (colors) and ISM equations of state qeos (symbol type),
for a subset of our simulations. Large symbols in each refer to the subset of ultra-high resolution simulations. Top: Total gas mass inside an inner radius of
10pc (left; from our intermediate-scale re-simulations) or 0:1pc (right; from our nuclear-scale re-simulations), near the peak of accretion in each simulation.
“Bulge” (B) refers to any spherical term in the potential (bulge+halo+black hole). Middle: Gas inflow rate through the relevant inner radius. Solid points denote
the time-averaged inflow rate over the full re-simulation while open points are the maximum seen in any individual snapshot (to give some sense of the duty
cycle). Bottom: Amplitude of the dominant non-axisymmetric mode in the stars (generally m = 2 at intermediate scales, m = 1 at nuclear scales), averaged
over the peak/strong inflow phase. Inflow properties scale broadly with B=T as expected, but the large scatter and sensitivity to other parameters makes a large
parameter survey critical.
that a self-gravitating disk becomes unstable to large-scale gravi-
tational modes roughly when the Ostriker-Peebles criterion is met,
i.e., when Trot=jW j > tcrit, where Trot Md V 2c is the kinetic energy
of rotation (in the rotationally supported – i.e. kinematically cold
– component), jW j  GM2=R is the potential energy, and tcrit is a
threshold value ' 0:15 0:26, depending on the nature of the in-
stability (e.g., Bardeen 1975; Aoki et al. 1979). If a mass fraction
fd = 1B=T is in the disk, the condition for instability generically
becomes
fd > fcrit(tcrit; a=R) (6)
where a is the bulge scale length, R is the disk scale length, and the
exact functional form of fcrit depends on the profiles of the disk and
bulge. For Kuzmin and Plummer profiles for the disk and bulge,
c
0000 RAS, MNRAS 000, 000–000

How Do Massive Black Holes Get Their Gas? 23
respectively, we find that this instability criterion becomes
Md[< R]
Mb + Md[< R]
>

3 tcrit
4[12 tcrit]
R
a
1=2
(R a) (7)

3 tcrit
16
[1 + (a=R)2]
R
a

(R a) (8)
where Md[< R] is the total (stars + gas) disk mass within R and Mb
is the total bulge mass.
If the disky material of interest has a scale-length compara-
ble to that of the bulge (generally true for the intermediate-scale
simulations here), equation (6) simply becomes B=T . 1 tcrit '
0:70:8. This is reasonably consistent with the fact that the inflow
rate and interior gas mass in the intermediate scale simulations in
Figure 13 (left column) asymptote for B=T & 0:60:7.
Shlosman et al. (1989) show that an initial gas disk must be
compressed to at least Rnew=Ri  (1:0 2:5) f 2gas=[1 + 0:2(B=T )
2]
to satisfy equation (6), again for systems in which the disk and
bulge initially have similar scale-lengths. For the case of major
mergers, we can compare this to Covington et al. (2008) and Hop-
kins et al. (2009a,c)’s estimates of how much angular momentum
the gas loses during the merger. These authors find that the gas of-
ten inflows until it becomes self-gravitating – as a result, secondary
bar instabilities should be ubiquitous. More quantitatively, they find
Rnew=Ri  fgas=(1+ fgas= f0) with f0  0:20:3. Using this, we es-
timate that mergers with reasonable gas fractions & 0:3 0:5 will
lead to secondary instabilities and rapid inflow. For isolated disk
galaxies with bars, the critical criterion is whether the inner Lin-
blad resonance lies inside or outside of Rnew; this will generally not
be true in low- fgas or high-B=T systems.
5.2 Nuclear (pc) Scales
Inside the radius where the BH begins to dominate the gravita-
tional potential, it becomes increasingly difficult for system to be
globally gravitationally unstable in the sense of Trot=jW j > tcrit.
Specifically, in the potential of the BH this criterion implies that
the disk would no longer be unstable inside a radius Rmin 
10pc(MBH=108 M)
1=2 (d=1011 M kpc
2)1=2. This confirms
our earlier claims (e.g., §3.3) that the character of the instabilities
responsible for angular momentum transport must change near the
BHs radius of influence. Indeed, we see numerically that the trans-
port on small scales is dominated by m = 1 modes (e.g., Fig. 4).
A full discussion of the origin of the m = 1 modes is beyond
the scope of this paper; we will present a detailed analytic analy-
sis (in simplified BH+disk models) of their structure, growth rates,
and pattern speeds in future work (in preparation). We can, how-
ever, present a simple analysis indicating how we believe the modes
arise. For m = 1 modes,
=m! 0 in a Keplerian potential. It is
well-known that this allows very low frequency (low pattern speed

p) modes to be present in a Keplerian potential (e.g., Tremaine
2001). Moreover, Tremaine (2001) showed that the only modes that
can be global and exert coherent strong torques over a large dy-
namic range in a quasi-Keplerian potential are these m = 1 modes.
We find that the modes in the simulations indeed have small
p,
i.e.,
p
and m = 1, deep in the potential well of the BH. Note
that for a collisionless particle, rather than the stellar+fluid disks of
interest here, this simply corresponds to a standing elliptical Ke-
plerian orbit. Small
p=
(and wavenumber m = 1) is important
because it implies near-resonance conditions at small radii; this al-
lows the stars to produce strong torques on the gas, which drives or-
bit crossings, shocks, and angular momentum loss, even when the
mass of the disk at small radii is much less than the BH mass (see
e.g. Chang et al. 2007). Following Tremaine (2001) and expanding
the equations of motion for a linear mode in a nearly-Keplerian po-
tential (in the WKB limit) it is straightforward to show that the
magnitude of the induced velocities (v=Vc) from a given mode
scale as v=Vc  (=)Md(< R)=MBH = jajMd(< R)=MBH for
all m 6= 1; for m = 1, however, v=Vc  jaj. In other words, the con-
ventional wisdom that torques from global modes are strongly sup-
pressed in the spherical potential of a BH is generally true, but the
resonance condition that  '
in Keplerian potentials allows for
large coherent eccentricities, self-gravitating collective effects, and
torques from m = 1 modes at arbitrarily small disk-to-BH masses.
As a result, although we certainly see higher-m modes in our small-
scale simulations, the global structure is almost always dominated
by a coherent global m = 1 eccentric disk.
However, Tremaine (2001) also showed that slow modes with
pattern speeds
p 
, are linearly stable in quasi-Keplerian po-
tentials. How, then, do these modes arise in the simulations? To
start, we consider the dispersion relation for linear waves in the
WKB limit (see e.g. Lau & Bertin 1978), which is given by
!m



2
= 1 +2 f˜dA

jkrj2 + m2
1=2
+ A

jkhj2 +
m2h2
r2

(9)
where
f˜d 
d R2
Menc(< R)

Md(< R)
Mtot(< R)
; (10)
A = 1 +
2m2(3)
(+ 1)(jkrj2 + m2)
; (11)
  1 + 2
@ lnVc
@ lnR

dlnMenc(< R)
dlnR
; (12)
and h = cs=
is the disk thickness. We now consider the limit of
equation (9) such that kr is modest. Of course, our analytic treat-
ment of such quasi-global gravitational modes should be consid-
ered with due caution, as the WKB approximation necessary to de-
rive the dispersion relation is not justified. But this nevertheless
points towards the conditions under which global gravitational in-
stability may be present.9 In this limit, equation (9) will have grow-
ing modes if the right-hand side is negative, which for a cold disk
(h r) requires
f˜d >
1
2m
(1 +)2
7
: (13)
As  ! 0, i.e., as the potential becomes Keplerian, equation 13
suggests that m = 1 modes can become unstable when f˜d & 0:1.
This is consistent with the fact that we find the m = 1 modes at
nearly all B=T (Fig. 13). Equation 13 also allows for the instability
of higher-m modes, but these cannot exert coherent strong torques
on the gas, nor can they propagate efficiently to small radii, so they
do not dominate the structure or inflows we see in our simulations.
Because the disk in general will have more mass at large radii,
any instability will grow outside-in; for massive disks, the most-
unstable point is where Md(< R) 0:51MBH ( f˜d  1=31=2).
9 Note also that equation 9 is derived with the WKB approximation, but in-
cluding terms up to second order in jk Rj1 (dropping out-of-phase terms),
which makes it somewhat more accurate for the large-scale modes of inter-
est here. Specifically, this leads to the additional enhancement of the growth
rate A  1 + sin i, where  @ ln
=@ lnR and i is the arm pitch angle;
this approximates (at this order) the effects of the swing amplifier, an im-
portant contribution to the instability criterion.
c
0000 RAS, MNRAS 000, 000–000

24 Hopkins and Quataert
Here modes can grow on a dynamical time. Equation 9 shows that
this growth is for fast modes: overstabilities with Re(!) m
. In
the simulations, the slow modes typically dominate the torques at
sub-pc scales. However, we find that the fast mode itself drives
the slow modes. Specifically, the m = 1 mode first appears around
the most unstable point Rcrit predicted above, where Md(< Rcrit)
0:5 1MBH. This rapidly leads to an eccentric disk at these radii.
As the instability causes gas to lose angular momentum, the gas
falls inwards, and begins turning into stars. The eccentric mass dis-
tribution at Rcrit then in turn efficiently induces a similar eccen-
tricity at smaller radii. As noted in Tremaine (2001), the m = 1
slow mode may be linearly stable, but it can be excited in the first
place – and will be long-lived – in response to an external asym-
metry in the potential. Our interpretation of the simulations is that
the unstable self-gravitating disk near  Rcrit provides the asym-
metry that generates the slow, low pattern speed, m = 1 mode at
smaller radii. The pattern speed
p =
(Rcrit) is conserved as the
mode moves in, but the circular speed
(R) is rapidly increasing
at smaller radii near the BH, so the mode at small radii is indeed
a slow mode – for the typical parameters here, the pattern speed is

p  15kms1 pc1, set by the circular speeds at 10100pc,
where Md(< R)MBH (for comparison,
& 1000kms1 pc1 at
R< 1pc). The slow mode at small radii responsible for the angular
momentum transport is therefore not a local small-scale instabil-
ity, but part of a global instability beginning at larger radii. In the
future, we plan to explore the origin of this global m = 1 mode
in more detail. Calculating the dynamics of the gas and stars, and
including self-gravity and star formation, are all important for cap-
turing the inward propagation of these modes.
Similar m = 1 modes have been studied previously as a mech-
anism for accretion in fluid disks (Adams et al. 1989; Shu et al.
1990; Ostriker et al. 1992; Bournaud et al. 2005), both in the galac-
tic context and in proto-planetary or proto-stellar disks (e.g., Baner-
jee et al. 2004; Boley et al. 2006; Vorobyov & Basu 2006, 2009;
Krumholz et al. 2007; Cai et al. 2008). In a purely gaseous disk,
wave packets can propagate through the OLR to r!1, eventually
becoming simple sound waves (Adams et al. 1989). Because the
waves can freely escape to the outer boundary, infinite pure gaseous
disks in nearly-Keplerian potentials do not support strong growing
m = 1 modes, and the mode growth is quite sensitive to the descrip-
tion of the disk boundary. However, for a stellar disk, the mode can-
not propagate beyond the OLR where
 2m(

p)2 = 0.
Refraction of the stellar waves at this boundary is important, and
implies that global m = 1 modes can grow even when the disk ex-
tends to R ROLR. As a result, we shall show in subsequent work
that the m = 1 modes in quasi-Keplerian stellar disks are relatively
insensitive to the outer boundary conditions. This is also implicit
in Tremaine (2001), where for many of the modes the outer radius
of the disk could be arbitrarily large. The ability of stellar m = 1
modes to grow more readily again highlights the importance of con-
sidering the two component stellar + gas system when considering
gas inflow in galactic nuclei.
6 ACCRETION VERSUS STAR FORMATION
Given the gas inflows modeled here, our results imply a correlation
between the accretion rate onto the central BH – and hence the
AGN luminosity – and the star formation rates on larger scales.
Figure 14 shows this correlation for a number of our simulations.
We take the inflow rate at 0:1 pc as the AGN accretion rate, which
corresponds to a bolometric luminosity of
Lbol; AGN = r _MBH c
2  5:71045
 r
0:1
 _MBH
M yr1

ergs1
(14)
where the radiative efficiency is r = 0:1. Figure 14 compares the
BH accretion rate to the total star formation rate inside a (pro-
jected) radius R, at a few representative radii: true nuclear scales
(R < 10pc), more readily observable nuclear scales (R < 100pc),
the smallest scales observable in most moderate to high-redshift
systems (R < 1kpc), and the total star formation rate, integrated
over all radii (R!1). For each point, in Figure 14, the projected
SFR(< R) is technically the median over 100 sightlines, but pro-
jection effects are negligible. For the larger radii, we measure the
SFR in an appropriate intermediate or galaxy-scale simulation. The
BH accretion rate is then determined from that in the appropriate
re-simulation or set of re-simulations matched to the galaxy scale
simulation central properties at the given time. Because a given dy-
namical time in a galaxy-scale simulation corresponds to many dy-
namical times (and simulation outputs) on very small scales, the
number of points is larger in the plots going out to larger scale.
Not surprisingly, Figure 14 shows that there is a correlation
between accretion rate and star formation rate, on all scales. There
is significant scatter in this correlation, although the scatter is less
when the SFR is evaluated inside a small nuclear annulus. If we
assume a linear proportionality between BH accretion rate and star
formation rate in some annulus, we find
_MBH  _M(R< 10pc) (15)
 0:1 _M(R< 100pc) (16)
 0:03 _M(R< 1kpc) (17)
 0:003 _M(total) (18)
for the median relation in the simulations, which corre-
sponds to a bolometric luminosity ratio of Lbol; SF=Lbol; AGN 
(0:01; 0:1; 0:3; 2:5) at (R < 10pc; R < 100pc; R < 1kpc; R <
1), respectively, or a ratio of the star formation-powered in-
frared luminosity to the hard X-ray luminosity from the AGN of
LIR; SF=LHX; AGN  (0:2; 2; 6; 60) in the same annuli.
In greater detail, the correlations are not precisely linear. Pa-
rameterizing the correlation between star formation and BH ac-
cretion as a power-law, _MBH / _M , we find that the correlations
above are approximately valid for _MBH  0:110M yr
1, but 
changes from near-linear at small scales (R < 10pc) to somewhat
super-linear (  1:5) on large scales. For the total star formation
rate, we find
_MBH  0:5M yr
1
 _M
50M yr1
1:5
: (19)
However, we caution that our simulations do not include the ap-
propriate physics to resolve the variety of processes that can pro-
duce low accretion rates 0:1M yr
1, nor do they sample a large
regime of galaxies with total star formation rates . 1M yr
1, so
the fits above should not be used in the limits of weak accretion
and star formation. The qualitative point is nonetheless important:
non-axisymmetric instabilities due to self-gravity become ineffi-
cient at low gas densities (low star formation rates), and so other
mechanisms must be important for fueling at least some of the low-
luminosity AGN population.
It is notable that even during the active phases modeled here,
the average relation between accretion rate and total SFR is such
that the systems are most often not AGN dominated in a bolometric
c
0000 RAS, MNRAS 000, 000–000

How Do Massive Black Holes Get Their Gas? 25
0.01 0.1 1 10 100

0.01
0.1
1
10
BH
AR
[M
O •
yr
-
1 ]
RSF < 10 pc
SF
R
SF
R/1
0
SF
R/1
00
SF
R/1
000
0.1 1 10 100 1000

0.01
0.1
1
10

RSF < 100 pc
1 10 100 1000
SFR(<R) [MO • yr-1]
0.01
0.1
1
10
BH
AR
[M
O •
yr
-
1 ]
RSF < 1 kpc
1 10 100 1000
SFR(<R) [MO • yr-1]
0.01
0.1
1
10

RSF < 50 kpc
Figure 14. The predicted relationship between the BH accretion rate (BHAR) and the star formation rate (SFR) inside a given radius. For each simulation,
we plot the SFR integrated within the indicated projected radius, at different times during the simulation, compared with the < 0:1pc BH accretion rate
from appropriate re-simulation. The characteristic timescale for variability at small radii – which sets the variability in the BHAR – is much shorter than the
corresponding timescale at large radii; at fixed SFR, each system can thus have a variety of BHARs. Nevertheless, the BHAR correlates with SFR, with the
least scatter for the nuclear star formation (upper left). The correlations seen here are a natural consequence of the fact that high gas densities, and thus high
star formation rates, are necessary to trigger secondary instabilities and drive accretion on small scales.
sense (although the two are certainly comparable, in these phases).
In other words, even the most active systems will be consistent with
SFR-dominated properties except for a small fraction of the time;
this is consistent with observations such as the FIR-radio correla-
tion (Walter et al. 2009; Riechers et al. 2009; Wang et al. 2008).
However, the scatter in accretion rates at fixed SFR is large, and
most of the growth in BH mass comes from the rare, short-lived
times when the systems scatter to higher LAGN=LSF. Roughly, we
can estimate these duty cycles with the scatter in Figure 14 – about
 30% of the systems in “peak” AGN phases have sufficient accre-
tion rates ( _MBH & 0:01 _M) so as to be AGN-dominated.
The predictions in Figure 14 can be compared to a number of
observations. Wu & Cao (2006) compile total star formation rate
estimators in  100 nearby AGN ranging in luminosity from LIN-
ERS ( _MBH < 0:01M yr
1) to PG quasars ( _MBH 15M yr
1).
The sample details are discussed in Satyapal et al. (2005), but the
authors find they are reasonably robust to the choice of SFR or BH
accretion rate indicator. Silverman et al. (2008) consider a simi-
lar compilation in more luminous systems in the COSMOS sur-
vey, spanning a redshift range z  0:3 1:1. At total SFR values
> 10M yr
1, these samples agree fairly well. Specifically, for
_M  101000M yr
1, these different compilations are consis-
tent with a median _MBH  0:0010:01 _M, with factor of  3 dis-
persion; this is quite similar to our results shown in Figure 14. On
the other hand, Ho (2005) used OII line emission to argue that the
star formation rate associated with luminous type 1 AGN is at best
comparable to the BH accretion rate, inconsistent with the median
integrated galaxy predictions in Figure 14. This may be because ob-
scuration makes OII an imperfect proxy for star formation in dense,
obscured, starbursts. In addition, luminous Type 1 AGN are prob-
ably, by selection, on the tail end of the BHAR/SFR distribution
(and could in principle be in an AGN feedback-dominated stage of
evolution not modeled here). Indeed, Kim et al. (2006) followed up
Ho’s work and found enhanced star formation in type 2 quasars,
similar to our predictions here.
Wang et al. (2007) consider spatially resolved estimates of
nuclear luminosity and emission-line (SFR) profiles in a range of
Seyferts and quasars (mostly at low redshift) at spatial scales from
 0:053kpc (see Imanishi 2002, 2003; Imanishi & Wada 2004).
Their results are similar to those in a variety of other studies (see
e.g. Kauffmann et al. 2003; Kewley et al. 2006; Evans et al. 2006;
Schawinski et al. 2009). At both R < 100pc (where the observa-
tions span a dynamic range in SFR(< 100pc) 0:110M yr
1)
and R < 1kpc (SFR(< 1kpc)  1 100M yr
1), our results are
reasonably consistent with these observations, in both normaliza-
tion and scatter. A small sample is studied at extremely high res-
olution in Davies et al. (2007), where star formation and AGN lu-
minosity are resolved at  1030pc scales. They find that, within
c
0000 RAS, MNRAS 000, 000–000

26 Hopkins and Quataert
these radii, the luminosity from star formation is  1 5% of that
from the BH, again consistent with our prediction here.
7 DISCUSSION
We have presented high resolution smoothed particle hydrody-
namic simulations of the dynamics of gas in the central regions of
galaxies, in order to study gas inflow from galactic scales  10kpc
down to the scales of a canonical AGN accretion disk ( 0:1pc).
We use the properties of galactic-scale simulations – both merg-
ers and isolated galaxy disks – to initialize new higher resolution
“re-simulations” that focus on the . 100 pc dynamics. This tech-
nique is then applied a second time to follow the gas to  0:1 pc.
By intention, our calculations do not employ particle splitting and
are not all exact realizations of the small-scale dynamics associated
with a single, specific larger-scale simulation. Our “re-simulations”
should not be taken as such an exact realization. Instead, since it is
not practical to directly solve the “full” problem (from  10kpc to
105 pc), we have focused on understanding and isolating the key
physics involved. To provide a more physical model, however, we
use initial conditions drawn from a variety of larger galaxy-scale
simulations, rather than an ad hoc set-up. This approach allows
us to systematically survey, for the first time, how a large variety
of initial conditions and galaxy properties affect gas dynamics in
galactic nuclei.
It is still not feasible, with any numerical scheme, to simul-
taneously model the small scales of an AGN accretion disk and
those of a galaxy for cosmological timescales. For this reason, di-
rect “zoom-in” approaches have thus far been fairly limited, usually
evolving each region for just a few local dynamical times. Our nu-
merical approach attempts to overcome this severe limitation. First,
the smallest regions of each simulation are simulated for many lo-
cal dynamical times. Even more importantly, we quantity the dy-
namics on small scales using  100 simulations of a large range
of initial conditions and galaxy properties This approach does not
require an exact mapping between one large-scale simulation and
another on smaller scales (although we often interpret our results in
this manner). For example, even if our nuclear-scale simulations of
the central 100pc of a galaxy are not self-consistently matched to
those of a specific larger-scale simulation, they still represent valid
simulations of a given set of initial conditions of a bulge+stellar
disk+gas disk+BH system, evolved for many dynamical times. As
such, we can still use them to understand the physics of angular
momentum transport and BH fueling at these radii (given these ini-
tial conditions). This methodology allows us to survey a wide range
of initial conditions and galaxy properties and isolate the physical
parameters that are most important for the physics of gas inflows
on small scales in galactic nuclei. We have also checked our "re-
simulation" approach by carrying out a small number of ultra-high
resolution simulations that bridge the gap between different resim-
ulated regions; we find in all cases that these ultra-high resolution
simulations validate our resimulation methodology.
Our calculations demonstrate that for sufficiently gas-rich,
disk-dominated systems, which have a large inflow of gas into
kpc scales, there is a cascade of non-axisymmetric gravitational
instabilities that ultimately produces BH accretion rates as high as
 10M yr
1 at . 0:1 pc (Fig. 13). Moreover, we have explicitly
shown that these conditions are satisfied during major mergers of
gas-rich galaxies, thus providing theoretical support for the con-
nection between mergers and BH growth. It is, however, also im-
portant to stress that our work only implies that galaxy mergers are
a sufficient condition for fueling a central quasar, not that they are
necessary: similar inflows can be obtained in isolated gas-rich, bar
or spiral wave-unstable galactic disks.
In broad terms, our results support the “bars within bars” sce-
nario proposed in Shlosman et al. (1989). However, we show that
the secondary instabilities on intermediate scales ( 10 100 pc)
exhibit a diverse range of morphologies, including standard nuclear
spirals, bars, rings, barred rings, one or three-armed spirals, and
irregular/clumpy disks and streams (Fig. 2 & 3). This is very im-
portant for comparing to real galaxies: observations have generally
found that nuclear bars are not ubiquitous in AGN (perhaps not
substantially more common than in non-active systems). However,
these same observations have found that some form of asymmet-
ric nuclear structure is ubiquitous, with the most common features
being nuclear spirals and rings very similar to those seen in Fig-
ure 3 (see Martini & Pogge 1999; Martini et al. 2003; Peletier et al.
1999; Knapen et al. 2000; Laine et al. 2002; Knapen et al. 2002;
Greene et al. 2009, and references therein). Our work shows that
these observations are not necessarily in conflict with the hypothe-
sis that gravitational torques dominate inflow on sub-kpc scales in
AGN. Rather there is a large diversity in the observational mani-
festation of these gravitational torques. “Bars within bars” should
not be taken too literally. Instead a more accurate characterization
would be: “it’s non-axisymmetric features all the way down” (or
“stuff within stuff”).
This diversity in observational appearance owes both to the
effects of global parameters such as gas fraction and bulge-to-
total ratio, and more subtle variations induced by the relative scale
lengths of disks and bulges, the precise profile shapes, and the
thickness/dispersion in the sub-kpc bulge, stellar disk, and gaseous
disk components. Observationally, these parameters are all seen to
vary considerably even in isolated galaxies, let alone in chaotic
merging systems; as a result, we expect the morphological diver-
sity exhibited in Figures 2 & 3 to be the norm.
Once the gas reaches the BH radius of influence ( 10pc for
a typical  L system), the potential becomes quasi-Keplerian and
spherical (dominated by the BH itself); the system is thus no longer
bar unstable. However, the gas is still locally self-gravitating and
prone to forming stars if it does not inflow rapidly. Indeed, this is
traditionally the range of radii at which it has been the most difficult
to produce the large accretion rates needed for luminous AGN and
quasars (Shlosman et al. 1990; Thompson et al. 2005). We find,
however, that new large-scale non-axisymmetric modes arise at 
pc-scales and robustly allow for continued gas transport to smaller
radii. Specifically, for sufficiently gas-rich, disky systems, we find
an eccentric/lopsided disk or a one-armed spiral instability (an m =
1 mode); see Figure 4. This leads to large inflow rates of  1
10M yr
1 into the central 0:1 pc, comparable to what is needed
to explain the brightest quasars observed. At radii . 0:010:1 pc,
the accretion flow is no longer self-gravitating and local viscous
heating is sufficient to maintain Q & 1 and transport gas to the BH
(Goodman 2003; Thompson et al. 2005); the disk is approaching a
canonical thin Keplerian accretion disk.
The m = 1 modes seen here have been studied previously as
a mechanism for accretion onto BHs and protostars (Adams et al.
1989; Shu et al. 1990; Ostriker et al. 1992), but have largely been
neglected in studies of fueling luminous AGN. This is in part be-
cause most previous simulations neglected either star formation or
the self-gravity of the gas and stars (owing to computational limits
or, of course, the fact that these complications are not present in the
case of a protostar). Moreover, previous calculations have shown
that although near-Keplerian potentials can support the “slow” pat-
tern speed m = 1 modes we find here, such modes are linearly
c
0000 RAS, MNRAS 000, 000–000

How Do Massive Black Holes Get Their Gas? 27
stable (Tremaine 2001). Why, then, are these modes ubiquitous in
our simulations? We find that the inclusion of stars and gas, self-
gravity, live star formation, and some model for stellar feedback are
all important for the behavior of these instabilities. In particular, we
typically find a complex interplay between the gas and stars (§4 &
5.2). The system often develops the m = 1 modes first near co-
rotation, where the disk and BH contribute comparably to the po-
tential; the mode growth is faster in the stellar distribution, which
can support self-crossing orbits and which is much less sensitive
to the outer properties of the disk. As gas streams move through
the lopsided stellar distribution, they are torqued, experience or-
bit crossing, and lose angular momentum and energy (as in, e.g.,
Chang et al. 2007). As the gas moves to smaller radii, the stars at
larger radii provide a less efficient angular momentum sink. How-
ever, star formation ensures that new stars are formed in situ (them-
selves in a lopsided distribution), allowing the gas to continue to
flow in. This process leads to the m = 1 mode propagating inwards
into the gravitational potential well of the BH, from the larger radii
where it originated. The torques and inflow rates are thus quite dif-
ferent in stellar+gas systems as compared to pure gas disks. In par-
ticular, orbit crossing induced by the stars can lead to much more
efficient gas inflow; we will study this analytically in Paper II.
Our calculations show that the inflow of gas on sub-kpc scales,
including the BH accretion rate at . 0:1 pc, is highly time vari-
able (Fig. 5). This variability is a result of gravitationally unstable
modes forming, damping, and, re-forming, and the fact that individ-
ual clumps or dense regions can dominate the inflow rate at a given
time. The “duty cycle” for large inflow rates is modest, even in sys-
tems with a large time averaged accretion rate. In one of our most
violent simulations, a massive nuclear ring forms at kpc; owing
to an inner Linblad-like resonance, the system has an outflow from
sub-kpc scales, and inflow from the merger at larger radii. Since
the outflow/inflow rates are large, the ring soon becomes strongly
self-gravitating and globally collapses, leading to one of the largest
net accretion rates into  10pc in any of our simulations. How-
ever, more than  95% of the time during the “active” phase, the
system has strong outflow from the sub-kpc region. Moreover, be-
cause the gravitational modes driving accretion at large scales (e.g.
the merger), at  100 pc (secondary instabilities), and at pc (the
lopsided nuclear disk) are physically separate, their variability on
short timescales is not fully coupled. This is clearly important when
relating observed torques on large scales to AGN accretion rates on
small scales. In many observed AGN or barred disks (with a proper
inner Linblad resonance), there are no strong inflows (and may even
be outflows) observable at kpc (Block et al. 2001; García-Burillo
et al. 2005; Haan et al. 2009; Stoklasová et al. 2009; Durbala et al.
2009). Our simulations demonstrate that this can in fact be the case
even in systems with large time-averaged accretion rates.
The key parameter determining the qualitative behavior at
each scale in our simulations is the ratio of the total mass in a ro-
tationally supported disk (gas and stars) to the pre-existing mass in
any spherical component (bulge, halo, or black hole); see Figure 13
and § 5. There is, of course, a considerable literature studying the
instability of self-gravitating disks in disk/bulge systems, which is
applicable to the intermediate scales ( 10 100 pc) in our simu-
lations (see e.g. Bardeen 1975; Athanassoula 2008; Narayan et al.
1987; Raha et al. 1991; Christodoulou et al. 1995; Earn & Sell-
wood 1995; Bournaud & Combes 2002; Dubinski et al. 2009). As
outlined in Shlosman et al. (1989), most of these criteria amount
to variations on the Ostriker-Peebles criterion that the kinetic en-
ergy of rotation inside some radius be & 20% of the potential en-
ergy. This requires increasing the surface density of disky material
within the pre-existing bulge/halo effective radius to a value com-
parable to, or equal to, that of the bulge/halo. For major, gas-rich
mergers ( fgas & 0:3), this condition is almost always met. However,
for gas-poor mergers, or for isolated barred galaxies, it is not clear
how often systems will in fact trigger secondary instabilities. We
find cases where weak large-scale bars trigger no significant subse-
quent activity on smaller scales (Figure 6).
In the gravitational potential of the central BH, the criterion
for further inflow is somewhat altered, although we believe that it is
qualitatively similar (§5.2): lopsided m = 1 modes are dynamically
important provided that the total mass in the disky component is
 0:11MBH in the central 10 parsecs (Figure 13; note the large
scatter even at fixed disk to BH mass). We find that this criterion
is generally satisfied, so long as the larger scales ( 10 100 pc)
are themselves unstable. The critical barrier, then, to triggering a
cascade of instabilities leading to significant BH accretion is the
presence of a large self-gravitating gas inflow at  a kpc.
There is significant evidence that galaxies in fact meet the cri-
teria outlined here for AGN fueling. Studies of ongoing gas-rich
mergers find that the gas densities reach and exceed the self-gravity
threshold, with gas densities similar to those predicted here (Fig-
ures 11-12; see e.g. Tacconi et al. 1999; Hibbard & Yun 1999; Laine
et al. 2003; Iono et al. 2007; Schweizer & Seitzer 2007). Similarly,
studies of early-type galaxies find that the stellar surface densities
can be separated into a distinct “dissipational” component – the
stars at small radii which must have formed in a dissipational star-
burst – and a “dissipationless” envelope of violently relaxed stars
(see Hopkins et al. 2008a). The starburst relic component is inferred
to have formed from self-gravitating gas; indeed, this appears to
define its size-mass relation (Hopkins et al. 2009a). In addition,
the criterion for significant gas inflow given a modest pre-existing
bulge corresponds to a surface density in the rotationally supported
component of 1011 M kpc
2 at 10pc (Figure 12), for an L
system. This is comparable to the characteristic stellar surface den-
sity observed at these radii in massive spheroids (see Hopkins et al.
2009e, and references therein).
One of the strong predictions of this work is that the link be-
tween high gas surface densities and inflow leads to a correlation
between BH accretion rate and star formation (§ 6). We have as-
sumed a star formation law in which _ / 1:5, which corresponds
to a fixed star formation efficiency per local dynamical time. This
is supported by models of turbulence-regulated star formation, in
which the efficiency per dynamical time is a weak function of
ambient conditions such as the Mach number of the turbulence
(Krumholz & McKee 2005). With this assumption, we predict a
correlation between BH accretion and the star formation at various
scales in galactic nuclei, albeit one with significant scatter (Fig-
ure 14). These predictions agree reasonably well with current ob-
servations (Wang et al. 2007; Wu & Cao 2006; Silverman et al.
2008). In the future, it will be interesting to explore the possibility
that the star formation efficiency in galactic nuclei is significantly
lower than the mean ISM value, as has been suggested by several
authors (Thompson et al. 2005; Begelman & Shlosman 2009; Lar-
son 2009).
If the asymmetric stellar structures found here are long-lived
(decay times of  Gyr), they should in principle be observable
around nearby massive BHs. One tantalizing possibility is that the
ubiquitous eccentric stellar and gaseous disks we find at  1 10
pc may explain the nuclear stellar disk seen in M31, which has
been the subject of considerable study (see e.g. Lauer et al. 1993;
Tremaine 1995; Salow & Statler 2001; Sambhus & Sridhar 2002;
Peiris & Tremaine 2003; Bender et al. 2005). Hopkins & Quataert
c
0000 RAS, MNRAS 000, 000–000

28 Hopkins and Quataert
(2010) examine this possibility in detail and show that the spa-
tial scales  1 10pc, moderate eccentricities, and stellar masses
 0:11MBH found in our simulations are all comparable to those
observed. Eccentric nuclear disks have also been observed (albeit
in less detail) around a number of other massive but inactive BHs,
including NGC4486b (Lauer et al. 1996), M83 (Thatte et al. 2000),
and VCC128 (Debattista et al. 2006). This strongly suggests that
we can probe the the instabilities that produced the growth of mas-
sive BHs in the “fossil” relics of that accretion.
The dominant uncertainties in the models presented here are
our treatment of the ISM physics, star formation, and feedback.
Lacking a full micro-physical understanding of these processes – or
the resolution to directly incorporate such a model – it is necessary
to adopt a sub-resolution prescription for their impact. Our model
is that stars form at all radii, in local, self-gravitating clumps, with
a fixed efficiency relative to the local dynamical time. In addition,
feedback from stars (supernovae, radiation, stellar winds/outflows,
etc.) contributes, potentially along with cosmic rays and magnetic
fields, to generating a turbulent velocity that is much larger than the
thermal sound speed of the ISM. It is this turbulent “sound speed”
that is included in our ISM equation of state (Figure 1). Observa-
tions suggest that these assumptions are at least plausible even at
the small scales we model here: a local Kennicutt-Schmidt relation
holds in nuclear starbursts and on small scales in galaxies (Kenni-
cutt 1998; Kennicutt et al. 2007; Bigiel et al. 2008) and star forma-
tion efficiencies per dynamical time do not appear to evolve even in
very dense regions (comparable to the highest gas densities mod-
eled here; see, e.g., Tan et al. 2006; Krumholz & Tan 2007). More-
over, AGN are observed to have associated nuclear star formation
on scales as small as a pc; the densities of stars formed in situ at
these scales (Davies et al. 2007; Hicks et al. 2009) are reasonably
consistent with our predictions. Large turbulent and/or non-thermal
gas velocity dispersions are also ubiquitous, from starburst nuclei
on kpc scales to the sub-pc scales around AGN probed by water
masers (Downes & Solomon 1998; Westmoquette et al. 2007; Tac-
coni et al. 1999; Iono et al. 2007; Kondratko et al. 2005); these ob-
servations motivate our choice of sub-resolution feedback models
(Figure 1). To study the impact of ISM uncertainties on our results,
we explicitly used different star formation efficiencies in some of
our simulations; although this does have a quantitative affect on the
resulting gas densities and inflow rates, the physical processes that
we have identified as driving accretion remain the same.
As noted in §2, the primary consequence of the large effective
sound speed in our models is to increase the Jeans and Toomre
masses, thus suppressing the formation of small-scale structure.
That is, we effectively smooth over the smaller scales on which
our treatment of the ISM physics is not appropriate. This implies
that in our model, most of the gas, most of the time, resides in rel-
atively diffuse structures, rather than being bound to very compact
clusters, as would be the case if we did not include stellar feed-
back and/or a sub-grid sound speed (see Appendix B). We suspect
that this approach is reasonable for the global gravitational dynam-
ics highlighted in this paper. In particular, our calculations indicate
that over a remarkably large dynamic range, from 0:1pc to 10kpc,
the torques and angular momentum transport are gravitational, and
scale simply with the magnitude of the asymmetry in the potential,
not with the sound speed of the gas (Fig. 9 and §4). In Paper II, we
will present a more detailed comparison between analytic models
and our numerical results that supports this conclusion.
In our calculations, we have neglected AGN feedback in order
to study how gas gets down to a BH in the first place. To the ex-
tent that feedback is important on these scales, our simulations may
better approximate the “early” stages of BH growth, when the BH
is still relatively small (compared to, e.g., the still-forming bulge),
and before it reaches a critical mass or luminosity at which feed-
back becomes dynamically important. When BH growth becomes
self-regulated, inflow and outflow are coupled, and the problem of
AGN fueling must be addressed with a model that includes AGN
feedback as a function of gas properties and spatial scale.
ACKNOWLEDGMENTS
We thank Phil Chang, Lars Hernquist, Norm Murray, and
Volker Springel for helpful discussions during the development of
this work, and the referee for a careful reading of the text and a
number of very useful suggestions. PFH and EQ were supported
by the Miller Institute for Basic Research in Science, University
of California Berkeley. EQ was also supported in part by NASA
grant NNG06GI68G and the David and Lucile Packard Foundation.
REFERENCES
Adams, F. C., Ruden, S. P., & Shu, F. H. 1989, ApJ, 347, 959
Allen, L., et al. 2007, in Protostars and Planets V, ed. B. Reipurth,
D. Jewitt, & K. Keil, 361–376
Aller, M. C., & Richstone, D. O. 2007, ApJ, 665, 120
Aoki, S., Noguchi, M., & Iye, M. 1979, PASJ, 31, 737
Athanassoula, E. 2008, MNRAS, 390, L69
Athanassoula, E., Bienayme, O., Martinet, L., & Pfenniger, D.
1983, A&A, 127, 349
Athanassoula, E., Lambert, J. C., & Dehnen, W. 2005, MNRAS,
363, 496
Balbus, S. A., & Hawley, J. F. 1998, Reviews of Modern Physics,
70, 1
Banerjee, R., Pudritz, R. E., & Holmes, L. 2004, MNRAS, 355,
248
Bardeen, J. M. 1975, in IAU Symposium, Vol. 69, Dynamics of
the Solar Systems, ed. A. Hayli, 297–+
Barnes, J. E. 1998, in Saas-Fee Advanced Course 26: Galax-
ies: Interactions and Induced Star Formation; Springer-Verlag
Berlin/Heidelberg, ed. R. C. Kennicutt, Jr., F. Schweizer, J. E.
Barnes, D. Friedli, L. Martinet, & D. Pfenniger, 275–+
Barnes, J. E., & Hernquist, L. 1996, ApJ, 471, 115
Barnes, J. E., & Hernquist, L. E. 1991, ApJL, 370, L65
Begelman, M., & Shlosman, I. 2009, ApJL, in press,
arXiv:0904.4247, 0904, 4247
Bell, E. F., & de Jong, R. S. 2001, ApJ, 550, 212
Bender, R., et al. 2005, ApJ, 631, 280
Berentzen, I., Shlosman, I., Martinez-Valpuesta, I., & Heller,
C. H. 2007, ApJ, 666, 189
Bigiel, F., Leroy, A., Walter, F., Brinks, E., de Blok, W. J. G.,
Madore, B., & Thornley, M. D. 2008, AJ, 136, 2846
Blitz, L. 1993, in Protostars and Planets III (University of Arizona
Press: Tucson), ed. E. H. Levy & J. I. Lunine, 125–161
Block, D. L., Puerari, I., Knapen, J. H., Elmegreen, B. G., Buta,
R., Stedman, S., & Elmegreen, D. M. 2001, Astronomy and As-
trophysics, 375, 761
Böker, T., Sarzi, M., McLaughlin, D. E., van der Marel, R. P., Rix,
H.-W., Ho, L. C., & Shields, J. C. 2004, AJ, 127, 105
Boley, A. C., Mejía, A. C., Durisen, R. H., Cai, K., Pickett, M. K.,
& D’Alessio, P. 2006, ApJ, 651, 517
Bonnell, I. A., Dobbs, C. L., Robitaille, T. P., & Pringle, J. E.
2006, MNRAS, 365, 37
c
0000 RAS, MNRAS 000, 000–000

How Do Massive Black Holes Get Their Gas? 29
Bouché, N., et al. 2007, ApJ, 671, 303
Bournaud, F., & Combes, F. 2002, A&A, 392, 83
Bournaud, F., Combes, F., Jog, C. J., & Puerari, I. 2005, A&A,
438, 507
Bryant, P. M., & Scoville, N. Z. 1999, AJ, 117, 2632
Bullock, J. S., et al. 2001, MNRAS, 321, 559
Cai, K., Durisen, R. H., Boley, A. C., Pickett, M. K., & Mejía,
A. C. 2008, ApJ, 673, 1138
Carlberg, R. G. 1986, ApJ, 310, 593
Chang, P., Murray-Clay, R., Chiang, E., & Quataert, E. 2007, The
Astrophysical Journal, 668, 236
Christodoulou, D. M., Shlosman, I., & Tohline, J. E. 1995, ApJ,
443, 551
Cossins, P., Lodato, G., & Clarke, C. J. 2009, MNRAS, 393, 1157
Côté, P., et al. 2006, ApJS, 165, 57
Covington, M., Dekel, A., Cox, T. J., Jonsson, P., & Primack, J. R.
2008, MNRAS, 384, 94
Cox, T. J., Dutta, S. N., Di Matteo, T., Hernquist, L., Hopkins,
P. F., Robertson, B., & Springel, V. 2006, ApJ, 650, 791
Cuadra, J., Armitage, P. J., Alexander, R. D., & Begelman, M. C.
2009, Monthly Notices of the Royal Astronomical Society, 393,
1423
Dasyra, K. M., et al. 2006, ApJ, 638, 745
—. 2007, ApJ, 657, 102
Davies, R. I., Sánchez, F. M., Genzel, R., Tacconi, L. J., Hicks,
E. K. S., Friedrich, S., & Sternberg, A. 2007, ApJ, 671, 1388
Debattista, V. P., Ferreras, I., Pasquali, A., Seth, A., De Rijcke, S.,
& Morelli, L. 2006, ApJL, 651, L97
DeBuhr, J., Quataert, E., Ma, C.-P., & Hopkins, P. 2009, MNRAS,
in press, arXiv:0909.2872
Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433,
604
Dotti, M., Ruszkowski, M., Paredi, L., Colpi, M., Volonteri, M.,
& Haardt, F. 2009, MNRAS, 396, 1640
Downes, D., & Solomon, P. M. 1998, ApJ, 507, 615
Dubinski, J., Berentzen, I., & Shlosman, I. 2009, ApJ, 697, 293
Durbala, A., Buta, R., Sulentic, J. W., & Verdes-Montenegro, L.
2009, MNRAS, in press [arXiv:0905.2340]
Earn, D. J. D., & Sellwood, J. A. 1995, ApJ, 451, 533
Elmegreen, B. G. 2007, ApJ, 668, 1064
Escala, A. 2007, ApJ, 671, 1264
Escala, A., Larson, R. B., Coppi, P. S., & Mardones, D. 2004, ApJ,
607, 765
Evans, A. S., Solomon, P. M., Tacconi, L. J., Vavilkin, T., &
Downes, D. 2006, AJ, 132, 2398
Ferrarese, L., & Merritt, D. 2000, ApJL, 539, L9
Ferrarese, L., et al. 2006a, ApJL, 644, L21
—. 2006b, ApJS, 164, 334
Förster Schreiber, N. M., et al. 2006, ApJ, 645, 1062
Gammie, C. F. 2001, ApJ, 553, 174
García-Burillo, S., Combes, F., Schinnerer, E., Boone, F., & Hunt,
L. K. 2005, Astronomy and Astrophysics, 441, 1011
Gebhardt, K., et al. 2000, ApJL, 539, L13
Goodman, J. 2003, MNRAS, 339, 937
Greene, J. E., Zakamska, N. L., Liu, X., Barth, A. J., & Ho, L. C.
2009, ApJ, in press, arXiv:0907.1086
Gunn, J. E. 1987, in Nearly Normal Galaxies. From the Planck
Time to the Present (New York: Springer-Verlag), ed. S. M.
Faber, 455–464
Haan, S., Schinnerer, E., Emsellem, E., García-Burillo, S.,
Combes, F., Mundell, C. G., & Rix, H.-W. 2009, The Astrophys-
ical Journal, 692, 1623
Heller, C., Shlosman, I., & Englmaier, P. 2001, The Astrophysical
Journal, 553, 661
Hernquist, L. 1989, Nature, 340, 687
—. 1990, ApJ, 356, 359
—. 1993, ApJ, 404, 717
Hernquist, L., & Ostriker, J. P. 1992, Astrophysical Journal, 386,
375
Hernquist, L., Spergel, D. N., & Heyl, J. S. 1993, ApJ, 416, 415
Hibbard, J. E., & Yun, M. S. 1999, ApJL, 522, L93
Hicks, E. K. S., Davies, R. I., Malkan, M. A., Genzel, R., Tacconi,
L. J., Sánchez, F. M., & Sternberg, A. 2009, ApJ, 696, 448
Ho, L. C. 2005, ApJ, 629, 680
Hopkins, P. F., Cox, T. J., Dutta, S. N., Hernquist, L., Kormendy,
J., & Lauer, T. R. 2009a, ApJS, 181, 135
Hopkins, P. F., Cox, T. J., & Hernquist, L. 2008a, ApJ, 689, 17
Hopkins, P. F., Cox, T. J., Younger, J. D., & Hernquist, L. 2009b,
ApJ, 691, 1168
Hopkins, P. F., & Hernquist, L. 2006, ApJS, 166, 1
—. 2009, MNRAS, in press, arXiv:0910.4582
Hopkins, P. F., Hernquist, L., Cox, T. J., Dutta, S. N., & Rothberg,
B. 2008b, ApJ, 679, 156
Hopkins, P. F., Hernquist, L., Cox, T. J., Kereš, D., & Wuyts, S.
2009c, ApJ, 691, 1424
Hopkins, P. F., Hernquist, L., Cox, T. J., Robertson, B., & Krause,
E. 2007, ApJ, 669, 67
Hopkins, P. F., Hernquist, L., Martini, P., Cox, T. J., Robertson,
B., Di Matteo, T., & Springel, V. 2005, ApJL, 625, L71
Hopkins, P. F., Lauer, T. R., Cox, T. J., Hernquist, L., & Kor-
mendy, J. 2009d, ApJS, 181, 486
Hopkins, P. F., Murray, N., Quataert, E., & Thompson, T. A.
2009e, MNRAS, in press, arXiv:0908.4088
Hopkins, P. F., & Quataert, E. 2010, MNRAS, L62+
Hopkins, P. F., Somerville, R. S., Cox, T. J., Hernquist, L., Jogee,
S., Kereš, D., Ma, C.-P., Robertson, B., & Stewart, K. 2009f,
MNRAS, 397, 802
Imanishi, M. 2002, ApJ, 569, 44
—. 2003, ApJ, 599, 918
Imanishi, M., & Wada, K. 2004, ApJ, 617, 214
Iono, D., et al. 2007, ApJ, 659, 283
Jogee, S. 2006, in Lecture Notes in Physics, Berlin Springer Ver-
lag, Vol. 693, Physics of Active Galactic Nuclei at all Scales, ed.
D. Alloin, 143–+
Jogee, S., Scoville, N., & Kenney, J. D. P. 2005, ApJ, 630, 837
Johansson, P. H., Naab, T., & Burkert, A. 2009a, ApJ, 690, 802
—. 2009b, ApJL, in press, arXiv:0910.2232
Kalnajs, A. J. 1971, ApJ, 166, 275
Kauffmann, G., et al. 2003, MNRAS, 346, 1055
Kawakatu, N., & Wada, K. 2008, ApJ, 681, 73
Kennicutt, Jr., R. C. 1998, ApJ, 498, 541
Kennicutt, Jr., R. C., et al. 2007, ApJ, 671, 333
Kewley, L. J., Groves, B., Kauffmann, G., & Heckman, T. 2006,
MNRAS, 372, 961
Kim, M., Ho, L. C., & Im, M. 2006, ApJ, 642, 702
Knapen, J. H., Pérez-Ramírez, D., & Laine, S. 2002, Monthly No-
tice of the Royal Astronomical Society, 337, 808
Knapen, J. H., Shlosman, I., & Peletier, R. F. 2000, ApJ, 529, 93
Kondratko, P. T., Greenhill, L. J., & Moran, J. M. 2005, ApJ, 618,
618
Kormendy, J. 1989, ApJL, 342, L63
Kormendy, J., Fisher, D. B., Cornell, M. E., & Bender, R. 2009,
ApJS, 182, 216
Kormendy, J., & Freeman, K. C. 2004, in IAU Symposium,
c
0000 RAS, MNRAS 000, 000–000

30 Hopkins and Quataert
Vol. 220, Dark Matter in Galaxies, ed. S. Ryder, D. Pisano,
M. Walker, & K. Freeman, 377–+
Kormendy, J., & Richstone, D. 1995, ARA&A, 33, 581
Krips, M., et al. 2007, A&A, 468, L63
Krumholz, M. R., Klein, R. I., & McKee, C. F. 2007, ApJ, 656,
959
Krumholz, M. R., Matzner, C. D., & McKee, C. F. 2006, ApJ,
653, 361
Krumholz, M. R., & McKee, C. F. 2005, ApJ, 630, 250
Krumholz, M. R., & Tan, J. C. 2007, ApJ, 654, 304
Laine, S., Shlosman, I., Knapen, J. H., & Peletier, R. F. 2002, The
Astrophysical Journal, 567, 97
Laine, S., van der Marel, R. P., Rossa, J., Hibbard, J. E., Mihos,
J. C., Böker, T., & Zabludoff, A. I. 2003, AJ, 126, 2717
Larson, R. B. 1981, MNRAS, 194, 809
—. 2009, Nature, in press, arXiv:0909.3110
Lau, Y. Y., & Bertin, G. 1978, ApJ, 226, 508
Lauer, T. R., et al. 1993, AJ, 106, 1436
—. 1996, ApJL, 471, L79+
—. 2007, ApJ, 664, 226
Lemoine-Busserolle, M., Bunker, A., Lamareille, F., & Kissler-
Patig, M. 2009, MNRAS, in press, arXiv:0909.1386
Levine, R., Gnedin, N. Y., Hamilton, A. J. S., & Kravtsov, A. V.
2008, ApJ, 678, 154
Lodato, G., & Rice, W. K. M. 2004, MNRAS, 351, 630
Lynden-Bell, D., & Kalnajs, A. J. 1972, MNRAS, 157, 1
Magorrian, J., et al. 1998, AJ, 115, 2285
Mandelbaum, R., Seljak, U., Kauffmann, G., Hirata, C. M., &
Brinkmann, J. 2006, MNRAS, 368, 715
Martini, P. 2004, in Coevolution of Black Holes and Galaxies
(Cambridge: Cambridge University Press), ed. L. C. Ho, 169
Martini, P., & Pogge, R. W. 1999, Astron. J., 118, 2646
Martini, P., Regan, M. W., Mulchaey, J. S., & Pogge, R. W. 2003,
ApJ, 589, 774
Mayer, L., Kazantzidis, S., Madau, P., Colpi, M., Quinn, T., &
Wadsley, J. 2007, Science, 316, 1874
McKee, C. F., & Ostriker, J. P. 1977, ApJ, 218, 148
Mihos, J. C., & Hernquist, L. 1994, ApJL, 431, L9
—. 1996, ApJ, 464, 641
Murray, N., Quataert, E., & Thompson, T. A. 2009, ApJ, in press
[arXiv:0906.5358]
Narayan, R., Goldreich, P., & Goodman, J. 1987, MNRAS, 228, 1
Nayakshin, S., Cuadra, J., & Springel, V. 2007, MNRAS, 379, 21
Nayakshin, S., & King, A. 2007, MNRAS, in press,
arXiv:0705.1686
Noguchi, M. 1988, A&A, 203, 259
O’Shea, B. W., Nagamine, K., Springel, V., Hernquist, L., & Nor-
man, M. L. 2005, ApJS, 160, 1
Ostriker, E. C., Shu, F. H., & Adams, F. C. 1992, ApJ, 399, 192
Ostriker, J. P. 1980, Comments on Astrophysics, 8, 177
Peiris, H. V., & Tremaine, S. 2003, ApJ, 599, 237
Peletier, R. F., et al. 1999, ApJS, 125, 363
Perets, H. B., & Alexander, T. 2008, ApJ, 677, 146
Perets, H. B., Hopman, C., & Alexander, T. 2007, ApJ, 656, 709
Quillen, A. C., Frogel, J. A., Kenney, J. D. P., Pogge, R. W., &
Depoy, D. L. 1995, ApJ, 441, 549
Raha, N., Sellwood, J. A., James, R. A., & Kahn, F. D. 1991,
Nature, 352, 411
Rice, W. K. M., Lodato, G., & Armitage, P. J. 2005, MNRAS,
364, L56
Riechers, D. A., Walter, F., Bertoldi, F., Carilli, C. L., Aravena,
M., Neri, R., Cox, P., Weiss, A., & Menten, K. M. 2009, ApJ, in
press, arXiv:0908.0018
Robertson, B., Bullock, J. S., Cox, T. J., Di Matteo, T., Hernquist,
L., Springel, V., & Yoshida, N. 2006a, ApJ, 645, 986
Robertson, B., Hernquist, L., Cox, T. J., Di Matteo, T., Hopkins,
P. F., Martini, P., & Springel, V. 2006b, ApJ, 641, 90
Sakamoto, K., Okumura, S. K., Ishizuki, S., & Scoville, N. Z.
1999, The Astrophysical Journal, 525, 691
Salow, R. M., & Statler, T. S. 2001, ApJL, 551, L49
Salucci, P., Szuszkiewicz, E., Monaco, P., & Danese, L. 1999,
MNRAS, 307, 637
Sambhus, N., & Sridhar, S. 2002, A&A, 388, 766
Sanders, D. B., & Mirabel, I. F. 1996, ARA&A, 34, 749
Sanders, D. B., Soifer, B. T., Elias, J. H., Madore, B. F., Matthews,
K., Neugebauer, G., & Scoville, N. Z. 1988a, ApJ, 325, 74
Sanders, D. B., Soifer, B. T., Elias, J. H., Neugebauer, G., &
Matthews, K. 1988b, ApJL, 328, L35
Satyapal, S., Dudik, R. P., O’Halloran, B., & Gliozzi, M. 2005,
ApJ, 633, 86
Schartmann, M., Meisenheimer, K., Klahr, H., Camenzind, M.,
Wolf, S., & Henning, T. 2009, MNRAS, 393, 759
Schawinski, K., et al. 2009, eprint arXiv, 0903, 3415
Schweizer, F. 1998, in Saas-Fee Advanced Course 26: Galax-
ies: Interactions and Induced Star Formation (Springer-Verlag
Berlin/Heidelberg), ed. R. C. Kennicutt, Jr., F. Schweizer, J. E.
Barnes, D. Friedli, L. Martinet, & D. Pfenniger, 105–+
Schweizer, F., & Seitzer, P. 2007, AJ, 133, 2132
Shankar, F., Salucci, P., Granato, G. L., De Zotti, G., & Danese,
L. 2004, MNRAS, 354, 1020
Shlosman, I., & Begelman, M. C. 1989, ApJ, 341, 685
Shlosman, I., Begelman, M. C., & Frank, J. 1990, Nature, 345,
679
Shlosman, I., Frank, J., & Begelman, M. C. 1989, Nature, 338, 45
Shu, F. H., Tremaine, S., Adams, F. C., & Ruden, S. P. 1990, ApJ,
358, 495
Silverman, J. D., et al. 2008, eprint arXiv, 0810, 3653
Soltan, A. 1982, MNRAS, 200, 115
Springel, V. 2005, MNRAS, 364, 1105
Springel, V., Di Matteo, T., & Hernquist, L. 2005a, ApJL, 620,
L79
—. 2005b, MNRAS, 361, 776
Springel, V., & Hernquist, L. 2002, MNRAS, 333, 649
—. 2003, MNRAS, 339, 289
—. 2005, ApJL, 622, L9
Stoklasová, I., Ferruit, P., Emsellem, E., Jungwiert, B., Pécontal,
E., & Sánchez, S. F. 2009, A&A, in press [arXiv:0905.3349]
Tacconi, L. J., Genzel, R., Tecza, M., Gallimore, J. F., Downes,
D., & Scoville, N. Z. 1999, ApJ, 524, 732
Tacconi, L. J., et al. 2006, ApJ, 640, 228
Tan, J. C., Krumholz, M. R., & McKee, C. F. 2006, ApJL, 641,
L121
Thatte, N., Tecza, M., & Genzel, R. 2000, A&A, 364, L47
Thompson, T. A., Quataert, E., & Murray, N. 2005, ApJ, 630, 167
Toomre, A. 1977, ARA&A, 15, 437
Tremaine, S. 1995, AJ, 110, 628
—. 2001, AJ, 121, 1776
Veilleux, S., et al. 2009, ApJ, 701, 587
Vorobyov, E. I., & Basu, S. 2006, ApJ, 650, 956
—. 2009, MNRAS, 393, 822
Wada, K., & Norman, C. A. 2002, ApJL, 566, L21
Wada, K., Papadopoulos, P., & Spaans, M. 2009, ApJ, in press,
arXiv:0906.5444
c
0000 RAS, MNRAS 000, 000–000

How Do Massive Black Holes Get Their Gas? 31
N = 5x105 N = 3x106 N = 2x107
Figure A1. Example resolution test for a nuclear scale simulation. The
same initial conditions re-simulated with the particle number shown (lower-
resolution and nominal cases of Nf8h2b1h; see Table 3). For & 0:5 106
particles, our results are converged.
Walter, F., Riechers, D., Cox, P., Neri, R., Carilli, C., Bertoldi, F.,
Weiss, A., & Maiolino, R. 2009, Nature, 457, 699
Wang, J.-M., Chen, Y.-M., Yan, C.-S., Hu, C., & Bian, W.-H.
2007, ApJL, 661, L143
Wang, R., et al. 2008, ApJ, 687, 848
Westmoquette, M. S., Smith, L. J., Gallagher, III, J. S., O’Connell,
R. W., Rosario, D. J., & de Grijs, R. 2007, ApJ, 671, 358
Wise, J. H., Turk, M. J., & Abel, T. 2007, ArXiv e-prints, 710
Woods, D. F., Geller, M. J., & Barton, E. J. 2006, AJ, 132, 197
Wu, Q., & Cao, X. 2006, The Publications of the Astronomical
Society of the Pacific, 118, 1098
Younger, J. D., Hopkins, P. F., Cox, T. J., & Hernquist, L. 2008,
ApJ, 686, 815
APPENDIX A: NUMERICAL TESTS
We have performed a number of tests of our methodology, in an
effort to ensure both that we have captured the most important
physics in our simulations, and to ensure that there are no artifi-
cial effects occurring due to numerical or resolution artifacts. We
describe some of the key tests here. For general tests of GADGET-
3, we refer the reader to Springel (2005); Springel & Hernquist
(2003); Springel et al. (2005b); Springel & Hernquist (2002).
A1 Resolution Tests
For the initial galaxy-scale simulations (Table 1), resolution tests
have been extensively described previously (Di Matteo et al. 2005;
Robertson et al. 2006b; Cox et al. 2006; Hopkins et al. 2009a);
these range from  105 particles to  107 particles, and spatial res-
olutions/softenings from  20pc to  100pc. For our purposes,
the quantities of interest (e.g. the gas inflows to  100pc scales)
are well-converged (see e.g. Appendix B of Hopkins et al. 2009a).
We have performed a series of analogous resolution tests on
our intermediate (Table 2) and small-scale simulations (Table 3),
covering a similar span in particle number. For the intermediate-
scale simulations, we have varied the SPH smoothing length and
gravitational softening length from < 1pc to  10pc. We find that
the results of particular importance for this paper, e.g., the total
gas mass that loses angular momentum and flows to smaller scales,
are well-converged for all the resolutions . 10pc used here. This
is not surprising, because for physical reasons (discussed in the
main text), the material reaches an angular momentum barrier at
 10pc. For our small-scale simulations, we find similar behav-
ior for all force resolutions 1pc (experimenting with a range of
softenings from 0:021pc), and with good convergence at higher
particle numbers & a few 105. Given the range of resolution ex-
plored, comparison with Figure 10 demonstrates that, especially
in our high-resolution studies, even the verticle structure of the
gaseous disks is well-resolved (in the highest-resolution cases, the
typical hsml=R< 0:01). Figure A1 shows an example resolution test
varying from 51052107 particles.
A2 Do the Results Depend on Large-scale Tidal Fields?
Because we are simulating the properties of large-scale non-
axisymmetric systems, one might worry that the tidal forces from
large radii could be important. This is certainly the case, e.g., in
cosmological simulations. In the present context, the question is
whether we need to include the matter distribution at  10kpc in
real-time to understand the dynamics of material at . 100pc.
We assess this directly as follows. For a given re-simulation
of the nuclear regions, we determine the matter distribution of the
larger-scale simulation from which the re-simulation is drawn and
analytically fit the non-axisymmetric contributions to the potential
at all radii.10 Figure A2 compares the results from a re-simulation
in which we include the large-scale non-axisymmetric potential at
& kpc with our standard approach, in which we do not include this
potential and in which the matter at large radii is not included in the
re-simulation (so the information is lost). This is simulation If9b5
in Table 2, a  0:01 1 kpc re-simulation of the galactic nucleus
during a galaxy-galaxy merger simulation at the peak of nuclear
activity, just after the coalescence of the two galactic nuclei. At
large scales, tidal tails are still present (these will take up to sev-
eral Gyr to fully relax), and the system is largely unrelaxed outside
of a few hundred pc, so this is a time when the large-scale effects
would a priori be the most important. In the simulations shown in
Figure A2 (left panel), we analytically fit the large-scale tidal field
in the galaxy-galaxy merger simulation at multiple times, and in-
terpolate it in time in our re-simulation.
Figure A2 shows that there is effectively no difference be-
tween the simulations with and without the large-scale tidal field, at
any of the radii of interest. Including large-scale tidal forces does
somewhat truncate the large spiral wave pattern that appears as a
secondary instability in the intermediate-scale simulation, but this
truncation occurs only on the largest radii & 500pc, radii that we
are not trying to model in detail in the re-simulation. Most impor-
tantly, the scales which the re-simulation is intended to accurately
represent, . 100pc, are indistinguishable in the two simulations.
The reason for this is ultimately that the mass profile is a relatively
steep function of radius. If a mass M is enclosed in an annulus at
a radius r, then the tidal force it produces is / jajM=r3, where jaj
is the magnitude of the non-axisymmetric component of the mass
distribution at that annulus. The contribution of a given large radius
to the tidal force at small radii is thus / (r), the local mass den-
sity. This, however, is always a decreasing function of (increasing)
radius, so that material at larger radii makes little contribution to
the tidal force. This can be compared with the role of tidal forces
in determining the angular momentum of large-scale dark matter
structures and ultimately (via these structures) dark matter halos;
since the large-scale matter field has constant on average, there
are significant contributions to tidal forces from all radii. This is
why including the information about large-scale tidal fields is crit-
ical for simulations of cooling and galaxy formation, but it is not
important when following the further inflows of material in indi-
vidual galaxies.
10 To do this, we adopt the radial basis expansion set proposed in Hernquist
& Ostriker (1992), with standard spherical harmonics. This is similar to the
density modeling described in § 2, and can be performed to arbitrary order,
but only the first few terms are not noise-dominated.
c
0000 RAS, MNRAS 000, 000–000

32 Hopkins and Quataert








With Large-Scale Fields
A








500 pc
Without Large-Scale Fields
B 10 100 1000R [pc]
9.0
9.5
10.0
10.5
11.0
log
(Σ ga
s) [M
O • k
pc-
2 ]
AB
Figure A2. Simulations testing the effects of large-scale tidal forces in our re-simulations. Left: Re-simulation of the central  10pc-kpc of a larger-scale
galaxy merger simulation, in which the full time-dependent potential from the large-scale simulation is analytically included in the re-simulation. We show the
gas density and star formation rate density as in Figure 3, three large-scale dynamical times after the beginning of the simulation. Center: Same re-simulation,
but without the additional potential from the large-scale material; i.e., the gravitational potential is only from the “live” material being re-simulated in the
central kpc. The image is at the same point in time. Right: Gas density profiles of both re-simulations at the peak of the resulting inflows. Because the density
declines rapidly with radius, the large-scale density/potential fields make no significant difference to the dynamics at the these scales.
We have performed similar experiments to those shown here
for other intermediate-scale simulations (both in mergers and iso-
lated galaxies), and for our smallest-scale simulations ( 0:1
10pc), and reach the same conclusions in each case.
A3 Do the Results Depend on Initial Conditions?
A reasonable question about our re-simulation technique is whether
it might be sensitive to the precise initial conditions in the nuclear
disk that is initialized at high resolution. In particular, there is lit-
tle detailed information about the structure or dynamics of the nu-
clear disk in the original larger-scale simulation, precisely because
it lacked the necessary resolution to study the small-scale nuclear
dynamics. Here we show that our results do not depend on the de-
tails of how we initialize the re-simulations in the nuclear region.
The fundamental reason for this is that the nuclear dynamics is gov-
erned by instabilities that are self-consistently generated on small
scales. These instabilities depend primarily on the gross structural
properties of the nucleus, and rapidly lose memory of the initial
conditions.
Figure A3 shows the results of five different re-simulations,
each with different initial conditions. These are all  100pc re-
simulations of the central regions from a galaxy-merger simula-
tion (intermediate-scale simulation If9b5 in Table 2), just after the
coalescence of the two galaxy nuclei. Panel A is a rigorous re-
simulation of the matter distribution in the central kpc of the large-
scale simulation, with the matter initialized “as-is” from the large-
scale simulation. Panel B is the same except we do not include
any of the non-axisymmetric modes in the density, pressure, po-
tential, etc. in the initial conditions; i.e., we azimuthally average
the initial conditions used in Panel A. Panel C is even more ide-
alized: the system is axisymmetric and with exponential disk and
Hernquist profiles for the “disky” and “spherical” components in
the re-simulation, respectively. The masses and effective radii of
these components match the original simulation, but the profiles
are these simple analytic approximations rather than the true pro-
file from the larger-scale simulation. To test whether the presence of
particular non-axisymmetries is important, Panel D includes non-
axisymmetric structure in the initial conditions, but the precise
modes are randomly chosen in phase and amplitude and so do not
match the actual non-axisymmetries from Panel A. Finally, all of
the above simulations include gas at small radii in the initial condi-
tions; this is because the finite smoothing length in the larger-scale
simulation spreads gas out over a region  100 pc even though the
dynamics below this scale is completely incorrect. In order to be as
conservative as possible, Panel E considers an initial “hole” in the
gas distribution at radii that are not resolved in the large scale simu-
lation; for numerical reasons, the “hole” is a sharp power-law cutoff
inside 100pc. This could represent a very hard angular momentum
barrier, such as a strong inner Linblad resonance.
The images and surface density plots in Figure A3 are 108 yr
after the beginning of the simulation. There is remarkably little dif-
ference in the mode structure visible in the images. More quantita-
tively, the surface density profiles are very similar, with significant
gas inflow to  10 pc. This demonstrates that the precise initial
conditions used in the re-simulations are not important for most
of our results. The fundamental reason for this is that the instabil-
ities that dominate the dynamics arise from the internal structure
of the material; they depend on the mass and scale-lengths of the
rotationally and dispersion supported components, but not on the
initial seed amplitudes of various modes. In the case of axisymmet-
ric initial conditions (Panel B), the initial perturbations take slightly
longer to appear because the initial power present to be amplified
is smaller (it is due to particle noise); and in the case of strong ini-
tial modes (Panel D), the instability is slightly more developed at
the fixed time shown here. However, these time shifts have little
long-term effect on the evolution of the system.
The robustness of our conclusions to resolution, tidal fields,
and initial conditions is also demonstrated by the fact that our ultra-
high resolution simulations, discussed at length in the text, give
identical results to our standard re-simulated runs. These ultra-high
resolution simulations allow us to bridge the boundaries of our
standard re-simulations, and so do not suffer the potential nega-
tive effects considered here. Figure A4 shows an illustrative ex-
ample of one such simulation about mid-way through the period
of peak activity. The simulation is initialized as one of our typical
intermediate-scale re-simulations, but with  107 particles, giving
a final resolution of . pc. On large scales, the dynamics is similar
c
0000 RAS, MNRAS 000, 000–000

How Do Massive Black Holes Get Their Gas? 33
Figure A3. Effects of initial conditions on our re-simulations of galactic nuclei. Panels A-E show the gas density and star formation rate density 108 yr after
the beginning of the simulation. Bottom right: Surface density profiles of all five simulations at the same time. A: Re-simulation of the central kpc from a
larger-scale merger simulation (as in Figure A2). The initial conditions are taken exactly from the larger-scale simulation, i.e., with inherited non-axisymmetric
modes. B: Same as A but in this case we initialize azimuthally averaged (i.e., axisymmetric) profiles. C: Even more idealized: we initialize an axisymmetric
bulge+disk+halo+BH system with analytic profile shapes; the mass and radius of each component is the same as that of the material in the original large-scale
simulation at these radii, but the profile shapes are not. D: Same as A, but with random seed non-axisymmetric modes having order unity random amplitudes
and phases E: Same as A, but the gas at . 100pc is removed so that there is initially a central “hole” in the gas distribution. Comparison of A-E demonstrates
that the precise initial conditions have little effect on the mode dynamics and gas transport. Morphologically, systems with stronger (weaker) initial seeds
develop modes slightly more quickly (slowly), but this time shift has little significant or lasting effect on the dynamics.
to our other intermediate-scale runs, with an m = 2 mode forming
rapidly (there are non-trivial mode amplitudes from m = 1 3).
The inflow leads to large gas masses at small radii; by the active
phase shown in Figure A4, the lopsided, eccentric nuclear disk is
evident. The properties of this disk are statistically indistinguish-
able from those in our nuclear-scale re-simulations, given similar
boundary conditions (compare e.g. Figure 2).
APPENDIX B: ISM PHYSICS
As discussed extensively in §2 & §7, the largest a priori uncer-
tainty in our modeling is the question of what physics should be
included to describe the behavior of the ISM. Our approach is to
include a large “turbulent” sound speed in the equation state, as a
subgrid model for the effects of feedback by stellar winds, radia-
tion pressure, and supernovae on the ISM. We choose the turbulent
velocity largely by comparison to observed systems (Fig. 1). The
resulting turbulent velocities  30 100 kms1 are also reason-
ably consistent with theoretical estimates of the feedback required
to maintain Q  1 and disrupt molecular clouds in starburst envi-
ronments (Thompson et al. 2005; Murray et al. 2009).
To illustrate the importance of including a subgrid model, in
this Appendix we compare the results of our standard models with
a simulation in which the ISM is isothermal at  104 K, i.e., cs '
10 kms1 (qeos = 0). We have experimented with a range of isother-
mal cooling floors from cs = 1 15kms1 (100 3 104 K), but
the results here are generic to this range. Figure B1 shows sev-
eral of the key results comparing the cs ' 10 kms1 simulation
with our fiducial qeos = 0:125 simulation, for a re-simulation on
 100pc scales. From the images it is clear that after just a few lo-
cal dynamical times in the intermediate-scale disk, the simulation
without feedback has violently fragmented into large clumps. It is
important to stress that in this simulation we fully resolve the Jeans
and Toomre lengths so the fragmentation is physical, not numerical.
The difficulty with the cs ' 10 kms1 (qeos = 0) simulation
is that, absent a full model for how feedback can either suppress
their formation or suppress their internal star formation, the only
possibility is that the clumps “run away” and turn efficiently into
stars. The runaway and final consequences of the fragmentation for
star formation and the IMF are sensitive to details of the cooling
rates and cooling function, but the end result of runaway star for-
mation is inevitable if the gas is allowed to be arbitrarily cold and
has a small cooling time (Gammie 2001; Nayakshin et al. 2007;
Thompson et al. 2005). This has several dramatic effects. First, the
global star formation efficiency is significantly enhanced – for the
same mean gas surface density, most of the gas will turn into stars
in  1 4 local dynamical times absent feedback. This is a fac-
tor of 1050 higher than implied by the Schmidt-Kennicutt laws.
For comparison, our fiducial qeos = 0:125 simulation has a star for-
mation efficiency per dynamical time of a few percent, much more
consistent with observations. Because of the rapid and efficient star
formation in the cs = 10 kms1 simulation, the clumps quickly
convert a large fraction of their mass to stars; the resulting mas-
c
0000 RAS, MNRAS 000, 000–000

34 Hopkins and Quataert








100 pc








100 pc










10 pc








10 pc







10 pc










10 pc
Figure A4. Illustration of one of our ultra-high resolution simulations (Inf28b2h). The gas density distribution is shown face-on, as in Figure A3, in most
panels, and edge-on at bottom right. The initial conditions are those of a typical intermediate-scale re-simulation, but with 107 particles and sub-pc resolution.
This obviates the need for a secondary “re-simulation” of the nuclear scales. The large-scale modes fuel gas inwards and here lead to a large nuclear gas mass,
which then forms a lopsided m = 1 mode around the BH and drives further inflow. The eccentric nuclear gas disk is clearly visible, and similar to those in our
nuclear-scale re-simulations; in fact there is no statistical difference, given similar boundary conditions.
sive stellar clusters sink and coalesce via clump-clump scattering
and/or dynamical friction, to the center of the galaxy. As the lower-
right panel of Figure B1 shows, most of the mass flowing into the
central  10 pc in the cs ' 10 kms1 simulation is in the form of
stars, rather than gas. This results in the formation of an extremely
massive nuclear star cluster, with a mass of  109 M inside the
central  10 20pc – a factor of at least  10 larger than ob-
served nuclear star clusters (Böker et al. 2004; Côté et al. 2006;
Ferrarese et al. 2006a). Indeed, the resulting stellar surface den-
sity at small scales (left panel of Fig. B1) significantly exceeds the
highest nuclear densities observed in any cusp ellipticals (Ferrarese
et al. 2006b; Kormendy et al. 2009; Lauer et al. 2007), or for that
matter any nuclear star clusters, globular clusters, nuclear disks, or
other high-density stellar systems (Hopkins et al. 2009e).
In contrast to the constant cs ' 10 kms1 simulation, our fidu-
cial simulation with significant subgrid feedback does not clump
up into many dense star clusters. This is also not completely phys-
ical since in reality we expect that the ISM should have signifi-
cant sub-structure on the scales we model (including, e.g., molec-
ular clouds and star clusters). However, when such self-gravitating
gaseous clumps form, observations strongly suggest that are likely
to be short-lived and/or to inefficiently form stars. Our subgrid
model assumes that most of the mass remains in a diffuse ISM,
rather than becoming bound in dense star clusters. Not only is the
latter physically implausible, but Figure B1 demonstrates that it
strongly violates observational constraints. Note that there are con-
ditions under which even an extremely cold gaseous disk could
avoid catastrophic fragmentation (see e.g. Lodato & Rice 2004;
Rice et al. 2005; Boley et al. 2006; Krumholz et al. 2007), partic-
ularly the case in which the cooling time is sufficiently long such
that continued cooling is offset by gravitational heating (giving rise
to sustained local structure such as tightly-wound spiral arms, but
not runaway fragmentation). However, the cooling rates under typ-
ical conditions in our simulations are much too rapid to reach this
regime; we could, of course, modify the cooling function to sup-
press the cooling rate and/or mimic some continuous heating term
– but this is effectively equivalent to adopting a new sub-grid feed-
back/microphysics model, and we find it has the same qualitative
effects.
Figure B2 shows the effects of varying the efficiency of feed-
back from supernovae and massive stars, encapsulated in the pa-
rameter qeos. Specifically, we show results of a set of nuclear scale
simulations (Nf8h2b4q in Table 3) having qeos = 0:02 0:40 (ef-
fective turbulent velocities  20 100kms1), roughly the lower
and upper limits allowed by observational constraints for the sys-
tems of interest in Figure 1. Figure B2 shows that the choice of
sub-resolution model has a significant effect on the amount of re-
solved sub-structure in the simulation. This is not surprising since
c
0000 RAS, MNRAS 000, 000–000

How Do Massive Black Holes Get Their Gas? 35
-0.2 0.0 0.2
x [kpc]
-0.2
0.0
0.2
y
[kp
c]
qeos=0
-0.2 0.0 0.2
x [kpc]
-0.2
0.0
0.2
y
[kp
c]
qeos=0.12
1 10 100
R [pc]
109
1010
1011
1012
1013
1014
Σ n
ew
s
ta
rs


[M
O •
kp
c-2
]
qeos=0 (100 K)qeos=0 (104 K)qeos=0.125Observed
∆t = 107 yr
0 2 4 6 8 10
time [Myr]
1
10
100
dM
/dt
(<1
0p
c)
[M
O •
yr
-
1 ] StarsGas
Figure B1. Examples of the problems that occur if systems are evolved without some “feedback” to prevent runaway cooling and fragmentation. Top Left:
Gas density in such an initially smooth re-simulation of the central kpc (unlike in the main text, here darker areas represent higher-density to highlight the
fragmentation). The system here can cool to 104 K and has no feedback (qeos = 0); the image is shown after  2 106 years (less than the global dynamical
time). Top Right: The same initial conditions but for our fiducial model – a moderate subgrid (“turbulent”) sound speed  30kms1 (qeos = 0:125) in dense
star-forming regions. In the no feedback qeos = 0 case, gas clumps at the Jeans scale and rapidly turns into stars – leading to a severe violation of the observed
global Kennicutt (1998) relation. Bottom Left: Stellar density profiles after 107 yr, for simulations with no feedback (two examples shown, with different
cooling floors as labeled) and our fiducial qeos = 0:125 model; the initial density is just 1010 M kpc
2 at small radii. We also show the observed stellar mass
density profiles of massive cusp ellipticals in Virgo from Kormendy et al. (2009), chosen to have the same mass at > 0:5 1kpc. In the absence of feedback,
runaway fragmentation and the inability of clumps to dissolve inevitably leads to the formation of an extremely massive nuclear stellar cluster with a mass
and surface density at least an order of magnitude larger than any observed. Bottom Right: The accretion rate into the central 10pc for both gas (dotted) and
already-formed stars (solid). In the qeos = 0 models, sinking clumps provide large gas accretion rates 10100M yr
1, but the stellar inflow rate typically
exceeds the gas inflow rate by a factor of  10; this is inconsistent with the observed stellar densities at small radii. Our fiducial, moderate-feedback case, by
contrast, drives primarily gas – not stars – to small radii.
larger turbulent velocities raise the Jeans mass/length, above which
gravity is the dominant force. The key point, however, is that all of
the simulations show a similar nuclear lopsided disk. In terms of the
properties of inflow and global instability discussed throughout this
paper, the differences produced by changing qeos are similar to the
differences produced by somewhat different galaxy properties. The
fundamental reason for the weak dependence on the subgrid model
is that the torques in our simulations are primarily determined by
gravity, not hydrodynamic forces or viscosity. The primary role of
the subgrid feedback model is simply to prevent catastrophic frag-
mentation of the galactic gas. Provided this can be accomplished,
our qualitative conclusions do not depend on the details of the feed-
back model.
c
0000 RAS, MNRAS 000, 000–000

36 Hopkins and Quataert
qeos = 0.02 qeos = 0.06 qeos = 0.12
qeos = 0.21 qeos = 0.30 qeos = 0.40
10 pc
Figure B2. Effects of varying the assumed efficiency of stellar feedback on the resulting nuclear eccentric disk, given identical inflow conditions at 100pc.
We show the gas surface density with color indicating specific SFR (as Figure A3), for several choices of qeos (from the survey Nf8h2b4q). The range shown
corresponds to sub-grid turbulent velocities ranging from  20 100kms1 in the high star-formation rate regions of the ISM. The efficiency of feedback
affects the level of substructure. However, the general character and appearance of the lopsided disk mode is similar over the entire plausible range. Provided
the catastrophic fragmentation in Figure B1 can be avoided, our qualitative conclusions are independent of qeos.
c
0000 RAS, MNRAS 000, 000–000