Cosmic Evolution of the H₂ Mass Density and the Epoch of Molecular Gas

The UKIDSS-UDS field was mapped at 850 μm as part of the SCUBA-2 Cosmology Legacy Survey (Geach et al. 2017). The full details of the data collection, reduction, and map properties are given in Geach et al. (2017). Briefly, the beam-convolved map spans approximately 1 deg² covering the bulk of the multiwavelength coverage of this field, with a uniform (instrumental) noise of σ₈₅₀ = 0.9 mJy beam⁻¹. Geach et al. (2017) estimated the SCUBA-2 confusion limit to be σ_conf = 0.8 mJy beam⁻¹. The beam FWHM is approximately 15″, with a full analytic description of the point-spread function (PSF) given by Geach et al. (2017).

The UDS Data Release 11 (DR11) 12-band matched catalog is K band selected with the 95% completeness limit estimated to be K_AB = 25 mag. The full details of this catalog will be comprehensively provided in O. Almaini et al. (2021, in preparation) and W. G. Hartley et al. (2021, in preparation), and only a summary is given here. The catalog provides photometry in 12 bands (U, B, V, R, $i^{\prime}$ , $z^{\prime}$ , Y, J, H, K, 3.6 and 4.5 μm), where available.

The J, H, and K photometry is taken from the DR11 release of UKIDSS-UDS. The UKIDSS project, described in Lawrence et al. (2007), utilizes the UKIRT Wide Field CAMera (WFCAM; Casali et al. 2007). The photometric system and calibration are outlined in Hewett et al. (2006) and Hodgkin et al. (2009), respectively, and the pipeline processing and science archive are described in M. J. Irwin et al. (2021, in preparation) and Hambly et al. (2008). UKIDSS-UDS covers an area of 0.8 deg², reaching median depths of J = 25.6, H = 25.1, and K = 25.3 (5σ, AB, estimated from 2″ apertures in source free areas; O. Almaini et al. 2021, in preparation).

The B, V, R, $i^{\prime}$ , and $z^{\prime}$ optical imaging is from the Subaru/XMM-Newton Deep Survey, which utilizes Suprime-Cam on the Subaru Telescope (Furusawa et al. 2008). U-band data are from the Canada–France–Hawaii Telescope Megacam instrument (O. Almaini et al. 2021, in preparation) and Y-band imaging is obtained from the VISTA Deep Extragalactic Observations survey (VISTA-VIDEO; Jarvis et al. 2013). The InfraRed Array Camera (IRAC) imaging at 3.6 and 4.5 μm is from the Spitzer UKIDSS Ultra Deep Survey (SpUDS; PI Dunlop), combined with deeper data from the Spitzer Extended Deep Survey (SEDS: Ashby et al. 2013). To expand the coverage to outer regions of the field, shallower data are also used from the SIRTF Wide-area InfraRed Extragalactic survey (SWIRE; Lonsdale et al. 2003).

UDS DR11 provides image masks, with masked regions corresponding to image boundaries, artefacts, and bright stars. We employ the UDS binary mask for “good” regions, which has an unmasked area of 0.64 deg². This binary mask combines the masked regions of the photometry images detailed above (not including the deeper SEDS or SpUDS IRAC images).

We also utilize the subsets feature of the UKIDSS-UDS catalog and for our galaxy sample chose the catalog-defined “good galaxy” subset, which comprises 217,429 sources. These sources have full 12-band photometry and lie within the corresponding “good mask regions, are not cross-talk sources (for which JHK photometry is likely compromised), and are not classified as stars.

UDS DR11 also includes photometric redshifts derived using the code eazy (Easy and Accurate Z_phot from Yale; Brammer et al. 2008). To estimate the photometric redshift for each source the 12-band broadband photometry was fit with a spectral energy distribution template producing a redshift probability distribution (W. G. Hartley et al. 2021, in preparation). We utilize both the maximum-likelihood photometric redshifts and redshift probability distributions provided with UDS DR11. eazy performs well compared to other commonly used photometric redshift codes (e.g., zphot, HyperZ, Rainbow; Dahlen et al. 2013), with the resulting normalized mean absolute deviation between eazy derived photometric redshifts and spectroscopic redshifts found to be only σ_nmad ≈ 0.02 (W. G. Hartley et al. 2021, in preparation).

In Table 1, we present the average observed 850 μm flux densities for our galaxy sample as a function of redshift. We quote the uncertainties to 1σ and include the additional uncertainty due to sample variance. We note that at z ≳ 6 the UDS redshifts are untested. However, as our sample is binned according to the ${ \mathcal P }(z)$ for each source, every galaxy effectively contributes to the flux in each redshift interval. Hence, we show the average 850 μm flux density estimates for our galaxy sample to z = 10.

