Introduction

Quasars residing in the so-called cosmic dawn (z 6, age of the Universe <1 Gyr) were first discovered at the turn of the century (e.g., Fan et al. 2000Fan et al. , 2003)). To date, we know ∼400 quasars at z > 5.5, including three at z > 7.5 (Bañados et al. 2018;Yang et al. 2020;Wang et al. 2021). Multiwavelength campaigns, in particular at submillimeter and millimeter (hereafter, mm) wavelengths, revealed that these early quasars reside in extremely active galaxies, experiencing intense starbursts (with star formation rates SFR = 100-1000 M yr -1 ; Bertoldi et al. 2003;Wang et al. 2008a,b;Drouart et al. 2014;Leipski et al. 2014;Venemans et al. 2018) and fast, radiatively efficient black hole accretion ( ṀBH > 10 M yr -foot_0 ; De Rosa et al. 2014;Schindler et al. 2020;Yang et al. 2021). This rapid gas consumption depletes the immense reservoirs of molecular gas (M H2 > 10 10 M ; Walter et al. 2003;Wang et al. 2010;Venemans et al. 2017a;Decarli et al. 2018) of the host galaxies, which can be refueled via accretion of cool gas from the circumgalactic medium (e.g., Farina et al. 2019;Drake et al. 2019) or via mergers with gas-rich galaxies (e.g., Trakhtenbrot et al. 2017;Decarli et al. 2017Decarli et al. , 2019;;Vito et al. 2019).

Both star formation and black hole accretion release huge amounts of energy which could affect the evolution of the host galaxies: Mechanical feedback (via winds or jets) can displace and even remove the gas reservoir; and radiative feedback can heat-up and excite the interstellar medium (ISM) thus slowing down or preventing cooling and gas fragmentation (see Fabian 2012;Heckman & Best 2014;King & Pounds 2015;Somerville & Davé 2015, for reviews on this topic). A quantitative description of the cold (T < 100 K) ISM, poised to constrain the mass, density, temperature, and excitation mechanism of the molecular gas in the quasar hosts, would shed light on the interplay between the ISM, star formation and nuclear activity.

In distant galaxies, the molecular gas mass, M H2 , is usually inferred based on photometry covering the Rayleigh-Jeans part of the dust continuum emission (see, e.g., Groves et al. 2015;Scoville et al. 2017a), or via far-infrared emission lines of species such as the carbon monoxide, 12 C 16 O (hereafter, CO; see, e.g., Walter et al. 2003;Wang et al. 2011, as well as the review in Carilli & Walter 2013), neutral carbon (in particular the fine-structure transitions [C i] 1-0 and [C i] 2-1 ; see, e.g., Popping et al. 2017;Valentino et al. 2018;Boogaard et al. 2020), or singly ionized carbon (via the bright [C ii] 3/2-1/2 line; see Venemans et al. 2017b;Zanella et al. 2018). All of these methods rely on a set of assumptions and conversion factors: for example, the gas-to-dust ratio, δ g/d (see, e.g., Berta et al. 2016) or the α CO factor that links the luminosity of the CO(1-0) ground transition with the associated molecular gas mass (see Bolatto et al. 2013, for a review). These factors carry implicit dependencies on the metallicity and other properties of the galaxy. Comparisons of M H2 estimates obtained via different prescriptions are pivotal to pin down systematics in the mass determinations. For instance, Dunne et al. (2021) used a combination of dust continuum, CO(1-0) and [C i] 1-0 line emission to gauge the molecular gas mass in a sample of massive galaxies at z = 0.35. Because four unknown quantities (molecular gas mass, carbon abundance, α CO , δ g/d ) are constrained with three observables, only relative values can be directly measured; however, by using additional information (e.g., from their sample variance), mean values for the calibration factors can be inferred. Similarly, Sommovigo et al. (2021Sommovigo et al. ( , 2022) ) proposed to use a combination of nested scaling relations to constrain the dust temperature, T d , from the [C ii] line luminosity.

The combination of different diagnostics of the cold ISM also reveal some precious information on the physical conditions of the gas. For instance, models of Photo-Dissociation Regions (PDRs) or X-ray Dominated Regions (XDRs) account for various chemical and radiative processes that take place in a gas cloud that is illuminated by an external radiation field in order to predict the emerging intensity of various lines. The former set of models assume that the excitation is driven by the UV radiation from young stars, and they are thus best suited to describe the emission of molecular clouds invested by the radiative output of an on-going starburst. The latter models assume the presence of a strong X-ray radiation field as the driver of the cold gas heating and excitation. These models are hence best suited to account for the quasar contribution to the excitation budget, although both PDR and XDR models have been extensively used in a wide range of astrophysical applications (e.g., Tielens & Hollenbach 1985;Maloney et al. 1996;Meijerink et al. 2007;Goicoechea et al. 2019;Vallini et al. 2019). By exploiting differences in the predicted line emission in PDRs versus XDRs, we can infer hints on the role of different radiative mechanisms (star formation or nuclear activity) that may drive the excitation of the cold ISM. Dust, [C ii], [C i], and CO(7-6) form an observationally convenient suite of diagnostics, as with only two frequency settings (one covering the [C ii], and one covering [C i] and CO(7-6), together with their respective underlying continuum) we can secure a rich collection of tracers of the cold interstellar medium.

