Stellar Feedback, AGN Accretion, & Obscuration

Stellar and Quasar Feedback in Concert: Effects on AGN Accretion, Obscuration, and Outflows

Philip F. Hopkins, Paul Torrey, Claude-André Faucher-Giguère, Eliot Quataert, & Norman Murray
TAPIR & The Walter Burke Institute for Theoretical Physics, Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA
Department of Physics, MIT, 77 Massachusetts Avenue, Cambridge, MA 02139, USA
Department of Physics and Astronomy and CIERA, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA
Department of Astronomy and Theoretical Astrophysics Center, University of California Berkeley, Berkeley, CA 94720
Canadian Institute for Theoretical Astrophysics, 60 St. George Street, University of Toronto, ON M5S 3H8, Canada
Canada Research Chair in Astrophysics

We study the interaction of feedback from active galactic nuclei (AGN) and a multi-phase interstellar medium (ISM), in simulations including explicit stellar feedback, multi-phase cooling, accretion-disk winds, and Compton heating. We examine radii pc around a black hole (BH), where the accretion rate onto the BH is determined and where AGN-powered winds and radiation couple to the ISM. We conclude: (1) The BH accretion rate is determined by exchange of angular momentum between gas and stars in gravitational instabilities. This produces accretion rates , sufficient to power luminous AGN. (2) The gas disk in the galactic nucleus undergoes an initial burst of star formation followed by several Myrs where stellar feedback suppresses the star formation rate (SFR). (3) AGN winds injected at small radii with momentum fluxes couple efficiently to the ISM and have dramatic effects on ISM properties within pc. AGN winds suppress the nuclear SFR by factors and BH accretion rate by factors . They increase the outflow rate from the nucleus by factors , consistent with observational evidence for galaxy-scale AGN-driven outflows. (4) With AGN feedback, the predicted column density distribution to the BH is consistent with observations. Absent AGN feedback, the BH is isotropically obscured and there are not enough optically-thin sightlines to explain Type-I AGN. A ‘torus-like’ geometry arises self-consistently as AGN feedback evacuates gas in polar regions.

galaxies: formation — galaxies: evolution — galaxies: active — star formation: general — cosmology: theory

1 Introduction

The masses of super-massive black holes (BHs) correlate with various host galaxy bulge properties (Magorrian et al. 1998; Ferrarese & Merritt 2000; Gebhardt et al. 2000; Hopkins et al. 2007; Aller & Richstone 2007; Feoli & Mancini 2009; Kormendy et al. 2011; for a review see Kormendy & Ho 2013). The small scatter in these correlations (relative to other galaxy properties; Hopkins et al., 2009b; Kormendy & Ho, 2013), together with constraints indicating that most BH mass is assembled in an optically bright quasar phase (Soltan, 1982; Salucci et al., 1999; Yu & Tremaine, 2002; Hopkins et al., 2006b), has led to the development of models where large-scale effects of feedback from accretion self-regulate BH growth at a critical mass (Silk & Rees, 1998; King, 2003; Di Matteo et al., 2005; Murray et al., 2005). Gas inflows triggered by some process fuel rapid BH growth in a nuclear starburst (Diamond-Stanic & Rieke, 2012; Mushotzky et al., 2014), until feedback begins to expel nearby gas and dust. This “blowout” results in a short-lived, bright optical quasar that, having expelled its fuel supply, fades and leaves a remnant on the observed BH-host correlations (Hopkins et al., 2005a, c). This general scenario has been able to explain many quasar observables, including luminosity functions, lifetimes, and BH mass functions (Hopkins et al., 2005b, 2006c, 2008b, 2009a; Volonteri et al., 2006; Menci et al., 2003; Somerville et al., 2008; Lapi et al., 2006; Tortora et al., 2009). It has also been speculated that this feedback might ultimately have a large impact throughout the AGN host galaxy, expelling or heating gas and explaining the rapid quenching of star formation in massive galaxies (Granato et al., 2004; Scannapieco & Oh, 2004; Croton et al., 2006; Hopkins et al., 2008a; Antonuccio-Delogu & Silk, 2008), and considerable observational evidence has emerged for this in recent years (see e.g. Nesvadba et al., 2010; Cimatti et al., 2013; LaMassa et al., 2013; Shimizu et al., 2015; Guillard et al., 2015; Alatalo et al., 2015).

High-velocity outflows can be driven from the BH accretion disk by a variety of physical processes including, e.g., radiation pressure on lines and dust, magnetic processes, or Compton heating (see e.g. Blandford & Payne, 1982; Begelman, 1985; Chang et al., 1987; Sanders et al., 1988; Konigl & Kartje, 1994; Murray et al., 1995; Elvis, 2000; Proga, 2000, 2007; Silk, 2005; Murray et al., 2005; Batcheldor et al., 2007; Tortora et al., 2009). These manifest themselves observationally as ultra-fast outflows (Tombesi et al., 2010, 2013, 2015, e.g.), the broad emission line regions and broad absorption line quasars (e.g. Weymann et al., 1981; de Kool et al., 2001; Gabel et al., 2006; Ganguly et al., 2007), more moderate velocity outflows () associated with the narrow line region (Laor et al., 1997; Crenshaw et al., 2000; Steenbrugge et al., 2005; Krongold et al., 2007), as well as quasar absorption and occultation systems (e.g. McKernan & Yaqoob, 1998; Turner et al., 2008; Miller et al., 2008). Observations on galaxy scales have also provided strong evidence for powerful molecular, atomic, and ionized outflows with velocities , outflow rates up to times the BH accretion rate, and spatial extents of  kpc (Tremonti et al., 2007; Prochaska & Hennawi, 2009; Moe et al., 2009; Wild et al., 2009; Fischer et al., 2010; Humphrey et al., 2010; Dunn et al., 2010; Bautista et al., 2010; Feruglio et al., 2010; Sturm et al., 2011; Rupke & Veilleux, 2011; Coil et al., 2011; Greene et al., 2011, 2012; Faucher-Giguère et al., 2012; Borguet et al., 2012; Cimatti et al., 2013; Cicone et al., 2014; Harrison et al., 2014, 2015; Zakamska & Greene, 2014; Zakamska et al., 2015). In some cases, however, it remains unclear to what extent these outflows are driven by AGN activity vs. star formation.

The physics of how AGN-powered outflows interact with the ISM and affect the fueling of the AGN itself – how inflow and outflow are governed on scales between the small-scale viscous accretion disk (pc) and the galaxy proper ( kpc) – remains highly uncertain. There have been many theoretical studies of different “modes” of AGN feedback (see references above, as well as Ostriker et al. 2010; Ciotti et al. 2010; Choi et al. 2014, 2015; Steinborn et al. 2015, and references therein). However, to date, most of these studies have treated the interstellar medium with relatively simple “sub-grid” prescriptions that ignore the additional complications introduced by stellar feedback, interstellar turbulence, and/or the small-scale phase structure of the medium around an AGN. In order to build on these models and model the interaction of AGN outflows and the ISM with greater fidelity, it is critical to include both a realistic description of the physics of the ISM, star formation, and stellar feedback, as well as a plausible description of AGN feedback mechanisms.

Towards this end, in this paper we use a suite of numerical simulations to study the interaction of quasar-driven winds and a multi-phase ISM. In a series of papers (Hopkins et al., 2011, 2012d) (hereafter Papers I & II, respectively), we have developed a new set of numerical methods to explicitly model some of the key processes that shape the multi-phase ISM; the simulations include physically motivated, but still subgrid, treatments of stellar radiation pressure, HII photoionization and photoelectric heating, and the heating, momentum, and mass deposition by supernovae (SNe) and stellar winds. The feedback is tied to the young stars with energetics and time-dependence taken directly from stellar evolution models – this is particularly important in galactic nuclei, since the dynamical times become shorter than stellar evolution timescales. In a series of papers (Hopkins et al., 2014; Muratov et al., 2015; Oñorbe et al., 2015; van de Voort et al., 2015; Ma et al., 2015, 2016; Faucher-Giguere et al., 2015), we showed that, on galactic scales, these models produce a quasi-steady ISM in which molecular clouds form and disperse rapidly, with phase structure, turbulent velocity dispersions, and disk and GMC properties in reasonable agreement with observations. Here, we combine these models with models for AGN accretion and feedback via both Compton heating and high-velocity winds from the AGN accretion disk, and examine how various forms of AGN feedback affect black hole accretion, AGN obscuration, and the generation of galaxy-scale outflows. We focus on scales of pc, where the accretion rate onto the black hole is determined, and where AGN-powered winds and Compton heating couple to the ISM. For comparison, the “radius of influence” of the BH, inside of which it dominates the potential, is pc in the case we study. This is the first in a new series of papers so we highlight a few of the key results but leave more detailed studies for future work.

The remainder of this paper is organized as follows. §2 summarizes our galaxy models and our treatment of radiative cooling, star formation and BH growth, and stellar and AGN feedback. §3 summarizes the results of simulations with stellar feedback only, while §4 compares these results to simulations that include AGN feedback. §5 summarizes and discusses our key results. A series of Appendices contain key technical results. Appendix A describes our implementation of BH feedback. Appendix B summarizes the effects of including short timescale variability in the assumed BH accretion rate. Appendix C describes convergence tests and the effects of using alternate numerical methods. Appendix D shows that in-shock cooling does not compromise our results.

Figure 1: Morphology of the gas in a standard simulation, in face-on (; left), side-on (; middle), and cylindrical (; right) projections. The time (Myr since the beginning of the simulation) is () orbital periods at pc ( pc). Brightness encodes projected gas density (increasing with density; logarithmically scaled with a dex stretch); color encodes gas temperature with blue material being K molecular gas, pink  K warm ionized gas, and yellow K hot gas. Top: Simulation with stellar, but no AGN feedback (no_BAL in Table 1). A multiphase disk forms; it is mostly molecular inside the central pc, with heating by HII regions very localized to small ionized “bubbles” and heating by SNe restricted to low-density regions where it can vent vertically. The central  pc develops a stellar+gas accretion disk dominated by modes. Bottom: Same, with broad-absorption line winds (v5000). The winds blow out a polar cavity and generate an expanding shell in-plane, with occasional dense clumps sinking through to the center. Feedback eventually evacuates the entire nuclear region.

[ caption=Simulations,center]lccccl \tnote[ ]Parameters describing the simulations in the text: Each employs a gas particle mass of and minimum SPH smoothing length of pc. Additional simulations for numerical tests are in Appendix C.
(1) Model name
(2) : Momentum-loading of BAL wind feedback ()
(3) : Energy-loading of BAL wind feedback ()
(4) : Mass-loading (determined by )
(5) : AGN wind launching velocity at the simulation resolution (in ; determined by ) Model & & & & & Notes
no_BAL & 0 & 0 & 0 & 0 & no AGN FB
v5000 & 1 & 0.008 & 6.0 & 5,000 & “default”
v5000_hiP & 10 & 0.08 & 60 & 5,000 & high-momentum
v5000_loP & 0.1 & 0.0008 & 0.6 & 5,000 & low-momentum
v30000 & 1 & 0.05 & 1.0 & 30,000 & high-energy
v500 & 1 & 0.0008 & 60 & 500 & low-energy
v5000_C & 1 & 0.008 & 6.0 & 5,000 & +Compton heating
v5000_iso & 1 & 0.008 & 6.0 & 5,000 & isotropic winds

2 The Simulations

