A high baryon fraction in massive haloes at z~3

TL;DR

This study examines the baryon content in the circumgalactic medium of massive haloes at z~3, using Lyα nebulae observations, and finds evidence for a high baryon fraction close to the cosmic value, challenging some feedback models.

Contribution

It provides new constraints on the baryon fraction in high-redshift galaxy haloes by modeling Lyα emission, considering different emission mechanisms and their implications.

Findings

CGM baryon fractions are close to the cosmic value, with cold gas occupying less than 1% of volume.

Reconciliation of low baryon fractions with observations requires more massive haloes or colder gas than standard models predict.

Scattered QSO broad line region photons imply a baryon fraction of at least 70%, indicating substantial baryonic content in these haloes.

Abstract

We investigate the baryon content of the circumgalactic medium (CGM) within the virial radius of $M_h \sim 10^{12} \; M_\odot$ haloes at z ~ 3, by modelling the surface brightness profile of the giant Ly$\alpha$ nebulae recently discovered by MUSE around bright QSOs at this redshift. We initially assume fluorescent emission from cold photo-ionized gas confined by the pressure of a hot halo. Acceptable CGM baryon fractions (equal or smaller than the cosmological value) require that the cold gas occupies $\lesssim$ 1% of the volume, but is about as massive as the hot gas. CGM baryon fractions as low as 30% of the cosmic value, as predicted by some strongly ejective feedback models at this redshift, are not easy to reconcile with observations, under our assumptions, unless both the QSO-hosting haloes at $z\sim3$ are more massive than recent BOSS estimates based on clustering and the photo-ionized gas is colder than expected in a standard QSO ionizing radiation field. We also consider the option that the emission is dominated by photons scattered from the QSO broad line region. In this scenario, a very stringent lower limit to the baryon fraction can be obtained under the extreme assumption of optically thin scattering. We infer in this case a baryon fraction of at least 70% of the cosmic value, for fiducial parameters. Lower values require halo masses or gas temperatures different than expected, or that some mechanism keeps the cold gas systematically over-pressured with respect to the ambient medium.

Figures

Equations

$$\Sigma =\int f_{V}B(T)n^2,$$

$$B(T)\equiv \frac{E_{\text{Ly}\alpha }}{4\pi }\left(1+2\frac{n_{\text{He}}}{n_{\text{H}}}\right){\alpha }_{\text{Ly}\alpha }^{\text{eff}}(T),$$

$${\alpha }_{\text{Ly}\alpha }^{\text{eff}}=1.6\times 10^{-13}{T}_4^{-0.96}{\text{cm}}^3{\text{s}}^{-1}$$

$$⟨n⟩=f_{V}n,$$

$$C=f_V^{-1}$$

$$\Sigma (R)=3.9\times 10^{-4}{\left(\frac{R}{10\text{kpc}}\right)}^{-1.5}\text{erg}{\text{s}}^{-1}{\text{cm}}^{-2}{\text{srd}}^{-1}.$$

$$n(r)=0.15f_V^{-1/2}{T}_4^{t/2}{\left(\frac{r}{10\text{kpc}}\right)}^{-1.25}{\text{cm}}^{-3},$$

$$\left.\begin{array}[]{c}n_{0}=f_{V}^{-1/2}\left(\dfrac{\Sigma_{0}}{2B(T)\chi(1+\beta)R_{0}}\right)^{1/2}\;,\\ \\ \gamma=\dfrac{1+\beta}{2}\;,\end{array}\right.$$

$$\chi (a)\equiv {\int }_0^{+\infty }{\left(1+x^2\right)}^{-a/2}dx.$$

$$M_{\text{cold}}={\int }_0^{r_{\text{vir}}}4\pi f_V\frac{m_p}{X}{r}^{2}n(r)dr,$$

$$r_{\text{vir}}=\frac{1}{1+z}{\left(\frac{G{M}_{\text{halo}}}{100{\Omega }_{\text{m}}{H}_0^2}\right)}^{1/3}=74M_{12}^{1/3}\text{kpc},$$

$$\left.\begin{array}[]{ll}f_{\textrm{CGM,cold}}&\equiv\dfrac{M_{\textrm{cold}}(r_{\textrm{vir}})}{(\Omega_{b}/\Omega_{m})M_{\textrm{halo}}}\\ &\\ &=f_{1}\;f_{V}^{1/2}T_{4}^{t/2}M_{12}^{-\gamma/3}\;,\end{array}\right.$$

$$\epsilon =\frac{n_{\text{hot}}}{n_{\text{cold}}}=\frac{T_{\text{cold}}}{T_{\text{hot}}},$$

$$n_{\text{hot}}(r)=\epsilon n_0{\left(\frac{r}{R_0}\right)}^{-\gamma }.$$

$$T_{\text{hot}}=\frac{1}{\gamma }\frac{\mu m_p}{k_B}\frac{G{M}_{\text{halo}}}{r_{\text{vir}}},$$

$$\begin{array}{rl}\epsilon & =\frac{\gamma k_B{T}_{\text{cold}}}{\mu m_p}{\left(1+z\right)}^{-1}{\left(100{\Omega }_{\text{m}}{H}_0^2\right)}^{-1/3}{\left(G{M}_{\text{halo}}\right)}^{-2/3}\\ & ={\epsilon }_1{T}_4{M}_{12}^{-2/3},\end{array}$$

$${\epsilon }_1=1.0\times 10^{-2}\frac{\gamma }{1+z}=3.0\times 10^{-3}$$

$$\frac{M_{\text{hot}}}{M_{\text{cold}}}=\epsilon \frac{1-f_V}{f_V}={\epsilon }_1{f}_V^{-1}{T}_4{M}_{12}^{-2/3},$$

$$\left.\begin{array}[]{ll}f_{\textrm{CGM,tot}}&\equiv\dfrac{M_{\textrm{cold}+\textrm{hot}}(r_{\textrm{vir}})}{(\Omega_{b}/\Omega_{m})M_{\textrm{halo}}}\\ &\\ &=f_{1}(1+\varepsilon_{1}f_{V}^{-1}T_{4}M_{12}^{-2/3})f_{V}^{1/2}T_{4}^{t/2}M_{12}^{-\gamma/3}\;,\end{array}\right.$$

$$T_{\text{eq}}=\frac{2h{\nu }_{\text{LL}}}{3(p+2)k_B}=2.8\times 10^4{\left(\frac{p+2}{3.7}\right)}^{-1}\text{K},$$

$$A(T)\equiv \frac{n_{He}}{n_H}\frac{{\alpha }_{\text{rec}}^{He}(T)}{{\alpha }_{\text{rec}}^H(T)}.$$

$$T_{\text{eq}}=5.6\times 10^4\text{K}.$$

$$f_{\text{min}}=2f_1{\epsilon }_1^{1/2}{T}_4^{(1+t)/2}{M}_{12}^{-(1+\gamma )/3}$$

$$f_{V,\text{eq}}={\epsilon }_1{T}_4{M}_{12}^{-2/3}.$$

$$\left.\begin{array}[]{ll}\Gamma_{\textrm{ion}}(r)&=\dfrac{1}{4\pi r^{2}}\displaystyle\int_{\nu_{\textrm{LL}}}^{+\infty}\dfrac{L_{\nu}}{h\nu}\sigma^{\textrm{ion}}(\nu)\;d\nu\\ &\\ &=\dfrac{1}{p+3}\dfrac{L_{\nu,\textrm{LL}}\;\sigma^{\textrm{ion}}_{\textrm{LL}}}{4\pi hr^{2}}\\ &\\ &=8.8\times 10^{-7}\xi\left(\dfrac{r}{10\;\textrm{kpc}}\right)^{-2}\;\textrm{s}^{-1}\;,\end{array}\right.$$

$$\xi \equiv \left(\frac{L_{\nu ,\text{LL}}}{5.2\times 10^{31}\text{erg}{\text{s}}^{-1}{\text{Hz}}^{-1}}\right){\left(\frac{p+3}{4.7}\right)}^{-1},$$

$$\left.\begin{array}[]{ll}\Gamma_{\textrm{rec}}&=\alpha_{\textrm{rec}}(T)\left(1+2\dfrac{n_{He}}{n_{H}}\right)n\\ &\\ &=4.8\times 10^{-13}\;T_{4}^{-0.75}\left(\dfrac{n}{\textrm{cm}^{-3}}\right)\;\textrm{s}^{-1}\;,\end{array}\right.$$

$$x_{HI}=7.9\times 10^{-8}{f}_V^{-1/2}{\xi }^{-1}{T}_4^{-0.27}{\left(\frac{r}{10\text{kpc}}\right)}^{0.75}.$$

$$\left.\begin{array}[]{l}n_{HI}=1.1\times 10^{-8}\;f_{V}^{-1}\;\xi^{-1}\;T_{4}^{0.21}\left(\dfrac{r}{10\;\textrm{kpc}}\right)^{-0.5}\;\textrm{cm}^{-3}\;,\end{array}\right.$$

$$⟨n_{HI}⟩=f_{V}n=1.1\times 10^{-8}{\xi }^{-1}{T}_4^{0.21}{\left(\frac{r}{10\text{kpc}}\right)}^{-0.5}{\text{cm}}^{-3},$$

Full text

††pagerange: A high baryon fraction in massive haloes at $z\sim 3$–References††pubyear: 2019

1 Introduction

Modern galaxy formation models rely on means to keep outside of the star-forming body of galaxies most of the baryons theoretically associated with their haloes. Different models achieve this result in different ways. In particular, the so-called missing baryons can be stored either relatively close to galaxies (in the circumgalactic medium, or CGM) or very far away from them (in the intergalactic medium, or IGM).111Here and afterwards, we define the boundary between the CGM and the IGM to be the virial radius. As a consequence, observational estimates of the (baryonic) mass of the CGM, as a function of halo mass and redshift, are of great importance to further constrain existing models and guide advancement of galaxy formation theory.

Shock heating as a consequence of gravitational collapse was first identified as a possible means to bring most of the baryons falling into the potential wells of forming galaxies to the virial temperature, leading to the formation of a low-density, hot and massive CGM (Rees & Ostriker 1977; White & Rees 1978). Whenever the dynamical time is long compared to the cooling time, however, the process is largely compensated by radiative cooling, which would lead to the collapse of most of the baryons towards the centre and to excessive star formation compared to observations (e.g. White & Frenk 1991). In some conditions, a virial shock does not even develop, or it is unstable, and filaments of cold gas can penetrate deep into the halo, possibly fuelling substantial star formation at the centre (e.g. Birnboim & Dekel 2003; Binney 2004; Kereš et al. 2005). The conditions for efficient radiative cooling are far from trivial, as they crucially depend on metallicity (e.g. Sutherland & Dopita 1993), on the intensity and shape of a meta-galactic ionizing background (e.g. Efstathiou 1992; Wiersma, Schaye & Smith 2009), on the time-scales needed to ensure collisional or photo-ionization equilibrium (Gnat 2017) and on whether local ionizing sources are consistently taken into account (Cantalupo 2010). Some of these effects have been shown to be important, especially at low masses (e.g. Kannan et al. 2016), while, at least at high masses, it is most of the times assumed that catastrophic cooling must be compensated by fresh injection of energy, which either provides additional heating to the CGM, or mechanically removes a substantial fraction of the CGM from galaxy haloes. The second option can be referred to as ejective feedback.

The most commonly advocated energy sources are supernovae (SNe). First proposed by Dekel & Silk (1986), SN-driven galactic winds can, if certain conditions are met, eject a large fraction of the baryons out of the potential wells of galaxies and thus solve the overcooling problem. For this reason, they have now become a key ingredient of most semi-analytical galaxy formation models (e.g. Henriques et al. 2013 and references therein), as well as hydrodynamical cosmological simulations (e.g. Springel & Hernquist 2003; Oppenheimer & Davé 2008; Dalla Vecchia & Schaye 2012). Expulsion of baryons due to SN feedback is expected to be particularly effective for low-mass galaxies, due to their low escape speed. Crucially, however, the process will impact the baryon content of high-mass haloes as well. Baryons ejected from the shallow potential wells of dwarf galaxies at very high redshift can reach considerably large distances. As a consequence, more massive haloes, formed later from the merging of the low-mass ones, would find themselves relatively deficient in baryons already at virialization. The baryon fraction in these haloes would then slowly increase towards the cosmological value, as the gas ejected at early times is gradually accreted back (‘re-incorporated’). With some dependence on implementation details, this scenario provides predictions for the baryon fraction within the virial radius as a function of halo mass and redshift (e.g. Davé 2009; Muratov et al. 2015; Liang et al. 2016; Mitchell et al. 2018), thus indicating a way to test or refine the paradigm observationally.

There are reasons to expect that reality may be more complicated than the simple scenario above. Models based on the described paradigm usually require that the energy produced by SN explosions (as well as by the stellar winds of young stars) must be transferred to the surrounding medium with an efficiency close to unity. Solutions must then be envisaged to stem the radiative losses that are expected on theoretical grounds (e.g. Weaver et al. 1977; Gatto et al. 2017; Naab & Ostriker 2017 and references therein). Even at maximum efficiency, it is widely believed that SN feedback is insufficient to explain the high-mass end of the galaxy stellar mass function, so that additional mechanisms are often advocated to regulate the evolution of galaxies in that regime, the most common being feedback from active galactic nuclei (AGN, e.g. Ciotti & Ostriker 2001; Brighenti & Mathews 2006; Beckmann et al. 2017). AGN feedback has been modelled in a variety of fashions: thermal, kinetic, radiative, continuous or intermittent, isotropic or anisotropic, with a variety of assumed efficiencies and a different impact on the baryon content of the CGM (e.g. Sijacki et al. 2007; Steinborn et al. 2015; Choi et al. 2015; Weinberger et al. 2017). In addition to the above, other mechanisms are being investigated which can influence accretion on to galaxies and star formation, while having a variety of impacts on the baryon content of haloes. These include, among others, various kinds of radiation pressure (e.g. Roškar et al. 2014; Kimm et al. 2018), cosmic rays (e.g. Samui, Subramanian & Srianand 2018; Jacob et al. 2018), thermal conduction (Nipoti & Binney 2007) and photo-ionization by local sources (Cantalupo 2010; Kannan et al. 2016).

