Feedback from Active Galactic Nuclei: Energy- versus momentum-driving

Feedback from Active Galactic Nuclei: Energy- versus momentum-driving

Tiago Costa111E-mail:, Debora Sijacki and Martin G. Haehnelt
Institute of Astronomy and Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB HA

We employ hydrodynamical simulations using the moving-mesh code AREPO to investigate the role of energy and momentum input from Active Galactic Nuclei (AGN) in driving large-scale galactic outflows. We start by reproducing analytic solutions for both energy- and momentum-driven outflowing shells in simulations of a spherical isolated dark matter potential with gas in hydrostatic equilibrium and with no radiative cooling. We confirm that for this simplified setup, galactic outflows driven by a momentum input rate of order can establish an relation with slope and normalisation similar to that observed. We show that momentum input at a rate of is however insufficient to drive efficient outflows once cooling and gas inflows as predicted by cosmological simulations at resolved scales are taken into account. We argue that observed large-scale AGN-driven outflows are instead likely to be energy-driven and show that such outflows can reach momentum fluxes exceeding within the innermost of the galaxy. The outflows are highly anisotropic, with outflow rates and a velocity structure found to be inadequately described by spherical outflow models. We verify that the hot energy-driven outflowing gas is expected to be strongly affected by metal-line cooling, leading to significant amounts () of entrained cold gas.

methods: numerical - black hole physics - cosmology: theory

1 Introduction

There is now mounting evidence that Active Galactic Nucleus (AGN) activity can drive powerful large-scale outflows. Observational signatures of outflows in the form of broadened (Doppler shifted) emission and absorption lines have been detected for a number of galaxies hosting an AGN. Detected outflows have characteristic speeds , spatial scales of and often appear to consist of a complex multi-phase medium. This typically comprises a hot ionised component that can travel at speeds as high as (Humphrey et al. 2010; Greene et al. 2011; Nesvadba et al. 2011; Westmoquette et al. 2012; Rupke & Veilleux 2013b; Liu et al. 2013; Mullaney et al. 2013; Förster Schreiber et al. 2013; Rodríguez Zaurín et al. 2013; Villar Martín et al. 2014; Arribas et al. 2014; Harrison et al. 2014; Genzel et al. 2014), a partially overlapping neutral atomic component at speeds not much exceeding (Krug et al. 2010; Sturm et al. 2011; Rupke & Veilleux 2013b) and a substantial portion of cold molecular gas, as revealed by spatially resolved CO-emission and as OH-emission/absorption features (Feruglio et al. 2010; Sturm et al. 2011; Alatalo et al. 2011; Aalto et al. 2012; Cicone et al. 2012; Veilleux et al. 2013; Rupke & Veilleux 2013a; Combes et al. 2013; Cicone et al. 2014; Sakamoto et al. 2014). Molecular outflows detected in ultra-luminous infrared galaxies (ULIRGs) and quasars (QSOs) possess speeds of the order of and high mass outflow rates up to (see e.g. Sturm et al. 2011; Cicone et al. 2012). Large-scale outflows at scales of and speeds have also been detected in luminous high redshift QSOs (Cano-Díaz et al. 2012; Maiolino et al. 2012), where the estimated outflow rates () exceed the inferred star formation rates of the host galaxies. While it is difficult to disentangle the combined effects of both supernova and AGN feedback effects in driving outflows, theoretical models of supernova-driven feedback exclude outflow speeds much exceeding on energetic grounds (Martin 2005; Sharma & Nath 2013). Most crucially, several studies indicate that the presence of an AGN boosts the properties of detected outflows. Using CO observations of a sample of local ULIRGs and QSOs, Cicone et al. (2014) have shown that the mass outflow rates, the kinetic luminosity, the momentum flux and the spatial extension of detected outflows all correlate with the fraction of the host galaxy’s bolometric luminosity attributed to a central AGN. The depletion times of molecular gas of estimated for outflows in Cicone et al. (2014) are lower than the depletion time for star formation, providing tentative evidence that AGN-driven outflows can quench the host galaxy.

Observed AGN-driven outflows thus appear to strongly affect their host galaxies, potentially realising the negative AGN feedback effect invoked by galaxy formation models in order to address a series of open questions concerning the properties of massive galaxies. AGN feedback may account for how massive galaxies move away from the ‘main-sequence’ of star forming galaxies to become ‘red and dead’ (Di Matteo et al. 2005; Springel et al. 2005a; Sijacki & Springel 2006; Okamoto et al. 2008; Hopkins et al. 2008; Dubois et al. 2013a; Martizzi et al. 2014), explaining the shape of the otherwise overpopulated high-mass end of the galactic stellar mass function (Scannapieco & Oh 2004; Churazov et al. 2005; Kawata & Gibson 2005; Croton et al. 2006; Bower et al. 2006; Sijacki et al. 2007; Somerville et al. 2008; Booth & Schaye 2009; Puchwein & Springel 2013; Vogelsberger et al. 2014). Other potential roles of AGN feedback include helping to drive the morphological transformation between spirals and ellipticals, leading to metal enrichment of the intracluster (ICM) and intergalactic medium (IGM) and thereby shaping their thermodynamic properties (Sijacki et al. 2007; Puchwein et al. 2008; McCarthy et al. 2010; Gaspari et al. 2011; Teyssier et al. 2011; Martizzi et al. 2012; Planelles et al. 2014). It has also been argued that AGN-driven outflows could provide a ‘positive feedback’ effect whereby star formation occurs in compressed outflowing gas (Silk 2005; Zubovas et al. 2013; Silk 2013) or even that it may drive the apparent size evolution of elliptical galaxies (Ishibashi et al. 2013; Ishibashi & Fabian 2014). Finally, AGN feedback also offers a compelling mechanism for the origin of the correlations observed between the mass of the supermassive black hole and the mass and velocity dispersion of the bulge of the host galaxy (e.g. Magorrian et al. 1998; Tremaine et al. 2002; Marconi & Hunt 2003; Ferrarese & Ford 2005; Gültekin et al. 2009; McConnell & Ma 2013; Kormendy & Ho 2013). Balancing the energy or momentum released by the AGN with the gravitational weight of its surrounding interstellar material has been argued to reproduce the scaling and normalisation of observed correlations (Silk & Rees 1998; Haehnelt et al. 1998; Fabian 1999; King 2003; Murray et al. 2005).

The efficiency of AGN feedback is ultimately dependent on the detailed hydrodynamics of AGN-driven outflows, which in turn is sensitive to a wide range of physical effects including the efficiency of (in-shock) cooling and star formation, the geometry, dynamical and thermodynamic state of the gaseous medium through which the outflow propagates. In Section 2, we review the current theoretical understanding of AGN-driven outflows as envisaged in models in which the energy and/or momentum deposited by the AGN on galactic scales couples hydrodynamically (as opposed to radiatively) with the interstellar medium (ISM) via a fast inner AGN wind. We clarify the nature of energy- and momentum-driven shells of outflowing material in spherically symmetric galactic potentials and discuss important distinctions between these two limiting cases. In Section 3, we discuss the implementation of energy- and momentum-driven outflow feedback into numerical simulations. In Section 4, using simulations of isolated dark matter potentials with gas in hydrostatic equilibrium and no radiative cooling, we reproduce analytical solutions for energy- and momentum-driven shells. Using simulations including radiative cooling, we show that under the right circumstances, large amounts of outflowing molecular gas are likely to form due to the cooling of shock-heated outflowing material. Finally, in Section 5, we implement the same AGN energy- and momentum-driven outflow models in fully cosmological hydrodynamic simulations and present several important departures from the idealised picture envisaged in analytic models. We place particular emphasis on the origin of an relation in the context of such simulations. We summarise our conclusions in Section 6. Tests of our numerical implementations of energy- and momentum-driven outflows against changes in numerical resolution and the energy/momentum injection procedure are discussed in Appendix A. Possible future research avenues in the context of AGN feedback are discussed in Appendix B.

2 Energy and Momentum-driven outflows

2.1 Accretion as source of energy

Energy and momentum deposited by AGN into the ISM of galaxies as well as into their surrounding IGM/ICM is initially released by accretion of baryonic matter on to supermassive black holes (Lynden-Bell 1969). A significant fraction of the accreted rest mass energy is radiated away close to the event horizon of the black hole, leading to an accretion luminosity given by


where is the black hole accretion rate, is the speed of light in vacuum and is the radiative efficiency of the black hole, which typically lies in the range depending on the spin of the black hole. Studies comparing observed AGN luminosities and the inferred black hole mass density suggest that, on average, (Soltan 1982; Fabian & Iwasawa 1999; Yu & Tremaine 2002), a value consistent with moderately spinning black holes (a recent study by Ueda et al. (2014) however gives ).

A point of much debate is the fraction of the luminous energy (Eq. 1) that can couple to surrounding interstellar gas and drive an outflow. Observed AGN show signatures of outflows down to the smallest resolvable scales. These inner outflows take the form of highly collimated jets directly observable by their radio emission as well as wide-angle accretion disc winds inferred from very broad absorption lines in so called ‘broad absorption line’ (BAL) QSOs (e.g. Pounds et al. 2003; Ganguly et al. 2007; Reeves et al. 2009; Tombesi et al. 2010). Inferred outflow rates are typically similar or higher than the black hole accretion rate. The luminosity of AGN is thereby often assumed to be limited to the Eddington luminosity, which is given by


where represents the opacity, the mass of the accreting black hole and the gravitational constant.

Here, we focus on the direct hydrodynamical coupling of an outflow launched from the inner parts of the accretion disc to interstellar gas at scales. In Section , we briefly discuss the alternative suggestion of a coupling of radiation escaping from the AGN to dusty interstellar material at galactic () scales. We next review simple analytical models addressing the interaction between the inner AGN wind and the ISM without cooling.

2.2 Spherical wind models

We start by investigating the simplified model proposed by King (2003, 2005) describing the interaction of a sub-relativistic accretion disc wind, from now on referred to as ‘inner wind’ for brevity, with the ISM surrounding the accreting black hole, which is assumed to lie at the centre of the galactic halo. In line with observed nuclear (ultra-)fast outflows (Pounds et al. 2003; Pounds & Reeves 2009; Reeves et al. 2009; Tombesi et al. 2010, 2011, 2012, 2013; Pounds & King 2013), King (2003, 2005) assumes the inner wind to have simple properties. It is taken to have a high covering fraction with opening angle such that and its mass outflow rate is taken to be equal to the black hole’s Eddington accretion rate . When combined, these two properties contribute to a high opacity in the inner wind and an optical depth to electron scattering of , supporting the view that such an inner wind can indeed be driven via radiation pressure from the AGN (King & Pounds 2003). Consequently, the inner wind initially transports momentum at a rate comparable with the momentum injection rate of the AGN as given by


where is the speed of the inner wind and is the mass outflow rate of the inner wind, which is, as already mentioned, assumed to be equal to the Eddington accretion rate. From Eq. 1, it thus follows that


Substitution of as given in Eq. 4 into Eq. 3 then fixes the velocity of the inner wind to


Mass continuity implies that the mass density of the inner wind must in turn be given by


Finally, the inner wind’s kinetic luminosity is given by


If for the radiative efficiency, the canonical is adopted, the inner wind speed is in agreement with observed blueshifted absorption lines in QSO spectra (see e.g. Pounds & King 2013) and the kinetic luminosity is with . For fixed momentum flux, note that Eqs. 3 and 7 however imply that a lower results in higher and lower , respectively.

With a speed of , the inner wind is highly supersonic and must drive a strong shock into the ISM. The structure of the shock pattern that ensues is analogous to that resulting from the interaction between a fast stellar wind and the surrounding interstellar gas (see e.g. Weaver et al. 1977; Dyson & Williams 1997). The inner wind drives a forward shock wave that thrusts into the ambient ISM, while the inner wind itself must decelerate strongly in an inner reverse shock facing towards the black hole. Thus, in order of increasing distance from the AGN, the flow pattern consists of four zones (see Fig. 1): (a) the unshocked highly supersonic inner wind; (b) the shocked inner wind material that has crossed the reverse shock, also often termed ‘wind shock’ in the literature (e.g. Zubovas & King 2012; Faucher-Giguère & Quataert 2012); (c) a shell of interstellar gas swept up by the forward shock and (d) the unperturbed ambient ISM.

