Unravelling the physics of multiphase AGN winds through emission line tracers

Observations of emission lines in Active Galactic Nuclei (AGN) often find fast (~1000 km s^-1) outflows extending to kiloparsec scales, seen in ionised, neutral atomic and molecular gas. In this work we present radiative transfer calculations of emission lines in hydrodynamic simulations of AGN outflows driven by a hot wind bubble, including non-equilibrium chemistry, to explore how these lines trace the physical properties of the multiphase outflow. We find that the hot bubble compresses the line-emitting gas, resulting in higher pressures than in the ambient ISM or that would be produced by the AGN radiation pressure. This implies that observed emission line ratios such as OIV 25 ${\mu}$m / NeII 12 ${\mu}$m , NeV 14 ${\mu}$m / NeII 12 ${\mu}$m and NIII 57 ${\mu}$m / NII 122 ${\mu}$m constrain the presence of the bubble and hence the outflow driving mechanism. However, the line-emitting gas is under-pressurised compared to the hot bubble itself, and much of the line emission arises from gas that is out of pressure, thermal and/or chemical equilibrium. Our results thus suggest that assuming equilibrium conditions, as commonly done in AGN line emission models, is not justified if a hot wind bubble is present. We also find that>50 per cent of the mass outflow rate, momentum flux and kinetic energy flux of the outflow are traced by lines such as NII 122 ${\mu}$m and NeIII 15 ${\mu}$m (produced in the 10^4 K phase) and CII 158 ${\mu}$m (produced in the transition from 10^4 K to 100 K).


INTRODUCTION
Galaxies that host active galactic nuclei (AGN) have been observed to contain fast (∼1000 km s −1 ) outflows of gas on kiloparsec scales (e.g. Cecil 1988;Veilleux et al. 2003), likely driven by the input of energy and momentum from accretion onto the central supermassive black hole (e.g. Feruglio et al. 2010;Cicone et al. 2014).
Given their high velocities, one might expect the outflowing material to be almost exclusively hot ( 10 7 K). However, these outflows have been observed in a wide range of emission and absorption lines, spanning molecular, neutral atomic and ionised gas phases (Fischer et al. 2010;Rupke & Veilleux 2011;Greene et al. 2012;Harrison et al. 2012Harrison et al. , 2014Aalto et al. 2015;Fiore et al. 2017;Feruglio et al. 2020).
Several possible mechanisms for the formation of the cool phase in outflows have been considered (applicable to both AGN-and star formation-driven winds), including the Email: alexander.j.richings@durham.ac.uk entrainment of cold gas clouds in a hot outflow (Scannapieco & Brüggen 2015;Gaspari & Sadowski 2017;Schneider & Robertson 2017;Gronke & Oh 2020), the direct acceleration of cold gas by cosmic rays or radiation pressure (Booth et al. 2013;Costa et al. 2018;Hopkins et al. 2020), or in-situ cooling of the hot outflowing material (Wang 1995;Silich et al. 2003;Zubovas & King 2014;Costa et al. 2015; Thompson et al. 2016). While the true origin of the cool phase in these powerful outflows remains uncertain, the observational evidence for such a multiphase structure is nonetheless overwhelming (see also Veilleux et al. 2020 for a recent review).
AGN-driven galactic outflows are important as they are likely to play a vital role in shaping the formation and evolution of galaxies. The energy and momentum injected by black hole winds can regulate the growth of the black holes and the stellar component of their host galaxies, giving rise to the observed scaling relations between the two (Silk & Rees 1998;King 2003;Murray et al. 2005;Zubovas & King 2012;Torrey et al. 2020). Energetic feedback from AGN can quench star formation in the most massive galaxies (Binney & Tabor 1995;Dubois et al. 2013;Bower et al. 2017), and is required by modern cosmological models of galaxy formation to reproduce observed 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é et al. 2019), the colours of massive elliptical galaxies (Springel et al. 2005;Feldmann et al. 2017;Trayford et al. 2016), and the stellar mass densities of massive galaxies (Choi et al. 2018;Wellons et al. 2020;Parsotan et al. 2020). Outflows can also enrich the circumgalactic medium (CGM) through the transport of metals from the galaxy (Hummels et al. 2013;Ford et al. 2013;Muratov et al. 2017;Tumlinson et al. 2017;Hafen et al. 2019), and can deplete CGM gas fractions by carrying baryons beyond the virial radius Oppenheimer et al. 2020).
To understand how such outflows can influence the surrounding environment, we need to measure their energetics (e.g. mass outflow rates, momentum fluxes, and energy fluxes), to determine whether they contain sufficient energy and momentum to have a significant impact on their host galaxy. This requires us to connect the emission and absorption line tracers in which the outflows have been observed to the physical properties of their constituent gas components.
Many studies have successfully employed photoionisation models to reproduce the emission line properties of AGN. Dopita et al. (2002) and Groves et al. (2004a) presented photoionisation models of dusty clouds dominated by radiation pressure in the Narrow Line Region (NLR) of AGN. The pressure and density structure of the clouds were determined by assuming that they are in hydrostatic balance with the radiation pressure exerted by the AGN, also known as Radiation Pressure Confinement (RPC; see also Draine 2011 andYeh &Matzner 2012 for the same effect in ionised gas around star forming regions). They then used the mappings iii photoionisation and shock code to calculate the intensities of emission lines on a grid of density, metallicity, ionisation parameter and the power-law index of the ionising spectrum. Observed line ratio diagnostics can then be compared to these models to deduce the physical properties of 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 a large range of distances from the nucleus. They showed that RPC can explain the observed optical emission line ratios and the overlap of extended X-ray and optical line emission in nearby Seyferts. Bianchi et al. (2019) later showed that these types of models can explain the emission measure distribution of X-ray emission lines. RPC models also explain observed Broad Line Region (BLR) emission line ratios over a range of 10 8 in AGN luminosity (Baskin et al. 2014) and the broad ionisation distribution in AGN outflows (Stern et al. 2014b). Stern et al. (2016) used these RPC models and observed emission line ratios in luminous quasars to constrain the relative importance of hot gas and radiation pressure on the dynamics of quasar outflows. They found no evidence for the compression of emission line gas expected in the presence of a hot bubble on any scale, and argued that a dynamically important hot bubble can be ruled out on spatial scales below 40 pc.
Such photoionisation models have had success in reproducing many of the observed emission line properties of AGN. However, they assume pressure, thermal and/or chemical equilibrium in the line-emitting gas. These assumptions need to be tested further.
In Richings & Faucher-Giguère (2018a) (hereafter RFG18), we ran a suite of hydro-chemical simulations of AGN winds to explore the origin of molecular outflows. We demonstrated that observed molecular outflows in AGN host galaxies can be produced by the in-situ formation of new molecules within the outflowing material. In that work, we focussed on the molecular phase. However, the chemical modelling in these simulations also includes the chemistry of the atomic and ionised phases. This enables us to make direct predictions for the emission lines from all phases of the outflow. Most previous hydrodynamic simulations of AGN outflows do not directly model the line emission, which is important as emission lines contain most of the observational constraints on the models. Our simulation suite therefore presents a unique tool with which to study the connection between the physical properties and energetics of multiphase AGN outflows and the observable emission line tracers, and to test the underlying assumptions that enter into alternative photoionisation models for AGN emission.
In this paper, we compute emission lines from the ionised, neutral atomic and molecular phases of the AGN outflows in the RFG18 simulations, which we use to explore the physical properties of the line-emitting gas. The remainder of this paper is organised as follows. In Section 2 we describe the simulations and the radiative transfer calculations used to model the line emission. In Section 3 we present predictions for the line emission from the simulations (Section 3.1), investigate the physical properties of the line-emitting gas (Section 3.2), and compare our model predictions to observations (Section 3.3). We explore how the emission line ratios computed from our simulations can be used to constrain the driving mechanisms of AGN outflows in Section 4, and in Section 5 we study the outflow energetics of the different gas phases traced by the various emission lines. Finally, we summarise our results in Section 6.