Observational determinations of the baryonic mass in the CGM around galaxies, as a function of halo mass and redshift, are of great help to confirm galaxy formation models or direct theoretical investigations further. At low redshift, a wealth of information comes from the X-ray bremsstrahlung of hot (close to virial temperature) gas in groups and clusters, as well as around galaxies (e.g. Laganá et al. 2013; Bogdán et al. 2017; Li et al. 2018). Some of these observations (e.g. Lovisari, Reiprich & Schellenberger 2015; Eckert et al. 2016) suggest that at $z=0$ the fraction of baryons (at least of those in the hot phase) found within haloes increases with halo mass, an input that has already proven useful to calibrate feedback models in some cosmological simulations (e.g. Le Brun et al. 2014; Pillepich et al. 2018). In addition to observations of the hot gas, an increasing amount of data is becoming available for the colder phase ($T\sim 10^{4}\;\textrm{K}$), entrained within hot haloes and mostly observable in absorption along the sightlines of background QSOs. Estimating the mass of this component is more difficult, because of the intrinsically serendipitous nature of absorption studies, as well as the need to resort to photo-ionization modelling to convert the observed column densities of the absorbing species into a total gas mass. It has however been suggested that the mass in the cold phase can in some cases also be significant (e.g. Werk et al. 2014; Stern et al. 2016; Zahedy et al. 2019).

The baryon content of the CGM of high $z$ galaxies is of special value to galaxy formation theory, as its determination would provide the most direct information on the nature of feedback mechanisms – and in general on the processes determining galaxy formation – in the earliest stages of their development. In particular, at early epochs, the amount of retained baryons is expected to change, with halo mass and redshift, in a way that is mostly determined by the efficiency of SN feedback (e.g. Henriques et al. 2019). Therefore, measurements of the CGM baryon fraction at high $z$ can, at least in principle, add useful constraints to the choice of the parameters, such as the fraction of SN energy that is converted into gas kinetic energy, that eventually modulate the ejection of baryons from galaxy haloes. Unfortunately, secure measurements of bremsstrahlung emission from hot haloes at high $z$ require very long integration times with the current instrumentation. Progress is expected from complementary probes, such as the Sunyaev-Zeldovich effect. This is unaffected by cosmological dimming and has already provided interesting indications at intermediate redshifts at the mass scale of galaxy clusters (e.g. Chiu et al. 2018 and references therein). The mass of the hot gas can also be estimated indirectly from that of the cold gas, when observations of the latter are available, based on the assumption of pressure balance between the two phases. The ansatz of pressure equilibrium was initially proposed by Spitzer (1956), who successfully applied it to the halo of the Milky Way, and was proven applicable up to the scale of galaxy clusters (e.g. Churazov et al. 2001; Zhuravleva et al. 2016). Absorption data provide some constraints on the cold gas around star-forming galaxies (Rudie et al. 2012) and QSOs (Prochaska, Hennawi & Simcoe 2013) out to $z\sim 2$ and suggest the presence of significant amounts of cold gas at projected distances comparable to virial radius, though a quantitative determination of the total gas mass is not straightforward for the reasons mentioned above. On the other hand, the cold phase of the CGM at $z\gtrsim 2$ is sufficiently bright in fluorescent Ly$\alpha$ emission to be detected with current facilities, as long as it is illuminated by an intense source of ionizing photons, such as a bright QSO. Indeed, numerous extended Ly$\alpha$ nebulae associated with AGN have been discovered around radio galaxies (e.g. Villar-Martín 2007 and references therein) and more recently around radio-quiet QSOs (e.g. Cantalupo et al. 2014; Hennawi et al. 2015).

Thanks to the Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010), Borisova et al. (2016) found 100% detection rate of giant Ly$\alpha$ nebulae around a sample of 17 radio-quiet and 2 radio-loud very luminous ($L\sim 10^{44}\;\textrm{erg}\;\textrm{s}^{-1}$) QSOs at $z\sim 3$. A similar result was also found by the MUSEUM survey (Arrigoni Battaia et al. 2019). Giant Ly$\alpha$ nebulae around luminous radio-quiet QSOs are the ideal laboratory where to investigate the CGM baryon fraction, for at least two reasons. First, bright QSOs are able to ionize the illuminated part of their own CGM almost completely, so that in this region contributions to Ly$\alpha$ emission from collisional excitation and scattering either are drastically reduced or their treatment is much simpler than in more ordinary environments. Second, both the typical QSO life time and the CGM light crossing time are much shorter than the CGM sound crossing time. This gives us some confidence that fluorescent nebulae around radio-quiet QSOs provide a genuine and unbiased snapshot of the state of CGM, before it might be significantly displaced by possible mechanical feedback from the current QSO activity event. 222The same does not necessarily apply to radio-loud QSOs, as relativistic jets might displace some portions of the CGM as rapidly as the ionizing photons make them visible through fluorescence. Note, nonetheless, that the two nebulae discovered by Borisova et al. (2016) around radio-loud QSOs have a surface brightness profile very similar to that of nebulae around radio-quiet objects.

The aim of this paper is to build upon the observations of Borisova et al. (2016) and provide an estimate of the CGM baryon fraction in relatively massive $(M_{\textrm{h}}\sim 10^{12}\;\textrm{M}_{\odot})$ haloes at $z\sim 3$. We find that a Ly$\alpha$ surface brightness as high as the observed one is difficult to explain with baryon fractions significantly smaller than the cosmological average. Lower values, typical of strong ejective feedback models, can be accounted for under specific conditions, which we discuss in detail and which can be tested by independent observations. Much of our analysis relies on the hypotheses that (i) Ly$\alpha$ emission is dominated by recombination radiation and (ii) that the cold phase and the hot phase are in pressure equilibrium with each other. We also discuss, however, the conditions and consequences of relaxing these assumptions.

In Section 2 we use the observations of Borisova et al. (2016) to infer the density profile and the total mass of the cold CGM gas in massive haloes at $z\sim 3$. In Section 3 we estimate the mass of the hot gas based on the assumption of pressure equilibrium. In Section 4 we discuss the total CGM mass, its dependence on model parameters and the conditions to reconcile it with the strong ejective feedback paradigm. In Section 5 we check the self-consistency of our model in terms of the ionization structure of the CGM. Section 6 addresses deviations from a pure-recombination scenario and in particular the impact of Ly$\alpha$ scattering. In Section 7 we discuss the possible origin and physical properties of the cold gas, including possible deviations from pressure equilibrium. A summary is given in Section 8.

In this paper we consistently adopt a flat cosmology with parameters $H_{0}=67.74\;\textrm{km}\;\textrm{s}^{-1}\;\textrm{Mpc}^{-1}$, $\Omega_{m}=0.3089$ and a cosmic baryon fraction $\Omega_{b}/\Omega_{m}=0.1573$ (Planck Collaboration 2016). Using cosmological parameters as in Borisova et al. (2016) has negligible impact on our results.

2 The cold CGM in MUSE QSO Ly$\alpha$ nebulae

2.1 Recombination radiation from the CGM

2.1.1 Ly$\alpha$ recombination surface brightness

If a galaxy hosts a sufficiently bright QSO, large portions of its CGM can be almost completely photo-ionized (see Section 5 below for more details). Recombinations in the medium would produce Ly$\alpha$ radiation, with a surface brightness given by the integral along the line of sight of the recombination Ly$\alpha$ emissivity

[TABLE]

where $f_{V}$ and $n$ are, respectively, the volume filling factor and the hydrogen number density of the Ly$\alpha$ emitting gas, while the coefficient $B$ is given by

[TABLE]

where $E_{\textrm{Ly}\alpha}$ is the energy of the Ly$\alpha$ line, $n_{\textrm{He}}/n_{\textrm{H}}=0.25Y/(1-Y)=1/12$ is the relative proportion of helium-to-hydrogen nuclei for a helium abundance by mass $Y=0.25$ and finally (posing $T_{4}\equiv T/10^{4}\;\textrm{K}$)

[TABLE]

is the effective Ly$\alpha$ recombination coefficient, as derived by fitting the calculations by Pengelly (1964) in the temperature range $5\times 10^{3}\;\textrm{K}<T<8\times 10^{4}\;\textrm{K}$. We adopt here case A coefficients. Note however that case A and B effective Ly$\alpha$ recombination coefficients differ by only $5\%$ (Pengelly 1964). They must not be confused with the total recombination coefficients, which differ in the two cases by a factor of $\sim 1.6$ (e.g. Osterbrock & Ferland 2006).

2.1.2 Filling factor

In a multiphase CGM, hot and diffuse gas at $T\sim 10^{6}-10^{7}\;\textrm{K}$ occupies the majority of the volume and compresses the relatively cold phase ($T\sim 1-5\times 10^{4}\;\textrm{K}$) to large densities and a small fraction of the volume. Given its dependence on density and temperature (equations 1 and 3), recombination radiation will be dominated by the cold and dense gas for any filling factor larger than $\sim 10^{-6}$. This justifies the approximation that Ly$\alpha$ emission is produced by the cold phase only and explains the importance of including the filling factor $f_{V}$ in equation (1).

We do not make assumptions here about the geometry or topology of the cold Ly$\alpha$ emitting regions (for instance, whether they are disconnected clouds or clumps, or a web of sheets and filaments), nor on their typical size, and we will loosely refer to them as cold gas structures. The possible nature and origin of the cold gas structure, as well as their typical sizes, are discussed further in Section 7.

If not stated otherwise, we interchangeably indicate with $n_{\textrm{cold}}$ or simply with $n$ the number density of hydrogen nuclei within the cold gas structures. This should not be confused with the volume-averaged density of cold gas, which is given by

[TABLE]

where the average is performed on volumes small with respect to the nebula but large with respect to the size of the emitting structures. Note that the term $f_{V}n^{2}$ in equation (1) could equivalently be replaced by $C\langle n\rangle^{2}$, where

[TABLE]

is the clumping factor of the emitting gas, as due to the fact that the latter only occupies a fraction of the total volume.

In this study, we neglect small-scale variations of the density within the structures themselves, i.e. their internal clumpiness, which would appear as an additional multiplicative factor $C_{\textrm{int}}$ on the right-hand side of equation (5) (though see a discussion in Section 7.3 below and see also Cantalupo et al. 2019). We do, however, allow for large-scale variations of the density $n$ with the physical distance $r$ from the centre of the nebula. This is both necessary – to explain the observed radially declining surface brightness profile (Section 2.2) – and expected, as the pressure of the confining hot gas is likely to decline with radius due to the gravitational pull of the halo (Section 3). The filling factor $f_{V}$ and the temperature of the cold gas are also in principle a function of radius, though in the following we assume them to be constant, for simplicity.

2.1.3 Main caveats

In the rest of this section, we apply equation (1) to the observed Ly$\alpha$ surface brightness profile of the MUSE QSO Giant Ly$\alpha$ nebulae, to infer the density and mass of the cold phase of the CGM around bright QSOs at $z\sim 3$. When doing so, one must be aware of a number of effects that could bias the result in one direction or the other. Dust absorption can attenuate the Ly$\alpha$ flux. We neglect dust, as we do not expect it to be present in the CGM or IGM in large amounts. Obviously, however, if the intrinsic Ly$\alpha$ surface brightness is larger than what we measure, our derived densities and masses should be interpreted as lower limits, strengthening our main finding of a large CGM baryon fraction. Similarly, if only a small fraction of the CGM is highly ionized (which may happen, for instance, because of a small aperture of the QSO ionization cone), then a significant portion of the existing cold CGM would not be visible in recombination radiation. This would also imply that our mass estimates are lower limits. On the other hand, mechanisms other than recombination, such as collisional excitation and scattering (or ‘photon pumping’) can contribute to the observed surface brightness. Scattering within the CGM can also mix photons belonging to different lines of sight, complicating the analysis. The most important of these points are discussed in greater detail in Sections 5 and 6. We argue there that they either do not affect our results, or can be taken into account by relatively simple modelling.

2.1.4 A note on ionization by stars

