# Outflows driven by quasars in high-redshift galaxies with radiation hydrodynamics

Rebekka A Bieri, Yohan A Dubois, Joakim A Rosdahl, Alexander A Wagner, Joseph A Silk, Gary A. Mamon

► To cite this version:

Rebekka A Bieri, Yohan A Dubois, Joakim A Rosdahl, Alexander A Wagner, Joseph A Silk, et al. Outflows driven by quasars in high-redshift galaxies with radiation hydrodynamics. Monthly Notices of the Royal Astronomical Society, 2017, 464 (2), pp.1854 - 1873. 10.1093/mnras/stw2380. hal-01484553

HAL Id: hal-01484553

https://hal.sorbonne-universite.fr/hal-01484553

Submitted on 7 Mar 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Outflows driven by quasars in high-redshift galaxies with radiation hydrodynamics

Rebekka Bieri,<sup>1★</sup> Yohan Dubois,<sup>1★</sup> Joakim Rosdahl,<sup>2</sup> Alexander Wagner,<sup>3</sup> Joseph Silk<sup>1,4,5,6</sup> and Gary A. Mamon<sup>1</sup>

Accepted 2016 September 16. Received 2016 September 12; in original form 2016 June 20

## ABSTRACT

The quasar mode of active galactic nuclei (AGN) in the high-redshift Universe is routinely observed in gas-rich galaxies together with large-scale AGN-driven winds. It is crucial to understand how photons emitted by the central AGN source couple to the ambient interstellar medium to trigger large-scale outflows. By means of radiation–hydrodynamical simulations of idealized galactic discs, we study the coupling of photons with the multiphase galactic gas, and how it varies with gas cloud sizes, and the radiation bands included in the simulations, which are ultraviolet, optical, and infrared (IR). We show how a quasar with a luminosity of $10^{46}$ erg s$^{-1}$ can drive large-scale winds with velocities of $10^2–10^3$ km s$^{-1}$ and mass outflow rates around $10^3$ M$_{\odot}$ yr$^{-1}$ for times of the order of a few million years. IR radiation is necessary to efficiently transfer momentum to the gas via multiscattering on dust in dense clouds. However, IR multiscattering, despite being extremely important at early times, quickly declines as the central gas cloud expands and breaks up, allowing the radiation to escape through low gas density channels. The typical number of multiscattering events for an IR photon is only about a quarter of the mean optical depth from the centre of the cloud. Our models account for the observed outflow rates of ~500–1000 M$_{\odot}$ yr$^{-1}$ and high velocities of ~$10^3$ km s$^{-1}$, favouring winds that are energy driven via extremely fast nuclear outflows, interpreted here as being IR radiatively driven winds.

**Key words:** methods: numerical – galaxies: active – galaxies: high-redshift – galaxies: ISM.

## 1 INTRODUCTION

Galaxy formation and evolution is highly non-linear and poses significant challenges to our current understanding of the Universe. One particular issue is the relation between the mass function of dark matter haloes, provided by theoretical models, and the luminosity function of galaxies, given by observations and usually fit by a Schechter (1976) function. By assuming that stellar mass follows halo mass, we are left with a theoretical prediction that leads to excessive numbers of galaxies at both the low-mass and the high-mass ends. Thus, some mechanisms have been advocated to regulate the baryon budget in galaxies. Feedback is likely an important mechanism, operating within galaxies and driving large-scale outflows, removing the star-forming gas and/or preventing

further infall. Feedback from stars, especially from supernovae, is assumed to impact the low end of the observational stellar luminosity function, by influencing star formation (e.g. Dekel & Silk 1986; Mac Low & Klessen 2004; McKee & Ostriker 2007) driving turbulence (e.g. Joung & Mac Low 2006; Agertz et al. 2009; Tamburro et al. 2009) and generating galactic scale outflows (e.g. Springel & Hernquist 2003; Oppenheimer & Davé 2006). Since the supernova-driven outflows become less effective as galaxies and haloes grow in mass (e.g. Dubois & Teyssier 2008), the growth of supermassive black holes (BHs) accelerates once the supernova-driven large-scale outflows are hindered (above $M_* > 10^9$ M$_{\odot}$; Dubois et al. 2015; Bower et al. 2016; Habouzit, Volonteri & Dubois 2016). Therefore, efficient feedback in high-mass galaxies is thought to be provided by active galactic nuclei (AGN) hosting supermassive BHs (e.g. Magorrian et al. 1998; Hu 2008; Kormendy, Bender & Cornell 2011). Gas accretion on to a BH leads to energy release capable of driving outflows that regulate star formation and the local baryonic

<footnote n="1"><sup>1</sup>Institut d'Astrophysique de Paris (UMR 7095: CNRS and UPMC - Sorbonne Universités), 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="2"><sup>2</sup>Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, the Netherlands</footnote>

<footnote n="3"><sup>3</sup>Center for Computational Sciences, University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8577, Japan</footnote>

<footnote n="4"><sup>4</sup>Laboratoire AIM-Paris-Saclay, CEA/DSM/IRFU, CNRS, Univ. Paris VII, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="5"><sup>5</sup>Department of Physics and Astronomy, The Johns Hopkins University Homewood Campus, Baltimore, MD 21218, USA</footnote>

<footnote n="6"><sup>6</sup>BIPAC, Department of Physics, University of Oxford, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote>* E-mail: bieri@iap.fr (RB); dubois@iap.fr (YD)</footnote>
content (Silk & Rees 1998), hence regulating the BH growth and that of the surrounding galaxy (Kormendy & Ho 2013).

The two main modes of AGN feedback identified so far are the so-called *radio*-mode (radiatively inefficient) powered by mechanical jets, and the radiatively efficient *quasar*-mode powered by photons that couple to the gas and transfer their momentum to it. Several studies have examined the effects of these different modes in cosmological hydrodynamical (HD) simulations (e.g. Di Matteo, Springel & Hernquist 2005; Sijacki et al. 2007; Di Matteo et al. 2008; Booth & Schaye 2009; Dubois et al. 2012) and also on smaller galaxy scales (e.g. Proga, Stone & Kallman 2000; Ostriker et al. 2010; Wagner & Bicknell 2011; Nayakshin & Zubovas 2012; Gabor & Bournaud 2014), and have highlighted the capacity of regulating the baryon content of galaxies with AGN feedback. The general findings of these HD simulations and also of semi-analytical models (e.g. Bower et al. 2006; Croton et al. 2006) are that AGN feedback suppresses star formation in massive galaxies, reproduces the observed high-end tail of the galaxy mass function, and is largely responsible for the morphological transformation of massive galaxies into ellipticals.

Although the picture of supermassive BHs exerting strong feedback on their host galaxies is very attractive, the details of the mechanism remain vague, primarily because the implementation of the BH feedback in these studies relies on subgrid recipes in pure HD simulations. The implementations in these studies vary slightly, but in most cases, quasar feedback is approximated by depositing thermal energy within the resolution element, with the efficiency of the radiation-gas coupling represented by a single parameter chosen to match global observations of BH mass-bulge velocity dispersion (Ferrarese & Merritt 2000). This is the quasar mode of AGN feedback implemented in recent ~100 Mpc state-of-the-art HD cosmological simulations such as Horizon-AGN (Dubois et al. 2014), Illustris (Vogelsberger et al. 2014), MassiveBlack-II (Khandai et al. 2015), or Eagle (Schaye et al. 2015).

In reality, this subgrid model of quasar feedback hides complex physics. There is a general consensus that photons emitted from a thin accretion disc surrounding the BH (Shakura & Sunyaev 1973) will effectively couple to the surrounding interstellar medium (ISM) and eventually drive a large-scale wind. These radiatively driven winds are not necessarily simply momentum-conserving with $\dot{p} \simeq L/c$ (as proposed by King 2003), where $\dot{p}$ is the momentum input rate, $L$ is the bolometric luminosity of the source, and $c$ is the speed of light, as opposed to energy conserving, where the momentum of the gas builds up from the pressure work of the gas, hence, $\dot{p} \gg L/c$ (Silk & Rees 1998; Faucher-Giguère & Quataert 2012; Zubovas & King 2012). For example, observations by Cicone et al. (2014) show outflows with mechanical advantages of a few tens, where the mechanical advantage is defined as the ratio between $\dot{p}$ and $L/c$. Recent observations of fast ($>10^3$ km s$^{-1}$) molecular outflows from AGN favour energy-conserving winds from the nuclear accretion disc (Feruglio et al. 2015; Tombesi et al. 2015) as do recent multiphase simulations discussed below (e.g. Costa, Sijacki & Haehnelt 2014).

Alternate approaches to simple internal energy input have been implemented, where the gas close to the BH is explicitly given a momentum input rate of $L/c$ multiplied by a factor of a few (Debuhr et al. 2010; Debuhr, Quataert & Ma 2011; Choi et al. 2012, 2014; Barai et al. 2014; Costa et al. 2014; Zubovas & Nayakshin 2014; Hopkins et al. 2016). However, simulations show contradictory results regarding the plausibility of regulating the BH growth and impact on the baryon content of the galaxy using momentum-driven winds. While Costa et al. (2014, see also Barai et al. 2014)

advocate an energy-conserving wind to significantly affect the galaxy and its surroundings, Choi et al. (2012) show that their momentum feedback drives faster winds than their energy feedback model. Since quasar-driven winds are powered by the complex coupling of photons with the gas, it becomes timely to improve our current understanding of the transfer of momentum and energy from AGN-emitted radiation to the ISM by means of radiation-hydrodynamical (RHD) simulations. In such RHD simulations, the emission, absorption, and propagation of photons and their interaction with the gas are self-consistently followed. Detailed RHD simulations should help to improve our understanding in how the quasar radiation couples to the gas, in particular via the coupling of infrared (IR) radiation to dust, and how powerful quasar winds are eventually driven. Therefore, RHD simulations can provide better subgrid models for large-scale cosmological simulations.

RHD simulations of quasar feedback on galactic scales have become feasible (e.g. Ciotti & Ostriker 2007, 2012; Kim et al. 2011; Novak, Ostriker & Ciotti 2012) and allow one to measure how much of the radiation momentum is transferred to the gas. These studies manage to resolve important physics involving the variability of the AGN (the *duty cycle*). RHD simulations also show how a strong radiation source at the centre of a galaxy affects the surrounding ISM. Ciotti & Ostriker (2007, 2012) and Novak et al. (2012) show that, for elliptical gas-poor galaxies, the momentum input rate very rarely exceeds $L/c$ due to the low opacity of dust to the re-radiated IR and the destruction of dust in high-temperature environments, which in turn limits the amount of momentum the radiation can transfer to the gas. Note that these studies investigate the effect of radiation from a BH at the centre of an elliptical galaxy. Thus, they do not directly study the problem of high-redshift galaxies but rather investigate the effect on the galaxy in a later phase at low redshift when the galaxy is passively evolving, and where the accretion on the central BH is fuelled by stellar mass losses. Kim et al. (2011) simulate the evolution of a high-redshift galaxy with radiation from the BH; however, the IR radiation has not been considered in this work, and, hence, the impact of radiation into the gas is probably underestimated. Similarly in Roos et al. (2015), the effect of the ultraviolet (UV) photoheating only from a central quasar (treated in post-processing) has very little impact on the evolution of a gas-rich isolated disc galaxy. The IR radiation is potentially important because in optically thick gas it is constantly absorbed by dust and re-emitted, and hence can give several times its momentum $L_{IR}/c$ to the dust up to the dust optical depth $\tau_{IR} = \Sigma_g \kappa_{IR}$, where $\Sigma_g$ is the gas surface density and $\kappa_{IR}$ is the dust opacity. It can create radiation pressure-driven winds where the momentum input rate exceeds the momentum flux of the source by one or two orders of magnitude, mimicking the effect of an energy-driven wind. Star clusters in starburst galaxies can easily reach values of $\tau_{IR} = 10-100$ (see fig. 3 of Agertz et al. 2013), possibly explaining the observed energy-conserving winds in high-redshift quasars.

High-resolution HD simulations of the AGN jet feedback (e.g. Bicknell et al. 2000; Sutherland & Bicknell 2007; Antonuccio-Delegu & Silk 2010; Wagner & Bicknell 2011; Gaibler et al. 2012) have shown that a clumpy interstellar structure results in interactions between the jet and the gas that differ from those in simulations with a homogeneous ISM. It is hence expected that multiphase ISM properties also affect the momentum transfer from the radiation to the gas. Modelling a clumpy interstellar structure might be even more important when considering photons, as low-density gas can trace escape paths in which the radiation-matter interaction is reduced because the radiation field is not fully trapped. Also, it is expected that as the radiation is sweeping up the dense gas, low-density channels
form and photons start to escape along these preferred directions, lowering the mechanical advantage (Krumholz & Thompson 2012, 2013; Davis et al. 2014; Rosdahl & Teyssier 2015).

