The ALMA Reionization Era Bright Emission Line Survey The molecular gas content of galaxies at z ∼ 7

A key to understanding the formation of the first galaxies is to quantify the content of the molecular gas as the fuel for star formation activity through the epoch of reionization. In this paper, we use the 158 µ m [C ii ] fine-structure emission line as a tracer of the molecular gas in the interstellar medium (ISM) in a sample of z = 6 . 5 − 7 . 5 galaxies recently unveiled by the Reionization Era Bright Line Emission Survey, REBELS, with the Atacama Large Millimeter / submillimeter Array. We find substantial amounts of molecular gas ( ∼ 10 10 . 5 M ⊙ ) comparable to those found in lower redshift galaxies for similar stellar masses ( ∼ 10 10 M ⊙ ). The REBELS galaxies appear to follow the standard scaling relations of molecular gas to stellar mass ratio ( µ mol ) and gas depletion timescale ( t dep ) with distance to the star-forming main-sequence expected from extrapolations of z ∼ 1 − 4 observations. We find median values at z ∼ 7 of µ mol = 2 . 6 4 . 1 − 1 . 4 and t dep = 0 . 5 + 0 . 26 − 0 . 14 Gyr, indicating that the baryonic content of these galaxies is gas-phase dominated and little evolution from z ∼ 7 to 4. Our measurements of the cosmic density of molecular gas, log( ρ mol / ( M ⊙ Mpc − 3 )) = 6 . 34 + 0 . 34 − 0 . 31 , indicate a steady increase by an order of magnitude from z ∼ 7 to 4.


Introduction
One of the most important advances in astrophysics in the past two decades has been the determination of the evolution of the cosmic star formation rate (SFR) density from cosmic dawn to the present time.This quantity was found to rapidly rise from z = 10 to z = 3 (a span of 3 Gyrs of cosmic time), reaching a plateau at z = 1 − 3, and then steadily declining at z < 1 (e.g., Madau & Dickinson 2014;Péroux & Howk 2020).The cosmic SFR activity appears dominated by star-forming galax-ies that grow secularly, forming a tight correlation in the SFR versus stellar mass plane, usually termed as the star-forming "main-sequence" (MS; Brinchmann et al. 2008;Daddi et al. 2007;Elbaz et al. 2007Elbaz et al. , 2011;;Noeske et al. 2007;Peng et al. 2010;Rodighiero et al. 2010;Whitaker et al. 2010Whitaker et al. , 2014;;Speagle et al. 2014;Schreiber et al. 2015;Iyer et al. 2018;Di Cesare et al. 2023).These galaxies sporadically will be involved in major galaxy interactions or mergers, leading to increased SFRs (e.g.; Kartaltepe et al. 2012), and eventually will halt star for-mation, possibly through exhaustion of their cold gas reservoirs that sustain star formation (Peng et al. 2010;Sargent et al. 2014;Spilker et al. 2018;Belli et al. 2021;Williams et al. 2021;Bezanson et al. 2019Bezanson et al. , 2022)).
This framework for galaxy growth has been built through intensive multiwavelength observational campaigns and simulations.Determination of the molecular gas content has been paramount, as the cold gas represents the fuel for star-forming activity.Since directly observing H 2 molecular gas is difficult, measurements of the molecular gas mass in distant galaxies have focused on the dust continuum and the CO and [C i] line emission (e.g., Heintz & Watson 2020).Dedicated targeted surveys and archival searches for such dust and CO/[C i] line observations have led to thousands of galaxies with measured molecular gas masses (e.g.; Daddi et al. 2008Daddi et al. , 2010;;Tacconi et al. 2013Tacconi et al. , 2018;;Dessauges-Zavadsky et al. 2015, 2017;Freundlich et al. 2019;Scoville et al. 2014Scoville et al. , 2017;;Fudamoto et al. 2017;Liu et al. 2019;Valentino et al. 2018Valentino et al. , 2020)).A complementary approach has been performing deep systematic observations of dust continuum and CO line emission in contiguous patches of the sky centered on cosmological deep fields (Carilli & Blain 2002;Decarli et al. 2014;Walter et al. 2014Walter et al. , 2016;;Aravena et al. 2016b;Dunlop et al. 2017;Pavesi et al. 2018;Riechers et al. 2019;Franco et al. 2020;Hatsukade et al. 2018;González-López et al. 2019, 2020;Decarli et al. 2020).Such deep field surveys have the advantage of avoiding pre-selection biases compared to targeted surveys, although by design they are restricted to fainter objects.Overall, these programs have enabled establishing a set of scaling relations that allow us to describe the galaxy growth in terms of the galaxies' stellar mass, specific SFRs, and redshift, out to z ∼ 4. In particular, current studies indicate that from z = 3 to 1, there is an increase in the molecular gas depletion timescales, a decrease in the molecular gas fractions, and a decrease of the specific SFR.Also, the molecular gas depletion timescales decrease with increasing specific SFR (e.g.; Tacconi et al. 2020).Of particular interest has been the determination of the cosmic density of molecular gas (ρ H2 ) as a function of redshift from various dust continuum and CO line surveys out to z ∼ 5 (Walter et al. 2014;Decarli et al. 2016Decarli et al. , 2019Decarli et al. , 2020;;Riechers et al. 2019;Magnelli et al. 2020).This quantity has been measured to closely follow the evolution of the cosmic density of SFR, supporting a relatively unchanged star formation efficiency on cosmic scales.Based on these results, it has also become evident that gas accretion from the intergalactic medium (IGM) is necessary to sustain the build-up of stellar mass in galaxies, at least since z = 1.5 (Walter et al. 2020).
Due to cosmic dimming, the decreasing metallicities at higher redshifts, which make the dust continuum and CO/[C i]line emission fainter, and the increase of the cosmic microwave background with redshift, it becomes increasingly difficult to measure the molecular gas reservoirs of galaxies at z > 4 (e.g.;Carilli & Walter 2013).Since the [C ii] line is among the brightest coolant lines in galaxies, and because it is redshifted into the millimeter atmospheric windows observable from the ground at z > 4, it has thus become a powerful tool to measure the interstellar medium (ISM) properties of galaxies.Even though its emission arises from various environments, including the neutral, ionized, and molecular gas phases (e.g.; Vallini et al. 2015;Pallottini et al. 2017;Lagache et al. 2018;Olsen et al. 2017), it has been reliably calibrated and used as a molecular gas tracer for massive galaxies at z < 0.2 (Hughes et al. 2017;Madden et al. 2020) and at z > 2 (Zanella et al. 2018;Vizgan et al. 2022).While [C ii] has been traditionally used as a tracer of star formation in galaxies, the close link between star formation and molecular gas content (i.e. the Schmidt-Kennicutt relation) yields naturally a relation between [C ii] and molecular gas.Current studies show that [C ii] would even trace the molecular gas better than CO for low-metallicity environments where the gas is CO-dark, which is the case for high-redshift galaxies (Madden et al. 2020;Vizgan et al. 2022).
The brightness of the [C ii] line has allowed for the molecular gas measurements of galaxies at z∼ 4 − 6 from the Atacama Large Millimeter/submillimeter Array (ALMA) Large Program to INvestigate C+ at Early Times (ALPINE) survey (Dessauges-Zavadsky et al. 2020).
Here, we aim to extend such exploration to even higher redshifts, by using data from the Reionization Era Bright Emission Line Survey (REBELS; Bouwens et al. 2022).REBELS is an ALMA large program aimed at confirming the redshifts of a sample of 40 UV-bright star-forming galaxy candidates at z ∼ 6 − 9 using [C ii] 158 µm and [O iii] 88µm emission line scans of the high-probability redshift range obtained from photometric redshift measurements.REBELS has enabled studies of the properties of these early galaxies through simultaneous studies of the dust continuum and [C ii] emission line, including: determination of the dust properties (Inami et al. 2022;Sommovigo et al. 2022;Ferrara et al. 2022;Dayal et al. 2022), the IR luminosity function (Barrufet et al. 2023) and obscured SFR density and fractions at z ∼ 7 (Algera et al. 2023b), measurements of the dust temperature and [OIII]]/[CII] line ratios in a few systems (Algera et al. 2023a), specific SFRs (Topping et al. 2022), the discovery of inconspicuous overdensities of dusty star-forming galaxies (Fudamoto et al. 2021), measurements of cosmic HI gas mass densities (Heintz et al. 2022), [CII] sizes (Fudamoto et al. 2022) identification of Lyman-α emission and velocity offsets (Endsley et al. 2022), and a resolved study of a massive starforming system at z = 7.3 (Hygate et al. 2023).
In this paper, we present molecular gas mass estimates for a sample of 28 [C ii]-detected star-forming galaxies at z ∼ 6 − 8 drawn from REBELS (Bouwens et al. 2022) and its pilot surveys (Smit et al. 2018;Schouws et al. 2022b,a).These measurements are thereby used to provide estimates of the evolution of the molecular gas depletion timescales, molecular gas ratios, and cosmic molecular gas density at z ∼ 7.In Section 2, we present details about the REBELS survey, the available multi-wavelength data, and comparison samples used throughout this work.In Section 3, we discuss the use of the [C ii]-based molecular gas calibration in the REBELS sample, comparing it to alternative estimates of the molecular gas content available.In Section 4, we present measurements of the molecular gas depletion timescales, gas fraction, and cosmic molecular gas density, and discuss them in the context of the established scaling relations and predicted evolution of these quantities out to z ∼ 7.In Section 5 we provide a summary of the main results of this work.Throughout the paper, we adopt a standard ΛCDM cosmology1 , and a Chabrier (2003) initial mass function (IMF), with stellar masses ranging from 0.1 − 300 M ⊙ .

