A high baryon fraction in massive haloes at z~3
Gabriele Pezzulli Sebastiano Cantalupo




firstpage–lastpage 2019

1ES 1927+654: An AGN Caught Changing Look on a Timescale of Months
    [
    Accepted 2019 March 22. Received 2019 March 18; in original form 2018 October 1
===================================================================================







We investigate the baryon content of the circumgalactic medium (CGM) within the virial radius of M_h∼ 10^12 M_⊙ haloes at z ∼ 3, by modelling the surface brightness profile of the giant Lyα 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 to or smaller than the cosmological value) require that the cold gas occupies ≲ 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∼3 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 overpressured with respect to the ambient medium.




galaxies: haloes – intergalactic medium – galaxies: formation – galaxies: evolution







§ 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).[Here 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  (; ). 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. ). 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. ; ; ). The conditions for efficient radiative cooling are far from trivial, as they crucially depend on metallicity (e.g. ), on the intensity and shape of a meta-galactic ionizing background (e.g. ; ), on the time-scales needed to ensure collisional or photo-ionization equilibrium () and on whether local ionizing sources are consistently taken into account (). Some of these effects have been shown to be important, especially at low masses (e.g. ), 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 <cit.>, 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.  and references therein), as well as hydrodynamical cosmological simulations (e.g. ; ; ). 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. ; ; ;  ), 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. ; ;  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. ; ; ). 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. ; ; ; ). 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. ; ), cosmic rays (e.g. ; ), thermal conduction () and photo-ionization by local sources (; ).

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. ; ; ). Some of these observations (e.g. ; ) 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. ; ). In addition to observations of the hot gas, an increasing amount of data is becoming available for the colder phase (T ∼ 10^4  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. ; ; ).

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. ). 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.  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 <cit.>, who successfully applied it to the halo of the Milky Way, and was proven applicable up to the scale of galaxy clusters (e.g. ; ). Absorption data provide some constraints on the cold gas around star-forming galaxies () and QSOs () out to z ∼ 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 ≳ 2 is sufficiently bright in fluorescent Lyα 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α nebulae associated with AGN have been discovered around radio galaxies (e.g.  and references therein) and more recently around radio-quiet QSOs (e.g. ; ).

Thanks to the Multi Unit Spectroscopic Explorer (MUSE; ), <cit.> found 100% detection rate of giant Lyα nebulae around a sample of 17 radio-quiet and 2 radio-loud very luminous (L ∼ 10^44 erg s^-1) QSOs at z ∼ 3. A similar result was also found by the MUSEUM survey (). Giant Lyα 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α 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. [The 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 <cit.> 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 <cit.> and provide an estimate of the CGM baryon fraction in relatively massive (M_h∼ 10^12 M_⊙) haloes at z∼3. We find that a Lyα 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α 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 <ref> we use the observations of <cit.> to infer the density profile and the total mass of the cold CGM gas in massive haloes at z ∼ 3. In Section <ref> we estimate the mass of the hot gas based on the assumption of pressure equilibrium. In Section <ref> 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 <ref> we check the self-consistency of our model in terms of the ionization structure of the CGM. Section <ref> addresses deviations from a pure-recombination scenario and in particular the impact of Lyα scattering. In Section <ref> we discuss the possible origin and physical properties of the cold gas, including possible deviations from pressure equilibrium. A summary is given in Section <ref>.

In this paper we consistently adopt a flat cosmology with parameters H_0 = 67.74  km s^-1 Mpc^-1, Ω_m = 0.3089 and a cosmic baryon fraction Ω_b/Ω_m = 0.1573 (). Using cosmological parameters as in <cit.> has negligible impact on our results. 






§ THE COLD CGM IN MUSE QSO LYΑ NEBULAE





 §.§ Recombination radiation from the CGM





  §.§.§ Lyα recombination surface brightness

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

    Σ = ∫ f_V B(T) n^2   ,

where f_V and n are, respectively, the volume filling factor and the hydrogen number density of the Lyα emitting gas, while the coefficient B is given by

    B(T) ≡E_Lyα/4 π( 1 + 2 n_He/n_H) α^eff_Lyα(T)   ,

where E_Lyα is the energy of the Lyα line, n_He/n_H = 0.25 Y/(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 ≡ T/10^4  K)

    α^eff_Lyα = 1.6 × 10^-13  T_4^-0.96 cm^3  s^-1

is the effective Lyα recombination coefficient, as derived by fitting the calculations by <cit.> in the temperature range 5 × 10^3  K < T < 8 × 10^4  K. We adopt here case A coefficients. Note however that case A and B effective Lyα recombination coefficients differ by only 5 % (). They must not be confused with the total recombination coefficients, which differ in the two cases by a factor of ∼ 1.6 (e.g. ).




  §.§.§ Filling factor

In a multiphase CGM, hot and diffuse gas at T ∼ 10^6 - 10^7  K occupies the majority of the volume and compresses the relatively cold phase (T ∼ 1 - 5 × 10^4  K) to large densities and a small fraction of the volume. Given its dependence on density and temperature (equations <ref> and <ref>), recombination radiation will be dominated by the cold and dense gas for any filling factor larger than ∼ 10^-6. This justifies the approximation that Lyα emission is produced by the cold phase only and explains the importance of including the filling factor f_V in equation (<ref>).

We do not make assumptions here about the geometry or topology of the cold Lyα 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 <ref>.

If not stated otherwise, we interchangeably indicate with n_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

    ⟨ n ⟩ = f_V n   ,

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 (<ref>) could equivalently be replaced by C⟨ n ⟩^2, where

    C = f_V^-1

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_int on the right-hand side of equation (<ref>) (though see a discussion in Section <ref> below and see also ). 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 <ref>) – 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 <ref>). 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.




  §.§.§ Main caveats