We sum ${ \mathcal P }(z)$ (which is completeness corrected) in each redshift bin, giving us the galaxy weighting for each redshift interval. The integral of the summed ${ \mathcal P }(z)$ across all redshift intervals should be approximately equal to the total number of galaxies in our sample. We calculate this to be 154,839, which is consistent with our galaxy sample of 153,399 sources (taking into account the completeness corrections). With the summed ${ \mathcal P }(z)$ and taking the area of our sample as the unmasked region of the SCUBA-2 850 μm map (which corresponds to the UDS binary mask for good galaxy regions), we calculate the number density of galaxies as a function of redshift. By combining this with our average flux density (see Table 1, Column 4 for CMB corrected values) we calculate the summed flux density for our galaxy sample in each redshift interval. In Figure 2, we present the number density and summed flux density for galaxies in our sample as a function of redshift. The distribution of the summed flux densities with redshift is broadly comparable to the redshift distribution found for SMGs, which peaks at z ≈ 2 (e.g., Blain et al. 2002; Chapman et al. 2005; Simpson et al. 2014; Miettinen et al. 2017; Zavala et al. 2018), while the number density distribution generally declines with increasing redshift as expected. The difference in the evolution of these distributions demonstrates that our derived flux densities are not biased by the number density of galaxies in each redshift bin.

Zoom In Zoom Out Reset image size

We adopt the approach of S16 and utilizing our flux density measurements, estimate the average molecular gas mass for our galaxy sample in each redshift interval. The full details of this approach are given in S16; however, we provide a brief description here.

The long-wavelength RJ tail of dust emission is nearly always optically thin (τ ≪ 1) and consequently this provides a direct probe of the total dust mass and hence the molecular gas mass. S16 utilize this to obtain an empirically calibrated RJ luminosity-to-gas mass ratio

$$\begin{eqnarray}&&\left(\displaystyle \frac{{M}_{\mathrm{mol}}}{{M}_{\odot }}\right)=\displaystyle \frac{1}{{\alpha }_{850}}\left(\displaystyle \frac{{L}_{850,\mathrm{rest}}}{\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}}\right)\mathrm{for}\ {\lambda }_{\mathrm{rest}}\gtrsim 250\,\mu {\rm{m}},\end{eqnarray} \tag{ 8 }$$

with ${\alpha }_{850}=(6.7\pm 1.7)\times {10}^{19}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}\,{{M}_{\odot }}^{-1}$. The restriction λ_rest ≳ 250 μm is required to ensure that at an observed wavelength of 850 μm the rest-frame emission stays on the RJ tail. S16 demonstrated that this luminosity-to-mass ratio is relatively constant for high-stellar mass (M_stellar = (2 − 40) × 10¹⁰ M_⊙) normal star-forming and star-bursting galaxies, both locally and at high z.

We estimate the average rest-frame 850 μm luminosity density of galaxies as a function of redshift using the average flux densities detailed in Table 1, assuming a mass-weighted dust temperature of 25 K ¹¹ and employing the relation from S16:

$$\begin{eqnarray}&&\begin{array}{l}{L}_{\nu 850}={S}_{\nu }[\mathrm{Jy}]\times \ 1.19\times {10}^{27}\times \ {\left(\displaystyle \frac{\nu (850\mu {\rm{m}})}{{\nu }_{\mathrm{obs}}(1+z)}\right)}^{3.8}\\ \times \ \displaystyle \frac{{\left({d}_{{\rm{L}}}[\mathrm{Mpc}]\right)}^{2}}{1\,+\,z}\times \ \displaystyle \frac{{{\rm{\Gamma }}}_{\mathrm{RJ}}(25,{\nu }_{850\mu {\rm{m}}},0)}{{{\rm{\Gamma }}}_{\mathrm{RJ}}(25,{\nu }_{\mathrm{obs}},z)}\ [\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}].\end{array}\end{eqnarray} \tag{ 9 }$$

The Γ_RJ term in Equation (10) corrects for departures from the RJ ν² dependence as the observed emission approaches the spectral energy distribution (SED) peak in the rest frame and where T_dust is the mass-weighted temperature characterizing the RJ dust emission

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{RJ}}({T}_{\mathrm{dust}},{\nu }_{\mathrm{obs}},z)=\displaystyle \frac{h{\nu }_{\mathrm{obs}}(1+z)/{{kT}}_{\mathrm{dust}}}{{e}^{h{\nu }_{\mathrm{obs}}(1+z)/{{kT}}_{\mathrm{dust}}}-1}.\end{eqnarray} \tag{ 10 }$$