The simulations were performed using the GIZMO code (Hopkins, 2015). GIZMO is a multi-method code which can be run with any of several hydro solvers; here we run the code in its smoothed-particle hydrodynamics (SPH) mode, specifically in the “pressure-entropy” (“P-SPH”) form which includes several improvements relative to older SPH implementations. Specifically, this is a heavily modified version of the parallel TreeSPH code GADGET-3 (Springel, 2005), in a fully conservative formulation (Springel & Hernquist, 2002) which is also density-independent in a manner that allows contact discontinuities and improved fluid mixing (Hopkins, 2013, see Appendix C). The artificial viscosity, adaptive timestepping, and smoothing kernel are updated following Hopkins (2013). The galaxy models and the treatment of star formation and stellar feedback are described in detail in Paper I (Sec. 2 & Tables 1-3) and Paper II (Sec. 2). We briefly summarize the salient properties here.

Figure 2: Face-on morphology of the gas in the additional simulations from Table 1, as Fig. 1, at the same time (and same scale). Top: Compton heating (v5000_C) and isotropic vs. disk-planar winds (v5000_iso) have no visible effects (the differences are consistent with stochastic variations run-to-run). Middle: Lowering (v5000_loP) or raising (v5000_hiP) the momentum-loading of the winds leads to smaller/larger bubbles after the initial accretion event, although these in turn alter the subsequent accretion rate (see § 4). Bottom: Lowering the wind velocity (at fixed momentum-loading; v500) has no significant effect (although the shocked gas is colder, as expected). Raising the wind velocity to (v30000) creates the most hot gas (as expected); however, this gas appears to mostly vent out from the central regions, so a nuclear disk similar to the no-feedback case re-forms within pc which drives a higher accretion rate at later times.

2.1 Initial Conditions

The initial conditions are a gas-rich nuclear disk in a massive galaxy, drawn from the large parameter survey of Hopkins & Quataert (2010a).We consider a BH (initial ) in a Hernquist (1990) stellar bulge (, isotropic orbits and scale-length kpc) and halo (, with virial radius, concentration, and velocity appropriate at ). The BH is surrounded by an exponential nuclear disk of gas and stars (scale-lengths pc and pc, and , respectively; stellar disk with vertical sech profile and dispersions such that , gas disk initially thermally supported with ). The initial surface densities of the gas and stellar disk are thus pc g cm.

This is chosen, based on the survey in Hopkins & Quataert (2010a), to provide a “best case” for extremely rapid BH growth and quasar-level fueling. It is motivated by large-scale simulations of major mergers which produce dense, torus-like structures and high accretion rates (Hopkins & Quataert, 2010b, 2011b; Hopkins & Quataert, 2011a; Hopkins et al., 2012a, b; Hopkins, 2012), as well as at least some observations indicating the presence of powerful AGN in nuclear starburst “cusps” even in galaxies which may not be experiencing extended star formation (Diamond-Stanic & Rieke, 2012; Mushotzky et al., 2014; Alatalo et al., 2015).

The initial gas disk contains particles; the initial gas particle mass is . We consider a limited resolution comparison in Appendix C. The force softening for the BH, gas, and star particles is set to  pc, with minimum SPH smoothing length times this. We note that all simulations employ the more sophisticated formulation of artificial viscosity described in Morris & Monaghan (1997), which greatly reduces numerical dissipation away from shocks relative to earlier implementations (see e.g. Rosswog et al., 2000; Price, 2008).

2.2 Cooling, Star Formation, & Stellar Feedback

Gas follows an atomic cooling curve with additional fine-structure cooling to  K. Metal-line cooling is followed species-by-species for 11 tracked species as in Wiersma et al. (2009a, b). The enrichment for each species is followed with the time dependent metal flux directly attached to the mass, momentum and energy flux from stellar winds and SNe Types Ia & II (see Hopkins et al., 2012c, 2013b).

Star formation is allowed only in dense, molecular, self-gravitating regions above . We follow Krumholz & Gnedin (2011) to calculate the molecular fraction in dense gas as a function of local column density and metallicity, and allow SF only from molecular gas. Following Hopkins et al. (2013c), we also restrict star formation to only gas which is locally self-gravitating, i.e. has (where the limit is taken as the “averaging radius” vanishes, allowing to be calculated in a resolution-independent manner only as a function of local properties). Gas which meets all of these criteria forms stars at a rate (i.e.  efficiency per free-fall time). As shown in Hopkins et al. (2013c), the molecular criterion is not especially important in galaxy centers since most of the dense gas is molecular already, but the self-gravity criterion is important for small scales around black holes, where any simple constant-density threshold for star formation fails to account for the radially-dependent tidal forces. Even in these regions however, the role of these criteria is primarily to determine where stars form; Hopkins et al. (2013c) and a number of other studies have shown that the total star formation rate, once fragmentation and stellar feedback are resolved, is set by stellar feedback, and is largely insensitive to details of both cooling and star formation prescriptions (see Saitoh et al., 2008; Hopkins et al., 2011, 2012d, 2012c, 2013b, 2013a, 2015; Agertz et al., 2013).

Once stars form, feedback is included in the form of radiation pressure (UV, optical, and IR, allowing for multiple-scattering), stellar winds (fast, young star winds and slow AGB winds), SNe (types Ia and II), photo-ionization and photo-electric heating. Every star particle is treated as a single stellar population with an age based on its formation time and metallicity and mass inherited from its parent gas particle. Feedback includes the relevant mass, metal (with separately tracked species), momentum, and energy injection to the neighboring gas; all of the relevant quantities (stellar luminosities, spectral shapes, SNe rates, wind mechanical luminosities, yields) for the mechanisms above are tabulated as a function of time directly from the stellar population models in STARBURST99, assuming a Kroupa (2002) IMF. For every SNe event (or every timestep for winds and single-scattering photon momentum), the relevant energy, momentum, mass, and metals are deposited into the nearest gas particles surrounding each star particle; long-range photo-heating and radiation pressure are treated in a simplified manner assuming spherically symmetric photon propagation from each star particle as an independent source. See Hopkins et al. (2011, 2012d) for details. The end result of this stellar feedback is a multiphase ISM with a broad range of densities and temperatures.

2.3 Black Hole Growth & Feedback

The simulations all include super-massive BHs. The BH is much more massive than the stellar/gas particles, so we do not need to artificially “force” the BH particle to stay in the center of the potential, but let it move freely. We cannot, however, directly resolve the viscous accretion disk of the BH on scales pc. We therefore simply assume that the BH immediately accretes any gas particle gravitationally bound to it (relative velocity less than escape) with apocentric radius (calculated from the particle position and velocity relative to the BH) (the minimum Keplerian distance). The rate of particle accretion is capped at the Eddington limit.

The BH radiates at a luminosity ( is assumed).111We also describe in Appendix B a model which imposes a spectrum of sub-grid time variability in the accretion rates; however this has no significant effects on the time-averaged results here. The explicit details of the BH feedback implementation are given in Appendix A; we briefly summarize them here. Since quasars are believed to have high-velocity, near-planar winds driven off the accretion disk (e.g., Murray et al. 1995), we assume that a fraction of the photon momentum drives a wind launched at the resolution scale around the BH from accreted gas. Specifically a fraction of any gas accreted is blown out as a wind with velocity , planar with the inflow (by launching particles directly at the accretion radius with this velocity). Two parameters define the wind, the mass-loading and velocity; this is equivalent to specifying the momentum-loading () and energy-loading () of the wind. Values for the simulation parameters are in Table 1.

We also include Compton heating & cooling from the radiation field. Following Sazonov et al. (2004), this can be approximated with a nearly obscuration-independent Compton temperature of K. We add the appropriate Compton rates to the standard cooling function (with a limiter following Faucher-Giguere & Quataert 2012 to account for rate-limiting by Coulomb collisions at the high temperatures that can obtain in strong shocks).

Figure 3: Top Panel: Star formation rate within 10 and 500 pc regions for simulations with (v5000) and without (no_BAL) AGN feedback. Middle and Bottom Panels: Location of the same simulations on the Kennicutt-Schmidt relations at different times. The star formation rate surface density and gas surface densities are averages within 10 pc and the rotation rate in the bottom panel is also measured at 10 pc. The observations in the middle panel (dashed line dex is from Narayanan et al. (2012)’s variable model) are based on a range of galaxies, not just galactic nuclei, but nonetheless provide a useful point of comparison. The star formation efficiency per dynamical time evolves significantly with time during the simulation, with a relatively high star formation efficiency in the burst of star formation at early times followed by a more prolonged period of lower star formation efficiency.
Figure 4: Time-averaged structural properties of the simulations. Top: mode amplitude in the cold molecular gas, as a function of radius. With no BAL feedback the large spiral modes are visible here; with BAL winds the order-unity asymmetries introduced by the AGN wind impacting the ISM dominate. Bottom: Gaussian disk scale-height () versus radius. With no BAL winds, a modest is supported by the combination of stellar feedback and gravitational instabilities. With BAL winds, is greatly enhanced because there is little gas and it is often dominated by escaping/venting polar winds. Even in the latter case, the scale-height of the cold rotating gas remains modest, similar to that in the non BAL wind simulation (see Fig. 1).
Figure 5: Top: BH accretion rate vs. time. To make systematic differences clear we smooth the rates in a yr window (see Appendix B for an un-smoothed case). AGN feedback suppresses the BH accretion rate relative to simulations without AGN feedback, with higher momentum-loading in the input AGN wind (higher ) leading to lower . Middle: Total momentum flux in the galaxy-scale outflow at pc in each model, vs. time. The models with BAL winds all equilibrate at broadly similar outflow momentum flux. This explains why higher- models adjust to have lower . Bottom: Corresponding mass-outflow rate of the wind at pc. Simulations with AGN feedback have dramatically larger outflow rates from the nuclear region relative to the simulation without AGN feedback.
Figure 6: Inflow/outflow rate vs. radius in the model with stellar feedback alone, averaged over time for a few dynamical times. Negative values (outflow) are dotted, positive values (inflow) are solid (absolute value plotted for the sake of a logarithmic projection). We compare several analytic accretion rate models (§ 3). The “gravitational torques” estimator (resonant exchange between gas+stellar gravitational instabilities) is accurate within a factor at all radii pc, and correctly predicts the major sign (inflow/outflow) changes. Spherical accretion models fare poorly: the “pure Bondi” estimator gives , too large to fit on the plot; the “modified Bondi-Hoyle” estimator over-predicts by dex. “Ballistic accretion” from turbulence fares poorly in the opposite manner (predicting at pc). “Gravito-turbulent viscosity” is dimensionally reasonable, but under-estimates by factors of near the BH radius of influence, where gravitational torques are most prominent, and does not capture the sign information. Around  pc, resolution effects from our discrete particle number become important (some predictions drop precisely because the BH is accreting particles).

3 Results with Stellar Feedback, but No Black Hole Feedback

Fig. 1 (top row) shows the morphology of the high-resolution no AGN feedback run at a typical time after a few orbital periods, when the system has reached an approximate statistical steady state (Fig. 2 compares the additional simulations in Table 1). Figure 3 shows the star formation rate as a function of time for simulations with and without AGN feedback (top panel) and two versions of the Kennicutt-Schmidt relation describing the star formation law for these nuclear-scale simulations (bottom two panels). Note that in the simulation with only stellar feedback, there is an initial burst of star formation but after a few Myr, the star formation rate settles into an approximate steady state at yr within kpc. The image in Fig. 1 is shown in the latter phase. Within pc, stellar feedback alone does clear most of the gas after a few Myr; this is recycled in a small-scale fountain on a similar timescale. The dynamics of these small-scale burst-quench cycles is explored in more detail in Torrey et al. (2016).

