Automated Mining of the ALMA Archive in the COSMOS Field (A³COSMOS). II. Cold Molecular Gas Evolution out to Redshift 6

In Paper I we presented the A³COSMOS project, which is an Automated ALMA Archive mining in the COSMOS field. We developed pipelines for producing continuum images using nearly all publicly available ALMA archival data in COSMOS (regardless of observing bands, but we discarded very high resolution (beam size $\lt 0\buildrel{\prime\prime}\over{.} 1$) data; see Paper I). We performed two major (sub)millimeter continuum photometry extractions, one prior-based and one blind extraction, to make sure the photometries are robust and outliers are identified (see below). Both photometries are verified by extensive Monte Carlo simulations and corrected for flux bias and uncertainty. Additional photometry tasks using apertures following S17 show good consistency for isolated sources ($\lt 20 \%$ difference on average, but they significantly differ for blended or merger-like sources for which aperture photometry is not suitable).

In order to obtain the most robust galaxy catalog from the initial (sub)millimeter continuum detections, we applied very strict criteria to select ALMA detections: a peak flux to rms noise ratio of 5.40 for blind extraction and 4.35 for prior photometry, which correspond to an expected spurious source fraction of $\sim 10 \%$ (according to our statistical analysis). These spurious sources are statistically unavoidable in the initial photometry catalogs, but we developed a series of assessments to identify the most reliable detections. We hence removed ALMA detections that (1) have inconsistent fluxes between blind- and prior-based (sub)millimeter photometry (identified by the Flag_inconsistent_flux in the A³COSMOS catalog, which are about 10 sources that likely are mergers or blended sources and exhibit a ≳0.5 dex difference between blind- and prior-photometry fluxes; see examples in Appendix B of Paper I); (2) have a peculiar counterpart association quality (Flag_outlier_CPA; which is likely due to chance alignment between a prior source and a noise peak); and/or (3) show an excess in ALMA flux relative to the galaxy SED (Flag_outlier_SED; which is likely due to inconsistent photometric redshift, blended sources, or noise). These criteria exclude sources that are either boosted by noise in the ALMA image or multiple galaxies coaligned, plus other less-clear situations. For more details, we refer the reader to Paper I.

After removing the spurious sources, our robust galaxy catalog from A³COSMOS (version 20180801) contains 669 galaxies (36% have spectroscopic redshifts mainly from the COSMOS spec-z catalog compiled by M. Salvato; see references in Paper I). Due to the strict additional selection criteria, the spurious fraction is reduced to close to zero according to our statistics in Paper I. Yet this implies that we miss a significant number of low ALMA signal-to-noise ratio (${\rm{S}}/{\rm{N}}$) sources that have a $\lt 50 \%$ chance of being real, faint galaxies. For comparison, S17 explored all ALMA Band 6 and 7 data in the ALMA public archive and selected sources with total flux of ${\rm{S}}/{\rm{N}}\gt 2$. Betti et al. (2019) analyzed ALMA continuum data for 101 galaxies and selected 68 as detections with an aperture-based total flux of ${\rm{S}}/{\rm{N}}\gt 2$ or peak flux of ${\rm{S}}/{\rm{N}}\gt 3$. The data used in Betti et al. (2019) were public in the ALMA archive before 2018 August and are therefore in our catalog. Ninety of their galaxies appear in our prior-fitting catalog without applying an ${\rm{S}}/{\rm{N}}$ selection; however, only eight galaxies have a peak flux ${\rm{S}}/{\rm{N}}\gt 4.35$, which is our selection criterion based on statistics (corresponding to a spurious fraction $\sim 10 \%$). This quick comparison demonstrates that our catalog has very strict constraints and only considers the statistically most robust ALMA detections. Lowering the selection criterion for A³COSMOS from a (peak flux) ${\rm{S}}/{\rm{N}}$ of 4.35 to 3.0 doubles the A³COSMOS galaxy sample, but 40% of the sample will be spurious based on our simulation statistics. Given this trade-off between increased sample size and decreased reliability, we resort to the original robust galaxy catalog containing only highly reliable sources from A³COSMOS.

Galaxy properties in the A³COSMOS galaxy catalog, including stellar mass (${M}_{\star }$), IR luminosity (${L}_{\mathrm{IR}}$), and dust mass (${M}_{\mathrm{dust}}$), are obtained from MAGPHYS (da Cunha et al. 2008, 2015) SED fitting to their optical-to-radio SEDs (see Paper I). We compute the dust-obscured SFR from IR luminosity following the Kennicutt (1998a) calibration and Chabrier (2003) IMF: $\mathrm{SFR}={L}_{\mathrm{IR}}/1\times {10}^{10}\ {M}_{\odot }\,{\mathrm{yr}}^{-1}$.

In addition, to understand whether using MAGPHYS SED fitting is biased by the built-in SED templates or the assumption of energy balance, we performed two more independent SED fittings for each galaxy to fit the stellar (up to IRAC ch2) and near-IR-to-radio data points separately, with the FAST (Kriek et al. 2009) and "super-deblended" (described in Liu et al. 2018; and also used by Jin et al. 2018) SED fitting tools, respectively. We find that the MAGPHYS-fitted stellar masses are systematically larger by about 0.25 dex than the FAST-fitted values (with a scatter of 0.30 dex), while the dust-obscured SFRs are fully consistent between MAGPHYS and the super-deblended SED fitting. The systematic discrepancy in stellar mass has also been found by Battisti et al. (2019) and reproduced in SED modeling with various nonparametric star formation histories (e.g., Leja et al. 2019). Since this is not yet fully understood, we still adopt the MAGPHYS SED fitting results. We tested that using FAST-fitted stellar masses will not change our main results, but only alter the coefficients in our equations (by $\lesssim 20 \%$).

We include 20 samples of galaxies with CO observations and well-constrained stellar mass and SFR properties from the literature as information complementary to our analysis. The full list is presented in Table 1 (starting from the second row). It encompasses most of the CO-observed samples analyzed by T18. Most of these samples are galaxies in the local universe, and the largest sample is the xCOLD GASS survey sample (Saintonge et al. 2017).

Table 1.  Galaxy Samples Used for This Study

Sample Name	z	${\mathrm{log}}_{10}({M}_{\star }/{M}_{\odot })$	${N}_{\det .}$ a	References
A3COSMOS (Paper I)	0.29–5.667b	∼10.0–12.0	658b	Liu et al. (2019, data set 20180801)
DGS	0–0.045	6.5–10.6	32	Rémy-Ruyer et al. (2014, 2015)
HRS	0.0034–0.006	8.7–11.7	99	Andreani et al. (2018, and refs. thereinc)
KINGFISH	0.0005–0.006	6.3–10.7	28	Groves et al. (2015, and refs. thereind)
Saintonge+2017 (xCOLD GASS)	0.01–0.05	9.0–11.4	330	Saintonge et al. (2017)
Cicone+2017 (ALLSMOG)	0.01–0.03	8.5–11.5	48	Cicone et al. (2017)
Lisenfeld+2017	0.01–0.07	9.0–11.3	41	Lisenfeld et al. (2017)
Cortzen+2019	0.03–0.29	8.9–11.5	51	Cortzen et al. (2019)
Villanueva+2017 (VALES)	0.03–0.33	10.1–11.3	49	Villanueva et al. (2017)
Bertemes+2018 (Stripe82)	0.03–0.20	10.0–11.3	78	Bertemes et al. (2018)
Kirkpatrick+2014 (5MUSE)	0.05–0.29	10.5–11.4	17	Kirkpatrick et al. (2014)
Bauermeister+2013 (EGNOG)	0.06–0.31	10.7–11.5	14	Bauermeister et al. (2013)
Lee+2017	0.27–0.62	10.0–11.1	20	Lee et al. (2017)
Spilker+2018	0.60–0.75	∼11.0	4	Spilker et al. (2018)
Combes+2013 (ULIRGs)	0.61–0.97	9.3–12.0	12	Combes et al. (2013)
Magdis+2012a (BzK)	0.51–1.60	10.5–11.0	9	Magdis et al. (2012a)
Tacconi+2018 (PHIBSS 1&2)	0.50–2.49	9.8–11.6	148	Tacconi et al. (2018)
Kaasinen+2019	1.78–2.93	10.6–11.7	10	Kaasinen et al. (2019)
Magdis+2017 (LBG)	2.8–2.9	11.28–11.38	1	Magdis et al. (2017)
Magdis+2012b (LBG)	2.9–3.2	11.0–11.3	1	Magdis et al. (2012b)
Tan+2014 (GN20)	4.05–4.06	10.6–11.0	3	Tan et al. (2014)

Notes.

^aWe only include sources with $\gt 3\sigma$ detections. ^bThe largest photometric redshift in the A³COSMOS catalog is 5.54 based on the prior redshift information from Laigle et al. (2016) or Davidzon et al. (2017), and it is 7.2 based on Jin et al. (2018). The largest spectroscopic redshift is 5.667 based on the prior information from Capak et al. (2015). Eleven sources have IR/mm photo-z = 5.7–7.2 only from Jin et al. (2018) and are very uncertain. However, our test in Section 5 shows that including or excluding them does not obviously alter our results. ^cWe used Table 1 of Andreani et al. (2018) for the SFR (IR luminosity), stellar mass, and molecular gas properties, and Table 3 of Hughes et al. (2013) for metallicity. The HRS survey is described in Boselli et al. (2010), with Herschel photometry originally from Ciesla et al. (2012, 2014) and gas properties originally from Boselli et al. (2014a, 2014b). ^dWe used Tables 1 and 2 of Groves et al. (2015) for the SFR, stellar mass, molecular gas, and metallicity properties. The KINGFISH survey is described in Kennicutt et al. (2011), with dust photometry originally from Dale et al. (2012) and molecular gas originally from Leroy et al. (2009).

Download table as:  ASCIITypeset image

Saintonge et al. (2017) applied a metallicity-dependent ${\alpha }_{\mathrm{CO}}$ according to Accurso et al. (2017) to convert their CO observations into molecular gas mass. A similar metallicity-dependent ${\alpha }_{\mathrm{CO}}$ is also adopted by the Bertemes et al. (2018) and T18 PHIBSS samples (with slightly different equations; see Appendix A.1). Most other complementary samples either assume only a single ${\alpha }_{\mathrm{CO}}$ value, that is, either a Galactic value or an ultraluminous infrared galaxy (ULIRG) value (see Appendix A.1), or bimodal values depending on the galaxy type (e.g., Villanueva et al. 2017).

To homogenize the complementary sample, we recalculated all molecular gas masses from the CO line luminosities by applying the metallicity-dependent ${\alpha }_{\mathrm{CO}}$ following T18. We use metallicity to calculate ${\alpha }_{\mathrm{CO}}$ when available (mostly for $z\lesssim 0.3$ galaxies; where metallicity is from optical emission lines using the Pettini & Pagel 2004 calibration or converted to that calibration following Kewley & Ellison 2008 where necessary). Otherwise we first estimate the metallicity using the mass–metallicity relation following Genzel et al. (2015, Equation 12(a); see also Appendix A.4), and then we calculate the ${\alpha }_{\mathrm{CO}}$. The recomputed molecular gas masses are within a factor of 2 ($\lesssim 0.36$ dex in logarithm) from their originally obtained values.

