APEX at the QSO MUSEUM: molecular gas reservoirs associated with $z\sim$3 quasars and their link to the extended Ly$\alpha$ emission

Cool gas (T$\sim$10$^{4}$~K) traced by hydrogen Ly$\alpha$ emission is now routinely detected around $z\sim3$ quasars, but little is known about their molecular gas reservoirs. Here, we present an APEX spectroscopic survey of the CO(6-5), CO(7-6) and [CI](2-1) emission lines for 9 quasars from the QSO MUSEUM survey which have similar UV luminosities, but very diverse Ly$\alpha$ nebulae. These observations ($\langle~\rm rms~\rangle=2.6$~mJy in 300~km~s$^{-1}$) detected three CO(6-5) lines with 3.4$\leq I_{\rm CO(6-5)} \leq$5.1~Jy~km~s$^{-1}$, 620$\leq$FWHM$\leq$707~km~s$^{-1}$, and three [CI](2-1) lines with 2.3$\leq I_{\rm [CI](2-1)} \leq$15.7~Jy~km~s$^{-1}$, 329$\leq$FWHM$\leq$943~km~s$^{-1}$. For the CO and [CI] detected sources, we constrain the molecular gas reservoirs to be $\rm M_{H_{2}} = (0.4-6.9) \times 10^{11} M_{\odot}$, while the non-detections imply $\rm M_{H_{2}}<1.1\times 10^{11} M_{\odot}$. We compare our observations with the extended Ly$\alpha$ properties to understand the link between the cool and the molecular gas phases. We find large velocity shifts between the bulk of Ly$\alpha$ and the molecular gas systemic redshift in five sources (from $\sim$-400 to $\sim+$1200~km~s$^{-1}$). The sources with the largest shifts have the largest Ly$\alpha$ line widths in the sample, suggesting more turbulent gas conditions and/or large-scale inflows/outflows around these quasars. We also find that the brightest ($I_{\rm[CI](2-1)}=15.7\pm3.7~\rm Jy~km~s^{-1}$) and the widest (FWHM$\sim$900~km~s$^{-1}$) lines are detected for the smallest and dimmest Ly$\alpha$ nebulae. From this, we speculate that host galaxy obscuration can play an important role in reducing the ionizing and Ly$\alpha$ photons able to escape to halo scales, and/or that these systems are hosted by more massive halos.