Figure 1: The shock pattern resulting from the collision of a fast inner AGN wind with the surrounding ISM (see also Weaver et al. 1977; Dyson & Williams 1997; Zubovas & King 2012; Faucher-Giguère & Quataert 2012). The flow is divided into four distinct regions, which are labelled with letters in order of increasing distance to the AGN. Region (a) is filled with the unshocked inner wind which flows out with a speed . Since the inner wind acts as a piston pushing into the ISM at highly supersonic speeds, a shock must naturally form. However, in the reference frame of the inner wind itself, it is the ISM that moves against it highly supersonically, giving rise to a reverse shock (region (b)). Neglecting relativistic effects, gas behind the reverse shock has a temperature (see text). On the outside, region (c) is bounded by a forward shock that sweeps up interstellar gas into an expanding shell. Regions (b) and (c) are separated by a contact discontinuity (dashed circle). Finally, region (d) is occupied by the undisturbed ISM. It is the cooling of region (b) that determines whether the outflowing shell in region (c) is energy- or momentum-driven.

If the forward shock is strong, the pressure from the ambient ISM (region d) can be neglected and the dynamics of the shell of swept up gas (region c) is fully determined by the difference between the outward pressure force exerted on it by the shocked inner wind (region b) and the inward pull due to gravity on the shell. The latter depends on the total mass (dark matter and black hole)***Note that the contribution of a stellar component is not explicitly treated in King (2003). enclosed by the radius of the shell and on the mass of the swept up shell itself , where is the gas mass fraction of the halo. Balancing the two forces gives the equation of motion of the outflowing shell as

Figure 2: Momentum and energy-driven AGN outflows. A momentum-driven outflow (shown on the left) occurs when the inner wind that shock-heats as it passes the reverse shock (region (b)) cools efficiently. The resulting (isothermal) shock consists of an initially adiabatic shock followed by a region where radiative cooling dominates shown in light blue. The loss of thermal energy results in a drop of pressure support and causes region (b) to be very thin. The post-shock pressure at the contact discontinuity in this case is . An energy-driven outflow (shown on the right) results when the energy injected by the inner wind (region (a)) is fully conserved throughout the outflow, i.e. no substantial cooling losses are sustained. Region (c) is driven out as the now hot and thick region (b) expands adiabatically. The separate question of how efficiently region (c) cools does not determine whether the outflow is energy- or momentum-driven. We revisit this question in Section 4.2 and for now represent region (c) as thick, while mindful it can collapse into a very thin shell if cooling is efficient.

The properties of the fluid right across the reverse shock at , where the subscript denotes ‘reverse shock’, can be obtained from the Rankine-Hugoniot jump conditions for a strong adiabatic shock. The number density of the shocked inner wind across the reverse shock is accordingly given by


using the result from Eq. 6 to substitute for .

Assuming the reverse shock to have a speed , the temperature just behind the shock is


Note that such high temperatures are only possible to achieve due to the high speed of of the inner (unshocked) AGN wind and by neglecting relativistic effects in the hot plasma. The thermal energy in the shocked inner wind can be either transferred to the ISM or lost via radiative cooling. The dynamics of the resulting (out)flow pattern depends critically on which of these is the fate of the thermal energy contained in the hot shocked inner wind.

2.3 Momentum-driven limit

If the reverse shock cools efficiently, shocked inner wind gas in region (b) traverses a radiative region with a size of the order of the cooling length until its temperature returns to its pre-shock value at a radius . This sequence of an initially adiabatic shock and a cooling region constitutes an isothermal shock. Since cooling is efficient, , it follows that the shock is approximately planar, i.e. region (b) is very thin (see Fig. 2). For a strong isothermal shock, the Rankine-Hugoniot jump conditions yield a post-shock pressure given by


The post-shock gas then drives an outer shock into the unperturbed ISM, thrusting with a pressure equal to the ram-pressure of the pre-shock inner wind, i.e. . This result becomes physically intuitive when recalling that region (b) is very thin. The inner wind can then be imagined to collide with the ISM directly, transferring its momentum fully in the interaction. In this limit, the outflow is said to be momentum-driven since the driving force is equal to the momentum flux of the inner wind. Eq. 8 can then be simplified to


For an isothermal halo with velocity dispersion (see King 2005; McQuillin & McLaughlin 2012; Ishibashi & Fabian 2012), this can be recast into the form


where . Neglecting the last term on the right hand side of Eq. 13 (valid if the launching speed of the shell is higher than the escape speed from the black hole’s gravitational potential), the differential equation admits unbound solutions only if . If AGN outflows are momentum-driven, they can in principle escape to arbitrarily high radii only when the central black hole grows beyond a critical mass, at which the outflowing shell becomes unbound. At this critical value, further accretion on to the black hole is expected to become inefficient and the black hole mass is limited to


where the superscript ‘’ denotes ‘momentum-driven’. Remarkably, the normalisation and the scaling are in close agreement with the observed relation (Ferrarese & Merritt 2000; Gebhardt et al. 2000; Tremaine et al. 2002; Gültekin et al. 2009; McConnell & Ma 2013; Kormendy & Ho 2013).

2.4 Energy-driven limit

If the reverse shock is instead unable to radiate away its thermal energy, the shell is driven by the adiabatic expansion of the hot shocked wind bubble (see Fig. 2). In the limit where the full energy of the inner wind is conserved, the resulting outflow pattern is termed energy-driven.

In an energy-driven outflow, the expansion rate of the shell of shocked interstellar gas is equal to the rate at which the shocked wind does ‘PdV’ work on its surroundings. By balancing energy losses due to ‘PdV’ work and work done for overcoming gravity with energy injection by the inner wind, we can write an energy equation that reads


where the first equality simply follows from the ideal gas law and is the adiabatic index of the gas. As mentioned, if . Using Eq. 8 to eliminate from Eq. 15, assuming an isothermal halo such as in King (2005), King et al. (2011) and McQuillin & McLaughlin (2013) and that the shell’s initial speed is sufficient for it to escape from the black hole potential yields the equation of motion for the shell


Eq. 16 also admits unbound solutions. Assuming a solution of the type (King et al. 2011) and taking into Eq. 16 gives


When this condition is satisfied, the energy-driven outflow can clear the halo off its gas and terminate black hole accretion. In this case, the black hole mass is fixed to a value given by


using ‘’ to denote ‘energy-driven’ (see also Silk & Rees 1998; Haehnelt et al. 1998; King 2010b; Fabian 2012; McQuillin & McLaughlin 2013). Energy-driven outflows therefore yield , somewhat steeper than the observed scaling (see e.g. Kormendy & Ho 2013).

2.5 Energy- versus momentum-driving and the role of Compton cooling/heating

2.5.1 The efficiency of energy- vs. momentum-driving

Using Eqs. 14 and  18, the ratio of the final black hole masses for energy- and momentum-driven outflows is given as


suggesting that lower black holes masses are required to drive a large-scale outflow for energy- than for momentum-driven solutions. In this section, we clarify the reason for the relative efficiency of energy- over momentum-driven outflows.

Whether an outflow is energy- or momentum-driven, its kinetic luminosity can be parametrised in terms of a kinetic conversion efficiency that quantifies the fraction of the AGN luminosity which is converted into kinetic energy of the outflow on galactic scales. The outflow’s kinetic luminosity is thus given by


The momentum flux of the outflow can, in turn, be written as


again, irrespective of whether it is energy- or momentum-driven.

If the outflow is energy-driven, then is related to via


where is the efficiency in converting the thermal energy injected by the AGN into kinetic energy of the outflow on galactic scales.

We now evaluate the ratio of the momentum flux of the outflow to that of the inner wind, a ratio often termed ‘momentum boost-factor’. Substituting Eq. 22 into Eq. 21 and using (Eq. 3) gives


Energy-driven shells can therefore acquire momentum fluxes substantially higher than the original momentum flux of the inner wind (see also Zubovas & King 2012; Faucher-Giguère & Quataert 2012).

Conversely, if the outflow is momentum-driven, momentum conservation across the isothermal shock implies that


and hence, by combining Eq. 21 and Eq. 24


As expected, the boost factor in the momentum-driven regime is given by


in contrast with the energy-driven case, where high boost factors can apply.

The reason for the higher relative efficiency of energy-driven outflows stems from the inequality , which applies for typical values for (). Energy-driven shells propagate more rapidly than momentum-driven outflows and their higher momentum flux means they can more vigorously remove gas from the halo.

2.5.2 The role of Bremsstrahlung and inverse Compton cooling/heating

In practice, if AGN-driven outflows arise due to an inner wind shock, they are likely to be a mixture of the energy- and momentum-driven limits discussed above. Whether, however, the outflows are more suitably approximated by their energy- or momentum-driven limit depends on the relative locations of the reverse shock radius and the cooling radius , since it is the cooling of the shocked wind that discriminates between the two cases. If , the hot inner wind shock cannot cool and the outflow falls towards its energy-driven limit. Conversely, if , the reverse shock cools efficiently and the outflow is more accurately represented by its momentum-driven limit.

For gas of temperatures in the range (see Eq. 10), the only significant cooling processes are free-free emission and inverse Compton scattering. We begin by estimating the timescale for free-free (thermal Bremsstrahlung) emission (see e.g. Mo et al. 2010) for the shocked inner wind gas


where we have taken the electron number density as (for a fully ionised primordial gas) and assumed the fiducial values given in Eq. Implicit in this calculation is the approximation that the density and temperature of the inner wind shock bubble is uniform (see Faucher-Giguère & Quataert 2012)..

The outflow time of the shocked inner wind is given by


Setting gives a cooling radius of , i.e. the shocked inner wind cannot cool via free-free emission.

A much more efficient cooling mechanism arises from the energy exchange that takes place as photons in the radiation field of the AGN and electrons in the hot shocked inner wind gas undergo Compton scattering. In this process, high energy photons can transfer energy to (non-relativistic) electrons if , thereby heating the gas, while low energy photons with gain energy, effectively cooling the gas. The temperature at which Compton heating and cooling balance out, the ‘Compton temperature’ , depends on the spectral energy density of the emitting quasar and will therefore vary from source to source and with accretion rate (see Gan et al. 2014, for a recent discussion on the AGN feedback effects of a varying Compton temperature). For the average QSO, Sazonov et al. (2004) estimate . Since (as given in Eq. 10), the hot shocked wind cools (see Bourne & Nayakshin 2013, for the potential observational imprint of inverse Compton cooling of the inner wind shocked bubble). Using the estimated cooling function for Compton heating/cooling of a gas exposed to a QSO continuum with provided in Sazonov et al. (2005) gives a cooling time


In the ultra-relativistic regime (for ), note that the inverse Compton cooling rate however scales as , while the gas thermal energy still scales as , resulting in even lower cooling times than suggested by Eq.  (which only holds for non-relativistic electrons) (see e.g. Pozdnyakov et al. 1983). For mildly relativistic electrons, as is the case for , the cooling time given in Eq. is however likely overestimated by a just a factor of a few. Setting gives the cooling radius for inverse Compton cooling


In the simulations presented in this study, we do not explicitly model Compton cooling/heating because neither the inner wind nor the reverse shock are resolved. Efficient Compton cooling is implicit in our momentum-driven model however. Note also that once gas cools below , Compton scattering can instead heat gas and can provide an important source of feedback in its own right as shown by Ciotti & Ostriker (1997, 2001) and Sazonov et al. (2005).

According to the estimate of Eq. 30, the shocked inner wind cools very efficiently below a radius . In this regime, the outflow will be momentum-driven (King 2003). Once the momentum-driven shell moves beyond , the outflow becomes energy-driven due to the relative inefficiency of inverse Compton cooling in the shocked inner wind. Note, however that if the electrons and protons in the hot shocked wind plasma are thermally decoupled, as is known to occur in very strong shocks, (e.g. in supernova remnants, Ghavamian et al. 2007), the cooling radius will be reduced to too close to the sphere of influence of the black hole for meaningful cooling to occur (Faucher-Giguère & Quataert 2012).

2.6 Radiation pressure-driven outflows