We restrict our estimates of the average rest-frame 850 μm luminosity to z ≲ 2.5 to ensure that rest-frame emission stays on the RJ tail. With the average rest-frame 850 μm luminosity density derived (detailed in Table 2) we use the RJ luminosity-to-gas mass ratio from Equation (8) to estimate the average molecular gas mass as a function of redshift to z ≈ 2.5. This calibration includes a factor of 1.36 to account for the associated mass of heavy elements (mostly helium at 8% by number), so we correct our results by a factor 1/1.36 (M_mol) to obtain ${M}_{{{\rm{H}}}_{2}}$.

Table 2. Redshift Intervals, Average Rest-frame 850μ m Luminosity, Average Molecular Gas Mass, Co-moving Molecular Gas Mass Density as a Function of z

Δz	〈L850,rest〉	$\langle {M}_{{{\rm{H}}}_{2}}\rangle$	${\rho }_{{{\rm{H}}}_{2}}$	${{\rm{\Omega }}}_{{{\rm{H}}}_{2}}$
	[1027 erg s−1 Hz−1 dz−1]	[107 M⊙ dz−1]	[106 M⊙ Mpc−3]	[10−7]
0.25	1.64 ± 0.02	1.80 ± 0.03	1.92 ± 1.80	115.66 ± 108.42
0.75	81.97 ± 3.89	89.95 ± 4.26	17.12 ± 4.39	569.51 ± 146.13
1.25	217.79 ± 18.32	239.01 ± 20.10	23.22 ± 4.09	430.70 ± 75.93
1.75	514.61 ± 64.72	564.77 ± 71.02	34.67 ± 5.32	380.54 ± 58.41
2.25	492.12 ± 85.19	540.08 ± 93.49	20.74 ± 4.60	143.48 ± 31.84

Note. Column 1: midpoint of each redshift interval (Δz); Columns 2 and 3: the average rest-frame 850 μm luminosity and average molecular gas mass for galaxies in each redshift interval; Column 4: molecular gas mass density as a function of redshift and Column 5: molecular gas mass density in terms of the critical mass density. We restrict our results to z ≲ 2.5 to ensure that the observed 850 μm dust emission is tracing the RJ tail.

Download table as:  ASCIITypeset image

Since the summed photometric redshift probability distributions inform us about the galaxy weighting in each redshift interval and the UDS binary mask for good regions gives us the unmasked area of the SCUBA-2 850 μm map, we combine this information with the average molecular gas mass and differential co-moving volume element to estimate the co-moving volume density of molecular gas

$$\begin{eqnarray}&&{\rho }_{{{\rm{H}}}_{2}}={\rm{\Omega }}{\int }_{z-{\rm{\Delta }}z/2}^{z+{\rm{\Delta }}z/2}N(z)\langle {M}_{{{\rm{H}}}_{2}}(z)\rangle \displaystyle \frac{{dV}}{{dzd}{\rm{\Omega }}}{dz}.\end{eqnarray} \tag{ 11 }$$

We present our values for ${\rho }_{{{\rm{H}}}_{2}}(z)$ in Table 2 as a function of redshift, also giving this in terms of the critical mass density ${{\rm{\Omega }}}_{{{\rm{H}}}_{2}}={\rho }_{{{\rm{H}}}_{2}}(z)/{\rho }_{\mathrm{crit}}(z)$. We use a Monte Carlo analysis to calculate the uncertainties for our values of L_850,rest, ${M}_{{{\rm{H}}}_{2}}$ and ${\rho }_{{{\rm{H}}}_{2}}(z)$, first drawing random values for S_ν from a Gaussian distribution where the mean is the average flux density and width the uncertainty on the average flux density, and then drawing values for a mass-weighted T_dust from a Gaussian distribution with a mean of 25 K (corresponding to the constant T_dust assumed by S16) and width 3 K. Observations (e.g., Planck Collaboration et al. 2011; S16) and simulations (e.g., Liang et al. 2018, 2019) find that a mass-weighted T_dust shows little variation with galaxy L₈₅₀ or redshift (see Behrens et al. 2018), and by utilizing a temperature distribution with σ = ±3 K, we recognize this minimal variance in our uncertainty calculations. We use these values to estimate L_850,rest, ${M}_{{{\rm{H}}}_{2}}$ and ${\rho }_{{{\rm{H}}}_{2}}(z)$ from Equations (9), (8), and (11), respectively, for 1000 runs, with the uncertainty being taken as the standard deviation across these trials.