In the rest of this section, we apply equation (<ref>) to the observed Lyα surface brightness profile of the MUSE QSO Giant Lyα nebulae, to infer the density and mass of the cold phase of the CGM around bright QSOs at z∼ 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α 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α 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 <ref> and <ref>. We argue there that they either do not affect our results, or can be taken into account by relatively simple modelling.




  §.§.§ 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 <ref> and Section <ref> 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α 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α luminosity. The first condition cannot be discarded, considering the large observational uncertainties (see e.g.  and references therein), although recent work suggests that the most luminous QSOs may have rather wide ionization cones (). 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_Lyα = 10^44 erg s^-1, one would need a star formation rate of 100  M_⊙ yr^-1 (e.g. ), assuming an implausibly large escape fraction of ionizing photons f_esc = 1, or 4 × 10^3  M_⊙ yr^-1 assuming a more realistic value f_esc≃ 2.5 % (e.g.  and references therein). For comparison, the observed stellar-to-halo mass relation (e.g. ; ) and main sequence of star-forming galaxies (e.g.  and references therein) imply that a halo with mass M_halo≃ 10^12 M_⊙ (close to the peak of star formation efficiency) at z∼3 would typically host a central galaxy with stellar mass M_⋆ = 10^10.4 M_⊙ and a star formation rate of 80  M_⊙ 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.





 §.§ De-projected cold gas density profile


The MUSE QSO Giant Lyα nebulae () have a median intrinsic (redshift-dimming corrected) surface brightness profile that is well described by a power-law function of the projected radius R:

    Σ (R) = 3.9 × 10^-4( R/10  kpc)^-1.5 erg s^-1 cm^-2 srd^-1  .

Note that the sample of <cit.> 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 (<ref>) is different from the slope reported in <cit.>. Our slope was obtained by fitting to the median of the individually dimming-corrected surface brightness profiles, improving on the rough first estimate in <cit.>. Our results are not sensitive to small changes in the slope of the surface brightness profile (see also Section <ref> 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 <ref>).

Assuming spherical symmetry (or as long as we focus on spherically averaged properties), we can deproject equation (<ref>) by means of equation (<ref>), to obtain the intrinsic radial profile of the density of the cold (Lyα-emitting) gas:

    n (r) = 0.15   f_V^-1/2 T_4^t/2( r/10  kpc)^-1.25 cm^-3  ,

where r is the spherical radius and t = 0.96 is the exponent in relation (<ref>). More in detail, the normalization and slope (n_0, γ) of the de-projected density profile (equation <ref>) are obtained from the normalization and slope (Σ_0, β) of the projected surface brightness profile (equation <ref>) by means of the relations

    .[ n_0 = f_V^-1/2( Σ_02B (T)χ(1+β)R_0)^1/2  ,;                                           ;                                 γ = 1+β2 , ].

where R_0 = 10  kpc and the dimensionless projection factor χ is defined (see also ) as

    χ(a) ≡∫_0^+∞( 1 + x^2 )^-a/2 d x   .


The derivation above, which is equivalent to a classic Abel de-projection, formally relies on the assumption that the surface brightness profile (equation <ref>) and the density profile (equation <ref>) extend to infinity. As this is not realistic, we show in Figure <ref> the effect, on the projected Lyα surface brightness profile, of introducing a truncation radius r_max in the density profile (<ref>). Predictions are shown for some values of r_max, equal to or larger than the virial radius of a 10^12 M_⊙ halo at z = 3 (r_vir≃ 74  kpc, see equation <ref> below). Truncation radii equal or larger than r_vir are all consistent with the current data. Recent results from the MUSEUM survey () possibly support a truncation radius not much larger than the virial radius. On the other hand (as already noted by ), there is at least one example of a Giant Lyα nebula whose surface brightness profile is identical to those of the MUSE nebulae and extends to considerably larger physical radii (the Slug nebula; ). 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_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 (<ref>) 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 (<ref>) on f_V and T_4 are readily understood: if a smaller fraction of the volume is occupied by Lyα 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 <ref>) and therefore again a larger density to give account of a fixed Lyα surface brightness. Note that the relation between inferred density and filling factor in equation (<ref>) is equivalent to the familiar relation between inferred average density and clumping factor ⟨ n ⟩ ∝  C^-1/2, as it follows from equations (<ref>) and (<ref>).





 §.§ The baryon fraction of the cold CGM


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

    M_cold = ∫_0^r_vir 4 π f_V m_p/X r^2 n(r) dr   ,

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:

    r_vir = 11+z( G M_halo100  Ω_m H_0^2)^1/3 = 74   M_12^1/3 kpc  ,

where M_12≡ M_halo/10^12 M_⊙ is the halo mass in units of 10^12 M_⊙ and z = 3.2 is the median redshift of the MUSE QSO Lyα nebulae. Finally, dividing M_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:

    . [                           f_CGM,cold       ≡M_cold(r_vir)(Ω_b/Ω_m) M_halo;                                     ;                                      = f_1   f_V^1/2 T_4^t/2 M_12^-γ/3  , ].

where f_1 = 7.3, t = 0.96 and γ = 1.25.  

As f_CGM,cold cannot exceed unity, equation (<ref>) 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 <ref>), in agreement with previous findings (e.g. ; ). 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α emission.





§ 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 <ref>. 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. ; ; ). 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 <ref>.

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

    ε = n_hot/n_cold = T_cold/T_hot  ,