Simulations
In RFG18 we presented a series of hydro-chemical simulations of AGN-driven galactic outflows. These simulations model an initially uniform ambient medium, into which we inject an isotropic AGN wind by spawning gas particles in the central parsec with an outward velocity of 30 000 km s −1 and a momentum injection rate determined by the AGN luminosity, LAGN. Each simulation follows the interaction of this wind with the ambient medium over 1 Myr, which corresponds to the typical flow times (r/v) of kiloparsecscale outflows that have been observed in luminous quasars (e.g. González-Alfonso et al. 2017). These simulations use the gravity+hydrodynamics code gizmo, with the Meshless Finite Mass (MFM) hydro solver (Hopkins 2015) and a fiducial resolution of 30 M per gas particle.
We model the time-dependent chemistry of the gas using the chimes non-equilibrium chemistry and cooling module 1 (Richings et al. 2014a,b), which follows the evolution of 157 ions and molecules that dominate the cooling rate from cold (∼10 K), molecular gas to hot (>10 9 K), highly ionised plasmas. The chimes chemical network contains various collisional reactions, including collisional ionisation, recombination, charge transfer reactions, and reactions on the surface of dust grains such as the formation of molecular hydrogen (assuming a constant dust-to-metals ratio), along with photoionisation and photodissociation reactions.
For the photochemistry, we use the average quasar spectrum from Sazonov et al. (2004) normalised according to the bolometric AGN luminosity and the distance of the gas particle from the black hole. This is an average between an obscured and an unobscured spectrum, and is characteristic of a typical quasar. This spectrum differs from that used in some other models of AGN emission lines. For example, Stern et al. (2016) use a power-law spectrum representative of an unobscured quasar, with a fiducial extreme UV (EUV) slope of −1.6 (along with variations in this slope). In particular, the UV to X-ray ratio is a factor ≈2 lower in the average Sazonov et al. (2004) spectrum than the fiducial power-law spectrum of Stern et al. (2016). We will consider the effects of the choice of spectrum below.
Self shielding is included using a Sobolev-like shielding length approximation based on the density gradient (see equation 5 of RFG18). The densities of individual species are multiplied by this shielding length to obtain their column densities, which are then used to attenuate the photoionisation, photodissociation and photoheating rates.
Further details of our simulations can be found in RFG18. For this paper, we focus on the low-luminosity (nH10 L45 Z1) and fiducial (nH10 L46 Z1) runs. Both use a uniform ambient medium with an initial hydrogen density nH = 10 cm −3 and solar metallicity 2 , but they differ in the bolometric AGN luminosity, which is 10 45 and 10 46 erg s −1 , respectively. For brevity, we will refer to these two runs as L45 and L46 for the remainder of this paper. The low-density run (nH = 1 cm −3 ) from RFG18 did not cool and form a multiphase wind before the end of the simulation after 1 Myr, while the low-metallicity run (Z = 0.1Z ) strongly under-predicts the molecular outflow rates compared to observations.

Radiative transfer calculations
To create maps of the emission lines from our simulations, we post-process the simulation snapshots using version 0.40 of the publicly available Monte Carlo radiative transfer code radmc-3d 3 (Dullemond et al. 2012), using the nonequilibrium ion and molecule abundances that were calculated during the simulations with the chimes chemistry module. While the full 3D spherical outflow is included in the simulations, only one octant of the simulation volume uses the highest resolution level, with the remainder of the volume using 8× lower mass resolution. We therefore only use the high-resolution octant to produce the emission line maps. This octant is mirrored in the line of sight direction, 2 Throughout this paper, we use the solar abundances listed in table 1 of Wiersma et al. (2009), for which the solar metallicity is Z = 0.0129. 3 http://www.ita.uni-heidelberg.de/~dullemond/software/ radmc-3d/ to capture both the receding and approaching sides of the outflow, as viewed by the observer. The resulting emission maps thus cover one quadrant of the outflow.
As radmc-3d is a grid-based code, we first project the gas particles from the simulations on to an Adaptive Mesh Refinement (AMR) grid, which is refined such that no cell contains more than 8 particles. The particles are smoothed using a cubic spline kernel with a smoothing length enclosing 32 neighbours, as used in the MFM hydro solver for these simulations. When projecting the gas temperatures and velocities on to the grid, we weight the contribution from each particle by the given ion or molecule abundance, to avoid unphysical effects from mixing particles with very different properties in the same cell (see the discussion in section 5 of RFG18).
We calculate the emissivities of Hα and Hβ including both recombination of Hii and collisional excitation of Hi, using the cascade matrix formalism described in Raga et al. (2015), with the atomic data and fits for the collision strengths and recombination coefficients that they present in their Appendix A. For all other ions and molecules, the level populations are computed by radmc-3d using an approximate non-LTE approach based on the Local Velocity Gradient (LVG) method. We utilised atomic data (energy levels, transition probabilities and collisional excitation rates) from version 7.1 of the chianti database 4 (Dere et al. 1997;Landi et al. 2013), supplemented with additional collisional excitation rate data from the lamda database 5 (Schöier et al. 2005) for Ci, Cii and Oi.
Dust absorption has little to no effect on the infrared emission lines, but those at optical and UV wavelengths are strongly affected by dust. The idealised setup of these simulations, in which the outflow is embedded within a dense ambient interstellar medium (ISM), is representative of an AGN in the buried QSO phase, before the outflow has broken out of the disk of its host galaxy. In observations of buried QSOs, the outflows are very difficult to detect in optical and UV lines due to the strong dust absorption from the surrounding ISM. Even in systems where the outflows have broken out of the dense medium, they are typically only seen in the blue wing of the emission line (i.e. from the near side of the outflow, moving towards the observer).
When we include the effects of dust in our radiative transfer calculations, we also find that the optical and UV lines are strongly suppressed, with the outflow only seen in the blue wings. However, in our idealised setup the host galaxy is represented by a box of initially uniform density gas up to 2.4 kpc across, which does not include the turbulent structures we would expect to see in the ISM nor the geometry of a galactic disc. We therefore cannot model the realistic dust attenuation from the host galaxy in this setup. Including dust absorption from the ambient medium also limits how much of the optical and UV emission we can study from the simulations.
For the radiative transfer calculations in this paper, we therefore include dust grains only in the outflow. For gas with a radial velocity 0 km s −1 (i.e. the ambient ISM, which is slowly inflowing by the end of the simulation due to the gravitational potential of the host galaxy), we set the dust density to zero. We also disable dust when the gas temperature is above 10 6 K, as grains will be rapidly destroyed by sputtering in this regime (e.g. Tsai & Mathews 1995). This approach has little impact on the infrared lines. However, we stress that the resulting optical and UV line luminosities do not include the effects of dust attenuation from the host galaxy. When comparing these lines to observations, it is therefore important that we only consider line ratios at similar wavelengths, for which the strength of any additional dust absorption on each line will be similar, thus leaving the ratios unaffected. By excluding dust from the ambient medium in this way, we can focus on the nature of the gas that produces each emission line.
For the dust in the outflow, we use a mixture of graphite and silicate grains with dust-to-gas ratios of 2.4 × 10 −3 and 4.0×10 −3 , respectively, at solar metallicity. We include dust absorption, scattering and thermal emission. To obtain the continuum emission, we repeat each radmc-3d calculation a second time with the emission lines disabled. The resulting continuum emission is then subtracted from the total emission. Fig. 1 shows the spectra of six infrared emission lines (top two rows) and three optical lines (bottom row) that are commonly observed in AGN host galaxies. These represent a subset of the lines that we study in this paper; the full sample of 24 spectra can be found in Appendix A. The dashed and solid curves in each panel are from the simulations L45 and L46, respectively. Each spectrum is only integrated over the spatial extent of the outflow, out to a radius of 0.57 kpc (L45) and 1.04 kpc (L46). The fluxes are normalised to the maximum flux in each spectrum.

