PHIBSS: Unified Scaling Relations of Gas Depletion Time and Molecular Gas Fractions^*

The MS evolution has been interpreted in terms of an "equilibrium-growth/gas-regulator model" (e.g., Bouché et al. 2010; Davé et al. 2011, 2012; Lilly et al. 2013; Dekel & Mandelker 2014; Peng & Maiolino 2014; Rathaus & Sternberg 2016). High SFRs and galaxy growth along the MS are sustained for several gigayears by a continuous supply of fresh gas from the cosmic web and through mostly minor mergers maintaining large gas reservoirs for star formation (e.g., Kereš et al. 2005; Dekel et al. 2009). At z ∼ 1–2.5 MS SFGs double their stellar mass on a typical timescale of ∼0.5–1 Gyr, but their growth appears to halt when they reach the Schechter mass, M_* ∼ 10^10.5–11 M_⊙, and they transition to the sequence of passive galaxies, in a process (or processes) termed "(mass) quenching" (Kauffmann et al. 2003; Conroy & Wechsler 2009; Peng et al. 2010). Simulations suggest that in parallel to the average growth along the MS, SFGs oscillate up and down in sSFR across the MS band (±0.6 dex) on a ∼1 Gyr timescale, owing to increased and decreased gas accretion rates and to internal gas transport and "compaction" events (Tacchella et al. 2016).

A number of studies in the past decade indicate that the MS line demarcates not only the location of the maximum number of galaxies as a function of δMS at constant M_* and z. It is also the location of disk galaxies in terms of structural (n_Sersic ∼ 1; e.g., Franx et al. 2008; Wuyts et al. 2011b; Lang et al. 2014) and kinematic (${v}_{\mathrm{rot}}/{\sigma }_{0}\gg 1$; e.g., Förster Schreiber et al. 2009; Wisnioski et al. 2015) properties, from z ∼ 0 to 2.5. In contrast, going up from the MS, dust temperatures increase at all redshifts (Elbaz et al. 2011; Nordon et al. 2012; Magnelli et al. 2014), the dense gas fraction increases (Gracia-Carpio et al. 2011; Lada et al. 2012), and the ratio of far-IR (FIR) cooling line to continuum luminosity drops ("FIR line deficit"), indicative of local volumetric changes in interstellar medium (ISM) properties (Stacey et al. 2010; Gracia-Carpio et al. 2011; Herrera-Camus et al. 2018).

The equilibrium-growth model predicts a close connection between specific star formation rates (sSFRs), gas fractions, and metallicities as a function of redshift, with only modest changes as compared to z = 0 MS star formation physics (e.g., Elmegreen et al. 2009; Elbaz et al. 2010, 2011; Krumholz & Dekel 2010; Gracia-Carpio et al. 2011; Nordon et al. 2012; Lilly et al. 2013; Peng & Maiolino 2014). In the Milky Way and nearby galaxies most and arguably all star formation occurs in massive (10^4...6 M_⊙), dense (n(H₂) ∼ 10^2...5 cm⁻³), cold (T_gas ∼ 10–30 K), gravitationally bound "giant molecular clouds" (Solomon et al. 1987; McKee & Ostriker 2007; Bolatto et al. 2008), and not in warm atomic gas (Bigiel et al. 2008; Leroy et al. 2008; Schruba et al. 2011). An important open issue is whether the depletion time for converting molecular gas to stars on galactic scales is set locally within clouds (Krumholz & McKee 2005) or on large galactic scales (Elmegreen 1997; Silk 1997). Another is how gas reservoirs change as a function of redshift, stellar mass, SFR, galaxy size/internal structure, gas motions, and environmental parameters (e.g., Bouché et al. 2010; Daddi et al. 2010a, 2010b; Genzel et al. 2010, 2015; Tacconi et al. 2010, 2013; Davé et al. 2011, 2012; Lagos et al. 2011, 2015a, 2015b; Fu et al. 2012; Lilly et al. 2013; Popping et al. 2015).

This work builds on the previous analysis by G15; we summarize here the salient points that motivate the present work. G15 (and references therein) took advantage of the finding of many observational studies during the past decade (see Sections 1.1 and 1.2) that in the physical framework of the evolution of SFGs along the MS, the most important galaxy-integrated, cool and dense ISM properties, namely, the molecular gas content relative to the stellar mass, $\mu ={M}_{\mathrm{molgas}}/{M}_{* }$, and the molecular gas depletion time, ${t}_{\mathrm{depl}}={M}_{\mathrm{molgas}}/\mathrm{SFR}$ (Gyr), mainly depend on cosmic time (or redshift) and on the location along and perpendicular to the MS line at a given redshift. Here M_molgas is the total molecular gas mass, including a 36% mass fraction of helium and a correction for the photodissociated surface layers of the molecular clouds that are fully molecular in H₂ but dissociated ("dark") in CO (see Equation (2); Sternberg & Dalgarno 1995; Wolfire et al. 2010; Bolatto et al. 2013). G15 show that the scaling relations for t_depl and μ_gas can then be written as products of functions depending on redshift, stellar mass, and offset from the MS line, $\delta \mathrm{MS}=\mathrm{sSFR}/\mathrm{sSFR}(\mathrm{MS},z,{M}_{* }$), and only indirectly on the absolute value of the SFR or sSFR. G15 show empirically that this separation of variables is justified, since the depletion time to first order does not depend on stellar mass, and the slope $d\mathrm{log}\,{t}_{\mathrm{depl}}/d\delta \mathrm{MS}$ does not depend significantly on z. We test and reestablish these fundamental assumptions in Section 4.2. With the scaling relations for t_depl established, ratios of molecular gas to stellar mass then follow straightforwardly by multiplying the depletion time with the sSFR of a galaxy, ${\mu }_{\mathrm{gas}}={M}_{\mathrm{molgas}}/{M}_{* }={t}_{\mathrm{depl}}\times (\mathrm{SFR}/{M}_{* }$). Separability is thus possible for the μ_gas dependence as well. G15 also show that good fits are obtained to the scaling relations with a product of power-law functions in the variables above, resulting in linear fitting functions in log–log space.

We caution that this empirical conclusion is not unique, and we further explore in this paper whether more complex fitting functions are required, as the quality of the data improves. We also caution that the parameterization in terms of offset from the MS, which in turn is a function of z and M_*, is well motivated by physical properties (Section 1.2), but is mathematically not unique, in part because MS recipes vary (Figure 1). G15 explored different MS recipes and also fitted directly in log sSFR, log(1 + z), and log(M_*) space. The main difference in fitting in $z,\mathrm{sSFR},{M}_{* }$ from fitting in $z,\delta \mathrm{MS},{M}_{* }$ space is the interpretation of the parameter values of the redshift scaling. While the slope $d\mathrm{log}\,{t}_{\mathrm{depl}}/d\mathrm{log}(1+z)$ at constant δMS = 0 obviously describes the redshift evolution of the entire population, the corresponding slope at constant sSFR is a redshift cut at that selected sSFR and as such does not have a well-defined meaning (G15). For completeness we include also in this paper fits in $z,\mathrm{sSFR},{M}_{* }$ space. In terms of ${{\chi }_{{\rm{r}}}}^{2}$ this fit is somewhat worse than that of the δMS fit.

