The cold molecular gas in z ≳ 6 quasar host galaxies

Context. Probing the molecular gas reservoirs of z ≳ 6 quasar (QSO) host galaxies is fundamental to understanding the coevolution of star formation and black hole growth in these extreme systems. Yet, there is still an inhomogeneous coverage of molecular gas tracers for z ≳ 6 QSO hosts. Aims. To measure the average excitation and mass of the molecular gas reservoirs in the brightest z > 6 . 5 QSO hosts, we combined new observations of CO(2–1) emission with existing observations of CO(6–5), CO(7–6), [C i ](2–1), [C ii ]158 µ m, and dust-continuum emission. Methods. We reduced and analysed observations of CO(2–1), taken with the Karl G. Jansky Very Large Array, in three z = 6 . 5 − 6 . 9 QSO hosts—the highest redshift observations of CO(2–1) to date. By combining these with the nine z = 5 . 7 − 6 . 4 QSO hosts for which CO(2–1) emission has already been observed, we studied the spread in molecular gas masses and CO excitation of z ≳ 6 QSOs. Results. Two of our three QSOs, P036 + 03 and J0305–3150, were not detected in CO(2–1), implying more highly excited CO than in the well-studied z = 6 . 4 QSO J1148 + 5251. However, we detected CO(2–1) emission at 5 . 1 σ for our highest-redshift target, J2348– 3054, yielding a molecular gas mass of (1 . 2 ± 0 . 2) × 10 10 M ⊙ , assuming α CO = 0 . 8(K km s − 1 pc 2 ) − 1 and r 2 , 1 = 1. This molecular gas mass is equivalent to the lower limit on the dynamical mass measured previously from resolved [C ii ]158 µ m observations, implying that there is little mass in stars or neutral gas within the [C ii ]-emitting region and that a low CO-to-H 2 conversion factor is applicable. On average, these z ≳ 6 QSO hosts have far higher CO(6–5)-, CO(7–6)-, and [C ii ]158 µ m versus CO(2–1) line ratios than the local gas-rich and IR-luminous galaxies that host active galactic nuclei, but with a large range of values, implying some variation in their interstellar medium conditions. We derived a mean CO(6–5)-to-CO(1–0) line luminosity ratio of r 6 , 1 = 0 . 9 ± 0 . 2. Conclusions. Our new CO(2–1) observations show that even at 780 Myr after the Big Bang, QSO host galaxies can already have molecular gas masses of 10 10 M ⊙ , consistent with a picture in which these z ≳ 6 QSOs reside in massive starbursts that are coevolving with the accreting supermassive black holes. Their high gas versus dynamical masses and extremely high line excitation imply the presence of extremely dense and warm molecular gas reservoirs illuminated by strong interstellar radiation fields.


Introduction
First discovered in optical surveys in the early 2000s, luminous z ≳ 6 quasars (QSOs) have challenged our understanding of how supermassive black holes (SMBHs) grow and coevolve with their host galaxies (see Fan et al. 2023, and references therein).Studies of the optical and near-infrared (NIR) spectra of these QSOs reveal that they host SMBHs with masses of ≳ 10 9 M ⊙ (e.g.Jiang et al. 2007;Kurk et al. 2007;De Rosa et al. 2014;Wu et al. 2015), implying the presence of enormous black hole seeds at z > 15 and/or early episodes of super-Eddington accretion (Latif & Ferrara 2016;Inayoshi et al. 2017;Trakhtenbrot 2021).To understand how these SMBHs grow and impact star formation, it is crucial to constrain the mass and properties of their gas reservoirs.Doing so will also help reveal the evolution-ary path of the gas-poor elliptical galaxies that host today's most massive SMBHs (M BH ≳ 10 10 M ⊙ , Kormendy & Ho 2013).
To discern the type of galaxies in which z ≳ 6 QSOs reside, they have been observed at radio and (sub)millimetre wavelengths, where the emission from the host galaxies' interstellar medium (ISM) dominates that of the QSO.Observations of their dust-continuum emission have revealed that many of the hosts are bright in the far-infrared (FIR), with high FIR luminosities of L FIR ∼ 10 12 − 10 13 L ⊙ , which have been interpreted as signatures of significant dust heating associated with high star formation rates (e.g.Bertoldi et al. 2003a;Beelen et al. 2006;Riechers et al. 2007;Decarli et al. 2017).The best studied of these appear to exhibit moderately high dust temperatures of 30 − 60 K (e.g.Leipski et al. 2013;Witstok et al. 2023) and dust masses of ≳ 10 9 M ⊙ (e.g.Venemans et al. 2017a;Shao et al. 2019;Meyer et al. 2022), implying large molecular gas masses of ≳ 10 10 M ⊙ (assuming similar or higher gas-to-dust mass ratios than in the Milky Way) and efficient supernova dust production (e.g.Michałowski et al. 2010;Gall et al. 2011;Kuo & Hirashita 2012).Observations of the dominant ISM coolant, [C ii] 158 µm, also show this emission to be bright, again indicating high SFRs (Walter et al. 2009;Venemans et al. 2012;Wang et al. 2013;Bañados et al. 2015;Venemans et al. 2016).Combined with the underlying continuum emission, these results yield [C ii]/FIR luminosity ratios similar to high-redshift starbursts and local Ultra-Luminous Infrared Galaxies (ULIRGS).Moreover, the recent wealth of resolved studies of [C ii] (e.g.Willott et al. 2017;Wang et al. 2019;Venemans et al. 2020;Novak et al. 2020;Neeleman et al. 2021), most recently pushing down to scales of ∼ 200 − 300 pc (Walter et al. 2022;Meyer et al. 2023), have revealed the presence of some rotating gas disks, in which most of the gas is concentrated in highly compact, central regions.
For some z ≳ 6 QSO hosts, observations of the warm and dense molecular gas component traced by  and/or CO(7-6) emission have also revealed the presence of massive molecular gas reservoirs of ∼ 10 10 M ⊙ (Bertoldi et al. 2003b;Carilli et al. 2007;Wang et al. 2010Wang et al. , 2011aWang et al. , 2013;;Venemans et al. 2017b).Yet, converting these observations to a molecular gas mass requires robust constraints on the CO excitation ladders.To test how much molecular gas is in the cold/diffuse state and hence correctly calibrate the CO excitation for these z ≳ 6 QSO hosts, it is critical to also observe their low-J (J = 1, 2) CO emission.But low-J CO emission is intrinsically much fainter than the mid-J CO and [C ii] emission.To date, nine QSOs at z = 5.7−6.4 have observations of CO(2-1) emission (Wang et al. 2010(Wang et al. , 2011b;;Stefan et al. 2015;Wang et al. 2016;Shao et al. 2019).In combination with their previously detected CO(5-4), , and/or CO(7-6) emission, these observations have revealed the presence of highly-excited CO, with what seems to be a dearth of diffuse cold gas.
To continue systematically exploring the molecular gas reservoirs of z ≳ 6 QSO host galaxies, we used the Karl G. Jansky Very Large Array (VLA) to observe the CO(2-1) emission of three QSOs with existing observations of CO(6-5), CO(7-6), [C i] and [C ii] emission (Bañados et al. 2015;Venemans et al. 2016Venemans et al. , 2017a)).By combining these lines, we robustly measured their molecular mass, anchored the CO excitation, and compared their ISM properties to those of highly star-forming galaxies and QSO hosts at lower redshift.
This paper is organised as follows.In Section 2, we present the VLA observations of our three QSO hosts and present the literature sample.In Section 3, we infer physical properties from these observations.We place these properties in the context of other galaxy samples in Section 4, and discuss the physical implications in Section 5.In Section 6, we present a short summary.Throughout this work, we adopt a Lambda-cold dark matter cosmology, with H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3, and Ω λ = 0.7.

