THE EVOLVING INTERSTELLAR MEDIUM OF STAR-FORMING GALAXIES SINCE z = 2 AS PROBED BY THEIR INFRARED SPECTRAL ENERGY DISTRIBUTIONS

We use deep 100 and 160 μm PACS and 250, 350, and 500 μm SPIRE observations from the GOODS-Herschel program. Details about the observations are given in Elbaz et al. (2011). Herschel fluxes are derived from point-spread function (PSF) fitting using galfit (Peng et al. 2002). A very extensive set of priors, including all galaxies detected in the ultra-deep Spitzer Multiband Imaging Photometer (MIPS) 24 μm imaging, is used for source extraction and photometry at 100, 160, and 250 μm, which effectively allow us to obtain robust flux estimates for relatively isolated sources, even beyond formal confusion limits at 250 μm. For 350 and 500 μm, this approach does not allow accurate measurements due to the increasingly large PSFs. Hence, we use a reduced set of priors based primarily on Very Large Array (VLA) radio detections, resulting in flux uncertainties consistent with the confusion noise at the IR wavelengths. Our measurements are in good agreement with the alternative catalogs used in Elbaz et al. (2011). The advantage of using galfit for PSF fitting is in its detailed treatment of the covariance matrix to estimate error bars in the flux measurements, which is crucial to eventually derive a reliable estimate of flux errors for the case of blended/neighboring sources. The effective flux errors in each band can vary substantially with position, depending on the local density of prior sources over areas comparable to the PSF. A detailed description of the flux measurements and Monte Carlo (MC) derivations of the uncertainties will be presented elsewhere (E. Daddi et al. 2012, in preparation). We correct the PACS photometry for a small flux bias introduced (due to source filtering for background subtraction) during data reduction (see H-GOODS public data release¹⁴; P. Popesso et al. 2012, in preparation).

Daddi et al. (2010a) presented Plateau de Bure Interferometer (PdBI) CO[2–1] emission line detections of six star-forming galaxies at z ∼ 1.5, originally selected by using the "star-forming BzK" color criterion (Daddi et al. 2004b). All six sources (BzK-21000, BzK-17999, BzK-4171, BzK-16000, BzK-12591, and BzK-22536) have spectroscopic redshifts (now confirmed by multiple CO detections; for optical redshifts see D. Stern et al. 2012, in preparation; Cowie et al. 2004) and robust PACS and/or SPIRE detections. Four of these sources also benefit from 1.3 millimeter continuum data derived as a by-product from the CO[5–4] emission line observations (H. Dannerbauer et al. 2012, in preparation). For BzK-21000, continuum detections and upper limits are also obtained at 1.1, 2.2, and 3.3 mm (Daddi et al. 2010a; Dannerbauer et al. 2009). In addition to IRAC and MIPS 24 μm data, the sources are seen in the 16 μm InfraRed Spectrograph peak-up image (Teplitz et al. 2011). Although we will revisit (and confirm) the SFRs of these sources based on the Herschel data, the existing UV, mid-IR, and radio SFR estimates, along with stellar mass measurements derived by fitting the Bruzual & Charlot (2003) model SEDs to their rest-frame UV to near-IR photometry (Daddi et al. 2010a), consistently place them in the SFR–M_* MS at z ∼ 1.5 (Daddi et al. 2010a). Finally, the UV rest-frame morphologies, the double-peaked CO profiles, the large spatial extent of their CO reservoirs, and the low gas excitation of the sources all provide strong evidence that they are large, clumpy, rotating disks (Daddi et al. 2010a).

We also consider three z ∼ 0.5 normal disks with Sérsic index n < 1.5 and spectroscopic redshifts, for which we have CO[2–1] emission line detections (Daddi et al. 2010b). Similar to the z ∼ 1.5 sample, the sources are detected in the PACS and/or SPIRE bands and are part of the SFR–M_* MS based on their 24 μm derived IR luminosities and SFRs, which are known to be robust in this redshift range (Elbaz et al. 2010, 2011). The Herschel and millimeter photometries of the MS galaxies considered here are summarized in Table 1.

As a comparison sample to our MS galaxies, we also include in our analysis a small sample of SMGs: GN 20, SMM J2135-0102, and HERMES J105751.1+573027. The selection of the targets was driven by the need for available multi-wavelength photometry, including Herschel and millimeter continuum observations, as well as CO[1–0] emission line detections.

GN 20 was already studied in Magdis et al. (2011). It is one of the best-studied SMGs to date, the most luminous, and also one of the most distant (z = 4.05; Pope et al. 2006; Daddi et al. 2009) in the GOODS-N field. Herschel photometry and the far-IR/millimeter properties of the source have already been presented in detail in Magdis et al. (2011). In brief, the source is detected in all Herschel bands (apart from 100 μm) and in the AzTEC 1.1 mm map (Perera et al. 2008), while continuum emission is also measured at 2.2, 3.3, and 6.6 mm (Carilli et al. 2011), as well as at 1.4 GHz with the VLA (Morrison et al. 2010). Furthermore, Carilli et al. (2010) reported the detection of the CO[1–0] and CO[2–1] lines with the VLA and CO[6–5] and CO[5–4] lines with the PdBI and the Combined Array for Research in Millimeter Astronomy (CARMA), respectively. A compilation of the photometric data is given in Table 1 of Magdis et al. (2011).

SMM J2135-0102 (SMM-J2135 hereafter) is a highly magnified SMG at z = 2.325, with an amplification factor of μ = 32.4 ± 4.5, serendipitously discovered by Swinbank et al. (2010) behind the cluster MACS J2135-01. The source has exquisite photometric coverage from rest-frame UV to radio including SPIRE broadband data (Swinbank et al. 2010, Table 1; Ivison et al. 2010, Table 1). Swinbank et al. (2010) also report the detection of CO[1–0] and CO[3–2] emission lines using the Zpectrometer on the Green Bank Telescope and PdBI observations, while Ivison et al. (2010), using the SPIRE Fourier Transform Spectrometer (FTS), presented the detection of the [C ii]158 μm cooling line.

Finally, HERMES J105751.1+573027 (HSLW-01 hereafter) is a Herschel/SPIRE-selected galaxy at z = 2.957, multiply lensed (magnification factor μ = 10.9 ± 0.7) by a foreground group of galaxies. The source was discovered in Science Demonstration Phase Herschel/SPIRE observations of the Lockman-SWIRE field as part of the Herschel Multi-tiered Extragalactic Survey (HerMES; S. Oliver et al. 2012). The optical to millimeter photometry of the source is presented in Table 1 of Conley et al. (2011), while Riechers et al. (2011a) report the detection of CO[5–4], CO[3–2], and CO[1–0] emission using the PdBI and CARMA and the Green Bank Telescope.

We note that due to the high magnification factor for SMM-J2135 and HSLW-01, these sources might not be representative of the bulk population of SMGs, traditionally selected as S₈₅₀ > 5 mJy. Also, a large number of non-lensed SMGs with CO[3–2] observations are available in the literature. However, we choose not to consider them in this study, as the uncertain gas excitation results in dubious CO[1–0] estimates (e.g., Ivison et al. 2010; Riechers et al. 2011b), which are essential for investigating the gas properties of the sources.

In addition to our samples of individually detected MS galaxies at z ∼ 1.5 and z ∼ 0.5, we also perform a stacking analysis to derive the average SED of MS galaxies at z ∼ 1 and z ∼ 2. Specifically, we use the z ∼ 1 GOODS samples of Salmi et al. (2012; see also Daddi et al. 2007; M. Pannella et al. 2012, in preparation). In order to remove bulge-dominated galaxies with low sSFR, we discard objects with a Sérsic index n > 1.5 (Salmi et al. 2012), as measured from the GOODS ACS images (Giavalisco et al. 2004), and consider only 24 μm detected galaxies from the GOODS ultra-deep Spitzer imaging. To remove SBs, we proceed as follows: (1) from M_* and redshift, we compute sSFR_MS, i.e., the fiducial mean sSFR expected for MS galaxies with a given mass and redshift, derived as described at the start of Section 2; (2) from Herschel data we derive L_IR and subsequently SFR estimates and exclude galaxies with measured sSFR/sSFR_MS > 3. Since not all SB systems are expected to be detected by Herschel, we also omit sources with sSFRs/sSFR_MS > 3 using the UV-based, corrected-for-extinction, SFR estimates by Daddi et al. (2007) and the Spitzer estimates of Salmi et al. (2012). The z ∼ 2 sources are drawn from the BzK-selected sample of Daddi et al. (2007; now extended to deeper K-band depths by M. Pannella et al. 2012, in preparation). Since this color-selection technique naturally excludes quiescent galaxies, we only exclude SBs identified using the same method as for the z ∼ 1 sample. The total number of z ∼ 1 and z ∼ 2 galaxies considered in the stacking analysis is 1569 and 3618, respectively. The average SEDs of these two samples are then measured by stacking, using various standard techniques depending on the wavelength:

• 1.  
  16 μm. We use the Teplitz et al. (2011) maps. The stacking is performed using the IAS stacking library (Béthermin et al. 2010) and PSF-fitting photometry.
• 2.  
  24 and 70 μm. We use the 24 μm and 70 μm (Frayer et al. 2006) maps of GOODS, the IAS stacking library, and aperture photometry with parameters similar to that in Béthermin et al. (2010).
• 3.  
  100 and 160 μm. We use the GOODS-Herschel PACS data (Elbaz et al. 2011), the IAS stacking library, and PSF fitting. We apply the appropriate flux correction for faint, non-masked sources to the PACS stacks (P. Popesso et al. 2012, in preparation).
• 4.  
  250, 350, and 500 μm. We use the GOODS-Herschel SPIRE (Elbaz et al. 2011) data and the mean of pixels centered on sources as in Béthermin et al. (2012b). The bias due to clustering is estimated to be about 20% at 500 μm and thus negligible compared with the statistical uncertainties. This bias is smaller at shorter wavelengths where the beam is narrower (Béthermin et al. 2012b).
• 5.  
  870 μm and 1100 μm. We used the Weiß et al. (2009) LABOCA map of eCDFS and the Perera et al. (2008) AzTEC map of GOODS-N. Contrary to SPIRE, (sub)millimeter data are noise limited; it is thus optimal to beam-smooth the map before stacking (J. D. Vieira et al. 2012, in preparation). We thus apply the Marsden et al. (2009) method to perform our stacking, which takes the mean of pixels centered on sources in the beam-smoothed map.

At all wavelengths, we used a bootstrap technique to estimate the uncertainties (Jauzac et al. 2011), and the mean fluxes measured in the two fields are combined quadratically to produce a mean SED at z ∼ 1 and z ∼ 2. For the z ∼ 2 normal galaxies, we also construct average SEDs in three stellar mass bins in the range of 9.7 < log (M_*/M_☉) < 11.2 with a bin size of Δlog (M_*/M_☉) = 0.5, as well as in four bins of sSFR/sSFR_MS over the range −0.4 < log (sSFR/sSFR_MS) < 0.4, with a bin size of Δlog (sSFR/sSFR_MS) = 0.2.

We employ the dust models of DL07, which constitute an update of those developed by Weingartner & Draine (2001) and Li & Draine (2001), and which were successfully applied to the integrated photometry of the Spitzer Nearby Galaxy Survey (SINGS) galaxies (Draine et al. 2007). These models describe the interstellar dust as a mixture of carbonaceous and amorphous silicate grains, whose size distributions are chosen to mimic the observed extinction law in the MW, the Large Magellanic Cloud (LMC), and the Small Magellanic Cloud (SMC) bar region. The properties of these grains are parameterized by the polycyclic aromatic hydrocarbon (PAH) index, q_PAH, defined as the fraction of the dust mass in the form of PAH grains. The majority of the dust is supposed to be located in the diffuse interstellar medium (ISM), heated by a radiation field with a constant intensity U_min. A smaller fraction γ of the dust is exposed to starlight with intensities ranging from U_min to U_max, representing the dust enclosed in photodissociation regions (PDRs). Although this PDR component contains only a small fraction of the total dust mass, in some galaxies it contributes a substantial fraction of the total power radiated by the dust. Then, according to DL07, the amount of dust dM_dust exposed to radiation intensities between U and U + dU can be expressed as a combination of a δ-function and a power law:

$$\begin{eqnarray} \frac{dM_{\rm dust}}{dU} &=& (1-\gamma) M_{\rm dust} \delta (U-U_{\rm min}) \nonumber\\ &&+ \gamma M_{\rm dust} {\alpha -1 \over U_{\rm min}^{1-\alpha } - U_{\rm max}^{1-\alpha }} U^{-\alpha }, \end{eqnarray} \tag{ 1 }$$

with (U_min ⩽ U_max, α ≠ 1), M_dust the total dust mass, γ the fraction of the dust mass that is associated with the power-law part of the starlight intensity distribution, and U_min, U_max, and α characterizing the distribution of starlight intensities in the high-intensity regions.

Following DL07, the spectrum of a galaxy can be described by a linear combination of one stellar component approximated by a blackbody with color temperature T_* = 5000 K and two dust components, one arising from dust in the diffuse ISM, heated by a minimum radiation field U_min ("diffuse ISM" component), and the other from dust heated by a power-law distribution of starlight, associated with the intense PDRs ("PDR" component). Then, the model emission spectrum of a galaxy at distance D is

$$\begin{eqnarray} f_{\nu }^{\rm model} &=& \Omega _* B_\nu (T_*) + {M_{\rm dust} \over 4 \pi D^2} \big[ (1-\gamma) p_\nu ^{(0)}(q_{\rm PAH},U_{\rm min}) \nonumber\\ && + \gamma p_\nu (q_{\rm PAH},U_{\rm min},U_{\rm max},\alpha) \big], \end{eqnarray} \tag{ 2 }$$

where Ω_* is the solid angle subtended by stellar photospheres, p⁽⁰⁾_ν(q_PAH, U_min) and p_ν(q_PAH, U_min, U_max, α) are the emitted power per unit frequency per unit dust mass for dust heated by a single starlight intensity U_min, and dust heated by a power-law distribution of starlight intensities dM/dU∝U^−α extending from U_min to U_max.

In principle, the dust models in their most general form are dictated by seven free parameters (Ω_*, q_PAH, U_min, U_max, α, γ, and M_d). However, Draine et al. (2007) showed that the overall fit is insensitive to the adopted dust model (MW, LMC, and SMC) and the precise values of α and U_max. In fact, they showed that fixed values of α = 2 and U_max = 10⁶ successfully described the SEDs of galaxies with a wide range of properties. They also favor the choice of MW dust models for which a set of models with q_PAH ranging from 0.4% to 4.6% is available. Furthermore, since small U_min values correspond to dust temperatures below ∼15 K that cannot be constrained by far-IR photometry alone, in the absence of rest-frame submillimeter data, they suggest using 0.7 ⩽U_min ⩽ 25. While this lower cutoff for U_min prevents the fit from converging to erroneously large amounts of cold dust heated by weak starlight (U_min < 0.7), the price to pay is a possible underestimate of the total dust mass if large amounts of cold dust are indeed present. However, Draine et al. (2007) concluded that omitting submillimeter data from the fit increases the scatter of the derived masses up to 50% but does not introduce a systematic bias in the derived total dust masses.

Under these assumptions, we fit the mid-IR to millimeter data points of each galaxy in our sample, searching for the best-fit model by χ² minimization and parameterizing the goodness of fit by the value of the reduced χ², χ²_ν ≡ χ²/df (where df is the number of degrees of freedom). The best-fit model yields a total dust mass (M_dust), U_min, γ, and q_PAH, while to derive L_IR¹⁵ estimates, we integrate the emerging SEDs from 8 to 1000 μm:

$$\begin{equation} L_{\rm IR} = \int _{8\,\mu {\rm m}}^{1000\,\mu {\rm m}} L_\nu (\lambda) \times \frac{c}{\lambda ^{2}} d\lambda. \end{equation} \tag{ 3 }$$

A by-product of the best-fit model is also the dust-weighted mean starlight intensity scale factor, 〈U〉, defined as