INTRODUCTION
Super-massive black holes are found at the centre of massive galaxies (e.g., Richstone et al. 1998;Ferrarese & Merritt 2000;Kauffmann & Haehnelt 2000;Kormendy & Ho 2013).They become visible as extremely luminous active galactic nuclei (AGN) through episodes of intense accretion across the history of the Universe (e.g., Schmidt 1963;Bañados et al. 2018;Lyke et al. 2020).Because of the large budget of rest-mass energy available, these objects can regulate their own growth and the evolution of their host galaxies, even if only a small fraction of their feedback couples efficiently to the surrounding material (e.g., Silk & Rees 1998;Di Matteo et al. 2005;Steinborn et al. 2015).
★ E-mail: nahir@mpa-garching.mpg.deWithin AGN, quasars are the most luminous sources where we can see the nuclear emission directly (e.g., Antonucci 1993;Elvis 2000).Clustering measurements suggest that quasars preferentially inhabit dark matter halos with masses  DM ∼ 10 12 −10 13 M (e.g., Porciani et al. 2004;Shen et al. 2007;White et al. 2012;Timlin et al. 2018 and references therein).This mass range should guarantee that a non-negligible fraction of cool ( ∼ 10 4 K) gas, inflowing from large intergalactic scales at redshifts  2 does not shock heat at the halo boundary, but accretes in cold form (e.g., Dekel & Birnboim 2006).Quasars at such epochs are therefore expected to sit in halos with both a cool and a warm/hot gas ( ∼ 10 5 − 10 7 K) phase (e.g., Kereš et al. 2005).Because the quasar number density peaks between  ∼ 2 and  ∼ 3 (e.g., Richards et al. 2006;Shen et al. 2020), these epochs (∼ 10.4 − 11.6 Gyr ago) are frequently targeted by observations to understand how quasars are triggered and which reservoirs, from halo to galaxy scales, sustain their central engines.
The halo gas, known as circumgalactic medium (CGM, e.g., Tumlinson et al. 2017), has been studied around quasars mostly targeting the cool phase both in absorption (e.g., Hennawi et al. 2006;Prochaska et al. 2013;Farina et al. 2013Farina et al. , 2014;;Lau et al. 2018) and in emission (e.g., Heckman et al. 1991;Bunker et al. 2003;Hennawi & Prochaska 2013;Farina et al. 2019;Fossati et al. 2021).While the absorption technique usually relies on only one background sightline per foreground halo to provide statistical information on the physical properties of the CGM of quasars, studies of the CGM in emission are currently able to map the quasar CGM around individual systems.At  ∼ 2 − 3, the study of projected quasar pairs has led to a number of new insights: (i) the measurement of the anisotropic clustering of H systems around quasars (Hennawi & Prochaska 2007;Jalan et al. 2019) suggested that their ionizing radiation escapes anisotropically or intermittently, (ii) the discovery of large reservoirs (> 10 10 M ) of cool and metal-enriched ( 0.1 Z ) halo gas (Prochaska et al. 2013(Prochaska et al. , 2014;;Lau et al. 2016), and (iii) the study of the kinematics of the halo, which seems to suggest that the gas is in virial equilibrium with the dark matter halo, though there is some evidence for outflowing gas (Prochaska et al. 2014;Lau et al. 2018).
In recent years, sensitive integral field unit spectrographs like the Multi-Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010), the Keck Cosmic Web Imager (KCWI; Morrissey et al. 2012) and the Palomar Cosmic Web Imager (PCWI; Matuszewski et al. 2010) revolutionized the study of CGM gas through emission lines by allowing deeper observations in reasonable amount of time.The seminal papers by Rees (1988) and Haiman & Rees (2001), predicted that gas surrounding quasars reprocesses the impinging strong UV radiation as Ly emission.Current studies routinely report extended Ly emission, with ∼ 200 quasars surveyed to date at 2 <  < 4 (e.g., Husband et al. 2015;Fumagalli et al. 2016;Borisova et al. 2016;Arrigoni Battaia et al. 2019a;Cai et al. 2019;O'Sullivan et al. 2020;Mackenzie et al. 2021;Fossati et al. 2021).The bulk of the extended emission traces gas on a few tens of kpc near the quasars, while large-scale structures extending to  80 kpc are seen at lower surface brightness (SB Ly ∼ 10 −18 erg s −1 cm −2 arcsec −2 ).
These studies reveal few extended structures over hundreds of kpc with SB Ly 10 −17 erg s −1 cm −2 arcsec −2 (Arrigoni Battaia et al. 2018a), likely pinpointing very dense environments (Hennawi et al. 2015, Nowotka et. al subm.).These rare and bright large-scale nebulae are also known as enormous Ly nebulae (ELAN; Cai et al. 2017).The Ly kinematics in the extended nebulae is consistent with gravitational motions in halos with masses consistent with typical quasar hosts (e.g., Arrigoni Battaia et al. 2019a;O'Sullivan et al. 2020), with a few exceptions with possible quasar winds extending over tens of kpc (Travascio et al. 2020).
While a large fraction of these Ly nebulae are likely powered by the quasars, the balance between different plausible mechanisms is still debated.Most previous work assumes that the Ly emission is due to recombination radiation following quasar photoionization (e.g., Heckman et al. 1991;Cantalupo et al. 2014).Resonant scattering of quasar Ly photons and active companions can, however, provide a non-negligible contribution on scales of tens of kpc near compact sources (e.g., Cantalupo et al. 2014;Husemann et al. 2018;Arrigoni Battaia et al. 2019b).On top of this, there are large uncertainties on the ionizing radiation that impinges on the surrounding gas, because quasars are expected to be anisotropic, intermittent sources with different degrees of obscuration.These uncertainties hamper the physical interpretation of properties of the emitting gas (e.g., density  H , metallicity; Fossati et al. 2021).
There is, however, evidence in few systems that gas at large projected distances ( 50 kpc) is not affected by resonant scattering effects, namely: (i) non-resonant lines follow the kinematics of the Ly emission (e.g., He , Cai et al. 2017), and (ii) there is no evidence for double-peaked line profiles at the current resolution of the observational data (e.g., Arrigoni Battaia et al. 2018a).Neglecting resonant scattering, photoionization models match the observed Ly and low HeII mission only if interstellar-medium-like densities ( h 1 cm −3 ) in small-scale structures (< 20 pc) are invoked (Cantalupo et al. 2014;Hennawi et al. 2015;Arrigoni Battaia et al. 2015;Borisova et al. 2016).This finding suggests the presence of dense CGM gas whose survival and entrainment in the warm/hot halo seem plausible from current high resolution "cloud-crushing" simulations(e.g., Gronke & Oh 2018, 2020;Kanjilal et al. 2021).Note that such processes are still largely unresolved by current cosmological simulations, even when attempts are made to resolve the CGM (e.g., Hummels et al. 2019;Peeples et al. 2019).
It is therefore of interest to ascertain observationally the maximum density the cool CGM gas is able to reach, and whether a fraction of the gas is able to transform into a molecular phase.The molecular gas around quasars can be best probed through different tracers depending on its physical properties (e.g., density, temperature) and those of the surrounding environment (e.g., radiation field) (e.g., Carilli & Walter 2013).Most previous works have focused on the rotational () transitions of carbon monoxide 12 C 16 O (hereafter CO), which is the most abundant molecule after H 2 .Low- CO transitions are good tracers of the total cold molecular gas due to their low excitation temperatures.The CO(=1-0) ground transition requires an excitation temperature of only ∼5.5K (e.g., Bolatto et al. 2013).Using observations of different CO transitions and radiative transfer models (e.g., large velocity gradient, LVG; e.g., van der Tak et al. 2007), it is possible to constrain the CO spectral line energy distribution (SLED), and probe the excitation conditions and physical properties of the gas, as the density and kinetic temperature (e.g., Weiß et al. 2007b;Riechers et al. 2009).
Currently, there is only tentative evidence for extended molecular gas reservoirs around individual quasars, but only few studies attempted long integrations.Riechers et al. (2006) presented CO(1-0) detections in three quasars at  ∼ 4. Using single component LVG models, they found that all the flux detected in CO(1-0) was associated with the molecular gas traced by higher CO transitions.An extended component up to 30% of the total CO(1-0) luminosity was allowed by the observations.The extended component could have larger mass if the CO conversion factor was taken to be higher on larger scales.Emonts et al. 2019 targeted the CO(1-0) transition from the MAMMOTH-I ELAN located at  ∼ 2 (Cai et al. 2017;Arrigoni Battaia et al. 2018b) and reported emission extended over tens of kpc, with roughly 50% of the CO(1-0) emission outside of galaxies.Finally, Decarli et al. (2021) targeted the CO(3-2) transition for two  ∼ 2 ELANe, the Slug (Cantalupo et al. 2014) and the Jackpot (Hennawi et al. 2015).Their NOEMA observations did not unveil any extended molecular reservoir in these objects down to molecular gas surface densities typical of starbursting systems (Σ H 2 < 12 − 68 M pc −2 ).
Fine structure lines of atomic carbon, for example [C ], are an additional tracer to probe the cold molecular phase (e.g., Papadopoulos & Greve 2004;Valentino et al. 2018).Observational studies in the local Universe have shown that CO and [C ] can coexist, suggesting that both transitions arise from the same regions (e.g., White et al. 1994;Ikeda et al. 2002;Israel & Baas 2002), though spatial variations could be present (e.g., Salak et al. 2019).Analysis using simultaneously [C ] and multi-transition CO observations at high redshifts found agreement between the H 2 masses determined through the two different tracers (e.g., Weiß et al. 2003;Alaghband-Zadeh et al. 2013), corroborating the assumption that [C ] and CO usually coexist (Carilli & Walter 2013).The carbon masses found in the literature for quasars at  ∼ 2 − 3 are typically of the order of 10 6 − 10 7  and do not differ significantly from those found for quasars at  > 3 (e.g., Weiß et al. 2003;Walter et al. 2011;Schumacher et al. 2012;Venemans et al. 2017;Banerji et al. 2018;Yang et al. 2019).Molecular masses are then usually obtained by assuming the same abundance of [C ] relative to H 2 as found in high-z quasars (e.g., Weiß et al. 2005).
In this framework, we targeted the CO(6-5) ( rest = 691.4731GHz), CO(7-6) ( rest = 806.6518GHz) and [C ] 3 P 2 -3 P 1 (hereafter [C ](2-1),  rest = 809.3420GHz) transitions with the SEPIA180 receiver (Belitsky et al. 2018b,a) on the Atacama Pathfinder Experiment (APEX) for a sample of nine z∼3 quasars, whose halo gas has been studied in the QSO MUSEUM survey (Arrigoni Battaia et al. 2019a).With these observations, we aim to (i) constrain the molecular phase around these massive systems and thus start characterizing the multiphase nature of the halo gas, and (ii) investigate the relation between the molecular gas content and the large-scale cool phase.The molecular line detections also pin down the systemic redshift of the quasar very accurately, allowing us to probe the kinematics of the halo gas.This work is structured as follows.In Section 2, we describe our sample, observations and data reduction.In Section 3 we present the observed line properties and refine the systemic redshift when possible.In Section 4, we describe the estimation of the molecular gas masses using different methods, and present results for these masses.In Section 5, we compare the derived molecular masses with the Ly properties of our sources.In Section 6 we discuss our main results, and explore the link between the molecular gas content and the large scale Ly emission.Finally, Section 7 summarizes our findings.
Throughout this paper, we adopt the cosmological parameter

SAMPLE AND OBSERVATIONAL DATA
Our sample is composed of nine quasars at  ∼3 selected from the QSO MUSEUM survey (Arrigoni Battaia et al. 2019a), which targeted 61 quasars with MUSE for the study of their CGM restframe UV line emission.The nine quasars were observed with the SEPIA180 1 receiver mounted on the APEX antenna, located in Llano de Chajnantor, Chile.The targets were selected using the following constraints: • Visibility from the telescope site and the presence of CO rotational transitions, CO(6-5) and CO(7-6), and the [C ](2-1) transition within the frequency range covered by the SEPIA180 instrument.
• The expected frequency of the targeted emission lines was required to be located far from the atmospheric 183 GHz waterabsorption feature to best exploit the sensitivity of the SEPIA180 instrument, even under high water vapor conditions.
• One of the targets was selected to be radio-loud, reflecting a similar fraction of such objects in the parent quasar sample (∼10%, Ivezić et al. 2002).
The CO(6-5) ( obs 2 = 167.3GHz) observations were carried out between October and December of 2018 under the ESO programme 0102.A-0394A (PI: F. Arrigoni Battaia), with a total of ∼133 hours of telescope time.The [C ](2-1) ( obs = 195.8GHz) and 1 The APEX/SEPIA180 dual polarisation 2SB receiver is a pre-production version of the ALMA Band 5 receiver, and covers the frequency range 159-211 GHz. 2 Observed frequency of the line transition at =3.133, the median redshift of the sample studied in this paper.(𝜈 obs = 195.2 GHz) observations were performed between May and December of 2019 under the ESO programme 0103.A-0306A (PI: F. Arrigoni Battaia), with a total of ∼140 hours of telescope time.The main beam full width half maximum (FWHM) of the SEPIA180 receiver is about ∼32" (∼249 kpc) for the CO(6-5) observations and ∼ 30" (∼234 kpc) for the [C ](2-1) -CO(7-6) observations.Figure 2 shows Ly images of the nine targets analysed in this work, with superimposed APEX beams shown as white circles.The acquired data will therefore provide an integrated spectrum of the emission within such beams.The median value of precipitable water vapor (PWV) was 1.4 and 1.5 mm for CO(6-5) and [C ](2-1) -CO(7-6), respectively.The full histograms of the PWV values for the observations are shown in Fig. 3.A summary of the sample and observational setup is shown in Table 1.Due to source visibility and weather constraints, we obtained CO(6-5) data for 7 sources and [C ](2-1) -CO(7-6) for 8 sources.
The data reduction was performed using the / 3 package version 1.1.For each source the data corresponding to a different date were processed separately, before combining them.In this procedure, for every target the noisy edges (∼3%) of the spectra were trimmed.Then, a velocity window4 was chosen to encompass the expected location of the emission line, according to the redshift of the source.First-degree polynomial baselines were computed neglecting the data within that window, and subtracted from the individual scans.All data for each source were then combined into one final spectrum, after visual inspection of individual scans.These spectra cover an average spectral window of ∼6000 km s −1 .
To further improve the root mean square (rms) of the final combined spectra, we applied the following procedure.For each target, we computed the rms for each used subscan in order to reject the nosiest data.We computed the median rms of the whole dataset (i.e., all dates for each source) and removed the data farthest away from this median.New reductions ignoring these subscans were performed following the steps described above, checking if the final rms of the dataset improved.We found that the removal of the noisiest subscans did not improve the final rms, because the decrease in exposure time compensates the improvement in rms, so we decided to keep all the data for the final reduction.The final rms for each tuning is reported in Table 1.
To illustrate the stability of the SEPIA180 instrument at the targeted frequencies, in Fig. A1 of Appendix A we show the rms as a function of the bin size (in km s −1 ), for the final combined spectra of each source, starting from the original resolution and up to 600 km s −1 .At the bottom of each panel is shown by how much the observed rms deviates from the expected value.At a bin size of 300 km s −1 , we found a mild median deviation of 12% and 14% for the CO(6-5) and [C ](2-1) -CO(7-6) observations, respectively.
As last step, we transformed the intensity units of the spectra, originally in temperature (K), to physical flux units.For this purpose, we assumed the conversion factor 38.4 ± 2.8 Jy K −1 , calculated for the SEPIA180 receiver (Belitsky et al. 2018b) 5 .

OBSERVATIONAL RESULTS
The resulting APEX spectra, reduced and converted to flux density units (mJy), are presented in Figs. 4 and 5 for the CO(6-5) and [C ](2-1)-CO(7-6) observations, respectively.For all sources, the left panel spectrum has a bin size between 150 and 200 km s −1 (depending on the depth of the data), and the right panel spectrum has a bin size of 300 km s −1 .These two different bin sizes are shown to highlight the reliability of the detections.In this work we report as detections the lines that fulfill the following conditions: i) have a peak emission at S/N> 2, at bin sizes of 300 km s −1 , ii) are also present at the resolution of 150 km s −1 but with lower significance than at 300 km s −1 , and consistent integrated fluxes, iii) have an integrated S/N> 3.In Appendix B we show that this detection algorithm is reliable, giving basically a zero-rate of false positive identifications in a negative and jack-knife tests.
Importantly, throughout this work we assume that the detected emission is due to the central quasar hosts, unless specified.We checked this assumption by computing the number of expected line-emitter companions for each transition within the APEX observations, down to the 3 limiting luminosities of  =3.13×10 10 , 3.02×10 10 and 2.65×10 10 K km s −1 pc 2 , respectively for CO(6-5), CO(7-6), and [C ](2-1).Specifically, we assumed (i) a cylindrical volume defined by the APEX primary beam and the covered velocity range of ∼ 6000 km s −1 , (ii) the luminosity function of CO(6-5), CO(7-6) and [C ](2-1) emission measured for similar redshifts ( ∼ 2; Decarli et al. 2020), and (iii) a deterministic bias model for the clustering of sources around quasars (e.g., García-Vergara et al. 2017, 2019).In this model, we assume a power-law shape for the clustering, with a fixed slope of  = 1.8, and we use the  ∼ 3 quasar clustering (Shen et al. 2007) and the clustering of Lyman-break galaxies (LBGs) at  = 4 (Ouchi et al. 2004), which are assumed to have similar clustering as CO and [C ] sources.
Following these assumptions, we found that the expected number of companions for the total number of observed fields per line (7, 8, 8) are 0.07, 0.01, 1.39, respectively for CO(6-5), CO(7-6), and [C ](2-1).We caution that the luminosity functions in Decarli et al. (2020) are still associated with large uncertainties, and only upper limits are reported at the bright-end sampled by our observations.We have thus extrapolated their measurements up to brighter fluxes, and therefore our estimations represent upper limits for the number of expected companions.Specifically, the worst case is the current [C ](2-1) luminosity function, which has a relatively flat shape (see Fig. 7 of Decarli et al. 2020).In this case, our extrapolation is flat and therefore represents a clear upper limit considering that the luminosity function is expected to steeply decrease for high luminosities.In summary, the only transition for which we may find a companion is [C ](2-1), with a conservative probability < 17% for each field.Follow-up high-resolution observations with interferometers (e.g., ALMA, NOEMA) are required to verify this assumption and assess whether any of the detected emission comes from companions and/or larger scales.
We also computed the probability that one detected line in one field is actually any CO or [C ] transition from an interloper galaxy at possible lower redshifts given our tunings, by assuming the luminosity functions in Decarli et al. (2019Decarli et al. ( , 2020)), and the comoving volume spanned by our observations.Once again, these estimates represent upper limits as our observations sample the bright-end of the luminosity functions.We found that the probability of observing a contaminant is <0.2% for any CO or [C ] line.Therefore, it is very likely that any detected emission in our observations is associated with the quasar or its environment.1, first column).The white dotted circles represent the APEX/SEPIA180 beam sizes for the CO(6-5) observations (∼32"), while the white dashed circles are the beam sizes for the [C ](2-1) -CO(7-6) observations (∼30").The cyan diamond indicates that the source is radio-loud.), which has an intrinsic uncertainties of 415 km s −1 .See Table 3 for updated systemic redshifts ( mol,sys ) calculated from the lines detected in this work.e Tuning frequency used for the observations.f Total ON-OFF exposure time per source.The total telescope time is roughly double this integration time.g RMS per bin of 300 km s −1 of the final combined spectrum, in antenna temperature  *  units.h Median PWV between all the observed dates for each source.

