Issue 
A&A
Volume 603, July 2017



Article Number  A69  
Number of page(s)  17  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/201730457  
Published online  10 July 2017 
Exponentially growing bubbles around early supermassive black holes
^{1} INAF–Osservatorio Astronomico di Bologna, via Gobetti 93/3, 40129 Bologna, Italy
email: roberto.gilli@oabo.inaf.it
^{2} Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA
^{3} Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
Received: 18 January 2017
Accepted: 13 March 2017
We address the as yet 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_{⊙} by z = 6 by Eddingtonlimited accretion and convert 5% of their bolometric output into a wind. We first considered the case of energydriven and momentumdriven 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: for an exponentially growing power with efolding time t_{Sal}, we find that the late time expansion of the outflow radius is also exponential, with efolding time of 5t_{Sal} and 4t_{Sal} in the energydriven and momentumdriven limit, respectively. We then considered energydriven outflows produced by quasistellar objects (QSOs) at the centre 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 × 10^{11}–10^{13}M_{⊙}. These radius and velocity scales compare well with those measured for the outflowing gas in the z = 6.4 QSO SDSS J1148+5251. By the time the QSO is observed, we found that the total thermal energy injected within the bubble in the case of an energydriven outflow is E_{th} ~ 5 × 10^{60} erg. This is in excellent agreement with the value of E_{th} = (6.2 ± 1.7) × 10^{60} erg measured through the detection of the thermal SunyaevZeldovich effect around a large population of luminous QSOs at lower redshifts. This suggests that QSO outflows are closer to the energydriven limit than to the momentumdriven limit. We investigated the stability of the expanding gas shell in the case of an energydriven supersonic outflow propagating within a dark matter halo with M_{h} = 3 × 10^{11}M_{⊙} at z = 6. We found that the shell is RayleighTaylor unstable already at early times 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 if this does not cool. Since the only effective cooling mechanism for such a gas is inverse Compton by the cosmic microwave background (CMB) photons (ICCMB), which is important only at z ≥ 6, we speculate that such shell fragments may be observed only around highz QSOs, where ICCMB cooling of the bubble gas can prevent their evaporation.
Key words: black hole physics / quasars: supermassive black holes / shock waves / galaxies: highredshift
© ESO, 2017
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 Xray spectra of nearby AGN revealed the presence of ultrafast outflows (UFOs; see e.g. Tombesi et al. 2012, 2013, and references therein) produced within a few millipc of the central black hole (BH). Ultrafast outflows are made of plasma moving at 0.1c on average, and up to ~ 0.3c (Tombesi et al. 2015; Nardini et al. 2015). Although it is difficult to constrain the opening angle of UFOs, in the powerful local quasistellar object (QSO) PDS456 the remarkable P Cygni profile of the FeXXVI Kα line indicates the presence of a subrelativistic wind expanding almost spherically. On much larger scales (~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 μm and HI 21 cm 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 largescale outflowing gas. The energy required to drive largescale 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 largescale molecular outflows suggest that AGN outflows are closer to the energydriven limit than to the momentumdriven limit (see Sect. 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 SunyaevZeldovich effect on the stacked farinfrared spectra of large populations of z ~ 2 QSOs (Ruan et al. 2015; Crichton et al. 2016)
AGNdriven 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 μm line with the PdBI interferometer have revealed the presence of largescale gas outflows extending out to 30 kpc from the QSO and moving with velocities of >1000 km s^{1} (Cicone et al. 2015). This largescale outflow extends well beyond the characteristic size of the host galaxies of z ~ 6 QSOs (~ 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; FaucherGiguè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, AGNdriven 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 BHs 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, BHs 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, supercritical growth of early BHs). 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 efolding times is one of the few ways of coping with the still unsolved issue of the formation of the first supermassive BHs (SMBHs), we note here that this is inconsistent with simulations of BH growth at high Eddington ratios. Preheating instabilities that halt BH fueling already arise at Eddington ratios of ~ 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 longterm 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 selfconsistently during the outflow expansion, in contrast to most previous works.
A concordance cosmology with H_{0} = 70 km s^{1} Mpc^{1}, Ω_{m} = 0.3, Ω_{Λ} = 0.7, in agreement within the errors with the Planck 2015 results (Planck Collaboration XIII 2016), is used throughout this paper.
2. Problem setup and definitions
We consider the case of an accreting supermassive BH radiating at a given fraction λ ≡ L_{bol}/L_{E} of its Eddington luminosity L_{E} = 4πGm_{p}cM/σ_{T} = Mc^{2}/t_{E} = 1.26 × 10^{38}M erg s^{1}, where M is the BH mass in units of M_{⊙} and t_{E} = σ_{T}c/ (4πGm_{p}) = 0.45 Gyr is the Eddington time.
If a fraction ϵ of the mass m_{acc} falling onto the BH is converted into radiation, then , with ṁ_{acc} = λL_{E}/ (ϵc^{2}) = λM/ (ϵt_{E}). The BH growth rate is then Ṁ = (1 − ϵ)ṁ_{acc} = (1 − ϵ)λM/ (ϵt_{E}), from which it follows that (1)\end{document}where (2)is the Salpeter (efolding) time, and M_{0} is the mass of the BH “seed”. Under these assumptions, the bolometric QSO luminosity grows as (3)where L_{bol,0} = λ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 Ω ≈ 4π. We then investigate the evolution of astrophysical bubbles produced by an exponentially growing input energy source, (4)where \begin{lxirformule}$\dot m_{w}$\end{lxirformule} and v_{w} are the wind mass rate and velocity, respectively.
Observations of ultrafast outflows in the Xray 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 ṁ_{w} ≈ ṁ_{acc}. We further assume that the QSO wind is launched with constant velocity^{1}. Therefore, the outflow rate grows as (5)where(6)In turn, the ejected wind mass grows as (7)
Fig. 1 Classic structure of a spherical outflow produced by a fast, supersonic wind colliding into an ambient medium (see also Weaver et al. 1977; Dyson & Williams 199713; Costa et al. 2014). The flow is divided into four main regions labelled a, b, c, d. Region a is occupied by the QSO wind launched at a velocity v_{w} and carrying a kinetic power L_{w}. A reverse shock between regions a and b shocks and heats the QSO wind. It is the cooling of the shocked wind in region b (the bubble) that determines whether the outflow is energydriven (no cooling) or momentumdriven (instantaneous cooling). Region c is the thin shell of shocked ambient gas in pressure equilibrium with the gas in region b through a contact discontinuity. Region d shows the still unperturbed ambient medium with density ρ_{gas} where the outflow expands. 
3. Energydriven vs. momentumdriven 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 momentumdriven. 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 “energydriven”^{2}. 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 “momentumdriven”, 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 ρ_{gas}. Their computations can be generalized to astrophysical bubbles powered by other energy sources. In the energydriven case, the following equations can be combined to obtain the full equation of motion of the expanding shell: (8)(9)(10)Equation (8) shows the relation between thermal energy E_{th} and pressure P for a monoatomic gas in a volume V = (4π/ 3)R^{3}.
Equation (9) describes the evolution of the momentum of the shell of sweptup gas (the lefthand side is the rate of change of momentum of the shell and the righthand side represents the force acting on it).
Equation (10) is the energy equation of the system: the rate at which the thermal energy of the hot gas in the bubble (lefthand 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 momentumdriven case, E_{th} = 0 and the pressure P in Eq. (9) is given by the wind ram pressure (11)where \begin{lxirformule}$L_w$\end{lxirformule} 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 (12)where is the mean baryon density at z = 0. For H_{0} = 70 km s^{1} Mpc^{1} and Ω_{b,0} = 0.045, ρ_{b,0} = 4 × 10^{31} g cm^{3}. At z = 6 the mean baryon density is therefore ρ_{gas} ≡ ρ_{b}(6) ~ 10^{28} g cm^{3}. For an ambient gas with uniform density ρ_{gas}, the mass of the shell of sweptup gas is equal to (4π/ 3)R^{3}ρ_{gas}.
4.1. Energydriven case
Following Dyson & Williams (199713), we combine Eqs. (8)–(10) to obtain the following equation of motion of the shell: (13)This is a nonlinear 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 (14)with(15)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 BH has grown to 10^{9}M_{⊙}, i.e. M(t_{9}) ≡ 10^{9}M_{⊙}. This is the typical mass observed for bright SDSS QSOs at z ~ 6. For this computation we assumed f_{w} = 0.05, ϵ = 0.1, λ = 1, M_{0} = 10^{4}M_{⊙}. 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 ~ 10^{2}M_{⊙}, as expected for the remnants of POPIII stars (Madau & Rees 2001), to ~ 10^{4}M_{⊙}, as expected from the collapse of stellar clusters (Devecchi & Volonteri 2009), and even up to 10^{6}M_{⊙}, 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 × 10^{40} erg s^{1}, and t_{9} ~ 580 Myr.
Fig. 2 Time evolution of the radius of a hot gas bubble inflated by a wind from an early SMBH accreting at its Eddington limit (energydriven case). The following parameters are assumed: f_{w} = 0.05, λ = 1, M_{0} = 10^{4}M_{⊙} (hence L_{w,0} = 6.3 × 10^{40} erg s^{1}), ϵ = 0.1, ρ_{gas} = 10^{28} g cm^{3}. The black solid curve shows the numerical solution obtained for the equation of motion (Eq. (13)) assuming an exponentially growing wind with L_{w}(t) = L_{w,0}e^{t/tSal} and assuming R,v → 0,0 for t → 0. The analytic approximation to the numerical solution (Eq. (16)), is shown by the solid green line. The purely exponential (but unphysical) solution given by Eq. (14) is shown by the black dotted curve. The black dotted line shows the solution obtained for a constant energy source L_{w}(t) = L_{w,0} (Eq. (17)). The vertical dotted lines mark the Salpeter time t_{Sal} and the time t_{9} needed to grow the SMBH to 10^{9}M_{⊙} as labelled. The bubble radius scales as t^{3 / 5} for t ≪ t_{Sal} and e^{t/ (5tSal)} for t ≫ t_{Sal}. The red points show the results of a lowresolution simulation run with the hydrodynamical code RAMSES assuming the same input parameters. The results of the simulation are in excellent agreement with the numerical solution to Eq. (13). 
Fig. 3 Effects of varying the initial seed BH mass M_{0} on the evolution of the bubble radius following Eqs. (15) and (16) (energydriven case). All other relevant parameters are the same as assumed in Fig. 2. The three green solid curves are computed for M_{0} = 10^{2},10^{4},10^{6}M_{⊙} as labelled. The BH reaches a final mass of 10^{9}M_{⊙} at different times t_{9} (shown by the dotted lines), but the bubble radius at that time, R(t_{9}), is the same whatever the initial seed mass is. For reference, the magenta line shows the time evolution of a bubble inflated by a constant wind power of L_{w} = f_{w}L_{E} = 6.3 × 10^{45} erg s^{1}, corresponding to a 10^{9}M_{⊙} BH radiating at its Eddington limit and with f_{w} = 0.05. Because of the scaling between BH mass and the normalization of the shell radius (; Eq. (15)), this curve has a (10^{9}/ 10^{4})^{1 / 5} = 10 times higher the normalization than the black dotted line shown in Fig. 2. 
Fig. 4 Time evolution of the bubble radius (green curve) and velocity (blue curve) for an exponentially growing luminosity source powered by an exponentially growing BH (black curve) for the energydriven case. The source parameters and vertical dotted lines are as in Fig. 2. The velocity of the bubble reaches a minimum around a few Salpeter times, after which the shell starts accelerating. RayleighTaylor instabilities arise at this stage (see Sect. 5). 
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 (), assuming a stellarsized seed with M_{0} = 10^{2}M_{⊙} 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 → 0 (and v ≡ Ṙ → 0) for t → 0 and obtained a simple analytic approximation to the numerical solution by manipulating Eq. (14) and subtracting its secondorder Taylor expansion from the exponential term (16)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 199713): (17)As a further check we ran a test lowresolution 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 Sect. 4.4).
The computational box has a volume of kpc)^{3} and is characterized by a maximum resolution of ~ 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 × 10^{40} erg s^{1} and matter at a rate given by Eq. (5), with ṁ_{w,0} ~ 2.2 × 10^{4}M_{⊙} yr^{1}. In order to avoid the occurrence of diamondshaped shock fronts which can sometimes be present in similar conditions in Cartesian gridbased simulations (Tasker et al. 2008), our source is nonpointlike but distributed over a spherical volume , with Δ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 Δt are and , respectively, where and . 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 discontinuitybased. 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 efolding 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 BH mass (e.g. 10^{9}M_{⊙}) 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_{⊙} to 10^{6}M_{⊙} are shown. The bubble radius when the BH has grown to 10^{9}M_{⊙} is R(t_{9}) = 267 kpc whatever is the initial seed mass. 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 × 10^{45} erg s^{1}, corresponding to a 10^{9}M_{⊙} 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_{⊙} 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 ≡ Ṙ, 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 ~6%) by the derivative of the analytic approximation to the bubble radius given in Eq. (16): (18)From Eq. (18) it can be seen that at late times (t ≫ t_{Sal}) the bubble velocity grows exponentially according to v(t) ~ R_{0}/ (5t_{Sal}) e^{t/ (5tSal)}, i.e. it approaches the derivative of Eq. (14). Instead, at early times (t ≪ t_{Sal}), the bubble velocity follows the relation (19)i.e. we recover the v(t) ∝ t^{− 2 / 5} dependence found by Weaver et al. (1977) and Dyson & Williams (199713) 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 ~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 RayleighTaylor unstable around a few Salpeter times, i.e. at 200–300 Myr after the BH starts accreting. The shell stability is discussed in Sect. 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π/ 3)R^{3}ρ_{gas} for a uniform gas density ρ_{gas}, we obtain (20)By integrating Eq. (20) twice in time and solving for R, we obtain (21)where (22)This solution is physically meaningful, as R → 0 for t → 0.
Fig. 5 Time evolution of the shell radius for a momentumdriven outflow produced by an early SMBH accreting at its Eddington limit. Input parameters are as in Fig. 2. The black curve shows the analytic solution described by Eq. (21). The vertical dotted lines are as in Fig. 2. The two asymptotic solutions, where the shell radius scales as t^{1 / 2} for t ≪ t_{Sal} and e^{t/ (4tSal)} for t ≫ t_{Sal}, are also shown as dotted lines. 
At early times, t ≪ t_{s}, Eq. (21) reduces to (23)The above asymptotic solution correctly recovers the dependence found by Steigman (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 (24)Equation (24) shows that the efolding 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: (25)In Fig. 6 (top and middle panels) we compare the time evolution of the shell radius and velocity in the energydriven vs. momentumdriven 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 ≳ 1 Gyr the mass ejected by the BH (growing as e^{t/tSal}) will be larger than the mass of the shell (which, for instance, grows as e^{0.75t/tSal} 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 Sect. 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 Ω_{Λ} ≠ 0 can be written as (e.g. Longair 2008) (26)where tan(θ) = (Ω_{m}/ Ω_{Λ})^{1 / 2}(1 + z)^{3 / 2}.
Fig. 6 Comparison between the time evolution of the shell radius (top panel), velocity (middle), and acceleration (bottom) for an outflow produced by an exponentially growing BH in the energydriven (black curves) and momentumdriven (red curves) limits. The source input parameters are as in Figs. 2 and 5. The ratio between the two radii is plotted as a dashed curve in the top panel. The vertical dotted lines are as in Fig. 2. At any given time, the energydriven outflow is about two times more extended and faster than the momentumdriven outflow (apart from very late times, t> 1 Gyr, where the ejected wind mass becomes comparable with the shell mass and our treatment is inaccurate). 
Equation (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 ~ 12.5.
From Eq. (12) it is then easy to see that the average baryon density of the Universe decreases by a factor of ~ 7 from z = 12.5 to z = 6. In the energydriven case, the normalization of the shell radius R_{0} scales as , so at early times the bubble radius can be expected to be a factor of ~ 7^{0.2} (~ 1.5) smaller than that shown in Fig. 2. In the momentumdriven limit, where , the bubble radius would decrease by a factor of ~ 7^{0.25} (~ 1.6) at most at early times.
4.4. Effects of the external pressure
In Sect. 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 reionization of the IGM, we can assume that the temperature of the IGM around a QSO is T ≈ 2 × 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} = (γk_{B}T/μm_{p})^{1 / 2} will then be ~ 22 km s^{1} at most (assuming γ = 5 / 3 and a mean molecular weight per particle of μ = 0.59, valid for primordial gas). This is much smaller than the minimum shell velocity computed in the previous section for the energydriven and the momentumdriven 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 Sect. 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 Sects. 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, Eqs. (9) and (10), turn into (27)and (28)respectively, where is the gravitational force exerted by the total mass M_{t}(<R) enclosed within the radius of the shell of sweptup gas on the shell mass itself M_{s}(R), and ℱ_{g}Ṙ is the rate of work done to lift the shell out of the gravitational potential. For a uniform density field M_{s}(R) = (4π/ 3)R^{3}ρ_{gas}, and, for sufficiently large radii where the BH attraction can be neglected, M_{t}( <R) = (4π/ 3)R^{3}(ρ_{gas}/f_{gas}), where f_{gas} = ρ_{gas}/ρ_{m} is the gas mass fraction of the density field. By combining Eq. (8) with Eqs. (27) and (28), after some algebra we get the following equation of motion of the shell: (29)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 ∝R^{4}Ṙ on the lefthand side.
We explored under which conditions the pull exerted by gravity significantly slows down the shell expansion compared to what was shown in Sect. 4. To make a comparison with the original solutions of Weaver et al. (1977) and Dyson & Williams (199713), we first considered a constant energy source L_{w}(t) = L_{w,0}, which, in the absence of gravity, produces a bubble expanding as . We then tried a powerlaw solution to Eq. (29) of the form R(t) = R_{0}t^{α}, which leads to the relationship (30)For the above equation becomes (31)which, for t → 0, returns , i.e. at early times the solution of Weaver et al. (1977) applies: (32)Instead, for the equation can be written as (33)which, for t → ∞, returns , and the shell equation of motion at late times becomes (34)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 (∝) than in the standard nogravity case of Dyson & Williams13 (199713; ∝).
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 Eqs. (32) and (34) cross each other, i.e. R_{early}(t_{c}) = R_{late}(t_{c}). This returns (35)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 (ρ_{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 ρ_{gas} = 10^{25}g cm^{3} on the bubble radius produced by a constant source with L_{w}(t) = L_{w,0} = 6.3 × 10^{40} erg s^{1} (the full solution to Eq. (29) was obtained numerically). In this case, t_{c} ~ 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/tSal}. Similarly to Eq. (13), even Eq. (29) admits a simple analytic (exponential) solution of the form , which differs from that obtained in the case of no gravity (Eq. (14)) only by its normalization. The relation between the two normalizations and R_{0} (Eq. (15)) is (36)where . For , , 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 g cm^{3}. Therefore, in the case of a QSO placed within a medium with density similar to the average at z = 6, ρ_{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 (), even by assuming a three dex higher background density, one would get ; 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).
Fig. 7 Time evolution of the bubble radius for a constant (magenta curves) and exponentially growing (green curves) source in the energydriven limit. The solid (dashed) lines show the radius evolution considering (neglecting) the effects of gravity. Solid curves have been obtained by numerically solving Eq. (29). The assumed parameters for the exponential input energy source L_{w}(t) = L_{w,0}e^{t/tSal} are the same as those in Fig. 2. The constant source has L_{w}(t) = L_{w,0}. The vertical dotted lines are as in Fig. 2. The density of the ambient medium is assumed to be ρ_{gas} = 10^{25} g cm^{3}, i.e. about 1000 times the average of the Universe at z = 6. The powerlaw asymptotic behaviour at early (R ∝ t^{3 / 5}, magenta dashed line) and at late (R ∝ t^{1 / 5}, magenta dotted line) times can be appreciated for the constant energy source. For the exponentially growing case, such a density field reduces the size of the bubble radius by ~ 25% at late times. 
Fig. 8 Time evolution of the gas density in the bubble ρ_{bubble} for the energydriven case (region b in Fig. 1). The source input parameters are as in Fig. 2. The black curve shows the analytic solution described by Eq. (37). The vertical dotted lines are as in Fig. 2. The two asymptotic solutions, where the gas density scales as t^{− 4 / 5} for t ≪ t_{Sal} and e^{2t/ (5tSal)} for t ≫ t_{Sal}, are also shown as dotted lines. 
5. Shell stability
The basic structure of winddriven bubbles was given in Sect. 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 BH 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:

1.
The RayleighTaylor instability, which has a classic growthrate given by , where a is the acceleration and λ 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.
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 nonlinear 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 RayleighTaylor instability in the following, noting clearly that the final nonlinear state will have to be calculated by highresolution simulations.
Our case is different from the classic behaviour of windblown bubbles as the gas density within the bubble ρ_{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), ṁ = const., m ∝ t and R ∝ t^{3 / 5}, from which it follows that ρ_{bubble} ∝ m/R^{3} ∝ t^{− 4 / 5}.
To compute ρ_{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} = ηR with η 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: (37)where the numerator (ejected wind mass) has been taken from Eq. (7). The time evolution of R is given by Eq. (16).
Assuming η ≪ 1, and considering Eqs. (14) and (17), Eq. (37) reduces to at late times (t ≫ t_{Sal}), and to at early times (t ≪ t_{Sal}). Here .
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 ρ_{shell} = (1 + ℳ^{2})ρ_{gas}, where ℳ is the Mach number. For instance, the minimum shell velocity of ~ 220 km s^{1} in the energydriven case considered here (see central panel of Fig. 6) would correspond to a minimum Mach number of ~ 10 and hence to a minimum shell density of ρ_{shell} ~ 10^{26} g cm^{3}. At early times, when the shell is decelerating, the bubble density follows the ∝t^{− 4 / 5} behaviour expected for winds outflowing at a constant mass rate. In this regime, the shell is RayleighTaylor stable. At later times (t> 100 Myr), when the exponential term dominates, the shell accelerates and is expected to be RayleighTaylor unstable. As shown in Fig. 6 (bottom panel), in the energydriven limit, the shell acceleration crosses zero at t ~ 2t_{Sal}, i.e. this is the time where RayleighTaylor instabilities start developing. Interestingly, for momentumdriven outflows, this instability time is shorter as the acceleration crosses zero at about 1.4t_{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_{⊙}) or large values of η (η> 0.9, i.e. the QSO wind freely streaming up to hundred kpc distances) ρ_{bubble} might be larger than ρ_{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.
Fig. 9 Time evolution of the halo mass (in units of 10^{10}M_{⊙}, green curve) and virial radius (in units of kpc; grey curve) as compared to the evolution of the BH mass (in units of 10^{4}M_{⊙}; red curve) during the QSO accretion history. By the time the SMBH has grown to 10^{9}M_{⊙}, the halo mass and virial radius have grown by a factor of ~ 150 and ~ 10, respectively. 
Fig. 10 Density profiles of dark matter (black dashed curves) and gas (solid blue curves) at z = 12.5 and z = 6 (as labelled) for a QSO hosting halo growing to 10^{13}M_{⊙} by z = 6. The gas density is up to five orders of magnitude higher than in the uniform density field case explored in Sect. 4. 
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 ~ 6 QSOs reside is largely unknown. Based on abundancematching arguments, early SMBHs should live in 10^{13}M_{⊙} 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 ~ 6 are hosted by halos with M_{h} ~ 3 × 10^{11}M_{⊙} 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 Nbody simulation encompassing a sufficiently large volume (~ 70 Gpc^{3}) 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 ~ 10^{13}M_{⊙} at z = 6 had a mass of ~ 6 × 10^{10}M_{⊙} at z = 12.5 and should evolve into a 10^{15}M_{⊙} halo at z = 0 (i.e. into a Comalike massive galaxy cluster). We define the halo mass and virial radius in the standard way, i.e. , where ρ_{crit} = 3H^{2}(z) / (8πG) is the critical density of the Universe. The evolution of the halo mass and virial radius compared to the evolution of the BH mass are shown in Fig. 9.
We assume that, at any redshift, the dark matter density in the halo follows the NavarroFrenkWhite profile (Navarro et al. 1997) and that as the halo grows its concentration follows the concentrationmass 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 ~ 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} = (μm_{p}/ 2k_{B}) × (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} (~ 10^{12}M_{⊙} 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_{⊙} 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.
Fig. 11 Time evolution of the bubble radius produced by a growing SMBH in a dark matter halo growing to 10^{13}M_{⊙} at z = 6 in the energydriven case (black curves). The dashed curve shows the evolution of the bubble radius when only the gravitational pull of the halo is considered. The solid curve also considers the external pressure P_{ext} by the halo gas. The mass and luminosity of the accreting SMBH are as in Fig. 2. The evolution of the halo virial radius is also shown (grey curve). For comparison, the evolution of the bubble radius in the case of a uniform density field as assumed in Fig. 2 is also plotted (green curve). The vertical dotted lines are as in Fig. 2. 
We assumed for the accreting BH the same input parameters as in the previous sections (i.e. M_{0} = 10^{4}M_{⊙}, L_{w0} = 6.3 × 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 ~ 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 sweptup gas, which now follows the t^{1 / 5} law found in Sect. 4.5. Finally, for t ≳ 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 ~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 lefthand side a term equal to − 4πR^{2}P_{ext}, where P_{ext} = ρ_{gas}k_{B}T_{gas}/μm_{p}. As shown by the solid curve in Fig. 11, the external pressure reduces the bubble radius by a factor of ~ 1.6 up to a few Salpeter times (~ 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_{⊙}), the bubble reaches radii of ~ 50 kpc by z = 6, with velocities of ~ 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 BH. 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 (). For halos with M(t_{9}) ≤ 3 × 10^{11}M_{⊙}, the expansion turns supersonic at all times (see Fig. 12 right) and our treatment provides a careful description of the shell motion.
Fig. 12 Left column: time evolution of the bubble radius (upper panel) and velocity (lower panel) produced by a growing SMBH in a dark matter halo growing to 10^{13}M_{⊙} by z = 6 in the energydriven case (black curves) as compared with the evolution of the halo virial radius and sound speed of the halo gas, respectively (grey curves). For such a massive halo the outflow expansion is always subsonic. Right column: as in the left column but for a smaller halo that grows to 3 × 10^{11}M_{⊙} by z = 6. In this case the outflow expansion is always supersonic. The meaning of the vertical dotted lines is as in the previous figures. 
7. Discussion
7.1. Continuous vs. intermittent accretion
The basic assumption in our computations is that the outflow is powered by a BH that has been accreting mass without interruption for many efolding times, corresponding to ~ 580 Myr, in the case presented before. This is necessary to produce the 10^{9}M_{⊙} BHs observed in z = 6 QSOs if the accretion rate is limited to the Eddington rate^{3} and if the radiative efficiency is that of standard geometrically thin, optically thick ShakuraSunyaev accretion discs, i.e. ϵ ~ 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 shutdown episodes can continue for many cycles, but it is clear that low duty cycles make the buildup of early SMBHs extremely challenging. A possible solution has been proposed by considering that when the accretion rate is close to Eddington, the thindisc solution no longer applies, and the accretion should occur through radiatively inefficient “slim” discs with ϵ< 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 subpopulations 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 efolding time of the BH growth may easily be much shorter than in the standard, radiatively efficient case, whereas the radiative output would be only mildly superEddington (Madau et al. 2014; Volonteri et al. 2015). As an example, based on the results of Madau et al. (2014), a nonspinning BH that is accreting at four times its Eddington rate should have a radiative efficiency of ϵ ~ 0.03 and a bolometric luminosity of L_{bol} ~ 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_{⊙} seed to a 10^{9}M_{⊙} 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, Eddingtonlimited case. By recalling that in the energydriven limit, the normalization of the bubble radius scales approximately as , where L_{w,0} is the initial wind power, we may expect that, in the case of SuperEddington accretion, the bubble radius is a factor of 2.6^{− 1 / 5} × 7^{3 / 5} ~ 2.6 smaller. By considering that the bubble will keep expanding inertially during the nonactive phases, this factor is likely an upper limit. We then conclude that intermittent superEddington 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 ~ T_{vir}. We recall that . For reference, for a halo growing to 3 × 10^{11}M_{⊙} (10^{13}M_{⊙}) by z = 6, the virial temperature increases from ~ 10^{5} K (~ 10^{6} K) at z = 12.5 to ~ 1.3 × 10^{6} K (~ 1.3 × 10^{7} K) at z = 6. We explored the effects of varying the gas temperature in the halo allowing the parameter V ≡ 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 BH 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} ≲ 3 × 10^{11}M_{⊙} 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} ≥ 5 × 10^{12}M_{⊙} 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 × 10^{11}M_{⊙} the bubble radius and velocity at t_{9} (z = 6) are ~ 1 Mpc and 6000 km s^{1} for T ~ T_{vir}/ 2, to be compared with ~ 200 kpc and ~ 1000 km s^{1} for T = 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_{⊙} in hosting halos whose mass may end up at z = 6 in a range between 10^{11} and 10^{13}M_{⊙}. We do not push our computations to seed masses smaller than 10^{3}M_{⊙} 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 (ϵ = 0.1,λ = 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 supersonictomixed and the mixedtosubsonic 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_{⊙}. Larger/smaller seeds simply mean shorter/longer accretion times to reach the same BH mass and bubble radii.
Fig. 13 Left: outflow regimes as a function of halo mass at z = 6 and halo gas temperature normalized to the virial temperature. The energydriven limit is considered. The seed BH mass (10^{4}M_{⊙}) and all the other AGN input parameters are kept as in the previous computations. The transitions from supersonic to mixed (i.e. the outflow is part of the time supersonic and part of the time subsonic), and from mixed to subsonic, occur at M_{h} ~ 3 × 10^{11}M_{⊙} and M_{h} ~ 5 × 10^{12}M_{⊙}, respectively (red curves), and are nearly independent on T_{gas}. Right: same as in the left panel but as a function of halo mass at z = 6 and initial seed BH mass. Here we assume T_{gas} = T_{vir} for simplicity. 
7.3. Detection of hot gas within QSO bubbles through the thermal SunyaevZeldovich 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. FaucherGiguè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 SunyaevZeldovich (tSZ) effect.
Recent works found evidence for the tSZ effect by stacking farIR and submm 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 HerschelSPIRE to build the average farIR spectral energy distribution of SDSS QSOs in different redshift bins, from z = 0.5 to z = 3.5. They found a 3–4σ evidence that the average farIR 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 yparameter (38)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 halflight 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 (39)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 (40)where μ_{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( ± 1.7) × 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} ≲ 5 × 10^{12}M_{⊙}, 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 largescale QSO feedback.
Although we are focusing on higherz 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 energydriven 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 (41)This value is remarkably similar to that measured by Crichton et al. (2016), hence suggesting that AGN outflows are closer to the energydriven limit than to the momentumdriven 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 largescale 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 energydriven 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 postshock wind temperature, that is the bubble temperature T_{b}, is (42)For such a hot gas, the only effective cooling mechanisms are freefree (thermal Bremsstrahlung) emission, inverse Comptonscattering, and, for sufficiently compact systems, pair production. The timescale for freefree cooling can be estimated as (Mo et al. 2010; Costa et al. 2014) (43)which is obviously much longer than the Hubble time for the temperatures and densities (n_{e} ~ 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 highenergy 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} ~ 2 × 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 (44)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} ~ 2.73(1 + z)K ≪ T_{b}, the hot electrons in the bubble might cool efficiently by inverse Compton scattering on CMB photons (ICCMB). Following Mo et al. (2010), for the cosmology adopted here and at z ≳ 2, the cooling time for the ICCMB process can be approximated with (45)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} ≳ 100 Myr. This means that at highz the ICCMB is an effective cooling mechanism for the hot gas in the bubble. This may produce significant departures from the pure energydriven limit discussed above, but we defer a detailed treatment of the outflow behaviour to future hydrodynamic simulations which include ICCMB cooling.
We finally discuss the possibility that the hightemperature plasma within the bubble (T_{b} ~ 10^{10} K) may cool down because of pair production. Indeed, if the system is sufficiently compact, highenergy photons may produce electronpositron 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 ≡ L_{b}σ_{T}/Rm_{e}c^{3} ≫ 1 (Lightman & Zdziarski 1987). Here L_{b} is the freefree luminosity (see e.g. Eq. (B.1) in FaucherGiguère & Quataert 2012) at ~MeV energies of the thermal plasma in the bubble and R is the bubble radius. For the values typical of our systems, l ~ 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_{⊙} at z = 12.5 placed within a dark matter halo growing to M_{h} = 3 × 10^{11}M_{⊙} at z = 6, for which the outflow is always supersonic (Fig. 12 right).
As opposed to the case treated in Sects. 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 RayleighTaylor (RT) instabilities. In particular, for a fluid with a radial acceleration , the growth rate of the perturbations at a given spatial scale λ is , where 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 ΔR ~ R/ (3 + 3ℳ^{2}), where R is the bubble radius and ℳ is the outflow Mach number (Weaver et al. 1977).
For the case discussed here, g′> 0 for t_{acc}> 5 Myr, and the efolding time of the smallscale perturbation growth t_{RT} ≡ 1 / Γ_{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 highresolution 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 energydriven 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 Δ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 ρ_{cl} and that of the hot tenuous gas within the bubble ρ_{b} ( χ ≡ ρ_{cl}/ρ_{b} ~ 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 falltime , 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 freefall time τ_{ff} is shorter than the cloud age τ ≡ 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} = ΔR/ 2 ~ 1.0–1.5 kpc are expected. Given the gas density within the cloud ρ_{cl} = ρ_{shell} = (1 + ℳ^{2})ρ_{gas} ~ 10^{26} g cm^{3}, this corresponds to gas masses of M_{cl} ~ 1 − 2 × 10^{6}M_{⊙} 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τ and R_{cl} = R_{b} + v_{s}τ − gτ^{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) (46)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 ζ relative to solar. The corresponding cooling times τ_{cool} ~ 1.5k_{B}T_{cl}/ (n_{cl}Λ) are then (47)For the system considered here, where T_{cl} = T_{gas} ~ 0.3 − 1.0 × 10^{6} K (assuming T_{gas} = T_{vir}) and n_{cl} ~ 0.01 cm^{3}, it follows that t_{cool} ≪ t_{acc} only if ζ ≳ 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 BH 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 ~ 10^{4} K gas, like the one probed by MUSE through the ubiquitous detection of giant Lyman α halos extending for a few tens of kpc around luminous QSOs at 3 <z< 4 (Borisova et al. 2016).
Fig. 14 Schematic view (to scale) of the outflow structure around a QSO at z = 6. Regions b, c, and d are as in Fig. 1. The host halo mass was fixed to 3 × 10^{11}M_{⊙}. As RayleighTaylor instabilities develop, the gas shell loses fragments (clouds) of the size of its thickness. Depending on when a cloud detaches from the shell, it will be seen as an outflowing or inflowing cloud, where “younger” clouds have larger distances from the QSO and higher receding velocities. The length of each arrow is proportional to the cloud velocity. The dashed curve at ~ 50 kpc shows the boundary where clouds are seen as infalling or outflowing. Clouds that have detached from the shell at t_{acc}< 186 Myr have already fallen back to the BH by the time the system is observed (t_{9}), and still, they may accumulate and produce large reservoirs of cold gas around the QSO. 
As the clouds move within the hot gas in the bubble they are subject to KelvinHelmoltz (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 λ ~ r_{cl}, for which the growth timescale of the KH instability is (Murray et al. 1993) (48)in the limit of a highdensity contrast between the density of the cloud and that of the bubble, χ ≫ 1, as applicable here.
Since the cloud loses mass at a rate , the characteristic stripping time for the cloud mass is (49)where is the mass of the cloud (Lin & Murray 2000).
For the observable clouds in the system discussed above, the stripping time τ_{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 (~ 100 km s^{1}), is close to the selfgravitating critical mass required to remain stable against KH instabilities (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 (50)where the dimensionless parameter σ_{0} is defined as (51)and (52)In our case, σ_{0} ≫ 1 and the evaporation time can be written as (53)Therefore, all clouds should quickly evaporate if T_{b} ~ 10^{10} K. However, as discussed in the previous section, for bubbles produced by highz QSOs, ICCMB cooling would be effective in bringing the plasma to lower temperatures. For T_{b} ~ a few ×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 ~ 6 and beyond. We defer a detailed treatment of the shell fragmentation and formation of gas clouds to future highresolution hydrodynamical simulations.
8. Conclusions
We investigated the physics and time evolution of largescale outflows produced by early QSOs powered by exponentially growing BHs. We assumed that these systems grow to M_{BH} = 10^{9}M_{⊙} 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 momentumdriven 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 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 QSOinflated bubble, and on the stability and structure of the expanding gas shell. Our main results can be summarized as follows:

For a SMBH/QSO growing in mass/power with an efolding (Salpeter) time t_{Sal}, the late time expansion of the bubble radius is also exponential, with an efolding time of 5t_{Sal} and 4t_{Sal} for an energydriven and a momentumdriven 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.

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.

We considered energydriven outflows produced by BHs growing from seeds with a mass range of 10^{3} − 10^{6}M_{⊙} and placed within growing dark matter halos spanning a mass range of 3 × 10^{11} − 10^{13}M_{⊙} 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_{⊙} 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 × 10^{11} − 10^{13}M_{⊙}. These radius and velocity scales compare well with those measured for the outflowing gas in the wellknown z = 6.4 QSO SDSS J1148+5251.

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} ≤ 3 × 10^{11}M_{⊙} 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} ≥ 5 × 10^{12}M_{⊙} 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 ×1000 km s^{1}, respectively.

In the case of an energydriven outflow, we computed a total thermal energy of ~ 5 × 10^{60} erg contained in the bubble around the QSO. This number is in excellent agreement with the value of (6.2 ± 1.7) × 10^{60} erg per QSO as derived from a large sample of QSOs through the detection of the thermal SunyaevZeldovich effect in their stacked farIR spectra. This suggests that QSO outflows are closer to the energydriven limit than to the momentumdriven limit.

We investigated the stability of the expanding gas shell in the case of an energydriven supersonic outflow propagating within a dark matter halo. We found that the shell is RayleighTaylor 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 (ICCMB), which is important only at z ≥ 6, we speculate that such shell fragments can be observed only around highz QSOs, where ICCMB cooling of the bubble gas can prevent their evaporation.

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 ~ 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 ~ 3.5 QSOs.
Under the above assumptions, Eq. (4) implies that 2f_{w}ϵ(c/v_{w})^{2} = 1, i.e. v_{w} = 0.1c for f_{w} = 0.05 and ϵ = 0.1 as assumed here, in agreement with the velocities measured for UFOs.
Here 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).
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 ≥ T_{vir}/ 4.
Acknowledgments
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 ASIINAF I/037/12/0.
References
 Agarwal, B., Dalla Vecchia, C., Johnson, J. L., Khochfar, S., & Paardekooper, J.P. 2014, MNRAS, 443, 648 [NASA ADS] [CrossRef] [Google Scholar]
 Angulo, R. E., Springel, V., White, S. D. M., et al. 2012, MNRAS, 425, 2722 [NASA ADS] [CrossRef] [Google Scholar]
 Begelman, M. C., Sikora, M., & Rees, M. J. 1987, ApJ, 313, 689 [NASA ADS] [CrossRef] [Google Scholar]
 Borisova, E., Cantalupo, S., Lilly, S. J., et al. 2016, ApJ, 831, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Brusa, M., Bongiorno, A., Cresci, G., et al. 2015, MNRAS, 446, 2394 [NASA ADS] [CrossRef] [Google Scholar]
 Calura, F., Few, C. G., Romano, D., & D’Ercole, A. 2015, ApJ, 814, L14 [NASA ADS] [CrossRef] [Google Scholar]
 Cicone, C., Maiolino, R., Sturm, E., et al. 2014, A&A, 562, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cicone, C., Maiolino, R., Gallerani, S., et al. 2015, A&A, 574, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cioffi, D. F., & Shull, J. M. 1991, ApJ, 367, 96 [NASA ADS] [CrossRef] [Google Scholar]
 Ciotti, L., & Ostriker, J. P. 2007, ApJ, 665, 1038 [NASA ADS] [CrossRef] [Google Scholar]
 Ciotti, L., & Ostriker, J. P. 2012, in Astrophys. Space Sci. Lib. 378, eds. D.W. Kim, & S. Pellegrini, 83 [Google Scholar]
 Correa, C. A., Wyithe, J. S. B., Schaye, J., & Duffy, A. R. 2015, MNRAS, 452, 1217 [NASA ADS] [CrossRef] [Google Scholar]
 Costa, T., Sijacki, D., & Haehnelt, M. G. 2014, MNRAS, 444, 2355 [NASA ADS] [CrossRef] [Google Scholar]
 Cowie, L. L., & McKee, C. F. 1977, ApJ, 211, 135 [NASA ADS] [CrossRef] [Google Scholar]
 Cowie, L. L., Ostriker, J. P., & Stark, A. A. 1978, ApJ, 226, 1041 [NASA ADS] [CrossRef] [Google Scholar]
 Crichton, D., Gralla, M. B., Hall, K., et al. 2016, MNRAS, 458, 1478 [NASA ADS] [Google Scholar]
 De Rosa, G., Venemans, B. P., Decarli, R., et al. 2014, ApJ, 790, 145 [NASA ADS] [CrossRef] [Google Scholar]
 Dekel, A., & Birnboim, Y. 2006, MNRAS, 368, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Devecchi, B., & Volonteri, M. 2009, ApJ, 694, 302 [NASA ADS] [CrossRef] [Google Scholar]
 Drazin, P. G. 2002, Introduction to Hydrodynamic Stability (Cambridge: Cambridge University Press), 276 [Google Scholar]
 Dubois, Y., Pichon, C., Devriendt, J., et al. 2013, MNRAS, 428, 2885 [NASA ADS] [CrossRef] [Google Scholar]
 Dyson, J., & Williams, D. 1997, The Physics of the Interstellar Medium, Second Edition, Series in Astronomy and Astrophysics (CRC Press) [Google Scholar]
 Eftekharzadeh, S., Myers, A. D., White, M., et al. 2015, MNRAS, 453, 2779 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, X., Strauss, M. A., Schneider, D. P., et al. 2003, AJ, 125, 1649 [NASA ADS] [CrossRef] [Google Scholar]
 Fanidakis, N., Macciò, A. V., Baugh, C. M., Lacey, C. G., & Frenk, C. S. 2013, MNRAS, 436, 315 [NASA ADS] [CrossRef] [Google Scholar]
 FaucherGiguère, C.A., & Quataert, E. 2012, MNRAS, 425, 605 [NASA ADS] [CrossRef] [Google Scholar]
 Ferrara, A., & Loeb, A. 2013, MNRAS, 431, 2826 [NASA ADS] [CrossRef] [Google Scholar]
 Ferrara, A., & Scannapieco, E. 2016, ApJ, 833, 46 [NASA ADS] [CrossRef] [Google Scholar]
 Ferrara, A., Pettini, M., & Shchekinov, Y. 2000, MNRAS, 319, 539 [NASA ADS] [CrossRef] [Google Scholar]
 Feruglio, C., Maiolino, R., Piconcelli, E., et al. 2010, A&A, 518, L155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Feruglio, C., Fiore, F., Carniani, S., et al. 2015, A&A, 583, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaspari, M., & Sądowski, A. 2017, ApJ, 837, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Gaspari, M., Temi, P., & Brighenti, F. 2017, MNRAS, 466, 677 [NASA ADS] [CrossRef] [Google Scholar]
 Germain, J., Barai, P., & Martel, H. 2009, ApJ, 704, 1002 [NASA ADS] [CrossRef] [Google Scholar]
 Harrison, C. M., Alexander, D. M., Mullaney, J. R., et al. 2016, MNRAS, 456, 1195 [NASA ADS] [CrossRef] [Google Scholar]
 Kakiichi, K., Graziani, L., Ciardi, B., et al. 2016, MNRAS, 468, 3718 [NASA ADS] [CrossRef] [Google Scholar]
 King, A. R. 2010, MNRAS, 402, 1516 [NASA ADS] [CrossRef] [Google Scholar]
 King, A. R., Zubovas, K., & Power, C. 2011, MNRAS, 415, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Koo, B.C., & McKee, C. F. 1992, ApJ, 388, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Krause, M. G. H., & Diehl, R. 2014, ApJ, 794, L21 [NASA ADS] [CrossRef] [Google Scholar]
 Krug, H. B., Rupke, D. S. N., & Veilleux, S. 2010, ApJ, 708, 1145 [NASA ADS] [CrossRef] [Google Scholar]
 Lapi, A., Cavaliere, A., & Menci, N. 2005, ApJ, 619, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Lightman, A. P., & Zdziarski, A. A. 1987, ApJ, 319, 643 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, D. N. C., & Murray, S. D. 2000, ApJ, 540, 170 [NASA ADS] [CrossRef] [Google Scholar]
 Longair, M. S. 2008, Galaxy Formation (Berlin: Springer) [Google Scholar]
 Mac Low, M.M., & Norman, M. L. 1993, ApJ, 407, 207 [NASA ADS] [CrossRef] [Google Scholar]
 Madau, P., & Rees, M. J. 2001, ApJ, 551, L27 [NASA ADS] [CrossRef] [Google Scholar]
 Madau, P., Haardt, F., & Dotti, M. 2014, ApJ, 784, L38 [NASA ADS] [CrossRef] [Google Scholar]
 Maiolino, R., Gallerani, S., Neri, R., et al. 2012, MNRAS, 425, L66 [NASA ADS] [CrossRef] [Google Scholar]
 Marcolini, A., Strickland, D. K., D’Ercole, A., Heckman, T. M., & Hoopes, C. G. 2005, MNRAS, 362, 626 [NASA ADS] [CrossRef] [Google Scholar]
 Marconi, A., Risaliti, G., Gilli, R., et al. 2004, MNRAS, 351, 169 [NASA ADS] [CrossRef] [Google Scholar]
 Michaut, C., Cavet, C., Bouquet, S. E., Roy, F., & Nguyen, H. C. 2012, ApJ, 759, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Mo, H., van den Bosch, F. C., & White, S. 2010, Galaxy Formation and Evolution (Cambridge: Cambridge University Press) [Google Scholar]
 Morganti, R., Veilleux, S., Oosterloo, T., Teng, S. H., & Rupke, D. 2016, A&A, 593, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Murray, S. D., White, S. D. M., Blondin, J. M., & Lin, D. N. C. 1993, ApJ, 407, 588 [NASA ADS] [CrossRef] [Google Scholar]
 Nardini, E., Reeves, J. N., Gofford, J., et al. 2015, Science, 347, 860 [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [NASA ADS] [CrossRef] [Google Scholar]
 Nayakshin, S., & Zubovas, K. 2012, MNRAS, 427, 372 [NASA ADS] [CrossRef] [Google Scholar]
 Novak, G. S., Ostriker, J. P., & Ciotti, L. 2011, ApJ, 737, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Ostriker, J. P., & McKee, C. F. 1988, Rev. Mod. Phys., 60, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Overzier, R. A. 2016, A&ARv, 24, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Pezzulli, E., Valiante, R., & Schneider, R. 2016, MNRAS, 458, 3047 [NASA ADS] [CrossRef] [Google Scholar]
 Pittard, J. M. 2013, MNRAS, 435, 3600 [NASA ADS] [CrossRef] [Google Scholar]
 Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ruan, J. J., McQuinn, M., & Anderson, S. F. 2015, ApJ, 802, 135 [NASA ADS] [CrossRef] [Google Scholar]
 Sazonov, S. Y., Ostriker, J. P., Ciotti, L., & Sunyaev, R. A. 2005, MNRAS, 358, 168 [NASA ADS] [CrossRef] [Google Scholar]
 Scannapieco, E., & Oh, S. P. 2004, ApJ, 608, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Shen, Y. 2016, ApJ, 817, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Shen, Y., McBride, C. K., White, M., et al. 2013, ApJ, 778, 98 [NASA ADS] [CrossRef] [Google Scholar]
 Soltan, A. 1982, MNRAS, 200, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Steigman, G., Strittmatter, P. A., & Williams, R. E. 1975, ApJ, 198, 575 [NASA ADS] [CrossRef] [Google Scholar]
 Tasker, E. J., Brunino, R., Mitchell, N. L., et al. 2008, MNRAS, 390, 1267 [NASA ADS] [CrossRef] [Google Scholar]
 Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tombesi, F., Cappi, M., Reeves, J. N., & Braito, V. 2012, MNRAS, 422, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Tombesi, F., Cappi, M., Reeves, J. N., et al. 2013, MNRAS, 430, 1102 [Google Scholar]
 Tombesi, F., Meléndez, M., Veilleux, S., et al. 2015, Nature, 519, 436 [Google Scholar]
 Trakhtenbrot, B., Volonteri, M., & Natarajan, P. 2017, ApJ, 836, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Valiante, R., Schneider, R., Volonteri, M., & Omukai, K. 2016, MNRAS, 457, 3356 [NASA ADS] [CrossRef] [Google Scholar]
 Venemans, B. P., Walter, F., Decarli, R., et al. 2017, ApJ, 837, 146 [NASA ADS] [CrossRef] [Google Scholar]
 Vishniac, E. T. 1983, ApJ, 274, 152 [NASA ADS] [CrossRef] [Google Scholar]
 Volonteri, M., & Bellovary, J. 2012, Rep. Prog. Phys., 75, 124901 [NASA ADS] [CrossRef] [Google Scholar]
 Volonteri, M., Silk, J., & Dubus, G. 2015, ApJ, 804, 148 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, R., Wagg, J., Carilli, C. L., et al. 2013, ApJ, 773, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, ApJ, 218, 377 [NASA ADS] [CrossRef] [Google Scholar]
 Willott, C. J., Albert, L., Arzoumanian, D., et al. 2010, AJ, 140, 546 [NASA ADS] [CrossRef] [Google Scholar]
 Yu, Q., & Tremaine, S. 2002, MNRAS, 335, 965 [NASA ADS] [CrossRef] [Google Scholar]
 Zakamska, N. L., Hamann, F., Pâris, I., et al. 2016, MNRAS, 459, 3144 [NASA ADS] [CrossRef] [Google Scholar]
 Zubovas, K., & King, A. 2012, ApJ, 745, L34 [Google Scholar]
All Figures
Fig. 1 Classic structure of a spherical outflow produced by a fast, supersonic wind colliding into an ambient medium (see also Weaver et al. 1977; Dyson & Williams 199713; Costa et al. 2014). The flow is divided into four main regions labelled a, b, c, d. Region a is occupied by the QSO wind launched at a velocity v_{w} and carrying a kinetic power L_{w}. A reverse shock between regions a and b shocks and heats the QSO wind. It is the cooling of the shocked wind in region b (the bubble) that determines whether the outflow is energydriven (no cooling) or momentumdriven (instantaneous cooling). Region c is the thin shell of shocked ambient gas in pressure equilibrium with the gas in region b through a contact discontinuity. Region d shows the still unperturbed ambient medium with density ρ_{gas} where the outflow expands. 

In the text 
Fig. 2 Time evolution of the radius of a hot gas bubble inflated by a wind from an early SMBH accreting at its Eddington limit (energydriven case). The following parameters are assumed: f_{w} = 0.05, λ = 1, M_{0} = 10^{4}M_{⊙} (hence L_{w,0} = 6.3 × 10^{40} erg s^{1}), ϵ = 0.1, ρ_{gas} = 10^{28} g cm^{3}. The black solid curve shows the numerical solution obtained for the equation of motion (Eq. (13)) assuming an exponentially growing wind with L_{w}(t) = L_{w,0}e^{t/tSal} and assuming R,v → 0,0 for t → 0. The analytic approximation to the numerical solution (Eq. (16)), is shown by the solid green line. The purely exponential (but unphysical) solution given by Eq. (14) is shown by the black dotted curve. The black dotted line shows the solution obtained for a constant energy source L_{w}(t) = L_{w,0} (Eq. (17)). The vertical dotted lines mark the Salpeter time t_{Sal} and the time t_{9} needed to grow the SMBH to 10^{9}M_{⊙} as labelled. The bubble radius scales as t^{3 / 5} for t ≪ t_{Sal} and e^{t/ (5tSal)} for t ≫ t_{Sal}. The red points show the results of a lowresolution simulation run with the hydrodynamical code RAMSES assuming the same input parameters. The results of the simulation are in excellent agreement with the numerical solution to Eq. (13). 

In the text 
Fig. 3 Effects of varying the initial seed BH mass M_{0} on the evolution of the bubble radius following Eqs. (15) and (16) (energydriven case). All other relevant parameters are the same as assumed in Fig. 2. The three green solid curves are computed for M_{0} = 10^{2},10^{4},10^{6}M_{⊙} as labelled. The BH reaches a final mass of 10^{9}M_{⊙} at different times t_{9} (shown by the dotted lines), but the bubble radius at that time, R(t_{9}), is the same whatever the initial seed mass is. For reference, the magenta line shows the time evolution of a bubble inflated by a constant wind power of L_{w} = f_{w}L_{E} = 6.3 × 10^{45} erg s^{1}, corresponding to a 10^{9}M_{⊙} BH radiating at its Eddington limit and with f_{w} = 0.05. Because of the scaling between BH mass and the normalization of the shell radius (; Eq. (15)), this curve has a (10^{9}/ 10^{4})^{1 / 5} = 10 times higher the normalization than the black dotted line shown in Fig. 2. 

In the text 
Fig. 4 Time evolution of the bubble radius (green curve) and velocity (blue curve) for an exponentially growing luminosity source powered by an exponentially growing BH (black curve) for the energydriven case. The source parameters and vertical dotted lines are as in Fig. 2. The velocity of the bubble reaches a minimum around a few Salpeter times, after which the shell starts accelerating. RayleighTaylor instabilities arise at this stage (see Sect. 5). 

In the text 
Fig. 5 Time evolution of the shell radius for a momentumdriven outflow produced by an early SMBH accreting at its Eddington limit. Input parameters are as in Fig. 2. The black curve shows the analytic solution described by Eq. (21). The vertical dotted lines are as in Fig. 2. The two asymptotic solutions, where the shell radius scales as t^{1 / 2} for t ≪ t_{Sal} and e^{t/ (4tSal)} for t ≫ t_{Sal}, are also shown as dotted lines. 

In the text 
Fig. 6 Comparison between the time evolution of the shell radius (top panel), velocity (middle), and acceleration (bottom) for an outflow produced by an exponentially growing BH in the energydriven (black curves) and momentumdriven (red curves) limits. The source input parameters are as in Figs. 2 and 5. The ratio between the two radii is plotted as a dashed curve in the top panel. The vertical dotted lines are as in Fig. 2. At any given time, the energydriven outflow is about two times more extended and faster than the momentumdriven outflow (apart from very late times, t> 1 Gyr, where the ejected wind mass becomes comparable with the shell mass and our treatment is inaccurate). 

In the text 
Fig. 7 Time evolution of the bubble radius for a constant (magenta curves) and exponentially growing (green curves) source in the energydriven limit. The solid (dashed) lines show the radius evolution considering (neglecting) the effects of gravity. Solid curves have been obtained by numerically solving Eq. (29). The assumed parameters for the exponential input energy source L_{w}(t) = L_{w,0}e^{t/tSal} are the same as those in Fig. 2. The constant source has L_{w}(t) = L_{w,0}. The vertical dotted lines are as in Fig. 2. The density of the ambient medium is assumed to be ρ_{gas} = 10^{25} g cm^{3}, i.e. about 1000 times the average of the Universe at z = 6. The powerlaw asymptotic behaviour at early (R ∝ t^{3 / 5}, magenta dashed line) and at late (R ∝ t^{1 / 5}, magenta dotted line) times can be appreciated for the constant energy source. For the exponentially growing case, such a density field reduces the size of the bubble radius by ~ 25% at late times. 

In the text 
Fig. 8 Time evolution of the gas density in the bubble ρ_{bubble} for the energydriven case (region b in Fig. 1). The source input parameters are as in Fig. 2. The black curve shows the analytic solution described by Eq. (37). The vertical dotted lines are as in Fig. 2. The two asymptotic solutions, where the gas density scales as t^{− 4 / 5} for t ≪ t_{Sal} and e^{2t/ (5tSal)} for t ≫ t_{Sal}, are also shown as dotted lines. 

In the text 
Fig. 9 Time evolution of the halo mass (in units of 10^{10}M_{⊙}, green curve) and virial radius (in units of kpc; grey curve) as compared to the evolution of the BH mass (in units of 10^{4}M_{⊙}; red curve) during the QSO accretion history. By the time the SMBH has grown to 10^{9}M_{⊙}, the halo mass and virial radius have grown by a factor of ~ 150 and ~ 10, respectively. 

In the text 
Fig. 10 Density profiles of dark matter (black dashed curves) and gas (solid blue curves) at z = 12.5 and z = 6 (as labelled) for a QSO hosting halo growing to 10^{13}M_{⊙} by z = 6. The gas density is up to five orders of magnitude higher than in the uniform density field case explored in Sect. 4. 

In the text 
Fig. 11 Time evolution of the bubble radius produced by a growing SMBH in a dark matter halo growing to 10^{13}M_{⊙} at z = 6 in the energydriven case (black curves). The dashed curve shows the evolution of the bubble radius when only the gravitational pull of the halo is considered. The solid curve also considers the external pressure P_{ext} by the halo gas. The mass and luminosity of the accreting SMBH are as in Fig. 2. The evolution of the halo virial radius is also shown (grey curve). For comparison, the evolution of the bubble radius in the case of a uniform density field as assumed in Fig. 2 is also plotted (green curve). The vertical dotted lines are as in Fig. 2. 

In the text 
Fig. 12 Left column: time evolution of the bubble radius (upper panel) and velocity (lower panel) produced by a growing SMBH in a dark matter halo growing to 10^{13}M_{⊙} by z = 6 in the energydriven case (black curves) as compared with the evolution of the halo virial radius and sound speed of the halo gas, respectively (grey curves). For such a massive halo the outflow expansion is always subsonic. Right column: as in the left column but for a smaller halo that grows to 3 × 10^{11}M_{⊙} by z = 6. In this case the outflow expansion is always supersonic. The meaning of the vertical dotted lines is as in the previous figures. 

In the text 
Fig. 13 Left: outflow regimes as a function of halo mass at z = 6 and halo gas temperature normalized to the virial temperature. The energydriven limit is considered. The seed BH mass (10^{4}M_{⊙}) and all the other AGN input parameters are kept as in the previous computations. The transitions from supersonic to mixed (i.e. the outflow is part of the time supersonic and part of the time subsonic), and from mixed to subsonic, occur at M_{h} ~ 3 × 10^{11}M_{⊙} and M_{h} ~ 5 × 10^{12}M_{⊙}, respectively (red curves), and are nearly independent on T_{gas}. Right: same as in the left panel but as a function of halo mass at z = 6 and initial seed BH mass. Here we assume T_{gas} = T_{vir} for simplicity. 

In the text 
Fig. 14 Schematic view (to scale) of the outflow structure around a QSO at z = 6. Regions b, c, and d are as in Fig. 1. The host halo mass was fixed to 3 × 10^{11}M_{⊙}. As RayleighTaylor instabilities develop, the gas shell loses fragments (clouds) of the size of its thickness. Depending on when a cloud detaches from the shell, it will be seen as an outflowing or inflowing cloud, where “younger” clouds have larger distances from the QSO and higher receding velocities. The length of each arrow is proportional to the cloud velocity. The dashed curve at ~ 50 kpc shows the boundary where clouds are seen as infalling or outflowing. Clouds that have detached from the shell at t_{acc}< 186 Myr have already fallen back to the BH by the time the system is observed (t_{9}), and still, they may accumulate and produce large reservoirs of cold gas around the QSO. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.