$$\begin{eqnarray} &&\langle U\rangle = \left\lbrace \begin{array}{@{} ll}\displaystyle \frac{L_{\rm dust}}{P_{0} M_{{\rm dust}}} \;{\rm or} \nonumber\\ \ms \displaystyle \left[(1-\gamma)U_{{\rm min}}+\frac{\gamma \ln (U_{{\rm max}}/U_{{\rm min}})}{U_{{\rm min}}^{-1}-U_{{\rm max}}^{-1}}\right], \;{\rm for}\;\alpha =2, \end{array} \right. \nonumber\\ \end{eqnarray} \tag{ 4 }$$

where P₀ is the power absorbed per unit dust mass in a radiation field with U = 1. Note that essentially 〈U〉 is proportional to L_IR/M_dust, and as we will discuss later for the definition of L_IR adopted here, i.e., L_{8–1000 μm}, our data suggest that P₀ ≈ 125. Uncertainties in L_IR and M_dust are quantified using MC simulations. To summarize, for each galaxy a Gaussian random number generator was used to create 1000 artificial flux sets from the original fluxes and measurement errors. These new data sets were then fitted in the same way, and the standard deviation in the new parameters was taken to represent the uncertainty in the parameters found from the real data set. Best-fit values along with their corresponding uncertainties are listed in Table 2 for all sources in our sample. To check for possible contamination of the submillimeter continuum broadband photometry from C+ (158 μm) emission (e.g., Smail et al. 2011), we repeated the fit excluding the affected bands (i.e., 350 μm and 500 μm for z ∼ 0.5 and z ∼ 1.5–2.0, respectively), without noticing any effect in the derived parameters. The best-fit models along with the observed SEDs of the z ∼ 1.5 and z ∼ 0.5 galaxies and SMGs are shown in Figures 1–3, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While Herschel data accurately probe the peak of the SED in the far-IR emission of distant galaxies, rest-frame submillimeter observations (λ_rest ⩾ 250 μm) are necessary to sample the Rayleigh–Jeans tail. Since the available far-IR photometry for five of our sources is restricted to Herschel observations, it is important to explore possible biases or systematics introduced by the lack of millimeter continuum data. Given the ever-increasing number of distant galaxies with Herschel photometry, this investigation will also provide guidance for similar studies in the future. To assess the significance of adding millimeter photometry in the derivation of the far-IR properties of the galaxies, and particularly their dust mass, we repeat the fitting procedure, this time excluding any data at wavelengths longer than rest-frame 200 μm, for the four z ∼ 1.5 BzK and the three SMGs with available millimeter continuum observations.

A comparison between the derived dust masses with and without the use of millimeter data is shown in Figure 4 (left). We find no evidence for a strong systematic bias, with an average ratio between the two estimates 〈M^mm_dust/M^nomm_dust〉 = 0.80 ± 0.2. Interestingly, the sources with the largest discrepancies are the SMGs, although the sample is too small to clearly demonstrate the existence of a systematic effect. However, the addition of millimeter data has a noticeable impact on the uncertainties of the derived M_dust estimates, which are reduced, on average, by a factor of ∼2. Studies in the local universe reach similar conclusions, with the addition that in the absence of rest-frame submillimeter data, dust masses tend to be underestimated for metal-poor galaxies. This could serve as an indirect indication that our sample mainly consists of metal-rich sources, something that we will also argue in Section 4.

Zoom In Zoom Out Reset image size

The derived γ and U_min are also in broad agreement between the two cases. However, we notice a weak systematic bias toward higher 〈U〉 values, i.e., stronger mean radiation fields, when millimeter data are considered in the fit, reflecting the fact that 〈U〉∝ L_IR/M_dust. We conclude that such trends in 〈U〉 and M_dust suggest that millimeter data can place better constraints on the diffuse ISM emission, as well as on the relative contribution of a PDR component to the total radiation field, and consequently on the M_dust of individual high-z galaxies. Finally, we note that the detailed sampling of the peak of the far-IR emission provided by Herschel data can result in robust L_IR estimates without the need of millimeter data.

Another method for deriving estimates for the dust masses and other far-IR properties, such as dust temperatures and dust emissivity indices (β), is to fit the far-IR to submillimeter SED of the galaxies with a single-temperature MBB, expressed as

$$\begin{equation} f_{\rm \nu } \propto \frac{\nu ^{3+\beta }}{e^{\frac{h\nu }{kT_{{\rm d}}}}-1}, \end{equation} \tag{ 5 }$$

where T_d is the effective dust temperature and β is the effective dust emissivity index. Then, from the best-fit model, one can estimate M_dust from the relation

$$\begin{eqnarray} &&M_{d} = \frac{S_{\nu }D^{2}_{L}}{(1+z)\,\kappa _{\rm rest}B_{\nu }(\lambda _{{\rm rest}},T_{{\rm d}})}, \;{\rm with}\;\kappa _{\rm rest}=\kappa _{\rm 0} \left(\frac{\lambda _{o}}{\lambda _{{\rm rest}}} \right)^{\beta }, \nonumber\\ \end{eqnarray} \tag{ 6 }$$

where S_ν is the observed flux density, D_L is the luminosity distance, and κ_rest is the rest-frame dust mass absorption coefficient at the observed wavelength. While this is a rather simplistic approach, mainly adopted due to the lack of sufficient sampling of the SED of distant galaxies, it has been one of the most widely used methods in the literature. Therefore, an analysis based on MBB models provides not only estimates of the effective dust temperature and dust emissivity of the galaxies in our sample but also a valuable comparison between dust masses inferred with the MBB and DL07 methods.

We fit the standard form of an MBB, considering observed data points with λ_rest > 60 μm, to avoid emission from very small grains that dominate at shorter wavelengths. For the cases where millimeter data are available we let β vary as a free parameter, while for the rest we have assumed a fixed value of β = 1.5, typical of disk-like, MS galaxies (Magdis et al. 2011; Elbaz et al. 2011). From the best-fit model, we then estimate the total M_dust through Equation (6). For consistency with the DL07 models, we adopt a value of κ_250 μm = 5.1 cm² g⁻¹ (Li & Draine 2001). To obtain the best-fit models and the corresponding uncertainties of the parameters, we followed the same procedure as for the DL07 models. The derived parameters are summarized in Table 2, and the best-fit models are shown in Figures 1–3.

A comparison between dust masses derived by DL07 and MBB models is shown in Figure 4 (right). We see that an MBB appears to infer dust masses that are lower than those derived based on DL07 models on average by a factor of ∼2 (〈M^DL07_dust/M^MBB_dust〉 = 1.96 ± 0.5). A similar trend was recently reported by Dale et al. (2012), who found that the discrepancy between the two dust mass estimates is smaller for sources with warmer S₇₀/S₁₆₀ colors. Interestingly, when we convolve the best-fit DL07 rest-frame SEDs of our galaxies with the 70 and 160 μm PACS filters, we find that the sources with the warmest S₇₀/S₁₆₀ are indeed those with the best agreement in the dust masses derived from the two methods. The reason behind this discrepancy has been addressed by Dunne et al. (2000). The grains of a given size and material are exposed to different intensities of the interstellar radiation field and thus attain different equilibrium temperatures, which will contribute differently to the SED. The single-temperature models cannot account for this range of T_d in the ISM and attempt to simultaneously fit both the Wien side of the graybody, which is dominated by warm dust, and the Rayleigh–Jean tail, sensitive to colder dust emission. The net effect is that the temperatures are driven toward higher values, consequently resulting in lower dust masses. While restricting the fit to λ ⩽ 100 μm does not provide an adequate solution to the problem (Dale et al. 2012), a two-temperature blackbody fit is more realistic and returns dust masses that are larger by a factor of ∼2 compared to those derived based on a single T_d MBB (Dunne & Eales 2001), in line with those inferred by the DL07 technique.

Despite the large uncertainties associated with the MBB analysis, it is worth noting that for all SMGs in our sample we measure an effective dust emissivity index β ∼ 2.0, in agreement with the findings of Magnelli et al. (2012). On the other hand, the far-IR to millimeter SEDs of all BzK galaxies are best described by β ∼ 1.5, similar to that found by Elbaz et al. (2011) for MS galaxies. However, we stress that the effective β does not necessarily reflect the intrinsic β of the dust grains due to the degeneracy between β and the temperature distribution of dust grains in the ISM. The β–T_d degeneracy also prevents a meaningful comparison between the derived T_d for the galaxies in our sample.