In this study, we present the first systematic campaign targeting CO(7-6) and [C i] in a sample of ten z ∼ 6 quasars. In addition, the unprecedented bandwidth capabilities of new mm receivers together with upgraded correlators enable to observe CO(6-5) as well in the same frequency setting as CO(7-6) and

for certain redshift windows. We used the Plateau de Bure Interferometer (PdBI), later upgraded to the NOrthern Extended Millimeter Array (NOEMA) of the Institut de Radioastronomie Millimétrique (IRAM) to secure CO(6-5), CO(7-6) and [C i] coverage in the 3mm atmospheric window. All of these quasars have been previously detected in the [C ii] and dust continuum emission. The paper is structured as follows: in Sect. 2, we introduce the sample, the observations, and the data reduction. In Sect. 3, we present our results, our estimates of M H2 based on various methods, and our constraints on the physical properties of the ISM. We draw our conclusions in Sect. 4. Throughout this paper we assume a ΛCDM cosmology, with H 0 = 70 km s -1 Mpc -1 , Ω m = 0.3, and Ω Λ = 0.7, and a Kroupa initial mass function to compute star formation rates (SFR) 1 . In this cosmological framework, the scale distance at z = 6.00 is 5.713 kpc arcsec -1 , the luminosity distance is D L = 57.742 Gpc, and the age of the Universe is 914 Myr. The notation ν refers to rest-frame frequencies.

Observations and data processing

Sample

Our goal is to capitalize on the combination of

CO, and the underlying dust continuum emission to study the molecular medium in quasar host galaxies at the dawn of cosmic time. (2019). We removed all of the sources previously observed in CO(7-6) from our selection. We also removed 12 additional quasars that reside at 5.9 z 6.1, as at these redshifts the CO(7-6) + [C i] lines are significantly affected by poorer atmospheric transmission. We also removed six quasars that have IR luminosity L IR < 10 12 L , that is, quasars are likely too faint at IR wavelengths for our survey. Finally, we drop 11 quasars that lie at declination < -10 • for the sake of observability from NOEMA. In addition to this sample, we include previously unpublished observations of the quasars SDSS J0338+0021 (z = 5.0267) and SDSS J2054-0005 (z = 6.0391). The resulting list of targets thus consists of ten quasars (see Table 1).

Observations and data reduction

The dataset consists of the collection of a main survey and a number of smaller programs conducted at NOEMA and previously at PdBI, for a total observing time of ∼70 h on source (eight antennas equivalent). The half power beam width of NOEMA's 15 m antennas is ≈48.6 at 100 GHz. Observations were mostly carried out under average or poor weather conditions with precipitable water vapor columns typically 5 mm and system temperatures of 80-200 K with the array in compact (C or D) configuration. We processed the data using clic from the GILDAS 2 suite. Table 2 provides a summary of our observations. J0338+0021 was observed on August 1, and September 28, 2008 (program ID: S052), using the narrow-band 2 mm receiver, and again in July 2013 (program ID: X04D) using WideX. The sources 3C454.3 and MWC349 were observed for bandpass and flux calibration. We observed J0406+066 as amplitude and phase calibrator. The final cube includes 7200 visibilities, corresponding to 9.0 h (five-antenna equivalent). J2054-0005 was observed in October-November 2013 (program ID: X04D). We used the WideX correlator and tuned the receiver to secure the [C i] and CO(7-6) lines in a single setup.

The sources 3C454.3 and MWC349 were observed for bandpass and flux calibration. For the amplitude and phase calibration, we observed 2059+034. The final cube includes 16881 visibilities, corresponding to 21.10 h (five-antenna equivalent).

PJ036+03 was observed in multiple short tracks in June-July 2015 (program ID: S15DA) and June-July 2017 (program ID: S17CD), in compact array configuration (5D and 6D in 2015; 7D and 8D in 2017) using the WideX correlator. The receiver was tuned to encompass the two CO(6-5) and CO(7-6) lines in two independent frequency settings. The calibration sources include the blazar 3C454.3 for bandpass calibration, the source MWC349 for absolute flux calibration, and the quasar 0215+015 for amplitude and phase calibration. The final cubes consist of a total of 22440 visibilities (13.33 h on source, assuming seven antennas) on the CO(7-6) line, and 22317 visibilities (13.28 h on source, assuming seven antennas) on CO(6-5).