The aim of this paper is to quantify the coupling of quasar radiation with a clumpy ISM and to study how photons (UV, optical, and IR) can drive powerful winds using RHD simulations. Since radiation-gas coupling depends on the clumpiness of the gas, we consider a two-phase (see McKee & Ostriker 2007) fractal ISM at pc-scale resolution. We investigate the momentum budget associated with the dispersion of the clouds in the galaxy and its dependence on different cloud sizes, filling factors, quasar luminosities, and energy bands of the quasar spectrum.

In Section 2, we describe our suite of RHD simulations. Our results are presented in Section 3. Section 4 is a discussion of the caveats of our methods and setup. Section 5 provides the final conclusions.

## 2 METHODS

We perform a suite of simulations using RAMSES-RT, an RHD extension of the adaptive mesh refinement (AMR) code RAMSES (Teyssier 2002). We model quasar-emitted radiation interacting with a surrounding multiphase ISM, in order to study how efficiently radiation couples to the gas within the galaxy.

### 2.1 Initial gas density distribution

We set up a gaseous disc with a two-phase ISM in pressure equilibrium, with a uniform hot phase with temperature $T \sim 10^6$ K, and a cold $T \sim 1-10^4$ K phase that is uniform on large scales, but very clumpy on small scales. Our ISM setup is described in detail in Wagner, Umemura & Bicknell (2013), which followed work done by Sutherland & Bicknell (2007) investigating the interaction of an AGN jet with a non-uniform ISM.

The density field for the cold phase is homogeneous on large scales, i.e. there is no radial density gradient, but it is very clumpy on small scales: it follows a single-point lognormal distribution and two-point fractal statistics. In a lognormal distribution, the logarithm of the ISM density is a Gaussian, with mean $m$ and variance $s^2$. With $P(\rho)$ being the lognormal probability distribution of the mass density $\rho$, one can write

$$P(\rho) = \frac{1}{\sqrt{2\pi s \rho}} \exp \left[ -\frac{(\ln \rho - m)^2}{2s^2} \right], \quad (1)$$

with

$$m = \ln \left( \frac{\mu^2}{\sqrt{\sigma^2 + \mu^2}} \right), \quad s = \sqrt{\ln \left( \frac{\sigma^2}{\mu^2} + 1 \right)}, \quad (2)$$

where $\mu$ and $\sigma^2$ are the mean and the variance of the linear density field. In the simulations presented here, we adopt $\mu = 1$ and $\sigma^2 = 5$, identically to what has been used by Wagner & Bicknell (2011). These values are in agreement with ranges found by Fischera, Dopita & Sutherland (2003) and Fischera & Dopita (2004), observing the column density distribution in a turbulent ISM. The variance in equation (2) gives a measure of how concentrated the mass is within the density cores, or, conversely, the fraction of the volume within the low-density regions. With these values, gas densities below the mean $\mu$ encompass one-quarter of the mass and occupy three-quarters of the volume of the simulated disc. For further discussion of the adopted values, we refer to Bicknell et al. (2000).

We define $F(k)$ to be the Fourier transform of the density $\rho(r)$, where $k$ and $r$ are the wavenumber and position vectors, respectively. The two-point structure of a homogeneous turbulent medium is characterized in Fourier space by an isotropic power spectrum $D(k)$ defined as

$$D(k) = \int k^2 F(k) F^*(k) d\Omega. \quad (3)$$

The power-spectrum is proportional to a power law with index $-5/3$ in order to reproduce the spectrum driven by Kolmogorov turbulence. The Fourier transform of the density $\rho$ is proportional to the turbulent field usually described by the velocity vector (Warhaft 2000). We follow Wagner & Bicknell (2011) and adopt a standard Kolmogorov power spectrum for our non-uniformly distributed gas within the disc.

We want to stress that our initial setup is stationary and hence does not capture the actual ISM turbulence (see Kritsuk, Norman & Wagner 2011 and references therein). We rather parametrize the non-uniform properties of a generic turbulent medium and focus on characteristics such as the variance of the gas density $\sigma^2$ and the two-point self-similar power-law structures, and hence rely on a range of previous experimental and theoretical results from the field of turbulence. The initial distribution of the ISM that we adopt should therefore be regarded as a physically motivated generalization of an inhomogeneous ISM, while it may not necessarily represent an accurate model of a turbulent ISM.

The density distribution is set up via an iterative process, described first by Lewis & Austin (2002), where we have adapted the PYFC package (Wagner & Bicknell 2011) to our needs (mainly parallelized part of the code in order to have a smaller minimum sampling wavenumber). Our density field simultaneously follows lognormal single-point statistics (see equation 1) and a Kolmogorov power-law self-similar structure in wavenumber with index $-5/3$ and minimum sampling wavenumber $k_{min}$. In real space, the minimum sampling wavenumber determines the scale of the largest fractal structure in the cube relative to the size of the cube. Effectively, it is the average number of clouds per dimension divided by 2. For example, in one setup of our simulation, we used $k_{min} = 5$ kpc$^{-1}$ for a cube mapped to a disc with a diameter of $d = 3$ kpc. Then the largest structures (clouds) extend to $R_{c,max} = 3/(2k_{min}) = 300$ pc.

Finally, we place the cube into the RAMSES simulation domain by reading each density value within the cube and placing it into a cell of the RAMSES grid. Here, the resolution of the generated cube does not necessarily need to be the resolution of the smallest cell within the RAMSES simulation box. In order to obtain a cylindrical shape resembling a galaxy, the density cube is filtered, in the xy-plane, by a symmetric flat mean density profile with mean cold phase density $\langle n_w \rangle$ and radius $r = 1.5$ kpc, and in the z-plane the density cube is filtered by a step function with height $h = 0.3$ kpc. The porosity of the ISM arises by imposing a temperature roof $T_{roof}$ above which the gas is defined to be in the hot phase. The mean density of the cold ISM phase is chosen such that the total mass of the cold phase is $\sim 2 \times 10^{10} M_{\odot}$ and is around 500 H cm$^{-3}$. The roof temperature is around 50 K for the different simulations.

In our simulations, clouds are initially in pressure equilibrium with the surrounding hot phase, where the pressure is set to be $P \sim 7 \times 10^{-12}$ Pa. For the hot phase, whose temperature is fixed at $T_h \sim 10^7$ K, the density is constant and set as $n_{H,h} = 0.01$ H cm$^{-3}$ for all the simulations, while the initial metallicity is set to zero. The cold gas is initialized with solar metallicity. With these values and with a disc radius of 1.5 kpc and thickness of 0.15 kpc, we simulate
a typical compact, gas-rich, high-redshift galaxy (e.g. Daddi et al. 2010; Genzel et al. 2010; Tacconi et al. 2010).

The density and pressure profiles of the hydrostatic environment in a massive gas-rich protogalaxy are fairly flat under the gravitational influence of the bulge and dark matter halo (assuming a Navarro, Frenk & White 1996 profile; e.g. Capelo, Natarajan & Coppi 2010). Hence this justifies the uniform hot phase distribution adopted in our simulation.

The filling factor of cold phase within the disc is given by

$$f_v = \frac{V_{\text{cold}}}{V_{\text{tot}}}, \qquad (4)$$

where $V_{\text{cold}}$ is the volume of the cold phase and $V_{\text{tot}}$ the total cylindrical volume of the region in which the density is distributed.

Since the minimum sampling wavenumber relates to the largest fractal structure in the cube, we will hereafter, respectively, refer to the simulations with $k_{\min} = 30 \text{ kpc}^{-1}$, $5 \text{ kpc}^{-1}$, and $1 \text{ kpc}^{-1}$ as smallC (smallest clouds), medC (medium clouds), and bigC (biggest clouds). For our simulations, we chose a mean cold phase density $\langle n_w \rangle$ of 508, 503, and 435 H cm$^{-3}$ for the smallC, medC, and bigC simulations, respectively. Additionally, we chose roof temperatures of 43, 70, and 70 K for the smallC, medC, and bigC simulations, respectively.

## 2.2 Radiation hydrodynamics

Among the possible RHD implementations that are both helpful for cosmological and galaxy-scale simulations (e.g. Petkova & Springel 2009; Krumholz, Klein & McKee 2011; Pawlik & Schaye 2011; Wise & Abel 2011; Jiang, Stone & Davis 2012; Skinner & Ostriker 2013), we chose RAMSES-RT (Rosdahl et al. 2013; Rosdahl & Teyssier 2015), implemented in the RAMSES (Teyssier 2002) AMR HD code, to model the interaction of radiation from the central BH with the galaxy’s interstellar gas. The evolution of the gas is computed using a second-order unsplit Godunov scheme for the Euler equations. We use the HLLC Riemann solver (Toro, Spruce & Speares 1994) with the MinMod total variation diminishing scheme to reconstruct the interpolated variables from their cell-centred values. The RAMSES-RT RHD extension to RAMSES self-consistently adds the propagation of photons and their on-the-fly interaction with hydrogen and helium via photoionization, heating, and momentum transfer, as well as their interaction with dust particles via momentum transfer. The advection of photons between grid cells is implemented with a first-order moment method, whereas the set of radiation transport equations is closed with the M1 relation for the Eddington tensor (Rosdahl et al. 2013). The M1 closure relation (Levermore 1984) can establish and retain bulk directionality of photon flows, and can to some degree model shadows behind opaque obstacles.

The radiation is split into different photon groups, defined by frequency bands. For each photon group, the radiation is described, in each grid cell, by the radiation energy density (energy per unit volume) and the bulk radiation flux (energy per unit area per unit time), which corresponds to the radiation intensity integrated over all solid angles. RAMSES-RT solves the non-equilibrium evolution of the ionization fractions of hydrogen and helium, along with photon fluxes and the gas temperature in each grid cell. For the lowest energy IR group, the photons can give momentum to the gas multiple times via absorption and re-emission. Higher energy radiation groups are absorbed by dust and re-emitted (conserving energy) into the IR group, but they can also interact with hydrogen and helium via photoionization. We use a subgrid scheme to account for the trapping of IR photons in regions where the mean free path is smaller than the grid spacing. This scheme recovers the proper

asymptotic limit in the radiation diffusion regime (see Rosdahl & Teyssier 2015 for a detailed discussion).

Since the Courant condition imposes that the time-step duration (and therefore the computational load) scales inversely with the speed of light c, we apply the so-called reduced speed of light approximation (RSLA; see also Gnedin & Abel 2001; Rosdahl et al. 2013). The rationale for the RSLA is that as long as the radiation travels faster than ionization fronts, the results of RHD simulations are more or less converged with respect to the (reduced) speed of light.

However, IR radiation is not photoionizing, so it is not obvious whether a reduced speed of light produces converging results, especially when IR trapping becomes important. For our simulations, we chose a reduced speed of light fraction $c_{\text{red}}/c = 0.2$, leading to $c_{\text{red}} \sim 6 \times 10^4 \text{ km s}^{-1}$. We test our chosen reduced speed of light by performing convergence tests in Appendix A but leave a detailed discussion to a forthcoming paper.

We implement outflow boundary conditions, such that any matter that leaves the simulation volume is lost to the system. We however choose a sufficiently large simulation box size $L_{\text{box}} = 96 \text{ kpc}$ to ensure negligible mass-loss. The simulation employs a coarse grid of cell size $L_{\text{box}}/2^9 = 187 \text{ pc}$ and allows up to five additional levels of refinement, so that the smallest cell size is $L_{\text{box}}/2^{14} = 5.9 \text{ pc}$. The refinement is triggered with a quasi-Lagrangian criterion ensuring that if the gas mass within a cell is larger than 80 M$_{\odot}$, a new refinement level is triggered. We ensure that the entire disc is maximally refined at the beginning of the simulation, leading to $\sim 10^7$ cells initially at the maximum level of refinement. For convergence studies, we have performed lower resolution runs with a spatial resolution of $\Delta x = 11.6 \text{ pc}$ using a maximum of six levels of additional refinement.