where n_hot and T_hot are the hydrogen number density and the temperature of the hot gas, respectively. Note that, in equation (<ref>), 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_cold and T_hot are independent on radius, then the density contrast ε 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 <ref>), with a downward normalization shift, due to the fact that the hot gas is more diffuse:

    n_hot(r) = ε n_0 ( r/R_0)^-γ .

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

    T_hot = 1/γμ m_p/k_BGM_halo/r_vir  ,

where m_p is the mass of the proton, k_B the Boltzmann constant and μ = 0.593 the mean molecular weight. This suggests that temperature gradients, if present, are not dramatic.

By combining equations (<ref>), (<ref>) and (<ref>), the density contrast can be written

    ε    = γ k_B T_cold/μ m_p( 1 + z )^-1( 100  Ω_m H_0^2 )^-1/3( G M_halo)^-2/3
        = ε_1 T_4 M_12^-2/3  ,

with

    ε_1 = 1.0 × 10^-2γ/1 + z = 3.0 × 10^-3

for z = 3.2 and γ = 1.25.

The total mass of the hot gas within the virial radius is readily obtained by integrating equation (<ref>). 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 

    M_hot/M_cold = ε1 - f_V/f_V = ε_1 f_V^-1 T_4 M_12^-2/3 ,

where we have used equation (<ref>) and taken into account that f_V ≪ 1 for every viable model (cf. Section <ref>). 