Observations and sample
Our study takes advantage of the [C ii] line emission observations in a sample of z ∼ 6 − 8 galaxies obtained by the ALMA REBELS large program (PID: 2019.1.01634.L).The REBELS survey details are provided in Bouwens et al. (2022) and Schouws et al. (in preparation).Here, we provide a brief description of the observations.
The REBELS program used line scans to search for redshifted [C ii] 158µm and [O iii]88µm lines from sources that had previously been selected as candidate z > 6 galaxies based on their UV luminosity (M UV < −21.5 mag) and photometric redshift probability distributions (6.5 < z < 9.5), following an approach similar to the pilot programs (Smit et al. 2018;Schouws et al. 2022b).
The REBELS survey targeted 40 galaxies selected from a total area of 7 deg 2 encompassing various deep fields (Sec.2.2).Line scans searching for the [C ii] and [O iii] lines were conducted for 36 galaxies at z ∼ 6.5 − 7.5 and four galaxies at z ∼ 8.0 − 9.5 using ALMA bands 6 and 8, respectively.For each target, the scanning range was chosen to cover 90% of the relevant frequency range following the photometric redshift probability distribution.
The ALMA observations used in this study were performed in the most compact configurations, C43-1 and C43-2, yielding synthesized beam sizes of 1.2 − 1.6 ′′ in band 6.The [CII] line scans were performed to a sensitivity of 2 × 10 8 L ⊙ at z = 7 (5σ; Bouwens et al. 2022).
Data calibration was done using the Common Astronomy Software Application (CASA) software following the observatory-provided pipeline.Datacubes were thus created, inverting the calibrated visibilities with the tclean task and using a natural weighting scheme to maximize the sensitivity.Continuum images were obtained by collapsing the datacubes and excluding the frequency ranges that contained the [C ii] lines (when detected) within 2×FWHM [CII] .The [C ii] line datacube and the continuum image were searched for emission at the location of the rest-frame UV source, following the procedures detailed in Schouws et al. (in preparation) and Inami et al. (2022).
We obtained 23 [C ii] line and 16 dust continuum detections, respectively, from the first year data.In this study, we thus focus on 23 galaxies with observations of [C ii] line and dust continuum emission at z ∼ 7, plus additional five galaxies with [C ii] detections (two of them with dust detections) from the pilot studies (Smit et al. 2018;Schouws et al. 2022a).We deliberately exclude [CII] non-detections from the analysis, as it is unclear whether they are [CII] faint galaxies or they are sources with redshifts falling outside the line scan covered ranges.The redshift distribution of the 28 confirmed [C ii] line emitters considered in this work is shown in Figure 1 (left).The mean redshift of these sources is z = 6.95, spanning a range of z = 6.50 − 7.68.
The available photometry and the [C ii]-derived spectroscopic redshifts were thus used to infer the galaxies' physical properties, including stellar masses as described by Topping et al. (2022).Estimates of total SFRs was done through restframe UV observations from HST (Bouwens et al. 2022) and estimates of the IR luminosities from the ALMA rest-frame 158µm continuum measurements (for details see: Inami et al. 2022).
Spectral Energy Distribution (SED) fitting for the REBELS sample was conducted using the Prospector code (Leja et al. 2019).Two different approaches were employed: a nonparametric prescription for the galaxies' star formation histories (SFHs) (SFHs; Johnson et al. 2021) and a constant SFH.In both methods, consistent physical assumptions were made, including sub-solar metallicities (0.2 Z ⊙ ), a Chabrier (2003) initial mass function (IMF) for stellar masses below 300 M ⊙ , and a Small Magellanic Cloud (SMC) dust extinction curve (see: Topping et al. 2022).
Additionally, SED fitting for the REBELS sample was performed using the Beagle code (Chevallard & Charlot 2016) with a constant SFH (Stefanon et al., in preparation) and identical physical assumptions, except for the fitting code and stellar population templates.Some of the derived properties obtained through this procedure have been reported in Bouwens et al. (2022).
As elaborated in detail by Topping et al. (2022), the use of Prospector with a constant SFH yields results similar to those obtained by Beagle.However, it is known that adopting constant SFHs can result in younger age estimates for starbursting galaxies, as the starburst can dominate the signal over older stellar populations (i.e. the 'outshining' problem: young stellar populations outshine the old ones).In such cases, non-parametric SFHs provide greater flexibility and can help mitigate this issue (Algera et al. 2023b).Nevertheless, these models tend to yield larger estimates of stellar masses, with an average increase of 0.43 dex compared to constant SFH models (Topping et al. 2022).For consistency with a related parallel study of the HI mass content of REBELS galaxies (Heintz et al. 2022), we adopt the results from Prospector with non-parametric SFHs for this work.Using this approach, the majority of the REBELS galaxies are consistent within the uncertainties with the standard star formation main-sequence extrapolated to z = 7 (Algera et al. 2023b).In Section 4.2.1, we discuss potential biases and the effect on our results introduced by the use of these non-parametric models.