The gas in our simulations follows the equation of state for an ideal monoatomic gas with an adiabatic index of $\gamma = 5/3$. To keep our setup as simple as possible and only probe the effect of radiation-matter coupling, we neglect gas cooling,<fnref n="1" /> star formation, feedback from stars, and gravitational forces, but these effects will be studied in future work. We discuss those caveats in Section 4.

## 2.3 Modelling the quasar

The spectrum of radiation from the central quasar is modelled using a maximum of five photon groups, defined by the photon energy bands listed in Table 1. The photon groups consist of one IR group (0.01–1 eV), one optical group (‘Opt’ from 1–13.5 eV), and three groups of ionizing UV photons, the first two (UV1 and UV2) bracketed by the ionization energies of H I, He I, and He II, and the third (UV3) extending from He II ionization to soft X-rays (1 keV). For a given quasar luminosity, we split the spectral energy distribution (SED) of a typical quasar, as calculated by Sazonov et al. (2004) and shown in Fig. 1, into these five photon groups, and the corresponding fractions of the quasar luminosity are given in Table 1. Sazonov et al. calculated the typical quasar SED using published AGN composite spectra in the optical, UV, and X-rays, also considering the cosmic X-ray background and the contribution of AGN to IR wavelengths, and the estimated local mass density of BHs. Note that we do not model the hard X-ray energy band ($E > 1$ keV), where photons can heat (or cool) the gas to an equilibrium temperature of $\sim 2 \times 10^7$ K through Compton (or inverse Compton) scattering of electrons, but this happens at a very small distance from

<footnote n="1">¹ Hence the only role of metals in our simulations is to set the dust opacities, as shown by equation (5).</footnote>
<table><caption>**Table 1.** Properties of the photon groups used in the simulations. Columns are name; minimum and maximum energies; cross-sections to ionization by H i, He i, and He ii; dust opacity; luminosity fraction (from quasar spectrum of Sazonov et al. 2004). The dust opacity $\bar{\kappa}$ for each group scales with the gas metallicity $\kappa_i = \bar{\kappa}_i Z/Z_\odot$. The luminosity fraction per photon group is used to calculate the designated energy from the quasar that goes into each corresponding photo group.</caption><thead><tr><th>Photon group</th><th>E<sub>min</sub> (eV)</th><th>E<sub>max</sub> (eV)</th><th>σ<sub>H1</sub> (cm<sup>2</sup>)</th><th>σ<sub>He1</sub> (cm<sup>2</sup>)</th><th>σ<sub>He1</sub> (cm<sup>2</sup>)</th><th>κ̃ (cm<sup>2</sup> g<sup>-1</sup>)</th><th>Energy fraction</th></tr></thead><tbody><tr><td>IR</td><td>0.01</td><td>1.</td><td>0</td><td>0</td><td>0</td><td>10</td><td>0.300</td></tr><tr><td>Opt</td><td>1.</td><td>13.5</td><td>0</td><td>0</td><td>0</td><td>10<sup>3</sup></td><td>0.250</td></tr><tr><td>UV1</td><td>13.5</td><td>24.6</td><td>3.1 × 10<sup>-18</sup></td><td>0</td><td>0</td><td>10<sup>3</sup></td><td>0.079</td></tr><tr><td>UV2</td><td>24.6</td><td>54.4</td><td>4.7 × 10<sup>-19</sup></td><td>4.2 × 10<sup>-18</sup></td><td>0</td><td>10<sup>3</sup></td><td>0.067</td></tr><tr><td>UV3</td><td>54.4</td><td>10<sup>3</sup></td><td>1.1 × 10<sup>-20</sup></td><td>2.3 × 10<sup>-19</sup></td><td>1.7 × 10<sup>-19</sup></td><td>10<sup>3</sup></td><td>0.085</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Broad-band spectrum of a typical quasar (adopted from Sazonov, Ostriker & Sunyaev 2004). The shaded areas show the energy range of each photon group: IR (red), optical (green), UV1 (blue), UV2 (yellow), and UV3 (pink). In the plot, we show the fraction of the total energy going into each photon group (per cent energy) as well as the energy per photon (Avg. energy) within the groups. Note that we do not model the hard X-ray energy band, which contributes 22 per cent to the total energy of the quasar.</figcaption></figure>

the source, and hence, these photons have little overall impact on the ISM (e.g. Ciotti & Ostriker 2012; Hopkins et al. 2016). The ionization cross-sections $\sigma(E)$ are taken from Verner et al. (1996). For each photon group, the cross-sections are luminosity-weighted averages over the energy interval, as described in Rosdahl et al. (2013).

Dust opacities are important for calculating the momentum transfer between the photons and the gas as it happens via scattering on dust. Interstellar dust grains are destroyed, mostly by sputtering, by shock-heated gas above temperatures of $T \ge 10^5$ K (Draine & Salpeter 1979, albeit the exact temperature and destruction time-scale depend on the dust grain size). For this work, we modified the calculation of the dust opacity in RAMSES-RT to include dust destruction by thermal sputtering when the gas temperature is above a cutoff temperature $T_{\text{cut}} = 10^5$ K. This is necessary when the source luminosity is sufficiently high to heat up the gas to very high temperatures. The cutoff temperature for sputtering is weakly dependent on density, an effect we ignore here. We did not include dust sublimation since RAMSES-RT does not follow the dust temperature (which sublimes above $10^3$ K) and leave this for future work.

In our simulations the opacities for the different photon groups are given by

$$ \kappa_i = \bar{\kappa}_i \frac{Z}{Z_\odot} \exp\left(-\frac{T}{T_{\text{cut}}}\right). \qquad (5) $$

We thus assume that the dust content simply scales with the metallicity of the gas below the cutoff temperature.

For the IR, we assume an opacity of $\kappa_i = \kappa_{\text{IR}} = 10(Z/Z_\odot) \text{cm}^2 \text{g}^{-1}$, whereas for the higher energy photons (optical and UV) we assume $\kappa_i = 1000(Z/Z_\odot) \text{cm}^2 \text{g}^{-1}$, i.e. one hundred times higher than of the IR. The chosen opacities are physically motivated by a combination of observations and dust-formation theory of the ISM and stellar nurseries (Semenov et al. 2003 for the IR photons, and Li & Draine 2001 for higher energy radiation). They are however uncertain by a factor of a few, due to model uncertainties as well as the temperature dependence on the opacity, which is ignored in our simulations. Past studies have used similar values (e.g. Hopkins, Quataert & Murray 2011; Agertz et al. 2013; Roškar et al. 2014). The usual IR opacities are in the range of $\kappa_{\text{IR}} = 5-10 \text{cm}^2 \text{g}^{-1}$ and hence our assumed IR opacity is at the high end of what is usually considered.

Table 1 lists the values for the photon group energies, where for each photon group, the energy intervals are between the lower bound $E_{\text{min}}$ and upper bound $E_{\text{max}}$. Also shown are the photoionization cross-sections ($\sigma_{\text{H1}}$, $\sigma_{\text{He1}}$, and $\sigma_{\text{He1}}$) for hydrogen and helium, calculated as described above, the dust-interaction opacities ($\bar{\kappa}_i$), and fractions of total quasar luminosity emitted into each photon group.

We have chosen two different quasar luminosities $10^{43}$ and $10^{46}$ erg s$^{-1}$ for our simulations. The bolometric quasar luminosity function (QLF) at redshift $z=3$, compiled by Hopkins, Richards & Hernquist (2007), shows that the chosen quasar luminosities are very common and nicely bracket the QLF.

In this work, we explore the effect of an AGN radiation source that is steady, isotropic, and located at the centre of the galaxy. In reality, the quasar luminosity is proportional to its accretion rate, which varies in time. Additionally, the BH can move around the potential well of the dark matter halo and galaxy and hence is not always exactly located at the centre of the galaxy. Usually, a lower (higher) density environment around the BH results in a lower (higher) luminosity of the quasar. Moreover, changing the location of the BH changes the optical depth around the source, as the encompassing density changes, which should alter the level of coupling between the radiation and matter. For this reason, we have ensured that the initial conditions for the density field are such that the BH for the $10^{46}$ erg s$^{-1}$ simulations is in the densest region of the entire density distribution, while the initial conditions place the lower luminosity source at an intermediate density. Aside from this, the statistical properties (as discussed above) of the simulations with both quasar luminosities are exactly the same. We have also performed a simulation with a $10^{46}$ erg s$^{-1}$ quasar surrounded by the exact same density distribution as for the $10^{43}$ erg s$^{-1}$ simulations to assess the role of quasar luminosity independently of the density around the quasar.
<table><caption>**Table 2.** Simulation parameters: quasar luminosity (*L*), largest possible cloud size (*R*<sub>c, max</sub>), and the gas density ('environment') around the quasar. If, for instance, the quasar environment is denoted with *max* the quasar is placed into the cell within a maximum density. Additionally, no suffix about the quasar position is used if the position of the quasar is in a maximum density environment for the L46 simulation, or in a medium density environment for the L43 simulation, respectively.</caption><thead><tr><th>Simulation name</th><th>log L (erg s<sup>-1</sup>)</th><th>R<sub>c, max</sub> (kpc)</th><th>Quasar environment</th></tr></thead><tbody><tr><td>L46_smallC</td><td>46</td><td>0.05</td><td>max</td></tr><tr><td>L46_medC</td><td>46</td><td>0.3</td><td>max</td></tr><tr><td>L46_bigC</td><td>46</td><td>1.5</td><td>max</td></tr><tr><td>L46_bigC_medρ<sub>Q</sub></td><td>46</td><td>1.5</td><td>med</td></tr><tr><td>L46_bigC_minρ<sub>Q</sub></td><td>46</td><td>1.5</td><td>min</td></tr><tr><td>L43_smallC</td><td>43</td><td>0.05</td><td>med</td></tr><tr><td>L43_medC</td><td>43</td><td>0.3</td><td>med</td></tr><tr><td>L43_bigC</td><td>43</td><td>1.5</td><td>med</td></tr><tr><td>L43_bigC_maxρ<sub>Q</sub></td><td>43</td><td>1.5</td><td>max</td></tr></tbody></table>

In summary, the key model parameters used in the different simulations are the quasar luminosity, the radius of the largest fractal structure, and the location of the quasar Table 2 summarizes the assigned and derived parameters that we used for our simulations and Fig. 2 shows an example of the *smallC*, *medC*, and *bigC* cloud distributions, where the volume filling factor is kept at 50 per cent. We have also performed simulations with a filling factor of 100 per cent, but the results do not significantly change from those shown in the paper.

## 3 RESULTS

We now present our simulation results and examine the interplay of the BH-emitted radiation with the surrounding gas, focusing on the momentum transferred from the radiation on to the gas. We start with a comparison of the effects of radiation on the galaxy for different cloud sizes for the more luminous $L = 10^{46}$ erg s$^{-1}$ source. Then we investigate how the wind becomes radiatively driven. We next compare the effects of different quasar positions relative to the clouds, and end with a comparison of different quasar luminosities. Additionally, we include in Appendix A a study of how the transferred momentum depends on the speed of light and conclude that this approximation is a crucial component in recovering the correct mechanical advantage from the radiation on to the gas for the first few Myr.

### 3.1 Effects of different cloud sizes

We explore the impact of radiation on the gas clouds and velocity evolution for galaxy discs with different cloud sizes. We focus here on the L46 simulations, where the source is embedded into the high-density environment of the clouds.

#### 3.1.1 Qualitative effects of cloud sizes

In Fig. 2, we show maps of the gas density in the disc at different times for the L46_smallC, L46_medC, and L46_bigC simulations. Comparing the three different cloud sizes, we observe that the outflow is less symmetric for less uniform initial conditions, i.e. larger clouds. For all the simulations, the radiation from the central source destroys the encompassing high-density cloud hosting the source before reaching a lower den-

sity environment. The ionized outflow generated by the radiation pushes the gas out into the circumgalactic medium through the lower density channels, which are more prominent with larger clouds.

Until the cloud is destroyed, the source is surrounded by a high-density environment where the optical depth is sufficiently high for the IR radiation to boost the momentum transfer to the gas. Once the radiation manages to destroy the central cloud, the optical depth drops to lower values, reducing the effect of momentum-boost from IR multiscattering. The time it takes for the radiation to break from the central source through the cloud is hence crucial and is much shorter for the L46_smallC simulation than for the larger-cloud L46_medC or L46_bigC simulations. This gives the photons in the L46_medC and L46_bigC simulations more time to transfer momentum to the gas via the trapped photons in the optically thick gas.