A separate class of models for AGN feedback envisages that AGN-driven outflows are driven by direct radiation pressure on dust grains embedded in the ISM of the galaxy (Fabian 1999; Murray et al. 2005; Debuhr et al. 2011). In this scenario, interstellar gas is pressurised by the radiation field of the AGN rather than by a shocked AGN inner wind as in King (2003, 2005). Observational support that has been invoked in favour of this picture is for instance the observed lack of AGN with simultaneously high luminosities and high hydrogen column densities (Raimundo et al. 2010; Fabian 2012), where radiation pressure would be most efficient at driving outflows. It has been argued that radiation pressure can lead to large-scale outflows and to momentum-boosts, as IR reprocessed photons are repeatedly absorbed (see e.g. Roth et al. 2012). Since their dynamics can be described by an equation essentially identical to Eq. 12 in the optically thick regime, radiation pressure-driven outflows have also been termed ‘momentum-driven’ in the literature. These outflows are however not to be confused with the momentum-driven outflows discussed in this study, which have a fundamentally different physical origin. We postpone a more detailed assessment of the plausibility of IR radiation-pressure wind models and their implementation in numerical simulations to future work. Note, however, that such radiation pressure driving requires large spatially coherent IR optical depth which is increasingly difficult to realise on galactic ( kpc) scales. See also Krumholz & Thompson (2013) and Davis et al. (2014) for recent discussions on the possibility of launching IR radiation pressure-driven winds.

3 Numerical implementation of the galactic wind models

3.1 Preliminaries

Models of AGN-driven outflows (see Section 2) involve a large disparity of physical scales; ranging from where energy is released from the accretion disc that feeds the supermassive black hole, all the way to galactic () scales at which the large-scale outflow sweeps across the galaxy and the surrounding halo. The ab initio treatment of the interaction between an AGN and its host galaxy in numerical simulations requires currently infeasible resolution and dynamical range. In the context of the inner wind model of AGN feedback (Section 2), lack of resolution means that it is currently neither possible to resolve the inner wind nor the reverse shock in numerical simulations. Numerical implementations of AGN feedback therefore always rely on ‘subgrid’ models that attempt to mimic the effects of unresolved outflow physics at resolved scales (see e.g. Wurster & Thacker 2013, for a comparison of several widely used AGN feedback models). The relation of popularly adopted ‘subgrid’ models with simple hydrodynamical models such as summarised in Section 2.2 however often remains unclear.

Our strategy in this study is to mimic the presence of a reverse shock and its behaviour in both the momentum- and energy-driven (limiting) regimes using ‘subgrid’ models (see also Nayakshin & Power 2010; Power et al. 2011; Nayakshin & Zubovas 2012). Uncertainties in the properties of the inner wind are bypassed in numerical simulations by assuming that either energy or momentum (scaled to the luminosity of the QSO) is injected at galactic scales (). Injection of energy is normally straightforward to implement and has the advantage that a fixed amount of energy is added to the wind. It realises efficient feedback as long as cooling is not overwhelming, e.g. at large distances from the dense galactic centre or at high energy input rates. It probably, however, underestimates the effect of AGN feedback when cooling times are short compared to outflow times and too little of the thermal energy is converted into kinetic energy. Injection of momentum, on the other hand, is usually motivated by arguing that cooling is efficient in the inner outflow, where the optical depth is expected to be high such that a momentum flux should be injected at large scales. This approach has the advantage that galactic winds have greater impact in rather dense environments and therefore smaller radii where cooling is important. It has, however, the disadvantage that the corresponding injected energy scales linearly with the outflow velocity of the gas into which the momentum is injected and thus depends strongly on the physical and dynamical state of the ISM. Like energy injection, it is also very sensitive to the resolved distance to the AGN. It is important to note that momentum or thermal energy injection in numerical simulations do not necessarily respectively translate to what is meant by energy- or momentum-driven outflows in the classical sense. Particles/cells receiving too large a momentum/velocity kick will shock and thermalise rapidly, driving an outflow via adiabatic expansion rather than by ram-pressure. Such an outflow rather resembles the energy-driven case and cannot be considered ‘momentum-driven’ even if it relies on momentum injection. Conversely, if thermal energy is added to gas which cools very efficiently so that most of the injected energy is radiated away, gas may be driven out due to residual bulk motion and ram-pressure rather than via adiabatic expansion. Such an outflow can no longer be described as ‘energy-driven’. Such uncertainties imply that the physical interpretation of any ‘subgrid’ model is dependent on the included physics (e.g. cooling) and the spatial resolution of the simulation. Bearing these difficulties in mind, we now describe implementations of energy- and momentum-driven outflows that successfully reproduce the predictions of the analytical models of King (2003, 2005) in their idealised setting.

3.2 Numerical setup

For our numerical simulations, we employ the moving-mesh code AREPO (Springel 2010), which uses a second-order accurate finite-volume method on an unstructured Voronoi mesh that is allowed to move with the flow in order to model gas hydrodynamics. For our simulations of an isolated halo, we adopt a static Hernquist dark matter potential (Hernquist 1990; Sijacki et al. 2012). Even though an isothermal profile is typically assumed in simple analytical models (Silk & Rees 1998; Fabian 1999; King 2003, 2005), the relevant equations of motion can be straightforwardly adapted for a Hernquist halo by taking the expression for the enclosed mass in Eq. 12 and Eq. 15 as


where is the total mass and is the scale length of the halo. The halo’s total mass density is given by

Figure 3: Projection of gas mass along a slab of thickness and length for the energy-driven outflow model (top row) and for the momentum-driven outflow model (bottom row) at different simulation times for a black hole of mass emitting at its Eddington luminosity. Injection of energy or momentum at the centre of the halo drives a shock that compresses almost all intervening gas into a propagating shell. For matching black hole mass, the energy-driven outflows propagate more rapidly than the momentum-driven outflows due to their higher efficiency in converting AGN luminosity to kinetic luminosity.

The Hernquist potential has the advantage that its total mass converges as (removing the need for a truncation at an arbitrary radius) while closely resembling the Navarro-Frenk-White profile (Navarro et al. 1997). Moreover, solutions to Eqs. 12 and 15 converge to a unique solution at large radii for a Hernquist halo (McQuillin & McLaughlin 2012), facilitating a comparison with analytical solutions. For our numerical experiments, we select a halo of total mass , which we populate with resolution elements, resulting in an initial gravitational softening length of . Note that in AREPO, gas has an adaptive gravitational softening set by the cell size. The minimum cell size in our simulations reaches . As in King (2003, 2005), we assume that a fraction of the total mass is in gas such that the gas density is . We choose in all our tests. Finally, we choose the halo to have a concentration of as expected for a halo at (e.g Macciò et al. 2008). This choice of parameters fixes the virial radius defined as the radius enclosing a mean density times that of the cosmic mean. to and the scale length to .

Initial conditions were generated by randomly sampling the halo’s density structure out to a very large radius of in order to ensure its density profile to be as close a match as possible to its analytical expression (Eq. 32). Starting from these initial conditions, the halo was then evolved non-radiatively for a period of to minimise halo relaxation effects. A black hole sink particle was then introduced at the centre of the halo and its position was fixed at this location for the entire duration of the various simulations. For the black hole mass, we explored values ranging from to , in order to probe solutions for black hole masses below, equal and above the critical value for the halo§§§For any halo with a peaked circular velocity () profile, it can be shown that (McQuillin & McLaughlin 2012) respectively. In order to establish a close comparison with the models of King (2003, 2005), it is assumed that the AGN constantly emits at its Eddington luminosity for all simulations. In Section 4.2, we relax this assumption.

3.3 Energy-driven limit

In the energy-driven regime, the hot shocked inner wind expands as it does ‘PdV’ work on its surroundings (see Section 2.4). Its thermal energy changes according to Eq. 15, where it was assumed that the kinetic energy of the inner wind is fully converted into thermal energy across the reverse shock. We reproduce this behaviour in our simulations by letting a fraction of the energy emitted by the AGN couple thermally with its neighbouring gas cells. The algorithm first identifies the nearest gas cell neighbours, so defined if they are located within an SPH-kernel of the black hole particle. Energy is then injected at a rate into the neighbouring cells in a kernel-weighted fashion, so that gas cells closer to the black hole receive a higher amount of thermal energy. For the simulations presented in Section 4.1, note that always. In Section 4.2, we explore varying the AGN lightcurve. In contrast with typical implementations of this widely used AGN model, we do not prevent gas cells from heating above a maximum temperature. This procedure is usually adopted in order to prevent gas particles/cells from heating to temperatures at which relativistic effects become relevant and for which our treatment of hydrodynamics breaks down. In the context of this study however, imposing a temperature ceiling would effectively limit the energy injected by the AGN and would introduce an obvious source of disagreement with analytical models. Our implementation is otherwise very close to previous implementations of AGN thermal feedback (Springel et al. 2005b; Di Matteo et al. 2005; Sijacki et al. 2007; Di Matteo et al. 2008; Sijacki et al. 2009; Booth & Schaye 2009; Dubois et al. 2013a, b). In Appendix A, we show that in the context of this study, this model is largely insensitive to the selected number of gas cell neighbours.

3.4 Momentum-driven limit

In order to reproduce the momentum-driven limit of AGN outflows, it is necessary to ensure that gas is accelerated by the impinging ram-pressure of the (unresolved) inner wind rather than ‘PdV’ work of the shocked inner wind gas, which now has negligible thermal pressure. Instead of injecting thermal energy, we impart momentum kicks to the nearest neighbour gas cells directly at rate in the same spirit as previous implementations of mechanical feedback (Ostriker et al. 2010; Nayakshin & Power 2010; Debuhr et al. 2011; Choi et al. 2012, 2013, 2014). Note that the magnitude of the effective velocity kick is kernel-weighted, in order to ensure that gas closest to the black hole is expelled most efficiently. We have also run a set of test simulations, in which we varied the number of gas cells into which momentum is injected (see Appendix A). For a large number of neighbours, we have found that the numerical and analytical solutions converge. The numerical results overshoot with respect to the analytical solution if the selected number of neighbours is too small. In this case, gas cells effectively receive larger velocity kicks. In this regime, agreement between numerical and analytical solutions can only be restored by increasing the time-stepping accuracy.

4 Comparison of hydrodynamic simulations with the analytical wind model

4.1 Simulations without radiative cooling

Figure 4: Time sequence of gas density profiles in the numerical models for an energy-driven (left hand plot) and momentum-driven (right hand plot) outflow for a black hole mass of emitting at its Eddington luminosity. Note that approximately matching radial distances are traversed in much less time by the energy-driven outflow than the momentum-driven outflow. The higher velocity of energy-driven outflows is a consequence of the higher conversion factor of AGN luminosity to kinetic luminosity of the outflowing shell. Note also that for the momentum-driven outflow, the speed of the outflowing shell quickly decelerates into the subsonic regime, generating a forward ripple travelling at the speed of sound.

A time sequence of the blast waves generated in the numerical models for energy and momentum-driven outflows detailed above is illustrated in Fig. 3, where the gas mass was projected along a thin slab of thickness . As in the theoretical picture proposed in Silk & Rees (1998), Fabian (1999) and King (2003), a thin shell propagates outwards, carrying with it almost all intervening gas. For matching black hole mass, the energy-driven outflow propagates more rapidly across the halo than a momentum-driven outflow in agreement with the considerations made in Section 2.5. Another feature that becomes clear in Fig. 3 is the development of Rayleigh-Taylor instabilities in the energy-driven outflow as hot gas rises buoyantly through under-densities in the shocked ISM component (see also King 2010a; Zubovas & King 2014). These instabilities ultimately lead to a pronounced departure from spherical symmetry in the outflowing shell.

The propagation of the shell is represented quantitatively in Fig. 4, where a time sequence of gas density profiles is shown for both energy- and momentum-driven models. At large radii, where gas has not yet been perturbed by the outflow, the profiles follow the Hernquist profile very closely. Strong deviations to the Hernquist profile occur at intermediate radii, where the gas density rises above the levels corresponding to a Hernquist potential at the same radius. This region is occupied by the shocked medium, swept up into a shell, and stretches all the way to the lowest radii, where the gas density then drops. In the momentum-driven outflow, note that the density profile has a very sharp peak in the innermost radii of the shell. This region is occupied by gas which is highly compressed after receiving a radial velocity boost and mimics the thin inner wind (isothermal) shock present in the momentum-driven outflow (see Fig. 2). For a black hole mass of , the thickness of the shell increases and becomes comparable with its distance to the black hole at . As shown later, this approximately matches the location at which the shell starts propagating subsonicallyNote that we have explicitly verified in our simulations that the shell always remains thin for solutions in which it always moves supersonically. (see Fig. 6). In this regime, the forward shock weakens and eventually propagates as a sound wave. At for instance, the forward disturbance satisfies , clearly indicating that it can no longer be described as a strong shock. Due to the shape of the Hernquist profile, note also that the density of the shell of shocked matter drops with distance from the black hole for both energy- and momentum-driven outflows. This is important when we later consider the cooling of the shocked medium, which, due to its dependence on , is particularly efficient in the central regions.