We note that the galaxy sample we use to derive the results in Tables 1 and 2 includes all galaxies in the good galaxy subset of the UDS DR11 catalog, regardless of the reliability of photometric redshifts for individual sources. If we apply a χ² cut to exclude galaxies with the least reliable redshifts (omitting galaxies with a reduced χ² value for the photometric redshift of >10) and repeat the process outlined in Sections 3 and 4 above, we find a less than 2% variation in our results with estimates consistent with those in Table 2 within the uncertainties.

We plot ${\rho }_{{{\rm{H}}}_{2}}(z)$ in Figure 3, compared to direct CO line estimates from ASPECS (Decarli et al. 2016b, 2020), COLDZ (Riechers et al. 2019), and VLASPECS (Riechers et al. 2020), as well as values derived using far-IR and UV photometry (Berta et al. 2013). We fit a function of the same form as the SFR density function presented in Madau & Dickinson (2014) to the log of our results, and derive the best-fit parameters for our data using emcee. This yields

$$\begin{eqnarray}&&\begin{array}{l}{\mathrm{log}}_{10}\left(\displaystyle \frac{{\rho }_{{{\rm{H}}}_{2}}}{{M}_{\odot }{\mathrm{Mpc}}^{-3}}\right)=(6.59\pm 0.30)\\ \,\times \,\displaystyle \frac{{\left(1+z\right)}^{0.38\pm 0.16}}{1+{\left[(1+z)/5.57\pm 1.69\right]}^{1.78\pm 0.61}},\end{array}\end{eqnarray} \tag{ 12 }$$

which we plot in Figure 3. Our results show a peak ${\rho }_{{{\rm{H}}}_{2}}(z)$ at z ≈ 2 mirroring existing constraints.

Zoom In Zoom Out Reset image size

Using an alternative approach we derive ${\rho }_{{{\rm{H}}}_{2}}(z)$ from first principles using the halo mass function from Murray et al. (2013) and assuming a constant halo mass range of 10^11.5–10¹⁵ M_⊙. We estimate the molecular gas mass density as a function of halo mass (for redshifts 0 ≤ z ≤ 7) using the stellar-halo mass ratio from Moster et al. (2013) and the ISM-stellar mass relation from Scoville et al. (2017). The ISM-stellar mass relation is calibrated using a sample of high mass galaxies (M_stellar ≳ 10¹⁰ M_⊙), therefore, we adopt a halo mass range for which the corresponding stellar masses are comparable with the Scoville et al. (2017) calibration sample. Integrating these estimates with respect to halo mass gives the total molecular gas density as a function of redshift, which we present in Figure 3.

We also estimate ${\rho }_{{{\rm{H}}}_{2}}(z)$ from the SFR density, ${\rho }_{\dot{{M}_{\star }}}(z)$, assuming a corresponding volume averaged star formation efficiency, $\eta (z)={\rho }_{{{\rm{H}}}_{2}}(z)/{\psi }_{\star }(z)$. We use the functional fit of Wilkins et al. (2019), a recalibration of the well-known Madau & Dickinson (2014) cosmic star formation history. We make the assumption that η(z) is constant and that the total molecular gas mass per galaxy can be related to ongoing star formation as $\xi {M}_{{{\rm{H}}}_{2}}=\mathrm{SFR}/\epsilon$ (Geach & Papadopoulos 2012). Here, ξ is the ratio of dense, actively star-forming molecular gas to the total molecular reservoir with ξ ≈ 0.04 for quiescent disks and ξ > 0.5 for starbursts (e.g., Papadopoulos & Geach 2012), while the factor describes the rate at which the dense molecular gas forms stars. Figure 3 shows the predicted ${\rho }_{{{\rm{H}}}_{2}}(z)$ inferred from the Wilkins et al. (2019) fit, assuming a constant “average” η(z) = 0.3 Gyr corresponding to ξ = 0.1 and = 37 Gyr⁻¹ (e.g., Geach & Papadopoulos 2012).

This value for η(z) is similar to the typical values of t_dep (${M}_{{{\rm{H}}}_{2}}/\mathrm{SFR}$) ≈1 Gyr (e.g., Tacconi et al. 2018) for main-sequence galaxies. Tacconi et al. (2018) found a relatively weak dependence of t_dep with redshift, t_dep ∝ (1 + z)^−0.57 to z < 2.5, consistent with our picture of a common mode of star formation in normal galaxies, at least out to the peak epoch. We use this relation from Tacconi et al. (2018) and the Wilkins et al. (2019) fit to derive an estimate of ${\rho }_{{{\rm{H}}}_{2}}(z)$, which incorporates a weakly evolving star formation efficiency. We present our predicted ${\rho }_{{{\rm{H}}}_{2}}(z)$ in Figure 3.