After the radiation has destroyed the encompassing cloud, the gas from the cloud continues to expand until it reaches the neighbouring overdensities. Since the clouds distributed within the L46_smallC simulation are all relatively small, they are rapidly destroyed by the radiation, and the cloud gas efficiently mixes with the background gas and fills up tunnels of lower density, quickly creating smooth isotropic shells of outflowing gas. Clouds in the L46_medC, and L46_bigC simulations are much bigger, and, at $t = 1$ Myr, the shell evolution is less spherically symmetric compared to the L46_smallC simulation. Additionally, the outflow generated by the radiation creates a shell of swept-up gas, which is most apparent in the L46_bigC and also at early times of the L46_medC simulation, while less dominant for the L46_smallC simulation. Indeed, the larger the clouds, the more gas mass they contain and the greater (and hence more apparent) the overdensity of the shell can become.

#### 3.1.2 Cloud evolution

Fig. 3 illustrates the evolution of the mass-weighted cloud mean density, temperature, and H II fraction for the different L46 simulations. The lines show the mean values of all gas with metallicity larger than 0.5 Solar. This picks out the gas originally belonging to dense clouds, since the gas within clouds is initialized with solar metallicity, while the diffuse gas outside them is initially metal free.

Looking at the mean density evolution for the three simulations, we see that the clouds start to disperse right from the beginning, causing the mean density of the clouds to drop. The decrease in mean density is similar for the different simulations. With the mean density of the clouds declining with time, the mean free path of the IR photons also increases and one therefore expects the IR photons to scatter less within the gas.

For the L46_medC and L46_bigC simulations, the mean temperature rises within a very short amount of time (~0.1 Myr) to $10^{4.3-4.5}$ K, whereas it takes ~10 Myr for the L46_smallC simulation to reach similar values. Additionally, the H II fraction of the clouds also increases with time, with a roughly 30 per cent ionization fraction for the gas within the clouds at the end of the simulation. Photoionization proceeds much more slowly for the small cloud than for the more massive clouds.

As seen in Fig. 2, the early evolution of the density evolution is sensitive to small-scale inhomogeneities. In the L46_smallC simulation, the photons manage to quickly destroy the encompassing cloud and then escape through lower density channels without efficiently heating, ionizing, and pushing the higher density gas. For
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Slices of the gas density for the *L46_smallC* simulations (top row), *L46_medC* simulations (middle row), and for the *L46_bigC* simulations (bottom row). The different columns show different times as labelled. The galaxies are shown both face-on (upper portion of rows) and edge-on (bottom portion of rows). The galaxy is destroyed by radiation, and a large-scale radiatively driven wind is generated. While in the *smallC* simulation the wave induced by the radiation expands almost uniformly in the *xy*-plane, this is not the case for the *medC* and *bigC* simulations where the outflow escapes via tunnels of low density. The red square overplotted at the 0 Myr plot for the *bigC* simulation corresponds to the zoomed-in region of Figs 4 and 5, whereas the red square in the *medC* simulation corresponds to the zoomed-in region in Fig. 8.</figcaption></figure>

the bigger cloud simulations, on the other hand, the radiation is trapped for longer within the encompassing cloud, leading to the photons to interact with the higher density gas. On long time-scales, the radiation has smoothed out the inhomogeneities, and hence the evolution of the mean temperature and ionization fraction is similar for all the cloud masses, and converges to the same values since the total masses are the same.

Fig. 4 shows a close-up of a slice of the gas density, gas temperature, H II fraction, and velocity from the *L46_bigC* simulation as indicated by the red square in Fig. 2. The position of the quasar is just outside the zoomed-in slice, on the right-hand side. Fig. 5 displays the time evolution of the radiation fluxes in two wavebands: IR and optical + UV, for the *L46_bigC* simulation in the same zoomed-in region as in Fig. 4.

At the start of the simulation (left-hand panels of Figs 4 and 5), the gas is in a fractal two-phase medium, with a cold, fully neutral, high-density component and a hot, fully ionized low-density component, with both components at rest.

At 0.1 Myr, the IR photons have already penetrated the whole cloud. Where the IR photon flux is high, the dense gas has already started moving outwards as indicated by the moderate velocities (10–50 km s⁻¹) on the right edge of the zoomed-in region. During this early stage, the gas in the high IR flux region is already at higher temperatures than the rest of the gas within the cloud (second panel of the second row of Fig. 5) but has not yet shock-heated to very high temperatures. Given the moderate temperatures well below dust destruction, the cloud is not yet transparent and the IR photons multiscatter within it, giving a boosted push from the inside-out. At
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Evolution of the mean density (top), temperature (middle), and H II fraction (bottom) of the clouds as a function of time for the L46SmallC (blue), L46MedC (green), and L46BigC (red) simulations. The lines show the mass-weighted mean value for the clouds, using only the cells with a metallicity above a cutoff value of 0.5 Solar (where the maximum and minimum metallicities are Solar and zero). The radiation manages to disperse the clouds, as can be seen by the decreasing mean densities. The mean temperature as well as ionization fraction increase with time. This lowers the influence of the radiation as time passes due to the gas being either transparent and/or ionized.</figcaption></figure>

0.1 Myr, the cloud is still fully neutral as the UV photons do not penetrate into the cloud (shown in Fig. 5) because the gas is too dense and UV shielded.

At 1 Myr, the gas has begun to move away from the radiation source. To guide the eye, we have marked the outflowing gas at 3 Myr with selected contours. The magenta contours at density 4000 H cm⁻³ mark the shock front. At 1 Myr, the shock front is still neutral with a temperature of roughly 10³ K and is travelling at a velocity of ~50 km s⁻¹. Given the mean density of the front at 1 Myr of ~4 × 10³ H cm⁻³, the mean free path, $1/(\kappa\rho)$, of IR photons is only ~1.8 pc, which indicates the IR radiation is trapped and multiscattering in this region allowing the IR radiation to efficiently transfer momentum on to the gas. In the regions where the gas is dense and the IR flux is high, the temperature of the gas is also higher compared to the rest of the cloud indicating a shock. Additionally, in the regions of high IR flux, the IR radiation is mixing the multiphase gas and creating a more uniform structure.

At the right edge of the cloud, i.e. in the direction of the source, the diffuse gas is heated up to high temperatures (~10⁶ K), is fully ionized, and has a velocity of ~10³ km s⁻¹. When looking at the

photon flux of the different radiation groups (see Fig. 5), we see that, unlike the IR photons, the UV and optical radiation still does not deeply penetrate the dense cloud. We will later see in more detail that the early evolution of dense clouds is clearly dominated by the IR radiation, whereas the UV radiation plays a more important role at later times.

When the radiation illuminates dense clouds, the gas is dispersed starting from the illuminated side, as seen in the density slice at 3 Myr. To guide the eye, we have again marked the outflowing gas with two selected contours. The magenta contour at density 4 × 10³ H cm⁻³ marks the shock front, while the yellow contour at density 2.5 H cm⁻³ marks the outflow tail before the density drops to the background density. The dispersed gas moves away from the source at a velocity of ~50–150 km s⁻¹ and has a temperature of around ~10⁵ K. The gas flowing away from the dispersed cloud still has only a marginally smaller density than the cloud, but is smoothed out by the IR radiation. A large portion of the outflow behind the density ridge (magenta contour) is fully ionized, whereas the cold cloud preceding the outflow is still neutral (see also Fig. 3). As seen in Fig. 5, the UV radiation has only just reached into the tail edge of the outflowing gas (yellow contour), started to ionize the gas from the outside and heat it via photoionization, and hence push it further from the back end. The rest of the gas is thus ionized via collisional ionization. We will see in Section 3.1.3 that gas is accelerated by the radiation pressure from the UV photons and by photoionization heating, but UV photons contribute mostly to the overall radiatively driven wind once they are reprocessed in the IR. The gas reaches temperatures of around ~10⁵ K and velocities up to 500 km s⁻¹ beside the neutral gas, but still within the smooth outflowing gas from the clouds (left from the yellow contour in Fig. 4). In the regions where the gas mixed with the background gas (right from the yellow contour), the gas reaches temperatures up to 10⁶.⁵ K and velocities of up to 1000 km s⁻¹. At t = 3 Myr, the cloud size with a mean density of ~10³ H cm⁻³ is now comparable to the mean free path of the IR radiation, ~30 pc, causing the IR flux, and thus the influence of the IR photons, to decrease.

In summary, we find that for smaller, more fragmented clouds encompassing the radiation source, the radiation has a tendency to escape, and hence it is less efficient at heating and ionizing the gas compared to larger and more coherent surrounding structures which efficiently trap the radiation. Focusing on a single cloud close to (but separated from) the source of radiation, we see that the IR efficiently penetrates the cloud and pushes from the inside-out, smoothing out inhomogeneities, while the UV (and optical) radiation cannot penetrate and acts more by pushing on and heating into the side of the cloud.

### 3.1.3 Velocity evolution

Fig. 6 shows mass-weighted velocity–density diagrams and velocity–temperature diagrams for the L46MedC simulation at different times (1, 3, and 7 Myr from left to right). Three regions can be highlighted in the velocity–density diagram. At low densities, there is a range of velocities in distinct intervals. In light of Fig. 4, we see that the low-density gas far from the source is not accelerated while that closer to the source is already accelerated close to 10³ km s⁻¹. At the highest densities, some of the gas has low to moderate velocities up to 100 km s⁻¹. Again, comparing to Fig. 4, the almost zero velocity gas is on the far side of the source and has not yet been accelerated. And finally, in the mid-density range (between 6 and 150 H cm⁻³) diagonal stripes of mass of
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Slices of the gas density (top row), temperature (second to top row), ionization fraction (second to bottom row) for a slice of a zoomed region in the L46_bigC simulation. The position of the zoomed-in region is marked by a red square in Fig. 2. To guide the eye, selected density contours are overplotted at 1 and 3 Myr. At 1 Myr, the magenta contour denotes a density of 4000 H cm⁻³. At 3 Myr, the magenta and yellow contours show densities of 4 × 10³ H cm⁻³ and 2.5 H cm⁻³, respectively. The radiation pushes the gas from the outside and disperses the outer region of the cloud. The dispersed gas moves at a speed of ~100 km s⁻¹ and is heated to a temperature of 10⁴⁻⁵ K.</figcaption></figure>

~10⁴ M☉ arise due to the dispersion of the fast-moving high-density gas to lower densities. Because of the dispersion and mixing with the background at rest, the dispersed gas also slows down building the diagonal stripes observed.

The highest velocity gas (tip of the velocity–density diagram) shows an anticorrelation with density. This corresponds to gas near the source, where the radiation from the source disperses the clouds to lower densities, which reach the highest velocities. Lower density gas is moved earlier by the outflow created by the photon–gas interaction, whereas it takes longer for high-density gas to reach the same velocities. The high-density gas eventually also reaches high velocities resulting in the whole cloud to move, which leads to the destruction of the whole disc. Comparing the velocity–density diagrams at different times, we see that at earlier times, gas at high densities has faster velocities than at later times, which is caused by several factors. First, the outflow front reaching an over-

density must decelerate while interacting with the gas from the cloud. Secondly, as we already have seen above, the high-density regions are dispersed, leading to lower densities and higher velocities as seen at 7 Myr. Finally, the flux of photons decreases with distance from the source, and late times correspond to clouds farther away from the central source, which receive a smaller flux of photons.

The mass-weighted temperature versus velocity evolution shows that the radiation accelerates the cold and dense gas that expands further into a more diffuse and broader wind. Then, an increasing amount of gas reaches higher temperatures, resulting in a high-mass fraction of high-speed hot gas at the end of the simulation. At later times, the dense and cold gas has lower velocities than at earlier times, as it corresponds to more distant clouds receiving a lower photon flux. The hot (T > 10⁵ K) gas is optically thin because it is fully ionized and the dust is destroyed in it. Therefore, the bulk of
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Slices of the IR (top row) and optical + UV (bottom row) radiation flux (over all directions) as a function of time for the same zoomed region as in Fig. 4 for the same (L46_bigC) simulation. Here the UV flux is the sum of the flux from the individual photon groups UV1, UV2, and UV3. The IR flux is computed by considering both the streaming and trapped contributions to the IR energy density. The same density contours as in Fig. 4 are overplotted to guide the eye. While the IR radiation is completely penetrating the cloud after $\le 0.1$ Myr, the column density of the cloud is too high for the optical and UV radiation to penetrate the overdensity of the cloud.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Top: mass-weighted velocity versus density for the L46_medC simulation. Bottom: mass-weighted velocity versus temperature for the L46_medC simulation. The points are coloured by the total mass within each 2D-histogram cell. The different columns show different times as labeled. The dense cold gas is accelerated by the radiation, expands and results in a broad, diffuse, and fast wind. The highest velocity for the dense gas is lower at late times, since it corresponds to more distant clouds receiving a lower photon flux. The fastest velocities are reached for gas with temperature below dust destruction.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Mass outflow rate as a function of radius for three different times 1, 3, and 7 Myr from top to bottom (same time as those used in Fig. 6) for the L46 simulations. The radiation causes the gas to move out of the galaxy, reaching mass outflow rates of up to 500–1200 M$_{\odot}$ yr$^{-1}$ after 7 Myr depending on the cloud encompassing the source. The outflow rate is larger the bigger the encompassing cloud because the radiation is trapped for longer within the big clouds and hence has longer time to impart momentum on to the gas.</figcaption></figure>