In order to establish a close comparison of the predictions of our numerical models with those presented by King (2003, 2005), we solved Eq. 12 and Eq. 15 numerically for and a Hernquist halo as described in Section 3.2. We then tracked the location of the contact discontinuity as a function of time in our simulated outflows and compared it with predictions based on the equations described in Section 2.2. The results of this comparison are shown in Fig. 5, where we plot the estimated location of the contact discontinuity in our various simulations as a function of time against the analytical predictions (solid lines) for a matching Hernquist halo and identical black hole masses. According to Fig. 5, the propagation of both energy- and momentum-driven outflows is in close agreement with analytical predictions throughout most of the simulation time. In energy-driven outflows, our simulations systematically underestimate the location of the contact discontinuity at later times when compared to analytical models. This occurs as a consequence of the formation of Rayleigh-Taylor instabilities at the interface between the shocked inner wind and the shocked medium components (see also King 2010a). Since the swept up material is no longer confined to a spherically symmetric shell, when tracking the contact discontinuity, our method preferentially picks up the location of overdense clumps that lag behind in the Rayleigh-Taylor unstable flow.

Figure 5: Radius of the expanding shell as a function of time for the energy- and momentum-driven models for three different black hole masses (filled circles). Solid lines represent the corresponding analytical solutions. The close match between analytical and numerical results indicate that the chosen sub-grid models provide an adequate approximation to the energy- and momentum-driven solutions. For the momentum-driven solution, the agreement is poorer for lower mass black holes. The source of disagreement is the back-reaction of the ambient medium on the outflow as it enters the subsonic regime, which is accounted for in our hydrodynamical simulations but neglected in the equations presented in Section 2, which assumes the forward shock is strong (see text for more details).

For the momentum-driven outflows shown on the right hand side in Fig. 5, agreement between numerical simulations and analytic solutions is less successful for the lower black hole masses, though close for . In order to establish the source of this disagreement, we show the outflow velocity as a function of radius according to the various analytical solutions in Fig. 6. We also show the sound speed profile in the simulated Hernquist halo as a dashed line. For the momentum-driven solutions with and , the speed of the outflowing shell becomes comparable to the speed of sound at and , respectively. The assumption of a strong forward shock is however implicit in Eq. 8, which should otherwise not neglect the confining pressure of the ambient gaseous mediumNote that Faucher-Giguère & Quataert (2012) include it in their analytic treatment. (see e.g. Weaver et al. 1977; Dyson & Williams 1997). As shown in Fig. 6, the assumption of a strong forward shock is valid for all energy-driven solutions for the black hole masses considered but is unsuitable in the momentum-driven limit for black hole masses . In this regime, the shock weakens (see Fig. 4) and the outflow becomes confined by the pressure of the ambient medium, causing it to decelerate. Note, however, that this finding does not alter the core of the conclusions of King (2003, 2005) qualitatively. Our results show that the shell of swept-up material still stalls when the black hole mass is lower than a critical value.

Figure 6: Analytical energy- (red lines) and momentum-driven (blue lines) solutions to Eq. 8. The local speed of sound is shown as a dashed line separating the sub- and supersonic regimes of the solutions. Momentum-driven solutions for black hole masses partially exist in the subsonic regime, where Eq. 8 is invalid. The behaviour of our simulated solutions deviates from the analytical predictions in this regime quantitatively, but do not change the key conclusion that there is a critical black hole mass, below which momentum-driven outflows stall and above which they can freely expand.

The structure of the shock in our energy-driven and momentum-driven models is illustrated in Fig. 7, which shows two-dimensional histograms of gas density against radial distance from the black hole. The colour coding reflects the local gas mass. Fig. 7 shows that the flow patterns can be clearly separated into different phases. For the energy-driven outflow, lower radii are sparsely populated by a diffuse hot component, which is generated in a ‘subgrid’ fashion in our simulations by directly injecting energy into gas cells surrounding the black hole. The adiabatic expansion of this component drives the forward shock containing the bulk of the outflow. Note that it is because the shell of dense gas lies on top of the very low density hot bubble that Rayleigh-Taylor instabilities arise in the outflow (see Fig. 3). The shock structure looks very different for the momentum-driven outflow shown on the right hand side of Fig. 7, where the mildly shocked gaseous medium is separated from the inner void by a thin region of highly compressed gas at which momentum is being injected. Because the shell lies on top of an even denser layer of gas (the ‘subgrid’ realisation of the cooled shocked inner wind), no Rayleigh-Taylor instabilities arise in this case******We have checked that this is true for the entire duration of all our simulations for momentum-driven solutions.. This finding is in good agreement with the analytical arguments of King (2010a), which suggest that energy-driven shells should always be Rayleigh-Taylor unstable, whereas momentum-driven shells for black hole masses close or much higher than should be Rayleigh-Taylor stable.

Temperature profiles for energy- and momentum-driven outflows are shown in Fig. 8. For the energy-driven outflow, three main regions can again be identified. In order of increasing radial distance from the AGN, these are the shocked wind phase at temperatures of , the shocked medium phase at temperatures of (in agreement with analytical estimates for an energy-driven outflow of speed as in Zubovas & King 2014) and the roughly isothermal ambient medium at the virial temperature . For the momentum-driven outflow, the temperature of the mildly shocked gaseous medium increases only slightly and drops towards the centre of the halo as it becomes increasingly poorly populated.

Figure 7: The structure of the outflow shown as a two dimensional histogram of the gas density as a function of radius for the energy-driven outflow model at (left) and for the momentum-driven outflow at (right) for a black hole mass of . The outflow can be separated into three distinct zones: (a) a shocked wind phase which is modelled in a ‘subgrid’ fashion by adding energy or momentum for energy- and momentum-driven outflows respectively, (b) a shocked medium phase, which is thick in the energy-driven outflow and (c) the ambient gas phase. In the innermost regions of the momentum-driven outflow, the shocked wind phase consists of a thin and highly compressed layer instead of a hot bubble, as in the energy-driven case. The arrow indicates the location of the forward (weak) shock.
Figure 8: Temperature profile of gas in the energy-driven outflow at (left) and in the momentum-driven outflow at (right) for a black hole with mass . The region directly affected by the ‘subgrid’ model is shaded in orange. For energy-driven outflows, our ‘subgrid’ model generates a very hot bubble that drives an outflow as it expands adiabatically, while for momentum-driven outflows, our ‘subgrid’ model evacuates the central regions entirely by imparting velocity kicks to the nearest neighbouring gas cells of the central black hole. The temperatures of the different regions are an indicator of the cooling processes that are likely to dominate in each phase if radiative cooling is taken into account. The shocked wind gas in the energy-driven outflow should cool efficiently via inverse Compton cooling close to the AGN (see text), whereas the shocked medium component should only cool via standard two-body processes such as Bremsstrahlung or metal-line cooling.

4.2 Cooling of the forward shock

Taking (assuming a fully ionised primordial gas) and for the shocked medium phase as in our simulated energy-driven outflow, Eq. 27 gives a free-free cooling timescale of . Only if the cooling time is comparable to the flow time , can the shell of shocked material fragment into clumps of colder material. Once the gas cools to and hydrogen is no longer ionised, further cooling must proceed through molecular and metal cooling processes, not included in our numerical simulations. The characteristic timescales for these processes are however expected to be very short (, see Zubovas & King 2014) and validate the approximation taken in this study that gas at traces molecular gas. Observed molecular outflows appear to have momentum fluxes (Sturm et al. 2011; Cicone et al. 2014), suggesting that they must be interpreted as energy-driven within the picture in which AGN energy/momentum and interstellar gas couple hydrodynamically (King et al. 2011; Faucher-Giguère & Quataert 2012; Zubovas & King 2014). We accordingly focus solely on simulations of energy-driven outflows in this section.

Our simulation setup is identical to that described in Section 3.2, but we now also include radiative cooling. We assume the cooling function presented in Katz et al. (1996), which is suitable for a primordial mixture of hydrogen and helium but probably underestimates the cooling of the shock-heated medium likely to be significantly enriched with metals and dust due to efficient star formation. The presence of metals and dust in the outflowing gas should greatly decrease the cooling times. In order to bracket the likely efficiency of radiative cooling, we therefore run an additional simulation with the same cooling function, increased by a factor of . For temperatures in the range of , the factor of increase results in a cooling function similar to that of gas enriched to solar metallicity (Sutherland & Dopita 1993).

In our various numerical tests, we also explore relaxing the assumption that the AGN constantly emits at its Eddington luminosity. The ratio is sensitive to the mass of the black hole and to the luminosity of the AGN that drives the outflow since . In order to account for the (degenerate) effects of lowering the black hole mass and/or AGN luminosity, we consider three different AGN light curves as shown in Fig. 9. The first case assumes the AGN to emit at a constant luminosity as before; in a second case we let the AGN luminosity rapidly rise from a value to and, in the third case, we adopt a variable light curve. In this last scenario, the chosen light curve is taken directly from a high resolution cosmological simulation that self-consistently follows black hole accretion (Costa et al. 2014). In the cosmological simulations, variability results from alternating accretion and outflow episodes, that respectively boost and suppress the AGN luminosity. For brevity, we shall name the simulations as ‘EnergyEdd’, ‘EnergyStep’, ‘EnergyVar’ and ‘EnergyEddCooling’ for simulations of energy-driven outflows adopting the identically named lightcurves in Fig. 9 respectively. EnergyEddCool is identical to EnergyEdd, but is simulated using the cooling function increased by a factor of . Note that in this section, we implicitly assume that an (unresolved) AGN inner wind is launched even if the Eddington ratio is low. We caution, however, that it is unclear whether an AGN inner wind can be launched if for instance if the driving mechanism is radiation pressure on UV lines (e.g. Proga & Kallman 2004; Proga 2005).

Figure 9: Eddington ratio as a function of time for the AGN light curves considered: () the AGN emits at its Eddington rate throughout the entire simulation (black), () the AGN initially emits at a hundredth of its Eddington rate, but rises sharply to its Eddington limit about half-way through the simulation and () the AGN luminosity curve varies with time as predicted by cosmological simulations (Costa et al. 2014).
Figure 10: The origin of multi-phase structure in outflowing material for energy-driven outflows with radiative cooling. The plot at the top shows the cumulative mass in outflowing gas as a function of temperature in the innermost of simulations EnergyEdd, EnergyStep, EnergyVar (primordial cooling) and EnergyEddCooling ( times more efficient cooling) for the light curves shown in Fig. 9 at a time where the bulk of the outflow is located at about from the AGN. Large quantities of cold gas () at form for EnergyStep, EnergyVar and EnergyEddCooling, where cooling times are comparable with the flow times. The plot at the bottom shows the flow time (estimated as ) (solid lines) and the cooling time as a function of temperature in the different simulations. Since they yield similar results, the region enclosed by the maximum and minimum cooling times and maximum and minimum flow times are shaded in blue and by red lines respectively for simulations EnergyEdd, EnergyStep and EnergyVar. Green lines show the cooling (solid) and flow (dashed) times for EnergyEddCooling.
Figure 11: Projection of gas mass along a slab of thickness and length for three different simulations using the energy-driven outflow model and the different light curves presented in Fig. 9. Note that due to the absence of any significant amount of outflowing cold gas in the simulation using standard cooling, here we show results for the run with the cooling function times increased. We show the output for snapshots corresponding to the time at which the outflow has reached radii of about in each simulation. Black contours enclose regions of outflowing cold gas with . In all cases, large amounts of cold gas form in the outflow.
Simulation time
Table 1: Properties of entrained cold gas () in simulations EnergyEddCooling, EnergyStep and EnergyVar at a time when the bulk of the outflow is at a distance of about from the AGN. From left to right, we provide values for the simulation time, the Eddington ratio of the AGN, followed by the total mass, the average speed, maximum speed, the outflow rate, the momentum flux normalised to of the AGN and the kinetic luminosity normalised to for outflowing cold gas. Results for simulation EnergyEdd are not included, because there is no cold gas present in this simulation at any time. The highest outflow rate, momentum boost and kinetic luminosity is obtained for simulation EnergyEddCooling, where cooling is most efficient, while AGN energy input is also very high.