In the SED-fitted dust mass method, we first obtain the dust mass (${M}_{\mathrm{dust}}$), dust mean temperature (${T}_{\mathrm{dust}}$), and dust emissivity (${\beta }_{\mathrm{dust}};$ describing the dust opacity κ's wavelength dependency; e.g., Li & Draine 2001) from optical-to-millimeter SED fitting; we then apply a gas-to-dust ratio, ${\delta }_{\mathrm{GDR}}\equiv {M}_{\mathrm{total}\mathrm{gas}}/{M}_{\mathrm{dust}}$, which relates total gas (molecular and atomic) mass to dust mass.

In the first step, different assumptions on dust grain models can lead to variations in the determined dust properties. Yet simulations (e.g., Hayward & Smith 2015) and observations of local galaxies (e.g., Hayward & Smith 2015; Hunt et al. 2019) indicate that SED fitting tools like MAGPHYS are able to reasonably recover galaxies' dust properties (at least for ${L}_{\mathrm{IR}}\gt {10}^{11}\ {L}_{\odot }$). For our work, we ignore the systematic uncertainty introduced by different dust grain models (i.e., different SED fitting tools). This might not be entirely correct, but further investigation of this topic requires a subsample of galaxies with well-sampled SEDs and accurate spectroscopic redshifts, which is beyond the scope of this paper.

In the second step, ${\delta }_{\mathrm{GDR}}$ is found to strongly depend on metallicity (e.g., Leroy et al. 2011; Rémy-Ruyer et al. 2014; De Vis et al. 2019; see more details in Appendix A.2), and the latter is correlated with the stellar mass (known as the mass–metallicity relation; see detailed discussion in Appendix A.4). Differences exist among the empirical scaling relations in the literature, whereas our ALMA continuum observations preferentially select intensely SFGs with ${M}_{\star }\gt 2\times {10}^{10}\ {M}_{\odot }$, which exhibit a close-to-solar metallicity based on the mass–metallicity relation at $z\sim 2.3$ of Erb et al. (2006). Our analysis is, therefore, only affected by the relatively small offset of 0.1–0.2 dex between the relations of Leroy et al. (2011) and Rémy-Ruyer et al. (2014) at $\gt 0.5$ solar metallicity.

As A³COSMOS galaxies do not have homogeneous metallicity measurements, we compute the metallicity based on redshift and stellar mass for each of the galaxies using the mass–metallicity relation of Genzel et al. (2015, Equation 12(a)) and compute the ${\delta }_{\mathrm{GDR}}$ using the Rémy-Ruyer et al. (2014) prescription. Our detailed comparison of various forms of the mass–metallicity relation and the "fundamental metallicity relation" (FMR; e.g., Mannucci et al. 2010, 2011; yet still debated) in Appendix A.4 shows that Equation 12(a) of Genzel et al. (2015, which is also used by T18) provides the most plausible predictions for the metallicity of high-redshift $z\gt 1$ galaxies. Here we adopt a slightly modified form of

$$\begin{eqnarray}&&\begin{array}{l}12+{\mathrm{log}}_{10}({\rm{O}}/{\rm{H}})\ =\\ \left\{\begin{array}{l}a,\qquad \qquad \qquad \ \ {\rm{if}}\ {\mathrm{log}}_{10}({M}_{\star }/{M}_{\odot })\geqslant b(z),\\ a-0.087\times {\left({\mathrm{log}}_{10}\left({M}_{\star }/{M}_{\odot }\right)-b(z)\right)}^{2},\qquad \ {\rm{else.}}\end{array}\right.\\ {\rm{where}}\ a=8.74\ {\rm{and}}\ b(z)=10.4+4.46\times {\mathrm{log}}_{10}(1+z)\\ -\ 1.78\times {\left({\mathrm{log}}_{10}(1+z)\right)}^{2}.\end{array}\end{eqnarray} \tag{ 6 }$$

The modification (under the ${\mathrm{log}}_{10}({M}_{\star }/{M}_{\odot })\geqslant b(z)$ condition) prevents a drop in metallicity for very massive galaxies at $z\lt 1$ (see Figure 18).

Finally, we consider that our A³COSMOS high-redshift galaxies are molecular-rich (same as assumed by T18 at $z\gt 0.4$); that is, the molecular-to-total-gas ratio ${f}_{\mathrm{mol}\mathrm{frac}}$ is unity. In this way, we obtain the dust-SED-based ${M}_{\mathrm{mol}\mathrm{gas}}$ by multiplying ${M}_{\mathrm{dust}}$ with the mass–metallicity-derived ${\delta }_{\mathrm{GDR}}$ and ignore the contribution from atomic gas. Hereafter we refer to this method as the "${\delta }_{\mathrm{GDR},Z}$" method.

We caution that, as discussed in Appendix A.3, observations of local galaxies actually indicate that ${f}_{\mathrm{mol}\mathrm{frac}}$ is usually below 50% even for a galaxy with ${M}_{\star }\sim 1\times {10}^{11}\ {M}_{\odot }$. Applying an actual ${f}_{\mathrm{mol}\mathrm{frac}}$, for example, based on the Krumholz et al. (2009) correlation with stellar mass or metallicity, will lead to a lower ${M}_{\mathrm{mol}\mathrm{gas}}$. Based on our next comparison of ${M}_{\mathrm{mol}\mathrm{gas}}$ calibrations (see Section 3.3), this will cause even a larger difference to the RJ-tail dust continuum methods where atomic gas is also not considered. Therefore, here we choose to not account for the atomic gas and leave the consideration of an actual ${f}_{\mathrm{mol}\mathrm{frac}}$ to future work.

Recent studies show that dust continuum luminosity at rest-frame RJ-tail wavelengths tightly correlates with gas mass or CO line luminosity across two orders of magnitude in local and high-redshift galaxies (e.g., Bourne et al. 2013; Scoville et al. 2014, 2016; Groves et al. 2015; Hughes et al. 2017; Saintonge et al. 2017; Bertemes et al. 2018; Kaasinen et al. 2019). Scoville et al. (2014) found a constant ratio between dust continuum luminosity and ${M}_{\mathrm{total}\mathrm{gas}}$,

$$\begin{eqnarray}&&{\alpha }_{{\rm{rj}},\mathrm{tot}}\equiv \displaystyle \frac{{L}_{{\nu }_{{\rm{rj}}}}/(\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1})}{{M}_{\mathrm{total}\mathrm{gas}}/({M}_{\odot })},\end{eqnarray} \tag{ 7 }$$

where they calibrated ${\alpha }_{{\rm{rj}},\mathrm{tot}}$ to be $1.0\pm 0.23\times {10}^{20}$ $(\mathrm{erg}\,{{\rm{s}}}^{-1}{\mathrm{Hz}}^{-1}\,{M}_{\odot }^{-1})$ at rest-frame 850 μm with a small sample of 12 local galaxies. Meanwhile, Groves et al. (2015) studied the atomic, molecular gas, and dust continuum at rest-frame 70, 100, 160, 250, 350, and 500 μm in 36 local spiral galaxies. They found a mean ${M}_{\mathrm{total}\mathrm{gas}}/\nu {L}_{\nu ,500}=28.5$ for near-solar-metallicity galaxies, corresponding to ${\alpha }_{500,\mathrm{tot}}=2.2\times {10}^{20}$ $(\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}\,{M}_{\odot }^{-1})$, and a factor of 10 lower values for much more metal-poor galaxies. According to Equation (9) of Scoville et al. (2014), ${\alpha }_{{\rm{rj}},\mathrm{tot}}$ is proportional to dust opacity ${\kappa }_{\nu }$, which scales with frequency by ${\kappa }_{\nu }\propto {\nu }^{1.7-2.0}$ (Li & Draine 2001), so it is expected that ${\alpha }_{{\rm{rj}},\mathrm{tot}}$ is a factor of 2.5 higher at 500 μm than at 850 μm. Therefore the calibrations are consistent between Groves et al. (2015) and Scoville et al. (2014). Meanwhile, the variation from metal-rich to metal-poor galaxies can also be explained by a dramatic change in ${\delta }_{\mathrm{GDR}}$ (Appendix A.2).

Focusing on molecular gas only, Scoville et al. (2016) calibrated the ratio between the RJ-tail dust continuum luminosity and ${M}_{\mathrm{mol}\mathrm{gas}}$,

$$\begin{eqnarray}&&{\alpha }_{{\rm{rj}},\mathrm{mol}}\equiv \displaystyle \frac{{L}_{{\nu }_{{\rm{rj}}}}/(\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1})}{{M}_{\mathrm{mol}\mathrm{gas}}/({M}_{\odot })},\end{eqnarray} \tag{ 8 }$$

to be $6.7\pm 1.7\times {10}^{19}$ $(\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}\,{M}_{\odot }^{-1})$ at rest-frame 850 μm for a few tens of local spirals, ULIRGs, and $z\sim 2$ submillimeter galaxies (SMGs). Later studies with larger samples of CO and RJ-tail continuum observations found slightly nonlinear correlations; that is, ${\alpha }_{{\rm{rj}},\mathrm{mol}}$ has a dependency on ${L}_{{\nu }_{{\rm{rj}}}}$ or ${L}_{\mathrm{CO}}^{{\prime} }$ (e.g., Hughes et al. 2017; Bertemes et al. 2018 and Saintonge et al. 2018). As their samples span a wide range of stellar mass from 10⁹ to ${10}^{12}\ {M}_{\odot }$ and CO $J=1\to 0$ line luminosity ${L}_{\mathrm{CO}(1-0)}^{{\prime} }$ from 10⁷ to ${10}^{12}\ {\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{pc}}^{2}$, galaxies have significantly varied metallicity, ${f}_{\mathrm{mol}\mathrm{frac}}$, and ${\delta }_{\mathrm{GDR}}$. A simple explanation for the variations is that ${\alpha }_{{\rm{rj}},\mathrm{mol}}$ scales with ${f}_{\mathrm{mol}\,\mathrm{frac}}^{-1}$ and ${\delta }_{\mathrm{GDR}}^{-1}$, which both relate to metallicity (see Appendix A.3 and A.2, respectively).

Since the literature on the calibration of ${\alpha }_{{\rm{rj}},\mathrm{mol}}$ is already very rich, we do not further discuss it here. In the following, we will adopt the three calibrations from S17, Hughes et al. (2017), and Groves et al. (2015), referred to as the "${\alpha }_{850,{\rm{S}}17}$," "${\alpha }_{850,{\rm{H}}17}$", and "${\alpha }_{160\mbox{--}500,{\rm{G}}15}$" method, respectively. For the "${\alpha }_{160\mbox{--}500,{\rm{G}}15}$" method, we use the calibration factors for the ${\mathrm{log}}_{10}({M}_{\star }/{M}_{\odot })\gt 9$ galaxies in Table 5 of Groves et al. (2015), and we assume that our A³COSMOS galaxies have negligible atomic gas contribution. These works directly calibrate the ratio between ${L}_{{\nu }_{{\rm{rj}}}}$ and ${M}_{\mathrm{mol}\mathrm{gas}}$ (and "${\alpha }_{850,{\rm{H}}17}$" and "${\alpha }_{160\mbox{--}500,{\rm{G}}15}$" include a luminosity dependency), so the need for a calibration of the underlying ${f}_{\mathrm{mol}\mathrm{frac}}$ and ${\delta }_{\mathrm{GDR}}$ is bypassed.