Observations
The CO(2-1) observations analysed here were taken as part of the VLA programmes, 16A-160 and 17A-336 (PI: Venemans).For programme 16A-160, 11 observational blocks were conducted between 2016/03/04 and 2017/05/17, resulting in 10 h of observations of J2348-3054, 2 h of J0305-3150 and no observations of P036+03.To complete the experiment, 20 further observational blocks were requested as part of the programme 17A-336, with 18 h to be spent on J2348-3054, 12 h on J0305-3150 and 10 h on P036+03.This supplementary programme was carried out in full between 2017/05/20 and 2017/08/08.Across the two programmes, almost all observation blocks were carried out in the C configuration, with one instead conducted in the B configuration.All observations were carried out in the Ka Band, spanning 26.5-40 GHz.

Data reduction
To calibrate the data from each two-hour observation set, we used the Common Astronomy Software Application (CASA Team et al. 2022) (CASA) pipeline version 5.6.2 as described in https://science.nrao.edu/facilities/vla/data-processing/pipeline#section-26 (that is, we included most of the standard calibration tasks and disabled the application of Hanning smoothing).After examining the calibrated data, we performed additional flagging, removing the RFI at 29.43 GHz remaining in the observation sets of J2348-3054 and the poor individual antenna data for short time windows in two observation sets (one for J0305-3150 and one set for J2348-3054).This leads to some of the gaps in data in the spectra of Fig. 1.
To test whether the CO(2-1) emission is detected, we made three-channel cubes, by integrating each channel over 1.2× the FWHM of the previously detected [C ii] 158 µm line emission (see Table 1) and centring the frequency range on the CO(6-5)-derived redshift.Assuming that the CO(2-1) emission has the same linewidth as [C ii] 158 µm and systemic velocity of the other CO transitions, this velocity range should encompass 84% (1.9 ± 0.9) × 10 10 (2.4 ± 1.0) × 10 10 (1.6 ± 0.9) × 10 10 M dust (M ⊙ ) (1.2 ± 0.7) × 10 8 (2.3 ± 0.7) × 10 9 (1.1 ± 0.6) × 10 8  2022), * from a Gaussian fit to the spectrum extracted from the NOEMA observations of programmes S15DA and W15FD (also used in Decarli et al. 2022) of the CO(2-1) emission.We applied a natural weighting scheme and cleaned down to 2σ within 3 ′′ . of the source centre.These cleaned maps are shown in Fig. 1, with the beam properties and rms values provided at the bottom of Table 1.Using the same central frequency, we also created image cubes of 100 km s −1 , cleaning down to 2σ and applying natural weighting.The spectra extracted from the single pixel corresponding to the peak in the moment-0 maps of these cubes, are shown in the right panels of Fig. 1.

Line fluxes and upper limits
We tested whether the CO(2-1) emission of the three QSO hosts was detected using the moment-0 maps described above.For our highest-redshift target, J2348-3054, we detect a 5.1σ peak in emission, consistent with the expected source position.We find a 3.3 σ peak offset by 0 ′′ .43 (2.4 kpc) from the source centre for P036+03; however, other peaks of emission exceeding this value are apparent in the integrated channel maps at the same frequency and in the channels on either side.Thus, we treat it as an upper limit.Likewise, we detect no CO(2-1) emission for J0305-3150, with a peak value near the source centre of 2.6 σ.
At a spectral resolution of 100 km s −1 , none of the three sources exceed S/N = 3 across the expected line width.Moreover, none of the three QSOs targeted show any evidence of continuum emission, nor is there any spectrally offset CO(2-1) emission.We also searched for evidence of the companions reported in Venemans et al. (2020), but find nothing exceeding ∼ 2.5 σ at the expected positions and frequencies.
Based on the moment-0 maps, we treated the CO(2-1) emission from J2348-3054 as a significant detection and derived 4σ upper limits for the other two sources, consistent with our detection threshold.We derived the line flux measurement for J2348-3054 by first measuring the value of the single pixel at the peak position (0 ′′ .3 offset from the source centre quoted in Table 1) in the central channel of the bottom row in Fig. 1, multiplying by the channel width, and accounting for the extra 16% of the line flux not included in 1.2 × FWHM.To derive the upper limits, The spectra are shown for channels of 100 km s −1 velocity width, the shaded region indicating the 1.2 FWHM region used to create the channel maps to the left.P036+03 and J0305-3150 are undetected (with a significant number of negative and positive peaks in the data cube exceeding the significance of the peaks in emission near the expected source centres), so we place 4σ upper limits on their CO(2-1) emission.J2348-3054 is detected at 5.1σ (central velocityintegrated map in bottom row).
we measured the root mean square (rms) of the moment-0 map, also multiplying this value by 1.2 × FWHM [CII] , accounting for the extra 16%, and multiplying by our detection threshold of 4σ.We also accounted for potentially extended emission underlying the noise, assuming that the CO(2-1) emission arises from the same size region as the [C ii] 158 µm emission.We deem it unlikely that the CO emission would be significantly more extended than the [C ii] 158 µm emission (see also Riechers et al. 2011a).To test how much flux we would miss by deriving an upper limit (line flux) from only the pixel rms (peak pixel value), we followed the approach of Carniani et al. (2020).That is, we simulated both a point source and extended (Gaussian) emission, using the [C ii] 158 µm size (2 r 1/2,[C II] , Table 1), position angle, and inclination measured for each galaxy.We then performed mock VLA observations under the same configuration used for the real observations and compared the flux extracted from the peak pixel for both the point source and extended source.We performed noise-free tests to accurately recover the difference between the flux of the point versus extended source (which we cannot recover from our data).Based on this modelling, we found that for P036+03, we would miss 15% of the flux if the CO(2-1) emission is as extended as the [C ii] 158 µm emission.
Similarly, for J0305-3150, we would miss 19% of the flux.We added these potentially missing extended components when deriving the 4σ upper limits, based on the rms.For J2348-3054, the peak pixel extraction for the extended source is equivalent to that of the point source (only 0.004% less), because the CO(2-1) beam (7.8 × 3.6 kpc), is significantly larger than twice the measured [C ii] 158 µm half-light radius (∼ 0.8 kpc) (Neeleman et al. 2021).This is also consistent with the fact that the CO(2-1) emission from J2348-3054 is well fit by a (convolved) point source and the flux extracted from this point source matches that extracted from a larger circular aperture of 2 ′′ . 5. Thus, we performed no correction for any missing flux for J2348-3054.
Given the low S/N of these observations, we derived the CO(2-1) line flux density and upper limits assuming that the CO(2-1) line width is the same as that measured for the much higher S/N [C ii] 158 µm emission; that is, we use the highest S/N prior we have.This seems to be a reasonable assumption, as for most of the z > 5 QSO host galaxies detected at > 5σ in both CO(2-1) and [C ii] 158 µm, the reported FWHM values of the two lines have been consistent within uncertainties (see Table A.1 for literature sample, values reported in Wang et al. 2013Wang et al. , 2016;;Shao et al. 2019;Khusanova et al. 2022;Meyer et al. 2022).Although we use the most recent measurements of the [C ii] 158 µm line FWHM (Venemans et al. 2020;Li et al. 2022) for the integrated CO(2-1) line flux density and upper limits, we also tested our results using the values derived in the original studies that reported the [C ii] 158 µm detections (Bañados et al. 2015;Venemans et al. 2016).The integrated CO(2-1) line flux densities (hereafter simply 'line fluxes') and upper limits remain consistent, within errors.However, the S/N of the detection of J2348-3054 increases with the wider FWHM measured in Venemans et al. (2020) versus Venemans et al. (2016).