For some aspects of our problem – in particular, to estimate the temperature of the emitting gas, which determines the value of the recombination coefficient (see equation 3 and Section 4.3 below) – it is important to determine what the main source of ionizing photons is. In the close vicinity of a QSO, the QSO itself is obviously the most likely dominant source of ionizing photons. However, other sources such as young stars could in principle provide ionizing photons as well. For young stars to dominate the photoionization budget of the MUSE Giant Ly$\alpha$ nebulae, both the following two conditions must be met: (i) that the ionization cone of the QSO is sufficiently narrow that most of the CGM is not illuminated by the QSO and (ii) that there is sufficient star formation (either in the central galaxy or in the CGM itself) to keep the CGM completely ionized and power the observed Ly$\alpha$ luminosity. The first condition cannot be discarded, considering the large observational uncertainties (see e.g. Ichikawa et al. 2019 and references therein), although recent work suggests that the most luminous QSOs may have rather wide ionization cones (Ricci et al. 2017). The second necessary condition, on the other hand, is disfavoured based on energetic grounds and observed scaling relations. To power the observed median luminosity $L_{\textrm{Ly}\alpha}=10^{44}\;\textrm{erg}\;\textrm{s}^{-1}$, one would need a star formation rate of $100\;\textrm{M}_{\odot}\;\textrm{yr}^{-1}$ (e.g. Dijkstra & Westra 2010), assuming an implausibly large escape fraction of ionizing photons $f_{\textrm{esc}}=1$, or $4\times 10^{3}\;\textrm{M}_{\odot}\;\textrm{yr}^{-1}$ assuming a more realistic value $f_{\textrm{esc}}\simeq 2.5\%$ (e.g. Haardt & Madau 2012 and references therein). For comparison, the observed stellar-to-halo mass relation (e.g. Behroozi et al. 2013; Moster, Naab & White 2013) and main sequence of star-forming galaxies (e.g. Speagle et al. 2014 and references therein) imply that a halo with mass $M_{\textrm{halo}}\simeq 10^{12}\;\textrm{M}_{\odot}$ (close to the peak of star formation efficiency) at $z\sim 3$ would typically host a central galaxy with stellar mass $M_{\star}=10^{10.4}\;\textrm{M}_{\odot}$ and a star formation rate of $80\;\textrm{M}_{\odot}\;\textrm{yr}^{-1}$, i.e. 50 times smaller than required. Similarly, diffuse star formation within the CGM itself is an unlikely explanation for the ionization state of the MUSE nebulae, as it would require that the CGM forms stars at a rate equal to, or larger than, that in the central galaxy.

2.2 De-projected cold gas density profile

The MUSE QSO Giant Ly$\alpha$ nebulae (Borisova et al. 2016) have a median intrinsic (redshift-dimming corrected) surface brightness profile that is well described by a power-law function of the projected radius $R$:

[TABLE]

Note that the sample of Borisova et al. (2016) contains 17 radio-quiet and 2 radio-loud QSOs. We used here the entire sample; excluding the two radio-loud QSOs has negligible impact on our results. Note that the slope in equation (6) is different from the slope reported in Borisova et al. (2016). Our slope was obtained by fitting to the median of the individually dimming-corrected surface brightness profiles, improving on the rough first estimate in Borisova et al. (2016). Our results are not sensitive to small changes in the slope of the surface brightness profile (see also Section 6.2 below). We did not attempt to model individual nebulae, mainly because a measurement of the baryon fraction requires an estimate of the halo mass and this is known only in a statistical sense (see Section 4.2).

Assuming spherical symmetry (or as long as we focus on spherically averaged properties), we can deproject equation (6) by means of equation (1), to obtain the intrinsic radial profile of the density of the cold (Ly$\alpha$-emitting) gas:

[TABLE]

where $r$ is the spherical radius and $t=0.96$ is the exponent in relation (3). More in detail, the normalization and slope $(n_{0},\gamma)$ of the de-projected density profile (equation 7) are obtained from the normalization and slope $(\Sigma_{0},\beta)$ of the projected surface brightness profile (equation 6) by means of the relations

[TABLE]

where $R_{0}=10\;\textrm{kpc}$ and the dimensionless projection factor $\chi$ is defined (see also Pezzulli, Fraternali & Binney 2017) as

[TABLE]

The derivation above, which is equivalent to a classic Abel de-projection, formally relies on the assumption that the surface brightness profile (equation 6) and the density profile (equation 7) extend to infinity. As this is not realistic, we show in Figure 1 the effect, on the projected Ly$\alpha$ surface brightness profile, of introducing a truncation radius $r_{\textrm{max}}$ in the density profile (7). Predictions are shown for some values of $r_{\textrm{max}}$, equal to or larger than the virial radius of a $10^{12}\;\textrm{M}_{\odot}$ halo at $z=3$ ($r_{\textrm{vir}}\simeq 74\;\textrm{kpc}$, see equation 11 below). Truncation radii equal or larger than $r_{\textrm{vir}}$ are all consistent with the current data. Recent results from the MUSEUM survey (Arrigoni Battaia et al. 2019) possibly support a truncation radius not much larger than the virial radius. On the other hand (as already noted by Borisova et al. 2016), there is at least one example of a Giant Ly$\alpha$ nebula whose surface brightness profile is identical to those of the MUSE nebulae and extends to considerably larger physical radii (the Slug nebula; Cantalupo et al. 2014). We note, however, that azimuthal symmetry is a poor description of the Slug nebula, as well as of the outer regions of some of those MUSE nebulae that can be traced at radii larger than $r_{\textrm{vir}}$, so that the interpretation of profiles extending beyond the virial radius is not obvious. In the following, we will conservatively assume that the profile (7) holds within the virial radius, which is what we need to estimate the CGM baryon fraction there. In any case, we emphasize that our analysis is not meant to be applied far out in the IGM, where both strong asymmetries and lack of pressure confinement from a hot atmosphere are expected.

The dependencies in equation (7) on $f_{V}$ and $T_{4}$ are readily understood: if a smaller fraction of the volume is occupied by Ly$\alpha$ emitting gas, individual structures must be denser in order to explain the observed emission; at the same time, a larger temperature implies a smaller effective recombination coefficient (equation 3) and therefore again a larger density to give account of a fixed Ly$\alpha$ surface brightness. Note that the relation between inferred density and filling factor in equation (7) is equivalent to the familiar relation between inferred average density and clumping factor $\langle n\rangle\;\propto\;C^{-1/2}$, as it follows from equations (4) and (5).

2.3 The baryon fraction of the cold CGM

The density profile (7) can be integrated over volume to obtain the total mass of cold gas within the virial radius:

[TABLE]

where $m_{p}$ is the mass of the proton and $X=0.75$ is the cosmic hydrogen abundance by mass, while, as customary, we approximate the virial radius with the radius $r_{200}$ enclosing a 200-fold matter overdensity:

[TABLE]

where $M_{12}\equiv M_{\textrm{halo}}/10^{12}\;\textrm{M}_{\odot}$ is the halo mass in units of $10^{12}\;\textrm{M}_{\odot}$ and $z=3.2$ is the median redshift of the MUSE QSO Ly$\alpha$ nebulae. Finally, dividing $M_{\textrm{cold}}$ by the cosmological baryon fraction gives the fraction, among all the baryons nominally associated with the halo, that are stored in the cold phase of the CGM:

[TABLE]

where $f_{1}=7.3$, $t=0.96$ and $\gamma=1.25$.

As $f_{\textrm{CGM,cold}}$ cannot exceed unity, equation (12) implies that the filling factor of the cold gas $f_{V}$ must be smaller than 1 (or, equivalently, that the clumping factor must be larger than 1, cfr. equation 5), in agreement with previous findings (e.g. Cantalupo et al. 2014; Arrigoni Battaia et al. 2015). We now consider the option that the physical reason why the cold gas has $f_{V}<1$ is pressure confinement by a hot corona, which is likely present in the halo, although invisible in Ly$\alpha$ emission.

3 The mass of the hot CGM

Here we calculate the mass of the hot phase of the CGM that is needed to confine the cold phase to the densities derived in Section 2. Our treatment is based on the assumption of pressure equilibrium between the two phases and is not dissimilar, in its principles, to other works on the subject (e.g. Haiman & Rees 2001; Dijkstra & Loeb 2009; Vedantham & Phinney 2019). The possible presence and impact of deviations from pressure equilibrium are briefly discussed at the end of this section and in more detail in Section 7.

Pressure equilibrium dictates that the ratio $\varepsilon$ of the densities of the two phases (i.e. the density contrast) equals the inverse ratio of their temperatures:

[TABLE]

where $n_{\textrm{hot}}$ and $T_{\textrm{hot}}$ are the hydrogen number density and the temperature of the hot gas, respectively. Note that, in equation (13), we do not need to account for different ionized fractions (and therefore a different mean molecular weight) for the cold gas and the hot gas, as both phases are assumed here to be virtually completely ionized (although by entirely different processes: collisional ionization and photo-ionization for the hot and cold phases, respectively).

If one assumes that both $T_{\textrm{cold}}$ and $T_{\textrm{hot}}$ are independent on radius, then the density contrast $\varepsilon$ is also a constant, with the consequence that the inferred density profile of the hot gas follows a similar power law as the cold gas (equation 7), with a downward normalization shift, due to the fact that the hot gas is more diffuse:

[TABLE]

Interestingly, a power-law profile like that in equation (14) is indeed consistent with isothermal hydrostatic equilibrium in a dark matter halo gravitational potential, for a gas with close-to-virial temperature

[TABLE]

where $m_{p}$ is the mass of the proton, $k_{B}$ the Boltzmann constant and $\mu=0.593$ the mean molecular weight. This suggests that temperature gradients, if present, are not dramatic.

By combining equations (11), (13) and (15), the density contrast can be written

[TABLE]

with

[TABLE]

for $z=3.2$ and $\gamma=1.25$.

The total mass of the hot gas within the virial radius is readily obtained by integrating equation (14). Equivalently, the ratio of hot to cold gas mass is given by the ratio of the densities times the ratio of the volumes occupied by the two phases

[TABLE]

where we have used equation (16) and taken into account that $f_{V}\ll 1$ for every viable model (cf. Section 2.3).

The quantitative implications of equations (12) and (18) are presented in Section 4. Before doing so, we briefly discuss here whether any significantly large errors on the mass of the hot gas are expected to be introduced by our simplifying assumptions.

It is easy to verify that the simple inference above is not affected by the detailed shape of the gravitational potential. In a realistic potential (e.g. Navarro, Frenk & White 1996), the hydrostatic equilibrium density profile of an isothermal hot gas would not be exactly described by a power law, but rather require an additional radial-dependent factor equal to $(V_{\textrm{circ}}(r)/V_{200})^{-2}$, where $V_{\textrm{circ}}(r)$ is the circular speed at radius $r$ and $V_{200}=V_{\textrm{circ}}(r_{\textrm{vir}})$. In the situation of interest, this correction is remarkably small (less than 10% for $0.2<r/r_{\textrm{vir}}<1$), mostly due to the relatively low concentration of dark matter haloes at this redshift ($c\sim 3.5$; e.g. Dutton & Macciò 2014). The impact on the total mass is even smaller (less than 1%), mainly because most of the mass is at large radii, where the approximation $V_{\textrm{circ}}\simeq V_{200}$ is the most accurate. Coherent rotation, which must be present as the CGM necessarily has non-negligible angular momentum, also affects the density profile by providing centrifugal support against gravity, but also in this case the impact is mostly limited to the very innermost regions, which contain a very small fraction of the mass (Pezzulli, Fraternali & Binney 2017; Oppenheimer 2018).

More than the detailed shape of the potential or the presence of coherent rotation, deviations from the simple scenario above can arise due to large-scale turbulence, or other sources of non-thermal pressure. Crucially, however, such deviations would contribute both to enhance the compression of the cold phase and to provide additional support to the hot gas against gravity, with the result that the density contrast $\varepsilon$ (which is essentially estimated by taking the ratio of the two effects) would remain largely unaffected. A formal way to see this is to reformulate equation (13) in terms of a ratio of sound speeds squared $\varepsilon=c^{2}_{s,cold}/c^{2}_{s,hot}$. If large-scale turbulence is in place, the cold gas would be subject to (and compressed by) additional ram pressure, which would be, at least to the first order, encapsulated by the replacement $c^{2}_{\textrm{s,hot}}\rightarrow c^{2}_{\textrm{s,hot}}+v^{2}_{\textrm{turb, hot}}$. An analogous reformulation and replacement would be needed in equation (15) as well, as the very same turbulent motions would also provide support to the hot gas against gravity.333viz., the equation $c_{\textrm{s, hot}}^{2}=GM_{\textrm{halo}}/\gamma r_{\textrm{vir}}$ (equivalent to equation 15) should be replaced by $(c_{\textrm{s, hot}}^{2}+v_{\textrm{turb,hot}}^{2})=GM_{\textrm{halo}}/\gamma r_{\textrm{vir}}$. Note that only here we are using the ‘isothermal sound speed’ $c_{s}\equiv\sqrt{P/\rho}$ to simplify the notation and without loss of generality. It is then easily seen that the final estimated contrast $\varepsilon$ remains unaffected.444The same cannot be said for the temperature of the hot gas and therefore, for instance, for predictions of the X-ray luminosity of the CGM. Although we cannot claim this result to be exact, physical insight leads us to believe that additional corrections, if present, should be second-order and sub-dominant. Similar arguments should hold essentially for any large-scale source of non-thermal pressure. Considerations about other kinds of pressure imbalance are discussed in Section 7.

4 The total CGM baryon fraction

By combining equations (12) and (18), we finally obtain the total (cold plus hot) CGM baryon fraction within the virial radius:

[TABLE]

where we recall that $f_{1}=7.3$, $\varepsilon_{1}=3.0\times 10^{-3}$, $t=0.96$ and $\gamma=1.25$ (equations 3, 8, 12 and 17).

The total CGM baryon fraction thus depends on three parameters: $f_{V}$, $M_{\textrm{halo}}$ and $T_{\textrm{cold}}$. We discuss these dependencies below, one at a time.

4.1 Dependence on the filling factor

In Figure 2 we show how the total CGM baryon fraction, inferred from MUSE observations and from our model, depends on the filling factor of the cold gas $f_{V}$. We adopt here a halo mass $M_{\textrm{halo}}=10^{12.3}\;\textrm{M}_{\odot}$ and a temperature of the cold gas $T_{\textrm{cold}}=10^{4}\;\textrm{K}$, but we stress that these values are fixed here for illustrative purposes only. The effects of varying $M_{\textrm{halo}}$ and $T_{\textrm{cold}}$ are discussed in detail in Sections 4.2 and 4.3, respectively. In Figure 2, we also show the individual contributions of the cold gas and the hot gas to the total baryon fraction of the CGM.