Spectra and images
All of these emission lines show a broad component from the outflowing shell, with velocities up to ≈300 and ≈700 km s −1 in L45 and L46, respectively. However, several of these lines show a strong narrow component arising from the ambient medium, which is seen either in emission (e.g. Cii158 µm) or absorption (e.g. Oi63 µm). For the infrared lines, we therefore only consider emission in the wings of each line for the rest of our analysis, with velocities |∆v|>100 km s −1 , to focus on the outflow component.
In the optical lines (and UV lines; not shown), the outflow is seen more strongly in the blue wing (approaching) than the red wing (receding) due to strong dust absorption at these wavelengths, as noted in Section 2.2. We again stress that we only include dust grains in the outflow for these radiative transfer calculations. This effect would be even more dramatic when we include the additional dust absorption from the ambient ISM in the host galaxy. For our further analysis of the optical and UV lines, we therefore only consider the blue wing, with velocities ∆v<−100 km s −1 . As the Sii lines are a doublet, we also need to exclude velocities <−450 km s −1 to avoid contamination of the 6731Å line by the 6716Å line. We apply this additional velocity cut to both of the sulphur lines. The velocity ranges that are excluded from the further analysis of each line are highlighted by the grey shaded bands in Fig. 1. Figure 1. Continuum-subtracted spectra of six infrared emission lines (top two rows) and three optical emission lines (bottom row) commonly observed in AGN, from the simulations L45 (dashed curves) and L46 (solid curves), normalised to the maximum flux in each spectrum. These represent a subset of the lines that we study in this paper; the spectra from the full sample of 24 lines can be found in Appendix A. The vertical dotted line in each panel indicates the line centre. The broad component of these lines arises from the outflowing shell, with velocities extending up to ≈300 and ≈700 km s −1 in L45 and L46, respectively. However, central velocities at |∆v| 100 km s −1 are typically dominated by narrow components from the ambient ISM. The grey shaded bands indicate velocity ranges that are excluded for the remainder of our analysis, to focus on the outflow component.
We can use these radiative transfer calculations to compare the spatial morphologies of each emission line from the outflow. In Figs. 2 and 3 we show velocity-integrated maps of the infrared and optical emission lines, where we include emission from the red and blue wings (infrared) or the blue wing only (optical) as discussed above. The panels are each 0.8 kpc and 1.2 kpc across in L45 and L46, respectively. We again only show a subset of the lines studied in this paper; the full sample can be found in Appendix A.
The H2, 17 µm line (top left panel of Figs. 2 and 3), which arises from the pure rotational S(1) transition of the H2 molecule, traces dense, clumpy structures that have condensed from the cooling material within the outflowing shell. Comparing to the maps of CO emission (see Appendix A), we find that the H2, 17 µm line traces the same molecular structures as the CO emission.
In the top middle panel, the infrared Oi63 µm line traces the same structures as the H2, 17 µm emission. In contrast, the Oi 6300Å optical line is somewhat more diffuse, even though both are produced by neutral atomic oxygen. We will see in Section 3.2 that the optical emission is only present in the warmer gas phase (∼10 4 K). This is unsurprising, as the 6300Å line requires a collisional excitation   energy E/kB∼10 4 K, while the 63 µm line requires 100× less energy and can thus be excited in colder gas.
The Cii158 µm emission (top right) is strong in these clumpy structures, but there is also a significant contribution in the more diffuse regions between the dense clumps. The intermediate ions (e.g. Nii122 µm and Neiii15 µm; lefthand and centre panels of the middle row) predominantly arise from the diffuse, volume-filling phase in the outflowing shell, while the highest ionisation states such as Nevi7 µm (right-hand panel of the middle row) are dominated by small, bright knots of emission. We will show in the next section that these high ionisation states are produced in gas that is undergoing a period of rapid cooling. As such, any one particular patch of gas spends a very short time producing these high ionisation lines. The small, bright knots that we see in Nevi7 µm are therefore short flashes of emission as these compact regions transition through a rapid cooling phase.