A key ingredient of this method is the metallicity of the galaxies, for which we have to rely on indirect measurements. We first derive stellar mass estimates by fitting the Bruzual & Charlot (2003) model SEDs to their rest-frame UV to near-IR spectrum. Then for the z ∼ 1.5 BzK galaxies and SMM J2135-0102 we use the M_*–Z relation at z ∼ 2 from Erb et al. (2006), while for the z ∼ 0.5 sources we adopt the fundamental metallicity relation (FMR) of Mannucci et al. (2010), which relates the SFR and the stellar mass to metallicity. To check for possible systematics, we also derive metallicity estimates for the z ∼ 1.5 sources with the FMR relation and find that the two methods provide very similar results. For GN 20 and HLSW-1, the situation is more complicated, as the FMR formula is only applicable up to z ∼ 2.5, and the Erb et al. (2006) relation saturates above M_* ∼ 10¹¹ M_☉. For these two sources, we adopt the line of reasoning of Magdis et al. (2011). Namely, in addition to the metallicity estimates based on the FMR relation, we also consider the extreme case where the huge SFR of the two sources (∼1500 M_☉ yr⁻¹) is due to a final burst of star formation triggered by a major merger that will eventually transform the galaxy into a massive elliptical. Once star formation ceases, the mass and metallicity of the resulting galaxy will not change further, and one might therefore apply the mass–metallicity relation of present-day elliptical galaxies (e.g., Calura et al. 2009). Combining the metallicity estimates based on this scenario with those derived based on the FMR relation, we estimate that the metallicities could fall in the ranges 12 + log [O/H] = 8.8–9.2 for GN 20 and 8.6–9.0 for HLSW-1. For the analysis, we will adopt 12 + log [O/H] = 9.0 ± 0.2 for GN 20 and 12 + log [O/H] = 8.8 ± 0.2 for HLSW-1. We also adopt a typical uncertainty of 0.2 for the rest of the sources. The assumed metallicities, all calibrated at the Pettini & Pagel (2004) (PP04) scale, along with the stellar masses, are given in Table 3.

To derive α_CO estimates for our sample, we first derive a δ_GDR–Z relation using the local sample presented by Leroy et al. (2011), after converting all metallicities to the PP04 scale (Figure 5, left). The data yield a tight correlation between the two quantities¹⁶ described as log δ_GDR = (10.54 ± 1.0) − (0.99 ± 0.12) × (12 + log (O/H)), with a scatter of 0.15 dex. Having obtained the M_dust and metallicity estimates, we then use the δ_GDR–Z relation to derive M_gas and subsequently estimate α_CO from the equation M_gas = α_CO × L'_CO. For the last step, we have assumed that at high-z, $M_{\rm H_{\rm 2}}\gg M_{{\rm H\,\mathsc{i}}}$ or equivalently that M_gas ≈ $M_{\rm H_{\rm 2}}$, as supported by both observational evidence (e.g., Daddi et al. 2010a; Tacconi et al. 2010; Geach et al. 2011) and theoretical arguments (Blitz & Rosolowsky 2006; Bigiel et al. 2008; Obreschkow et al. 2009). All values and their corresponding uncertainties, which take into account both the dispersion of the δ_GDR–Z relation and the uncertainties in Z and M_dust, are summarized in Table 3. In Figure 5 (right), we plot the derived α_CO values as a function of metallicity, along with the estimates of Leroy et al. (2011) for a sample of local normal galaxies, as well as for a sample of local ULIRGs from Downes & Solomon (1998), for which we were able to compute their metallicities on the PP04 scale. Fitting high-z MS galaxies along with the local sample of Leroy et al. (2011) yields

$$\begin{eqnarray} \log (\alpha _{\rm CO}) &=& -(1.39\,{\pm}\, 0.3)\,{\times}\, (12\,{+}\,\log [{\rm O/H}])_{\rm PP04}\nonumber\\ &&+12.8\,{\pm}\, 2.2, \end{eqnarray} \tag{ 8 }$$

with α_CO decreasing for higher metallicities, in agreement with previous studies (e.g., Wilson 1995; Israel 1997; Schruba et al. 2012, and references within). It is evident that MS galaxies at any redshift have higher α_CO values than those of local starbursting ULIRGs and more similar to those of local normal galaxies. In particular, for high-z MS galaxies, we find a mean α_CO of 5.5 ± 0.4, very close to the MW value of 4.6 (e.g., Strong & Mattox 1996; Dame et al. 2001). Interestingly, SMM-J2135 falls close to the locus of MS galaxies with α_CO ∼ 3, while for the remaining SMGs we find α_CO ∼ 1, close to the average value of local ULIRGs (α_CO ∼ 0.8; e.g., Solomon et al. 1997; Tacconi et al. 2008). Finally, simply as a consistency check, we derive and plot in Figure 5 (left) the M_gas/M_dust values of our targets, using α_CO estimates as obtained from Equation (8).

Zoom In Zoom Out Reset image size

Before trying to interpret our findings regarding the α_CO values, it is important to discuss possible caveats and limitations of the method. A key ingredient for the derivation of α_CO is the M_dust estimate. These estimates rely heavily on the adopted model of the grain-size distribution. Dust masses based on different grain-size compositions and opacities in the literature can vary even by a factor of ∼3, indicating that the absolute values should be treated with caution (e.g., Galliano et al. 2011). However, the relative values of dust masses derived based on the same assumed dust model should be correct and provide a meaningful comparison, as long as the dust has similar properties in all galaxies. Therefore, any trends arising from the derived M_dust estimates are essentially insensitive to the assumed dust model.

One important assumption of our technique is that the δ_GDR–Z relation defined by local galaxies does not evolve substantially with cosmic time. These quantities are very intimately related, as dust is ultimately made of metals. For example, solar metallicity of Z_☉ = 0.02 means that 2% of the gas is in metals. The local δ_GDR–Z relation suggests an average trend where 1% of the gas is in dust. This means that about 50% of the metals are locked into dust. This is already quite an efficient process of metal condensations into dust, and it is hard to think how this could be even more efficient at high z (where, e.g., less time is available for the evolution of low-mass stars into the asymptotic giant branch phase), with a much higher proportion of metals locked into the dust phase, a situation that would lead to an overestimate of M_gas with our technique. Similarly, chemical evolution studies suggest a small evolution of dust-to-metal ratio, i.e., by a maximum factor of ∼2 from 2 Gyr after the big bang up to the present time (e.g., Edmunds 2001; Inoue 2003; Calura et al. 2008). Furthermore, if we plot the ratio of the direct observables, L'_CO/M_dust, as a function of metallicity, both local and high-z galaxies follow the same trend, suggesting that the local δ_GDR–Z relation is valid at high z (Figure 6). Finally, while the Leroy et al. (2011) sample yields a slope of −1.0 in the δ_GDR–Z relation, Muñoz-Mateos et al. (2009) find a rather steeper slope of −2.45. However, adopting the Muñoz-Mateos relation (after rescaling it to PP04 metallicity) returns very similar α_CO values, as the two relations intercept at metallicities very close to the average metallicity of our sample (12 + log [O/H] ≈ 8.6). We note though that for GN 20 and HLSW (the sources with the highest Z estimates), the Muñoz–Mateos relation would reduce α_CO values by a factor of ∼3.

Zoom In Zoom Out Reset image size

With robust L_IR measurements in hand, we are in a position to derive accurate estimates of the SFR and sSFR for galaxies in our sample. We first convert the derived L_IR to SFRs by using the Kennicutt (1998) relation, scaled for a Chabrier IMF for consistency with our stellar mass estimates. We then infer the sSFR of each source and compare it to the characteristic sSFR_MS of the MS galaxies at the corresponding redshift of the source, using the methodology described in Section 2 up to z = 2.5. For the two SMGs at z ∼ 3 and z ∼ 4, we use the SFR–M_* relations of Magdis et al. (2010b) and Daddi et al. (2009), respectively. In order to conservatively identify MS galaxies, we classify SBs as galaxies with an excess sSFR relative to that of the MS galaxies by at least a factor of three, i.e., with sSFR/sSFR_MS > 3. We note, however, that the sSFR/sSFR_MS indicator is only a statistical measure of the star formation mode of the galaxies, and there is no rigid limit that separates MS from SB galaxies.