As it is clear from the figure and from equation (12), the inferred masses of the cold gas and the hot gas have opposite dependencies on $f_{V}$. The estimated mass of the cold CGM increases with increasing $f_{V}$. This is a consequence of the dependence of the recombination emissivity on the density squared (see Section 2). On the other hand, the estimated mass of the hot gas increases with decreasing $f_{V}$. This is mostly because a larger amount of hot gas is needed to compress the cold gas to the larger densities that are associated with smaller filling factors (Section 3).

These opposite dependencies have two obvious consequences. First, there is at most a finite range of filling factors that are consistent with a physically meaningful baryon fraction (equal to or smaller than the cosmological value, shown as a horizontal black dashed line). Second, for any given choice of the other two parameters ($M_{\textrm{halo}}$ and $T_{\textrm{cold}}$), there is a minimum value of the total CGM baryon, below which no model based on recombination radiation and pressure balance would be able to explain the observed Ly$\alpha$ surface brightness of the MUSE nebulae. We will use this minimum in Section 4.4 to put a lower limit to the total CGM baryon fraction, as inferred from the MUSE observations within our assumptions. We emphasize that the minimum CGM baryon fraction is obtained when the two phases contribute equally to the total mass: if either of the two phases dominates in mass, the total baryon fraction will be larger than this lower limit.

4.2 Dependence on halo mass

Figure 3 shows the dependence of the inferred CGM baryon fraction on the assumed mass of the halo. Different lines correspond to the total (cold plus hot) CGM baryon fraction, for different values of the assumed halo mass and for one fixed value of the temperature of the cold gas $T_{\textrm{cold}}=10^{4}\;\textrm{K}$.

The inferred CGM baryon fraction decreases with increasing assumed halo mass. This is partly due to the obvious fact that increasing the halo mass also increases the amount of baryons theoretically associated with the halo, so that a lower fraction of those is needed to explain the observed Ly$\alpha$ surface brightness. Note however that the virial radius and the virial temperature are also changing accordingly, which explains why the dependence on halo mass is not trivially linear (equation 19).

With decreasing halo mass, it becomes increasingly difficult to accommodate a total CGM baryon fraction significantly smaller than the cosmological value, until eventually no physically acceptable solution is found. This makes the knowledge of the typical mass of QSO-hosting haloes at $z\sim 3$ crucial to draw conclusions from our results.

Probably the most reliable estimates of QSO-hosting halo masses are those based on QSO clustering, or on the QSO-galaxy cross-correlation function. Several studies of this kind have been attempted at a redshift $z\sim 3$ relevant to our application, finding a broad range of results, from $M_{\textrm{halo}}\sim 10^{12}\;\textrm{M}_{\odot}$ to $M_{\textrm{halo}}\sim 10^{13}\;\textrm{M}_{\odot}$ (see e.g. Shen et al. 2007 and references therein). Big uncertainties are to a large extent due to low-number statistics, inherent to studies of rare objects such as QSOs, but possibly also to differences in selection techniques and especially in the methods used to determine the redshift of the sources, which is crucial for clustering measurements. Taking these aspects into account, we chose to use, as a reference, the study by Eftekharzadeh et al. (2015). Their sample in fact is – by far – the one with the cleanest selection and the better statistics to our knowledge at the redshift of interest, with more than 20 000 QSOs with spectroscopic redshifts in the range $2.64<z<3.4$, as determined from the final sample of the Baryon Oscillation Spectroscopic Survey (BOSS; Dawson et al. 2013). These authors find a typical halo mass $M_{\textrm{halo}}=10^{12}\;\textrm{M}_{\odot}$, which is on the low side of the range of masses mentioned above. We refer to Eftekharzadeh et al. (2015) for a discussion of why larger values had been found by the previous work of Shen et al. (2007) on a smaller sample.

Despite the lack of new studies with sample sizes comparable to that of Eftekharzadeh et al. (2015), the debate on the mass of QSO-hosting haloes at $z\sim 3$ is still ongoing in the recent literature. For instance, Allevato et al. (2016) inferred a large mass $M_{\textrm{halo}}\sim 10^{12.9}\;\textrm{M}_{\odot}$ based on a sample of $\sim$ 400 X-ray selected AGN in a range $2.9<z<5.5$, $\sim 100$ of which with spectroscopic redshift. Also, Timlin et al. (2018) determined halo masses with a large uncertainty $M_{\textrm{halo}}=1.70-9.83\times 10^{12}\;\textrm{h}^{-1}\;\textrm{M}_{\odot}$ based on a sample of 1378 objects, identified as QSO candidates and assigned a photometric redshift in the range $2.9<z<5.1$ by a machine-learning algorithm. The results by Allevato et al. (2016) and Timlin et al. (2018) emphasize the importance of further investigations to attenuate uncertainties in this field. However, given the much smaller samples of these studies with respect to Eftekharzadeh et al. (2015), as well as the less robust selection techniques (especially concerning the determination of redshift), we think that the BOSS result should, for the time being, be considered as the reference.

Even if the BOSS result is (at least formally) the best measurement available to date, it is still not obvious that it can be safely applied to our study. The sample of Eftekharzadeh et al. (2015) in the relevant redshift range contains QSOs as luminous as $\log(L_{\textrm{bol}}/\textrm{erg}\;\textrm{s}^{-1})=47.72$. This overlaps with the luminosity range of the QSOs in the sample of Borisova et al. (2016) ($47.29<\log(L_{\textrm{bol}}/\textrm{erg}\;\textrm{s}^{-1})<48.14$, using equation 1 in Shen et al. 2009), but the latter also contains objects up to 2.6 times more luminous than the maximum luminosity probed by Eftekharzadeh et al. (2015). If, as some theoretical models predict, the halo mass correlates strongly with QSO luminosities (see e.g. Hopkins et al. 2007; Conroy & White 2013), it is possible that the QSOs of Borisova et al. (2016) live in haloes with a larger mass than the QSOs in the BOSS sample. On the other hand, a number of studies have recently shown that statistical measurements of QSO host halo masses correlate surprisingly weakly with QSO luminosities (e.g. Shen et al. 2009; Padmanabhan et al. 2009; Eftekharzadeh et al. 2015; Uchiyama et al. 2018; He et al. 2018), so that this option is not very clearly supported by the current data.

Furthermore, we notice that the MUSEUM survey (Arrigoni Battaia et al. 2019) has investigated Giant Ly$\alpha$ nebulae around a sample of QSOs less luminous than those of Borisova et al. (2016), reaching a maximum luminosity of $\log(L_{\textrm{bol}}/\textrm{erg}\;\textrm{s}^{-1})=47.87$. This is very similar to the maximum luminosity probed by Eftekharzadeh et al. (2015), possibly making comparisons less prone to potential bias. The Giant Ly$\alpha$ nebulae in the MUSEUM sample have a surface brightness profile that is undistinguishable from those of Borisova et al. (2016) in the radial range of interest for this study, possibly suggesting that a bias in QSO luminosity is not the main explanation of the large surface brightnesses observed in all these objects, although studying Ly$\alpha$ nebulae around QSOs with luminosities even lower than those in the MUSEUM survey is needed to achieve a more complete picture on this matter.

Based on the discussion above, in the following, we will sometimes use as fiducial the halo mass $M_{\textrm{halo}}=10^{12}\;\textrm{M}_{\odot}$ found by Eftekharzadeh et al. (2015). However, we will also show and discuss results obtained in a broad range of masses $M_{\textrm{halo}}=10^{12.3\pm 0.5}\;\textrm{M}_{\odot}$ (Trainor & Steidel 2012), which we chose as representative of the range of uncertainty in the value of the halo mass.

If a low value $M_{\textrm{halo}}\sim 10^{12}\;\textrm{M}_{\odot}$ is confirmed by future studies, it would imply that it is not easy to explain the surface brightness of MUSE Giant Ly$\alpha$ nebulae with a baryon fraction significantly smaller than cosmological, at least within the assumptions discussed so far.

On the other hand, a significantly larger value for the halo mass would not necessarily imply that the observed Ly$\alpha$ surface brightness is trivial to reproduce. Most models in fact agree that at sufficiently large halo masses most of the gas is expected to be in the hot phase. This would violate the first condition to maximize the predicted Ly$\alpha$ surface brightness and specifically that the two phases contribute equally to the total mass of the CGM (Section 4.1). Equivalently, the relevant curve in Figure 3 would need to be read at some (model-dependent) displacement from the minimum. The quantitative details of the transition can be rather different in different models, depending on both numerical aspects (e.g. Nelson et al. 2013) and on the adopted feedback recipe (e.g. Correa et al. 2018).

4.3 Dependence on the temperature of the cold gas

In Figure 4 we show the dependence of the inferred total CGM baryon fraction on the assumed temperature of the cold photo-ionized gas $T_{\textrm{cold}}$, for a fixed halo mass $M_{\textrm{halo}}=10^{12.3}\;\textrm{M}_{\odot}$. The inferred CGM baryon fraction increases with increasing the assumed $T_{\textrm{cold}}$. The reason is two-fold. First, as the recombination coefficient is a declining function of temperature (equation 3), a larger $T_{\textrm{cold}}$ requires a larger amount of cold gas to explain the same observed Ly$\alpha$ emission. Second, a larger $T_{\textrm{cold}}$ also implies a larger pressure of the cold gas and therefore a larger amount of hot gas needed to keep it confined. The magnitude of the effect is significant and it is therefore important to make physically motivated assumptions concerning the value of $T_{\textrm{cold}}$.

In an almost fully photo-ionized gas, the temperature is mainly determined by the equilibrium between photo-electric heating and recombination cooling and therefore, ultimately, by the shape of the ionizing spectrum (e.g. Osterbrock & Ferland 2006). In our case, the ionizing spectrum is completely dominated by the QSO (we will see in Section 5 that the UVB can be neglected by orders of magnitude). In a pure hydrogen nebula, the resultant equilibrium temperature would be

[TABLE]

where $h$ is the Planck’s constant, $\nu_{\textrm{LL}}$ is the Lyman-limit frequency and $p$ is the power-law slope of the hydrogen ionizing continuum, for which we take $p=1.7$ as a fiducial value (Lusso et al. 2015, see also Section 5 below).

In a QSO ionizing field, however, the photo-electric effect due to the second ionization of helium is also a very important heating channel. Taking this into account and assuming that a single power-law spectrum can be extrapolated to $E>4\;\textrm{Ry}$, the equilibrium temperature increases by a factor $(1+4A(T))/(1+A(T))$, where

[TABLE]

By using the temperature-dependent recombination coefficients as in Meiksin (2009), we then find an equilibrium temperature

[TABLE]

This is large enough that no physically acceptable solution can be found for a halo mass $M_{\textrm{halo}}=10^{12.3}\;\textrm{M}_{\odot}$, as it is clear for instance from Figure 4. Note that the problem would only be exacerbated for $M_{\textrm{halo}}=10^{12}\;\textrm{M}_{\odot}$, as suggested by clustering measurements (Section 4.2).

Deviations from the simple estimate above would mostly come from (i) line cooling from residual non-completely ionized hydrogen and helium and (ii) traces of metals, which would constitute additional line cooling channels, but also additional targets to photo-electric heating. Other potentially relevant effects are those whereby the energy of a photon can be distributed to more than one electron and primarily (i) ionization of hydrogen by high-energy HeII recombination photons and (ii) secondary ionization (of hydrogen or other species) by high-energy electrons produced in a previous ionization event (Osterbrock & Ferland 2006). To assess the role of all these effects, we performed more detailed calculations of the equilibrium temperature using the photo-ionization code Cloudy (version 13.04, last described by Ferland et al. 2013). We modelled the hydrogen and helium ionizing spectrum as a single power law with slope $p=1.7$ . We also included a hardening at X-ray frequencies (constant $\nu L_{\nu}$ in the energy range $2\;\textrm{keV}<E<100\;\textrm{keV}$) and a cut-off for $E>100\;\textrm{keV}$, though we found the details of the X-ray spectrum to be unimportant. We adopted a specific luminosity at the Lyman limit $\tilde{L}_{\nu,\textrm{LL}}=5.2\times 10^{31}\;\textrm{erg}\;\textrm{s}^{-1}\;\textrm{Hz}^{-1}$ (see Section 5.1 below) and a distance from the source $d=50\;\textrm{kpc}$. Results are shown in Figure 5 for different values of gas density and metallicity.

We see that the results are in broad agreement with the simple analytical estimate above (equation 22) especially for intermediate densities and metallicities. As also expected, the equilibrium temperature decreases with increasing density, due to collisionally excited line cooling, the effect being larger for larger metallicity. Note that at low densities the effect of non-primordial metallicity is reversed, as photo-electric heating on highly ionized metal species (as well as a small contribution of Compton heating) dominates over additional line cooling, so that the equilibrium temperature is even larger than the analytic estimate.

An equilibrium temperature equal or smaller than $T\sim 2\times 10^{4}\;\textrm{K}$, as desirable from inspection of Figure 4, can be reached only for a rather high density $n=10\;\textrm{cm}^{-3}$ and solar metallicity. Note that at a distance $d=50\;\textrm{kpc}$ and a temperature $T=2\times 10^{4}\;$, a density $n=10\;\textrm{cm}^{-3}$ corresponds to a filling factor $f_{V}=7.3\times 10^{-6}$ (equation 7), which would give an absurdly large baryon fraction even at low temperature (Figure 4). Furthermore, it is not obvious that a metallicity as high as solar can be easily reached in the CGM at $z=3$. Doing the calculation at smaller distances $d<50\;\textrm{kpc}$ would not help in obtaining lower temperatures, as the ionization parameter increases with decreasing radius for our density slope $\gamma=1.25<2$.