The S16 RJ luminosity-to-gas mass ratio has been shown to provide molecular gas mass estimates accurate to within a factor of around 2 when compared with measurements made via direct CO ( J = 1 → 0) line observations (e.g., Scoville et al. 2017; Kaasinen et al. 2019), with variations in the dust emissivity index, temperature, and gas-to-dust ratios being accountable for the deviations. This factor of 2 accuracy is based on samples of galaxies with high-stellar masses M_stellar = (2–40) × 10¹⁰ M_⊙ as these galaxies are likely to have near-solar metallicity (Tremonti et al. 2004). This avoids probing low metallicity sources for which the dust-to-gas abundance ratio is likely to drop or the CO gas fraction is low (Bolatto et al. 2013). In Figure 4, we show the marginalized stellar mass estimates for UDS galaxies (O. Almaini et al. 2021, in preparation) and the corresponding 95% stellar mass completeness (derived using the method of Pozzetti et al. 2010, A. Wilkinson et al. 2021, in preparation) for the UDS catalog. As can be seen in Figure 4, our galaxy sample includes a proportion of galaxies with stellar masses lower than those used to derive the S16 RJ luminosity-to-gas mass ratio, with these sources being more abundant in lower redshift bins.

Zoom In Zoom Out Reset image size

The dust-to-gas relation has been found to be relatively consistent for nearby galaxies with M_stellar > 10⁹ M_⊙ (e.g., Groves et al. 2015), but drops for galaxies with lower stellar masses (hence lower metallicities). Cosmological galaxy formation simulations have shown that deviations from this relation become significant (≳0.5 dex) at L₈₅₀ ≲ 10²⁸ erg s⁻¹ Hz⁻¹ in the redshift range of 0 < z < 9.5 (e.g., Privon et al. 2018). As shown in Figure 4, at z ≳ 1 the majority of our sample are likely to have M_stellar > 10⁹ M_⊙ and Table 2 shows that the mean rest-frame 850 μm luminosity for our sample in all but the lowest redshift interval (Δz = 0.25) is L_{850 μm} ≳ 10²⁸ erg s⁻¹ Hz⁻¹. Therefore, while the RJ luminosity-to-gas mass ratio has been calibrated on high-stellar mass galaxies (M_stellar = (2 − 40) × 10¹⁰ M_⊙), we make the assumption that applying this calibration to our sample at redshifts z ≳ 0.5 is likely to result in comparable uncertainties (i.e., a factor of 2). However, in the redshift bin Δz = 0.25 the mean rest-frame luminosity is 〈L_{850 μm}〉 = 1.6 × 10²⁷ erg s⁻¹ Hz⁻¹, and as such, our results are likely to be underpredicted by ≳0.5 dex (e.g., Privon et al. 2018) due to the abundance of lower mass (low metallicity) galaxies in this redshift interval.

In Figure 3, we compare our results, which are revised to account for the influence of the CMB, to those from direct CO line surveys (e.g., Decarli et al. 2016b, 2020; Riechers et al. 2019, 2020). We limit our discussion to the results from these surveys at z ≲ 2.5 as we are restricted to this redshift range due to the wavelength of our observations (λ_obs = 850 μm). While the results from ASPECS/ASPECS LP (Decarli et al. 2016b, 2020) and COLDz (Riechers et al. 2019, 2020) were not corrected for the influence of the CMB, we note that this is not necessary at z ≲ 4.5 as the effect of the CMB on measurements of molecular gas mass density from direct CO line observations is minimal (≲15%, e.g., Decarli et al. 2019) and as such this does not impact on our analysis here.

We find that our results are broadly consistent with the estimates from direct CO line surveys within uncertainties and show notably good agreement with results obtained through observations of the ground-state CO line (ASPECS LP at Δz ≃ 0.25 and COLDZ at Δz ≃ 2.25). Albeit, we caution that our results at Δz ≃ 0.25 are likely underestimated due to the abundance of low stellar mass galaxies in this redshift bin.

