COMBINED CO AND DUST SCALING RELATIONS OF DEPLETION TIME AND MOLECULAR GAS FRACTIONS WITH COSMIC TIME, SPECIFIC STAR-FORMATION RATE, AND STELLAR MASS^*

To explore the cold molecular gas in SFGs covering the entire redshift range from z = 0 to 4, the stellar mass range of M_* = 10^9.8–10^11.8 M_☉, and at a given redshift and stellar mass, SFRs from about 10⁻¹–10² times the main-sequence SFR, we collected 500 CO detections of SFGs near, below, and above the main sequence from a number of concurrent molecular surveys with CO 1–0, 2–1, 3–2 (and in two cases 4–3) rotational line emission (Table 1). We include:

• 1.  
  Two hundred and sixteen detections and one stack detection (much below the main sequence) of CO 1–0 emission above and below the main sequence between z = 0.025 and 0.05 from the final COLDGASS survey with the IRAM 30 m telescope (Saintonge et al. 2011a, 2011b; A. Saintonge et al., in preparation). We note that the SFRs in that survey have been updated from earlier UV-/optical spectral energy distribution (SED) fitting (Saintonge et al. 2011a) with mid-IR SFRs from WISE (A. Saintonge et al., in preparation; Huang & Kauffmann 2014);
• 2.  
  Ninety CO 1–0 detections with the IRAM 30 m of z = 0.002–0.09 luminous and ultra-luminous 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, J. Gracia-Carpio et al. (2014, private communication), and Garcia-Burillo et al. (2012);
• 3.  
  Thirty-one CO 1–0 or 3–2 detections of (above main-sequence) SFGs between z = 0.06 and 0.5 with the CARMA millimeter array from the EGNOG survey (Bauermeister et al. 2013);
• 4.  
  Fourteen 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-main-sequence) ULIRGs with the IRAM 30 m telescope from Combes et al. (2011, 2013);
• 5.  
  Eleven CO 2–1 or 3–2 detections of near main-sequence SFGs between z = 0.5 and 3.2 from Daddi et al. (2010a) and Magdis et al. (2012b), obtained with the IRAM PdBI;
• 6.  
  Six CO 2–1 detections of z = 1–1.2 main-sequence SFGs selected from the Herschel–PACS Evolutionary Probe (PEP) survey (Lutz et al. 2011), obtained with the IRAM PdBI (Magnelli et al. 2012a);
• 7.  
  Fifty-two detections of CO 3–2 emission in main-sequence SFGs in two redshift slices at z = 1–1.5 (38) and z = 2–2.5 (14), as part of the PHIBSS1 survey with the IRAM PdBI (Tacconi et al. 2010, 2013);
• 8.  
  Thirty-one detections (and two upper limits) of CO 2–1 or 3–2 in main-sequence SFGs between z = 0.5 and 1, and 3 at z ∼ 2, as part of the PHIBSS2 survey with the IRAM PdBI (L. J. Tacconi et al., in preparation);
• 9.  
  Nineteen CO 2–1, 3–2, or 4–3 detections of above main-sequence 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.  
  Eight CO 3–2 detections of z = 1.4–3.2 lensed main-sequence SFGs obtained with the IRAM PdBI (Saintonge et al. 2013, and references therein).

The redshift–sSFR coverage of this sample is shown in Figure 1, with the different symbols denoting the various surveys mentioned in our listing. The overall distribution in z–sSFR space is biased to SFGs above the main sequence because of the sensitivity limits. However, the more recent extensive surveys at the IRAM telescopes at z ∼ 0.03 (COLDGASS), z ∼ 0.7 (PHIBSS2), z ∼ 1.2 (PHIBSS1+2), and z = 2.2 (PHIBSS1) have begun to establish a decent coverage of massive SFGs above and below the main-sequence line. Most of these data are benchmark sub-samples of large UV/optical/IR/radio imaging surveys with spectroscopic redshifts, and well-established and relatively homogeneous stellar and star-formation properties. The COLDGASS sample is drawn from SDSS. PHIBSS1+2 and the data of Daddi et al. (2010a), Magdis et al. (2012b), and Magnelli et al. (2012a) are selected from deep rest-frame UV/optical imaging surveys in the Extended Groth Strip (Davis et al. 2007; Newman et al. 2013; Cooper et al. 2012), GOODS N (Giavalisco et al. 2004; Berta et al. 2010), and COSMOS (Scoville et al. 2007; Lilly et al. 2007, 2009), including the recent CANDELS J- and H-band Hubble Space Telescope (HST) imaging (Grogin et al. 2011; Koekemoer et al. 2011) and 3D-HST grism spectroscopy (Brammer et al. 2012; Skelton et al. 2014), as well as D3a (Kong et al. 2006) and the BX/BM UGR color selected samples of Steidel et al. (2004) and Adelberger et al. (2004).

Zoom In Zoom Out Reset image size

We have binned the 500 SFGs of our CO sample into 6 redshift bins (Table 1). The number of SFGs in each of the five highest z bins is comparable (28–49). The highest four bins have a good coverage of the main-sequence population, while the lowest of these non-local bins (z = 0.05–0.45) contains mostly above main-sequence, starburst outliers. There are few galaxies significantly below the main sequence, for the obvious reason of detectability. The lowest redshift bin (mostly COLDGASS) naturally contains the largest number of galaxies (296 of the 500 galaxies). This imbalance needs to be carefully taken into account when considering the scaling relations.

We emphasize that the majority of our final sample of ∼500 galaxies are near-main sequence SFGs Δlog(sSFR/sSFR(ms)) = ±0.6 (dashed horizontal lines in Figure 1), with a few below main sequence (mainly from the COLDGASS sample at z = 0), and ∼130 (26%) above main-sequence starburst outliers. The focus of this paper is on the near-main-sequence population.

2.1.1. Derivation of Molecular Gas Masses

Observations of GMCs in the Milky Way and nearby galaxies have established that the integrated line flux of ¹²CO millimeter rotational lines can be used to infer molecular gas masses, although the CO molecule only makes up a small fraction of the entire gas mass, and its lower rotational lines (1–0, 2–1, 3–2) are almost always very optically thick (Dickman et al. 1986; Solomon et al. 1987; Bolatto et al. 2013). This is because the CO emission comes from moderately dense (volume average densities 〈n(H₂)〉 ∼ 200 cm⁻³, column densities N(H₂) ∼ 10²² cm⁻²), self-gravitating GMCs of kinetic temperature 10–50 K. Dickman et al. (1986) and Solomon et al. (1987) have shown that in this virial regime—or if the emission comes from an ensemble of similar mass, near-virialized clouds spread in velocity by galactic rotation—the integrated line CO line luminosity $L^{\prime} _{{\rm CO}} = \int_{{\rm source}} \int_{{\rm line}} {T_R }\, {dv}\, {dA}$ (in K km s⁻¹ pc², T_R is the Rayleigh–Jeans source brightness temperature as a function of Doppler velocity v) is proportional to the total gas mass in the cloud/galaxy. In this cloud counting technique the total molecular gas mass (including a 36% mass correction for helium) then depends on the observed CO J → J − 1 line flux F_{CO J}, source luminosity distance D_L, redshift z, and observed line wavelength λ_{obs J} = λ_rest J (1 + z) as (Solomon et al. 1997)

$$\begin{eqnarray} M_{{\rm mol\ gas}}/M_ \odot &=& \alpha _{{\rm CO }\,1} \times L^{\prime}_{{\rm CO }\,1} \nonumber\\ &=& 1.75 \times 10^9 \left({\frac{{\alpha _{{\rm CO}\,1} \times R_{1J} }}{{\alpha _{{\rm MW}} }}} \right) \times \left({\frac{{F_{{{\rm CO} \, J}} }}{{{\rm Jy\, km}\,{\rm s^{-1}}}}} \right)\nonumber\\ &&\times ({1 + z})^{ - 3} \times \left({\frac{{\lambda _{{{\rm obs}\, J}} }}{{{\rm mm}}}} \right)^2 \times \left({\frac{{D_L }}{{{\rm Gpc}}}} \right)^2. \end{eqnarray} \tag{ 4 }$$

Here α_{CO 1} is an empirical conversion factor to transform the observed quantity (CO luminosity in the 1–0 transition) to the inferred physical quantity (molecular gas mass), and R_1J is the ratio of the 1–0 to the J – (J – 1) CO line luminosity, R_1J = L'_{CO 1-0}/L'_{CO J – (J-1)}.

From conversion factor to conversion function. From theoretical considerations the CO conversion factor in Equation (4) is expected to be a function of several physical parameters (Narayanan et al. 2011, 2012; Feldmann et al. 2012a, 2012b). In the virial/cloud counting model, α depends on the ratio of the square root of the average cloud density 〈n(H₂)〉 and the equivalent Rayleigh–Jeans brightness temperature T_RJ of the CO transition J → J − 1. It also increases with the inverse of the metallicity Z (see Leroy et al. 2011; Genzel et al. 2012; and Bolatto et al. 2013 for more detailed discussions of the observational evidence):

$$\begin{equation} \alpha _{{{\rm CO}\, J}} = \zeta \left({\frac{{\left({\left\langle {n(H_2)} \right\rangle } \right)^{1/2} }}{{T_{R{ J}} }}} \right) \times \chi (Z). \end{equation} \tag{ 5 }$$

In the Milky Way and nearby SFGs with near solar metallicity, as well as in dense star-forming clumps of lower mass, lower metallicity galaxies, the empirical CO 1–0 conversion factor α_{CO 1} determined with dynamical, dust, and γ-ray calibrations are broadly consistent with a single value of α_{CO 1} = α_MW = 4.36 ± 0.9 (M_☉ (K km s⁻¹ pc²)⁻¹), which is equivalent to X_CO = N(H₂)/(T_{RJ = 1}Δv) = 2 × 10²⁰ (cm⁻² (K km s⁻¹)⁻¹; Strong & Mattox 1996; Dame et al. 2001; Grenier et al. 2005; Bolatto et al. 2008, 2013; Leroy et al. 2011; Abdo et al. 2010; Ostriker et al. 2010).

