Exponentially growing bubbles around early super massive black holes

TL;DR

This paper models the growth and impact of outflows driven by early supermassive black holes, revealing exponential expansion dynamics and matching observed thermal energies around quasars at high redshift.

Contribution

It provides analytic solutions for outflow evolution driven by exponentially growing black holes and compares predictions with observational data at z=6.

Findings

Outflow radius grows exponentially with time, with specific e-folding times.

Final bubble sizes depend on halo mass, not seed black hole mass.

Predicted thermal energies align with Sunyaev-Zeldovich observations.

Abstract

We addressed the so far unexplored issue of outflows induced by exponentially growing power sources, focusing on early supermassive black holes (BHs). We assumed that these objects grow to $10^9\;M_{\odot}$ by z=6 by Eddington-limited accretion and convert 5% of their bolometric output into a wind. We first considered the case of energy-driven and momentum-driven outflows expanding in a region where the gas and total mass densities are uniform and equal to the average values in the Universe at $z>6$. We derived analytic solutions for the evolution of the outflow, finding that, for an exponentially growing power with e-folding time $t_{Sal}$, the late time expansion of the outflow radius is also exponential, with e-folding time of $5t_{Sal}$ and $4t_{Sal}$ in the energy-driven and momentum-driven limit, respectively. We then considered energy-driven outflows produced by QSOs at the center of early dark matter halos of different masses and powered by BHs growing from different seeds. We followed the evolution of the source power and of the gas and dark matter density profiles in the halos from the beginning of the accretion until $z=6$. The final bubble radius and velocity do not depend on the seed BH mass but are instead smaller for larger halo masses. At z=6, bubble radii in the range 50-180 kpc and velocities in the range 400-1000 km s$^{-1}$ are expected for QSOs hosted by halos in the mass range $3\times10^{11}-10^{13}\;M_{\odot}$. By the time the QSO is observed, we found that the total thermal energy injected within the bubble in the case of an energy-driven outflow is $E_{th}\sim5 \times 10^{60}$ erg. This is in excellent agreement with the value of $E_{th}=(6.2\pm 1.7)\times 10^{60}$ erg measured through the detection of the thermal Sunyaev-Zeldovich effect around a large population of luminous QSOs at lower redshift. [abridged]

Figures

Equations

$$M(t)=M_0{e}^{t/t_{Sal}},$$

$$t_{Sal}=\frac{ϵ}{1-ϵ}\frac{t_E}{\lambda }=50Myr\left(\frac{9ϵ}{1-ϵ}\right){\lambda }^{-1}$$

$$L_{bol}(t)=\lambda L_E=\lambda \frac{M(t)c^2}{t_E}=L_{bol,0}{e}^{t/t_{Sal}},$$

$$L_w(t)=\frac{1}{2}{\dot{m}}_w{v}_w^2=L_{w,0}{e}^{t/t_{Sal}}=f_w{L}_{bol,0}{e}^{t/t_{Sal}},$$

$${\dot{m}}_w(t)={\dot{m}}_{w,0}{e}^{t/t_{Sal}},$$

$${\dot{m}}_{w,0}=\frac{\lambda M_0}{ϵt_E}=2.2\times 10^{-4}{M}_{\odot }y{r}^{-1}\lambda \left(\frac{0.1}{ϵ}\right)\left(\frac{M_0}{10^4{M}_{\odot }}\right).$$

$$m_w(t)={\int }_0^t{\dot{m}}_{w}d{t}^{\prime }={\dot{m}}_{w,0}{t}_{Sal}(e^{t/t_{Sal}}-1).$$

$$E_{th}=\frac{3}{2}P\frac{4\pi }{3}{R}^3=2\pi P{R}^3,$$

$$\frac{d}{dt}[M(R)\dot{R}]=4\pi R^{2}P,$$

$${\dot{E}}_{th}=L_w-4\pi R^{2}P\dot{R}.$$

$$P=L_w/(2\pi v_w{R}^2),$$

$${\rho }_b(z)={\rho }_{b,0}(1+z{)}^3,$$

$$15R^2{\dot{R}}^3+12R^3\dot{R}\ddot{R}+R^4\stackrel{...}{R}=\frac{3}{2\pi }\frac{L_w(t)}{{\rho }_{gas}}.$$

$$R(t)=R_0{e}^{t/(5t_{Sal})},$$

$$R_0={\left(\frac{375}{56\pi }\right)}^{\frac{1}{5}}{t}_{Sal}^{\frac{3}{5}}{\left(\frac{L_{w,0}}{{\rho }_{gas}}\right)}^{\frac{1}{5}}$$

$$=28.5kpc{\left(\frac{9ϵ}{1-ϵ}\right)}^{\frac{3}{5}}{\lambda }^{-\frac{2}{5}}{\left(\frac{f_w}{0.05}\right)}^{\frac{1}{5}}{\left(\frac{M_0}{10^4{M}_{\odot }}\right)}^{\frac{1}{5}}{\left(\frac{{\rho }_{gas}}{10^{-28}gc{m}^{-3}}\right)}^{-\frac{1}{5}}.$$

$$R(t)\simeq {\left(\frac{8}{11}\right)}^{\frac{1}{5}}{R}_0{\left[e^{t/t_{Sal}}-1-\frac{t}{t_{Sal}}-\frac{1}{2}{\left(\frac{t}{t_{sal}}\right)}^2\right]}^{\frac{1}{5}}.$$

$$R(t)\simeq {\left(\frac{4}{33}\right)}^{\frac{1}{5}}{R}_0{\left(\frac{t}{t_{Sal}}\right)}^{\frac{3}{5}}\sim 0.65R_0{\left(\frac{t}{t_{Sal}}\right)}^{\frac{3}{5}}.$$

$$v(t)\simeq \frac{8}{55}\frac{R_0}{t_{Sal}}{\left(\frac{R_0}{R(t)}\right)}^4\left(e^{t/t_{Sal}}-1-\frac{t}{t_{Sal}}\right).$$

$$v(t)\simeq {\left(\frac{4}{33}\right)}^{\frac{1}{5}}\frac{3}{5}\frac{R_0}{t_{Sal}}{\left(\frac{t}{t_{Sal}}\right)}^{-\frac{2}{5}},$$

$$\frac{d}{dt}[\frac{4\pi }{3}{\rho }_{gas}{R}^3\dot{R}]=\frac{2L_{0,w}}{v_w}{e}^{t/t_{Sal}}.$$

$$R(t)=R_0{\left(e^{t/t_s}-1-\frac{t}{t_s}\right)}^{\frac{1}{4}},$$

$$R_0={\left(\frac{6L_{w,0}}{\pi v_w{\rho }_{gas}}\right)}^{\frac{1}{4}}{t}_s^{\frac{1}{2}}$$

$$=10.5kpc{\left(\frac{9ϵ}{1-ϵ}\right)}^{\frac{1}{2}}{\lambda }^{-\frac{1}{2}}{\left(\frac{f_w}{0.05}\right)}^{\frac{1}{4}}{\left(\frac{M_0}{10^4{M}_{\odot }}\right)}^{\frac{1}{4}}{\left(\frac{{\rho }_{gas}}{10^{-28}gc{m}^{-3}}\right)}^{-\frac{1}{4}}.$$

$$R(t)\simeq {\left(\frac{3L_{w,0}}{\pi v_w{\rho }_{gas}}\right)}^{\frac{1}{4}}{t}^{\frac{1}{2}}.$$

$$R(t)\simeq {\left(\frac{6L_{w,0}}{\pi v_w{\rho }_{gas}}\right)}^{\frac{1}{4}}{t}_s^{\frac{1}{2}}{e}^{t/(4t_s)}.$$

$$v(t)=\frac{R_0}{4t_{Sal}}{\left(\frac{R_0}{R(t)}\right)}^3(e^{t/t_{Sal}}-1).$$

$$t_U(z)=\frac{2}{3H_0{\Omega }_{\Lambda }^{1/2}}ln\left(\frac{1+cos\theta }{sin\theta }\right),$$

$$\frac{d}{dt}[M_s(R)\dot{R}]=4\pi R^{2}P-{\mathcal{F}}_g,$$

$${\dot{E}}_{th}=L_w-4\pi R^{2}P\dot{R}-{\mathcal{F}}_g\dot{R},$$

$$15R^2{\dot{R}}^3+12R^3\dot{R}\ddot{R}+R^4\stackrel{...}{R}+12\pi G\frac{{\rho }_{gas}}{f_{gas}}{R}^4\dot{R}=\frac{3}{2\pi }\frac{L_w(t)}{{\rho }_{gas}}.$$

Full text

1 Introduction

Active galactic nuclei (AGN) may drive powerful gas outflows out to distances much greater than kpc. Observations have shown that the outflowing gas may be in a broad range of ionization states and may cover large solid angles. For instance, the detection of blueshifted metal absorption lines in the X-ray spectra of nearby AGN revealed the presence of ultra-fast outflows (UFOs; see e.g. Tombesi et al. 2012, 2013 and references therein) produced within a few milli-pc of the central black hole (BH). Ultra-fast outflows are made of plasma moving at 0.1c on average, and up to $\sim 0.3$c (Tombesi et al., 2015; Nardini et al., 2015). Although it is difficult to constrain the opening angle of UFOs, in the powerful local quasi-stellar object (QSO) PDS456 the remarkable P Cygni profile of the FeXXVI $K\alpha$ line indicates the presence of a subrelativistic wind expanding almost spherically. On much larger scales ($\sim$kpc), outflows of less ionized gas have been traced through observations of asymmetric emission lines in the optical regime, such as the [OIII]5007Å line (Brusa et al., 2015; Harrison et al., 2016; Shen, 2016; Zakamska et al., 2016). Outflows of neutral atomic and molecular gas have also been observed. The neutral atomic component has been probed by observations of the Na I D absorption doublet (Krug et al., 2010) or by the [CII]158$\mu$m and HI 21cm fine structure emission lines (Maiolino et al., 2012; Morganti et al., 2016). The molecular component is mostly probed by observations of CO transitions (Feruglio et al., 2010; Cicone et al., 2014). Velocities higher than 1000 km s*-1* are commonly observed in all of these gas phases. The molecular component generally dominates the total mass budget of the large-scale outflowing gas. The energy required to drive large-scale outflows in AGN is so large that in most cases it can only be provided by nuclear activity rather than by star formation in the host (Cicone et al., 2014). In at least two QSOs, namely Mrk231 (Feruglio et al., 2015) and IRAS F11119 (Tombesi et al., 2015), a UFO and a molecular outflow have both been observed. The comparison between the energy of the inner UFOs and that of the large-scale molecular outflows suggest that AGN outflows are closer to the energy-driven limit than to the momentum-driven limit (see Section 2); that is, the UFOs are able to inflate a bubble of hot gas that does not cool efficiently and whose thermal pressure can push layers of colder gas to large distances. Very recently, further support for the existence of hot gas bubbles around AGN came from the detection of the thermal Sunyaev-Zeldovich effect on the stacked far-infrared spectra of large populations of $z\sim 2$ QSOs (Ruan et al., 2015; Crichton et al., 2016)