The quasar PJ338+29 was observed in July 2018 (program ID: S18DM). The expanded bandwidth offered by the Poly-Fix correlator enables to encompass the CO(6-5) and CO(7-6) lines in our targets using a single frequency setting. The blazar 3C454.3 served as bandpass calibrator, MWC349 as absolute flux calibrator, and the quasars 2234+282 and J2217+243 as amplitude and phase calibrators. The final cubes comprises a total of 8670 visibilities (3.87 h on source, eight-antenna equivalent).

We observed the quasar J2219+0102 as part of two programs, one covering the CO(5-4) and CO(6-5) lines (August-September 2018, program ID: S18DM), and one covering the CO(6-5), CO(7-6) and [C i] lines (June 2019, program ID: S19DM). The bright sources 3C454.3 and MWC349 acted as bandpass and flux calibrators, while the quasars 2216-038 and 2223-052 were observed as amplitude and phase calibrators. We reduced the two settings independently, and then merged the visibilities in the overlapping frequency range. The 2018 and 2019 observations comprise 13380 and 8640 visibilities, respectively, corresponding to 5.97 and 3.86 h on source (eight-antenna equivalent).

We observed the remaining targets between June 2019 and January 2020 (program ID: S19DM). The observations of PJ159-02 relied on 0851+202, LKHA101, and J1028-0236 as bandpass, flux, and phase and amplitude calibrators respectively. The final cube consists of 7170 visibilities (2.0 h on source, ten-antenna equivalent). Observations of J1319+0950 were performed in excellent weather conditions (precipitable water vapor <2 mm, system temperature T sys = 55-70 K) and used on the calibrators 3C273 (bandpass), MWC349 (flux), and 1307+121 (phase, amplitude, pointing, focus). The final cubes consists of 8310 visibilities (2.88 h on source, nine-antenna equivalent). Observations of the quasar J1048-0109 were performed using 3C84, 3C273, and LKHA101 as bandpass and flux calibrators, and 1055+018 as phase and amplitude calibrator. The final cube consists of 10170 visibilities (3.53 h on source, nine-antenna equivalent). We calibrated the observations of the quasar J1104+2134 using MWC349, 0851+202, 3C273, 0923+392 for flux and bandpass calibration; and 1040+244 for phase and amplitude. The final cube consists of 12660 visibilities (4.40 h on source, nine-antenna equivalent). Finally, the observations of PJ359-06 capitalized on 3C454.3, MWC349, and 0003-066 for bandpass, flux, and amplitude and phase cali-bration. The final cube comprises 3.46 h on source (nine-antenna equivalent).

Imaging

We imaged the cubes using the GILDAS suite mapping. Natural weighting was adopted in all the cases. The beam FWHM ranges between 3 and 7 (see Table 2). We resampled the spectral axis in 50 km s -1 wide channels (computed at the CO(7-6) frequency). In creating the continuum images, we first masked out the channels where strong emission lines are found, based on visual inspection of the spectra. When no obvious lines were detected, we masked out a range of ∼600 km s -1 around the expected frequency of the lines. We created continuum images from the lower and upper side bands independently using the unmasked channels, and use them to subtract from the cubes using the task uv_subtract. Finally, we created continuumsubtracted line maps using the channels within the line full width at zero intensity. In the absence of an obvious line detection, we used a 6 channels (300 km s -1 ) window around the expected frequency of the transition, based on the [C ii] redshift. For all the maps, we defined cleaning boxes around the quasar, and cleaned down to 1.5-σ per channel using Hogbom cleaning (via the task clean). Finally, we extracted spectra in single-pixel extraction at the formal optical and near-infrared coordinates of the targets (i.e., we did not optimize the signal extraction based on the observed maps). All astrometric uncertainties are negligible with respect to the beam size of the presented observations. Figures 1 and2 show the extracted spectra and the maps of all the sources in our sample.

Results

In the following, we present our results on the spectral analysis of the new CO(7-6) + [C i] and CO(6-5) observations (Sect. 3.1). We then infer molecular gas mass estimates (Sect. 3.2), and constrain the properties of the ISM based on the available lines (Sect. 3.3).

Spectral analysis

We fit the 3 mm spectra from our sample and the 1 mm photometric measurements from the literature (see Table 1) with a Gaussian curve for each line: CO(6-5), CO(7-6), and [C i]; and a dust continuum described as:

where Ω is the source angular size in steradiants; z is the source redshift; B ν (T ) is the black body emissivity spectrum; T d is the dust temperature; T CMB is the CMB temperature at the source redshift, T CMB = 2.725 (1 + z) K. Finally, τ ν is the optical depth:

where κ d (ν) is the frequency-dependent dust absorption coefficient; Σ d is the dust surface density; M d is the total dust mass within the angular size defined by Ω; D A is the angular diameter distance, D A = D L /(1 + z) 2 ; N H2 is the column density of molecular gas; m H2 is the mass of the H 2 molecule; and δ g/d is the gasto-dust mass ratio. We interpolate the κ d (ν) values reported in Table 5 of Draine (2003) throughout the IR regime, and extrapolate assuming a power-law with slope β = 1.9 at ν < 350 GHz (see also da Cunha et al. 2021). At the CO(7-6) frequency, Eq. ( 2) implies that τ ν ∼ 1 for a molecular gas column density of ∼1.8 × 10 25 cm -2 (assuming δ g/d = 100, see, e.g., Berta et al. 2016). Such high column densities have been observed in the nuclei of some ULIRGs (e.g., Arp 220; see Scoville et al. 2017b), but are not expected to apply to the galaxy-averaged regimes sampled in our spatially unresolved observations (see also Sect. 3.3). At [C ii] frequencies, the same optical depth is reached at ≈6.3× lower column densities. We perform the spectral fits using our custom Monte Carlo Markov chain algorithm smc. We fit the redshift of the quasars, informed by the [C ii] redshift; the width of the emission lines; the integrated fluxes of the lines; the dust continuum temperature, T d ; and the dust mass within the beam, M d . Due to the modest signal-to-noise ratio of some of the lines in our survey, we impose that all of the fitted CO and [C i] lines of a given source have the same velocity width. As priors, we adopt Maxwellian distributions for the line width and the dust temperature (using 300 km s -1 and 50 K as scale lengths of the distributions); and loose Gaussian priors for the line fluxes and the logarithm of the dust mass. In Eqs. ( 1) and (2), we assume Ω = Ω b f , where Ω b is the resolution element in steradiants, and f is the (unconstrained) filling factor. Operatively, we assume f = 1 in the fits, and stress that, in the optically thin scenario, Ω has no impact on S ν and M d as it disappears from Eq. ( 1). The best fit and 1-σ errors on the fitted parameters are derived from the 50%, 14%, and 86% quantiles of the marginalized posterior distributions. The results of the fits are listed in Table 3. In the following, we treat as upper limits measurements of line fluxes and continuum flux densities for which the posterior flux estimates are below their upper 3-σ confidence level. All of the quasars are detected in the CO(7-6) emission line. Conversely, only 6/10 and 6/8 of the targeted quasars are detected in their [C i] and CO(6-5) emission, respectively.

In Fig. 3, we plot the continuum flux density at 1 mm (at the rest frequency of [C ii], as reported in the literature; see Table 1 for references) as a function of the observed flux densities at 3 mm [at the reference frequency of CO(7-6), from this work]. We find a typical 1 mm/3 mm flux density ratio of ∼19. The observed ratios favor relatively low optical depth values, τ 1900 GHz 1 against the optically thick scenario (τ 1900 GHz 5) for the quasars in our sample, unless extremely high temperatures (T d 100 K) are invoked. This suggests that the bulk of the dust is optically thin at 3 mm in these galaxy-wide observations, in contrast with what has been observed at the very center of some quasar host galaxies (e.g., Walter et al. 2022).

For most of the quasars in our sample, only the two photometric points shown in Fig. 3 are available in the rest-frame far-infrared, thus making our estimates of the IR luminosities quite uncertain. In particular, the lack of constraints on the dust spectral energy distribution at wavelengths shorter than the expected peak hinders any informative constraints on the dust temperature, T d . For instance, the observed F ν (1 mm)/F ν (3 mm) ratios in our sample are equally fit with (T d ,τ 1900 GHz ) = (60 K, 1) or (47 K, 0.2). With this caveat in mind, in the remainder of the analysis we infer IR luminosities following Eq. ( 2) under the assumption of τ 1900 GHz = 0.2 and T d = 47 K, scaled to match the observed 1 mm flux density. These values are in line with the dust temperature reported for other high-z quasars that have been extensively studied in their dust spectral A60, page 5 of 13 energy distribution (e.g., Beelen et al. 2006;Drouart et al. 2014;Leipski et al. 2014;Meyer et al. 2022). For the quasars J0338+0021, J2054-0005, and J1319+0950, we instead adopt the values published in Wang et al. (2013) and Leipski et al. (2014), based on the fits of their Spectral Energy Distributions. In all cases, the IR luminosities are computed by integrating the dust continuum between 8-1000 µm.

We convert line fluxes into luminosities, following, for instance, Carilli & Walter (2013):

where F line is the line integrated flux.