Comparison samples
Our prime comparison sample is based on the ALPINE survey (Le Fèvre et al. 2020).This survey obtained [C ii] emission line and dust continuum observations for a sample of 118 star-forming galaxies in the redshift range z = 4 − 6 (Béthermin et al. 2020).These galaxies were selected based on their rest-frame UV properties and the availability of rest-UV spectroscopic redshift measurements.The derived stellar masses, in the range 10 8.4 −10 11 M ⊙ , and SFRs, in the range 3−630 M ⊙ yr −1 , follow the expected trend of the star-forming main sequence at z = 4 − 6.The ALPINE observations yielded [CII] line and dust continuum detections for 75 and 23 sources, respectively, at a signal-to-noise level of 3.5.While the range of stellar masses and SFRs covered by ALPINE is similar to that of REBELS (see Fig. 1), by design, the latter did not require the existence of previous rest-UV (e.g., Ly-α) redshift determinations.
Figure 1 (middle panels) shows the distribution in stellar mass and SFR for the REBELS galaxies, compared with that of the ALPINE sample.Our current REBELS sample covers wide ranges in both parameters, with M stars ∼ 10 8.5 − 10 10.5 M ⊙ and SFR∼ 30 − 300 M ⊙ yr −1 , matching the range covered by ALPINE.Figure 2 shows the location of the REBELS sources in the SFR(UV+IR)-stellar mass plane, compared to the expected location of the star-forming main sequence (MS) at z ∼ 5 − 7, based on various prescriptions (Schreiber et al. 2015;Iyer et al. 2018;Khusanova et al. 2021).For clarity, we split our sample into sources detected and not in 158µm dust continuum emission.Most galaxies in our sample of [CII] detected sources fall within the boundaries of the MS at these redshifts (±0.3 dex) compared to the widely used Schreiber et al. (2015) prescription for galaxies at z = 7.However, they appear to be above the MS compared with the Iyer et al. (2018) prescription at z = 7, although still within ±0.5 dex.
Based on this comparison, we find that REBELS galaxies have parameter distributions similar to the ones from the ALPINE sample.We note that the REBELS galaxies were not confirmed spectroscopically through rest-UV lines (Lyman-α) as was the ALPINE sample, although both samples were selected based on their rest-UV brightness.Furthermore, we find that most REBELS galaxies are consistent with being in the MS and thus represent examples of massive yet "typical" star-forming galaxies at high redshift.A similar conclusion is found by a recent study by Di Cesare et al. (2023).
At lower redshifts (z < 4), we adopt the compilation of CO line and dust continuum measurements obtained by the IRAM Plateau de Bureau HIgh-z Blue Sequence Survey (PHIBSS; Tacconi et al. 2013Tacconi et al. , 2018;;Freundlich et al. 2019).
The PHIBSS is a targeted survey that provides stellar masses, SFRs, and molecular gas masses (from dust and CO observations) for large samples of 1444 massive star-forming galaxies at z = 0 − 4.
We complement this with similar data from the ALMA Spectroscopic Survey in the Hubble UDF (ASPECS; Walter et al. 2016;Aravena et al. 2019Aravena et al. , 2020;;Boogaard et al. 2020) and other surveys (Papovich et al. 2016;Magdis et al. 2017;Gowardhan et al. 2019;Kaasinen et al. 2019;Molina et al. 2019;Bourne et al. 2019;Pavesi et al. 2018;Cassata et al. 2020).We remark that since the [CII] emission line is difficult to access at these lower redshifts, the molecular gas mass estimates in the comparison samples are based on observations of CO and dust continuum emission.