AGN-driven outflows are obviously best studied in nearby objects, but they are observed in AGN at all redshifts (Brusa et al., 2015; Zakamska et al., 2016), even beyond redshift 6 (Maiolino et al., 2012), i.e. in the most distant QSOs known. A remarkable example is given by the famous QSO SDSSJ1148+5251 (Fan et al., 2003), one of the first $z>6$ QSOs discovered by the SDSS survey, where observations of the [CII]158$\mu$m line with the PdBI interferometer have revealed the presence of large-scale gas outflows extending out to 30 kpc from the QSO and moving with velocities of $>1000$ km s*-1* (Cicone et al., 2015). This large-scale outflow extends well beyond the characteristic size of the host galaxies of $z\sim 6$ QSOs ($\sim 1-3$ kpc Wang et al. 2013; Venemans et al. 2017), and then expands on scales comparable with that of the hosting dark matter halo, possibly affecting the QSO local environment.

The theory of QSO driven outflows has been presented in many works, either using an analytic/numerical approach (King, 2010; King et al., 2011; Zubovas & King, 2012; Faucher-Giguère & Quataert, 2012), detailed hydrodynamical simulations (Nayakshin & Zubovas, 2012; Costa et al., 2014), or a combined approach linking simulations at different scales through analytic recipes (Gaspari & Sa̧dowski, 2017). A common assumption in most of these works is that the source is constant in power and that the outflow expands within a dark matter halo potential and gas density profile which are constant. In general, despite the vast literature on outflows produced by astrophysical sources (supernovae- , star- , AGN-driven winds), only energy sources that are impulsive, constant, or evolve as power laws have been considered. In this work we explore the issue of outflows produced by exponentially growing power sources, such as should be the case in early QSOs. In fact, to explain the billion solar mass black holes observed in QSOs at $z>6$ (Willott et al., 2010; De Rosa et al., 2014), i.e. when the Universe was less than 1 Gyr old, black holes must have grown exponentially by accreting continuously at the Eddington limit (but see e.g. Madau et al. 2014; Volonteri et al. 2015; Pezzulli et al. 2016 for intermittent, super-critical growth of early black holes). If a fixed fraction of the QSO power is transferred to an inner UFO, then the power of this wind must also grow exponentially. Although assuming continuous and efficient BH accretion for several e-folding times is one of the few ways of coping with the still unsolved issue of the formation of the first supermassive black holes (SMBHs), we note here that this is inconsistent with simulations of BH growth at high Eddington ratios. Pre-heating instabilities that halt BH fueling already arise at Eddington ratios of $\sim 0.01$ in the case of spherical accretion (Cowie et al., 1978). Hydrodynamic simulations in 2D have shown that this also occurs when removing the constraint of spherical accretion (Novak et al., 2011; Ciotti & Ostriker, 2012). At high Eddington ratios an additional instability effect is expected to take place, due to a recurrent thickening of the accretion disc and the consequent increase of its covering factor. Because of these arguments, uninterrupted BH accretion is therefore unlikely, yet it can be regarded as a useful approximation of the long-term behaviour of a more erratic growth. Finally, we note that at high redshifts the mass growth of the hosting dark matter halo and the change in the gas profile within the halo are not negligible during the accretion time of the QSO, so we followed them self-consistently during the outflow expansion, in contrast to most previous works.

A concordance cosmology with $H_{0}=70$ km s*-1* Mpc*-1* , $\Omega_{m}=0.3$ , $\Omega_{\Lambda}=0.7$, in agreement within the errors with the Planck 2015 results (Planck Collaboration et al., 2016), is used throughout this paper.

2 Problem set-up and definitions

We consider the case of an accreting supermassive black hole radiating at a given fraction $\lambda\equiv L_{bol}/L_{E}$ of its Eddington luminosity $L_{E}=4\pi Gm_{p}cM/\sigma_{T}=Mc^{2}/t_{E}=1.26\times 10^{38}M$ erg s*-1*, where $M$ is the black hole mass in units of $M_{\odot}$ and $t_{E}=\sigma_{T}c/(4\pi Gm_{p})=0.45$ Gyr is the Eddington time.

If a fraction $\epsilon$ of the mass $m_{acc}$ falling onto the black hole is converted into radiation, then $L_{bol}=\epsilon\dot{m}_{acc}c^{2}$, with $\dot{m}_{acc}=\lambda L_{E}/(\epsilon c^{2})=\lambda M/(\epsilon t_{E})$. The black hole growth rate is then $\dot{M}=(1-\epsilon)\dot{m}_{acc}=(1-\epsilon)\lambda M/(\epsilon t_{E})$ , from which it follows that

[TABLE]

where

[TABLE]

is the Salpeter (e-folding) time, and $M_{0}$ is the mass of the black hole “seed”. Under these assumptions, the bolometric QSO luminosity grows as

[TABLE]

where $L_{bol,0}=\lambda c^{2}M_{0}/t_{E}$. It is assumed that a constant fraction $f_{w}$ of the AGN bolometric luminosity powers a supersonic wind pushing on the surrounding ambient medium. This fraction is typically considered to be $f_{w}$=0.05 (Scannapieco & Oh, 2004; Lapi et al., 2005; Germain et al., 2009). The outflow is also assumed to be almost spherical (i.e. a bubble) with covering factor $\Omega\approx 4\pi$. We then investigate the evolution of astrophysical bubbles produced by an exponentially growing input energy source,

[TABLE]

where $\dot{m}_{w}$ and $v_{w}$ are the wind mass rate and velocity, respectively.

Observations of ultra-fast outflows in the X-ray spectra of QSOs (e.g. Tombesi et al. 2013) suggest that the gas outflow rate is of the order of the gas accretion rate onto the BH, so we assume $\dot{m}_{w}\approx\dot{m}_{acc}$. We further assume that the QSO wind is launched with constant velocity.111Under the above assumptions, Eq. 4 implies that $2f_{w}\epsilon(c/v_{w})^{2}=1$, i.e. $v_{w}=0.1c$ for $f_{w}=0.05$ and $\epsilon=0.1$ as assumed here, in agreement with the velocities measured for UFOs. Therefore, the outflow rate grows as

[TABLE]

where

[TABLE]

In turn, the ejected wind mass grows as

[TABLE]

3 Energy-driven vs. momentum-driven outflows and outflow structure

The classic theory of astrophysical bubbles driven by fast winds (see Weaver et al. 1977 and Ostriker & McKee 1988 for reviews) shows that the structure of the outflow is characterized by a forward and a reverse shock, and is divided into four main regions (see Fig. 1). The forward shock propagates in the ambient medium (region $d$ in Fig. 1), enhancing its density in a shell (region $c$) that progressively sweeps up more gas. The outflowing wind (region $a$) instead encounters the innermost reverse shock, and is therefore heated. The region $b$ of hot shocked wind is usually termed a “bubble”. The cooling properties of the shocked wind define whether an outflow is energy- or momentum-driven. In the limit in which this gas does not cool, all the energy originally provided by the wind is retained within the system, and we refer to the outflow as “energy-driven” 222Here we neglect the cooling of the gas within the shell. This would, in any case, affect the estimates of the bubble radius and velocity by less than 15% (Weaver et al., 1977).. It is the thermal pressure of the hot gas within the bubble that accelerates the shell of gas into the ambient medium. On the contrary, if all thermal energy within the bubble is radiated away instantaneously, energy is not conserved, whereas all the momentum carried by the wind is transferred directly to the shell of shocked ambient gas. The outflow is therefore referred to as “momentum-driven”, which effectively corresponds to maximal cooling. Realistic outflows obviously fall between these two limits. It is, however, instructive to treat them separately.

We follow an approach similar to that of Weaver et al. (1977), who provided the basics for treating the structure and evolution of bubbles inflated by fast stellar winds within a medium of constant density $\rho_{gas}$. Their computations can be generalized to astrophysical bubbles powered by other energy sources. In the energy-driven case, the following equations can be combined to obtain the full equation of motion of the expanding shell:

[TABLE]

[TABLE]

[TABLE]

Equation 8 shows the relation between thermal energy $E_{th}$ and pressure $P$ for a monoatomic gas in a volume $V=(4\pi/3)R^{3}$.

Eq. 9 describes the evolution of the momentum of the shell of swept-up gas (the left-hand side is the rate of change of momentum of the shell and the right-hand side represents the force acting on it).

Eq. 10 is the energy equation of the system: the rate at which the thermal energy of the hot gas in the bubble (left-hand side) changes is equal to the kinetic luminosity of the wind minus the rate at which the hot gas does work on the surrounding medium. In the momentum-driven case, $E_{th}=0$ and the pressure $P$ in Eq. 9 is given by the wind ram pressure

[TABLE]

where $L_{w}$ and $v_{w}$ are the wind kinetic power and velocity, respectively.

Weaver et al. (1977) considered the case of a constant source of energy. Energy inputs evolving in time as power laws were considered by e.g. Koo & McKee (1992). Here we consider sources with exponentially growing luminosities.

In our equations we have assumed that the effects of the pressure of the ambient gas and of gravity are both negligible. We explore the effects of relaxing these assumptions in the next sections.

4 QSO in a uniform density field

We first assume that the BH forms and grows, and hence the outflow propagates, within a region with uniform density equal to the mean matter density of the Universe. We recall that the cosmic evolution of the baryon density is given by