Metallicity dependence of the conversion factor. For galaxies of sub-solar gas phase metallicity, the conversion factor and metallicity are inversely correlated, as the result of an increasing fraction of the molecular hydrogen gas column that is photo-dissociated, resulting in molecular gas that is deficient (dark) in CO (Wilson 1995; Arimoto et al. 1996; Israel 2000; Wolfire et al. 2010; Leroy et al. 2011; Genzel et al. 2012; Bolatto et al. 2013; Sternberg et al. 2014). Motivated by the theoretical work of Wolfire et al. (2010) on the photo-dissociation of clouds with a range of hydrogen densities and UV radiation field intensities, but with a constant hydrogen column, Bolatto et al. (2013) proposed the following fitting function for χ(Z):

$$\begin{equation} \chi (Z) = 0.67 \times \exp \left(0.36 \times 10^{ - (12 + \log ({\rm O}/{\rm H}) - 8.67)}\right), \end{equation} \tag{ 6 }$$

where 12+log(O/H) is the gas phase oxygen abundance in the galaxy on the Pettini & Pagel (2004) calibration scale, with the solar abundance of 8.67 (Asplund et al. 2004). Equation (6) assumes an average GMC hydrogen column density of 100 M_☉ pc⁻², or 9 × 10²¹ cm⁻². Genzel et al. (2012) combined the local (Leroy et al. 2011) and high-z empirical evidence for a second fitting function,

$$\begin{equation} \chi ({\rm Z}) = 10^{ - 1.27 \times (12 + \log ({\rm O}/{\rm H}) - 8.67)}. \end{equation} \tag{ 7 }$$

For the near-solar metallicities typical for most SFGs in our overall sample (96% of the 1012 SFGs are between 12 + log(O/H) = 8.55 and 8.75 on the PP04 scale), Equations (6) and (7) yield values for χ(Z) within ±0.12 dex of each other. We thus took the geometric mean of Equations (6) and (7) in estimating the gas masses from CO in this paper. Note that this approach is not applicable for significantly sub-solar metallicity galaxies. Between 12+log(O/H) = 7.9 and 8.4, Equation (7) implies a correction 0.22–0.32 dex greater than Equation (6).

CO ladder excitation dependence of the conversion factor. To convert the CO 2–1 and 3–2 luminosities in the near-main-sequence SFGs (at all redshifts) to an equivalent CO 1–0 luminosity, we apply a correction factor of R_1J = 1.3 and 2 to correct for the lower Rayleigh–Jeans brightness temperature of the J – (J−1) relative to the 1–0 transition. This excitation correction entails a combination of the Planck correction (for a finite rotational temperature), as well as a correction for a sub-thermal population in the upper rotational levels. For above main-sequence SMGs and ULIRGs (sSFR/sSFR(ms, z, M_*) > 4) we take R_1J = 1.2, 1.9, and 2.4 for the 2–1, 3–2, and 4–3 transitions. These correction factors are empirically motivated by recent CO ladder observations in low- and high-z SFGs (Weiss et al. 2007; Dannerbauer et al. 2009; Ivison et al. 2011; Riechers et al. 2010; Combes et al. 2013; Bauermeister et al. 2013; Bothwell et al. 2013; Aravena et al. 2014; Daddi et al. 2015). While these correction factors undoubtedly vary from galaxy to galaxy, their scatter is unlikely to be greater than ±0.1 dex, as judged from the recent data sets.

Density–temperature dependence of the conversion factor. This leaves the correction factor/function ζ((〈n(H₂)〉)^1/2/T_RJ), which is correlated with the SFR at a given mass and redshift, that is, the vertical location in the stellar mass–SFR plane (Elbaz et al. 2011; Gracia-Carpio et al. 2011; Nordon et al. 2012; Lada et al. 2012). Our initial assumption is that this function is a constant of order unity. As we show in Section 4.1, this assumption is justified for the z = 0–3 near-main-sequence population, which is at the focus of this paper. However, this assumption is very likely not appropriate for outliers or starbursts high above the main sequence (Daddi et al. 2010b; Genzel et al. 2010, and references therein; Magdis et al. 2012a; Sargent et al. 2012, 2014; Tan et al. 2014). For instance, for extreme z = 0 ULIRGs there is good evidence from dynamical arguments that α_CO is 0.8–1.5, or 0.46–0.74 dex smaller than the Milky Way value, perhaps implying the presence of a second, starburst star-formation mode with ∼5 times greater star-formation efficiency (Scoville et al. 1997; Downes & Solomon 1998). Daddi et al. (2010b), Genzel et al. (2010), and Sargent et al. (2014) incorporated this information for the input assumptions of Equation (5). However, this approach has since come under criticism (e.g., Ostriker & Shetty 2011; Kennicutt & Evans 2012; Krumholz et al. 2012) on the grounds that the resulting smaller depletion timescales then immediately introduce a bi-modal, and thus probably unphysical gas-star-formation relation. In this paper we do not a priori include such a correction factor for the outlier population, but then derive in Section 4.1 quantitative constraints on the scaling of α_CO with sSFR/sSFR(ms, z, M_­*) from the comparison of the CO data to the dust data (not affected by the conversion factor issue). From this comparison we will show that for the near-main population that is the focus of this paper, this simplified approach delivers a good description of the conversion factor.

Our starting point in this paper thus is to use

$$\begin{equation} \alpha _{0J} = \alpha _{{\rm MW}} \times \chi (Z) \times R_{1{ J}} \end{equation} \tag{ 8 }$$

to derive molecular gas masses from CO observations for all 500 SFGs.

As part of the Herschel–PEP (Lutz et al. 2011) and Herschel-HERMES (Oliver et al. 2012) far-IR continuum surveys, Magnelli et al. (2014) established 100–500 μm far-IR SEDs by stacking PACS and SPIRE photometry in 8846, 4753, and 254,749 K- and I-selected SFGs in the GOODS-N, GOODS-S, and COSMOS fields, respectively. For details of the methodology we refer to Magnelli et al. (2014). Briefly, SFRs are calibrated onto the Wuyts et al. (2011a) ladder of UV-, mid-IR, and far-IR based indicators. Since the far-IR detection rate drops with increasing redshift and decreasing SFR and M_*, it is necessary to average many individual data points to determine good far-IR SEDs as a function of z, SFR, and M_*. For this purpose, Magnelli et al. (2014) binned the data onto a three-dimensional grid in z, SFR, and M_*, and then stacked the photometry in each bin. Next Magnelli et al. determined the dust temperature for each resulting SED by fitting to model SEDs from the library of Dale & Helou (2002), for which dust temperatures were established from single optically thin, modified blackbody fits with emission index β = 1.5 (Table A.1 in Magnelli et al. 2014). To ensure constrained SED shapes, we only use bins where the stacked photometry is detected at >3σ in at least three bands that encompass the fitted SED peak, and the reduced χ² of the fit is less than two. In the median, detections in our stacked SEDs reach out to a rest wavelength of 223 μm, and only ∼10% stop at ⩽160 μm.

From these stacks we computed dust masses using Draine & Li (2007) models. We fitted the models following the procedure prescribed by these authors, and widely used in the literature, adopting the Li & Draine (2001) values of dust opacity as a function of wavelength. We limited the parameter space to the range suggested by Draine et al. (2007) for galaxies missing submillimeter data, based on their analysis of local SINGS galaxies. Draine et al. (2007) compared dust masses of local SINGS galaxies obtained using SEDs sampled out to rest-frame ∼160 μm and SEDs that also include the submillimeter. They found that in the absence of submillimeter data, dust masses are a factor ⩽2 more uncertain but exhibit no net bias. Magdis et al. (2012a) confirm similar results for a small sample of high-redshift galaxies detected in the submillimeter, excluding any data points at wavelength >200 μm. Based on a Monte Carlo analysis, the errors on the stacked photometry correspond to a median uncertainty of 0.14 dex for our dust masses. S. Berta et al. (in preparation) present a more comprehensive analysis of uncertainties in Herschel-based dust masses. In total, we obtain Draine & Li (2007) dust masses and Magnelli et al. (2014) dust temperatures for 512 bins in z, SFR, and M_*

As recognized by several authors, there can be systematic differences between Draine & Li (2007)-based dust masses and the typically smaller ones derived using single-temperature modified blackbody models, with details depending on the treatment of dust opacities and emissivity indices (e.g., Magnelli et al. 2012a; Magdis et al. 2012a; Bianchi 2013; S. Berta et al. in preparation), as well as between dust temperatures that are based on different conventions. We defer a detailed discussion to S. Berta et al. (in preparation), but note that Equation (9) below and our subsequent results refer to dust masses from the Draine & Li (2007) method, dust temperatures as in Magnelli et al. (2014), and SFRs based on the Wuyts et al. (2011a) ladder. For all bins, these SFRs agree within ±0.3 dex with the IR luminosities from the stacks (Magnelli et al. 2014).

If dust is a calorimeter of the incident UV, radiating at an average temperature and optically thin in the far-IR, a simple scaling is expected between M_dust, T_dust, and SFR. For our adopted scales and an emissivity index β = 1.5, this relation takes the form

$$\begin{equation} \left({\frac{{M_{\rm dust} }}{{M_ \odot }}} \right) = 1.2 \times 10^{15} \left({\frac{{\rm SFR}}{{M_ \odot\, {\rm yr}^{ - 1} }}} \right) \times \left({\frac{{T_{\rm dust} ({\rm MBB})}}{K}} \right)^{ - 5.5}. \end{equation} \tag{ 9 }$$

where the constant has been calibrated on the data for the 512 bins.

The conversion to gas masses requires the application of a metallicity dependent dust-to-gas ratio correction, which also enters the redshift evolution through the redshift dependence of the mass–metallicity relation below (e.g., Bethermin et al. 2015). Following Magdis et al. (2012a) and Magnelli et al. (2012a) we converted the Draine & Li (2007) model dust masses to (molecular) gas masses by applying the metallicity dependent dust-to-gas ratio fitting function for z ∼ 0 SFGs found by Leroy et al. (2011):