the high-velocity gas ($v > 100$ km s$^{-1}$) corresponds to intermediate values of temperature of a few $10^4$ K.

We have measured the velocities of the gas for the two other simulations (not shown here), and they have a very similar evolution. The main difference is in the intermediate density (1–100 H cm$^{-3}$) and temperature (10$^3$–10$^5$ K) range, where the gas velocity is higher with increasing cloud size around the source: 100–500, 100–600, and 200–1000 km s$^{-1}$ for the smallC, medC, and bigC simulations, respectively. A similar relationship between velocity and density of outflowing gas is seen in kinetic feedback simulations (e.g. Wagner et al. 2013; Mukherjee et al. 2016).

We now measure the mass outflow rate with

$$ \dot{M}_{\text{gas}} = \oint \rho \mathbf{v} \cdot \hat{\mathbf{r}} \, dS = \sum_{i \in \text{shell}} m_i \mathbf{v}_i \cdot \hat{\mathbf{r}}_i \frac{S}{V}, \qquad (6) $$

by considering only the outward flow (cells with $v_i \cdot \hat{r}_i > 0$) across a spherical shell of radius $r$, where $i$ denotes the index of a cell within the spherical shell of surface $S$ and volume $V$. Finally, $\hat{r}_i$ is the unit vector of the cell with velocity $v_i$. Here, we adopt a shell of thickness 0.25 kpc. Fig. 7 shows the mass outflow rate as a function of radius for the L46 simulations measured at different times. The radiation pressure causes the gas to move out, reaching outflow rates of up to 500 to 2400 M$_{\odot}$ yr$^{-1}$ depending on the initial cloud setup. The mass outflow rate is always larger the bigger the encompassing cloud around the source, due to the radiation being trapped for longer within the large clouds. The values of the mass outflow rates of 1000–2000 M$_{\odot}$ yr$^{-1}$ measured in the largest cloud simulation at $t \ge 3$ Myr as well as the high velocities of ~1000 km s$^{-1}$ are close to those measured by Tombesi et al. (2015). As already stated

above, the collapse time of the cloud encompassing the source is ~4 Myr, which hence dynamically influences the mass outflow of the galaxy, at least towards the end of the simulation. Thus, our predictions of the mass outflow rates are optimistic and have to be tested with simulations including gravity. We leave this to future work.

## 3.2 Effects of different photon groups on the cloud evolution

To better determine the specific contribution of each photon group, we have performed simulations of the same density distribution and quasar luminosity for the medium cloud size simulation (L46_medC) where we excluded certain photon groups. The density, temperature, H II fraction, and velocity maps of these simulations at 5 Myr can be seen in Fig. 8. The position of the quasar source, at the coordinate origin, is at the top-right corner of the images. In the top row, all the photon groups are included. In the middle row, the IR radiation is excluded, and in the bottom row only the IR radiation is included (i.e. the UV groups and optical are excluded).

Comparing the different rows in Fig. 8, we see that each photon group contributes to dispersing the dense gas, but the IR contribution dominates as the outflow is clearly more advanced in the IR-only run than in the run with only UV and optical. However, even if the IR photons are most important in the gas dispersion, the effect of the UV radiation is non-negligible.

Comparing the middle row with the bottom row of Fig. 8, we see that the main difference between the effects of IR and UV + optical photons is that the IR radiation plays the role of smoothing out the dense gas, especially the regions which the UV and optical radiation cannot reach. When not including the IR photons, the cloud structure is maintained with a similar multiphase state to that of the initial conditions. Generally, the UV photons only manage to disperse and ionize the gas at the shock front and do not as efficiently mix and disperse the gas of the cloud with that of the background, but instead push the overdense gas by direct radiation pressure and photoionization. The inability of UV photons in efficiently mixing the multiphase gas is most apparent when looking at the temperature (and velocity) structure.

With only the IR radiation included, the gas of the cloud is much more mixed with the background gas creating a more uniform density structure at the shock front. This arises because the IR photons are isotropically pushing the gas from the inside of the clouds and are, thus, responsible for the smoothing of the multiphase density distribution and the puffing-up of the clouds.

Comparing the temperature slices of the rows, we observe that when the UV photons are included, the temperature of the smoothed gas with densities of ~100 H cm$^{-3}$ and fully ionized is higher, due to the photoionization heating and extra momentum input. We also see that the increased velocity in the simulation incorporating all photon groups is due to the combined contribution of all these photon groups. As expected, the ionization front is more advanced when both the UV (+ optical) and IR photons are included compared to when only the IR radiation is included in the simulation.

Looking at the velocity maps from the different simulations at 5 Myr, we see that when the IR radiation is excluded, only the ionized hot gas (10$^4$–10$^{6.5}$ K) moves with large velocities ranging from 100 to 1000 km s$^{-1}$. However, with the IR photons, the neutral gas is also moving with a velocity of up to 100 km s$^{-1}$, where the gas is warm (10$^4$ K), and ~10 km s$^{-1}$, where the gas is cold (<10$^2$ K). The IR photons are, hence, capable of moving the dense, neutral gas, which the UV and optical radiation cannot reach due to the high optical depth of the cloud. We see that when all photon groups
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Evolution of the total momentum for simulations where the contribution of different photon groups are included: IR only (red lines), UV only (dark blue lines), UV and optical (light blue lines), and all groups (green lines). For the simulations, the exact same initial conditions as for the L46_medC simulation, with all the photon groups included (solid green line), are used. The dashed lines show (L<sub>group</sub>/c) t, where t is the time elapsed since the start of the simulations and L<sub>group</sub> is the luminosity in the photon group bands used in the simulation with the same colour. Through multiple scatterings on the dust, the IR radiation imparts many times a momentum L<sub>group</sub>/c on to the gas, thus greatly boosting the total momentum transferred to the gas. The main contribution to the total momentum from the UV and optical photons comes from the reprocessed UV photons into IR photons that then multiscatter and impart a momentum boost on to the gas. Photoionization heating has a small but non-negligible effect. Finally, the optical photons give a small contribution to the momentum.</figcaption></figure>

of the simulation, and Δt<sub>N</sub> the time interval between two snapshots. As shown in Fig. 1, ~80 per cent of the total bolometric luminosity is covered by the photon groups used in the simulations (recall that we do not cover the hard X-ray band). A mechanical advantage above unity occurs when the photons boost the momentum transfer between the radiation and the gas. This can happen when, the IR photons are multiply scattered, as well as by photoionization heating and subsequent shock formation. The reason we used the bolometric momentum in this calculation is because it is consistent with subgrid models of BH feedback in the literature (e.g. Debuhr et al. 2011; Barai et al. 2014; Choi et al. 2014; Costa et al. 2014; Zubovas & Nayakshin 2014; Hopkins et al. 2016).

Fig. 10 displays the evolution of the mechanical advantage for the L46_smallC, L46_medC, and L46_bigC simulations. For all the simulations, the mechanical advantage decreases with time. However, the magnitude of the mechanical advantage is not the same for the three different cloud size simulations, in particular in the early stages. The efficiency of the momentum transfer for the L46_smallC, L46_medC, and L46_bigC simulations is bigger, the larger the cloud encompassing the source of radiation. Because it takes more time to destroy larger clouds, the photons are trapped and scatter for a longer time. Indeed, a bigger encompassing region around the quasar results in larger the momentum boost from the radiation. For all the simulations, the radiation carves, right from the beginning, a hole into the centre of the galaxy, through which the radiation can escape (Fig. 2). This causes the mechanical advantage to decrease with time as the photons injected by the source are more likely to escape the system without scattering. We see that the optical depth (or cloud size) plays a very important factor in the evolution of the mechanical advantage. In addition, dust destruction in hot gas can play a role in suppressing the amount of momentum passed from radiation to gas. We have seen in Section 3.1.2, and especially

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Evolution of the mechanical advantage (momentum input rate *p* over L/c, where L is bolometric luminosity) as a function of time for the L46 simulations. The mechanical advantage is above unity until 10 Myr, it is larger the bigger the encompassing cloud, and it decreases with time. For the medC and bigC simulations, the mechanical advantage starts decreasing before the central cloud is fully destroyed as the efficiency of the momentum transfer is already less efficient once the radiation manages to carve a hole into the centre of the galaxy.</figcaption></figure>

in Fig. 3, that the gas is quickly heated to high temperatures for the L46_medC and L46_bigC simulations whereas it takes longer for the L46_smallC gas to reach similarly high temperatures. It leads us to the conclusion that dust destruction has more influence for the bigger cloud simulations compared to the L46_smallC simulation. However, as we will see in Section 3.4, this effect is not as important as the effect of the cloud size and gas density surrounding the source.

## 3.4 Evolution of the optical depth

For a better understanding of the mechanical advantage from the radiation and the efficiency of photon-gas coupling, we measure the optical depth through the disc for the quasar IR radiation as a function of time for the L46_smallC, L46_medC, and L46_bigC simulations. The IR optical depth is defined as

$$ \tau_{IR} = \int \rho \kappa_{IR} dl = \sum_{i \in \text{LOS}} \rho_i \kappa_{IR}(T_i, Z_i) \Delta l_i, \quad (7) $$

where the opacity κ<sub>IR</sub>, function of temperature and metallicity, is given in equation (5), *l* is the line-of-sight (LOS) coordinate, and where ρ<sub>i</sub>, T<sub>i</sub>, and Z<sub>i</sub> are, respectively, the density, temperature, and metallicity of the cell of index *i* along the LOS, while Δl<sub>i</sub> is the length of the LOS through the *i*th cell.

Fig. 11 shows the evolution of the mean optical depth calculated over 500 randomly selected LOS, uniformly sampling a sphere from the centre of the disc, where the quasar source is located, up to a distance of 1.5 kpc.

At the beginning of the simulation, the optical depth calculated over the central cloud is ~0.8 times that of the optical depth calculated over the whole disc for L46_medC and L46_bigC, while it is 0.1 that of the disc for L46_smallC. Once the radiation carves a hole into the central cloud, the optical depth calculated over the cloud drops to zero. This happens later for bigger encompassing clouds. Note that at 10 Myr, the mean optical depths from the source for the three simulations converge to the same value of τ<sub>IR</sub> ≃ 4.

As discussed above, there are two important phases for the IR radiation. First, the radiation is trapped within an optically thick region of the central cloud and imparts momentum on to the gas,
<figure><img src="image_12.png" /><figcaption>Figure 11. Left: evolution of the optical depth $\tau_{IR}$ as a function of time for the L46_smallC, L46_medC, and L46_bigC simulations. The lines show the mean values of the optical depth calculated by sampling the sphere along 500 different lines of sight, over the galaxy disc (up to 1.5 kpc, solid lines) or over the central cloud using a radius corresponding to the largest fractal structure within the respective simulation box (see Table 2, dashed lines). The shaded areas show the $\pm 1 \sigma$ distribution of the optical depth over the galaxy. The mean optical depth of the three simulations depends on the cloud size and is larger the bigger the clouds within the disc. Right: fraction of solid angle over the sphere covered by an optical depth greater than a threshold optical depth $\tau_{thres}$ of 0.1 (solid), 1 (dashed), or 10 (dotted). The optical depth is calculated over the whole sphere. The figure shows that with increasingly bigger clouds, the covering fraction of $\tau_{IR} > 1$ is enhanced at given times and decreases at later times.</figcaption></figure>

which dominates the evolution of the outflow. Once the radiation has destroyed the central cloud, the IR radiation is only trapped within the densest regions of the gas. We have seen that the mean density of the cloud decreases with time, which, in turn, reduces the optical depth of the gas and hence reduces the influence of the trapped photons at later stages. Therefore, in this second phase, the IR radiation plays a less dominant role in pushing the gas out of the disc. These two phases are clear in Fig. 11: the cloud starts are high optical depth with little decrease, and then there is a rapid drop.

