Evidence for Cold-stream to Hot-accretion Transition as Traced by Ly{\alpha} Emission from Groups and Clusters at 2

We present Keck Cosmic Web Imager (KCWI) observations of giant Lya halos surrounding 9 galaxy groups and clusters at 2<z<3.3, including five new detections and one upper limit. We find observational evidence for the cold-stream to hot-accretion transition predicted by theory by measuring a decrease in the ratio between the spatially extended Lya luminosity and the expected baryonic accretion rate (BAR), with increasing elongation above the transition mass Mstream). This implies a modulation of the share of BAR that remains cold diminishing quasi-linearly (logarithmic slope of 0.97+-0.19, 5 sigma significance) with the halo to Mstream mass ratio. The integrated star-formation rates (SFRs) and AGN bolometric luminosities display a potentially consistent decrease, albeit significant only at 2.6 sigma and 1.3 sigma, respectively. The higher scatter in these tracers suggests the Lya emission might be mostly a direct product of cold accretion in these structures rather than indirect, mediated by outflows and photo-ionization from SFR and AGNs; this is also supported by energetics considerations. Below Mstream (cold-stream regime) we measure LLya/BAR=10^{40.51+-0.16}~erg/s/Msun*yr, consistent with predictions, and SFR/BAR=10^{-0.54+-0.23}: on average 30_{-10}^{+20}% of the cold streams go into stars. Above Mstream (hot-accretion regime), LLya is set by Mstream (within 0.2~dex scatter in our sample), independent of the halo mass but rising tenfold from z=2 to 3.