Our analysis suggests that high-z MS galaxies, associated with a secular star formation mode, exhibit higher α_CO values at any redshift compared to those of merger-driven local ULIRGs. Since the latter are known to be strong outliers from the local SFR–M_* relation, we expect a dependence of the α_CO value with sSFR. However, the functional dependence between the two quantities is not straightforward. In the local universe, there is evidence for a bimodality in the α_CO value, linked with the two known star formation modes: normal disks (for Z ∼ Z_☉) are associated with α_CO ∼ 4.4 and merger-driven ULIRGs with α_CO ∼ 0.8. As we will discuss later, if the α_CO of normal galaxies depends only on metallicity, then we expect only a weak increase of α_CO as a function of sSFR/sSFR_MS within the MS. Consequently, in order to incorporate the observed steep decrease of α_CO in SB systems, α_CO and the relative offset from the MS (i.e., sSFR/sSFR_MS) will be related with a step function. On the other hand, if α_CO also depends on other parameters than metallicity, such as the compactness and/or the clumpiness of the ISM of normal galaxies, and if the ISM conditions of MS galaxies do vary significantly, then α_CO could smoothly decrease as we depart from the MS and move toward the SB regime (e.g., Narayanan et al. 2012).

The scenario of a continuous variation of α_CO with sSFR/sSFR_MS is depicted in Figure 7(a). In addition to our galaxies, we also include a sample of local normal galaxies for which we have robust M_* and SFR estimates from Leroy et al. (2008) and α_CO measurements from Schruba et al. (2012). A linear regression fit yields the following relation:

$$\begin{equation} \alpha _{\rm CO} = (5.8\pm 2.0) \times [{{\rm sSFR}/{\rm sSFR}_{{\rm MS}}}]^{(-0.67 \pm 0.25)}, \end{equation} \tag{ 9 }$$

with a scatter of 0.3 dex. We stress that, given the lack of a sufficient number of SBs in our sample, this relation is poorly constrained in the SB regime. Nevertheless, the locus of local ULIRGs seems to be consistent with the emerging trend. Note that since we lack a sample of local ULIRGs for which both α_CO and sSFR are accurately determined, the position of local ULIRGs in Figure 7 indicates the average α_CO and sSFR of the population, as derived by Downes & Solomon (1998) and da Cunha et al. (2010b), respectively. To further check the robustness of the fit, we also perturbed the original data within the errors and repeated the fit for 1000 realizations. The mean slope derived with this method is very close to the one describing the original data. A Spearman's rank correlation test yields a correlation coefficient of ρ = −0.51, with a p-value of 0.03, suggesting a moderately significant correlation between sSFR/sSFR_MS and α_CO. However, repeating the Spearman's test, this time excluding starbursting systems, does not yield a statistically significant correlation, with ρ = −0.15 and a p-value of 0.3, indicating a very small dependence, if any, between α_CO and sSFR within the MS.

Zoom In Zoom Out Reset image size

It is therefore unclear whether α_CO varies inside the MS or if instead there is a bifurcation between the mean values appropriate for MS versus SB galaxies. If the α_CO conversion factor primarily depends on metallicity, one would expect little or no variations with sSFR at constant mass within the MS. Detailed computations for this scenario are presented in M. Sargent et al. (2012, in preparation), in which burst-like activity is assumed to be characterized by a constant, low conversion factor, chosen here to be α^SB_CO = 0.8. To predict the variations of α_CO as a function of sSFR/sSFR_MS according to this scenario, Sargent et al. (2012) compute the relative importance of MS and SB activities at a given position within the stellar mass versus SFR plane using the decomposition of the sSFR distribution at fixed stellar mass derived in Sargent et al. (2012; based on the data of Rodighiero et al. 2011). Metallicities, which form the basis for computing α_CO for an ISM experiencing "normal" star formation activity, vary smoothly as a function of stellar mass and SFR according to the calibration of the fundamental metallicity plane given in Mannucci et al. (2010; see also Lara-López et al. 2010). To compute the corresponding value of the conversion factor, M. Sargent et al. (2012, in preparation) assume a variation of α_CO∝Z⁻¹ as they find that this relation between metallicity and α_CO reproduces the faint-end slope of the z = 0 CO LF of Keres et al. (2003) best, under the set of assumptions just described.

The colored tracks in Figure 7(b) delineate the variations of α_CO with sSFR/〈sSFR〉_MS expected for galaxies in three bins of total stellar mass. Three regions can be clearly distinguished: (1) the MS locus, characterized by a gradual increase of α_CO with sSFR that is caused by the slight rise in metallicity predicted by the FMR; (2) the SB region at high sSFR/〈sSFR〉_MS ≫ 4, where the conversion factor, by construction, assumes a mass- and redshift-independent, constant value; and (3) a narrow transition region (spanning roughly sSFR/〈sSFR〉_MS ∈ [2, 4]), where α_CO drops from an approximately MW-like conversion factor to the SB value. Note that, in contrast to the previous case, α_CO varies only slightly among MS galaxies. This weak dependence of α_CO on sSFR/sSFR_MS introduces a step as we move from MS to SB systems. In the emerging picture, the majority of the star-forming population is dominated by either the MS or SB mode, and only a small fraction consists of "composite" star-forming galaxies hosting both normal and SB activities. For composite sources, the α_CO should be interpreted as a mass-weighted, average conversion factor that reflects the relative amount of the molecular gas reservoir that fuels star-forming sites experiencing burst-like and secular star formation events, respectively. Note also that variations with redshift (indicated by different colored lines in Figure 7(b)) are small for the limited range of stellar mass covered by our sample and would likely be indistinguishable within the natural scatter about the median trends plotted in the figure.

Both scenarios appear to be consistent with the data, leaving open the question of whether the transition from normal to SB galaxies is followed by a smooth or a step-like variation of α_CO. However, both scenarios agree in that α_CO should not vary much within the MS. Indeed, even for the case of continuous variation, the decrease of α_CO becomes statistically significant only for strongly starbursting systems, as Equation (9) indicates that α_CO would vary at most by a factor of ∼2 within the full range of the MS, well within the observed scatter of the relation. Therefore, the average value α_CO ∼ 5, as derived in this study, can be regarded as representative for the whole population of MS galaxies at any redshift, with stellar masses ⩾10¹⁰ M_☉. We remind the reader that because sSFR_MS increases with redshift as (1 + z)^2.95, at least out to z ≈ 2.5, galaxies with (U)LIRG-like luminosities can enter the MS at higher redshifts or at large stellar masses (e.g., Sargent et al. 2012). We stress that our analysis suggests that the α_CO of these systems, i.e., high-z MS ULIRGs, is on average a factor of ∼5 larger than that of local ULIRGs.

The determination of the α_CO value has also been the focus of several theoretical studies. In particular, Narayanan et al. (2011) investigated the dependence of α_CO on the galactic environment in numerical simulations of disk galaxies and galaxy mergers and reported a relationship between α_CO and the CO surface brightness of a galaxy. Here, we compare our observationally constrained α_CO values with those predicted by their theoretical approach. According to Narayanan et al. (2011),

$$\begin{equation} \alpha _{\rm CO} = \frac{6.75 \times 10^{20} \langle W_{\rm CO} \rangle ^{-0.32}}{6.3 \times 10^{19} \times (Z/Z_{\odot })^{0.65}}, \end{equation} \tag{ 10 }$$

where 〈W_CO〉 is the CO surface brightness of the galaxy in K Km s⁻¹. According to Solomon & Vanden Bout (2005),

$$\begin{equation} L^{\prime }_{\rm CO} = 23.5 \times \Omega _{{\rm B*S}} D^{2}_{L} I^{^{\prime }}_{\rm CO} (1 + z)^{-3}, \end{equation} \tag{ 11 }$$

where L'_CO is measured in K km s⁻¹ pc², Ω_B*S is the solid angle of the source convolved with the telescope beam measured in arcsec², and $I^{^{\prime }}_{\rm CO}$ is the observed integrated line intensity in K km s⁻¹. In the case of an infinitely good resolution, as assumed in Narayanan's simulations, this becomes the solid angle subtended by the source. The CO line luminosity can also be expressed for a source of any size in terms of the total line flux:

$$\begin{equation} L^{^{\prime }}_{\rm CO} = 3.25 \times 10^{7} S_{\rm CO} \Delta v \nu ^{-2}_{{\rm obs}} D^{2}_{L} (1 + z)^{-3}, \end{equation} \tag{ 12 }$$

where S_COΔv is the velocity integrated flux in Jy km s⁻¹. From Equations (11) and (12), and for the case of an infinitely good resolution, we derive

$$\begin{equation} W_{\rm CO} = 1.38\times 10^{6} \frac{4 \times \nu _{{\rm obs}}^{-2}}{\pi (1+z) \theta ^{-2}} \times S_{\rm CO}\Delta v, \end{equation} \tag{ 13 }$$

where ν_obs is the observed frequency and θ is the size of the source in arcsec.¹⁷ For our sources we use, where possible, the CO sizes, and UV sizes when CO observations are too noisy to reliably measure the source extent. Using Equations (10) and (13), we then get an estimate of α_CO for each source in our sample according to the prescription of Narayanan et al. (2011), and in Figure 8 we compare them to the values derived based on our method. Although there seems to be a small systematic offset (on average by a factor of ∼1.4) toward lower values for the theoretical approach, the overall agreement of the α_CO values is very good when one considers the large uncertainties and assumptions of the two independent methods. We stress that the derived 〈W_CO〉 and therefore the estimated α_CO values are very sensitive to the choice of the size indicator.

Zoom In Zoom Out Reset image size

While the sSFR probes the star formation mode of a galaxy only indirectly and from a statistical point of view, accurate measurements of its star formation efficiency, SFE = SFR/M_gas, can provide more direct and reliable answers. With the robust α_CO estimates derived in this study, we can convert the CO measurements into M_gas and subsequently infer the SFE of galaxies in our sample. In what follows we use SFE = SFR/M_gas and SFE = L_IR/M_gas interchangeably, since SFR is linearly related to L_IR through SFR = L_IR × 10⁻¹⁰ (i.e., through the Kennicutt (1998) relation, calibrated for a Chabrier IMF).

We have so far established that even though there appears to be only small variation of α_CO among MS galaxies, the α_CO of high-z MS galaxies is similar to that of local disks with lower L_IR and is ∼5 times larger than the average α_CO value observed for local ULIRGs. The direct implication of this finding is that for a given L'_CO, one could expect a large variation in the amount of molecular gas for galaxies between MS and SB galaxies. This also suggests that for a given L_IR (or equally, SFR), galaxies in the MS have lower SFEs, as compared to starbursting systems. Indeed, using the inferred α_CO values, we derive SFE estimates for each source in our sample and plot them against sSFR/sSFR_MS in Figure 9(a). We supplement our sample with a compilation of local disks (Leroy et al. 2008), as well as with a sample of local ULIRGs (Solomon et al. 1997; Rodríguez Zaurín et al. 2010).

Zoom In Zoom Out Reset image size

Although the actual dependence of SFE on sSFR/sSFR_MS within the MS and the transition from MS to the SB regime will be discussed in detail in the next section, here we can derive some crucial results about the SFE of the two populations (MS and SB galaxies) by employing a simple, empirical relation between SFE and sSFR/sSFR_MS. Indeed, a linear regression fit and a Spearman's test to our high-z data (including SBs) yield a strong and statistically significant correlation (ρ = 0.71, p-value = 0.0012) between SFE and sSFR/sSFR_MS, with

$$\begin{equation} {\rm SFE} = (8.1\pm 1.2) \times [{{\rm sSFR}/{\rm sSFR}_{{\rm MS}}}]^{1.34\pm 0.13}, \end{equation} \tag{ 14 }$$

suggesting substantially higher SFEs for galaxies with enhanced sSFR (Figure 9(a)). For MS galaxies, we find an average 〈SFE〉 = (14  ±  2) L_☉/M_☉, corresponding to an average gas consumption timescale (τ_gas = M_gas/SFR) of ∼0.7 Gyr, indicative of long-lasting star formation activity. This is in direct contrast to the short-lived, merger-driven SB episodes observed in local ULIRGs, with an average 〈SFE〉 = 200 L_☉/M_☉, ∼12 times higher than that of MS galaxies. Since the majority of galaxies at any redshift are MS galaxies (e.g., Elbaz et al. 2011) and SB galaxies seem to play a minor role in the star formation density throughout the cosmic time (e.g., Rodighiero et al. 2011; Sargent et al. 2012), our results come to support the recently emerging picture of a secular, long-lasting star formation as the dominant mode of star formation in the history of the universe.

To ensure that our result is not an artifact of the assumed α_CO values, in Figure 10 we also plot the ratio of the direct observables, L_IR/L'_CO versus the offset from the MS for the same set of objects, omitting any assumptions for α_CO. We find the two values to strongly correlate with ρ = 0.79 and a p-value of 3.4 × 10⁻⁶ and a functional relationship of

$$\begin{equation} L_{\rm IR}/L^{\prime }_{\rm CO} = 34.6(\pm 5.0) \times [{{\rm sSFR}/{\rm sSFR}_{{\rm MS}}}]^{0.80\pm 0.1}, \end{equation} \tag{ 15 }$$

with MS galaxies exhibiting lower L_IR/L'_CO values by a factor of ∼4. We note, however, that the absence of a well-defined sample of high-z SBs is striking. SMGs, which are frequently regarded as prototypical high-z SBs, have recently turned out to be more likely a mixed ensemble of objects, with a significant fraction of them exhibiting SFRs and SFEs typical of MS galaxies (Ivison et al. 2011; Rodighiero et al. 2011). Indeed, while we have only three SMGs in this study, we seem to face the same situation. Of the three SMGs considered here, only HSLW-01 appears to be a strong SB, with sSFR/sSFR_MS ∼ 10 and SFE ∼ 300 L_☉/M_☉. GN 20 is only marginally outside the MS regime, with sSFR/sSFR_MS ∼ 4, although with high SFE ∼120 L_☉/M_☉, while SMM-J2135 behaves like an MS galaxy with sSFR/sSFR_MS ∼ 2.5 and SFE ∼ 30 L_☉/M_☉, much lower than the SFE observed for local ULIRGs (>100 L_☉/M_☉). Clearly a larger, objectively selected sample of high-z SBs is essential to improve this investigation.

Zoom In Zoom Out Reset image size

5.1.1. Scenario I: A Global SF Law for MS Galaxies; f_gas as a Key Parameter

In this scenario, the physical parameter that drives the sSFR of an MS galaxy is the gas fraction, while the SFE remains roughly constant. In this case we have an SF law in the form

$$\begin{equation} \log M_{\rm gas}= \xi _{\rm gas} \times \log {\rm SFR} + C_{\rm 1}, \end{equation} \tag{ 16 }$$

where ξ_gas = 1 would imply M_gas/SFR = constant. However, there is observational evidence that ξ_gas ≈ 0.8, i.e., slightly lower than unity (e.g., Daddi et al. 2010a; Genzel et al. 2010; M. Sargent et al. 2012, in preparation), suggesting a very mild dependence of SFE from SFR (or from sSFR/sSFR_MS, for fixed stellar mass), with a slope of ≈0.2. On the other hand, dividing Equation (16) by stellar mass indicates that M_gas/M_* would vary as a function of sSFR (or equally as a function of SFR at fixed stellar mass) with

$$\begin{equation} \frac{M_{\rm gas}}{M_{\ast }} \propto \left[\frac{{\rm sSFR}}{{\rm sSFR}_{{\rm MS}}}\right]^{0.8}. \end{equation} \tag{ 17 }$$