Physical properties of the line-emitting gas
In this section we explore the physical properties of the gas probed by each emission line. Before we break down the gas properties by the separate emission line tracers, it is useful to look at the properties of the total gas distribution in the simulations. Fig. 4 shows the temperature versus density of all gas in the high-resolution octant of the L46 simulation, where the colour scale indicates the distribution of total gas mass in each pixel. We have highlighted several regions in this plot to guide the discussion. We only show L46 here, but the L45 simulation exhibits the same structures in temperature-density space (see also Fig. 5 of RFG18). The evolution of gas through the temperature-density plane proceeds as follows: 1) The ultra-fast outflow (30 000 km s −1 ) injected in the central parsec encounters the reverse shock (also called the wind shock) and is heated to ∼10 10 K, creating a hot bubble. The pressure of this hot bubble is P hot /2.3kB = nHT = 2.0 × 10 8 cm −3 K. This pressure accelerates the outflow, boosting the momentum beyond that expected in a momentum-driven outflow (Faucher-Giguère & Quataert 2012). 2) The ambient ISM is intially in thermal equilibrium at a density nH = 10 cm −3 before it has been hit by the outflow. 3) When the ambient ISM is swept up by the forward shock driven by the outflow, it is compressed and shock heated to ∼10 6.5−7 K. This corresponds to the post-shock temperature for a shock velocity of ≈500 − 1000 km s −1 , according to the Rankine-Hugoniot jump conditions for a strong shock, which is consistent with the 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 with the hot bubble, but remains in pressure equilibrium. 5) The outflowing shell of gas swept up from the ambient medium initially remains close to the post-shock temperature. 6) Once the swept up gas reaches the peak in the cooling curve, the cooling time becomes very short and it rapidly cools. We can estimate the cooling time in this regime based on the analytic models presented in Richings & Faucher-Giguère (2018b). In Fig. 6 of that work, we presented the temperature evolution of the swept up shell of an isotropic AGN wind in a 1D analytic model based on Faucher-Giguère & Quataert (2012). The fiducial model (black curves) uses the same AGN luminosity, ambient ISM density and metallicity as our L46 simulation. In this model, the swept up shell cools from 10 6 K to 10 4 K in a time t cool = 0.002 Myr. For comparison, the sound crossing time of the shell, with a thickness of ≈25 pc (from the simulations) and using the sound speed at 10 6 K, is t sound = 0.2 Myr. Hence the cooling time is 100× less than the sound crossing time in this regime, and so it cools isochorically, as the pressure of the surrounding hot medium has insufficient time to compress the cooling gas and maintain pressure equilibrium during this phase. 7) Gas that has just undergone rapid cooling and is approaching thermal equilibrium is now under-pressurised compared to the surrounding hot medium. 8) Once the gas reaches thermal equilibrium at a temperature ∼10 4 K, it is subsequently compressed, increasing in density and pressure. 9) At higher densities, the gas undergoes a thermal instability and cools down to the cold phase at temperatures ∼100 K. The cooling in this regime proceeds isobarically, as the 10 4 K and 100 K phases remain approximately in pressure equilibrium with one another, albeit with a large scatter. As we demonstrated in RFG18, this cold phase is conducive to the formation of new molecules.
Figs. 5 and 6 show the temperature-density distribution of the line-emitting gas for the millimetre/infrared lines and the optical/UV lines, respectively. In these figures we show the full sample of 24 emission lines, not just a subset as in previous figures. The colour scale indicates the fraction of the total line luminosity that is produced from gas in each 2D temperature-density bin, normalised by the size of the bin in logarithmic temperature-density space. Plots showing the temperature and density of gas producing the millimetre and infrared lines from the L46 simulation, arranged in order of increasing ionisation energy (or dissociation energy for molecules). The colour scale indicates the fraction of the total line luminosity that originates from gas in each 2D temperature-density bin. The dashed, dot-dashed and dotted 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 the radius of the outflowing shell (1.04 kpc), and the pressure of the ambient ISM (P ISM ). The grey contours enclose 90 per cent of the total emission. The top axis shows the ionisation parameter, evaluated at 1.04 kpc. The hot bubble has compressed the lineemitting gas, resulting in higher pressures than in the ambient ISM or that would be produced from radiation pressure alone. However, this gas is under-pressurised compared to the hot bubble itself, especially for the highest ionisation states, which are located in gas that has recently undergone a rapid cooling phase.
These are calculated using the individual line luminosities and ion/molecule-weighted temperatures and densities from each cell in the AMR grid used for the radiative transfer calculations (see Section 2.2). We only show the L46 simulation here, but L45 exhibits the same trends. The grey contours in each panel enclose 90 per cent of the total line emission.
In Fig. 5, we see in the top row that the CO2.6 mm emission (from the J=1−0 transition) and the H2, 17 µm emis- Like the millimetre and infrared lines, emission at optical and UV wavelengths also arises from gas that is over-pressurised compared to the ambient ISM and the radiation pressure alone, but underpressurised compared to the hot phase.
sion (from the pure rotational S(1) transition) trace cold gas, at temperatures 10 3 K. In contrast, the H2, 2.2 µm line, from the v1−0 S(0) ro-vibrational transition, traces warmer molecular gas up to several thousand K. This is unsurprising given that the vibrational energy levels require a higher collisional excitation energy than the purely rotational transitions, and has also been noted many times elsewhere in the literature (e.g. Veilleux et al. 2020). In observations of the Seyfert galaxy NGC 2110, Rosario et al. (2019) demonstrate that the warm molecular gas traced by the ro-vibrational v1−0 S(1) transition is spatially anti-correlated with the CO J=2−1 transition. This further supports the picture that the vibrationally excited H2 and the CO trace different phases of molecular gas.
The infrared lines of the low-ionisation states (e.g. Oi and Cii) arise from gas that is cooling from the warm (10 4 K) to the cold (100 K) phase, corresponding to region (9) in Fig. 4. Intermediate ions such as Nii and Oiii originate from gas that lies on the thermal equilibrium branch at 10 4 K, i.e. region (8) in Fig. 4. Finally, the highest ionisation states such as Nev and Nevi are produced by gas that is coming to the end of the rapid cooling phase and is approaching the thermal equilibrium temperature.
In Fig. 6 we have separated the hydrogen Balmer line emission between that produced by the recombination of Hii ('rec'; left panels of the top two rows) and that driven by the collisional excitation of Hi ('col'; central panels of the top two rows). The recombination radiation arises from gas undergoing the rapid cooling phase, similarly to the highionisation states. In contrast, collisional excitation radiation arises only from neutral gas on the 10 4 K branch. Raga et al. (2015) demonstrated that, at 10 5 K, the emission coefficients of Hα and Hβ from collisional excitation are up to five orders of magnitude higher than those due to recombination. This can lead to significant contributions from collisional excitations at these high temperatures even for fairly low Hi fractions. The lack of such collisional excitation radiation at 10 5 K in our simulations indicates that the Hi fractions in the rapid cooling phase are very low, <10 −5 . In L46, collisional excitation contributes 54 and 19 per cent to the total Hα and Hβ emission, respectively. In L45 (not shown), it contributes 43 and 10 per cent to Hα and Hβ, respectively.
The 6300Å emission from Oi arises predominantly from the 10 4 K phase. This contrasts with the infrared lines from Oi, 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, as the 6300Å line requires a higher excitation energy. The optical and UV emission from Sii, Nii, Oiii and Neiii also trace the 10 4 K phase, while the Nev UV emission includes gas towards the end of the rapid cooling phase.
To understand the role of photoionisation in the lineemitting gas, it is instructive to look at the ionisation parameter, Uion. This is defined as the ratio of the densities of hydrogen-ionising photons and hydrogen nuclei (nγ and nH, respectively), and can be calculated as follows: where L 13.6 eV−1 keV ion is the ionising luminosity from 13.6 eV to 1 keV, r is the distance from the AGN, c is the speed of light, and hν 13.6 eV−1 keV is the mean energy of ionising photons in the 13.6 eV − 1 keV band. We include ionising photons only up to 1 keV here, but the ionisation parameter would be almost unaffected if we instead included all photons >13.6 eV, as there are relatively few photons in the X-ray band. For the average quasar spectrum from Sazonov et al. (2004) that we use in our simulations, hν 13.6 eV−1 keV = 32.3 eV, and L 13.6 eV−1 keV ion /LAGN = 0.10, where LAGN 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 r = 1.04 kpc (in L46) from the AGN. The ionising flux, and hence the density of ionising photons, is therefore approximately uniform for this gas, and thus Uion depends simply on the inverse of the gas density. We therefore show Uion at 1.04 kpc on the top axis of each panel in Figs. 5 and 6. In general, emission from the highest ionisation states arises from gas with 10 −2 Uion 10 −1 , while intermediate ions have 10 −4 Uion 10 −2 . Gas traced by molecules and low ionisation states extends to even lower ionisation parameters, down to 10 −5 .
The dashed, dotted and dot-dashed lines in Figs. 5 and 6 indicate the hot gas pressure (P hot , defined from region 1 in Fig. 4), the initial pressure of the ambient ISM (PISM) and the radiation pressure from the AGN at a radius of 1.04 kpc corresponding to the outflowing shell (P rad ), respectively.
The radiation pressure can be calculated from the ionising luminosity as follows: where we again include ionising radiation only up to 1 keV, as the X-rays will penetrate beyond the ionised layer and therefore 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 photons to the radiation pressure. They show that β≈1 in dust-free gas, and β≈2 in dusty gas. We showed in RFG18 that we require dust to be present in the outflowing shell to produce molecules consistent with observations. We therefore take β = 2 when calculating the radiation pressure. Note that, while photoionisation and photodissociation from the AGN radiation field are included in the thermochemistry in these simulations, direct radiation pressure itself is not included. This would tend to compress gas at thermal pressures less than P rad (e.g. Dopita et al. 2002;Draine 2011). We see in Figs. 5 and 6 that some of the emission lines arise from gas at pressures below P rad . We would therefore expect this gas to undergo additional compression if we were to include radiation pressure in the simulations, but the vast majority of emission in our simulations arises from gas with Pgas P rad , and thus we do not expect the inclusion of radiation pressure in our simulations to significantly affect our results.
We see that the hot bubble has compressed the line emitting gas to become over-pressurised compared to the ISM. For most emission lines, it also reaches higher pressures than would be achieved through direct radiation pressure alone. However, as this gas has cooled below the initial post-shock temperature of the outflowing shell, it is no longer in pressure equilibrium with the hot phase. The high ionisation species (e.g. Nevi) trace gas that is most strongly under-pressurised compared to the hot medium, as it has not had time to be compressed and so is still at relatively low densities where higher ionisation states are prevalent. In contrast, the low-and intermediate-ionisation states trace gas at higher densities that have had longer to be compressed and are therefore closer to the pressure of the hot bubble, although they still remain under-pressurised.
Photoionisation models of AGN emission lines such as Groves et al. (2004a) and Stern et al. (2014a) also assume thermal equilibrium in the line-emitting gas. Our simulations support this assumption for the intermediate lines such as Nii122 µm, Nii 6585Å , Sii 6716, 6731Å and Oi 6300Å . In Figs. 5 and 6, these lines trace a narrow region at the thermal equilibrium temperature ∼10 4 K. However, high ionisation states such as Nev and Nevi arise from a broad range of temperatures at a given density (spanning up to ∼0.5−1 dex) as they are coming to the end of the rapid cooling phase, and they have not yet reached thermal equilibrium. Emission from molecules and low-ionisation states such as Oi63, 145µm and Cii158 µm that trace the 100−10 4 K phase also cover a broad range of temperatures at each density, rather than simply tracking the thermal equilibrium temperature as a function of density. This is likely due to turbulent shock heating. We find that, in our simulations, the dense (nH>100 cm −3 ) gas clouds that condense within the outflowing shell after it cools typically have velocity dispersions up to a few tens of km s −1 . The shocks resulting from this turbulence can heat the gas up to a few thousand K.
The photoionisation models also assume that the ions and molecules are in chemical equilibrium. To test this assumption in our simulations, we computed the equilibrium abundances of each gas particle at fixed density and temperature, and then repeated the radiative transfer calculations using the equilibrium abundances. The most notable difference is in the Oiii 5007Å line, which is ≈7−8 times higher when we set the abundances to equilibrium than when we use non-equilibrium abundances. Other emission lines also show modest deviations when we assume equilibrium abundances, in particular Nev14 µm (up to 80 per cent), Nii122 µm (up to 62 per cent), Oiv25 µm (up to 60 per cent), Niii57 µm (up to 44 per cent), Siii18 µm (up to 39 per cent), Nii 6585Å (up to 35 per cent), H2; 2.2 µm (up to 31 per cent), H2; 17 µm (up to 29 per cent), and CO2.6 mm (up to 26 per cent). The remaining lines differ by <20 per cent when we use equilibrium abundances.