Fig. 10 illustrates the impact of radiative cooling in our four numerical experiments. At the top of the panel, we show the cumulative mass of outflowing gas in the innermost of the halo as a function of temperature, at a time when the outflowing shell is located at a distance of about from the AGN. For simulation EnergyEdd (thin black line), the lowest temperature in outflowing gas is about , two orders of magnitude above the minimum temperature achievable due to primordial cooling. In this particular case, it does not seem possible to cool the shocked medium gas to temperatures at which molecules can form and the outflow is solely composed of hot ionised gas. However, in all other numerical tests, large amounts () of gas with form. We here term as ‘cold’ any gas with temperatures , though none of our conclusions is very sensitive to this threshold as most of the ‘cold gas’ has in fact reached the temperature floor of . As we have seen, gas at such temperatures would likely cool rapidly to much lower temperatures. Due to the very high cooling rates for gas in the shocked medium phase, simulation EnergyEddCooling produces of gas in the outflow. Large amounts of cold material () also form in simulation EnergyStep, where less efficient energy input from the AGN leads to a lower outflow speed and hence a lower . Finally, simulation EnergyVar, which has an identical cooling function as simulations EnergyEdd and EnergyStep, but an AGN luminosity intermediate between these two, correspondingly contains less cold material than EnergyStep but more than EnergyEdd. Note also that for every curve shown at the top of Fig. 10, there is a characteristic temperature at which the cumulative mass rises very sharply. In order of increasing temperature, this occurs first for run EnergyVar, followed by EnergyEddCooling, EnergyStep and finally EnergyEdd. This sharp rise in mass traces the uncooled component of the shell of swept-up interstellar gas. Its temperature is sensitive to the outflow speed, which determines the temperature of the gas that passes the forward shock, and on the efficiency of cooling.

At the bottom of Fig. 10, we explore the efficiency of cooling by plotting both and as function of outflowing gas temperature for all four simulations. We shade the region enclosed by the maximum and minimum cooling times in blue and the region enclosed by the maximum and minimum flow times in red for simulations EnergyEdd, EnergyStep and EnergyVar. Results for simulation EnergyEddCooling are shown as green lines as labelled on the plot. Note that for all simulations, cooling times rise extremely rapidly with increasing temperature. In all cases, for gas with . The hot bubble that drives the outflowing shell is therefore largely unaffected by cooling and the interpretation of the ‘subgrid’ model as realising the effects of a shocked wind phase remains valid. Cooling is efficient whenever and this can be seen to occur for simulations with standard cooling for outflowing gas at temperatures of . For simulation EnergyEdd, both timescales become comparable at , just above the virial temperature, while most outflowing gas has (see Fig. 10). Note also the slight bump in the flow time for all simulations at for EnergyEdd, EnergyStep and EnergyVar and at for EnergyEddCooling. In all cases, this traces the gas in the ambient medium which is currently passing through the forward shock. Whereas in the simulations with only primordial cooling, this gas is still at the virial temperature (see Fig. 8) and the bump therefore occurs at , in EnergyEddCooling, cooling is so efficient that it affects even the less dense ambient medium. The temperature of the gas that passes the forward shock in this case is therefore slightly lower, explaining the behaviour seen in Fig. 10.

Fig. 11 illustrates that the distribution of cold gas in the outflow is clumpy. In this panel, the location of cold clouds is shown as black contours plotted against the outflowing gas mass projected along a thin slab of thickness for simulations EnergyEddCooling, EnergyStep, EnergyVar at a time when the outflow has reached approximately in spatial extent. Gas preferentially cools in the over-densities seeded by Poisson noise in the initial conditions of the gaseous halo, which leads to the formation of the observed clumpy structure. Note that a departure from spherical symmetry in the outflow is present in all cases. This is particularly evident in EnergyStep, where the very low luminosity of the AGN has not been able to drive gas to more than away from the black hole, while the shocked medium has expanded adiabatically, leading to the diffuse appearance of the hot outflowing component. The presence of Rayleigh-Taylor instabilities, now accentuated by cooling in the shocked medium phase, is also clearly visible in Fig. 11 as hot and cold gas mix in the interface between shocked medium and shocked wind. Fig. 11 also shows the mass of cold gas found in the projected slab. In agreement with our previous discussion as well as Fig. 10, the largest quantities of cold outflowing gas are found in EnergyEddCooling, followed by EnergyStep and EnergyVar.

We quantify various properties of the simulated energy-driven (cold) outflows in Table 1. Even though average outflow speeds are in the range , maximum speeds can reach for sufficiently high AGN power, as is the case for EnergyEddCooling. In this case, outflow rates are very high (), momentum-boost factors reach (from Eq. 23, this value is similar to what would be expected by taking ) and kinetic luminosities are , in reasonable agreement with observational estimates for outflows driven by AGN powered by black holes of similar mass (e.g. Cicone et al. 2014). Note than in order to obtain simultaneously high momentum boost factors and high kinetic luminosities, both high outflow rates and high outflow speeds are required. In simulations EnergyStep and EnergyVar, outflow speeds are high (), but outflow rates are much lower than in EnergyEddCooling (only ). Both the momentum-boost and the kinetic luminosity are comparatively low in these cases ( and ). In particular, the relative inefficiency of the energy-driven outflow in EnergyStep despite the high AGN power () highlights the importance of efficient cooling. In this particular case and in general, metal-line cooling is probably crucial in generating quantities of cold material entrained in the outflow as large as observed. More efficient cooling leads to higher cold outflow rates and to higher momentum-boosts and kinetic luminosities in improved agreement with observed outflows.

5 Cosmological simulations

5.1 Feeding black holes from the cosmic web

In order to achieve efficient growth, supermassive black holes require a steady inflow of gas into their haloes at a high rate. Large amounts of gas can be brought in to the central galaxy as a result of a merger with another galaxy or as a consequence of filamentary inflow from the cosmic web (Di Matteo et al. 2008; Sijacki et al. 2009; Di Matteo et al. 2012; Khandai et al. 2012; Dubois et al. 2012, 2013b; Costa et al. 2014; Feng et al. 2014). Inflowing material slows down the propagation of outflows (see e.g. Nayakshin & Power 2010) as additional work to overcome the inflow is required and will reduce the efficiency of AGN feedback. Note also that a large portion of the ISM with which the luminous energy from the accreting black hole couples is cold unlike in the simple cases considered in Sections 2 and 3 and is distributed anisotropically around the AGN. Cold gas cools radiatively more efficiently and can lead to lower conversion efficiency between AGN luminosity and outflow kinetic energy, whereas an anisotropic mass distribution will cause the efficiency of feedback to be direction dependent. A full description of the growth and evolution of supermassive black holes and the propagation of outflows launched by AGN activity therefore ultimately requires a detailed account of their cosmological context. In this section, we do this by incorporating the ‘subgrid’ prescriptions that successfully reproduce the King (2003, 2005) models of energy- and momentum-driven winds in their idealised settings as presented in Section 3 into fully cosmological simulations. In order to highlight complications introduced by the cosmological environment of AGN, we compare the outflows driven in the cosmological simulations with those launched in an isolated halo with a gravitational potential equivalent to the spherically averaged (cosmological) halo. We assume a CDM cosmology with a set of parameters identical to Costa et al. (2014), i.e. , , , and . Due to the cosmological nature of the simulations, units are given in comoving coordinates in this section. Units of distance are accordingly prefixed with ‘c’ for ‘comoving’.

5.2 Numerical Setup

The numerical setup of the cosmological simulations presented in this study is very similar to that of Sijacki et al. (2009) and Costa et al. (2014) and here we only briefly review its main features. We employ output from the Millennium simulation (Springel et al. 2005c), which follows the collisionless dynamics of dark matter in a cubic cosmological volume of side length , in order to generate initial conditions for a QSO hosting halo at . The selected halo††††††The chosen halo corresponds to region ‘O’ in Costa et al. (2014). grows to a virial mass by and is thus a suitable host for a high redshift QSO (Sijacki et al. 2009; Angulo et al. 2012; Costa et al. 2014). The zoom-in technique (see Sijacki et al. 2009, for a more detailed description) is adopted in order to simulate a small region centred on the selected halo at high resolution. Resolution is gradually degraded with radial distance from the selected halo in order to accurately represent the large-scale cosmological tidal field, while maintaining computational efficiency. The high resolution region is then populated by gas cells and is evolved from a starting redshift down to .

The hydrodynamics of gas is followed using AREPO. Gas in our simulations is assumed to consist of an optically thin primordial mixture of hydrogen () and helium (). The N-body dynamics of dark matter is now followed self-consistently using a TreePM algorithm. In order to isolate the effects of AGN-driven outflows from other feedback effects, we deliberately limit the physical processes included in our simulations to radiative cooling (according to Katz et al. 1996) and star formation (following Springel & Hernquist 2003). Star-forming gas is thereby assumed to lie on an effective equation of state that accounts for the (self-regulated) balance of cold cloud formation via a thermal instability and subsequent evaporation in a hot component pressurised by supernova explosions. The resulting pressure gradients lead to dynamical stabilisation of the star forming gas, but do not trigger outflows. Note further that supernova-driven outflows are not incorporated in our simulations. Supernova-driven outflows have nevertheless been found not to launch efficient outflows from the very deep potential well of the quasar host halo. Instead, as they suppress star formation in less massive galaxies surrounding the QSO, supernova-driven outflows generate a surplus of infalling material that causes more efficient black hole growth and hence more dominant AGN-driven outflows in the QSO host halo (see Costa et al. 2014). In order to avoid resorting to further ‘subgrid’ models, black hole accretion is also not followed explicitly. Instead we seed a black hole sink particle at the centre of the most massive halo when this reaches a mass of at . Our choice of black hole mass is in close agreement with that predicted by self-consistently following black hole growth from a seed black hole in cosmological simulations for the same halo and at the same redshift (Sijacki et al. 2009; Costa et al. 2014). Note that in a halo at this redshift, gas cooling is very efficient and the bulk of the baryons cannot form a quasi-static ‘atmosphere’. At the time of black hole seeding, the galaxy thus has a stellar mass of and a star formation rate of (estimated within the galaxy’s half mass radius). We estimate the D velocity dispersion of star particles within the half mass radius of the galaxy hosting the black hole to be . For this velocity dispersion, a black hole mass of agrees well with the expected value if the galaxy were to lie on the observed relation (, Kormendy & Ho 2013). The black hole mass is then kept fixed at this value for the entire duration of the simulations. We set the dynamical mass of the black hole to a value hundred times the mass of the dark matter particle () in order to prevent it from unphysically oscillating around the centre of the halo. Note that the position of the black hole particle is then simply used as a tracer of the location at which AGN feedback energy or momentum is to be injected. We then carried out separate simulations employing either energy- or momentum-driven outflow models as described in Section 3.

Figure 12: The formation and evolution of an energy- (top row) and momentum-driven (two bottom rows) outflows in cosmological simulations. The AGN is powered by a black hole with ( in the bottom row) and lies at the centre of a halo. The background rendering illustrates the gas density projected along a slab with thickness of . White lines show contours for the radial velocity of outflowing gas in the cosmological simulations, whereas yellow dashed lines show identical contours for the outflow in an identical spherically averaged halo simulated in isolation. Contour levels are , and for the top, middle and bottom panels respectively. Outflows driven in cosmological simulations for energy-driving with and momentum-driving with have a comparable spatial extent and are highly anisotropic when compared to the spherical outflows in the isolated haloes. Depending on direction, cosmological outflows propagate to lower or larger distances from the AGN than isolated halo outflows. In the case of momentum-driving with a black hole, no outflow is observed in the cosmological simulation, whereas clearly present in the isolated halo simulation, clearly showing that a momentum input is required to launch large-scale momentum-driven outflows in cosmological environments.
Figure 13: Mass outflow rate profiles per logarithmic unit distance for the AGN in isolated halo simulations (dashed lines) and in cosmological simulations (solid lines) at different output times (colour-coded as labelled). Note that all the quantities are given in comoving units. Outflows are confined to an expanding shell in the isolated halo for both the energy- and momentum-driven cases. For the energy-driven outflow in the cosmological simulation, the outflow is spatially extended, but roughly located at similar distance. Note that no such outflow is present in the cosmological simulation of a momentum-driven outflow powered by a black hole, in contrast with the isolated halo case. Substantially higher masses (here shown for as a black solid line) are required to launch large-scale momentum-driven outflows with . The enhancement in outflow rates seen for the momentum-driven outflow at is due to an accretion shock and occurs simultaneously at the same spatial scales in simulations with and without AGN feedback.