$$\begin{equation} \delta _{dg} = \frac{{M_{{\rm dust}} }}{{M_{{\rm mol\ gas}} }} = 10^{(- 2 + 0.85 \times (12 + \log ({\rm O}/{\rm H}) - 8.67))}, \end{equation} \tag{ 10 }$$

where 12+log(O/H) is the gas phase oxygen abundance (see also Draine et al. 2007 for dust to gas with metallicity scalings of the SINGS nearby galaxy sample, and Galametz et al. (2011) or Rémy-Ruyer et al. (2014) for lower metallicity galaxies down to 12+log(O/H) = 8.0). We note that the metallicity dependence in Equation (10) is within a few percent of that found in the last section from averaging Equations (6) and (7). This means that the metallicity (and hence, mass) corrections we choose in this paper for the dust and CO data are very similar.

As in the case of the CO sample, the 512 stacks are grouped into six redshift bins comparable to those of the CO sample. These 512 stacks provide complete and unbiased estimates of the mean far-IR/submillimeter properties of all SFGs with 0.16 < z < 2, M_* > 10¹⁰ M_☉, and log(sSFR/sSFR(ms, z, M_*)) > −0.3 (see Figures 4 and 5 of Magnelli et al. 2014). In contrast to the CO sample, the dust sample has comparable numbers of SFGs (83–191) in the middle four redshift bins, while the number in the lowest and highest bins are significantly smaller (∼30 each), introducing substantially greater uncertainties in the redshift scaling relation for the dust data as compared to the CO data. This difference actually turns out to be advantageous in the discussion of the scaling relations below, as the dust sample is obviously not dominated in number by the lowest redshift bin.

For the few SFGs in this paper with estimates of gas phase metallicities from strong line rest-frame optical line ratios, we determine individual estimates of logZ = 12 + log (O/H). For instance, if the λ6583 [N ii]/λ6563 Hα line flux ratio is measured, the Pettini & Pagel (2004) indicator yields

$$\begin{eqnarray} &&12 + \log ({\rm O/H})_{{\rm PP04}} = 8.9 + 0.57\nonumber\\ &&\quad \times \log (F(6583[{\rm N}\,{\scriptsize{II}}])/F(6563{\rm H\alpha})). \end{eqnarray} \tag{ 11 }$$

The scatter in the above relation is ±0.18 dex. However, for the large majority of the SFGs in our CO and dust samples, such line ratios are not available and it is necessary, for the metallicity corrections discussed above, to refer to the mass–metallicity relation. Following Maiolino et al. (2008) we combined the mass–metallicity relations at different redshifts presented by Erb et al. (2006), Maiolino et al. (2008), Zahid et al. (2014), and Wuyts et al. (2014) in the following fitting function:

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

Mannucci et al. (2010) presented evidence for a dependence of metallicity on the SFR for z ∼ 0 SDSS galaxies at a given stellar mass (the fundamental metallicity relation), yielding an alternative version of Equation (12a):

$$\begin{eqnarray} &&\hspace{-21pt}\Delta Z_{{\rm M08}} = 12 + \log ({\rm O/H})_{{\rm M08}} - 8.69\nonumber\\ &&\hspace{-21pt}\,\qquad= 0.21 + 0.39 \times x - 0.2 \times x^2\nonumber\\ &&\hspace{-21pt}\qquad\quad\,- 0.077 \times x^3 + 0.064 \times x^4, \nonumber\\ &&\hspace{-21pt}{\rm with }\quad x = \log M_* - 0.32 \times \log SFR - 10,\quad {\rm and} \nonumber\\ &&\hspace{-21pt} 12 + \log ({\rm O/H})_{{\rm PP04}} - 8.9 = - 0.4408 + 0.7044\nonumber\\ &&\hspace{-21pt}\quad\qquad\times \Delta Z_{{\rm M08}} - 0.1602 \times (\Delta Z_{{\rm M08}})^2 \nonumber\\ &&\hspace{-21pt}\quad\qquad- 0.4105 \times (\Delta Z_{{\rm M08}})^3 - 0.1898 \times (\Delta Z_{{\rm M08}})^4, \end{eqnarray} \tag{ 12b }$$

where stellar masses are in solar masses and SFRs in solar masses per year. The Mannucci et al. relation (their Equation (4)) is on the Maiolino et al. (2008, M08) scale, which is then converted to the PP04 scale using the coefficients in their Table 4.

This relation implies that at constant stellar mass and within ±0.6 dex of the main-sequence, the PP04 metallicity changes by ±0.04 dex near solar metallicity, which is a second-order correction for Equations (6) and (7), given the uncertainties in metallicity and SFRs (even for SFGs far from the main sequence the correction is only ∼−0.1 dex). At high-z several recent studies of strong emission line indicators also indicate a weak dependence of metallicity on SFR for massive (log(M_*/M_☉) > 10) galaxies (Steidel et al. 2014; Wuyts et al. 2014; Sanders et al. 2015). Wuyts et al. (2014) found that the fundamental metallicity relation in Equation (12b) does broadly trace the redshift evolution of the mass–metallicity relation in Equation (12a). At constant z Wuyts et al. find that for a 0.6 dex change in SFR, the implied metallicity does not change more than 0.08 dex, at 〈z〉 = 0.9 and 2.3. While the application of the strong emission line indicators to metallicity determinations at high z is currently under debate (e.g., Steidel et al. 2014), at face value these results imply a negligible correction in Equations (6) and (7) when applying Equation (12b). For these reasons our default assumption in the following is Equation (12a). We will discuss in Section 4.2 how replacing Equation (12a) by Equation (12b) quantitatively affects the scaling relations.

For the COLDGASS sample the stellar masses and luminosities are calibrated in the frame of SDSS, Galaxy Evolution Explorer, and WISE (Saintonge et al. 2011a, 2012). For the various LIRG/ULIRG samples at z ∼ 0–1 we refer to the original papers for a discussion of the stellar masses, which were converted, if necessary, to the Chabrier IMF adopted here. Infrared luminosities were obtained from the far-IR (30–300 μm) SEDs, assuming L_IR = 1.3 × L_FIR, and SFRs were estimated from Kennicutt (1998) with a correction to the Chabrier IMF adopted here, SFR (M_☉ yr⁻¹) = 10⁻¹⁰ × L_IR (L_☉) (see discussion in Genzel et al. 2010).

At z ⩾ 0.1 global stellar properties for all optically/UV-selected SFGs (both for the CO and dust samples) were derived following the ladder of indicators procedure as outlined by Wuyts et al. (2011a). In brief, stellar masses were obtained from fitting the rest-UV to near-IR SEDs with Bruzual & Charlot (2003) population synthesis models, the Calzetti et al. (2000) reddening law, a solar metallicity, and a range of star-formation histories (in particular including constant SFR, as well as exponentially declining or increasing SFRs with varying e-folding timescales). SFRs were obtained from rest-UV+IR luminosities through the Herschel–Spitzer-calibrated ladder of SFR indicators of Wuyts et al. (2011a) or, if not available, from the UV–optical SED fits. For the main-sequence population (with near-constant star-formation histories) we adopt uncertainties of ±0.15 dex for the stellar masses, and ±0.2 dex for the SFRs, although somewhat smaller uncertainties may be appropriate for SFGs with measurements of individual far-IR luminosities (Wuyts et al. 2011a).

For SMGs we adopted the stellar masses and luminosities of Magnelli et al. (2012b, 2014) and B. Magnelli et al. (2014, private communication), the latter being derived from PACS/SPIRE Herschel SEDs and converted to SFRs with the modified Kennicutt (1998) conversion as given above. The stellar masses and SFRs of most above main-sequence outliers (ULIRGs, SMGs) are more uncertain than those of the main-sequence populations (±0.3 dex). The outliers are more dusty (e.g., Wuyts et al. 2011b), making population synthesis analysis of the UV/optical rest-frame SEDs more uncertain. The bursty nature of the star-formation histories adds substantial additional uncertainties to the average SFRs (e.g., Figure 6 in Genzel et al. 2010), as well as to the inferred stellar masses, as is demonstrated by the up to one order of magnitude varying estimates of stellar masses in SMGs in different publications (Hainline et al. 2009, 2011; Michalowski et al. 2012, 2014; Davé et al. 2010; Bussmann et al. 2012). In the specific case of the stellar masses for the GOALS LIRG/ULIRG sample used in this paper (Howell et al. 2010), this may result in an overestimate of the intrinsic stellar masses. Active galactic nucleus may contribute to the bolometric luminosity for the extreme starburst population (e.g., for z ∼ 0 ULIRGs; Genzel et al. 1998); this means that SFRs in these systems may be overestimated (by 0.1–0.4 dex; see Genzel et al. 2010, and references therein).

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.

In the left panel of Figure 2 we plot for the six redshift bins (different colored symbols) the CO-based depletion timescale as a function of normalized sSFR offset from the main-sequence line at a given redshift (sSFR(ms, z, M_*), Equation (1)). In this log–log presentation log(t_depl) appears to scale linearly with log(sSFR/sSFR(ms)) over more than three orders of magnitude in sSFR, from more than a factor of 10 below to two orders of magnitude above the main sequence, and even in the regime of the extreme outliers, such as z ∼ 0–0.5 ULIRGs and some SMGs. This means that the dependence of t_depl on sSFR is fit by a single s law to within the uncertainties dictated by the scatter of the relation. We note that the conclusion of a single power law is strictly correct only if α_CO does not vary significantly with log(sSFR/sSFR(ms)), as has been assumed implicitly in Equation (8). We explore the justification of this assumption for the near-main-sequence SFGs quantitatively in Section 4.1. For the extreme above main-sequence outlier/starburst ULIRG population at z ∼ 0 (log(sSFR/sSFR(ms)) > 1) there is persuasive evidence from mass modeling that α_CO is four to five times lower (Scoville et al. 1997; Downes & Solomon 1998; Daddi et al. 2010b; Genzel et al. 2010). If this correction is applied to the data in Figure 2, the resulting log(t_depl)–log(sSFR/sSFR(ms)) distribution would show a downward kink and no longer be fit by a single power law. This may imply a second starburst mode of star formation (Sargent et al. 2012, 2014; Magdis et al. 2012a).

Zoom In Zoom Out Reset image size