Comparison with observations
We now compare the predictions from our simulations to observations of outflows in AGN host galaxies. Fig. 7 compares the line luminosity (L line ) divided by the bolometric AGN luminosity (LAGN) plotted against LAGN for a subset of the infrared lines (the full sample can be found in Appendix A). We divide by LAGN on the y-axis to reduce the dynamic range; a simple linear scaling between L line and LAGN would thus manifest itself as a horizontal trend in these plots. For the simulations (RFG18; grey symbols), our radiative transfer calculations include only one quadrant of the spherical outflow (see Section 2.2). We therefore multiply the resulting luminosities by four to obtain the luminosity of the whole outflow.
In the left-hand column of Fig. 7 we compare mid-IR lines from our simulations to a sample of observed AGN host galaxies from Veilleux et al. (2009) (V09; squares). We divide the V09 sample according to the AGN classification as a Palomar-Green Quasar (PG QSO), UltraLuminous In-fraRed Galaxy Seyfert 1 (ULIRG Sy1), ULIRG Sy2, or an ULIRG Low Ionisation Nuclear Emission Regions (LINER). We exclude galaxies classified as Hii-like. The different classifications are denoted by the colour. In the right-hand column we compare far-IR lines from our simulations to the observations from Herrera-Camus et al. (2018) (HC18; diamonds), divided as a Sy1, Sy2, LIRG Sy2 and LIRG LINER, again indicated by the colour. We exclude those classified as starburst galaxies in HC18. For the V09 sample, we use the bolometric luminosities listed in their work. However, this information is not included in HC18. For this latter sample we therefore use the bolometric luminosity given in V09 where available, otherwise we derive it from the far-IR (40 µm−120 µm) luminosity using the bolometric correction from Table 10 of V09 and assuming a 100 per cent AGN fraction, i.e. LAGN = 1.22LFIR. For many of the lines shown in Fig. 7, the two simulations overlap with the observations. A notable exception is the Oi63 µm line. While there are two LIRG LINER systems that are close to, or higher than, the Oi63 µm luminosity predicted from L45, most of the observations lie at least an order of magnitude below the simulations. Comparing to the individual AGN classifications, we see that in the H2, 17 µm and Neiii15 µm lines the L45 simulation is closest to the ULIRG LINER systems, while L46 most closely overlaps with the ULIRG Sy1/Sy2 galaxies. In Nevi7 µm L46 is again closest to the ULIRG Sy1/Sy2, while L45 is closer to the PG QSOs. In the far-IR lines Cii158 µm and Nii122 µm, it is more difficult to associate the simulations with a particular AGN classification. The observed LIRG LINERs bracket the line luminosities predicted by the simulations. However, in the observations these systems show a trend of decreasing L line /LAGN with increasing LAGN, which is not replicated in the simulations. Most of the observed Sy1 and Sy2 galaxies have lower (higher) luminosities of Cii158 µm (Nii122 µm) than the simulations, typically by a fac- While many of the emission line luminosities predicted from the simulations overlap with the observations in Fig. 7, the large dynamic covered by the observations limits the constraining power in comparing the luminosities of individual lines. However, emission line ratios are often confined to a narrower range in observed systems, which provides tighter constraints on the models. In Fig. 8 we therefore compare infrared line ratios from the simulations (grey symbols) and the observed samples from V09 (squares) and HC18 (diamonds). In the top-left panel of Fig. 8, the simulations lie on the observed correlation between these two line ratios. Veilleux et al. (2009) also demonstrated this correlation in the observations, which represents an excitation sequence extending from the ULIRG LINERs in the bottom left, through the ULIRG Sy2 and Sy1 galaxies and on to the PG QSOs in the top right. They note that this trend can be understood if Nev14 µm and Oiv25 µm are produced by the AGN while Neii12 µm comes from starburst activity. For fixed Nev14 µm / Oiv25 µm, varying the relative contribution of the starburst will move the ratios along the observed correlation, with the PG QSOs exhibiting the smallest contribution from the starburst, while the ULIRG LINERs contain the highest starburst contribution. Our simulations lie closest to the observed ULIRG Sy2 systems, although we do not include contributions from a starburst component, Figure 9. BPT diagrams from our simulations (black symbols), star-forming galaxies and type 2 AGN from SDSS (Kewley et al. 2006;grey 2D histogram), and type 1 AGN from SDSS with bolometric luminosities >10 45 erg s −1 (Stern & Laor 2013; blue contours). The black solid and dotted curves show the boundaries between star-forming and active galaxies from Kewley et al. (2001) and Kauffmann et al. (2003), respectively, while the black dashed lines show the boundary between LINERs and Seyferts from Kewley et al. (2006). The Sii/Hα and Oi/Hα line ratios are up to an order of magnitude higher in the simulations than is seen in the observations, due to strong X-ray heating in the assumed partially-obscured incident quasar spectrum and unresolved ionised layers.
which would tend to move the simulations down and to the left in this plot.
In the two panels showing ratios including Nevi (topright and bottom-left), the simulations lie above the observations by ≈0.5 dex. This suggests that the Nevi7 µm line may be too strong in our simulations. The predicted Nevi7 µm / Neii12 µm and Nevi7 µm / Neiii15 µm ratios from the simulations in these two panels overlap with the observations, however the observations cover almost two orders of magnitude and so would still be consistent with reducing the Nevi7 µm emission in the simulations by ≈0.5 dex.
In the bottom right panel of Fig. 8, the simulations again deviate from the observations. While the Niii57 µm / Nii122 µm ratios from the simulations and observations overlap, Oi63 µm / Cii158 µm is an order of magnitude too high in the simulations. As we saw in fig. 7, this is due to the Oi 63 µm line, which is 10× too high in the simulations.
As noted in Section 2.2, our radiative transfer calculations do not include dust attenuation from the ambient ISM of the host galaxy. While this has little effect on the infrared lines, the resulting luminosities of the optical and UV lines are much stronger than we would expect were such attenuation to be included. Nevertheless, we can compare the 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) of optical emission line ratios, which are commonly used to distinguish star-forming and AGN-dominated galaxies. The panels show the ratios of Nii 6585Å / Hα 6563Å (top left), Sii 6716+6731Å / Hα 6563Å (top right) and Oi 6300Å / Hα 6563Å (bottom left) on the x-axis, plotted against Oiii 5007Å / Hβ 4861Å on the y-axis in each case. The simulations are shown by the black symbols, while the grey 2D histograms show the distribution of star-forming galaxies and type 2 AGN in the Sloan Digital Sky Survey (SDSS) from Kewley et al. (2006) (updated to SDSS-DR8). The blue contours show observed type 1 AGN in SDSS from Stern & Laor (2013), for which we only include those with a bolometric luminosity >10 45 erg s −1 to select the most powerful AGN in this sample. The black solid and dotted curves show the boundaries between star-forming and active galaxies from Kewley et al. (2001) and Kauffmann et al. (2003), respectively, while the black dashed lines show the boundary between LINERs and Seyferts from Kewley et al. (2006).
In the top left panel panel of Fig. 9, the L45 simulation overlaps with the observed AGN from SDSS and Stern & Laor (2013), although it lies close to the boundary between AGN and star-forming galaxies from Kewley et al. (2001). L46 would also be classified as an AGN from this panel, although the Oiii 5007Å / Hβ 4861Å is somewhat higher than is seen in even the most extreme AGN in the observations, by up to a factor ≈2.
In the top right and bottom left panels, the simulations lie on the LINER side of the classification boundaries. However, the Sii 6716+6731Å / Hα 6563Å and Oi 6300Å / Hα 6563Å ratios are an order of magnitude higher than the observations. The Hα emission traces recombinations driven by photoionisation, while the Sii and Oi emission traces cooling in the Warm Neutral Medium (WNM) phase at ∼10 4 K which is radiated via these metal lines. The high line ratios in the BPT diagrams predicted by the simulations therefore suggest that the ratio of the WNM to the photoionised phases is likely to be too high in our simulations, by up to an order of magnitude.
We find there are two effects that contribute to these anomalous ratios. Firstly, the WNM in our simulations is supported by X-ray heating. If we recompute the thermal equilibrium temperature without the X-ray component of the incident radiation field from the quasar, we find that much of the gas in the WNM cools from 10 4 K to a few thousand K. As noted in Section 2.1, the average incident quasar spectrum that we use from Sazonov et al. (2004) has a lower UV to X-ray ratio, by a factor ≈2, than typical unobscured quasar spectra used in other AGN emission line models that successfully reproduce the BPT diagram (e.g. Stern et al. 2016). This means there will be relatively stronger X-ray heating, compared to photoionisation driven by the UV band, which will increase the ratio of the WNM to photoionised phases and will thus contribute to the high line ratios that we find in the BPT diagrams.
Secondly, the transition from ionised to atomic hydrogen occurs at a column density NH, tot≈10 20.5 cm −2 in our simulations. However, at a density of nH=100 cm −3 and our fiducial mass resolution of mgas=30 M , this is comparable to the column density of a single gas particle, which suggests that the photoionised layer in our simulations is barely resolved. To estimate the extent to which we may be underestimating the photoionised component due to limited resolution, we can estimate the expected Hα luminosity from recombinations of Hii given the incident ionising luminosity of the quasar. For the recombination emission coefficient that we use from Raga et al. (2015), together with the rate coefficient for case B recombination used in chimes, we expect that each recombination of Hii will produce 0.27 Hα photons at 10 4 K. If we assume that every incident ionising photon is absorbed and thus leads to one recombination, the resulting Hα luminosity that we would expect given the ionising luminosity is ≈3−4 times higher than that found in the simulations (before radiative transfer effects such as dust attenuation are included). This suggests that we underestimate the photoionised component by a factor ≈3−4.
We therefore conclude that these two effects of strong Xray heating in the assumed partially obscured quasar spectrum and unresolved photoionised layers lead to the anomalously high line ratios in the BPT diagrams. Our simulations therefore do not provide reliable predictions for the BPT diagram. However, we caution that other emission line predictions presented above, in particular those comparing the photoionised and WNM phases, are also likely to be influenced by these effects. For example, the enhancement of the WNM due to strong X-ray heating will contribute to the anomalously high Oi63 µm luminosity seen in Fig. 7. Stern et al. (2016) argued that, in AGN outflows driven by the thermal pressure of a hot bubble, the line-emitting gas is compressed to higher pressures than we would expect in an outflow driven purely by the radiation pressure exerted by the AGN. They used idealised photoionisation calculations assuming hydrostatic balance between the hot phase and the line-emitting gas to show that this compression by the hot bubble leads to a dependence of the emission line ratios on the ratio of P hot / P rad . By comparing these predictions to observations of quasar outflows, they concluded that, on large scales ( 1 kpc, comparable to the scales that we probe with our simulations), models with P hot /P rad > 6 are incompatible with the observed line ratios. This limits the dynamical role that hot gas pressure can play in driving galactic outflows in quasars. This comparison considered emission lines at optical and UV wavelengths, which selects unobscured quasars for which there are clear paths to the central AGN. As noted by Stern et al. (2016), this selection might limit the effects of hot gas pressure as the hot bubble can begin to vent out and is no longer constrained by the dense ISM of the host galaxy (e.g. Torrey et al. 2020). These results therefore do not rule out an earlier buried phase in which the hot gas pressure may dominate. While our simulations do include the photoionisation and photodissociation from the AGN radiation, we do not include the direct radiation pressure. However, we saw in Figs. 5 and 6 that the radiation pressure is more than an order of magnitude less than the hot gas pressure in our simulations. We can therefore use our simulations to test the assumptions in the Stern et al. (2016) models in the P hot P rad regime. However, note that in our simulations Figure 10. Infrared line ratios versus the ratio of hot-gas to radiation pressure (P hot /P rad ) in the simulations (black symbols) and an idealised cloudy photoionisation model in hydrostatic equilibrium based on Stern et al. (2016) (blue dotted curves). The horizontal blue lines indicate the P hot P rad (i.e. RPC) limit of the cloudy model. In the top five panels, the line ratios from the simulations are close to the RPC limit. However, in the bottom three panels the simulations lie more than 1 dex below the RPC limit predictions. These three line ratios in the bottom panels therefore present the best opportunity to observationally distinguish between the hot gas pressure scenario and the RPC limit. Figure 11. The line ratios Oiii 5007Å / Hβ 4861Å (left) and Nev 3426Å / Neiii 3869Å (right) plotted against the ratio of hot gas pressure to radiation pressure (P hot /P rad ). The black symbols show our simulations, while the blue dotted curves show the line ratios in an idealised cloudy photoionisation model in hydrostatic equilibrium based on Stern et al. (2016) using the Sazonov et al. (2004) average quasar spectrum. The line ratios in the P hot P rad (i.e. RPC) limit of the cloudy model are indicated by the horizontal blue lines. The line ratios from the simulations bracket the RPC limit of the cloudy model, which means that we cannot use these particular ratios to distinguish between the hot gas pressure dominated limit and RPC. the outflow remains constrained by the dense ISM, while Stern et al. (2016) focussed on comparisons with unobscured quasars in which the hot bubble may begin to vent out.