As a sanity check, we also repeated our calculation including only the (QSO-dominated) ionizing background at $z=3.2$, a metallicity $Z=0.1\;Z_{\odot}$ and a density equal to 10 times the critical density at $z=3.2$, in order to mimic the conditions in the IGM far away from an individual QSO. We found in this case an equilibrium temperature $T_{\textrm{eq}}=3\times 10^{4}\;\textrm{K}$, in agreement with models of HeII re-ionization (e.g. McQuinn et al. 2009; Upton Sanderbeck, D’Aloisio & McQuinn 2016). The relatively low equilibrium temperature in this case is due to the relatively low ionization parameter at large distances from a single bright QSO. This emphasizes that temperatures as high as those shown in Figure 5 are expected, for a gas photo-ionized by a QSO, only in the very close vicinity to the source, such as in the CGM of the host galaxy, which is the main focus of this work.

Another possibility to achieve lower equilibrium temperatures at the small distances of interest could be to change the shape of the ionizing spectrum, especially at energies $E>4\;\textrm{Ry}$, which contribute most of the heating and on which observational constraints are very scant. We therefore also considered a modified spectrum, identical to the one above but for being double as steep at $E>4\;\textrm{Ry}$ than in the range $1\;\textrm{Ry}<E<4\;\textrm{Ry}$. However, both analytic estimate and Cloudy modelling showed in this case only a modest $\sim 25\%$ decrease in $T_{\textrm{eq}}$. Virtually, only an exponential cut-off in the EUV spectrum would be effective in beating the temperature down as required. We note that this is in fact the theoretical expectation for the emission spectrum of a classical accretion disc around a supermassive black hole (see e.g. Shakura & Sunyaev 1973; Laor & Davis 2011). The predictions of this kind of model are illustrated in Figure 6, where we assumed the spectrum of a classical (geometrically thin, optically thick) accretion disc (Shakura & Sunyaev 1973) with Eddington accretion around a supermassive black hole with mass $M_{BH}=10^{8}\;\textrm{M}_{\odot}$, ignoring in particular any comptonization effect. As expected, this model easily predicts an equilibrium temperature $T_{\textrm{eq}}\simeq 2\times 10^{4}\;\textrm{K}$ for a metallicity $Z=0.1\;Z_{\odot}$ and a broad range of densities and even lower temperatures ($T_{\textrm{eq}}\simeq 10^{4}\;\textrm{K}$) for solar metallicity. This kind of solution is formally appealing and it emphasizes the crucial role of assumptions on the shape of the ionizing spectrum of QSOs. Nonetheless, we refrain from considering this option as our fiducial case, as it would raise significant issues with a large body of independent observations including (i) widespread evidence of comptonization in AGN (e.g. Done et al. 2012 and references therein) and (ii) indirect evidence based on HeII re-ionization (e.g. Worseck et al. 2016 and references therein).

Motivated by the discussion above, in the following we will sometimes adopt as fiducial the value $T_{\textrm{cold}}=5.6\times 10^{4}\;\textrm{K}$ (cf. equation 22 and Figure 5). We will nonetheless also present our results for a large range of temperatures to allow for more generality.

4.4 The minimum CGM baryon fraction

A useful way to summarize the results of the previous sections is to compute, for each value of the parameters $M_{\textrm{halo}}$ and $T_{\textrm{cold}}$, the minimum total CGM baryon fraction. As it is readily found from equation (19), this is equal to

[TABLE]

and is achieved for the filling factor $f_{V,\textrm{eq}}$ such that cold gas and hot gas contribute equally to the total mass:

[TABLE]

Equation (23) gives the lower limit to the total CGM baryon fraction needed to explain the observations within our assumptions. It is shown in Figure 7, as a function of the two parameters $M_{\textrm{halo}}$ and $T_{\textrm{cold}}$. The solid black line indicates the cosmological baryon fraction. The blank region at low masses and large temperatures corresponds to unphysical solutions (baryon fraction larger than cosmological), while progressively smaller baryon fractions are inferred at large halo masses and low temperatures.

The couple of fiducial values $M_{\textrm{halo}}\simeq 10^{12}\;\textrm{M}_{\odot}$ and $T_{\textrm{cold}}\simeq 5.6\times 10^{4}\;\textrm{K}$ (determined in Sections 4.2 and 4.3, respectively) fall in the unphysical region. The simplest explanation for this fact is that at least one of these values is incorrect. We do not have strong reasons to prefer a revision of either $M_{\textrm{halo}}$ or $T_{\textrm{cold}}$. We stress, however, that both parameters can be constrained by independent observations and in particular by further studies of clustering and spectral properties of luminous QSOs (Sections 4.2 and 4.3). Another option is that both our fiducial values are correct, but that some of our other assumptions (recombination-dominated emission and pressure equilibrium) is wrong. These possibilities are discussed in Sections 6 and 7, but we anticipate that they are not trivially achieved either.

Figure 7 can also be compared to cosmological hydrodynamical simulations or to semi-analytic models, as long as these include predictions for the CGM baryon fraction in massive haloes at high redshift.

Cosmological hydrodynamical simulations with strong ejective feedback (see Section 1) typically predict a CGM baryon fraction equal to $\sim 30-40\%$ of the cosmic average at the relevant halo mass and redshift (e.g. Davé 2009; Liang et al. 2016; Kulier et al. 2019; Mitchell, private communication). Unfortunately, not all simulation projects have published suitable predictions, so we cannot make an exhaustive comparison here. On the other hand, we hope that our results can encourage more of these predictions in the future. Based on the references above and on Figure 7, we expect that simulations with strong ejective feedback should be able to reproduce the observed Ly$\alpha$ surface brightness, under the assumption of fluorescent emission (recombination), only if the analysis is performed on rather massive haloes ($M_{\textrm{halo}}\gtrsim 10^{12.5}\;\textrm{M}_{\odot}$) and the cold CGM temperature is assumed to be sufficiently low ($T_{\textrm{cold}}\lesssim 2\times 10^{4}\;\textrm{K}$), despite photo-heating by the central QSO. Furthermore, it is required that the CGM mass is split almost equally among the cold phase and the hot phase, so that the former is both massive and dense enough to give rise to sufficiently bright recombination radiation. In Section 4.5 we discuss in some detail the only experiment published so far that can be directly compared to our expectations.

Among semi-analytic models with strong ejective feedback, the one of Henriques et al. (2019) predicts that a typical halo with mass $M_{\textrm{halo}}=10^{12}\;\textrm{M}_{\odot}$ at $z=3$ would have a median of 53% of its baryons in the CGM. This is slightly larger than what found by the simulations mentioned above and therefore potentially closer to what is needed to explain the observations. Note however that the CGM in the model of Henriques et al. (2019) is all hot and therefore invisible in Ly$\alpha$, unless a large fraction of it is somehow assumed to be converted into cold CGM.

4.5 Comparison with Gronke & Bird 2017

A special mention is needed to the study of Gronke & Bird (2017). These authors analysed massive haloes at $z=2$ in the Illustris simulation (Vogelsberger et al. 2014) and calculated the expected recombination Ly$\alpha$ surface brightness, under the assumption that the CGM is highly ionized. They found agreement with the observed surface brightness profile of the Slug nebula and of MUSE Giant Ly$\alpha$ nebulae (note that the Borisova et al. 2016 profile in their Figure 9 should be shifted upwards by 0.58 dex due to differential redshift-dimming between the observed $z=3.2$ and the simulated $z=2$). As they assumed recombination radiation from a highly ionized medium, their experiment is directly comparable to our analytic calculations. As discussed in Section 4.4, we expect simulations with strong ejective feedback to be able to reproduce a Ly$\alpha$ surface brightness as high as observed, as long as the assumed halo mass and gas temperature are sufficiently large and small, respectively. Indeed, Gronke & Bird (2017) adopted a range of halo masses centred on $M_{\textrm{halo}}=10^{12.5}\;\textrm{M}_{\odot}$, which is on the high side of those discussed in Section 4.2. Furthermore, their post-processing analysis neglected photo-heating from the QSO, with the result that they most likely underestimated the temperature of the photo-ionized CGM. Note that Illustris does include a recipe for radiative feedback from QSOs (Vogelsberger et al. 2013). Gronke & Bird (2017), however, did not select their haloes on the requirement that the central super massive black hole be in an active phase and therefore a QSO would be found therein with a probability that can be approximated by the QSO duty cycle $f_{\textrm{duty}}\ll 1$ (e.g. Haiman, Ciotti & Ostriker 2004; Eftekharzadeh et al. 2015). Note also that any imprint of past QSO activity on photo-ionized gas temperatures would die out on roughly one recombination time-scale, or less than $1\;\textrm{Myr}$ in our case.

We therefore think that the formal success found by this work is not surprising and actually in line with our analytic expectations, but it may be driven by a particularly favourable choice of the parameters. In particular, our analysis shows that neglecting in post-processing photo-heating of the cold CGM by the central QSO can significantly bias the recombination rate and therefore the predicted Ly$\alpha$ surface brightness in the direction needed to make the comparison with observations easier.

Gronke & Bird (2017) also considered the option that most of the CGM is ionized by star formation rather than by the central QSO, in which case the temperature of the photo-ionized gas can be lower than our estimates. This option, however, is strongly disfavoured for the MUSE Giant Ly$\alpha$ nebulae, as we have argued already in Section 2.1.4.

5 Ionization structure

Here we verify the self-consistency of our model. We start from the key requirement that the CGM can be almost completely photo-ionized by the central QSO (Section 5.1), before discussing the ionization of helium (Section 5.2) and the possible contribution of collisional excitation to the Ly$\alpha$ emission (Section 5.3). We obviously refer here only to the fraction of the solid angle that is illuminated by the QSO. We recall that considering any substantial QSO obscuration (by a dusty torus or any other intervening material) would only go in the direction of reinforcing our conclusions (see Section 2.1.3).

5.1 Hydrogen

The QSOs in the sample of Borisova et al. (2016) have a median specific luminosity at the Lyman limit frequency $\nu_{\textrm{LL}}$ of $\tilde{L}_{\nu,\textrm{LL}}=5.2\times 10^{31}\;\textrm{erg}\;\textrm{s}^{-1}\;\textrm{Hz}^{-1}$, which we estimated following Arrigoni Battaia et al. (2015) and using the absolute i-band magnitudes of the QSOs in the Borisova et al. (2016) sample (see their Table 2). The corresponding hydrogen ionization rate at distance $r$ from the QSO is

[TABLE]

where $\sigma^{\textrm{ion}}(\nu)\simeq\sigma^{\textrm{ion}}_{\textrm{LL}}\;(\nu/\nu_{\textrm{LL}})^{-3}$, with $\sigma^{\textrm{ion}}_{\textrm{LL}}=6.3\times 10^{-18}\;\textrm{cm}^{2}$, is the hydrogen photo-ionization cross section and $L_{\nu}\simeq L_{\nu,\textrm{LL}}(\nu/\nu_{\textrm{LL}})^{-p}$ is a power-law approximation of the ionizing continuum of the QSO, with a typical slope $p=1.7$ (e.g. Lusso et al. 2015), while the factor

[TABLE]

incorporates possible deviations of the normalization and slope of the QSO spectrum from our fiducial values. Note that adopting a more precise double power-law approximation to the frequency dependence of the hydrogen ionization cross section (see e.g. Meiksin 2009) is equivalent to replacing in the expressions above $\frac{1}{p+3}\rightarrow\frac{p+s+\beta}{(p+s)(p+s+1)}$, with $s=2.99$ and $\beta=1.34$. This gives a correction of less than $7\%$ for any $p>1$ and is negligible to our purposes.

The photo-ionization rate given in equation (25) should be compared to the (case A) total hydrogen recombination rate (per proton) for a fully ionized gas:

[TABLE]

where again we adopted a power-law approximation to a more complex temperature dependence (Meiksin 2009) over the temperature range of interest and we accounted for electrons from doubly ionized helium. It is readily seen that $\Gamma_{\textrm{rec}}\ll\Gamma_{\textrm{ion}}$ for any plausible density for CGM conditions (the equality would require $n\simeq 2\times 10^{4}\;\textrm{cm}^{-3}$ at $r=100\;\textrm{kpc}$ and even larger densities at smaller radii). Note also that $\Gamma_{\textrm{ion}}$ exceeds the photo-ionization rate due to the cosmic UV background (e.g. Haardt & Madau 2012) by several orders of magnitude at any distance plausibly associated with the CGM.

Again following Meiksin (2009), the equilibrium hydrogen neutral fraction under these conditions is achieved on very short time-scale $t_{\textrm{ion}}=\Gamma_{\textrm{ion}}^{-1}$ (everywhere smaller than a few years in our case) and very well approximated by $x_{HI}=\Gamma_{\textrm{rec}}/\Gamma_{\textrm{ion}}$, which gives

[TABLE]

Combining equations (7) and (28), we obtain the neutral hydrogen number density $n_{HI}=x_{HI}n$, which reads

[TABLE]

and therefore the volume-averaged neutral hydrogen number density

[TABLE]

which is independent on the filling factor. This is not surprising, if one considers that both recombination radiation (on which our inferred density profile of equation 7 is based) and the neutral fraction are determined by the same physical process (obviously, recombination). It is similarly easy to show that the inferred volume-averaged neutral hydrogen number density $\langle n_{HI}\rangle$ is actually independent also on assumptions on the ionization state of helium (which enters in both equations 2 and 27) and on possible internal clumpiness $C_{\textrm{int}}$ within the cold gas structures (cf. Section 2.1.2), as all these factors equally contribute to the recombination emissivity and the recombination rate. The small dependence on temperature is due to a non-constant efficiency of Ly$\alpha$ photons production per recombination. The main dependence is, of course, on the intensity and spectrum of the ionizing source, as parametrized here by the factor $\xi$ (equation 26).