When compared with the ASPECS LP survey (Decarli et al. 2020) our results generally trace the lower boundaries of their estimates between 0.75 ≲ z ≲ 1.75. The S16 RJ luminosity-to-gas mass ratio is calibrated using the ground-state CO( J = 1 → 0) line, whereas at z > 0.75 the ASPECS LP results are derived from observations of higher state excitation CO lines. Therefore, this offset could be explained by the uncertainties associated with translating higher excitation CO lines observations to ground-state CO( J = 1 → 0) luminosities. However, in the redshift interval Δz = 1.25, even if the extreme case of thermalized gas is assumed, our upper estimate (taking into account 1σ uncertainties) falls a factor of 1.42 below the lower boundary of the ASPECS LP survey (e.g., 38.90 × 10⁶ M_⊙ Mpc⁻³, Table A3, Decarli et al. 2019). As such the uncertainties in CO line ratios do not fully account for the offset we see.

Building on previous studies, Liu et al. (2019) derived the molecular gas mass density using a data set comprised of ≈700 ALMA continuum detected galaxies and ≈1000 galaxies with CO observations (taken from the literature). To derive molecular gas masses for the continuum detected galaxies Liu et al. (2019) employed the Hughes et al. (2017) luminosity-to-gas mass calibration. Liu et al. (2019) estimated an SMF (stellar mass function) integrated molecular gas mass density based on the SMF integrated to M_stellar = 10⁹ M_⊙ and using a gas fraction function derived from their composite sample of ≈1700 galaxies. Their results trace the upper boundaries of the molecular gas mass density derived from the most recent blank field CO line surveys (e.g., Decarli et al. 2019; Riechers et al. 2019) and are ≈1 dex higher than our estimates. This offset with our estimates of molecular gas mass density could be in part due to assumptions made in Liu et al. (2019) to derive an SMF integrated molecular gas mass density (i.e., that all star-forming galaxies are on the main sequence) or potentially differences in sample selection. The majority (≈800) of the CO detected sources in the Liu et al. (2019) composite sample are in the local universe (i.e z < 0.3), so at z > 0.3 their data set is dominated by ALMA continuum detected galaxies, which are preferentially massive and dust rich (and hence, using a dust-to-gas mass conversion, gas rich). In contrast our galaxy sample is NIR selected and as such our selection is less likely to sample these luminous dust-rich SMGs, with previous studies finding that ≈20% of SMGs are missed in optical/new-IR surveys (e.g., Dudzevičiūtė et al. 2020). Magnelli et al. (2020) use a stacking approach to measure the co-moving gas mass density of a sample of 555 NIR selected galaxies, with galaxies split into bins of z and M_stellar. Their stacking method accounts for the metallicity of galaxies in these bins (inferred using the stellar mass–metallicity from Tacconi et al. 2018) and is equivalent to the S16 calibration at solar metallicity. Our results trace the lower boundaries of their estimates, but are inconsistent (within 1σ uncertainties) at z > 1. This discrepancy may be in part due to our method not accounting for the metallicity of low mass galaxies, resulting in an underestimation of the molecular gas mass for galaxies with M_stellar < 10⁹ M_⊙. However, if this was the sole reason for the difference in our results, we would expect this to have more of an impact at z < 1 where this effect will be more prominent. We caution that our results also rely solely on photometric redshifts, which despite the high-quality 12-band photometry of the UDS catalog, cannot compete with the accuracy of redshifts derived via spectroscopic surveys. By utilizing the redshift probability distributions in our 3D stacking approach, we aim to provide mitigation against these uncertainties. However, while our estimates rely exclusively on the use of photometric redshifts, the results obtained in both ASPECS LP (Decarli et al. 2020; Magnelli et al. 2020) and Liu et al. (2019) benefit from the inclusion of sources with more reliable spectroscopic redshifts. This may also play a part in deviations seen when we compare our estimates with these previous surveys.

In order to present the most complete view of the evolution of the molecular gas mass density, we stack all sources in our NIR selected sample, including counterparts, to the bright submillimeter sources in the SCUBA-2 UDS map. To test the contribution of these galaxies to our results, we repeat our stacking analysis with the SCUBA-2 UDS source subtracted map. As expected excluding the ∼1000 UDS submillimeter sources reduces the average observed 850 μm flux in our redshift intervals, which propagates to our estimate of the co-moving molecular gas mass density. At z < 1.5 the exclusion of the UDS submillimeter sources has a minimal impact on our estimates of ${\rho }_{{{\rm{H}}}_{2}}(z)$ and these remain consistent within the 1σ uncertainties. However, at z > 1.5 our estimates of ${\rho }_{{{\rm{H}}}_{2}}(z)$ drop by a factor of 2.05 and 2.33 in the redshift intervals Δz = 1.75 and Δz = 2.25, respectively. This coincides with the peak number density of SMGs at z ≈ 2. This indicates that approximately 50% of the molecular gas mass density at the peak of the SFR density is locked in dust-rich SMGs. We note that our inferred contribution of SMGs is also likely underestimated as we expect that approximately 20% of SMGs are undetected in our NIR selected sample (e.g., Dudzevičiūtė et al. 2020). Our finding is in keeping with that of Zavala et al. (2021) who found that bright SMGS (L_IR > 10¹² L_⊙) dominate the obscured SFR density at z ≈ 2 and also Magnelli et al. (2020) who found that the bulk of dust and gas in galaxies is locked in massive star-forming galaxies.