3.2.1. Band Conversion from Observed-frame to Rest-frame RJ Tail

The good agreement between RJ-tail dust continuum-to-gas mass calibrations and the overwhelming observational efficiency compared to (sub)millimeter line observations make the RJ-tail dust method very favorable and promising for large surveys at high redshift.

Our high-redshift galaxies are most commonly observed in ALMA Band 6 and 7, which correspond to rest-frame $\lesssim 250\,\mu {\rm{m}}$ and ≲160 μm, respectively, for galaxies at $z\gtrsim 4$. In Figure 2, we show the longest rest-frame wavelengths of the available ALMA data (${\lambda }_{\mathrm{rest}}$) for each galaxy in our sample. Some 85% of our sources have ${\lambda }_{\mathrm{rest}}\geqslant 250\,\mu {\rm{m}}$, while the rest only probe shorter-wavelength dust continua.

Zoom In Zoom Out Reset image size

In order to apply the ${\alpha }_{{\rm{rj}},\mathrm{mol}}$ conversion from dust continuum to molecular gas mass, a "band conversion"¹⁹ is needed to obtain the corresponding flux density at the calibrated rest-frame wavelength, that is, rest-frame 850 μm for applying the ${\alpha }_{850,{\rm{S}}17}$ and ${\alpha }_{850,{\rm{H}}17}$ methods, and either 160, 250, 350, or 500 μm for applying the ${\alpha }_{160\mbox{--}500,{\rm{G}}15}$ method.

We use our MAGPHYS SED fitting for the band conversion, that is, predicting a longer-wavelength flux density with an SED covering only shorter wavelengths. MAGPHYS fits the dust SED with two dust components, one associated with actual star-forming birth clouds and the other exposed to the ambient interstellar radiation field. The former dust usually has a high temperature and dominates the short-wavelength (e.g., ${\lambda }_{\mathrm{rest}}\lt 60\,\mu {\rm{m}}$) flux density, while the latter dust is constrained to have a temperature in the range of 15–25 K (da Cunha et al. 2008, 2015) and dominates the long-wavelength flux density. A similar idea of composite dust models is also adopted by Draine & Li (2007) and used in fitting local SFGs (e.g., Draine et al. 2007; Aniano et al. 2012) and high-redshift galaxies (e.g., Magdis et al. 2012a, 2014, 2017).

Using such composite-model SED fitting for band conversion has a large advantage over using a single-temperature modified blackbody, as it is much less biased toward the luminosity-weighted dust temperature. Privon et al. (2018) studied the systematic bias of the band conversion using their zoomed-in cosmological simulations, finding that assuming a single-temperature modified blackbody SED for conversion leads to a more than 0.5 dex overestimation in ${L}_{{\nu }_{850\mu {\rm{m}},\mathrm{rest}}}$ when the true dust temperature is a factor of two different than assumed (see their Figure 5).

Whereas MAGPHYS performs well in fitting the dust SED shape, the sampling of the dust SED is usually limited by the available data for $z\gt 4$ sources (as shown in Figure 2). In Appendix B we perform a test to estimate the bias of lacking long-wavelength data in predicting longer-wavelength flux density. We find that when having only ${\lambda }_{\mathrm{rest}}\leqslant 160\,\mu {\rm{m}}$ data points, MAGPHYS underpredicts the rest-frame 850 μm flux density by up to 0.8 dex (on average 0.4 dex) when the dust continua photometries have a quadratic-added mean ${\rm{S}}/{\rm{N}}\lesssim 15$. Meanwhile, the worse case of having only the rest-frame 160 μm data point available over the 8 μm to 3 mm range causes a similar bias by MAGPHYS.

To apply the band conversion, we first compute the ratio between the SED-predicted flux densities at $850\times (1+z)\mu {\rm{m}}$ and the observed wavelength,

$$\begin{eqnarray}&&{{\rm{\Gamma }}}^{\mathrm{SED}}\equiv {S}_{{\nu }_{850\times (1+{\rm{z}})\mu {\rm{m}}}}^{\mathrm{SED}}/{S}_{{\nu }_{\mathrm{obs}}}^{\mathrm{SED}},\end{eqnarray} \tag{ 9 }$$

and then we scale the observed ALMA flux density by ${{\rm{\Gamma }}}^{\mathrm{SED}}$ and compute the luminosity:

$$\begin{eqnarray}&&{L}_{{\nu }_{850\mu {\rm{m}},\mathrm{rest}}}=4\pi {d}_{{\rm{L}}}^{2}\times {S}_{{\nu }_{\mathrm{obs}}}^{\mathrm{ALMA}}\times {{\rm{\Gamma }}}^{\mathrm{SED}}/(1+z).\end{eqnarray} \tag{ 10 }$$

In principle, we can also directly take the SED-predicted rest-frame 850 μm flux density ${S}_{{\nu }_{850\times (1+{\rm{z}})\mu {\rm{m}}}}^{\mathrm{SED}}$, but this would lead to underpredicted scatter in our analysis due to the degeneracy within SED models.

Finally, we divide the luminosity ${L}_{{\nu }_{850\mu {\rm{m}},\mathrm{rest}}}$ by ${\alpha }_{850}$ derived from S17 and Hughes et al. (2017) to obtain the "${\alpha }_{850,{\rm{S}}17}$" and "${\alpha }_{850,{\rm{H}}17}$" molecular gas masses. In the "${\alpha }_{160\mbox{--}500,{\rm{G}}15}$" method, because Groves et al. (2015) provided calibrations at six calibration wavelengths (70, 100, 160, 250, 350, and 500 μm), we perform the band conversion from the longest-wavelength ALMA data to its nearest calibration wavelength.

In Figure 3 we compare the molecular gas masses estimated from the above-mentioned "${\delta }_{\mathrm{GDR},Z}$," "${\alpha }_{850,{\rm{S}}17}$," "${\alpha }_{850,{\rm{H}}17}$", and "${\alpha }_{160\mbox{--}500,{\rm{G}}15}$" methods. As shown in the bottom row of the figure, "${\delta }_{\mathrm{GDR},Z}$" leads to systematically lower gas masses than the other three RJ-tail continuum methods. The bias is stronger for sources that do not have long-wavelength (${\lambda }_{\mathrm{rest}}\gt 250\,\mu {\rm{m}}$) coverage. This is closely related to the MAGPHYS SED fitting feature, where missing long-wavelength data seems to lead to an underestimation of the cold, ambient dust that dominates the total dust mass (consistent with the tests in Appendix B).

Zoom In Zoom Out Reset image size

In the top panel, the "${\alpha }_{850,{\rm{S}}17}$" and "${\alpha }_{850,{\rm{H}}17}$" methods agree within 0.1 dex for sources with long-wavelength coverage (but up to about 0.3 dex for sources lacking $\gt 160\,\mu {\rm{m}}$ data). However, a systematic offset of about 0.1 dex exists, which is likely because "${\alpha }_{850,{\rm{S}}17}$" uses a single conversion factor while "${\alpha }_{850,{\rm{H}}17}$" uses a luminosity-dependent conversion factor. The latter has been confirmed by many other works (e.g., Bertemes et al. 2018; Saintonge et al. 2018) and therefore is more reliable.

For panels in the second row, the gas masses based on the "${\alpha }_{850,{\rm{S}}17}$" and "${\alpha }_{850,{\rm{H}}17}$" methods are compared to those using the "${\alpha }_{160\mbox{--}500,{\rm{G}}15}$" method. The "${\alpha }_{160\mbox{--}500,{\rm{G}}15}$" method leads to 0.25 dex lower molecular gas masses than "${\alpha }_{850,{\rm{S}}17}$," or 0.15 dex lower than "${\alpha }_{850,{\rm{H}}17}$" for the majority of sources. A small number of sources with poor long-wavelength coverage, however, have smaller differences. This is probably due to the smaller ${M}_{\star }\gt {10}^{9}\ {M}_{\odot }$ sample in Groves et al. (2015) and the intrinsic variation in ${L}_{{\nu }_{\mathrm{RJ},\mathrm{rest}}}$/${M}_{\mathrm{mol}\mathrm{gas}}$.

To summarize, we find that the gas mass calibrations are as follows: ${{M}_{\mathrm{mol}\mathrm{gas}}}_{({\delta }_{\mathrm{GDR},Z})}$ ≲ ${{M}_{\mathrm{mol}\mathrm{gas}}}_{({\alpha }_{160\mbox{--}500,{\rm{G}}15})}\lt {{M}_{\mathrm{mol}\mathrm{gas}}}_{({\alpha }_{850,{\rm{H}}17})}$ $\lt {{M}_{\mathrm{mol}\mathrm{gas}}}_{({\alpha }_{850,{\rm{S}}17})}$. The systematic offsets are about 0.15–0.25 dex, but are comparable to the scatter of the data. Considering the relatively better agreement of the "${\alpha }_{850,{\rm{H}}17}$" method to other methods as well as recent observations (Bertemes et al. 2018; Saintonge et al. 2018), we choose the "${\alpha }_{850,{\rm{H}}17}$" method as our final gas mass calculation for the A³COSMOS galaxies. We also tested our full analysis with other gas mass calibrations in Section 5, finding that our results are not obviously altered.

In Figure 4 we show the specific star formation rate ($\mathrm{sSFR}$) versus redshift distribution, where sSFR is normalized by the MS sSFR of Speagle et al. (2014) at each redshift. The left and right panels have the same data points but have different X-axis scales and color schemes: the X-axis (redshift) is logarithmic in the left panel and only the complementary samples are color coded, while the redshift is linear in the right panel and only A³COSMOS data points are color coded. Whereas our A³COSMOS sample primely populates the $z\gt 1$ regime, the complementary samples provide coverage at $z\lt 1$. However, we do notice that the MS is not well sampled at $0.1\lesssim z\lesssim 0.5$ and $z\gt 1$. Only the most massive A³COSMOS galaxies (${\mathrm{log}}_{10}({M}_{\star }/{M}_{\odot })\gtrsim 11.5$) sample well the MS, while less-massive ones lie above.

Zoom In Zoom Out Reset image size

The majority of $1\lesssim z\lesssim 2$ complementary sample sources are from the PHIBSS 1&2 surveys (T18) and Kaasinen et al. (2019). Compared to the A³COSMOS galaxies, they are slightly less massive (see Table 1), so T18 sources are able to represent the ${\mathrm{log}}_{10}({M}_{\star }/{M}_{\odot })\sim 10\mbox{--}11$ MS, while the A³COSMOS galaxies are probing the ${\mathrm{log}}_{10}({M}_{\star }/{M}_{\odot })\sim 11\mbox{--}12$ MS at $z\gt 1$.

We also notice that only very low redshift ($z\lesssim 0.03$) galaxies cover the ${\mathrm{log}}_{10}({M}_{\star }/{M}_{\odot })\sim 9\mbox{--}10$ parameter space. In this low-stellar-mass range, the metallicity-dependent ${\alpha }_{\mathrm{CO}}$ might be more uncertain, and so are the estimated molecular gas masses. However, this regime is important in understanding molecular gas scaling relations as shown in later sections, so here we still fit these galaxies from the complementary samples.