Initial conditions for the isolated halo were obtained by considering the spherically averaged profiles of dark matter, stars and gas in the cosmological halo at the time the black hole particle is seeded in its centre. Mass density profiles were then separately produced for the gas component and for the collisionless components (dark matter and stars). The collisionless component was treated as a static potential, as in the static halo models discussed in Sections 3 and 4. The gaseous component, whose distribution was obtained by sampling the (non-analytical) density profile, was further assumed to be in hydrostatic equilibrium. Subsequently, a black hole was seeded at the centre of the isolated halo and its position and mass were fixed for the duration of the simulations. Note that since the density profiles used to construct initial conditions for the isolated halo were produced using the comoving units adopted in cosmological simulations, all quantities output in our isolated halo simulations must therefore also be interpreted as being given in comoving units. For the gas particle mass and gravitational softening length in our isolated halo simulations, we adopt values identical to those taken for the cosmological simulations (see Costa et al. 2014), i.e. and respectively. In AREPO, the gravitational softening length is adaptive and is set by the cell size. The typical minimum cell size in our cosmological simulations is . Simulations including radiative cooling and either energy- or momentum AGN feedback were then performed for this halo separately. Note that our simulations for an isolated halo here differ from the simulations presented in Section 4 only in resolution, the non-analytical nature of the potential and the fact that gas no longer exactly follows dark matter.

5.3 Comparison with spherical isolated haloes

Fig. 12 shows a time sequence of the density field for our cosmological simulations for energy- (top panel) and momentum-driven (two bottom panels) outflow models. In contrast with the spherically symmetric galactic halo, the density field is very anisotropic and most matter is concentrated in thin filaments with an overall configuration that changes little over the time span covered in Fig. 12. Each plot is centred on the galactic disc‡‡‡‡‡‡This is only true at early times, since the outflow completely disrupts the disc after some time. hosting the black hole, which is seen edge-on at as indicated by the projected angular momentum vector shown in the top left corner of the top panel. The host galaxy is located at a prominent region of the density field at about the centre of the over-density and is continuously fed by filamentary infall of gas. It is also in close proximity to another massive galaxy (about south-east at ) with which it eventually merges. The location of the AGN host galaxy in a collapsing over-density favours the presence of copious amounts of gas that can fuel black hole accretion and is what ultimately enables the formation of a very massive black hole (Sijacki et al. 2009; Di Matteo et al. 2012; Costa et al. 2014).

Figure 14: Profile of total enclosed mass for energy-driven (top row) and momentum-driven outflows (bottom row). The left column shows results for isolated halo simulations and the right column for cosmological simulations. For energy-driven outflows, the innermost regions of the halo are emptied with similar efficiency in isolated halo and cosmological simulations, except at late times, when the outflow is more efficient in the isolated halo. Momentum-driving for a black hole mass of yields very different results for isolated halo and cosmological simulations. For the latter, the total gas mass in the central regions grows due to efficient inflows, whereas a shell is efficiently driven outwards in the isolated halo. A black hole mass and ten times higher yields results similar to energy-driven outflows for at matching times. The absence of radiative cooling also enables an outflow to efficiently propagate (dotted line). Note that all quantities are given in comoving units.

Overplotted on the density fields are radial velocity contours (see figure caption for the contour levels), which trace the geometry and spatial extent of the AGN-driven outflows. Yellow contours refer to the outflows driven in the simulation for the isolated halo, while white contours refer to the cosmological simulations. For the latter, energy-driven outflows are highly anisotropic and in striking contrast with the outflowing spherical shell driven in the isolated halo simulation. Initially, the cosmological outflow propagates in directions parallel to the angular momentum vector of the galaxy hosting the luminous AGN because it is confined by the galactic disc (see also Zubovas et al. 2011). At later times, the central regions of the halo are disrupted and the outflow becomes more spherical, though still more pronounced in regions of low density and less prominent in the direction of the filaments (see also Costa et al. 2014). Note that whether the outflow driven in cosmological simulations has travelled to higher or lower distances from the black hole than in the isolated halo now depends on direction. As we shall see, the fact that the outflow driven in cosmological simulations travels more inefficiently than for the isolated halo along certain directions means that outflows driven in cosmological simulations are generally less efficient at clearing gas away from the centre of the halo.

For the equivalent cosmological simulation with momentum injection for (second row), there is no discernible outflow. Note that in the isolated halo case, the same black hole mass however gives rise to a momentum-driven shell that efficiently evacuates the central regions of the halo, indicating that a momentum input is required to drive a large-scale outflow in a cosmological setting (cf. Silk & Nusser 2010). For momentum injection at a rate , a larger black hole mass is required, as shown in the bottom row of Fig. 12, where results for the cosmological simulation are shown for times matching those of the energy-driven outflow, but for a black hole mass of . Results for the isolated halo are given as in the first row, i.e. for a black hole mass of and for the energy-driven outflow model. Note that the momentum-driven outflow now extends over spatial scales comparable to the energy-driven outflow (at matching times for a black hole mass ten times lower), suggesting that both outflows are now similarly efficient. In this case, the momentum-driven outflow is also highly anisotropic and, as in the energy-driven outflow, inefficient in directions of filamentary infall.

Fig. 13 shows the outflow rate as a function of radius for simulations of energy- (left) and momentum-driven (right) outflows for a black hole in the cosmological (solid lines) and isolated halo (dashed lines) simulations. Here, all gas with positive radial velocity is taken into account. While for the isolated halo, this velocity cut ensures that only gas which flows out due to AGN feedback is being traced, for the cosmological simulations, results must interpreted as upper limits because motion of gas which is not associated with AGN feedback will also be inevitably included. Consistent with Fig. 12, Fig. 13 shows that in cosmological simulations, the energy-driven outflow is spatially extended and anisotropic and not adequately described by a shell-like geometry. The outflow rate levels are significantly above those occurring in the absence of an AGN (dotted line), but at any location lower than in a shell-like configuration. Fig. 13 also shows that the peak outflow rates in cosmological simulations systematically trail behind the outflowing shell in the isolated case. This inefficiency of energy-driven outflows in cosmological simulations compared to the isolated halo case is moderate, however, and does not prevent the outflow from effectively clearing the centre of the halo.

The right hand side of the panel in Fig. 13 shows the momentum-driven outflow case. In the isolated halo, a momentum input is clearly sufficient to drive a large-scale outflow, since a shell can be seen to propagate outwards. Note that the free expansion of a momentum-driven shell is expected in this case since the black hole and its host galaxy lie on the relation (see Section 2). Results for the cosmological simulations confirm that, as seen in Fig. 12, there is no significant bulk outflow and the gas instead piles up in the central regions of the halo. The associated increase in outflow rates in the central regions is due to an increase in the speed of random motions of gas in these regions. On the same plot, we show the outflow rates resulting from a momentum-driven outflow for a black hole mass of at . The outflow that forms in this case has reached a distance of from the black hole. The peak outflow rate is at , similar to results for an energy-driven outflow due to a black hole mass of at a similar simulation time of .

Fig. 14 shows the total enclosed mass as a function of radius for the isolated halo (left column) and cosmological (right column) simulations. Results for the energy-driven outflow are given in the first row and for the momentum-driven outflow in the second row. Consistent with results presented in Fig. 13, the innermost regions of the AGN host galaxy are emptied with comparable efficiency in cosmological and isolated halo simulations for energy-driven outflows for a black hole mass of at early times. We verified that the star formation rate of the AGN host galaxy accordingly drops to zero. For the momentum-driven outflows, the total mass of gas in the central regions increases with time, leading to an extremely high star formation rate of by in the AGN host galaxy. For a black hole mass of however, the now efficient outflow is in fact comparable to the energy-driven outflow for a black hole mass ten times lower at matching times and the central regions of the galaxy are cleared.

From Fig. 12, we had seen that even when able to propagate out to large distances, outflows are always far more (if not totally) inefficient along directions of filamentary infall. Thus, a fundamental requirement for efficient AGN feedback is the outflow’s ability in preventing the central regions from being replenished by these inflows. We next verified that in the absence of such (cold) inflows, a momentum input rate of is sufficient to clear the innermost regions of the galaxy. For this purpose we performed another cosmological simulation with momentum injection from a black hole mass of , but artificially suppressed the cooling function from the time at which the black hole is seeded. This procedure guarantees that the density distribution of this and the previous simulation with standard primordial cooling is identical at the time at which momentum injection starts taking place. Suppressing cooling equates to reducing the efficiency at which the galaxy is replenished with fresh material. The dotted black line in Fig. 14 shows that momentum-driving now efficiently expels gas from the inner regions of the galaxy.

The results presented in this section indicate that if the AGN energy/momentum couples to the ISM of the host galaxy hydrodynamically via an inner wind, an efficient large-scale outflow must be energy-driven already at scales of (minimum cell size). In the cosmological simulations, a momentum flux of is required to efficiently disrupt the central disc and revert the filamentary inflows from which it is fed, while is sufficient to drive an unbound shell for the spherically averaged isolated halo. In Fig. 15, we verify that the momentum fluxes and kinetic luminosities of energy-driven outflows resulting from a black hole in our cosmological simulations are comparable to observed values. Both quantities can be seen to rise with time as thermal energy is converted into kinetic energy of the outflow. The net momentum fluxes and kinetic luminosities, estimated by integrating over the curves shown in Fig. 15 are of the order of and , respectively, at scales in good agreement with observed outflows (e.g. Maiolino et al. 2012; Cicone et al. 2014).

Figure 15: Momentum flux and kinetic luminosity per logarithmic unit distance for an energy-driven outflow in the cosmological simulations at different times. Integrated quantities are labelled next to the corresponding curves. Note that, since units are given in comoving coordinates, these results imply that an outflow with momentum flux and kinetic luminosity builds up at scales of . The magnitude of the momentum fluxes, kinetic luminosities and spatial scales agree well with observed AGN-driven outflows and indicate that the energy-driven model is a viable mechanism for AGN feedback in cosmological simulations.

5.4 Implications for the origin of the relation

We have seen that in a realistic cosmological context, a much larger amount of momentum input is required to suppress the infall of (cold) gas than is available in the momentum driven model suggested by King (2003) based on spherical static haloes. This is not surprising as not only is the gas not static, but a significant fraction of the gas falls in along well defined filamentary structures with a small covering factor which does not offer a large working surface for the outward streaming AGN-driven outflow. A significantly larger momentum input than is therefore required to even temporarily interrupt the inflow. Most of the injected momentum and energy escapes to large radii with little effect on the majority of the inflowing gas. Note that similar conclusions are drawn by Nayakshin (2014); Zubovas & Nayakshin (2014) and Bourne et al. (2014) for a clumpy gaseous halo. Our cosmological simulations thereby suggest that the self-regulation of the black holes is not due to a complete ejection of the gas from the surrounding galactic haloes, but due to a temporary suppression/reduction of the inflow rate of cold gas. It is thereby not obvious at what radius the feedback loop limiting the black hole mass is closed. If cooling is efficient, the momentum inflow rate of infalling gas and the energy input rate required to eject the inflowing matter scale with and , respectively. The requirement of stopping the inflow and partially ejecting the gas should thus lead to a self-regulated relation with a slope somewhere in between these two cases, as is indeed observed. As discussed in considerable detail by Kormendy & Ho (2013), such a self-regulated relation is then almost certainly modulated by the late-time merging of then gas-poor galaxies, especially at the high-mass end. Note also, that in real galaxies is only weakly dependent on distance to the centre. The scale at which the feedback loop is closed should therefore have only a modest effect on an relation resulting from such self-regulation. This may explain why a wide range of feedback implementations in simulations of widely varying resolution appear to be able to reproduce at least approximately the slope of the observed relation. The main conclusion from our study regarding the relation is, however, that the observed AGN-driven winds must become energy-driven at sufficiently small galactic scales such that the outflow driven by the inner AGN wind is still sufficiently fast to convert about of the initially liberated energy into kinetic energy of the outflow. As we have explicitly shown there is then energetically no problem to realise momentum fluxes of in the outflow as observed. Note that the picture would be the same if a similar amount of energy was injected by a different means than via thermalisation of a fast inner AGN wind. A possibility to be explored further here is Compton heating in phases in which the AGN spectral energy distribution is particularly hard. Implications of our findings on future AGN feedback models are discussed in Appendix B.