The quantitative implications of equations (<ref>) and (<ref>) are presented in Section <ref>. 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. ), 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_circ(r)/V_200)^-2, where V_circ (r) is the circular speed at radius r and V_200 = V_circ (r_vir). In the situation of interest, this correction is remarkably small (less than 10% for 0.2 < r/r_vir < 1), mostly due to the relatively low concentration of dark matter haloes at this redshift (c ∼ 3.5; e.g. ). 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_circ≃ 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 (; ).

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 ε (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 (<ref>) in terms of a ratio of sound speeds squared ε = 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_s,hot→ c^2_s,hot + v^2_turb, hot. An analogous reformulation and replacement would be needed in equation (<ref>) as well, as the very same turbulent motions would also provide support to the hot gas against gravity.[viz., the equation c_s, hot^2 = G M_halo/γ r_vir (equivalent to equation <ref>) should be replaced by (c_s, hot^2 + v_turb,hot^2) = G M_halo/γ r_vir. Note that only here we are using the `isothermal sound speed' c_s ≡√(P/ρ) to simplify the notation and without loss of generality.] It is then easily seen that the final estimated contrast ε remains unaffected.[The 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 <ref>.








§ THE TOTAL CGM BARYON FRACTION


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

    .[                                                           f_CGM,tot                                  ≡M_cold+hot(r_vir)(Ω_b/Ω_m) M_halo;                                                                    ;                                                                     = f_1 ( 1 + ε_1 f_V^-1 T_4 M_12^-2/3 ) f_V^1/2 T_4^t/2 M_12^-γ/3  , ].

where we recall that f_1 = 7.3, ε_1 = 3.0 × 10^-3, t = 0.96 and γ = 1.25 (equations <ref>, <ref>, <ref> and <ref>).

The total CGM baryon fraction thus depends on three parameters: f_V, M_halo and T_cold. We discuss these dependencies below, one at a time.




 §.§ Dependence on the filling factor

In Figure <ref> 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_halo = 10^12.3 M_⊙ and a temperature of the cold gas T_cold = 10^4  K, but we stress that these values are fixed here for illustrative purposes only. The effects of varying M_halo and T_cold are discussed in detail in Sections <ref> and <ref>, respectively. In Figure <ref>, 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 (<ref>), 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 <ref>). 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 <ref>).

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_halo and T_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α surface brightness of the MUSE nebulae. We will use this minimum in Section <ref> 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.





 §.§ Dependence on halo mass

Figure <ref> 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_cold = 10^4  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α 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 <ref>).

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∼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∼3 relevant to our application, finding a broad range of results, from M_halo∼ 10^12 M_⊙ to M_halo∼ 10^13 M_⊙ (see e.g.  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 <cit.>. 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; ). These authors find a typical halo mass M_halo = 10^12 M_⊙, which is on the low side of the range of masses mentioned above. We refer to <cit.> for a discussion of why larger values had been found by the previous work of <cit.> on a smaller sample.

Despite the lack of new studies with sample sizes comparable to that of <cit.>, the debate on the mass of QSO-hosting haloes at z∼3 is still ongoing in the recent literature. For instance, <cit.> inferred a large mass M_halo∼ 10^12.9 M_⊙ based on a sample of ∼ 400 X-ray selected AGN in a range 2.9 < z < 5.5, ∼ 100 of which with spectroscopic redshift. Also, <cit.> determined halo masses with a large uncertainty M_halo = 1.70-9.83 × 10^12 h^-1 M_⊙ 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 <cit.> and <cit.> emphasize the importance of further investigations to attenuate uncertainties in this field. However, given the much smaller samples of these studies with respect to <cit.>, 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 <cit.> in the relevant redshift range contains QSOs as luminous as log(L_bol / erg s^-1) = 47.72. This overlaps with the luminosity range of the QSOs in the sample of <cit.> (47.29 < log(L_bol / erg s^-1) < 48.14, using equation 1 in ), but the latter also contains objects up to 2.6 times more luminous than the maximum luminosity probed by <cit.>. If, as some theoretical models predict, the halo mass correlates strongly with QSO luminosities (see e.g. ; ), it is possible that the QSOs of <cit.> 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. ; ; ; ; ), so that this option is not very clearly supported by the current data.

Furthermore, we notice that the MUSEUM survey () has investigated Giant Lyα nebulae around a sample of QSOs less luminous than those of <cit.>, reaching a maximum luminosity of log(L_bol / erg s^-1) = 47.87. This is very similar to the maximum luminosity probed by <cit.>, possibly making comparisons less prone to potential bias. The Giant Lyα nebulae in the MUSEUM sample have a surface brightness profile that is undistinguishable from those of <cit.> 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α 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_halo = 10^12 M_⊙ found by <cit.>. However, we will also show and discuss results obtained in a broad range of masses M_halo = 10^12.3 ± 0.5 M_⊙ (), which we chose as representative of the range of uncertainty in the value of the halo mass.

If a low value M_halo∼ 10^12 M_⊙ is confirmed by future studies, it would imply that it is not easy to explain the surface brightness of MUSE Giant Lyα 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α 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α surface brightness and specifically that the two phases contribute equally to the total mass of the CGM (Section <ref>). Equivalently, the relevant curve in Figure <ref> 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. ) and on the adopted feedback recipe (e.g. ).





 §.§ Dependence on the temperature of the cold gas


In Figure <ref> we show the dependence of the inferred total CGM baryon fraction on the assumed temperature of the cold photo-ionized gas T_cold, for a fixed halo mass M_halo = 10^12.3 M_⊙. The inferred CGM baryon fraction increases with increasing the assumed T_cold. The reason is two-fold. First, as the recombination coefficient is a declining function of temperature (equation <ref>), a larger T_cold requires a larger amount of cold gas to explain the same observed Lyα emission. Second, a larger T_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_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. ). In our case, the ionizing spectrum is completely dominated by the QSO (we will see in Section <ref> that the UVB can be neglected by orders of magnitude). In a pure hydrogen nebula, the resultant equilibrium temperature would be

    T_eq = 2 h ν_LL/3 (p+2) k_B = 2.8 × 10^4 ( p+2/3.7)^-1 K  ,

where h is the Planck's constant, ν_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 (, see also Section <ref> 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  Ry, the equilibrium temperature increases by a factor (1 + 4 A(T))/(1 + A(T)), where 

    A(T) ≡n_He/n_Hα^He_rec(T)/α^H_rec(T)  .

By using the temperature-dependent recombination coefficients as in <cit.>, we then find an equilibrium temperature

    T_eq = 5.6 × 10^4  K  .

This is large enough that no physically acceptable solution can be found for a halo mass M_halo = 10^12.3 M_⊙, as it is clear for instance from Figure <ref>. Note that the problem would only be exacerbated for M_halo = 10^12 M_⊙, as suggested by clustering measurements (Section <ref>).

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 (). 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 ). 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 ν L_ν in the energy range 2  keV < E < 100  keV) and a cut-off for E > 100  keV, though we found the details of the X-ray spectrum to be unimportant. We adopted a specific luminosity at the Lyman limit  L̃_ν, LL = 5.2 × 10^31 erg s^-1 Hz^-1 (see Section <ref> below) and a distance from the source d = 50  kpc. Results are shown in Figure <ref> for different values of gas density and metallicity. 



We see that the results are in broad agreement with the simple analytical estimate above (equation <ref>) 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 ∼ 2 × 10^4  K, as desirable from inspection of Figure <ref>, can be reached only for a rather high density n = 10  cm^-3 and solar metallicity. Note that at a distance d = 50  kpc and a temperature T = 2 × 10^4, a density n = 10  cm^-3 corresponds to a filling factor f_V = 7.3 × 10^-6 (equation <ref>), which would give an absurdly large baryon fraction even at low temperature (Figure <ref>). 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  kpc would not help in obtaining lower temperatures, as the ionization parameter increases with decreasing radius for our density slope γ = 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_⊙ 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_eq = 3 × 10^4  K, in agreement with models of HeII re-ionization (e.g. ; ). 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 <ref> 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  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  Ry than in the range 1  Ry < E < 4  Ry. However, both analytic estimate and Cloudy modelling showed in this case only a modest ∼ 25 % decrease in T_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. ; ). The predictions of this kind of model are illustrated in Figure <ref>, where we assumed the spectrum of a classical (geometrically thin, optically thick) accretion disc () with Eddington accretion around a supermassive black hole with mass M_BH = 10^8  M_⊙, ignoring in particular any comptonization effect. As expected, this model easily predicts an equilibrium temperature T_eq≃ 2 × 10^4  K for a metallicity Z = 0.1   Z_⊙ and a broad range of densities and even lower temperatures (T_eq≃ 10^4  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.  and references therein) and (ii) indirect evidence based on HeII re-ionization (e.g.  and references therein).

Motivated by the discussion above, in the following we will sometimes adopt as fiducial the value T_cold = 5.6 × 10^4  K (cf. equation <ref> and Figure <ref>). We will nonetheless also present our results for a large range of temperatures to allow for more generality.





 §.§ 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_halo and T_cold, the minimum total CGM baryon fraction. As it is readily found from equation (<ref>), this is equal to

    f_min = 2 f_1 ε_1^1/2 T_4^(1+t)/2 M_12^-(1+γ)/3

and is achieved for the filling factor f_V,eq such that cold gas and hot gas contribute equally to the total mass:

    f_V, eq = ε_1   T_4 M_12^-2/3  .

Equation (<ref>) gives the lower limit to the total CGM baryon fraction needed to explain the observations within our assumptions. It is shown in Figure <ref>, as a function of the two parameters M_halo and T_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_halo≃ 10^12 M_⊙ and T_cold≃ 5.6 × 10^4  K (determined in Sections <ref> and <ref>, 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_halo or T_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 <ref> and <ref>). 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 <ref> and <ref>, but we anticipate that they are not trivially achieved either.

Figure <ref> 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 <ref>) typically predict a CGM baryon fraction equal to ∼ 30-40 % of the cosmic average at the relevant halo mass and redshift (e.g. ; ; ; 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 <ref>, we expect that simulations with strong ejective feedback should be able to reproduce the observed Lyα surface brightness, under the assumption of fluorescent emission (recombination), only if the analysis is performed on rather massive haloes (M_halo≳ 10^12.5 M_⊙) and the cold CGM temperature is assumed to be sufficiently low (T_cold≲ 2 × 10^4  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 <ref> 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 <cit.> predicts that a typical halo with mass M_halo = 10^12 M_⊙ 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 <cit.> is all hot and therefore invisible in Lyα, unless a large fraction of it is somehow assumed to be converted into cold CGM.





 §.§ Comparison with 

A special mention is needed to the study of <cit.>. These authors analysed massive haloes at z = 2 in the Illustris simulation () and calculated the expected recombination Lyα 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α nebulae (note that the  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 <ref>, we expect simulations with strong ejective feedback to be able to reproduce a Lyα surface brightness as high as observed, as long as the assumed halo mass and gas temperature are sufficiently large and small, respectively. Indeed, <cit.> adopted a range of halo masses centred on M_halo = 10^12.5 M_⊙, which is on the high side of those discussed in Section <ref>. 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 (). <cit.>, 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_duty≪ 1 (e.g. ; ). 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  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α surface brightness in the direction needed to make the comparison with observations easier.

<cit.> 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α nebulae, as we have argued already in Section <ref>.





§ 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 <ref>), before discussing the ionization of helium (Section <ref>) and the possible contribution of collisional excitation to the Lyα emission (Section <ref>). 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 <ref>).




 §.§ Hydrogen

The QSOs in the sample of <cit.> have a median specific luminosity at the Lyman limit frequency ν_LL of L̃_ν, LL = 5.2 × 10^31 erg s^-1 Hz^-1, which we estimated following <cit.> and using the absolute i-band magnitudes of the QSOs in the <cit.> sample (see their Table 2). The corresponding hydrogen ionization rate at distance r from the QSO is

    . [                               Γ_ion (r) = 14 π r^2∫_ν_LL^+∞L_νh νσ^ion(ν)   d ν;                                        ;                                                 = 1p+3L_ν, LL σ^ion_LL4 π h r^2;                                        ;                                            = 8.8 × 10^-7ξ( r10  kpc)^-2 s^-1  , ].

where σ^ion(ν) ≃σ^ion_LL  (ν/ν_LL)^-3, with  σ^ion_LL = 6.3 × 10^-18 cm^2, is the hydrogen photo-ionization cross section and L_ν≃  L_ν, LL  (ν/ν_LL)^-p is a power-law approximation of the ionizing continuum of the QSO, with a typical slope p = 1.7 (e.g. ), while the factor

    ξ≡( L_ν, LL/5.2 × 10^31 erg s^-1 Hz^-1) ( p+3/4.7)^-1  ,

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. ) is equivalent to replacing in the expressions above 1/p+3→p + s + β/(p + s)(p+s+1), with s = 2.99 and β = 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 (<ref>) should be compared to the (case A) total hydrogen recombination rate (per proton) for a fully ionized gas:

    .[                                       Γ_rec                = α_rec(T) ( 1 + 2n_Hen_H) n;                                                                                        ;                                             = 4.8 × 10^-13  T_4^-0.75( ncm^-3)  s^-1  , ].

where again we adopted a power-law approximation to a more complex temperature dependence () over the temperature range of interest and we accounted for electrons from doubly ionized helium. It is readily seen that Γ_rec≪Γ_ion for any plausible density for CGM conditions (the equality would require n ≃ 2 × 10^4  cm^-3 at r = 100  kpc and even larger densities at smaller radii). Note also that Γ_ion exceeds the photo-ionization rate  due to the cosmic UV background (e.g. ) by several orders of magnitude at any distance plausibly associated with the CGM.

Again following <cit.>, the equilibrium hydrogen neutral fraction under these conditions is achieved on very short time-scale t_ion = Γ_ion^-1 (everywhere smaller than a few years in our case) and very well approximated by x_HI = Γ_rec / Γ_ion, which gives

    x_HI = 7.9 × 10^-8 f_V^-1/2ξ^-1 T_4^-0.27( r10  kpc)^0.75  .

Combining equations (<ref>) and (<ref>), we obtain the neutral hydrogen number density n_HI = x_HIn, which reads

    .[ n_HI = 1.1 × 10^-8  f_V^-1 ξ^-1  T_4^0.21( r10  kpc)^-0.5 cm^-3  , ].

and therefore the volume-averaged neutral hydrogen number density

    ⟨ n_HI⟩ = f_V n  = 1.1 × 10^-8 ξ^-1  T_4^0.21( r10  kpc)^-0.5 cm^-3  ,

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 <ref> 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 ⟨ n_HI⟩ is actually independent also on assumptions on the ionization state of helium (which enters in both equations <ref> and <ref>) and on possible internal clumpiness C_int within the cold gas structures (cf. Section <ref>), 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α 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 ξ (equation <ref>).

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

    N_HI (r) = 7.1 × 10^14ξ^-1 T_4^0.21( r10  kpc)^0.5 cm^-2  ,

which is everywhere much smaller than the Lyman limit column density N_H, LL = σ_ion(ν_LL)^-1 = 10^17.2 cm^-2, for any plausible value of ξ 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α scattering, on the other hand, are more complicated and are addressed in detail in Section <ref>.





 §.§ Helium

Similarly, the column density of HeII is

    N_HeII = 1.3 × 10^16ξ_He^-1 T_4^0.27( r10  kpc)^0.5 cm^-2  ,

with

    ξ_He≡ 4^-(p_1-1.7) p_1+3/p_2+3 ξ  ,

where we used again atomic coefficients as in <cit.> 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  Ry < E < 4  Ry and E > 4  Ry, respectively. If for simplicity we assume p_2 = p_1 = 1.7 (i.e. if the spectrum by  can be extrapolated above 4 Ry), we see that the column density (<ref>) is still safely below the limit for self-shielding N_HeII, crit = 6.3 × 10^17 cm^-2, implying that the entire nebula should be transparent to the helium-ionizing radiation as well. 

Note, however, that <cit.> estimated very low HeII/Lyα 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α emission is contributed by processes other than pure recombination (see Sections <ref> and <ref> below; see also ), 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. ). 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 <ref> below; see also ). 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.





 §.§ Collisional excitation


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

The local Lyα emissivity due to collisional excitation and recombination are, respectively:

    j_coll = E_Lyα/4 πα^eff_coll n_HI n_e

and

    j_rec = E_Lyα/4 πα^eff_rec n_HII n_e

(we omitted `Lyα' super-scripts or sub-scripts j's and α's for clarity) and their ratio is

    j_coll/j_rec = α^eff_coll/α^eff_recx_HI/1 - x_HI  ,

where we take x_HI from equation (<ref>) and the ratio of recombination to collisional excitation coefficients as in <cit.>. The ratio α^eff_coll / α^eff_rec increases rapidly with temperature (much more rapidly than the decrease with temperature of x_HI, see equation <ref>). To obtain a conservative upper limit to the contribution of collisional excitation, we therefore consider a fairly high temperature T = 5 × 10^4  K and also a fairly large distance from the centre r = 100  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:

    j_coll/j_rec < 1.4 × 10^-2ξ^-1 f_V^-1/2

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).






§ THE IMPACT OF SCATTERING


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




 §.§ 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α 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α 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 <ref>) into a line-centre optical depth to Lyα scattering (following e.g. ). We find

    τ_Lyα, QSO (r) = 0.99  ξ^-1 T_4^0.21σ_300^-1( r/10  kpc)^0.5  ,