Line luminosities and dust properties
For the three galaxies studied here, and for the comparison samples, we computed the line luminosities as where D L is the luminosity distance in Mpc, ν obs is the observed central frequency of the line in GHz and S line ∆v is the line flux in Jy km s −1 .At the redshift of these QSOs, the high temperature of the cosmic microwave background (CMB) may have a pronounced effect on the observed CO line fluxes, becoming a strong background against which lines are measured as well as potentially boosting the higher-versus lower-J transitions (da Cunha et al. 2013a;Zhang et al. 2016;Tunnard & Greve 2016).We therefore derived both CMB-corrected and uncorrected values for the CO line luminosities.As shown in da Cunha et al. (2013a), the ratio of the observed versus intrinsic line flux depends on the temperature and density of the gas.To perform the CMB corrections, we adopted their non-LTE, high-density (n H 2 = 10 4.2 cm −3 ) and high-temperature (T kin = 40 K) scenario.In their model, they assume that the kinetic temperature of the gas is coupled to that of the dust: T k = T d .Thus, we partly motivate the choice of this scenario by the high inferred dust temperatures (∼40 K), as well as the compact gas reservoirs and high observed line ratios, which imply high densities and temperatures (see Sec. 5).Adopting this correction, the fraction of CO(2-1) line flux observed against the CMB, f CMB , for P036+03, J0305-3150, and J2348-3054 is 0.56, 0.55, and 0.53 (respectively), whereas the fraction of CO(6-5) observed against the CMB is 0.67, 0.65, and 0.64 (respectively).This additional correction results in higher intrinsic line fluxes, and hence luminosities, than without the CMB correction.
We also refit the dust SEDs, ensuring that we take into account the effects of the CMB by adopting a similar approach to da Cunha et al. ( 2021) and Meyer et al. (2022) (see Appendix B).We fit the dust SEDs with a modified blackbody, assuming optically thin dust emission at λ > 40 µm.Our method accounts for the effects of the CMB, correcting for both the heating by and contrast against the CMB (as prescribed by da Cunha et al. 2013a) and includes the most up-to-date continuum data.In this way, we obtained new values for the FIR luminosities, dust masses, dust temperatures, and emissivity indices of all z ≳ 6 QSO hosts, including our three and the literature sample (Table B.1).We note that in most cases, the dust temperatures and emissivity indices remain unconstrained due to the limited number of continuum data points.For these, the errors on the dust masses and FIR luminosities (L FIR ) are higher, taking into account the full prior range on T dust and β.We provide two sets of L FIR values to ease the literature comparison, integrating between 40-400 µm and 42.5-122.5µm (e.g. as in Liu et al. 2015versus De Breuck et al. 2011).
If the QSO host galaxies have high dust columns, the dust may actually be optically thick at λ > 40 µm, contrary to our assumption.Given the large dust masses of 10 8 − 10 9 M ⊙ inferred here (albeit for the optically thin assumption) and small galaxy sizes of r half ∼ 1 kpc (as in the measurements for the highestquality [C ii] 158 µm observations of J2348+3054, Walter et al. 2022) this is certainly plausible.As shown by da Cunha et al. ( 2021), accounting for the wavelength up to which dust remains optically thick (λ thick = 60 − 140 µm for their sample) has the effect of shifting the best-fit dust temperatures of their sample 10 K higher, on average, and dust masses to a few times lower values.However, the derived FIR luminosities and β remain practically the same.Thus, our L FIR should not change significantly in the case of optically thick dust emission but our dust temperatures and masses may be biased to slightly lower and higher values, respectively.

Molecular gas masses
We derived the molecular gas masses from the CO(2-1) emission assuming that the low-J emission is thermalised, L ′ CO(2−1) = L ′ CO(1−0) , as in low-redshift QSOs (e.g.Molyneux et al. 2023).To convert L ′ CO(1−0) to a molecular gas mass, we assumed an intrinsic ULIRG-like α CO = 0.8 M ⊙ (K km s −1 pc 2 ) −1 , as derived in Downes & Solomon (1998) (for ULIRG centres) based on resolved CO(2-1) or CO(1-0) observations of 10 local ULIRGs.These low values were found to be necessary in order for the molecular gas masses not to exceed the inferred dynamical masses.Similarly, early comparisons of dynamical mass estimates and CO-derived molecular gas masses in z ≳ 6 QSO hosts motivated the choice of such low α CO also in these systems (e.g.Walter et al. 2003Walter et al. , 2004)).Following the majority of the literature on z ≳ 6 QSOs (see the review of Carilli & Walter 2013), we adopted the same value but it is important to note that this may need to be adjusted in future, once robust constraints on α CO can be determined for each QSO.For example, by modelling the multi-line and continuum emission, Harrington et al. (2021) show that for systems in which the CO gas column density is mostly traced by higher-J lines, α CO is shifted to higher values.We also inferred the molecular gas masses based on the CMBcorrected CO(2-1) emission (see Sec. 3.1).This additional correction results in a higher intrinsic molecular gas masses than without the CMB correction (although for lower assumed temperatures and densities these corrections would be much larger).In effect, this additional CMB correction leads to an increased α CO , that is α CO,tot.= α CO,intr./ f CMB .Thus, Milky-Way-like values of α CO ∼ 4 M ⊙ (K km s −1 pc 2 ) −1 are feasible for host galaxies with lower average molecular gas temperatures or densities than in the model we have assumed here.
Fig. 2: Comparison of molecular gas masses inferred from different sets of observations for our three QSO hosts (top three) and the literature sample of z ≳ 6 QSO hosts (bottom nine).For the CO-based values, we have assumed a ULIRG CO-to-H 2 conversion factor of α CO = 0.8 M ⊙ (K km s −1 pc 2 ) −1 .We show two masses for CO(2-1), one derived without any CMB correction being applied (filled blue circles) and one derived by applying a correction assuming the high-density (n H 2 = 10 4.2 cm −3 ) and high-temperature (T kin = 40 K), non-LTE scenario presented in da Cunha et al. (2013a) (blue outlined circles).We also show molecular gas masses derived from [C i] (2-1) emission (green squares) assuming the [C i] is optically thin, in LTE and has the same abundance as the sample of lensed starbursts studied in Harrington et al. (2021), namely [C i]/H 2 = (6.82± 3.04) × 10 −5 .We also compare to molecular gas masses derived from the dust masses fit here, assuming the local gas-to-dust ratio of 100 (red triangles), as well as molecular gas masses derived from the [C ii] 158 µm line luminosities, assuming the mean Zanella et al. (2018).
Dust emission can also be used as a molecular gas tracer, as long as the dust temperature, emissivity index, and dust-to-gas ratios can be well-constrained or safely assumed (e.g.Groves et al. 2015;Scoville et al. 2016;Liang et al. 2018).In this work, we derived the molecular gas masses from the dust masses fit here (see Appendix B).To derive molecular gas masses, we scaled the dust masses by the local gas-to-dust ratio of 100 (e.g.Draine et al. 2007), which has been found to be consistent with some z ≳ 6 QSO hosts (e.g.Wang et al. 2016;Decarli et al. 2022;Feruglio et al. 2023).
Fig. 3: Dynamical masses (black squares), dust masses (red triangles), and black hole masses (green stars) compared to the CO(2-1)-derived molecular gas masses assuming a ULIRG α CO , with a CMB correction (light shaded symbols) and without a CMB correction (dark shaded symbols).For the CO(2-1)derived values, we assume that the low-J CO emission is ther- ).Most molecular gas masses are high, spanning 0.1 − 1× the dynamical mass.
In recent years, atomic carbon, [C i], has also garnered significant attention as a potentially reliable and efficient molecular gas tracer for z > 2 star-forming galaxies (e.g.Walter et al. 2011;Alaghband-Zadeh et al. 2013;Bothwell et al. 2017;Yang 2017;Andreani et al. 2018;Valentino et al. 2018;Nesvadba et al. 2019).At these high redshifts, the two [C i] transitions are shifted to millimetre wavelengths and can be detected in less time than CO(2-1).Unlike for 12 CO, [C i] emission has mostly been found to be optically thin, (e.g.Ikeda et al. 2002;Weiß et al. 2003;Harrington et al. 2021), potentially making it a more straightforward tracer.Moreover, [C i] emission is widespread within clouds, stemming from both the cloud interiors and surfaces (Ikeda et al. 2002;Kramer et al. 2008;Bisbas et al. 2015Bisbas et al. , 2017;;Glover et al. 2015;Clark et al. 2019).[C i] has therefore been proposed to be a more reliable molecular gas tracer than low-J CO emission in diffuse and/or metal-poor gas (e.g.Papadopoulos et al. 2004;Glover & Clark 2016).
To derive molecular gas masses from the observed [C i] (2-1) emission, we assumed the [C i] emission is optically thin, and in LTE (following Weiß et al. 2005).In this case, the mass of neutral carbon is given by, where Q ex = 1 + 3 exp(−23.6/Tex ) + 5 exp(−62.5/Tex ) is the partition function and T ex is the excitation temperature.We assumed the atomic carbon is in thermal equilibrium with the dust, with the average dust temperature typically quoted in the literature, T ex = 47 K (Beelen et al. 2006;Leipski et al. 2014).If we instead assume T ex = 30 K, the [C i]-derived molecular gas masses would be 1.5× the values presented here, whereas if we assume T ex = 100 K the molecular gas masses would be 0.8× our values, therefore making little difference to the derived values.
To convert the atomic carbon mass to a molecular gas mass, we applied the mean atomic carbon abundance of ⟨[C i]/H 2 ⟩ = (6.82± 3.04) × 10 −5 , measured by Harrington et al. (2021) for their sample of 18, z = 2 − 3 lensed starburst galaxies with [C i] observations.This value is consistent with the implied high value of 7 × 10 −5 found by Bothwell et al. (2017) for their sample of 13 lensed z ∼ 4 dusty star-forming galaxies, when assuming the same low α CO as in this work.However, these values are higher, on average, than those measured for less IR-luminous, main-sequence star-forming galaxies at z = 1.2 − 2.5.For example, Boogaard et al. (2020) derive a value of [C i]/H 2 = (1.9 ± 0.4) × 10 −5 based on the [C i] detections of six main-sequence star-forming galaxies, and Valentino et al. (2018) find a value of [C i]/H 2 = (1.6 ± 0.8) × 10 −5 for their sample of z = 1.2 main-sequence galaxies.These lower values may be attributed to the difference in physical conditions (e.g.cosmic ray, X-ray, and/or mechanical heating effects) and/or the calibration methods. 1 We take the value of Harrington et al. (2021) here since their sample are comprised of starbursts, and thus are more similar to the z ≳ 6 QSO hosts studied to date (e.g.Wang et al. 2011aWang et al. , 2016;;Shao et al. 2019).Moreover, they do not assume an α CO , but instead model the α CO and [C i] abundance simultaneously.
[C ii] 158 µm has been proposed as an even more efficient molecular gas tracer for high-redshift and metal-poor galaxies (e.g.Zanella et al. 2018;Madden et al. 2020).However, as a multi-phase gas tracer, [C ii] 158 µm emission is an indirect probe of the bulk of the cold, molecular gas and has not yet been systematically tested against more common molecular gas tracers at z > 6.To infer the molecular gas masses from [C ii], we applied the mean conversion factor from the calibration of Zanella et al. (2018), α [CII] = 31 M ⊙ /L ⊙ , folding in the 0.3 dex standard deviation around this value as the uncertainty.Zanella et al. (2018) derive this value based on a comparison of molecular gas masses versus [C ii] 158 µm emission.For most of the z = 1 − 5 main-sequence and starburst galaxies underlying their empirical calibration, the molecular gas masses were derived from CO observations.However, for the local starbursts underlying this calibration, the molecular gas masses were instead mostly derived using relations between the depletion timescale (t depl = M mol /S FR) and specific star formation rate, sSFR (see Sec. 4.1).