Emission line measurements
We measured the molecular velocity-integrated emission line fluxes  CO(6−5) ,  [C ] (2−1) and I CO(7−6) , by fitting Gaussian to the detected lines.For those spectra that presented two peaks, a double Gaussian fit was applied.The uncertainties of the measured fluxes include the aforementioned error on the flux conversion factor.The Gaussian fits also provided an estimate of the full width at half maximum (FWHM) of the emission lines and their respective uncertainties.
When an emission line was not detected, we derived 3 upper limits using the rms noise within the same velocity range for the emission lines that were detected for that target, and within a velocity width of 300 km s −1 (expected average line width for quasar hosts, e.g., Weiß et al. 2003;Weiss et al. 2007a;Walter et al. 2011) when no line was detected.In Table 2 we tabulate these fluxes (or upper limits), FWHM and their respective uncertainties for each source.
For the quasars with IDs 39 ([C ](2-1) line) and 43 (CO(6-5) line), it was not possible to obtain a good Gaussian fit at bin sizes of 300 km s −1 (see Figs. 4 and 5).Therefore, their integrated fluxes were first estimated by adding the area covered by each bin of 300 km s −1 contained within the emission line and verified against a Gaussian fit at bin sizes between 150 and 200 km s −1 .Their integrated fluxes computed from both bin sizes are consistent, as per our detection criteria.The FWHMs and fluxes given in Table 2 are those estimated with the Gaussian fit.For the quasar with ID 3 (or J 0525-233), all line properties are listed using bin sizes of 150 km s −1 , because at this resolution we obtained a better fit (see Fig. 5) than at bin sizes of 300 km s −1 .
We started the analysis from the CO(6-5) tuning.From the spectra shown in Fig. 4, we found CO(6-5) emission line detections (represented by a green shaded area) for three sources: quasars with IDs 8, 43, and 50.Two of these detections are closer to the redshift of the Ly nebula (IDs 43 and 50; purple vertical line) than the uncertain systemic redshift computed from C (orange vertical line).They present fluxes of 3.4≤  CO(6−5) <5.1 Jy km s −1 .The other targeted sources, quasars with IDs 13, 18, 29 and 39, do not present any CO(6-5) emission down to the current rms.We therefore report 3 upper limits for these sources (see Table 2).
From the spectra in Fig. 5 GHz (or +3966 km s −1 , integrated S/N = 2.7,  line ∼2.8Jy km s −1 ) and 193.4 GHz (or -957 km s −1 , integrated S/N = 2.1,  line ∼2.5 Jy km s −1 ), respectively.The redshift of these lines is not consistent with the CO(6-5) detections for these sources, but they might be an associated object.We note that sub-millimeter galaxies (SMG) have often being found in the surroundings of quasars (e.g., Silva et al. 2015), however the CO(6-5) transition (detected for IDs 43 and 50) is expected to be stronger in quasars than in SMG (Carilli & Walter 2013).Also, we stress that the detections in the CO(6-5) tuning must come from the quasar, due to the low probability of finding a companion in this transition in our observations (see Section 3).The rest of the sources do not have any detected emission lines at the rms of the current observations.It is important to stress that we considered our detected lines as [C ](2-1) and not CO(7-6) because they are closer to the expected position of [C ](2-1) based on the systemic redshift of the sources (within 1, IDs 3 and 39), or because if they were CO(7-6) emission, CO(6-5) detections would be also expected unless unphysical line ratios are assumed (quasar with ID 29, see next paragraph).We acknowledge that the used systemic redshifts (based on C ) are uncertain (see Section 3.2) and further data are needed to verify our line identifications.
The spectrum of the quasar with ID 3 (or J 0525-233) shows the [C ](2-1) emission line (at 196.8 GHz) and a S/N = 1.4 bump (binning of 300 km s −1 ) at the expected location of CO(7-6) (196.2GHz), which we estimate would contribute ∼ 20% of the [C ](2-1) flux when using the 300 km s −1 binning.In contrast, the quasar with ID 29 (or Q-0115-30)6 presents two emission lines (at 191.9 and 191.4 GHz) consistent with the redshift of its Ly  .The seven APEX CO(6-5) spectral line observations.Each observed spectrum is shown with bin sizes between 150 and 200 km s −1 (left column) and 300 km s −1 right (column).In each panel, the vertical black error bar represents the uncertainty (rms) per bin, while the ID number of each quasar is shown in the top right corner.The orange vertical lines at the velocity zero point represent the systemic redshift of the sources before our APEX observations (from C , see section 3.2), and the shaded regions correspond to its uncertainty (∼415 km s −1 , Arrigoni Battaia et al. 2019a).The purple vertical lines represent the nebulae Ly redshift from MUSE extracted within the APEX beam (see Fig. 2), with their corresponding uncertainty (magenta shaded regions).The green shaded area indicates the detected CO(6-5) emission line, while the dark green curve represents the Gaussian fit applied to this line, when possible.nebula.Both these lines likely correspond to [C ](2-1).To confirm this, we first assumed that the line at frequency 191.4 GHz was CO(7-6) and estimated the observed I CO(7−6) /I CO(6−5) ratio (> 3.73).We then compared this observed ratio to theoretical predictions from large velocity gradient (LVG) models (Goldreich & Kwan 1974) with different physically plausible kinetic temperatures and densities (see Section 4.2 for a detailed explanation of these models and input parameters).This test did not allowed us to find any modelled I CO(7−6) /I CO(6−5) ratio as high as that observed, indicating that the line peak at lower frequencies also corresponds to [C ](2-1) 7 .We stress that this identification needs further confirmation, e.g. by targeting additional transitions from this source. 7All the LVG models explored in this work have I CO(7−6) /I CO(6−5) < 1.34 when considering kinetic temperatures of 20 − 200 K, and densities of 10 3 − 10 5 cm −3 , as frequently used in the literature (e.g., Riechers et al. 2006;Weiss et al. 2007a).These LVG models allow a maximum I CO(7−6) =3.48 Jy km s −1 at a kinetic temperature of 200 K, which represent the most extreme SLED found in the literature (Weiss et al. 2007a).Even in this extreme case, a fraction of 63% of the observed second line peak has to correspond to [C ](2-1).Quasar ID 29, with a total  [C ] (2−1) = 15.7 ± 3.7 Jy km s −1 , is therefore the brightest [C ](2-1) detection in this sample.Furthermore, it is ∼3 times brighter than the lensed Cloverleaf quasar ( [C ] (2−1) ∼ 5.2 Jy km s −1 , Weiß et al. 2003), which is as far as we know the brightest [C ](2-1)  ∼ 3 quasar detected in the literature.We further discuss the nature of this double-peak emission of quasar with ID 29 (or Q-0115-30) in Section 6.1.2.Overall, our observations suggest that [C ](2-1) is stronger than CO(7-6) in these  ∼ 3 quasars.This is in contrast with other observations of the same transitions in different  ∼ 4 − 7 quasars (Riechers et al. 2009;Venemans et al. 2017;Lu et al. 2018;Wang et al. 2019;Yang et al. 2019) and  ∼ 2.5-3 quasars (Weiß et al. 2003;Walter et al. 2011;Schumacher et al. 2012), likely indicating that the excitation conditions in the molecular gas are different.Also, Banerji et al. (2018) observed these transitions in two quasars and their companions at  ∼ 2.5, finding a [C ](2-1) line slightly stronger than CO(7-6) in one quasar and one companion.To our knowledge, the studies mentioned above are the only available observations of these lines in the literature.Therefore, if our observations are confirmed in deeper datasets, current literature could be affected by low number statistics.We further note that [C ](2-1) is stronger than CO(7-6) in high-redshift radio galaxies (∼5 times stronger, e.g., The integrated fluxes, S/N and FWHM have been calculated using a bin size of 300 km s −1 .For the cases with non-detections, 3 upper limits on the fluxes are provided.* Feature added in the Gaussian fit of quasar with ID 3 (or J 0525-233), to not overestimate the integrated flux of the [C ](2-1) line at bin size of 300 km s −1 (see Fig. 5).
† These values have been calculated using a bin size between 150 and 200 km s −1 .‡ These values have been calculated as the sum of the two Gaussian components (see Fig. 5).For each component, the individual values of I CI(2−1) are 6.1±2.9 km s −1 and 9.5±2.3Jy km s −1 , with FWHM CI(2−1) values of 471±224 km s −1 and 471±108 km s −1 , respectively.Gullberg et al. 2016), which can be explained as enhancement of atomic carbon in cosmic ray dominated regions (e.g., Bisbas et al. 2017).

Redshifts from molecular lines
A precise estimate of the systemic redshifts of quasars plays a fundamental role in understanding the physical processes and kinematics of each system.For instance, an accurate systemic redshift would allow us to better constrain the cool gas kinematics mapped on large scales by the Ly emission, and to compare it to cosmological simulations (Arrigoni Battaia et al. 2019a).
The uncertain systemic redshifts for our sample, estimated by Arrigoni Battaia et al. (2019a), are shown in the fifth column of Table 1.These redshifts were determined from the peak of the C line, after correcting from the expected luminosity-dependent blue-shift (Shen et al. 2016), and have an intrinsic uncertainty of 415 km s −1 .This large uncertainty, comparable to outflow/inflow velocities expected in quasar halos, hampers any kinematical study of these systems.However, the molecular emission lines should provide a more robust measure of the quasar's systemic redshift (e.g., Banerji et al. 2017;Bischetti et al. 2021).We therefore derived new systemic redshifts for the objects with detected molecular lines.For the case of ID 29, the redshift was estimated as an average between the centroids of the two peaks.These new redshifts and their uncertainties are listed in Table 3.They have on average an uncertainty of 74.8 km s −1 and a difference of +1045 km s −1 (or Δ = +0.014)with respect to the systemic redshifts from C (see also Section 5.2).Hereafter we will assume these new values as systemic redshifts.

MOLECULAR MASS ESTIMATES
In this section we estimate molecular gas masses M H 2 by using different methods.Specifically, we compute i) carbon masses M C and derive the respective M H 2 by assuming a carbon abundance relative to H 2 , ii) M H 2 for sources that have a measured constraint (upper limit) on the CO(7-6)/CO(6-5) ratio, iii) M H 2 by combining the two previously obtained mass ranges from [C ] and the CO ratio, and iv) M H 2 for sources with no clear constraint on the CO ratio (i.e.sources with non-detections).