[TABLE]

where $\rho_{b,0}=3H_{0}^{2}\Omega_{b,0}/(8\pi G)$ is the mean baryon density at $z=0$. For $H_{0}=70$ km s*-1* Mpc*-1* and $\Omega_{b,0}=0.045$, $\rho_{b,0}=4\times 10^{-31}$ g cm*-3*. At $z=6$ the mean baryon density is therefore $\rho_{gas}\equiv\rho_{b}(6)\sim 10^{-28}$ g cm*-3*. For an ambient gas with uniform density $\rho_{gas}$, the mass of the shell of swept-up gas $M_{s}(R)=\int_{0}^{R}\rho_{gas}(R)dV$ is equal to $(4\pi/3)R^{3}\rho_{gas}$.

4.1 Energy-driven case

Following Dyson & Williams (1997), we combine Eqs. 8, 9, and 10 to obtain the following equation of motion of the shell:

[TABLE]

This is a non-linear differential equation that admits more than one solution, and where not all of the solutions can be obtained by linear combination of the others. For an exponential input energy source described by Eq. 4, the simplest (but unphysical, see below) solution to Eq.13 is analytic and takes the form

[TABLE]

with

[TABLE]

In Fig. 2 we show the time evolution of the bubble radius $R(t)$ up to the accretion time $t_{9}$, defined as the time at which the black hole has grown to $10^{9}\;M_{\odot}$, i.e. $M(t_{9})\equiv 10^{9}\;M_{\odot}$. This is the typical mass observed for bright SDSS QSOs at $z\sim 6$. For this computation we assumed

$f_{w}=0.05$ ,

$\epsilon=0.1$ ,

$\lambda=1$ ,

$M_{0}=10^{4}\;M_{\odot}$ .

The adopted seed mass $M_{0}$ is an intermediate value within the broad range predicted by different models. This range spans the whole interval from $\sim 10^{2}\;M_{\odot}$, as expected for the remnants of POPIII stars (Madau & Rees, 2001), to $\sim 10^{4}\;M_{\odot}$, as expected from the collapse of stellar clusters (Devecchi & Volonteri, 2009), and even up to $10^{6}\;M_{\odot}$, as postulated for large seeds forming from the direct collapse of large gas clouds under specific environmental conditions (Agarwal et al. 2014; see also Volonteri & Bellovary 2012 for a review and Valiante et al. 2016 for a model of seed formation at different mass scales). With the above assumptions the initial power of the wind is $L_{w,0}=6.3\times 10^{40}$ erg s*-1*, and $t_{9}\sim 580$ Myr.

The simple solution given in Eq. 14 would imply that at $t=0$ the radius of the bubble is $R=28.5$ kpc, i.e. the bubble “jump starts” on scales that are far larger than the sphere of influence of the seed BH and even larger than its putative host galaxy. Since $R_{0}$ scales weakly with the seed mass ($R_{0}\propto M_{0}^{0.2}$), assuming a stellar-sized seed with $M_{0}=10^{2}\;M_{\odot}$ would only reduce $R_{0}$ by a factor of 2.5. The pure exponential solution therefore appears unphysical.

We therefore solved Eq. 13 numerically by imposing $R\rightarrow 0$ (and $v\equiv\dot{R}\rightarrow 0$) for $t\rightarrow 0$ and obtained a simple analytic approximation to the numerical solution by manipulating Eq. 14 and subtracting its second-order Taylor expansion from the exponential term

[TABLE]

As shown in Fig. 2, Eq. 16 represents an excellent approximation (within 5%) of the bubble radius over the entire accretion history of the BH. At late times ($t>>t_{Sal}$) this solution approaches the simple analytic solution given by Eq. 14. At early times ($t<<t_{Sal}$), the numerical solution approaches the analytic solution that is valid for a constant source of power $L_{w}(t)=L_{w,0}$ (as assumed by Weaver et al. 1977 or Dyson & Williams 1997):

[TABLE]

As a further check we ran a test low-resolution hydrodynamical simulation in which a source embedded in a uniform density field injects energy and matter at an exponential rate. The simulation was performed using a customized version of the RAMSES code (Teyssier, 2002) as described by Calura et al. (2015). The initial conditions are represented by a uniform and homogeneous distribution of gas characterized by a density of $10^{-28}$ g cm*-3* and temperature $10^{4}$ K (see Section 4.4).

The computational box has a volume of $L_{box}^{3}=(650$ kpc$)^{3}$ and is characterized by a maximum resolution of $\sim 130$ pc. At $t=0$, a source located in the origin of the computational box starts injecting energy at a rate given by Eq. 4, with $L_{w,0}=6.3\leavevmode\nobreak\ \times\leavevmode\nobreak\ 10^{40}$ erg s*-1* and matter at a rate given by Eq. 5, with $\dot{m}_{w,0}\sim 2.2\leavevmode\nobreak\ \times\leavevmode\nobreak\ 10^{-4}\;M_{\odot}$ yr*-1*. In order to avoid the occurrence of diamond-shaped shock fronts which can sometimes be present in similar conditions in cartesian grid-based simulations (Tasker et al., 2008), our source is non-point-like but distributed over a spherical volume $V=\frac{4}{3}\pi\Delta R^{3}$, with $\Delta R=200$ pc. The total thermal energy and mass injected by the source in the surrounding environment per unit volume and in the time step $\Delta t$ are $\Delta\epsilon=\dot{\epsilon}\Delta\,t$ and $\Delta\rho=\dot{\rho}\Delta\,t$, respectively, where $\dot{\rho}=\frac{\dot{m}_{w}}{V}$ and $\dot{\epsilon}=\frac{L_{w}}{V}$.

At each time step, a number of cells with the highest refinement level are created at the position of the source. In the remainder of the computational domain, the refinement strategy is both geometry- and discontinuity-based. Outflow boundary conditions are used. The simulation is adiabatic, i.e. the effects of radiative cooling are assumed to be negligible. The effects of gravity are also neglected. As shown in Fig. 2, the expansion of the shell from the simulation is in excellent agreement with the numerical solution to Eq. 13.

At late times, Eq. 16 is well approximated by Eq. 14 and the growth of the bubble radius is purely exponential. Some interesting properties can be inferred by looking at Eq. 14. For instance, the e-folding time of the bubble radius is five times the Salpeter time. Also, it can be shown that the bubble radius (and velocity) at a given black hole mass (e.g. $10^{9}\;M_{\odot}$) does not depend on the initial seed mass: for smaller seeds, it will simply take longer to accrete that mass and blow the bubble to that radius. This is visible in Fig. 3, where the effects of varying the seed BH mass from $M_{0}=10^{2}M_{\odot}$ to $10^{6}M_{\odot}$ are shown. The bubble radius when the BH has grown to $10^{9}M_{\odot}$ is $R(t_{9})=267$ kpc whatever the initial seed mass is. This value is three times smaller than the radius reached in the same time span by a bubble inflated by a constant wind power of $L_{w}=f_{w}L_{bol}=f_{w}L_{E}=6.3\times 10^{45}$ erg s*-1*, corresponding to a $10^{9}M_{\odot}$ BH radiating at its Eddington limit and with $f_{w}=0.05$. This is shown by the magenta curve in Fig. 3: as expected, the normalization of this curve is a factor of $(10^{9}/10^{4})^{1/5}$ = 10 higher than the black dotted line in Fig. 2, which was computed for a $10^{4}M_{\odot}$ BH (again radiating at its Eddington limit and with $f_{w}=0.05$). As a matter of fact, previous works on QSO outflows do not consider the increase in the QSO luminosity and wind power as the BH grows, but rather assume a constant release of energy.

By solving numerically Eq. 13 we also derived a solution for the time evolution of the bubble velocity $v\equiv\dot{R}$, which is compared in Fig. 4 with the time evolution of the bubble radius and BH growth. The numerical solution for the bubble velocity can be approximated (within $\sim$6%) by the derivative of the analytic approximation to the bubble radius given in Eq. 16:

[TABLE]

From Eq. 18 it can be seen that at late times ($t>>t_{Sal}$) the bubble velocity grows exponentially according to $v(t)\sim R_{0}/(5t_{Sal})\,e^{t/(5t_{Sal})}$, i.e. it approaches the derivative of Eq. 14. Instead, at early times ($t<<t_{Sal}$), the bubble velocity follows the relation

[TABLE]

i.e. we recover the $v(t)\propto t^{-2/5}$ dependence found by Weaver et al. (1977) and Dyson & Williams (1997) for a constant source of power (Eq. 19 can be obtained by taking the derivative of Eq. 17).

As shown in Fig. 4, the most notable feature in the velocity curve is that it reaches a minimum of $\sim$220 km s*-1* at a few Salpeter times, and then starts increasing exponentially. This means that the shell undergoes an acceleration and it becomes Rayleigh-Taylor unstable around a few Salpeter times, i.e. at 200-300 Myr after the BH starts accreting. The shell stability is discussed in Section 5.

4.2 Momentum driven case (maximal cooling)

In momentum driven outflows, the shocked wind cools immediately and the bubble (region $b$ in Fig. 1) collapses, so the wind ram pressure is directly pushing on the gas shell (see e.g. Zubovas & King 2012, Costa et al. 2014). By combining Eqs. 9 and  11 and recalling that the mass of the shell $M(R)=(4\pi/3)R^{3}\rho_{gas}$ for a uniform gas density $\rho_{gas}$, we obtain

[TABLE]

By integrating Eq. 20 twice in time and solving for $R$, we obtain

[TABLE]

where

[TABLE]

This solution is physically meaningful, as $R\rightarrow 0$ for $t\rightarrow 0$.

At early times, $t<<t_{s}$, Eq. 21 reduces to

[TABLE]

The above asymptotic solution correctly recovers the $R\propto t^{\frac{1}{2}}$ dependence found by Steigman et al. (1975) for momentum driven outflows in the case of constant power sources (at early times the exponential term can indeed be approximated with a constant). At late times, $t>>t_{s}$, the exponential term takes over and the equation of motion collapses to

[TABLE]