Molecular gas mass comparison
Based on the molecular gas mass comparison shown in Fig. 2, we find that the CO(2-1)-based molecular gas masses of the three, z = 6.6 − 6.9 QSO hosts studied here (top three) are consistent with the z = 5.7 − 6.4 QSO hosts with existing CO(2-1) measurements (bottom nine in Fig. 2, Wang et al. 2010Wang et al. , 2011b;;Stefan et al. 2015;Wang et al. 2016;Shao et al. 2019), with M mol ∼ (0.8 − 4.5) × 10 10 M ⊙ .Even our newly CO(2-1)-detected QSO host, J2348-3054 at z = 6.9, already has a molecular gas mass of ∼ 1 × 10 10 M ⊙ .The fact that these extremely IR-bright z ≳ 6 QSO hosts all have similarly massive molecular gas reservoirs implies that there is something about these systems that enables massive molecular gas reservoirs to form very rapidly.It also implies that the radiation from the rapidly accreting SMBHs is not dissociating significant amounts of CO.
In Fig. 2, we compare the CO(2-1)-derived masses to those inferred from [C i] (2-1), [C ii] 158 µm and dust-continuum emission.We find that the molecular gas masses inferred from the dust mass measurements are mostly consistent with those inferred from CO(2-1), implying that, on average, the gas-to-dust ratios are consistent with, or only somewhat higher than, the assumed Milky Way value of 100 (also shown in Fig. 3).The only two exceptions, J0305-3150 and J2310+1855, are systems for which the best-fit dust temperatures are low (≲ 30 K).It may be that for these galaxies, the gas-to-dust ratios are actually lower (implying more plentiful dust), or that the best-fit solution is inaccurate due to optical depth effects or, poor sampling of the dust SED for J0305-3150.
For our our three z > 6.5 QSO hosts, and J1148+5251 (Riechers et al. 2009;Stefan et al. 2015), we find that the [C i] (2-1)-derived molecular gas masses are consistent with the CO(2-1)-derived values, when assuming both the ULIRG α CO value and the high atomic carbon abundance measured in Harrington et al. (2021).If we instead applied a Milky-Way-like α CO , the [C i]-derived molecular gas masses would be ∼ 5× lower in comparison.In that case, the lower carbon abundances measured by Boogaard et al. (2020) and Valentino et al. (2018) would be appropriate.If the assumed α CO is correct, this comparison indicates that the gas in these QSO hosts may be highly carbon enriched and/or subject to several sources of strong heating.We discuss the implications, including additional effects that complicate this comparison, in Sec. 5.
In contrast to the other tracers, we find that the [C ii] 158 µmbased molecular gas masses are systematically higher than those derived from the CMB-uncorrected CO(2-1), by a factor of 6.8 ± 1.9 (mean and standard deviation for the whole sample of z ≳ 6 QSO hosts).For J2348-3054, the [C ii]-derived value is 4.6× that of the CO(2-1) derived value.Assuming the other values are accurate, this implies that the mean conversion factor derived for local and high-redshift galaxies, α [CII] = 31 M ⊙ /L ⊙ , is too high for these z ≳ 6 QSO host galaxies.This may be because α [CII] is still uncertain for such highly star-forming galaxies.At high sSFR, the calibration from Zanella et al. (2018) mostly relies on local starbursts taken from Díaz-Santos et al. (2017), which were without CO observations.To derive the molecular gas masses of this sample (and hence calibrate these against the [C ii] 158 µm luminosity), Zanella et al. (2018) 2017) only (see Fig. 9 of Zanella et al. 2018).When also using the Sargent et al. (2014) relation (which is based on a metallicity-dependent α CO ) they find α [CII] values of ≤ 10 at the highest sSFRs (sSFR/sSFR MS > 10).Thus, they also derive some values consistent with what we find for the z ≳ 6 QSO host galaxies, depending on the assumed scaling relation.