3.1 Black Hole Accretion

Fig. 1 shows that the gas disk exhibits strong non-linear spiral wave and eccentric/lopsided disk modes, which are visible in spite of the inhomogeneous structure of the ISM. Using simulations on similar spatial scales but with a much less realistic model of the ISM, Hopkins & Quataert (2010a) showed that non-linear modes generated by stellar-gas interactions dominate the angular momentum transport in galactic nuclei at and inside the BH sphere of influence. However, that study was limited by the assumption that the ISM gas was described by a simple “effective equation of state” (i.e. no resolved phase structure, winds, or turbulence, simply single-phase gas with a barytropic non-thermal pressure following Springel & Hernquist 2003). We confirm their result here with a much more realistic ISM model.

Fig. 4 (top panel) plots the mode amplitudes222Mode amplitudes are measured in the gas surface density as

versus radius for the simulations with and without AGN feedback. For the simulation without AGN feedback, the mode amplitude found here is similar to that found in Hopkins & Quataert (2010a)’s simulations with gas fractions . This suggests that the mode excitation and saturation physics is at least broadly similar in spite of the more dynamic multi-phase ISM present in our simulations.

Fig. 5 shows the black hole accretion rate,333To highlight the systematic differences between runs, we boxcar-smooth each curve in a time window of yr. In Appendix B we show that there is variability on all resolved timescales in the simulations (as also seen by Novak et al., 2011; Gan et al., 2014; Dubois et al., 2014), and we even consider a model for un-resolved time variability. However, we caution that some of the resolved small-timescale variability is almost certainly artificial here (owing to the assumption that individual gas particles are accreted discretely) – there are such accretion events per simulation. The BH accretion rate would be likely smoothed out if these were accreted into a viscous disk, which then accretes onto the BH. the outflow rate from the galactic nucleus, and the total momentum flux in the outflow as a function of time.444The net inflow rate at radius is given by in an annulus. The outflow rate is the same integral, but only over where . The rates are time-averaged in each annulus (which also removes the spurious radial velocity contribution from e.g. stationary modes). Because of finite bin-widths the inflow rate can change sign discretely from bin-to-bin. The simulations clearly find large inflows up to to the central pc. Hopkins & Quataert (2011a) derive an analytic approximation for the inflow rate through each annulus for inflows driven by strong gravitational torques and resonant angular momentum exchange between gas and stars. For modes with complex potential and pattern speed , this is:


with a phase function and an order-unity amplitude correction derived in Hopkins & Quataert (2011a), which can be measured directly in the simulations. For an mode in a quasi-Keplerian potential, this is approximately (with the mode amplitude).

Fig. 6 compares equation 2 to the simulation inflow/outflow rates; the agreement is reasonable, particularly given that the analytic result was derived under the assumption of smooth (non-turbulent) gas flows. Note that divergence between models at pc is primarily a consequence of our discrete accretion model actually removing particles at this radius, and should be considered with caution.

Fig. 6 also compares the inflow rate in our simulations to four alternative proposed accretion rate estimators, none of which does as good a job of reproducing the simulation results. (1) Bondi: . This over-predicts the accretion rate by an enormous factor as most of the gas is cold and molecular, supported not by pressure but by angular momentum. (2) Modified Bondi-Hoyle: . This allows for the fact that mass outside the BH itself (e.g. the bulge or nuclear cluster) should, when viewed from gas at large enough distance, act as a point mass in the same way; it also allows for super-sonic relative motion of the gas and BH. This is dimensionally closer to what we measure than [1], but given the low , it amounts to assuming all gas is in free-fall (neglects angular momentum) – i.e. it is quite similar to simply taking – and it over-predicts by factors . Note that replacing , as in the standard Bondi-Hoyle formulation (used in Springel et al., 2005; Hopkins et al., 2006a, 2005a; Di Matteo et al., 2008; Croft et al., 2009) only slightly decreases the discrepancy. (3) Ballistic Accretion: (this corresponds to accretion of the randomly-populated low-angular momentum “tail” of highly turbulent flows from Hobbs et al., 2011, we generalize their formulae for accretion through each annulus). This disagrees with the simulations as well; dimensionally it gives but with a “reduction factor” , which for found here is very small, so that there is very little ballistic accretion. (4) Gravito-turbulent viscosity: where is the (cooling function-dependent) effective turbulent viscosity for a disk (Gammie, 2001; Thompson et al., 2005; Debuhr et al., 2010, 2011).555In Fig. 6, we adopt , close to the predictions from Gammie (2001) given the measured Mach number in the diffuse gas. Changing the value of will systematically shift the normalization of the predicted curve. This is dimensionally similar to the gravitational torques scaling but with free-fall slowed by a term instead of ; over some radii the two are comparable but the former decreases rapidly inside the BH radius of influence (implying accretion would be “throttled”) while can remain order-unity all the way to the true accretion disk (see Tremaine, 1995; Bacon et al., 2001; Hopkins & Quataert, 2011a; Hopkins & Quataert, 2010b; Hopkins, 2010).

The comparisons in this section are based on simulations without AGN feedback. In the presence of feedback, the net accretion rate onto the BH is determined by a competition between the inflow rate from large scales set by gravitational torques and the efficiency of AGN feedback at suppressing this inflow in the galactic nucleus. Our simulations explicitly resolve this competition and produce accretion rates a factor of lower than in simulations without AGN feedback (Fig. 5). For lower resolution galaxy-scale or cosmological simulations it is unclear what the best time averaged accretion rate estimator is to capture this competition between inflow by gravitational torques and AGN feedback; this merits further study in future work.

3.2 Star Formation and Vertical Disk Structure

Fig. 4 (top panel) shows the vertical scale height of the gas disk as a function of radius. The disk is in vertical equilibrium but the dispersions are turbulent (much larger than thermal). As shown in Paper II on larger scales, stars form roughly until feedback can maintain Toomre () and offset further collapse. At large radii this gives ; at pc this is (dispersions ). The cylindrical image in Figure 1 highlights the modest thickening of the disk at larger radii that is qualitatively analogous to that required in AGN ”torus” obscuration models. As we describe in §4 this effect is much more dramatic in simulations with AGN feedback because feedback efficiently evacuates the polar region of gas.

The bottom panels of Figure 3 shows our simulations in two common versions of the Kennicutt-Schmidt relation. The star formation rate surface density and gas surface densities are averages within 10 pc and the rotation rate in the bottom panel is also measured at 10 pc. The observations in Figure 3 are best fits from Narayanan et al. (2012) based on a variable factor. They are shown to provide a point of comparison, but include a range of galaxies, not just galactic nuclei. The time averaged star formation efficiency in Figure 3 is broadly consistent with observations. The efficiencies evolve significantly with time, however, with a relatively high star formation efficiency in the burst of star formation at early times followed by a more prolonged period of lower star formation efficiency. Perhaps most striking is that the star formation efficiency per dynamical time decreases by nearly a factor of during the course of the simulations without AGN feedback. Thus the decline in the star formation rate is not simply due to gas depletion but is also due to the decreasing star formation efficiency. Note that the duration of the simulation is comparable to the lifetime of massive stars. Thus the stellar feedback that is effective for most of the duration of the simulation is that due to stellar radiation and stellar winds, since supernovae only start after Myr and have not had significant time to operate. In addition, because the local dynamical time is short compared to the lifetimes of massive stars, the efficiency of stellar feedback depends primarily on the surface density of young stars, rather than the star formation rate. We explore the consequences of this for the “burstiness” of nuclear star formation and origins of the nuclear-scale Kennicutt-Schmidt relation in a companion paper (Torrey et al., 2016).

4 Results with Black Hole Feedback

Figure 7: Column density distribution on sightlines to the BH in each simulation. We integrate the column along sightlines to each BH at each time following Hopkins et al. (2005d), uniformly sampling the sky in solid angle, and show the distribution over all sightlines and times in each simulation. Stellar feedback alone produces a relatively narrow range of very large columns. Simulations with BAL winds have evacuated polar regions with column densities cm and overall broader obscuring column density distributions. The “clumpy torus” thus naturally arises from AGN feedback interacting with the large-scale ISM at pc.
Figure 8: Median column density as a function of polar angle (averaged over azimuthal angle ), calculated over all times as Fig. 7, for the simulations in Table 1. The disk mid-plane () features Compton-thick columns, as expected, which decline towards the poles (). By clearing the central regions, feedback suppresses the column densities at all , but proportionally more towards the poles (where the smaller initial column makes complete evacuation easier) – this is needed to create optically thin sightlines. In all the simulations, the scatter in column density at a fixed is dex, owing to sub-structure in the gas (Fig. 1). Thus the inclusion range for v5000 (shaded) is larger than the systematic offset between it and any other simulation with feedback (but not the no-feedback case).

We now consider the results of simulations with AGN feedback, focusing on the fiducial v5000 run in which the AGN wind at small radii is injected with and . However, we consider variations in both the mass and momentum-loading as well, as outlined in Table 1.

Fig. 1 (bottom) shows the gas morphology at a few Myr in our v5000 simulation; there is a clear dramatic impact of feedback on the gas, with the central 30 pc relatively evacuated of gas by the 3 Myr time of these images. Qualitatively, our other simulations with feedback (Fig. 2) resemble this case, but with subtle differences discussed below.

The v5000 outflows are launched in the dense disk mid-plane. This drives an expanding shell in the disk plane, with gas piled up in a narrow ring/shell at the outer (radiative) shock where the winds are encountering the ISM. This is similar to Faucher-Giguere & Quataert (2012)’s models for galaxy-scale winds driven by AGN, though it is not clear if those models quantitatively apply because the hot shocked gas created by the AGN wind is not well-confined – this may be a limitation of the small scales we are simulating (there is no full galaxy and halo into which the winds can propagate in these simulations). Indeed, out of the midplane, the entrained mass is modest so outflows coast or are accelerated by hot gas pressure filling the growing central cavity in the disk.

The large impact of the AGN wind on the ambient gas has three closely related effects. First, it strongly suppresses the star formation in the galactic nucleus, by a factor of (Fig. 3). Secondly, it increases the net outflow rate from the galactic nucleus by a factor of , to yr (Fig. 5). Finally, on longer timescales the BH feedback roughly regulates the BH accretion rate. Specifically, the feedback momentum flux scales as ; balancing infall with feedback therefore implies a critical value of , so in equilibrium . This scaling provides a reasonable approximation over a sufficient time average (Fig. 5), but the evacuation of the central regions clearly leads to very large-amplitude variability on  yr timescales.

It is useful to directly compare the momentum flux in the galaxy scale winds (Fig. 5) with those injected at small radii in the AGN wind. These need not be the same if AGN feedback produces a bubble of hot gas that does work on the surrounding material, increasing the momentum flux in the wind (e.g., Faucher-Giguere & Quataert 2012). For low input wind velocities the momentum fluxes in Fig. 5 are comparable to that injected in the AGN wind at small radii, while for higher input wind velocities (in particular, the v30000 simulation), there is a factor of few boost in the AGN wind momentum flux – although we inject a momentum flux , the outflow momentum fluxes reach , comparable to observed winds on larger-scales which have been decelerated to (Borguet et al., 2012; Cimatti et al., 2013; Cicone et al., 2014; Harrison et al., 2014; Zakamska & Greene, 2014; Zakamska et al., 2015). The modest boosts found here are because gas shocked heated by the AGN wind is able to escape relatively easily along the polar direction. In a more self-consistent calculation, it is possible that the existence of large warps between the disk axis on small scales and that on large scales (e.g., Hopkins et al. 2012b) might act to better confine the outflow.