Here we study the scaling relations for the two most important molecular gas properties: molecular gas depletion time, ${\tau }_{\mathrm{depl}}\equiv {M}_{\mathrm{mol}\mathrm{gas}}/\mathrm{SFR}$, and molecular gas to stellar mass ratio, ${\mu }_{\mathrm{molgas}}\equiv {M}_{\mathrm{mol}\mathrm{gas}}/{M}_{\star }$. In Figures 5 and 6, we show their distributions versus other galaxy properties: redshift, cosmic age, offset to the MS (${\rm{\Delta }}\mathrm{MS};$ using the MS of Speagle et al. 2014), ${M}_{\star }$, SFR, and sSFR. Two diagrams are shown for each distribution: a scatter plot (upper panels) and a contour plot (lower panels). In the contour plot, we show three sets of contours representing the data densities of A³COSMOS galaxies (orange), PHIBSS 1&2 $0.5\lesssim z\lesssim 2$ galaxies (green), and all other local/low-redshift galaxies (gray), respectively.

The molecular gas depletion time ${\tau }_{\mathrm{depl}}$ spans about one order of magnitude in the high-redshift range from $z\sim 1$ to 6, but has more than two orders of magnitude variation at $z\sim 0$. The latter can be due to the strong correlations with either ${\rm{\Delta }}\mathrm{MS}$, SFR, or sSFR, as shown in the corresponding scatter plots in Figure 5.

Zoom In Zoom Out Reset image size

However, because the SFR and sSFR have redshift dependency and the ${M}_{\star }$ distribution is biased differently from low to high redshift, it is unclear from just this figure which galaxy property mostly determines ${\tau }_{\mathrm{depl}}$. There is even a break or turnover feature in the cosmic age and SFR versus ${\tau }_{\mathrm{depl}}$ panels, which is likely caused by the selection bias at $z\gt 3$, where we only cover the most massive galaxies (${\mathrm{log}}_{10}{M}_{\star }\sim 11\mbox{--}12$). In the intermediate-redshift range ($1\lesssim z\lesssim 3$), the A³COSMOS and PHIBSS 1&2 surveys' galaxies have very similar distributions, as can be seen in the contour plots.

Similar plots are shown in Figure 6 for the molecular gas to stellar mass ratio, ${\mu }_{\mathrm{molgas}}$. Compared to gas depletion time distributions, ${\mu }_{\mathrm{molgas}}$ has a nearly three orders of magnitude variation from local to high redshift, and even at $z\gt 1$ the variation is still as large as two orders of magnitude. From the first two panels, we see a moderate redshift evolution and a strong dependency on stellar mass, respectively. Note that ${\mu }_{\mathrm{molgas}}$ also exhibits a strong dependency on ${\rm{\Delta }}\mathrm{MS}$, but local galaxies are systematically offset from the high-redshift ones by nearly one dex below.

Zoom In Zoom Out Reset image size

As shown in the ${M}_{\star }$–${\mu }_{\mathrm{molgas}}$ panel, local galaxies and high-redshift galaxies seem to follow different distributions: local galaxies have similar ${\mu }_{\mathrm{molgas}}$ across different stellar masses, while high-redshift ones exhibit a steep slope. However, we caution that this is likely an artifact of the high-redshift sample selection using submillimeter data, as the submillimeter selection is similar to an SFR selection or a dust-mass selection, picking up massive MS galaxies and less massive but starbursty galaxies (see Figure 4).

The last two columns of Figure 6 show clear and tight correlations between ${\mu }_{\mathrm{molgas}}$ and SFR and sSFR. Local, PHIBSS 1&2 intermediate-redshift and A³COSMOS higher-redshift galaxies form a contiguous distribution from SFR ∼0.1 to >1000 $\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ and sSFR ∼0.01 to $\sim 50\ {\mathrm{Gyr}}^{-1}$. This is likely a combined effect of the SFR or sSFR evolution and the evolution of molecular gas content.

The last scatter plot, when canceling out the ${M}_{\star }$ term in both axes, is equivalent to the ${M}_{\mathrm{mol}\mathrm{gas}}$ versus SFR correlation, that is, the Kennicutt–Schmidt law (Schmidt 1959; Kennicutt 1998a). A log–log space linear fitting gives a slope of ∼0.7 with a scatter of 0.3 dex, consistent with Kennicutt (1998a), as well as the slope of ∼0.8 as measured by Sargent et al. (2014) for MS and strong starbursts separately (see also Section 4.4).

Again, while Figures 5 and 6 show that observational data from a variety of samples span a wide range of parameter space and are consistent where they overlap, we caution that not all parameters are independent in these plots, and the apparent correlations have degeneracies. Therefore, a high-dimensional-space fitting to the data is important to characterize the relative contribution of each key parameter to the observed gas fraction and depletion time. In Figure 7 (online only), we show our data points in three-dimensional (3D) space to better illustrate the complexity, and we perform such a high-dimensional-space function fitting in Section 5.

We show in Figure 8 the composite view, named as the "spindle" diagram, of the distributions of the four parameters: gas fraction, redshift, stellar mass, and SFR. The first-level information in the figure is the sSFR evolution of our galaxies binned by redshift and stellar mass (curves are the Speagle et al. 2014 MS). The second-level information is that in each redshift and stellar mass bin (the boxes in the figure), the horizontal span represents the gas fraction ${f}_{\mathrm{mol}\mathrm{gas}}$ (Equation (1); 0%–100% from bin center to edges; shown symmetric for illustration purposes), and the Y position is still sSFR as indicated by the global Y-axis. We can see that in the high-redshift bins, a higher sSFR corresponds to a larger horizontal span inside each box, which means a higher gas fraction.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (5 panels)

Tarball of figure set images

Since all data in these boxes share the same Y-axis, the sSFR of the data can be directly read from the figure and compared to the MS curves. The inhomogeneity of our sample is obvious in the less-massive galaxy bins (i.e., blue and yellow boxes) extending one to two dex above their corresponding MS, while the most massive galaxies (i.e., the red boxes) merely extend more than one dex above the MS.

To summarize, with this "spindle" diagram, we can more clearly see the following:

• (a)  
  At a fixed redshift, more-massive galaxies have both lower sSFR and gas fraction than less-massive ones.
• (b)  
  At a fixed redshift, galaxies that lie farther above the MS exhibit higher gas fractions (${f}_{\mathrm{mol}\mathrm{gas}}$ approaching 100% for the ones with lowest mass and highest sSFR). (Yet such a trend is debated for individual or small (∼10) samples of strong starburst galaxies (${\rm{\Delta }}\mathrm{MS}\gt 0.6$), e.g., Silverman et al. 2015, 2018.)
• (c)  
  Similar to the sSFR evolution, gas fraction evolves with redshift: galaxies of similar stellar mass and distance to the MS but at an earlier cosmic time tend to have a higher gas fraction. (These trends have been known anecdotally for more than 10 years, e.g., Daddi et al. 2008, 2010a; Tacconi et al. 2008, 2010, 2013; Genzel et al. 2010; Magdis et al. 2012b, to name a few, but they have only been quantified recently with sufficiently large samples as presented here.)

The star formation law (or Kennicutt–Schmidt law; Schmidt 1959; Kennicutt 1998a; hereafter SF law) describes the correlation between molecular gas mass and star formation rate and has an empirical form of $\mathrm{SFR}=A\times {M}_{\mathrm{mol}\,\mathrm{gas}}^{N}$, with a slope $N\approx 1.4$ in the log–log space (Kennicutt 1998a; Gao & Solomon 2004). It is physically motivated by the fact that star formation is fueled by molecular gas. However, galaxies show a large scatter in the ${M}_{\mathrm{mol}\mathrm{gas}}$–$\mathrm{SFR}$ plane, and some galaxies like local ULIRGs (e.g., Sanders et al. 2003) and bright, high-redshift SMGs (e.g., Smail et al. 1997; Blain et al. 2002) are more than a one dex offset from "normal" SFGs. This is also referred to as the bimodal SF law (e.g., see Daddi et al. 2010b; Genzel et al. 2010 for a strictly bimodal scenario; and Sargent et al. 2014 for a continuous dichotomy between "normal" star-forming and starburst galaxies). However, why these galaxies are offset from the normal star-forming SF law and whether they are also starbursts in the MS relations are still poorly explored. Given the popular assumption (or intense debate) on the MS/starburst dichotomy and bimodality of SF laws, for example in analytic galaxy modeling (Béthermin et al. 2012, 2017; Sargent et al. 2012, 2014) and observational studies (Daddi et al. 2010b; Genzel et al. 2010; Silverman et al. 2015, 2018; Elbaz et al. 2018; Cibinel et al. 2019), we investigate these two topics with our large A³COSMOS and compiled sample and present how current models are fitting the data.

Figure 9 shows the correlation between SFR and ${M}_{\mathrm{mol}\mathrm{gas}}$ for all galaxies in this work. Data points are color coded by ${\rm{\Delta }}\mathrm{MS}$, and the A³COSMOS and complementary samples are distinguished by different symbols (circle and cross, respectively). For A³COSMOS galaxies with large SFR ($\sim 100\mbox{--}3000\ {M}_{\odot }\,{\mathrm{yr}}^{-1}$) and ${M}_{\mathrm{mol}\mathrm{gas}}$ ($\sim 5\times {10}^{10}\mbox{--}5\times {10}^{11}\ {M}_{\odot }$), a higher ${\rm{\Delta }}\mathrm{MS}$ means more deviation from the normal "star-forming SF law" (the blue line in Figure 9; adopted from Sargent et al. 2014). The strongest starbursts with more than one dex offset from the MS show a 0.43 dex (median) offset from the "star-forming SF law," while  MS galaxies (${\rm{\Delta }}\mathrm{MS}\lt 0.5$) exhibit a small deviation of 0.12 dex (median). Considering that the A³COSMOS sample does not sample well the below-MS region, the 0.12 dex offset does not prevent us from drawing the conclusion that MS galaxies also follow the normal star-forming SF law.

Zoom In Zoom Out Reset image size

However, the starbursts lying significantly above the MS (${\rm{\Delta }}\mathrm{MS}\sim 1$) seem to behave differently between the high-redshift and low-redshift/local samples. High-redshift starbursts do not show large enough offsets to reach the "starburst SF law," as indicated by the red line in Figure 9 (also from Sargent et al. 2014), which is offset by about one dex from the "star-forming SF law." This is also recently found by CO observations of a small sample of 12 strong starbursts at $z\sim 1.5$ by Silverman et al. (2015, 2018). Meanwhile, some low-redshift/local starbursts with ${\rm{\Delta }}\mathrm{MS}\sim 0.5$ are able to reach the "starburst SF law," and the positive correlation between ${\rm{\Delta }}\mathrm{MS}$ and the offset to the "star-forming SF law" is more clear there.

