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

Daizhong Liu, E. Schinnerer, B. Groves, B. Magnelli, P. Lang, S. Leslie, E.

Jimenez-Andrade, D.A. Riechers, G. Popping, Georgios E. Magdis, et al.

► To cite this version:

Daizhong Liu, E. Schinnerer, B. Groves, B. Magnelli, P. Lang, et al. Automated Mining of the ALMA Archive in the COSMOS Field (A3COSMOS): II. Cold Molecular Gas Evolution out to Redshift 6. The Astrophysical Journal, 2019, 887 (2), pp.235. 10.3847/1538-4357/ab578d hal-02433107

HAL Id: hal-02433107

https://hal.science/hal-02433107v1

Submitted on 13 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

Daizhong Liu<fnref n="1" />, E. Schinnerer<fnref n="1" />, B. Groves<fnref n="2" />,<fnref n="3" />, B. Magnelli<fnref n="4" />, P. Lang<fnref n="1" />, S. Leslie<fnref n="1" />, E. Jiménez-Andrade<fnref n="4" />,<fnref n="5" />,
D. A. Riechers<fnref n="1" />,<fnref n="6" />, G. Popping<fnref n="1" />,<fnref n="7" />, Georgios E. Magdis<fnref n="8" />,<fnref n="9" />,<fnref n="10" />,<fnref n="11" />, E. Daddi<fnref n="12" />, M. Sargent<fnref n="13" />, Yu Gao<fnref n="14" />,
Y. Fudamoto<fnref n="15" />, P. A. Oesch<fnref n="15" />, and F. Bertoldi<fnref n="4" />

Received 2019 September 11; revised 2019 November 11; accepted 2019 November 12; published 2019 December 23

## Abstract

We present new measurements of the cosmic cold molecular gas evolution out to redshift 6 based on systematic mining of the Atacama Large Millimeter/submillimeter Array (ALMA) public archive in the COSMOS deep field (A³COSMOS). Our A³COSMOS data set contains ~700 galaxies ($0.3 \lesssim z \lesssim 6$) with high-confidence ALMA detections in the (sub)millimeter continuum and multiwavelength spectral energy distributions. Multiple gas mass calibration methods are compared, and biases in band conversions (from observed ALMA wavelength to rest-frame Rayleigh-Jeans tail continuum) have been tested. Combining our A³COSMOS sample with ~1000 CO-observed galaxies at $0 \lesssim z \lesssim 4$ (75% at $z < 0.1$), we parameterize galaxies' molecular gas depletion time ($\tau_{depl}$) and molecular gas to stellar mass ratio ($\mu_{molgas}$) each as a function of the stellar mass ($M_*$), offset from the star-forming main sequence (ΔMS) and cosmic age (or redshift). Our proposed functional form provides a statistically better fit to current data (than functional forms in the literature) and implies a “downsizing” effect (i.e., more-massive galaxies evolve earlier than less-massive ones) and “mass quenching” (gas consumption slows down with cosmic time for massive galaxies but speeds up for low-mass ones). Adopting galaxy stellar mass functions and applying our $\mu_{molgas}$ function for gas mass calculation, we for the first time infer the cosmic cold molecular gas density evolution out to redshift 6 and find agreement with CO blind surveys as well as semianalytic modeling. These together provide a coherent picture of cold molecular gas, star formation rate, and stellar mass evolution in galaxies across cosmic time.

Unified Astronomy Thesaurus concepts: Galaxy evolution (594); High-redshift galaxies (734); Interstellar medium (847); Submillimeter astronomy (1647)

Supporting material: figure set, interactive figure

## 1. Introduction

The interstellar medium (ISM), especially the cold molecular gas, is the fuel of star formation activity in galaxies. In recent years, our knowledge of the cosmic evolution of star formation and stellar mass growth has been obtained out to redshift ~5 (e.g., see the latest reviews by Lutz 2014 and Madau & Dickinson 2014; see also Davidzon et al. 2017; Liu et al. 2018; to name a few). However, the cosmic evolution of the cold molecular gas is much less well constrained, and the validity of different gas tracers is debated (e.g., Magdis et al. 2012a; Santini et al. 2014; Genzel et al. 2015; Tacconi et al. 2018; Decarli et al. 2019; Riechers et al. 2019).

There are several widely used tracers of the molecular gas content in galaxies, including the commonly used carbon monoxide (CO) rotational transition lines, dust masses from the

dust spectral energy distribution (SED), and the cold dust continua at the Rayleigh-Jeans (RJ) tail of the dust SED. We introduce each case below.