We have added further valuable constraints to this picture of cosmic molecular gas evolution using two alternative approaches.

We estimate ${\rho }_{{{\rm{H}}}_{2}}(z)$ from the halo mass function (Murray et al. 2013) assuming a constant halo mass range of 10^11.5–10¹⁵ M_⊙, and using the stellar–halo mass ratio from Moster et al. (2013) and the ISM–stellar mass relation from Scoville et al. (2017). For the latter relation we make the assumption that all galaxies are on the star-forming main sequence (e.g., sSFR/sSFR_MS = 1). The evolution of our halo mass derived ${\rho }_{{{\rm{H}}}_{2}}(z)$ (shown in Figure 3) follows a similar shape to the SFR density, rising to a peak at 1 ≲ z ≲ 3 and decreasing to the present day. The minimum halo mass we assume corresponds to stellar masses of ≈10^9.5 M_⊙ (e.g., Moster et al. 2013), equivalent to the lowest stellar masses probed in the ASPECS LP survey (Decarli et al. 2020). Our ${\rho }_{{{\rm{H}}}_{2}}(z)$ estimates show good agreement with the ASPECs/ASPECS LP surveys (Decarli et al. 2016a, 2020) at z ≳ 0.6. However, as shown in Figure 3 at z ≲ 0.6 our estimate of ${\rho }_{{{\rm{H}}}_{2}}(z)$ lies above the lowest redshift bins from the ASPECS LP survey (Decarli et al. 2020) and is ≈1 dex higher than our estimate of ${\rho }_{{{\rm{H}}}_{2}}(z)$ derived from measurements of observed 850 μm flux. To obtain an estimate of ${\rho }_{{{\rm{H}}}_{2}}(z)$ from the halo mass function we make the assumption that all galaxies are star-forming. As such our estimate of ${\rho }_{{{\rm{H}}}_{2}}(z)$ derived from the halo mass function can be seen as an upper limit ${\rho }_{{{\rm{H}}}_{2}}(z)$ for the stellar mass range sampled. It follows that we see a more significant deviation between our estimate and observationally derived results at lower redshifts as the fraction of passive galaxies is higher at later epochs.

We also estimate ${\rho }_{{{\rm{H}}}_{2}}(z)$ from the SFR density (Wilkins et al. 2019), assuming a constant (Geach & Papadopoulos 2012) and weakly evolving (Tacconi et al. 2018) star formation efficiency. These constant efficiency models predict a co-moving molecular gas mass density in good agreement with both measurements of molecular gas mass via observations of direct CO line emission (Decarli et al. 2016b, 2020; Riechers et al. 2019, 2020) and our results derived from measurements of the long-wavelength dust emission, out to a peak at z ≈ 2. A simple conclusion is that the peak epoch of star formation at z ≈ 2 is not driven by significantly more efficient (or starburst-like) star formation in galaxies, but by a higher abundance of molecular fuel in galaxies. We note that the estimate derived from weakly evolving star formation is ≈1 dex higher than our results at Δz ≃ 0.25. This is likely a consequence of the latter being underestimated due to the abundance of low stellar mass galaxies in this redshift bin.

We recognize that our assumption of a constant efficiency model is at odds with Scoville et al. (2017), who argued that while cold molecular gas reservoirs increase with z (as (1 + z)^1.84), the SFR increases more rapidly (as (1 + z)^2.9), indicating that the peak of star formation is a consequence of both increased molecular gas content in galaxies and higher star formation efficiency. We also note that at z ≳ 1 early-type galaxies have been shown to be more compact for a given stellar mass than their local counterparts (e.g., Daddi et al. 2005; Cappellari et al. 2009), which taken in combination with the “Kennicutt–Schmidt” relation (a power-law relation between SFR and gas surface densities, Schmidt 1959; Kennicutt 1998) implies that star formation may be more efficient at z ≳ 1.