In Figure 10, we more clearly illustrate the correlation between galaxies' offsets to the MS and the SF law. The X-axis, $\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}}$, represents the offset to the MS, with sSFR_MS computed following Speagle et al. (2014). The Y-axis, ${\tau }_{\mathrm{mol}\mathrm{gas}}/{\tau }_{\mathrm{mol}\mathrm{gas},\mathrm{MS}\mathrm{SF}\mathrm{law}}$, represents the offset to the "star-forming SF law," where ${\tau }_{\mathrm{mol}\mathrm{gas},\mathrm{MS}\mathrm{SF}\mathrm{law}}$ ≡ ${M}_{\mathrm{mol}\mathrm{gas},\mathrm{SF}\mathrm{law}}/{\mathrm{SFR}}_{\mathrm{MS}}$ = $\alpha \times {\mathrm{SFR}}^{\beta }/{\mathrm{SFR}}_{\mathrm{MS}}$, and the α and β coefficients are taken from Sargent et al. (2014). Galaxies are binned into four panels by their stellar masses in Figure 10. Data points are color coded by SFR. Model-predicted curves from Sargent et al. (2014) are shown for comparison. Their model, named the two-star formation-mode (2-SFM) model, assumes that galaxies have two modes of star formation: an MS mode and a starburst mode. MS galaxies (e.g., $\mathrm{sSFR}\lesssim 3\times {\mathrm{sSFR}}_{\mathrm{MS}}$) obey the SF law with a Galactic-like ${\alpha }_{\mathrm{CO}}$, while starbursts with sSFR above the MS (e.g., $\mathrm{sSFR}\gtrsim 3\times {\mathrm{sSFR}}_{\mathrm{MS}}$) are shifted toward the "starburst SF law," and they also have a much lower ${\alpha }_{\mathrm{CO}}$. The shift in the SF-law plane happens most rapidly when the sSFR increases from ∼3 to $\sim 4\times$ the MS's sSFR (see Figure 9 of Sargent et al. 2014), thus causing the steep model turnover seen in Figure 10.

Zoom In Zoom Out Reset image size

The data are more complicated than what the 2-SFM model predicts. Galaxies in the lowest-mass bin (${\mathrm{log}}_{10}{M}_{\star }\lt 9.8$) are below the model-predicted curve, while in the mid-stellar-mass bins (${\mathrm{log}}_{10}{M}_{\star }\sim 9.8\mbox{--}11.2$), some galaxies are above it. The turnover is likely seen in the two higher-mass bins (${\mathrm{log}}_{10}{M}_{\star }\,\gt 10.5$), whereas it is less obvious in the two lower-mass bins. The difference cannot be explained by the calibration of the MS because of the reasonably good agreement between MS calibrations (see Figure 10 caption). The molecular gas masses for the lowest-mass galaxies, which are mostly from complementary samples, are calculated via a metallicity-dependent ${\alpha }_{\mathrm{CO}}$ (see Section 2.2); therefore, ${\alpha }_{\mathrm{CO}}$ seems to be not strongly biased. While their SFRs are derived using optical photometry and lack of far-IR data, they intrinsically have a low metallicity and are dust poor, so the lack of far-IR/millimeter should not introduce a significant bias. Unfortunately, observational evidence is still scarce. Coogan et al. (2019) presented CO nondetections for five low-mass ($\langle {\mathrm{log}}_{10}{M}_{\star }\rangle =9.8$) galaxies at $z\sim 2$, resulting in an upper limit on their gas depletion times of $\lt 0.8\,$ Gyr, or ${\tau }_{\mathrm{mol}\mathrm{gas}}/{\tau }_{\mathrm{MS},\mathrm{SF}\mathrm{law}}\lt 0.6$ (assuming a Galactic ${\alpha }_{\mathrm{CO}}$), in agreement with our findings. If the difference between data and model is truly significant, then it implies that low-mass (${\mathrm{log}}_{10}{M}_{\star }\lesssim 10.0$) MS galaxies might follow a different SF law with $3\times$ faster molecular gas depletion than higher-mass MS galaxies, but this has yet to be confirmed with more observations.

In the other three higher-mass bins, from MS to starburst regime, we find good agreements between the data and model for galaxies close to and below the MS, meaning again that MS galaxies also obey the star-forming SF law. Nevertheless, a number of strong starburst galaxies show slower gas depletion (longer gas depletion times) than they should have according to the model. The majority of these strongest starburst outliers with long gas depletion times are from the complementary samples (e.g., Combes et al. 2013; Villanueva et al. 2017) with CO observations but without metallicity information. In this case, their ${\alpha }_{\mathrm{CO}}$ values are indirectly inferred from their stellar masses and SFRs (see Appendix A.4). Silverman et al. (2015, 2018) found that a different choice of ${\alpha }_{\mathrm{CO}}$ alters the ${\tau }_{\mathrm{mol}\mathrm{gas}}/{\tau }_{\mathrm{MS},\mathrm{SF}\mathrm{law}}$ ratio from close to one to 0.2 for a starburst galaxy with ${\rm{\Delta }}\mathrm{MS}\sim 1$ dex (see their Figure 8). Therefore, it is still unclear how well the molecular gas masses (or stellar mass) can be constrained in these strongest starbursts (more detailed multiline gas studies are needed, e.g., with RJ-tail dust continuum plus multi-J CO (e.g., Liu et al. 2015) plus other tracers, to settle this issue).

We also caution that our high-redshift, intermediate-mass (${\mathrm{log}}_{10}{M}_{\star }\sim 10\mbox{--}11$) sample has a strong bias toward a higher ${\rm{\Delta }}\mathrm{MS}$, so we sample better the region above the model curve than below it. This sample bias is less significant for the most massive bin (${\mathrm{log}}_{10}{M}_{\star }\gtrsim 11$), where we most clearly see the turnover.

To summarize the link between the MS and the SF law, we find that (a) the massive (${\mathrm{log}}_{10}{M}_{\star }\sim 10\mbox{--}12$) MS galaxies obey the SFGs' SF law; (b) from the MS to ${\rm{\Delta }}\mathrm{MS}\gtrsim 1$, galaxies start to deviate from the SFGs' SF law toward the starbursts' SF law, with a rapid change at ${\rm{\Delta }}\mathrm{MS}\sim 0.4\mbox{--}0.6$ dex roughly in agreement with the 2-SFM prediction; (c) low-mass (${\mathrm{log}}_{10}{M}_{\star }\lesssim 10$) galaxies appear to have systematically shorter gas depletion times, and even the MS ones do not obey the SFGs' SF law.

These details will likely stimulate further refinement of the popular models and observing strategies.

In the previous sections, we illustrated the wide dynamical range of our sample. Such a data set is the largest sample for the study of the gas scaling relation and its evolution to date, and it will grow with future processing of the ALMA archive in the COSMOS deep field under A³COSMOS. The distribution of our sample in the (z, ${M}_{\star }$, ${\rm{\Delta }}\mathrm{MS}$) high-dimensional space is roughly contiguous from $(z,{\mathrm{log}}_{10}{M}_{\star },{\rm{\Delta }}\mathrm{MS})\sim (0.0,9.5,-1)$ to $\sim (5.0,12.0,+1)$. The gas mass calibrations for the CO and dust subsamples are in good agreement where they overlap in the parameter space.

Nevertheless, the comprehensive presentation of our data set in previous figures also reveals that the sample is nonuniformly distributed and only partially covers the full parameter space. Our sample is biased toward submillimeter-detected (i.e., IR-bright), massive high-redshift galaxies, as well as CO-detected (gas-rich), massive local/low-redshift galaxies. The impact of such sample biases on the results is hard to quantify with the current data set. Stacking ALMA data for sufficient numbers of faint galaxies with similar properties can help to cover additional portions of the parameter space and will be presented in future work. Meanwhile, low-J CO and RJ-tail dust observations toward samples covering the less-probed areas of the parameter space hold the key to further improving such studies.

Here we discuss the cosmic evolution of the molecular gas depletion time ${\tau }_{\mathrm{depl}}$ as predicted by our best-fit function (Equation (11)) for galaxies in bins of stellar mass and MS offset. In Figure 12, we bin all our 1663 galaxies into 4 × 3 panels, with the ${\mathrm{log}}_{10}{M}_{\star }$ bin center ranging from 9.0 to 12.0 (bin width 1.0) and the ${\rm{\Delta }}\mathrm{MS}$ bin center ranging from −0.5 to $+0.5$ (bin width 0.5). The predictions of our best-fit function are shown as solid lines, while the predictions from the T18 and S17 functions (with their fitting) are shown as long- and short-dashed lines, respectively, for comparison. Galaxies in each panel are also binned in small redshift interval so as to show the mean and scatter at each redshift.

Zoom In Zoom Out Reset image size

From the figure, we can see that our function behaves differently than the other two functions. The evolution of ${\tau }_{\mathrm{depl}}$ exhibits a much stronger dependency on stellar mass in our function. For very massive (${\mathrm{log}}_{10}{M}_{\star }\sim 12.0$) galaxies, our function predicts a factor of about 20 increase in ${\tau }_{\mathrm{depl}}$ from very early cosmic time to the present, while the T18 and S17 functions predict only a factor of 5–8 increase. Data from this work favors our function in these bins. Meanwhile, for low-mass (${\mathrm{log}}_{10}{M}_{\star }\sim 9.0$) galaxies, our function predicts a reversed evolutionary trend compared to the T18 and S17 functions. That means a galaxy with a stellar mass as low as ${\mathrm{log}}_{10}{M}_{\star }\sim 9.0$ has a longer depletion time at an earlier cosmic time, and its star formation speeds up with cosmic age. Current data in these bins are not sufficient to clearly distinguish which function is better. The few CO observations available for local dwarf galaxies (e.g., Bolatto et al. 2011; Cormier et al. 2014) show that ${\tau }_{\mathrm{depl}}$ ranges from 0.1 to $1.0\ \mathrm{Gyr}$ (with $\mathrm{SFR}\sim 20\mbox{--}30$ to $0.04\ {M}_{\odot }\,{\mathrm{yr}}^{-1}$, i.e., from high to low ${\rm{\Delta }}\mathrm{MS}$, respectively). These observations still agree with the predictions of our function.

We caution that this figure does not track the evolution of individual galaxies, as they grow in stellar mass and may have rapidly changed ${\rm{\Delta }}\mathrm{MS}$ with time. Thus, for example, the flat ${\tau }_{\mathrm{depl}}$ versus redshift trend for less-massive (${\mathrm{log}}_{10}{M}_{\star }\sim 10.0$) galaxies seen in the middle columns of the figure does not imply a constant ${\tau }_{\mathrm{depl}}$ for an individual galaxy across its evolution history: its stellar mass growth will move it into a higher stellar mass ${\tau }_{\mathrm{depl}}$ evolution track. In the ${\mathrm{log}}_{10}{M}_{\star }\sim 10.0$ and ${\rm{\Delta }}\mathrm{MS}\sim 0.5$ bin, our function does not fit well the $z\sim 3\mbox{--}5$ galaxies, while the T18 and S17 functions do. This is mainly driven by the small number of low-mass starburst galaxies in this redshift range. As already discussed in Section 4.4, our sample within this range is sparse and biased, and the statistics are expected to be less significant.

If only looking at the function predictions, our function actually provides a coherent picture of galaxy "downsizing" (e.g., Cowie et al. 1996; Thomas et al. 2005); that is, more-massive galaxies (possibly in more-massive dark matter halos) evolve earlier than less-massive galaxies. Meanwhile, the star formation in the most massive galaxies quickly slows down at redshift 2–3, which probably points to the "mass-quenching" effect (e.g., Peng et al. 2010).