Measurements of molecular gas
Several far-infrared to radio emission line and dust continuum measurements have been used as standard tracers of the molecular gas content in galaxies (e.g., Carilli & Walter 2013;Hodge & da Cunha 2020).These include the CO line, particularly in the lowest energy transitions, the [C i] line, and dust continuum emission in the Rayleigh-Jeans regime.These tracers have provided molecular gas mass measurements in galaxies out to z ∼ 3 − 4.However, due to their faintness, these tracers become harder to detect at high redshifts with current facilities.Furthermore, at z > 6, the Cosmic Microwave Background (CMB) temperature becomes comparable to the excitation temperature of the most standard tracers like CO J < 3 and dust emission in the Rayleigh-Jeans tail (da Cunha et al. 2013).The radiation from these tracers thus becomes almost impossible to detect against the CMB at these redshifts.
The [C ii] 158µm line, due to its brightness, high excitation temperature, and accessibility from the ground at the observed redshifted frequencies (due to the atmosphere transparency), has appeared as a top ISM tracer in galaxies at z > 4. Zanella et al. (2018) found a correlation between the [C ii] luminosity and the molecular gas mass with a mean absolute deviation of 0.3 dex, and without evident systematic variations.They find a [CII] luminosity to molecular gas mass conversion factor α [CII] ≃ 30 (M ⊙ /L ⊙ ).As indicated by Zanella et al. (2018), the [C ii]-to-H2 conversion factor, α [CII] , appears to be largely independent of the molecular gas depletion time, metallicity, and redshift for galaxies in the star-formation main sequence.Such correlation appears to hold for galaxies on and above the MS, in the redshift range 0 < z < 6, and with metallicities 12+log(O/H)=7.9 to 8.8.A larger scatter is observed for galaxies above the MS at lower metallicities.
There is vast literature indicating that the [C ii] emission in galaxies arises from various ISM phases, including ionized, neutral, and molecular gas (for a detailed discussion see: Lagache et al. 2018;Ferrara et al. 2019;Pallottini et al. 2019;Dessauges-Zavadsky et al. 2020).Numerical simulations and observations have shown evidence that > 60% of the [C ii] emission is dominated by molecular clouds and photon-dominated regions (PDRs), with a minor contribution of 20 − 25% each from diffuse ionized and neutral gas (Vallini et al. 2015;Olsen et al. 2017;Katz et al. 2017Katz et al. , 2019;;Pallottini et al. 2017Pallottini et al. , 2022;;Vizgan et al. 2022).As such, [C ii] can be readily calibrated as a tracer of these various phases, including star formation (e.g., De Looze et al. 2014), neutral gas (Heintz et al. 2021) and molecular gas (Zanella et al. 2018;Madden et al. 2020).These calibrations are a likely function of stellar mass (and/or metallicities), evolutionary stage, and cosmic time.The main caveat to using [C ii] as a molecular gas tracer is the unknown fraction of [C ii] emission that arise from PDRs.However, the current results indicate that [C ii] is a tracer as good as others, such as CO and dust for galaxies with stellar masses ∼ 10 10 M ⊙ , particularly considering the need to assume an excitation ladder and a metallicity-dependent α CO conversion factor for CO measurements, and given the necessary assumption of a gas-to-dust conversion (or similar calibrations) for dust continuum based measurements.
Based on this, we measured molecular gas masses for the REBELS galaxies in our sample using Equation (2) from Zanella et al. (2018), and following the approach adopted for the ALPINE survey by Dessauges-Zavadsky et al. (2020), We thus adopt α . Figure 1 (right) shows the distribution of molecular gas masses obtained in this way.The masses range between 10 9 − 10 11 M ⊙ , comparable to the range obtained for the ALPINE sample.Evidently, this is a consequence of the constant conversion factor and the same stellar mass and SFR ranges.
Before analyzing further the implications of these molecular gas measurements, in the rest of this section, we focus on checking the [C ii]-based molecular gas mass estimates against other derivations provided by the available data.

IR luminosity as a proxy of molecular gas mass
One of the best-established relations in extragalactic astronomy is that between the IR luminosity (L IR ) and the CO(1-0) luminosity (L ′ CO ).This has long been known to hold for a variety of galaxies and environments, from local spirals (Leroy et al. 2008) to the most extreme starburst galaxies known through a wide range of redshifts, out to z ∼ 6 (e.g., Aravena et al. 2016c).Following Dessauges-Zavadsky et al. ( 2020), we convert the [C ii]-derived molecular gas masses (M mol, [CII] ) into the expected CO(1-0) luminosities under the assumption of typical α CO values, with L ′ CO = M mol,[CII] /α CO .The α CO is known to depend on the metallicity, although with significant scatter (e.g., Bolatto et al. 2013), and thus a range of typical values is considered between 1.0 to 4.36 M ⊙ (K km s −1 ) −1 .
Figure 3 (left) shows the L IR vs. L ′ CO diagram for the REBELS targets compared to standard curves representing the ranges occupied by star-forming galaxies, for reference.This includes the universal fit to the compilation of CO observations for star-forming galaxies by Dessauges-Zavadsky et al. (2015), along with the same line with a factor ±0.3 dex added, representing the typical scatter around the fit.
We find that assuming a standard value of α CO = 4.36 M ⊙ (K km s −1 ) −1 , similar to what is found in the Milky Way and mainsequence galaxies at z ∼ 2, the REBELS galaxies fall right on top of the Dessauges-Zavadsky et al. ( 2015) curve.Variations of α CO between 2.4−6.4M ⊙ (K km s −1 ) −1 would still yield a correlation within the 0.3 dex of the L IR vs. L ′ CO relation.Conversely, using a lower α CO of 1.0 (same units), comparable to what is found in local starbursts (Downes & Solomon 1998), would lead to significantly larger CO luminosities being inconsistent with the L IR vs. L ′ CO relations.Using values of α CO > 6.4 M ⊙ (K km s −1 ) −1 , as expected for sub-solar metallicities for high redshift galaxies, would also lead to L ′ CO values incompatible with the L IR vs. L ′ CO relation.We note that the values of L IR appear to have a threshold at ∼ 10 11.5 L ⊙ , which is due in part to the depth of the REBELS program.However, we caution that the REBELS L IR values were computed using a single dust continuum measurements and realistic assumptions of the dust temperature (46 K).Warmer (or colder) dust temperatures will yield larger or lower L IR values.
While it is not surprising to find a linear scaling between both quantities, as the [C ii] luminosities are known to scale with SFR, the fact that the [C ii]-based CO values are consistent with the known L IR -L ′ CO relation supports the idea of using [C ii] as a molecular gas tracer in the REBELS sample.