The optical depths displayed in the left-hand panel of Fig. 11 show considerable scatter, so that much radiation can escape even when the mean optical depth is high. The right-hand panel of Fig. 11, the fraction of the solid angle covered with $\tau_{IR}$ larger than a threshold optical depth $\tau_{thres}$ is shown. The fraction of the solid angle covered with $\tau_{IR} > 1$ gives an indication of the trapping of the IR photons, i.e. their multiple scattering, which in turns gives an indication of the efficiency of photon to gas coupling.

The general evolution of the optical depth distribution for the L46_smallC, L46_medC, and L46_bigC simulation helps to understand the evolution of the mechanical advantage. At the start of the simulations, before the central cloud is destroyed, the optical depth distributions calculated over the cloud for the L46_medC and L46_bigC simulations overlap, with mean optical depths of $\langle \tau_{IR} \rangle = 75$ and 56 for the L46_bigC and L46_medC simulations, respectively (left-hand panel). However, the L46_bigC simulation has more scatter in the optical depth distribution compared to that of the L46_medC simulation. Additionally, the two simulations have the same unity fraction of solid angle over the sphere covered by $\tau_{IR}$ larger than unity for the first ~1 Myr (right-hand panel). The larger optical depth for the L46_bigC simulation compared with the L46_medC simulation leads to a larger mechanical advantage at the start of the simulation. The mean optical depth of the two simulations decreases with time and they converge after ~10 Myr. The mean optical depth at the beginning of the L46_smallC simulation, on the other hand, is much smaller ($\langle \tau_{IR} \rangle \sim 11$). Moreover, the fraction of solid angle around the quasar with $\tau_{IR} > 1$ is at 80 per cent from the start for L46_smallC. Thus, it explains the small mechanical advantage for the small cloud simulation when

compared to the larger cloud simulations since a significant fraction (20 per cent) of all possible lines of sights are optically thin.

Indeed, the optical depth around the source is the important factor in understanding how much momentum from the photons can be transferred to the gas. However, as we have seen, channels of optically thin gas can form within the optically thick layers of gas, reducing the expected amount of momentum that is effectively transferred to the gas. Following Hopkins et al. (2011), we write the expected momentum boost as

$$ \dot{p} = (1 + \eta \tau_{IR}) \frac{L}{c}, \qquad (8) $$

where the factor of $\eta \tau_{IR} L/c$ accounts for the momentum boost passed on to the gas by the total number of IR scattering events. Equation (8) includes a dimensionless ad hoc reduction factor $\eta$ that accounts for extra sources of momentum (e.g. UV photoheating, X-ray Compton scattering off the electrons, dust photoelectric heating increasing $\eta$ above unity) and inhomogeneities in the gas (decreasing $\eta$ below unity). Hence, the reduction factor accounts for the fraction of expected multiple scattering that effectively happen. Note that Hopkins et al. (2011, and others) also use the mean optical depth when defining $\eta$.

Given the change in momentum estimated from the simulation, we can determine the reduction factor $\eta$ from equation (8), and this is shown in Fig. 12. The solid lines show the reduction factor calculated with the mean optical depth within the galaxy ($\eta_{gal}$), whereas the dashed line show the reduction factor calculated using the mean optical depth over the central cloud ($\eta_{cloud}$).

For all the three simulations, $\eta_{gal}$ starts below unity starting at ~0.2 for the L46_smallC and L46_medC simulation and ~0.3 for the L46_bigC simulation and then decreases slowly with time. On the other hand, $\eta_{cloud}$ starts with a similar value as that of $\eta_{galaxy}$ for the L46_medC and L46_bigC simulations, whereas $\eta_{cloud}$ for the L46_smallC simulation already starts at a higher value due to the reduced optical depth over the cloud for this simulation. After ~1 Myr, $\eta_{cloud}$ rises steeply in the L46_medC simulation due to the optical depth dropping to zero around that time. The two reduction factors for the L46_bigC simulation evolve very similarly, again
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Evolution of the reduction factor $\eta$, providing a measure of the coupling efficiency of IR photons (see equation 8) as a function of time for the L46 simulations. The solid and dashed lines show the reduction factor calculated with the mean optical depth within the galaxy ($\eta_{\text{gal}}$) and within the central cloud ($\eta_{\text{cloud}}$). The values of the reduction factor become unphysical when $\tau_{\text{IR}}$ falls below values around unity (cloud destruction).</figcaption></figure>

due to the similar behaviour of the optical depth calculated either over the whole galaxy or the encompassing cloud.

It is important to keep in mind that equation (8) does not hold when the IR radiation is not trapped anymore and escapes without multiscattering, i.e. when $\langle \tau_{\text{IR}} \rangle < 1$, which explains why the value of $\eta$ diverges in this regime. This is mostly the case for $\eta_{\text{cloud}}$ when the central cloud is destroyed. Our measurements show that, for the cases where $\langle \tau_{\text{IR}} \rangle \ge 1$, the measured mechanical advantage is well below unity even when the source is placed in a region surrounded by a large optical depth. It shows that the non-uniform structure of the ISM and the subsequent building of low-density channels as well as the building of a hole in the centre of the galaxy has a great influence in setting the reduction factor.

This is in line with the argumentation by Krumholz & Matzner (2009) and Krumholz & Thompson (2012, 2013), who use analytical arguments to show that the effective optical depth can never be larger than a few, due to the tendency of photons to escape through lower density channels, even though if the average optical depth

is much larger than unity. Their findings have been challenged by Davis et al. (2014), who showed that the results highly depend on the solver used for the radiative transfer. The flux-limited diffusion (FLD) method used in Krumholz & Thompson (2012) and Krumholz & Thompson (2013) leads to a significantly lower efficiency in accelerating the gas via radiation pressure than in the more accurate method (variable Eddington tensor; VET) used by Davis et al. that showed a stronger acceleration of the gas. Rosdahl & Teyssier (2015) showed that the M1 method used in our simulations lies somewhat between those of FLD and VET. Using an Implicit Monte Carlo radiation transfer scheme, Tsang & Milosavljević (2015) find results consistent with Davis et al., demonstrating again the importance of accurate radiative transfer in simulations of radiative feedback.

We have already seen that the initial stage, when the radiation is trapped within the central cloud, plays a crucial role in the early acceleration of the gas. We thus expect the IR radiation to be an important driver of the outflow at early times. However, the large IR luminosity rapidly destroys the cloud encompassing the source, which leads to a decrease in the optical depth. This, in turn, shortens the time the radiation is trapped within the central cloud, where it can scatter sufficiently to impart a larger momentum boost on to the gas.

## 3.5 Effects of the quasar position

As radiatively driven AGN winds are stronger when the source is embedded within more massive clouds, we expect that the location of the source relative to the cloud should also influence the momentum given to the gas. We have re-simulated the L46_bigC simulation, changing the underlying density at the position of the quasar from a high-density environment (L46_bigC simulation shown before, hereby referred to as $\text{max}\rho_{\text{Q}}$ simulation) to a medium-density region ($\text{med}\rho_{\text{Q}}$) and to a low-density environment ($\text{min}\rho_{\text{Q}}$). The regions are chosen such that the average density, calculated over a region including the direct neighbouring cells, is maximum (~8000 H cm⁻³), around the mean (~320 H cm⁻³), and minimum (~0.7 H cm⁻³), respectively. The initial density distribution of the three different simulations is shown in Fig. 13.

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Initial density distribution of the different environments for the quasar with L46_bigC_minρQ (left-hand panel), L46_bigC_medρQ (middle panel), and L46_bigC_maxρQ (right-hand panel). The quasar position is marked by a red cross at the centre of each image.</figcaption></figure>
The top panel of Fig. 14 displays the evolution of the mechanical advantage for the three different quasar positions. This figure indicates that the maximum mechanical advantage reached in the simulation increases by up to a factor of 10 when the gas density at the position of the quasar is increased from ~0.7 H cm⁻³ to ~8000 H cm⁻³.

The middle panel of Fig. 14 shows that the quasar position has an important effect on the mean IR optical depth, which increases by a factor of ~30 with increasing density at the quasar location. According to the lower panel of Fig. 14, at the start of the simulation, the fractions of solid angle covered with $\tau_{IR} > 1$ increases from 40 per cent for the low-gas density around the source to 100 per cent for the run with the high gas density at the source location.

When the quasar is placed within the minimum density environment, the mean optical depth is initially at $\langle\tau_{IR}\rangle \sim 3$, but the solid angle covered by an optical depth larger than $\tau_{IR} = 1$ is only ~30 per cent, and ~60 per cent of the IR photons have $\tau_{IR} < 0.1$. A significant fraction of the photons can hence free-stream out of the disc via optically thin channels and multiscattering is negligible. The mechanical advantage is at roughly unity for the first few 0.1 Myr before dropping off. Thus, there is no overall momentum boost in the beginning, but since a fraction of the radiation escapes freely without interactions, another fraction of the radiation must give a boosted momentum to the gas.

With the quasar placed at an intermediate density, ~80 per cent of the solid angle around the source initially has $\tau_{IR} > 1$, but a non-negligible fraction of ~15 per cent of lines of sight have $\tau_{IR} < 0.1$ and hence are more or less free-streaming channels. In this setup, the mechanical advantage is intermediate between the two other cases, starting at ~3, steadily dropping, and going below unity after about 1 Myr. Thus, the probability of the photons scattering more than once is much smaller when the local density around the quasar is intermediate compared to when it has the largest value. The environment of the source determines the optical depth around the quasar, which influences the amount of momentum the photons can impart on to the gas. If the source is already in an environment where the covering fraction for $\tau_{IR} > 1$ is small, the mechanical advantage is much smaller than when the source is in an environment with a large optical depth and the covering fraction for $\tau_{IR} > 1$ is higher.

## 3.6 Comparison between different luminosities

We now study the effect of the quasar luminosity on the mechanical advantage. We compare four simulations, using two different quasar luminosities, $L = 10^{43}$ and $10^{46}$ erg s⁻¹ (the L43 and L46 simulations, respectively), and two different quasar locations: intermediate and high gas density around the source, for the largest cloud size simulation.

The top panel of Fig. 15 shows the mechanical advantage of the four simulations as a function of time. For a given density at the quasar location, the mechanical advantage at early times (0.1 Myr) has similar values for the two different luminosities. However, in the higher luminosity case it peaks and drops sooner than in the lower luminosity case.

The reason for the higher mechanical advantage in the low-luminosity simulations is that the time it takes for the L46 source to disperse the encompassing cloud and to carve a hole into the centre of the disc is shorter than for the L43 source. This is clear when comparing the evolution of the optical depth for low and high luminosity (middle panel of Fig. 15). The optical depth decreases faster with the high-luminosity quasar than with the lower luminosity one. In-

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Top: evolution of the mechanical advantage for the L46_bigC simulations, where the position of the quasar within the density field has been varied. The labels max, mean, and min stand for the quasar position within the environment of maximum, mean and minimum density, respectively. Middle: evolution of the mean galaxy optical depth $\tau_{IR}$ with the one $\pm\sigma$ standard deviation (shaded areas). Bottom: fractions of solid angle over the sphere covered with $\tau_{IR}$ larger than a threshold optical depth $\tau_{thres}$. We see that the environment of the luminous source has a strong effect on the maximum boost gained by the photons.</figcaption></figure>

deed, with the lower luminosity quasar, the radiation has thus more time during which the surrounding cloud is intact and the radiation multiscatters inside it. On the other hand, with a more luminous source, the surrounding cloud is more quickly dispersed, allowing the photons to escape, and giving a lower mechanical advantage.

It is important to stress that the mechanical advantage is defined as the fraction of the momentum change over the total bolo-metric momentum, where the use of the bolometric momentum is consistent with subgrid models of BH feedback in the literature (e.g. Costa et al. 2014; Zubovas & Nayakshin 2014; Hopkins et al. 2016). Hence, even though the mechanical advantage is higher for the lower luminosity source, the total momentum at ~1 Myr for the high-luminosity source is $82 M_{\odot} \text{ km s}^{-1} (7.8 M_{\odot} \text{ km s}^{-1})$ and
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Top: evolution of the mechanical advantage for low and high quasar luminosities, respectively, $L = 10^{43}$ erg s$^{-1}$ (L43) and $L = 10^{46}$ erg s$^{-1}$ (L46), and two different local gas densities (see Fig. 13). Middle: evolution of the mean IR optical depth with the $\pm 1\sigma$ standard deviation (shaded areas for L46 simulations, hatched areas for L43 simulations). Bottom: evolution of the reduction factor $\eta$. The mechanical advantage for the lower luminosity simulation is higher than that of the higher luminosity simulation with a larger value of $\eta$. The higher luminosity leads to a faster destruction of the central cloud compared to the lower luminosity, leaving less time for the radiation to impart its momentum to the gas.</figcaption></figure>