Atomic carbon mass
The mass of the atomic carbon can be estimated from the [C ](2-1) line luminosity through the formulation presented in Weiß et al. (2003Weiß et al. ( , 2005)), under the assumption that this [C ] transition is optically thin: where ( ex ) = 1+3 − 1 / ex +5 − 2 / ex corresponds to the C partition function,  ex is the excitation temperature, and  1 = 23.6K and  2 = 62.5 K correspond to the energies above the ground state.We used  ex =30 K as frequently found in high-redshift quasars (e.g., Weiß et al. 2003;Walter et al. 2011).The [C ] line luminosity,  [CI] , can be estimated via (Solomon et al. 1992) where  line is the velocity-integrated line flux in units of Jy km s −1 ,  obs is the observed frequency in units of GHz, and  L corresponds to the luminosity distance in Mpc.The final units of  line are K km s −1 pc 2 .
As noted in Section 3.1, the [C ](2-1) emission line is detected for three sources of our sample: quasars with IDs 3, 29 and 39.Using equation 1, we found atomic carbon masses in the range of ∼4.0×10 7 M -3.1×10 8 M .For the five sources with upper limits on the [C ](2-1) transition (quasars with IDs 7, 8, 13, 43 and 50), we found atomic carbon masses < 8.0×10 7 M .The detected sources show higher values than usually reported for quasars in the literature.Indeed, it is common to find values of the order 10 6 -10 7 M , when the same excitation temperature of 30 K is assumed (e.g., Walter et al. 2011;Venemans et al. 2017;Yang et al. 2019).The value for each source is listed in the third column of Table 4.It is important to note that the assumption of a higher  ex in equation 1, for instance  ex =50 K, implies a  C ∼38 % lower.
The atomic carbon mass can be used to determine the molecular gas mass using the atomic carbon abundance relative to H 2 : In this work we assumed  [C ]/X[H 2 ]= (8.4 ± 3.5) ×10 −5 as usually done in the literature for high-redshift quasars (e.g., Walter et al. 2011;Venemans et al. 2017).We note that this value is higher than what has been recently found for  ∼ 1 − 3 main sequence galaxies (∼1.7×10 −5 , Valentino et al. 2018;Boogaard et al. 2020).
Using this method, we constrained molecular gas masses for quasars with IDs 3, 29 and 39, which span from ∼8.0×10 10 M to ∼6.1×10 11 M .For quasars with IDs 7, 8, 13, 43 and 50, we found molecular gas masses < 1.5×10 11 M .Most of these molecular gas masses are higher than the typical range (∼109 -10 11 M ) found in the literature for other high- quasars (e.g., Weiß et al. 2003;Walter et al. 2011;Anh et al. 2013;Hill et al. 2019).We list all these values in the fourth column of Table 4.
Importantly, we are assuming that all the [C ](2-1) emission comes from the quasar hosts.If part of this emission comes from larger areas, the  ex might break into different values (e.g., lower/higher on larger/smaller scales), affecting the molecular masses estimated in this work.To constraint  ex around these quasars, resolved [C ](2-1) and [C ](1-0) observations are needed (e.g., ALMA, NOEMA).

LVG Models and CO constraint
The molecular gas mass can also be derived from CO following the equation with  being the CO luminosity-to-gas mass conversion factor, and  CO(1−0) is the luminosity of the CO(1-0) that can be estimated from equation 2. In this work, we assume a value of  = 0.8 M (K km s −1 pc 2 ) −1 (Downes & Solomon 1998), which has been estimated for local ultra-luminous infrared galaxies and has been typically adopted to calculate molecular masses in high-redshift quasars (e.g., Wang et al. 2010;Venemans et al. 2017).
To estimate the molecular gas masses from the CO(6-5) and CO(7-6) transitions, we need to assume a CO spectral line energy distribution (hereafter SLED) to find the CO(1-0) line intensity.We did not find in the literature any quasar characterized by a SLED 8 that agrees with the observed constraints on the CO ratios.For this reason, we modelled the CO SLED using the large velocity gradient (hereafter LVG) method, which has been widely applied to high- 8 There are only a handful of z∼3 quasars with a well-characterized SLED (see Fig. 8).
quasars by several authors (e.g., Weiß et al. 2007b;Riechers et al. 2009;Schumacher et al. 2012).These models include a velocity gradient (that indicates the change in the line of sight velocity in the turbulent medium) that is considerably larger than local thermal velocities of the gas, leading to photons being able to escape due to the different velocities along the cloud, following a photon escape probability.This method finds the populations of the molecular energy levels excited by collisions with H 2 (main collision partner for CO) for certain given parameters as CO abundance, kinetic temperature (T kin ), H 2 density (n H 2 ) and velocity gradient (dv/dr).It is then possible to identify the physical parameters that best describe the conditions of the gas through the comparison of the model predictions to the observed line ratios (e.g., Carilli & Walter 2013).
In this work, we have used the radiative transfer code 9 (van der Tak et al. 2007), considering a spherical and single component LVG model.For the calculations, we adopted an H2 orthoto-para ratio of 3 and collision rates from Yang et al. (2010).We set the following input parameters and explore a grid of different models in which we vary some of these parameters: T kin equal to 30 K, n H 2 range of 10 3 -10 5 cm −3 (typical for quasar host galaxies, Carilli & Walter 2013), background temperature of ∼11 K (cosmic microwave background at  ∼3), and a column density of the molecular gas given by: where Δ turb corresponds to the turbulence line width fixed here at 100 km s −1 , and [CO]/dv/dr = 1×10 5 pc (km s −1 ) −1 , following procedures commonly adopted in the literature (e.g., Weiß et al. 2007b;Schumacher et al. 2012).LVG models are intrinsically degenerate in the parameters T kin and n H 2 , meaning that different combinations of these parameters can give the same SLED.We focused only on models with T kin = 30 K because the value of the kinetic temperature is expected to be comparable to the excitation temperature of neutral carbon (Israel & Baas 2002), which we assumed to be ∼30 K (see Section 4.1).In total, 31 different CO ladders were modelled, which are shown in Fig. 6.In this figure, the CO line flux normalized to the CO(1-0) line is shown as a function of rotational quantum number J, and the colour bar represents the different values of n H 2 .The Note: Upper limits of masses provided in each column are due to non-detections of [C ] and/or CO (see Section 4). a Atomic carbon mass assuming an excitation temperature of 30 K (see Section 4.1).b Molecular gas mass derived from the atomic carbon mass (see Section 4.1, equation 3).c Molecular gas mass derived from CO (see Section 4.2).d Molecular gas mass derived from applying both CO and [C ] constraints (see Section 4.3).e Inclination-dependent dynamical masses (see Section 6.1.1).The reported uncertainties include only the errors on the FWHM values.* In the case the detected line for ID 3 is CO(7-6), we estimate 1.2 <  H 2 /[10 11 M ] < 4.5 from the final CO SLEDs (see Fig. 8).† In the case the detected line for ID 39 is CO(7-6), we estimate 0.2 <  H 2 /[10 11 M ] < 0.6 from the final CO SLEDs (see Fig. 8).‡ This value correspond to the dynamical mass for the disk scenario (see Section 6.1.2).This value correspond to the dynamical mass for the merger scenario (see Section 6.1.2). peak of the modelled ladders varies from J∼3 to J∼8, for n H 2 =10 3 cm −3 and 10 5 cm −3 , respectively.We note that according to current observations, quasars are expected to show the SLED peak between J∼6 and J∼8 (e.g., Riechers et al. 2009Riechers et al. , 2011;;Carilli & Walter 2013;Banerji et al. 2018;Wang et al. 2019;Bischetti et al. 2021).Importantly, the emitted CO flux is proportional to the source solid angle (e.g., Weiß et al. 2007b).As our observations are not able to constrain the size of the emitting source, we focus on line ratios in this paper.These LVG models will therefore represent only solutions for the component of highly excited gas, which likely coexist with a less excited component emitting the [C ](2-1) emission that we observe in some of our targets.Indeed, it has been shown that emission from species with very different critical densities likely originate from gas at different densities (e.g., Harrington et al. 2021), and also on different scales (e.g., Emonts et al. 2014;Casey et al. 2018;Spingola et al. 2020).
We first find models that reproduce the observed upper limits on the CO(7-6)/CO(6-5) line ratios.In our dataset there are only three sources for which it was possible to estimate an observed  CO(7−6) / CO(6−5) upper limit: the quasars with IDs 8, 43 and 50 (see Table 2).For these sources, we selected LVG models constrained by their observed  CO(7−6) / CO(6−5) ratio and found densities in the range n H 2 ∼10 3 -10 4.4 cm −3 .These models constrained the peak of the SLED to be at J∼3-7, and the molecular gas masses to be (0.5 -64.1)×10 11 M when considering the predicted  CO(1−0) from the SLEDs.These results are shown in the fifth column of Table 4.
It is noteworthy that using a higher T kin in our LVG models, for instance T kin = 50 K (found in some high- quasars, e.g., Weiss et al. 2007a;Riechers et al. 2011), implies an upper limit in the molecular masses ∼32% lower, i.e. this difference does not alter our results significantly.In contrast, using a lower T kin , for example T kin = 20 K, implies non-physical values for the upper limit in the molecular masses of up to ∼10 13 M .