Dust masses
The dust content in galaxies is related to the molecular gas mass through the gas-to-dust ratio (δ GDR ).This parameter depends strongly on the galaxies' metallicities and thus on the stellar mass through the mass-to-metallicity (MZ) relation.
Dust mass (M d ) estimates for the REBELS galaxies have recently been obtained by two independent studies by Ferrara et al. (2022) and Sommovigo et al. (2022).It is thus interesting to check how well the dust masses derived by these studies relate to the [C ii]-based molecular gas measurements in the REBELS sample.The former constructed a semi-empirical model for dust reprocessing that uses the UV spectral slope, the observed UV continuum flux at 1500Å and the observed FIR continuum at 158µm as input parameters, and considers various extinction curves and dust geometries.The latter adopted the Sommovigo et al. (2021) method to derive the dust properties of galaxies based on a combination of the [C ii] line emission and the underlying 158µm continuum, using standard functional forms with varying normalization for the star formation law ("Schmidt-Kennicutt" relation Ferrara et al. 2019;Pallottini et al. 2019) and standard assumptions for the dust properties (Sommovigo et al. 2021(Sommovigo et al. , 2022)).2014), we find typical values of δ GDR = 1000 − 3000 for the stellar masses and redshifts obtained for the REBELS galaxies (median ∼ 10 9.7 M ⊙ , z ∼ 7).These δ GDR values are consistent with recent predictions for these stellar masses and redshifts, based on hydrodynamical simulations (DustyGadget; see Fig. 4 from Graziani et al. 2020).
Figure 3 (middle) shows the comparison of both M d estimates and M mol, [CII] , in the available galaxies.We find that both dust mass estimates are consistent and linearly related to the [C ii]-based molecular gas mass estimates within the scatter (∼ 30%), suggesting δ GDR ∼ 1500.While it is not surprise that our molecular gas estimates correlate well with the dust mass estimates from Sommovigo et al. (2022), as both use [CII] line information, the comparison suggest an overall consistency among the different methods and measurements, including the independent method used by Ferrara et al. (2022).

Estimates from dynamical mass
An alternative way to compute the molecular gas mass is provided by considering the mass budget in each galaxy.The dynamical mass within half-light radius will be given by M dyn (< r e ) = 0.5(M stars + M gas ) + M DM (< r e ) (e.g., Daddi et al. 2010).Assuming that the mass budget within r e will be dominated by baryons, thus the dark matter fraction is negligible and that the ISM is dominated by the molecular gas (M gas ≃ M mol ), we can use the dynamical and stellar masses to obtain M mol,dyn .
Unfortunately, our [C ii] observations can spatially resolve only mildly some of the sources, and therefore we do not have information about their resolved [C ii] kinematics.
For simplicity, we estimate the dynamical masses of the REBELS sources using equations 10 and 11 in Bothwell et al. (2013).In the first case, the gas is distributed in a virialized spherical system, e.g. a compact starburst, with M dyn = 1.56 × 10 6 σ2 R, where σ is the velocity dispersion (σ = ∆v FWHM /2.35) and R is the source radius.In the second case, the gas is distributed in a disk with M dyn = 4 × 10 4 ∆ 2 FWHM R/sin(i) 2 , where i is the inclination angle of the disk.In both cases, we estimate the source radius from the circularized estimate R = √ a maj × a min , where a maj and a min are the semi-major and semi-minor axes, respectively.We restrict the estimates on R by considering only sources that have been resolved at least marginally in a given axis with a significance of 2, or a/δa > 2. In unresolved cases, we consider the source to be small and simply fix their radius to 2 kpc, based on previous size measurements in galaxies at similar redshifts (Fujimoto et al. 2020;Herrera-Camus et al. 2021) and for the REBELS sample (Fudamoto et al. 2022).This translates into 0.36 ′′ at z = 7 (or a source diameter of ∼ 0.7 ′′ ), corresponding to ∼ 25% of the observed beam size (1.6 ′′ ).In the disk case, we adopt a mean inclination angle sin(i) = π/4 ≈ 0.79.The virialized spherical geometry dynamical mass estimator will yield a mass 4.4 larger than the one obtained by assuming a disk-like gas distribution for the same source size, linewidth, and an average inclination parameter.Figure 3 (right) shows the comparison between the [C ii]based molecular gas mass (M mol, [CII] ) and the one estimated from the dynamical and stellar masses (M mol,dyn−stars ).From Fig. 1 we note that the stellar masses are systematically smaller than the derived molecular gas masses and thus have little effect in this comparison.Despite the number of assumptions and uncertainties in the parameters, we find overall agreement between both molecular gas estimates in the case of the disk geometry, with most of the galaxies falling within ±0.3 dex from a linear relation.In the case of virialized geometry, we find a systematic offset between the [C ii] and dynamical estimates and large scatter.The reason for such large scatter is likely due to the simple geometry assumption for very complex star-forming systems at these redshifts, where most of them are not expected necessarily to be either rotating disks or fully spherical systems.
In summary, we find that all the independent checks obtained by comparing the [C ii]-based molecular gas mass estimates with the IR luminosity, the dust masses, and the dynamical masses are consistent and support the application of the Zanella et al. (2018) calibration to estimate the molecular gas mass, assuming that the molecular gas-phase dominates the baryon budget at z > 6 (which might not be the case; Heintz et al. 2022;Vizgan et al. 2022).

Analysis and discussion
The relationship between the SFR and M mol in galaxies is usually termed the global "Schmidt-Kennicutt" (SK) relation or "star formation" law.This relationship probes how star formation activity is linked to the molecular gas supply.The ratio between the two quantities is usually called the molecular gas depletion timescale, t dep = M mol /SFR, which is seen as the amount of time left before the galaxy runs out of molecular gas at the current rate of star formation.Per definition, this quantity does not account for possible feedback processes in the galaxy (outflows) and interaction with the IGM (inflows).The fraction of baryonic mass in the galaxy that is in the form of molecular gas is usually described by the molecular gas fraction, f mol = M mol /(M mol + M stars ) or as the molecular gas to stellar mass ratio, µ mol = M mol /M stars .
Several studies in the last decade have found that these quantities are intrinsically linked to other galaxy properties such as their stellar mass, specific SFR (sSFR), and/or their distance to the star-forming main sequence at a given redshift, ∆ MS = sSFR/sSFR MS (z), yielding the differentiation of different modes of star-formation as a function of ∆ MS (e.g., Sargent et al. 2014).Furthermore, t dep and µ mol were found to evolve with redshift (and stellar mass).While mild evolution is observed for the average t dep for galaxies at z < 3, strong evolution is seen for the average µ mol (Tacconi et al. 2018;Liu et al. 2019, e.g.,).Recent determinations of the molecular gas mass in galaxies from the ALPINE survey permitted extending these results to the redshift range 4 − 6, indicating little evolution on the average values of both t dep and µ mol from z ∼ 4 − 6 to z ∼ 3 (Dessauges-Zavadsky et al. 2020), being therefore interesting to check if similar results are found in an independent sample, and at higher redshift.
In the following, we analyze and put in context the ISM properties of the REBELS galaxies compared to previous observations of distant galaxies.The comparison samples used are briefly described in Section 2.3 and in the following analysis were restricted to the stellar mass range of the REBELS sample.This filter is important for the lower redshift samples which contain a significant number of more massive galaxies.