Note that the outflow rates of that we find are relatively modest, compared to many observations at  kpc scales (see references above). This, of course, owes in part to the limited region we are simulating (pc) – there simply is not much material in this volume for the winds to “sweep up.” On these scales directly, there are actually relatively few constraints on AGN outflow rates and velocities (most of the constraints above apply either to AGN accretion disk scales, or spatially-resolved outflows on kpc scales). However, if we assume the outflow remains momentum-conserving, then for typical galaxy mass profiles we would easily expect it to entrain an order or magnitude or more mass as it propagates from  kpc to a few kpc.

Figs. 7-8 plots the column density distribution and dependence on viewing angle (both averaged over time and in various time intervals) for each of our simulations, both with and without AGN feedback. The model without AGN feedback predicts virtually no systems with column densities below cm even along polar sightlines, in stark contrast to observations. Even though the ISM is highly inhomogeneous on larger scales, a small, dense thick-disk or “halo” component surrounding the BH in the central  pc is sufficient to produce these extremely high column densities even in the polar direction. The BAL winds have, however, an enormous impact on the column density distribution. This is not surprising given their impact on the nuclear gas morphology. The polar regions are completely evacuated, giving a large fraction of sightlines that are fully un-obscured. The remaining sightlines follow a broad column density distribution, driven in part by the fragmentation and asymmetries seen in the expanding equatorial shells. The evacuation of the central regions out to some radius where gives a canonical “torus-like” global morphology. This is particularly clear in the cylindrical image shown in the lower right panel of Figure 1.

4.1 Dependence on the Strength and Form of AGN Feedback

In Table 1, we outline a series of runs changing the energy and momentum loading of the AGN-driven winds. Here we compare their properties.

First, we add Compton heating/cooling to our “standard” case (run v5000_C). Consistent with previous studies (Ostriker et al., 2010; Ciotti et al., 2010; Choi et al., 2012, 2014, 2015; Park et al., 2014) and observational constraints (Stern et al., 2015; Chatterjee et al., 2015), this has little or no effect on the properties we measure (see Figs. 2-7). This should not be surprising: all but the lowest-density gas in the simulations has cooling times much shorter than the Compton heating time. To the extent that Compton heating/cooling is important, it has been speculated that it may be important for cooling in the reverse shock of the BAL winds; however, in the cases we simulate this is not expected to dominate even for a spherical blastwave (see Faucher-Giguere & Quataert, 2012), and “venting” rather than cooling dominates the escape of energy in the hot gas.

We also consider a case (v5000_iso) where the winds from the BH are directed isotropically from the BH, as opposed to in the accretion plane. This is discussed in § A; Figs. 2-7 demonstrate this produces only very small changes (within the range produced by purely stochastic effects) relative to our standard v5000. This is also consistent with previous studies (Debuhr et al., 2012), and follows simply from the fact that the distribution of accretion directions/angles (the fractional scale-height of accreting gas) is relatively large.

Our v5000_hiP and v5000_loP runs keep the wind velocity fixed but increase/decrease the mass (and momentum) loading by an order of magnitude, respectively, relative to the BH accretion rate. Not surprisingly, a lower (higher) mass-loading produces initially weaker (stronger) outflows: the “bubble” in Fig. 2 is smaller (larger), and this produces slightly higher (lower) column density sightlines in Fig. 7. This has relatively weak effects on the gas properties outside the region being evacuated (Fig. 4). The SFR and Kennicutt-Schmidt relation are essentially indistinguishable from those in Fig. 3 for our v5000 run, except higher/lower mass-loading translates to earlier/later suppression of the SFR in the immediate vicinity of the BH. The initially stronger outflow (for the same accretion rate) in v5000_hiP suppresses the BH accretion rate substantially in Fig. 5, while the weaker outflow in v5000_loP allows more rapid growth of the accretion rate compared to v5000. Interestingly, this ends up producing a very similar total momentum flux and outflow rate in the winds at Myr. This suggests that indeed we are seeing self-regulation when the BH injects sufficient momentum into the medium to drive outflows that clear its vicinity. This critical asymptotic value is reached eventually in the different runs, and appears to only weakly depend on the feedback parameters. However, how quickly the critical point is reached, and (correspondingly) how much the BH is able to grow, is strongly dependent on the initial feedback mass-loading.

In our v500 and v30000 runs, we keep the momentum-loading fixed but vary the initial wind velocity from . Decreasing the velocity to leads to slightly higher average accretion rates in Fig. 5, but the effect on the circum-nuclear structure (Figs. 2-4) and column density distribution (Fig. 7), and outflow rates (Fig. 5) are weak. This is not surprising since in this limit, the winds are primarily momentum-conserving and so the absolute wind speed (at fixed momentum) does not qualitatively change the dynamics. With winds, however, we see interesting differences. The faster velocity produces more shock-heated high-temperature gas (unsurprisingly); the lower initial wind mass-loading, however, also appears to lead to more efficient “venting” of the wind. This means there is less of a circum-nuclear “bubble” carved out in the cold/neutral gas in Fig. 2, but rather more pronounced hot gas channels escaping. This in turn allows a thin nuclear disk to re-form (visible morphologically in Fig. 2, but also manifest in higher typical column densities in Fig. 7), which then produces higher BH accretion rates seen in Fig. 5. These higher inflow rates rates lead to a larger net wind momentum injection and outflow rate. Some of these differences are similar to the behavior seen in the “thermal energy deposition” models of Choi et al. (2014, 2015); however, those simulations did not include the physics driving the multiphase structure of the ISM, and so could not capture the full magnitude of “venting” effects we see here.

5 Discussion

We have used simulations with pc resolution to study BH accretion and feedback in gas-rich nuclear disks around massive BHs accreting at quasar-like luminosities. Our calculations include an explicit treatment of star formation and stellar feedback, which produce a self-consistently inhomogeneous ISM. We model AGN feedback via Compton heating/cooling and high-speed accretion disk winds injected at small radii.

5.1 The Role of Stellar Feedback

Absent AGN feedback the properties of the gas disk inside pc are as follows. Gas cools efficiently and collapses in a mini-starburst until sufficient young stars are formed to maintain (mostly via radiation pressure-driven turbulence), leading to dispersions in a cold nuclear molecular disk. As in previous simulations which adopt highly simplified sub-grid models of the ISM (Hopkins & Quataert, 2010a), the disk develops large-amplitude modes in gas and stars, and resonant angular momentum transfer between gas and stellar disks drives rapid inflow of gas, with accretion rates of yr at pc. This agrees well with the analytic (Hopkins & Quataert, 2011a) predictions for “gravitational torque”-driven accretion. In contrast, the Bondi-Hoyle, viscous, or ballistic accretion rate estimators fail to capture the simulation results and are not appropriate for the regimes simulated here, in which much of the gas resides in a rotationally supported (albeit geometrically quite thick) disk (Fig. 6).

Stellar feedback does operate somewhat differently in galactic nuclei, as opposed to larger galactic radii, because the local dynamical time is  Myr. Young massive stars are sheared into an un-clustered mass distribution (e.g. executing hundreds of orbits at pc) before they explode. Rather than local, Jeans-scale clouds evolving independently, we should think of the disk as a coherently evolving, disky “star cluster” (see e.g. Thompson et al., 2005). On longer timescales, this leads to episodic “burst-quench” cycles on small scales, studied in Torrey et al. (2016).

5.2 The Role of AGN Feedback

Our calculations demonstrate that high-velocity winds from the central pc with momentum fluxes suggested by observations (e.g., Cimatti et al. 2013; Tombesi et al. 2015; Harrison et al. 2014; Zakamska et al. 2015) have a dramatic effect on the circum-BH ISM. In particular, such winds can evacuate gas from the circum-BH disk (see Figs. 1-2). This suppresses the star formation rate and black hole accretion rate in the galactic nucleus by a factor of and enhances the gas outflow rate at pc by a comparable factor (Figs. 3 & 5), also similar to observations in at least some systems with powerful on-going outflows (Shimizu et al., 2015; Guillard et al., 2015; Alatalo et al., 2015). As expected, the amount of BH growth required to produce this level of feedback and evacuate gas from the central regions depends inversely on the mass and momentum-loading of the BH accretion-disk winds. The amount of hot gas generated and its “venting” depend on the initial wind velocities. Our simulations thus provide support for models in which luminous AGN significantly disrupt the ISM of their host galaxies, at least on scales pc. Our simulations also specifically support the hypothesis that luminous AGN may play a key role in driving galaxy-scale outflows from gas-rich galactic nuclei.

In the plane of the circum-BH disk, the AGN winds decelerate as material is entrained into expanding rings/shells. In the polar direction, however, the galaxy-scale outflows powered by the AGN retain high velocities () as they reach kpc scales; although not isotropic, the opening angle for the high-velocity outflow is large ( of the sky). A detailed comparison with observations is outside the scope of this work, because we consider only one initial condition, and do not model the large-scale galaxy beyond pc on which many AGN-driven outflows are observed (although a follow-up study designed to compare in detail with these observations is in progress). However, the broad range of velocities present simultaneously within the same system is consistent with outflows observed outside of accretion disk scales in molecular and ionized gas (see e.g. Tombesi et al., 2015; Harrison et al., 2014; Zakamska & Greene, 2014; Alatalo, 2015). Similarly, we identify outflowing material across a broad range of temperatures and spatial scales, from molecular to K gas, and from scales of pc. Much of the material in the disk plane is accelerated to low velocities and will not escape, but entrains a large mass (comparable to or larger than starburst-driven winds). At least some energy of shocked AGN winds is converted into work, allowing the outflow momentum to reach in some cases (again, qualitatively consistent with observations; see Borguet et al., 2012; Cimatti et al., 2013; Cicone et al., 2014; Harrison et al., 2014; Zakamska & Greene, 2014). This is dictated largely by the geometry of the surrounding ISM, rather than the AGN wind at small radii. In particular, simulations with isotropically directed AGN winds on small scales give similar results to our default calculations that utilize primarily planar winds (this highlights that once the AGN wind shocks the gas follows the “path of least resistance” in the polar direction independent of exactly how the wind is initially directed). Understanding whether the outflows we find will be confined or halted by the galactic ISM or will continue to escape out of the galaxy will require galaxy-scale simulations. It is important to stress that the present calculations are not well-suited for addressing this question because our idealized initial conditions do not have, e.g., a gaseous halo or the nuclear warps/mis-alignments seen in both simulations and observations of galactic nuclei.

In our calculations, AGN-driven outflows also have a dramatic impact on obscuration of the AGN itself. AGN winds evacuate the polar region to allow a fully un-obscured view of the BH. AGN winds thus self-consistently produce a torus-like morphology (see, in particular, the lower right panel of Fig. 1). Quantitatively, we find a broad column density distribution from , in reasonable agreement with observations (e.g. Malizia et al., 2009; Treister et al., 2009; Risaliti et al., 1999; Burlon et al., 2011, and references therein). The inhomogeneous nature of the ISM also inevitably introduces large (dex) variation in obscuring columns even at roughly fixed polar angle – similar to observational suggestions of “clumpy” torii (Risaliti et al., 2002; Mason et al., 2006; Sánchez et al., 2006; Nenkova et al., 2008; Ramos Almeida et al., 2009; Hönig & Kishimoto, 2010; Deo et al., 2011).