Fitting a power law to each of the t_depl–(sSFR/sSFR(ms, z, M_*)) z-bins shows no significant redshift evolution of the slope ξ_g1(z) = d(log t_depl(z)/dlog (sSFR/sSFR(ms, z, M_*)) (Table 2, Figure 3). A weighted fit to the CO slope data in Figure 3 (filled blue circles) yields dξ_g1(z)/dlog(1+z)) = −0.08 (±0.13, 1σ). We note that there is no evidence for a z-dependence of the dust-based depletion timescales (black filled circles in Figure 3) discussed below (Section 3.2), for which a weighted fit to the six redshift bins yields dξ_g1(z)/dlog(1+z) = +0.13 (±0.18, 1σ).

Zoom In Zoom Out Reset image size

A similar analysis in the logt_depl–log M_* shows that a linear function (a power law in the original variables) with slope 0 (±0.1) can account for the data in each of the redshift bins.

These are important constraints. If a function of three independent variables (z, sSFR/sSFR(ms, z, M_*), M_*) is a power law in each of these variables, and each has a slope that does not depend on the other variable, then the function can be written as a product of power-law functions, with each dependent only on one variable. That means that the variables can be separated:

$$\begin{eqnarray} && t_{\rm depl} (z, {\rm sSFR}, M_*)|_{\alpha = \alpha _{0j} } = f_1 (z)|_{{\rm sSFR} = {\rm sSFR} ({\rm ms},z,M_*)} \nonumber\\ && \quad \times g_1 ({{\rm sSFR/sSFR}} ({\rm ms},z,M_*)) \times h_1 (M_*){\rm.} \end{eqnarray} \tag{ 13 }$$