Below we also compare the predictions of our functions with other works in the literature. Our formula predicts that, for local galaxies with stellar mass $3\times {10}^{9}$, $3\times {10}^{10}$, $3\times {10}^{11}$, and $3\,\times {10}^{12}\,{M}_{\odot }$, their ${\tau }_{\mathrm{depl}}=0.7$, 1.3, 2.6, and 5.0 Gyr, respectively. In comparison, Huang & Kauffmann (2014, 2015) studied about 600 local galaxies from the HERACLES (Leroy et al. 2009), ATLAS3D (Cappellari et al. 2011; Alatalo et al. 2013), and COLD GASS (Saintonge et al. 2011a, 2011b) surveys, and they found ${\tau }_{\mathrm{depl}}=-0.36{\mathrm{log}}_{10}\mathrm{sSFR}-0.14{\mathrm{log}}_{10}({{\rm{\Sigma }}}_{\star })+5.87$.²³ This translates into ${\tau }_{\mathrm{depl}}=1.2$, 1.5, 2.0, and 2.9 Gyr for the four aforementioned stellar masses, assuming that the galaxy size follows the Fernández Lorenzo et al. (2013) and Shen et al. (2003) size–mass relation. Thus the predictions agree within 20% for the two intermediate stellar mass ranges, or $\sim 50 \%$ for all ranges. We note that the $3\times {10}^{12}\,{M}_{\odot }$ case is an extrapolation of their function as their data only probe galaxies with ${10}^{10}\lt {M}_{\star }/{M}_{\odot }\lt {10}^{11.5}$.

New observations are needed to clearly distinguish which function is better and confirm whether our function can reproduce "downsizing" and "mass quenching." Such observations should prioritize low-mass galaxies at high redshift (with enough sensitivity and integration time), as well as the highest-mass but below-MS galaxies at the early cosmic time (though such galaxies are still rarely found).

Similar to the previous section, we show in Figure 13 the binned view of the evolution of ${\mu }_{\mathrm{molgas}}$ as predicted by our best-fit function Equation (11). The three evolution functions in Figure 13, that is, from our Equation (11), T18, and S17, consistently show that ${\mu }_{\mathrm{molgas}}$ has a strong dependency on stellar mass. More-massive galaxies have a lower gas fraction at the same redshift. These functions are also very close to each other for ${\mathrm{log}}_{10}{M}_{\star }\gtrsim 11$ galaxies at all redshifts below ∼3. For lower-mass galaxies, our function is located between the T18 and S17 ones. The S17 function does not fit well local galaxies because they do not include local samples in their analysis. But at high redshift, our function in this work predicts similarly high gas fractions as the S17 function, which are a factor of 10 higher than those expected from the T18 function (for ${\mathrm{log}}_{10}{M}_{\star }\sim 9$ galaxies).

Zoom In Zoom Out Reset image size

The dependency of ${\mu }_{\mathrm{molgas}}$ or ${f}_{\mathrm{gas}}$ on stellar mass was also found much earlier for local galaxies (e.g., Young & Scoville 1991; McGaugh & de Blok 1997; Kennicutt 1998b; Schombert et al. 2001). Young & Scoville (1991) reported an increase in gas fraction by two orders of magnitude from early-type to late-type galaxies (along the Hubble sequence) in the local universe. This is equivalent to similar orders of magnitude increase in their IR to H-band luminosity ratio, that is, ∝ sSFR (Kennicutt 1998b). McGaugh & de Blok (1997) and Schombert et al. (2001) also found strong decreases in the gas fraction with brighter B-band magnitude (higher stellar mass) and higher stellar surface density including low-surface-brightness local galaxies. This is in agreement with our function's prediction.

More recently, Jiang et al. (2015) reported a similarly strong decrease in gas fraction versus stellar mass as reported here, down to a stellar mass of ${10}^{8.5\mbox{--}9.0}\ {M}_{\odot }$ (see also Cao et al. 2017; Saintonge et al. 2017) with a nonlinear behavior. In their sample, ${\mu }_{\mathrm{molgas}}$ is about 0.08–0.3 for ${\mathrm{log}}_{10}{M}_{\star }\sim 9\mbox{--}10$ galaxies, then decreases to about 0.02–0.1 for ${\mathrm{log}}_{10}{M}_{\star }\sim 10\mbox{--}11$ galaxies, which is slightly below this work at $z\sim 0$ and is likely caused by their use of a constant ${\alpha }_{\mathrm{CO}}$, while the true ${\alpha }_{\mathrm{CO}}$ might be higher for low-mass, metal-poor galaxies.

In summary, the predictions from this work, S17, and T18 only obviously differ in those regimes where not much data are currently available, that is, at low stellar masses across cosmic time and for all stellar masses at $z\gt 5\mbox{--}6$. This work's predictions agree with other individual observations in the literature, and our evolution function has physical implications of "downsizing" and "mass quenching" in galaxy evolution. In general, our analysis also raises the need for future CO and RJ-dust observations of below-MS or less-massive (${\mathrm{log}}_{10}{M}_{\star }\lesssim 10$) galaxy samples.

In recent years, deep HST, Spitzer, and ground-based near-IR observations in deep fields have pushed the accurate measurements of the SFG SMFs out to $z\sim 4\mbox{--}5$ (e.g., Marchesini et al. 2009; Peng et al. 2010; Baldry et al. 2012; Santini et al. 2012; Ilbert et al. 2013; Moustakas et al. 2013; Muzzin et al. 2013; Grazian et al. 2015; Song et al. 2016; Davidzon et al. 2017; Wright et al. 2017, 2018). Similarly, deep Herschel far-infrared/submillimeter and ground-based submillimeter surveys pushed the accurate measurements of cosmic SFR density (CSFRD) out to $z\sim 3\mbox{--}4$ as well (e.g., Madau & Dickinson 2014; Liu et al. 2018; and references therein). The CSFRD represents the SFR at each cosmic epoch, so by integrating the SFR across all of the previous cosmic time, we will be able to obtain the total stellar mass density at that time. Meanwhile, the integration of the (SFG) SMF at that cosmic time should in principle equal the total stellar mass integrated from the CSFRD. In Appendix D, we verify that they are in good agreement for the redshift bins where empirical SMFs are available.

Note that the CSFRD has been described as a function of redshift (double power law; Madau & Dickinson 2014), while it is still difficult to characterize the SMF as a contiguous function of redshift. Wright et al. (2018) provide such a functional form, but their function exhibits certain deviations from direct measurements (see Appendix D). Therefore we construct our SMFs by adopting the z-evolving shape of the SMFs in the literature and normalize them according to the integrated CSFRD. The full description of this procedure and its verification on observational data can be found in Appendix D. Thus our assumed (SFG) SMFs and CSFRDs are consistent with each other at each redshift.

Based on the assumption that "all" SFGs exactly follow our gas fraction function (Equation (11)), and their number density obeys the SMF at each redshift, we can compute the molecular gas mass density by integrating the product of gas fraction, stellar mass, and SMF in each stellar mass bin at each redshift:

$$\begin{eqnarray}&&\begin{array}{l}{\rho }_{\mathrm{molgas}}(z)\,=\\ \displaystyle \sum _{{M}_{\star }\ \mathrm{bins}}{{\rm{\Phi }}}_{\mathrm{SMF}}(z,{M}_{\star })\times {M}_{\star }\times {\mu }_{\mathrm{molgas}}(z,{M}_{\star },{\rm{\Delta }}\mathrm{MS}).\end{array}\end{eqnarray} \tag{ 13 }$$

In Figure 14 we present the integrated cosmic cold molecular gas mass density versus redshift, using three different gas fraction functions ${\mu }_{\mathrm{molgas}}(z,{M}_{\star },{\rm{\Delta }}\mathrm{MS})$, our Equation (11) (orange solid line), T18 (green long-dashed line), and S17 (pink short-dashed line). The same SMFs are used for the three gas fraction functions.

Zoom In Zoom Out Reset image size

Note that the result is sensitive to the lower stellar mass limit down to which the integration is performed. Davidzon et al. (2017) adopt a lower limit of ${M}_{\star }={10}^{8.0}\ {M}_{\odot }$ when integrating SMFs to compute the cosmic stellar mass density. To match the CO blind deep-field data (e.g., Decarli et al. 2019; Riechers et al. 2019), we integrate only down to ${M}_{\star }={10}^{9.0}\ {M}_{\odot }$, that is, an order of magnitude shallower.

In Figure 14, we compare results from three recent CO blind deep-field surveys (Decarli et al. 2016, 2019 and Riechers et al. 2019 from the ASPECS-pilot, COLDz, and ASPECS-LP surveys, respectively) to the gas evolution curves derived from our, T18, and S17 functions. The form of the function ${\mu }_{\mathrm{molgas}}(z,{M}_{\star },{\rm{\Delta }}\mathrm{MS})$ significantly impacts the resulting cosmic cold gas mass density curve. Both our and the T18 functions provide very reasonable fits to the data without any tuning (except for the integration limit). Because the observational CO luminosity detection limit varies with redshift and sample (or excitation "correction") and is in general higher than the integration limits we chose, the currently available data cannot sufficiently constrain these functions.

The drawback of the SMF $\times {\mu }_{\mathrm{molgas}}$ integration in the previous section is that it only accounts for galaxies located exactly on the MS. In order to account for starburst galaxies as well as the scatter of the MS, we adopt here an alternative approach to derive the cosmic cold molecular gas mass density: we calculate for each mock galaxy (simulated under the 2-SFM framework by Béthermin et al. 2017) the cold molecular gas mass using our ${\mu }_{\mathrm{molgas}}$ function before summing them up within each redshift bin.

The "SIDES" simulation (Simulated Infrared Dusty Extragalactic Sky²⁴ ; Béthermin et al. 2017) generated 1,489,629 mock galaxies within a 2 deg² light cone from redshift 0.02 to 9.95. Different sets of SMFs were adopted according to redshift (Kelvin et al. 2014 for local galaxies; Moutard et al. 2016 at $z\lt 1.5$; Davidzon et al. 2017 at $1.5\lt z\lt 4$; and Grazian et al. 2015 at $z\gt 4$). Stellar masses were assigned to dark matter halos via abundance matching, and a certain recipe for the SFG fraction was assumed at each redshift. The modeling also accounts for the scatter of the star-forming MS coming from both the MS population itself and starburst galaxies, so it reasonably reproduces true galaxy distributions.

We use the SIDES mock galaxy catalog and select ${\mathrm{log}}_{10}{M}_{\star }\gt 9.0$ SFGs (as in the previous section to match CO luminosity function studies), and then we apply Equation (11) (as well as the T18 and S17 functions) to each galaxy to obtain its gas fraction, and hence to derive its molecular gas mass. We integrate the molecular gas mass for all galaxies in a given redshift bin, then divide it by the corresponding comoving volume to obtain the cosmic cold gas mass density ${\rho }_{\mathrm{mol}\mathrm{gas}}$. We sample the redshift range from 0 to 15 with 500 bins (i.e., bin size ∼0.00253 in ${\mathrm{log}}_{10}(1+z);$ as in the previous section).