Eq. 24 shows that the e-folding time of the radius of the expanding shell is 4 times the Salpeter time, i.e. the shell expansion is slower than that of the SMBH and of the power it releases. The expansion of the shell is shown in Fig. 5, where its asymptotic behaviour is also highlighted. By taking the time derivative of Eq. 21 we also derived an expression for the shell velocity in the momentum driven regime as follows:

[TABLE]

In Fig. 6 ($top$ and $middle$ panels) we compare the time evolution of the shell radius and velocity in the energy-driven vs. momentum-driven scenarios: the shell radius in the energy driven case is always larger – by about a factor of two – than in the momentum driven case, apart from very late times when the faster exponential growth of the momentum driven case takes over. Similarly, the shell velocity in the energy driven regime is always higher than in the momentum driven regime except at very late times. It is worth stressing that for $t\gtrsim 1$ Gyr the mass ejected by the BH (growing as $e^{t/t_{Sal}}$) will be larger than the mass of the shell (which, for instance, grows as $e^{0.75t/t_{Sal}}$ in the momentum driven regime). In such conditions our system of equations, which is designed to describe the motion of a dense shell enclosing an empty bubble, is no longer valid. In the $bottom$ panel of Fig. 6 we also show the behaviour of the shell acceleration in the two regimes: in the energy driven limit, the acceleration crosses zero at later times. We return to this in the stability analysis in Section 5.

4.3 Effects of decreasing gas density because of cosmological expansion

We discuss here the effects of the Hubble expansion on the evolution of the bubble radius. Since we placed the accreting BH within an unbound region, the net effect of cosmic expansion is a reduction in the (proper) average background gas density. The relation between cosmic time and redshift for a flat cosmological model with $\Omega_{\Lambda}\neq 0$ can be written as (e.g. Longair 2008)

[TABLE]

where $tan(\theta)=(\Omega_{m}/\Omega_{\Lambda})^{1/2}(1+z)^{3/2}$.

Eq. 26 shows that at z=6 the age of the Universe is $t_{U}=914$ Myr and that for a total accretion time $t_{9}$ of 575 Myr, the epoch at which accretion is supposed to start ($t_{U}=914-575=339$ Myr) corresponds to $z\sim 12.5$.

From Eq. 12 it is then easy to see that the average baryon density of the Universe decreases by a factor of $\sim 7$ from $z=12.5$ to $z=6$. In the energy-driven case, the normalization of the shell radius $R_{0}$ scales as $\rho_{gas}^{-1/5}$, so at early times the bubble radius can be expected to be a factor of $\sim 7^{0.2}$ ($\sim 1.5$) smaller than that shown in Fig. 2. In the momentum-driven limit, where $R_{0}\propto\rho_{gas}^{-1/4}$, the bubble radius would decrease by a factor of $\sim 7^{0.25}$ ($\sim 1.6$) at most at early times.

4.4 Effects of the external pressure

In Section 3 we derived the equation of motion of the shell neglecting the effects of the pressure of the ambient gas which in principle may slow down the expansion of the bubble significantly and alter the physical properties of the shell. In particular, if the expansion velocity of the bubble is comparable to the sound speed in the ambient gas, no shock will be generated and consequently no shell of dense gas will form. In the redshift range considered here (z=6-12.5), i.e. before a full re-ionization of the IGM, we can assume that the temperature of the IGM around a QSO is $T\approx 2\times 10^{4}$ K at most if either the QSO radiation or that of stars in its host or in nearby galaxies is heating the IGM through photoionization (Ferrara et al., 2000; Mo et al., 2010; Kakiichi et al., 2016). The sound speed in the IGM, $c_{s}=(\gamma k_{B}T/\mu m_{p})^{1/2}$ will then be $\sim 22$ km s*-1* at most (assuming $\gamma=5/3$ and a mean molecular weight per particle of $\mu=0.59$, valid for primordial gas). This is much smaller than the minimum shell velocity computed in the previous section for the energy-driven and the momentum-driven driven scenarios. The bubble expansion is therefore always supersonic and, as opposed to the case of a constant power source where the shell velocity decreases and the ambient pressure exerts a significant resistance at late times (Weaver et al., 1977), in the exponential case the bubble velocity grows at late times, making the effect of the ambient pressure negligible. As discussed in Section 6, the effects of the external pressure are instead significant even for an exponentially increasing source power if the QSO is placed at the centre of a large dark matter halo.

4.5 Effects of gravity

In Sections 3 and 4 we derived the equation of motion of the shell neglecting the gravitational pull exerted on it by the total mass $M_{t}(<R)$ enclosed within the expanding shell radius. When including gravity, the energy and momentum equations, Eq. 9 and Eq. 10, turn into

[TABLE]

and

[TABLE]

respectively, where $\mathcal{F}_{g}=\frac{GM_{s}(R)M_{t}(<R)}{R^{2}}$ is the gravitational force exerted by the total mass $M_{t}(<R)$ enclosed within the radius of the shell of swept-up gas on the shell mass itself $M_{s}(R)$, and $\mathcal{F}_{g}\dot{R}$ is the rate of work done to lift the shell out of the gravitational potential. For a uniform density field $M_{s}(R)=(4\pi/3)R^{3}\rho_{gas}$, and, for sufficiently large radii where the black hole attraction can be neglected, $M_{t}(<R)=(4\pi/3)R^{3}(\rho_{gas}/f_{gas})$, where $f_{gas}=\rho_{gas}/\rho_{m}$ is the gas mass fraction of the density field. By combining Eq. 8 with Eq. 27 and Eq. 28, after some algebra we get the following equation of motion of the shell:

[TABLE]

The above simple generalization of Eq. 13 accounts for the effects of gravity on a gas shell expanding in a constant density field. It is indeed identical to Eq. 13, but for the addition of the gravitational term $\propto R^{4}\dot{R}$ on the left-hand side.

We explored under which conditions the pull exerted by gravity significantly slows down the shell expansion compared to what was shown in Section 4. To make a comparison with the original solutions of Weaver et al. (1977) and Dyson & Williams (1997), we first considered a constant energy source $L_{w}(t)=L_{w,0}$, which, in the absence of gravity, produces a bubble expanding as $R(t)\propto t^{\frac{3}{5}}$. We then tried a power-law solution to Eq. 29 of the form $R(t)=R_{0}t^{\alpha}$, which leads to the relationship

[TABLE]

For $\alpha=\frac{3}{5}$ the above equation becomes

[TABLE]

which, for $t\rightarrow 0$, returns $R_{0}=\left(\frac{125}{154\pi}\frac{L_{w,0}}{\rho_{gas}}\right)^{\frac{1}{5}}\equiv R_{early,0}$, i.e. at early times the solution of Weaver et al. (1977) applies:

[TABLE]

Instead, for $\alpha=\frac{1}{5}$ the equation can be written as

[TABLE]

which, for $t\rightarrow\infty$, returns $R_{0}=\left(\frac{5}{8\pi^{2}G}\frac{L_{w,0}f_{gas}}{\rho_{gas}^{2}}\right)^{\frac{1}{5}}\equiv R_{late,0}$, and the shell equation of motion at late times becomes

[TABLE]

Therefore, at late times, when the enclosed mass becomes sufficiently large, gravity is effective, and the bubble radius $R(t)$ expands at a lower rate ($\propto t^{\frac{1}{5}}$) than in the standard no-gravity case of Dyson & Williams (1997) ($\propto t^{\frac{3}{5}}$).

A characteristic transition time $t_{c}$ between the two regimes can be defined as the time where the two asymptotic limits on the bubble radius described by Eq. 32 and Eq. 34 cross each other, i.e. $R_{early}(t_{c})=R_{late}(t_{c})$. This returns

[TABLE]

From Eq. 35 it is easy to see that gravity effects are negligible for a bubble inflated by a QSO expanding in a region with density equal to the mean matter density of the Universe at z=6 ($\rho_{gas}=10^{-28}$g cm*-3*). Instead, if the expansion happens within highly overdense regions, or at very early times when the Universe was smaller and denser, gravity effects may become relevant. In Fig. 7 we show the dragging effect exerted by a density field with $\rho_{gas}=10^{-25}$g cm*-3* on the bubble radius produced by a constant source with $L_{w}(t)=L_{w,0}=6.3\times 10^{40}$ erg s*-1* (the full solution to Eq. 29 was obtained numerically). In this case, $t_{c}\sim 77$ Myr, and the asymptotic behaviour of the bubble radius both at early and late times can be appreciated.

We then considered the case of an exponentially growing source of energy $L_{w}(t)=L_{w,0}\;e^{t/t_{Sal}}$. Similarly to Eq. 13, even Eq. 29 admits a simple analytic (exponential) solution of the form $R(t)=R^{G}_{0}\;e^{t/(5t_{Sal})}$, which differs from that obtained in the case of no gravity (Eq. 14) only by its normalization. The relation between the two normalizations $R^{G}_{0}$ and $R_{0}$ (Eq. 15) is

[TABLE]

where $\rho^{G}_{gas}\equiv\frac{7}{75\pi G}t_{Sal}^{-2}f_{gas}$. For $\rho_{gas}<<\rho^{G}_{gas}$, $R^{G}_{0}\sim R_{0}$, and the effects of gravity are negligible. Assuming a universal gas fraction of $f_{gas}=0.16$ and the same Salpeter time used in the previous section ($t_{Sal}=$50 Myr), we get $\rho^{G}_{gas}=3\times 10^{-26}$ g cm*-3*. Therefore, in the case of a QSO placed within a medium with density similar to the average at z=6, $\rho_{gas}=10^{-28}$ g cm -3, the effects of gravity can safely be neglected. We note that, because of the weak dependence of the bubble radius on the density ($R\approx\rho_{g}^{-0.2}$), even by assuming a three dex higher background density, one would get $R^{G}_{0}\sim 0.75R_{0}$; in other words, for an exponentially growing source of power the gravity would reduce the bubble radius by only 25% at late times (see the numerical solutions in Fig. 7).

5 Shell stability

The basic structure of wind-driven bubbles was given in Section 2. We now analyse the stability of the shells. For our terminology we refer to to the ambient medium that has passed through the outermost bubble shock as the “shell” (region $c$ in Fig 1). The wind material that has passed through the reverse shock, where the driving wind from the central black hole and its accretion disc impinges on the bubble, is called the shocked wind medium or simply the “bubble” (region $b$ in Fig 1). There are two main instabilities:

The Rayleigh-Taylor instability, which has a classic growth rate given by $\sqrt{2\pi a/\lambda}$ , where $a$ is the acceleration and $\lambda$ is the spatial scale, and where the classic instability pattern of bubble and spikes occurs as a result of the instability. This is highly disruptive for shells in the linear analysis. 2. 2.

The Vishniac instability (Vishniac, 1983), which occurs for thin shells bounded by thermal pressure on one side and ram pressure on the other side. If the shell is perturbed, the ram pressure has an oblique component on it, whereas the thermal pressure is always normal to the shell surface. As a result, material is transported along the shell causing it to oscillate and producing an “overstability”.

In the initial phases, i.e. at $t<100$ Myr, when the bubble growth is similar to that of a constant energy source, the decelerated shell may indeed be prone to the Vishniac overstability (Pittard, 2013; Krause & Diehl, 2014). However, numerical results in the non-linear regime have shown that these oscillations saturate and do not lead to any real disruption (Mac Low & Norman, 1993; Michaut et al., 2012). We therefore only discuss the Rayleigh-Taylor instability in the following, noting clearly that the final non-linear state will have to be calculated by high-resolution simulations.

Our case is different from the classic behaviour of wind-blown bubbles as the gas density within the bubble $\rho_{bubble}$ is expected to increase with time once the exponential term dominates, whereas the density of bubbles inflated by constant winds is expected to decrease with time. For instance, for the stellar winds studied by Weaver et al. (1977), $\dot{m}=const.$, $m\propto t$ and $R\propto t^{3/5}$, from which it follows that $\rho_{bubble}\propto m/R^{3}\propto t^{-4/5}$.

To compute $\rho_{bubble}$ we assume that the wind shock boundary $R_{w}$ (i.e. the boundary between region a and b in Fig.1) is proportional to the shell radius $R$ so that the shells move together as a piston with $R_{w}=\eta R$ with $\eta$ a slowly varying constant of the order of unity. We then obtain a simple formula for the bubble density, assuming for simplicity that this is uniform across the bubble:

[TABLE]

where the numerator (ejected wind mass) has been taken from Eq. 7. The time evolution of $R$ is given by Eq. 16.

Assuming $\eta<<1$, and considering Eqs. 14 and  17, Eq. 37 reduces to $\rho_{bubble}^{late}\simeq\rho_{o}e^{\frac{2}{5}\frac{t}{t_{Sal}}}$ at late times ($t>>t_{Sal}$), and to $\rho_{bubble}^{early}\simeq\rho_{o}(33/4)^{3/5}(t/t_{Sal})^{-4/5}$ at early times ($t<<t_{Sal}$). Here $\rho_{o}=\frac{3}{4\pi}\frac{\dot{m}_{w,0}t_{Sal}}{R_{0}^{3}}$.

The full behaviour of the bubble density and its asymptotic limits are shown in Fig. 8. For the parameters assumed here, the density in the bubble is always much lower than that in the shell $\rho_{shell}=(1+\mathcal{M}^{2})\rho_{gas}$, where $\mathcal{M}$ is the Mach number. For instance, the minimum shell velocity of $\sim 220$ km s*-1* in the energy-driven case considered here (see central panel of Fig. 6) would correspond to a minimum Mach number of $\sim 10$ and hence to a minimum shell density of $\rho_{shell}\sim 10^{-26}$ g cm*-3*. At early times, when the shell is decelerating, the bubble density follows the $\propto t^{-4/5}$ behaviour expected for winds outflowing at a constant mass rate. In this regime, the shell is Rayleigh-Taylor stable. At later times ($t>100$ Myr), when the exponential term dominates, the shell accelerates and is expected to be Rayleigh-Taylor unstable. As shown in Fig. 6 (bottom panel), in the energy-driven limit, the shell acceleration crosses zero at $t\sim 2t_{Sal}$, i.e. this is the time where Rayleigh-Taylor instabilities start developing. Interestingly, for momentum-driven outflows, this instability time is shorter as the acceleration crosses zero at about 1.4$t_{Sal}$. We note that in the fully exponential regime the ejected mass grows more rapidly than the bubble volume, producing an increasing gas density within the bubble. However, only for accretion times longer than 1 Gyr (i.e. BH masses $>>10^{10}\;M_{\odot}$) or large values of $\eta$ ($\eta>0.9$, i.e. the QSO wind freely streaming up to hundred kpc distances) $\rho_{bubble}$ might be larger than $\rho_{shell}$; however, these two situations appear unphysical. Furthermore, as specified in the previous section, for $t>1$ Gyr, the ejected mass becomes comparable with the shell mass and our set of equations is no longer valid.

6 QSO in a growing dark matter halo

We explore here the effects of placing the accreting BH within a spherically symmetric massive dark matter halo. In particular, we follow the evolution of the bubble as the dark matter halo grows between z=12.5 and z=6 and so does its gas content. In fact, the mass of early dark matter halos can grow by more than two orders of magnitude in this redshift range (see e.g. Correa et al. 2015). Furthermore, the gas content of the halo is expected to change its radial profile (e.g. Ferrara & Loeb 2013). All these effects, which are generally neglected in simple analytic models, regulate the expansion rate of the outflows produced by early QSOs.

The typical mass of the halos in which $z\sim 6$ QSOs reside is largely unknown. Based on abundance-matching arguments, early SMBHs should live in $10^{13}\;M_{\odot}$ halos if their duty cycle is of the order of unity. On the other hand, for smaller duty cycles, or if a large, hidden population of early obscured QSOs exists (which is likely, based on what is known at lower redshifts), early SMBHs should be hosted by smaller, more abundant halos. Furthermore, some simulations suggest that QSO feedback may preferentially inhibit gas accretion in large halos, and that early SMBHs at $z\sim 6$ are hosted by halos with $M_{h}\sim 3\times 10^{11}M_{\odot}$ on average (Fanidakis et al., 2013).

To trace the evolution of massive DMHs at early times up to the largest mass scales, we considered the results of the Millennium XXL simulation (MXXL, Angulo et al. 2012). The MXXL is indeed the only N-body simulation encompassing a sufficiently large volume ($\sim 70$ Gpc3) to be able to follow the growth of the rarest and most massive halos at any time.

Angulo et al. (2012) report an analytic fit to the average evolution of the most massive halos from z=6 to z=0. We extrapolated that fit back in time by considering the average halo accretion rates as quoted in Angulo et al. (2012) and derived the following expression for the time evolution of massive dark matter halos from z=12.5 to z=0: $M_{h}(z)=M_{h}(0)e^{-0.783z}$ , where $M_{h}(0)$ is the halo mass at redshift zero. It is then easy to see that, on average, a halo of $\sim 10^{13}\;M_{\odot}$ at z=6 had a mass of $\sim 6\times 10^{10}\;M_{\odot}$ at z=12.5 and should evolve into a $10^{15}\;M_{\odot}$ halo at z=0 (i.e. into a Coma-like massive galaxy cluster). We define the halo mass and virial radius in the standard way, i.e. $M_{h}=(4/3)\pi R_{vir}^{3}200\rho_{crit}$, where $\rho_{crit}=3H^{2}(z)/(8\pi G)$ is the critical density of the Universe. The evolution of the halo mass and virial radius compared to the evolution of the black hole mass are shown in Fig. 9.

We assume that, at any redshift, the dark matter density in the halo follows the Navarro-Frenk-White profile (Navarro et al., 1997) and that as the halo grows its concentration follows the concentration-mass relation derived by Correa et al. (2015) (see their Eq.20). For the massive halos in the redshift range z=6-12.5 considered here, the concentration parameter is nearly constant, $c\sim 2.8$.

For the gas within the halo, we consider the radial profile reported by Ferrara & Loeb (2013), who extended the work by Makino et al. (1998) on the profile of gas in hydrostatic equilibrium within halos with NFW dark matter density profiles, and allowed for gas temperatures lower than the virial temperature $T_{vir}=(\mu m_{p}/2k_{B})\times(GM_{h}/R_{vir})$. Indeed, it has been suggested that gas accreting onto dark matter halos may not be shock heated to $T_{vir}$ if the halo is smaller than a critical mass $M_{shock}$ ($\sim 10^{12}M_{\odot}$ at all redshifts; Dekel & Birnboim 2006). At high redshifts most halos have masses below $M_{shock}$, hence gas accretion onto them might occur primarily through streams of gas that remain cold, and a large fraction of cold gas may be even present in larger halos (see e.g. Fig. 11 in Overzier 2016).

We investigated the energy driven case and modified Eqs. 27 and 28 to account for the time evolution of the mass and profile of the dark matter and gas. We first considered a QSO within a halo whose mass reaches $10^{13}\;M_{\odot}$ at z=6, and assumed that the ambient gas is always at $T_{gas}=T_{vir}$ as the halo grows. In Fig. 10 we show the halo DM and gas profiles at z=12.5 and at z=6. The QSO outflow now has to sweep gas densities that are up to five orders of magnitude higher than in the uniform density field case explored in the previous sections.

We assumed for the accreting black hole the same input parameters as in the previous Sections (i.e. $M_{0}=10^{4}\;M_{\odot}$, $L_{w0}=6.3\times 10^{40}$ erg s*-1*). The expansion of the bubble radius was obtained numerically and is shown by the dashed line in Fig. 11. Based on the results obtained in the previous sections we are now in a position to understand all the phases of the bubble expansion. At early times / small radii (e.g. $<$1 Myr, or $<0.2$ kpc) the exponential source power can be approximated as a constant and the gradient in the gas profile is mild (see Fig. 10); therefore, the bubble radius expands following the classic $t^{3/5}$ law valid for constant sources within a uniform density field. At $t\sim 10$ Myr, the source power can still be approximated as a constant and the mass enclosed within the bubble radius has become sufficiently large to slow down the expansion of the shell of swept-up gas, which now follows the $t^{1/5}$ law found in Section 4.5. Finally, for $t\gtrsim t_{Sal}$, the exponential source power takes over gravity, and the expansion of the bubble is also exponential. At $t_{9}$ the bubble radius is 60 kpc. That is, the outflow has reached scales comparable with those observed in the $z=6.4$ QSO SDSS J1148 (Cicone et al., 2015). As expected, the bubble size is smaller (by a factor of $\sim 4-5$) than what was found in the case of a uniform density field since now the outflow is crossing a gas density that is on average a few dex higher. As for the velocity of the bubble expansion, this reaches $v=450$ km s*-1* at $t=t_{9}$, whereas velocities up to 1000 km s*-1* have been observed for SDSS J1148.