In Fig. 4, we compare the luminosity of the lines under examination with the infrared luminosity. Our sample spans a factor ∼10 in IR luminosity, but only a factor of 3 in the luminosity of any of the detected lines. Both the measured CO(6-5) and CO(7-6) luminosities are in line with empirical scaling relations based on local ULIRGs (see, e.g., Greve et al. 2014). This, together with the widespread detections of CO(7-6), and the similar luminosities found for CO(6-5) and CO(7-6), suggests that the CO excitation is relatively high up to J up ≈ 7. We measure [C i] luminosities of (3 -10) × 10 7 L , as expected based on the empirical [C i]-IR relation from Valentino et al. (2020). Finally, the [C ii]/IR luminosity ratios in our sample span a factor ∼6, from 1.4 × 10 -3 (J2219+0102) to 2.4 × 10 -4 (J2054-0005), with lower [C ii]/IR ratios typically observed for quasars with L IR > 3 × 10 12 L (see, e.g., Stacey et al. 2010;Gullberg et al. 2015). For comparison, the observed ratios in low and 5)-( 7) Fitted line integrated fluxes. Limits are at 3-σ significance. ( 8)-( 10) Inferred line luminosities. Capitalizing on the self-similarity of the line luminosities within our sample, we compute the average spectrum by shifting all of the 3 mm spectra to z = 6, and averaging the observed flux densities with weights set by the channel variance. The noise spectrum is derived as the sum of the inverse variance of the spectra of each source. We do not rescale the spectra based on the source luminosity. The stacked spectrum is then fitted in the similar fashion as individual sources. The stacked spectrum and the fitted values are reported in Fig. 2 and Table 3. In the following, when comparing the luminosity of the stacked lines with the ), derived as described in the text. (6-8) Observed line luminosity ratios. (a) For the line ratios of the stacked spectrum, we adopt the average value L [CII] = 2.5 × 10 9 L based on the median value of the literature measurements in Table 1.

one of [C ii], for the stack we adopt the median luminosity of [C ii] in our sample, L [CII] = 2.5 × 10 9 L .

Molecular gas masses

Here we infer molecular mass estimates via four recipes, based on the dust continuum emission and different molecular gas tracers. This allows us to test the consistency of the various mass estimators and to study the impact of working assumptions. First, we derive the dust mass from the dust continuum emission. As discussed in the previous section, τ ν 1 at 3 mm in our galaxy-scale observations. In this regime, the dust is optically thin and Eq. ( 1) yields a nearly linear dependence of the observed flux density on τ ν , and hence, on the dust mass. In the following, we adopt the template used for the estimate of the IR luminosiy (T d = 47 K, τ 1900 GHz = 0.2). Because the optical depth is fixed, Eq. ( 1) implies a constraint on the size of the emitting region, Ω, which thus sets the dust mass via Eq. ( 2). In this framework, an observed 3 mm flux density of 0.1 mJy for a source at z = 6 corresponds to a size of 0.097 arcsec 2 (∼3.2 kpc 2 ), and a dust mass of M d = 2.8 × 10 8 M . Our CO(7-6) observations have a median beam area of 16.8 arcsec 2 , yielding a typical filling factor f = 5.8 × 10 -foot_3 . If instead we leave τ ν free, the dependence on the angular size is lost in the optically thin regime, and the flux density is simply a function of temperature and dust mass. As discussed in Sect. 3.1, however, we cannot constrain the dust temperature in the majority of our sample. If we refer to the bestfit values from our continuum fits, T d = 30-110 K, the inferred dust masses differ from the fiducial ones by an average factor ∼1.6. We convert the dust mass into the associated molecular gas mass via a gas-to-dust ratio, δ g/d = 100 (see, e.g., Bolatto et al. 2013;Sandstrom et al. 2013;Genzel et al. 2015;Berta et al. 2016). For comparison, the analysis by Dunne et al. (2021) led to a mean value of δ g/d = 128 ± 16, in broad agreement with the round value assumed here. The inferred molecular gas masses are listed in Table 4.

Molecular gas mass estimates in high-z galaxies typically rely on CO observations (see, e.g., reviews in Carilli & Walter 2013;Combes 2018;Tacconi et al. 2020;Hodge & da Cunha 2020). While low-J transitions should be preferred, as they are less sensitive to uncertainties on the CO excitation, the molecular mass estimates in z > 6 quasars tend to rely on intermediate (J up = 5-7) transitions (e.g., Wang et al. 2010;Venemans et al. 2017a;Yang et al. 2019) which appear to be at the peak of the CO Spectral Line Energy Distribution for these quasars (e.g., Li et al. 2020). In this work, we pin our estimates to the CO(7-6) emission 3 . Molecular gas mass is derived from CO as:

where we adopt a CO-to-H 2 conversion factor α CO = 0.8 M (K km s -1 pc 2 ) -1 typical for ULIRG and quasars (see discussions in Bolatto et al. 2013;Carilli & Walter 2013), and r 71 is the CO(7-6)/CO(1-0) luminosity ratio, which we assume to be r 71 = 0.17 based on the high-z CO excitation template in Boogaard et al. (2020). For comparison, the wellstudied quasar J1148+5251 shows r 71 > 0.12, and r 73 = 0.58, suggesting very high CO excitation (Bertoldi et al. 2003;Riechers et al. 2009). Table 4 lists the inferred values.