CONSTRAINTS ON THE DRIVING MECHANISMS OF AGN OUTFLOWS
As discussed in Section 3.2, the hot bubble in our simulations does compress the line-emitting gas, as suggested in Stern et al. (2016). However, the rapid cooling phase leaves this gas under-pressurised compared to P hot . We therefore cannot assume hydrostatic balance between the line-emitting gas and the hot phase. We explore the implications that this has for the emission line ratio constraints on 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 ratios versus P hot /P rad . The black symbols show our simulations, and the blue dotted curves show an idealised cloudy photoionisation model in hydrostatic equilibrium, including dust grains, based on the Stern et al. (2016) models but using the Sazonov et al. (2004) average quasar spectrum as in the hydrodynamic simulations. The horizontal blue lines indicate the Radiation Pressure Confinement (RPC) limit of the cloudy model, i.e. for P hot P rad .
The strongest distinction between our simulations and the RPC limit can be seen in the bottom three panels of Fig. 10, for the infrared line ratios of Oiv25 µm / Neii12 µm, Nev14 µm / Neii12 µm and Niii57 µm / Nii122 µm. The predictions for these ratios from the simulation are more than 1 dex lower than the RPC limit of the cloudy model. We therefore conclude that these ratios present the strongest constraints between the hot gas pressure-driven and radiation pressure-driven scenarios.
However, not all line ratios can be used to constrain the driving mechanism. In both panels of Fig. 11 the line ratio predictions from our simulations bracket the RPC limit of the idealised cloudy model, with L45 and L46 producing lower and higher ratios than in RPC, respectively. If we ran further simulations covering a wider range of the parameter space, we would expect some of these to produce intermediate ratios between the two simulations shown here, which could thus overlap with the RPC predictions. In the top five panels of Fig. 10, the line ratios predicted by the simulations are also close to the RPC limit. These particular line ratios therefore cannot distinguish between the hot gas pressuredriven scenario and the RPC regime. This contrasts with the predictions of the cloudy model in the P hot P rad limit, in which these line ratios become much lower than in the RPC limit. This discrepancy between the simulations and the hot gas pressure-dominated limit of the cloudy model is due to the lower pressures that we find in the simulations compared to the assumption of hydrostatic balance with the hot phase.
In general, the three infrared ratios most sensitive to the driving mechanism (i.e. the bottom three panels of Fig. 10) compare high ionisation and low ionisation gas. We saw in Fig. 5 that the high ionisation states are found in gas at lower densities that have had less time to undergo compression, while the low ionisation states trace gas that exhibit stronger compression and are thus at higher densities. Since the compression by the hot bubble produces the dense lowionisation gas, the ratios between high and low ionisation states are most sensitive to the existence of the hot bubble and hence to P hot /P rad . However, line ratios such as Nevi7 µm / Neiii15 µm in the simulations are consistent with the RPC limit.