If radially integrated, equation (29) gives the neutral hydrogen column density between the QSO and a point in the nebula at distance $r$ from the centre:

[TABLE]

which is everywhere much smaller than the Lyman limit column density $N_{H,LL}=\sigma_{\textrm{ion}}(\nu_{LL})^{-1}=10^{17.2}\;\textrm{cm}^{-2}$, for any plausible value of $\xi$ and $T_{4}$. This confirms that our main assumption that the entire nebula is almost completely ionized and transparent to hydrogen ionizing radiation is self-consistent. The implications for Ly$\alpha$ scattering, on the other hand, are more complicated and are addressed in detail in Section 6.

5.2 Helium

Similarly, the column density of HeII is

[TABLE]

with

[TABLE]

where we used again atomic coefficients as in Meiksin (2009) and we further assumed that the filling factor $f_{V}>10^{-6}$, while $p_{1}$ and $p_{2}$ are the slopes of the QSO ionizing continuum in the ranges $1\;\textrm{Ry}<E<4\;\textrm{Ry}$ and $E>4\;\textrm{Ry}$, respectively. If for simplicity we assume $p_{2}=p_{1}=1.7$ (i.e. if the spectrum by Lusso et al. 2015 can be extrapolated above 4 Ry), we see that the column density (32) is still safely below the limit for self-shielding $N_{HeII,\textrm{crit}}=6.3\times 10^{17}\;\textrm{cm}^{-2}$, implying that the entire nebula should be transparent to the helium-ionizing radiation as well.

Note, however, that Borisova et al. (2016) estimated very low HeII/Ly$\alpha$ line ratios, inconsistent with recombination radiation from fully ionized hydrogen and helium. Though the estimate above suffers from several approximations and extrapolations (most notably, the shape of the ionizing spectrum and the extrapolation of the modelled density and filling factor to radii smaller than 10 kpc), we found that it is not easy to accommodate HeII column density that are high enough to guarantee self-shielding. If this is true, and if the very low limits on the line ratios are confirmed by deeper observations, this could suggest that either some of the Ly$\alpha$ emission is contributed by processes other than pure recombination (see Sections 5.3 and 6 below; see also Humphrey et al. 2019), or that the majority of the emission comes from regions that are dense enough that the helium is only partially ionized despite being fully illuminated (e.g. Arrigoni Battaia et al. 2015). This in our case would imply a very low filling factor $f_{V}<10^{-6}$, or additional internal clumpiness within the cold gas structures (see Section 7.3 below; see also Cantalupo et al. 2019). It is of course also possible that the explanation is a combination of those mentioned above and it may be different for different nebulae. We plan to explore these aspects in more detail in future work.

5.3 Collisional excitation

We finally verify that the contribution of collisional excitation to the Ly$\alpha$ emission can be consistently neglected in our model.

The local Ly$\alpha$ emissivity due to collisional excitation and recombination are, respectively:

[TABLE]

and

[TABLE]

(we omitted ‘Ly$\alpha$’ super-scripts or sub-scripts j’s and $\alpha$’s for clarity) and their ratio is

[TABLE]

where we take $x_{H}I$ from equation (28) and the ratio of recombination to collisional excitation coefficients as in Cantalupo, Porciani & Lilly (2008). The ratio $\alpha^{\textrm{eff}}_{\textrm{coll}}/\alpha^{\textrm{eff}}_{\textrm{rec}}$ increases rapidly with temperature (much more rapidly than the decrease with temperature of $x_{HI}$, see equation 28). To obtain a conservative upper limit to the contribution of collisional excitation, we therefore consider a fairly high temperature $T=5\times 10^{4}\;\textrm{K}$ and also a fairly large distance from the centre $r=100\;\textrm{kpc}$ (to maximize the neutral fraction). Note that choosing even higher temperatures would not help, as collisional ionization would start to dominate over photo-ionization driving an additional drop of the neutral fraction and therefore of the collisional excitation rate. In these conservative conditions, we obtain as an upper limit:

[TABLE]

so that we can safely ignore collisional excitation as a powering mechanism for filling factors larger than $10^{-4}$ (a limit that goes further down for smaller distances and temperatures).

6 The impact of scattering

Here we discuss the possible impact on our results of Ly$\alpha$ scattering. We distinguish between the addition of Ly$\alpha$ photons, which may be produced by the central QSO and scattered into the line of sight (Section 6.1), from the redistribution of Ly$\alpha$ photons due to scattering within the CGM (Sections 6.2 and 6.3).

6.1 Scattering from the broad-line region

Even if the CGM is almost completely ionized, some of the very rare neutral hydrogen atoms can still intercept some of the Ly$\alpha$ photons coming from the broad-line region (BLR) of the central QSO and scatter them in to the line of sight. These photons would contribute to the observed Ly$\alpha$ surface brightness of the nebulae, in addition to the recombination radiation we have focused on so far.

To quantify this effect, we must convert the neutral hydrogen column density as seen by the QSO (equation 31) into a line-centre optical depth to Ly$\alpha$ scattering (following e.g. Meiksin 2009). We find

[TABLE]

where $\sigma_{300}$ is the 1D velocity dispersion of the cold CGM, as seen by the QSO, in units of 300 $\textrm{km}\;\textrm{s}^{-1}$.

We have therefore found a model-dependent optical depth of order unity. This is the most difficult regime for analytic treatment, as both the single scattering and the diffusion limit are prone to be inadequate. On the other hand, a brute-force numerical approach would not give a clear-cut solution either, as the limiting factor here is not the calculation itself, but rather that the result is bound to be strongly dependent on the assumptions on the intrinsic kinematics of the cold gas, which is largely unknown.

The reference value for $\sigma$ that we used in equation (38) corresponds to the typical width of the Ly$\alpha$ line as observed in emission (Borisova et al. 2016). Note however that this can be a biased estimate of the true value of $\sigma$ for at least two reasons. First, if radiative transfer effects were indeed important (which is what we would like to establish here), the observed line width would likely be larger than the intrinsic one, with the result that the fiducial $\sigma$ above would be overestimated and the optical depth underestimated. Second, the velocity dispersion as seen by the QSO will be different from that along the line of sight unless the kinematics of the cold CGM is isotropic. For instance, if the cold CGM were dominated by coherent radial flows (such as inflows or outflows) the dispersion seen by the QSO (as determined by the acceleration or deceleration of the cold gas flows) would again be smaller than that along the line of sight (which would instead be dominated by projection effects). Both effects would go in the direction of favouring optical depths larger than unity. As we show below, a large optical depth would imply a small contribution of BLR photons to the observed emission and therefore that the model considered so far is self-consistent in this respect (although not necessarily unique).

Our ignorance on the intrinsic kinematics of the CGM prevents us from calculating the true contribution of scattered BLR photons to the observed emission. We can none the less easily put a very conservative upper limit to this contribution by calculating the maximum possible number of BLR Ly$\alpha$ scattered photons per fluorescent Ly$\alpha$ photon produced in a recombination event, based on fundamental atomic physics, as described below (we have also verified the validity of our reasoning with Cloudy simulations).

Recombination Ly$\alpha$ photons are produced, with some probability $\eta$, every time an electron recombines with a proton to form an atom of neutral hydrogen. Statistically, the electron will be freed again (and therefore become again eligible to recombination and Ly$\alpha$ emission) in a time equal to the inverse of the photo-ionization rate $t_{\textrm{ion}}=\Gamma_{\textrm{ion}}^{-1}$. During this time, the neutral atom is available to undergo a number of scattering events, statistically equal to $\Gamma_{\textrm{scatt}}t_{\textrm{ion}}=\Gamma_{\textrm{scatt}}/\Gamma_{\textrm{ion}}$, where $\Gamma_{\textrm{scatt}}$ is the scattering probability per unit time

[TABLE]

where $\bar{\sigma}=\int_{0}^{+\infty}\sigma_{\textrm{scatt}}(\nu)d\nu=1.1\times 10^{-2}\;\textrm{cm}^{2}\;\textrm{Hz}$ is the frequency-integrated cross-section to Ly$\alpha$ scattering (e.g. Rybicki & Lightman 1979; Meiksin 2009) and

[TABLE]

is the specific flux density of Ly$\alpha$ photons propagating directly from the BLR to distance $r$ from the QSO. The cumulative Ly$\alpha$ optical depth $\tau_{\textrm{Ly}\alpha,\textrm{QSO}}(r)$ describes attenuation by scattering at smaller radii.555Obviously attenuation here is not synonym to absorption. The fate of intercepted BLR photons after first scattering is addressed in Section 6.2. Note that we do not need to introduce an attenuation factor for the Lyman limit photons, as the optical depth to photo-ionization is negligible irrespective of assumptions on kinematics, at the column densities of interest (Section 5.1).

The local ratio of the scattering to recombination Ly$\alpha$ emissivity is given by

[TABLE]

where the case A Ly$\alpha$ production efficiency per recombination is well approximated by

[TABLE]

It is clear from equations (40) and (41) that a very conservative upper limit for the ratio of scattering to recombination emissivity is found in the Ly$\alpha$ optically thin limit $\tau_{\textrm{Ly}\alpha,\textrm{QSO}}\ll 1$. Furthermore, in this limit, the ratio $j_{\textrm{scatt}}/j_{\textrm{rec}}$ is just a constant, determined by the intrinsic spectrum of the QSO. This is not surprising if one considers that, in an almost completely ionized gas, one single physical mechanism (recombination) is responsible both for Ly$\alpha$ fluorescence and for the creation of scattering targets (neutral atoms) and that in both cases the energy source is the same (the QSO). Assuming the typical QSO spectrum from Lusso et al. (2015) and a UV slope $p=1.7$ we then find

[TABLE]

We stress that the value above is a very strict upper limit on the true ratio of the scattering-to-recombination emissivity, as it assumes that the whole CGM around the QSO sees the Ly$\alpha$ emission from the BLR entirely and unattenuated. Having clarified this, in such a limit one would have that the intrinsic recombination radiation surface brightness is $(1+\mathcal{R}_{\textrm{max}}(T))$ smaller than the observed one. As a consequence, all the densities and the masses discussed in the previous sections would need to be rescaled downwards by a factor $\sqrt{1+\mathcal{R}_{\textrm{max}}(T)}$. For instance, for a distance $r=50\;\textrm{kpc}$, a temperature $T_{\textrm{cold}}=5.6\times 10^{4}\;\textrm{K}$ and a volume filling factor $f_{V}=10^{-3}$, the inferred density of the cold gas decreases from $n=1.4\;\textrm{cm}^{-3}$, in the recombination-dominated case, to $n=0.23\;\textrm{cm}^{-3}$, assuming maximum scattering. The neutral hydrogen column density and optical depth should also be rescaled, this time by a factor $1+\mathcal{R}_{\textrm{max}}(T)$, due to the extra dependence of the neutral fraction of density (cf. Section 5.1).

The minimum CGM baryon fraction assuming maximum possible contribution from scattering is given in Figure 8, as a function of the assumed halo mass $M_{\textrm{halo}}$ and temperature of the cold gas $T_{\textrm{cold}}$. In this case, the fiducial values $M_{\textrm{halo}}\simeq 10^{12}\;\textrm{M}_{\odot}$ and $T_{\textrm{cold}}\simeq 5.6\times 10^{4}\;\textrm{K}$ (discussed in Sections 4.2 and 4.3, respectively) lead to a physically acceptable solution, with a minimum CGM baryon fraction equal to 70% of the cosmological value. This is still larger than the median baryon fraction predicted by many strong ejective feedback models (see references in Sections 1 and 4.4). None the less, allowing for a few $10\%$ of baryons in stars and ISM, this model is consistent with some moderate fraction of the baryons being expelled from haloes at higher redshift and not re-accreted yet. We recall that the room for expelled material gets lower if any of the two phases (cold or hot) dominates in mass over the other, as the inferred CGM baryon fraction is larger than the lower limit above if any of the two phases dominate over the other in mass (see Figure 2). On the other hand, smaller CGM baryon fractions, similar to those predicted by strong ejective feedback models, can be obtained if just one of the two parameters $M_{\textrm{halo}}$ and $T_{\textrm{cold}}$ has a value much different from fiducial (rather than both of them, as in the recombination-dominated scenario; see Section 4.4).

We end this section by discussing the self-consistency and the plausibility of a scattering-dominated scenario. As mentioned above, the neutral hydrogen column density in this case needs to be rescaled, with respect to a recombination-dominated scenario, by a factor $1+\mathcal{R}_{\textrm{max}}$. Equations (31) and (38) therefore become

[TABLE]

and

[TABLE]

respectively, where we omitted for clarity a negligible explicit dependence on temperature. Equation (45) in particular implies that a Ly$\alpha$ optically thin scenario is self-consistent out to a distance $r=100\;\textrm{kpc}$ for an intrinsic velocity dispersion as low as $\sigma=35\;\textrm{km}\;\textrm{s}^{-1}$ (assuming $\xi=1$).

We have therefore found that the recombination-dominated scenario and the scattering-dominated scenario are both formally self-consistent. As, however, they have distinctly different implications for the CGM baryon fraction (see Figures 7 and 8), it would be desirable to find additional means to discriminate between these two options.

One possible way to verify the viability of the scattering-dominated scenario could be to look for proximate absorption signatures in the QSO Ly$\alpha$ spectra. We have estimated that, in order to extract their luminosities from the BLR of the central QSOs, the MUSE Giant Ly$\alpha$ nebulae should offer to the central source, on average, a rest-frame equivalent width

[TABLE]