The results are presented in Figure 15. The wiggling at the low-redshift end is likely due to the cosmic variance. At higher redshifts ($z\gt 0.5$), our curve coincidentally agrees with the semianalytic model (SAM) simulation by Popping et al. (2019). Other simulations, Obreschkow & Rawlings (2009) and Lagos et al. (2011), can be seen in Figure 5 of Riechers et al. (2019): at z ∼ 5, the Popping et al. (2019) simulation exhibits a 0.2 dex lower ${\rho }_{\mathrm{mol}\mathrm{gas}}$ than that of Lagos et al. (2011), and the Obreschkow & Rawlings (2009) ${\rho }_{\mathrm{mol}\mathrm{gas}}$ is 0.1 dex lower than Lagos et al. (2011), while the three are reversed at z ∼ 0.5, but still within 0.2 dex. Thus, in general, the simulations and the predictions with the functional form derived here are in good agreement.

Zoom In Zoom Out Reset image size

When using the T18 and S17 functions for the computation, the corresponding cold molecular gas mass density curves show a large difference. The S17 function leads to a much higher cold molecular gas mass density at all redshifts, which is likely because their function predicts significantly higher ${\mu }_{\mathrm{molgas}}$ (see Figure 13).²⁵ The T18 function results in fully (marginally) consistent cold molecular gas mass densities as our function at z ≲ 1 (z ∼ 2–3); however, it predicts 0.4–0.9 dex lower values at $z\gt 4$. This is mainly driven by the downturn of their ${\mu }_{\mathrm{molgas}}$ function at z ≳ 4 (as mentioned in Section 5) and probably also affected by their systematically lower ${\mu }_{\mathrm{molgas}}$ for low-mass galaxies (see Figure 13). Nevertheless, due to the large uncertainties in the CO blind deep-field data, it is still hard to distinguish whether our function is statistically better than the T18 function. We will further investigate this with simulated galaxies (Popping et al. 2019) in future work.

Meanwhile, all three gas evolution function based cosmic gas mass densities are consistent with the result of Klitsch et al. (2019), who obtained relatively shallow upper limits of about $1.6\times {10}^{8}\,{M}_{\odot }\,{\mathrm{Mpc}}^{-3}$ on the molecular gas mass densities at redshift z ∼ 0–1.7 using the quasar absorption line method. Deeper such observations in the future will also better verify the validity of our function.

In summary, the above comparisons indicate that our knowledge of the galaxy star-forming MS (e.g., the two-star formation model, 2-SFM; Sargent et al. 2014), SMFs (see references in Béthermin et al. 2017), and molecular gas fraction parameterization (using our functional form of ${\mu }_{\mathrm{molgas}}$ in Equation (11)) is moving toward a coherent picture.

The CO-to-H₂ conversion factor, ${\alpha }_{\mathrm{CO}}$, is an empirical ratio converting CO line luminosity to total molecular gas mass. It has been found to be relatively constant in the inner Galactic giant molecular clouds (GMCs), being around $4.6\ {M}_{\odot }{({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1}{\mathrm{pc}}^{2})}^{-1}$ (${X}_{\mathrm{CO}}=2.1\times {10}^{20}\ {\mathrm{cm}}^{-2}{({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1}{\mathrm{pc}}^{2})}^{-1}$; Solomon et al. 1987; Solomon & Barrett 1991; Solomon & Vanden Bout 2005), or $6.5\ {M}_{\odot }\ {({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1}{\mathrm{pc}}^{2})}^{-1}$ (${X}_{\mathrm{CO}}=3.0\,\times {10}^{20}\ {\mathrm{cm}}^{-2}\,{({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1}{\mathrm{pc}}^{2})}^{-1}$) when including heavy elements, which are mostly helium, while it is as low as $0.8\ {M}_{\odot }{({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1}{\mathrm{pc}}^{2})}^{-1}$ (${X}_{\mathrm{CO}}=2.0\,\times {10}^{20}\ {\mathrm{cm}}^{-2}\ {({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1}{\mathrm{pc}}^{2})}^{-1};$ $\pm 0.5\,\mathrm{dex}$) in local ULIRGs (Solomon et al. 1997; Downes & Solomon 1998; Solomon & Vanden Bout 2005). The calibration of ${\alpha }_{\mathrm{CO}}$ relies on a number of other (molecular) gas mass tracers, including virial mass, optically thin CO isotopologues, dust extinction, dust emission (via the gas-to-dust ratio, e.g., Section A.2), and diffuse γ-ray radiation. More details are given in the recent review by Bolatto et al. (2013). In this work, we only focus on the established ${\alpha }_{\mathrm{CO}}$–metallicity relations presented in Genzel et al. (2015) and T18, as shown in the left panel of Figure 16, to homogenize the molecular gas mass calculation for our complementary samples.

Zoom In Zoom Out Reset image size

The gas-to-dust mass ratio, ${\delta }_{\mathrm{GDR}}\equiv {M}_{\mathrm{total}\mathrm{gas}}/{M}_{\mathrm{dust}}$, describes the correlation between the total amount of gas (molecular plus atomic, compositing almost all of the ISM) and dust. As dust grains are usually assumed to be well mixed within the ISM, ${\delta }_{\mathrm{GDR}}$ should be predictable by ISM chemical models (e.g., see recent review by Galliano et al. 2018). We will skip the physical mechanism behind this and refer the reader to Galliano et al. (2018). Here we aim to understand how ${\delta }_{\mathrm{GDR}}$ can be applied for molecular gas mass estimation for high-redshift galaxies.

The calibration of ${\delta }_{\mathrm{GDR}}$ is usually based on observations of CO and H i emission lines plus multiwavelength photometry to which SED fitting is performed (e.g., Lisenfeld et al. 2000; Leroy et al. 2011; Magdis et al. 2011; Sandstrom et al. 2013; Rémy-Ruyer et al. 2014; Tan et al. 2014). These works found that ${\delta }_{\mathrm{GDR}}$ is correlated with galaxies' gas phase metallicity, as illustrated by data from our sample compilation (see Table 1) in Figure 16 (right panel). Note that ${\delta }_{\mathrm{GDR}}$ is around 100 for galaxies with solar and supersolar metallicity, while it increases nonlinearly toward lower metallicity, reaching over 1000 in extremely metal-poor ($\lt 10 \% \,{Z}_{\odot }$) galaxies (e.g., Elmegreen et al. 2013; Shi et al. 2014, 2016). The difference between the derived relations of Leroy et al. (2011) and Rémy-Ruyer et al. (2014) is about 0.1 dex in the supersolar-metallicity regime, increases to 0.2 dex at 0.2 solar metallicity, and then quickly becomes much larger at even lower metallicity.

As ${\delta }_{\mathrm{GDR}}$ is calibrated with total gas mass instead of molecular gas mass, a molecular-to-total gas mass ratio, ${f}_{\mathrm{mol}\mathrm{frac}}\equiv {M}_{\mathrm{mol}\mathrm{gas}}/{M}_{\mathrm{total}\mathrm{gas}}$, needs to be considered. It is discussed in the next section (Appendix A.3).

The molecular hydrogen fraction, ${f}_{\mathrm{mol}\mathrm{frac}}\equiv {M}_{\mathrm{mol}\mathrm{gas}}/{M}_{\mathrm{total}\mathrm{gas}}$, is the ratio between molecular gas and molecular+atomic gas. In the following, we use ${\mu }_{\mathrm{mol}\mathrm{frac}}\equiv {M}_{\mathrm{mol}\mathrm{gas}}/{M}_{\mathrm{atomic}\mathrm{gas}}$ for the molecular-to-atomic gas mass ratio.

Note that ${f}_{\mathrm{mol}\mathrm{frac}}$ (or ${\mu }_{\mathrm{mol}\mathrm{frac}}$) also correlates with metallicity; for example, the amount of dust grains, as hydrogen molecules form mainly on the surface of dust grains (e.g., Hollenbach & Salpeter 1971), and the abundance of dust grains depend on the metal enrichment by recent star formation activities (e.g., Draine 2003). The correlations between ${f}_{\mathrm{mol}\mathrm{frac}}$ and $12\,+{\mathrm{log}}_{10}({\rm{O}}/{\rm{H}})$ and ${M}_{\star }$ are illustrated in Figure 17, where theoretical models from Krumholz et al. (2009), Popping et al. (2014), and Davé et al. (2016) are compared to the data.

Zoom In Zoom Out Reset image size

In Figure 17, we show ${\mu }_{\mathrm{mol}\mathrm{frac}}$ versus metallicity and stellar mass with a large compilation of 524 galaxies from the literature (see labels and figure caption). All galaxies have ${M}_{\mathrm{mol}\mathrm{gas}}$ from CO observations, ${M}_{\mathrm{atomic}\mathrm{gas}}$ from H i observations, $12+{\mathrm{log}}_{10}({\rm{O}}/{\rm{H}})$ from optical spectroscopy, and ${M}_{\star }$ from multiwavelength optical/near-infrared data. The data points exhibit a large scatter in both panels, which is probably caused by the uncertainties in ${M}_{\mathrm{mol}\mathrm{gas}}$, ${M}_{\mathrm{atomic}\mathrm{gas}}$, and metallicity. The metallicity-dependent CO-to-H₂ conversion factor has on average a ∼38% uncertainty in the Saintonge et al. (2017) catalog, where the conversion factor is computed based on Accurso et al. (2017), and the observed CO line flux has a ∼5%–30% uncertainty. The H i line flux has a ∼2%–20% uncertainty in their catalog, and in addition the conversion from H i line flux to ${M}_{\mathrm{atomic}\mathrm{gas}}$ may have a 30% or higher uncertainty, due to the assumption of optically thin H i (e.g., Fukui et al. 2018). These uncertainties add up to at least ∼50%–60% uncertainty for the Y-axis.

Three theoretical models from Popping et al. (2014), Davé et al. (2016), and Krumholz et al. (2009) are overlaid as colored lines. Comparing with the data, the Krumholz et al. (2009) model provides the best fit at the low-metallicity end. At high metallicities, it seems the data are not statistically meaningful, and all three models provide reasonable predictions.

The figure shows that for local galaxies with solar-abundance metallicity and ${M}_{\star }\gt {10}^{10}\ {M}_{\odot }$, molecular gas nearly dominates the total gas mass ($\left\langle {f}_{\mathrm{mol}\mathrm{frac}}\right\rangle \gt 50 \%$). At higher redshifts, however, there is no observational constraint. We can only assume such scaling relations are still valid at higher redshifts. In principle, higher-redshift galaxies have higher SFRs and gas density at the same stellar mass (see Appendix A.5), and ${f}_{\mathrm{mol}\mathrm{frac}}$ should be at least as high as local analogs of similar stellar mass and metallicity. Thus it is common to ignore the atomic gas contribution in high-redshift galaxies with ${M}_{\star }\gt {10}^{10}\ {M}_{\odot }$ (e.g., T18).