ENERGETICS OF DIFFERENT PHASES IN THE AGN OUTFLOW
Emission lines are vital tools for measuring the energetics of galactic outflows. This is important for understanding how the AGN can influence its host galaxy, for example whether it can provide sufficient energy to unbind the gas in the galaxy and hence quench star formation (e.g. Sturm et al. 2011;Cicone et al. 2014). However, we have seen in Section 3.2 that different lines trace different phases of the outflow. If we only measure the outflow energetics (e.g. mass outflow rate, momentum flux and energy flux) from a single line, we therefore would not capture the total energetics of the outflow. In this section, we investigate what fraction of the total energetics are captured by each emission line. We calculate this from the simulations as follows. First, we calculate the mass outflow rate (ṁpart), momentum flux (ṗpart) and kinetic energy flux (Ė kin, part ) of each gas particle in the highresolution octant of the simulation volume. We define these quantities analogously to how they would be calculated in observations as follows: where mpart, v los, part and rpart are the particle mass, line-ofsight velocity and distance from the black hole, respectively. For each of the ion and molecule species, we then project these quantities on to the AMR grid used for the radiative transfer calculations weighted by the abundance of the given species. Next, we calculate the emissivity ( ) of each emission line in every cell of the AMR grid (i.e. the line luminosity per unit volume). We then sort the cells in order of increasing , and we find the emissivity 90 for which 90 per cent of the total line luminosity orginates from cells with > 90. Finally, we calculate theṁ,ṗ andĖ kin summed over cells with > 90 for each line, which we compare to the same quantities summed over all cells to find the fraction of the total energetics that are traced by each line.
Defined in this way, we are measuring the fraction of the energetics in cells that are 'bright' in a given emission line (such that 90 per cent of the total emission is included). Note that a given cell can be bright in multiple lines, even from different ionisation states or molecules of the same element. For example, a cell containing a mixture of CO and Cii could be bright in both and thus would be accounted for in the energetics of each line; we do not divide a cell between the relative contributions to each line. Thus the fractions that we present here should not be considered as a division of the total energetics between each separate species. Rather, they tell us how much of the energetics can be captured by looking at each individual line, compared to how much arises from gas that is dark in the given line and thus would not be observable.
The results of this analysis are presented in Fig. 12, which shows these fractions for each emission line. The species 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 (ṁ), momentum flux (ṗ) and kinetic energy flux (Ė kin ), denoted by the background colour as light-, mid-and dark-blue, respectively. Results from the simulations L45 and L46 are shown by the grey and black lines, respectively.
We again caution that our predictions are likely to be influenced by uncertainties in the relative contributions of the photoionised and WNM phases, due to strong X-ray heating in the assumed partially obscured incident quasar spectrum and unresolved ionised layers, as discussed in Section 3.3. These results should therefore be considered as a qualitative indication of the relative importance of the different phases traced by each line, but the exact values remain uncertain.
In general, the relative contributions of each line to the total momentum flux and kinetic energy flux are similar to the fraction of the total mass outflow rate that they trace. This is unsurprising given the idealised setup of these simulations, as most of the outflowing mass is located in a thin shell moving radially outwards with similar velocities across the shell. There are some exceptions particularly iṅ E kin /Ė kin, tot , which can differ fromṁ/ṁtot andṗ/ṗtot by factors of a few in some cases (e.g. Cii158 µm, Nev14 µm and Nevi7 µm). The kinetic energy flux has the strongest dependence on velocity, so variations in the velocity distributions of different species will have a larger impact onĖ kin /Ė kin, tot than the other quantities. Nevertheless, the trends that we see here are primarily driven by the fraction of the mass that is located in cells of the AMR grid that are bright in the given line.
The Cii158 µm line (produced in the transition from 10 4 K to 100 K) and the Nii122 µm and Neiii15 µm lines (pro-duced in the 10 4 K phase) trace a large fraction of the mass outflow rate and momentum flux (≈50−70 per cent). The latter two lines also capture most of the kinetic energy flux, although this is somewhat lower in Cii158 µm. It may seem surprising that Neiii traces such a high fraction of the energetics, given that it has an ionisation energy close to species such as Oiii, which trace much lower fractions of the energetics. Comparing these emission lines in Figs. 5 and 6, we see this is because Neiii emission extends to higher densities than the Oiii lines. Photoionisation models of RPC clouds using cloudy also exhibit significant Neiii abundances extending into the neutral region (see e.g. Fig. A1 in Stern et al. 2014a).
Emission from the high ionisation states, such as Oiv25 µm, Nev14 µm and Nevi7 µm, only trace low fractions of the total energetics (typically 10 per cent). As we saw in Section 3.2, these species arise from gas in a transitionary phase as it reaches the end of a period of rapid cooling, before it is subsequently compressed to higher densities due to the external pressure exerted by the hot medium. Thus gas spends a relatively short period of time evolving through this phase.
We again stress that these results do not include dust attenuation from the host galaxy. This will not affect the millimetre and infrared lines, but we would expect those at optical and UV wavelengths to be strongly absorbed. The fractions presented in Fig. 12 for the optical and UV lines should therefore be interpreted as the fraction of the energetics in gas that emits these lines, but might not necessarily be observed if they are strongly attenuated by the host galaxy.