Finally, the fraction of total cold gas mass to total baryonic mass of a galaxy, the gas fraction, is ${f}_{\mathrm{gas}}\,=({M}_{\mathrm{molgas}}+{M}_{{\rm{H}}{\rm{I}}})/({M}_{* }+{M}_{\mathrm{molgas}}+{M}_{{\rm{H}}{\rm{I}}}$). Here M_{H i} is the integrated atomic hydrogen mass of a galaxy. In massive z ∼ 0 SFGs the integrated atomic molecular hydrogen content dominates the total mass of the neutral ISM, M_{H i} ∼ (2–3)M_molgas (e.g., Saintonge et al. 2011a, and references therein; Catinella et al. 2010, 2013). The evolution of the atomic gas content of galaxies with redshift is relatively poorly known, since the H i 21 cm emission line cannot be detected outside the local universe with available technology. Bauermeister et al. (2010, and references therein) summarize the results coming from UV damped Lyα absorbers toward high-z QSOs and conclude that H i columns likely do not vary strongly with redshift, while the molecular component strongly evolves with redshift (this paper; Daddi et al. 2010b; Tacconi et al. 2010, 2013; G15; Lagos et al. 2015a). For these reasons we make the approximation ${\mu }_{\mathrm{gas}}\sim {\mu }_{\mathrm{molgas}}$ and ${f}_{\mathrm{gas}}\sim {f}_{\mathrm{molgas}}$, which is valid at z > 0.4 (see also Section 4.2.1).

There are currently three main avenues to obtain molecular gas masses:

• 1.  
  The most common and well-established method is the observation of a low-lying CO emission line (CO 1–0, 2–1, 3–2), and using its integrated line luminosity L_CO and a conversion factor (or function) ${\alpha }_{\mathrm{CO}}$ to convert L_CO to molecular gas mass, in the regime where CO comes from optically thick, virialized clouds (e.g., Dickman et al. 1986; Solomon et al. 1987; Bolatto et al. 2013);
• 2.  
  More recently, high-quality FIR/submillimeter dust emission SEDs have become available with the Herschel mission. From these SEDs, dust masses are inferred by fitting dust emissivity models (e.g., Draine & Li 2007). Then molecular gas masses are estimated assuming a gas-to-dust mass ratio (which can be a function of metallicity; see below; e.g., Leroy et al. 2011; Magdis et al. 2011, 2012a; Eales et al. 2012; Magnelli et al. 2012b; Berta et al. 2013, 2016; Rémy-Ruyer et al. 2014; Santini et al. 2014; Sargent et al. 2014; Béthermin et al. 2015; Schreiber et al. 2015; G15). Detailed explanations of this method are given in the above papers, and especially in G15 and Berta et al. (2016).
• 3.  
  A single-band measurement of the dust emission flux on the Rayleigh–Jeans side of the SED (at ∼1 mm) can be used to infer a dust/gas mass, if a single, constant dust temperature is a sufficiently good approximation. The long-wavelength dust emissivity can be estimated from a model, or calibrated from observations of sources in which gas masses are known from method 1 (e.g., Scoville et al. 2014, 2016, 2017).

All three methods have strengths and weaknesses (e.g., Bolatto et al. 2013; G15; Scoville et al. 2016, A. Weiss et al. 2018, in preparation). The "CO method" (method 1) is very well established and calibrated from observations in the local universe (Bolatto et al. 2013), but at high z the CO 1–0 line typically has to be replaced by a higher-energy state line, which requires a calibration of the temperature- and density-dependent ratio ${R}_{1J}={T}_{1}/{T}_{J}$, where T_J is the beam-averaged, Rayleigh–Jeans brightness temperature of the line J $\to$J–1. Observations in low- and high-z SFGs suggest R₁₂ ∼ 1.16–1.3, R₁₃ ∼ 1.8, and R₁₄ ∼ 2.4 (e.g., Weiss et al. 2007; Dannerbauer et al. 2009; Bothwell et al. 2013; Bolatto et al. 2015; Daddi et al. 2015).

With Herschel the FIR dust SED method (method 2) has been widely used for large samples between z = 0 and 1 and for stacks between z = 1 and 2.5. However, the FIR data covering the emission peak are luminosity weighted toward the warm, star-forming dust component and less sensitive to colder dust between star formation regions (Scoville et al. 2016; Carleton et al. 2017). This bias can lead to an underestimate of the ISM mass, especially at low z, where only a fraction of the cold ISM is actively star-forming. The effect will be less at high z, where the entire galaxy is globally unstable to star formation (Genzel et al. 2008, 2011; Elbaz et al. 2011). The Herschel-based references used in this study address this concern by fitting DL07 dust models that implement a distribution of dust temperatures. The rest wavelength coverage needed to constrain the colder dust has been studied by Draine et al. (2007), Magdis et al. (2012a), and Berta et al. (2016). Another concern is that reaching down to the MS line at z > 1 cannot be done with single source detection photometry but requires stacking. We refer to the above papers for discussions of how stacking impacts the deduced masses. Analysis of the same or similar high-redshift data sets can lead to significant differences (0.02–0.4 dex) in inferred dust masses, depending exactly on the assumptions and methodology (Santini et al. 2014; Béthermin et al. 2015; G15; Berta et al. 2016). The details of the different dust SED modeling methods are given in these papers.

Finally the "1 mm" method (method 3) is becoming increasingly popular, as it is much more efficient than the CO method (factors of ∼5–10 in observing time at ALMA), if a constant dust temperature for the emitting dust grains at 1 mm can be assumed (T_dust = 25 K: see Scoville et al. 2016, 2017, A. Weiss et al. 2018, in preparation, for a full discussion). Alternatively, a second photometric measurement at shorter wavelength can constrain the dust temperature, but this can be costly in observing time (G15).

All three methods assume that zero-point calibrations established at z = 0 are valid at higher redshifts without much change. Masses inferred from both CO and dust emission are sensitive to metallicity. In the case of CO, the conversion factor increases with decreasing metallicity Z, because CO is photodissociated to a larger depth in each cloud (Wolfire et al. 2010; Bolatto et al. 2013). Various metallicity-dependent conversion functions have been proposed in the literature. As in G15, we adopt the geometric mean of the ${\alpha }_{\mathrm{CO}}$(Z) recipes of Bolatto et al. (2013) and Genzel et al. (2012) (Equations (6) and (7) in G15),