A third path to molecular gas estimates capitalizes on the neutral carbon emission, which arises from the outer layers of molecular clouds as well as the diffuse molecular gas in the interstellar medium (see theoretical predictions in, e.g., Glover & Clark 2016;Clark et al. 2019;Bisbas et al. 2019;Galactic observations in, e.g., Frerking et al. 1989;Kramer et al. 2004;Cubick et al. 2008;Beuther et al. 2014; and searches at high redshift by, e.g., Walter et al. 2011;Alaghband-Zadeh et al. 2013;Bothwell et al. 2017;Popping et al. 2017;Papadopoulos et al. 2018;Valentino et al. 2018). In the optically thin assumption, following a three-level model of the carbon emission, the mass of neutral carbon can be inferred from the observed line luminosity using:

where Q ex = 1 + 3 exp(-23.6/T ex ) + 5 exp(-62.5/T ex ) is the partition function and T ex is the excitation temperature in Kelvin (see Weiß et al. 2003Weiß et al. , 2005)). We assume that the carbon emission is in thermal equilibrium with the dust, and thus adopt T ex = 47 K as measured on average in high-z quasars (Beelen et al. 2006;Leipski et al. 2014). The mass estimate in Eq. ( 6) is not very sensitive to T ex : any temperature in the range T ex = 30-100 K yields a mass estimate within ∼0.15 dex from our fiducial value.

The resulting masses can be converted into H 2 masses through the neutral carbon to H 2 abundance, [C]/[H 2 ], and accounting for the mass ratio of carbon atoms and H 2 molecules. Following Boogaard et al. (2020), here we adopt a carbon abundance value of (1.9 ± 0.4) × 10 -5 (under the assumption that CO is in thermal equilibrium; see Bothwell et al. 2017;Boogaard et al. 2020). This value is consistent with the

where mates show the largest scatter in our work, as they spread over ∼1 dex; on the contrary, the CO-based estimates display a much smaller range (<0.5 dex). The plot also highlights the impact of our operative assumptions: Our M H2 estimates depend linearly on α CO , δ g/d , and α [CII] , and antilinearly on r 71 , [C]/[H 2 ], and (to first order) T d . Different assumptions for these parameters could mitigate or increase discrepancies among these mass estimates. For instance, by adopting a higher CO-to-H 2 value, α CO = 3.6 M (K km s -1 pc 2 ) -1 (see, e.g., Daddi et al. 2010), the CO-based estimates of M H2 would increase by a factor 4.5, thus narrowing the gap with respect to the [C ii]-based ones, but at the same time separating them from the [C i]-and dust-based measurements. On the other hand, if we adopt a higher COexcitation, for example, r 71 = 0.38 from the quasar template in Carilli & Walter (2013) instead of the fiducial value r 71 = 0.17 from Boogaard et al. (2020), we would reduce the CO-based dust mass estimates, thus widening the gap with respect to the [C ii]-based estimates. Similarly, assuming a higher carbon abundance, [C]/[H 2 ] = 5×10 -5 as estimated by Weiß et al. (2005) in a sample of z ∼ 2.5 submillimeter galaxies, would yield a decrease in the [C i]-based mass of molecular gas. Along the same lines, one could adopt the method proposed by Sommovigo might be suppressed by thermalization, changes in the photoelectric efficiency of dust grains, optical depth, and other processes; see e.g. Sutter et al. 2021). While a deeper understanding of the physical properties of the line and dust emitting ISM is required in order to directly constrain these unknown quantities, our analysis demonstrates that the fiducial values lead to consistent molecular gas mass estimates (within a factor 3).

Physical properties of the cold ISM

In order to further explore the interplay among CO, [C i], [C ii], and dust in high-z quasars, we assume that radiative processes (photodissociation, photoelectric heating, etc.) dominate over other nonradiative mechanisms altering the energy budget in the molecular gas (shocks, turbulence, cosmic rays, etc). This assumption allows us to directly compare our observations with the line and continuum luminosity predictions of Photo-Dissociation Regions (PDRs) and X-ray Dominated Regions Fig. 6. Luminosity line ratios predicted in the context of PDR and XDR models, as a function of column density, gas density, and strength of the incident radiation field (see text for details). The contours show the ratios measured in each source in our sample, color coded as shown in the legend. XDR models can reproduce the observed high [C ii]/[C i] ratios only at low column densities, N H 10 22 cm -2 , whereas PDR models are compatible with the observed constraints over a wider range of parameters. The relatively low [C ii]/CO(7-6) ratio observed in the targeted quasars is suggestive of typically large column densities, N H 10 23 cm -2 for both PDR and XDR environments. Finally, the low observed [C i]/CO(7-6) ratios point to a high density n H 10 3 (10 5 ) cm -3 in PDRs (XDRs). Overall, a dense (n H > 10 3 cm -3 ) gas cloud with column densities of N H ∼ 10 23 cm -2 impinged by UV radiation from newborn stars appears to best describe the observed line ratios. However, we stress that no single model accurately reproduces simultaneously all of the observed line ratios, thus suggesting that a single cloud is too simplistic a description of the observed phenomenology.