Scaling relations
Figure 4 shows the relationships formed by t dep and µ mol with ∆ MS for the REBELS sources, compared to the samples and prescriptions in the literature for the expected scaling relations at z = 7 (Tacconi et al. 2018;Liu et al. 2019;Scoville et al. 2017), computed at the median stellar mass of the REBELS sample, 10 9.5 M ⊙ at z = 7.We note that the Scoville et al. (2017) measurements provide the total ISM masses (HI and H 2 ), and thus, by construction, yield larger masses compared to molecular gas fractions only.For the MS, we use the prescription from Schreiber et al. (2015).We note that these prescriptions for µ mol , ∆ MS and the MS correspond to extrapolations based largely on CO and dust observations of galaxy samples z < 4, as no similar observations are available for galaxies at z > 6.
The REBELS sample does not show a strong dependency of t dep on ∆ MS as opposed to lower redshift samples and expectations from scaling relations.This could be partly due to IR luminosity limits (see Fig. 3).
We find that all galaxies are consistent within the scatter with the expected low values for t dep at these redshifts.A similar trend is followed by the ALPINE sample, indicating that this trend is not associated with sample selection.Thus, the lack of dependency on ∆ MS is likely related to the fact that most sources fall in the MS (∆ MS ∼ 0) and to the scatter in this relationship expected at high-redshift.Moreover, the location of both REBELS and ALPINE samples in the t dep vs. ∆ MS plane agrees well with the expectation from the Tacconi et al. (2018) and Scoville et al. (2017) relations at these redshifts.
A significant difference in the dependency of µ mol with ∆ MS is seen between the high and low redshift samples and among the various prescriptions for the scaling relations.We find that for the REBELS and ALPINE samples, µ mol consistently increases with distance from the MS, as expected (i.e.galaxies that have greater SFR per unit stellar mass have larger M mol per unit stellar mass).Furthermore, µ mol values are systematically larger for the REBELS and ALPINE samples compared to galaxies at z ∼ 1−3 (i.e.galaxies at higher redshifts are richer in molecular gas) for a given stellar mass.However, a significant difference is observed when comparing the various extrapolations for scaling relations at these redshifts, being the Tacconi et al. (2018) prescription the one that predicts values closer to our measurements.This is similar to previous findings from Dessauges-Zavadsky et al. (2020).
Overall, the REBELS (and ALPINE) galaxies mimic the scaling relation trends expected based on lower redshift samples.However, we find large discrepancies among the various models at these redshifts.This is particularly evident for the dependence of µ mol on ∆ MS , where model predictions differ by nearly two orders of magnitude for MS galaxies.As mentioned earlier, an important source of uncertainty, which affects the values of ∆ MS and µ mol , is the lack of sensitive constraints on the rest-frame optical SED.This will affect the stellar mass estimates and their derived uncertainties.

Evolution of t dep and µ mol
Figure 5 shows the measured dependency of t dep and µ mol with redshift for the REBELS sources, compared with the literature samples and prescriptions for the expected scaling relations (Tacconi et al. 2018;Liu et al. 2019;Scoville et al. 2017) for MS galaxies with a stellar mass of 10 9.75 M ⊙ .To check for evolutionary trends, we divided the sample into ranges of redshift, grouping the ALPINE and REBELS samples into two bins each (4.3 < z < 4.7 and 5.0 < z < 6.0 for ALPINE; 6.4 < z < 7.0 and 7.0 < z < 7.8 for REBELS) and dividing the literature galaxy samples at lower redshifts into six bins.With this, we computed the median among all galaxies in each bin.Uncertainties in each median point taken are obtained by computing the quartile values for the subsamples in each bin.These values are found to be in excellent agreement with weighted average values for the subsamples in each bin, and with previously reported values for these mass ranges (Tacconi et al. 2018;Dessauges-Zavadsky et al. 2020, e.g.,).
We find that the REBELS galaxies exhibit a mild (or no) evolution of t dep from z = 7 to z ∼ 4 − 6 (and at most redshifts), as expected from previous observations and models, yielding a median t dep = 0.5 +0.26  −0.14 Gyr at z = 7, where the uncertainty corresponds to the interquartile range.Little scatter is seen among the REBELS sample, and there is an overall consistency with the Tacconi et al. (2018) and Scoville et al. (2017) prescriptions.The short t dep values found here require that there is additional  gas accretion for the next 3 Gyr to sustain the high star formation activity until z ∼ 2.
Similarly, we find a median µ mol = 2.6 +4.1 −1.4 at z = 7, equivalent to a molecular gas fraction f mol = 0.73 +0.17 −0.11 .The scatter seen here is mostly due to the wide range of stellar masses of galaxies per redshift bin, since there is a strong dependence of µ mol on M stars (see Fig. 8 bottom panel from Dessauges-Zavadsky et al. 2020).This scatter could also be caused by the fact that the REBELS galaxies are likely younger than lower redshift objects at matched stellar mass and thus caught at earlier stages in their assembly process, or simply because of the uncertainties in stellar mass measurements, given the lack of rest-frame optical photometry.
In this case, the scaling relation model prescriptions diverge by one order of magnitude, with the Tacconi et al. (2018) model being the one that agrees best with the data.Since the Scoville et al. (2017) measurements yield the total ISM masses (HI and H2), they yield larger gas fractions when compared to molecular gas fractions only.