Mass budget
In Fig. 3, we compare the molecular gas masses derived here against the other main mass components of the QSO host galaxies-the dust, black hole, and total (dynamical) mass.We use the dust masses derived here (Sec.3.1 and Appendix B).For the black hole masses of our three QSOs, we take the measured values from Mazzucchelli et al. (2023) and Farina et al. (2022), which were derived from Mg ii following the method from Shen et al. (2011).Likewise for the literature sample, we take the Mg ii-based values from Shen et al. (2019), Farina et al. (2022), andMazzucchelli et al. (2023).For the dynamical mass estimates, we take the values derived by Neeleman et al. (2021), based on 0 ′′ .25 observations of [C ii] 158 µm emission.These dynamical masses were estimated from the best-fit [C ii] 158 µm half-light radius, rotation velocity, and velocity dispersion (where these values were fit using the parametric kinematic fitting code, qubefit; Neeleman 2021).
Assuming a ULIRG-like α CO , the CO(2-1)-derived gas masses (opaque filled symbols with errorbars in Fig. 3) are 3-11 times greater than the black hole masses, 20 − 190× greater than the dust masses, and 0.1 − 1× the dynamical masses.In order for the CO(2-1)-derived molecular gas masses (and upper limits) of the three QSOs studied here not to exceed the estimated dynamical masses, the applied CO-to-H 2 conversion factors need to be low, similar to the ULIRG value we have adopted.This also means that only a minimal CMB correction (as for the non-LTE, high temperature and density case assumed here) is required.In particular for J2348-3054, the molecular gas mass is almost equivalent (0.98×) to the lower limit on the dynamical mass of ≥ 1.2 × 10 10 M ⊙ , whereas the black hole mass is (3.3 ± 1) × 10 9 M ⊙ .If the dynamical masses are reliable, the CO(2-1)-derived molecular gas masses also imply that the molecular gas in most of these systems, especially in J2348-3054, dominates the mass budget, leaving barely any mass in neutral gas or stars within the region traced by [C ii] 158 µm.
Our finding that the molecular gas comprises a large fraction of the total mass, at least in the central regions, is consistent with earlier findings for some of the literature sample and for J2348-3054 based on dynamical modelling of [C ii] 158 µm versus CO(6-5)-or [C ii] 158 µm-derived gas masses (Yue et al. 2021;Izumi et al. 2021;Neeleman et al. 2021;Walter et al. 2022).However, it appears to be in tension with the mean stellar masses of ≥ 10 10 M ⊙ recently inferred for six z ≳ 6 QSOs by modelling the spectral energy distribution of their JWST photometry (Yue et al. 2023).The only QSO host with both a dynamical mass (from Neeleman et al. 2021) and stellar mass (from their study) is J0100+2802, for which the dynamical mass is 1.2 +0.6 −0.7 × 10 10 , the CO(2-1)-derived molecular gas mass is (0.95 ± 0.32) × 10 10 and the stellar mass derived by (Yue et al. 2023) is < 30 × 10 10 .For this source, which is in both our samples, the upper limit could still yield a consistent mass budget.Better photometry are needed for the full sample of QSO hosts studied here to accurately characterise their stellar masses and compare these to the other baryonic components.
The above argument that the molecular gas masses dominate the mass budget, relies on a big 'if'; in reality, the measured dynamical masses are highly uncertain and could be systematically underestimated if the measured sizes and/or circular velocities are too small, or if the assumed geometry is incorrect (see discussion in Neeleman et al. 2021).If the gas is not in a planar stucture or if the dynamics are not supported by gravity (but instead by large outflows or turbulence), then the models assumed in fitting the dynamics are inappropriate.The dynamical masses may also be underestimated if these z ≳ 6 QSO host systems are not in dynamical equilibrium.Indeed, all evidence points to these systems being in a state of rapid evolution, implying that the assumptions used in the modelling are not appropriate for these systems and leaving room for higher α CO factors.
An additional complication for the dynamical versus gas mass comparison, is that current [C ii] 158 µm observations may be missing the much fainter, extended component, as shown Novak et al. (2020) based on the stacked [C ii] 158 µm emission of 20 z ≳ 6 QSO hosts (including those studied here).The much larger extent of this faint component would mean that the measured rotation curves do not probe far enough into the galaxy 'disk' to sample the maximum rotation velocity, nor would the measured sizes then represent the whole disk.This is why dynamical masses are quoted within certain radii.In the case of the QSO hosts studied here, we quote values inferred from within twice the [C ii] 158 µm scale length (2r d, [CII] ).In order for the molecular gas traced by CO(2-1) to make up a small fraction of the total mass within this region, it would have to be significantly more extended than the [C ii] 158 µm emission.Although we cannot test this scenario at the resolution of our data, it seems unlikely based on low-redshift galaxies (e.g. de Blok et al. 2016;Bigiel et al. 2020).Accurately testing how much of the total mass is in the form of molecular gas, will require deep, highresolution observations of the host galaxy (and its environment), as well as better sampling of the CO excitation ladders, enabling robust constraints on α CO .

Depletion timescales
Combining the CO(2-1)-derived molecular gas masses and IRinferred SFRs, we find that these z ≳ 6 QSO hosts have extremely short gas depletion timescales.Our highest-redshift target, J2348-3054, has a high molecular gas mass of ∼ 1 × 10 10 M ⊙ .Assuming that 25-60% of its IR luminosity is attributable to star formation (as for the QSO sample of Leipski et al. 2013) and following the SFR-IR calibration of Hao et al. (2011), this equates to a SFR of 500 − 1300 M ⊙ yr −1 .Combined with the CO(2-1)-derived molecular gas mass, this implies a short depletion timescale (t depl = M mol /SFR) of only 8 − 20 Myr.This value is consistent with the full literature sample, which span SFRs of 400 − 7000 M ⊙ yr −1 and depletion timescales of 4 − 50 Myr.These depletion timescales are significantly shorter than for most ULIRGS and high-redshift star-forming galaxies, which have t depl ∼ 10 7 − 10 8 Myr (e.g, Solomon & Vanden Bout 2005;Tacconi et al. 2020), but are comparable to a handful of hyperluminous starbursts at z > 2 (e.g.Combes et al. 2012;Aravena et al. 2016;Ciesla et al. 2020).These results imply that IRluminous z ≳ 6 QSO hosts are undergoing an extreme starburst phase.However, better sampling of the dust-continuum emission is still required to disentangle the AGN and SFR contributions in these systems.

CO line excitation
The fact that we do not detect the CO(2-1) emission of P036+03 and J0305-3150 implies that their average CO excitation is more extreme than that of the z = 6.4 QSO host J1148+5251, which is the best-studied CO line excitation ladder of any z ≳ 6 QSO studied to date.To investigate this further, we compare the CO line excitation ladders to those of the local AGN and QSOs investigated by Kamenetzky et al. (2016) and Rosenberg et al. (2015) in Fig. 4. Unlike for the local AGN and QSOs, many of the z ≳ 6 QSO hosts have thermalised or even suprathermal mid-to-low J line ratios (left panel) (as also shown in Riechers et al. 2009;Wang et al. 2010Wang et al. , 2011b)) Fig. 4: CO excitation ladders for the three z > 6.5 QSO host galaxies studied here (filled squares), compared to other z ≳ 6 QSO hosts with CO(2-1) measurements (semi-transparent lines connected by small filled circles) and local AGN and QSO hosts from Rosenberg et al. (2015) and Kamenetzky et al. (2016) (grey line and shading indicating the median and 16th to 84th percentiles).
The colour-coding indicates the redshift, with the legend at the right providing the name of the corresponding QSO.Left panel: CO line luminosity ratios relative to CO(2-1).Right panel: CO line luminosities in units of L ⊙ .Most z ≳ 6 QSO host galaxies show significantly (factor of ∼3×) higher mid-J versus CO(2-1) emission than the local AGN and QSO hosts (near/exceeding the thermalised value).The mid-J CO line emission for most z ≳ 6 QSO hosts exceeds that of the local AGN and QSO hosts by a greater factor than the CO(2-1) emission, implying that a greater fraction of these more massive molecular gas reservoirs is in more highly excited states.this is caused by significantly lower CO(2-1) or significantly higher mid-J CO emission we also show the CO excitation ladder in units of L ⊙ (right hand panel).From this comparison, it seems that the difference between the samples is much more pronounced at J = 6 and 7 than for J = 2, implying that there is an even greater fraction of molecular gas in these, more highly excited states in the z ≳ 6 versus z ∼ 0 QSOs.We discuss this further in Sec. 5. Based on our results, in combination with the literature z > 5.7 QSO hosts, we caution that when inferring molecular gas masses from mid-J observations of these extreme systems, it is important to apply a higher excitation factor than what is measured for local AGN and ULIRGs (r 6,1 = 0.23) or z = 2 − 3 star-forming galaxies (r 6,1 = 0.28, Boogaard et al. 2020).For the full sample of z ≳ 6 QSOs with observations of both the CO(6-5) and CO(2-1) transitions, we measure a mean and standard deviation on the CO(6-5)-to-CO(2-1) line luminosity ratio of r 6,2 = 0.9 ± 0.2.Assuming that the first three levels are thermalised, this also equates to a CO(6-5)-to-CO(1-0) line luminosity ratio of r 6,1 = 0.9 ± 0.2.