(XDRs). The models are based on the analysis presented in Pensabene et al. (2021) and briefly summarized here. We use the CLOUDY radiative-transfer code (version c.17.01;Ferland et al. 2017) to predict the line emission of a single plane-parallel semi-infinite cloud impinged by a radiation field in both the PDR and XDR regimes. We created grids of 270 PDR and XDR models with total hydrogen density in the range log(n H [cm -3 ]) = [2, 6] (15 values in steps of ∼0.29 dex) and strength of the incident radiation field (18 models) in the range log G 0 [Habing units] = [1, 6] (for PDRs), and log F X [erg s -1 cm -2 ] = [-2.0, 2.0] (for XDRs). We consider three cases of total hydrogen column density log N H [cm -2 ]={22, 23, 24} (see Fig. 6). We assume that all the line emission arises from the neutral or molecular gas phase, in particular the contribution of Hii regions to the [C ii] emission is neglected (see, e.g., Díaz-Santos et al. 2017;Pensabene et al. 2021).

The relative strength of all of the targeted lines is sensitive to the gas density, as well as to the intensity and hardness of the incident radiation field. For instance, the [C i]/CO(7-6) ratio is sensitive to the gas density in the regime n H < 10 4 cm -3 , and only marginally affected by other parameters for N H 10 23 cm -2 . The [C ii]/[C i] typically shows values >10 in PDRs, while it is always <10 in XDRs, unless N H 10 22 cm -2 . This is because X-rays penetrate deeper into the clouds than UV photons, thus heating the cloud cores and enhancing [C i]. For the same reason, CO molecules receive additional energy and thus the CO spectral line energy distribution shows stronger emission in the high-J transitions. This translates into a suppressed [C ii]/CO(7-6) luminosity ratio in XDRs compared with predictions for PDRs, at least for high-density clouds (n H > 10 4 cm -3 ). At low column densities, N H 10 22 cm -2 , predictions for this set of observables based on PDR and XDR models tend to converge. The expected [C ii]/[C i] ratios decrease in both models. Most notably, CO(7-6) emission is largely suppressed, thus leading to high [C ii]/CO(7-6) and [C i]/CO(7-6) luminosity ratios.

As discussed in the previous sections, the line luminosities in our sample span a small range (<0.5 dex). This implies that also the probed line ratios show similarly homogeneous values within the sample. We argue that the selection of IR-luminous (L IR > 10 12 L ) sources might bias our sample toward a specific subclass of host galaxies hosting a compact starburst. In Fig. 6, we compare the predicted line ratios with the values observed for the quasars in our sample. We find that PDR models successfully predict the observed [C ii]/[C i] ratios over a wide range of column densities. That is not the case for XDR models, which require low column densities in order to account for the high [C ii]/[C i] luminosity ratios in our data. All models predict a wide range of [C ii]/CO(7-6) luminosity ratios, although at low column densities it becomes increasingly difficult to explain the low values observed in our study, in particular if compared to predictions from PDR models. Finally, the also displays a range of CO(7-6)/IR ratios comparable with other high-z sources, but a factor ∼ 2-3× higher than the bulk of the local galaxies.

Right: predictions of the [C ii]/[C i] ratio as a function of the CO(7-6)/IR luminosity ratio for various gas densities log n H [cm -3 ] = (3, 4, 5) for a PDR with incident radiation field log G 0 = (2, 3, 4) and an XDR with log F X = (-1, 0, +1). Our sources (shown in filled circled, with the same color scheme as in the left-hand panel) are best described by a PDR-like environment with N H ∼ 10 23 cm -2 , n H 10 4 cm -3 , and G 0 ∼ 10 3 . observed [C i]/CO(7-6) luminosity ratios point to intermediate to high densities (n H 10 3 cm -3 in PDRs, n H 10 5 cm -3 in XDRs).