Potential effects of our adoption of stellar masses
As discussed in Section 2.2, we have adopted stellar masses obtained through SED fitting using non-parametric SFHs as this approach helps overcome the outshining problem and for consistency with a parallel study (Heintz et al. 2022).For the sample considered in this study, these non-parametric stellar mass estimates are on average ∼0.35 dex higher than those computed using parametric constant SFHs.
This selection has only a minor effect on the main conclusions drawn in this and in the previous section.The choice of non-parametric stellar masses make the REBELS galaxies mostly consistent with the star-forming main-sequence, whereas using constant SFHs yields a slightly higher fraction of starburst galaxies, or above 0.5 dex from the expected star-forming main sequence at z ∼ 7. We find that only two sources are above this threshold for non-parametric stellar masses versus 9 sources for parametric ones.
The choice of non-parametric stellar masses has little effect on the computed ∆ MS .The median ∆ MS for non-parametric SFHs is only 0.1 dex lower than that computed with parametric constant SFHs.Similarly, our choice of stellar mass measurements does not affect the t dep estimates; however, it has a more important effect on µ mol .The median µ mol obtained for non-parametric stellar masses is a factor 2.4 lower than that obtained with parametric ones.Choosing the parametric stellar masses would still put the REBELS galaxies within the expected scaling relations in the log(µ mol ) vs. ∆ MS plot (Fig. 4).However, they would now be more consistent with the Liu et al. ( 2019) prescription for µ mol at z = 7.This would be in contradiction to the consistently low value of t dep at z ∼ 7, which is consistent with the prescription of Tacconi et al. (2020).

Cosmic density of molecular gas
Great observational efforts have been made in recent years to measure the CO luminosity function (LF) at different redshifts and, thereby, the cosmic density of molecular gas ρ mol .Initial attempts to measure the CO LF at high redshift were done with the Very Large Array (VLA) in an overdense field at z ∼ 1.5 (Aravena et al. 2012).However, the first dedicated survey for this purpose used the IRAM Plateau de Bureau Interferometer (PdBI) to perform a spectral scan over frequency and thus conduct a blind search for CO line emission in the Hubble Deep Field North (< 1 arcmin 2 ; Walter et al. 2014;Decarli et al. 2014).Later on, CO line and dust continuum surveys with ALMA were performed to this end in the Hubble UDF, including the ASPECS pilot (∼ 1 arcmin 2 , Walter et al. 2016;Aravena et al. 2016b,a;Decarli et al. 2016) and large programs (∼ 5 arcmin 2 , Decarli et al. 2019;González-López et al. 2020;Aravena et al. 2019Aravena et al. , 2020;;Decarli et al. 2020;Boogaard et al. 2020).Similar efforts were done with the VLA through the CO Luminosity Density at High-Redshift (COLDz) survey (Riechers et al. 2019) covering larger areas at shallower depths for lower-J CO transitions.A parallel approach has taken advantage of targeted CO and dust continuum surveys, as each independently observed pointing would essentially yield a portion of a cosmological deep field, and combining it with statistical methods to assess the completeness and compute the cosmological volumes covered.This yielded estimates for the CO LF from the PHIBSS (Lenkić et al. 2020), A3COSMOS (Liu et al. 2019), and the ALMACAL surveys (Klitsch et al. 2019).Similar estimates have been obtained from dust continuum observations (Scoville et al. 2017;Magnelli et al. 2020), and most recently using the ALPINE [C ii] survey (Yan et al. 2020).
It is interesting to compare these measurements with the value of ρ mol at z = 6.5 − 7.5 from the REBELS survey.Following previous studies (Yan et al. 2020;Heintz et al. 2021Heintz et al. , 2022)), we compute ρ mol from the [C ii] luminosity density, L, obtained for galaxies at z ≈ 7 from the REBELS survey by Oesch et al. (in preparation).In this study, the [C ii] luminosity density is found by integrating the [C ii] luminosity function down to a limit of log(L The latter is obtained from the UV LF (Bouwens et al. 2021), using an empirical L [CII] − L UV relation, derived for galaxies at the same redshift (for details, see Oesch et al., in prep.; and also Barrufet et al. 2023).With this, we obtain a [C ii] luminosity density log(L/[L ⊙ Mpc −3 ]) = 4.85 +0.25  −0.20 for the REBELS sample.Now, this can be converted to a molecular gas mass density using the calibration described above, α Zanella et al. 2018), yielding log(ρ mol /(M ⊙ Mpc −3 )) = 6.34 +0.34  −0.31 at z = 6.5 − 7.5.This result includes the uncertainties associated with the α [CII] calibration (Zanella et al. 2018).For the ALPINE data, we use the updated calculations of the [C ii] luminosity density from Heintz et al. (2022).
Figure 6 shows the result of this comparison.The measurements obtained from the different methods and tracers agree that ρ mol resembles the evolution of the cosmic SFR density (scaled by a typical molecular gas depletion timescale of 0.5 Gyr, shown by the solid green line in Fig. 6; e.g., Decarli et al. 2020), although the studies by A3COSMOS (Liu et al. 2019) and Scoville et al. (2017) find different normalizations.
However, the few measurements obtained at z ∼ 4 − 6, including the PHIBSS (Lenkić et al. 2020) and COLDz (Riechers et al. 2019) CO-based and the ALPINE [C ii] based estimates (Yan et al. 2020;Heintz et al. 2022), appear to agree, indicating a continuous decline of ρ mol at increasing redshifts.The REBELS measurements agree with the trend of ρ mol declining at higher redshifts.These results indicate that ρ mol increases by an order of magnitude from z ∼ 7 to z ∼ 4, and that beyond z = 6, the values of ρ mol reach levels similar and below the current ρ mol (at z = 0) of log(ρ mol /(M ⊙ Mpc −3 )) ∼ 7.