The FMR (the correlation between metallicity, stellar mass, and SFR; Mannucci et al. 2010, 2011) and the mass–metallicity relation (MZR; the correlation between metallicity and stellar mass; e.g., Kewley & Ellison 2008) are usually used to infer the metallicity and metallicity-related properties (e.g., ${\alpha }_{\mathrm{CO}}$, ${\delta }_{\mathrm{GDR}}$) of high-redshift galaxies when no sufficient optical nebular emission line information is available. A number of FMRs and MZRs exist in the literature (see below), with metallicity ($12+{\mathrm{log}}_{10}({\rm{O}}/{\rm{H}})$) parameterized as a function of ${M}_{\star }$ and/or z or SFR. However, whether the FMR and MZRs are valid across cosmic time or within a given stellar mass range is less studied.

Here we take the following seven FMRs and MZRs most widely used for high-redshift studies and compared them in Figure 18 so that their validities can be more clearly seen in bins of redshift:

• 1.  
  Mannucci et al. (2010) presented their FMR in their Equations (2) and (4) for local galaxies. Their metallicity values are originally calibrated from optical emission lines following the Maiolino et al. (2008, M08) prescription instead of Pettini & Pagel (2004, PP04), so we converted their derived metallicity to the PP04 calibration by solving both Mannucci et al. (2010) Equation (1) (the polynomial form, instead of their Equation (2)) and Maiolino et al. (2008) Equation (1) (with their Table 4's second row of coefficients). Both their Equations (2) and (4) are shown in Figure 18. The caveat of their Equation (2) includes that (a) at a given redshift and SFR, it first increases with stellar mass and then drops quickly when ${\mathrm{log}}_{10}{M}_{\star }\gt 11.2;$ and (b) it predicts the lowest metallicity for starburst galaxies at $z\gt 1$. And the caveat of their Equation (4) is the nonphysical extrapolation for (a) MS galaxies at all redshift with ${\mathrm{log}}_{10}{M}_{\star }\gt 11.2$, and (b) $z\gtrsim 1$ MS galaxies with ${\mathrm{log}}_{10}{M}_{\star }\lt 9.5$.
• 2.  
  Mannucci et al. (2011) provided an updated version of their FMR in Mannucci et al. (2010) for lower-mass galaxies. The update is only for the ${\mathrm{log}}_{10}{M}_{\star }-0.32\times {\mathrm{log}}_{10}\mathrm{SFR}\lt 9.5$ case, that is, low mass or high SFR. Therefore, we show this equation only in the bottom panels for starburst galaxies. Note that there is a nonphysical jump in metallicity at ${\mathrm{log}}_{10}{M}_{\star }\sim 10.2\mbox{--}10.5$.
• 3.  
  Equation 12(a) in Genzel et al. (2015) is originally from Wuyts et al. (2014) and was also adopted by T18 in the identical form. This formula considers both redshift and stellar mass as the parameters determining metallicity. It predicts reasonable metallicities except at ${\mathrm{log}}_{10}{M}_{\star }\gtrsim 11$ at local (or ${\mathrm{log}}_{10}{M}_{\star }\gtrsim 11.5\mbox{--}12$ at $z\gtrsim 1$). This motivates our modification of this equation as described in Equation (6).
• 4.  
  Equation 12(b) in Genzel et al. (2015) is based on Equation (4) in Mannucci et al. (2010), with a M08-to-PP04 conversion applied by solving both Mannucci et al. (2010) Equation (2) and Maiolino et al. (2008) Equation (1) (see their Table 4's second row of coefficients). Therefore, the curve of this equation is very similar to the Mannucci et al. (2010) Equation (4) curve, yet we caution that the M08-to-PP04 conversion is different in their work than here, and our conversion (by solving Mannucci et al. (2010) Equation (2) and Maiolino et al. 2008 Equation (1)) should be more precise.
• 5.  
  Equation (5) in Magnelli et al. (2012) uses the Denicoló et al. (2002) calibration. Thus we convert this calibration to the PP04 N2 calibration following Kewley & Ellison (2008). They assumed two different MZRs distinguished by redshift at 1.5. We caution that it is much lower at low redshift ($z\lesssim 1$). It also always predicts subsolar metallicity for galaxies at $z\gt 1$.
• 6.  
  Kewley & Ellison (2008) use PP04 O3N2 and N2 MZRs as listed in their Table 2.²⁶ Their equation only depends on stellar mass and has no redshift evolution, so we only show their curve in the first panel. It predicts a too-high metallicity for high-redshift galaxies with ${\mathrm{log}}_{10}{M}_{\star }\,\lt 10.5$, and like Mannucci et al. (2010) Equation (2), it also has a nonphysical drop with increasing stellar mass when ${\mathrm{log}}_{10}{M}_{\star }\gt 11.2$.
• 7.  
  See Maiolino et al. (2008) Equation (2) with the coefficients listed in their Table 5. They fitted five different MZRs at five redshifts they analyzed. Here we linearly interpolate their coefficients in redshift so as to plot their curves in Figure 18. The equation seems reasonable at low z ($z\lt 2\mbox{--}3$) but predicts significantly subsolar metallicity at $z\gtrsim 5$ even for starbursts.

Zoom In Zoom Out Reset image size

Figure 18 shows that most of these formulae are consistent (within 0.2 dex) only for MS galaxies at $z\lesssim 2$ and with ${\mathrm{log}}_{10}{M}_{\star }\sim 10\mbox{--}11.2$. Subtle differences exist among these curves, and the reader should consider the proper choice of MZR or FMR to use. A small difference, such as a 0.1 dex lower/higher metallicity, could translate into a factor of 1.6 higher/lower ${\delta }_{\mathrm{GDR}}$ when assuming the ${\delta }_{\mathrm{GDR}}$–metallicity relation in Figure 16. Also note that the prescription for deriving $12+{\mathrm{log}}_{10}({\rm{O}}/{\rm{H}})$ from optical emission lines is important, as can be seen by comparing the first and third formulae at the high-mass end. Finally, we caution that these formulae do not agree well at the low-mass regime (${M}_{\star }\lt {10}^{10.0}\ {M}_{\odot }$). But for the study in this work, although with such a large ALMA sample, we still do not probe such low-mass galaxies. Therefore these discrepancies are currently not an issue.

In Section 5 we mentioned that we adopt the Speagle et al. (2014) MS (the #49 fitting in their Table 7) with cosmic time as the variable. In Figure 19 we compare a number of MS relations in the literature (see the labels therein). The Whitaker et al. (2014), Lee et al. (2015), and Tomczak et al. (2016) MS are only valid at $z\lesssim 3$. The Sargent et al. (2014) MS predicts the lowest sSFR_MS at $z\gt 5$, while the Béthermin et al. (2015) MS predicts a factor of 2 higher sSFR_MS than average at $z\gt 5$ and is in general higher for ${\mathrm{log}}_{10}{M}_{\star }\gt {10}^{11}$ galaxies. The Pearson et al. (2018) MS is a factor of $\lt 2$ lower than others at z ∼ 2–3 and in general lower for ${\mathrm{log}}_{10}{M}_{\star }\lt {10}^{10}$ galaxies. These MS calibrations have large scatter in the ${\mathrm{log}}_{10}{M}_{\star }\lt {10}^{9}$ and ${\mathrm{log}}_{10}{M}_{\star }\gtrsim {10}^{12}$ regimes, which lack observational data. The Speagle et al. (2014) MS (with cosmic age) is closer to the average of all MS analyzed, so we adopt it for our work.

Zoom In Zoom Out Reset image size

The Leslie et al. (2019) MS is potentially an alternative choice for a most reasonable MS to use. They derived the MS correlation from redshift ∼0.2 to ∼6 by stacking the VLA 3 GHz large program data (Smolčić et al. 2017; covering 2 sq. deg. COSMOS field with a sensitivity of $1\,\sigma \sim 2.3\,\mu \mathrm{Jy}/\mathrm{beam}$ at a spatial resolution of 075) using a large sample of ∼300,000 galaxies from the Laigle et al. (2016) and Davidzon et al. (2017) catalogs. The largest difference between the Leslie et al. (2019) MS and the Speagle et al. (2014) one is that the former exhibits a flattening for a higher stellar mass, while the latter is a straight line at each redshift. Such a flattening, yet debated, has also been reported by other stacking studies (e.g., Lee et al. 2015; Schreiber et al. 2015; Tomczak et al. 2016). We refer the reader to these papers and Leslie et al. (2019) for details on the shapes of different MS functions. Here we investigate how an MS with flattening affects our functional form by adopting the Leslie et al. (2019) MS and repeating our analysis from Sections 5 to 6.2. In Figure 20, we compare the obtained ${\mu }_{\mathrm{molgas}}$ evolution curves using the Leslie et al. (2019) MS to those using the Speagle et al. (2014) MS.

Zoom In Zoom Out Reset image size

Figure 20 shows that adopting the Leslie et al. (2019) MS results in an at most factor of two higher gas fraction for low-mass, high-redshift galaxies (${\mathrm{log}}_{10}{M}_{\star }\lesssim 9.2$, z ∼ 2–6), while being indistinguishable from using the Speagle et al. 2014 MS for galaxies with ${\mathrm{log}}_{10}{M}_{\star }\gtrsim 10\mbox{--}11$. The difference at the low-mass end is likely because the Leslie et al. (2019) MS predicts two times higher sSFRs for low-mass galaxies at $z\sim 0$ while a factor of two lower sSFR at $z\sim 3$, as shown in the left panel of Figure 19. This leads to a systematically lower ${\rm{\Delta }}\mathrm{MS}$ for low-mass galaxies at $z\sim 0$, altering the slope of ${\mu }_{\mathrm{molgas}}$ versus ${\rm{\Delta }}\mathrm{MS}$ to be shallower (by a small change of 0.04 in the coefficient ${\mathsf{a}}$ in Equation (11)), meanwhile steepening the slope of ${\mu }_{\mathrm{molgas}}$ versus ${\mathrm{log}}_{10}{M}_{\star }$ (by a change of 0.14 in the coefficient ${\mathsf{b}}$ in Equation (11)). Thus it results in a $\lt 2\times$ higher extrapolation for the gas fraction at the low-mass end. We note that the fits are indistinguishable where data are rich; that is, using either the Leslie et al. (2019) or Speagle et al. (2014) MS makes no obvious difference for ${\mathrm{log}}_{10}{M}_{\star }\gtrsim 10\mbox{--}11$ galaxies at all redshifts (z ∼ 0–6).

In Figure 19, we additionally show the MS relations used in S17 following the equations in their Sections 2.1 and 2.2. Although their equation is not aimed at extrapolating out to $z\gt 3\mbox{--}4$ (and exhibits a large excess compared to others), here we show their curve and use their MS to compute ${\rm{\Delta }}\mathrm{MS}$ for our data out to $z\sim 6$ for the sole purpose of evaluating the validity of their MS at these redshifts. In Figure 20 we also repeated the fitting with our Equation (11) functional form and used the S17 MS for ${\rm{\Delta }}\mathrm{MS}$ normalization. The implied gas fraction evolution curve is $\sim 2\times$ ($\sim 4\times$) higher than that using the Leslie et al. (2019) MS (Speagle et al. 2014 MS) at the low-mass end, meanwhile it also exhibits a $\sim 2.5\times$ lower gas fraction at the massive end. This means that the difference in MS can indeed explain about half of the difference between our and the S17 molecular gas mass density curves seen in Figure 13 (the shape of the functional form is likely responsible for the other half of the difference).