$$\begin{eqnarray}{\alpha }_{\mathrm{CO}J} & = & 4.36\times {R}_{J1}\\ & & \times \,\sqrt{\begin{array}{l}0.67\times \exp (0.36\times {10}^{-(12+\mathrm{log}({\rm{O}}/{\rm{H}})-8.67)}\\ \quad \times {10}^{-1.27\times (12+\mathrm{log}({\rm{O}}/{\rm{H}})-8.67)}\end{array}}\\ & & ({M}_{\odot }/({\rm{K}}\,\mathrm{km}/{\rm{s}}\ {\mathrm{pc}}^{2})).\end{eqnarray} \tag{ 2 }$$

Here log Z = 12+log(O/H) is the metallicity on the Pettini & Pagel (2004) scale. Molecular gas masses are then computed from Equation (4) in G15. The recipes of Bolatto et al. and Genzel et al. are similar near solar metallicity but then deviate from each other below ∼0.5 Z_⊙. In the subsolar regime the exponential dependence of the Bolatto et al. recipe drives a much steeper increase of α than the power law in the Genzel et al. recipe. As a compromise we took an average of the recipes (the harmonic mean corresponds to the average in log space). We discuss the impact of choosing individual prescriptions in the low-metallicity regime in Section 4.2.3.

In the case of the dust methods 2 and 3 we assume that the dust-to-gas ratio is nearly linearly correlated with metallicity, at least for metallicities 12 + log(O/H) > ∼8, as found in Leroy et al. (2011) and Rémy-Ruyer et al. (2014). We adopt the ratio of molecular gas to dust mass as

$$\begin{eqnarray}&&{\delta }_{{gd}}=\displaystyle \frac{{M}_{\mathrm{molgas}}}{{M}_{\mathrm{dust}}}={10}^{(+2-0.85\times (12+\mathrm{log}({\rm{O}}/{\rm{H}})-8.67))}.\end{eqnarray} \tag{ 3 }$$

Note that we deviate here from the assumption ${\delta }_{{gd}}=\mathrm{const}$. in Scoville et al. (2016) but otherwise use their Equation (16).²⁴

For the few SFGs in this paper with estimates of gas phase metallicities from rest-frame optical strong line ratios, we determine individual estimates of log Z = 12 + log(O/H), adopting the Pettini & Pagel (2004) scale (e.g., Kewley & Ellison 2008, for a detailed discussion). However, for most of the SFGs in the CO and dust samples, such line ratios are not available and it is necessary to use the mass–metallicity relation for the metallicity corrections discussed above. Following G15, we adopt

$$\begin{eqnarray}&&12+\mathrm{log}{({\rm{O}}/{\rm{H}})}_{{PP}04}=a-0.087\times {(\mathrm{log}{M}_{* }-b)}^{2},\mathrm{with}\\ &&\quad a=8.74(0.06),\ \mathrm{and}\\ &&\quad b=10.4(0.05)+4.46(0.3)\times \mathrm{log}(1+z)-1.78(0.4)\\ &&\qquad \qquad \quad \times \,{(\mathrm{log}(1+z))}^{2}.\end{eqnarray} \tag{ 4 }$$

To set the scene, we summarize in Table 1 previous work on the molecular gas scaling relations, including several theoretical papers describing the results from hydro-simulations and semianalytic work.

• 1.  
  Redshift dependence. There is broad qualitative agreement in the literature that the depletion timescale is about 1 Gyr, and dropping by a factor of 2–4 between z = 0 and 2.5 (Bigiel et al. 2008; Daddi et al. 2010b; Genzel et al. 2010; Tacconi et al. 2010, 2013; Saintonge et al. 2011b, 2013, 2016; Leroy et al. 2013; Santini et al. 2014; Sargent et al. 2014, G15). In comparing different estimates, those with a greater redshift range are naturally preferable (see the Appendix, where we investigate the effects that limited parameter space and source statistics have on the inferred scaling relations). Since ${\mu }_{\mathrm{gas}}=\mathrm{sSFR}\times {t}_{\mathrm{depl}}$, a depletion time that is slowly varying with redshifts means that molecular-to-stellar-mass ratios nearly track sSFR(z), and gas fractions increase strongly with redshift. The theoretical work also finds strong redshift evolution of increasing molecular (not atomic!) gas fractions with redshift (e.g., Lagos et al. 2012, 2015a; Genel et al. 2014; Popping et al. 2015).
• 2.  
  Depletion time along and perpendicular to MS. There is also broad qualitative agreement that t_depl decreases as one steps upward in sSFR perpendicular to the MS line at a given z and M_*, as long as one considers a large enough range in δMS (Saintonge et al. 2011b, 2012; Leroy et al. 2013; Tacconi et al. 2013; Huang & Kauffmann 2014; Sargent et al. 2014; G15). The theoretical work is in agreement with these findings (e.g., Lagos et al. 2012, 2015a; Genel et al. 2014; Popping et al. 2015). Depletion times are constant or increase slowly, stepping along the MS line. Since ${\mu }_{\mathrm{gas}}=\mathrm{sSFR}\times {t}_{\mathrm{depl}}$, this implies that gas fractions track the mass dependence of sSFR (Saintonge et al. 2011b).
• 3.  
  Other parameters. Studies of scalings on kiloparsec scales within galaxies are so far only available at z ∼ 0, in particular through the HERACLES survey (Bigiel et al. 2008, 2011; Leroy et al. 2008, 2013). The HERACLES data do not exhibit strong intragalactic parameter dependencies, with the exception of a significant drop of t_depl at low galaxy masses, which likely reflects the impact of UV photodissociation in low-metallicity ISM clouds and a resulting change in ${\alpha }_{\mathrm{CO}}$ (Leroy et al. 2013). The HERACLES data do not show a significant dependence of t_depl on (gas, stellar) density or galactic radius, with the exception of a drop in the circumnuclear regions. There is also no significant change in t_depl (or its inverse, the star formation efficiency) between arm and interarm regions in M51, NGC 628, and NGC 6946 (Foyle et al. 2010). In contrast, Huang & Kauffmann (2015) find a significant dependence on star formation and stellar surface density (log t_depl ∼ −0.36 × log Σ_SFR − 0.5 × log Σ_*) from a different analysis of HERACLES, in combination with COLD GASS. From COLD GASS alone Saintonge et al. (2012, 2016) find that the depletion time may increase with stellar surface density, as the mass fraction of quenched bulges/spheroids increases. Expanding a subset of the HERACLES sample with HCN 1–0/CO 1–0 line ratios (as dense gas tracer), Usero et al. (2015) find that the L(IR)-to-L(HCN) ratio, thought to be closely related to the star formation efficiency of dense molecular gas, decreases systematically with these same parameters and is much lower near galaxy centers than in the outer regions of the galaxy disks. For fixed conversion factors, these results are incompatible with a simple model in which star formation depends only on the amount of gas mass above some density threshold (Lada et al. 2012).

Table 1.  Summary of Scaling Relations in the Current Literature

Method	z Range	$d\mathrm{log}\,{t}_{\mathrm{depl}}/d\mathrm{log}(1+z)$	$d\mathrm{log}{\mu }_{\mathrm{gas}}/d\mathrm{log}(1+z)$	$d\mathrm{log}\,{t}_{\mathrm{depl}}/d\mathrm{log}({\boldsymbol{\delta }}\mathrm{MS})$	$d\mathrm{log}\,{t}_{\mathrm{depl}}/d\mathrm{log}({M}_{* })$	$d\mathrm{log}{{\boldsymbol{\mu }}}_{\mathrm{gas}}/d\mathrm{log}({M}_{* })$	Scaling with Other Parameters	${dSFR}/{{dM}}_{\mathrm{molgas}}d{{\rm{\Sigma }}}_{\mathrm{SFR}}/d{{\rm{\Sigma }}}_{\mathrm{molgas}}$	Reference
CO COLD GASS galaxy integrated	0.025–0.05(N = 440 incl. nondetects.)			$-{0.5}_{0.1}$	$+{0.36}_{0.1}$	$-{0.24}_{0.1}$	dlog tdepl/dlog ${{\rm{\Sigma }}}_{* }$ ∼ −0.5, dlog tdepl/dlog ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ ∼ −0.36 no bars/interactions no or weak AGNs	${1.18}_{0.24}$	Saintonge+11a, b, 12, 13, 116, Huang & Kauffmann 14, 15 Accurso+17
CO HERACLES galaxy 1 kpc resolved	6-35e-4 (N = 30)			$-{0.25}_{0.25}$	$+{0.2}_{0.15}$		weak with radius, spiral arms, surface density, except nucleus	${1}_{0.15}$ (>1.5 kpc)	Bigiel +08, 11 Leroy+08, 13, Schruba +11
CO PHIBSS1	1–2.5 (N = 52 detects.)	$-{0.7}_{0.3}$	$+{2.7}_{0.4}$	$-{0.32}_{0.2}$		$-{0.6}_{0.15}$		${1.1}_{0.2}$	Tacconi+10, 13 Genzel+10
CO + dust FIR/submillimeter SED (plus literature)	0–4 (N = 131 CO detects. 15 FIR stacks)	$-{0.8}_{0.3}$	$+{2.6}_{0.4}$	$-{0.2}_{0.1}$			"two-mode" SF, with strong increase of bursts for δMS > 4	${1.25}_{0.15}$	Daddi +10b, Magdis+12a, Sargent+14, Bethermin+15
Dust FIR SED	0.3–2 (N = 121 stacks)	$-{1.5}_{0.4}$	$+{1.6}_{0.5}$	$-{0.3}_{0.1}$					Santini+14
CO PHIBSS1+2 +dust FIR SED (plus literature)	0–2.5 (N = 500 CO > 4σ detects., 512 FIR stacks)	$-{0.34}_{0.3}$	$+{2.7}_{0.2}$	$-{0.49}_{0.05}$	$+{0.01}_{0.1}$	$-{0.37}_{0.1}$		${1.1}_{0.15}$ increases with δMS	Genzel+15 (including data from Magnelli +14, cf. Berta+16)
Dust 1 mm	1–4.4 (N = 51 > 4σ detects., 35 stacks)	$-{1.2}_{0.4}$	+${1.8}_{0.5}$	$-{0.55}_{0.1}$	$+{0.23}_{0.2}$	$-{0.02}_{0.2}$		${0.9}_{0.1}$	Scoville+16
Dust 1 mm	0.3–4.5 (N = 708)	$-{1.05}_{0.05}$	$+{1.8}_{0.14}$	$-{0.7}_{0.02}$	$-{0.01}_{0.01}$	$-{0.7}_{0.04}$		1.0	Scoville+17
Semianalytic model	All		$+{2.5}_{0.3}$						Popping+15
Illustris hydro sim	All		$+{1.8}_{0.4}$			$-{0.6}_{0.2}$			Genel+14
Eagle hydro sim	All	$-{1.3}_{0.4}$	$+{1.5}_{0.4}$						Lagos+15a

Note. The subscripts in columns (3)–(9) are the formal uncertainties as listed in the individual studies.

Download table as:  ASCIITypeset image

Table 1 shows that despite the qualitative agreement mentioned above, uncertainties are large enough and methodologies sufficiently different to result in ambiguous and sometimes contradictory conclusions. The different methodologies for computing SFRs and stellar masses alone can lead to significant differences. At z = 0, for instance, Bigiel et al. (2008) and Leroy et al. (2008, 2013) find t_depl(MS) ∼ 2–2.5 Gyr versus 1–1.5 Gyr from Saintonge et al. (2011b, 2012), emphasizing the importance of homogenized definitions and calibrations (benchmarking). At z = 0, Bigiel et al. (2008, 2011) and Leroy et al. (2008, 2013) find a constant depletion time, with the exception of galactic nuclei, while Saintonge et al. (2011b, 2012) and Huang & Kauffmann (2014, 2015) emphasize that depletion time is correlated inversely with sSFR. The different conclusions could be caused in part by the narrower range of sSFRs covered in the HERACLES sample analyzed in Bigiel et al. (2008, 2011) and Leroy et al. (2008, 2013), relative to COLD GASS. The derived values for the slope dlog t_depl/dlog δMS vary from −0.2 to −0.7 (Saintonge et al. 2012, 2013; Santini et al. 2014; Sargent et al. 2014; G15; Scoville et al. 2016, 2017). Simulations and semianalytic work continue to find somewhat lower (by a factor of 1.5–2) SFRs and gas fractions at high z than implied by most observations (e.g., Davé et al. 2011, 2012; Genel et al. 2014; Lagos et al. 2015a).

Another important example is the interpretation of the MS itself. A massive (M_* ∼ 10¹¹ M_⊙) galaxy on the midplane of the MS at z ∼ 2.3 has an SFR of ∼200 M_⊙yr⁻¹. The same SFR is reached at z = 0 only for (ultra)luminous infrared galaxies ((U)LIRGs), which are placed well above the MS by major-merger-triggered starbursts (Sanders & Mirabel 1996). This does not mean that all z = 2.3 MS SFGs are mergers, however. On the contrary, the advent of Herschel SEDs (see Elbaz et al. 2011; Nordon et al. 2012) and the firm establishment of the MS picture discussed in Section 1 favor an explanation where the increase in SFRs on the MS in Equation (1) implies that normal star-forming disk galaxies at high z are more gas-rich. While the former ("starburst") interpretation had dominated earlier work, the more recent results now favor the "equilibrium-growth model" (Section 1.1; see Elbaz et al. 2011; Nordon et al. 2012). Santini et al. (2014) and Scoville et al. (2016) find only a moderate slope of μ_gas with z, suggesting the need for a more efficient star formation process at high z, perhaps driven by the increased merger rates. In contrast, G15 find that μ_gas changes rapidly with redshift, tracking sSFR(z) and favoring a single dominant star formation process on the MS at all redshifts between z = 0 and 2.5. Daddi et al. (2010a, 2010b), Genzel et al. (2010), Magdis et al. (2012a), and Sargent et al. (2014) all applied a Galactic CO conversion factor near the MS, but a much smaller one above the MS, motivated by observations of local ULIRGs (Downes & Solomon 1998). This resulted in the proposal that galactic star formation is "bi-modal," with large depletion timescales (low efficiency) at δMS ∼ 1 and short depletion timescales (high efficiency) at $\delta \mathrm{MS}\gg 1$, with a fairly sudden transition in between. G15 did not see much difference in the dependence of t_depl on δMS between CO and dust-based data, suggesting that strong ${\alpha }_{\mathrm{CO}}(\delta \mathrm{MS}$) changes are probably not justified, at least at high z. The work of Scoville et al. (2016, 2017) independently reaches the same conclusion.

As we show below, many of these differences are caused by output parameter estimation of often modest data sets across limited input variable ranges, in addition to the systematic differences in calibrations of the input variables entering the analyses. By substantially increasing these variable ranges and the total number of data points and including different calibration methods, we expect significantly more robust results.

We collected 667 CO detections of SFGs from a number of molecular surveys with CO 1–0, 2–1, and 3–2 (and in two cases 4–3) rotational line emission. These data cover the redshift range from z = 0 to 4.0, the stellar mass range of M_* = 10^9.0 to 10^11.8 M_⊙ (M_*  <  10¹⁰ M_⊙ for z = 0 only), and, at a given redshift and stellar mass, SFRs from about 10⁻¹ to 10² times the MS SFR. We include the following:

• 1.  
  216 detections of CO 1–0 emission above and below the MS between z = 0.025 and 0.05 from the final xCOLD GASS survey with the IRAM 30 m telescope (Saintonge et al. 2011a, 2011b, 2016, 2017), and the single COLD GASS stack detection of galaxies much below the MS. We also include 89 detections of the low-mass extension of XCOLD GASS (log(M_*/M_⊙) = 9.0–10.0; Saintonge et al. 2017). We note that the SFRs in COLD GASS have been updated from earlier UV/optical SED fitting (Saintonge et al. 2011a) with mid-IR (MIR) SFRs from WISE (Huang & Kauffmann 2014; Saintonge et al. 2016).
• 2.  
  90 CO 1–0 detections with the IRAM 30 m of z = 0.002−0.09 luminous and ultraluminous IR galaxies (LIRGs and ULIRGs) from the GOALS survey (Armus et al. 2009), from the work of Gao & Solomon (2004), Gracia-Carpio et al. (2008, 2009, 2018, private communication), and Garcia-Burillo et al. (2012).
• 3.  
  31 CO 1–0 or 3–2 detections of above MS SFGs between z = 0.06 and 0.5 with the CARMA millimeter array from the EGNOG survey (Bauermeister et al. 2013).
• 4.  
  14 CO 2–1 or 3–2 detections at z = 0.6–0.9 and 18 CO 1–0 detections at z = 0.2–0.58 (significantly above MS) of ULIRGs with the IRAM 30 m telescope from Combes et al. (2011, 2013).
• 5.  
  51 detections of CO 3–2 emission in MS SFGs in two redshift slices at z = 1–1.5 and z = 2–2.5 as part of the PHIBSS1 survey with the IRAM PdBI (now NOEMA; Tacconi et al. 2010, 2013).
• 6.  
  97 detections of CO 2–1 or 3–2 in MS SFGs between z = 0.5 and 2.7, as part of the PHIBSS2 survey with the updated IRAM NOEMA interferometer (G15; Freundlich et al. 2018; PHIBSS2 team 2018, in preparation).
• 7.  
  9 CO 2–1 or 3–2 detections of near MS SFGs between z = 0.5 and 3.2 from Daddi et al. (2010a) and Magdis et al. (2012b), obtained with the IRAM PdBI.
• 8.  
  6 CO 2–1 detections of z = 1–1.2 MS SFGs selected from the Herschel-PEP survey (Lutz et al. 2011), obtained with the IRAM PdBI (Magnelli et al. 2012a).
• 9.  
  19 CO 2–1, 3–2, or 4–3 detections of above-MS submillimeter galaxies (SMGs) between z = 1.2 and 3.4, obtained with the IRAM PdBI by Greve et al. (2005), Tacconi et al. (2006, 2008), and Bothwell et al. (2013).
• 10.  
  8 CO 3–2 detections of z = 1.4–3.2 lensed MS SFGs obtained with the IRAM PdBI (Saintonge et al. 2013, and references therein).
• 11.  
  10 ALMA CO 3–2 and 2–1 detections between z = 1 and 2.5 (DeCarli et al. 2016; R. Genzel et al. 2018, in preparation).
• 12.  
  7 CO 2–1 and 3–2 detections of z = 1.4–2.2 "outliers" above the MS, with ALMA and NOEMA (Silverman et al. 2015).

We take from the literature thermal continuum dust observations from Herschel and ALMA, which have molecular gas masses estimated from methods 2 and 3 described above. We include the following:

1. 512 (Magnelli et al. 2014; Berta et al. 2016, G15) and 121 (Santini et al. 2014) stacks of z = 0.1–1.9 SFGs with deep Herschel PACS/SPIRE spectrophotometry in the COSMOS and GOODS (N/S) fields as part of the PEP (Lutz et al. 2011) and HerMES (Oliver et al. 2012) surveys, and 15 stacks of z = 0.4–3.8 SFGs in COSMOS, again with deep PACS/SPIRE photometry, plus additional short- and long-wavelength coverage with Spitzer, LABOCA, and AzTEC (Béthermin et al. 2015). In all cases we adopted dust masses from these references and converted to gas mass as described in Section 2.1. For the final analysis, we removed 136 stacks at z < 0.4 that are likely significantly affected by dust in H i gas, as discussed above.

2. 102 detections and 21 stacks of z = 0.9–4.4 of 850–1300 μm continuum fluxes from recent ALMA observations (Tadaki et al. 2015, 2017; Barro et al. 2016; DeCarli et al. 2016; Scoville et al. 2016; Dunlop et al. 2017; S. Lilly et al. 2018, in preparation), using the methodology and calibration proposed by Scoville et al. (2016) with T_dust = 25 K = const., but correcting for the metallicity dependence of gas-to-dust ratios as discussed above. In the case of Scoville et al. (2016) we included the 72 individual detections with significance $\geqslant 4\sigma$ and redshifts verified by additional data (S. Wuyts 2018, private communication). For the 21 stacks in the z ∼ 1 and ∼2 bands presented by Scoville et al. (2016), 18 (0.86) have a significance >4σ and were included. We also included the average of 1 mm photometry detections of 45 SFGs between z = 2.8 and 3.8 as published in Schinnerer et al. (2016). Finally, we included six 1 mm dust continuum detections between z = 1.2 and 2.3 with the IRAM NOEMA interferometer (from the PHIBSS2 survey), for a total of 108 individual detections and 22 stacks in the 1 mm photometry technique.

Our basic approach is that the core, near-MS, CO, or dust data sets are "benchmark" subsamples of large panchromatic (UV/optical/infrared/radio) imaging surveys, preferably with spectroscopic redshifts, and with well-established and relatively homogeneous stellar and star formation properties. We note that we eliminated some of the z ∼ 4 data points of Scoville et al. (2016), since the grism redshifts from 3D-HST and the literature were discrepant with the photometric redshifts used in that paper. The xCOLD GASS sample is mass selected from the Sloan Digital Sky Survey (SDSS; Saintonge et al. 2011a, 2011b, 2016, 2017). PHIBSS1 and 2 and the Herschel and ALMA dust samples are selected from deep rest-frame UV/optical imaging surveys in the EGS (Davis et al. 2007; Cooper et al. 2012; Newman et al. 2013), GOODS N/S (Giavalisco et al. 2004; Berta et al. 2010), and COSMOS fields (Lilly et al. 2007, 2009; Scoville et al. 2007), including the recent CANDELS J- and H-band HST imaging (Grogin et al. 2011; Koekemoer et al. 2011) and 3D-HST grism spectroscopy (Brammer et al. 2012; Skelton et al. 2014; Momcheva et al. 2016). We also include galaxies from the Deep-3a survey (Kong et al. 2006) and the BX/BM surveys of Steidel et al. (2004) and Adelberger et al. (2004). We have supplemented these core samples with smaller data sets addressing outliers above the MS, mainly starburst sources (LIRGS, ULIRGs, submillimeter galaxies, etc.) as described in Sections 3.1 and 3.2.

With these selections it is possible to place the basic galaxy parameters—stellar masses, SFRs, effective radii in the rest-frame optical—on a common "ladder" system (see Wuyts et al. 2011a, 2011b, and Saintonge et al. 2011a, 2016, for details), where SFRs are based on FIR emission, MIR emission, and UV to NIR SED fitting in decreasing preference. The impact of the availability or absence of MIR/FIR photometry is quite important, as the comparison of the left and right panels of Figure 1 shows. SFRs based only on optical/UV SED analysis tend to underestimate the total SFRs, which is especially relevant at high redshifts, at low stellar masses, and below the MS. Wherever necessary and possible, we adjusted the stellar masses and SFRs from the literature to the same assumptions. Typical fractional stellar mass uncertainties (including systematic errors) are ±0.13 dex on the MS and ±0.2 dex for outliers, SFR uncertainties are ±0.2 dex for Herschel/Spitzer-detected galaxies and ±0.25 dex for SED-inferred SFRs, or starbursts, and gas mass uncertainties are ±0.23 dex (G15). Note that throughout the paper we define stellar mass as the "observed" mass ("live" stars plus remnants), after mass loss from stars. This is about 0.15–0.2 dex smaller than the integral of the SFR over time. The redshift–sSFR–M_* coverage of the various samples is shown in Figure 2, with the different symbols denoting the various surveys mentioned in our listing above.

Zoom In Zoom Out Reset image size

Figures 2 and 3 give the distribution of our data in the log sSFR–log(1 + z) plane separately for the three methods, in the log sSFR–log M_* (at z = 0 for the CO method) and in the $\mathrm{log}{R}_{e}$(5000 Å)–log(1 + z) plane. We remind the reader that it is not realistic to construct an unbiased gas sample, whose distribution function in these planes is proportional to the distribution function of the parent sample. Rather, the question is how broad and unbiased the coverage in each of the relevant parameters is, from which we then can attempt to determine scaling relations.

Zoom In Zoom Out Reset image size

Overall the parameter coverage of the combined sample is quite broad: redshift z = 0–4.4, stellar mass from log(M_*/M_⊙) = 9 to 11.9, SFRs from 10⁻² to 10² times those on the MS reference line, and sizes 0.25–2 times the typical size at a given mass and redshift. The best coverage occurs at z = 0 with the xCOLD GASS and LIRG/ULIRG CO surveys. Unfortunately, we have no access to equivalent surveys in the dust tracers of the molecular ISM, due to the strong contribution from dust in the atomic medium (e.g., Draine et al. 2007; Kennicutt et al. 2011; Dale et al. 2012; Eales et al. 2012; Sandstrom et al. 2013). The low-mass galaxy coverage from log(M_*/M_⊙) = 9 to 10 is also exclusively from the xCOLD GASS CO survey. In the mid-z range (z = 0.5–2.5) we have the best comparison of the three independent methods, and with several independent analysis methods of the SED technique. Owing to sensitivity limits, the overall distribution in z-sSFR space at all z is somewhat biased to SFGs on and above the MS and at higher stellar masses ($\langle \mathrm{log}\,\delta \mathrm{MS}\rangle =0.2\mbox{--}0.34$). However, the recent extensive surveys at the IRAM telescopes at z ∼ 0.03 (xCOLD GASS), z ∼ 0.7 (PHIBSS2), z ∼ 1.2 (PHIBSS1+2), and z = 2.2 (PHIBSS1 + 2) now establish good coverage of massive SFGs above and below the MS line (Figure 2). In comparison, the other large survey of the molecular gas properties in high-z SFGs from 1 mm ALMA dust photometry (Scoville et al. 2016, 2017) has $\langle \mathrm{log}\,\delta \mathrm{MS}\rangle =0.5\mbox{--}1$ between z = 1 and 3 and thus is even more strongly biased to above MS galaxies.

As mentioned above, with the exception of xCOLD GASS (Saintonge et al. 2017), none of the other samples reach below M_* ≤ 10¹⁰ M_⊙. This is in part because of sensitivity limitations for the higher-redshift surveys, but more importantly it is also per design, as in this mass range the (line or continuum) luminosity per ISM mass decreases rapidly owing to the metallicity dependence of ${\alpha }_{\mathrm{CO}}$ and ${\delta }_{{gd}}$ (Section 2.1). Given these limitations, we discuss in Sections 4.2.3 and 4.3 what we can infer from xCOLD GASS for the stellar mass dependence of the scaling relations for M_* ≤ 10¹⁰ M_⊙, but we caution the reader that the relations at the low-mass end are not represented by galaxies with z > 0.05.

In the framework of the MS prescription (Equation (1)) and following the analysis in G15 (Section 2) and other papers, our Ansatz is to separate the parameter dependencies of t_depl as products of power laws, first in redshift (1 + z), next at a given redshift above and below the MS at a fixed stellar mass ($\delta \mathrm{MS}=\mathrm{sSFR}/\mathrm{sSFR}(\mathrm{MS},z,{M}_{* }$)), and then along the MS (δM_* = M_*/5 × 10¹⁰M_⊙, corrected to fiducial stellar mass of 5 × 10¹⁰ M_⊙). Finally, we investigate the residuals as a function of effective radius (half-light radius in rest-frame optical [5000 Å] R_e), relative to the average radius of the star-forming population, ${R}_{e0}=8.9$ kpc (1 + z)^−0.75 (M_*/5 × 10¹⁰M_⊙)^0.23 (van der Wel et al. 2014), such that $\delta R={R}_{e}/{R}_{e0}$. This means that

$$\begin{eqnarray}&&\mathrm{log}({t}_{\mathrm{depl}}(z,\mathrm{sSFR},{M}_{* },{R}_{e}))=\,{A}_{t}+{B}_{t}\times \mathrm{log}(1+z)\\ &&\quad +\,{C}_{t}\times \mathrm{log}(\delta \mathrm{MS})+\,{D}_{t}\times \mathrm{log}\delta {M}_{* }\\ &&\quad +\,{E}_{t}\times \mathrm{log}(\delta R).\end{eqnarray} \tag{ 5 }$$

Separation of variables requires that the parameters C_t, D_t, and E_t should not depend significantly on redshift. We show in Sections 4.2.3 and 4.2.4 that D_t and E_t are indeed close to zero and can be neglected to first order.

To explore the redshift dependence of C_t, we first split the independent data sets of each of the three methods (CO, dust FIR, dust 1 mm photometry) into six redshift bins (z = 0–0.1, 0.1–0.5, 0.5–0.9, 0.9–1.6, 1.6–2.5, and 2.5–4.4). In each of the redshift bins and separately for each of the three methods, we fit for ${A}_{{\rm{t}}}^{{\prime} }$ = A_t + B_t × log(1 + z) and C_t. In the literature there are three independent analyses of the FIR/submillimeter dust SEDs from Herschel (+Spitzer, ground-based: Berta et al. 2016 and G15 [using fluxes from Magnelli et al. 2014], Santini et al. 2014, and Béthermin et al. 2015), in which t_depl and μ_gas were determined in stacks, as a function of z, δMS, and δM_*. We analyzed each of these data sets separately. The right panel of Figure 4 summarizes the inferred slopes C_t. There is no overall significant redshift trend of C_t. The distribution of individual values of C_t around the best-fit, error-weighted average (−0.44) has a scatter of ±0.22 dex, somewhat larger than the median uncertainty of the individual data points (±0.15 dex). If data points of Béthermin et al. (2015) are not considered, that scatter further decreases to ±0.16 dex, suggesting that the data can be described mostly with scatter around a constant, redshift-independent slope, in excellent agreement with the Separation Ansatz. The C_t values inferred from the Magnelli et al. (2014) and Berta et al. (2016) points (with very similar input data) differ on average by ΔC_t ∼ 0.36 dex. This suggests that the scatter is significantly affected by systematics in the analyses. We adopt the overall best slope of C_t = −0.44, which includes all the data sets listed in Section 3 (see Table 3(a)).

Zoom In Zoom Out Reset image size

Table 3.  (a) Fit Parameters for Equation (5) Obtained from Error-weighted, Multiparameter Regression. (b) Fit Parameters for Equation (6) Obtained from Error-weighted, Multiparameter Regression

(a) log tdepl (Gyr) = At + Bt*log(1 + z) + Ct*log(sSFR/sSFR(MS, z, M*)) + ${D}_{{\rm{t}}}^{* }$(log M*−10.7) + Et*log(Re/Re0(z, M*))
Data	Parameter	N	${{\chi }^{2}}_{r}$	${A}_{t}$ a	${B}_{t}$	${C}_{t}$	${D}_{t}$	${E}_{t}$
All (error weighted) (after removal of z < 0.4 dust points)	tdepl (Gyr) S14	1309	0.86	+0.090.03	−0.620.08	−0.440.03	+0.0950.03	−
All (err.w.) including z = 0.1–0.4 dust points	tdepl (Gyr) S14	1444	0.85	+0.120.03	−0.660.1	−0.450.03	+0.090.03	⋯
SFGs with 5000 Å Re (error weighted)	tdepl (Gyr) S14	512	0.85	+0.0060.02				+0.110.08
All (equal weight per z bin)	tdepl (Gyr) S14	1309	0.73	+0.090.03	$-{0.52}_{0.1}$	−0.450.03	$+{0.02}_{0.04}$	⋯
CO(error w.)	tdepl (Gyr) S14	667	0.89	+0.060.03	$-{0.44}_{0.13}$	−0.430.03	+0.170.04	⋯
Dust FIR (error w.)	tdepl (Gyr) S14	512	0.56	$+{0.25}_{0.13}$	$-{1.0}_{0.5}$	$-{0.53}_{0.07}$	$-{0.07}_{0.07}$	⋯
Dust 1 mm (error w.)	tdepl (Gyr) S14	130	1.3	$+{0.42}_{0.45}$	−0.70.9	$-{0.64}_{0.15}$	$-{0.27}_{0.2}$	⋯
Best (with bootstrap errors)	tdepl (Gyr) S14	1309		$+{\bf{0}}.{{\bf{09}}}_{{\bf{0.05}}}$	$-{\bf{0}}.{{\bf{62}}}_{{\bf{0.13}}}$	−0.440.04	+0.090.05	$+{\bf{0}}.{{\bf{11}}}_{{\bf{0.12}}}$
Best (with bootstrap errors)	tdepl (Gyr) W14	1309		+0.0020.04	−0.370.08	−0.400.04	+0.170.05	⋯
All (error weighted) (after removal of z < 0.4 dust points)	tdepl (Gyr) S14 fit sSFR instead of δMS	1309	0.86	−0.530.04	${0.95}_{0.15}$	−0.450.04	−0.080.05

(b) log Mmolgas/M* = Aμ + Bμ*(log(1+z)–Fμ)β + Cμ*log(sSFR/sSFR(MS, z, M*)) + Dμ*(log M*-10.7) + Eμ*log(Re/Re0(z, M*))b
Data	Parameter	N	${{\chi }^{2}}_{r}$	${A}_{\mu }$ a	${B}_{\mu }$	${F}_{\mu }$ c	${C}_{\mu }$	${D}_{\mu }$	${E}_{\mu }$
All (error w.)	μ = Mmolgas/M*	1309	0.7	+0.070.15	−3.80.4	+0.630.1	+0.530.03	−0.330.03	⋯
β = 2 in log(1+z)	S14
β = 0			0.95	−1.350.03	+2.950.08	⋯	+0.540.03	−0.310.03
All with optical Re (error w.)	Mmolgas/M*	512	0.85	−0.0120.03	⋯	⋯	⋯	⋯	+0.110.07
	S14
All (equal weight per z bin)	Mmolgas/M*	1309	0.7	+0.170.15	−3.250.4	+0.70.1	+0.530.03	−0.380.04	⋯
	S14
CO (ew.) β = 2	Mmolgas/M*	667	0.8	+0.190.24	−3.380.8	+0.70.17	+0.560.03	−0.300.04	⋯
β = 0	S14		0.89	−1.40.03	+3.10.12	⋯	+0.560.03	−0.280.04
Dust FIR (error w.)	Mmolgas/M*	512	0.6	−0.90.1	+1.90.3	⋯	+0.460.15	−0.400.07	⋯
β = 0	S14
Dust 1 mm (error w.)	Mmolgas/M*	130	1.2	−0.440.5	+1.00.8	⋯	+0.360.15	−0.520.2	⋯
β = 0	S14
All (error w.)	Mmolgas/M*	1309	0.7	+0.070.15	−3.80.4	+0.630.1	+0.530.03	−0.330.03	⋯
β = 2 in log(1+z)	S14
β = 0			0.95	−1.350.03	+2.950.08	⋯	+0.540.03	−0.310.03
Bestc β = 2	Mmolgas/M*	1309		+0.120.15	−3.620.4	+0.660.1	+0.530.03	−0.350.03	+0.110.1
β = 0	S14			−1.190.04	+2.490.2	⋯	+0.520.03	−0.360.03	+0.110.1
Bestc β = 2	Mmolgas/M*	1309		+0.160.15	−3.690.4	+0.650.1	+0.520..03	−0.360.03	⋯
β = 0	W14			−1.250.03	+2.60.25		+0.530.03	−0.360.03
All (error w.)	μ = Mmolgas/M*	1309	0.7	+0.070.15	−3.80.4	+0.630.1	+0.530.03	−0.330.03	⋯
β = 2 in log(1+z)	S14
β = 0	fit sSFR instead of δMS		0.95	−1.350.03	+2.950.08	⋯	+0.540.03	−0.310.03

Notes.

^aAfter introduction of zero-points for each method: zero(CO) = +0.03 dex, zero(FIR Berta/Magnelli) = −0.22, zero(FIR Santini) = +0.21 zero, (FIR Bethermin) = −0.023, zero(1 mm) = −0.003. ^bR_e0 is the mean effective radius of the star-forming population as a function of z and M_*, as derived by van der Wel et al. (2014) based on the HST CANDELS data: ${R}_{e0}=8.9\,\mathrm{kpc}$ (1 + z)^−0.75(M_*/5 × 10¹⁰M_⊙)^0.22. ^cZero-point offset for β = 2. If the data are fit with $\beta =0\to F=0$, such that log M_molgas/M_* = A + B*log(1 + z) + C*log(sSFR/sSFR(MS, z, M_*)) + D*(log M_*–10.7) +E*log(${R}_{e}/\langle {R}_{e}(z,{M}_{* })\rangle$).

Download table as:  ASCIITypeset image

Once the scaling relations for t_depl are determined, it is in principle straightforward to determine the equivalent relations for the ratio of molecular gas to stellar mass, μ_gas, from the combination of Equations (1) and (5),

$$\begin{eqnarray}&&\mathrm{log}(\mu (z,\mathrm{sSFR},{M}_{* },{R}_{e}))=\mathrm{log}({M}_{\mathrm{molgas}}/{M}_{* })\\ &&=\,\mathrm{log}\,{t}_{\mathrm{depl}}+\mathrm{log}(\mathrm{sSFR})\\ &&=\,\mathrm{log}\,{t}_{\mathrm{depl}}+\mathrm{log}(\mathrm{sSFR}(\mathrm{MS},z,{M}_{* }))\\ &&\quad +\,\mathrm{log}\left(\displaystyle \frac{\mathrm{sSFR}}{\mathrm{sSFR}(\mathrm{MS},z,{M}_{* })}\right)\\ &&=\,{A}_{\mu }+{B}_{\mu }\times {(\mathrm{log}(1+z)-{F}_{\mu })}^{\beta }+{C}_{\mu }\times \mathrm{log}(\delta \mathrm{MS})\\ &&\quad +\,{D}_{\mu }\times \mathrm{log}(\delta {M}_{* })+{E}_{\mu }\times \mathrm{log}(\delta R).\end{eqnarray} \tag{ 6 }$$

Unfortunately, the slope of the MS line (dsSFR(MS)/dlog(1 + z)) in S14 and also in W14 varies quite strongly with redshift (from +3.6 at z ∼ 0–1 to +1.2 at z > 2.2, at log M_* ∼ 10.8) and mass. Linear functions in log(1 + z) and log δM_*, based on the scaling relations obtained in Equations (1) and (5), thus are not sufficient. To capture the slope variations, we introduced two more parameters, ${F}_{\mu }$ and β (=2), as shown in Equation (6). Because of the slope variations of sSFR(MS) in z, fitting of the data is quite sensitive to the range and distribution of data points (especially for the pure power-law case β = 0 and ${F}_{\mu }=0$). For this reason we also fit the data by first binning in z and then giving all z bins equal weight for the determination of the parameters ${A}_{\mu }$, ${B}_{\mu }$, ${C}_{\mu }$, and ${F}_{\mu }$.

In the following we use two approaches in parallel. First, to visualize specific trends in Figures 4 and 5, we use data averages/medians of typically 20–100 individual measurements, separately for each technique. Binning in this way elucidates possible deviations from the assumed basic power-law parameterizations introduced in Equations (5) and (6). Second, for quantitative fitting of the data, we employ multiparameter linear regression fitting in 3-, 4-, or 5-space of logarithmic variables (log(1 + z), log²(1 + z), log δMS, log δM_*, log δR_e) using all data points individually, weighted by the inverse square of their uncertainties. For the specific fitting of log(1+z) versus log t_depl or log ${\mu }_{\mathrm{gas}}$, we explored how the fits changed when we gave equal weight to each of the six redshift bins used in the fit. This is important for establishing the best overall log(1 + z)–log ${\mu }_{\mathrm{gas}}$ scaling relation (Equation (6)) in the presence of the nonlinear fitting function (β = 2). The equal bias for different z bins removes the otherwise overly strong weight of the large number of z ∼ 0 CO data points. The results of the fitting exercises are reported in Table 3(a), including a recommended, overall "best" set of fit parameters (boldface). We determined the uncertainties in Table 3(a) by splitting the sample randomly into two halves, fitting the parameters, and then repeating the splitting and fitting to establish the range of uncertainties by bootstrapping. As a second check, we also did the fitting procedures after eliminating one or more of the smaller data sets (jackknifing), thus checking for systematics of the individual data sets (see also the Appendix).

Zoom In Zoom Out Reset image size

4.2.1. Redshift Dependence of Depletion Timescale

4.2.2. Variations of Depletion Timescale above and below the MS

4.2.3. Mass Dependence

4.2.4. Size Dependence

We now repeat the same exercise for the ratio of molecular gas to stellar mass, ${\mu }_{\mathrm{gas}}={M}_{\mathrm{molgas}}/{M}_{* }$, using Equation (6) as the basis of our fitting. Figure 6 and Table 3(b) show the results. We first fit for the dependence of μ_gas as a function of redshift normalized to the MS reference line of S14, applying the zero-point corrections to the different data sets. The result is shown in the bottom left panel of Figure 6. The decline in μ_gas over cosmic time is consistent with the results found previously in the literature based on fewer data points (see Table 1). We find that fitting the redshift dependence of ${\mu }_{\mathrm{gas}}$ across the entire range requires a curved, second-order function in log(1 + z) –log μ_gas space, and the best fit with β = 2 is given in Table 3(b) and clearly improves the quality of the fit over the simple power law. This second-order fitting function should not be surprising, as it merely reflects the shape of sSFR in Equation (1) and demonstrates that sSFR and μ_gas of MS galaxies track each other, as found previously by Tacconi et al. (2010, 2013) and others. Because of this curvature of μ_gas with local slope getting flatter the higher the mean z, fitting subsets of data with only partial redshift coverage will reflect this systematic flattening, for instance, when only fitting the FIR dust or 1 mm dust data that do not have data at z < 0.4. We also show the mean total gas fraction (H i + H₂) at z = 0 for a mean stellar mass of log(M_*/M_⊙) = 10.7, taken from Saintonge et al. (2011a). As mentioned above, the dust data of Santini et al. (2014), G15, and Berta et al. (2016) in the lowest redshift bin (z ∼ 0.1) are located above best-fit relation, suggesting that dust from the atomic medium in the z = 0 data of Saintonge et al. (2011a, 2017) could be becoming significant (Sections 3 and 3.2). In this interpretation z ∼ 0.1–0.3 would mark the transition between a primarily atomic ISM and a primarily molecular ISM. The mass contribution of the atomic component is approximately constant with redshift (Bauermeister et al. 2010), while the molecular component is strongly evolving with redshift (see also Lagos et al. 2011, 2015a).

Zoom In Zoom Out Reset image size

With the redshift evolution established, we now fit for the dependence with sSFR (δMS) perpendicular to the MS. The top left panel of Figure 6 confirms the tight correlation found previously (e.g., Saintonge et al. 2012, 2016; Tacconi et al. 2013; Sargent et al. 2014, G15) that gas fractions increase with $\delta \mathrm{MS}$ (${\mu }_{\mathrm{gas}}\sim \delta {\mathrm{MS}}^{0.52}$). The increase in sSFR above the MS is then a combination of an increase in the available gas for star formation and an increase in star formation efficiency (decrease in depletion time as described in the previous section).

Next, we fit for the mass dependence after subtracting the dependences on redshift and sSFR and show the result in the bottom right panel of Figure 6. The plot clearly shows the drop in μ_gas at high M_* as found by us and others, at both low (Saintonge et al. 2011b) and high (e.g., Magdis et al. 2012b; Tacconi et al. 2013; G15; Scoville et al. 2017) redshift. Driven by the xCOLD GASS data at log(M_*/M_⊙) < 10 SFGs, there is a tendency for the slope to steepen at the high-mass tail, and requiring a second-order fitting function, log μ_gas = 0.35_0.15 − 0.085_0.05 × (log(M_*/M_⊙) − 8.7₁). As in the discussion of the depletion time dependence on mass, conclusions on the distribution at low mass strongly depend on the metallicity correction of ${\alpha }_{\mathrm{CO}}$ applied. This drop reflects the drop seen in the MS reference curve (W14; Schreiber et al. 2015) and is most likely correlated with the quenching of galaxies beyond the Schechter mass at all redshifts, the so-called "mass quenching" and the formation of bulges (e.g., Peng et al. 2010). Finally, the top right panel of Figure 6 reflects the weak dependence μ_gas on size, as for t_depl, once one has marginalized over the other three parameters. Potential reasons for this lack of a size dependence are discussed in the previous section.