The 3D stacking approach we use derives the average properties for galaxies in our sample as a function of redshift, and thus we do not measure the molecular gas mass and SFRs for individual sources. While the UDS DR11 catalog does include M_stellar estimates (which are evaluated at the peak maximum-likelihood redshift) for individual galaxies, our 3D stacking method bins galaxies according to the discretized redshift probability distribution (${ \mathcal P }(z)$), and as such each galaxy in our sample effectively contributes to the flux in all redshift intervals. Hence, using this 3D stacking method precludes a M_stellar selection relative to our redshift bins. Therefore, we are not able to repeat the analysis from Scoville et al. (2017) to test their assertion of an evolving star formation efficiency.

Although we cannot quantify the contribution of higher star formation efficiencies to the peak of SFR density at z ≈ 2, the symmetry between our constant efficiency models with our statistically derived ${\rho }_{{{\rm{H}}}_{2}}(z)$ indicates a star formation history that is predominantly driven by an increased supply of molecular gas in galaxies, rather than a significant evolution in star formation efficiency (consistent with the findings of Decarli et al. 2020; Magnelli et al. 2020). With this in mind we now turn to the formation of H₂ itself. Cazaux & Spaans (2004) combined a microscopic model for the relative rates of gas-phase and dust H₂ production with a cosmological model to show the more efficient dust-phase production becomes the dominant route to H₂ formation at z ≈ 3–6 for reasonable assumptions about the conditions of the ISM of early galaxies. Therefore, there is a perfect storm for massive galaxy growth at z ≈ 2: not only is the cosmic accretion rate at its peak, massive halos have had time to grow, galaxies have increased gas densities, and previous generations of stars in the progenitors of these systems have provided the metal enrichment that accelerates the formation of H₂, which as the fuel for star formation, drives galaxy growth; this could be described as the epoch of molecular gas.

S16 intentionally restricted their calibration sample to galaxies at z ≤ 3 to ensure observed 850 μm emission was from the rest-frame long-wavelength RJ tail, where dust is optically thin and emission is dominated by the contribution of cold dust (which is well represented by a mass-weighted T_dust = 25 K). In Figure 3 we have shown an average of our ${\rho }_{{{\rm{H}}}_{2}}(z)$ estimates at 2.5 ≲ z ≲ 6 (the UDS redshifts are untested at earlier epochs as there are no UDS galaxies with spectroscopic redshifts at z ≳ 6.5), but note that at these redshifts estimates are less reliable due to large uncertainties in the RJ correction (see Equation (10)) as rest-frame emission approaches the peak of the SED.

In the optically thick regime (as rest-frame dust emission moves off the long-wavelength RJ tail) the rest-frame emission no longer correlates with the total dust mass of a galaxy and probes only the surface dust, which using the approach of S16 would result in an underestimation of L₈₅₀ and hence the molecular gas mass. However, as the rest-frame emission approaches the peak of the SED we are increasingly sensitive to the dense, warm dust component, which significantly boosts the luminosity (with only a small mass fraction) and dominates the emission close to the SED peak. Consequently, rest-frame dust emission at high z is not well represented by a mass-weighted T_dust = 25 K, which would result in an overestimate of the dust and gas mass.

In addition to these competing effects, we are also likely to be missing a significant population of lower mass galaxies at z > 2.5. As shown in Figure 4, the 95% stellar mass completeness at z ≈ 2.5 is predicted to be ≃10^9.5 M_⊙, so we are are simply not sensitive to the majority of low mass galaxies at the highest redshifts. In addition, although relatively rare (with number counts N( > 3.5 mJy) ≃ 3000 deg⁻²; Geach et al. 2017) SMGs are dust rich (M_dust ∼ 10⁹ M_⊙; e.g., da Cunha et al. 2015; Magnelli et al. 2019) and about 20% are undetected in optical/NIR surveys (e.g., Dudzevičiūtė et al. 2020). This non-detection of SMGs is unlikely to have a significant impact on our estimates at low z. However, at z ≳ 2.5 since we are significantly under-sampling the galaxy population and as the number of galaxies in our redshift bins fall the non-detection of dust-rich SMGs becomes more statistically significant, further contributing to an underestimation of the molecular gas mass density at z ≳ 2.5.

The overall impact of the above is difficult to quantify. However, as shown by Figure 3 our results at z ≳ 2.5 are systematically lower than the estimates obtained via direct CO line emission, which suggests that the use of this method past z ≈ 2.5 (when λ_obs = 850 μm no longer probes the rest-frame RJ tail) results in an underestimation of the molecular gas mass density. In consequence, while our results are highly uncertain at z > 2.5, we suggest that to z ≲ 6 these can be seen as providing a lower limit to the molecular gas mass density.