where σ_300 is the 1D velocity dispersion of the cold CGM, as seen by the QSO, in units of 300 km 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 σ that we used in equation (<ref>) corresponds to the typical width of the Lyα line as observed in emission (). Note however that this can be a biased estimate of the true value of σ 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 σ 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α scattered photons per fluorescent Lyα 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α photons are produced, with some probability η, 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α emission) in a time equal to the inverse of the photo-ionization rate t_ion = Γ_ion^-1. During this time, the neutral atom is available to undergo a number of scattering events, statistically equal to Γ_scatt t_ion = Γ_scatt/Γ_ion, where Γ_scatt is the scattering probability per unit time

    Γ_scatt = ∫_0^+∞F_ν (ν)/h νσ_scatt (ν) d ν = F_ν, Lyα/hν_Lyασ̅  ,

where σ̅ = ∫_0^+∞σ_scatt (ν) d ν = 1.1 × 10^-2 cm^2  Hz is the frequency-integrated cross-section to Lyα scattering (e.g. ; ) and

    F_ν, Lyα = L_ν, Lyα/4 π r^2 e^-τ_Lyα, QSO (r)

is the specific flux density of Lyα photons propagating directly from the BLR to distance r from the QSO. The cumulative Lyα optical depth τ_Lyα, QSO (r) describes attenuation by scattering at smaller radii.[Obviously attenuation here is not synonym to absorption. The fate of intercepted BLR photons after first scattering is addressed in Section <ref>.] 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 <ref>). 