We find that Compton heating/cooling from the AGN produces weak effects on these scales, consistent with previous studies (Choi et al., 2012, 2014, 2015; Park et al., 2014).

5.3 Future Work: Other Scales & Forms of Feedback

These simulations are a first exploration of the interaction between AGN and stellar feedback on scales between the BH accretion disk and galaxy. We focused on these scales because they are relatively unexplored and yet critical for understanding BH growth and the impact of AGN winds and radiation on the ISM. It is, however, also clearly important to extend our models to cover a broader range of spatial scales. On smaller scales, understanding the origin of AGN winds and radiation and their “escape” from the accretion disk is critical for setting the magnitude and geometry of AGN feedback on pc scales. On galactic scales, we need to understand how the outflows found here interact with the galaxy ISM on long timescales: in particular how this changes galaxy star formation histories and regulates future episodes of BH inflow. This is necessary to determine the effects of feedback on BH-host galaxy correlations. Also, many observations of galaxy-scale, AGN-driven winds suggest large momentum-loading in the winds, with up to (see Sturm et al., 2011; Faucher-Giguere & Quataert, 2012; Cicone et al., 2014). It will be particularly interesting to see whether these observations are consistent with a model in which this large momentum flux is generated on accretion disk scales (our v5000_hiP model), e.g., by super-Eddington accretion, or whether they require additional large-scale effects such as confinement (and buildup of a pressure-driven bubble) of radiation or hot shocked gas in the galaxies ISM. This could, e.g., be produced by misalignment between the nuclear scale disk and the galaxy ISM as a whole.

The two AGN feedback mechanisms we have studied here (fast AGN winds and Compton heating), are by no means exhaustive. In future work, we will extend this to include radiation pressure on both narrow lines and dust, photo-heating, and the effects of relativistic jets, all of which can act directly on gas both on the scales we model here but also on much larger scales up to kpc. We will also study the effects of different initial conditions (e.g. gas fraction and disk-to-BH mass ratio).