Here f₁(z) tracks the dependence of t_depl on redshift at the main-sequence line (Equation (1)), g₁(sSFR/sSFR(ms, z, M_*)) describes the dependence of t_depl on sSFR relative to the main-sequence line, and h₁(M_*) delineates the stellar mass dependence. Now at first glance, the way we have written Equation (13) (and analyzed the data) might seem a contradiction of the statement on variable separation above, since g₁(sSFR/sSFR(ms, z, M_*) contains Equation (1) in the denominator, which is a function of both z and M_*. However, Equation (1) is a product of power laws, so the dependence of sSFR(ms, z, M_*) can be easily pulled out of g₁. Equation (13) can then be resorted into a product of power laws of the individual variables (z, sSFR, M_*), as needed for separation. The slope of g₁ remains unaffected by this renormalization. It turns out that because of the shallow mass dependence of Equation (13), the mass dependence also does not change much when resorting g₁ as a function of sSFR only. The only function strongly affected is f₁(z) because that now acquires a strong redshift dependence from Equation (1). This discussion already suggests that the separation of variables and the general conclusions on the quality of fits are independent of the choice of the main-sequence prescription, which we will discuss in more detail in Section 4.2. The depletion timescale may also depend on other parameters, such as bulge mass, gas volume, surface density, and environmental density around the galaxy, but these cannot be explored with the current data sets.

If the parameter dependencies of t_depl can be separated according to Equation (13), Equation (3) shows that the parameter dependencies of molecular gas fractions can be separated as well:

$$\begin{eqnarray} && \frac{{\rm M_{\rm mol\ gas} }}{{M_* }}(z, {\rm sSFR}, M_*)|_{\alpha = \alpha _{0J} } = t_{\rm depl} \times {\rm sSFR} \nonumber\\ && \quad = f_2 (z)|_{{\rm sSFR} = {\rm sSFR} ({\rm ms},z)} \times g_2 ({ {\rm sSFR/sSFR}} ({\rm ms},z,M_*)) \nonumber\\ && \qquad \times h_2 (M_*){\rm.} \end{eqnarray} \tag{ 14 }$$

We thus assume log (g₁(sSFR/sSFR(ms, z, M_*)) = a_g1 + ξ_g1 × log (sSFR/sSFR(ms, z, M_*)) (cf. Saintonge et al. 2011b, 2012). In the first iteration we took the best-fit constant slope from Figure 3 (ξ_g1 = −0.46) and fitted the data in each redshift bin for the zero offset a_g1(z). The resulting z-distribution of the zero offsets is shown in the right panel of Figure 2. Their z-dependence can be well described by a power law, log (f₁(z)) = a_f1 + ξ_f1 × log (1 + z). In the second iteration we removed the fitted zero-point offsets as a function of redshift by dividing the original data by f₁(z) and fit the sSFR dependences with a power-law function, as above, but this time for all 500 data points simultaneously. This has the obvious advantage of giving a much more robust estimate of the slope ξ_g1. It also strengthens our assumption that the data set can be fit by a single power law. The result is shown in the right panel of Figure 2 and the left panel of Figure 4. We find that the best-fitting parameters are a_f1 = −0.043 (±0.01), ξ_f1 = −0.16 (±0.04), and a_g1 = 0, ξ_g1 = −0.46 (±0.03). The numbers in the parentheses are the statistical 1σ fit uncertainties only; they do not include uncertainties due to systematics and cross-terms in the co-variance matrix.

Zoom In Zoom Out Reset image size

The redshift dependence is shallower than that found by Tacconi et al. (2013; ξ_f1 = −0.7 to −1). This can be explained by the fact that we now quote the redshift dependence of the depletion time on the main-sequence line. Tacconi et al. (2013) instead took an average of the COLDGASS data (〈t_depl〉_COLDGASS (z = 0) = 1.5 Gyr) and the z = 1.2 PHIBSS data (〈t_depl〉_PHIBSS (z = 1.2) = 0.6–0.7 Gyr). In that case the logarithmic slope is −1. The COLDGASS sample has many massive galaxies below the main-sequence line, whereas PHIBSS1 has mostly galaxies above, such that the average depletion times are biased high (at z = 0) and low (at z = 1.2) because the negative slope in the left inset of Figure 2 results in an overestimate of the redshift dependence. The slope in the log t_depl–log sSFR/sSFR(ms, z, M_*) plane inferred above is in agreement with Saintonge et al. (2012) and Huang & Kauffmann (2014) from COLDGASS data alone.

One might be concerned that the observed correlation between sSFR/sSFR(ms, z, M_*) and t_depl is caused (or at least affected) by the fact that the sSFR is proportional, and the depletion time is inversely proportional to SFR, introducing an artificial negative correlation if the SFR has a substantial uncertainty. To explore the impact of this effect, we created Monte Carlo mock data sets. We started with the scaling relations in Table 3 (t_depl as a function of sSFR/sSFR(ms, z, M_*)) for a mock data set spanning an order of magnitude in stellar mass and centered around some redshift, and then added ±0.2 (respectively ±0.3 dex for outliers) scatter in SFR. We find that the artificial anti-correlation between sSFR and t_depl is only significant if the data have near constant sSFR. For data with a range in sSFR comparable to our observed sample (2.5–3.5 dex), the effect merely leads to a very small increase in scatter, but does not change the slope of the intrinsic relation.

The dispersion of the data in the left panel of Figure 4 around the best-fitting power-law function is ±0.24 dex. This is quite tight given the uncertainties of stellar masses (±0.15–0.3 dex), SFRs (±0.2–0.3 dex), and molecular gas masses (>±0.2 dex), and considering the possibility of substantial variations of the CO conversion factor across the more than three orders of magnitude sSFR/sSFR(ms, z, M_*) variation spanned by the data in Figure 4. There is a tendency for the log(t_depl/f₁) residuals as a function of log(sSFR/sSFR(ms, z, M_*)) to exhibit an excess of negative values for log(sSFR/sSFR(ms, z, M_*)) > 0.6. This may indicate that the depletion timescale above the main sequence, in the starburst-outlier regime, drops faster than that captured by the power-law fit above.

The right panel of Figure 4 explores whether the residuals log(t_depl/(f₁ × g₁)) depend on the remaining internal parameter, stellar mass. There appears to be no significant trend, which is in excellent agreement with the findings of Huang & Kauffmann (2014) at z ∼ 0.

Next we repeat the same exercise for the dust-based depletion time estimates from Herschel. The left panel of Figure 5 again shows the depletion time measurements as a function of sSFR/sSFR(ms, z, M_*) in the six redshift bins. As for CO, the dust-based depletion timescales do not exhibit a significant redshift evolution of ξ_g1(z) = d(log t_depl(z)/dlog (sSFR/sSFR(ms, z, M_*)) (black filled circles in Figure 3), so a constant slope ξ_g1 = −0.59 is an adequate description of the current dust data. The dependences of t_depl on redshift, sSFR, and stellar mass can again be written as a product of power laws, as in Equations (13) and (14). Proceeding as before, we determine the zero points of these power-law functions in each bin, and plot them as a function of redshift in the right panel of Figure 5 (black filled circles), along with the zero points of the CO-based determinations from Figure 2 (blue filled circles).

Zoom In Zoom Out Reset image size

These totally independent estimates are remarkably close, especially given the possibly hidden systematic uncertainties—in CO conversion factor on the one hand and dust modeling and conversion from dust to gas on the other. The dust-based depletion time appears to vary faster with redshift (ξ_f1 = −0.77 (±0.19) than the CO data (ξ_f1 = −0.16 (±0.04)). The difference is only moderately significant (∼3σ), because the dust samples in the lowest and highest redshift bins only contain two dozen data points and the z-coverage is smaller than in the CO data (0.16 < z < 2). The extrapolated (z = 0) zero point of the dust data (a_f1 = 0.34 (±0.07)) is larger than that of the CO data (a_f1 = −0.04 (±0.01)). However, on average between z = 0 and 2.5 the CO- and dust-based depletion time estimates do not differ by more than ∼30%. We note that an even better and more straightforward comparison would be the direct comparison of gas masses determined in the same galaxies from each of the two techniques. This route is not possible, however, because the two samples do not significantly overlap, and the dust technique involves stacking a number of galaxies, rather than observing individual galaxies.

The similarity between CO- and dust-based scaling relations is further strengthened when comparing the redshift-corrected dependences of depletion timescale on sSFR offset, for CO (blue circles) and dust data (red squares) in the left and right panels of Figure 6. The slope of the dust data in the sSFR scaling relation is somewhat steeper than that of the CO data (dust: ξ_g1 = −0.59 (±0.05), CO: ξ_g1 = −0.46 (±0.03)), but the overall distributions overlap over almost three orders of magnitude in sSFR/sSFR(ms, z, M_*), from below the main sequence to the starburst outliers. The lack of a trend as a function of stellar mass also agrees in dust and gas (right panel of Figure 6). This excellent agreement of the separated scaling relations is also particularly relevant because of the very different redshift distributions of the CO and dust data in terms of numbers of SFGs per redshift bin. While the CO data in Figure 6 are heavily weighted toward the z ∼ 0 COLDGASS measurements, the dust data are strongly weighted to z ∼ 1. The agreement in scalings with sSFR and M_* thus cannot be an artifact of biased redshift distributions.

Zoom In Zoom Out Reset image size

Table 3 summarizes the fit parameters for the power-law fitting functions for the CO and dust data. It also gives the parameters for an average between the gas and dust relations that might currently be considered the best description of the molecular gas depletion timescaling relations.

3.2.1. Global Fits and Error Estimates

Next we determined the equivalent relations for the molecular gas fractions. As discussed in the Introduction, once the scaling relations for t_depl are established, those for M_{mol gas}/M_* follow straightforwardly from Equations (1) and (3). Given the slopes of the best-fitting power laws in Table 3, one would then expect the following for the equivalent power-law M_{mol gas}/M_*-fitting function from Equations (1), (3), and (14):

$$\begin{eqnarray} &&\hspace{-19pt} {\rm log}(M_{\rm mol\ gas} /M_* (z,{\rm sSFR},M_*)|_{\alpha = \alpha _{0J} }) \nonumber\\ &&\hspace{-19pt}\quad {\rm } = \log (f_2 (z)|_{sSFR = sSFR({\rm ms},z,M_*)}) \nonumber\\ &&\hspace{-19pt}\qquad + \log (g_2 ({\rm sSFR/sSFR}({\rm ms},z,M_*))) {+} \log (h_2 (M_*)) \nonumber\\ &&\hspace{-19pt}\quad {\rm } = a{}_{f2} + {\rm }\xi _{f2}\times \log (1 + z) + \xi _{g2} \nonumber\\ &&\hspace{-19pt}\qquad\times \log ({\rm sSFR/sSFR}({\rm ms},z,M_*)) + \xi _{h2} \times \log (M_*), \end{eqnarray} \tag{ 15 }$$

with slopes ξ_f2 = ξ_f1 + 3.2 ∼ + 2.9, ξ_g2 = ξ_g1 + 1 ∼ +0.5, and ξ_h2 = ξ_h1 − 0.3 ∼ −0.3.

To check for systematic effects, we decided to check the consistency of the scaling relations (together with Equation (1)) for the CO and dust samples by computing M_{mol gas}/M_* for each galaxy (or galaxy stack) and then establishing the scaling relations. Figures 7–10 show this fitting carried out for the CO data (Figures 7 and 8) and the dust data (Figures 9 and 10).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Consistent with the results of the last section, the agreement between the CO and dust data is very good. We fitted the data with the global fit method described in Section 3.2.1., establishing that both the best-fit values and their uncertainties are robust and well captured by the most probable errors (statistical plus systematic) of the underlying parameters. The parameters of the best-fitting power-law functions are summarized in Table 4. We find ξ_f2 = 2.7, ξ_g2 = 0.5, and ξ_h2 = −0.4 for both the binned and the global fitting methods.

These slopes (and also the zero points) are within ∼0.1 dex of the expectations from Equation (3), and give a measure of the internal systematic uncertainties. We will come back to this topic in Section 4.2.

The depletion times and molecular gas to stellar mass ratios were derived from two independent and very different techniques (CO and dust). The ∼500 measurements for each of the two methods covered the redshift range from 0 to 2–3, which is the range in the sSFR from below to much above the main-sequence line at each redshift, and stellar mass from 10¹⁰ to several 10¹¹ M_☉ yielding similar zero points and, to within the uncertainties, the same scaling indices. A more direct cross-check of the two techniques through a direct, galaxy-by-galaxy comparison would be highly desirable, but is currently not possible because of the lack of sample overlap. The good agreement is better than we would have expected (but see Magdis et al. 2012a). Given the multiple parameter dependences of the CO and dust techniques, one might easily have predicted offsets and trends between these techniques of 0.2–0.5 dex. We only corrected the CO mass estimates for excitation and metallicity dependences, and the dust-to-gas-mass ratios for metallicity dependence. Given the systematic parameter dependencies and uncertainties of the calibrations used in each of the two techniques, the good agreement (in zero point and logarithmic scaling indices) is gratifying.

On this basis, our analysis yields the following main results:

• 1.  
  To first order, the dependences on redshift, sSFR, and stellar mass can be well separated into a product of three power-law functions depending on the individual parameters.
• 2.  
  The depletion time on the main-sequence line is smaller than the Hubble time at all z, and changes only slowly with redshift (t_depl∝ (1+z)^−0.3±0.15, by a factor of 0.7 from z = 0 to z = 2.5), which is in broad agreement with our earlier study in Tacconi et al. (2010, 2013) and less rapid than found by Santini et al. (2014). The molecular gas to stellar mass ratios and the sSFRs as a function of redshift track each other fairly closely (∝ (1+z)³), again in broad agreement with Tacconi et al. (2013), Magdis et al. (2012a), Saintonge et al. (2013), Santini et al. (2014), and Sargent et al. (2014). This finding in turn suggests that the factor of ∼20–30 increase in galactic SFRs between the local universe and the peak of cosmic SFR at z ∼ 1–3 is mainly driven by the increased supply rate of fresh gas, rather than changes in galaxy scale star-formation efficiency (in starbursts with small t_depl). This is consistent with the "gas-regulator" model (Bouché et al. 2010; Tacconi et al. 2010, 2013; Daddi et al. 2010a; Davé et al. 2012; Lilly et al. 2013).
• 3.  
  Changes in the sSFR at constant z and M_* are due to a combination of variations in gas fraction and depletion timescale throughout the redshift range probed, which is in agreement with the z ∼ 0 COLDGASS results of Saintonge et al. (2012), Magdis et al. (2012a), and Huang & Kauffmann (2014). Galaxies above the main sequence have larger gas fractions but also smaller depletion timescales (or equivalently, higher star-formation efficiency), in approximately equal measure, than galaxies at or below the main sequence. The dependence on gas fraction may reflect the time variation in the average gas supply rate from the cosmic web. The increase in star-formation efficiency with sSFR (by a factor of 20 from the lower envelope of the main sequence to the star-bursting outliers above the main sequence) may be driven by the internal properties of the star-forming interstellar medium (ISM), such as the dense gas fraction (Gao & Solomon 2004; Lada et al. 2012; Juneau et al. 2009; Gracia-Carpio et al. 2011; Elbaz et al. 2011) in the more compressed, cuspier SFGs above the main sequence (Wuyts et al. 2011b; Elbaz et al. 2011). The increasing depletion times below the main sequence, especially at high masses, may also be a signature of suppression of the gravitational instability by large shear velocities driving up the Toomre Q-parameter of the ISM above the critical value (morphological quenching; Martig et al. 2009; Genzel et al. 2014; Tacchella et al. 2015).
• 4.  
  Throughout the redshift range probed, the molecular gas to stellar mass ratios decrease as a function of stellar mass (M_gas/M_* ∝ $M_*^{-0.4}$) but the depletion timescale does not vary with stellar mass, which is in agreement with Tacconi et al. (2013), Santini et al. (2014), Huang & Kauffmann (2014), and Sargent et al. (2014). This is in contrast to the ideal gas regular model (d(sSFR)/d(logM_*)|_{main sequence} ∼ 0; Lilly et al. 2013), for which no or only little dependence of gas fractions on stellar mass would be expected. This means that the dropping gas fractions of the observed SFGs at and above the Schechter mass (log(M_S/M_☉) ∼ 10.9) are a direct consequence of the fact that the observed SFGs on the main sequence have lower sSFRs than expected for an ideal gas regulator (sSFR = constant). We interpret these findings as empirical evidence for the expected quenching process(es) that are theoretically expected to happen near and above the Schechter mass.

Our results so far are based on the tenet that the ratio of molecular gas mass to CO luminosity only needs a correction for metallicity and ladder excitation (Equation (8)); α_CO is assumed to not vary with z, sSFR, and M_*. This approach was explained in Section 2.1.1. We will now return to the issue of variations of the CO conversion factor/function as a function of these parameters.

In this section we will take advantage of the fact that the dust data are not dependent on a conversion factor, although they do depend on other uncertain prescriptions for the metallicity dependence of the gas-to-dust ratio, as well as for the mass–metallicity relation and the Draine & Li (2007) dust emissivity modeling. As such, we can plausibly assume that the dust data can be treated as ground truth for the impact of α_CO variations.

A sensitive, straightforward, and robust test of the dependence of α_CO on sSFR, z, and M_* can be derived from the well-determined dependence of dust temperature on sSFR/sSFR(ms, z, M_*) in the Herschel dust data of Magnelli et al. (2014; see also Magdis et al. 2012a; Santini et al. 2014). Of course, such a test can also be made by comparing the dependence of depletion time estimates as a function of sSFR/sSFR(ms, z, M_*) in Figure 6 (with very similar results), but the test on the dust temperature distribution is more robust, as the latter does not depend implicitly on the specific dust model (i.e., Draine & Li 2007), or the prescriptions for the dust-to-gas ratio and stellar mass as function of metallicity. As we will see, the dust-to-gas ratio prescription appears explicitly in the equations, and the dust modeling enters only in the zero point of the dust mass estimates (the constant in Equation (9)).

Following Magnelli et al. (2014) and Magdis et al. (2012a) the relation between the molecular gas depletion timescale and dust temperature in the limit of optically thin far-IR dust emission (and optically thick dust absorption in the UV/optical, dust as a calorimeter) is given by (see Equation (9))

$$\begin{eqnarray} \frac{{L_{{\rm IR}} }}{{M_{\rm dust} }} &=& c_1 \times \frac{{\rm SFR}}{{M_{\rm mol\ gas} \times \delta _{dg} (Z)}} = c_1 \times \frac{1}{{t_{\rm depl} \times \delta _{\rm dg} (Z)}} \nonumber\\ & = & c_2 \times T_{\rm dust}^{4 + \beta }, \end{eqnarray} \tag{ 16 }$$

where the metallicity dependent dust-to-molecular gas ratio δ_dg is given in Equation (10), β is the logarithmic scaling index of the frequency dependence of dust emissivity (β ∼ 1.5; Magnelli et al. 2014), and c₁ and c₂ (as well as c₃, c₄, c₅ below) are constants. We now expand the CO-conversion function (for J = 1) as a function of the four main parameters to the linear terms in the log:

$$\begin{eqnarray} &&\log (\alpha _{\rm CO 1}) = \log (\alpha _{\rm MW}) + \log (\alpha _{{\rm CO}{\rm }1{\rm }0} /\alpha _{\rm MW}) \nonumber\\ && \hphantom{\log (\alpha _{\rm CO 1}) =} + \varepsilon _z \times \log (1 + z) + \varepsilon _Z \times (\log (Z) - \log (Z_ \odot)) \nonumber\\ &&\hphantom{\log (\alpha _{\rm CO 1}) =} + \varepsilon _{\rm sSFR} \times \log ({{\rm sSFR/sSFR}}({\rm ms},z,M_*)){\rm } \nonumber\\ && \hphantom{\log (\alpha _{\rm CO 1}) =} + \varepsilon _* \times (\log M_* - 10.5), \end{eqnarray} \tag{ 17 }$$

where ε_z, ε_Z, ε_sSFR, and ε_* are the logarithmic slopes of the variations of α_CO as a function of redshift, metallicity, sSFR relative to that at the main sequence, and stellar mass offset from logM_* = 10.5. We define α_{CO 10} as the zero point of the conversion factor, that is, its value at z = 0, sSFR = sSFR(ms, z, M_*), logM_* = 10.5 and Z = Z₈. The metallicity dependence is assumed to follow Equations (6) and (7) with ε_Z ∼ −1.2. Next we can write

$$\begin{eqnarray} &&\log t_{\rm depl} = \log \left({\frac{{\alpha _{\rm CO} \times L^{\prime}_{\rm CO} }}{{\rm SFR}}} \right) = \log t_{\rm depl} |_{\alpha _{\rm CO{\rm }1} = \alpha _{\rm MW} } \nonumber\\ &&\hphantom{\log t_{\rm depl} =} + {\rm }\log (\alpha _{{\rm CO}{\rm }1{\rm }0} /\alpha _{\rm MW}) \nonumber\\ &&\hphantom{\log t_{\rm depl} =} + \varepsilon _z \times \log (1 + z) + \varepsilon _Z \times (\log (Z) - \log (Z_ \odot)) \nonumber\\ &&\hphantom{\log t_{\rm depl} =} + \varepsilon _{\rm sSFR} \times \log ({\rm sSFR/sSFR}({\rm ms},z,M_*)) \nonumber\\ &&\hphantom{\log t_{\rm depl} =} + \varepsilon _* \times (\log M_* - 10.5), \end{eqnarray} \tag{ 18 }$$

where the left side of the equation is the "true" depletion timescale as given in Equation (16), while the first term on the right side is the "observed" depletion timescale under the assumption that α_CO = α_MW = constant, as in Table 3 and Equation (13):

$$\begin{eqnarray} &&\hspace{-25pt} \log t_{\rm depl} |_{\alpha _{{\rm CO}1} = \alpha _{\rm MW} } = a_{f1} + \xi _{f1} \times \log (1 + z) + \xi _{g1} \nonumber\\ &&\hspace{-25pt}\hphantom{\log t_{\rm depl} |_{\alpha _{{\rm CO}1} = \alpha _{\rm MW} } = } \times \log ({\rm sSFR/sSFR}({\rm ms},z,M_*))\nonumber\\ &&\hspace{-25pt}\hphantom{\log t_{\rm depl} |_{\alpha _{{\rm CO}1} = \alpha _{\rm MW} } = } + \xi _{h1} \times (\log (M_*) - 10.5). \end{eqnarray} \tag{ 19 }$$

Proceeding as before for the depletion timescale and gas to stellar mass ratio, the dependence of the observed dust temperature on redshift, sSFR offset, and stellar mass can also be separated into the product of three power-law functions (as for t_depl and M_{mol gas}/M_*) to yield

$$\begin{eqnarray} && \log T_{\rm dust} = a_{f3} + \xi _{f3} \, \times \, \log (1 + z) + \xi _{g3} \nonumber\\ &&\hphantom{\log T_{\rm dust} =}\times \log ({\rm sSFR/sSFR}({\rm ms},z,M_*))\nonumber\\ &&\hphantom{\log T_{\rm dust} =} + \xi _{h3} \times (\log (M_*) - 10.5). \end{eqnarray} \tag{ 20 }$$

Figures 10 and 11 show the scaling relations of the dust temperature in the same 512 galaxy stacks as in Figures 5 and 6 and Figures 9 and 10, respectively, and Table 5 summarizes the best-fit values for the power-law fitting function in Equation (20).

Zoom In Zoom Out Reset image size

Equations (16) and (18)–(20) can now be combined to yield

$$\begin{eqnarray} && \log (T_{\rm dust}) = c_3 - \frac{1}{{4 + \beta }} \nonumber\\ && \begin{array}{l} \times \left[ \begin{array}{@{} l @{}} (\xi _{\rm dg} + \varepsilon _Z) \times (\log Z - \log Z_ \odot) + (\xi _{f1} + \varepsilon _z) \times \log (1 + z) \nonumber\\ \quad+ (\xi _{g1} + \varepsilon _{\rm sSFR}) \times \log ({\rm sSFR/sSFR}({\rm ms},z,M_*)) \nonumber\\ \quad+ (\xi _{h1} + \varepsilon _*) \times (\log M_* - 10.5) \nonumber\\ \end{array} \right] \nonumber\\ {\rm } = c_4 + \xi _{f3} \times \log (1 + z) + \xi _{g3} \times \log ({\rm sSFR/sSFR}({\rm ms},z,M_*)) \nonumber\\ \quad {\rm } + \xi _{h3} \times (\log M_* - 10.5), \nonumber\\ \end{array}\hspace{-8pt} \nonumber\\ \end{eqnarray} \tag{ 21 }$$

which after sorting finally results in

$$\begin{eqnarray} 0 &=& {\rm } ({0.35(\pm 0.14){\rm } - {\rm }\log (\alpha _{{\rm CO}{\rm }1{\rm }0} /\alpha _{\rm MW} }) \nonumber\\ && + {\rm } ({ - \xi _{dg} - \varepsilon _Z }) \times (\log Z - 8.67) \nonumber\\ && + ({ - \xi _{f1} - (4 + \beta) \times \xi _{f3} - \varepsilon _z }) \times \log (1 + z) \nonumber\\ && +{\rm } ({ - \xi _{g1} - (4 + \beta) \times \xi _{g3} - \varepsilon _{\rm sSFR} }) \nonumber\\ && \times \log ({\rm sSFR/sSFR} ({\rm ms},z,M_*)) \nonumber\\ && + ({ - \xi _{h1} - (4 + \beta) \times \xi _{h3} - \varepsilon _* }) \times (\log M_* - 10.5). \qquad \end{eqnarray} \tag{ 22 }$$

This shows that the logarithmic slopes of the dependence of α_CO on Z, z, sSFR/sSFR(ms, z, M_*), and M_* can be uniquely constrained from the scaling relations of t_depl and T_dust. If Equation (22) is fulfilled everywhere in the sampled parameter space and the different variables on the right are independent, each of the coefficients in front of the variables on the right must be zero.

Constraints for main-sequence population. With the fit results in Tables 3 and 5, and the assumption that the Herschel dust observations provide ground truth, we find for the near-main-sequence population

$$\begin{eqnarray} \alpha _{{\rm CO}{\rm }1{\rm }0} &=& 9.8\,{\rm }(\pm 3.2) \nonumber\\ \varepsilon _Z &=& - 0.9{\rm }(\pm 0.3) \nonumber\\ \varepsilon _z &=& - 0.4{\rm }(\pm 0.13) \nonumber\\ \varepsilon _{\rm sSFR} &=& 0{\rm }(\pm 0.03) \nonumber\\ \varepsilon _* &=& 0.07{\rm }(\pm 0.02). \end{eqnarray} \tag{ 23 }$$

The inferred metallicity dependence of the CO-conversion factor is broadly consistent with Equations (6) and (7) (ε_Z ∼ −1.2 ± 0.3). The derived redshift dependence suggests that the conversion factor at z ∼ 2.2 is ∼0.7 times that at z ∼ 0, which would then imply a steeper gradient (ξ_f1 ∼ −0.5), which is in agreement with the dust depletion time evolution in Figure 5. The fact that the inferred zero point of the conversion factor is twice α_MW is also consistent with the z = 0 shift between dust- and CO-depletion times in Figure 5. This last conclusion is strongly dependent on the dust modeling assumptions in Draine & Li (2007). For instance, for a simple modified blackbody, modeling the dust masses would come down by a factor of ∼0.5–0.7 dex (Magnelli et al. 2012a; S. Berta et al. in preparation), in that case resulting in a zero point of α_{CO 10} = 2–3.

In our opinion the most important result is that the CO conversion factor near the main sequence depends little on sSFR. This is consistent with the earlier studies of Magdis et al. (2012b) and Magnelli et al. (2012a), but with much improved statistical confidence. Across Δlog(sSFR/sSFR(ms, z, M_*)) = ±0.6 the dust temperature measurements set an upper limit to a change in the CO conversion factor of 6% for the mean value, and 25% if the 2σ uncertainties are included. This is a strong constraint that applies as long as the dust measurements indeed provide a ground-truth estimate.

Constraints for above main-sequence outliers. The dust temperature gradient in Figure 12 appears to steepen (ξ_g3 increases from 0.086 to ∼0.135) for the outliers above the main sequence (up-bend of the distribution of points in the left panel for Δlog(sSFR/sSFR(ms, z, M_*)) ⩾ +0.9), implying greater radiation field densities than on the main sequence (see Magdis et al. 2012a). At the same time, the CO- and dust-depletion times in this regime are comparable (no significant change in ξ_g1 (CO) ∼ −0.46), albeit with increased scatter (left panel of Figure 6). Equation (22) then implies a drop in the conversion factor (ε_sSFR = −ξ_g1 − (4 + β) × ξ_g3). At Δlog(sSFR/sSFR(ms, z, M_*)) ∼ +1–1.3 the residuals from the slope 0.086 power-law fit (red-dotted line in Figure 12) are about +0.05 dex, which implies a decrease of the conversion factor to α_CO ∼ 2 (ε_sSFR ∼ −0.3). Within the 0.3 dex uncertainties this is consistent with α_CO ∼ 0.8–1.5, the value empirically estimated by Downes & Solomon (1998) and Scoville et al. (1997) for z ∼ 0 ULIRGs (see also Bolatto et al. 2013). If such a change of α_CO is applied, depletion timescales for the extreme outliers decrease by 0.3–0.7 dex, suggesting the presence of a second, starburst mode, as discussed by Daddi et al. (2010b), Genzel et al. (2010), Magdis et al. (2012a, 2012b), and Sargent et al. (2014). The main quantitative difference to Magdis et al. (2012a, 2012b) and Sargent et al. (2014) is that our data perhaps suggest that the deviations in α_CO become significant only about one order of magnitude, and not, as in these papers, already a factor of four above the main-sequence line. The statistics of our data in Figures 6 and 12, however, are not sufficient to distinguish a continuous from an abrupt change in α_CO, and one needs to keep in mind the large uncertainties of stellar masses and SFRs for this obscured, bursty population (Section 2.4).

Zoom In Zoom Out Reset image size

Finally the last constraint on the stellar mass dependence implies a change of α_CO of less than 7% from 10¹⁰ to 10¹¹ M_☉ in stellar mass.

4.2.1. Final Global Fits

4.2.2. Dependence of the Results on the Prescription of the Main Sequence

4.2.3. Fitting in sSFR Space

4.2.4. Impact of Metallicity Descriptions

4.2.5. Systematic Uncertainties and Missing Parameters

From a theoretical perspective, the slow change of the molecular gas depletion time with cosmic epoch (t_depl ∝ (1+z)^{−0.34±0.15}) is somewhat surprising (e.g., Kauffmann et al. 1993; Cole et al. 1994; Elmegreen 1997; Silk 1997). Considering the definition of the depletion time in the context of the Kolmogorov–Smirnov (K-S) relation (Equation (2)), t_depl might naturally be thought to be proportional to the galaxy's dynamical time, with the proportionality being the inverse of the galaxy's star-formation efficiency η (Kauffmann et al. 1993; Cole et al. 1994; Silk 1997; Elmegreen 1997; Kennicutt 1998; Genzel et al. 2010). In the Mo et al. (1998) framework of disk formation in a dark matter dominated universe, the disk's dynamical time t_dyn (R_d) (expressed in terms of its scale length R_d and maximum rotation velocity v_d) is tied to the properties of the dark matter halo, such that

$$\begin{equation} t_{\rm depl} = \frac{{t_{\rm dyn} (R_d)}}{\eta } = \frac{1}{\eta } \times \frac{{R_d }}{{{\rm v}_d }} = \frac{{\lambda R_h }}{{\eta \times C_h v_h }} = \frac{\lambda }{{\eta C_h }} \times 0.1H(z)^{ - 1}. \end{equation} \tag{ 24 }$$

Here R_h and v_h are the halo size and circular velocity, and C_h is the ratio of the disk's rotation velocity to the halo circular velocity, which depends on the halo concentration (Bullock et al. 2001b). The parameter λ is the angular momentum parameter of the baryons (〈λ〉_{dark matter} ∼0.04; Bullock et al. 2001a) and H(z) = H₀ × (Ω_Λ + Ω₀ × (1+z)³)^1/2 is the Hubble parameter. In a matter dominated universe (applicable at high-z) the depletion time should then be proportional to (1+z)^−3/2 (see Davé et al. 2011), which is inconsistent with our results. A more careful evaluation requires two corrections. First the Hubble parameter in a ΛCDM universe changes more slowly at late times. If one approximates H(z)∝(1+z)^β_, the average β for the redshift range z = 0–2.5 is −0.98. Second, the concentration parameter of dark matter halos was smaller at high-z than at z = 0 (Bullock et al. 2001b), such that C_h ∼ 1.025 at z ∼ 2.5 and 1.24 at z ∼ 0 (see Somerville et al. 2008). Taken together, these two effects change the effective redshift dependence in Equation (16) to (1+z)^β with β = −0.83 (again between z = 0 and 2.5). In a recent evaluation of star-forming disk sizes in CANDELS/3D HST between z = 0 and 3, van der Wel et al. (2014) empirically find β = −0.75 (±0.05), which is smaller than β = −0.83, but comes close to it. Additional baryonic processes connected to the processing and dissipation of the angular momentum from the scale of the cosmic web to the inner, star-forming disk and feedback processes inside the disk can both increase and decrease λ of the baryonic component (Dutton et al. 2011; Danovich et al. 2014). Empirically the average observed λ of the star-forming gas at z ∼0.8–2.5 is about 0.035, which is (fortuitously) close to the dark matter λ (A. Burkert et al. in preparation).

The difference between the values of van der Wel's −0.75 (±0.05) and our average estimate of −0.34 (±0.15) for the slope in the redshift dependence of the depletion time is not highly significant. If it were, the shallow slope may suggest that the depletion timescale is not set by the global galactic dynamical time but by local processes. There are good reasons for this view. Cloud collapse and star formation are local processes that proceed on the free-fall time τ_ff, which depends on the inverse square root of the local gas density (e.g., Krumholz & McKee 2005; Krumholz et al. 2012). Empirically Leroy et al. (2013) find from spatially resolved observations of the K-S relation in 30 z ∼ 0 disk galaxies that the depletion time is near constant and independent of the local or global orbital times. However, the high-z disk galaxies in our sample appear to be globally marginally stable systems with a Toomre parameter Q ⩽ Q_crit ∼ 1. In this case it is easy to show that the depletion timescale in the K-S relation on large scales is locked to the average galactic orbital time, even if in principal the local volumetric SFR density is tied to the local free-fall timescale (e.g., Genzel et al. 2010; Krumholz et al. 2012). The dependence of the depletion timescale on a sSFR (t_depl∝(sSFR/sSFR(ms, z, M_*))^−1/2) may be considered to be another argument in favor of a local origin of the depletion timescale.

The scaling relations in Tables 3 and 4 can be used to predict the form of the galaxy integrated molecular gas–star formation relations at different redshifts and under varying selection functions. Because of the strong dependence of the depletion time on sSFR (t_depl ∝ (sSFR/sSFR(ms, z, M_*)^−0.5), the relation between M_{mol gas} and SFR (the molecular K-S relation; Kennicutt 1998) becomes super-linear for a sample of SFGs with a spread in sSFR, although intrinsically the relation at constant sSFR is linear.

We explored the impact of this dependence more quantitatively by creating mock galaxy samples at different redshifts and with varying spread in sSFR/sSFR(ms, z, M_*), and establishing the resulting M_{mol gas}–SFR from the relations in Equation (1) and Table 3. For this purpose we repeatedly drew samples of 10–100 SFGs at a given z (the upper value is an upper limit to the samples available now and in the near future, at least for z > 0.5). At each redshift we varied the redshifts by ±10%–30% around the mean redshift and varied stellar masses from log(M_*/M₈) = 10.3 to 11.3. We also varied SFR and M_{mol gas} estimates by ±0.2 dex to reflect measurement uncertainties, and Δ(log(sSFR/sSFR(ms, z, M_*)) to sample the vertical direction in the M_*–SFR plane. The final mock sample gives approximately uniform coverage in all these parameters.

At all z, the resulting integrated K-S slope N = dlog(SFR)/d(M_{mol gas}) in these mock data sets approaches ∼1 for large N and a small range in sSFR. However, even for a pure main-sequence sample Δ(log(sSFR/sSFR(ms, z, M_*)) = ±0.3, the slope becomes super linear, N ∼ 1.1–1.2, and with a substantial scatter of ±0.07 dex resulting from the variation in SFG parameters. The scatter naturally increases with the decreasing number n_G of galaxies in the sample, from ±0.06 for n_G = 50–100, to ±0.13 for n_G = 20 and ±0.2 for n_G = 10. The slope increases with the range of sSFR, from N = 1.08 for Δ(log(sSFR/sSFR(ms, z, M_*)) = ±0.3, to N = 1.33 for Δ(log(sSFR/sSFR(ms, z, M_*)) = ±0.6 and N = 1.6 for Δ(log(sSFR/sSFR(ms, z, M_*)) = ±1. Finally, the slope for fixed Δ(log(sSFR/sSFR(ms, z, M_*)) increases by about 0.2 from z = 0 to z = 2.

These findings are in excellent agreement with the integrated galaxy K-S slopes in the literature, both at low-z (Kennicutt 1998; Bigiel et al. 2008; Leroy et al. 2013; Saintonge et al. 2012; Sargent et al. 2014), as well as z ⩾ 1 (Daddi et al. 2010b; Genzel et al. 2010; Tacconi et al. 2013; Sargent et al. 2014). In particular, the dependence of the slope on Δ(log(sSFR/sSFR(ms, z, M_*)) is in excellent agreement with the findings of Kennicutt (1998), Daddi et al. (2010b), Genzel et al. (2010), Magdis et al. (2012b), Sargent et al. (2014), and Tacconi et al. (2013). Note that for galaxy samples with extreme selection functions, such as a combination of a near-main-sequence sample with a sample of starburst outliers (such as those studied by Daddi et al. 2010b and Genzel et al. 2010), the continuous change in t_depl with sSFR would then appear like the superposition of two separate gas–SFR relations, as proposed by Sargent et al. (2014).

Due to the lack of size information, we cannot study the more common K-S surface density relation Σ_{mol gas}–Σ_SFR. Whether or not extra factors come into play in the surface density relation as compared to the integrated quantities depends on the underlying physical reason of the sSFR dependence on t_depl, and the probably connected inverse scaling of sizes with sSFR (Wuyts et al. 2011b; Elbaz et al. 2011). Krumholz et al. (2012) have pointed out that if the K-S relation is intrinsically a volumetric relation between the gas and star-formation volume densities, ρ_SFR∝ρ_gas/τ_ff(ρ_gas), the scale height h_gas of the gas comes in as an additional factor, such that the galaxy-averaged surface density relation becomes

$$\begin{eqnarray} \Sigma _{\rm SFR} &=& {\rm } \langle \rho _{\rm SFR} \times h_{\rm gas} \rangle {\rm } \propto {\rm }h_{\rm gas}^{ - 1/2} {\rm } \times {\rm } \langle \rho _{\rm gas} \times h_{\rm gas} \rangle ^{3/2} \nonumber\\ & \propto & {\rm } h_{\rm gas}^{ - 1/2} \times \Sigma _{\rm gas}^{3/2}, \end{eqnarray} \tag{ 25 }$$

since the local free-fall timescale $\tau _{ff} (\rho _{\rm gas}) \propto \rho _{\rm gas}^{ - 1/2}$. If the scaling t_depl ∝ sSFR^−1/2 is physically connected with average local gas densities increasing with sSFR, as indicated by the findings of Wuyts et al. (2011b) and Elbaz et al. (2011), Equation (25) would be in qualitative agreement with our scaling relation, with the addition of the scale height factor.

To summarize these results qualitatively, an intrinsically linear K-S relation at a constant vertical position in the sSFR–M_* plane transforms into a super-linear relation once there is a range of sSFRs, either in sample of galaxy integrated measurements, or in a spatially resolved sample within a galaxy for which the initial K-S relation is applicable. This is because at greater sSFR (or SFR or Σ_SFR) the gas fraction increases (pushing a data point to the right in the K-S plane) and the depletion time decreases (pushing the data point up in the K-S plane), resulting in a combined shift to the upper right along the diagonal. When data points with a range of sSFR (or SFR or Σ_SFR) are combined, this appears as a super-linear relation with scatter.

We have shown in this paper that the derivation of molecular gas masses in SFGs either from a low-J CO rotational line emission (using a CO conversion function that is dependent on metallicity and contains a simple excitation correction) or from full far-IR SEDs (with a conversion to gas mass that uses the Draine & Li 2007 dust model, a metallicity dependent gas to dust ratio from Leroy et al. 2011, and the mass–metallicity relation), yield consistent results across a wide range of z, sSFR/sSFR(ms, z, M_*), and M_*. We argued that the combined scaling relations in Tables 3–5 capture the most important variations of the molecular gas depletion timescale, molecular to stellar mass ratios, and average dust temperature, within the ±0.24 dex, ±0.24 dex, and ±0.033 dex scatter of the three relations, and to an absolute level of ±25%.

Extrapolating to future work in the measurement of cold gas masses in galaxies, it is likely that the focus will be on the expansion of the statistics and parameter range, including studies on the dependence on parameters that we have not been able to explore (such as galaxy morphology and environment). On the other hand, spatially resolved measurements are likely to play an increasingly important role, especially at higher redshift, in order to explore the dependence of the depletion timescale and gas fraction on internal galaxy structure, such as bulge-to-disk mass ratio, molecular volume and surface densities, clumpiness, and internal galaxy kinematics, including galactic turbulence. The latter measurements will require high-resolution kinematic data, which come for free with molecular line observations of CO, HCN, etc.

The former goal primarily requires efficient measurements of galaxy integrated gas masses. To this end, it is well known that the detection of a given gas mass from broadband detection of its submillimeter dust emission is substantially faster than from CO 2–1 or 3–2 data. In the case of ALMA, the detection of a molecular gas mass of 10¹⁰ M_☉ from band 7 (350 GHz) broadband data requires only a few minutes at z ∼ 1–2 (for a 5σ detection), whereas a CO-based measurement (again at 5σ) requires more than 1 hr at z ∼ 0.6–1 and several hours at z ∼ 2 (Scoville 2013; Scoville et al. 2014).²⁴

For these reasons, Scoville (2013) and Scoville et al. (2014) proposed that a single frequency, broadband measurement in the Rayleigh–Jeans tail of the dust SED (for instance at 345 GHz) is sufficient to establish dust and gas masses. Scoville and his colleagues argue that the variation of dust temperature on galactic scales is sufficiently small, such that the assumption of a constant dust temperature, T₀ ∼ 25 K, is justified. Qualitatively this assumption is in very good agreement with the slow changes of average T_dust with the redshift and sSFR in the stacked Herschel data (Magnelli et al. 2014) we presented in Figures 11 and 12. Based on SCUBA observations of a subset of a sample of z ∼ 0 disks from Draine et al. (2007), Scoville et al. also argue that the dust-to-gas ratio does not vary significantly with metallicity (δ_dg ∼ 0.0067 ∼ constant). This assumption is obviously at tension with the Leroy et al. (2011) scaling relation in Equation (10), which predicts a fairly strong metallicity dependence of δ_dg. This tension needs to be resolved in future studies.

The Rayleigh–Jeans tail method thus is based on a single broad-band flux measurement, using calibrations of the 350 GHz dust emissivity to dust/gas mass from Planck observations in the Milky Way (Planck Collaboration et al. 2011a, 2011b) and of the nearby SINGS galaxies (Draine et al. 2007), yielding

$$\begin{eqnarray} \left({\frac{{M_{\rm mo lg as} }}{{1 \times 10^{10} M_ \odot }}} \right) &=& \left({\frac{{S_{350} D_L^2 }}{{{\rm mJy} \times Gpc^2 }}} \right) \times (1 + z)^{ - (3 + \beta)} \nonumber\\ && \times \left({\frac{{\nu _{\rm obs} }}{{350 \; {\rm GHz} }}} \right)^{ - (2 + \beta)}, \end{eqnarray} \tag{ 26 }$$

for T_dust = 25 K, the above calibration of κ_dust (350 GHz) (Scoville et al. 2014) and β = 1.8.

It is instructive to compare this method to more detailed measurements that would include the temperature variations in Table 5, the full Planck formula of the dust modified gray-body emission, and the Leroy et al. (2011) recipe for δ_dg(Z) (Equation (10)).

For this purpose we built a simulation to verify the performance of Equation (26) on simulated data points with the above assumptions. We use the scaling relations in Tables 4 and 5 to define a grid of modified blackbody SEDs (MBBs) in the M_*–sSFR–z space, using the binning by Magnelli et al. (2014). Each point of the grid is characterized by M_*, sSFR, z, T_dust, and M_{mol gas}. Based on the scatter of our scaling relations, we adopted σ(logM_{mol gas}) = 0.23 and σ(logT_dust) = 0.033 to reflect variations in the average properties of a bin. We convolved the MBBs with a box filter centered on ALMA's band 7 (350 GHz), and computed the resulting flux density (in mJy). We then added Gaussian noise as computed for a given integration time from the ALMA sensitivity calculator assuming 34 12 m antennas, which were given as the default in the cycle two time estimator. The resulting observed flux density is then converted to M_{mol gas} using Equation (26). Figure 13 (left panel) presents the results of the simulation and compares input to output gas masses. This figure shows that the Rayleigh–Jeans approximation in Equation (26), assuming no metallicity dependent dust-to-gas ratio, leads to a significant underestimation of all inferred gas masses (on average −0.3 dex), and to artificial systematic trends throughout the probed parameter space (±0.4 dex), especially at high-z.

Zoom In Zoom Out Reset image size

A first-order improvement comes from introducing a constant Planck correction for all galaxies in any sample. For T_dust = T₀ = 25 K = constant one multiplies Equation (26) with

$$\begin{equation} \Gamma _{\rm Pl} = \frac{{\exp (h\nu _{\rm obs} (1 + z)/kT_0) - 1}}{{h\nu _{\rm obs} (1 + z)/kT_0 }},\end{equation} \tag{ 27 }$$

where ν_obs = 350 GHz (Scoville et al. 2014). The central panel in Figure 12 shows that with this global Planck correction, the overall underestimate of gas masses is corrected and in fact over-compensated, mainly because the adopted value of T₀ is below the actual mean dust temperature near the main sequence. The large (±0.35 dex) systematic trends remain because of the intrinsic variation in T_dust in Table 5 and the dust to gas ratio variations as a function of metallicity/mass. Correcting for the latter effect with the mass–metallicity relation in Equation (10) improves the estimates further (right panel of Figure 13), but the temperature variations still cause significant systematic trends that should be accounted for in future studies attempting to measure fairly accurate relative trends in gas mass or depletion time. Obviously if the scaling relations in Table 5 are applicable to the galaxy sample, a full Planck correction with a variable T_dust can be executed with Equation (27), which should then give the correct output gas masses, making the single-band technique an efficient method for determining gas masses for large samples.

If one does not wish to rely on the scaling relations in Table 5, a dust temperature must be established for every individual galaxy. This can be done using measurements in two ALMA bands, exploiting the frequency dependence of the Planck correction in Equation (27). We simulated the performance of this two-band technique, starting with 10σ photometry in either the band 6 (240 GHz) and band 7 (350 GHz) combination (Figure 14), or the band 7 and band 9 (670 GHz) combination (Figure 15). The band 6/7 combination has the advantage of less demanding observing conditions, but it delivers less accurate output dust temperatures (left panels in Figure 14) and gas masses (middle panels in Figure 14). For 10σ photometry, the resulting fractional precision of gas masses is not better than 0.7. The band 7/9 combination performs better in this respect (by about a factor of two), but band 9 observations require much better, and thus somewhat rarer atmospheric conditions. For the quoted fractional precisions, the total observing times for an input gas mass of ∼10¹⁰ M_☉ are >100 minutes at z ∼ 0.6–2.2. At least for z > 1.5 this is still somewhat more favorable than CO detections (which are usually done in the 3 mm atmospheric window) but the requirement of band 9 measurements for many galaxies will likely be quite restrictive in terms of sample size.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size