Conclusions
We present molecular gas mass estimates and properties for the sample of galaxies at z ∼ 7 unveiled by the REBELS large program survey, using the [C ii] line emission as a proxy for the molecular ISM.The main conclusions of this study are as follows.
1. Comparison of the [C ii]-based molecular gas mass estimates with independent estimates and proxies, including the IR luminosity (as a proxy for the molecular gas mass through the L IR − L ′ CO relation), dust masses, and dynamical mass estimates yield reasonable agreement when applied to the REBELS sample at z ∼ 7. Furthermore, the matched stellar mass and SFR distributions of REBELS and ALPINE surveys suggest that the calibrations to obtain M mol from L [CII] apply for massive main-sequence galaxies at z ∼ 7 and z ∼ 4 − 6. 2. The molecular to stellar mass ratios (µ mol ) and molecular gas depletion timescales (t dep ) appear to follow the standard scaling relations with distance to the main sequence ∆ MS for REBELS galaxies at z ∼ 7.While µ mol increases toward The green dotted curve shows the cosmic SFR density (Madau & Dickinson 2014), scaled by a typical molecular gas depletion timescale of 0.5 Gyr (Tacconi et al. 2020).REBELS measurements extending out to z ∼ 7 indicate a sustained decline of ρ mol at earlier epochs larger ∆ MS , t dep is observed to stay relatively constant.The location of the REBELS galaxies at z ∼ 7 in these plots agrees with the regions occupied by ALPINE galaxies and the distant end of lower redshift samples at matched stellar masses.3. The median t dep and µ mol stay almost constant in the range z ∼ 2 − 7 for galaxies with stellar masses 10 8.5−10.5 M ⊙ .4. We measure a value of the cosmic density of molecular gas at z = 6.5−7.5 of log(ρ mol /(M ⊙ Mpc −3 )) = 6.34 +0.34 −0.31 , indicating an increase of an order of magnitude from z ∼ 7 to z ∼ 4.
The advent of infrared facilities such as the JWST and Roman will allow us to explore samples of hundreds to thousands of galaxies at z > 6, opening the study of the formation processes of galaxies through the end of the epoch of reionization.Investigating the cold ISM of such early galaxies will require observations with submillimeter facilities as ALMA.However, due to the dimming of traditional ISM tracers such as CO line and dust continuum emission due to cosmological distances, lower metallicities, and increasingly low-contrast to the cosmic microwave background at higher redshift (da Cunha et al. 2013), measurements of the ISM gas content (a key ingredient to study the galaxy assembly) will be increasingly difficult.The [C ii] line emission has therefore come to save the day, enabling relatively accurate measurements of M mol for massive systems in the early universe.Calibration of the [C ii] line for larger samples of galaxies, as well as other fine-structure lines, will be essential for future galaxy formation studies.

Fig. 1 .
Fig. 1.Distribution of properties of the REBELS galaxies (yellow filled histogram) compared to the ALPINE galaxies (blue empty histogram).From left to right, we show the distribution of redshift, stellar mass, SFR (UV+IR), and [C ii]-based molecular gas masses.Stellar masses are derived from non-parametric SED fitting as described in the text (see also, Topping et al. 2022).The histogram for the ALPINE sample has been normalized to the maximum number of REBELS galaxies.No comparison of the redshift distributions is provided, since, by construction, the samples cover different redshift ranges.

Fig. 2 .
Fig. 2. SFR(UV+IR) vs. stellar mass diagram for the REBELS sources.Blue triangles show sources detected in [C ii] emission without IR continuum detection.Red color circles show sources detected in [C ii] and IR emission.The dotted and dashed lines represent models of the starforming main sequence at z = 7(Schreiber et al. 2015;Iyer et al. 2018) and at z ∼ 5.5(Khusanova et al. 2021).Most of the REBELS galaxies are consistent with being within ±0.3 dex of the expected main sequence at z = 7 (seeTopping et al. 2022).

Fig. 3 .
Fig. 3. (Left:) L IR versus L ′ CO(1−0) for the REBELS galaxies.The filled orange circles show the result of assuming a "Milky Way" conversion factor, α CO = 4.6, while the open circles show the case of α CO = 1.0.The linear fit to CO(1-0) observations of local spiral and disk galaxies at high redshift compiled by Dessauges-Zavadsky et al. (2015) is shown as a dashed blue line (Daddi et al. 2010; Aravena et al. 2016c, see also,).The dotted line shows the same line, with a factor of ±0.3 dex, representing the typical scatter.(Middle:) Comparison of the [C ii]-based molecular gas masses for the REBELS galaxies with the dust mass estimates from the models from Ferrara et al. (2022) and Sommovigo et al. (2022) are shown as blue diamonds and orange circles, respectively.Dotted lines highlight curves of constant gas-to-dust ratios.(Right:) Comparison of the [C ii]-based molecular gas masses (M mol,[CII] ) with estimates obtained from the dynamical and stellar masses (M mol,dyn−stars ).Derivations of the dynamical mass based on an assumed spherical (virial) and disk geometries are shown as orange-filled and empty circles, respectively.The dotted lines represent the 1:1 relation and the range ±0.3 dex.

Fig. 4 .
Fig. 4. ISM scaling relations for the REBELS galaxies compared to literature.(Left:) Molecular gas depletion timescale (t dep ) as a function of offset to the MS (∆ MS ).(Right:) Molecular gas to stellar mass ratio (µ mol ) as a function of offset to the MS (∆ MS ) .In both panels, orange-filled circles show the REBELS galaxies z = 6 − 8 with molecular gas masses derived from [CII] measurements.Green triangles show [C ii] measurements in ALPINE galaxies at z = 4 − 6 (Béthermin et al. 2020).Gray squares and blue crosses show CO and dust measurements from PHIBSS at z < 3 (Tacconi et al. 2018) and additional sources from the literature (see text), respectively.All samples have been restricted to galaxies with M stars = 10 8.5−10.5 M ⊙ to match the range spanned by REBELS galaxies.The solid blue, dashed magenta and black dot-dashed lines represent the prescriptions for the scaling relations computed for the median mass of the REBELS sample, M stars = 10 9.75 M ⊙ and z = 7, from Scoville et al. (2017), Tacconi et al. (2018) and Liu et al. (2019), respectively.

Fig. 5 .
Fig. 5. Evolution of the molecular ISM out to z = 7.5.(Left:) Molecular gas depletion timescale (t dep ) as a function of redshift.(Right:) Molecular gas to stellar mass ratio (µ mol ) as a function of redshift.Symbols are the same as in Fig. 4. Similarly, literature comparison samples are trimmed to have the stellar mass range of the REBELS galaxies, M stars = 10 8.5−10.5 M ⊙ .The black open squares represent the averages in each relevant redshift range.