The local ratio of the scattering to recombination Lyα emissivity is given by

    j_scatt/j_rec = Γ_scatt/ηΓ_ion = p+3/ησ̅/σ^ion_LLν_LyαF_ν, Lyα/F_ν, LL = 0.71 p+3/ηF_ν, Lyα/F_ν, LL  ,

where the case A Lyα production efficiency per recombination is well approximated by

    η(T) = 0.41   T_4^-0.23  .

It is clear from equations (<ref>) and (<ref>) that a very conservative upper limit for the ratio of scattering to recombination emissivity is found in the Lyα optically thin limit τ_Lyα, QSO≪ 1. Furthermore, in this limit, the ratio j_scatt/j_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α 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 <cit.> and a UV slope p = 1.7 we then find

    ℛ≡j_scatt/j_rec≤ℛ_max (T) = 26   T_4^0.23  .

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α 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 + ℛ_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 √(1 + ℛ_max(T)). For instance, for a distance r = 50  kpc, a temperature T_cold = 5.6 × 10^4  K and a volume filling factor f_V = 10^-3, the inferred density of the cold gas decreases from n = 1.4  cm^-3, in the recombination-dominated case, to n = 0.23  cm^-3, assuming maximum scattering. The neutral hydrogen column density and optical depth should also be rescaled, this time by a factor 1 + ℛ_max(T), due to the extra dependence of the neutral fraction of density (cf. Section <ref>).

The minimum CGM baryon fraction assuming maximum possible contribution from scattering is given in Figure <ref>, as a function of the assumed halo mass M_halo and temperature of the cold gas T_cold. In this case, the fiducial values M_halo≃ 10^12 M_⊙ and T_cold≃ 5.6 × 10^4  K (discussed in Sections <ref> and <ref>, 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 <ref> and <ref>). 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 <ref>). 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_halo and T_cold has a value much different from fiducial (rather than both of them, as in the recombination-dominated scenario; see Section <ref>).



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 + ℛ_max. Equations (<ref>) and (<ref>) therefore become

    N_HI (r) = 2.7 × 10^13ξ^-1( r10  kpc)^0.5 cm^-2  ,