The computations above neglect the effects of the external pressure on the bubble evolution, which, unlike the case of a QSO placed in the field, are now expected to be relevant. We then modified the momentum equation (Eq. 27) by adding on the left-hand side a term equal to $-4\pi R^{2}P_{ext}$, where $P_{ext}=\rho_{gas}k_{B}T_{gas}/\mu m_{p}$. As shown by the solid curve in Fig.11, the external pressure reduces the bubble radius by a factor of $\sim 1.6$ up to a few Salpeter times ($\sim 1.2$ at late times, in the fully exponential regime of the QSO power).

For a given QSO power, the evolution of the bubble expansion strongly depends on the halo mass in which the BH resides. In the case of a very massive halo ($M(t_{9})=10^{13}\;M_{\odot}$), the bubble reaches radii of $\sim 50$ kpc by z=6, with velocities of $\sim 400$ km s*-1*. This is shown in Fig. 12 (left), where the comparison between the velocity of the bubble expansion and the sound speed of the gas in the halo is also shown. For massive halos, the bubble expansion is subsonic during the entire accretion history of the black hole. This means that the QSO wind is not able to create a shock and then no shell of dense gas propagates outwards. In this case the treatment provided by our equations can still be used to approximately describe the expansion of the contact discontinuity between the bubble and the ambient gas (Koo & McKee, 1992). In smaller halos, the gas density encountered by the outflow is on average smaller, hence the outflow velocity is higher, and the sound speed in the ambient gas is lower ($c_{s}\propto M_{h}^{1/3}$). For halos with $M(t_{9})\leq 3\times 10^{11}\;M_{\odot}$, the expansion turns supersonic at all times (see Fig. 12 right) and our treatment provides a careful description of the shell motion.

7 Discussion

7.1 Continuos vs. intermittent accretion

The basic assumption in our computations is that the outflow is powered by a black hole that has been accreting mass without interruption for many e-folding times, corresponding to $\sim 580$ Myr, in the case presented before. This is necessary to produce the $10^{9}\;M_{\odot}$ black holes observed in z = 6 QSOs if the accretion rate is limited to the Eddington rate333Defined as the critical accretion rate producing an Eddington luminosity for a radiatively efficient thin disc around a non-rotating black hole, i.e. $\dot{m}_{E}=17L_{E}/c^{2}$. and if the radiative efficiency is that of standard geometrically thin, optically thick Shakura-Sunyaev accretion discs, i.e. $\epsilon\sim 0.06-0.3$. Recent works have shown that in fact most QSOs at z = 6 may grow in this way given their measured radiative efficiencies and Eddington ratios (Trakhtenbrot et al., 2017). However, the likelihood of such a long, uninterrupted accretion at high rates is highly debated. Hydrodynamical simulations (Ciotti & Ostriker, 2007; Dubois et al., 2013) have shown that strong feedback effects occur already above a fraction of the Eddington ratio (although this is based on spherical rather than disc accretion models). These feedback effects rapidly shut down the accretion flow, which can possibly restart after a quiescence time longer than the burst of activity. The accretion and shut-down episodes can continue for many cycles, but it is clear that low duty cycles make the build-up of early SMBHs extremely challenging. A possible solution has been proposed by considering that when the accretion rate is close to Eddington, the thin-disc solution no longer applies, and the accretion should occur through radiatively inefficient “slim” discs with $\epsilon<0.05$ (Madau et al., 2014; Volonteri et al., 2015). We recall that the Soltan (1982) argument indicates that the bulk of the BH growth in the Universe must occur through radiatively efficient accretion (see also Yu & Tremaine 2002; Marconi et al. 2004). However, radiatively inefficient accretion can still power sub-populations of AGN that provide little contribution to the total mass density of local SMBHs, like e.g. QSOs at very high redshifts. In this case, the e-folding time of the black hole growth may easily be much shorter than in the standard, radiatively efficient case, whereas the radiative output would be only mildly super-Eddington (Madau et al., 2014; Volonteri et al., 2015). As an example, based on the results of Madau et al. (2014), a non-spinning black hole that is accreting at four times its Eddington rate should have a radiative efficiency of $\epsilon\sim 0.03$ and a bolometric luminosity of $L_{bol}\sim 2.6L_{E}$. In this case, the Salpeter time is a factor of 7 smaller than what we assumed for continuous accretion, as is the duration of the accretion needed to grow a $10^{4}\;M_{\odot}$ seed to a $10^{9}\;M_{\odot}$ SMBH. Since the wind power is assumed to be proportional to the radiative (bolometric) QSO output, this is 2.6 times more powerful than what we assumed for the continuous, Eddington-limited case. By recalling that in the energy-driven limit, the normalization of the bubble radius scales approximately as $R\propto L_{w,0}^{1/5}t_{Sal}^{3/5}$, where $L_{w,0}$ is the initial wind power, we may expect that, in the case of Super-Eddington accretion, the bubble radius is a factor of $2.6^{-1/5}\times 7^{3/5}\sim 2.6$ smaller. By considering that the bubble will keep expanding inertially during the non-active phases, this factor is likely an upper limit. We then conclude that intermittent super-Eddington accretion should produce bubble sizes comparable within a factor of 2 to what has been discussed in the previous sections. Interestingly, the low duty cycles implied by intermittent accretion would imply that early SMBHs are much more abundant than the active QSOs detected at z=6, and then that they should be hosted on average by smaller halos which may favour the development of supersonic outflows (see next section).

7.2 Exploring the parameter space: outflow properties as a function of halo mass, gas temperature, and BH seed

Recent hydrodynamical cosmological simulations of early BH formation suggest that part of the gas in the hosting halo may be at temperatures significantly lower than $T_{vir}$ when it is crossed by the QSO outflow (Dubois et al., 2013). Cold gas ($T<10^{5}$ K) in the form of filaments accreted by the halo itself may indeed constitute most of the gas mass in the halo. However, this cold gas is expected to fill only a minor fraction of the halo volume, which should instead be mostly occupied by a hot diffuse medium with $T\sim T_{vir}$. We recall that $T_{vir}\propto M_{h}^{2/3}$. For reference, for a halo growing to $3\times 10^{11}\;M_{\odot}$ ($10^{13}\;M_{\odot}$) by z=6, the virial temperature increases from $\sim 10^{5}$ K ($\sim 10^{6}$ K) at z=12.5 to $\sim 1.3\times 10^{6}$ K ($\sim 1.3\times 10^{7}$ K) at z=6. We explored the effects of varying the gas temperature in the halo allowing the parameter $V\equiv T_{gas}/T_{vir}$ to vary in the gas profile equations of Ferrara & Loeb (2013). As expected, for $T<T_{vir}$ the gas has higher central densities and steeper profiles. We explored the gas temperature variations in halos of different mass and the results are shown in Fig. 13 (left). The seed mass of the black hole and all the other source parameters were kept as in the previous computations. We find that the onset of the subsonic and supersonic regimes of the outflow are nearly independent on the gas temperature. This can be understood by considering that, on the one hand, the sound speed decreases with decreasing gas temperature, on the other hand, the gas central density increases, slowing down the outflow expansion. For halos with masses $M_{h}\lesssim 3\times 10^{11}\;M_{\odot}$ the outflow is always supersonic. As more massive halos are considered, the outflow spends a larger fraction of time at subsonic velocities. For $M_{h}\geq 5\times 10^{12}\;M_{\odot}$ the outflow is subsonic at any time above $t_{acc}>1$ Myr (we simply refer to this as always subsonic). When all the other parameters are fixed, the bubble radius and velocity of supersonic outflows can reach very large values for $T<T_{vir}$. As an example, for a halo with $M_{h}=3\times 10^{11}\;M_{\odot}$ the bubble radius and velocity at $t_{9}$ (z=6) are $\sim 1$ Mpc and $6000$ km s*-1* for $T\sim T_{vir}/2$, to be compared with $\sim 200$ kpc and $\sim 1000$ km s*-1* for $T=T_{vir}$. 444 As noted by Ferrara & Loeb (2013), the central gas densities may reach implausibly high values (and extremely steep profiles) for $T<<T_{vir}$. In our computations, this would also turn into implausibly high values of $v(t_{9})$ and $R(t_{9})$, so we limited our computation to $T\geq T_{vir}/4$. We also investigated what happens if the seed BH mass is left free to vary in the range $10^{3}$ to $10^{6}\;M_{\odot}$ in hosting halos whose mass may end up at z=6 in a range between $10^{11}$ and $10^{13}\;M_{\odot}$. We do not push our computations to seed masses smaller than $10^{3}\;M_{\odot}$ because this would mean starting the accretion at redshifts greater than $z=17$, where i) our extrapolations for the evolution of the host halo mass are highly uncertain and ii) the central densities of these early halos as modelled here reach implausibly high values (see Ferrara & Loeb 2013). For simplicity, we kept $T_{gas}=T_{vir}$ in the following computations. Also, we kept all the input AGN parameters as in the previous computations ($\epsilon=0.1,\lambda=1,f_{w}=0.05$) so that the QSO wind power scales linearly with the seed mass. The results are shown in Fig. 13 $(right)$. The supersonic-to-mixed and the mixed-to-subsonic outflow transitions occur at increasing halo masses as the seed BH mass increases because varying the seed mass means varying the initial conditions for the outflow expansion: larger seeds imply more powerful initial winds and a later start of the accretion (i.e. at a lower redshift) when the halo mass has grown and the central gas density is lower. This has the effect of increasing the speed of the outflow, which then needs bigger halos to be slowed down to subsonic velocities. As for the bubble “final” radius and velocity at $t_{9}$, the same conclusions reached in Fig. 3 apply: as the Salpeter time does not depend on the seed mass, the bubble radius and velocity only depend on the final BH mass, which is fixed to $10^{9}\;M_{\odot}$. Larger/smaller seeds simply mean shorter/longer accretion times to reach the same BH mass and bubble radii.