6 Main conclusions and summary

  • We have confirmed that simplified analytical models of energy- and momentum-driven AGN outflows in a spherically symmetric isolated halo without cooling can be accurately reproduced in numerical simulations for similar assumptions. Our simulations, based on energy and momentum injection into the black hole’s nearest neighbour cells, give rise to expanding shells, with speeds and structure in close agreement with analytical solutions.

  • The velocity of the propagating shells in the hydrodynamical simulations are found to fall below those of the analytical solution when the shells move (sub)sonically through the halo’s gaseous medium. For a Hernquist halo with total mass , this occurs for momentum-driven outflows for black holes with mass at a distance of a few from the centre of the halo. This behaviour does however not alter the key conclusion of King (2003) that there is a critical mass below which the momentum input from the AGN cannot fully expel material from the centre of the halo. The necessary black hole mass to launch an unbound shell in the simulated halo was found to be that which is required to place the black hole and its host halo on the observed relation.

  • We have performed simulations of energy-driven outflows including radiative cooling in order to investigate the formation of cold material entrained in hot shocked gas in the outflowing shell. If the cooling time is comparable with the outflow time of the shell, large amounts of gas cool out of the expanding shell of shocked material. The radius at which this occurs depends sensitively on the outflow rate as well as the cooling properties of the gas. By artificially varying the AGN light curve (previously assumed to always correspond to the Eddington luminosity of the AGN) and the efficiency of cooling, we found that the amount of cold material can vary from no material able to cool, e.g. if the AGN emits at its Eddington limit but cooling is inefficient, to if either the AGN luminosity is sub-Eddington or if cooling is very efficient. In all cases, cold material entrained in the (energy-driven) outflows attains a specific kinetic energy close to the available specific thermal energy.

  • In cases in which only small quantities of cold gas form, the total momentum flux in cold gas is only a factor of a few times higher than the momentum flux of the AGN radiation field. In order to generate cold outflows with very high momentum fluxes and kinetic luminosities as observed, efficient cooling and high AGN luminosities are simultaneously required. Only in a simulation including a cooling rate ten times higher than that for primordial cooling and an Eddington limited AGN, do we find high momentum-boost factors and high kinetic luminosities . This suggests that efficient cooling, e.g. via metal-lines, is likely to play an important role in generating sufficient amounts of cold material in the otherwise hot outflows.

  • When included into cosmological simulations that account for the environment of rapidly growing supermassive black holes at , the adopted ‘subgrid’ prescriptions for AGN (energy and momentum) feedback give rise to highly anisotropic large-scale outflows. The simulated outflows propagate most efficiently along paths of least resistance. At small scales, the outflows preferentially move along directions perpendicular to the host galactic disc and at large-scales, they avoid regions of filamentary gas infall.

  • We have compared the efficiency of the AGN-driven outflows in cosmological simulations to our idealised spherically symmetric simulations by adopting the spherically averaged gravitational potential of the galactic halo in which the black hole is seeded in our cosmological runs. For matching black hole mass, energy-driven outflows remove gas from the innermost regions of the AGN host halo with comparable efficiency in cosmological and isolated halo simulations. At large radii (i.e. late times) however, energy-driven outflows are unable to disrupt the dense filaments that continuously transport infalling cold gas into the main halo.

  • The efficiency of momentum-driven outflows in cosmological simulations is drastically lower than in isolated halo simulations. In cosmological simulations, a momentum input of is unable to revert the inflow of fresh material into the centre of the galaxy, where the total gas mass instead increases with time, while the same input is sufficient to expell a momentum-driven shell for the isolated halo simulation. A momentum input of at least is required to prevent the central (resolved) regions of the AGN host galaxy to be replenished by inflows. For this to be energetically feasible, the outflows have to become energy-driven at sufficiently small galactic scales that the outflow driven by the shocked inner AGN wind rapidly converts about of the energy released by the AGN into kinetic energy of the outflow. Models in which AGN energy and momentum couple to the ISM radiatively would require very high and probably unrealistic effective infrared optical depths of at scales in order to produce efficient feedback.

By comparing idealised numerical simulations of galactic outflows expanding in spherical gravitational potentials, carefully calibrated against analytical solutions, we have shown that in realistic cosmological environments, about a factor ten larger momentum flow is required for AGN feedback to efficiently affect cooling-driven inflows and to self-regulate black hole growth. We have further demonstrated that such large momentum flow rates occur naturally if about 5% of the expected accretion luminosity of AGN is thermalised in the galaxy host on scales of ). The resulting energy-driven outflows contain substantial amounts of entrained cold gas which is cooling from the mostly very hot outflow in encouraging agreement with observed molecular outflows. The most plausible mechanism for the input of thermal energy at these distances from the central engine is the thermalisation of the kinetic energy of the ultra-fast nuclear outflows observed in many AGN. However, direct Compton heating during accretion phases resulting in exceptionally hard SEDs with Compton temperature also merits further study. With the dynamic range and resolution of cosmological simulations of AGN-driven outflows further improving, there is the exciting prospect of modelling the transition from ultra-fast nuclear to energy-driven galactic outflows and the effect of Compton cooling as well as Compton heating self-consistently in such simulations in the not too distant future.

7 Acknowledgments