In Fig. 7, we compare the observed [C ii]/[C i] and CO(7-6)/IR luminosity ratios in individual sources from our sample with other high-redshift quasars and SMGs (from Weiß et al. 2003Weiß et al. , 2013;;Riechers et al. 2009;Bradford et al. 2009;Ivison et al. 2010;Cox et al. 2011;De Breuck et al. 2011;Danielson et al. 2011;Wagg et al. 2012;Salomé et al. 2012;Walter et al. 2012;Riechers et al. 2013;Rawle et al. 2014;Leipski et al. 2014;Decarli et al. 2014;Venemans et al. 2017a;Yang et al. 2019;Zhao et al. 2020;Pensabene et al. 2021) as well local galaxies (based on the compilation in Rosenberg et al. 2015) and with the predictions from the models presented in Fig. 6. The quasars in our study appear to populate a rather narrow parameter space in this diagram: They all have [C ii]/[C i] 20, consistent with most of the high-z sources in our comparison, as well as the bulk of the local sample from Rosenberg et al. (2015). Notably, only the local Seyfert 2 galaxy NGC 1068 falls (although, only marginally) in the [C ii]/[C i] < 10 regime that is often adopted as a bona fide XDR domain. In terms of CO(7-6)/IR, we find values of 10 -5 -10 -4 , in line with most of the high-z comparison sample, and (on average) a factor 2-3 higher than most of the local galaxies in Rosenberg et al. (2015).

In Fig. 7 we also show how the individual sources presented in our study compare with respect to the predicted [C ii]/[C i]

and CO(7-6)/IR luminosity ratios from our PDR and XDR models. As already mentioned, the high values observed for the [C ii]/[C i] ratio point to a relatively low gas column density, N H 10 23 cm -2 for PDRs, 10 22 cm -2 for XDRs. However, at very low column densities, N H ∼ 10 22 cm -2 , the CO(7-6) luminosity is suppressed, thus models fail to reproduce the relatively high CO(7-6)/IR ratios observed in our quasars. The quasar population in our study appears to be best described by PDR models with column densities of N H ∼ 10 23 cm -2 , densities of n H 10 4 cm -3 , impinged with a radiation field with G 0 ∼ 10 3 . These modest column density estimates are in agreement with the estimates based on the dust continuum slope discussed in Sect. 3.1. At face value, these constraints imply a lineof-sight thickness of the medium of only a few parsecs. However, we stress that the simplified models adopted here (in particular, because of the single-cloud treatment at fixed gas density, and because we neglect nonradiative ingredients contributing to the gas excitation, such as turbulence and cosmic rays) likely fail to capture the complex interplay of different phases in the ISM, thus these findings need to be taken with caution.

Conclusions

We present a survey of [C i], CO(7-6), CO(6-5) and their underlying dust continuum (observed at 3 mm) in a sample of ten quasar host galaxies at z ∼ 6. All of the quasars in our sample had previously been detected in their [C ii] and dust continuum at ∼1 mm (in the observer's frame). Our main findings are:

(i) We detect CO(7-6) in all of the targeted quasars, and [C i] in 6 of them. Out of the 8 quasar host galaxies for which we have CO(6-5) coverage, we also report six detections in this emission line. The underlying 3 mm dust continuum is detected in 4/10 of the sample. This work doubles the number of z ∼ 6 quasars with CO(7-6) measurements in the literature. (ii) The observed flux density of the dust continuum emission at 3 mm and 1 mm is consistent with a modified black body emission with temperature T d = 47 K and an optical depth at the frequency of [C ii] of τ 1900 GHz = 0.2.

(iii) Our targets span a factor ∼10 in IR luminosity, but only a factor 3 in the luminosity of any line detected in our study. (iv) We derive molecular gas masses using four independent methods, based on the dust continuum and the CO, [C i], and [C ii] line emission. To first order, all the methods point to similar values of M H2 = 10 10-11 M . This confirms earlier results suggesting that immense gaseous reservoirs power the intense star formation and nuclear activity of the first quasar host galaxies. By combining the four independent measurements, we infer relative constraints on the unknown scaling factors (gas-to-dust, α CO , carbon abundance, α [CII] , etc.). Our estimates for these conversion factors agree with the values inferred in a similar study by Dunne et al. (2021) the expectations from radiative transfer models, we find that the best (albeit imperfect) description of the observations requires a PDR-like environment with dense (n H > 10 4 cm -3 ) clouds impinged by a radiation field of G 0 ∼ 10 3 (in Habing units). Column densities of N H ∼ 10 23 cm -2 appear favored by our modeling, thus yielding a typical size of the clouds of only a few parsec. This study confirms and further constrains the presence of massive reservoirs of molecular gas in the host galaxies of z ∼ 6 quasars, and suggests that star formation is the main mechanism responsible for the gas excitation at galactic scales. This work demonstrates the diagnostic power of combined 3 mm and 1 mm observations of very high-redshift sources, which can reveal a wide range of properties of the cold ISM in massive galaxies at cosmic dawn with only two frequency setups. Further NOEMA and ALMA observations will allow us to expand this analysis to a broader range of IR luminosities, and to extend the investigation to other sets of diagnostics (e.g., dense gas tracers, metallicity tracers, etc.).