We thank Todd Thompson for helpful discussions, and the anonymous referee for helpful suggestions. Support for PFH was provided by the Gordon and Betty Moore Foundation through Grant #776 to the Caltech Moore Center for Theoretical Cosmology and Physics, an Alfred P. Sloan Research Fellowship, NASA ATP Grant NNX14AH35G, and NSF Collaborative Research Grant #1411920. Numerical calculations were run on the Caltech compute cluster “Zwicky” (NSF MRI award #PHY-0960291) and allocation TG-AST130039 granted by the Extreme Science and Engineering Discovery Environment (XSEDE) supported by the NSF. EQ is supported in part by NASA ATP Grant 12-ATP12-0183, the David and Lucile Packard Foundation, and a Simons Investigator Award from the Simons Foundation. CAFG was supported by a fellowship from the Miller Institute for Basic Research in Science, by NASA through Einstein Postdoctoral Fellowship Award PF3-140106 and grant 10-ATP10-0187, by NSF through grant AST-1412836, and by Northwestern University funds.


  • Agertz et al. (2013) Agertz, O., Kravtsov, A. V., Leitner, S. N., & Gnedin, N. Y. 2013, ApJ, 770, 25
  • Agertz et al. (2007) Agertz, O., et al. 2007, MNRAS, 380, 963
  • Alatalo (2015) Alatalo, K. 2015, ApJL, 801, L17
  • Alatalo et al. (2015) Alatalo, K., et al. 2015, ApJ, 798, 31
  • Aller & Richstone (2007) Aller, M. C., & Richstone, D. O. 2007, ApJ, 665, 120
  • Antonuccio-Delogu & Silk (2008) Antonuccio-Delogu, V., & Silk, J. 2008, MNRAS, 389, 1750
  • Bacon et al. (2001) Bacon, R., Emsellem, E., Combes, F., Copin, Y., Monnet, G., & Martin, P. 2001, A&A, 371, 409
  • Basko et al. (1974) Basko, M. M., Sunyaev, R. A., & Titarchuk, L. G. 1974, A&A, 31, 249
  • Batcheldor et al. (2007) Batcheldor, D., Tadhunter, C., Holt, J., Morganti, R., O’Dea, C. P., Axon, D. J., & Koekemoer, A. 2007, ApJ, 661, 70
  • Bauer & Springel (2012) Bauer, A., & Springel, V. 2012, MNRAS, 423, 3102
  • Bautista et al. (2010) Bautista, M. A., Dunn, J. P., Arav, N., Korista, K. T., Moe, M., & Benn, C. 2010, ApJ, 713, 25
  • Begelman (1985) Begelman, M. C. 1985, ApJ, 297, 492
  • Blandford & Payne (1982) Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883
  • Borguet et al. (2012) Borguet, B. C. J., Arav, N., Edmonds, D., Chamberlain, C., & Benn, C. 2012, ApJ, in press, arXiv:1211.6250
  • Burlon et al. (2011) Burlon, D., Ajello, M., Greiner, J., Comastri, A., Merloni, A., & Gehrels, N. 2011, ApJ, 728, 58
  • Chang et al. (1987) Chang, C. A., Schiano, A. V. R., & Wolfe, A. M. 1987, ApJ, 322, 180
  • Chatterjee et al. (2015) Chatterjee, S., Newman, J. A., Jeltema, T., Myers, A. D., Aird, J., Coil, A. L., Cooper, M., Finoguenov, A., Laird, E., Montero-Dorta, A., Nandra, K., Willmer, C., & Yan, R. 2015, PASP, 127, 716
  • Choi et al. (2014) Choi, E., Naab, T., Ostriker, J. P., Johansson, P. H., & Moster, B. P. 2014, MNRAS, 442, 440
  • Choi et al. (2012) Choi, E., Ostriker, J. P., Naab, T., & Johansson, P. H. 2012, ApJ, in press, arXiv:1205.2082
  • Choi et al. (2015) Choi, E., Ostriker, J. P., Naab, T., Oser, L., & Moster, B. P. 2015, MNRAS, 449, 4105
  • Cicone et al. (2014) Cicone, C., et al. 2014, A&A, 562, A21
  • Cimatti et al. (2013) Cimatti, A., Brusa, M., Talia, M., Mignoli, M., Rodighiero, G., Kurk, J., Cassata, P., Halliday, C., Renzini, A., & Daddi, E. 2013, ApJL, 779, L13
  • Ciotti et al. (2010) Ciotti, L., Ostriker, J. P., & Proga, D. 2010, ApJ, 717, 708
  • Coil et al. (2011) Coil, A. L., Weiner, B. J., Holz, D. E., Cooper, M. C., Yan, R., & Aird, J. 2011, ApJ, 743, 46
  • Creasey et al. (2011) Creasey, P., Theuns, T., Bower, R. G., & Lacey, C. G. 2011, MNRAS, 415, 3706
  • Crenshaw et al. (2000) Crenshaw, D. M., et al. 2000, AJ, 120, 1731
  • Croft et al. (2009) Croft, R. A. C., Di Matteo, T., Springel, V., & Hernquist, L. 2009, MNRAS, 400, 43
  • Croton et al. (2006) Croton, D. J., et al. 2006, MNRAS, 365, 11
  • Cullen & Dehnen (2010) Cullen, L., & Dehnen, W. 2010, MNRAS, 408, 669
  • de Kool et al. (2001) de Kool, M., et al. 2001, ApJ, 548, 609
  • Debuhr et al. (2011) Debuhr, J., Quataert, E., & Ma, C. 2011, MNRAS, 412, 1341
  • Debuhr et al. (2010) Debuhr, J., Quataert, E., Ma, C., & Hopkins, P. 2010, MNRAS, 406, L55
  • Debuhr et al. (2012) Debuhr, J., Quataert, E., & Ma, C.-P. 2012, MNRAS, 420, 2221
  • Dehnen & Aly (2012) Dehnen, W., & Aly, H. 2012, MNRAS, 425, 1068
  • Deo et al. (2011) Deo, R. P., Richards, G. T., Nikutta, R., Elitzur, M., Gallagher, S. C., Ivezić, Ž., & Hines, D. 2011, ApJ, 729, 108
  • Di Matteo et al. (2008) Di Matteo, T., Colberg, J., Springel, V., Hernquist, L., & Sijacki, D. 2008, ApJ, 676, 33
  • Di Matteo et al. (2005) Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604
  • Diamond-Stanic & Rieke (2012) Diamond-Stanic, A. M., & Rieke, G. H. 2012, ApJ, 746, 168
  • Dubois et al. (2014) Dubois, Y., Volonteri, M., & Silk, J. 2014, MNRAS, 440, 1590
  • Dunn et al. (2010) Dunn, J. P., et al. 2010, ApJ, 709, 611
  • Durier & Dalla Vecchia (2012) Durier, F., & Dalla Vecchia, C. 2012, MNRAS, 419, 465
  • Elvis (2000) Elvis, M. 2000, ApJ, 545, 63
  • Faucher-Giguere et al. (2015) Faucher-Giguere, C.-A., Hopkins, P. F., Keres, D., Muratov, A. L., Quataert, E., & Murray, N. 2015, MNRAS, 449, 987
  • Faucher-Giguere & Quataert (2012) Faucher-Giguere, C.-A., & Quataert, E. 2012, MNRAS, in press, arXiv:1204.2547
  • Faucher-Giguère et al. (2012) Faucher-Giguère, C.-A., Quataert, E., & Murray, N. 2012, MNRAS, 420, 1347
  • Feoli & Mancini (2009) Feoli, A., & Mancini, L. 2009, ApJ, 703, 1502
  • Ferrarese & Merritt (2000) Ferrarese, L., & Merritt, D. 2000, ApJL, 539, L9
  • Feruglio et al. (2010) Feruglio, C., Maiolino, R., Piconcelli, E., Menci, N., Aussel, H., Lamastra, A., & Fiore, F. 2010, A&A, 518, L155
  • Fischer et al. (2010) Fischer, J., et al. 2010, A&A, 518, L41
  • Gabel et al. (2006) Gabel, J. R., Arav, N., & Kim, T.-S. 2006, ApJ, 646, 742
  • Gammie (2001) Gammie, C. F. 2001, ApJ, 553, 174
  • Gan et al. (2014) Gan, Z., Yuan, F., Ostriker, J. P., Ciotti, L., & Novak, G. S. 2014, ApJ, 789, 150
  • Ganguly et al. (2007) Ganguly, R., Brotherton, M. S., Cales, S., Scoggins, B., Shang, Z., & Vestergaard, M. 2007, ApJ, 665, 990
  • Gebhardt et al. (2000) Gebhardt, K., et al. 2000, ApJL, 539, L13
  • Gingold & Monaghan (1983) Gingold, R. A., & Monaghan, J. J. 1983, MNRAS, 204, 715
  • Granato et al. (2004) Granato, G. L., De Zotti, G., Silva, L., Bressan, A., & Danese, L. 2004, ApJ, 600, 580
  • Greene et al. (2011) Greene, J. E., Zakamska, N. L., Ho, L. C., & Barth, A. J. 2011, ApJ, 732, 9
  • Greene et al. (2012) Greene, J. E., Zakamska, N. L., & Smith, P. S. 2012, ApJ, 746, 86
  • Guillard et al. (2015) Guillard, P., Boulanger, F., Lehnert, M. D., Pineau des Forêts, G., Combes, F., Falgarone, E., & Bernard-Salas, J. 2015, A&A, 574, A32
  • Hamann et al. (2011) Hamann, F., Kanekar, N., Prochaska, J. X., Murphy, M. T., Ellison, S., Malec, A. L., Milutinovic, N., & Ubachs, W. 2011, MNRAS, 410, 1957
  • Harrison et al. (2014) Harrison, C. M., Alexander, D. M., Mullaney, J. R., & Swinbank, A. M. 2014, MNRAS, 441, 3306
  • Harrison et al. (2015) Harrison, C. M., Thomson, A. P., Alexander, D. M., Bauer, F. E., Edge, A. C., Hogan, M. T., Mullaney, J. R., & Swinbank, A. M. 2015, ApJ, 800, 45
  • Hernquist (1990) Hernquist, L. 1990, ApJ, 356, 359
  • Hobbs et al. (2011) Hobbs, A., Nayakshin, S., Power, C., & King, A. 2011, MNRAS, 413, 2633
  • Hönig & Kishimoto (2010) Hönig, S. F., & Kishimoto, M. 2010, A&A, 523, A27
  • Hopkins (2010) Hopkins, P. F. 2010, MNRAS, in press, arXiv:1009.4702 [astro-ph]
  • Hopkins (2012) —. 2012, MNRAS, 420, L8
  • Hopkins (2013) —. 2013, MNRAS, 428, 2840
  • Hopkins (2015) —. 2015, MNRAS, 450, 53
  • Hopkins et al. (2013a) Hopkins, P. F., Cox, T. J., Hernquist, L., Narayanan, D., Hayward, C. C., & Murray, N. 2013a, MNRAS, 430, 1901
  • Hopkins et al. (2008a) Hopkins, P. F., Cox, T. J., Kereš, D., & Hernquist, L. 2008a, ApJS, 175, 390
  • Hopkins et al. (2012a) Hopkins, P. F., Hayward, C. C., Narayanan, D., & Hernquist, L. 2012a, MNRAS, 420, 320
  • Hopkins et al. (2005a) Hopkins, P. F., Hernquist, L., Cox, T. J., Di Matteo, T., Martini, P., Robertson, B., & Springel, V. 2005a, ApJ, 630, 705
  • Hopkins et al. (2005b) Hopkins, P. F., Hernquist, L., Cox, T. J., Di Matteo, T., Robertson, B., & Springel, V. 2005b, ApJ, 630, 716
  • Hopkins et al. (2005c) —. 2005c, ApJ, 632, 81
  • Hopkins et al. (2006a) —. 2006a, ApJS, 163, 1
  • Hopkins et al. (2008b) Hopkins, P. F., Hernquist, L., Cox, T. J., & Kereš, D. 2008b, ApJS, 175, 356
  • Hopkins et al. (2007) Hopkins, P. F., Hernquist, L., Cox, T. J., Robertson, B., & Krause, E. 2007, ApJ, 669, 67
  • Hopkins et al. (2012b) Hopkins, P. F., Hernquist, L., Hayward, C. C., & Narayanan, D. 2012b, MNRAS, 425, 1121
  • Hopkins et al. (2005d) Hopkins, P. F., Hernquist, L., Martini, P., Cox, T. J., Robertson, B., Di Matteo, T., & Springel, V. 2005d, ApJL, 625, L71
  • Hopkins et al. (2009a) Hopkins, P. F., Hickox, R., Quataert, E., & Hernquist, L. 2009a, MNRAS, 398, 333
  • Hopkins et al. (2014) Hopkins, P. F., Keres, D., Onorbe, J., Faucher-Giguere, C.-A., Quataert, E., Murray, N., & Bullock, J. S. 2014, MNRAS, 445, 581
  • Hopkins et al. (2013b) Hopkins, P. F., Kereš, D., Murray, N., Hernquist, L., Narayanan, D., & Hayward, C. C. 2013b, MNRAS, 433, 78
  • Hopkins et al. (2009b) Hopkins, P. F., Murray, N., & Thompson, T. A. 2009b, MNRAS, 398, 303
  • Hopkins et al. (2006b) Hopkins, P. F., Narayan, R., & Hernquist, L. 2006b, ApJ, 643, 641
  • Hopkins et al. (2013c) Hopkins, P. F., Narayanan, D., & Murray, N. 2013c, MNRAS, 432, 2647
  • Hopkins & Quataert (2010a) Hopkins, P. F., & Quataert, E. 2010a, MNRAS, 407, 1529
  • Hopkins & Quataert (2010b) —. 2010b, MNRAS, 405, L41
  • Hopkins & Quataert (2011a) —. 2011a, MNRAS, 415, 1027
  • Hopkins & Quataert (2011b) —. 2011b, MNRAS, 411, L61
  • Hopkins et al. (2011) Hopkins, P. F., Quataert, E., & Murray, N. 2011, MNRAS, 417, 950
  • Hopkins et al. (2012c) —. 2012c, MNRAS, 421, 3522
  • Hopkins et al. (2012d) —. 2012d, MNRAS, 421, 3488
  • Hopkins et al. (2006c) Hopkins, P. F., Somerville, R. S., Hernquist, L., Cox, T. J., Robertson, B., & Li, Y. 2006c, ApJ, 652, 864
  • Hopkins et al. (2015) Hopkins, P. F., Torrey, P., Faucher-Giguere, C.-A., Quataert, E., & Murray, N. 2015, MNRAS, submitted, arXiv:1504.05209
  • Humphrey et al. (2010) Humphrey, A., et al. 2010, MNRAS, 408, L1
  • Hutchings & Thomas (2000) Hutchings, R. M., & Thomas, P. A. 2000, MNRAS, 319, 721
  • King (2003) King, A. 2003, ApJL, 596, L27
  • Konigl & Kartje (1994) Konigl, A., & Kartje, J. F. 1994, ApJ, 434, 446
  • Kormendy et al. (2011) Kormendy, J., Bender, R., & Cornell, M. E. 2011, Nature, 469, 374
  • Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511
  • Krongold et al. (2007) Krongold, Y., Nicastro, F., Elvis, M., Brickhouse, N., Binette, L., Mathur, S., & Jiménez-Bailón, E. 2007, ApJ, 659, 1022
  • Kroupa (2002) Kroupa, P. 2002, Science, 295, 82
  • Krumholz & Gnedin (2011) Krumholz, M. R., & Gnedin, N. Y. 2011, ApJ, 729, 36
  • LaMassa et al. (2013) LaMassa, S. M., Heckman, T. M., Ptak, A., & Urry, C. M. 2013, ApJL, 765, L33
  • Laor et al. (1997) Laor, A., Fiore, F., Elvis, M., Wilkes, B. J., & McDowell, J. C. 1997, ApJ, 477, 93
  • Lapi et al. (2006) Lapi, A., Shankar, F., Mao, J., Granato, G. L., Silva, L., De Zotti, G., & Danese, L. 2006, ApJ, 650, 42
  • Ma et al. (2016) Ma, X., Hopkins, P. F., Faucher-Giguère, C.-A., Zolman, N., Muratov, A. L., Kereš, D., & Quataert, E. 2016, MNRAS, 456, 2140
  • Ma et al. (2015) Ma, X., Kasen, D., Hopkins, P. F., Faucher-Giguère, C.-A., Quataert, E., Kereš, D., & Murray, N. 2015, MNRAS, 453, 960
  • Magorrian et al. (1998) Magorrian, J., et al. 1998, AJ, 115, 2285
  • Malizia et al. (2009) Malizia, A., Stephen, J. B., Bassani, L., Bird, A. J., Panessa, F., & Ubertini, P. 2009, MNRAS, 399, 944
  • Mason et al. (2006) Mason, R. E., Geballe, T. R., Packham, C., Levenson, N. A., Elitzur, M., Fisher, R. S., & Perlman, E. 2006, ApJ, 640, 612
  • McKernan & Yaqoob (1998) McKernan, B., & Yaqoob, T. 1998, ApJL, 501, L29+
  • Menci et al. (2003) Menci, N., Cavaliere, A., Fontana, A., Giallongo, E., Poli, F., & Vittorini, V. 2003, ApJL, 587, L63
  • Miller et al. (2008) Miller, L., Turner, T. J., & Reeves, J. N. 2008, A&A, 483, 437
  • Moe et al. (2009) Moe, M., Arav, N., Bautista, M. A., & Korista, K. T. 2009, ApJ, 706, 525
  • Morris & Monaghan (1997) Morris, J. P., & Monaghan, J. J. 1997, Journal of Computational Physics, 136, 41
  • Muratov et al. (2015) Muratov, A. L., Kereš, D., Faucher-Giguère, C.-A., Hopkins, P. F., Quataert, E., & Murray, N. 2015, MNRAS, 454, 2691
  • Murray et al. (1995) Murray, N., Chiang, J., Grossman, S. A., & Voit, G. M. 1995, ApJ, 451, 498
  • Murray et al. (2005) Murray, N., Quataert, E., & Thompson, T. A. 2005, ApJ, 618, 569
  • Mushotzky et al. (2014) Mushotzky, R. F., Shimizu, T. T., Meléndez, M., & Koss, M. 2014, ApJL, 781, L34
  • Narayan & Yi (1995) Narayan, R., & Yi, I. 1995, ApJ, 444, 231
  • Narayanan et al. (2012) Narayanan, D., Krumholz, M. R., Ostriker, E. C., & Hernquist, L. 2012, MNRAS, 421, 3127
  • Nenkova et al. (2008) Nenkova, M., Sirocky, M. M., Nikutta, R., Ivezić, Ž., & Elitzur, M. 2008, ApJ, 685, 160
  • Nesvadba et al. (2010) Nesvadba, N. P. H., Boulanger, F., Salomé, P., Guillard, P., Lehnert, M. D., Ogle, P., Appleton, P., Falgarone, E., & Pineau Des Forets, G. 2010, A&A, 521, A65
  • Novak et al. (2011) Novak, G. S., Ostriker, J. P., & Ciotti, L. 2011, ApJ, 737, 26
  • Oñorbe et al. (2015) Oñorbe, J., Boylan-Kolchin, M., Bullock, J. S., Hopkins, P. F., Kereš, D., Faucher-Giguère, C.-A., Quataert, E., & Murray, N. 2015, MNRAS, 454, 2092
  • Ogle et al. (1999) Ogle, P. M., Cohen, M. H., Miller, J. S., Tran, H. D., Goodrich, R. W., & Martel, A. R. 1999, ApJS, 125, 1
  • Ostriker et al. (2010) Ostriker, J. P., Choi, E., Ciotti, L., Novak, G. S., & Proga, D. 2010, ApJ, 722, 642
  • Park et al. (2014) Park, K., Ricotti, M., Di Matteo, T., & Reynolds, C. S. 2014, MNRAS, 445, 2325
  • Price (2008) Price, D. J. 2008, Journal of Computational Physics, 227, 10040
  • Price (2012) —. 2012, MNRAS, 420, L33
  • Price & Federrath (2010) Price, D. J., & Federrath, C. 2010, MNRAS, 406, 1659
  • Prochaska & Hennawi (2009) Prochaska, J. X., & Hennawi, J. F. 2009, ApJ, 690, 1558
  • Proga (2000) Proga, D. 2000, ApJ, 538, 684
  • Proga (2007) —. 2007, ApJ, 661, 693
  • Ramos Almeida et al. (2009) Ramos Almeida, C., et al. 2009, ApJ, 702, 1127
  • Read & Hayfield (2012) Read, J. I., & Hayfield, T. 2012, MNRAS, 422, 3037
  • Risaliti et al. (2002) Risaliti, G., Elvis, M., & Nicastro, F. 2002, ApJ, 571, 234
  • Risaliti et al. (1999) Risaliti, G., Maiolino, R., & Salvati, M. 1999, ApJ, 522, 157
  • Rosswog et al. (2000) Rosswog, S., Davies, M. B., Thielemann, F.-K., & Piran, T. 2000, A&A, 360, 171
  • Rupke & Veilleux (2011) Rupke, D. S. N., & Veilleux, S. 2011, ApJL, 729, L27+
  • Saitoh et al. (2008) Saitoh, T. R., Daisaka, H., Kokubo, E., Makino, J., Okamoto, T., Tomisaka, K., Wada, K., & Yoshida, N. 2008, PASJ, 60, 667
  • Saitoh & Makino (2009) Saitoh, T. R., & Makino, J. 2009, ApJL, 697, L99
  • Salucci et al. (1999) Salucci, P., Szuszkiewicz, E., Monaco, P., & Danese, L. 1999, MNRAS, 307, 637
  • Sánchez et al. (2006) Sánchez, F. M., Davies, R. I., Eisenhauer, F., Tacconi, L. J., Genzel, R., & Sternberg, A. 2006, A&A, 454, 481
  • Sanders et al. (1988) Sanders, D. B., Soifer, B. T., Elias, J. H., Madore, B. F., Matthews, K., Neugebauer, G., & Scoville, N. Z. 1988, ApJ, 325, 74
  • Sazonov et al. (2005) Sazonov, S. Y., Ostriker, J. P., Ciotti, L., & Sunyaev, R. A. 2005, MNRAS, 358, 168
  • Sazonov et al. (2004) Sazonov, S. Y., Ostriker, J. P., & Sunyaev, R. A. 2004, MNRAS, 347, 144
  • Scannapieco & Oh (2004) Scannapieco, E., & Oh, S. P. 2004, ApJ, 608, 62
  • Schmidt & Hines (1999) Schmidt, G. D., & Hines, D. C. 1999, ApJ, 512, 125
  • Shimizu et al. (2015) Shimizu, T. T., Mushotzky, R. F., Meléndez, M., Koss, M., & Rosario, D. J. 2015, MNRAS, 452, 1841
  • Sijacki et al. (2012) Sijacki, D., Vogelsberger, M., Keres, D., Springel, V., & Hernquist, L. 2012, MNRAS, 424, 2999
  • Silk (2005) Silk, J. 2005, MNRAS, 364, 1337
  • Silk & Rees (1998) Silk, J., & Rees, M. J. 1998, A&A, 331, L1
  • Soltan (1982) Soltan, A. 1982, MNRAS, 200, 115
  • Somerville et al. (2008) Somerville, R. S., Hopkins, P. F., Cox, T. J., Robertson, B. E., & Hernquist, L. 2008, MNRAS, 391, 481
  • Spitzer (1962) Spitzer, L. 1962, Physics of Fully Ionized Gases; New York: Interscience (New York: Interscience)
  • Springel (2005) Springel, V. 2005, MNRAS, 364, 1105
  • Springel et al. (2005) Springel, V., Di Matteo, T., & Hernquist, L. 2005, MNRAS, 361, 776
  • Springel & Hernquist (2002) Springel, V., & Hernquist, L. 2002, MNRAS, 333, 649
  • Springel & Hernquist (2003) —. 2003, MNRAS, 339, 289
  • Steenbrugge et al. (2005) Steenbrugge, K. C., et al. 2005, A&A, 434, 569
  • Steinborn et al. (2015) Steinborn, L. K., Dolag, K., Hirschmann, M., Prieto, M. A., & Remus, R.-S. 2015, MNRAS, 448, 1504
  • Stern et al. (2015) Stern, J., Faucher-Giguere, C.-A., Zakamska, N. L., & Hennawi, J. F. 2015, MNRAS, in press, arxiv:1510.07690
  • Sturm et al. (2011) Sturm, E., et al. 2011, ApJL, 733, L16+
  • Sunyaev & Churazov (1996) Sunyaev, R. A., & Churazov, E. M. 1996, Astronomy Letters, 22, 648
  • Thompson et al. (2005) Thompson, T. A., Quataert, E., & Murray, N. 2005, ApJ, 630, 167
  • Tombesi et al. (2013) Tombesi, F., Cappi, M., Reeves, J. N., Nemmen, R. S., Braito, V., Gaspari, M., & Reynolds, C. S. 2013, MNRAS, 430, 1102
  • Tombesi et al. (2010) Tombesi, F., Cappi, M., Reeves, J. N., Palumbo, G. G. C., Yaqoob, T., Braito, V., & Dadina, M. 2010, A&A, 521, A57
  • Tombesi et al. (2015) Tombesi, F., Meléndez, M., Veilleux, S., Reeves, J. N., González-Alfonso, E., & Reynolds, C. S. 2015, Nature, 519, 436
  • Torrey et al. (2016) Torrey, P., Hopkins, P. F., Faucher-Giguère, C.-A., Vogelsberger, M., Quataert, E., Kereš, D., & Murray, N. 2016, MNRAS, in press, arXiv:1601.07186
  • Tortora et al. (2009) Tortora, C., Antonuccio-Delogu, V., Kaviraj, S., Silk, J., Romeo, A. D., & Becciani, U. 2009, MNRAS, 396, 61
  • Treister et al. (2009) Treister, E., Urry, C. M., & Virani, S. 2009, The Astrophysical Journal, 696, 110
  • Tremaine (1995) Tremaine, S. 1995, AJ, 110, 628
  • Tremonti et al. (2007) Tremonti, C. A., Moustakas, J., & Diamond-Stanic, A. M. 2007, ApJL, 663, L77
  • Turner et al. (2008) Turner, T. J., Reeves, J. N., Kraemer, S. B., & Miller, L. 2008, A&A, 483, 161
  • van de Voort et al. (2015) van de Voort, F., Quataert, E., Hopkins, P. F., Kereš, D., & Faucher-Giguère, C.-A. 2015, MNRAS, 447, 140
  • Volonteri et al. (2006) Volonteri, M., Salvaterra, R., & Haardt, F. 2006, MNRAS, 373, 121
  • Weymann et al. (1981) Weymann, R. J., Carswell, R. F., & Smith, M. G. 1981, ARA&A, 19, 41
  • Wiersma et al. (2009a) Wiersma, R. P. C., Schaye, J., & Smith, B. D. 2009a, MNRAS, 393, 99
  • Wiersma et al. (2009b) Wiersma, R. P. C., Schaye, J., Theuns, T., Dalla Vecchia, C., & Tornatore, L. 2009b, MNRAS, 399, 574
  • Wild et al. (2009) Wild, V., Walcher, C. J., Johansson, P. H., Tresse, L., Charlot, S., Pollo, A., Le Fèvre, O., & de Ravel, L. 2009, MNRAS, 395, 144
  • Yu & Tremaine (2002) Yu, Q., & Tremaine, S. 2002, MNRAS, 335, 965
  • Zakamska & Greene (2014) Zakamska, N. L., & Greene, J. E. 2014, MNRAS, 442, 784
  • Zakamska et al. (2015) Zakamska, N. L., Hamann, F., Pâris, I., Brandt, W. N., Greene, J. E., Strauss, M. A., Villforth, C., Wylezalek, D., Alexandroff, R. M., & Ross, N. P. 2015, MNRAS, in press, arxiv:1512.02642