with this ~210 (~360) times higher than the total momentum in the low-luminosity quasar simulation, with the quasar placed in both simulations in a maximum (medium) density environment.

The reduction factor $\eta$ starts between 0.1 and 0.3 for the simulations with the bright quasar and is roughly twice as great in the simulations with the faint quasar (bottom panel Fig. 15). The higher reduction factor for lower luminosities is not only caused by the higher mechanical advantage for the low-luminosity simulation, but is also a consequence of the different distributions of the optical depth. The mean and lower error bound of the optical depth distribution of the high quasar luminosity simulation reaches much lower values than the corresponding values in the low quasar luminosity run. Hence, while the high quasar luminosity run produces a rapidly evolving outflow, with low-density channels that arise, in the low-

luminosity run the gas has more time to mix with the surrounding gas, leading to a smoother outflow, with no low-density channels arising. Therefore, a larger fraction of photons manage to escape through lower density channels for the higher luminosity simulation that then in turn lowers the reduction factor.

In summary, we find that the mechanical advantage is higher for lower luminosities because the radiation manages to destroy the encompassing cloud faster if the source is more luminous. The reduction factor starts in a similar range (0.1–0.3) for the two different luminosities, but then reaches highest values for the simulation with low luminosity and quasar placed in a maximum density environment. The reason for the higher reduction factor for lower luminosities is that in the high-luminosity run the radiation causes the formation of low-density tunnels that are not formed in the lower luminosity simulation due to the slower propagation of the wind. A larger amount of low-density tunnels lead to a smaller efficiency of the momentum transfer between the radiation and the gas and hence a smaller reduction factor for the higher luminosity simulation. However, the amount of momentum gained still scales with the luminosity of the source.

## 4 DISCUSSION

While our simulations are well suited for studying the interactions between radiation and gas, they come with caveats. We now discuss the shortcomings of our simulations and possible implications and changes when taking them into account.

First of all, including gravity in our simulations could have different, possibly opposing, effects. The first effect is that gravity may slow down the propagation of the outflows, as the gravitational potential will make it harder for the gas to escape the galaxy, which in turn would decrease the mass outflow rate out of the galaxy. However, part of the outflow reaches velocities up to $1000 \text{ km s}^{-1}$, which is fast enough for the gas to escape the central encompassing cloud and even escape out of the massive halo. The escape velocity out of the largest encompassing cloud of mean density of $\sim 100 \text{ H cm}^{-3}$ and radius $\sim 1.5 \text{ kpc}$ is $\sim 450 \text{ km s}^{-1}$. The cloud collapse time, on the other hand, is $\sim 9 \text{ Myr}$. Since the outflow velocities of $1000 \text{ km s}^{-1}$ are greater than the escape velocity out of the cloud and are reached on a faster time-scale ($\le 5 \text{ Myr}$) than the free-fall time of the central cloud, it is reasonable to assume that the central cloud indeed gets destroyed and that the outflow progresses further out of the galaxy.

Within the disc, the influence of gravity is less clear. In our simulation, radiation is capable of carving a hole in the central regions of the galaxy. In reality, part of the gas displaced by the radiation (that displaced vertically above the now rotating disc) would likely fall back to the centre, which would enhance the optical depth around the BH and thus could help to boost the momentum transfer from the radiation for a longer time than observed in our simulations.

With the inclusion of radiative cooling, in addition to gravity, the collapse of massive clouds would be enhanced, leading to a larger cloud mass of smaller size. This would likely strengthen the matter–radiation coupling due to the larger optical depths of the clouds. However, radiative cooling (in combination with gravity) could also lead to a fragmentation of the massive clouds (by reducing the Jeans length) and thus speed up the formation of lower density channels through which the radiation would escape.

In addition, our simulations explore the effect of a steady, uniform quasar, located at the centre of the galaxy. Actual BHs radiate with a luminosity that is proportional to their accretion rate, and hence their luminosities are not steady. Once the radiation manages to
carve a hole in the central regions of the galaxy, the luminosity of the quasar should drop to lower values hence decreasing the amount of momentum transferred to the gas.

The interplay between the effect of gravity and the changing luminosity of the quasar naturally leads to a self-regulating feedback cycle, where AGN activity pushes the gas out of the central regions of the galaxy leading to starvation of the BH, shutdown of radiative emission, and subsequent fall-back of the gas on to the BH. Furthermore, for gas-rich galaxies, the ISM gas is multiphase and filamentary, with dense star-forming clouds that can reach and feed the BH, which certainly leads to bursty and chaotic accretion (see e.g. Ciotti & Ostriker 2012, for a review). For groups and clusters, Gaspari, Ruszkowski & Oh (2013) argue that the combination of cooling and thermal instability, induced by BH feedback, causes the cold gas to collapse and decouple from the hot gas, leading to cold and chaotic BH accretion. In addition, simulations have shown that anisotropic thermal conduction enhances the formation of cold filamentary gas along magnetic field lines (e.g. Lecoanet, Parrish & Quataert 2012; Sharma et al. 2012). Nevertheless, since conductivity scales steeply with temperature (to the power 2.5), this mechanism that may be relevant to the atmosphere of galaxy clusters likely plays a minor role, if any, at the lower temperatures of the ISM.

The inclusion of gravity, cooling, (chaotic) accretion, the time variability of the quasar luminosity, as well as thermal conduction, increases the complexity of the non-linear interplay between the radiation and the gas, and we defer this study to future work.

Finally, it is still under debate whether an outflow driven by radiation from the BH may also lead to a (local) enhancement of star formation due to the compression of the clouds and the subsequent formation of more stars. We have neglected the formation of stars and stellar feedback in our simulations. Stars generally form in dense, cold gas regions. Removing gas mass from the galaxy is thought to negatively impact the formation of stars. As deduced in observations (e.g. Cresci et al. 2015), it may also be plausible that the compression of the gas during a burst of quasar activity also triggers star formation. Such a positive feedback effect has been shown for the non-radiative modes of AGN feedback (Gaibler et al. 2012; Zubovas et al. 2013; Bieri et al. 2015, 2016) and has been discussed for the radiative mode of AGN feedback by Ciotti & Ostriker (2007), who argue that AGN feedback can in fact trigger star formation in the central regions of galaxies. We defer to future work a deeper discussion of the possibility of triggered star formation due to quasar feedback.

## 5 CONCLUSIONS

Quasar-driven winds are powered by complex interactions between radiation and gas. In most recent state-of-the-art HD cosmological simulations (but also in HD simulations of isolated disc galaxies), quasar feedback is approximated by depositing thermal energy within the resolution element, where the efficiency of the radiation–gas coupling is represented by a single parameter chosen to match global observations such as the BH mass–bulge velocity dispersion (Ferrarese & Merritt 2000). Additionally, there is no consensus from these simulations on whether these AGN winds are momentum conserving or energy conserving, although recent observations favour energy-conserving winds from the nuclear accretion disc (Feruglio et al. 2015; Tombesi et al. 2015). Generally we refer to the outflow as energy conserving if the radiation is capable of efficiently boosting the momentum transfer between the photons and the gas (i.e. $\dot{p} \gg L/c$), and to momentum conserving if the amount of mo-

mentum transferred to the gas is similar to the momentum flux of photons provided by the source (i.e. $\dot{p} \approx L/c$).

Given the likely importance of AGN feedback in the evolution of massive galaxies and the increasing amount of radio galaxy observations – e.g. from the Low-Frequency Array, the Australian Square Kilometer Array Pathfinder, as well the Atacama Large Millimeter Array, which promise to provide us a better view of the ISM properties in these galaxies – it is time to improve our theoretical understanding of the mechanisms that drive the momentum and energy transfer from the photons to the ISM in order to properly quantify the role of AGN feedback in the evolution of galaxies.

We have performed idealized galactic disc RHD simulations to study the coupling of photons with the highly multiphase galactic gas using different cloud sizes. In the simulations, the emission, absorption, and propagation of photons and their interaction with the gas via photoionization, momentum transfer, and absorption/scattering on dust, is followed self-consistently.

We find that radiation from a bright quasar is capable of driving a powerful wind, with multiscattering IR photons playing the dominant role, especially at early times. The UV and optical photons never manage to reach the high-density regions of clouds due to the large optical depths. Therefore, the UV and optical photons exert their direct push to the clouds from outside. Some of the UV and optical radiation is dust absorbed and reprocessed into IR radiation that can propagate much deeper into the dense gas. This reprocessed IR radiation is the most significant contribution of the higher energy radiation to the total momentum of the gas (see Fig. 9).

The mass outflow rates as well as the velocity reached by the galactic wind depend on the structure of the ISM (and luminosity of the source). Larger clouds have greater IR optical depths, hence longer time during which IR photons are trapped within and scatter multiple times. This first phase plays a crucial role in the early acceleration of the gas. Once the radiation has destroyed the central encompassing cloud, the IR radiation can stream out of the galaxy through low-density regions, decreasing the efficiency of the momentum transfer from the photons to the gas.

Both the outflow rates of ~500 to 1000 $M_{\odot}$ yr$^{-1}$ and high velocities of ~1000 km s$^{-1}$ that we measured in our simulations with a large encompassing cloud around the quasar are close to those observed by Tombesi et al. (2015) and Feruglio et al. (2015). This agreement favours winds that appear as energy driven, which in our case are actually radiatively driven. But as mentioned above, our measured outflow rates are optimistic predictions as we did not include gravity in our simulations.

The mechanical advantage, defined as the ratio of the imparted momentum rate $\dot{p}$ and $L/c$, decreases with time, from ~3 to 30 at the beginning of the simulation to unity after ~10 Myr, thanks to the decreasing efficiency of the photon–gas coupling. This mechanical advantage largely depends on the size of the encompassing cloud as well as the position of the quasar (and the quasar luminosity). It varies by a factor of 10, depending on the size of the encompassing cloud (50 pc to 1.5 kpc) used in our simulations, with the IR optical depths ranging from $\tau_{IR} = 10$ to 100. The position of the quasar plays an important role in setting the amount of momentum transferred from the photons to the gas, where the mechanical advantage changes by a factor of 50, depending on whether the source is buried within the cloud or illuminating it from the outside. Thus, the IR optical depth $\tau_{IR}$ (i.e. the cloud size and mass) around the quasar is the most important criterion in determining the momentum injection rate at a given quasar luminosity.

The reduction factor $\eta = [\dot{p}/(L/c) - 1]/\tau_{IR}$ is an empirical estimate that accounts for extra sources of momentum (e.g. UV
photoionization heating, X-ray Compton scattering off the elec-
trons, dust photoelectric heating, etc.) and inhomogeneities in the
gas. The measured reduction factor never reaches values of unity,
showing that the mechanical advantage never reaches values as high
as τIR. The non-uniform structure of the ISM and the formation
of low-density channels is responsible in setting this low value of
η = 0.2–0.3, that decreases significantly once the most central cloud
starts breaking up and becomes optically thin. Thus, the number of
scattering events of IR photons is roughly one quarter of the IR
optical depth.

The radiation emitted by the central quasar destroys the cloud encompassing the source, which leads to a rapid decrease in the optical depth and reduction factor. We demonstrated that the destruction time is shorter for larger luminosities, which leads to a smaller mechanical advantage. Additionally, the fast evolution of the wind in the high-luminosity simulation leads to the formation of low-density channels leading to a smaller efficiency of the radiation–gas coupling. For the low-luminosity simulation, on the other hand, the evolution is slower which leaves the gas more time with to mix with the surrounding gas and thus leads to a smaller amount of optically thin gas. This in turn leads to a higher reduction factor for the lower luminosity simulation. Nonetheless, large luminosities are required to obtain radiatively driven quasar winds with fast velocities and high-mass outflow rates.

In the future, this study will be extended to AGN radiation within
self-regulated turbulent discs including gravity, gas cooling, star
formation, stellar feedback, and a model for dust creation and
destruction.

ACKNOWLEDGEMENTS