INTRODUCTION
It has long been understood from theory that galaxies in dark matter halos below M shock ≈ 10 12 M are fed by cold accretion, delivering gas ready to form stars (White & Frenk 1991;Birnboim & Dekel 2003;Keres et al. 2005;Dekel & Birnboim 2006;DB06 hereafter), driving high SFRs in distant galaxies (e.g. Genel et al 2008). Above M shock , cooling times are longer than dynamical times and shocks can efficiently heat incoming baryons.
However, numerical simulations and analytical work from Dekel & Birnboim (2006) first showed that cold accretion continues to penetrate at high redshifts in the form of cold streams even above M shock , in massive halos located at the intersection of multiple and dense filaments, narrower than the halos they accrete onto. This is crucial for feeding even more massive galaxies (e.g., Daddi et al. 2007) residing in massive halos at high redshifts (e.g., Bethermin et al. 2014). Subsequent numerical and analytical modelling developed this theory with respect to stream stability (Nelson et al. 2015;Mandelker et al. 2019), multi-phase properties (Cornuault et al. 2018), additional inner-halo cooling (Mandelker et al. 2019(Mandelker et al. , 2020Zinger et al. 2018), angular momentum transfer to galaxies (Danovich et al. 2015). This cold stream mode should be effective upto an evolving halo mass M stream (z), expected of order of 10 12.5 M at z = 2 and growing to 10 13.5 M at z = 3 (DB06), and rapidly diluted and disappearing at even higher masses.
Observational confirmation is still lacking for cold streams, and for evidence that the M stream transition affects gas accretion observables. Cold streams are predicted to be best detectable via their collision powered Lyα emission (Dijkstra et al 2009;Goerdt et al. 2010;Rosdhal & Blaizot 2012). However, deep Lyα observations of distant massive groups and clusters are still scarce. In this letter we present results from our ongoing KCWI survey targeting 9 massive galaxy structures at 2 < z < 3.5, providing Lyα-based evidence of the predicted dilution of cold streams across M stream . We adopt concordance cosmology (0.3; 0.7; 70) and a Chabrier IMF.

DATA AND MEASUREMENTS
We describe our sample selection and characterization of some important aspects of the structures (Table 1). A complete description of the fields will be given elsewhere (E. Daddi et al. in preparation North is up and East is left). Cl-1449 is from HST, the rest is ground-based. See D21 for a similar RO-1001 image. The blue soft layer shows Lyα emission, with contours displayed in log steps from -18.5 to -17.5 erg s −1 arcsec −2 as labeled. The dotted lines show the KCWI field. The orientation of the RO-fields was chosen to maximise overlap with the radio detections. RO-1001(Daddi et al 2021D21 hereafter), CC-0958 (Strazzullo et al. 2015), and XLSSC122 (Mantz et al. 2018). We present here two new radio overdensities (RO-0959 and RO-0958), selected following Daddi et al. (2017). We also include two Lyα blobs, SXDS-N-LAB1 (Matsuda et al. 2011;Subaru narrow-band imaging), and FVX-LAB (from our own narrow-band imaging in COSMOS).

KCWI observations and redshift identification
All the targets were observed with KCWI during observing runs in January 2018 and February 2019. The data reduction and analysis, including diffuse Lyα identification and characterisation via adaptive smoothing over a 3σ threshold follows D21's work for RO-1001. We detect giant Lyα nebulae in all structures (> 100 kpc; Figure 1), except the most massive/evolved XLSSC122 where we determine a conservative 3σ upper limit over 1000 km s −1 and 200 arcsec 2 . We confirm known nebulae: the Lyα luminosity of Cl-1449 agrees with Valentino et al. (2016), while SXDS-N-LAB1 is ×2.4 brighter than in Matsuda et al. (2011), consistent with its emission being partly redshifted out of their narrow-band filter. Fig. 1 shows three-band color-images of 8 structures newly observed in Lyα (see D21 for RO-1001). This includes a new giant Lyα halo discovery inside Cl-1001. The photometric redshift of RO-0958 and RO-0959 were confirmed by Lyα detections implying z = 3.29 and 3.09, respectively. These have been subsequently confirmed with ALMA from the detection of multiple CO lines (E. Daddi et al. in preparation). CC-0958 had tentative z ∼ 2.18 (Strazzullo et al. 2015) but KCWI revealed its giant Lyα halo at z = 2.51, still consistent with its photometric redshift.
The Lyα luminosities in Table 1 are integrated above a redshift dependent surface brightness (SB) of 2 × 10 −18 × [(1 + 2.78)/(1 + z)] 4 erg s −1 cm −2 arcsec −2 to account for SB dimming. We need (small) positive luminosity corrections only for the three z > 3 structures (Table 1), where observed SB limits are shallower than this threshold. Corrections were estimated using the other 5 profiles as a guide, with uncertainties < 0.05 dex.
The Lyα luminosities in Table 1 refer to the extended, diffuse emission only: for all structures we identified Lyα components arising from galaxies (generally AGNs, see below), and removed their contribution modelling them with the PSF as all remain unresolved at our resolution (typically 0.6-0.8 ). This correction is ∼ 15% for FVX-LAB and RO-0958, and much smaller elsewhere.

Estimates of host halo masses, SFRs and AGN content
Host halo mass (M DM ) estimates were already presented elsewhere for several structures (XLSSC122; Cl-1449; Cl-1001; RO-1001), where derivations based on their stellar mass (M * ) content were confirmed via X-ray luminosities and in two cases from the Sunyaev-Zel'dovich effect (SZ; XLSSC122 in Mantz et al. 2018;Cl-1449in Gobat et al. 2019. For the other structures we derive estimates from the M * following D21. We consider spectroscopic members, and those with consistent photometric redshifts (Laigle et al. 2016;Muzzin et al. 2013;Mehta et al. 2018), spatially coincident within the area of the structure, as gauged by the Lyα halo and self-consistently with the implied virial radius. The integrated M * above the completeness limit for the redshift is corrected to total using the mass functions in Muzzin et al. (2013), and converted to M DM using van der Burg et al. (2014). For CC-0958 we find a marginal 2.4σ detection in Chandra, fully consistent with the M * -based estimate. For the higher mass systems with M DM 4 × 10 13 M , we estimate M DM uncertainties at the level of 0.2 dex or better, at least in relative terms, by comparing the M * , X-ray and SZ derivations.
For the lower-mass systems (M DM 1×10 13 M ) we check that consistent estimates are derived using the brightest group galaxy, applying the relations in van der Burg et al. (2014) and Behroozi et al. (2013). We expect these estimates to be uncertain at the 0.3-0.4 dex level (e.g., Looser et al. 2021).
Integrated bolometric IR luminosities were presented elsewhere for Cl-1449, Cl-1001 and RO-1001 (refs in Sect. 2.1). We derived them for the other structures using Herschel PACS and SPIRE (plus other submm) observations, as described in D21 for RO-1001. Their uncertainties are small, below 0.2 dex. We set a conservative 3σ upper limit from nondetections in XLSSC122. We convert IR luminosities into SFRs following Daddi et al. (2007).
We used ancillary Chandra X-ray catalogs ( and from Fig.7 in DB06 we use M shock = 6 × 10 11 M and approximate M stream as: Fig. 2-left shows the location of our sample in the DB06 diagram, spreading across the M stream boundary. Fig.2-right shows the ratio of L Lyα to BAR for our sample as a function of M stream /M DM ratio. According to theory, cold streams should feed galaxies and halos for M DM < M stream at z > 1.4, where cold accretion is equal to total accretion as in Eq. 1, while for M DM > M stream one can expect a smooth transition where an increasingly lower fraction of accretion will be cold (DB06). We define: We model observables that are expected to be dependent on the availability of cold fuel as: where C Lyα is a constant and the slope from Eq. 3 becomes α Lyα . Similar relations are considered for integrated SFRs and L bol, AGN , with their modulation slopes α SFR , α AGN , and constants C SFR and C AGN . Fig. 2-right shows a behavior quite consistent with Eq.4. A linear fit to the data attempting to constrain its two free parameters returns α Lyα = 0.97 ± 0.19 and log C Lyα = 40.51 ± 0.16 with a scatter of 0.30 dex. Paired bootstrap (with replacements) implies similar uncertainties. The modulation of decreasing Lyα luminosity to accretion ratio α Lyα , when halo mass is larger than M stream , is hence detected at 5σ. Note that we are assuming the flattening in the right side of Fig. 2-right, not measuring it. However, if we were to extrapolate the linear fit in Fig. 2right above M stream /M DM > 1 (dashed line), the observed L Lyα there would deviate by 4.8σ in one case (RO-0958) and 1.5σ in the other two, all the three weaker than predicted by the extrapolation. This is unlikely to happen by chance, supporting the assumed flattening.
For the ratio SFR/BAR we find a consistent behavior (Fig. 3-left), but less significant (2.6σ): α SFR = 0.78 ± 0.28 and log C SFR = −0.54 ± 0.23, with a 0.45 dex scatter. Below M stream (cold-stream regime) some 20-50% of the cold accretion goes into SFR, on average (with a scatter of ×3). These fractions appear reasonable (Dekel et al. 2009), given that some reduced efficiency seems inevitable as not all the cold gas will be rapidly consumed. Above M stream (hot regime), SFR/BAR in z ∼ 2-2.5 structures is higher by ×3-10 than predicted by Behroozi et al. (2013).

Reliability of the detection
From our admittedly small and inhomogeneous sample we find observational evidence in support of the cold-stream to hot-accretion transition predicted by theory. It is crucial to assess its validity. Fig. 2 contains relations among various quantities including estimates based on observables (L Lyα and M DM ) and calculated from theory (M stream and BAR; Eqs. 1 and 2) that depend in turn on z and M DM . Redshift and Lyα luminosity errors are negligible with respect to the scatter observed in the fit. Hence, M DM is the most critical quantity in our analysis. For M DM < M stream (cold-stream regime), L Lyα /BAR is set to constant (C Lyα ; Eq. 1), because cold BAR equals total BAR. The uncertainty in C Lyα is 0.16 dex, consistent with M DM error propagation (0.3/ √ 5 dex), indirectly confirming the error estimates. In the hot regime (M DM > M stream ), the three structures with highest M DM /M stream are the ones with the best halo mass determination, estimated from X-ray and SZ measurements. By themselves, they set a slope α Lyα ≈ 1 respect to the rest of the sample. This simplifies to first order Eqs. 1 and 2 to L Lyα ∝ M stream , and we find indeed L Lyα /M stream = 10 30.75±0.20 erg s −1 M −1 with just 0.20 dex scatter, implying that L Lyα grows by 1 dex over ∆z = 1, similarly to M stream . Hence in the hot-regime (and for α ≈ 1) a precise determination of M DM is not required to measure α. Our conclusions are thus not critically dependent on M DM uncertainties.

Physical interpretation
It is important to question whether the correlation shown in Fig. 2-right arises from a direct link between cold accretion rate and Lyα luminosity (Fig. 4-right), or an indirect one, with accretion regulating SFR and L bol,AGN that in turn determine Lyα emission, e.g. by photo-ionization and subsequent recombinations (see D21 for more discussions). If the latter, it would be difficult to explain how Lyα can define a tighter relation to accretion than SFR and AGN if Lyα was a byproduct of these quantities. D21 already suggested cold accretion as the main source of Lyα powering in RO-1001, rather than AGN or SFR. However, it is worth reconsidering the matter here, with a larger sample of structures.
Starting from SFRs, basically in all cases, the most highly star-forming members are heavily dust-extinguished (Fig. 1); their contribution to Lyα photo-ionization from UVunattenuated SFR appears negligible, as in D21.
AGN photoionization should be more carefully considered as a plausible source for powering Lyα, at least in some of our structures. We estimate Lyα photo-ionisation rates from L bol,AGN (measurements or upper limits) using bolometric corrections to ultraviolet (Trakhtenbrot & Netzer 2012) and the Type1 QSO average spectrum (Lusso et al. 2015). Fig. 4-left shows that the maximum theorethical AGN ionizing radiation is potentially sufficient to power L Lyα only in 4 of the 9 structures. When considering that at these luminosities the Type 1 opening angle is expected ∼ 30% (Simpson et al. 2005; possibly too large for our sample where no unobscured AGN is found from the 9 structures, Fig. 4), and the Lyman Continuum escape ∼ 30% (Smith et al. 2020; also likely an overestimate given that our AGNs are embedded in high dust optical dephts based on ALMA detections), only two structures remain marginally viable to be fully AGN photo-ionized, neglecting further geometrical effects (Valentino et al. 2016).
Our sample is compared in Fig. 4 to Lyα nebulae selected around QSOs, from Borisova et al. (2016) and Mackenzie et al. (2021). Their Lyα photo-ionization rates to L Lyα ratios are 1-2 orders of magnitude larger, on average, with respect to structures in our sample (Fig. 4-left). Also, QSOs are sig- nificantly overluminous in their average ratio of Lyα emission to cold accretion rates (Fig. 4-right). It is thus possible that our structures' Lyα powering is not coming from photoionization, as for the QSOs.
Although these figures show scatter, and the impact of AGNs (and SFR) might vary and be somewhat larger in individual cases, the favored scenario is currently that their contribution is not dominant to Lyα in massive groups and clusters at 2 < z < 3.3, at least for those in the cold-stream regime.
In the hot-regime, if α SFR and α AGN were to be truly flatter than α Lyα , this could suggest the contributions from SFR and AGN to Lyα being increasingly important with growing M DM /M stream . Such a behaviour would be physically motivated by the longer timescales required for SFR quenching and residual gas consumption (of order 100-300 Myr; AGNs just reflecting the SFR with increased stochasticity), with respect to Lyα from accretion that is likely a more instantaneous measure. Determining α SFR (α AGN ) to sufficient precision, e.g. 5σ, would require larger samples of ∼ 50 (∼ 200) massive groups and clusters.

Predicting Lyα luminosities for dark matter halos
The good fit of Eqs. 3-4 to Fig. 2-right is encouraging in terms of using Lyα to trace accretion, hence ultimately halo masses, e.g. for unveiling dark matter halo locations in protoclusters. Using the measured C Lyα , Eq. 1-4 can be re-written in the cold-stream regime as: log L Lyα /erg s −1 ∼ 43.6+log M DM 10 13 M +2.25 log 1 + z 1 + 3 , depending (quasi)linearly on M DM . L Lyα can increase only until M DM reaches M stream . From that point on, in the hot-accretion regime, L Lyα is roughly constant depending only on redshift, regardless of how large M DM can reach (because α Lyα ∼ 1). The data (although using only 6 objects) allow us to fix the numerical parameters as: L Lyα /erg s −1 ∼ 10 42.6 ( 1 + z 1 + 1.4 ) ∼7 for M DM > M stream (6) implying that the typical Lyα nebula in the hot-regime (and saturation level for the cold-regime) is ≈ 10 43.3 erg s −1 at z = 2 and ≈ 10 44 erg s −1 at z = 3. The exponent recovered is lower than what fixed by theory (Eqs. 1 and 2), perhaps suggesting a less steep M stream redshift dependence (see discussion in Dekel et al. 2009 appendixes). These predictions can be tested with larger samples.
In conclusion, we report widespread giant Lyα nebulae in massive groups/clusters at 2 < z < 3.3, the only non-detection being the most evolved cluster. The Lyα luminosity is a smoothly decreasing fraction of the total baryonic accretion onto these massive halos, for the range where models predict that cold-streams should progressively cease feeding halos, thus supporting these models.
We thank Dawn Erb for sharing calibrations, Sebastiano Cantalupo for his CubEx code and discussions, and the referee for a constructive report. RMR acknowledges GO-15910.002 from the Space Telescope Science Institute. Data were obtained at the W. M. Keck Observatory, operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration, made possible by the generous financial support of the W. M. Keck Foundation. The authors also acknowledge the indigenous Hawaiian community and are grateful for the opportunity to collect data from the sum-mit of Maunakea.