FIR versus CO line luminosities
We test whether these z ≳ 6 QSO hosts follow empirical relations between the CO line luminosities and FIR luminosities, and hence whether they are consistent with trends found for local (U)LIRGs and z = 1−4 star-forming galaxies, in Fig. 5.The FIR luminosities of the z ≳ 6 QSO hosts are corrected for CMB effects as described in Sec.3.1.However, for the CO line luminosities, we show both the measured values and those calculated using a CMB correction consistent with the high-density and temperature, non-LTE case of da Cunha et al. (2013a).To compare FIR luminosities derived for different integration widths, we apply two scaling factors.For the sample of Valentino et al. (2020), we use the provided scaling L TIR,8−1000µm = 1.2LFIR,40−400µm .We also assumed L TIR,8−1000µm = 1.45LFIR,42.5−122.5µm(De Breuck et al. 2011).We compare to three sets of L FIR versus L ′ CO relations derived by Greve et al. (2014), Liu et al. (2015) and Kamenetzky et al. (2016).We adopt the FIR wavelength range convention of L FIR,40−400µm also used in Liu et al. (2015) and scale the Greve et al. (2014) relations using the mean value of their local sample, L TIR,8−1000µm /L FIR,30−500µm = 1.77.By converting in this manner, we are introducing an additional level of uncertainty.However, the conversion is crucial to avoid systematic offsets.
We find that the measured CO(2-1) line luminosities of the z ≳ 6 QSOs are offset, on average, to 2-6× lower values that predicted from the relations with L FIR (although consistent with the scatter of lower-redshift sources), whereas the higher-J transitions are consistent with these literature relations.This minor systematic offset may be the result of various factors: 1) the FIR luminosities of QSOs may be the result of both the SFRs and a significant contribution from the AGN, such that this comparison no longer reflects the correlation between the molecular gas mass and SFR, 2) the combined effects of the heating by, and contrast against the CMB serves to decrease the observed CO(2-1) line flux more than the J = 6 or 7 transitions, or 3) previous calibrations do not accurately constrain the slope of the relation as they did not sample the extremely high L FIR regime well enough.
To test whether the dust is being significantly heated by the AGN in these QSOs (explanation 1), would require consistently well-sampled dust SEDs, with which the different dust components can be disentangled, as in Leipski et al. (2013).By modelling the contributions to the FIR luminosities of 11, z > 5 QSO hosts, they find that the AGN contribute 40-75% to the FIR luminosity.Subtracting the mean AGN contribution of their sample, shifts the z ≳ 6 QSO hosts studied here close to the L ′ CO(2−1) versus L FIR relation (systematically above but consistent with the scatter), and systematically below but within the scatter of the L ′ CO(6−5) and L ′ CO(7−6) versus L FIR relations.In Fig. 2, we also show that most of the z ∼ 1.2 AGN from Valentino et al. (2020) lie towards lower CO(2-1) line luminosities than predicted from the total FIR luminosities.However, when subtracting the AGN contribution (which they model in their SED fitting) this systematic offset also disappears.The adopted CMB correction (faint outlined squares and circles) also brings the values more in line with the CO(2-1) relation, while still keeping the high-J lines consistent with these relations.However, it may well be that no correction is needed, due to the extremely high gas densities and temperatures in these QSO hosts.Isolating the CMB effects from the heating and turbulence generated through star formation and the central QSO will require dedicated simulations of these systems.