It is our pleasure to thank Romain Teyssier, Marta Volonteri, Pierre
Guillard, and Jonathan Coles for valuable discussions. We thank
the referee, Luca Ciotti, for his constructive comments and sug-
gestions that improved this paper. This work was granted access to
the HPC resources of CINES under the allocations c2015047421
made by GENCI. RB has been supported by the Institute Lagrange
de Paris. JR was funded by the European Research Council under
the European Union's Seventh Framework Programme (FP7/2007-
2013)/ERC Grant agreement 278594-GasAroundGalaxies, and the
Marie Curie Training Network CosmoComp (PITN-GA-2009-
238356). JS acknowledges support from project 267117 (DARK)
hosted by UPMC – Sorbonne Universités, from JHU by National
Science Foundation grant OIA-1124403. This work has been par-
tially supported by grant Spin(e) ANR-13-BS05-0005 of the French
ANR. This work has made use of the Horizon cluster, hosted by the
Institut d'Astrophysique de Paris. We warmly thank S. Rouberol
for running it smoothly.

REFERENCES

Agertz O., Lake G., Teyssier R., Moore B., Mayer L., Romeo A. B., 2009, MNRAS, 392, 294

Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, ApJ, 770, 25

Antonuccio-Delegu V., Silk J., 2010, MNRAS, 405, 1303

Barai P., Viel M., Murante G., Gaspari M., Borgani S., 2014, MNRAS, 437, 1456

Bicknell G. V., Sutherland R. S., van Breugel W. J. M., Dopita M. A., Dey
A., Miley G. K., 2000, ApJ, 540, 678

Bieri R., Dubois Y., Silk J., Mamon G. A., 2015, ApJ, 812, L36

Bieri R., Dubois Y., Silk J., Mamon G. A., Gaibler V., 2016, MNRAS, 455, 4166

Booth C. M., Schaye J., 2009, MNRAS, 398, 53

Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S., Baugh C.
M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645

Bower R., Schaye J., Frenk C. S., Theuns T., Schaller M., Crain R. A., McAlpine S., 2016, preprint (arXiv:e-prints)

Capelo P. R., Natarajan P., Coppi P. S., 2010, MNRAS, 407, 1148

Choi E., Ostriker J. P., Naab T., Johansson P. H., 2012, ApJ, 754, 125

Choi E., Naab T., Ostriker J. P., Johansson P. H., Moster B. P., 2014, MNRAS, 442, 440

Cicone C. et al., 2014, A&A, 562, A21

Ciotti L., Ostriker J. P., 2007, ApJ, 665, 1038

Ciotti L., Ostriker J. P., 2012, in Kim D.-W., Pellegrini S., eds, Astrophysics and Space Science Library, Vol. 378, Hot Interstellar Matter in Elliptical Galaxies. Springer-Verlag, Berlin, p. 83

Costa T., Sijacki D., Haehnelt M. G., 2014, MNRAS, 444, 2355

Cresci G. et al., 2015, A&A, 582, A63

Croton D. J. et al., 2006, MNRAS, 365, 11

Daddi E. et al., 2010, ApJ, 713, 686

Davis S. W., Jiang Y.-F., Stone J. M., Murray N., 2014, ApJ, 796, 107

Debuhr J., Quataert E., Ma C.-P., Hopkins P., 2010, MNRAS, 406, L55

Debuhr J., Quataert E., Ma C.-P., 2011, MNRAS, 412, 1341

Dekel A., Silk J., 1986, ApJ, 303, 39

Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604

Di Matteo T., Colberg J., Springel V., Hernquist L., Sijacki D., 2008, ApJ,
676, 33

Draine B. T., Salpeter E. E., 1979, ApJ, 231, 77

Dubois Y., Teyssier R., 2008, A&A, 477, 79

Dubois Y., Devriendt J., Slyz A., Teyssier R., 2012, MNRAS, 420, 2662

Dubois Y. et al., 2014, MNRAS, 444, 1453

Dubois Y., Volonteri M., Silk J., Devriendt J., Slyz A., Teyssier R., 2015, MNRAS, 452, 1502

Faucher-Giguère C.-A., Quataert E., 2012, MNRAS, 425, 605

Ferrarese L., Merritt D., 2000, ApJ, 539, L9

Feruglio C. et al., 2015, A&A, 583, A99

Fischera J., Dopita M. A., 2004, ApJ, 611, 919

Fischera J., Dopita M. A., Sutherland R. S., 2003, ApJ, 599, L21

Gabor J. M., Bournaud F., 2014, MNRAS, 441, 1615

Gaibler V., Khochfar S., Krause M., Silk J., 2012, MNRAS, 425, 438

Gaspari M., Ruszkowski M., Oh S. P., 2013, MNRAS, 432, 3401

Genzel R. et al., 2010, MNRAS, 407, 2091

Gnedin N. Y., Abel T., 2001, New Astron., 6, 437

Habouzit M., Volonteri M., Dubois Y., 2016, preprint (arXiv:e-prints)

Hopkins P. F., Richards G. T., Hernquist L., 2007, ApJ, 654, 731

Hopkins P. F., Quataert E., Murray N., 2011, MNRAS, 417, 950

Hopkins P. F., Torrey P., Faucher-Giguère C.-A., Quataert E., Murray N., 2016, MNRAS, 458, 816

Hu J., 2008, MNRAS, 386, 2242

Jiang Y.-F., Stone J. M., Davis S. W., 2012, ApJS, 199, 14

Joung M. K. R., Mac Low M.-M., 2006, ApJ, 653, 1266

Khandai N., Di Matteo T., Croft R., Wilkins S., Feng Y., Tucker E., DeGraf
C., Liu M.-S., 2015, MNRAS, 450, 1349

Kim J.-h., Wise J. H., Alvarez M. A., Abel T., 2011, ApJ, 738, 54

King A., 2003, ApJ, 596, L27

Kormendy J., Ho L. C., 2013, ARA&A, 51, 511

Kormendy J., Bender R., Cornell M. E., 2011, Nature, 469, 374

Kritsuk A. G., Norman M. L., Wagner R., 2011, ApJ, 727, L20

Krumholz M. R., Matzner C. D., 2009, ApJ, 703, 1352

Krumholz M. R., Thompson T. A., 2012, ApJ, 760, 155

Krumholz M. R., Thompson T. A., 2013, MNRAS, 434, 2329

Krumholz M. R., Klein R. I., McKee C. F., 2011, ApJ, 740, 74

Lecoanet D., Parrish I. J., Quataert E., 2012, MNRAS, 423, 1866

Levermore C. D., 1984, J. Quant. Spectrosc. Radiat. Transfer, 31, 149

Lewis G. M., Austin P. H., 2002, Proc. of the 11th Conf. on Atmospheric Radiation, An Iterative Method for Generating Scaling Log-normal Simulations

Li A., Draine B. T., 2001, ApJ, 554, 778

Mac Low M.-M., Klessen R. S., 2004, Rev. Mod. Phys, 76, 125

McKee C. F., Ostriker E. C., 2007, ARA&A, 45, 565
Magorrian J. et al., 1998, AJ, 115, 2285

Mukherjee D., Bicknell G. V., Sutherland R., Wagner A., 2016, MNRAS, 461, 967

Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563

Nayakshin S., Zubovas K., 2012, MNRAS, 427, 372

Novak G. S., Ostriker J. P., Ciotti L., 2012, MNRAS, 427, 2734

Oppenheimer B. D., Davé R., 2006, MNRAS, 373, 1265

Ostriker J. P., Choi E., Ciotti L., Novak G. S., Proga D., 2010, ApJ, 722, 642

Pawlik A. H., Schaye J., 2011, MNRAS, 412, 1943

Petkova M., Springel V., 2009, MNRAS, 396, 1383

Proga D., Stone J. M., Kallman T. R., 2000, ApJ, 543, 686

Roos O., Juneau S., Bournaud F., Gabor J. M., 2015, ApJ, 800, 19

Rosdahl J., Teyssier R., 2015, MNRAS, 449, 4380

Rosdahl J., Blaizot J., Aubert D., Stranex T., Teyssier R., 2013, MNRAS, 436, 2188

Roškar R., Teyssier R., Agertz O., Wetzstein M., Moore B., 2014, MNRAS, 444, 2837

Sazonov S. Y., Ostriker J. P., Sunyaev R. A., 2004, MNRAS, 347, 144

Schaye J. et al., 2015, MNRAS, 446, 521

Schechter P., 1976, ApJ, 203, 297

Semenov D., Henning T., Helling C., Ilgner M., Sedlmayr E., 2003, A&A, 410, 611

Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337

Sharma P., McCourt M., Quataert E., Parrish I. J., 2012, MNRAS, 420, 3174

Sijacki D., Springel V., Di Matteo T., Hernquist L., 2007, MNRAS, 380, 877

Silk J., Rees M. J., 1998, A&A, 331, L1

Skinner M. A., Ostriker E. C., 2013, ApJS, 206, 21

Springel V., Hernquist L., 2003, MNRAS, 339, 289

Sutherland R. S., Bicknell G. V., 2007, ApJS, 173, 37

Tacconi L. J. et al., 2010, Nature, 463, 781

Tamburro D., Rix H.-W., Leroy A. K., Mac Low M.-M., Walter F., Kennicutt R. C., Brinks E., de Blok W. J. G., 2009, AJ, 137, 4424

Teyssier R., 2002, A&A, 385, 337

Tombesi F., Meléndez M., Veilleux S., Reeves J. N., González-Alfonso E., Reynolds C. S., 2015, Nature, 519, 436

Toro E. F., Spruce M., Speares W., 1994, Shock Waves, 4, 25

Tsang B. T.-H., Milosavljević M., 2015, MNRAS, 453, 1108

Verner D. A., Ferland G. J., Korista K. T., Yakovlev D. G., 1996, ApJ, 465, 487

Vogelsberger M. et al., 2014, MNRAS, 444, 1518

Wagner A. Y., Bicknell G. V., 2011, ApJ, 728, 29

Wagner A. Y., Umemura M., Bicknell G. V., 2013, ApJ, 763, L18

Warhaft Z., 2000, Annu. Rev. Fluid Mech., 32, 203

Wise J. H., Abel T., 2011, MNRAS, 414, 3458

Zubovas K., King A., 2012, ApJ, 745, L34

Zubovas K., Nayakshin S., 2014, MNRAS, 440, 2625

Zubovas K., Nayakshin S., King A., Wilkinson M., 2013, MNRAS, 433, 3079

## APPENDIX A: EFFECT OF DIFFERENT SPEED OF LIGHT APPROXIMATIONS

With the explicit M1 scheme that we use for the transport of radiation between grid cells, the RHD timestep is limited by the speed of

<figure><img src="image_17.png" /><figcaption>**Figure A1.** Convergence test of the mechanical advantage for the lower resolution L46_medC_medρQ simulation using different values for the reduced speed of light $c_{\text{red}}$, whose default value is $0.2\ c$ in our simulations. Convergence occurs for $c_{\text{red}} \ge 0.2\ c$.</figcaption></figure>

light, since the solver breaks down if radiation is allowed to travel more than one cell width in one timestep. We thus use the RSLA (Gnedin & Abel 2001), with a fiducial light speed of $c_{\text{red}} = 0.2\ c$ in this work.

To understand how much our reduced speed of light affects our results, we have run lower resolution equivalents to L46_medC_medρQ with a spatial resolution of $\Delta x = 11.6\ pc$ (as opposed to fiducial resolution of $\Delta x = 5.8\ pc$ used in the rest of the paper), and where we adopt different values for the reduced speed of light $c_{\text{red}}$. The evolution of the mechanical advantage, shown in Fig. A1, indicates that the choice of $c_{\text{red}}$ has a significant effect on the mechanical advantage, especially at the beginning of the simulation. In fact, the mechanical advantage can increase by an order of magnitude, at early times (0.1 Myr), if the reduced speed of light fraction is increased from $c_{\text{red}} = 0.005$ to unity. Yet, this increase is very limited when changing from $c_{\text{red}} = 0.2\ c$ (as in our study) to $c$. After ~2 Myr, when the radiation has managed to carve a hole through which it escapes, all runs converge towards the same momentum input rate.

We conclude that our results regarding the mechanical advantage are not significantly altered by the reduced speed of light to $0.2\ c$ and that our results are thus well converged. In a forthcoming paper (Bieri et al., in preparation), we will discuss in more detail the measured reduction factor and its dependence on the reduced speed of light.

This paper has been typeset from a TeX/LA TeX file prepared by the author.