where we have used the median observed luminosity of the nebulae $L_{\textrm{Ly}\alpha}=10^{44}\;\textrm{erg}\;\textrm{s}^{-1}$ (Borisova et al. 2016), a median QSO Ly$\alpha$ specific luminosity $L_{\lambda,\textrm{Ly}\alpha}^{\textrm{QSO}}=3.3\times 10^{44}\;\textrm{erg}\;\textrm{s}^{-1}\;\textrm{\r{A}}^{-1}$ (again assuming the typical QSO spectrum of Lusso et al. 2015) and we finally assumed that the nebulae surround their central QSOs completely, so that $\Omega=4\pi$. We have looked for associated absorption signatures in the spectra of the radio-quiet QSO of Borisova et al. (2016) within the wavelength range defined by the FWHM of the nebular spectra and we have found an average rest-frame equivalent width $EW_{0}=0.27\;\textrm{\r{A}}$, fairly close to the expected value. We emphasize however that this kind of measurements should be taken with great caution, for a number of important reasons, including (i) the inherent difficulty of modelling absorption features on top of an emission line (here we simply modelled the ‘continuum’ with a local linear interpolation), (ii) we cannot securely distinguish absorption features arising in the CGM from those associated with the BLR itself, the ISM of the host galaxy or even the IGM, so that the measurement should probably be interpreted as an upper limit and (iii) the down-the-barrel sightline offers a very limited and possibly biased view, which may not be representative of the average. We speculate however that studies of associated absorbers on larger samples and with a more detailed modelling of kinematics and radiative transfer could possibly help to distinguish between scenarios.

Finally, we want to emphasize that there is at least one example of a Giant Ly$\alpha$ nebula (the Slug nebula; Cantalupo et al. 2014) where a scattering-dominated scenario can be excluded at very high confidence. If applied to the Slug, a calculation similar to that in equation (46) leads to a very different result than for the MUSE nebulae. This is due to (i) the relatively high luminosity of the Slug ($L_{\textrm{Ly}\alpha}=2\times 10^{44}\;\textrm{erg}\;\textrm{s}^{-1}$; Cantalupo et al. 2014) 666Note that the illuminating QSO is instead very similar to those of the sample of Borisova et al. (2016); see also Arrigoni Battaia et al. (2015). and most importantly (ii) to the fact that the Slug nebula is strongly asymmetric. Most of the luminosity arises from a region with size $l\simeq 80\;\textrm{kpc}$ at a projected distance $R\simeq 80\;\textrm{kpc}$ (Cantalupo et al. 2014). The solid angle subtended by the nebula as seen by the QSO can therefore be estimated as $\Omega\simeq 0.79\;\textrm{srd}$, or $\Omega/4\pi=0.06$ (assuming that the true distance and the projected distance coincide) or even smaller (if the true distance is larger than the projected one, as argued for instance in Cantalupo et al. 2019, though note that we have little knowledge of the structure of this nebula along the line of sight). As a consequence, the rest-frame equivalent width necessary to power the nebula with scattering is exceptionally high: $EW_{0}=9.6\;\textrm{\r{A}}$, or possibly even larger that that. Note that, to be consistent with the optically thin assumption, this would require a velocity dispersion in excess of $1000\;\textrm{km}\;\textrm{s}^{-1}$ (e.g. Meiksin 2009). We can therefore conclude that scattering cannot be the dominant powering mechanism of the Slug nebula, as also found by the radiative transfer simulations of Cantalupo et al. (2014). This is also in agreement with the recent detection, in the same object, of the non-resonant H$\alpha$ line, with a $\textrm{Ly}\alpha/\textrm{H}\alpha$ ratio consistent with recombination radiation (Leibler et al. 2018). Future H$\alpha$ observations, possible with the James Webb Space Telescope (JWST), will help to determine whether this is also the case for the MUSE Giant Ly$\alpha$ nebulae. Finally, polarization studies of Giant Ly$\alpha$ nebulae would also help constraining the contribution of scattering from a central source to the total observed emission (e.g. Hayes, Scarlata & Siana 2011).

6.2 Scattering within the nebula

Some of the Ly$\alpha$ photons, coming either from recombinations or from the first scattering of a BLR photon, can be scattered (or further scattered) within the nebula itself, thus affecting the surface brightness profile, while keeping the luminosity constant.777We neglect dust attenuation or any sort of net absorption, as they would only reinforce our conclusions; see Section 2.1.3. A difference between the intrinsic and observed surface brightness profiles would imply a difference between the intrinsic density distribution of the cold gas and the one derived neglecting internal scattering (Section 2.2). Although the total luminosity is conserved in the process, it is not obvious that the inferred total mass (which is what we are interested in here) would remain unaffected as well. Here we show that this is indeed not the case, but that the required correction is small enough to have negligible impact to our conclusions.

Suppose that the ‘intrinsic’ surface brightness distribution due to recombination is different from the observed one:

[TABLE]

[TABLE]

with $\Sigma_{0,\textrm{intr}}\neq\Sigma_{0,\textrm{oss}}$ and $\beta_{\textrm{intr}}\neq\beta_{\textrm{oss}}$. Then, using the notation of Section 2.2, the intrinsic density profile would be described by some intrinsic normalization and slope $(n_{0,\textrm{intr}},\gamma_{\textrm{intr}})$ that are different from those that we derived from the observed profile, which we can call here $(n_{0,\textrm{oss}},\gamma_{\textrm{oss}})$. We can integrate equations (47) and (48) to obtain the total luminosity and the luminosity due to recombinations, respectively:

[TABLE]

[TABLE]

where

[TABLE]

Note that in the integration above we use, as a lower limit, $R_{0}=10\;\textrm{kpc}$, which is a reasonable estimate for the CGM/ISM transition, while $R_{\textrm{max}}$ is left free, but can be assumed to be of the order of the virial radius (equation 11).

To derive the necessary corrections to the density profile, we then impose that the total luminosity be

[TABLE]

where $f_{\textrm{BLR}}$ is the global average of the ratio defined in equation (41) and will therefore vary in the range $0<f_{\textrm{BLR}}<\mathcal{R}_{\textrm{max}}$ (see equation 43). Combining equations (49), (50) and (52) with equations (8), we find

[TABLE]

The total cold gas mass is then found integrating over radius (as in Section 2.3), thence

[TABLE]

where the $\gamma$’s are as in equation (8).

The final correction factor on the cold gas mass, which then propagates to a correction to the total (cold plus hot) inferred mass of the CGM, is readily obtained by combining equations (53) and (54). The resulting expression (which we do not write explicitly as it is as long as trivial) contains a factor $(1+f_{\textrm{BLR}})^{-1/2}\leq(1+\mathcal{R}_{\textrm{max}})^{-1/2}$, which is just the same that was discussed already in Section 6.1. The remaining part of the expression gives any additional bias in the determination of the CGM mass due to internal scattering within the nebula. This is visualized in Figure 9, as a function of the intrinsic slope $\beta_{\textrm{intr}}$, for various values of the radial range parameter $x$, chosen to comprise the range of virial radii plausible for our haloes. Though the internal scattering of recombination radiation would only act in making the profile shallower, we consider here both the options $\beta_{\textrm{intr}}<\beta_{\textrm{oss}}$ and $\beta_{\textrm{intr}}>\beta_{\textrm{oss}}$, as the first scattering from the BLR can instead result in a centrally concentrated contribution in the high optical depth regime (large $\tau_{\textrm{Ly}\alpha,\textrm{QSO}}$, which however corresponds to small $\mathcal{R}$ and $f_{\textrm{BLR}}$, as discussed in Section 6.1). It is clear from the figure that the correction is of the order of a few per cent at most and therefore has no significant impact on our results.

6.3 Scattering off a neutral screen

An additional redistribution effect, which can be partly treated with the same formalism of Section 6.2, would also arise if some of the Ly$\alpha$ photons originated in the CGM and travelling away from the observer are scattered back into the line of sight after interacting with a layer or region of material with a relatively high neutral fraction and neutral hydrogen column density. This can happen, for instance, if the central QSO is obscured along some directions away from the line of sight, so that some regions of the CGM are less ionized than others and offer higher neutral hydrogen column density, as suggested, for instance, by absorption studies of the CGM (Prochaska et al. 2013) around QSOs of similar redshift as (although about 1 dex less luminous than) those of Borisova et al. (2016).

This process can change the shape of the surface brightness profile – an effect that we addressed already in Section 6.2 – but, in addition to that, it could also introduce a global normalization bias, as light is being added along some lines of sight. At the same time, however, this potential normalization enhancement would be counteracted by the fact that the whole ‘far’ side of the nebula would be largely obscured by the intervening screen and therefore not contribute to the observed Ly$\alpha$ flux. Insight suggests that the latter effect should in general be stronger than the former, so that the combination of the two is unlikely to bias the total flux high and therefore to weaken our conclusions. The idea can be intuitively illustrated as follows.

Suppose that the QSO is able to ionize only a fraction $f_{\textrm{ion}}<1$ of the entire solid angle. We include in this both the volume ionized by the QSO directly and the volume ionized indirectly by ionizing photons produced in recombination events (e.g. Osterbrock & Ferland 2006). Then only a fraction $0.5f_{\textrm{ion}}$ of the CGM volume would be bright in recombination radiation (or capable to directly see and scatter photons from the BLR) and also be directly visible by the observer. At the same time, some fraction $f_{\textrm{screen}}$ of the photons produced in this region would intercept the screen, where they can be scattered back (and possibly reach the observer) with some probability $f_{\textrm{back}}$, or find their way through the screen with a probability $1-f_{\textrm{back}}$. Also, with a probability $f_{\textrm{screen}}(1-f_{\textrm{back}})$, photons from the ‘far’ side of the CGM would intercept and traverse the screen and also contribute to the observed flux. The total observed Ly$\alpha$ flux would be different from that of an ideal case with perfect illumination, by a factor

[TABLE]

Although highly idealized, the inequality above suggests that, if anything, the presence of a neutral screen is more likely to introduce a bias in the direction of strengthening our results, rather than viceversa. We plan to demonstrate this point more thoroughly with radiative transfer simulations in future work.

7 Discussion

We discuss here the following interrelated questions: (i) What are the typical size and the number of the cold gas structures that produce the emission in MUSE Giant Ly$\alpha$ nebulae? (ii) What is the origin of the cold gas? (iii) How well justified is our assumption of pressure equilibrium with the hot gas? When relevant, we also try to connect our discussion to resolution or other numerical effects that may be important for the interpretation of cosmological hydrodynamical simulations of the CGM.

7.1 The typical size of the cold gas structures

So far, we have not addressed the question of what is the typical size and the number of the cold gas structures in the MUSE Giant Ly$\alpha$ nebulae. This is because the main parameter in our calculations is the volume filling factor $f_{V}$ and any given volume filling factor can be obtained either with a large number of small clouds or with a smaller number of larger ones. One can however attempt to break the degeneracy, at least qualitatively, by considering the area covering factor $f_{C}$, which is more easily accessible to observations than $f_{V}$ and is given by

[TABLE]

where the integral is along the line of sight and $R_{c}$ is the typical size of the individual cold gas structures along the same direction. Assuming that the intrinsic orientation of the structures is independent on our line of sight, we can take $R_{c}$ as the typical physical scale of the emitting structures, averaged over all directions; we stress that $R_{c}$ does not necessarily represent the radius of spherical clouds, as we did not make any assumptions on the geometry and not even on the topology of the emitting structures.

As an order-of-magnitude estimate, equation (56) provides

[TABLE]

where we assumed that the typical length of the line of sight through the CGM is of the order of the virial radius (equation 11) and we also assumed $f_{C}\gtrsim 1$. While we do not know what the true covering factor of the cold CGM is (due to the seeing-limited spatial resolution of MUSE), the continuity and smoothness of the emission suggest that it is not much smaller than unity. Arrigoni Battaia et al. (2015) have shown with synthetic observations that, even accounting for finite resolution, the observed morphology of extended Ly$\alpha$ emission requires an intrinsic covering factor $f_{C}>0.5$. Large covering factors are also supported by absorption studies, which are more sparse but have virtually infinite angular resolution (e.g. Prochaska et al. 2013). Note that $f_{C}>1$ is possible if multiple structures overlap in projection along the line of sight, which explains why equation (57) is best interpreted as an upper limit. We stress that the one above is just an order-of-magnitude estimate and it is only useful as a guidance to discuss which scenarios for the nature and origin of the cold gas should be considered more plausible than others.

7.2 Most of the cold gas is not in gravitationally bound substructures

Using the first equality in (57) and the definition of volume filling factor, one can also estimate (again by order of magnitude) what number of cold structures with size $R_{c}$ can be found within the virial radius:

[TABLE]

assuming $f_{V}\lesssim 10^{-2}$ and $f_{C}\gtrsim 1$. This simple estimate disfavours a scenario in which the cold gas structures are individually associated with dark matter sub-haloes (and therefore suggests that the Ly$\alpha$ emission cannot originate, for instance, from ISM bound to undetected satellite galaxies), as the latter are only a few hundreds at most (e.g. Moore et al. 1999; Rodríguez-Puebla et al. 2016). We recall that larger values of the filling factor $f_{V}$ would require an implausibly large mass of cold gas to explain the observed surface brightness (see Section 2.3); note also that this conclusion is independent on assumptions about pressure confinement.

The size estimate (57) also helps determining whether the cold gas structures are collapsed or on the verge of collapsing under their own self-gravity. This is an intriguing option, as it could imply that some of the Ly$\alpha$ emitting structures could be related to precursors of globular clusters. To this aim, we should compare $R_{c}$ to the Jeans length

[TABLE]

where again $T_{4}\equiv T/10^{4}\;\textrm{K}$, $c_{s}$ is the sound speed of an ionized plasma at temperature $T$ and we used equation (7), evaluated at a reference radius of 50 kpc. The negligible dependence on temperature originates from a combination of the temperature dependencies of the sound speed and the effective Ly$\alpha$ recombination coefficient. Combining equations (57) and (59) gives the ratio