In Figure 9(b), we explore the variations of SFE with normalized sSFR for three representative bins of mass and redshift, as predicted by M. Sargent et al. (2012, in preparation), based on the framework described above. The two factors that contribute to the evolutionary tracks are (1) the observation that at fixed SFR (or L_IR) SBs display more than an order of magnitude higher SFE than that of "normal" galaxies (e.g., Daddi et al. 2010b; Genzel et al. 2010; this study), and (2) the fact that the dependence of SFR on M_gas is slightly supra-linear, such that SFE increases with L_IR even when the MS galaxies and SBs are considered individually (i.e., Equation (16)). Point (2) is reflected in the weak increase of SFE throughout the MS and SB regimes (sSFR/〈sSFR〉_MS < 2 and sSFR/〈sSFR〉_MS > 4, respectively), while the rapid rise in SFE in the transition region is the manifestation of point (1) above. Note that the "jump" in SFE, from the normal to the SB galaxies, is not forced by the assumption of a bimodal star formation activity but naturally results from the weak increase of SFE with sSFR/sSFR_MS within the MS. Similar tracks are presented in Figure 10, considering L_IR/L'_CO instead of L_IR/M_gas.

5.1.2. Scenario II: A Varying SF Law for MS Galaxies; SFE as a Key Parameter

This scenario assumes that f_gas for MS galaxies remains constant and that the physical parameter that dictates the position of the galaxy on the SFR–M_* plane is the SFE. Since f_gas = const, then expressing SFE as a function of sSFR/sSFR_MS, i.e., at fixed stellar mass, we have

$$\begin{equation} \log {\rm SFE} = \log ({\rm SFR}) +C_{\rm 2}, \end{equation} \tag{ 18 }$$

suggesting a linear dependence of SFE on sSFR/sSFR_MS, in reasonable agreement with the slope derived by a linear fit (ξ_SFE = 1.34 ± 0.13, Equation (14)). Although a similar slope is found when the fit is redone for MS galaxies alone, excluding the SBs, a Spearman's test suggests that the correlation is statistically insignificant (ρ = 0.58, p-value = 0.13), leaving open the question of whether the trend is valid within the MS, or if it is driven solely by the SBs. A similar picture emerges when, instead of M_gas, we consider the directly observable quantity L'_CO (Figure 10). Again, in this case, excluding the SB data results in a statistically insignificant correlation. If indeed the SFE of MS galaxies follows such a steep increase as a function of sSFR/sSFR_MS, aside from the existence of parallel M_gas–L_IR SF laws, it would also imply a smooth, continuous transition from MS to SB galaxies.

Both scenarios considered here appear to be consistent with our data, and we cannot formally distinguish between them based on these individual objects. However, as we will see in the next section, the shape of the SED, both for individually elected sources and for populations averaged via stacking, can provide more definite answers.

For the case where f_gas varies within the MS (scenario I, Section 5.1.1), we have Equation (16) still holding. Also, using the z ∼ 2 stellar mass—metallicity relation of Erb et al. (2006), we get

$$\begin{equation} \log Z= \xi _{z} \times \log M_{\ast } + C_{\it 3}, \end{equation} \tag{ 20 }$$

with ξ_z ∼ 0.15 for a stellar mass range of 10 ⩽ log (M_*/M_☉) ⩽ 11. Finally, from the SFR–M_* relation, we have

$$\begin{equation} \vspace*{-2pt} \log {\rm SFR} = \xi _{\rm ms} \times \log M_{\ast } + C_{\rm 4},\vspace*{-2pt} \end{equation} \tag{ 21 }$$

with ξ_ms ≈ 0.80. Combining the above relations yields

$$\begin{equation} \vspace*{-2pt} \log \frac{{\rm SFE}}{Z} = (\xi _{\rm ms}-\xi _{\rm gas}*\xi _{\rm ms}-\xi _{\rm z}) \times \log M_{\ast }+C_{\rm 5}.\vspace*{-2pt} \end{equation} \tag{ 22 }$$

Substituting the values quoted above for the various ξ coefficients very nearly cancels out any dependence on log M_*, giving

$$\begin{equation} \log \frac{L_{\rm IR}}{M_{\rm dust}} \propto \log \frac{{\rm SFE}}{Z} \approx C_{\rm 5}, \end{equation} \tag{ 23 }$$

i.e., a negligible dependence of the dust-mass-weighted luminosity L_IR/M_dust with M_*. Expressing SFE as a function of sSFR, i.e., as a function of SFR for fixed stellar mass (and therefore metallicity), yields

$$\begin{equation} \log \frac{L_{\rm IR}}{M_{\rm dust}} = \log {\rm SFE} = (1-\xi _{\rm gas}) \times \log {\rm SFR} +C_{\rm 1}, \end{equation} \tag{ 24 }$$

suggesting a mild dependence of L_IR/M_dust with sSFR, with a slope of ∼0.2.

In the case where f_gas is constant among MS galaxies (scenario II, Section 5.1.2), where the thickness of the MS is the result of a strongly varying SFE (or equally, varying SF law), we have

$$\begin{equation} \log \frac{L_{\rm IR}}{M_{\rm dust}} = \log {\rm SFE} + C_{\rm 6} = \log {\rm SFR} + C_{\rm 7} \end{equation} \tag{ 25 }$$

since M_gas is constant at a fixed stellar mass and f_gas. This suggests a strong dependence of L_IR/M_duston sSFR/sSFR_MS. Note that L_IR/M_dust is not expected to vary with M_*, as was also the case for the previous scenario (Equation (23)).

Summarizing the analytical predictions, the two scenarios result in two distinct behaviors of L_IR/M_dust and f_gas as a function of sSFR/sSFR_MS. In the case of a weak increase of SFE within the MS (SFE ≈ const, or equally, of a global SF law for MS galaxies), we expect a variation of f_gas (Equation (17)), while L_IR/M_dust remains roughly constant within the MS (Equation (25)). The opposite trends are expected for the case where f_gas remains constant (or equally, of a varying SF law), i.e., a steep increase of SFE within the MS and a linear increase of L_IR/M_dust with sSFR/sSFR_MS (Equation (26)). We note that both scenarios predict only weak dependence of L_IR/M_dust with M_*. Furthermore, the physical meaning of L_IR/M_dust is the luminosity emitted per unit of dust mass, or the strength of the mean radiation field heating the dust, and could serve as a rough proxy of the effective dust temperature. Indeed, since, L_IR∝σT⁴_d, two sources with the same L_IR/M_dust should have similar T_d. Similarly, for a given M_dust, higher L_IR indicates higher T_d. Therefore, the above analysis also suggests a strong variation of the SED shape of the galaxies within the MS for the scenario where f_gas ≈ constant but only very little, if any, change in the shape of the SED for the case where SFE ≈ constant.

Using the stacked samples of z ∼ 2 MS galaxies in three stellar mass and four sSFR bins, described in Section 2, we can investigate what trends, if any, are present between L_IR/M_dust–M_* and L_IR/M_dust–sSFR/sSFR_MS. To derive the far-IR properties of the stacked samples, we followed the same SED fitting procedure used for the individually detected sources. However, this time we take into account the redshift distribution of the stacked sources, to account for artificial SED broadening in the far-IR and smearing of the spectral features in the mid-IR regime. Namely, we construct the whole set of DL07 models at various redshifts and create an average SED at the median redshift of the sources considered in the stack by co-adding the model SEDs at each redshift, weighted by the redshift distribution of stacked sources. The two methods return far-IR luminosities that are in good agreement, although not accounting for the redshift distribution of the stacked sources results in higher L_IR/M_dust values. This is expected as in this case the artificial broadening drives the fit, erroneously, to higher γ values, i.e., higher contribution of the PDR component, that result in lower M_dust for a given L_IR. We also validate the reliability of the inferred dust masses of the stacked samples against a possible artificial broadening of the SEDs from the dispersion of the shape of the SED of the sources included in the stacking. Namely, based on the value of 〈U〉 derived from the best fit to the stacked data, we generated 1000 artificial SEDs, adopting a scatter of 0.2 dex in 〈U〉, and assigned to each of them a redshift based on the redshift distribution of the original sample. Then, we produced an average SED that we fit in the same manner as the original data. We find that the derived 〈U〉 value for the simulated data is in very good agreement with the one obtained for the real data. The best-fit models for the various sSFR and stellar mass bins are shown in Figure 11, along with the stacked photometric points, while the inferred parameters are presented in Table 4. We note that for every sSFR bin, the derived L_IR implies an average SFR, and subsequently an average sSFR, that is very close to the initial corrected-for-extinction, UV-based sSFR estimate. This is also the first time that Herschel data provide evidence that the thickness of the MS at z = 2 is real and not due to scatter arising from uncertainties in the derivation of SFR.