CONCLUSIONS
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 a hot gas bubble generated by a central AGN. These simulations include an on-the-fly treatment for the non-equilibrium chemistry of ions and molecules coupled to the hydrodynamics and radiative cooling. The resulting ion and molecule abundances are used together with a radiative transfer code in post-processing to make predictions for the emission lines. We use these calculations to study how this emission traces the physical conditions and energetics of the outflow, and we compare the predicted emission lines from our models to observations. Our main results are as follows: (i) We find that molecules (CO, H2) and low-ionisation states (Cii and the infrared lines of Oi) trace clumpy structures that have condensed within the outflowing shell as it cooled, while the intermediate species (e.g. Nii, Siii, Neiii) arise from a more diffuse phase throughout the shell. Emission from the highest ionisation states (e.g. Nev, Nevi) is concentrated in small, bright knots, which are produced by regions that are passing through a period of rapid cooling. As such, any one particular region of gas spends a relatively short period of time in this phase, and thus this high ionisation emission appears as short, bright flashes concentrated in small regions (see Figs. 2 and 3). (ii) Comparing the temperature-density distribution traced by each emission line, we find that emission from high ionisation states (e.g. Nev, Nevi) is produced by gas that is coming to the end of a period of rapid cooling. Intermediate ions such as Nii and Oiii arise from gas at the thermal equilibrium temperature of ∼10 4 K, while low-ionisation states (e.g. Cii and the infrared lines of Oi) and molecules trace the transition from the warm (∼10 4 K) to the cold (∼100 K) phase (Figs. 5 and 6). (iii) The hot bubble compresses the line-emitting gas beyond the initial pressure of the ambient ISM by 1−2 orders of magnitude. For many emission lines, it also reaches higher pressures than would be achieved in an outflow driven and compressed by radiation pressure (Figs. 5 and 6). (iv) The gas is under-pressurised compared to the pressure of the hot medium by more than an order of magnitude. Also, while the intermediate ions are in thermal equilibrium at ∼10 4 K, the high ionisation states are up to ∼1 dex higher than the thermal equilibrium temperature, as they are still cooling. The molecules and low-ionisation states trace a broad range of temperatures above the thermal equilibrium, due to turbulent shock heating. This has implications for photoionisation models that compute AGN emission line intensities assuming thermal and/or pressure equilibrium (Figs. 5 and 6). (v) While the luminosities of many of the infrared lines predicted from our simulations overlap with AGN observations, there are notable discrepancies, in particular in the Oi63 µm line, which is an order of magnitude too high in the simulations (Fig. 7). Some of the predicted line ratios at infrared wavelengths are also inconsistent with observations. For example, Nevi7 µm / Oiv25 µm and Nevi7 µm / Nev14 µm are 3× too high, while Oi63 µm / Cii158 µm is 10× too high (Fig. 8). (vi) At optical wavelengths, the BPT diagnostic diagrams show that the Sii 6716+6731Å / Hα 6563Å and Oi 6300Å / Hα 6563Å ratios are ≈1 dex higher in the simulations than the observations (Fig. 9). These anomalous ratios are due to two effects. Firstly, our assumed incident spectrum for a partially obscurred quasar has a ≈2× higher X-ray to UV ratio than a more typical unobscured spectrum, which enhances the Warm Neutral Medium (WNM) supported by X-ray heating. Secondly, we do not fully resolve photoionised layers, which results in an underestimate of the photoionised phase by a factor ≈3−4. Our simulations therefore do not provide a reliable prediction for the BPT diagrams. We also caution that these uncertainties may also affect other emission line predictions from our simulations, in particular those that compare the WNM and photoionised phases. (vii) As the line-emitting gas is compressed by the hot bubble, we find that certain line ratios are sensitive to the ratio of the hot gas pressure to radiation pressure, P hot /P rad , providing a constraint on the driving mechanism of AGN outflows, as suggested by the photoionisation models of Stern et al. (2016) that assumed thermal, chemical and hydrostatic equilibrium and an idealised slab geometry. We find that the ratios of Oiv25 µm / Neii12 µm, Nev14 µm / Neii12 µm and Niii122 µm / Nii122 µm show the strongest distinction between our simulations (for which P hot P rad ) and the radiation pressure-dominated limit of cloudy pho-toionisation models in hydrostatic equilibrium (Fig. 10). We therefore suggest that these line ratios provide the strongest constraints on the relative dynamical importance of radiation pressure and hot gas pressure on the outflow. (viii) We quantify the fraction of the total mass outflow rate, momentum flux and kinetic energy flux of the outflow that is traced by each emission line. We find that the Nii122 µm and Neiii15 µm lines (arising from the 10 4 K phase) trace ≈50−70 per cent of the totals in all three quantities. Cii158 µm (arising from the transition from 10 4 K to 100 K) also traces ≈60−70 per cent of the mass outflow rate and momentum flux, but only ≈30−40 per cent of the kinetic energy flux. Meanwhile, the high ionisation states such as Oiv25 µm, Nev14 µm and Nevi7 µm (produced in a transitionary phase as the gas undergoes rapid cooling) trace 10 per cent of the energetics (Fig. 12).
Our simulations demonstrate that, in outflows driven by the thermal pressure of a hot gas bubble (as can be produced 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, which allows us to constrain the dynamical importance of the hot gas versus radiation pressure, provided that we consider line ratios that are most sensitive to this effect. We also find that much of the emission arises from gas in highly transitionary phases, in particular for high ionisation states such as Nev and Nevi that are produced towards the end of a rapid cooling phase, which leads to the line-emitting gas being out of thermal-, pressure-and chemical-equilibrium in many cases. This highlights the importance of simulations to capture these dynamical effects and their impact on emission line predictions.
Some of the emission line ratios predicted by our simulations are inconsistent with observations of all AGN subtypes, as noted above. However, in this study we only have two simulations, which cover a very limited range of the parameter space. These simulations also use an idealised setup, for example with an initially uniform ambient ISM that lacks turbulent gas structures or a realistic host galaxy geometry. The outflowing shell also remains constrained by the ambient ISM throughout the simulations, and so we cannot probe the regime after the outflow breaks out of the dense regions of the host galaxy. The tensions between our simulations and the observations are therefore insufficient evidence to rule out the hot gas pressure-driven scenario for kiloparsecscale AGN outflows, as we cannot determine whether these tensions are a limitation of the model or if they are simply due to the limited range and idealised nature of the simulations.
To conclusively distinguish between possible driving mechanisms of AGN outflows, we would therefore need to run a more extensive suite of simulations covering a wide range of the parameter space and different physical conditions. The resulting emission line predictions can then be compared to observations to constrain the models, similar to how grids of photoionisation models have been employed to constrain 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 expect this to be achievable in the near future. Figure A1. Continuum-subtracted spectra of the millimetre and infrared emission lines (top panels) and the optical and UV lines (bottom panels) in the full sample, from simulations L45 (dashed curves) and L46 (solid curves). The grey shaded bands indicate velocity ranges that are excluded from further analysis, to focus on the outflow component. Figure A2. Maps of the velocity-integrated continuumsubtracted emission from all millimetre, infrared, optical and UV emission lines in the full sample from simulation L45. Figure A3. As Fig. A2, but for the L46 simulation.