and

    τ_Lyα, QSO (r) = 3.7 × 10^-2 ξ^-1σ_300^-1( r/10  kpc)^0.5  ,

respectively, where we omitted for clarity a negligible explicit dependence on temperature. Equation (<ref>) in particular implies that a Lyα optically thin scenario is self-consistent out to a distance r = 100  kpc for an intrinsic velocity dispersion as low as σ = 35  km s^-1 (assuming ξ = 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 <ref> and <ref>), 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α spectra. We have estimated that, in order to extract their luminosities from the BLR of the central QSOs, the MUSE Giant Lyα nebulae should offer to the central source, on average, a rest-frame equivalent width 

    EW_0 = L_Lyα/L_λ, Lyα^QSO( Ω/4 π)^-1 =  0.3  Å  ,

where we have used the median observed luminosity of the nebulae L_Lyα = 10^44 erg s^-1 (), a median QSO Lyα specific luminosity L_λ, Lyα^QSO = 3.3 × 10^44 erg s^-1 Å^-1 (again assuming the typical QSO spectrum of ) and we finally assumed that the nebulae surround their central QSOs completely, so that Ω = 4 π. We have looked for associated absorption signatures in the spectra of the radio-quiet QSO of <cit.> 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  Å, 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α nebula (the Slug nebula; ) where a scattering-dominated scenario can be excluded at very high confidence. If applied to the Slug, a calculation similar to that in equation (<ref>) leads to a very different result than for the MUSE nebulae. This is due to (i) the relatively high luminosity of the Slug (L_Lyα = 2 × 10^44 erg s^-1; ) [Note that the illuminating QSO is instead very similar to those of the sample of <cit.>; see also <cit.>.] 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 ≃ 80  kpc at a projected distance R ≃ 80  kpc (). The solid angle subtended by the nebula as seen by the QSO can therefore be estimated as Ω≃ 0.79  srd, or Ω/4 π = 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 , 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  Å, 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  km s^-1 (e.g. ). 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 <cit.>. This is also in agreement with the recent detection, in the same object, of the non-resonant Hα line, with a Lyα/Hα ratio consistent with recombination radiation (). Future Hα observations, possible with the James Webb Space Telescope (JWST), will help to determine whether this is also the case for the MUSE Giant Lyα nebulae. Finally, polarization studies of Giant Lyα nebulae would also help constraining the contribution of scattering from a central source to the total observed emission (e.g. ).




 §.§ Scattering within the nebula


Some of the Lyα 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.[We neglect dust attenuation or any sort of net absorption, as they would only reinforce our conclusions; see Section <ref>.] 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 <ref>). 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:

    Σ_oss (R) = Σ_0, oss( R/R_0)^-β_oss


    Σ_int (R) = Σ_0, intr( R/R_0)^-β_intr

with Σ_0, intr≠Σ_0, oss and β_intr≠β_oss. Then, using the notation of Section <ref>, the intrinsic density profile would be described by some intrinsic normalization and slope (n_0, intr, γ_intr) that are different from those that we derived from the observed profile, which we can call here (n_0, oss, γ_oss). 
We can integrate equations (<ref>) and (<ref>) to obtain the total luminosity and the luminosity due to recombinations, respectively:

    L_tot = 8 π^2 Σ_0, oss R_0^2/2 - β_oss( x^2 - β_oss - 1 )   ,


    L_rec = 8 π^2 Σ_0, int R_0^2/2 - β_int( x^2 - β_int - 1 )   ,

where

    x ≡R_max/R_0  .

Note that in the integration above we use, as a lower limit, R_0 = 10  kpc, which is a reasonable estimate for the CGM/ISM transition, while R_max is left free, but can be assumed to be of the order of the virial radius (equation <ref>).

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

    L_tot = (1 + f_BLR) L_rec  ,

where f_BLR is the global average of the ratio defined in equation (<ref>) and will therefore vary in the range 0 < f_BLR < ℛ_max (see equation <ref>). Combining equations (<ref>), (<ref>) and (<ref>) with equations (<ref>), we find

    n_0, intr/n_0, oss = 1/√(1+f_BLR) ( χ(1+β_oss)/χ(1+β_intr) 2 - β_intr/2 - β_oss x^2-β_oss-1/x^2-β_intr-1)^1/2  .

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

    M_intr/M_oss = 3 - γ_oss/3 - γ_intr n_0, intr/n_0, oss x^γ_oss-γ_intr  ,

where the γ's are as in equation (<ref>).

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 (<ref>) and (<ref>). The resulting expression (which we do not write explicitly as it is as long as trivial) contains a factor (1 + f_BLR)^-1/2≤ (1 + ℛ_max)^-1/2, which is just the same that was discussed already in Section <ref>. 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 <ref>, as a function of the intrinsic slope β_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 β_intr < β_oss and β_intr > β_oss, as the first scattering from the BLR can instead result in a centrally concentrated contribution in the high optical depth regime (large τ_Lyα, QSO, which however corresponds to small ℛ and f_BLR, as discussed in Section <ref>). 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.







 §.§ Scattering off a neutral screen


An additional redistribution effect, which can be partly treated with the same formalism of Section <ref>, would also arise if some of the Lyα 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 () around QSOs of similar redshift as (although about 1 dex less luminous than) those of <cit.>.

