Unravelling the physics of multiphase AGN winds through emission line tracers
Alexander J. Richings, Claude-Andre Faucher-Giguere, Jonathan Stern
MMNRAS , 1–17 (2021) Preprint 24 February 2021 Compiled using MNRAS L A TEX style file v3.0
Unravelling the physics of multiphase AGN winds throughemission line tracers
Alexander J. Richings (cid:63) , Claude-Andr´e Faucher-Gigu`ere and Jonathan Stern , Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham, DH1 3LE, United Kingdom Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA) and Department of Physics and Astronomy,Northwestern University, 1800 Sherman Ave, Evanston, IL 60201, USA School of Physics and Astronomy, Tel Aviv University, Tel Aviv 69978, Israel
Accepted 2021 February 20. Received 2021 February 19; in original form 2020 December 11.
ABSTRACT
Observations of emission lines in Active Galactic Nuclei (AGN) often find fast( ∼ − ) outflows extending to kiloparsec scales, seen in ionised, neutral atomicand molecular gas. In this work we present radiative transfer calculations of emissionlines in hydrodynamic simulations of AGN outflows driven by a hot wind bubble,including non-equilibrium chemistry, to explore how these lines trace the physicalproperties of the multiphase outflow. We find that the hot bubble compresses theline-emitting gas, resulting in higher pressures than in the ambient ISM or that wouldbe produced by the AGN radiation pressure. This implies that observed emissionline ratios such as [O iv ] µ m / [Ne ii ] µ m , [Ne v ] µ m / [Ne ii ] µ m and [N iii ] µ m /[N ii ] µ m constrain the presence of the bubble and hence the outflow driving mecha-nism. However, the line-emitting gas is under-pressurised compared to the hot bubbleitself, and much of the line emission arises from gas that is out of pressure, thermaland/or chemical equilibrium. Our results thus suggest that assuming equilibrium con-ditions, as commonly done in AGN line emission models, is not justified if a hot windbubble is present. We also find that (cid:38)
50 per cent of the mass outflow rate, momentumflux and kinetic energy flux of the outflow are traced by lines such as [N ii ] µ m and[Ne iii ] µ m (produced in the 10 K phase) and [C ii ] µ m (produced in the transitionfrom 10 K to 100 K).
Key words: astrochemistry - galaxies: active - quasars: emission lines - quasars:general
Galaxies that host active galactic nuclei (AGN) have beenobserved to contain fast ( ∼ − ) outflows of gas onkiloparsec scales (e.g. Cecil 1988; Veilleux et al. 2003), likelydriven by the input of energy and momentum from accretiononto the central supermassive black hole (e.g. Feruglio et al.2010; Cicone et al. 2014).Given their high velocities, one might expect the out-flowing material to be almost exclusively hot ( (cid:38) K). How-ever, these outflows have been observed in a wide range ofemission and absorption lines, spanning molecular, neutralatomic and ionised gas phases (Fischer et al. 2010; Rupke& Veilleux 2011; Greene et al. 2012; Harrison et al. 2012,2014; Aalto et al. 2015; Fiore et al. 2017; Fluetsch et al.2019; Feruglio et al. 2020; Lutz et al. 2020).Several possible mechanisms for the formation of thecool phase in outflows have been considered (applicable to (cid:63)
Email: [email protected] both AGN- and star formation-driven winds), including theentrainment of cold gas clouds in a hot outflow (Scanna-pieco & Br¨uggen 2015; Gaspari & Sadowski 2017; Schneider& Robertson 2017; Gronke & Oh 2020), the direct accelera-tion of cold gas by cosmic rays or radiation pressure (Boothet al. 2013; Costa et al. 2018; Hopkins et al. 2020), or in-situcooling of the hot outflowing material (Wang 1995; Silich etal. 2003; Zubovas & King 2014; Costa et al. 2015; Thomp-son et al. 2016). While the true origin of the cool phase inthese powerful outflows remains uncertain, the observationalevidence for such a multiphase structure is nonetheless over-whelming (see also Veilleux et al. 2020 for a recent review).AGN-driven galactic outflows are important as theyare likely to play a vital role in shaping the formationand evolution of galaxies. The energy and momentum in-jected by black hole winds can regulate the growth ofthe black holes and the stellar component of their hostgalaxies, giving rise to the observed scaling relations be-tween the two (Silk & Rees 1998; King 2003; Murray etal. 2005; Zubovas & King 2012; Torrey et al. 2020). En- © a r X i v : . [ a s t r o - ph . GA ] F e b A. J. Richings, C.-A. Faucher-Gigu`ere and J. Stern ergetic feedback from AGN can quench star formation inthe most massive galaxies (Binney & Tabor 1995; Duboiset al. 2013; Bower et al. 2017), and is required by moderncosmological models of galaxy formation to reproduce ob-served galaxy populations in terms of their stellar masses(Bower et al. 2012; Crain et al. 2015; Tremmel et al. 2017;Weinberger et al. 2018; Dav´e et al. 2019), the colours ofmassive elliptical galaxies (Springel et al. 2005; Feldmann etal. 2017; Trayford et al. 2016), and the stellar mass densi-ties of massive galaxies (Choi et al. 2018; Wellons et al. 2020;Parsotan et al. 2020). Outflows can also enrich the circum-galactic medium (CGM) through the transport of metalsfrom the galaxy (Hummels et al. 2013; Ford et al. 2013; Mu-ratov et al. 2017; Tumlinson et al. 2017; Hafen et al. 2019),and can deplete CGM gas fractions by carrying baryons be-yond the virial radius (Davies et al. 2019; Oppenheimer etal. 2020).To understand how such outflows can influence the sur-rounding environment, we need to measure their energet-ics (e.g. mass outflow rates, momentum fluxes, and energyfluxes), to determine whether they contain sufficient energyand momentum to have a significant impact on their hostgalaxy. This requires us to connect the emission and absorp-tion line tracers in which the outflows have been observed tothe physical properties of their constituent gas components.Many studies have successfully employed photoionisa-tion models to reproduce the emission line properties ofAGN. Dopita et al. (2002) and Groves et al. (2004a) pre-sented photoionisation models of dusty clouds dominatedby radiation pressure in the Narrow Line Region (NLR) ofAGN. The pressure and density structure of the clouds weredetermined by assuming that they are in hydrostatic balancewith the radiation pressure exerted by the AGN, also knownas Radiation Pressure Confinement (RPC; see also Draine2011 and Yeh & Matzner 2012 for the same effect in ionisedgas around star forming regions). They then used the map-pings iii photoionisation and shock code to calculate theintensities of emission lines on a grid of density, metallicity,ionisation parameter and the power-law index of the ionis-ing spectrum. Observed line ratio diagnostics can then becompared to these models to deduce the physical propertiesof the line-emitting clouds (Groves et al. 2004b).Stern et al. (2014a) used the photoionisation code cloudy to model the emission from RPC clouds spanning alarge range of distances from the nucleus. They showed thatRPC can explain the observed optical emission line ratiosand the overlap of extended X-ray and optical line emissionin nearby Seyferts. Bianchi et al. (2019) later showed thatthese types of models can explain the emission measure dis-tribution of X-ray emission lines. RPC models also explainobserved Broad Line Region (BLR) emission line ratios overa range of 10 in AGN luminosity (Baskin et al. 2014) andthe broad ionisation distribution in AGN outflows (Sternet al. 2014b). Stern et al. (2016) used these RPC modelsand observed emission line ratios in luminous quasars toconstrain the relative importance of hot gas and radiationpressure on the dynamics of quasar outflows. They found noevidence for the compression of emission line gas expectedin the presence of a hot bubble on any scale, and arguedthat a dynamically important hot bubble can be ruled outon spatial scales below 40 pc.Such photoionisation models have had success in re- producing many of the observed emission line properties ofAGN. However, they assume pressure, thermal and/or chem-ical equilibrium in the line-emitting gas. These assumptionsneed to be tested further.In Richings & Faucher-Gigu`ere (2018a) (hereafterRFG18), we ran a suite of hydro-chemical simulations ofAGN winds to explore the origin of molecular outflows.We demonstrated that observed molecular outflows in AGNhost galaxies can be produced by the in-situ formation ofnew molecules within the outflowing material. In that work,we focussed on the molecular phase. However, the chemicalmodelling in these simulations also includes the chemistry ofthe atomic and ionised phases. This enables us to make di-rect predictions for the emission lines from all phases of theoutflow. Most previous hydrodynamic simulations of AGNoutflows do not directly model the line emission, which is im-portant as emission lines contain most of the observationalconstraints on the models. Our simulation suite thereforepresents a unique tool with which to study the connectionbetween the physical properties and energetics of multiphaseAGN outflows and the observable emission line tracers, andto test the underlying assumptions that enter into alterna-tive photoionisation models for AGN emission.In this paper, we compute emission lines from theionised, neutral atomic and molecular phases of the AGNoutflows in the RFG18 simulations, which we use to ex-plore the physical properties of the line-emitting gas. Theremainder of this paper is organised as follows. In Section 2we describe the simulations and the radiative transfer cal-culations used to model the line emission. In Section 3 wepresent predictions for the line emission from the simula-tions (Section 3.1), investigate the physical properties of theline-emitting gas (Section 3.2), and compare our model pre-dictions to observations (Section 3.3). We explore how theemission line ratios computed from our simulations can beused to constrain the driving mechanisms of AGN outflowsin Section 4, and in Section 5 we study the outflow energet-ics of the different gas phases traced by the various emissionlines. Finally, we summarise our results in Section 6. In RFG18 we presented a series of hydro-chemical simu-lations of AGN-driven galactic outflows. These simulationsmodel an initially uniform ambient medium, into which weinject an isotropic AGN wind by spawning gas particles inthe central parsec with an outward velocity of 30 000 km s − and a momentum injection rate determined by the AGNluminosity, L AGN . Each simulation follows the interactionof this wind with the ambient medium over 1 Myr, whichcorresponds to the typical flow times ( r/v ) of kiloparsec-scale outflows that have been observed in luminous quasars(e.g. Gonz´alez-Alfonso et al. 2017). These simulations usethe gravity+hydrodynamics code gizmo , with the MeshlessFinite Mass (MFM) hydro solver (Hopkins 2015) and a fidu-cial resolution of 30 M (cid:12) per gas particle.We model the time-dependent chemistry of the gas us-ing the chimes non-equilibrium chemistry and cooling mod-
MNRAS , 1–17 (2021) ultiphase AGN winds ule (Richings et al. 2014a,b), which follows the evolutionof 157 ions and molecules that dominate the cooling ratefrom cold ( ∼
10 K), molecular gas to hot ( > K), highlyionised plasmas. The chimes chemical network contains var-ious collisional reactions, including collisional ionisation, re-combination, charge transfer reactions, and reactions on thesurface of dust grains such as the formation of molecularhydrogen (assuming a constant dust-to-metals ratio), alongwith photoionisation and photodissociation reactions.For the photochemistry, we use the average quasar spec-trum from Sazonov et al. (2004) normalised according to thebolometric AGN luminosity and the distance of the gas par-ticle from the black hole. This is an average between an ob-scured and an unobscured spectrum, and is characteristic ofa typical quasar. We chose this spectrum as our simulationsfocus on outflowing material at kiloparsec scales, but do notexplicitly model the small-scale structures around the AGNthat are likely to contribute to the obscuration of the AGNradiation. By using an average spectrum in this way, ouraim was to capture the partial obscuration at small scales.This spectrum differs from that used in some other mod-els of AGN emission lines. For example, Stern et al. (2016)use a power-law spectrum representative of an unobscuredquasar, with a fiducial extreme UV (EUV) slope of − . ≈ n H = 10 cm − and solar metallicity , but they differ in thebolometric AGN luminosity, which is 10 and 10 erg s − ,respectively. For brevity, we will refer to these two runs asL45 and L46 for the remainder of this paper. The low-densityrun ( n H = 1 cm − ) from RFG18 did not cool and forma multiphase wind before the end of the simulation after1 Myr, while the low-metallicity run ( Z = 0 . Z (cid:12) ) stronglyunder-predicts the molecular outflow rates compared to ob-servations. To create maps of the emission lines from our simulations,we post-process the simulation snapshots using version 0.40of the publicly available Monte Carlo radiative transfercode radmc-3d (Dullemond et al. 2012), using the non- https://richings.bitbucket.io/chimes/home.html Throughout this paper, we use the solar abundances listed intable 1 of Wiersma et al. (2009), for which the solar metallicityis Z (cid:12) = 0 . equilibrium ion and molecule abundances that were cal-culated during the simulations with the chimes chemistrymodule. While the full 3D spherical outflow is included inthe simulations, only one octant of the simulation volumeuses the highest resolution level, with the remainder of thevolume using 8 × lower mass resolution. We therefore onlyuse the high-resolution octant to produce the emission linemaps. This octant is mirrored in the line of sight direction,to capture both the receding and approaching sides of theoutflow, as viewed by the observer. The resulting emissionmaps thus cover one quadrant of the outflow.As radmc-3d is a grid-based code, we first project thegas particles from the simulations on to an Adaptive MeshRefinement (AMR) grid, which is refined such that no cellcontains more than 8 particles. The particles are smoothedusing a cubic spline kernel with a smoothing length enclosing32 neighbours, as used in the MFM hydro solver for thesesimulations. When projecting the gas temperatures and ve-locities on to the grid, we weight the contribution from eachparticle by the given ion or molecule abundance, to avoidunphysical effects from mixing particles with very differentproperties in the same cell (see the discussion in section 5 ofRFG18).We calculate the emissivities of H α and H β includingboth recombination of H ii and collisional excitation of H i ,using the cascade matrix formalism described in Raga etal. (2015), with the atomic data and fits for the collisionstrengths and recombination coefficients that they present intheir Appendix A. For all other ions and molecules, the levelpopulations are computed by radmc-3d using an approxi-mate non-LTE approach based on the Local Velocity Gra-dient (LVG) method. We utilised atomic data (energy lev-els, transition probabilities and collisional excitation rates)from version 7.1 of the chianti database (Dere et al. 1997;Landi et al. 2013), supplemented with additional collisionalexcitation rate data from the lamda database (Sch¨oier etal. 2005) for C i , C ii and O i .Dust absorption has little to no effect on the infraredemission lines, but those at optical and UV wavelengths arestrongly affected by dust. The idealised setup of these sim-ulations, in which the outflow is embedded within a denseambient interstellar medium (ISM), is representative of anAGN in the buried QSO phase, before the outflow has bro-ken out of the disk of its host galaxy. In observations ofburied QSOs, the outflows are very difficult to detect in op-tical and UV lines due to the strong dust absorption fromthe surrounding ISM. Even in systems where the outflowshave broken out of the dense medium, they are typicallyonly seen in the blue wing of the emission line (i.e. from thenear side of the outflow, moving towards the observer).When we include the effects of dust in our radiativetransfer calculations, we also find that the optical and UVlines are strongly suppressed, with the outflow only seenin the blue wings. However, in our idealised setup the hostgalaxy is represented by a box of initially uniform densitygas up to 2.4 kpc across, which does not include the tur-bulent structures we would expect to see in the ISM northe geometry of a galactic disc. We therefore cannot model https://home.strw.leidenuniv.nl/~moldata/ MNRAS000
10 K), molecular gas to hot ( > K), highlyionised plasmas. The chimes chemical network contains var-ious collisional reactions, including collisional ionisation, re-combination, charge transfer reactions, and reactions on thesurface of dust grains such as the formation of molecularhydrogen (assuming a constant dust-to-metals ratio), alongwith photoionisation and photodissociation reactions.For the photochemistry, we use the average quasar spec-trum from Sazonov et al. (2004) normalised according to thebolometric AGN luminosity and the distance of the gas par-ticle from the black hole. This is an average between an ob-scured and an unobscured spectrum, and is characteristic ofa typical quasar. We chose this spectrum as our simulationsfocus on outflowing material at kiloparsec scales, but do notexplicitly model the small-scale structures around the AGNthat are likely to contribute to the obscuration of the AGNradiation. By using an average spectrum in this way, ouraim was to capture the partial obscuration at small scales.This spectrum differs from that used in some other mod-els of AGN emission lines. For example, Stern et al. (2016)use a power-law spectrum representative of an unobscuredquasar, with a fiducial extreme UV (EUV) slope of − . ≈ n H = 10 cm − and solar metallicity , but they differ in thebolometric AGN luminosity, which is 10 and 10 erg s − ,respectively. For brevity, we will refer to these two runs asL45 and L46 for the remainder of this paper. The low-densityrun ( n H = 1 cm − ) from RFG18 did not cool and forma multiphase wind before the end of the simulation after1 Myr, while the low-metallicity run ( Z = 0 . Z (cid:12) ) stronglyunder-predicts the molecular outflow rates compared to ob-servations. To create maps of the emission lines from our simulations,we post-process the simulation snapshots using version 0.40of the publicly available Monte Carlo radiative transfercode radmc-3d (Dullemond et al. 2012), using the non- https://richings.bitbucket.io/chimes/home.html Throughout this paper, we use the solar abundances listed intable 1 of Wiersma et al. (2009), for which the solar metallicityis Z (cid:12) = 0 . equilibrium ion and molecule abundances that were cal-culated during the simulations with the chimes chemistrymodule. While the full 3D spherical outflow is included inthe simulations, only one octant of the simulation volumeuses the highest resolution level, with the remainder of thevolume using 8 × lower mass resolution. We therefore onlyuse the high-resolution octant to produce the emission linemaps. This octant is mirrored in the line of sight direction,to capture both the receding and approaching sides of theoutflow, as viewed by the observer. The resulting emissionmaps thus cover one quadrant of the outflow.As radmc-3d is a grid-based code, we first project thegas particles from the simulations on to an Adaptive MeshRefinement (AMR) grid, which is refined such that no cellcontains more than 8 particles. The particles are smoothedusing a cubic spline kernel with a smoothing length enclosing32 neighbours, as used in the MFM hydro solver for thesesimulations. When projecting the gas temperatures and ve-locities on to the grid, we weight the contribution from eachparticle by the given ion or molecule abundance, to avoidunphysical effects from mixing particles with very differentproperties in the same cell (see the discussion in section 5 ofRFG18).We calculate the emissivities of H α and H β includingboth recombination of H ii and collisional excitation of H i ,using the cascade matrix formalism described in Raga etal. (2015), with the atomic data and fits for the collisionstrengths and recombination coefficients that they present intheir Appendix A. For all other ions and molecules, the levelpopulations are computed by radmc-3d using an approxi-mate non-LTE approach based on the Local Velocity Gra-dient (LVG) method. We utilised atomic data (energy lev-els, transition probabilities and collisional excitation rates)from version 7.1 of the chianti database (Dere et al. 1997;Landi et al. 2013), supplemented with additional collisionalexcitation rate data from the lamda database (Sch¨oier etal. 2005) for C i , C ii and O i .Dust absorption has little to no effect on the infraredemission lines, but those at optical and UV wavelengths arestrongly affected by dust. The idealised setup of these sim-ulations, in which the outflow is embedded within a denseambient interstellar medium (ISM), is representative of anAGN in the buried QSO phase, before the outflow has bro-ken out of the disk of its host galaxy. In observations ofburied QSOs, the outflows are very difficult to detect in op-tical and UV lines due to the strong dust absorption fromthe surrounding ISM. Even in systems where the outflowshave broken out of the dense medium, they are typicallyonly seen in the blue wing of the emission line (i.e. from thenear side of the outflow, moving towards the observer).When we include the effects of dust in our radiativetransfer calculations, we also find that the optical and UVlines are strongly suppressed, with the outflow only seenin the blue wings. However, in our idealised setup the hostgalaxy is represented by a box of initially uniform densitygas up to 2.4 kpc across, which does not include the tur-bulent structures we would expect to see in the ISM northe geometry of a galactic disc. We therefore cannot model https://home.strw.leidenuniv.nl/~moldata/ MNRAS000 , 1–17 (2021)
A. J. Richings, C.-A. Faucher-Gigu`ere and J. Stern the realistic dust attenuation from the host galaxy in thissetup. Including dust absorption from the ambient mediumalso limits how much of the optical and UV emission we canstudy from the simulations.For the radiative transfer calculations in this paper, wetherefore include dust grains only in the outflow. For gaswith a radial velocity (cid:54) − (i.e. the ambient ISM,which is slowly inflowing by the end of the simulation dueto the gravitational potential of the host galaxy), we setthe dust density to zero. We also disable dust when the gastemperature is above 10 K, as grains will be rapidly de-stroyed by sputtering in this regime (e.g. Tsai & Mathews1995). This approach has little impact on the infrared lines.However, we stress that the resulting optical and UV lineluminosities do not include the effects of dust attenuationfrom the host galaxy. When comparing these lines to obser-vations, it is therefore important that we only consider lineratios at similar wavelengths, for which the strength of anyadditional dust absorption on each line will be similar, thusleaving the ratios unaffected. By excluding dust from theambient medium in this way, we can focus on the nature ofthe gas that produces each emission line.For the dust in the outflow, we use a mixture of graphiteand silicate grains with dust-to-gas ratios of 2 . × − and4 . × − , respectively, at solar metallicity. We include dustabsorption, scattering and thermal emission. To obtain thecontinuum emission, we repeat each radmc-3d calculationa second time with the emission lines disabled. The result-ing continuum emission is then subtracted from the totalemission. Fig. 1 shows the spectra of six infrared emission lines (toptwo rows) and three optical lines (bottom row) that are com-monly observed in AGN host galaxies. These represent asubset of the lines that we study in this paper; the full sam-ple of 24 spectra can be found in Appendix A. The dashedand solid curves in each panel are from the simulations L45and L46, respectively. Each spectrum is only integrated overthe spatial extent of the outflow, out to a radius of 0.57 kpc(L45) and 1.04 kpc (L46). The spatial extent in each caseis determined by the radius of the outflowing shell after 1Myr, when the simulation ends. Thus the L46 simulationreaches a larger spatial extent due to the higher outflow ve-locity in this case. We focus our analysis at 1 Myr becausethis corresponds to the typical flow times, t flow = r/v , ofoutflows observed in luminous AGN (e.g. Gonz´alez-Alfonsoet al. 2017). The fluxes are normalised to the maximum fluxin each spectrum.All of these emission lines show a broad componentfrom the outflowing shell, with velocities up to ≈
300 and ≈
700 km s − in L45 and L46, respectively. However, sev-eral of these lines show a strong narrow component aris-ing from the ambient medium, which is seen either in emis-sion (e.g. [C ii ] µ m ) or absorption (e.g. [O i ] µ m ). For theinfrared lines, we therefore only consider emission in thewings of each line for the rest of our analysis, with velocities | ∆ v | >
100 km s − , to focus on the outflow component.In the optical lines (and UV lines; not shown), the out- Figure 1.
Continuum-subtracted spectra of six infrared emissionlines (top two rows) and three optical emission lines (bottom row)commonly observed in AGN, from the simulations L45 (dashedcurves) and L46 (solid curves), normalised to the maximum fluxin each spectrum. These represent a subset of the lines that westudy in this paper; the spectra from the full sample of 24 lines canbe found in Appendix A. The vertical dotted line in each panel in-dicates the line centre. The broad component of these lines arisesfrom the outflowing shell, with velocities extending up to ≈ ≈
700 km s − in L45 and L46, respectively. However, centralvelocities at | ∆ v | (cid:54)
100 km s − are typically dominated by narrowcomponents from the ambient ISM. The grey shaded bands indi-cate velocity ranges that are excluded for the remainder of ouranalysis, to focus on the outflow component. flow is seen more strongly in the blue wing (approaching)than the red wing (receding) due to strong dust absorptionat these wavelengths, as noted in Section 2.2. We again stressthat we only include dust grains in the outflow for these ra-diative transfer calculations. This effect would be even moredramatic when we include the additional dust absorptionfrom the ambient ISM in the host galaxy. For our furtheranalysis of the optical and UV lines, we therefore only con-sider the blue wing, with velocities ∆ v< −
100 km s − . As the[S ii ] lines are a doublet, we also need to exclude velocities < −
450 km s − to avoid contamination of the 6731˚A line bythe 6716˚A line. We apply this additional velocity cut to bothof the sulphur lines. The velocity ranges that are excludedfrom the further analysis of each line are highlighted by thegrey shaded bands in Fig. 1.We can use these radiative transfer calculations to com-pare the spatial morphologies of each emission line from theoutflow. In Figs. 2 and 3 we show velocity-integrated mapsof the infrared and optical emission lines, where we includeemission from the red and blue wings (infrared) or the bluewing only (optical) as discussed above. The panels are each0.8 kpc and 1.2 kpc across in L45 and L46, respectively. We MNRAS , 1–17 (2021) ultiphase AGN winds Figure 2.
Velocity-integrated maps of the continuum-subtractedinfrared and optical emission lines from simulation L45, excludingthe velocity ranges indicated by the grey bands in Fig. 1 to avoidemission from the ambient ISM. Only a subset of lines is shownhere; the full sample can be found in Appendix A. Each panel is0.8 kpc across. The fluxes have been normalised to a distance of184 Mpc, i.e. the luminosity distance to the quasar Mrk 231. Ingeneral, molecular and low-ionisation species trace clumpy struc-tures in the outflowing shell, while intermediate ions trace a morediffuse phase, and the highest ionisation states arise from smallregions of bright emission.
Figure 3.
As Fig. 2, but for simulation L46. Each panel is 1.2 kpcacross. We again see the same trends in the spatial morphologiesof the different species as for L45.
Figure 4.
Temperature versus density of all gas in the high-resolution octant of the L46 simulation after 1 Myr. The colourscale indicates the distribution of the total gas mass in this plane.We have highlighted several regions in this plot to guide the dis-cussion. again only show a subset of the lines studied in this paper;the full sample can be found in Appendix A.The H , µ m line (top left panel of Figs. 2 and 3), whicharises from the pure rotational S(1) transition of the H molecule, traces dense, clumpy structures that have con-densed from the cooling material within the outflowing shell.Comparing to the maps of CO emission (see Appendix A),we find that the H , µ m line traces the same molecularstructures as the CO emission.In the top middle panel, the infrared [O i ] µ m line tracesthe same structures as the H , µ m emission. In contrast,the [O i ] optical line is somewhat more diffuse, eventhough both are produced by neutral atomic oxygen. Wewill see in Section 3.2 that the optical emission is onlypresent in the warmer gas phase ( ∼ K). This is unsur-prising, as the 6300 ˚A line requires a collisional excitationenergy
E/k B ∼ K, while the 63 µ m line requires 100 × lessenergy and can thus be excited in colder gas.The [C ii ] µ m emission (top right) is strong in theseclumpy structures, but there is also a significant contributionin the more diffuse regions between the dense clumps. Theintermediate ions (e.g. [N ii ] µ m and [Ne iii ] µ m ; left-handand centre panels of the middle row) predominantly arisefrom the diffuse, volume-filling phase in the outflowing shell,while the highest ionisation states such as [Ne vi ] µ m (right-hand panel of the middle row) are dominated by small,bright knots of emission. We will show in the next sectionthat these high ionisation states are produced in gas thatis undergoing a period of rapid cooling. As such, any oneparticular patch of gas spends a very short time produc-ing these high ionisation lines. The small, bright knots thatwe see in [Ne vi ] µ m are therefore short flashes of emissionas these compact regions transition through a rapid coolingphase. In this section we explore the physical properties of the gasprobed by each emission line. Before we break down the gasproperties by the separate emission line tracers, it is use-
MNRAS , 1–17 (2021)
A. J. Richings, C.-A. Faucher-Gigu`ere and J. Stern ful to look at the properties of the total gas distribution inthe simulations. Fig. 4 shows the temperature versus den-sity of all gas in the high-resolution octant of the L46 sim-ulation, where the colour scale indicates the distribution oftotal gas mass in each pixel. We have highlighted several re-gions in this plot to guide the discussion. We only show L46here, but the L45 simulation exhibits the same structuresin temperature-density space (see also Fig. 5 of RFG18).The evolution of gas through the temperature-density planeproceeds as follows:1) The ultra-fast outflow (30 000 km s − ) injected in thecentral parsec encounters the reverse shock (also calledthe wind shock) and is heated to ∼ K, creating a hotbubble. The pressure of this hot bubble is P hot / . k B = n H T = 2 . × cm − K. This pressure accelerates theoutflow, boosting the momentum beyond that expectedin a momentum-driven outflow (Faucher-Gigu`ere &Quataert 2012).2) The ambient ISM is intially in thermal equilibrium ata density n H = 10 cm − before it has been hit by theoutflow.3) When the ambient ISM is swept up by the forwardshock driven by the outflow, it is compressed andshock heated to ∼ . − K. This corresponds to thepost-shock temperature for a shock velocity of ≈ − − , according to the Rankine-Hugoniot jumpconditions for a strong shock, which is consistent withthe outflow velocities that we find in the simulations.This gas is in pressure equilibrium with the hot bubble.4) Some of the material in the swept up shell mixes withthe hot bubble, but remains in pressure equilibrium.5) The outflowing shell of gas swept up from the ambientmedium initially remains close to the post-shock tem-perature.6) Once the swept up gas reaches the peak in the cool-ing curve, the cooling time becomes very short and itrapidly cools. We can estimate the cooling time in thisregime based on the analytic models presented in Rich-ings & Faucher-Gigu`ere (2018b). In Fig. 6 of that work,we presented the temperature evolution of the sweptup shell of an isotropic AGN wind in a 1D analyticmodel based on Faucher-Gigu`ere & Quataert (2012).The fiducial model (black curves) uses the same AGNluminosity, ambient ISM density and metallicity as ourL46 simulation. In this model, the swept up shell coolsfrom 10 K to 10 K in a time t cool = 0 .
002 Myr. Forcomparison, the sound crossing time of the shell, witha thickness of ≈
25 pc (from the simulations) and usingthe sound speed at 10 K, is t sound = 0 . × less than the sound crossingtime in this regime, and so it cools isochorically, as thepressure of the surrounding hot medium has insufficienttime to compress the cooling gas and maintain pressureequilibrium during this phase.7) Gas that has just undergone rapid cooling and is ap-proaching thermal equilibrium is now under-pressurisedcompared to the surrounding hot medium.8) Once the gas reaches thermal equilibrium at a temper-ature ∼ K, it is subsequently compressed, increasingin density and pressure.9) At higher densities, the gas undergoes a thermal insta-
Figure 5.
Plots showing the temperature and density of gas pro-ducing the millimetre and infrared lines from the L46 simula-tion, arranged in order of increasing ionisation energy (or dis-sociation energy for molecules). The colour scale indicates thefraction of the total line luminosity that originates from gas ineach 2D temperature-density bin. The dashed, dot-dashed anddotted lines show the pressure of the hot bubble ( P hot , i.e. region(1) in Fig. 4), the radiation pressure ( P rad ; see equation 3) at theradius of the outflowing shell (1.04 kpc), and the pressure of theambient ISM ( P ISM ). The grey contours enclose 90 per cent ofthe total emission. The top axis shows the ionisation parameter,evaluated at 1.04 kpc. The hot bubble has compressed the line-emitting gas, resulting in higher pressures than in the ambientISM or that would be produced from radiation pressure alone.However, this gas is under-pressurised compared to the hot bub-ble itself, especially for the highest ionisation states, which arelocated in gas that has recently undergone a rapid cooling phase. bility and cools down to the cold phase at temperatures ∼
100 K. The cooling in this regime proceeds isobari-cally, as the 10 K and 100 K phases remain approxi-mately in pressure equilibrium with one another, albeitwith a large scatter. As we demonstrated in RFG18,this cold phase is conducive to the formation of newmolecules.
MNRAS , 1–17 (2021) ultiphase AGN winds Figure 6.
As Fig. 5, but for the optical and UV emission lines.Like the millimetre and infrared lines, emission at optical and UVwavelengths also arises from gas that is over-pressurised comparedto the ambient ISM and the radiation pressure alone, but under-pressurised compared to the hot phase.
Figs. 5 and 6 show the temperature-density distribu-tion of the line-emitting gas for the millimetre/infrared linesand the optical/UV lines, respectively. In these figures weshow the full sample of 24 emission lines, not just a sub-set as in previous figures. The colour scale indicates thefraction of the total line luminosity that is produced fromgas in each 2D temperature-density bin, normalised by thesize of the bin in logarithmic temperature-density space.These are calculated using the individual line luminositiesand ion/molecule-weighted temperatures and densities fromeach cell in the AMR grid used for the radiative transfer cal-culations (see Section 2.2). We only show the L46 simulationhere, but L45 exhibits the same trends. The grey contoursin each panel enclose 90 per cent of the total line emission.In Fig. 5, we see in the top row that the CO . emis-sion (from the J =1 − , µ m emis-sion (from the pure rotational S(1) transition) trace coldgas, at temperatures (cid:46) K. In contrast, the H , . µ m line,from the v1 − − J =2 − and the CO trace different phasesof molecular gas.The infrared lines of the low-ionisation states (e.g. O i and C ii ) arise from gas that is cooling from the warm (10 K)to the cold (100 K) phase, corresponding to region (9) inFig. 4. Intermediate ions such as N ii and O iii originate fromgas that lies on the thermal equilibrium branch at 10 K,i.e. region (8) in Fig. 4. Finally, the highest ionisation statessuch as Ne v and Ne vi are produced by gas that is comingto the end of the rapid cooling phase and is approaching thethermal equilibrium temperature.In Fig. 6 we have separated the hydrogen Balmer lineemission between that produced by the recombination of H ii (‘rec’; left panels of the top two rows) and that driven bythe collisional excitation of H i (‘col’; central panels of thetop two rows). The recombination radiation arises from gasundergoing the rapid cooling phase, similarly to the high-ionisation states. In contrast, collisional excitation radiationarises only from neutral gas on the 10 K branch. Raga etal. (2015) demonstrated that, at 10 K, the emission coeffi-cients of H α and H β from collisional excitation are up to fiveorders of magnitude higher than those due to recombination.This can lead to significant contributions from collisional ex-citations at these high temperatures even for fairly low H i fractions. The lack of such collisional excitation radiation at10 K in our simulations indicates that the H i fractions inthe rapid cooling phase are very low, < − . In L46, colli-sional excitation contributes 54 and 19 per cent to the totalH α and H β emission, respectively. In L45 (not shown), itcontributes 43 and 10 per cent to H α and H β , respectively.The 6300 ˚A emission from O i arises predominantly fromthe 10 K phase. This contrasts with the infrared lines fromO i , which we saw in Fig. 5 extends down to the cold phase.As we noted in Section 3.1, this is as we would expect, asthe 6300 ˚A line requires a higher excitation energy. Theoptical and UV emission from S ii , N ii , O iii and Ne iii alsotrace the 10 K phase, while the Ne v UV emission includesgas towards the end of the rapid cooling phase.To understand the role of photoionisation in the line-emitting gas, it is instructive to look at the ionisation pa-rameter, U ion . This is defined as the ratio of the densities ofhydrogen-ionising photons and hydrogen nuclei ( n γ and n H ,respectively), and can be calculated as follows: U ion ≡ n γ n H (1)= L . − πr cn H (cid:104) hν (cid:105) . − , (2)where L . − is the ionising luminosity from 13.6 eVto 1 keV, r is the distance from the AGN, c is the speedof light, and (cid:104) hν (cid:105) . − is the mean energy of ionis-ing photons in the 13 . − > (cid:104) hν (cid:105) . − = 32 . L . − /L AGN = 0 . L AGN is the bolometric AGN luminosity.Since the line-emitting gas is located in a thin ( ≈
25 pc)spherical shell, it is all at approximately the same distance of
MNRAS000
MNRAS000 , 1–17 (2021)
A. J. Richings, C.-A. Faucher-Gigu`ere and J. Stern r = 1 .
04 kpc (in L46) from the AGN. The ionising flux, andhence the density of ionising photons, is therefore approxi-mately uniform for this gas, and thus U ion depends simplyon the inverse of the gas density. We therefore show U ion at 1.04 kpc on the top axis of each panel in Figs. 5 and6. In general, emission from the highest ionisation statesarises from gas with 10 − (cid:46) U ion (cid:46) − , while intermediateions have 10 − (cid:46) U ion (cid:46) − . Gas traced by molecules andlow ionisation states extends to even lower ionisation pa-rameters, down to 10 − .The dashed, dotted and dot-dashed lines in Figs. 5 and6 indicate the hot gas pressure ( P hot , defined from region 1in Fig. 4), the initial pressure of the ambient ISM ( P ISM ) andthe radiation pressure from the AGN at a radius of 1.04 kpccorresponding to the outflowing shell ( P rad ), respectively.The radiation pressure can be calculated from the ionisingluminosity as follows: P rad = β L . − ion πr c , (3)where we again include ionising radiation only up to 1 keV,as the X-rays will penetrate beyond the ionised layer andtherefore do not contribute to the radiation pressure here.The correction factor β , introduced by Stern et al. (2014a),takes into account the contribution of non-ionising photonsto the radiation pressure. They show that β ≈ β ≈ β = 2 when calculating the radiation pressure.Note that, while photoionisation and photodissociationfrom the AGN radiation field are included in the thermo-chemistry in these simulations, direct radiation pressure it-self is not included. This would tend to compress gas at ther-mal pressures less than P rad (e.g. Dopita et al. 2002; Draine2011). We see in Figs. 5 and 6 that some of the emission linesarise from gas at pressures below P rad . We would thereforeexpect this gas to undergo additional compression if we wereto include radiation pressure in the simulations, but the vastmajority of emission in our simulations arises from gas with P gas (cid:38) P rad , and thus we do not expect the inclusion of radi-ation pressure in our simulations to significantly affect ourresults.We see that the hot bubble has compressed the lineemitting gas to become over-pressurised compared to theISM. For most emission lines, it also reaches higher pres-sures than would be achieved through direct radiation pres-sure alone. However, as this gas has cooled below the ini-tial post-shock temperature of the outflowing shell, it is nolonger in pressure equilibrium with the hot phase. The highionisation species (e.g. Ne vi ) trace gas that is most stronglyunder-pressurised compared to the hot medium, as it has nothad time to be compressed and so is still at relatively lowdensities where higher ionisation states are prevalent. In con-trast, the low- and intermediate-ionisation states trace gasat higher densities that have had longer to be compressedand are therefore closer to the pressure of the hot bubble,although they still remain under-pressurised.Photoionisation models of AGN emission lines such asGroves et al. (2004a) and Stern et al. (2014a) also assumethermal equilibrium in the line-emitting gas. Our simula- tions support this assumption for the intermediate lines suchas [N ii ] µ m , [N ii ] , [S ii ] , and [O i ] . InFigs. 5 and 6, these lines trace a narrow region at the ther-mal equilibrium temperature ∼ K. However, high ionisa-tion states such as Ne v and Ne vi arise from a broad rangeof temperatures at a given density (spanning up to ∼ − i ] , µ m and [C ii ] µ m that trace the 100 − K phasealso cover a broad range of temperatures at each density,rather than simply tracking the thermal equilibrium tem-perature as a function of density. This is likely due to tur-bulent shock heating. We find that, in our simulations, thedense ( n H >
100 cm − ) gas clouds that condense within theoutflowing shell after it cools typically have velocity disper-sions up to a few tens of km s − . The shocks resulting fromthis turbulence can heat the gas up to a few thousand K.The photoionisation models also assume that the ionsand molecules are in chemical equilibrium. To test this as-sumption in our simulations, we computed the equilibriumabundances of each gas particle at fixed density and tem-perature, and then repeated the radiative transfer calcula-tions using the equilibrium abundances. The most notabledifference is in the [O iii ] line, which is ≈ − v ] µ m (up to 80 percent), [N ii ] µ m (up to 62 per cent), [O iv ] µ m (up to 60per cent), [N iii ] µ m (up to 44 per cent), [S iii ] µ m (up to39 per cent), [N ii ] (up to 35 per cent), H
2; 2 . µ m (upto 31 per cent), H
2; 17 µ m (up to 29 per cent), and CO . (up to 26 per cent). The remaining lines differ by <
20 percent when we use equilibrium abundances.
We now compare the predictions from our simulations toobservations of outflows in AGN host galaxies. Fig. 7 com-pares the line luminosity ( L line ) divided by the bolometricAGN luminosity ( L AGN ) plotted against L AGN for a subsetof the infrared lines (the full sample can be found in Ap-pendix A). We divide by L AGN on the y-axis to reduce thedynamic range; a simple linear scaling between L line and L AGN would thus manifest itself as a horizontal trend inthese plots. For the simulations (RFG18; grey symbols), ourradiative transfer calculations include only one quadrant ofthe spherical outflow (see Section 2.2). We therefore multiplythe resulting luminosities by four to obtain the luminosityof the whole outflow.In the left-hand column of Fig. 7 we compare mid-IRlines from our simulations to a sample of observed AGNhost galaxies from Veilleux et al. (2009) (V09; squares). Wedivide the V09 sample according to the AGN classificationas a Palomar-Green Quasar (PG QSO), UltraLuminous In-fraRed Galaxy Seyfert 1 (ULIRG Sy1), ULIRG Sy2, or anULIRG Low Ionisation Nuclear Emission Regions (LINER).We exclude galaxies classified as H ii -like. The different clas-sifications are denoted by the colour. In the right-hand col-umn we compare far-IR lines from our simulations to theobservations from Herrera-Camus et al. (2018) (HC18; dia-monds), divided as a Sy1, Sy2, LIRG Sy2 and LIRG LINER, MNRAS , 1–17 (2021) ultiphase AGN winds Figure 7.
Line luminosity ( L line ) divided by bolometric AGNluminosity ( L AGN ) plotted against L AGN . Our simulations fromRFG18 are shown by the grey symbols. In the left column, show-ing a subset of mid-IR lines, we compare to the observations inVeilleux et al. (2009) (V09; squares), which we divide accord-ing to the AGN classification as PG QSO, ULIRG Sy1, ULIRGSy2 and ULIRG LINER, as indicated by the colour. In the rightcolumn, showing a subset of far-IR lines, we compare to the ob-servations in Herrera-Camus et al. (2018) (HC18; diamonds), di-vided according to classification as Sy1, Sy2, LIRG Sy2 and LIRGLINER (by colour). The full sample of all lines can be found inAppendix A. While many of these lines from the simulations over-lap with the observations, a notable exception is [O i ] µ m (topright), for which the simulations overpredict the luminosity by anorder of magnitude. again indicated by the colour. We exclude those classifiedas starburst galaxies in HC18. For the V09 sample, we usethe bolometric luminosities listed in their work. However,this information is not included in HC18. For this lattersample we therefore use the bolometric luminosity given inV09 where available, otherwise we derive it from the far-IR(40 µ m − µ m) luminosity using the bolometric correctionfrom Table 10 of V09 and assuming a 100 per cent AGNfraction, i.e. L AGN = 1 . L FIR .For many of the lines shown in Fig. 7, the two simula-tions overlap with the observations. A notable exception isthe [O i ] µ m line. While there are two LIRG LINER systemsthat are close to, or higher than, the [O i ] µ m luminosity Figure 8.
Infrared line ratios from our simulations (RFG18; greysymbols) and observed AGN host galaxies from Veilleux et al.(2009) (V09; square symbols) and Herrera-Camus et al. (2018)(HC18; diamond symbols). The observed samples are again di-vided according to the AGN classification, as in Fig. 7. Upperand lower limits are denoted by arrows. The dotted lines showidealised cloudy photoionisation models in hydrostatic equilib-rium based on Stern et al. (2016), using the Sazonov et al. (2004)average quasar spectrum (S04; dark blue) and the fiducial Sternet al. (2016) model with an ionising slope of − . P hot (cid:28) P rad (i.e. RPC) limit of these models. Inthe top-left panel, the simulations are consistent with Sy2-likeline ratios. However, the [Ne vi ] µ m / [O iv ] µ m and [Ne vi ] µ m / [Ne v ] µ m ratios are ≈ × higher in the simulations, while the[O i ] µ m / [C ii ] µ m ratios are ≈ × higher. predicted from L45, most of the observations lie at least anorder of magnitude below the simulations.Comparing to the individual AGN classifications, wesee that in the H , µ m and [Ne iii ] µ m lines the L45 sim-ulation is closest to the ULIRG LINER systems, while L46most closely overlaps with the ULIRG Sy1/Sy2 galaxies.In [Ne vi ] µ m L46 is again closest to the ULIRG Sy1/Sy2,while L45 is closer to the PG QSOs. In the far-IR lines[C ii ] µ m and [N ii ] µ m , it is more difficult to associatethe simulations with a particular AGN classification. Theobserved LIRG LINERs bracket the line luminosities pre-dicted by the simulations. However, in the observations thesesystems show a trend of decreasing L line /L AGN with increas-ing L AGN , which is not replicated in the simulations. Mostof the observed Sy1 and Sy2 galaxies have lower (higher)luminosities of [C ii ] µ m ([N ii ] µ m ) than the simulations,typically by a factor ≈
3, although the most extreme exam-ples of Sy1 and Sy2 galaxies in the observations bracket thesimulations.While many of the emission line luminosities predictedfrom the simulations overlap with the observations in Fig. 7,the large dynamic covered by the observations limits the con-
MNRAS , 1–17 (2021) A. J. Richings, C.-A. Faucher-Gigu`ere and J. Stern straining power in comparing the luminosities of individuallines. However, emission line ratios are often confined to anarrower range in observed systems, which provides tighterconstraints on the models. In Fig. 8 we therefore compareinfrared line ratios from the simulations (grey symbols) andthe observed samples from V09 (squares) and HC18 (dia-monds). The dotted curves show idealised cloudy photoion-isation models in hydrostatic equilibrium where we vary theratio of P hot /P rad , based on the models of Stern et al. (2016).The dark blue curves show a model using the same Sazonovet al. (2004) average quasar spectrum as used in our simu-lations, while the cyan curves show the fiducial Stern et al.(2016) model with an ionising slope of − .
6, typical of anunobscured AGN spectrum. The blue stars in Fig. 8 indicatethe P hot (cid:28) P rad limit of these models. As the cloudy calcula-tions end at a temperature of 100 K, they do not capture thefull [C ii ] and [O i ] emission in the cloud, so we do not includethe cloudy model predictions in the bottom-left panel ofFig. 8. We discuss these models further in Section 4.In the top-left panel of Fig. 8, the simulations lie on theobserved correlation between these two line ratios. Veilleuxet al. (2009) also demonstrated this correlation in the ob-servations, which represents an excitation sequence extend-ing from the ULIRG LINERs in the bottom left, throughthe ULIRG Sy2 and Sy1 galaxies and on to the PG QSOsin the top right. They note that this trend can be under-stood if [Ne v ] µ m and [O iv ] µ m are produced by the AGNwhile [Ne ii ] µ m comes from starburst activity. For fixed[Ne v ] µ m / [O iv ] µ m , varying the relative contribution ofthe starburst will move the ratios along the observed correla-tion, with the PG QSOs exhibiting the smallest contributionfrom the starburst, while the ULIRG LINERs contain thehighest starburst contribution. Our simulations lie closestto the observed ULIRG Sy2 systems, although we do notinclude contributions from a starburst component, whichwould tend to move the simulations down and to the leftin this plot.In the two panels showing ratios including [Ne vi ] (top-right and bottom-left), the simulations lie above the observa-tions by ≈ vi ] µ m line maybe too strong in our simulations. The predicted [Ne vi ] µ m /[Ne ii ] µ m and [Ne vi ] µ m / [Ne iii ] µ m ratios from the sim-ulations in these two panels overlap with the observations,however the observations cover almost two orders of mag-nitude and so would still be consistent with reducing the[Ne vi ] µ m emission in the simulations by ≈ iii ] µ m / [N ii ] µ m ratios from the simulations and observationsoverlap, [O i ] µ m / [C ii ] µ m is an order of magnitude toohigh in the simulations. As we saw in fig. 7, this is due tothe [O i ] µ m line, which is 10 × too high in the simulations.As noted in Section 2.2, our radiative transfer calcu-lations do not include dust attenuation from the ambientISM of the host galaxy. While this has little effect on theinfrared lines, the resulting luminosities of the optical andUV lines are much stronger than we would expect were suchattenuation to be included. Nevertheless, we can comparethe optical and UV line ratios to observations, as any ad-ditional dust attenuation would leave the ratios unaffected,provided that the wavelengths of the two lines are similar.In Fig. 9 we show BPT diagrams (Baldwin et al. 1981) Figure 9.
BPT diagrams from our simulations (black symbols),star-forming galaxies and type 2 AGN from SDSS (Kewley etal. 2006; grey 2D histogram), and type 1 AGN from SDSS withbolometric luminosities > erg s − (Stern & Laor 2013; bluecontours). The black solid and dotted curves show the bound-aries between star-forming and active galaxies from Kewley et al.(2001) and Kauffmann et al. (2003), respectively, while the blackdashed lines show the boundary between LINERs and Seyfertsfrom Kewley et al. (2006). The [S ii ] /H α and[O i ] /H α line ratios are up to an order of magnitudehigher in the simulations than is seen in the observations, due tostrong X-ray heating in the assumed partially-obscured incidentquasar spectrum and unresolved ionised layers. of optical emission line ratios, which are commonly usedto distinguish star-forming and AGN-dominated galaxies.The panels show the ratios of [N ii ] / H α (topleft), [S ii ] / H α (top right) and [O i ] / H α (bottom left) on the x-axis, plotted against[O iii ] / H β on the y-axis in each case. The sim-ulations are shown by the black symbols, while the grey2D histograms show the distribution of star-forming galax-ies and type 2 AGN in the Sloan Digital Sky Survey (SDSS)from Kewley et al. (2006) (updated to SDSS-DR8). The bluecontours show observed type 1 AGN in SDSS from Stern &Laor (2013), for which we only include those with a bolo-metric luminosity > erg s − to select the most power-ful AGN in this sample. The black solid and dotted curvesshow the boundaries between star-forming and active galax-ies from Kewley et al. (2001) and Kauffmann et al. (2003),respectively, while the black dashed lines show the boundarybetween LINERs and Seyferts from Kewley et al. (2006).In the top left panel panel of Fig. 9, the L45 simulationoverlaps with the observed AGN from SDSS and Stern &Laor (2013), although it lies close to the boundary betweenAGN and star-forming galaxies from Kewley et al. (2001).L46 would also be classified as an AGN from this panel,although the [O iii ] / H β is somewhat higher than MNRAS , 1–17 (2021) ultiphase AGN winds is seen in even the most extreme AGN in the observations,by up to a factor ≈ ii ] / H α and [O i ] / H α ratios are an order of magnitude higher than the observa-tions. The H α emission traces recombinations driven by pho-toionisation, while the [S ii ] and [O i ] emission traces cool-ing in the Warm Neutral Medium (WNM) phase at ∼ Kwhich is radiated via these metal lines. The high line ratiosin the BPT diagrams predicted by the simulations there-fore suggest that the ratio of the WNM to the photoionisedphases is likely to be too high in our simulations, by up toan order of magnitude.We find there are two effects that contribute to theseanomalous ratios. Firstly, the WNM in our simulations issupported by X-ray heating. If we recompute the thermalequilibrium temperature without the X-ray component ofthe incident radiation field from the quasar, we find thatmuch of the gas in the WNM cools from 10 K to a fewthousand K. As noted in Section 2.1, the average incidentquasar spectrum that we use from Sazonov et al. (2004)has a lower UV to X-ray ratio, by a factor ≈
2, than typi-cal unobscured quasar spectra used in other AGN emissionline models that successfully reproduce the BPT diagram(e.g. Stern et al. 2016). This means there will be relativelystronger X-ray heating, compared to photoionisation drivenby the UV band, which will increase the ratio of the WNMto photoionised phases and will thus contribute to the highline ratios that we find in the BPT diagrams.Secondly, the transition from ionised to atomic hydro-gen occurs at a column density N H , tot ≈ . cm − in oursimulations. However, at a density of n H =100 cm − and ourfiducial mass resolution of m gas =30 M (cid:12) , this is comparableto the column density of a single gas particle, which suggeststhat the photoionised layer in our simulations is barely re-solved. To estimate the extent to which we may be under-estimating the photoionised component due to limited res-olution, we can estimate the expected H α luminosity fromrecombinations of H ii given the incident ionising luminos-ity of the quasar. For the recombination emission coefficientthat we use from Raga et al. (2015), together with the ratecoefficient for case B recombination used in chimes , we ex-pect that each recombination of H ii will produce 0.27 H α photons at 10 K. If we assume that every incident ionisingphoton is absorbed and thus leads to one recombination,the resulting H α luminosity that we would expect given theionising luminosity is ≈ − ≈ − i ] µ m luminosity seen in Fig. 7. Figure 10.
Infrared line ratios versus the ratio of hot-gas toradiation pressure ( P hot /P rad ) in the simulations (black symbols)and an idealised cloudy photoionisation model in hydrostaticequilibrium based on Stern et al. (2016) (blue dotted curves).The horizontal blue lines indicate the P hot (cid:28) P rad (i.e. RPC) limitof the cloudy model. In the top four panels, the line ratiosfrom the simulations are close to the RPC limit. However, in thebottom three panels the simulations lie more than 1 dex belowthe RPC limit predictions. These three line ratios in the bottompanels therefore present the best opportunity to observationallydistinguish between the hot gas pressure scenario and the RPClimit. Stern et al. (2016) argued that, in AGN outflows driven bythe thermal pressure of a hot bubble, the line-emitting gasis compressed to higher pressures than we would expect inan outflow driven purely by the radiation pressure exertedby the AGN. They used idealised photoionisation calcula-tions assuming hydrostatic balance between the hot phase
MNRAS000
MNRAS000 , 1–17 (2021) A. J. Richings, C.-A. Faucher-Gigu`ere and J. Stern
Figure 11.
The line ratios [O iii ] / H β (left) andNe v / Ne iii (right) plotted against the ratio of hotgas pressure to radiation pressure ( P hot /P rad ). The black sym-bols show our simulations, while the blue dotted curves showthe line ratios in an idealised cloudy photoionisation model inhydrostatic equilibrium based on Stern et al. (2016) using theSazonov et al. (2004) average quasar spectrum. The line ratios inthe P hot (cid:28) P rad (i.e. RPC) limit of the cloudy model are indi-cated by the horizontal blue lines. The line ratios from the simu-lations bracket the RPC limit of the cloudy model, which meansthat we cannot use these particular ratios to distinguish betweenthe hot gas pressure dominated limit and RPC. and the line-emitting gas to show that this compression bythe hot bubble leads to a dependence of the emission lineratios on the ratio of P hot / P rad . By comparing these pre-dictions to observations of quasar outflows, they concludedthat, on large scales ( (cid:38) P hot /P rad > P hot (cid:29) P rad regime. However, note that in our simulationsthe outflow remains constrained by the dense ISM, whileStern et al. (2016) focussed on comparisons with unobscuredquasars in which the hot bubble may begin to vent out.As discussed in Section 3.2, the hot bubble in our sim-ulations does compress the line-emitting gas, as suggestedin Stern et al. (2016). However, the rapid cooling phaseleaves this gas under-pressurised compared to P hot . Wetherefore cannot assume hydrostatic balance between theline-emitting gas and the hot phase. We explore the impli-cations that this has for the emission line ratio constraintson P hot /P rad below.Fig. 10 shows infrared line ratios plotted against P hot /P rad , while Fig. 11 shows optical and UV line ra-tios versus P hot /P rad . The black symbols show our simula- tions, and the blue dotted curves show an idealised cloudy photoionisation model in hydrostatic equilibrium, includingdust grains, based on the Stern et al. (2016) models but us-ing the Sazonov et al. (2004) average quasar spectrum asin the hydrodynamic simulations. The horizontal blue linesindicate the Radiation Pressure Confinement (RPC) limitof the cloudy model, i.e. for P hot (cid:28) P rad .The strongest distinction between our simulations andthe RPC limit can be seen in the bottom three panels ofFig. 10, for the infrared line ratios of [O iv ] µ m / [Ne ii ] µ m ,[Ne v ] µ m / [Ne ii ] µ m and [N iii ] µ m / [N ii ] µ m . The pre-dictions for these ratios from the simulation are more than1 dex lower than the RPC limit of the cloudy model. Wetherefore conclude that these ratios present the strongestconstraints between the hot gas pressure-driven and radia-tion pressure-driven scenarios.However, not all line ratios can be used to constrain thedriving mechanism. In both panels of Fig. 11 the line ratiopredictions from our simulations bracket the RPC limit ofthe idealised cloudy model, with L45 and L46 producinglower and higher ratios than in RPC, respectively. If we ranfurther simulations covering a wider range of the parameterspace, we would expect some of these to produce interme-diate ratios between the two simulations shown here, whichcould thus overlap with the RPC predictions. In the top fourpanels of Fig. 10, the line ratios predicted by the simulationsare also close to the RPC limit. These particular line ratiostherefore cannot distinguish between the hot gas pressure-driven scenario and the RPC regime. This contrasts with thepredictions of the cloudy model in the P hot (cid:29) P rad limit, inwhich these line ratios become much lower than in the RPClimit. This discrepancy between the simulations and the hotgas pressure-dominated limit of the cloudy model is dueto the lower pressures that we find in the simulations com-pared to the assumption of hydrostatic balance with the hotphase.In general, the three infrared ratios most sensitive to thedriving mechanism (i.e. the bottom three panels of Fig. 10)compare high ionisation and low ionisation gas. We saw inFig. 5 that the high ionisation states are found in gas atlower densities that have had less time to undergo compres-sion, while the low ionisation states trace gas that exhibitstronger compression and are thus at higher densities. Sincethe compression by the hot bubble produces the dense low-ionisation gas, the ratios between high and low ionisationstates are most sensitive to the existence of the hot bub-ble and hence to P hot /P rad . However, line ratios such as[Ne vi ] µ m / [Ne iii ] µ m in the simulations are consistentwith the RPC limit.In Fig. 8 we compared the infrared line ratios predictedby the cloudy photoionisation models in hydrostatic equi-librium to observations. The dark blue curves show the samemodels as in Fig. 10, using the average quasar spectrumfrom Sazonov et al. (2004) as used in the simulations. Wealso show the fiducial model from Stern et al. (2016), withan ionising slope of − .
6, which is typical for an unobscuredquasar spectrum (cyan curves). The stars in Fig. 8 show theRPC limit of the cloudy models. In the top two panelsof Fig. 8, we see that the RPC limit predictions lie closeto the Palomar-Green Quasars (PG QSO). However, in thebottom-left panel the RPC limit predictions are offset fromthe observations by ≈ cloudy calculations MNRAS , 1–17 (2021) ultiphase AGN winds end at a temperature of 100 K, they do not fully capture the[C ii ] and [O ii ] emission, so we do not show these models inthe bottom-right panel. Comparing the dark blue and cyancurves, we see that the assumed quasar spectrum can affectthe line ratio predictions by up to a factor ≈ Emission lines are vital tools for measuring the energetics ofgalactic outflows. This is important for understanding howthe AGN can influence its host galaxy, for example whetherit can provide sufficient energy to unbind the gas in thegalaxy and hence quench star formation (e.g. Sturm et al.2011; Cicone et al. 2014; Fluetsch et al. 2019; Lutz et al.2020). However, we have seen in Section 3.2 that differentlines trace different phases of the outflow. If we only measurethe outflow energetics (e.g. mass outflow rate, momentumflux and energy flux) from a single line, we therefore wouldnot capture the total energetics of the outflow.In this section, we investigate what fraction of the totalenergetics are captured by each emission line. We calculatethis from the simulations as follows. First, we calculate themass outflow rate ( ˙ m part ), momentum flux ( ˙ p part ) and ki-netic energy flux ( ˙ E kin , part ) of each gas particle in the high-resolution octant of the simulation volume. We define thesequantities analogously to how they would be calculated inobservations as follows:˙ m part = m part v los , part r part . (4)˙ p part = ˙ m part v los , part = m part v , part r part . (5)˙ E kin , part = 12 ˙ m part v , part = 12 m part v , part r part , (6)where m part , v los , part and r part are the particle mass, line-of-sight velocity and distance from the black hole, respectively.For each of the ion and molecule species, we then projectthese quantities on to the AMR grid used for the radiativetransfer calculations weighted by the abundance of the givenspecies. Next, we calculate the emissivity ( (cid:15) ) of each emissionline in every cell of the AMR grid (i.e. the line luminosity perunit volume). We then sort the cells in order of increasing (cid:15) , and we find the emissivity (cid:15) for which 90 per cent ofthe total line luminosity orginates from cells with (cid:15)>(cid:15) .Finally, we calculate the ˙ m , ˙ p and ˙ E kin summed over cellswith (cid:15)>(cid:15) for each line, which we compare to the samequantities summed over all cells to find the fraction of thetotal energetics that are traced by each line.Defined in this way, we are measuring the fraction of theenergetics in cells that are ‘bright’ in a given emission line(such that 90 per cent of the total emission is included). Notethat a given cell can be bright in multiple lines, even fromdifferent ionisation states or molecules of the same element.For example, a cell containing a mixture of CO and C ii could be bright in both and thus would be accounted for in theenergetics of each line; we do not divide a cell between therelative contributions to each line. Thus the fractions thatwe present here should not be considered as a division ofthe total energetics between each separate species. Rather,they tell us how much of the energetics can be captured bylooking at each individual line, compared to how much arisesfrom gas that is dark in the given line and thus would notbe observable.The results of this analysis are presented in Fig. 12,which shows these fractions for each emission line. Thespecies are arranged in order of increasing ionisation energy(ions) or dissociation energy (molecules) along the x-axis.For each species, we plot the fractions traced by the line-emitting gas for the mass outflow rate ( ˙ m ), momentum flux( ˙ p ) and kinetic energy flux ( ˙ E kin ), denoted by the back-ground colour as light-, mid- and dark-blue, respectively.Results from the simulations L45 and L46 are shown by thegrey and black lines, respectively.We again caution that our predictions are likely to beinfluenced by uncertainties in the relative contributions ofthe photoionised and WNM phases, due to strong X-rayheating in the assumed partially obscured incident quasarspectrum and unresolved ionised layers, as discussed in Sec-tion 3.3. These results should therefore be considered as aqualitative indication of the relative importance of the differ-ent phases traced by each line, but the exact values remainuncertain.In general, the relative contributions of each line tothe total momentum flux and kinetic energy flux are simi-lar to the fraction of the total mass outflow rate that theytrace. This is unsurprising given the idealised setup of thesesimulations, as most of the outflowing mass is located in athin shell moving radially outwards with similar velocitiesacross the shell. There are some exceptions particularly in˙ E kin / ˙ E kin , tot , which can differ from ˙ m/ ˙ m tot and ˙ p/ ˙ p tot byfactors of a few in some cases (e.g. [C ii ] µ m , [Ne v ] µ m and[Ne vi ] µ m ). The kinetic energy flux has the strongest depen-dence on velocity, so variations in the velocity distributionsof different species will have a larger impact on ˙ E kin / ˙ E kin , tot than the other quantities. Nevertheless, the trends that wesee here are primarily driven by the fraction of the massthat is located in cells of the AMR grid that are bright inthe given line.The [C ii ] µ m line (produced in the transition from10 K to 100 K) and the [N ii ] µ m and [Ne iii ] µ m lines(produced in the 10 K phase) trace a large fraction of themass outflow rate and momentum flux ( ≈ −
70 per cent).The latter two lines also capture most of the kinetic en-ergy flux, although this is somewhat lower in [C ii ] µ m . Itmay seem surprising that [Ne iii ] traces such a high frac-tion of the energetics, given that it has an ionisation energyclose to species such as [O iii ], which trace much lower frac-tions of the energetics. Comparing these emission lines inFigs. 5 and 6, we see this is because [Ne iii ] emission ex-tends to higher densities than the [O iii ] lines. Photoionisa-tion models of RPC clouds using cloudy also exhibit sig-nificant Ne iii abundances extending into the neutral region(see e.g. Fig. A1 in Stern et al. 2014a).Emission from the high ionisation states, such as[O iv ] µ m , [Ne v ] µ m and [Ne vi ] µ m , only trace low frac-tions of the total energetics (typically (cid:46)
10 per cent). As we
MNRAS , 1–17 (2021) A. J. Richings, C.-A. Faucher-Gigu`ere and J. Stern
Figure 12.
Fractions of the total mass outflow rate ( ˙ m , light blue), momentum flux ( ˙ p , mid-blue) and kinetic energy flux ( ˙ E kin , dark blue)in gas traced by each emission line from the L45 (grey line) and L46 (black line) simulations, arranged in order of increasing ionisation(or dissociation) energy. The [C ii ] µ m , [N ii ] µ m and [Ne iii ] µ m lines each trace a significant fraction of the total energetics ( (cid:38) iv ] µ m , [Ne v ] µ m and [Ne vi ] µ m traces much less ( (cid:46)
10 per cent). saw in Section 3.2, these species arise from gas in a transi-tionary phase as it reaches the end of a period of rapid cool-ing, before it is subsequently compressed to higher densitiesdue to the external pressure exerted by the hot medium.Thus gas spends a relatively short period of time evolvingthrough this phase.We again stress that these results do not include dustattenuation from the host galaxy. This will not affect themillimetre and infrared lines, but we would expect those atoptical and UV wavelengths to be strongly absorbed. Thefractions presented in Fig. 12 for the optical and UV linesshould therefore be interpreted as the fraction of the energet-ics in gas that emits these lines, but might not necessarily beobserved if they are strongly attenuated by the host galaxy.Observations of AGN- and star formation-driven out-flows find prominent [C ii ] µ m wings, both locally and athigh redshift (e.g. Maiolino et al. 2012; Cicone et al. 2015;Janssen et al. 2016; Bischetti et al. 2019; Herrera-Camus etal. 2021). This supports our prediction that the [C ii ] µ m emission traces a large fraction of the outflow. Fluetsch etal. (2020) recently compared the relative contributions ofthe molecular, neutral atomic and ionised phases to the to-tal mass outflow rate and energetics in local ULIRGs. Thefractions that we present in Fig. 12 do not directly indicatethe relative contributions of each species to the total. In-stead, they show whether the emission from a given speciesis widely spread out over the whole outflow or concentratedin small regions. To compare our simulations to the datafrom Fluetsch et al. (2020), we therefore also measured therelative contribution of the three phases from the fraction ofthe total outflowing mass in H , H i and H ii . We thus foundthat, in simulation L45 (L46), the relative mass fractions were as follows: molecular – 25% (17%); neutral atomic –69% (71%); ionised – 5% (12%). Fluetsch et al. (2020) alsofound a negligible contribution from the ionised phase, sim-ilar to our simulations. However, the outflows in Fluetsch etal. (2020) were dominated by the molecular phase, with onaverage more than 60% of the mass outflow rate in H , andincreasing even higher in AGN-dominated systems. In con-trast, our simulations predict that the neutral atomic phasedominates. In RFG18 we found that the H fraction in thesimulations increases with increasing resolution, so it is pos-sible that the simulations may approach the observed H fractions at higher resolution. In this paper we explore the line emission from molecular,neutral atomic and ionised gas in simulations of multiphase,kiloparsec-scale outflows driven by the thermal pressure of ahot gas bubble generated by a central AGN. These simula-tions include an on-the-fly treatment for the non-equilibriumchemistry of ions and molecules coupled to the hydrodynam-ics and radiative cooling. The resulting ion and moleculeabundances are used together with a radiative transfer codein post-processing to make predictions for the emission lines.We use these calculations to study how this emission tracesthe physical conditions and energetics of the outflow, andwe compare the predicted emission lines from our models toobservations. Our main results are as follows:(i) We find that molecules (CO, H ) and low-ionisationstates (C ii and the infrared lines of O i ) trace clumpystructures that have condensed within the outflowingshell as it cooled, while the intermediate species (e.g. MNRAS , 1–17 (2021) ultiphase AGN winds N ii , S iii , Ne iii ) arise from a more diffuse phase through-out the shell. Emission from the highest ionisationstates (e.g. Ne v , Ne vi ) is concentrated in small, brightknots, which are produced by regions that are passingthrough a period of rapid cooling. As such, any oneparticular region of gas spends a relatively short pe-riod of time in this phase, and thus this high ionisationemission appears as short, bright flashes concentratedin small regions (see Figs. 2 and 3).(ii) Comparing the temperature-density distributiontraced by each emission line, we find that emissionfrom high ionisation states (e.g. Ne v , Ne vi ) is producedby gas that is coming to the end of a period of rapidcooling. Intermediate ions such as N ii and O iii arisefrom gas at the thermal equilibrium temperature of ∼ K, while low-ionisation states (e.g. C ii and theinfrared lines of O i ) and molecules trace the transitionfrom the warm ( ∼ K) to the cold ( ∼
100 K) phase(Figs. 5 and 6).(iii) The hot bubble compresses the line-emitting gas be-yond the initial pressure of the ambient ISM by 1 − ∼ K, the high ionisation states areup to ∼ i ] µ m line, which is an order of magnitudetoo high in the simulations (Fig. 7). Some of the pre-dicted line ratios at infrared wavelengths are also in-consistent with observations. For example, [Ne vi ] µ m /[O iv ] µ m and [Ne vi ] µ m / [Ne v ] µ m are 3 × too high,while [O i ] µ m / [C ii ] µ m is 10 × too high (Fig. 8).(vi) At optical wavelengths, the BPT diagnostic diagramsshow that the [S ii ] / H α and [O i ] / H α ratios are ≈ ≈ × higher X-ray to UV ratio than a more typical un-obscured spectrum, which enhances the Warm NeutralMedium (WNM) supported by X-ray heating. Secondly,we do not fully resolve photoionised layers, which resultsin an underestimate of the photoionised phase by a fac-tor ≈ −
4. Our simulations therefore do not provide areliable prediction for the BPT diagrams. We also cau-tion that these uncertainties may also affect other emis-sion line predictions from our simulations, in particularthose that compare the WNM and photoionised phases. (vii) As the line-emitting gas is compressed by the hot bub-ble, we find that certain line ratios are sensitive tothe ratio of the hot gas pressure to radiation pressure, P hot /P rad , providing a constraint on the driving mecha-nism of AGN outflows, as suggested by the photoionisa-tion models of Stern et al. (2016) that assumed thermal,chemical and hydrostatic equilibrium and an idealisedslab geometry. We find that the ratios of [O iv ] µ m /[Ne ii ] µ m , [Ne v ] µ m / [Ne ii ] µ m and [N iii ] µ m /[N ii ] µ m show the strongest distinction between oursimulations (for which P hot (cid:29) P rad ) and the radiationpressure-dominated limit of cloudy photoionisationmodels in hydrostatic equilibrium (Fig. 10). We there-fore suggest that these line ratios provide the strongestconstraints on the relative dynamical importance of ra-diation pressure and hot gas pressure on the outflow.(viii) We quantify the fraction of the total mass outflowrate, momentum flux and kinetic energy flux of the out-flow that is traced by each emission line. We find thatthe [N ii ] µ m and [Ne iii ] µ m lines (arising from the10 K phase) trace ≈ −
70 per cent of the totals in allthree quantities. [C ii ] µ m (arising from the transitionfrom 10 K to 100 K) also traces ≈ −
70 per cent ofthe mass outflow rate and momentum flux, but only ≈ −
40 per cent of the kinetic energy flux. Meanwhile,the high ionisation states such as [O iv ] µ m , [Ne v ] µ m and [Ne vi ] µ m (produced in a transitionary phase asthe gas undergoes rapid cooling) trace (cid:46)
10 per cent ofthe energetics (Fig. 12).Our simulations demonstrate that, in outflows drivenby the thermal pressure of a hot gas bubble (as can be pro-duced by the shock heating of a fast accretion disk wind,for example identified observationally as a BAL or UFO),the line-emitting gas is compressed by the hot phase, whichallows us to constrain the dynamical importance of the hotgas versus radiation pressure, provided that we consider lineratios that are most sensitive to this effect. We also findthat much of the emission arises from gas in highly transi-tionary phases, in particular for high ionisation states suchas Ne v and Ne vi that are produced towards the end of arapid cooling phase, which leads to the line-emitting gas be-ing out of thermal-, pressure- and chemical-equilibrium inmany cases. This highlights the importance of simulationsto capture these dynamical effects and their impact on emis-sion line predictions.Some of the emission line ratios predicted by our sim-ulations are inconsistent with observations of all AGN sub-types, as noted above. However, in this study we only havetwo simulations, which cover a very limited range of the pa-rameter space. These simulations also use an idealised setup,for example with an initially uniform ambient ISM that lacksturbulent gas structures or a realistic host galaxy geometry.The outflowing shell also remains constrained by the ambi-ent ISM throughout the simulations, and so we cannot probethe regime after the outflow breaks out of the dense regionsof the host galaxy. The tensions between our simulationsand the observations are therefore insufficient evidence torule out the hot gas pressure-driven scenario for kiloparsec-scale AGN outflows, as we cannot determine whether thesetensions are a limitation of the model or if they are simply MNRAS , 1–17 (2021) A. J. Richings, C.-A. Faucher-Gigu`ere and J. Stern due to the limited range and idealised nature of the simula-tions.To conclusively distinguish between possible drivingmechanisms of AGN outflows, we would therefore need torun a more extensive suite of simulations covering a widerange of the parameter space and different physical condi-tions. The resulting emission line predictions can then becompared to observations to constrain the models, similar tohow grids of photoionisation models have been employed toconstrain such models in the past (e.g. Groves et al. 2004b;Stern et al. 2016). With ongoing improvements to the com-putational codes used to run these simulations, we expectthis to be achievable in the near future.
ACKNOWLEDGEMENTS
DATA AVAILABILITY
The data underlying this article will be shared on reasonablerequest to the corresponding author. A public version ofthe gizmo code is available at . A public version of the chimes code is available at https://richings.bitbucket.io/chimes/home.html . REFERENCES
Aalto S. et al., 2015, A&A, 574, 85Baldwin J. A., Phillips M. M., Terlevich R., 1981, PASP, 93, 5 Baskin A., Laor A., Stern J., 2014, MNRAS, 438, 604Bianchi S., Guainazzi M., Laor A., Stern J., Behar E., 2019, MN-RAS, 485, 416Bischetti M., Maiolino R., Carniani S., Fiore F., Piconcelli E.,Fluetsch A., 2019, A&A, 630, 59Binney J., Tabor G., 1995, MNRAS, 276, 663Booth C. M., Agertz O., Kravtsov A. V., Gnedin N. Y., 2013,ApJ, 777, L16Bower R. G., Benson A. J., Crain R. A., 2012, MNRAS, 422, 2816Bower R. G., Schaye J., Frenk C. S., Theuns T., Schaller M.,Crain R. A., McAlpine S., 2017, MNRAS, 465, 32Cecil G., 1988, ApJ, 329, 38Choi E., Somerville R. S., Ostriker J. P., Naab T., HirschmannM., 2018, ApJ, 866, 91Cicone C. et al., 2014, A&A, 562, 21Cicone C. et al., 2015, A&A, 574, 14Costa T., Sijacki D., Haehnelt M. G., 2015, MNRAS, 448, L30Costa T., Rosdahl J., Sijacki D., Haehnelt M. G., 2018, MNRAS,479, 2079Crain R. A. et al., 2015, MNRAS, 450, 1937Dav´e R., Angl´es-Alc´azar D., Narayanan D., Li Q., RafieferantsoaM. H., Appleby S., 2019, MNRAS, 486, 2827Davies J. J., Crain R. A., McCarthy I. G., Oppenheimer B. D.,Schaye J., Schaller M., McAlpine S., 2019, MNRAS, 485, 3783Dere K. P., Landi E., Mason H. E., Monsignori Fossi B. C., YoungP. R., 1997, A&AS, 125, 149Dopita M. A., Groves B. A., Sutherland R. S., Binette L., CecilG., 2002, ApJ, 572, 753Draine B. T., 2011, ApJ, 732, 100Dubois Y., Gavazzi R., Peirani S., Silk J.. 2013, MNRAS, 433,3297Dullemond C. P., Juhasz A., Pohl A., Sereshti F., Shetty R., Pe-ters T., Commercon B., Flock M., 2012, Astrophysics SourceCode Library, record ascl:1202.015Faucher-Gigu`ere C-A., Quataert E., 2012, MNRAS, 425, 605Feldmann R., Quataert E., Hopkins P. F., Faucher-Gigu`ere C-A.,Kereˇs D., 2017, MNRAS, 470, 1050Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., King-don J. B., Verner E. M., 1998, PASP, 110, 761Feruglio C., Maiolino R., Piconcelli E., Menci N., Aussel H.,Lamastra A., Fiore F., 2010, A&A, 518, L155Feruglio C., Fabbiano G., Bischetti M., Elvis M., Travascio A.,Fiore F., 2020, ApJ, 890, 29Fiore F. et al., 2017, A&A, 601, 143Fischer J. et al., 2010, A&A, 518, L41Fluetsch A. et al., 2019, MNRAS, 483, 4586Fluetsch A. et al., 2020, arXiv:2006.13232Ford A. B., Oppenheimer B. D., Dav´e R., Katz N., Kollmeier J.A., Weinberg D. H., 2013, MNRAS, 432, 89Gaspari M., Sadowski A., 2017, ApJ, 837, 149Gonz´alez-Alfonso E. et al., 2017, ApJ, 836, 11Greene J. E., Zakamska N. L., Smith P. S., 2012, ApJ, 746, 86Gronke M., Oh S. P., 2020, MNRAS, 492, 1970Groves B. A., Dopita M. A., Sutherland R. S., 2004a, ApJS, 153,9Groves B. A., Dopita M. A., Sutherland R. S., 2004b, ApJS, 153,75Hafen Z. et al., 2019, MNRAS, 488, 1248Harrison C. M. et al., 2012, MNRAS, 426, 1073Harrison C. M., Alexander D. M., Mullaney J. R., Swinbank A.M., 2014, MNRAS, 441, 3306Herrera-Camus R. et al., 2018, ApJ, 861, 94Herrera-Camus R. et al., 2021, arXiv:2101.05279Hopkins P. F., 2015, MNRAS, 450, 53Hopkins P. F., Chan T. K., Ji S., Hummels C., Kereˇs D., QuataertE., Faucher-Gigu`ere C-A., 2020, arXiv:2002.02462Hummels C. B., Bryan G. L., Smith B. D., Turk M. J., 2013, 430,1548 MNRAS , 1–17 (2021) ultiphase AGN winds Janssen A. W. et al., 2016, ApJ, 822, 43Kauffmann G. et al., 2003, MNRAS, 346, 1055Kewley L. J., Dopita M. A., Sutherland R. S., Heisler C. A.,Trevena J., 2001, ApJ, 556, 121Kewley L. J., Groves B., Kauffmann G., Heckman T., 2006, MN-RAS, 372, 961King A., 2003, ApJ, 596, 27Lamperti I. et al., 2017, MNRAS, 467, 540Landi E., Young P. R., Dere K. P., Del Zanna G., Mason H. E.,2013, ApJ, 763, 86Landt H., Bentz M. C., Ward M. J., Elvis M., Peterson B. M.,Korista K. T., Karovska M., 2008, ApJS, 174, 282Liu G., Zakamska N. L., Greene J. E., Nesvadba N. P. H., Liu X.,2013, MNRAS, 430, 2327Lutz D. et al., 2020, A&A, 633, 134Maiolino R. et al., 2012, MNRAS, 425, L66Muratov A. L. et al., 2017, MNRAS, 468, 4170Murray N., Quataert E., Thompson T. A., 2005, ApJ, 618, 569Oppenheimer B. D. et al., 2020, MNRAS, 491, 2939Parsotan T., Cochrane R., Hayward C. C., Angl´es-Alc´azar D.,Feldmann R., Faucher-Gigu`ere C.-A., Wellons S., Hopkins P.F., 2020, preprint (arXiv:2009.10161)Raga A. C., Castellanos-Ram´ırez, Esquivel A., Rodr´ıguez-Gonz´alez A., Vel´azquez P. F., 2015, RMxAA, 51, 231Richings A. J., Schaye J., Oppenheimer B. D., 2014a, MNRAS,440, 3349Richings A. J., Schaye J., Oppenheimer B. D., 2014b, MNRAS,442, 2780Richings A. J., Faucher-Gigu`ere C.-A., 2018a, MNRAS, 474, 3673(RFG18)Richings A. J., Faucher-Gigu`ere C.-A., 2018b, MNRAS, 478, 3100Rosario D. J., Togi A., Burtscher L., Davies R. I., Shimizu T. T.,Lutz D., 2019, ApJL, 875, L8Rupke D. S. N., Veilleux S., 2011, ApJ, 729, 27Sazonov S. Y., Ostriker J. P., Sunyaev R. A., 2004, MNRAS, 347,144Scannapieco E., Br¨uggen M., 2015, ApJ, 805, 158Sch¨oier F. L., van der Tak F. F. S., van Dishoeck E. F., Black J.H., 2005, A&A, 432, 369Schneider E. E., Robertson B. E., 2017, ApJ, 834, 144Silich S., Tenorio-Tagle G., Mu˜noz-Tu˜n´on C., 2003, ApJ, 590, 791Silk J., Rees M. J., 1998, A&A, 331, L1Springel V., Di Matteo T., Hernquist L., 2005, ApJ, 620, L79Stern J. & Laor A., 2013, MNRAS, 431, 836Stern J., Laor A. & Baskin A., 2014a, MNRAS, 438, 901Stern J., Behar E., Laor A., Baskin A., Holczer T., 2014b, MN-RAS, 445, 3011Stern J., Faucher-Gigu`ere C.-A., Zakamska N. L., Hennawi J. F.,2016, ApJ, 819, 130Sturm E. et al., 2011, ApJ, 733, L16Sutherland R., Dopita M., Binette L., Groves B., 2013, Astro-physics Source Code Library, record ascl:1306.008Thompson T. A., Quataert E., Zhang D., Weinberg D. H., 2016,MNRAS, 455, 1830Torrey P. et al., 2020, MNRAS, 497, 5292Tumlinson J., Peeples M. S., Werk J. W., 2017, ARA&A, 55, 389Trayford J. W., Theuns T., Bower R. G., Crain R. A., Lagos C.d. P., Schaller M., Schaye J., 2016, MNRAS, 460, 3925Tremmel M., Karcher M., Governato F., Volonteri M., Quinn T.R., Pontzen A., Anderson L., Bellovary J., 2017, MNRAS,470, 1121Tsai J. C., Mathews W. G., 1995, ApJ, 448, 84Veilleux S., Shopbell P. L., Rupke D. S., Bland-Hawthorn J., CecilG., 2003, AJ, 126, 2185Veilleux S. et al., 2009, ApJS, 182, 628Veilleux S., Maiolino R., Bolatto A. D., Aalto S., 2020, A&ARv,28, 2Wang B., 1995, ApJ, 444, 590 Weinberger R. et al., 2018, MNRAS, 479, 4056Wellons S., Faucher-Gigu`ere C.-A., Angl´es-Alc´azar D., HaywardC. C., Feldmann R., Hopkins P. F., Kereˇs D., 2020, MNRAS,497, 4051Wiersma R. P. C., Schaye J., Smith B. D., 2009, MNRAS, 393,99Yeh S. C. C., Matzner C. D., 2012, ApJ, 757, 108Zubovas K., King A. R., 2012, MNRAS, 426, 2751Zubovas K., King A. R., 2014, MNRAS, 439, 400
APPENDIX A: FULL SAMPLE OF EMISSIONLINE SPECTRA
In the main text of this manuscript we presented only a sub-set of the emission line spectra (Fig. 1), velocity-integratedemission maps (Figs. 2 and 3) and comparisons of the in-frared line luminosities with observations (Fig. 7), for thesake of brevity. In this appendix we present these figures forthe full sample of all emission lines considered in this study.
MNRAS000
MNRAS000 , 1–17 (2021) A. J. Richings, C.-A. Faucher-Gigu`ere and J. Stern
Figure A1.
Continuum-subtracted spectra of the millimetre andinfrared emission lines (top panels) and the optical and UV lines(bottom panels) in the full sample, from simulations L45 (dashedcurves) and L46 (solid curves). The grey shaded bands indicatevelocity ranges that are excluded from further analysis, to focuson the outflow component.
Figure A2.
Maps of the velocity-integrated continuum-subtracted emission from all millimetre, infrared, optical and UVemission lines in the full sample from simulation L45.MNRAS , 1–17 (2021) ultiphase AGN winds Figure A3.
As Fig. A2, but for the L46 simulation.
Figure A4.
Line luminosity ( L line ) divided by bolometric AGNluminosity ( L AGN ) versus L AGN from the simulations (RFG18;grey symbols) and observed samples of AGN host galaxies fromCicone et al. (2014) (C14; stars), Veilleux et al. (2009) (V09;squares), Herrera-Camus et al. (2018) (HC18; diamonds), Landtet al. (2008) (L08; left triangles) and Lamperti et al. (2017) (L17;right triangles). The observations are divided according to AGNclassification, as indicated by the colour. Upper limits are denotedby arrows.MNRAS000