Molecular mass constraints from a joint CO and [C ] analysis
After applying the CO constraint explained above, we set another condition based on the molecular gas masses already estimated from the [C ](2-1) transition (see Section 4.1).As  is fixed, this condition is identical to impose a constraint on the CO(1-0) luminosity.We caution that this approach may introduce a bias in our calculation, which depends on the goodness of the assumed parameters to model the molecular mass from [C ], and on  itself with respect to the physical conditions in each individual source.Observations of the CO(1-0) transition for these sources are definitely needed to confirm this methodology.
For each source, this step excluded some of the LVG models at the lowest densities selected in Section 4.2.In this way, we obtained the final CO SLEDs from the union of the CO and [C ] constraints.We show these final constrained CO SLEDs in Fig. 7 for quasars with IDs 8, 43 and 50.Their ranges of molecular masses using these joint constraints are (0.7-1.5)×10 11 M , (0.5-0.9)×10 11 M and (0.5-1.5)×10 11 M , respectively.These values are also listed in the second to last column of Table 4.
Fig. 8 shows the CO SLEDs obtained in this work (discontinuous lines) in comparison to three  ∼3 quasars with well sampled SLEDs (grey continuous lines), Cloverleaf ( = 2.6, Weiss et al. 2007a;Bradford et al. 2009), MG 0751+2716 ( = 3.2, Weiss et al. 2007a) and PSS1409 ( = 2.6, Weiss et al. 2007a).From this figure, we see that none of the SLEDs obtained in this work are similar in shape to those three obtained previously, indicating that our quasars have different physical properties.Our results yield SLEDs with J peak somewhat lower than the expected range for quasars (J∼6-8), varying from J∼5 to J∼7.
The six sources with IDs 3, 7, 13, 18, 29 and 39, have only upper limits for the targeted CO transitions.Our observations are therefore not able to put any constraint on the SLED of these sources.To compute the upper limits for the CO(1-0) luminosity, we assumed that the models selected for the previous three sources with IDs 8, 43 and 50, apply to the full observed sample.
We first selected the CO SLEDs found for the quasars with IDs 8, 43 and 50 (shown in Fig. 8).Depending on which emission line upper limit flux was measured for each of the six sources mentioned above, we used the CO(1-0)/CO(6-5) and/or CO(1-0)/CO(7-6) ratios.Then, we multiplied these CO ratios by the observed upper limit flux for each of them.After this, we estimated upper limits for the molecular gas masses of the six quasars with IDs 3, 7, 13, 18, 29 and 39.For the sources that have flux upper limits for both CO lines (quasars with IDs 13, 29 and 39), two upper limits in molecular gas mass were derived and the higher value was selected.
To summarize: the resulting molecular masses of our sample, obtained from applying only the CO constraints, are in the range of (0.5 -64.1)×10 11 M for the sources with IDs 8, 43 and 50, and are <6.9×10 11M for the objects with IDs 3, 7, 13, 18, 29 and 39 10 .We note that the range mentioned includes masses of the order of 10 12 M , which are implausible given the expected dark-matter halo mass for these quasars ( DM∼ 10 12.5 M , e.g., White et al. 2012).These values are also much higher than the values found in the literature for other high- quasars (∼10 9 -10 11 M , e.g., Weiß et al. 2003;Walter et al. 2011;Anh et al. 2013;Hill et al. 2019).These high masses are caused by models with densities of n H 2 ∼10 3 -10 3.8 cm −3 .This disagreement with the literature values suggests that such low densities are not plausible in explaining our current constraints on the emission from the high-J CO transitions.We note that the molecular gas masses obtained from applying the CO and [C ] constraints jointly, remove such models.
We emphasize that the final selected densities (n H 2 ∼10 3.8 -10 We note that the latter upper limit is based on the final CO SLEDs obtained for the sources with IDs 8, 43, and 50 (Figure 8).  10 4.4 cm −3 , see Fig. 8) represent those values needed by the high-J CO lines to match the molecular masses estimated from [C ](2-1).These densities should not be regarded as representative of those emitting the [C ](2-1) line, unless this emission is more extended than the high-J CO transitions. 11 The molecular masses that are finally derived after applying a more restricted set of models are in the range of 0.4 -6.9×10 11 M for quasars with IDs 3, 8, 29, 39, 43 and 50, and <1.1×10 11 M for quasars with IDs 7, 13 and 18.The brightest of our detected sources still has an estimated molecular gas mass that is higher than typically found in the literature.We will discuss this further in the next sections.
We finally check the molecular gas mass estimate for ID 3 and ID 39, for which we identify [C ](2-1) based on the uncertain C systemic redshift (section 3.1).In the case the detected lines are CO(7-6) instead of [C ](2-1), the final models shown in Figure 8 would imply molecular gas masses of (1.2 -4.5)×10 11 M and (0.2 -0.6)×10 11 M , respectively.These values would still confirm our findings, with ID 3 having a large molecular mass.

COMPARISON WITH THE LARGE-SCALE LY𝛼 EMISSION
To investigate any relation between the molecular gas phase and the cool halo gas, in this section we compare the APEX observations with the Ly line properties of the associated large-scale nebulae discovered with MUSE (Arrigoni Battaia et al. 2019a).First, we focus on finding any trend in the properties of each phase when 11 For completeness, we report here the predictions of the [C ](2-1)/CO(6-5) ratio using our 10 finally selected LVG models.These models with T kin = 30 K are characterized by [C ](2-1)/CO(6-5)<1.33.In the case of ID 29, this value is lower than the observed [C ](2-1)/CO(6-5)>4.36.This apparent discrepancy could be solved by invoking a [C ](2-1) solid angle > 3× the high-J CO solid angle (which corresponds to > 1.7× the high-J CO source radius) and/or with multi-component models (e.g., Weiß et al. 2003;Harrington et al. 2021).
compared with the other.Secondly, we study the line emission locations and profiles.

Molecular lines versus Ly𝛼 line properties
The FWHM of molecular lines from high-redshift galaxies are frequently used as a dynamical mass estimator (assuming a rotating disk geometry, e.g., Walter et al. 2003;Narayanan et al. 2009;Wang et al. 2013), and hence their stellar mass and therefore their halo mass can be determined (see Section 6.From this figure, we find that there is no clear correlation between FWHM of the molecular lines and Ly properties.Also, we note that all the quasars show similar values of FWHM (in the range of 329 -943 km s −1 , average FWHM 647±129 km s −1 ), considering the uncertainties.In the literature, the molecular linewidths (CO and [C ] lines) found for high- quasars have values between ∼150 -450 km s −1 (e.g., Weiß et al. 2003Weiß et al. , 2007b;;Walter et al. 2011;Schumacher et al. 2012;Venemans et al. 2017).The values found for our quasars are larger on average.For example, quasars with IDs 29 and ID 50 have a FWHM of 943±249 km s −1 ([C ](2-1) emission) and 707±100 km s −1 (CO(6-5) emission) respectively, suggesting that they have a different kinematics and/or physical properties compared to the quasars studied in the literature.
As molecular gas is expected to form from the cooling of 10 4 K gas (e.g., Dalgarno & McCray 1972), one may naively expect to observe the largest molecular reservoirs in sources with the largest and brightest Ly nebulae.As a next step, we focus on the molecular gas masses obtained from the joint constraints from CO and [C ] (see Section 4.3), and compare them to L Ly , SB Ly and the area of the large-scale nebulae.
The results of this analysis are presented in Figure 10, where the colours represent different sources, the vertical dashed lines are ranges and the arrows upper limits of the estimated molecular gas masses.The horizontal black lines represent the expected molecular mass for the quasar host galaxies on the main sequence (MS) of star-forming galaxies (log( H 2 [ ]) = 10.98 +0.32  −0.51 ), on 4×MS (log( H 2 [ ]) =11.29) and 10×MS (log( H 2 [ ]) =11.51).These values were computed under the following assumptions: i) a quasar halo mass of log( halo [ ]) = 12.68 0.81 −0.67 (for z∼3 quasars), estimated in the study of Kim & Croft (2008) from Ly forest statistics (this value encompasses estimates from quasar clustering, e.g., Shen et al. 2007;White et al. 2012;Timlin et al. 2018) mass of the objects in our sample, and iii) a molecular gas fraction expected for objects on MS, 4×MS and 10×MS of star-forming galaxies, as defined in the empirical relation of Tacconi et al. (2018).The grey shaded region represents the large uncertainties in the calculation for 1×MS.We note that quasars hosts are estimated to have star formation rates higher or similar to MS galaxies at the same redshift (e.g., Zhang et al. 2016;Shangguan et al. 2020;Jarvis et al. 2020;Circosta et al. 2021).The histogram shows the  H 2 distribution for our sample (red) and for a sample of  ∼2.5 -3 quasars (blue) extracted from the literature as explained in the figure caption.The molecular masses reported for these quasars are overall below the computed predictions (see histogram in Fig. 10), possibly indicating some gas depletion with respect to similarly massive star-forming galaxies.
It is clear that the molecular gas masses of the quasars with APEX detections are in agreement with the expected values for MS galaxies, with the exception of quasars with IDs 3 and 29, which have masses well above the MS.Intriguingly, these two objects characterised by the highest molecular gas masses are associated with the Ly nebulae with the lowest values of L Ly and SB Ly .The third panel of Fig. 10 shows that the most massive molecular reservoirs are associated with some of the smallest nebulae.They also have the highest FWHM [C ](2-1) emission lines (see Fig. 9).We will discuss this result in more detail in Section 6.3.

Molecular tracers vs Ly𝛼 line profiles
Now that we have a more robust estimate of the systemic redshift  mol,sys for some of our sources, we can compare its location with the redshift of the discovered Ly nebulae, assessing if the Ly emission spans similar velocity ranges with respect to the molecular phase.For this purpose, we compare the molecular and Ly line profiles extracted within the same aperture (the APEX beam; see Fig. 2), and also compare those with the quasar spectra.
Figures 11 and 12 show the normalized APEX CO(6-5) and [C ](2-1) detections (blue) in comparison to their normalized Ly nebula spectrum (upper row) and quasar spectrum (bottom row).The quasar spectra have been extracted within circular apertures of 1.5 arcsec radius.The vertical red dashed line represents the current systemic redshift estimated from the molecular lines and the grey shaded area corresponds to its uncertainty (see Section 3.2).
From both figures, we note that some of the peaks of the Ly nebulae present a significative shift from the current systemic redshift.For each object, we listed these velocity shifts Δ LyNeb−mol between the quasar systemics  mol,sys and the Ly of the nebulosities in the last column of Table 3.The values of Δ LyNeb−mol are in the range of -423 to 1236 km s −1 , with the quasars with ID 8 (or UM683) and ID 43 (or CTSH22.05)having the bluest and the reddest shift, respectively.The velocity shifts do not show any trend with respect to the Ly physical properties explored in this paper.Considering that the Ly photons experiment changes in frequency due to the scattering processes, these shifts can be an indication of bulk inflows (blueshift) or outflows (redshift) (e.g., Prescott et al. 2015;Dĳkstra 2017).
Figure 13 shows the velocity shifts between  sys and  mol,sys (blue), and redshift of the Ly nebulae and  mol,sys (orange).The now obtained systemic redshifts from molecular tracers are found to be, in most of the cases, more consistent with the Ly nebular redshifts, showing an average shift of 176±39 km s −1 (red vertical line).We note that this shift is larger than the average shift found at  ∼ 6, 69±36 km s −1 (Farina et al. 2019), obtained by comparing [C ] redshifts and the nebular redshifts for nine sources (green vertical line).The small number statistics hampers any conclusion from this comparison.For example, if we remove from our sample the radio-loud object, ID 3, we would get a smaller average shift of 46±35 km s −1 , which is consistent with the work by Farina et al. (2019).
This analysis shows that the systemic redshifts obtained from C are not reliable in 4 out of 6 cases (∼ 70%).This poor reliability was already indicated by the large velocity shifts of the extended nebular Ly emission and the C redshifts in surveys targeting  ∼ 2 − 3 quasars, reporting values as high as ∼ 6000 km s −1 (Arrigoni Battaia et al. 2019a;Cai et al. 2019).The same works found that the peak of the nebular emission has smaller velocity shifts with respect to the observed peak of the Ly emission of the quasar.Figures 11 and 12 indicate that for 3 out of 6 quasars (ID 3, 29, and 50), this was the case because of that peak being at the real systemic.We discuss the most significant shifts between  mol,sys and the redshift of the Ly nebulae in Section 6.2.4).The horizontal lines show the expected location of massive star-forming main sequence MS, 4×MS and 10×MS galaxies using the gas fractions in Tacconi et al. (2018) and the halo mass-stellar mass relation in Moster et al. (2018).The grey area is the uncertainty on this expectation, for 1×MS.The histogram show the M H2 distribution for our sample considering the average values of the mass ranges (red), and for a sample of  ∼2.5 -3 quasars (blue) extracted from the literature (Barvainis et al. 2002;Weiß et al. 2003;Beelen et al. 2004;Walter et al. 2011;Schumacher et al. 2012;Hill et al. 2019;Bischetti et al. 2021) 2).The APEX spectra are shown with a bin size of 300 km s −1 .Each quasar ID number is shown in the top right corner of each panel.The red vertical lines with grey shaded area correspond to the systemic redshift estimated in this work with its uncertainty.The dark green curves and green shaded areas follow the same notation as in Figure 4. To guide the eye, the purple curves represent a Gaussian fit applied to the Ly emission lines.Bottom panels: same as above, but comparing with the quasar spectra extracted within a circular aperture of 1.5 arcsec radius.