This process can change the shape of the surface brightness profile – an effect that we addressed already in Section <ref> – 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α 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_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. ). Then only a fraction 0.5 f_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_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_back, or find their way through the screen with a probability 1 - f_back. Also, with a probability f_screen(1- f_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α flux would be different from that of an ideal case with perfect illumination, by a factor

    . [                                                      F_ossF_ideal = 12 f_ion (1 + f_screen f_back) + 12 f_ion f_screen (1 - f_back);                                                                                                                                    ;                                                                                                 = 12 f_ion ( 1 + f_screen ) ≤ 1   . ].

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.






§ 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α 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.




 §.§ 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α 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

    f_C = ∫f_V/R_c dl   ,

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 (<ref>) provides

    R_c ≃ f_V r_vir f_C^-1≲ 74 ( f_V/10^-3)  pc  ,

where we assumed that the typical length of the line of sight through the CGM is of the order of the virial radius (equation <ref>) and we also assumed f_C ≳ 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. <cit.> have shown with synthetic observations that, even accounting for finite resolution, the observed morphology of extended Lyα 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. ). Note that f_C > 1 is possible if multiple structures overlap in projection along the line of sight, which explains why equation (<ref>) 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.





 §.§ Most of the cold gas is not in gravitationally bound substructures

Using the first equality in (<ref>) 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:

    N ≃ f_V^-2 f_C^3≳ 10^4   ,

assuming f_V ≲ 10^-2 and f_C ≳ 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α emission cannot originate, for instance, from ISM bound to undetected satellite galaxies), as the latter are only a few hundreds at most (e.g. ; ). 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 <ref>); note also that this conclusion is independent on assumptions about pressure confinement.

The size estimate (<ref>) 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α emitting structures could be related to precursors of globular clusters. To this aim, we should compare R_c to the Jeans length

    λ_J = √(π c_s^2/G ρ) = 16   T_4^0.02 f_V^1/4 kpc  ,

where again T_4 ≡ T / 10^4  K, c_s is the sound speed of an ionized plasma at temperature T and we used equation (<ref>), 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α recombination coefficient. Combining equations (<ref>) and (<ref>) gives the ratio

    r_J ≡R_c/λ_J≲ 4.6   T_4^-0.02 f_V^3/4  .

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 <ref>), but note that our conclusion would only be reinforced if scattering is non-negligible, as the estimated densities decrease in that case (Section <ref>), further enhancing stability against gravitational collapse.





 §.§ 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 <ref>), 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 <ref>. 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

    t_sound = R_c/c_s = 2.0  ( T_cold/5.6 × 10^4  K)^-1/2( R_c/74  pc)  Myr  ,

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

One possibility to explain a large Lyα 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 ∼ 2  Myr ago. There is one simple way to achieve this. If the QSO lifetime is significantly smaller than ∼ 2  Myr, the cold gas that we see bright in Lyα should be modelled as having been suddenly heated to the equilibrium temperature (see Section <ref>) and may still find itself temporarily overpressured with respect to the ambient medium, by an amount T_cold/T_in, where T_in is its temperature before being illuminated by the QSO. Quantitatively, this sudden photo-heating scenario is equivalent to rescaling the factor ε in equation (<ref>) by T_in/T_cold. As a result, equation (<ref>) for the minimum CGM baryon fraction reads in this case

    f_min = 2 f_1 ε_1^1/2 T_in, 4^1/2 T_4^t/2 M_12^-(1+γ)/3  ,

with T_in, 4≡ T_in / 10^4  K and the other quantities as in Section <ref>. The minimum CGM baryon fraction for this model is shown in Figure <ref> assuming an initial cold gas temperature T_in = 10^4  K and recombination-dominated emission. Comparison to Figure <ref> 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 <ref>), although relatively large halo masses M_halo≳ 10^12.4 M_⊙ 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.  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.  and references therein), leading to an additional internal clumping factor (see Section <ref>), which we did not account for in our calculations, though there are some recent observational indications pointing in this direction (). Significant internal clumpiness would have an impact on the estimated CGM baryon fraction of the cold gas. For instance, following <cit.>, internal supersonic turbulence with a Mach number M = 2.5 would induce in the cold gas a lognormal density distribution with dispersion σ_ln = 0.7. This would imply an internal clumping factor C_int  = e^σ_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 + Γ M^2) larger than the thermal component alone (with Γ = 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 (<ref>), by a factor of ∼ 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α 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 <ref>.





 §.§ 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 ). 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

    t_sound, fil = R_fil/c_s = 640   T_4^-1/2( R_fil/10  kpc)  Myr

and its halo crossing time

    t_cross, fil = r_vir/V_fil = 980 ( r_vir/100  kpc)  ( V_fil/100  km s^-1)^-1Myr  .

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α 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.[This 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. ; ; 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. ). <cit.> 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.[We 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 (). 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 <ref>), 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. ) 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 ).

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.





§ 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_halo∼ 10^12 M_⊙) haloes at z∼3, based on the observed surface brightness profile of Giant Lyα nebulae discovered with MUSE around luminous QSOs at this redshift (). 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α 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α 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α 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 ∼ 3 is larger than the current best estimates based on clustering (M_halo∼ 10^12 M_⊙; Section <ref>);

  * The equilibrium temperature of the photo-ionized gas is significantly smaller than what expected in a `standard' QSO ionizing spectrum (T_eq∼ 5 × 10^4  K; Section <ref>);

  * The intrinsic radial velocity dispersion of the CGM is large enough that the photo-ionized gas is optically thin to Lyα scattering (Section <ref>);

  * 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 <ref>).

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α 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 ∼ 2  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α 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∼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.

mn2e