We thank the anonymous referee for thoroughly reading the manuscript and for useful suggestions on how to improve it. We also thank Volker Springel, Cathie Clarke, Andrew King, Sergei Nayakshin, Andy Fabian, Martin Rees, Claudia Cicone and Roberto Maiolino for insightful comments and helpful discussions. TC thanks Matt Young, Mike Curtis and Laura Keating for useful discussions and suggestions and Amanda Smith for assistance with designing Fig. 1. All simulations presented in this study were performed using the Darwin High Performance Computer facilities based at the University of Cambridge, UK and at the DiRAC Complexity Cluster based at the University of Leicester, UK. TC is supported by an STFC studentship. MGH and TC acknowledge support from the FP7 ERC Advanced Grant Emergence-320596.


  • Aalto et al. (2012) Aalto S., Garcia-Burillo S., Muller S., Winters J. M., van der Werf P., Henkel C., Costagliola F., Neri R., 2012, A&A, 537, A44
  • Alatalo et al. (2011) Alatalo K. et al., 2011, ApJ, 735, 88
  • Angulo et al. (2012) Angulo R. E., Springel V., White S. D. M., Cole S., Jenkins A., Baugh C. M., Frenk C. S., 2012, MNRAS, 425, 2722
  • Arribas et al. (2014) Arribas S., Colina L., Bellocchi E., Maiolino R., Villar-Martín M., 2014, A&A, 568, A14
  • Booth & Schaye (2009) Booth C. M., Schaye J., 2009, MNRAS, 398, 53
  • Bourne & Nayakshin (2013) Bourne M. A., Nayakshin S., 2013, MNRAS, 436, 2346
  • Bourne et al. (2014) Bourne M. A., Nayakshin S., Hobbs A., 2014, MNRAS, 441, 3055
  • Bower et al. (2006) Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645
  • Cano-Díaz et al. (2012) Cano-Díaz M., Maiolino R., Marconi A., Netzer H., Shemmer O., Cresci G., 2012, A&A, 537, L8
  • Choi et al. (2013) Choi E., Naab T., Ostriker J. P., Johansson P. H., Moster B. P., 2013, MNRAS
  • Choi et al. (2012) Choi E., Ostriker J. P., Naab T., Johansson P. H., 2012, ApJ, 754, 125
  • Choi et al. (2014) Choi E., Ostriker J. P., Naab T., Oser L., Moster B. P., 2014, MNRAS
  • Churazov et al. (2005) Churazov E., Sazonov S., Sunyaev R., Forman W., Jones C., Böhringer H., 2005, MNRAS, 363, L91
  • Cicone et al. (2012) Cicone C., Feruglio C., Maiolino R., Fiore F., Piconcelli E., Menci N., Aussel H., Sturm E., 2012, A&A, 543, A99
  • Cicone et al. (2014) Cicone C. et al., 2014, A&A, 562, A21
  • Ciotti & Ostriker (1997) Ciotti L., Ostriker J. P., 1997, ApJ, 487, L105
  • Ciotti & Ostriker (2001) Ciotti L., Ostriker J. P., 2001, ApJ, 551, 131
  • Combes et al. (2013) Combes F. et al., 2013, A&A, 558, A124
  • Costa et al. (2014) Costa T., Sijacki D., Trenti M., Haehnelt M. G., 2014, MNRAS
  • Croton et al. (2006) Croton D. J. et al., 2006, MNRAS, 365, 11
  • Davis et al. (2014) Davis S. W., Jiang Y.-F., Stone J. M., Murray N., 2014, ApJ
  • Debuhr et al. (2011) Debuhr J., Quataert E., Ma C.-P., 2011, MNRAS, 412, 1341
  • Di Matteo et al. (2008) Di Matteo T., Colberg J., Springel V., Hernquist L., Sijacki D., 2008, ApJ, 676, 33
  • Di Matteo et al. (2012) Di Matteo T., Khandai N., DeGraf C., Feng Y., Croft R. A. C., Lopez J., Springel V., 2012, ApJ, 745, L29
  • Di Matteo et al. (2005) Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604
  • Dubois et al. (2013a) Dubois Y., Gavazzi R., Peirani S., Silk J., 2013a, MNRAS, 433, 3297
  • Dubois et al. (2013b) Dubois Y., Pichon C., Devriendt J., Silk J., Haehnelt M., Kimm T., Slyz A., 2013b, MNRAS, 428, 2885
  • Dubois et al. (2012) Dubois Y., Pichon C., Haehnelt M., Kimm T., Slyz A., Devriendt J., Pogosyan D., 2012, MNRAS, 423, 3616
  • Dyson & Williams (1997) Dyson J. E., Williams D. A., 1997, The physics of the interstellar medium
  • Fabian (1999) Fabian A. C., 1999, MNRAS, 308, L39
  • Fabian (2012) Fabian A. C., 2012, ARA&A, 50, 455
  • Fabian & Iwasawa (1999) Fabian A. C., Iwasawa K., 1999, MNRAS, 303, L34
  • Faucher-Giguère & Quataert (2012) Faucher-Giguère C.-A., Quataert E., 2012, MNRAS, 425, 605
  • Feng et al. (2014) Feng Y., Di Matteo T., Croft R., Khandai N., 2014, MNRAS, 440, 1865
  • Ferrarese & Ford (2005) Ferrarese L., Ford H., 2005, Space Sci. Rev., 116, 523
  • Ferrarese & Merritt (2000) Ferrarese L., Merritt D., 2000, ApJ, 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
  • Förster Schreiber et al. (2013) Förster Schreiber N. M. et al., 2013, ApJ
  • Gan et al. (2014) Gan Z., Yuan F., Ostriker J. P., Ciotti L., Novak G. S., 2014, ApJ
  • Ganguly et al. (2007) Ganguly R., Brotherton M. S., Cales S., Scoggins B., Shang Z., Vestergaard M., 2007, ApJ, 665, 990
  • Gaspari et al. (2011) Gaspari M., Melioli C., Brighenti F., D’Ercole A., 2011, MNRAS, 411, 349
  • Gebhardt et al. (2000) Gebhardt K. et al., 2000, ApJ, 543, L5
  • Genzel et al. (2014) Genzel R. et al., 2014, ApJ
  • Ghavamian et al. (2007) Ghavamian P., Laming J. M., Rakowski C. E., 2007, ApJ, 654, L69
  • Greene et al. (2011) Greene J. E., Zakamska N. L., Ho L. C., Barth A. J., 2011, ApJ, 732, 9
  • Gültekin et al. (2009) Gültekin K. et al., 2009, ApJ, 698, 198
  • Haehnelt et al. (1998) Haehnelt M. G., Natarajan P., Rees M. J., 1998, MNRAS, 300, 817
  • Harrison et al. (2014) Harrison C. M., Alexander D. M., Mullaney J. R., Swinbank A. M., 2014, MNRAS, 441, 3306
  • Hernquist (1990) Hernquist L., 1990, ApJ, 356, 359
  • Hopkins et al. (2008) Hopkins P. F., Cox T. J., Kereš D., Hernquist L., 2008, ApJS, 175, 390
  • Humphrey et al. (2010) Humphrey A., Villar-Martín M., Sánchez S. F., Martínez-Sansigre A., Delgado R. G., Pérez E., Tadhunter C., Pérez-Torres M. A., 2010, MNRAS, 408, L1
  • Ishibashi & Fabian (2012) Ishibashi W., Fabian A. C., 2012, MNRAS, 427, 2998
  • Ishibashi & Fabian (2014) Ishibashi W., Fabian A. C., 2014, MNRAS, 441, 1474
  • Ishibashi et al. (2013) Ishibashi W., Fabian A. C., Canning R. E. A., 2013, MNRAS, 431, 2350
  • Katz et al. (1996) Katz N., Weinberg D. H., Hernquist L., 1996, ApJS, 105, 19
  • Kawata & Gibson (2005) Kawata D., Gibson B. K., 2005, MNRAS, 358, L16
  • Khandai et al. (2012) Khandai N., Feng Y., DeGraf C., Di Matteo T., Croft R. A. C., 2012, MNRAS, 423, 2397
  • King (2003) King A., 2003, ApJ, 596, L27
  • King (2005) King A., 2005, ApJ, 635, L121
  • King (2010a) King A. R., 2010a, MNRAS, 408, L95
  • King (2010b) King A. R., 2010b, MNRAS, 402, 1516
  • King & Pounds (2003) King A. R., Pounds K. A., 2003, MNRAS, 345, 657
  • King et al. (2011) King A. R., Zubovas K., Power C., 2011, MNRAS, 415, L6
  • Kormendy & Ho (2013) Kormendy J., Ho L. C., 2013, ARA&A, 51, 511
  • Krug et al. (2010) Krug H. B., Rupke D. S. N., Veilleux S., 2010, ApJ, 708, 1145
  • Krumholz & Thompson (2013) Krumholz M. R., Thompson T. A., 2013, MNRAS, 434, 2329
  • Liu et al. (2013) Liu G., Zakamska N. L., Greene J. E., Nesvadba N. P. H., Liu X., 2013, MNRAS, 436, 2576
  • Lynden-Bell (1969) Lynden-Bell D., 1969, Nature, 223, 690
  • Macciò et al. (2008) Macciò A. V., Dutton A. A., van den Bosch F. C., 2008, MNRAS, 391, 1940
  • Magorrian et al. (1998) Magorrian J. et al., 1998, AJ, 115, 2285
  • Maiolino et al. (2012) Maiolino R. et al., 2012, MNRAS, 425, L66
  • Marconi & Hunt (2003) Marconi A., Hunt L. K., 2003, ApJ, 589, L21
  • Martin (2005) Martin C. L., 2005, ApJ, 621, 227
  • Martizzi et al. (2014) Martizzi D., Jimmy, Teyssier R., Moore B., 2014, MNRAS, 443, 1500
  • Martizzi et al. (2012) Martizzi D., Teyssier R., Moore B., 2012, MNRAS, 420, 2859
  • McCarthy et al. (2010) McCarthy I. G. et al., 2010, MNRAS, 406, 822
  • McConnell & Ma (2013) McConnell N. J., Ma C.-P., 2013, ApJ, 764, 184
  • McQuillin & McLaughlin (2012) McQuillin R. C., McLaughlin D. E., 2012, MNRAS, 423, 2162
  • McQuillin & McLaughlin (2013) McQuillin R. C., McLaughlin D. E., 2013, MNRAS, 434, 1332
  • Mo et al. (2010) Mo H., van den Bosch F. C., White S., 2010, Galaxy Formation and Evolution
  • Mullaney et al. (2013) Mullaney J. R., Alexander D. M., Fine S., Goulding A. D., Harrison C. M., Hickox R. C., 2013, MNRAS, 433, 622
  • Murray et al. (2005) Murray N., Quataert E., Thompson T. A., 2005, ApJ, 618, 569
  • Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  • Nayakshin (2014) Nayakshin S., 2014, MNRAS, 437, 2404
  • Nayakshin & Power (2010) Nayakshin S., Power C., 2010, MNRAS, 402, 789
  • Nayakshin & Zubovas (2012) Nayakshin S., Zubovas K., 2012, MNRAS, 427, 372
  • Nesvadba et al. (2011) Nesvadba N. P. H., Polletta M., Lehnert M. D., Bergeron J., De Breuck C., Lagache G., Omont A., 2011, MNRAS, 415, 2359
  • Okamoto et al. (2008) Okamoto T., Nemmen R. S., Bower R. G., 2008, MNRAS, 385, 161
  • Ostriker et al. (2010) Ostriker J. P., Choi E., Ciotti L., Novak G. S., Proga D., 2010, ApJ, 722, 642
  • Planelles et al. (2014) Planelles S., Borgani S., Fabjan D., Killedar M., Murante G., Granato G. L., Ragone-Figueroa C., Dolag K., 2014, MNRAS, 438, 195
  • Pounds & King (2013) Pounds K. A., King A. R., 2013, MNRAS, 433, 1369
  • Pounds & Reeves (2009) Pounds K. A., Reeves J. N., 2009, MNRAS, 397, 249
  • Pounds et al. (2003) Pounds K. A., Reeves J. N., King A. R., Page K. L., O’Brien P. T., Turner M. J. L., 2003, MNRAS, 345, 705
  • Power et al. (2011) Power C., Nayakshin S., King A., 2011, MNRAS, 412, 269
  • Pozdnyakov et al. (1983) Pozdnyakov L. A., Sobol I. M., Syunyaev R. A., 1983, Astrophysics and Space Physics Reviews, 2, 189
  • Proga (2005) Proga D., 2005, ApJ, 630, L9
  • Proga & Kallman (2004) Proga D., Kallman T. R., 2004, ApJ, 616, 688
  • Puchwein et al. (2008) Puchwein E., Sijacki D., Springel V., 2008, ApJ, 687, L53
  • Puchwein & Springel (2013) Puchwein E., Springel V., 2013, MNRAS, 428, 2966
  • Raimundo et al. (2010) Raimundo S. I., Fabian A. C., Bauer F. E., Alexander D. M., Brandt W. N., Luo B., Vasudevan R. V., Xue Y. Q., 2010, MNRAS, 408, 1714
  • Reeves et al. (2009) Reeves J. N. et al., 2009, ApJ, 701, 493
  • Rodríguez Zaurín et al. (2013) Rodríguez Zaurín J., Tadhunter C. N., Rose M., Holt J., 2013, MNRAS, 432, 138
  • Roth et al. (2012) Roth N., Kasen D., Hopkins P. F., Quataert E., 2012, ApJ, 759, 36
  • Rupke & Veilleux (2013a) Rupke D. S. N., Veilleux S., 2013a, ApJ, 775, L15
  • Rupke & Veilleux (2013b) Rupke D. S. N., Veilleux S., 2013b, ApJ, 768, 75
  • Sakamoto et al. (2014) Sakamoto K., Aalto S., Combes F., Evans A., Peck A., 2014, ApJ, 797, 90
  • 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
  • Sharma & Nath (2013) Sharma M., Nath B. B., 2013, ApJ, 763, 17
  • Sijacki & Springel (2006) Sijacki D., Springel V., 2006, MNRAS, 366, 397
  • Sijacki et al. (2007) Sijacki D., Springel V., Di Matteo T., Hernquist L., 2007, MNRAS, 380, 877
  • Sijacki et al. (2009) Sijacki D., Springel V., Haehnelt M. G., 2009, MNRAS, 400, 100
  • Sijacki et al. (2012) Sijacki D., Vogelsberger M., Kereš D., Springel V., Hernquist L., 2012, MNRAS, 424, 2999
  • Silk (2005) Silk J., 2005, MNRAS, 364, 1337
  • Silk (2013) Silk J., 2013, ApJ, 772, 112
  • Silk & Nusser (2010) Silk J., Nusser A., 2010, ApJ, 725, 556
  • 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
  • Springel (2010) Springel V., 2010, MNRAS, 401, 791
  • Springel et al. (2005a) Springel V., Di Matteo T., Hernquist L., 2005a, ApJ, 620, L79
  • Springel et al. (2005b) Springel V., Di Matteo T., Hernquist L., 2005b, MNRAS, 361, 776
  • Springel & Hernquist (2003) Springel V., Hernquist L., 2003, MNRAS, 339, 289
  • Springel et al. (2005c) Springel V. et al., 2005c, Nature, 435, 629
  • Sturm et al. (2011) Sturm E. et al., 2011, ApJ, 733, L16
  • Sutherland & Dopita (1993) Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253
  • Teyssier et al. (2011) Teyssier R., Moore B., Martizzi D., Dubois Y., Mayer L., 2011, MNRAS, 414, 195
  • Tombesi et al. (2012) Tombesi F., Cappi M., Reeves J. N., Braito V., 2012, MNRAS, 422, L1
  • 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. (2011) Tombesi F., Cappi M., Reeves J. N., Palumbo G. G. C., Braito V., Dadina M., 2011, ApJ, 742, 44
  • 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
  • Tremaine et al. (2002) Tremaine S. et al., 2002, ApJ, 574, 740
  • Ueda et al. (2014) Ueda Y., Akiyama M., Hasinger G., Miyaji T., Watson M. G., 2014, ApJ, 786, 104
  • Veilleux et al. (2013) Veilleux S. et al., 2013, ApJ, 776, 27
  • Villar Martín et al. (2014) Villar Martín M., Emonts B., Humphrey A., Cabrera Lavers A., Binette L., 2014, MNRAS, 440, 3202
  • Vogelsberger et al. (2014) Vogelsberger M. et al., 2014, MNRAS, 444, 1518
  • Weaver et al. (1977) Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377
  • Westmoquette et al. (2012) Westmoquette M. S., Clements D. L., Bendo G. J., Khan S. A., 2012, MNRAS, 424, 416
  • Wurster & Thacker (2013) Wurster J., Thacker R. J., 2013, MNRAS, 431, 2513
  • Yu & Tremaine (2002) Yu Q., Tremaine S., 2002, MNRAS, 335, 965
  • Zubovas & King (2012) Zubovas K., King A., 2012, ApJ, 745, L34
  • Zubovas & King (2014) Zubovas K., King A. R., 2014, MNRAS
  • Zubovas et al. (2011) Zubovas K., King A. R., Nayakshin S., 2011, MNRAS, 415, L21
  • Zubovas & Nayakshin (2014) Zubovas K., Nayakshin S., 2014, MNRAS, 440, 2625
  • Zubovas et al. (2013) Zubovas K., Nayakshin S., Sazonov S., Sunyaev R., 2013, MNRAS, 431, 793

Appendix A Convergence properties of outflow models

All results presented in this paper rely on the usage of two ‘subgrid’ models, whose goal is to realise energy- and momentum-driven AGN outflows at scales that can be accurately resolved in our simulations. In this Appendix, we summarise results from several numerical tests performed on our ‘subgrid’ models in order to test their robustness against changes in various numerical parameters.

Figure 16: Effect of varying the numerical resolution of simulations for energy- (left) and momentum-driven (right) outflows in a static Hernquist potential with total mass and black hole mass . Our adopted ‘subgrid’ models yield numerically converged solutions even for resolution elements. We conclude that both energy- and momentum-driven outflows are well converged in this study.
Figure 17: Effect of varying the number of cell neighbours into which the black hole deposits energy (left) and momentum (right). Here, we use resolution elements. The numerical solution resulting from our energy-driven model is robust against changes in the number of neighbour cells to the black hole and in each case agrees with the analytical solution. Our momentum-driven model is more sensitive to the number of cell neighbours and diverges from the analytical solution for lower number of neighbours. In this regime, gas cells are kicked with very high speeds and a higher time-stepping precision is required in order to bring such numerical solutions to agreement with the analytical predictions, as shown with blue triangles for a simulation with cell neighbours but a maximum time-step five times lower.

In Fig. 16, we examine the numerical convergence properties of our simulations. For this purpose, we performed three simulations for energy- and momentum-driven outflows for a static Hernquist potential with total mass (i.e. the same halo as discussed in the main body of this study) and a black hole mass of . In this set of simulations, we varied the number of resolution elements with which the gaseous halo is sampled, keeping all other paramaters fixed. We investigated simulations with ,