DISCUSSION
In this section, we discuss different scenarios regarding the quasar host properties found in this work, and we try to interpret the link between the molecular gas content and the large scale Ly emission.Future studies are essential to confirm or reject what we propose in the following sections.

Dynamical masses
For high- galaxies, it is possible to derive the dynamical mass of the system using the formula commonly used in the literature (e.g., Walter et al. 2003;Wang et al. 2013;Venemans et al. 2016), assuming rotational support: where  is the disk radius from the observed line and  circ is the maximum circular velocity of the gas disk, which is calculated as  circ = 0.75 × FWHM/sin(), i.e. it depends on the FWHM of the observed molecular line and inclination  of the galaxy (with respect to the plane of the sky).From our molecular data it is not possible to estimate the parameters  and , for this reason we proceed with a set of assumptions.Arrigoni Battaia et al. in preparation studied the quasar with ID 13 (or PK-1017+109) through ALMA observations, and found a radius of ∼3 kpc for this object.This is on the high side of the typical disk radii range found in literature (from molecular gas) for high-z quasar hosts (∼0.5 -2.5 kpc, e.g., Carilli et al. 2003;Walter et al. 2004;Aravena et al. 2008;Riechers et al. 2009;Polletta et al. 2011;Schumacher et al. 2012;Molina et al. 2021;Stacey et  2021).Here we adopt  = 3 kpc, which is in agreement with our ALMA observations of one of these systems.
Applying the above to five sources in our sample with FWHM estimations, excluding quasar ID 29 that is discussed in Section 6.1.2),we find inclination-dependent dynamical masses in the range of  dyn sin 2 () ∼ (0.4 − 1.9) × 10 11 M for quasars with IDs 3, 8, 39, 43 and 50.The quasar with ID 39 (or SDSS J0100+2105) is the least massive and the quasar with ID 50 (or SDSS J0819+0823) the most massive.The specific values for these dynamical masses are tabulated in the last colum of Table 4, with uncertainties considering only the errors from their FWHMs.Even though these dynamical masses are highly uncertain, we can get a rough idea on the inclination.Indeed, assuming that the amount of dark matter is negligible on host scales,  dyn ≈  * +  H 2 , where  * is the stellar mass.Assuming that the stellar mass is at least as high as the highest allowed molecular mass from our calculations, we find that two sources, IDs 3 and 39, require relatively small angles (∼ 20 • ), while the other four sources, IDs 8, 29, 43 and 50, require an inclination (∼ 50 • − 60 • ) close to the mean expected value of sin() = 0.79 ( = 57.3• , Law et al. 2009).
One caveat is that the estimates are subject to the different molecular lines used.We computed that the CO lines reported in the literature for  ∼ 3 quasars are on average larger than [C ] lines by 112±43 km s −1 (e.g., Weiß et al. 2003;Schumacher et al. 2012;Banerji et al. 2018;Yang et al. 2019).As these transitions have different critical densities, this could be an indication that the lines trace different gas components, which can be the case if [C ] and the high- CO transitions (as CO(6-5)) traces gas also on different scales, with [C ] extending to outer scales (Schumacher et al. 2012).
In Section 3.1, we showed that the quasar with ID 29 (or Q-0115-30) has a [C ](2-1) line emission consisting of two Gaussian components separated by ∼800 km s −1 .While this line identification needs further confirmation, in this section we discuss their possible origins.This double-peaked feature is unique within our sample and could be produced by a rotating disk or a merger between two molecular gas reservoirs (e.g., Neri et al. 2003;Narayanan et al. 2006;Greve & Sommer-Larsen 2008;Polletta et al. 2011).In Figure 14 we show this double-peaked line emission at the highest resolution possible (bin size of 80 km s −1 ) for which there is still a decent peak S/N of ∼2.5 and ∼3.5.
We now calculate the dynamical masses, under the assumption that the system is virialized.For the case of the gas distributed in a rotating disk, we estimate the dynamical mass  disk dyn as explained in Section 6.1.1,but considering the velocity difference between the two [C ](2-1) peaks Δ, following the formalism of Neri et al. (2003).For the disk radius, we assume a value of  = 3 kpc.We obtain  disk dyn  2 () ∼ 4.4 × 10 11 M .This quasar has the highest dynamical mass of the sample, which is somewhat expected given its larger molecular mass.A larger molecular mass and dynamical mass with respect to the other quasars could also imply a larger dark-matter halo mass.
The combination of an almost edge-on disk and a larger halo mass, could explain a low level of extended Ly emission as observed for this system.In a unification scenario for AGN (e.g., Urry & Padovani 1995) radiation from the quasar escapes through two ionization cones determined by obscuration from the dust distribution on small scales.Since we observe this object as a quasar, our line-of-sight is inside its ionizing cones.A large fraction of the quasar emission will pass through its massive host molecular disk, and will be absorbed before reaching the CGM.In other words, the obscuration from such a misaligned host galaxy decreases the budget of ionizing photons able to boost the Ly emission on CGM scales, resulting in a dimmer and less extended Ly nebula.Any Ly emission boosted on the scales of the host galaxy likely experience a large number of scatterings and final absorption due to dust, expected to be present in large quantities in such a large molecular reservoir (e.g., Venemans et al. 2017; see Section 6.3 for a rough indication of the dust mass).
On the other hand, a large halo mass implies a smaller fraction of cool gas able to survive the shock-heating process and penetrate the halo down to the host galaxy (e.g., Dekel & Birnboim 2006).In turn, warmer halo and a smaller cool CGM reservoir would result in a lower level of Ly emission.We therefore conclude that in this scenario, the mass of the halo together with the geometry of the system and Ly radiative transfer effects could determine the small extent and low SB level of the Ly nebula surrounding this quasar.Radiative transfer calculations possibly coupled to cosmological hydrodinamical simulations with dust implementations are needed to confirm this interpretation.
For the merger scenario, we can estimate the dynamical mass  merger dyn following Greve & Sommer-Larsen (2008): where FWHM is the full width half maximum of the observed [C ](2-1) emission lines, and  is the half projected distance between the two merging objects.We used the same  value assumed in the disk scenario (to obtain a lower limit in  merger dyn ) 12 , and FWHM = 943 ± 249 kms −1 (see Table 2), which was calculated as the sum of each single Gaussian fit corresponding to each [C ](2-1) observed line.We obtain a mass of  merger dyn = 1.2 × 10 12 M .Assuming that the amount of dark matter is negligible on host scales, the stellar mass is given by 12 In this scenario the companion object could be located anywhere inside the APEX beam, with the 800 km s −1 velocity shift representing a large peculiar velocity or a distance of at maximum ∼ 1.6 Mpc within the Hubble flow.Given the expected virial radius of such massive objects, the latter scenario will represent a merger in its very early phases or two objects still separated in the Hubble flow.We also estimated the individual stellar masses of each component of the merging system.The molecular gas masses for the two different objects computed from the atomic carbon masses are  ∼ 3.6 × 10 11 M .The current MUSE data are not able to verify whether such a massive companion object exists, but this field is characterized by nine faint continuum objects with unknown redshifts within the beam of our APEX observations (Figure 15).These objects do not have any emission lines or absorption features that allow us to verify their redshifts.In particular, we note the presence of one continuum source at ∼ 3 (or ∼ 23 projected kpc) from the quasar.In summary, given the brightness of the quasar and the seeing-limited MUSE observations, we are not able to assess whether there is an ongoing merger on very small scales.A strongly dust-obscured and massive companion could be missed by our optical observations (e.g., Omont et al. 1996).Finally.we note that within the APEX beam we also find three low- interlopers (cyan circles in Figure 15), whose redshifts can be clearly determined ( ∼ 0.52, 1.05, 0.07).These redshifts ensure that the molecular emission from our targeted object is not contaminated by these sources.
Our consideration of ID 29 shows that its dynamical mass is very large, independent of the origin of the double-peaked [C ](2-1) line emission.Observing this system with an interferometer (e.g., ALMA, NOEMA) would allow us to map the [C ](2-1) emission and ascertain whether ID 29 is associated with an exceptionally massive molecular disk or is merging with a similarly massive companion.
Such observations, by probing the ∼ 1 − 2 mm continuum, would be in turn able to verify the reason why the Ly nebula around ID 29 is dimmer and less extended than similarly bright quasars.