7.3 Detection of hot gas within QSO bubbles through the thermal Sunyaev-Zeldovich effect

From an observational point of view, if would be difficult to obtain a direct detection of the hot and tenuous gas filling the bubble, as this is expected to radiate inefficiently (e.g. Faucher-Giguère & Quataert 2012, Costa et al. 2014; see also next section). A promising technique is instead searching for signatures of the inverse Compton effect produced by the hot electrons in the bubble on the cosmic microwave background (CMB) photons, i.e. the thermal Sunyaev-Zeldovich (tSZ) effect.

Recent works found evidence for the tSZ effect by stacking far-IR and sub-mm data around large samples of QSOs (e.g. Ruan et al. 2015 , Crichton et al. 2016). In particular, Crichton et al. (2016) used data from both ACT and Herschel-SPIRE to build the average far-IR spectral energy distribution of SDSS QSOs in different redshift bins, from z=0.5 to z=3.5. They found a 3-4$\sigma$ evidence that the average far-IR SEDs of QSOs deviate from what is expected from pure dust emission, and that these deviations can be explained through the tSZ effect.

The spectral distortions produced by the tSZ effect are normally parameterized by the Compton y-parameter

[TABLE]

where $n_{e}$ and $T_{e}$ are the density and temperature of the hot electrons, respectively, and the integral is performed along the line of sight.

Because of the limitations in the angular resolution of ACT (FWHM = 1 arcmin, corresponding to a half-light radius of 250 kpc at the median redshift of their sample, z=1.85), Crichton et al. (2016) could only measure an integrated Compton parameter over the source solid angle, defined as

[TABLE]

where $d_{A}(z)$ is the angular diameter distance, i.e. the ratio between the object’s physical transverse size and its angular size; $P_{e}$ is the electron thermal pressure; and $E_{e}$ the electron thermal energy. For protons and electrons in thermal equilibrium the total thermal energy in the ionized gas is then

[TABLE]

where $\mu_{e}$ is the mean particle weight per electron (1.17 for primeval gas). Crichton et al. (2016) measured an average amount of thermal energy per source of $E_{th}=6.2(\pm 1.7)\times 10^{60}$ erg. This is about one dex higher than the thermal energy of the gas heated by the gravitational collapse of their dark matter halos (having $M_{h}\lesssim 5\times 10^{12}\;M_{\odot}$, as derived from clustering measurements; Shen et al. 2013; Eftekharzadeh et al. 2015), and was interpreted as evidence for thermal energy deposition from QSO outflows in their circumgalactic media, i.e. as evidence for large-scale QSO feedback.

Although we are focusing on higher-z QSOs, we are in the position of checking the amount of thermal energy deposited in our exponentially growing bubbles of hot gas at $t_{9}$, i.e. the evolutionary time where we expect to observe the source.

By considering that in energy-driven winds about half of the wind kinetic energy goes into bubble thermal energy (e.g. Weaver et al. 1977) we can integrate Eq. 4 deriving

[TABLE]

This value is remarkably similar to that measured by Crichton et al. (2016), hence suggesting that AGN outflows are closer to the energy-driven limit than to the momentum-driven limit. This conclusion is in agreement with the findings of Tombesi et al. (2015) and Feruglio et al. (2015) based on the comparison between the energy of the inner wind and that of the large-scale molecular outflow measured in two nearby QSOs.

7.4 Cooling of the hot gas in the bubble

The above computation was performed under the hypothesis that the shocked wind within the bubble does not cool and the outflow is then in the fully energy-driven limit. We will explore here under which conditions this assumption is satisfied. Assuming that the QSO wind velocity $v_{w}$ is much higher than the velocity at which the reverse shock is moving, the post-shock wind temperature, that is the bubble temperature $T_{b}$, is

[TABLE]

For such a hot gas, the only effective cooling mechanisms are free-free (thermal Bremsstrahlung) emission, inverse Compton-scattering, and, for sufficiently compact systems, pair production. The timescale for free-free cooling can be estimated as (Mo et al., 2010; Costa et al., 2014)

[TABLE]

which is obviously much longer than the Hubble time for the temperatures and densities ($n_{e}\sim 10^{-5}-10^{-7}$ cm*-3*) of our systems.

As the hot gas in the bubble is exposed to the radiation field of the QSO, the high-energy electrons in the bubble can exchange energy with the QSO photons through Compton scattering. If the temperature of the system is higher/lower than the QSO Compton temperature ($T_{C}\sim 2\times 10^{4}$ K, e.g. Sazonov et al. 2005), which depends only on the QSO spectral shape, the gas can be cooled/heated. The relevant timescale for Compton cooling/heating can be written as

[TABLE]

For our system, $T_{b}>>T_{C}$ and the bubble might cool because of inverse Compton scattering. However, when putting in the above equation the values of the bubble radius $R(t)$, and BH mass $M_{BH}(t)$ at any given accretion time, the timescales for Compton cooling are from 30 to 1000 times larger than the flow time for the whole halo mass range explored in the previous sections. We therefore conclude that Compton scattering of the QSO photons is an inefficient mechanism to cool down the bubble.

At high redshift, however, the CMB provides a large reservoir of photons available for Compton scattering because its energy density scales as $(1+z)^{4}$. Since the CMB temperature $T_{CMB}\sim 2.73(1+z)\;\rm{K}<<T_{b}$, the hot electrons in the bubble might cool efficiently by inverse Compton scattering on CMB photons (IC-CMB). Following Mo et al. (2010), for the cosmology adopted here and at $z\gtrsim 2$, the cooling time for the IC-CMB process can be approximated with

[TABLE]

where $t_{U}(z)$ is the age of the Universe at redshift $z$. By considering that $t_{U}(z)=t_{acc}+t_{U}(12.5)$ and inserting it in the above equation, we find that $t_{IC-CMB}$ is always larger than $t_{acc}$, but they are of the same order for $t_{acc}\gtrsim 100$ Myr. This means that at high-z the IC-CMB is an effective cooling mechanism for the hot gas in the bubble. This may produce significant departures from the pure energy-driven limit discussed above, but we defer a detailed treatment of the outflow behaviour to future hydrodynamic simulations which include IC-CMB cooling.

We finally discuss the possibility that the high-temperature plasma within the bubble ($T_{b}\sim 10^{10}$ K) may cool down because of pair production. Indeed, if the system is sufficiently compact, high-energy photons may produce electron-positron pairs that can slow down the fast protons in the plasma through Coulomb interactions, effectively cooling down the bubble (Begelman et al., 1987). This would happen if the system is optically thick to pair production, i.e. when the dimensionless compactness parameter $l\equiv L_{b}\sigma_{T}/Rm_{e}c^{3}>>1$ (Lightman & Zdziarski, 1987). Here $L_{b}$ is the free-free luminosity (see e.g. Eq. B1 in Faucher-Giguère & Quataert 2012) at $\sim$MeV energies of the thermal plasma in the bubble and $R$ is the bubble radius. For the values typical of our systems, $l\sim 10^{-13}$ at any given accretion time, making pair production an inefficient cooling mechanism.

7.5 Stability analysis and shell fragmentation

We provide here considerations on the stability and fragmentation of the gas shell pushed by the QSO. We consider the case of a BH starting from a seed mass of $10^{4}\;M_{\odot}$ at $z=12.5$ placed within a dark matter halo growing to $M_{h}=3\times 10^{11}\;M_{\odot}$ at $z=6$, for which the outflow is always supersonic (Fig. 12 right).

As opposed to the case treated in Sections 4 and 5 of a QSO bubble expanding in the IGM where the effects of gravity are negligible, for an outflow expanding within a dark matter halo gravity combines with the acceleration of the expanding shell in generating Rayleigh-Taylor (RT) instabilities. In particular, for a fluid with a radial acceleration $\ddot{R}$, the growth rate of the perturbations at a given spatial scale $\lambda$ is $\Gamma_{RT}=\sqrt{2\pi g^{\prime}/\lambda}$, where $g^{\prime}=\ddot{R}+g$ is the apparent gravitational, or net radial, acceleration of the system (Drazin, 2002). The growth rate of the perturbations is faster at the smallest scales, i.e. those comparable with the shell thickness $\Delta R\sim R/(3+3\mathcal{M}^{2})$, where $R$ is the bubble radius and $\mathcal{M}$ is the outflow Mach number (Weaver et al., 1977).

For the case discussed here, $g^{\prime}>0$ for $t_{acc}>5$ Myr, and the e-folding time of the small-scale perturbation growth $t_{RT}\equiv 1/\Gamma_{RT}$ becomes shorter than $t_{acc}$ soon afterwards. Therefore, RT instabilities have enough time to develop and alter the whole structure of the expanding shell. The detailed structure and fragmentation history of the shell can be studied only with high-resolution simulations. We assume here that the overall spherical structure of the shell is preserved even at late accretion times, and hence its expansion history presented in the previous section holds. The hydrodynamics simulations performed by Costa et al. (2014) assuming a constant source of energy and mass show that, after developing RT instabilities, the global structure of an expanding shell in energy-driven QSO outflows is preserved even at late times and can be described with good accuracy by simple numerical and analytic methods.