Zoom In Zoom Out Reset image size

As a reference point, we also consider the large compilation of normal galaxies in the nearby universe (〈z〉 ≈ 0.05) and local ULIRGs by da Cunha et al. (2010a, 2010b). The dust mass estimates in these studies are derived based on a two-component fit for warm and cold dust, known to give results consistent with those from the DL07 models (Magrini et al. 2011). Therefore, a direct comparison with our sample is meaningful after applying a correction factor of ∼1.7 to the dust masses of the da Cunha et al. (2010a, 2010b) sample, to account for the different κ value adopted in their study.

Figure 11 reveals a remarkable similarity in the average SEDs of z ∼ 2 MS galaxies for different sSFR bins, indicating a very small variation of the SED shape of the galaxies within the MS. This is further manifested by the weak dependence of L_IR/M_dust with sSFR/sSFR_MS (inset panel), both for the z ∼ 2 sample and also for the local, normal galaxies (slope of 0.11 ± 0.04). We reach similar conclusions when considering the SEDs for different stellar mass bins: the shape of the SED and the value of L_IR/M_dust for both z ∼ 2 and z ∼ 0 MS galaxies appear to be independent of the stellar mass. The observed trends between L_IR/M_dust versus sSFR/sSFR_MS and L_IR/M_dust versus M_* are in striking agreement with those derived by our analytical approach, favoring a small variation of SFE within the MS, and therefore (1) the existence of a global SF law for MS galaxies and (2) a step-like dependence of SFE with sSFR/sSFR_MS. Indeed, the scenario of continuously increasing SFE would imply a considerable increase of L_IR/M_dust (shown in the inset panel of Figure 11 (top) with a dashed line) and a noticeable change of SED shape of galaxies within the MS, something that is not supported by our data.

Our data can also be used to infer the dependence of f_gas on sSFR and M_*. Namely, from the derived M_dust and Z values, we can use the δ_GDR–Z relation to estimate M_gas and subsequently f_gas or M_gas/M_*. In Figure 12, we plot the derived M_gas/M_* as a function of sSFR/sSFR_MS and M_* for the normal galaxies at z ∼ 2.0 and z ∼ 0. We find a clear trend of increasing M_gas/M_* with increasing sSFR for both samples, with

$$\begin{equation} \frac{M_{\rm gas}}{M_{\ast }} =(2.05\pm 0.32) \times \left[\frac{{\rm sSFR}}{{\rm sSFR}_{\rm MS}}\right]^{0.87\pm 0.15}, \end{equation} \tag{ 26 }$$

very close to the slope derived in Equation (17). This result further supports that it is variations of f_gas, rather than variations of SFE, that are responsible for the thickness of the SFR–M_* relation at any redshift. We also reveal a clear trend of decreasing f_gas with increasing stellar mass, with

$$\begin{eqnarray} \log \left(\frac{M_{\rm gas}}{M_{\ast }}\right) &=&(5.63\pm 0.39) - (0.51\pm 0.10)\nonumber\\ && \times \log (M_{\ast }), [M_{\ast }\;{\rm in}\;M_{\odot }], \end{eqnarray} \tag{ 27 }$$

in agreement with various studies that predict a similar behavior (e.g., Daddi et al. 2010a; Saintonge et al. 2011; Davé et al. 2012; Popping et al. 2012; Fu et al. 2012). Note that these trends are due to the sublinear slopes of the SFR versus M_* and M_gas versus SFR relations (see Daddi et al. 2010a).

Zoom In Zoom Out Reset image size

We recall that the M_gas derived with the δ_GDR method traces the total gas of a galaxy, i.e., $M_{\rm H_{\rm 2}}$ + $M_{{\rm H\,\mathsc{i}}}$. As we have discussed above, in our analysis we have assumed that H₂ dominates the gas mass of relatively massive (M_* > 10¹⁰ M_☉) high-z galaxies. An indirect way to investigate this assumption is to compare the inferred M_gas and L_IR of the z ∼ 2 sample to the L_IR–$M_{\rm H_{\rm 2}}$ SF law for MS galaxies as derived by Daddi et al. (2010a), based on the dynamical properties of the individually detected BzK galaxies considered in our study. If $M_{{\rm H\,\mathsc{i}}}$ contributed significantly to the total gas mass, then the galaxies would appear to show an excess in $M_{\rm H_{\rm 2}}$ with respect to L_IR. This comparison, for various subsets of our z ∼ 2 sample, is shown in Figure 13. The very good agreement between the estimated M_gas and the SF law seems to validate our assumption (i.e., negligible contribution from H i), at least for the stellar mass range considered in this study. Indeed, our M_dust measurements are (obviously) luminosity weighted and thus effectively SFR weighted. Similarly, the metallicity estimates that we adopt from the literature are luminosity/SFR weighted as well, so it is reasonable that our final estimates are most sensitive to the gas connected to the ongoing SFRs, i.e., the molecular gas.

Zoom In Zoom Out Reset image size

We have seen that MS galaxies at a given redshift tend to have uniform radiation fields, parameterized by L_IR/M_dust, and an overall small variation in shape of their SED. We also presented evidence that the mean radiation field of MS galaxies increases with redshift as ∼(1 + z)^1.15. This indicates that we could build template SEDs of MS galaxies at different redshifts. In the DL07 models, the mean radiation field 〈U〉 is an interplay between γ, which is the fraction of dust emission arising from PDRs, and U_min, which is the radiation field heating the diffuse ISM. We stress that the same 〈U〉 could arise from a different combination of γ and U_min and it does not monotonically define the shape of the SED. Although our sample is not sufficiently large to break this degeneracy, we can use the mean 〈U〉 values at each redshift to construct template SEDs of MS galaxies, under the assumption that γ does not change significantly with time. We adopt a constant γ = 0.02, similar to that reported for local galaxies by Dale et al. (2012) and Draine et al. (2007). Indeed, if we were to reproduce the same 〈U〉 with higher γ values, then the shape of the SED would become unphysically flat close to the peak. Also, the galaxies in our sample have a mean γ very close to that value, so it seems like a valid simplification. Then, using the relation between 〈U〉 and redshift derived in the previous section, we built template SEDs of MS galaxies for a grid of redshifts up to z = 2.5 using the DL07 models. Beyond z > 2.5, we assume a flattening of the evolution of 〈U〉, similar to that observed for the sSFR (e.g., González et al. 2010). For the remaining free parameter, q_PAH, we adopt the mean value of 3.19% up to redshift z ∼ 1.5 and q_PAH = 2.50% at z > 1.5 as indicated by our data. The resulting MS SEDs are shown in Figure 16.

Zoom In Zoom Out Reset image size

We note that the overall philosophy is very similar to that of Elbaz et al. (2011), who presented a single template for MS galaxies, with a unique far-IR-to-mid-IR luminosity ratio, IR8, at all redshifts up to z > 2 and a single template for starbursting systems with a significantly higher value of IR8. Here, we extend this approach but introduce some evolution of the SEDs of MS galaxies, following the increase of 〈U〉 with redshift. Our templates result in IR8 values that also mildly increase from IR8 ∼ 4.0 at z = 0 to IR8 ∼ 6.0 at z > 2, in agreement with Elbaz et al. (2011; D. Elbaz et al. 2012, in preparation) and with Reddy et al. (2012). Similar to Elbaz et al. (2011), for SB galaxies, we show here a unique template, using the best-fit SED of GN 20. Despite the simple and straightforward nature of these templates, we note that Béthermin et al. (2012a) find that they manage to reproduce recent Herschel source counts (e.g., Berta et al. 2011; Béthermin et al. 2012a), including counts measured in redshift slices. We note that to investigate the contribution of an active galactic nucleus (AGN) in our templates, we also repeated the SED fitting for the stacked samples, excluding this time photometric points that are likely to be "polluted" by an AGN activity (i.e., λ_rest < 60 μm). The inferred 〈U〉 values are very close to the values derived when considering the whole set of photometric points, in agreement with various studies that find a negligible effect of an AGN in the far-IR colors of a galaxy (e.g., Hatziminaoglou et al. 2010; Mullaney et al. 2012). We conclude that our templates should be representative for the bulk of MS galaxies, even for those that contain an average AGN contribution in the mid-IR.