Ly𝛼 nebulae kinematics with respect to molecular gas
The use of Ly emission as a tracer of gas kinematics is a complex task due to its resonant nature (e.g., Neufeld 1990).Ly photons are expected to interact several times with Hydrogen gas before escaping most astrophysical systems (e.g., Dĳkstra 2017).During this process, scattered photons can experience both large changes in their frequency and a large displacement in space, possibly washing out any information on gas kinematics.Also, the larger the number of scatterings in a medium, the higher the probability for a Ly photon to be absorbed by dust.If the dust distribution is not homogeneous, the Ly line shape could be affected, possibly hiding information on the kinematics of the system.
In Section 5.2 we compared the molecular emission lines observed with APEX, with the Ly emission obtained with MUSE and integrated within the APEX beam.The main result is that we find velocity shifts between the Ly emission and the systemic redshift  sys,mol in the range -423 to 1236 km s −1 .The presence of a velocity shift between this Ly emission line integrated on halos scales and the systemic redshift of the quasar host galaxy, could be associated with a variety of physical processes occurring in the halo, e.g.substructures or gas infalling onto the central quasar, large-scale outflows, rotating structures or projection effects along the line-ofsight (see Arrigoni Battaia et al. 2018a, and references therein).In the following discussion we focus on the possibility that these shifts represent bulk inflows or outflows.
In particular, as indicated by Ly radiative transfer modelling (e.g, Verhamme et al. 2006;Laursen et al. 2009), photons scattering through outflowing or inflowing gas should appear redshifted or blueshifted, respectively, from the systemic redshift.The quasar with ID 29 has a shift consistent with  sys,mol .We cannot therefore draw any firm conclusion for this source with the current dataset.It could be that Ly radiative transfer effects and a balanced interplay of outflows and inflows result in the absence of a strong line shift or that very fast outflows or inflows bring the gas out of resonance, allowing the observation of Ly at  mol,sys .For the remaining four sources, IDs 3, 8, 39, 43 and 50, we find significant offsets.In particular, IDs 3 and 43 present the largest shifts and also show the largest values of FWHM of the Ly emission lines in the sample (with FWHM Ly > 1000 km s −1 ), which is typically characterized by relatively quiescent kinematics (median FWHM Ly = 915 km s −1 ) 13 .In the following, we review these four sources.
The quasars with IDs 8, 39 and 50 show negative shifts of -423 ± 76, -152 ± 64 and -377 ± 44 km s −1 , respectively, which could indicate an overall inflow signature from CGM scales onto the central quasar.The quasar with ID 39 (or SDSS J0100+2105) is surrounded by one of the less extended Ly nebulae in our sample (see Fig. 10; Ly extends out to ∼ 50 projected kpc).The observed blueshift for this source is lower than the inflow velocities ( in ≈ 0.9 vir ∼ 320 − 460 km s −1 ; Goerdt & Ceverino 2015) expected 13 The FWHM Ly values reported in this work are from Gaussian fits of the Ly emission integrated within the APEX beam.They therefore differ from the 2D first moment analysis on resolved maps presented in Arrigoni Battaia et al. (2019a).The FWHM Ly values are not corrected for the instrument resolution.
for the cool gas in quasar host halos (10 12 − 10 13 M , e.g., White et al. 2012;Timlin et al. 2018 and references therein).This could be due to projection effects or to the presence of other violent motions (e.g.winds, turbulences) that could wash out the inflow signature.
The observed blueshift for the quasars with ID 8 (or UM683) and ID 50 (or SDSS J0819+0823) is comparable with the inflow velocities expected for the cool gas in quasar host halos.The quasar with ID 8 has an average molecular reservoir and a Ly nebula with intermediate surface brightness, extending out to ∼ 145 kpc (Arrigoni Battaia et al. 2019a.)Interestingly, ID 50 is surrounded by one of the brightest and more extended (in area) Ly nebulae in the QSO museum sample (see Fig. 10; Ly extends out to ∼ 130 projected kpc), with a clearly asymmetric morphology.The widths of their integrated Ly lines, FWHM Ly = 915 ± 46 km s −1 and FWHM Ly = 967 ± 38 km s −1 (for IDs 8 and 50, respectively), may therefore be explained by gravitational motions in such massive halos (FWHM DM 1D ∼ 580 − 870 km s −1 ; Arrigoni Battaia et al. 2019a) and by the presence of large turbulence within the cool gas reservoir.
The quasar with ID 3 is radio-loud, characterized by a large reservoir of molecular gas and has one of the Ly nebulae with the lowest surface brightness, extending out to a maximum of ∼ 70 projected kpc (Arrigoni Battaia et al. 2019a).Its integrated Ly emission shows the most active kinematics of the sample (FWHM Ly = 1542 ± 196 km s −1 ), which together with the observed large Ly shift likely indicates that small-scale radio feedback is affecting the surrounding gas distribution in this system.The quasar with ID 43 has an average molecular reservoir and a Ly nebula with intermediate surface brightness, extending out to ∼ 118 kpc (Arrigoni Battaia et al. 2019a).Its radio nature is not clear because there are no radio observations of this target in the literature.Its relatively active Ly kinematics (FWHM Ly = 1068 ± 142km s −1 ), together with the large observed positive shift of the Ly line with respect to CO(6-5) and the asymmetric morphology of the nebula with a northern bright region, might all be hints of outflowing gas in this system.

Relation between large molecular reservoirs and Ly𝛼 nebulae
In section 5.1, we found a tentative trend indicating that the two quasars with the lowest surface brightness for their Ly nebulae are associated with the highest molecular gas masses.As the molecular gas content is also closely linked to the amount of dust 14 , this 14 To have an idea of the dust masses of our sample, one could assume a gas-to-dust ratio of 70 (e.g., Riechers et al. 2013;Venemans et al. 2017) and considered that the ∼75% of the total gas mass corresponds to molecular gas (e.g., Wang et al. 2016).Under these very simple assumptions, our quasars should have a dust mass in the range of  dust = 0.8 -12.9×10 9 M for the sources with IDs 3, 8, 29, 39, 43 and 50, and are <2.1×10 9 M for the objects could imply that the radiation of these quasars is more obscured by their host galaxies, therefore reducing the reprocessed Ly emission from the surrounding gas distribution.Two effects could then be in play in the radiative transfer of these systems: first, scattered Ly photons on small scales (tens of kpc) encounter a very high column density and dusty medium resulting in absorption of the emission, and second, most ionizing photons are not able to escape the massive host galaxy towards the CGM, limiting the detectability of cool halo gas in emission.On top of this, the halo mass could play an important role in determining the reservoir of cool gas on CGM scales, with more massive halos having a smaller fraction of their gas in this phase.Taking also this aspect into consideration, there are at least two possible interpretations for this tentative trend: 1) IDs 3 and 29 (which present the most massive molecular reservoirs together with dimmer nebulae) occupy similar massive dark matter halos as the rest of quasars in the sample, but their molecular gas fraction (and consequently dust) are higher, resulting in smaller escape of ionizing and Ly photons, 2) IDs 3 and 29 live in more massive dark matter halos, which are warmer and characterized by a smaller cool gas fraction, hence lower levels of Ly.
Unfortunately, we cannot unambiguously infer the halo masses with the available observations, nor we can firmly quantify obscuration in these systems.Future analysis based on high-resolution interferometric observations and spectra from other CO transitions (e.g., with ALMA, NOEMA) are needed to better constrain the nature of these sources.Such follow-up studies together with constraints on the far-infrared continuum will be able to accurately determine the dust and molecular content of these systems, allowing a better description of the radiative transfer around these quasars.

SUMMARY
With the purpose of gaining a new leverage on the gas cycle around  ∼ 3 quasars, we exploit the frequency window opened-up by SEPIA180 at APEX, targeting the CO(6-5), CO(7-6) and [C ](2-1) transitions for a sample of 9  ∼ 3 quasars already observed with MUSE/VLT.Specifically, these 9 targets (quasars with IDs 3,8,13,18,29,39,43 and 50) cover the most of the parameter space of Ly nebulae properties in the QSO MUSEUM survey (Arrigoni Battaia et al. 2019a).These data (average rms of 2.6 mJy) allowed us to investigate any relation between the molecular gas content and the large-scale atomic phase (∼10 4 K) as traced by the Ly emission.The main results of this work are summarized as follows.
• The brightest and widest detection in the sample is the [C ](2-1) line (I [C ] (2−1) ∼15.7 Jy km s −1 , FWHM∼943 km s −1 ) associated to ID 29, which also presents the dimmest and smallest Ly nebula in the sample.
• We obtained molecular gas masses  H 2 from the CO and with IDs 7 and 18.Under these assumptions objects with ID 3 and 29 with the most massive molecular reservoirs also present the highest dust masses.
[C ](2-1) constraints (Sections 4.1 and 4.2).Applying the two estimates together, the molecular masses for sources with detections are in the range M H 2 = (0.4 − 6.9) × 10 11 M , while for quasars with non-detections, the upper limits are M H 2 < 1.1 × 10 11 M .These masses are on the high side of the typical ranges found in the literature (∼10 9 -10 11 M ) for other  ∼ 3 quasars, with ID 29 and ID 3 being outliers.
• The systemic redshift of 6 quasars is refined using the CO(6-5) and [C ](2-1) detections.These new systemic redshifts have an average uncertainty of 74.8 km s −1 .
• We found significant offsets between the extended Ly nebulae (Ly emission extracted within the APEX beam) and the newly estimated systemic redshifts in five sources (quasars with IDs 3, 8, 39, 43 and 50).Two of these sources also show the largest values of FWHM Ly in our sample (> 1000 km s −1 ), which is otherwise characterized by relatively quiescent kinematics (median FWHM Ly = 915 km s −1 ).In Section 6.2 we discuss how these velocity shifts could be signatures of bulk inflows/outflows.For the quasars with IDs 8, 39 and 50, the nebular Ly line shows a blueshift of -423 ± 76, -152 ± 64 and -377 ± 44 km s −1 , respectively, which could indicate a large-scale inflow.In contrast, quasars with IDs 3 and 43 have the bulk of the Ly line redshifted by 786 ± 109 km s −1 and 1236 ± 133 km s −1 , respectively, which can indicate that the bulk of the large-scale gas is outflowing in both cases.For ID 3, this scenario is further strengthened by the fact that this quasar is radio-loud.
• We found that the two most massive molecular reservoirs in our sample (quasars with IDs 3 and 29, M H 2 = (2.3 − 6.9) × 10 11 M and M H 2 = (3.2− 4.4) × 10 11 M , respectively) are associated with the dimmest and smallest Ly nebulae.This suggests that the quasar host galaxy properties are key in understanding the powering and physics of surrounding Ly nebulae.Obscuration from the host galaxy, due to physical properties (e.g., higher molecular gas fraction, higher dust fraction) or geometry, could reduce the escape of ionizing and Ly photons emitted by the quasar, ultimately reducing the emission from the cool CGM.This interpretation is further strengthened by the detection of a double-peaked [C ](2-1) line emission for ID 29, which could be due to an inclined massive molecular disk, likely misaligned with the ionization cones of ID 29 (see Section 6.1.2).Another possibility is that these quasars with more massive molecular reservoirs could be hosted by more massive dark-matter halos.Such halos are expected to be warmer and to show a smaller fraction of cool CGM gas, which in turn would result in lower values of Ly emission.
These APEX observations therefore suggest that the quasar host-galaxy properties could affect the Ly radiative transfer, thus regulating the amount of detectable emission on large scales around quasars.Overall, this work stresses the importance of simultaneously studying the physical properties of the multiphase gas reservoirs surrounding quasars.By characterizing the quasar host galaxy properties and system geometry (e.g., systemic redshift, inflow/outflow), future studies will be able to assess the radiative transfer effects affecting the quasar radiation and ultimately better constrain the powering mechanisms for extended gas reservoirs.Deep high-resolution interferometric observations (e.g., with ALMA, NOEMA) are indispensable to map the extent and geometry of the molecular reservoirs, ultimately linking its physics to the other CGM phases.