At any given time, the first fragments (clouds) separating from the expanding shell will be those whose sizes are similar to the shell thickness $\Delta R$. Once detached from the outflowing shell, the clouds are subject to the gravitational field of the halo. Because of the high contrast between the gas density within the cloud $\rho_{cl}$ and that of the hot tenuous gas within the bubble $\rho_{b}$ ( $\chi\equiv\rho_{cl}/\rho_{b}\sim 10^{4-5}$), the ram pressure exerted by the hot gas is not sufficient to slow down the cloud, which eventually falls back towards the BH in a free fall-time $\tau_{ff}=(v_{s}+\sqrt{v_{s}^{2}+2gR})/g$, where $v_{s}$, $R$, and $g$ are the shell velocity, the bubble radius and the acceleration of gravity on the shell, respectively, at the time of the cloud detachment $t_{det}$. If the free-fall time $\tau_{ff}$ is shorter than the cloud age $\tau\equiv t_{9}-t_{det}$, by the time we observe the system the cloud has already fallen back to the BH, otherwise it may be observed. For the system considered here, this means that only clouds detaching from the shell at $t_{det}>186$ Myr can be observed. At these late times, the shell thickness varies only between 2 and 3 kpc, so clouds with radii $r_{cl}=\Delta R/2\sim 1.0-1.5$ kpc are expected. Given the gas density within the cloud $\rho_{cl}=\rho_{shell}=(1+\mathcal{M}^{2})\rho_{gas}\sim 10^{-26}$ g cm*-3*, this corresponds to gas masses of $M_{cl}\sim 1-2\times 10^{6}\;M_{\odot}$ per cloud. Also, each cloud detaching at $t_{det}>186$ Myr should have the following velocity and radius by $t_{9}$: $v_{cl}(t_{9})=v_{s}(t_{det})-g\tau$ and $R_{cl}=R_{b}+v_{s}\tau-g\tau^{2}/2$. Clouds that detach from the shell at late times, should be seen close to the shell itself and with large outflowing velocities (see Fig. 14). Clouds that detach from the shell at times just after 186 Myr, should be seen close to the BH and with infalling velocities. For the case considered here, the transition between the regions populated by outflowing and infalling clouds occurs at $R_{cl}=51$ kpc, where clouds at zero velocities should be observed (see Fig. 14).

Individual clouds are expected to suffer from radiative losses and cool down following the cooling function (Cioffi & Shull, 1991)

[TABLE]

where the first relation is valid for gas made only of hydrogen and helium and the second is valid for a gas with $10^{5}<T<10^{7}$ K and metallicity $\zeta$ relative to solar. The corresponding cooling times $\tau_{cool}\sim 1.5k_{B}T_{cl}/(n_{cl}\Lambda)$ are then

[TABLE]

For the system considered here, where $T_{cl}=T_{gas}\sim 0.3-1.0\times 10^{6}$ K (assuming $T_{gas}=T_{vir}$) and $n_{cl}\sim 0.01$ cm*-3*, it follows that $t_{cool}<<t_{acc}$ only if $\zeta\gtrsim 0.05$, i.e. the clouds would cool significantly by the time we can observe them only if a significant fraction of metals is already present in the ambient gas. We speculate here that those clouds that detach early from the shell and fall back towards the black hole may both provide fuel for further BH accretion, in analogy with chaotic cold accretion models (Gaspari et al., 2017), and a reservoir of “cold”, $T\sim 10^{4}$ K gas, like the one probed by MUSE through the ubiquitous detection of giant Lyman $\alpha$ halos extending for a few tens of kpc around luminous QSOs at $3<z<4$ (Borisova et al., 2016).

As the clouds move within the hot gas in the bubble they are subject to Kelvin-Helmoltz (KH) instabilities, which can effectively remove material from the surface of the clouds until their disruption (see also Ferrara & Scannapieco 2016). The mass loss is particularly relevant at scales $\lambda\sim r_{cl}$, for which the growth timescale of the KH instability is (Murray et al., 1993)

[TABLE]

in the limit of a high-density contrast between the density of the cloud and that of the bubble, $\chi>>1$, as applicable here.

Since the cloud loses mass at a rate $\dot{M}_{cl}\sim 4\pi\rho_{cl}r_{cl}^{3}/\tau_{KH}$, the characteristic stripping time for the cloud mass is

[TABLE]

where $M_{cl}=4\pi\rho_{cl}r_{cl}^{3}/3$ is the mass of the cloud (Lin & Murray, 2000).

For the observable clouds in the system discussed above, the stripping time $\tau_{strip}$ is of the same order of the cloud age at the time of the observation. Therefore, individual clouds may survive the KH instabilities and be observed. This conclusion is reinforced by considering that the mass of those clouds that detach from the shell soon after $t_{det}$, and hence that move at relatively low velocities ($\sim 100$ km s*-1*), is close to the self-gravitating critical mass required to remain stable against KH instabilities $M_{cr}\sim(6^{1/2}\pi v_{cl}^{3})/(G^{3/2}\chi^{2}\rho_{b}^{1/2})\sim 7\times 10^{6}\;M_{\odot}$ (Murray et al., 1993). Because of the high temperature in the bubble, the clouds are subject to significant evaporation (Cowie & McKee, 1977; Marcolini et al., 2005). The mass evaporation rate can be written as

[TABLE]

where the dimensionless parameter $\sigma_{0}$ is defined as

[TABLE]

and

[TABLE]

In our case, $\sigma_{0}>>1$ and the evaporation time can be written as

[TABLE]

Therefore, all clouds should quickly evaporate if $T_{b}\sim 10^{10}$ K. However, as discussed in the previous section, for bubbles produced by high-z QSOs, IC-CMB cooling would be effective in bringing the plasma to lower temperatures. For $T_{b}\sim$ a few $\times 10^{7}K$, the evaporation time would become as long as the cloud age at $t_{9}$. Therefore, an intriguing prediction of this simple model is that a distribution of clouds around QSOs like the one shown in Fig. 14 can be only observed in systems at $z\sim 6$ and beyond. We defer a detailed treatment of the shell fragmentation and formation of gas clouds to future high-resolution hydrodynamical simulations.

8 Conclusions

We investigated the physics and time evolution of large-scale outflows produced by early QSOs powered by exponentially growing black holes. We assumed that these systems grow to $M_{BH}=10^{9}\;M_{\odot}$ by $z=6$ by accreting at the Eddington limit and converting a fixed fraction of their bolometric output into a wind. This means that the outflow source power is also growing exponentially. We first considered the cases of energy- and momentum-driven outflows expanding in a region where the gas and total mass densities are uniform and equal to the average values in the Universe at $z\geq 6$. We then extended our computations to the case of QSOs placed at the centre of early dark matter halos of different masses and starting from different seed BH masses. We made considerations on the energetics of the outflow, on the cooling of the hot gas in the QSO-inflated bubble, and on the stability and structure of the expanding gas shell. Our main results can be summarized as follows:

$\bullet$ For a SMBH/QSO growing in mass/power with an e-folding (Salpeter) time $t_{Sal}$, the late time expansion of the bubble radius is also exponential, with an e-folding time of $5t_{Sal}$ and $4t_{Sal}$ for an energy-driven and a momentum-driven outflow, respectively. In the case of a QSO expanding within a uniform density field we provided analytic solutions to the time evolution of the bubble radius.

$\bullet$ For a QSO outflow expanding within a field where the gas and total mass densities are uniform and equal to the average field values at $z>6$, the expansion of the bubble is only affected by the gas density, whereas the gravitational drag exerted by dark matter is negligible. The latter is instead relevant for outflows produced by QSOs at the centre of large dark matter halos.

$\bullet$ We considered energy-driven outflows produced by black holes growing from seeds with a mass range of $10^{3}-10^{6}M_{\odot}$ and placed within growing dark matter halos spanning a mass range of $3\times 10^{11}-10^{13}\;M_{\odot}$ at $z=6$. We followed the evolution of the source power and of the gas and dark matter density profile in the halos from the beginning of the accretion until $z=6$. For a given final BH mass ($10^{9}\;M_{\odot}$ in our case), the bubble radius and velocity at z=6 do not depend on the initial seed mass: a bubble inflated by a smaller (larger) seed simply takes more (less) time to grow to the same final value $R(t_{9})$. The final bubble radius and velocity are instead smaller for larger halo masses. At z=6, bubble radii in the range 50-180 kpc and velocities in the range 400-1000 km s*-1* are expected for QSOs hosted by halos in the mass range $3\times 10^{11}-10^{13}\;M_{\odot}$. These radius and velocity scales compare well with those measured for the outflowing gas in the well-known z=6.4 QSO SDSS J1148+5251.

$\bullet$ We assumed that the gas in the halo is at the virial temperature. For large enough halos, where the gas temperature and sound speed are higher, the expansion of the bubble may become subsonic in a given time interval. We find that for halos with $M_{h}\leq 3\times 10^{11}\;M_{\odot}$ at z=6, the outflow is always supersonic. The fraction of time spent at subsonic velocities increases for larger masses, until it is always subsonic for $M_{h}\geq 5\times 10^{12}\;M_{\odot}$ at z=6. We also explored the effects of assuming a lower ambient gas temperature, down to $T_{vir}/4$. We found that the halo mass thresholds for fully subsonic and fully supersonic outflows do not strongly depend on the assumed gas temperature and corresponding density profile. For lower temperatures and steeper profiles, the bubble radii and velocities can reach values up to 1 Mpc and a few $\times 1000$ km s*-1*, respectively.

$\bullet$ In the case of an energy-driven outflow, we computed a total thermal energy of $\sim 5\times 10^{60}$ erg contained in the bubble around the QSO. This number is in excellent agreement with the value of $(6.2\pm 1.7)\times 10^{60}$ erg per QSO as derived from a large sample of QSOs through the detection of the thermal Sunyaev-Zeldovich effect in their stacked far-IR spectra. This suggests that QSO outflows are closer to the energy-driven limit than to the momentum-driven limit.

$\bullet$ We investigated the stability of the expanding gas shell in the case of an energy-driven supersonic outflow propagating within a dark matter halo. We found that the shell is Rayleigh-Taylor unstable already at $t>5$ Myr and, by means of a simple model, we investigated the fate of the fragments detaching from the shell. We found that these fragments should rapidly evaporate because of the extremely high temperature of the hot gas bubble, unless this cools effectively. Since the only effective cooling mechanism for such a gas is inverse Compton by the CMB photons (IC-CMB), which is important only at $z\geq 6$, we speculate that such shell fragments can be observed only around high-z QSOs, where IC-CMB cooling of the bubble gas can prevent their evaporation.

$\bullet$ We finally propose that those shell fragments that have already fallen back towards the centre of the dark matter halo by the time we observe the QSO may accumulate and constitute a reservoir of relatively cold gas ($T\sim 10^{4}$ K) on scales of up to a few tens of kpc. This mechanism could explain the ubiquitous presence of such a gas observed by MUSE around $z\sim 3.5$ QSOs.

Acknowledgements.

We acknowledge stimulating discussions with M. Gaspari, L. Ciotti, M. Meneghetti, S. Ettori, and A. Negri, and a clear and constructive report from the referee. We acknowledge financial contribution from the agreement ASI-INAF I/037/12/0.