Appendix A Black Hole Feedback Implementation

a.1 Broad Absorption Line Quasar Winds

Bright quasars often have BAL winds with velocities of . We model these in the most direct manner possible: when a gas particle is accreted, a fraction is actually assumed to accrete on the BH while the remaining is blown out as a BAL wind with velocity . There is both observational (Schmidt & Hines, 1999; Ogle et al., 1999) and theoretical (Murray et al., 1995) evidence that BAL winds are approximately planar (in or slightly out of the accretion disk plane). Assuming that the angular momentum vector of the small-scale accretion disk is correlated with that of the sub-pc accreting material, then this corresponds to directing along the radial vector from the BH. On an accretion event we therefore take (for the accreted gas particle) , apply the “kick” , and hold (internal energy per unit mass) constant. In model v5000_iso we instead assign the wind direction randomly. Previous work has shown that randomly directed winds yield results similar to planar winds, but require somewhat larger wind momentum fluxes to achieve the same feedback on the ambient gas (Debuhr et al., 2012).

Two parameters must be chosen: the initial outflow mass loading and velocity ; observations and theoretical models suggest values of order , (see e.g. Moe et al., 2009; Dunn et al., 2010; Hamann et al., 2011; Borguet et al., 2012, and references therein)

Equivalently, we can translate these parameters into the wind momentum and energy-loading. Since BAL winds are believed to be driven by line radiation pressure in the accretion disk, the available momentum flux is , where is the luminosity (with the radiative efficiency). The “initial” wind momentum and energy are and respectively; thus the energy and momentum-loading are


Note for , we recover the canonical adopted in previous simulations with purely thermal AGN feedback (e.g. Di Matteo et al., 2005; Hopkins et al., 2005d).

a.2 Compton Heating/Cooling

The radiation field of the BH will also Compton heat/cool gas in its vicinity. As discussed in Sazonov et al. (2004, 2005), this effect is nearly independent of obscuration: Compton heating is entirely dominated by photons with energies keV (for which we can usually safely ignore obscuration) and Compton cooling by the bolometric luminosity in lower-energy photons (re-distributed, but not, in integral, altered by obscuration). As such even Compton-thick columns result in factor changes in the heating/cooling rates. We therefore neglect obscuration and assume the radiation field is isotropic, so that the X-ray/bolometric flux from the AGN on all particles is given by , with Compton temperature  K as calculated in Sazonov et al. (2004) for a broad range of observed QSO SED shapes.666We propagate this flux through the gravity tree, since it follows an inverse-square law when we can neglect obscuration. This makes it trivial to apply the appropriate flux to arbitrary particle numbers, geometries, and numbers of black holes. In the cooling function, we add the appropriate Compton heating and cooling terms.777As is standard, cooling is solved implicitly within this function in the regime where the heating/cooling times are short compared to the particle timesteps. Although Compton cooling depends explicitly on the free electron fraction, for the photon energies dominating heating (much greater than the ionization energy of hydrogen), we can safely approximate Compton heating of bound electrons as identical to free electrons (see e.g. Basko et al., 1974; Sunyaev & Churazov, 1996).

Finally, as shown in Faucher-Giguere & Quataert (2012), some care is needed at the highest temperatures: if the timescale for Coulomb collisions to transfer energy from ions to electrons is longer than the Compton or free-free cooling time of the electrons, this is the rate-limiting process and a two-temperature plasma develops. We therefore do not allow the Compton+free-free cooling rate to exceed the Coulomb energy transfer rate between ions and electrons calculated for an ion temperature in the limit where the electrons are efficiently cooling (see Spitzer, 1962; Narayan & Yi, 1995). It is important to note that AGN wind-shocked electrons are generally non-relativistic: either immediately post-shock (where most energy is in protons, with electron temperature ), or in later stages when competition between Compton cooling and Coulomb heating regulates the temperature.