APPENDIX A: STABILITY OF DATA
We demonstrated the stability of the SEPIA180 receiver by checking if the noise scales down properly when using larger velocity bins.In Figure A1 we show the rms of the final reduced spectra as a function of the bin size, for each of our sources in the different observed tunings.The smallest bin size corresponds to the original resolution of the spectra (i.e., the unbinned data).The solid lines represent the expected variation of the rms (rms expected ) with respect to the bin size, which was estimated as: rms expected = rms original × √︃ binsize original /binsize observed (A1) where rms original (binsize original ) is the rms (bin size) at the original resolution, and binsize observed is the bin size used in the final combined spectrum.Each panel has a lower subpanel, which shows the deviation of the observed rms from the predicted value, estimated as (rms observed − rms expected )/rms expected , where rms observed is the final rms obtained in the final binned spectrum.At a bin size of 300 km s −1 , we obtained a median deviation of only 12% and 14% for CO(6-5) and [C ](2-1), respectively.Such a small deviation ensures that we can preform the analysis discussed in the main text.
Further, as a proof of the stability over time of the calibration of the receiver, in Fig. A2 we show the rms for each scan versus its PWV for the source with ID 39 (J 0100+2105).Each colour indicates a different observation date, spanning several months of APEX operations.These values correspond to both the CO(6-5) and [C ](2-1) -CO(7-6) observations.A linear fit shown by the grey solid line has been applied to the points, with a resulting squared correlation coefficient of  2 = 0.93 and a residual standard error of  = 0.009 K.The dashed lines enclose 95 per cent of The negative spectrum has one feature at +1200 km s −1 with peak S/N>2, but its integrated S/N is only S/N int = 2.3.c The negative spectrum has one feature at +3500 km s −1 with peak S/N>2, but its integrated S/N is only S/N int = 1.9.d The negative spectrum has one feature at +2000 km s −1 with peak S/N>2, but its integrated S/N is only S/N int = 2.5.e This jackknife spectrum has one feature with peak S/N>2, but its integrated S/N is only S/N int = 2.8.

Figure 1 .
Figure 1.The peak specific luminosity at the Ly line of the quasar versus  1450 for the QSO MUSEUM sample (Arrigoni Battaia et al. 2019a).The symbols sizes are proportional to the area of the Ly nebulae, while the colours indicate their average Ly surface brightness.We highlight the nine sources observed with APEX with black circles and their ID numbers.

Figure 2 .
Figure 2. 40 × 40 (or ∼ 312 kpc × 312 kpc) MUSE Ly surface brightness maps of the nebulae surrounding the 9 targets of this study (adapted from Arrigoni Battaia et al. 2019a).The numbers in the bottom right corner are the IDs from the QSO MUSEUM survey (see Table1, first column).The white dotted circles represent the APEX/SEPIA180 beam sizes for the CO(6-5) observations (∼32"), while the white dashed circles are the beam sizes for the [C ](2-1) -CO(7-6) observations (∼30").The cyan diamond indicates that the source is radio-loud.
, we found [C ](2-1) emission line detections in quasars with IDs 3, 29, and 39.They have fluxes of 2.3≤  [C ] (2−1) ≤15.7 Jy km s −1 .The quasar with ID 8 (or UM683) has a feature at 196.1 GHz (or -328 km s −1 , integrated S/N = 3.1,  line ∼3 Jy km s −1 , which could be a spurious line given the results of the jack-knife test for this spectrum (see Appendix B).For completeness, we also report that the quasars with ID 43 (or CTSH22.05) and ID 50 (or SDSS J0819+0823) have tentative features at 195.4

Figure 4
Figure 4.The seven APEX CO(6-5) spectral line observations.Each observed spectrum is shown with bin sizes between 150 and 200 km s −1 (left column) and 300 km s −1 right (column).In each panel, the vertical black error bar represents the uncertainty (rms) per bin, while the ID number of each quasar is shown in the top right corner.The orange vertical lines at the velocity zero point represent the systemic redshift of the sources before our APEX observations (from C , see section 3.2), and the shaded regions correspond to its uncertainty (∼415 km s −1 , ArrigoniBattaia et al. 2019a).The purple vertical lines represent the nebulae Ly redshift from MUSE extracted within the APEX beam (see Fig.2), with their corresponding uncertainty (magenta shaded regions).The green shaded area indicates the detected CO(6-5) emission line, while the dark green curve represents the Gaussian fit applied to this line, when possible.

Figure 5 .
Figure 5. Same as Figure 4, but for the eight APEX CO(7-6)-[C ](2-1) spectral line observations.Here, the zero velocities are placed at the expected [C ](2-1) transition, considering the systemic redshift from C (orange vertical line) with its uncertainty (shaded region) (see section 3.2).The grey shaded areas indicate a likely spurious source (ID8, see section 3.1 and Appendix B) or tentative features (IDs 43 and 50) that do not belong to the quasar hosts.For sources with IDs 3 and 29, we applied a double Gaussian fit, where the black dashed lines represent the single Gaussian components.For ID 3, a second component was added not to overestimate the integrated flux of the [C ](2-1) line.The emission lines of quasar with ID 29 are both due to [C ](2-1) (see section 3.1), and its systemic redshift (orange line) lies outside the plotted window.

Figure 6 .
Figure 6.CO SLEDs obtained from LVG modelling, for T kin = 30 K, using RADEX (van der Tak et al. 2007).The CO line flux normalized to the CO(1-0) line is plotted as a function of the upper rotational quantum number J. Different values in H 2 density are shown by the colour bar.

Figure 7 .
Figure 7. CO SLEDs (at T kin = 30 K) obtained with the joint CO and [C ] constraints, for the sources with IDs 8, 43 and 50 (top, middle and bottom panels, respectively), the SLEDS that can fit the joint constraints are represented by the coloured curves.The ladders show the CO line flux normalized to the CO(6-5) line as a function of rotational quantum number J, and the vertical black arrow represents the observed upper limit on the CO(7-6)/CO(6-5) ratio for each source.The different values in H 2 density are shown by the colour bar.The grey curves in the background show the discarded CO SLEDs (from Fig. 6).

Figure 8 .
Figure 8. CO SLEDs for the sources with IDs 8, 43 and 50 (discontinuous lines; see legend) compared to the SLEDs for the quasars Cloverleaf, MG0751 and PSS14, obtained by Weiss et al. (2007a) (continuous grey lines).The ladders correspond to the CO line flux normalized to the CO(1-0) line as a function of rotational quantum number J.
1.1 for further discussion).As halos of different masses are expected to be characterised by different fractions of cool and hot gas (e.g.,Dekel & Birnboim 2006), the FWHM of molecular lines could show important trends with the Ly properties.For instance, more massive halos should in principle show larger FWHM of the molecular lines and smaller fractions of cool gas, with consequently smaller Ly luminosities (L Ly ) and Ly areas compared to less massive halos.We start by comparing the FWHM of the observed molecular lines with the total L Ly , the average Ly surface brightness SB Ly , and the area encompassed by the Ly emission by the nebulae surrounding the quasars in our sample.Figure9shows this comparison for the different sources in our sample (represented by different colours).Specifically, the left panel shows the FWHM of the molecular lines as a function of SB Ly (corrected for the cosmological dimming), the central panel shows the FWHM versus L Ly , and the right panel shows the FWHM versus the area of the Ly nebulae.The legend in each panel indicates the different markers used for the different molecular lines (CO(6-5) and [C ](2-1)).

−Figure 10 .
Figure 10.Molecular gas mass M H2 derived from the joint CO and [C ] constraints versus Ly surface brightness corrected for cosmological dimming (first panel), the total Ly luminosity (second panel), and physical area (third panel) for the large-scale Ly nebulae discovered around the quasars in our sample (Arrigoni Battaia et al. 2019a).The ID numbers of the sources are shown in the first panel.The vertical dashed lines represent ranges of M H2 , and the arrows indicate upper limits (see Table4).The horizontal lines show the expected location of massive star-forming main sequence MS, 4×MS and 10×MS galaxies using the gas fractions inTacconi et al. (2018) and the halo mass-stellar mass relation inMoster et al. (2018).The grey area is the uncertainty on this expectation, for 1×MS.The histogram show the M H2 distribution for our sample considering the average values of the mass ranges (red), and for a sample of  ∼2.5 -3 quasars (blue) extracted from the literature(Barvainis et al. 2002;Weiß et al. 2003;Beelen et al. 2004;Walter et al. 2011;Schumacher et al. 2012;Hill et al. 2019;Bischetti et al. 2021), whose M H2 have been estimated with different molecular gas tracers and have a typical uncertainty of ∼ 0.1 dex.

Figure 11 .
Figure 10.Molecular gas mass M H2 derived from the joint CO and [C ] constraints versus Ly surface brightness corrected for cosmological dimming (first panel), the total Ly luminosity (second panel), and physical area (third panel) for the large-scale Ly nebulae discovered around the quasars in our sample (Arrigoni Battaia et al. 2019a).The ID numbers of the sources are shown in the first panel.The vertical dashed lines represent ranges of M H2 , and the arrows indicate upper limits (see Table4).The horizontal lines show the expected location of massive star-forming main sequence MS, 4×MS and 10×MS galaxies using the gas fractions inTacconi et al. (2018) and the halo mass-stellar mass relation inMoster et al. (2018).The grey area is the uncertainty on this expectation, for 1×MS.The histogram show the M H2 distribution for our sample considering the average values of the mass ranges (red), and for a sample of  ∼2.5 -3 quasars (blue) extracted from the literature(Barvainis et al. 2002;Weiß et al. 2003;Beelen et al. 2004;Walter et al. 2011;Schumacher et al. 2012;Hill et al. 2019;Bischetti et al. 2021), whose M H2 have been estimated with different molecular gas tracers and have a typical uncertainty of ∼ 0.1 dex.

Figure 14 .
Figure 14.APEX [C ](2-1) spectral line observation for Q-0115-30 (quasar with ID 29).The spectrum is shown with a bin size of 80 km s −1 , the highest resolution possible.The red vertical line represent the systemic redshift estimated from the molecular lines.The curves and symbol follow the same notation as Figure 5.

Figure 15 .
Figure 15.MUSE white-light image for the field around ID 29 (or Q-0115-30, at the center) highlighting the location of continuum sources within the APEX beam for the [C ](2-1) observations (white dashed circle).Low- interlopers (cyan circles) and sources with unknown redshift (green circles) are indicated (see Section 6.1.2for details).

Figure A2 .
Figure A2.RMS values (for each scan) versus PWV relation for the source SDSS J0100+2105 (ID 39), corresponding to the CO(6-5) and [C ](2-1) -CO(7-6) observations.Each colour indicates a different observation date.The grey solid line represents a linear fit to the points, and the two dashed lines represent the ±2 (±0.018K) deviation from the fit.

Table 1 .
Battaia et al. 2019a setup and log.Quasar systemic redshift taken from the QSO MUSEUM survey (ArrigoniBattaia et al. 2019a d

Table 3 .
Redshifts and central frequencies obtained from the targeted molecular lines, and velocity shifts from the large-scale Ly emission within the same APEX aperture.Redshift obtained from the observed CO(6-5) emission line.bRedshiftobtainedfrom the observed [C ](2-1) emission line.cRedshiftobtainedfrom the observed CO(7-6) emission line.dQuasarsystemicredshiftestimatedusing the observed molecular lines.eVelocityshiftbetweentheLy line peak of the nebulosities (extracted within the APEX beam, see Fig.2) and the systemic redshift  mol,sys .‡Thesevalues have been calculated as the average between the two Gaussian components (see Fig.5).For each component, the individual values of  C (2−1) are 191.91±0.11GHz and 191.39±0.06GHz. a

Table 4 .
Molecular, carbon and dynamical mass estimations.

Table B1 .
Numbers of spurious lines in the negative and jack-knife tests.