[TABLE]

We see that the condition $r_{J}<1$ for stability is safely met for any filling factor of interest. We conclude that the cold gas structures are most likely stable against gravitational collapse. In the calculation above, we used our density estimate based on a pure recombination scenario (Section 2), but note that our conclusion would only be reinforced if scattering is non-negligible, as the estimated densities decrease in that case (Section 6.1), further enhancing stability against gravitational collapse.

7.3 Can the cold gas be overpressured?

As the gravitational pull (by either dark matter sub-haloes or the gas itself) can be neglected (see Section 7.2), we are confident that the cold gas structures are most likely confined by the pressure of the surrounding medium, which justifies the main assumption of our Section 3. We also remind that such external pressure can include non-thermal contributions, such as those due to large-scale turbulence, and that the precise amount of thermal and non-thermal pressure has little impact on our results. With this specification in mind, any residual pressure imbalance should be smeared out roughly on one sound crossing time-scale

[TABLE]

where we scaled the result to the reference temperature and size estimated by means of equations (22) and (57), respectively.

One possibility to explain a large Ly$\alpha$ surface brightness without resorting to large baryon fractions could therefore be that the majority of the cold gas structures were either created, or they acquired their current pressure, less than $\sim 2\;\textrm{Myr}$ ago. There is one simple way to achieve this. If the QSO lifetime is significantly smaller than $\sim 2\;\textrm{Myr}$, the cold gas that we see bright in Ly$\alpha$ should be modelled as having been suddenly heated to the equilibrium temperature (see Section 4.3) and may still find itself temporarily overpressured with respect to the ambient medium, by an amount $T_{\textrm{cold}}/T_{\textrm{in}}$, where $T_{\textrm{in}}$ is its temperature before being illuminated by the QSO. Quantitatively, this sudden photo-heating scenario is equivalent to rescaling the factor $\varepsilon$ in equation (16) by $T_{\textrm{in}}/T_{\textrm{cold}}$. As a result, equation (23) for the minimum CGM baryon fraction reads in this case

[TABLE]

with $T_{\textrm{in},4}\equiv T_{\textrm{in}}/10^{4}\;\textrm{K}$ and the other quantities as in Section 4.4. The minimum CGM baryon fraction for this model is shown in Figure 10 assuming an initial cold gas temperature $T_{\textrm{in}}=10^{4}\;\textrm{K}$ and recombination-dominated emission. Comparison to Figure 7 shows that sudden photo-heating facilitates the achievement of a physically acceptable solution, without requiring uncomfortable assumptions on the QSO spectrum or CGM metallicity (Section 4.3), although relatively large halo masses $M_{\textrm{halo}}\gtrsim 10^{12.4}\;\textrm{M}_{\odot}$ are still required.

Short-lived overpressures could also originate, for instance, if some of the cold gas is ballistically ejected from the ISM of a galaxy, without significantly interacting with the surrounding medium before reaching a considerable distance within the CGM. A discussion of the conditions under which this can happen would go beyond the scope of this paper. We note, however, that schemes corresponding to this description are sometimes implemented in some recipes for feedback in cosmological hydrodynamical simulations (e.g. Vogelsberger et al. 2013 and references therein). We therefore expect some amount of overpressured gas to be found in such simulations, but we warn that it can be, at least to some extent, an unphysical consequence of the adopted hydrodynamic/feedback scheme and should therefore be considered with caution.

One should finally consider the option that the cold gas structures are internally turbulent. Super-sonic turbulence would induce a lognormal density distribution within the cold gas structures (e.g. Molina et al. 2012 and references therein), leading to an additional internal clumping factor (see Section 2.1.2), which we did not account for in our calculations, though there are some recent observational indications pointing in this direction (Cantalupo et al. 2019). Significant internal clumpiness would have an impact on the estimated CGM baryon fraction of the cold gas. For instance, following Nolan, Federrath & Sutherland (2015), internal supersonic turbulence with a Mach number $M=2.5$ would induce in the cold gas a lognormal density distribution with dispersion $\sigma_{\textrm{ln}}=0.7$. This would imply an internal clumping factor $C_{\textrm{int}}=e^{\sigma_{\textrm{ln}}^{2}}=2$ and drive our mass estimates of the cold gas down by a factor of 2. Two potential problems are worth being pointed out concerning this scenario. First, internal turbulence would imply a total (thermal plus turbulent) internal pressure that is about $(1+\Gamma M^{2})$ larger than the thermal component alone (with $\Gamma=5/3$ for a monoatomic gas). This would require an equal enhancement of the external confining pressure and imply an increase of the estimated mass of the hot gas (e.g. by a factor of 11 for $M=2.5$), more than compensating the decrease in the estimated mass of the cold gas. The last conclusion does not hold, of course, if the phenomenon is short-lived and set to be smeared out on a turbulent crossing time. Note that for supersonic turbulence the turbulent crossing time would be even smaller than the sound crossing time in equation (61), by a factor of $\sim M$, so that the emitting structures should have become turbulent, in this scenario, roughly 1 Myr ago. Also note, in contrast to the ‘sudden photo-heating’ scenario discussed above, that it is not clear why supersonic turbulence should have been set on so recently. Second, supersonic turbulence would quickly heat the cold gas, through repeated shocks, to high temperatures, where the Ly$\alpha$ emissivity drops, further lowering the survivability of this situation, by an amount that we did not try to quantify. We cannot exclude that a large internal clumpiness can be due to processes other than turbulence. We are not aware, however, of plausible alternatives. We recall that self-gravity (which would in principle be an appealing candidate) can be neglected here, as we have argued already in Section 7.2.

7.4 Progressive fragmentation and incomplete pressurization

Although persistent overpressures are challenging to explain, an underpressured cold CGM can instead be achieved under somewhat easier circumstances, which we briefly discuss below.

Cold gas accreting onto a halo from IGM filaments can find itself, at the moment of entering the halo, with a lower pressure with respect to the hot gas already present there (see e.g. the phase diagrams in Nelson et al. 2016). The hot gas will then compress the incoming gas to higher densities, but the process is not necessarily well approximated as being instantaneous. It is therefore important to establish whether the incoming gas has the time to be pressurized before penetrating deep into the halo. One way to do this is to compare the sound-crossing time of a filament

[TABLE]

and its halo crossing time

[TABLE]

It is readily seen that sufficiently thick or fast filaments could in principle be able to travel significant distances within the CGM before being completely pressurized by the hot gas. Part of this gas would therefore be underpressured and emit less Ly$\alpha$ than otherwise.

It is possible and even likely, however, that pressurization proceeds much faster then the simple estimate above. Unless perfect symmetry conditions are met, in fact, an underpressured structure would be progressively deformed, during compression, and possibly fragment in smaller units. Fragments (or highly deformed portions of the original structure) will have a smaller sound-crossing time and therefore be prone to even faster compression and further fragmentation. The process has a run-away nature and results in very rapid progressive fragmentation in structures with smaller and smaller size, until complete pressurization is reached, or some other process (e.g. thermal conduction) takes over, either halting further fragmentation or destroying the fragments altogether.888This process should not be confused with the interaction between the hot gas and an inflowing filament assumed to be in pressure equilibrium since the beginning and which can then be disrupted by the Kelvin-Helmholtz instability; this can result in moderate local deviations from pressure equilibrium (e.g. Mandelker et al. 2016; Padnos et al. 2018; Vossberg et al. in prep.)

It is also possible that some of the cold gas originates from the condensation of the hot gas itself. Progressive fragmentation is indeed usually described in the context of thermal instability, in which radiative cooling is the means to create local underpressures and initiate fragmentation. The final size of the cold gas structures in this case is roughly given by the acoustic length, for which the cooling time equals the sound crossing time (e.g. Inoue & Omukai 2015). McCourt et al. (2018) recently described the phenomenon under the name of shattering and clarified that fragmentation would stop at the cooling length (the minimum acoustic length achieved during the cooling process), unless this is smaller than the Field length, below which thermal conduction would lead the cold gas to evaporate. We also notice that, if the radial profiles of heating and cooling meet certain conditions (specifically, if the heating rate declines with radius more steeply than the cooling rate), some cold gas can condensate in a shell of finite radius, possibly meeting the requirements for the onset of the Rayleigh-Taylor instability and eventually lead to the creation of filamentary structures within the CGM.999We thank A. Negri for bringing this option to our attention.

Crucially, high resolution is required for the processes mentioned above to be captured in a hydrodynamical simulation. The cooling length, in particular, can be as small as a fraction of a pc (McCourt et al. 2018). It is therefore possible that underpressured gas may be present in some cosmological simulations, as a consequence of non-negligible nominal compression time-scales (equation 63), but it is also possible that some of these features may partially be artefacts, due to the numerical resolution being insufficient to follow the shattering process and the resultant accelerated compression. Note that this last consideration should mostly apply to grid codes. SPH codes could suffer from different issues, including on one side spurious clumpiness (e.g. Torrey et al. 2012) and on the other side artificial underpressures associated to the non-trivial nature of density evaluations in SPH codes (see e.g. Appendix D in Oppenheimer et al. 2018).

We conclude by emphasizing that a systematically underpressured cold CGM, if real, would imply that the hot halo is even more massive than what inferred based on pressure balance, and would therefore go in the direction of strengthening our results.

8 Summary and conclusions

Galaxy formation models rely on mechanisms to keep most of the baryons outside of galaxies, but different scenarios make different predictions about how far away from galaxies the gas is expected to be found. Measuring what fraction of baryons is retained within the virial radius of haloes – especially at high redshift, when most of the relevant processes are expected to take place – is therefore crucial to inform and constrain theoretical models and especially different flavours of feedback.

In this work, we have tried to determine the total baryonic mass stored in the CGM of relatively massive ($M_{\textrm{halo}}\sim 10^{12}\;\textrm{M}_{\odot}$) haloes at $z\sim 3$, based on the observed surface brightness profile of Giant Ly$\alpha$ nebulae discovered with MUSE around luminous QSOs at this redshift (Borisova et al. 2016). To this aim, we modelled the CGM as a multiphase medium, where a hot (collisionally ionized) plasma confines, due to its pressure, a colder (photo-ionized) phase into a small fraction of the volume and to large densities, with correspondingly large Ly$\alpha$ emissivity.

We have found that a CGM dominated in mass by any of the two phases would fall short in reproducing the observed surface brightness, either because the cold emitting gas is not enough, or because it has too low density (and therefore too low emissivity) due to insufficient compression from the ambient medium. On the other hand, by assuming that the two phases provide similar contributions to the total mass (i.e. by maximizing the theoretical radiative output of the model), one can put a lower limit to the total mass of the CGM required to explain the observed surface brightness. Our CGM mass estimates are lower limits also because we have neglected a number of effects (dust extinction, IGM attenuation, partial illumination due to a narrow QSO ionization cone), which would diminish the predicted Ly$\alpha$ luminosity and therefore increase the amount of gas needed to explain the observations.

We have thereby found that a large total (hot plus cold) baryonic mass of the CGM – close to the total amount of baryons nominally associated with these haloes – is required to give account for a Ly$\alpha$ surface brightness as high as observed. The physical constraint that the baryon fraction be equal to or smaller than the cosmological value requires at least one of the following four conditions to be met:

The mass of QSO-hosting dark matter haloes at $z\sim 3$ is larger than the current best estimates based on clustering ($M_{\textrm{halo}}\sim 10^{12}\;\textrm{M}_{\odot}$; Section 4.2); 2. 2.

The equilibrium temperature of the photo-ionized gas is significantly smaller than what expected in a ‘standard’ QSO ionizing spectrum ($T_{\textrm{eq}}\sim 5\times 10^{4}\;\textrm{K}$; Section 4.3); 3. 3.

The intrinsic radial velocity dispersion of the CGM is large enough that the photo-ionized gas is optically thin to Ly$\alpha$ scattering (Section 6.1); 4. 4.

Most of the cold gas is in a short-lived overpressured phase, which we observe before it had the time to expand in the surrounding medium (Section 7.3).

All these options can potentially be constrained by independent observations. Particular help would come from further constraints on clustering and spectral properties of high-$z$ luminous QSOs, as well as from follow-up studies on CGM non-resonant emission lines and/or polarimetry of the Ly$\alpha$ emission. The last possibility (short-lived overpressures) is the most difficult to test, but we estimated that it could be fulfilled if the majority of bright QSOs at $z=3$ have a life-time shorter than $\sim 2\;\textrm{Myr}$.

At least two of the conditions above must be met, if the observed surface brightness is to be reconciled with the relatively small CGM baryon fractions ($30-40\%$ of the cosmic average) predicted at the relevant halo mass and redshift by the extremely strong ejective feedback sometimes adopted in galaxy formation models. More comparisons with models and simulations would be desirable to derive more robust conclusions. Regarding the comparison of QSO Giant Ly$\alpha$ nebulae with simulations, we stress in particular the importance of explicitly including photo-electric heating due to photo-ionization of HeII, as this significantly impacts the temperature, recombination coefficient and ultimately the emissivity of the gas.

In conclusion, our results seem to indicate a preference for a relatively high baryon fraction (close to cosmological) in the CGM of massive haloes at $z\sim 3$, although more observations are needed to alleviate the existing uncertainties. If confirmed, our findings could help adding constraints on fundamental parameters of galaxy formation models, such as the efficiency with which SN feedback can expel baryons from galactic haloes in the early Universe.

Acknowledgements

We thank Arif Babul, Andrea Ferrara, Max Gronke, Andrea Negri and Benny Trakhtenbrot for useful discussion, Raffaella Anna Marino for help and support and Peter Mitchell for sharing some of his analysis of the EAGLE simulations. This work was supported by the Swiss National Science Foundation grant PP00P2_163824.