Observationally, CO lines at rest-frame millimeter wavelengths have been established as the most commonly used tracers of total molecular gas content in galaxies near and far since the 1970s (e.g., see the latest reviews by Carilli & Walter 2013; Combes 2018). At high redshift, this method relies on galaxy samples with accurate spectroscopic redshifts and usually has uncertainties from the CO-to-H₂ conversion factor ($\alpha_{CO} \equiv M_{molgas}/L'_{CO}$) and CO excitation. With this method, Genzel et al. (2010) and Tacconi et al. (2013, 2018) conducted the largest survey for individual galaxies (named PHIBSS) by observing hundreds of star-forming galaxies (SFGs) at $z \sim 1-3$ to study the molecular gas scaling relation

<footnote n="1">¹Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany; dzliu@mpia.de</footnote>

<footnote n="2">²Research School of Astronomy and Astrophysics, Australian National University, Canberra ACT, 2611, Australia</footnote>

<footnote n="3">³International Centre for Radio Astronomy Research, University of Western Australia, Crawley, Perth, Western Australia, 6009, Australia</footnote>

<footnote n="4">⁴Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, D-53121 Bonn, Germany</footnote>

<footnote n="5">⁵National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA</footnote>

<footnote n="6">⁶Department of Astronomy, Cornell University, Space Sciences Building, Ithaca, NY 14853, USA</footnote>

<footnote n="7">⁷European Southern Observatory, Karl-Schwarzschild-Strasse 2, D-85748, Garching, Germany</footnote>

<footnote n="8">⁸Cosmic Dawn Center at the Niels Bohr Institute, University of Copenhagen and DTU-Space, Technical University of Denmark, Denmark</footnote>

<footnote n="9">⁹DTU Space, National Space Institute, Technical University of Denmark, Elektrovej 327, DK-2800 Kgs. Lyngby, Denmark</footnote>

<footnote n="10">¹⁰Niels Bohr Institute, University of Copenhagen, DK-2100 Copenhagen Ø, Denmark</footnote>

<footnote n="11">¹¹Institute for Astronomy, Astrophysics, Space Applications and Remote Sensing, National Observatory of Athens, 15236, Athens, Greece</footnote>

<footnote n="12">¹²CEA Saclay, Laboratoire AIM-CNRS-Université Paris Diderot, Irfu/SAP, Orme des Merisiers, F-91191, Gif-sur-Yvette, France</footnote>

<footnote n="13">¹³Astronomy Centre, Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH, UK</footnote>

<footnote n="14">¹⁴Purple Mountain Observatory & Key Laboratory for Radio Astronomy, Chinese Academy of Sciences, 10 Yuanhua Road, Nanjing 210033,
People's Republic of China</footnote>

<footnote n="15">¹⁵Department of Astronomy, Université de Genève, Chemin des Maillettes 51, 1290 Versoix, Switzerland</footnote>
and evolution. Meanwhile, Walter et al. (2016) and Decarli et al. (2016, 2019) have been conducting the largest blank-field survey (named ASPECS) by scanning a range of millimeter spectra within a fixed sky area to determine the CO luminosity function and thereby study the molecular gas mass density evolution.

Alternatively, in the past few years, emission from dust grains located in the star-forming regions of galaxies has also been widely used as a proxy for the ISM. These dust grains absorb rest-frame ultraviolet (UV) photons from massive young stars and re-emit thermal radiation in the infrared (IR) to millimeter wavelengths. By fitting a galaxy's full dust SED with models, such as modified blackbody models or multi-component physical models (e.g., Draine & Li 2007), the dust mass and dust temperature (or interstellar radiation field) can be obtained (e.g., Santini et al. 2010, 2014; Magdis et al. 2011, 2012a; Magnelli et al. 2012, 2014; Saintonge et al. 2013; Sandstrom et al. 2013; Tan et al. 2014; Béthermin et al. 2015; Berta et al. 2016; Hunt et al. 2019). The dust mass can then be converted to gas mass via the application of empirical gas-to-dust ratios ($\delta_{GDR}$).

However, a galaxy's full dust SED is a composite of a variety of dust components with different temperatures. Warmer dust exposed to strong radiation fields (e.g., photodominated regions; Dale et al. 2001; Draine & Li 2007) globally outshines the colder dust at shorter wavelengths of the SED, but the former is much less abundant (e.g., <10% in mass) and does not represent the bulk of dust in a galaxy. Thus, obtaining a reliable dust mass usually requires longer wavelength coverage that includes the RJ tail (e.g., $\lambda_{rest} \gtrsim 250$ µm). Also, different dust SED models can result in strong and not easily predictable systematic effects (Berta et al. 2016). Therefore, instead of fitting the SED, the RJ-tail method has been proposed by Scoville (2013) and Scoville et al. (2014), which directly uses the RJ-tail dust continuum to trace gas (yet the underlying physics of using dust mass to trace gas mass is the same as in the dust SED method).

The RJ-tail method has recently been proven to be as reliable as the CO method (e.g., Scoville et al. 2014, 2016; Groves et al. 2015; Hughes et al. 2017; Bertemes et al. 2018; Saintonge et al. 2018; Kaasinen et al. 2019; and theoretical works, e.g., Liang et al. 2018, 2019; Privon et al. 2018) and is much more efficient in surveying large galaxy samples at high redshift. This method relies on the assumption that the dust grains providing most of the dust mass in galaxies are cold and mixed within the ISM. Their temperatures are likely always as cold as $T_{dust} \approx 25$ K (see Scoville et al. 2014, 2016), and hence they can trace the total gas content via a relatively stable gas-to-dust ratio ($\delta_{GDR}$; e.g., Leroy et al. 2011; Rémy-Ruyer et al. 2014). Yet we bear in mind that metallicity, true dust temperature, and mass distributions are all unsolved issues.

These studies have led to a rough picture of dust and gas evolution from redshift 3 to present, where (a) the fraction of molecular gas mass to the total of molecular gas and stellar masses,

$$ 
\begin{array}{l}
f_{\text{mol gas}} \equiv M_{\text{mol gas}}/(M_{\text{mol gas}} + M_{*}), \\
\text{or } \mu_{\text{molgas}} \equiv M_{\text{mol gas}}/M_{*},
\end{array}
\qquad (1) $$

decreases with cosmic age from $z \sim 3$ to $z \sim 0$ and depends on star formation rate (SFR) and stellar mass; and (b) the molecular gas depletion time,

$$ \tau_{\text{depl}} \equiv M_{\text{mol gas}}/\text{SFR}, \qquad (2) $$

increases from $z \sim 3$ to present and is significantly different between typical SFGs (which follow a tight $M_* - \text{SFR}$ main sequence (MS) at each redshift; e.g., Brinchmann et al. 2004; Elbaz et al. 2007; Daddi et al. 2007; Noeske et al. 2007) and starbursts (i.e., located significantly above the MS; e.g., Rodighiero et al. 2011, 2014).

Genzel et al. (2015) first compiled a large sample of local and high-redshift ($0 < z < 3$) galaxies with both CO (500 galaxies) and dust SED (512 galaxies) methods. They studied the gas scaling relations by characterizing $f_{\text{mol gas}}$ or $\mu_{\text{molgas}}$ and $\tau_{\text{depl}}$ as functions of $M_*$, SFR, and redshift. More specifically, they found that gas fraction and depletion time are more strongly correlated with the SFR offset to the MS,

$$ 
\begin{array}{l}
\delta\text{MS} \equiv \text{SFR}/\text{SFR}_{\text{MS}}, \\
\text{or } \Delta\text{MS} \equiv \log_{10}(\text{SFR}/\text{SFR}_{\text{MS}}),
\end{array}
\qquad (3) $$

rather than the absolute SFR.

Utilizing the RJ-tail dust continuum method (at rest-frame 850 µm), Scoville et al. (2017, hereafter S17) studied the gas ($f_{\text{mol gas}}$ or $\mu_{\text{molgas}}$, and $\tau_{\text{depl}}$) scaling relations with a large sample of 708 high-redshift Herschel far-IR-selected galaxies ($0.3 < z < 4.5$), including a large amount of public data in the Atacama Large Millimeter/submillimeter Array (ALMA) archive (at a 2.5–3σ detection threshold), and characterized the $\mu_{\text{molgas}}$ and $\tau_{\text{depl}}$ functional forms:

$$ 
\begin{array}{l}
\mu_{\text{molgas}} = 0.71 \times (\delta\text{MS})^{+0.32} \times M_{*,10}^{-0.70} \times (1+z)^{+1.84}, \\
\tau_{\text{depl}} = 3.23 \text{ Gyr} \times (\delta\text{MS})^{-0.70} \times M_{*,10}^{-0.01} \times (1+z)^{-1.04},
\end{array}
\qquad (4) $$

where $M_{*,10}$ is $M_*/(10^{10} M_\odot)$. With the same method but at rest-frame 250–500 µm, Schinnerer et al. (2016) studied a smaller sample of optically selected galaxies at $z = 2.8–3.6$. However, discrepancies exist due to the slightly different methods and samples.

Tacconi et al. (2018, hereafter T18) expanded the work of Genzel et al. (2015) by obtaining nearly 100 new CO detections in the PHIBSS2 survey and compiling more samples of local to high-redshift galaxies in the literature. They used all three methods for obtaining molecular gas measurements for 1444 galaxies at $0 < z < 4$ and fitted them all together to derive the $\mu_{\text{molgas}}$ and $\tau_{\text{depl}}$ functions:

$$ 
\begin{array}{l}
\mu_{\text{molgas}} = 2.32 \times (\delta\text{MS})^{+0.53} \times M_{*,10}^{-0.35} \\
\qquad \times 10^{-3.62 \times (\log_{10}(1+z)-0.66)^2}, \\
\tau_{\text{depl}} = 1.06 \text{ Gyr} \times (\delta\text{MS})^{-0.44} \times M_{*,10}^{+0.09} \\
\qquad \times (1+z)^{-0.62},
\end{array}
\qquad (5) $$

where we adopted their $\beta = 2$ best fit with the Speagle et al. (2014) MS and expressed their stellar mass in $M_{*,10}$ to match Equation (4).

Comparing Equations (4) and (5) at redshift 3 and $M_* = 5 \times 10^{10} M_\odot$ reveals a factor of 2.3 difference in $\mu_{\text{molgas}}$ and a factor of 1.5 in $\tau_{\text{depl}}$. Such noticeable differences exist for other parameter values as well, raising concerns on the validity of the $\mu_{\text{molgas}}$ and $\tau_{\text{depl}}$ functions and the predictability of $\mu_{\text{molgas}}$ and $\tau_{\text{depl}}$ from a galaxy's redshift, stellar mass, and SFR properties. In addition, previous works have constraints only for $z \lesssim 3–4$.

To solve the discrepancies and understand systematic biases especially for the latest RJ-tail dust method, a large, robust
galaxy sample from local to high redshift is needed to carry out the comprehensive analysis. Therefore, in this work, we present an independent study on the characterization of the molecular gas fraction (μ<sub>molgas</sub>) and depletion time (τ<sub>depl</sub>) functional forms utilizing a large (~700), robust galaxy sample at 0.3 ≲ z ≲ 6 in the 2 deg<sup>2</sup> COSMOS field (Scoville et al. 2007) from the A<sup>3</sup>COSMOS project,<sup>16</sup> together with ~1000 CO-detected galaxies at 0 ≲ z ≲ 4 (75% at z < 0.1) from recent large surveys in the literature. All A<sup>3</sup>COSMOS galaxies have robust (sub)millimeter continuum detections from public ALMA archival data (release date up to 2018 August 1) with an expected spurious fraction close to zero and the flux bias corrected statistically (Liu et al. 2019; hereafter Paper I).

With such a combined large sample, we provide new molecular gas fraction (μ<sub>molgas</sub>) and depletion time (τ<sub>depl</sub>) functional forms that are valid from redshift 0 to 6. We adopt galaxies’ stellar mass functions (SMFs) or realistic galaxy modeling to analytically derive the cosmic molecular gas mass density evolution for the first time with such a large data set out to redshift 6. The result supports a coherent picture of the evolution of galaxies’ stellar mass, star formation, and cold molecular gas.

This paper is organized as follows. Galaxy samples are presented in Section 2, with the A<sup>3</sup>COSMOS high-redshift sample in Section 2.1 and complementary local-to-high-redshift samples from the literature in Section 2.2. Molecular gas mass calculation and comparison are presented in Section 3 (dust SED method in Section 3.1, RJ-tail method in Section 3.2, and comparison in Section 3.3). The complexity and apparent correlations between μ<sub>molgas</sub>, τ<sub>depl</sub>, and galaxies’ redshifts, stellar masses, and SFRs are discussed in Section 4. The characterizations of the functional forms for μ<sub>molgas</sub> and τ<sub>depl</sub> are presented in Section 5, and their implications are discussed in Section 6. Finally, the cosmic evolution of molecular gas mass density is analytically obtained in Section 7, followed by the summary in Section 8.

In the appendices, we thoroughly compare several important correlations related to our analysis: CO-to-H<sub>2</sub> conversion factor versus metallicity in Appendix A.1; gas-to-dust ratio versus metallicity in Appendix A.2; molecular to total gas fraction versus stellar mass or metallicity in Appendix A.3; and stellar mass-metallicity relation in Appendix A.4. These comparisons give useful insights into how different correlations impact the results presented in this work, as well as supporting our fiducial model used in this work.

We adopt a flat ΛCDM cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω<sub>M</sub> = 0.3, and Λ₀ = 0.7,<sup>17</sup> and a Chabrier (2003) initial mass function (IMF).

## 2. Sample and Data

### 2.1. The A<sup>3</sup>COSMOS Galaxy Sample

In Paper I we presented the A<sup>3</sup>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 <0.′′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 (<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 ~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<sup>3</sup>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<sup>3</sup>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 (S/N) sources that have a <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 S/N > 2. Betti et al. (2019) analyzed ALMA continuum data for 101 galaxies and selected 68 as detections with an aperture-based total flux of S/N > 2 or peak flux of S/N > 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 S/N selection; however, only eight galaxies have a peak flux S/N > 4.35, which is our selection criterion based on statistics (corresponding to a spurious fraction ~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<sup>3</sup>COSMOS from a (peak flux) S/N of 4.35 to 3.0 doubles the A<sup>3</sup>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<sup>3</sup>COSMOS.

<footnote n="16"><sup>16</sup> https://sites.google.com/view/a3cosmos</footnote>

<footnote n="17"><sup>17</sup> Same as those adopted by T18.</footnote>
Galaxy properties in the $A^3\text{COSMOS}$ galaxy catalog, including stellar mass ($M_*$), IR luminosity ($L_{IR}$), and dust mass ($M_{\text{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: $\text{SFR} = L_{IR}/1 \times 10^{10} M_\odot \text{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\%$).

## 2.2. Complementary Local-to-high-redshift Galaxy Samples

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).

Saintonge et al. (2017) applied a metallicity-dependent $\alpha_{CO}$ according to Accurso et al. (2017) to convert their CO observations into molecular gas mass. A similar metallicity-dependent $\alpha_{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_{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_{CO}$ following T18. We use metallicity to calculate $\alpha_{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_{CO}$. The recomputed molecular gas masses are within a factor of 2 ($\lesssim 0.36$ dex in logarithm) from their originally obtained values.

# 3. Molecular Gas Mass Calculation

We summarize the three most commonly used molecular gas mass calibration methods for high-redshift galaxies in Figure 1.

As mentioned in the introduction, they are (a) CO lines, (b) SED-fitted dust mass, and (c) RJ-tail dust continuum.<fnref n="18" /> The CO method infers the molecular gas mass via the $\alpha_{CO}$ conversion factor, which relates CO luminosity to $H_2$ gas mass and is correlated with metallicity (see details in Appendix A.1). When the observed CO line is not the ground transition ($J = 1 \rightarrow 0$), an excitation ladder is needed to convert the higher-$J$ line luminosity to the $J = 1 \rightarrow 0$ one (e.g., Carilli & Walter 2013).

Given that the CO and dust RJ-tail 850 µm based gas mass calibrations have been extensively verified to be tightly correlated in a number of recent works at local and high redshift up to $z \sim 2-3$ (e.g., Scoville et al. 2014, 2017; Hughes et al. 2017; Bertemes et al. 2018; Saintonge et al. 2018; Kaasinen et al. 2019), we do not further discuss the CO method here, but focus on popular dust-based methods. In Section 3.1 we describe the use of SED-fitted dust mass to compute molecular gas mass, and in Section 3.2 we describe the use of the RJ-tail dust continuum for molecular gas mass calculations. There are multiple choices for calibration factors and wavelengths, so we compare these methods thoroughly in Section 3.3.

Later we will combine CO- and dust-based samples for our data fitting analysis (in Sections 4 and 5). We assume that the consistency between CO- and (our adopted) dust-based gas mass calibration extends to all galaxies in our combined sample, which is at least supported by the aforementioned CO and dust calibration studies (but we also discuss the current caveats at the end of Section 5).

## 3.1. Molecular Gas Mass from SED-fitted Dust Mass

In the SED-fitted dust mass method, we first obtain the dust mass ($M_{\text{dust}}$), dust mean temperature ($T_{\text{dust}}$), and dust emissivity ($\dot{\beta}_{\text{dust}}$; describing the dust opacity $\kappa$'s wavelength dependency; e.g., Li & Draine 2001) from optical-to-millimeter SED fitting; we then apply a gas-to-dust ratio, $\delta_{\text{GDR}} \equiv M_{\text{total gas}}/M_{\text{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_{IR} > 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_{\text{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_* > 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

<footnote n="18">¹⁸ The (b) and (c) methods have the same underlying physics, which uses dust mass to trace the total gas mass. Here we separate them because they have different technical steps and assumptions. See details in Sections 3.1 and 3.2.</footnote>
<table><caption>**Table 1**
Galaxy Samples Used for This Study</caption><thead><tr><th>Sample Name</th><th>z</th><th>log<sub>10</sub>(M<sub>*</sub>/M<sub>☉</sub>)</th><th>N<sub>det.</sub><sup>a</sup></th><th>References</th></tr></thead><tbody><tr><td>A<sup>3</sup>COSMOS (Paper I)</td><td>0.29–5.667<sup>b</sup></td><td>~10.0–12.0</td><td>658<sup>b</sup></td><td>Liu et al. (2019, data set 20180801)</td></tr><tr><td>DGS</td><td>0–0.045</td><td>6.5–10.6</td><td>32</td><td>Rémy-Ruyer et al. (2014, 2015)</td></tr><tr><td>HRS</td><td>0.0034–0.006</td><td>8.7–11.7</td><td>99</td><td>Andreani et al. (2018, and refs. therein<sup>c</sup>)</td></tr><tr><td>KINGFISH</td><td>0.0005–0.006</td><td>6.3–10.7</td><td>28</td><td>Groves et al. (2015, and refs. therein<sup>d</sup>)</td></tr><tr><td>Saintonge+2017 (xCOLD GASS)</td><td>0.01–0.05</td><td>9.0–11.4</td><td>330</td><td>Saintonge et al. (2017)</td></tr><tr><td>Cicone+2017 (ALLSMOG)</td><td>0.01–0.03</td><td>8.5–11.5</td><td>48</td><td>Cicone et al. (2017)</td></tr><tr><td>Lisenfeld+2017</td><td>0.01–0.07</td><td>9.0–11.3</td><td>41</td><td>Lisenfeld et al. (2017)</td></tr><tr><td>Cortzen+2019</td><td>0.03–0.29</td><td>8.9–11.5</td><td>51</td><td>Cortzen et al. (2019)</td></tr><tr><td>Villanueva+2017 (VALES)</td><td>0.03–0.33</td><td>10.1–11.3</td><td>49</td><td>Villanueva et al. (2017)</td></tr><tr><td>Bertemes+2018 (Stripe82)</td><td>0.03–0.20</td><td>10.0–11.3</td><td>78</td><td>Bertemes et al. (2018)</td></tr><tr><td>Kirkpatrick+2014 (5MUSE)</td><td>0.05–0.29</td><td>10.5–11.4</td><td>17</td><td>Kirkpatrick et al. (2014)</td></tr><tr><td>Bauermeister+2013 (EGNOG)</td><td>0.06–0.31</td><td>10.7–11.5</td><td>14</td><td>Bauermeister et al. (2013)</td></tr><tr><td>Lee+2017</td><td>0.27–0.62</td><td>10.0–11.1</td><td>20</td><td>Lee et al. (2017)</td></tr><tr><td>Spilker+2018</td><td>0.60–0.75</td><td>~11.0</td><td>4</td><td>Spilker et al. (2018)</td></tr><tr><td>Combes+2013 (ULIRGs)</td><td>0.61–0.97</td><td>9.3–12.0</td><td>12</td><td>Combes et al. (2013)</td></tr><tr><td>Magdis+2012a (BzK)</td><td>0.51–1.60</td><td>10.5–11.0</td><td>9</td><td>Magdis et al. (2012a)</td></tr><tr><td>Tacconi+2018 (PHIBSS 1&2)</td><td>0.50–2.49</td><td>9.8–11.6</td><td>148</td><td>Tacconi et al. (2018)</td></tr><tr><td>Kaasinen+2019</td><td>1.78–2.93</td><td>10.6–11.7</td><td>10</td><td>Kaasinen et al. (2019)</td></tr><tr><td>Magdis+2017 (LBG)</td><td>2.8–2.9</td><td>11.28–11.38</td><td>1</td><td>Magdis et al. (2017)</td></tr><tr><td>Magdis+2012b (LBG)</td><td>2.9–3.2</td><td>11.0–11.3</td><td>1</td><td>Magdis et al. (2012b)</td></tr><tr><td>Tan+2014 (GN20)</td><td>4.05–4.06</td><td>10.6–11.0</td><td>3</td><td>Tan et al. (2014)</td></tr></tbody></table>

**Notes.**

<sup>a</sup> We only include sources with >3σ detections.

<sup>b</sup> The largest photometric redshift in the A<sup>3</sup>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.

<sup>c</sup> We 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).

<sup>d</sup> We 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).

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 >0.5 solar metallicity.

As A<sup>3</sup>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 δ<sub>GDR</sub> 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 > 1 galaxies. Here we adopt a slightly modified form of

$$12 + \log_{10}(\text{O/H}) = \begin{cases} a, & \text{if } \log_{10}(M_{\star}/M_{\odot}) \ge b(z), \\ a - 0.087 \times (\log_{10}(M_{\star}/M_{\odot}) - b(z))^2, & \text{else.} \end{cases} \\ \text{where } a = 8.74 \text{ and } b(z) = 10.4 + 4.46 \times \log_{10}(1 + z) \\ - 1.78 \times (\log_{10}(1 + z))^2. \tag{6}$$

The modification (under the $\log_{10}(M_{\star}/M_{\odot}) \ge b(z)$ condition) prevents a drop in metallicity for very massive galaxies at $z < 1$ (see Figure 18).

Finally, we consider that our A<sup>3</sup>COSMOS high-redshift galaxies are molecular-rich (same as assumed by T18 at $z > 0.4$); that is, the molecular-to-total-gas ratio $f_{\text{mol frac}}$ is unity. In this way, we obtain the dust-SED-based $M_{\text{mol gas}}$ by multiplying $M_{\text{dust}}$ with the mass–metallicity-derived $\delta_{\text{GDR}}$ and ignore the contribution from atomic gas. Hereafter we refer to this method as the “$\delta_{\text{GDR},z}$” method.

We caution that, as discussed in Appendix A.3, observations of local galaxies actually indicate that $f_{\text{mol frac}}$ is usually below 50% even for a galaxy with $M_{\star} \sim 1 \times 10^{11} M_{\odot}$. Applying an actual $f_{\text{mol frac}}$, for example, based on the Krumholz et al. (2009) correlation with stellar mass or metallicity, will lead to a lower $M_{\text{mol gas}}$. Based on our next comparison of $M_{\text{mol 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_{\text{mol frac}}$ to future work.

## 3.2. Molecular Gas Mass from RJ-tail Dust Continuum

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).
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Overview of popular molecular gas mass ($M_{mol gas}$) calibration methods for high-redshift (i.e., $z > 1$) galaxies sorted into three categories: (a) (sub)millimeter emission line observations, (b) galaxy SED-fitted dust mass, and (c) RJ-tail dust continuum. See description in Section 3. The corresponding sections in this paper are labeled in the flow chart. The references for the $\alpha_{850, mol}$ conversion factors are Scoville et al. (2017; S17) and Hughes et al. (2017; H17); and for $\alpha_{160-500, mol}$, Groves et al. (2015; G15).</figcaption></figure>

Scoville et al. (2014) found a constant ratio between dust continuum luminosity and $M_{total gas}$,

$$ \alpha_{rj,tot} \equiv \frac{L_{\nu_{ij}}/(\text{erg s}^{-1} \text{Hz}^{-1})}{M_{\text{total gas}}/(M_{\odot})}, \qquad (7) $$

where they calibrated $\alpha_{rj,tot}$ to be $1.0 \pm 0.23 \times 10^{20}$ (erg s$^{-1}$ 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_{total gas}/\nu L_{\nu,500} = 28.5$ for near-solar-metallicity galaxies, corresponding to $\alpha_{500,tot} = 2.2 \times 10^{20}$ (erg s$^{-1}$ 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_{rj,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_{rj,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_{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_{mol gas}$,

$$ \alpha_{rj, mol} \equiv \frac{L_{\nu_{ij}}/(\text{erg s}^{-1} \text{Hz}^{-1})}{M_{\text{mol gas}}/(M_{\odot})}, \qquad (8) $$

to be $6.7 \pm 1.7 \times 10^{19}$ (erg s$^{-1}$ 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_{rj, mol}$ has a dependency

on $L_{\nu_{ij}}$ or $L'_{CO}$ (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^9$ to $10^{12} M_{\odot}$ and CO $J = 1 \rightarrow 0$ line luminosity $L'_{CO(1-0)}$ from $10^7$ to $10^{12} \text{ K km s}^{-1} \text{ pc}^2$, galaxies have significantly varied metallicity, $f_{\text{mol frac}}$, and $\delta_{\text{GDR}}$. A simple explanation for the variations is that $\alpha_{rj, mol}$ scales with $f_{\text{mol frac}}^{-1}$ and $\delta_{\text{GDR}}^{-1}$, which both relate to metallicity (see Appendix A.3 and A.2, respectively).

Since the literature on the calibration of $\alpha_{rj, 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, S17}$,” “$\alpha_{850, H17}$”, and “$\alpha_{160-500, G15}$” method, respectively. For the “$\alpha_{160-500, G15}$” method, we use the calibration factors for the $\log_{10}(M_{\star}/M_{\odot}) > 9$ galaxies in Table 5 of Groves et al. (2015), and we assume that our $A^3\text{COSMOS}$ galaxies have negligible atomic gas contribution. These works directly calibrate the ratio between $L_{\nu_{ij}}$ and $M_{\text{mol gas}}$ (and “$\alpha_{850, H17}$” and “$\alpha_{160-500, G15}$” include a luminosity dependency), so the need for a calibration of the underlying $f_{\text{mol frac}}$ and $\delta_{\text{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\text{m}$ and $\lesssim 160 \mu\text{m}$, respectively, for galaxies at $z \gtrsim 4$. In Figure 2, we show the longest rest-frame wavelengths of the available ALMA data ($\lambda_{\text{rest}}$) for each galaxy in our sample.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Upper panel: redshift versus the longest rest-frame wavelengths ($\lambda_{rest}$) of the available ALMA data for each galaxy in our sample. Color indicates whether the source has a spectroscopic redshift (spec-z; red) or a photometric redshift (photo-z; blue). The dark-gray and light-gray shading represent $\lambda_{rest} < 160$ µm and $160 \le \lambda_{rest} < 250$ µm, respectively. Lower panel: histograms of the longest ALMA $\lambda_{rest}$. Color indicates the same subsample as above. The three labels are the percentages of sources with $\lambda_{rest} < 160$ µm, $160 \le \lambda_{rest} < 250$ µm, and $\lambda_{rest} \ge 250$ µm, respectively.</figcaption></figure>

Some 85% of our sources have $\lambda_{rest} \ge 250$ µm, while the rest only probe shorter-wavelength dust continua.

In order to apply the $\alpha_{rf,mol}$ conversion from dust continuum to molecular gas mass, a “band conversion”<fnref n="19" /> 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,S17}$ and $\alpha_{850,H17}$ methods, and either 160, 250, 350, or 500 µm for applying the $\alpha_{160-500,G15}$ 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_{rest} < 60$ µ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 m,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 > 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_{rest} \le 160$ µ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 $S/N \le 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 × (1 + z) µm and the observed wavelength,

$$ \Gamma^{\text{SED}} \equiv S_{\nu/850 \times (1+z)\mu m}^{\text{SED}} / S_{\nu_{\text{obs}}}^{\text{SED}}, \quad (9) $$

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

$$ L_{\nu/850 \mu m,rest} = 4\pi d_L^2 \times S_{\nu_{\text{obs}}}^{\text{ALMA}} \times \Gamma^{\text{SED}}/(1+z). \quad (10) $$

In principle, we can also directly take the SED-predicted rest-frame 850 µm flux density $S_{\nu/850 \times (1+z)\mu m}^{\text{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 m,rest}$ by $\alpha_{850}$ derived from S17 and Hughes et al. (2017) to obtain the “$\alpha_{850,S17}$” and “$\alpha_{850,H17}$” molecular gas masses. In the “$\alpha_{160-500,G15}$” 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.

### 3.3. Comparing Gas Mass Calibrations

In Figure 3 we compare the molecular gas masses estimated from the above-mentioned “$\delta_{GDR,Z}$,” “$\alpha_{850,S17}$,” “$\alpha_{850,H17}$,” and “$\alpha_{160-500,G15}$” methods. As shown in the bottom row of the figure, “$\delta_{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_{rest} > 250$ µ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).

In the top panel, the “$\alpha_{850,S17}$” and “$\alpha_{850,H17}$” methods agree within 0.1 dex for sources with long-wavelength coverage (but up to about 0.3 dex for sources lacking >160 µm data). However, a systematic offset of about 0.1 dex exists, which is likely because “$\alpha_{850,S17}$” uses a single conversion factor while “$\alpha_{850,H17}$” 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,S17}$” and “$\alpha_{850,H17}$” methods are compared to those using the “$\alpha_{160-500,G15}$” method. The “$\alpha_{160-500,G15}$” method leads to 0.25 dex lower molecular gas masses than “$\alpha_{850,S17}$,” or 0.15 dex lower than “$\alpha_{850,H17}$” for the majority of sources. A small number of sources with poor long-wavelength coverage,

<footnote n="19">¹⁹ This means first applying the K-correction (Humason et al. 1956; Oke & Sandage 1968) to the best-fit SED, then interpolating or extrapolating to certain calibration wavelengths, and then scaling the observed ALMA flux accordingly; see details afterwards.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Comparisons between four methods of gas mass calibration based on dust SED or RJ-tail continuum as presented in Section 3 and labeled at top right. We divide the sources into three categories based on their longest available rest-frame ALMA wavelength (denoted as $\lambda_{\text{rest}}$): $\lambda_{\text{rest}} \le 160~\mu\text{m}$ (red), $160 < \lambda_{\text{rest}} \le 250~\mu\text{m}$ (orange), and $\lambda_{\text{rest}} > 250~\mu\text{m}$ (blue), which also correspond to the three same color-shaded areas in Figure 2, respectively. In each panel, the dashed line is a one-to-one line, and the parallel dotted lines indicate a factor of two variation. The embedded histogram plotted in each Y-versus-X scatter plot is the normalized distribution of $\log_{10}(Y/X)$. See discussion in Section 3.3.</figcaption></figure>

however, have smaller differences. This is probably due to the smaller $M_{\star} > 10^9~M_{\odot}$ sample in Groves et al. (2015) and the intrinsic variation in $L_{\nu_{\text{RJ,rest}}}/M_{\text{mol gas}}$.

To summarize, we find that the gas mass calibrations are as follows: $M_{\text{mol gas}(\delta_{\text{GDR,Z}}} \lesssim M_{\text{mol gas}(\alpha_{160-500,G15}) < M_{\text{mol gas}(\alpha_{850,H17}) < M_{\text{mol gas}(\alpha_{850,S17})}$. 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,H17}$” method to other methods as well as recent observations (Bertemes et al. 2018; Saintonge et al. 2018), we choose the “$\alpha_{850,H17}$” method as our final gas mass calculation for the $A^3\text{COSMOS}$ galaxies. We also tested our full analysis with other gas mass calibrations in Section 5, finding that our results are not obviously altered.

## 4. Galaxy Molecular Gas and Star Formation Properties

After the calculation of molecular gas mass for $A^3\text{COSMOS}$ galaxies, we combine them with our complementary galaxy samples listed in Table 1, allowing us to study galaxy molecular gas and star formation scaling relations and gas evolution in the following sections. In total, we have 1663 galaxies with redshift, SFR, stellar mass, and molecular gas mass measurements. All complementary galaxies are selected to have CO detections, and their molecular gas masses are homogenized with metallicity-dependent $\alpha_{\text{CO}}$, as detailed in Section 2.2. Such a combined sample is the largest, most robust individually detected sample so far, yet it still exhibits a certain incompleteness in the parameter space of redshift, stellar mass, and star formation, due to sample selection biases. Therefore,

before analyzing the gas scaling relations and the resulting gas evolution, we first provide detailed inspections to constrain potential sample selection biases.

### 4.1. Sample Distribution across the MS

In Figure 4 we show the specific star formation rate (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^3\text{COSMOS}$ data points are color coded. Whereas our $A^3\text{COSMOS}$ sample primarily populates the $z > 1$ regime, the complementary samples provide coverage at $z < 1$. However, we do notice that the MS is not well sampled at $0.1 \lesssim z \lesssim 0.5$ and $z > 1$. Only the most massive $A^3\text{COSMOS}$ galaxies ($\log_{10}(M_{\star}/M_{\odot}) \gtrsim 11.5$) sample well the MS, while less-massive ones lie above.

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^3\text{COSMOS}$ galaxies, they are slightly less massive (see Table 1), so T18 sources are able to represent the $\log_{10}(M_{\star}/M_{\odot}) \sim 10-11$ MS, while the $A^3\text{COSMOS}$ galaxies are probing the $\log_{10}(M_{\star}/M_{\odot}) \sim 11-12$ MS at $z > 1$.

We also notice that only very low redshift ($z \lesssim 0.03$) galaxies cover the $\log_{10}(M_{\star}/M_{\odot}) \sim 9-10$ parameter space. In this low-stellar-mass range, the metallicity-dependent $\alpha_{\text{CO}}$ might be more uncertain, and so are the estimated molecular
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Specific star formation rate (sSFR ≡ SFR/M<sub>⋆</sub>) versus redshift distribution of our galaxies (see Table 1) normalized by the star-forming MS SFR at each redshift. The Y = 0 dashed line means exactly on the Speagle et al. (2014) MS. Data points are the same in the left and right panels, except that the redshift axis is logarithmic in the left panel to better illustrate the complementary samples (with A<sup>3</sup>COSMOS data points in gray) and is linear in the right panel to illustrate our A<sup>3</sup>COSMOS sample (with all complementary samples shown in gray).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Molecular gas depletion time $\tau_{\text{depl}}$ versus various galaxy properties (from left to right): redshift, cosmic age, $\Delta\text{MS}$, $M_\star$, SFR, and sSFR, respectively (see Section 4). For each distribution, we show two vertically adjacent panels: the upper panel shows the scatter plot, and the lower panel provides density contours. In the scatter plots, data points are color coded by $M_\star$ in the first two top panels and by redshift (in $\log_{10}(1+z)$ scale, but tick labels are $z$) in the remaining panels. In the contour plots, we divide the sample into three main subsamples: orange solid contours represent the A$^3$COSMOS galaxies ($0.5 \lesssim z \lesssim 6$), green dashed contours the PHIBSS 1&2 galaxies ($0.5 \lesssim z \lesssim 2$), and gray dotted contours are all other local/low-$z$ galaxies.</figcaption></figure>

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.

## 4.2. Correlating Molecular Gas Fraction and Depletion Time to Galaxy Stellar Mass and Star Formation Properties

Here we study the scaling relations for the two most important molecular gas properties: molecular gas depletion time, $\tau_{\text{depl}} \equiv M_{\text{mol gas}}/\text{SFR}$, and molecular gas to stellar mass ratio, $\mu_{\text{molgas}} \equiv M_{\text{mol gas}}/M_\star$. In Figures 5 and 6, we show their distributions versus other galaxy properties: redshift, cosmic age, offset to the MS ($\Delta\text{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$^3$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_{\text{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 $\Delta\text{MS}$, SFR, or sSFR, as shown in the corresponding scatter plots in Figure 5.

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_{\text{depl}}$. There is even a break or turnover feature in the cosmic age and SFR versus $\tau_{\text{depl}}$ panels, which is likely caused by the selection bias at $z > 3$, where we only cover the most massive galaxies ($\log_{10} M_\star \sim 11-12$). In the intermediate-redshift range ($1 \lesssim z \lesssim 3$), the A$^3$COSMOS
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Molecular gas to stellar mass ratio $\mu_{\text{molgas}}$ as a function of various galaxy properties. See Figure 5 caption for the description of data points and contours.</figcaption></figure>

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_{\text{molgas}}$. Compared to gas depletion time distributions, $\mu_{\text{molgas}}$ has a nearly three orders of magnitude variation from local to high redshift, and even at $z > 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_{\text{molgas}}$ also exhibits a strong dependency on $\Delta\text{MS}$, but local galaxies are systematically offset from the high-redshift ones by nearly one dex below.

As shown in the $M_{*}-\mu_{\text{molgas}}$ panel, local galaxies and high-redshift galaxies seem to follow different distributions: local galaxies have similar $\mu_{\text{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_{\text{molgas}}$ and SFR and sSFR. Local, PHIBSS 1&2 intermediate-redshift and $A^3\text{COSMOS}$ higher-redshift galaxies form a contiguous distribution from SFR $\sim 0.1$ to $>1000\ M_{\odot}\ \text{yr}^{-1}$ and sSFR $\sim 0.01$ to $\sim 50\ \text{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_{*}$ term in both axes, is equivalent to the $M_{\text{mol gas}}$ versus SFR correlation, that is, the Kennicutt–Schmidt law (Schmidt 1959; Kennicutt 1998a). A log–log space linear fitting gives a slope of $\sim 0.7$ with a scatter of 0.3 dex, consistent with Kennicutt (1998a), as well as the slope of $\sim 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

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Three-dimensional view of the molecular gas to stellar mass ratio $\mu_{\text{molgas}}$ as a function of lookback time and stellar mass. Data points are colored by their offsets from the MS ($\Delta\text{MS}$); purple is below the MS, green is on the MS, and yellow is above the MS. (This interactive figure is only available in the online journal, where the reader can rotate/zoom/shift the view to see how the data points are distributed in 3D.)</figcaption></figure>

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.

## 4.3. Composite View of Galaxy Gas Fraction and MS Evolution

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_{\text{mol 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
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Galaxies' redshift and specific star formation rate (sSFR ≡ SFR/M<sub>*</sub>) and gas fraction ($f_{mol\ gas}$; Equation (1)). We divide our sample into four stellar mass ranges as labeled at the bottom-right corner. The colored curves represent the evolution of the galaxy MS sSFR (using Speagle et al. 2014 MS) for each subsample. Data points of each subsample are binned by redshift, and in each bin (i.e., each box in the figure) the sSFR versus $f_{mol\ gas}$ distribution is shown in a symmetric vertical histogram style (i.e., a “spindle” diagram). The horizontal axis indicates $f_{mol\ gas}$ in the following way: $f_{mol\ gas}$ = 0% when the width of the histogram is zero (i.e., a thin line at the bin center), and $f_{mol\ gas}$ = 100% when the width equals the box width. The vertical position in each box shares the same Y-axis of the whole figure, that is, corresponds to sSFR. See discussion in Section 4.3. The complete figure set (five images) is available in the online journal, where each stellar mass bin is shown individually for better readability. (The complete figure set (five images) is available.)</figcaption></figure>

high-redshift bins, a higher sSFR corresponds to a larger horizontal span inside each box, which means a higher gas fraction.

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_{mol\ 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 ($\Delta$MS > 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.)

### 4.4. Linking to the Galaxy Star Formation Law

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 SFR = $A \times M_{mol\ 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_{mol\ gas}$–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^3$COSMOS and compiled sample and present how current models are fitting the data.

Figure 9 shows the correlation between SFR and $M_{mol\ gas}$ for all galaxies in this work. Data points are color coded by $\Delta$MS, and the $A^3$COSMOS and complementary samples are distinguished by different symbols (circle and cross, respectively). For $A^3$COSMOS galaxies with large SFR (~100–3000 $M_{\odot}$ yr$^{-1}$) and $M_{mol\ gas}$ (~5 × $10^{10}$–5 × $10^{11} M_{\odot}$), a higher $\Delta$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 ($\Delta$MS < 0.5) exhibit a small deviation of 0.12 dex (median). Considering that the $A^3$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.

However, the starbursts lying significantly above the MS ($\Delta$MS ~ 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 $\Delta$MS ~ 0.5 are able to reach the “starburst SF law,” and the positive correlation between $\Delta$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, sSFR/sSFR<sub>MS</sub>, represents the offset to the MS, with sSFR<sub>MS</sub> computed following Speagle et al. (2014). The Y-axis, $\tau_{mol\ gas}/\tau_{mol\ gas, MS\ SF\ law}$, represents the offset to the “star-forming SF law,” where $\tau_{mol\ gas, MS\ SF\ law} \equiv M_{mol\ gas, SF\ law}/SFR_{MS} = \alpha \times SFR^3/SFR_{MS}$, and the $\alpha$ and $\beta$ coefficients are
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Galaxy $M_{\text{mol gas}}$–SFR relation, that is, star formation law. Data points are color coded by galaxies’ offsets from the $M_*$–SFR MS ($\Delta$MS; using Speagle et al. 2014 MS). Solid circles are $\text{A}^3\text{COSMOS}$ galaxies, and crosses represent all other galaxies from the complementary samples. The blue and red lines represent, respectively, the star formation laws from Sargent et al. (2014) for normal star-forming galaxies and extreme starbursts with a factor of $\sim 15$ offset. See discussion in Section 4.4.</figcaption></figure>

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., sSFR $\lesssim 3 \times \text{sSFR}_{\text{MS}}$) obey the SF law with a Galactic-like $\alpha_{\text{CO}}$, while starbursts with sSFR above the MS (e.g., sSFR $\gtrsim 3 \times \text{sSFR}_{\text{MS}}$) are shifted toward the “starburst SF law,” and they also have a much lower $\alpha_{\text{CO}}$. The shift in the SF- law plane happens most rapidly when the sSFR increases from $\sim 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.

The data are more complicated than what the 2-SFM model predicts. Galaxies in the lowest-mass bin ($\log_{10} M_* < 9.8$) are below the model-predicted curve, while in the mid-stellar-mass bins ($\log_{10} M_* \sim 9.8–11.2$), some galaxies are above it. The turnover is likely seen in the two higher-mass bins ($\log_{10} M_* > 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_{\text{CO}}$ (see Section 2.2); therefore, $\alpha_{\text{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 ($\log_{10} M_* = 9.8$) galaxies at $z \sim 2$, resulting in an upper limit on their gas depletion times of $<0.8$ Gyr, or $\tau_{\text{mol gas}}/\tau_{\text{MS,SF law}} < 0.6$ (assuming a Galactic $\alpha_{\text{CO}}$), in agreement with our findings. If the difference between data and model is truly significant, then it implies that low-mass ($\log_{10} M_* \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_{\text{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_{\text{CO}}$ alters the $\tau_{\text{mol gas}}/\tau_{\text{MS,SF law}}$ ratio from close to one to 0.2 for a starburst galaxy with $\Delta$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 ($\log_{10} M_* \sim 10–11$) sample has a strong bias toward a higher $\Delta$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 ($\log_{10} M_* \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 ($\log_{10} M_* \sim 10–12$) MS galaxies obey the SFGs’ SF law; (b) from the MS to $\Delta$MS $\gtrsim 1$, galaxies start to deviate from the SFGs’ SF law toward the starbursts’ SF law, with a rapid change at $\Delta$MS $\sim 0.4–0.6$ dex roughly in agreement with the 2-SFM prediction; (c) low-mass ($\log_{10} M_* \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.

## 4.5. Intermediate Summary on the Advantage and Caveats of This Sample

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 $\text{A}^3\text{COSMOS}$. The distribution of our sample in the (z, $M_*$, $\Delta$MS) high-dimensional space is roughly contiguous from (z, $\log_{10} M_*$, $\Delta$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
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Correlation between galaxies' offsets from the star formation law (as indicated by the Y-axis $\tau_{mol\ gas}/\tau_{mol\ gas, MS SF\ law}$) and from the galaxy MS (as indicated by the X-axis sSFR/sSFR<sub>MS</sub>) for four stellar mass bins as labeled. We adopt the Sargent et al. (2014) SF law and Speagle et al. (2014) MS (adopting a different MS from Sargent et al. (2014) will only result in a small shift of -0.16, -0.24, -0.17, and -0.05 dex in the X direction (and smaller offsets in the Y direction) for the four panels, respectively). Symbol shapes are the same as in Figure 9 and are color coded by log<sub>10</sub> SFR. The black line is the predicted median trend from the 2-SFM framework by Sargent et al. (2014). The gray shaded area indicates a 0.3 dex scatter in both X and Y.</figcaption></figure>

of the parameter space hold the key to further improving such studies.

5. Characterizing the Galaxy Molecular Gas Scaling Relation via Functional Fitting

Through our previous discussion of galaxy properties (molecular gas to stellar mass ratio $\mu_{molgas}$, molecular gas depletion time $\tau_{depl}$, stellar mass $M_*$, SFR, and redshift $z$ or the corresponding cosmic age $t_{cosmic\ age}$), we can already see the complexity inherent in their scaling relations. In this section, we provide high-dimensional functional fittings to simultaneously quantify the underlying dependencies of $\mu_{molgas}$ and $\tau_{depl}$ on $z$ (or $t_{cosmic\ age}$), $M_*$, and SFR.

We propose a new functional form that accounts for the different behaviors of galaxies that are due to their stellar masses seen in the previous figures:

$$ \log_{10} \tau_{\text{depl}} = (a + ak \times \log_{10}(M_*/10^{10})) \times \Delta\text{MS} \\ + b \times \log_{10}(M_*/10^{10}) \\ + (c + ck \times \log_{10}(M_*/10^{10})) \times t_{\text{cosmic age}} \\ + d $$

(best-fit: $a = -0.5724$, $b = -0.5174$,

$c = -0.002997$, $d = +0.02964$,

$ak = +0.1120$ and $ck = +0.0568$)

$$ \log_{10} \mu_{\text{molgas}} = (a + ak \times \log_{10}(M_*/10^{10})) \times \Delta\text{MS} \\ + b \times \log_{10}(M_*/10^{10}) \\ + (c + ck \times \log_{10}(M_*/10^{10})) \times t_{\text{cosmic age}} \\ + d $$

(best-fit $a = +0.4195$, $b = -0.6906$,

$c = -0.1543$, $d = +0.9339$,

$ak = +0.1195$ and $ck = +0.0320$)

(11)

where $t_{cosmic\ age}$ and $M_*$ are in units of Gyr and $M_\odot$, respectively.<fnref n="20" />

Here we adopt the MS function from the #49 fitting of Table 7 of Speagle et al. (2014), which is their preferred fit (see their abstract and Table 9). This functional form also uses cosmic age as a free parameter (as in our function).<fnref n="21" /> We compared various MS in the literature (e.g., Sargent et al. 2014; Whitaker et al. 2014; Béthermin et al. 2015; Lee et al. 2015; Schreiber et al. 2015; Tomczak et al. 2016; Pearson et al. 2018), finding that the Speagle et al. (2014) cosmic-age MS function provides the most reasonable fitting (see Appendix A.5). It can be rewritten in the same style as the above functions as follows:

$$ \log_{10} \text{SFR}_{\text{MS}} = b \times \log_{10}(M_*/10^{10}) \\ + (c + ck \times \log_{10}(M_*/10^{10})) \times t_{\text{cosmic age}} \\ + d \\ \text{where } b = +0.84, c = -0.15, d = +1.89, \text{ and} \\ ck = -0.026. \qquad (12) $$

We fit our new functional form to the combined sample in this work and refitted both the S17 and T18 functional forms, as described in Equations (4) and (5), respectively. We use the Python packages `pymc3` and `scipy.optimization` for the fitting.<fnref n="22" /> The former package performs Markov chain Monte Carlo (MCMC) fitting to calculate the probability distribution of the fitting, while the latter one performs least-squared minimization to find the best fit. The two algorithms agree very well, and the former one provides better uncertainty estimates for the fitted parameters.

We list our best-fit parameters in Table 2 as well as in Equation (11). The parameters fitted by S17 and T18 for their own functional forms are also provided in Table 2 for comparison.

<footnote n="20">²⁰ See Appendix C for the probability distributions of the fitted coefficients. We also provide a Python package named `a3cosmos-gas-evolution` for the calculation with our functions: https://ascl.net/code/v/2377 (Liu & A3COSMOS Team 2019).</footnote>

<footnote n="21">²¹ We remind the reader that in Speagle et al. (2014) the functional fitting with redshift ($\log_{10}(1+z)$) in their Table 8 is not suitable for use at high redshifts, such as $z \gtrsim 4$. Those fits are very different from the functions with cosmic age. For example, with the same #49 fitting set, the two MS functions agree only at $z \lesssim 1.5$, while the difference can be 0.5 dex at $z \sim 4.3$ (with the MS function with $\log_{10}(1+z)$ being higher).</footnote>

<footnote n="22">²² `pymc3` documentation: https://docs.pymc.io/; `scipy.optimization.curve_fit`: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimization.curve_fit.html.</footnote>
<table><caption>**Table 2**

Best-fit Coefficients for the Molecular Gas to Stellar Mass Ratio and Molecular Gas Depletion Time Functions</caption><thead><tr><th>log<sub>10</sub>(&mu;<sub>molgas</sub>)<br>(&equiv;log<sub>10</sub>(M<sub>mol gas</sub>/M<sub>*</sub>))</th><th>a (+ak)<br>(&Delta;MS)</th><th>b<br>(log<sub>10</sub>M<sub>*</sub>,10)<sup>a</sup></th><th>c (+ck)<br>(log<sub>10</sub>(1 + z) or t)</th><th>d<br>(norm.)</th></tr></thead><tbody><tr><td>This work, Equation (11)</td><td>+0.4195 + 0.1195 &times; log<sub>10</sub>M<sub>*</sub>,10</td><td>-0.6907</td><td>(-0.1543 + 0.0320 &times; log<sub>10</sub>M<sub>*</sub>,10) &times; t</td><td>+0.9339</td></tr><tr><td>This work, Equation (5)</td><td>+0.54</td><td>-0.39</td><td>-4.42 &times; (log<sub>10</sub>(1 + z) - 0.58)<sup>2</sup></td><td>+0.309</td></tr><tr><td>T18 best-fit, Equation (5)</td><td>+0.53</td><td>-0.35</td><td>-3.62 &times; (log<sub>10</sub>(1 + z) - 0.66)<sup>2</sup></td><td>+0.365</td></tr><tr><td>T18 best-fit, Equation (5)</td><td>+0.53</td><td>-0.35</td><td>-3.62 &times; (log<sub>10</sub>(1 + z) - 0.66)<sup>2</sup></td><td>+0.120</td></tr><tr><td>This work, Equation (4)</td><td>+0.57</td><td>-0.28</td><td>+2.27 &times; log<sub>10</sub>(1 + z)</td><td>-1.106</td></tr><tr><td>S17 best-fit, Equation (4)</td><td>+0.32</td><td>-0.70</td><td>+1.84 &times; log<sub>10</sub>(1 + z)</td><td>-0.149</td></tr><tr><td>log<sub>10</sub>(&tau;<sub>depl</sub>/(Gyr))<br>(&equiv;log<sub>10</sub>(M<sub>mol gas</sub>/SFR))</td><td></td><td></td><td></td><td></td></tr><tr><td>This work, Equation (11)</td><td>-0.5724 + 0.1120 &times; log<sub>10</sub>M<sub>*</sub>,10</td><td>-0.5174</td><td>(-0.0030 + 0.0568 &times; log<sub>10</sub>M<sub>*</sub>,10) &times; t</td><td>+0.0269</td></tr><tr><td>This work, Equation (5)</td><td>-0.40</td><td>+0.08</td><td>-0.91 &times; log<sub>10</sub>(1 + z)</td><td>+0.022</td></tr><tr><td>T18 best-fit, Equation (5)</td><td>-0.44</td><td>+0.09</td><td>-0.62 &times; log<sub>10</sub>(1 + z)</td><td>+0.027</td></tr><tr><td>This work, Equation (4)</td><td>-0.39</td><td>+0.08</td><td>-0.91 &times; log<sub>10</sub>(1 + z)</td><td>+0.021</td></tr><tr><td>S17 best-fit, Equation (4)</td><td>-0.70</td><td>-0.01</td><td>-1.04 &times; log<sub>10</sub>(1 + z)</td><td>+0.509</td></tr></tbody></table>

**Note.**

<sup>a</sup> log<sub>10</sub> M<sub>*</sub>,10 represents log<sub>10</sub>(M<sub>*</sub>/10<sup>10</sup> M<sub>&odot;</sub>).

Our refitting of the T18 function agrees with their original fitting: only the redshift coefficient is slightly changed by about 10%, which implies that the two fittings are consistent (<30%) at z < 2 and slightly discrepant at z ~ 3–5, where our fitting predicts about 30%–50% lower gas fractions and shorter depletion times, mainly driven by the new data coverage from this work (their data only covers z ~ 0–3).

For our new function, the fitted dependencies of gas fraction and depletion time on each parameter are presented in Figure 11. We show in each panel the best-fit function curve and the data points with a rescaling to remove the dependencies on other parameters than the current one presented by the X-axis of that panel. This rescaling uses our best-fit result. For example, for the rescaling in the first panel, the gas-to-stellar mass ratio $\mu_{molgas}$ of a galaxy with $\Delta MS = 1$ and $log_{10} M_{*} = 10.5$ will be shifted along the Y-axis by $-0.4195 \times \Delta MS$ dex (see the a coefficient in the $\mu_{molgas}$ formula in Equation (11)), bringing it down to the MS galaxy level. In this way, each panel only indicates the dependency of our function fitting on the parameter presented by the X-axis.

We also show the original best fits of T18 and S17 (to their own functional forms, i.e., Equations (5) and (4), respectively) in Figure 11. In comparison, our new functional form has a log-linear dependency on cosmic age, so $\mu_{molgas}$ and $\tau_{depl}$ almost flatten beyond redshift ~4 for the same stellar mass and $\Delta MS$ galaxies. The T18 best-fit function predicts a drop at z ~ 4 in $\mu_{molgas}$, while the S17 best-fit function predicts $\mu_{molgas}$ to continue increasing with redshift. In the $\Delta MS$ and $M_{*}$ panels, we also see certain differences, but our function is between the T18 and S17 ones.

The current fitting still has some minor caveats. For example, in the redshift panel, the limited number of z > 4.5 data points is mostly below our best-fit function. But as we discussed in Section 3.2.1, the band conversion for the rest-frame RJ-tail dust continuum has a large uncertainty when there are no long-wavelength data, which is the case for z > 4 galaxies. The test in Appendix B shows that our SED fitting tends to underestimate their true RJ-tail dust continuum by a factor of 2–6.

In addition, we tested the stability of our fitting for subsamples of galaxies: (a) only z < 4 data, and (b) without z > 1 CO (which are mostly from the PHIBSS 1&2 surveys from T18), that is, only using $A^3\text{COSMOS}$ dust-based data at z > 1. The tests show that the z > 4 data and z > 1 CO data do not statistically bias our fitting results, likely because their numbers are not large enough compared to the full sample. The $\chi^2$ information of these test fittings is listed in Table 3, which shows that our proposed functional form in Equation (11) gives statistically better fits to the data in this work than both the T18 (Equation (5)) and S17 (Equation (4)) functions, and that the T18 one is better than the S17 one. This is likely because our function (Equation (11)) has one more free parameter than the T18 function, which further has one more degree of freedom than the S17 function.

Moreover, we have run our full fitting process for other gas mass calibration methods. We find that using the S17 gas mass calibration, which slightly overpredicts gas masses compared to the H17 calibration, leads to ~11% changes in the coefficients in Equation (11) and results in a 11% shallower $\mu_{molgas}$ versus stellar mass (negative) dependency. This in turn increases the prediction of $\mu_{molgas}$ for MS, $log_{10} M_{*} \sim 10.5$ galaxies by a small amount of about 20% at z ~ 6. On the other hand, using the $\delta_{GDR,Z}$ gas mass calibration, which tends to underestimate the gas masses, we find the coefficients change by ~40%. This results in a ~40% steeper $\mu_{molgas}$ versus stellar mass (negative) dependency, and consequently lower $\mu_{molgas}$ for MS, $log_{10} M_{*} \sim 10.5$ galaxies at z ~ 6 by about 50%. The scatters in the diagnostic plots similar to those in Figure 11 are also larger by about 0.06 dex (e.g., the scatters around the best-fit lines in Figure 11 are about 0.28–0.30 dex with the $\alpha_{850,H17}$ calibration, while they are 0.33–0.36 dex with the $\delta_{GDR,Z}$ calibration). We also note that the slope of the $\mu_{molgas}$ versus $\Delta MS$ correlation is much less obviously affected (a ~ 0.39–0.41), and the fits close to z ~ 0 are not obviously affected, due to the large number of local/low-z galaxies in our sample with CO-based gas masses. These tests show that the choice of gas mass calibration method is not significantly altering our result,
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Characterizing molecular gas depletion time $\tau_{\text{depl}}$ (upper panels) and molecular gas to stellar mass ratio $\mu_{\text{molgas}}$ (lower panels) in the functional form of Equation (11). From left to right, we show $\tau_{\text{depl}}$ versus redshift, $\tau_{\text{cosmic age}}$, $\Delta\text{MS}$, and $M_\star$, respectively. Data points in each panel are rescaled using the best-fit function to remove the dependency on other parameters and leave only the correlation with the current X-axis parameter (with coefficient(s) labeled at the bottom of each panel). Orange data points are from $\Lambda^3\text{COSMOS}$, while green ones are from the PHIBSS 1&2 surveys (T18), and gray ones are from the literature as listed in Table 1 and at the top. We distinguish these samples by different symbols in order to better reveal outliers and sample biases against each parameter after removing other parameter dependencies. Our best-fit function is shown as the orange solid line in each panel, while the functions from T18 (see Equation (5)) and S17 (see Equation (4)) are shown as green dashed and pink dotted lines, respectively.</figcaption></figure>

by at most a factor of two at $z \sim 6$ and $\log_{10} M_\star \sim 10.5$, and less at lower redshifts.

Finally, we emphasize that, despite the fact that nearly all gas masses for our high-redshift ($z > 2.5$) galaxies are dust-based and similarly those for local/low-redshift ($z < 1$) galaxies are CO-based, we verified that our results are not significantly biased. Specifically we have excluded all CO-based galaxies and all dust-based galaxies from our fitting. We find that the slopes of both the $\mu_{\text{molgas}}$ versus $\Delta\text{MS}$ and $\mu_{\text{molgas}}$ versus $\log_{10} M_\star$ correlations are quite stable within 20%. However, the trend of the time evolution is significantly driven by the lack of constraining data at either low or high redshift: excluding all CO-based galaxies leads to a factor of 10 higher gas fraction at $z \sim 0$, because there is basically no constraint at $z < 1$, while excluding all dust-based galaxies gives a factor of 2 higher gas fraction at $z \sim 4-6$, due to little constraint at $z > 3$. Taking together the good consistency when fitting the non-redshift-dependent correlations and the good agreement between CO- and dust-based gas mass calibrations in the literature (see beginning of Section 3), it is not only very reasonable to combine the CO- and dust-based samples but also necessary to achieve sensible results.

# 6. Predictions from Fitted Gas Evolution Functions

## 6.1. Evolution of Molecular Gas Depletion Time

Here we discuss the cosmic evolution of the molecular gas depletion time $\tau_{\text{depl}}$ as predicted by our best-fit function (Equation (11)) for galaxies in bins of stellar mass and $M_\star$ offset. In Figure 12, we bin all our 1663 galaxies into $4 \times 3$ panels, with the $\log_{10} M_\star$ bin center ranging from 9.0 to 12.0

(bin width 1.0) and the $\Delta\text{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.

From the figure, we can see that our function behaves differently than the other two functions. The evolution of $\tau_{\text{depl}}$ exhibits a much stronger dependency on stellar mass in our function. For very massive ($\log_{10} M_\star \sim 12.0$) galaxies, our function predicts a factor of about 20 increase in $\tau_{\text{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 ($\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 $\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_{\text{depl}}$ ranges from 0.1 to 1.0 Gyr (with SFR $\sim$ 20–30 to $0.04 M_\odot \text{yr}^{-1}$, i.e., from high to low $\Delta\text{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 $\Delta\text{MS}$ with time. Thus, for example, the flat $\tau_{\text{depl}}$ versus redshift trend for less-massive ($\log_{10} M_\star \sim 10.0$)
<table><caption>**Table 3**

Statistics of Fitting the Functions of Molecular Gas to Stellar Mass Ratio ($\mu_{molgas} \equiv M_{mol gas}/M_{\star}$) with Subsamples</caption><thead><tr><th rowspan="2">Fitting Function<br>(For μ<sub>molgas</sub>)</th><th colspan="3">All Data Points</th><th colspan="3">Without z > 4 Data</th><th colspan="3">Without z > 1 CO</th></tr><tr><th>N</th><th>χ<sup>2</sup></th><th>χ<sup>2</sup><sub>redu.</sub></th><th>N</th><th>χ<sup>2</sup></th><th>χ<sup>2</sup><sub>redu.</sub></th><th>N</th><th>χ<sup>2</sup></th><th>χ<sup>2</sup><sub>redu.</sub></th></tr></thead><tbody><tr><td>Equation (11) (This work)</td><td>1663</td><td>1426.98</td><td>0.86</td><td>1617</td><td>1390.16</td><td>0.86</td><td>1554</td><td>1349.56</td><td>0.87</td></tr><tr><td>Equation (5) (This work)</td><td>1663</td><td>1444.83</td><td>0.87</td><td>1617</td><td>1404.74</td><td>0.87</td><td>1554</td><td>1369.32</td><td>0.88</td></tr><tr><td>Equation (5) (T18s fit)</td><td>1663</td><td>1503.94</td><td>0.91</td><td>1617</td><td>1478.11</td><td>0.92</td><td>1554</td><td>1427.91</td><td>0.92</td></tr><tr><td>Equation (4) (This work)</td><td>1663</td><td>2086.55</td><td>1.26</td><td>1617</td><td>1863.42</td><td>1.16</td><td>1554</td><td>1921.61</td><td>1.24</td></tr><tr><td>Equation (4) (S17s fit)</td><td>1663</td><td>10230.35</td><td>6.17</td><td>1617</td><td>10029.8</td><td>6.22</td><td>1554</td><td>10141.08</td><td>6.54</td></tr></tbody></table>

galaxies seen in the middle columns of the figure does not imply a constant $\tau_{depl}$ for an individual galaxy across its evolution history: its stellar mass growth will move it into a higher stellar mass $\tau_{depl}$ evolution track. In the $\log_{10} M_{\star} \sim 10.0$ and $\Delta\text{MS} \sim 0.5$ bin, our function does not fit well the $z \sim 3-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_{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_{depl} = -0.36 \log_{10} \text{sSFR} - 0.14 \log_{10}(\Sigma_{\star}) + 5.87^{23}$. This translates into $\tau_{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 ~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} < M_{\star}/M_{\odot} < 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).

## 6.2. Evolution of Molecular Gas Fraction

Similar to the previous section, we show in Figure 13 the binned view of the evolution of $\mu_{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_{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 $\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 $\log_{10} M_{\star} \sim 9$ galaxies).

The dependency of $\mu_{molgas}$ or $f_{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, $\propto \text{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-9.0} M_{\odot}$ (see also Cao et al. 2017; Saintonge et al. 2017) with a nonlinear behavior. In their sample, $\mu_{molgas}$ is about $0.08-0.3$ for $\log_{10} M_{\star} \sim 9-10$ galaxies, then decreases to about $0.02-0.1$ for $\log_{10} M_{\star} \sim 10-11$ galaxies, which is slightly below this work at $z \sim 0$ and is likely caused by their use of a constant $\alpha_{CO}$, while the true $\alpha_{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 > 5-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 ($\log_{10} M_{\star} \lesssim 10$) galaxy samples.

## 7. Implication for the Cosmic Evolution of Cold Molecular Gas Density

In this section, we study the implication of our cold molecular gas fraction function (Equation (11)) for the cosmic molecular gas mass density evolution. This requires us to know

<footnote n="23">23 Note that we are using their sSFR function instead of the $\Sigma_{\text{sSFR}}$ function in their abstract.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Evolution of the molecular gas depletion time $\tau_{\text{depl}} \equiv M_{\text{mol gas}}/\text{SFR}$ (in units of Gyr) with redshift in $4 \times 3$ bins of $\Delta\text{MS}$ and $\log_{10} M_*$. From left to right, $\log_{10} M_*$ increases from 9.0 to 12.0 with a step of 1.0 and bin width of 1.0, and from bottom to top $\Delta\text{MS}$ increases from -0.5 to +0.5 with a step of 0.5 and bin width of 0.5. We show the evolution function Equation (11) from this work and those from T18 and S17 (i.e., their best fits to Equations (5) and (4), respectively) in each panel (see the labels at the bottom). These functions are calculated with the mean $\Delta\text{MS}$ and $\log_{10} M_*$ of the subsample data available within each bin. Blue rectangles represent the mean($\tau_{\text{depl}}$) $\pm$ 1$\sigma$ ranges of all galaxies from A$^3$COSMOS and the literature in bins of redshift in each panel. We caution that this figure does not show the quality of data fitting because data still have variations in $\Delta\text{MS}$ and $\log_{10} M_*$ even within each panel. See Figure 11 for the fitting quality, and see Section 6.1 for the discussion of this figure.</figcaption></figure>

(a) the number density of SFGs at each redshift, (b) their stellar mass distribution at each redshift, and (c) their SFRs or $\Delta$MS.

The number density and stellar mass distribution evolution of SFGs has been reasonably well measured through SFG SMF studies. We discuss them in detail in Section 7.1.

Then, by either simply assuming that all SFGs are MS galaxies (Section 7.2), or more realistically adopting the aforementioned 2-SFM galaxy model (Sargent et al. 2014; Béthermin et al. 2017) as we do in Section 7.3, we obtain an SFR for each galaxy corresponding to its stellar mass and redshift. With the SFR and $\Delta$MS, the stellar mass is further converted to gas mass by applying our gas fraction function. Finally, by integrating over all SFGs, we obtain the cosmic molecular gas mass density at each redshift as presented in Sections 7.2 and 7.3.

Such a method is also used by Maeda et al. (2017), who fitted molecular gas fraction versus stellar mass correlations at two redshift bins ($z \sim 0$ and $z = 1-1.5$), and then integrated the cosmic molecular gas mass density using SMFs. Other earlier works (Sargent et al. 2013; see also Carilli & Walter 2013) instead fitted a molecular gas mass versus SFR correlation (i.e., SF law; independent of redshift and stellar mass) to infer gas mass and integrate over SMFs to obtain the cosmic molecular gas mass density.

## 7.1. Adopting the SMFs

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-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-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.

## 7.2. Integrating Cosmic Molecular Gas Mass Density

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,
<figure><img src="image_15.png" /><figcaption>Figure 13. Analogous to Figure 12, but for the evolution of the molecular gas to stellar mass ratio $\mu_{molgas}$. The parameters for each panel are the same as in Figure 12. See Figure 12 caption for the description and Section 6.2 for the discussion of this figure.</figcaption></figure>

stellar mass, and SMF in each stellar mass bin at each redshift:

$$ \rho_{molgas}(z) = \sum_{M_{\star} \text{ bins}} \Phi_{\text{SMF}}(z, M_{\star}) \times M_{\star} \times \mu_{\text{molgas}}(z, M_{\star}, \Delta\text{MS}). \quad (13) $$

In Figure 14 we present the integrated cosmic cold molecular gas mass density versus redshift, using three different gas fraction functions $\mu_{molgas}(z, M_{\star}, \Delta\text{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.

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_{molgas}(z, M_{\star}, \Delta\text{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.

### 7.3. Alternative Method: Cosmic Molecular Gas Mass Density with Mock Galaxy Models

The drawback of the SMF $\times \mu_{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_{molgas}$ function before summing them up within each redshift bin.

The “SIDES” simulation (Simulated Infrared Dusty Extra-galactic Sky$^{24}$; Béthermin et al. 2017) generated 1,489,629 mock galaxies within a 2 deg$^2$ 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 < 1.5$; Davidzon et al. 2017 at $1.5 < z < 4$; and Grazian et al. 2015 at $z > 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 $\log_{10} M_{\star} > 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_{molgas}$. We sample the redshift range from 0 to 15 with 500 bins (i.e., bin size $\sim 0.00253$ in $\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 > 0.5$), our curve coincidentally agrees with the semianalytic model (SAM) simulation by Popping et al. (2019). Other simulations, Obreschkow & Rawlings (2009)

<footnote>$^{24}$ Available at http://cesam.lam.fr/sides.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Cosmic evolution of the cold molecular gas mass density. Results from high-z CO blind deep-field studies from Decarli et al. (2016, 2019) and Riechers et al. (2019) are shown as green, blue, and red boxes, with the X-sides (Y-sides) indicating the observed redshift range (5th and 95th percentiles). The $z \sim 6$ arrow is an upper limit from Riechers et al. (2019). The orange solid, green long-dashed, and pink short-dashed lines are the SMF-integrated molecular gas mass density based on the SMFs presented in Section 7.1 (see also Appendix D; integrated down to $M_* = 10^{9.0} M_\odot$) and the gas fraction function from this work (Equation (11)), T18, and S17, respectively. The black dashed-dotted line is from the semianalytic model (SAM) simulation of Popping et al. (2019; based on Popping et al. 2014).</figcaption></figure>

and Lagos et al. (2011), can be seen in Figure 5 of Riechers et al. (2019): at $z \sim 5$, the Popping et al. (2019) simulation exhibits a 0.2 dex lower $\rho_{\text{mol gas}}$ than that of Lagos et al. (2011), and the Obreschkow & Rawlings (2009) $\rho_{\text{mol gas}}$ is 0.1 dex lower than Lagos et al. (2011), while the three are reversed at $z \sim 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.

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_{\text{molgas}}$ (see Figure 13).<fnref n="25" /> The T18 function results in fully (marginally) consistent cold molecular gas mass densities as our function at $z \lesssim 1$ ($z \sim 2-3$); however, it predicts 0.4–0.9 dex lower values at $z > 4$. This is mainly driven by the downturn of their $\mu_{\text{molgas}}$ function at $z \gtrsim 4$ (as mentioned in Section 5) and probably also affected by their systematically lower $\mu_{\text{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 \text{ Mpc}^{-3}$ on the molecular gas mass densities at redshift $z \sim 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_{\text{molgas}}$ in Equation (11)) is moving toward a coherent picture.

## 8. Summary

In this work, we present a comprehensive analysis of galaxy molecular gas scaling relations and their evolution using a robust ALMA-detected galaxy catalog from our Paper I (A³COSMOS). Each galaxy in the catalog has a redshift, stellar mass, SFR, and dust mass from far-infrared SED fitting including the ALMA data and rich multiwavelength data from the literature (see Paper I for details). We compared four

<footnote n="25">²⁵ We caution that S17 used a different MS function than this work and T18. Our test in Appendix A.5 shows that their MS can explain half of the discrepancy seen in Figures 14 and 15. The other major contributor to the discrepancy is the functional form. As shown in Figures 13, their gas fraction’s functional form is too high at both low and high redshift ($z < 1$ and $z > 4$) and for less-massive ($\log_{10} M_* < 10$) galaxies. While integrating all galaxies to compute the cosmic gas mass density, such a difference in the functional forms causes a large discrepancy.</footnote>
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Analogous to Figure 14, now using the “SIDES” mock galaxy catalog based on the 2-SFM galaxy model (Béthermin et al. 2017) to derive the molecular gas mass density. See Section 7.3 for details. See Figure 14 caption for symbols, curves, and labels. (The curves in this figure are also provided in the Python a3cosmos-gas-evolution package; Liu & A3COSMOS Team 2019.)</figcaption></figure>

methods of molecular gas mass calibration using SED-fitted dust mass or RJ-tail dust continuum (Section 3.3), from which we determine that the RJ-dust continuum method (with the Hughes et al. 2017 luminosity-dependent calibration) better infers the gas mass. Meanwhile, we also comprehensively discuss several related topics in the gas mass calibration, that is, $\alpha_{CO}$, $\delta_{GDR}$, molecular-to-atomic fraction, and metallicity, and their biases to this work in Appendices A.1–A.4.

Due to the sample inhomogeneity, higher-redshift (e.g., $z > 4$) galaxies do not always have RJ-tail wavelength coverage. Thus, we investigated the effect of band conversion with MAGPHYS high-z SED fitting for galaxies whose longest-wavelength ALMA data do not cover RJ-tail wavelengths. We found that it potentially results in a factor of 2–6 underestimation of gas mass at $z > 4$ (see Section 3.2.1 and Appendix B).

We combine our $A^3\text{COSMOS}$ sample with 20 complementary samples in the literature from local to high redshift (see Table 1) to study the scaling relations and cosmic evolution of molecular gas depletion time $\tau_{\text{depl}}$ and molecular gas to stellar mass ratio $\mu_{\text{molgas}}$. We parameterize the $\tau_{\text{depl}}$ and $\mu_{\text{molgas}}$ as functions of galaxy’s cosmic age, stellar mass, and SFR. We tested both Tacconi et al. (2018, T18) and Scoville et al. (2017, S17) functions (shown in Equations (5) and (4) respectively), and we meanwhile propose a new functional form of Equation (11) that accounts for the galaxies’ different evolution paths driven by their stellar masses. Then, by applying the gas fraction scaling relation to galaxies’ SMFs and integrating over all stellar masses, we obtain the evolution of cosmic cold molecular gas mass density, which is in a

coherent picture with the known CSFRD evolution and the semianalytic modeling of galaxies in the cosmological simulations (Figures 14 and 15 respectively).

Furthermore, we emphasize the following points:

1. The distributions of our sample’s redshifts, stellar masses, and SFRs are consistent with previous studies where they overlap in the parameter space (e.g., see contours in Figures 5 and 6). Given our total sample of 1663 galaxies, we see that the composite sample selection is biased to strong starbursts with $\Delta\text{MS} \sim 0.5–1.5$ and $\log_{10}M_\star \sim 10–11$ at $z \sim 0.08–1.0$ (see Figure 4), and biased to the most massive galaxies with $\Delta\text{MS} \sim 0.0$ and $\log_{10}M_\star \sim 12$ at $z > 3$ (see Figures 7 and 8). In particular, at $z > 4$ the dust continuum observations are mainly probing rest-frame wavelengths shorter than 250 µm, for which the SED-fitting-extrapolated RJ-tail flux might be underpredicted. However, they do not statistically affect our functional fitting because of their low number.

2. The parameterizations of $\tau_{\text{depl}}$ and $\mu_{\text{molgas}}$ with the functions in this work, T18, and S17 are roughly consistent where the data are commonly sampled in the parameter space, that is, $z \sim 1–3$, $\Delta\text{MS} > 0$, and $\log_{10}M_\star > 10.5$ (see Figures 12 and 13). They differ significantly for low-mass and MS or below-MS galaxies, which, however, could not be verified with the current data set. The chi-squared statistics for these parameterizations show that our new functional form and the T18 are similarly good and are better than the S17 functional form, which has one (two) fewer free parameter(s) than the T18 one (ours). We
emphasize that our new functional form implicitly leads to a “downsizing” in galaxy evolution and probably a “mass-quenching” effect. Although further data are needed to verify these effects (see Sections 6.1 and 6.2 as well as Figures 12 and 13), the results are promising for building a most comprehensive picture of gas evolution.

3. The integration of the galaxies' SMF with the application of the gas fraction scaling relation involves many assumptions. Noticeable differences are found between the simpler assumption that all SFGs exactly follow the MS (Figure 14) and the more realistic 2-SFM galaxy modeling (Figure 15), which accounts for the starburst/ MS dichotomy and uses different SMFs than in this work (Appendix D). The realistic galaxy modeling has a better agreement with semianalytic models (Popping et al. 2014, 2019). Among the three functional forms discussed in this work, only our new functional form (Equation (11)) of the gas fraction scaling relation could achieve such a high consistency.

4. Compared to CO blind deep-field surveys, our analytically derived cold molecular gas mass densities agree within their upper boundary. This is understandable as the current CO surveys usually could not sample well enough the faint end of the CO line luminosity function, so the integration of CO luminosity functions is usually down to only $log_{10}(L'_{CO}/(\text{K km s}^{-1}\text{pc}^2)) \sim 9.5$ (to avoid extrapolating the faint end; see, e.g., Decarli et al. 2019; Riechers et al. 2019).

5. Finally, our large, robust data set strongly supports a coherent picture of the evolution of galaxies' gas, stellar mass, and SFR, which can be parameterized by the MS functions (e.g., Speagle et al. 2014; Leslie et al. 2019; Appendix A.5), SMFs (e.g., Davidzon et al. 2017; Appendix D), and gas scaling functions (Equation (11)). The integration of SMF times the MS function (over stellar mass at each redshift) gives the CSFRD, and the integration of SMF times the molecular gas fraction function (over stellar mass at each redshift) results in the cosmic molecular gas mass density. Integrating the CSFRD curve (across cosmic time) further leads to the cosmic stellar mass density growth curve, which in turn is consistent with the integration of SMFs across cosmic time.

D.L., E.S., and P.L. acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 694343). B.M., E.F.J.A., and F.B. acknowledge support of the Collaborative Research Center 956, subprojects A1 and C4, funded by the Deutsche Forschungsgemeinschaft (DFG), project ID 184018867. S.L. acknowledges funding from Deutsche Forschungsgemeinschaft (DFG) Grant SCH 536/9-1. D.R. acknowledges support from the National Science Foundation under grant Nos. AST-1614213 and AST-1910107 and from the Alexander von Humboldt Foundation through a Humboldt Research Fellowship for Experienced Researchers. G.E.M. acknowledges support from the Villum Fonden research grant 13160 "Gas to stars, stars to dust: Tracing star formation across cosmic time," the Cosmic Dawn Center of Excellence funded by the Danish National Research Foundation, and the ERC Consolidator Grant funding scheme (project ConTExt, grant No. 648179). Y.G.'s research is supported by National Key Basic Research and Development Program of

China (grant No. 2017YFA0402704), National Natural Science Foundation of China (grant Nos. 11861131007, 11420101002), and Chinese Academy of Sciences Key Research Program of Frontier Sciences (grant No. QYZDJSSW-SLH008). We thank the anonymous referee for constructive comments. We thank Matthieu Béthermin for a very helpful discussion on galaxy modeling.

This paper makes use of the following ALMA data: ADS/JAO. ALMA#2011.0.00064.S, ADS/JAO.ALMA#2011.0.00097.S, ADS/JAO.ALMA#2011.0.00539.S, ADS/JAO.ALMA#2011.0.00742.S, ADS/JAO.ALMA#2012.1.00076.S, ADS/JAO.ALMA#2012.1.00323.S, ADS/JAO.ALMA#2012.1.00523.S, ADS/JAO.ALMA#2012.1.00536.S, ADS/JAO.ALMA#2012.1.00919.S, ADS/JAO.ALMA#2012.1.00952.S, ADS/JAO.ALMA#2012.1.00978.S, ADS/JAO.ALMA#2013.1.00034.S, ADS/JAO.ALMA#2013.1.00092.S, ADS/JAO.ALMA#2013.1.00118.S, ADS/JAO.ALMA#2013.1.00151.S, ADS/JAO.ALMA#2013.1.00171.S, ADS/JAO.ALMA#2013.1.00208.S, ADS/JAO.ALMA#2013.1.00276.S, ADS/JAO.ALMA#2013.1.00668.S, ADS/JAO.ALMA#2013.1.00815.S, ADS/JAO.ALMA#2013.1.00884.S, ADS/JAO.ALMA#2013.1.00914.S, ADS/JAO.ALMA#2013.1.01258.S, ADS/JAO.ALMA#2013.1.01292.S, ADS/JAO.ALMA#2015.1.00026.S, ADS/JAO.ALMA#2015.1.00055.S, ADS/JAO.ALMA#2015.1.00122.S, ADS/JAO.ALMA#2015.1.00137.S, ADS/JAO.ALMA#2015.1.00260.S, ADS/JAO.ALMA#2015.1.00299.S, ADS/JAO.ALMA#2015.1.00379.S, ADS/JAO.ALMA#2015.1.00388.S, ADS/JAO.ALMA#2015.1.00540.S, ADS/JAO.ALMA#2015.1.00568.S, ADS/JAO.ALMA#2015.1.00664.S, ADS/JAO.ALMA#2015.1.00704.S, ADS/JAO.ALMA#2015.1.00853.S, ADS/JAO.ALMA#2015.1.00861.S, ADS/JAO.ALMA#2015.1.00862.S, ADS/JAO.ALMA#2015.1.00928.S, ADS/JAO.ALMA#2015.1.01074.S, ADS/JAO.ALMA#2015.1.01105.S, ADS/JAO.ALMA#2015.1.01111.S, ADS/JAO.ALMA#2015.1.01171.S, ADS/JAO.ALMA#2015.1.01212.S, ADS/JAO.ALMA#2015.1.01495.S, ADS/JAO.ALMA#2015.1.01590.S, ADS/JAO.ALMA#2015.A.00026.S, ADS/JAO.ALMA#2016.1.00478.S, ADS/JAO.ALMA#2016.1.00624.S, ADS/JAO.ALMA#2016.1.00735.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ.

Facility: ALMA.

## Appendix A

### Empirical Galaxy Scaling Relations Used in This Work

Scaling relations describe how galaxy properties correlate with each other and are important for understanding galaxy populations and their evolution over cosmic time. As this work studies the molecular gas evolution in galaxies, the four scaling relations below are relevant and sometimes needed for our analysis. Calibrations of these correlations are widely studied in the literature, but their validity for different types of galaxies (i.e., different z, $M_*$, and SFR) is rarely studied. Here we compare a number of empirical calibrations and discuss their biases. This comparison guides our choice of the most suitable correlations to use in the analysis described in the main body of the paper.
<figure><img src="image_18.png" /></figure>

**Figure 16.** Left panel: metallicity dependence of the CO-to-H₂ conversion factor α<sub>CO</sub>. Symbols and lines are from Wilson (1995), Genzel et al. (2015), Bolatto et al. (2013), Accurso et al. (2017), and T18 as labeled. The horizontal blue-, green-, and orange-shaded regions correspond to α<sub>CO</sub> = 6.5, 4.3, and 0.8 M<sub>☉</sub> (K km s<sup>-1</sup> pc<sup>2</sup>)<sup>-1</sup>, representing Galactic GMCs, the inner Galactic disk, and ULIRGs, respectively. See Appendix A.1 for more details. The vertical dashed line indicates solar metallicity (12 + log<sub>10</sub>(O/H)<sub>☉</sub> = 8.69; Asplund et al. 2009). Right panel: metallicity dependence of the gas-to-dust mass ratio δ<sub>GDR</sub>. Symbols and lines are from Leroy et al. (2011), Magdis et al. (2012a), and Rémy-Ruyer et al. (2014) as labeled. See Appendix A.2 for more details. See Table 1 for references of the data points in both panels.

<figure><img src="image_19.png" /></figure>

### A.1. CO-to-H₂ Conversion Factor (α<sub>CO</sub>) versus Metallicity

The CO-to-H₂ conversion factor, α<sub>CO</sub>, 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<sub>☉</sub> (K km s<sup>-1</sup> pc<sup>2</sup>)<sup>-1</sup> (X<sub>CO</sub> = 2.1 × 10<sup>20</sup> cm<sup>-2</sup>(K km s<sup>-1</sup> pc<sup>2</sup>)<sup>-1</sup>; Solomon et al. 1987; Solomon & Barrett 1991; Solomon & Vanden Bout 2005), or 6.5 M<sub>☉</sub> (K km s<sup>-1</sup> pc<sup>2</sup>)<sup>-1</sup> (X<sub>CO</sub> = 3.0 × 10<sup>20</sup> cm<sup>-2</sup>(K km s<sup>-1</sup> pc<sup>2</sup>)<sup>-1</sup>) when including heavy elements, which are mostly helium, while it is as low as 0.8 M<sub>☉</sub> (K km s<sup>-1</sup> pc<sup>2</sup>)<sup>-1</sup> (X<sub>CO</sub> = 2.0 × 10<sup>20</sup> cm<sup>-2</sup>(K km s<sup>-1</sup> pc<sup>2</sup>)<sup>-1</sup>; ±0.5 dex) in local ULIRGs (Solomon et al. 1997; Downes & Solomon 1998; Solomon & Vanden Bout 2005). The calibration of α<sub>CO</sub> 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 α<sub>CO</sub>-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.

### A.2. Gas-to-dust Ratio (δ<sub>GDR</sub>) versus Metallicity

The gas-to-dust mass ratio, δ<sub>GDR</sub> ≡ M<sub>total gas</sub>/M<sub>dust</sub>, 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, δ<sub>GDR</sub> 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 δ<sub>GDR</sub> can be applied for molecular gas mass estimation for high-redshift galaxies.

The calibration of δ<sub>GDR</sub> 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 δ<sub>GDR</sub> 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 δ<sub>GDR</sub> is around 100 for galaxies with solar and supersolar metallicity, while it increases nonlinearly toward lower metallicity, reaching over 1000 in extremely metal-poor (<10% Z<sub>☉</sub>) 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 δ<sub>GDR</sub> is calibrated with total gas mass instead of molecular gas mass, a molecular-to-total gas mass ratio, f<sub>mol frac</sub> ≡ M<sub>mol gas</sub>/M<sub>total gas</sub>, needs to be considered. It is discussed in the next section (Appendix A.3).

### A.3. Molecular Hydrogen Fraction (f<sub>mol frac</sub>) versus Metallicity

The molecular hydrogen fraction, f<sub>mol frac</sub> ≡ M<sub>mol gas</sub>/M<sub>total gas</sub>, is the ratio between molecular gas and molecular+atomic gas. In the following, we use μ<sub>mol frac</sub> ≡ M<sub>mol gas</sub>/M<sub>atomic gas</sub> for the molecular-to-atomic gas mass ratio.

Note that f<sub>mol frac</sub> (or μ<sub>mol frac</sub>) 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<sub>mol frac</sub> and 12 + log<sub>10</sub>(O/H) and M<sub>*</sub> 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.

In Figure 17, we show μ<sub>mol frac</sub> versus metallicity and stellar mass with a large compilation of 524 galaxies from the literature (see labels and figure caption). All galaxies have M<sub>mol gas</sub> from CO observations, M<sub>atomic gas</sub> from H I observations, 12 + log<sub>10</sub>(O/H) from optical spectroscopy, and M<sub>*</sub> from multiwavelength optical/near-infrared data. The data points exhibit a large scatter in both panels, which is probably
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Molecular-to-atomic gas mass ratio $\mu_{\text{mol frac}} \equiv M_{\text{mol gas}}/M_{\text{atomic gas}}$ versus gas-phase metallicity ($12 + \log_{10}(\text{O/H})$; left panel; using the Pettini & Pagel 2004, PP04 calibration) and stellar mass (right panel; using a Chabrier 2003 initial mass function). Data points are compiled from the literature: see Table 1 for the references for the DGS, HRS, Stripe82, and KINGFISH surveys; in addition, we used the atomic gas mass $M_{\text{H}\text{I}}$ from Catinella et al. (2018) for the xCOLD GASS survey (Saintonge et al. 2017). Theoretical models from Popping et al. (2014, their Equation (8)), Davé et al. (2016, their Equations (1) and (2)), and Krumholz et al. (2009, their Equation (2)) are overlaid as colored lines.</figcaption></figure>

caused by the uncertainties in $M_{\text{mol gas}}$, $M_{\text{atomic gas}}$, and metallicity. The metallicity-dependent CO-to-H$_2$ 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_{\text{atomic 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_* > 10^{10} M_\odot$, molecular gas nearly dominates the total gas mass ($\langle f_{\text{mol frac}} \rangle > 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_{\text{mol 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_* > 10^{10} M_\odot$ (e.g., T18).

### A.4. Mass-Metallicity Relation (MZR)

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_{\text{CO}}$, $\delta_{\text{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 + \log_{10}(\text{O/H})$) parameterized as a function of $M_*$ 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 $\log_{10} M_* > 11.2$; and (b) it predicts the lowest metallicity for starburst galaxies at $z > 1$. And the caveat of their Equation (4) is the nonphysical extrapolation for (a) MS galaxies at all redshift with $\log_{10} M_* > 11.2$, and (b) $z \gtrsim 1$ MS galaxies with $\log_{10} M_* < 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 $\log_{10} M_* - 0.32 \times \log_{10}$ SFR < 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 $\log_{10} M_* \sim 10.2-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 $\log_{10} M_* \gtrsim 11$ at local (or $\log_{10} M_* \gtrsim 11.5-12$ at $z \gtrsim 1$). This motivates our modification of this equation as described in Equation (6).
<figure><img src="image_21.png" /><figcaption>**Figure 18.** Comparison of different metallicity relations in the literature for main-sequence galaxies in six redshift bins (top two rows) and for starburst galaxies in three redshift bins (bottom row). Metallicity relations are equations that compute metallicity from stellar mass, SFR, and/or redshift. The curves in each panel are the metallicity computed with each equation in bins of stellar mass, but at a given redshift and SFR as labeled in the top right corner. In each panel, equations are as labeled at the bottom right and described in Appendix A.4, while data points are taken from observations in the literature as labeled in the bottom left. The horizontal black dashed line indicates solar metallicity.</figcaption></figure>

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 > 1$.

6. Kewley & Ellison (2008) use PP04 O3N2 and N2 MZRs as listed in their Table 2.<fnref n="26" /> 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 $\log_{10} M_* < 10.5$, and like Mannucci et al. (2010) Equation (2), it also has a nonphysical drop with increasing stellar mass when $\log_{10} M_* > 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 < 2-3$) but predicts significantly subsolar metallicity at $z \gtrsim 5$ even for starbursts.

Figure 18 shows that most of these formulae are consistent (within 0.2 dex) only for MS galaxies at $z \lesssim 2$ and with $\log_{10} M_* \sim 10-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_{\text{GDR}}$ when assuming the $\delta_{\text{GDR}}$-metallicity relation in Figure 16. Also note that the prescription for deriving $12 + \log_{10}(\text{O/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

<footnote n="26">²⁶ Note that their equation in the Table 2 caption should be $y = a + bx + cx^2 + dx^3$.</footnote>
<figure><img src="image_22.png" /><figcaption>**Figure 19.** Comparison of galaxy star-forming main-sequence functions in the literature as labeled at the top: Speagle et al. (2014), Sargent et al. (2014), Whitaker et al. (2014), Béthermin et al. (2015), Schreiber et al. (2015), Lee et al. (2015), Tomczak et al. (2016), Pearson et al. (2018), S17, and Leslie et al. (2019). We show the MS evolution curves predicted by these functions at four representative stellar masses (from left to right): $log_{10}(M_*/M_\odot) = 9.0, 10.0, 11.0, \text{ and } 12.0$. See discussion in Appendix A.5.</figcaption></figure>

formulae do not agree well at the low-mass regime ($M_* < 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.

### A.5. Stellar Mass–SFR MS

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 $\text{sSFR}_\text{MS}$ at $z > 5$, while the Béthermin et al. (2015) MS predicts a factor of 2 higher $\text{sSFR}_\text{MS}$ than average at $z > 5$ and is in general higher for $\log_{10} M_* > 10^{11}$ galaxies. The Pearson et al. (2018) MS is a factor of <2 lower than others at $z \sim 2–3$ and in general lower for $\log_{10} M_* < 10^{10}$ galaxies. These MS calibrations have large scatter in the $\log_{10} M_* < 10^9$ and $\log_{10} M_* \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.

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 ~0.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\text{Jy/beam}$ at a spatial resolution of 0.75" 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_{\text{molgas}}$ evolution curves using the Leslie et al. (2019) MS to those using the Speagle et al. (2014) MS.

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 ($\log_{10} M_* \lesssim 9.2, z \sim 2–6$), while being indistinguishable from using the Speagle et al. 2014 MS for galaxies with $\log_{10} M_* \gtrsim 10–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 $\Delta\text{MS}$ for low-mass galaxies at $z \sim 0$, altering the slope of $\mu_{\text{molgas}}$ versus $\Delta\text{MS}$ to be shallower (by a small change of 0.04 in the coefficient a in Equation (11)), meanwhile steepening the slope of $\mu_{\text{molgas}}$ versus $\log_{10} M_*$ (by a change of 0.14 in the coefficient b in Equation (11)). Thus it results in a <2× 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 $\log_{10} M_* \gtrsim 10–11$ galaxies at all redshifts ($z \sim 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 > 3–4$ (and exhibits a large excess compared to others), here we show their curve and use their MS to compute $\Delta\text{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 $\Delta\text{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).
<figure><img src="image_23.png" /><figcaption>**Figure 20.** Similar to Figure 13, comparing the evolution of gas fraction $\mu_{molgas}$ curves obtained by using the Leslie et al. (2019) MS (green long-dashed line), S17 MS (pink short-dashed line), and the Speagle et al. (2014) MS (orange solid line; same as in Figure 13). See description of the binning scheme and data boxes in the captions for Figures 12 and 13, and see discussion in Appendix A.5.</figcaption></figure>

# Appendix B

## Biases in Band Conversion with MAGPHYS SED Fitting

To verify the potential bias in using MAGPHYS SED fitting to predict the rest-frame RJ-tail (i.e., 850 µm) dust continuum, we have done some tests using the multiwavelength data from UV to submillimeter for the JINGLE survey galaxy sample (Saintonge et al. 2018; Smith et al. 2019). JINGLE galaxies are MS SFGs in the local universe ($z < 0.05$) with well-sampled SEDs, including JCMT/SCUBA2 850 µm (Smith et al. 2019), Herschel 70–500 µm (Griffin et al. 2010; Pilbratt et al. 2010; Poglitsch et al. 2010), Spitzer 3.6–24 µm (Werner et al. 2004), Wide-field Infrared Survey Explorer 3.4–22 µm (Wright et al. 2010), VISTA 1.2–2.2 µm (Sutherland et al. 2015), 2MASS J, H and K (Skrutskie et al. 2006), SDSS optical (York et al. 2000; Eisenstein et al. 2011), and Galaxy Evolution Explorer UV (Morrissey et al. 2007).

We run the following MAGPHYS fitting tests: (a) fitting all data points at $\lambda \le 250$ µm, mimicking the $\lambda_{rest} \le 250$ µm cases in Figure 2 in the main text; (b) fitting all $\lambda < 8$ µm photometry data points plus only one $\lambda = 160$ µm data point, mimicking the cases where we have only one ALMA data point for fitting the whole dust SED (see Figure 2).

We compare the SED-predicted 850 µm fluxes from both tests to the true observed 850 µm fluxes in Figure 21 (left and right panels, respectively). As also mentioned in Section 3.2.1, the SED-predicted fluxes tend to be lower than the true fluxes, and the accuracy of predicting the 850 µm flux seems to depend on the S/N of all data points fitted for the dust component SED. When dust SED data points have an S/N ~ 20, the rest-frame 850 µm fluxes tend to be underestimated by ~0.5 dex. If S/N > 30, it seems some galaxies have no bias, while some still are underpredicted. For S/N < 20, the 850 µm fluxes are significantly underestimated by 0.3–0.8 dex.
<figure><img src="image_24.png" /><figcaption>**Figure 21.** Tests of MAGPHYS SED fitting with the JINGLE survey local ($z < 0.05$) galaxy sample (Saintonge et al. 2018) and multiwavelength data (from UV to SCUBA2 850 µm; Smith et al. 2019), which show the biases of MAGPHYS SED fitting when (left) fitting all photometry data points at $\lambda \le 250$ µm while no data point at longer wavelength is used, and (right) fitting photometry data points up to $\lambda < 8$ µm plus only one $\lambda = 160$ µm data point for the dust component. Galaxies with observed SCUBA2 850 µm flux $S/N \ge 3$ ($S/N < 3$) are shown as blue circles with error bars indicating the observational errors (red arrows representing $3\sigma$ upper limits in the observed fluxes). See Appendix B for details.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 22.** Probability distributions of the coefficients in the $\tau_{\text{depl}}$ and $\mu_{\text{molgas}}$ functional forms in Equation (11) as fitted by our MCMC fitting to all our data. See the fitting in Section 5. The left panel is for the $\tau_{\text{depl}}$ function, and the right panel is for the $\mu_{\text{molgas}}$ function.</figcaption></figure>

## Appendix C
MCMC Fitting

Figure 22 shows the probability distributions of the coefficients in Equation (11) obtained from our MCMC fitting in Section 5. In the fitting, we allow the coefficients to vary within a relatively large range of (−10, 10). The probability distribution of each coefficient is shown to include the most probable areas as automatically determined by the Python package corner. All coefficients have a clear peak in their probability distribution with a small width (uncertainty) of about 10%. The $\tau_{\text{depl}}$ function's coefficients seem to have some second peaks that are probably due to the nonuniform, complicated sample biases (see Section 2). The overall constraint of our fitting (to Equation (11) and other fittings) is tight.

## Appendix D
Stellar Mass Function

The galaxies' SMF at each redshift should in principle be consistent with the evolution of the CSFRD. As mentioned in Section 7.1, at a given cosmic time, the integration of the CSFRD over previous cosmic times should be equal to the integration of the SMF at that cosmic time over all stellar masses. Note that the SMF is usually divided into two galaxy types, SFGs and quiescent galaxies (QGs), and the integration of the SMF should be the sum of both SFGs and QGs.

Therefore we adopt SMFs for this work by adjusting known SMFs according to the integration of the CSFRD. For example, at $z < 0.085$, we adopt an SMF with a shape the same as for
<figure><img src="image_26.png" /><figcaption>**Figure 23.** Comparison of SMFs in nine redshift bins. The top nine panels are SMFs of star-forming galaxies (SFG); the middle nine panels are SMFs of quiescent galaxies (QG); and the bottom nine panels are the sum of SFG+QG, that is, total galaxy SMFs. Labels represent the following references: Davidzon et al. (2017), Peng et al. (2010), and Wright et al. (2018). The shape of the SMFs adopted in this work is from Peng et al. (2010) if *z* < 0.08, or the linear interpolation of Davidzon et al. (2017) SMFs at intermediate redshifts, or the shape of the Davidzon et al. (2017) SMFs at *z* = 4 if *z* > 4. The normalization of the SMFs adopted in this work is set to be consistent with the integration of the cosmic SFR density (Madau & Dickinson 2014). See Appendix D for details.</figcaption></figure>

the SMF from Peng et al. (2010), for both SFG and QG types, and with the normalization of SFG+QG adjusted to the CSFRD-integrated total stellar mass at that redshift, while keeping the SFG and QG SMFs' relative normalization the same as in Peng et al. (2010). Similarly, at higher redshifts, we adopt the shape of our SMF from an interpolation of the Davidzon et al. (2017) SMFs, as their SMFs are measured over multiple redshift bins ($0.085 < z < 0.2$ and $z > 4.0$), we adopt their $z = 0.2$ and $z = 4.0$ SMF shapes, respectively. The normalization is also adjusted such that the SFG+QG SMFs' integrated total stellar mass equals the CSFRD-integrated total stellar mass.

Note that during the integration of CSFRD over cosmic time, we have considered the loss of mass due to stellar evolution following Conroy & Wechsler (2009, see their Equation (11)). The choice of the mass-losing timescale can be different; for example, Ilbert et al. (2013) adopt 3 Myr and Behroozi & Silk (2015) adopt 1.4 Myr. Compared to the Conroy & Wechsler (2009) timescale, adopting 3 Myr would lead to a 0.09 dex higher integrated total stellar mass at $z = 0$.

In Figure 23 we compare our adjusted SMFs with the measured SMFs from Davidzon et al. (2017) and Peng et al. (2010), which are recent measurements with the deepest data available at high redshift and in the local universe, respectively. The very good agreement between our SMFs and theirs supports our knowledge of galaxy evolution characterized by SMFs and CSFRD. We also compare the SMFs from Wright et al. (2018), who compiled a large amount of data and performed a function fitting to characterize the evolution of the SMF. We show both their single- and double-Schechter function fitting in Figure 23; however, their SMFs are too high at the massive end, while changing rapidly in shape at $z > 6$. This perhaps shows the difficulty in obtaining a best-fitting function and is the reason that we adopt the CSFRD-adjusted SMFs rather than the function-characterized ones.

ORCID iDs

Daizhong Liu https://orcid.org/0000-0001-9773-7479

E. Schinnerer https://orcid.org/0000-0002-3933-7677

B. Groves https://orcid.org/0000-0002-9768-0246
B. Magnelli https://orcid.org/0000-0002-6777-6490

P. Lang https://orcid.org/0000-0002-5681-3575

S. Leslie https://orcid.org/0000-0002-4826-8642

D. A. Riechers https://orcid.org/0000-0001-9585-1462

G. Popping https://orcid.org/0000-0003-1151-4659

Georgios E. Magdis https://orcid.org/0000-0002-4872-2294

E. Daddi https://orcid.org/0000-0002-3331-9590

M. Sargent https://orcid.org/0000-0003-1033-9684

Yu Gao https://orcid.org/0000-0003-0007-2197

Y. Fudamoto https://orcid.org/0000-0001-7440-8832

P. A. Oesch https://orcid.org/0000-0001-5851-6649

F. Bertoldi https://orcid.org/0000-0002-1707-1775

## References

Accurso, G., Saintonge, A., Catinella, B., et al. 2017, MNRAS, 470, 4750

Alatalo, K., Davis, T. A., Bureau, M., et al. 2013, MNRAS, 432, 1796

Andreani, P., Boselli, A., Ciesla, L., et al. 2018, A&A, 617, A33

Aniano, G., Draine, B. T., Calzetti, D., et al. 2012, ApJ, 756, 138

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

Baldry, I. K., Driver, S. P., Loveday, J., et al. 2012, MNRAS, 421, 621

Battisti, A. J., da Cunha, E., Grasha, K., et al. 2019, ApJ, 882, 61

Bauermeister, A., Blitz, L., Bolatto, A., et al. 2013, ApJ, 763, 64

Behroozi, P. S., & Silk, J. 2015, ApJ, 799, 32

Berta, S., Lutz, D., Genzel, R., Förster-Schreiber, N. M., & Tacconi, L. J. 2016, A&A, 587, A73

Bertemes, C., Wuyts, S., Lutz, D., et al. 2018, MNRAS, 478, 1442

Béthermin, M., Daddi, E., Magdis, G., et al. 2015, A&A, 573, A113

Béthermin, M., Le Floc'h, E., Ilbert, O., et al. 2012, A&A, 542, A58

Béthermin, M., Wu, H.-Y., Lagache, G., et al. 2017, A&A, 607, A89

Betti, S. K., Pope, A., Scoville, N., et al. 2019, ApJ, 874, 53

Blain, A. W., Smail, I., Ivison, R. J., Kneib, J.-P., & Frayer, D. T. 2002, PhR, 369, 111

Bolatto, A. D., Leroy, A. K., Jameson, K., et al. 2011, ApJ, 741, 12

Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207

Boselli, A., Cortese, L., & Boquien, M. 2014a, A&A, 564, A65

Boselli, A., Cortese, L., Boquien, M., et al. 2014b, A&A, 564, A66

Boselli, A., Eales, S., Cortese, L., et al. 2010, PASP, 122, 261

Bourne, N., Dunne, L., Bendo, G. J., et al. 2013, MNRAS, 436, 479

Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151

Cao, T.-W., Wu, H., Du, W., et al. 2017, AJ, 154, 116

Capak, P. L., Carilli, C., Jones, G., et al. 2015, Natur, 522, 455

Cappellari, M., Emsellem, E., Krajnović, D., et al. 2011, MNRAS, 413, 813

Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105

Catinella, B., Saintonge, A., Janowiecki, S., et al. 2018, MNRAS, 476, 875

Chabrier, G. 2003, PASP, 115, 763

Cibinel, A., Daddi, E., Sargent, M. T., et al. 2019, MNRAS, 485, 5631

Cicone, C., Bothwell, M., Wagg, E., et al. 2017, A&A, 604, A53

Ciesla, L., Boquien, M., Boselli, A., et al. 2014, A&A, 565, A128

Ciesla, L., Boselli, A., Smith, M. W. L., et al. 2012, A&A, 543, A161

Combes, F. 2018, A&ARv, 26, 5

Combes, F., García-Burillo, S., Braine, J., et al. 2013, A&A, 550, A41

Conroy, C., & Wechsler, R. H. 2009, ApJ, 696, 620

Coogan, R. T., Sargent, M. T., Daddi, E., et al. 2019, MNRAS, 485, 2092

Cormier, D., Madden, S. C., Lebouteiller, V., et al. 2014, A&A, 564, A121

Cortzen, I., Garrett, J., Magdis, G., et al. 2019, MNRAS, 482, 1618

Cowie, L. L., Songaila, A., Hu, E. M., & Cohen, J. G. 1996, AJ, 112, 839

da Cunha, E., Charlot, S., & Elbaz, D. 2008, MNRAS, 388, 1595

da Cunha, E., Walter, F., Smail, I. R., et al. 2015, ApJ, 806, 110

Daddi, E., Bournaud, F., Walter, F., et al. 2010a, ApJ, 713, 686

Daddi, E., Dannerbauer, H., Elbaz, D., et al. 2008, ApJL, 673, L21

Daddi, E., Dickinson, M., Morrison, G., et al. 2007, ApJ, 670, 156

Daddi, E., Elbaz, D., Walter, F., et al. 2010b, ApJL, 714, L118

Dale, D. A., Aniano, G., Engelbracht, C. W., et al. 2012, ApJ, 745, 95

Dale, D. A., Helou, G., Contursi, A., Silbermann, N. A., & Kolhatkar, S. 2001, ApJ, 549, 215

Davé, R., Thompson, R., & Hopkins, P. F. 2016, MNRAS, 462, 3265

Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70

De Vis, P., Jones, A., Viaene, S., et al. 2019, A&A, 623, A5

Decarli, R., Walter, F., Aravena, M., et al. 2016, ApJ, 833, 69

Decarli, R., Walter, F., Gónzalez-López, J., et al. 2019, ApJ, 882, 138

Denicoló, G., Terlevich, R., & Terlevich, E. 2002, MNRAS, 330, 69

Downes, D., & Solomon, P. M. 1998, ApJ, 507, 615

Draine, B. T. 2003, ARA&A, 41, 241

Draine, B. T., Dale, D. A., Bendo, G., et al. 2007, ApJ, 663, 866

Draine, B. T., & Li, A. 2007, ApJ, 657, 810

Eisenstein, D. J., Weinberg, D. H., Agol, E., et al. 2011, AJ, 142, 72

Elbaz, D., Daddi, E., Le Borgne, D., et al. 2007, A&A, 468, 33

Elbaz, D., Leiton, R., Nagar, N., et al. 2018, A&A, 616, A110

Elmegreen, B. G., Rubio, M., Hunter, D. A., et al. 2013, Natur, 495, 487

Erb, D. K., Shapley, A. E., Pettini, M., et al. 2006, ApJ, 644, 813

Fernández Lorenzo, M., Sulentic, J., Verdes-Montenegro, L., & Argudo-Fernández, M. 2013, MNRAS, 434, 325

Fukui, Y., Hayakawa, T., Inoue, T., et al. 2018, ApJ, 860, 33

Galliano, F., Galametz, M., & Jones, A. P. 2018, ARA&A, 56, 673

Gao, Y., & Solomon, P. M. 2004, ApJ, 606, 271

Genzel, R., Tacconi, L. J., Gracia-Carpio, J., et al. 2010, MNRAS, 407, 2091

Genzel, R., Tacconi, L. J., Lutz, D., et al. 2015, ApJ, 800, 20

Grazian, A., Fontana, A., Santini, P., et al. 2015, A&A, 575, A96

Griffin, M. J., Abergel, A., Abreu, A., et al. 2010, A&A, 518, L3

Groves, B. A., Schinnerer, E., Leroy, A., et al. 2015, ApJ, 799, 96

Hayward, C. C., & Smith, D. J. B. 2015, MNRAS, 446, 1512

Hollenbach, D., & Salpeter, E. E. 1971, ApJ, 163, 155

Huang, M.-L., & Kauffmann, G. 2014, MNRAS, 443, 1329

Huang, M.-L., & Kauffmann, G. 2015, MNRAS, 450, 1375

Hughes, T. M., Cortese, L., Boselli, A., Gavazzi, G., & Davies, J. I. 2013, A&A, 550, A115

Hughes, T. M., Ibar, E., Villanueva, V., et al. 2017, MNRAS, 468, L103

Humason, M. L., Mayall, N. U., & Sandage, A. R. 1956, AJ, 61, 97

Hunt, L. K., De Looze, I., Boquien, M., et al. 2019, A&A, 621, A51

Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55

Jiang, X.-J., Wang, Z., Gu, Q., Wang, J., & Zhang, Z.-Y. 2015, ApJ, 799, 92

Jin, S., Daddi, E., Liu, D., et al. 2018, ApJ, 864, 56

Kaasinen, M., Scoville, N. Z., Walter, F., et al. 2019, ApJ, 880, 15

Kelvin, L. S., Driver, S. P., Robotham, A. S. G., et al. 2014, MNRAS, 444, 1647

Knicicutt, R. C., Calzetti, D., Aniano, G., et al. 2011, PASP, 123, 1347

Knicicutt, R. C., Jr. 1998a, ApJ, 498, 541

Knicicutt, R. C., Jr. 1998b, ARA&A, 36, 189

Kewley, L. J., & Ellison, S. L. 2008, ApJ, 681, 1183

Kirkpatrick, A., Pope, A., Aretxaga, I., et al. 2014, ApJ, 796, 135

Klitsch, A. Péroux, C., Zwaan, M. A., et al. 2019, MNRAS, 490, 1220

Kriek, M., van Dokkum, P. G., Labbé, I., et al. 2009, ApJ, 700, 221

Krumholz, M. R., McKee, C. F., & Tumlinson, J. 2009, ApJ, 699, 850

Lagos, C. D. P., Baugh, C. M., Lacey, C. G., et al. 2011, MNRAS, 418, 1649

Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24

Lee, N., Sanders, D. B., Casey, C. M., et al. 2015, ApJ, 801, 80

Lee, N., Sheth, K., Scott, K. S., et al. 2017, MNRAS, 471, 2124

Leja, J., Carnall, A. C., Johnson, B. D., Conroy, C., & Speagle, J. S. 2019, ApJ, 876, 3

Leroy, A. K., Bolatto, A., Gordon, K., et al. 2011, ApJ, 737, 12

Leroy, A. K., Walter, F., Bigiel, F., et al. 2009, AJ, 137, 4670

Leslie, S. K., Schinnerer, E., Magnelli, B., et al. 2019, ApJ, submitted

Li, A., & Draine, B. T. 2001, ApJ, 554, 778

Liang, L., Feldmann, R., Faucher-Giguère, C.-A., et al. 2018, MNRAS, 478, L83

Liang, L., Feldmann, R., Kereš, D., et al. 2019, MNRAS, 489, 1397

Lisenfeld, U., Alatalo, K., Zucker, C., et al. 2017, A&A, 607, A110

Lisenfeld, U., Isaak, K. G., & Hills, R. 2000, MNRAS, 312, 433

Liu, D. & A3COSMOS Team 2019, a3cosmos-gas-evolution: Galaxy Cold Molecular Gas Evolution Functions, Astrophysics Source Code Library, ascl:1910.003

Liu, D., Daddi, E., Dickinson, M., et al. 2018, ApJ, 853, 172

Liu, D., Gao, Y., Isaak, K., et al. 2015, ApJL, 810, L14

Liu, D., Lang, P., Magnelli, B., et al. 2019, ApJS, 244, 40

Lutz, D. 2014, ARA&A, 52, 373

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Maeda, F., Ohta, K., & Seko, A. 2017, ApJ, 835, 120

Magdis, G. E., Daddi, E., Béthermin, M., et al. 2012a, ApJ, 760, 6

Magdis, G. E., Daddi, E., Elbaz, D., et al. 2011, ApJL, 740, L15

Magdis, G. E., Daddi, E., Sargent, M., et al. 2012b, ApJL, 758, L9

Magdis, G. E., Rigopoulou, D., Daddi, E., et al. 2017, A&A, 603, A93

Magdis, G. E., Rigopoulou, D., Hopwood, R., et al. 2014, ApJ, 796, 63

Magnelli, B., Lutz, D., Saintonge, A., et al. 2014, A&A, 561, A86

Magnelli, B., Saintonge, A., Lutz, D., et al. 2012, A&A, 548, A22

Maiolino, R., Nagao, T., Grazian, A., et al. 2008, A&A, 488, 463

Mannucci, F., Cresci, G., Maiolino, R., Marconi, A., & Gnerucci, A. 2010, MNRAS, 408, 2115
Mannucci, F., Salvaterra, R., & Campisi, M. A. 2011, MNRAS, 414, 1263

Marchesini, D., van Dokkum, P. G., Förster Schreiber, N. M., et al. 2009, ApJ, 701, 1765

McGaugh, S. S., & de Blok, W. J. G. 1997, ApJ, 481, 689

Morrissey, P., Conrow, T., Barlow, T. A., et al. 2007, ApJS, 173, 682

Moustakas, J., Coil, A. L., Aird, J., et al. 2013, ApJ, 767, 50

Moutard, T., Arnouts, S., Ilbert, O., et al. 2016, A&A, 590, A103

Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18

Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJL, 660, L43

Obreschkow, D., & Rawlings, S. 2009, MNRAS, 394, 1857

Oke, J. B., & Sandage, A. 1968, ApJ, 154, 21

Pearson, W. J., Wang, L., Hurley, P. D., et al. 2018, A&A, 615, A146

Peng, Y.-j., Lilly, S. J., Kovač, K., et al. 2010, ApJ, 721, 193

Pettini, M., & Pagel, B. E. J. 2004, MNRAS, 348, L59

Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1

Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2

Popping, G., Pillepich, A., Somerville, R. S., et al. 2019, ApJ, 882, 137

Popping, G., Somerville, R. S., & Trager, S. C. 2014, MNRAS, 442, 2398

Privon, G. C., Narayanan, D., & Davé, R. 2018, ApJ, 867, 102

Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2014, A&A, 563, A31

Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2015, A&A, 582, A121

Riechers, D. A., Pavesi, R., Sharon, C. E., et al. 2019, ApJ, 872, 7

Rodighiero, G., Daddi, E., Baronchelli, I., et al. 2011, ApJL, 739, L40

Rodighiero, G., Renzini, A., Daddi, E., et al. 2014, MNRAS, 443, 19

Saintonge, A., Catinella, B., Tacconi, L. J., et al. 2017, ApJS, 233, 22

Saintonge, A., Kauffmann, G., Kramer, C., et al. 2011a, MNRAS, 415, 32

Saintonge, A., Kauffmann, G., Wang, J., et al. 2011b, MNRAS, 415, 61

Saintonge, A., Lutz, D., Genzel, R., et al. 2013, ApJ, 778, 2

Saintonge, A., Wilson, C. D., Xiao, T., et al. 2018, MNRAS, 481, 3497

Sanders, D. B., Mazzarella, J. M., Kim, D.-C., Surace, J. A., & Soifer, B. T. 2003, AJ, 126, 1607

Sandstrom, K. M., Leroy, A. K., Walter, F., et al. 2013, ApJ, 777, 5

Santini, P., Fontana, A., Grazian, A., et al. 2012, A&A, 538, A33

Santini, P., Maiolino, R., Magnelli, B., et al. 2010, A&A, 518, L154

Santini, P., Maiolino, R., Magnelli, B., et al. 2014, A&A, 562, A30

Sargent, M. T., Béthermin, M., Daddi, E., & Elbaz, D. 2012, ApJL, 747, L31

Sargent, M. T., Daddi, E., Béthermin, M., et al. 2014, ApJ, 793, 19

Sargent, M. T., Daddi, E., Béthermin, M., & Elbaz, D. 2013, AAABS, 4, 116

Schinnerer, E., Groves, B., Sargent, M. T., et al. 2016, ApJ, 833, 112

Schmidt, M. 1959, ApJ, 129, 243

Schombert, J. M., McGaugh, S. S., & Eder, J. A. 2001, AJ, 121, 2420

Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74

Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1

Scoville, N., Aussel, H., Sheth, K., et al. 2014, ApJ, 783, 84

Scoville, N., Lee, N., Vanden Bout, P., et al. 2017, ApJ, 837, 150

Scoville, N., Sheth, K., Aussel, H., et al. 2016, ApJ, 820, 83

Scoville, N. Z. 2013, in Evolution of Star Formation and Gas, ed. J. Falcón-Barroso & J. H. Knapen (Cambridge: Cambridge Univ. Press), 491

Shen, S., Mo, H. J., White, S. D. M., et al. 2003, MNRAS, 343, 978

Shi, Y., Armus, L., Helou, G., et al. 2014, Natur, 514, 335

Shi, Y., Wang, J., Zhang, Z.-Y., et al. 2016, NatCo, 7, 13789

Silverman, J. D., Daddi, E., Rodighiero, G., et al. 2015, ApJL, 812, L23

Silverman, J. D., Rujopakarn, W., Daddi, E., et al. 2018, ApJ, 867, 92

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Smail, I., Ivison, R. J., & Blain, A. W. 1997, ApJL, 490, L5

Smith, M. W. L., Clark, C. J. R., De Looze, I., et al. 2019, MNRAS, 486, 4166

Smolčić, V., Novak, M., Bondi, M., et al. 2017, A&A, 602, A1

Solomon, P. M., & Barrett, J. W. 1991, in IAU Symp. 146, Dynamics of Galaxies and Their Molecular Cloud Distributions, ed. F. Combes & F. Casoli (Berlin: Springer), 235

Solomon, P. M., Downes, D., Radford, S. J. E., & Barrett, J. W. 1997, ApJ, 478, 144

Solomon, P. M., Rivolo, A. R., Barrett, J., & Yahil, A. 1987, ApJ, 319, 730

Solomon, P. M., & Vanden Bout, P. A. 2005, ARA&A, 43, 677

Song, M., Finkelstein, S. L., Ashby, M. L. N., et al. 2016, ApJ, 825, 5

Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15

Spilker, J., Bezanson, R., Barišić, I., et al. 2018, ApJ, 860, 103

Sutherland, W., Emerson, J., Dalton, G., et al. 2015, A&A, 575, 25

Tacconi, L. J., Genzel, R., Neri, R., et al. 2010, Natur, 463, 781

Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179

Tacconi, L. J., Genzel, R., Smail, I., et al. 2008, ApJ, 680, 246

Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74

Tan, Q., Daddi, E., Magdis, G., et al. 2014, A&A, 569, A98

Thomas, D., Maraston, C., Bender, R., & Mendes de Oliveira, C. 2005, ApJ, 621, 673

Tomczak, A. R., Quadri, R. F., Tran, K.-V. H., et al. 2016, ApJ, 817, 118

Villanueva, V., Ibar, E., Hughes, T. M., et al. 2017, MNRAS, 470, 3775

Walter, F., Decarli, R., Aravena, M., et al. 2016, ApJ, 833, 67

Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, ApJS, 154, 1

Whitaker, K. E., Franx, M., Leja, J., et al. 2014, ApJ, 795, 104

Wilson, C. D. 1995, ApJL, 448, L97

Wright, A. H., Driver, S. P., & Robotham, A. S. G. 2018, MNRAS, 480, 3491

Wright, A. H., Robotham, A. S. G., Driver, S. P., et al. 2017, MNRAS, 470, 283

Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868

Wuyts, E., Kurk, J., Förster Schreiber, N. M., et al. 2014, ApJL, 789, L40

York, D. G., Adelman, J., Anderson, J. E., Jr., et al. 2000, AJ, 120, 1579

Young, J. S., & Scoville, N. Z. 1991, ARA&A, 29, 581