Discussion
The high line ratios and implied low α CO of these z ≳ 6 QSO hosts indicate the presence of more extreme molecular ISM conditions than in the less-luminous, and mostly gas-poorer, local AGN and QSO hosts and ULIRGs we have compared them to.These are likely caused by a combination of both strong ISRFs (and high X-ray fluxes), and the presence of compact, dense gas reservoirs.Our qualitative argument for this can be summarised as follows.
To keep a significant component of molecular gas (without dissociating it) this gas must be extremely dense and/or well shielded.Both appear to be the case, in that the z ≳ 6 QSO hosts galaxies studied here already have a high implied gas-todust ratios and small measured sizes; with [C ii] 158 µm halflight radi of 0.8 to 4.5 kpc (Neeleman et al. 2021;Walter et al. 2022;Meyer et al. 2022).The extent of this [C ii] 158 µm emission is significantly smaller (by at least a factor of five) than the scale of the emission in the local AGN and QSO hosts we are comparing to.Compounding this, the molecular gas masses inferred via the same α CO conversion factors are a few times larger on average (or of the same order as) the local comparison sample; the molecular gas masses of the z ≳ 6 QSO hosts are ∼few ×10 10 M ⊙ versus the ∼few ×10 9 M ⊙ measured for gasrich local AGN and QSO hosts (e.g.Papadopoulos et al. 2012b;Rosenberg et al. 2015;Kamenetzky et al. 2016).Indeed, Shao et al. (2019) infer densities of a 10 5 − 10 6 cm −3 for three of the literature QSO hosts included here.Not only do the z ≳ 6 QSO hosts have more molecular gas per unit volume than the local AGN and QSO hosts, they also have an order of magnitude higher (on average) FIR luminosities than the local AGN and QSO hosts.The implied SFRs of 400 − 7000 M ⊙ yr −1 , which stem from a smaller volume than in the local galaxies, imply higher ISRF strengths (with G 0 ≳ 10 3 in Habing units for some of the z ≳ 6 QSO hosts included here, according to the modelling performed by Shao et al. 2019 andDecarli et al. 2022).In combination, this is consistent with a picture in which the molecular clouds in these z ≳ 6 QSO hosts are more closely spaced and/or much denser than in local AGN and QSO hosts and are illuminated by stronger radiation fields.
The high CO(7-6)-and CO(6-5)-to-CO(2-1) line ratios (consistent with thermalised and even supra-thermal CO) imply the presence of dense, warm, and highly turbulent gas.Similarly high line ratios have been found in the local photodissociation region (PDR) known as the Orion Bar.Modelling the far more extensive set of molecular and atomic emission lines (including isotopologues) observed for this region, several groups (Goicoechea et al. 2015(Goicoechea et al. , 2016(Goicoechea et al. , 2019;;Andree-Labsch et al. 2017;Parikka et al. 2018) have found evidence for clumps of extremely high-density (n H ∼ 10 5 − 10 6 cm −3 ), and high-temperature (T ∼ 150 − 300s K) molecular gas.To produce such high line ratios, the gas temperatures may exceed the average dust temperature (see also Krumholz 2014;Harrington et al. 2021).Several sources of heating could elevate the kinetic temperatures of the gas.High star formation rate densities should give rise to more cosmic rays (compared to local ULIRGS/AGN), providing an additional, strong source of heating that increases the gas temperatures and causes the CO excitation ladders to peak at higher-J (e.g.Papadopoulos et al. 2012b;Bisbas et al. 2023).If large amounts of energy (e.g. from star formation and/or the QSO) are being injected into such a small volume of gas, then this gas is likely also more turbulent, providing another source of heating and implying a lower intrinsic α CO (e.g.Narayanan et al. 2011;Papadopoulos et al. 2012a).In combination with the high densities, this increase in velocity dispersion would serve to counteract any increase in α CO caused by a higher ISRF or higher cosmic ray ionisation rate (e.g.Clark & Glover 2015).
The high carbon abundances implied by our results are also consistent with the presence of highly turbulent gas.Mixing caused by strong turbulence may result in elevated [C i]/H 2 in the cloud interiors, thereby giving rise to bright [C i] emission (Xie et al. 1995;Glover et al. 2015).Moreover, strong contributions from cosmic rays, X-rays and/or shocks (Israel et al. 2014(Israel et al. , 2015) ) would enhance the brightness of [C i] compared to low-J CO emission (Israel et al. 2014(Israel et al. , 2015;;Papadopoulos et al. 2018Papadopoulos et al. , 2022)).These additional heating sources may drive atomic carbon out of the state of LTE asumed here.Indeed, Harrington et al. (2021) and Papadopoulos et al. (2022) show that the [C i] excitation of high-redshift star-forming galaxies and local (U)LIRGS is best described by non-LTE radiative transfer models; with the LTE assumption easily accounting for a factor of a few too high M [CI] .In that case, a lower atomic carbon abundance than what we assumed may still be appropriate.Testing the impact of turbulence, cosmic rays, and X-rays will require better-sampled CO excitation ladders, [C i] transitions (including the addition of the [C i](1-0) line) and dust spectral energy distributions than are currently available for these z ≳ 6 QSO hosts.
Given our data, it is unclear to what extent the extreme conditions within the molecular gas of these z ≳ 6 QSO host galaxies are driven by star formation versus the luminous accreting black holes.In Fig. 6, we show that the higher line ratios of the z ≳ 6 versus z ∼ 0 QSO (and AGN) hosts correlate with proxies for both the AGN power and star formation activity.The left column illustrates that such luminous QSO hosts have not been observed in CO in the local Universe; instead the most powerful AGN in the local Universe are hosted in elliptical galaxies with little molecular gas.However, we compare to these z ∼ 0 systems to span a wider range of bolometric and FIR luminosities.We use the bolometric luminosity L bol (middle column) as a proxy of the AGN power2 and the FIR luminosity, L FIR (right-most column) as a proxy for the SFR.The fact that these three line ratios correlate with the FIR luminosity is expected, e.g.high star formation rate densities lead to greater dust heating (and hence higher L FIR ) and also higher CO excitation (e.g.Narayanan & Krumholz 2014).We also find a clear correlation with the bolometric luminosity, albeit with a large scatter.This scatter can be attributed to the different, indirect methods used to infer the L bol : scaling the 3000Å luminosity following Shen et al. (2011) 3 for the z ≳ 6 QSO hosts, scaling the luminosity at 14-150 keV following Ricci et al. (2017) for the local sample (BASS catalogue of Koss et al. 2022) and scaling the [O iii]λ5007Å line luminosity4 for the sample of Molyneux et al. (2023).Thus, neither the tightness nor slope of these correlations allow us to disentangle the impact of star formation versus AGN activity on the molecular gas.
The lack of high-J CO emission and inhomogeneous X-ray coverage also do not allow us to disentangle the impact of star formation versus AGN activity.These QSOs may be heavily obscured and have inhomogeneous X-ray coverage.Two of our sample (J0305-3150 and J2348-3054) have only upper limits on the X-ray fluxes from Chandra, with net counts at 0.5-7 keV of ≤ 15 and ≤ 21 (Wang et al. 2021), whereas P036+03 and two of the literature sample, J0100+2802 and J1429+5447, are detected with XMM Newton (Medvedev et al. 2021;Zappacosta et al. 2023).For this sample of five, the data reveal no correlation between the CO line ratios probed here and the X-ray luminosity (e.g. at 2-10 keV).Future X-ray observations and observations of high-J CO transitions, up to CO(11-10), will enable tests against predicted trends with the X-ray flux (e.g.Vallini et al. 2019).
Overall, we find little evidence for any loss of signal due to the contrast against the CMB.The lack of a large systematic offset from empirical relations between the CO line versus CMBcorrected FIR luminosities indicates that little correction for the CMB contrast is needed for the CO lines.Likewise, the low α CO value needed to ensure the molecular gas masses do not exceed the dynamical masses also indicates that minimal correction for the CMB contrast is needed.However, the CMB can have multiple effects, which we cannot disentangle from variations in the ISRF and/or ISM density: it is a stronger background against which the lines are observed, but it is also an additional source of heating, thereby changing the temperature and density structure of the molecular ISM.As shown by da Cunha et al. (2013a), heating by the CMB can serve to make the CO excitation ladders peak at higher J, similar to effects of heating through strong ISRFs, cosmic rays, or X-rays.To disentangle these effects it is crucial to observe CO, [C i], and dust-continuum emission from samples of QSOs across different epochs, with the same high IR luminosities and molecular gas masses.
We check whether the compactness (proxy for density) impacts the molecular gas conditions as traced by the relative line strengths.To this end, high-resolution observations of the [C ii] 158 µm and FIR emission are available for our set of QSOs (P036+03, J0305-3150 and J2348-3054) and four of the literature sample (J0129-0035, J2054-0005, J0100+2802, and J1148+5251) (Venemans et al. 2020;Neeleman et al. 2021;Meyer et al. 2022;Tripodi et al. 2023).Measuring sizes has so far proven challenging, with different assumed surface brightness profiles and modelling tools being applied.Nonetheless, by any metric, J2348-3054 is the most compact of the seven galaxies (r eff,[CII] ∼ 0.46 kpc).J0129-0035 and J2054-0005 are the next most compact (r eff,[CII] ∼ 0.8 kpc).We find marginal evidence that the compactness is related to the CO excitation.Whereas J2348-3054 and J0129-0035 exhibit some of the highest CO(6-5)-to-CO(2-1) line ratios, J1148+5251, the most extended of these resolved z ≳ 6 QSOs, has line ratios consistent with the sample mean (neither at the low nor high end).The fact that the two most compact sources seem to have the highest excitation is consistent with the picture of a denser and warmer ISM in galaxies with higher UV flux per unit volume (originating either from the central black hole or star formation).However, the sample size with both multi-line observations and resolved [C ii] 158 µm emission remains too small to draw any strong conclusions.
We also attempt to test whether the presence of companions (and hence interactions) affects the relative line strengths or molecular gas masses.J2348-3054 and J2054-0005 have no confirmed companions, whereas J0305-3150 and J0100+2802 do (Venemans et al. 2019;Tripodi et al. 2023;Wang et al. 2023).J1429+5447 is a known merger (e.g.Khusanova et al. 2022).Of these, J0305-3150 has the highest mid-to-low-J CO line ratios, considering the upper limit on CO(2-1).J0100+2802 and J1429+5447 are consistent with the mean excitation of the full sample.Within this small sample, we find no clear connection yet between the presence of companions and the molecular gas masses or excitation.However, future JWST observations, as in Wang et al. (2023), will help shed light on this issue.

Conclusions
In this paper, we present the VLA Ka Band observations of CO(2-1) in three QSO host galaxies at 6.5 < z < 6.9, thereby probing their cold molecular gas reservoirs.These three QSO hosts were already detected in [C ii] 158 µm emission and the underlying FIR continuum (Bañados et al. 2015;Venemans et al. 2016) and in CO(6-5), CO(7-6) and [C i] (2-1) emission (Venemans et al. 2017a).For the highest-redshift QSO host, J2348-3054 at z = 6.9, we robustly detected the CO(2-1) emission, revealing a molecular gas reservoir of ∼ 1 × 10 10 M ⊙ already 780 Myr after the Big Bang.However, for the other two QSO hosts, P036+03 and J0305-3150, the CO(2-1) observations are consistent with noise, yielding 4σ upper limits on the molecular gas masses of (1.0 − 1.5) × 10 10 M ⊙ .We detect no underlying continuum emission for these three targets, nor any CO or continuum emission from their known companions.We combine these three QSO hosts with the full set of nine z > 5 QSO host galaxies that already have CO(2-1) observations (Table A.1), thereby performing a broader study of trends in the molecular gas masses and CO excitation.To this end, we consistently re-derive the FIR luminosities and dust masses, taking CMB effects into account.Based on the CO(2-1) observations, and the assumptions of a ULIRG α CO = 0.8 M ⊙ (K km s −1 pc 2 ) −1 (see Carilli & Walter 2013) and thermalised low-J transitions (r 2,1 = 1), the molecular gas masses of the three QSO hosts studied here are consistent with the literature sample, which spans (0.8 − 5) × 10 10 M ⊙ .These values are mostly consistent with those inferred from the dust masses, assuming a gas-to-dust ratio of 100.However, the molecular gas masses estimated from the [C ii] 158 µm line luminosity using the mean α [C II] of Zanella et al. (2018) are a factor of 6.8 ± 1.9 higher than the CO(2-1)-based values, reflecting a potential decrease in α [C II] with sSFR.In order for these CO(2-1)-derived values to also be consistent with those inferred from [C i] (2-1), atomic carbon must be already be abundant ([C i]/H 2 ∼ 9 × 10 −5 ) and/or be undergoing significant heating through cosmic rays, X-rays, turbulence, and/or shocks.
For P036+03 and J0305-3150, the gas excitation is significantly higher than for the z = 6.4 QSO J1148+5251 (Riechers et al. 2009;Stefan et al. 2015), implying more extreme physical conditions in these systems.Despite the small range of molecular gas masses, the combined set of z ≳ 6 QSO hosts spans a wide range in CO excitation.Nonetheless, the line ratios studied here are systematically higher than for local AGN and QSO host galaxies with observations of the same transitions; the mean CO(6-5)-to-CO(2-1) line luminosity ratio is r 6,2 = 0.9 ± 0.2 (2.5× higher than in local AGN and QSOs) and the mean [C ii] 158 µm-to-CO(2-1) line luminosity ratios is 1.0 ± 0.4 (2× higher than in local AGN and QSOs).Assuming that the CO(2-1) emission is thermalised, the mean CO(6-5)-to-CO(1-0) line luminosity ratio is also r 6,1 = 0.9 ± 0.2.These higher line ratios should be considered when inferring the molecular gas masses of z > 5 QSO hosts from only their mid-J CO lines.
Comparing the CO(2-1) observations to the existing dynamical mass measurements, we find that the mass budget of the three QSOs studied here is likely dominated by the molecular gas.In particular, J2348-3054 has a molecular gas mass consistent with the lower limit on the dynamical mass, leaving little mass in neutral gas or stars within the [C ii] 158 µm-emitting region.If the dynamical masses based on the resolved [C ii] observations are correct, these CO(2-1) observations imply that the CO-to-H 2 conversion factor must be low-similar to the ULIRG value-to yield molecular gas masses that are lower than the dynamical masses.Thus, little correction for the contrast against the CMB is needed for these bright, compact systems.Combined with the high CO excitation, these results imply the presence of warm and dense molecular gas reservoirs.In combination with the literature sample, we find that the CO(2-1)-derived molecular gas masses are 20 − 190× higher than the dust masses, implying the presence of significant, enriched reservoirs of dense molecular gas and dust.Moreover, the CO(2-1)-derived molecular gas masses are 3 − 11× higher than the black hole masses.
Our new CO(2-1) observations support a picture in which these IR-luminous z ≳ 6 QSO host galaxies have rapidly assembled copious amounts of molecular gas (either through accretion or mergers).Combining the previously detected CO(6-5), CO(7-6), and [C i] (2-1) emission with the new CO(2-1) constraints, indicates that all three QSO hosts studied here harbour significant amounts of warm and dense, metal-enriched gasindicating that these are already highly evolved galaxies.In particular, the newly CO(2-1)-detected z = 6.9 QSO host J2348-3054, has properties consistent with a massive starburst that is rapidly coevolving with the luminous central AGN (as found previously for some of the literature sample).
With these highest-redshift CO(2-1) observations, we confirm that it is possible to detect the cold, molecular gas traced by CO(2-1) emission at z > 6.5 within ∼ 15 hours of VLA time.These observations are also now feasible with ALMA Band 1 (although it remains less efficient than the VLA at z > 5).Observing CO(3-2) is also now feasible up z = 8.8 with ALMA's Band 1, whereas [C i](1-0) can soon be observed with ALMA's Band 2, up to z = 6.3.Complementing the existing plentiful [C ii] 158 µm observations of z ≳ 6 QSO hosts with observations of the cold, diffuse molecular gas is therefore achievable, albeit time-intensive.Next-generation facilities, such as the Next-Generation VLA (Eric Murphy and the ngVLA Science Advisory Council 2018), will enable such observations to be conducted even more efficiently.The next step is to systematically model the average conditions within these massive molecular gas reservoirs-constraining the turbulence within the molecular gas, the carbon abundance and α CO .To achieve this, it is now crucial to obtain QSO samples with a homogeneous coverage of CO (up to J > 8), [C i], and dust emission.The compilation of z ≳ 6 QSO hosts studied here will serve as a crucial basis for this work.

Fig. 1 :
Fig.1: CO(2-1) observations of P036+03, J0305-3150, and J2348-3054 (from top to bottom, respectively).Left: Cleaned channel maps with the middle panel depicting the moment-0 map centred on the [C ii] 158 µm redshift and integrated over 1.2× full width at half maximum (FWHM) of the [C ii] 158 µm line, which should encompass 84% of the CO(2-1) emission assuming there is no spatial or velocity offset.Right: Observed CO(2-1) spectra extracted from the pixel at the peak position of the moment 0 maps (blue line) and 1σ rms per channel (orange line).The spectra are shown for channels of 100 km s −1 velocity width, the shaded region indicating the 1.2 FWHM region used to create the channel maps to the left.P036+03 and J0305-3150 are undetected (with a significant number of negative and positive peaks in the data cube exceeding the significance of the peaks in emission near the expected source centres), so we place 4σ upper limits on their CO(2-1) emission.J2348-3054 is detected at 5.1σ (central velocityintegrated map in bottom row).
rely on two literature scaling relations between the depletion time and sSFR-those of Sargent et al. (2014) and Scoville et al. (2017)comparing the calibration derived from the average of the two, and from Scoville et al. ( Fig. 5: FIR luminosities versus CO line luminosities of our three QSOs (coloured squares), the literature sample of z ≳ 6 QSOs (filled, coloured circles), the z > 3 SFGs in Table A.2 (orange circles), the z ∼ 1.2 star-forming galaxies and AGN from Valentino et al. (2020) (orange and blue triangles), and the local (U)LIRGS, AGN, and QSOs from Rosenberg et al. (2015) and Kamenetzky et al. (2016) (grey, right-pointing triangles).The z ≳ 6 QSOs follow the literature relations for the CO(6-5) and CO(7-6) transitions but are shifted to 2-6× lower CO(2-1) line luminosities than predicted by these relations.

Fig. 6 :
Fig.6: Line luminosity ratios versus redshift (left column), L bol (middle), and L FIR,40−400µ m (right).The data for the three QSO hosts studied here are depicted via filled squares, whereas the previously studied z ≳ 6 QSO hosts are depicted by filled circles.The colour-coding is the same as for Fig.4.We also show the local AGN and QSO hosts fromKamenetzky et al. (2016) andMolyneux et al. (2023) (grey and black, respectively) and a few z = 1 − 4 QSO hosts fromRiechers et al. (2008Riechers et al. ( , 2009));Aravena et al. (2008), corrected for lensing where appropriate (green circles).The bolometric luminosities of the different samples were inferred in different ways: scaling the 3000Å luminosity followingShen et al. (2011) for the z ≳ 6 QSO hosts, scaling the luminosity at 14-150 keV followingRicci et al. (2017) for the local sample (BASS catalogue ofKoss et al. 2022) and scaling the [O iii]λ5007Å line luminosity for the sample ofMolyneux et al. (2023).For the z = 1 − 4 QSO hosts we use values determined from the [O iii]λ5007Å line, 6 µm emission, or the full spectral fitting.The FIR luminosities of the z < 4 samples were taken from the same studies, based on their SED fitting, whereas for the z ≳ 6 QSOs, we use the values derived in this work.

Fig. B. 1 :
Fig. B.1: Outputs of the Bayesian dust SED fitting for our three QSO hosts.For each source (labelled at top), the top right panel shows the best-fit dust models (solid lines) and the observed fluxes (blue and purple points).The data that were included in the fits are shown in light blue, those not included in purple.The remaining panels from top left to bottom right show the marginalised likelihood distributions for the three fit parameters: dust temperature (T dust ), mass (M dust ), and emissivity index (β).

Table 1 :
QSO properties and data.