Appendix B Variability on Un-Resolved Timescales

Figure 9: BH accretion rate vs. time as Fig. 5, for our “no BAL” simulation. We show the rate smoothed in a rolling yr window as in Fig. 5 (thick dashed), and smoothed in a  yr window (solid). There is variability on all resolved timescales; however, some of this owes to numerical effects, namely the assumption of instantaneous accretion of discrete gas particles. In Appendix B, we describe a model allowing variability even on unresolved (arbitrarily small) timescales.

We are able in these simulations to follow inflows to sub-pc scales. This, coupled to the fact that our accretion model is by definition discrete (i.e. particles are swallowed individually and instantaneously accreted onto the BH) leads to large variability on very small timescales, illustrated in Fig. 9. However, there are still several orders of magnitude between these scales and the BH event horizon, spanned by the Shakura-Sunyaev accretion () disk. Empirically, AGN exhibit variability on all observed timescales, corresponding to these unresolved spatial scales. Although we cannot resolve these scales, we can make a crude estimate of the effects of this variability by including a sub-grid power-spectrum of luminosity fluctuations and integrating over this to obtain the (modified) momentum flux in every resolved simulation timestep. We quantitatively implement this following the prescription in Hopkins & Quataert (2011a), integrating over a power spectrum with equal logarithmic power per logarithmic time interval, from the minimum resolved timestep down to the orbital time at the innermost stable circular orbit for a non-rotating BH. Performing such an experiment, we find almost no effect on our conclusions. Given the resolved dynamic range in the simulations, this additional variability occurs on extremely small timescales compared to the dynamical times of the outflow – the timescale over which feedback determines the equilibrium accretion rate. As such, other than adding the chosen random variance to the lightcurve on small timescales, this introduces (relatively) little dynamical effect.

Appendix C Numerical Tests

Figure 10: Disk structure as Fig. 4, but comparing different numerical methods (see § C: our “pressure-entropy” formulation of SPH developed in Hopkins (2013), which performs well in tests of fluid mixing instabilities while maintaining good conservation; a “density-entropy” formation (the “standard” GADGET SPH method); a run with a simplified artificial viscosity prescription and SPH smoothing kernel; and a run with the “pressure-entropy” formulation but ten times as many particles. The results appear robust to these variations.
Figure 11: Inflow and outflow properties as Fig. 5, comparing different numerical methods as Fig. 10. Again, the results are consistent. Convergence to the asymptotic level of outflow occurs more rapidly at high resolution and more slowly with the modified viscosity and kernel, but this is a result of the non-equilibrium initial conditions. The variations in the late-time accretion rates owe to the chaotic bursts as individual cold gas clumps sink to the center, so we only expect statistical convergence in this stage.

We now consider some tests of the robustness of the numerical methods used here. Figs. 10-11 repeat Figs. 4-5, but with varied numerical prescriptions. Our default simulations use the “pressure-entropy” SPH formulation described in Hopkins (2013), which is shown there to give dramatically improved results on in situations with fluid mixing around contact discontinuities (e.g. the Kelvin-Helmholtz and Rayleigh-Taylor instabilities) while retaining excellent conservation properties, and includes a number of additional improvements to the treatment of artificial viscosity (see Cullen & Dehnen, 2010), SPH smoothing kernel accuracy (Dehnen & Aly, 2012), and timestep communication relevant for treating extremely high Mach-number shocks (Saitoh & Makino, 2009; Durier & Dalla Vecchia, 2012).

To test whether these subtleties may be strongly influencing our results, we re-run our standard v5000 simulation instead using a “density-entropy” SPH formulation, as in Springel & Hernquist (2002) (the “standard” GADGET formulation of the SPH equations-of-motion). 888To ensure the simulations are otherwise exactly identical, we have had to re-run the v5000 simulation in the “pressure-entropy” case with a number of small modifications to the algorithm, and on an identical node configuration with pre-set values for certain random number calls. This is done for all tests in this section. For convenience we run the test cases at the particle number in the text. This produces a “surface tension” term at contact discontinuities that suppresses some fluid mixing instabilities, which has been the subject of much discussion in the literature (see Agertz et al., 2007; Read & Hayfield, 2012, and references therein). We also re-run the simulations with the pressure-entropy formulation, but adopting the much simpler and more numerically dissipative constant form of artificial viscosity from Gingold & Monaghan (1983) (which can significantly alter the behavior in sub-sonic turbulence; see Price 2012), and a greatly reduced-accuracy SPH smoothing kernel (a 32-neighbor cubic spline, as opposed to our standard 128-neighbor quintic spline). Together these variations produce the range of numerical effects which span the major SPH-grid code differences often discussed in the literature (see references above and Price & Federrath, 2010; Bauer & Springel, 2012; Sijacki et al., 2012). We also run a standard resolution test, increasing the number of particles by a factor of .

We see very little difference in the results in Figs. 10-11. Likewise there is relatively little difference in the column density distributions and gas morphologies. There are some slight differences in the phase diagrams, which correspond to the degree of fluid mixing along phase boundaries (the quantity most affected by these differences), but it is mostly at low temperatures where it does not have significant dynamical effects. This probably relates to the fact that in the numerical comparison studies discussed above, it is generally well-established that different methods agree well in the regimes of super-sonic turbulence and/or systems with dominant external forces. Moreover all of these changes preserve good energy and linear and angular momentum conservation, so to the extent that the outflow is primarily a simple momentum or energy-conserving “piston,” and the steady-state stellar feedback is the result of momentum input balancing runaway collapse (see Hopkins et al., 2011, 2012d), our conclusions should be robust.

In our resolution tests, we see quite good agreement in the BH accretion rate and wind momentum, up to the level of stochastic effects (random differences between simulations). In fact re-running our standard model with different random number seeds for the placement of the initial particles in the disk leads to comparable variations. In the outflow rates, we see the largest differences between runs, at times Myr. These also vary the most dramatically owing to purely random effects; at these times (order one dynamical time at pc), the outflow rate measured at pc has just begun to develop, so even small differences in the dynamics can appear significant. The high-resolution case develops a wind more quickly because fragmentation and star formation in the disk at these radii, being better resolved, proceeds faster, making the disk thicker and putting more material into a stellar-feedback driven fountain or wind which the BH wind is then able to entrain. Reassuringly, however, the other simulations soon “catch up” as their outer-disk star formation proceeds and more material is entrained, until the long-timescale outflow rates agree well.

Appendix D Resolving In-Shock Cooling: Numerical & Resolution Requirements

Figure 12: Shock tube tests designed to verify that numerical in-shock cooling does not significantly affect our results at our typical simulation resolution. See Appendix D for details. From left to right, three idealized shock tube problems are shown which represent a BAL wind encountering a cold, dense ISM at radii , respectively, and evolved to kyrs. Different colors correspond to the particle masses (as labeled); the highest-resolution case is comparable to our simulations in the main text. At lower resolutions there is noticeable numerical shock-broadening, however the post-shock temperatures are still well-converged (i.e. they are not affected by in-shock cooling, which would systematically change the post-shock temperatures at different resolutions).

The SPH hydro solver employed in this work captures the jump conditions associated with shocks over a finite width of order the kernel softening length. Gas particles passing through a numerically broadened shock can radiate away energy through traditional cooling channels. The post-shock gas temperature can therefore be numerically reduced via this “in-shock cooling” effect when shocks are broadened (e.g., Hutchings & Thomas, 2000; Creasey et al., 2011). In-shock cooling can become significant when the cooling timescale for gas moving through the shock is comparable to the resolution-dependent shock crossing timescale. For the case of BAL winds with input velocities of km s, the post-shock gas is expected to be heated to of order K where it will cool inefficiently owning to two-temperature plasma effects (Faucher-Giguere & Quataert, 2012). In the absence of efficient cooling channels for the post-shock gas, the outflows will remain energy conserving, efficiently driving outflows via PdV work on the ambient ISM. As a result, any significant amount of in-shock cooling can impact the post-shock gas temperature, and thus numerically modulate the quasar feedback efficiency studied in this paper. In this section, we investigate the magnitude of in-shock cooling effects via idealized numerical experiments and find that in-shock cooling should be minimal for the appropriate physical conditions and resolutions used throughout this paper.

In the full feedback simulations presented in this paper, BAL winds are implemented by imparting kicks of 5,000 km s to particles near the central black hole based on the accretion rate. The wind particles shock when they reach the ambient static ISM, thermalizing their kinetic energy, and giving rise to a physical situation similar to that shown in Figure 1 of Faucher-Giguere & Quataert (2012). To properly capture the full impact of BAL wind injection on quasar outflows, it is important that the thermalization of the BAL wind kinetic energy at the reverse shock does not suffer substantial from in-shock cooling. Gas densities, temperatures, and velocities for the reverse shock are set by the pre-shock BAL wind material, which is assumed to be free streaming. The non-homogenous density structure of the ISM and variability of the AGN radiation field make identifying the impact of in-shock cooling difficult in the full simulations directly. We instead construct idealized shock tube tests to recreate these conditions in a setting where resolution-dependent in-shock cooling can be directly identified.

We use a three-dimensional shock tube to explore the reverse shock density and temperature profile as a function of physical conditions (i.e. pre-shock density) and numerical resolution. The idealized initial conditions for the reverse shock include a fast moving medium (imitating the BAL wind material) moving into a static medium (imitating the ambient ISM). The BAL wind material is given an initial temperature of K, however this can change rapidly at the onset of the simulation due to Compton heating off of the AGN radiation if the gas density is sufficiently low. The BAL material is given a velocity of km s. The ambient ISM material is given an initial temperature of K and is initially static. The initial density for the BAL wind, initial density for the ambient ISM, and incident flux of AGN radiation are dependent on the location of the shock. We approximate the density of the pre-shock free-streaming BAL wind material to be given by

the density of the ambient ISM as

and the incident AGN radiation flux as

We assume fiducial values of , , and and run tests for pc. The shock tube uses periodic boundary conditions in a rectangular prism of dimension kpc.

Figure 12 shows the gas density and temperature profiles across the idealized shock at kyrs. The three panels show different values for the ambient gas density and AGN radiation flux (corresponding to pc, as described in the previous paragraph) tests with the legend indicating the gas particle mass resolution in each test. We find that the lowest resolution test (black line) blurs the location of the reverse shock substantially in the and pc tests. The two higher resolution tests present with less blurring the of the reverse shock, however there is still an offset present in the location of the reverse shock owing to the low particle number in the pre shock low density BAL wind material. In terms of in-shock cooling, we find that the post shock gas forms a stable and nearly flat temperature profile which shows little variation as we change the mass/particle resolution for our highest two resolution tests. Although some shock broadening is present, the post-shock gas temperatures are not strongly (if at all) impacted by the increasing resolution. If in-shock cooling were present, we would instead expect the post-shock gas temperature to decrease with lower mass resolution. Since the high resolution tests explored here () have resolution comparable to that used in the simulations presented in this paper and show little indication in-shock cooling, we conclude that in-shock cooling should not significantly impact the post-shock gas temperatures, and therefor should not have a significant impact on the results presented in this paper.

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description