© ESO, 2014

1. Introduction

Dust is an important component for understanding the galaxy formation and evolution paradigm. Dust abundance is directly connected with galaxy growth through the formation of new stars. Indeed, dust is composed of metals produced by stellar nucleosynthesis, and then expelled into the interstellar medium (ISM) via stellar winds and supernovae explosions. A fraction of these metals mixes with the gas phase, while about 30−50% (Draine et al. 2007) of them condenses into dust grains. Therefore, dust represents a consistent fraction of the total mass of metals and can be considered as a proxy for the gas metallicity. While dust is produced by the past star formation history, it also affects subsequent star formation, since it enhances the formation of molecules, hence allowing the formation of molecular clouds, out of which stars are born. Moreover, dust may affect the shape of the initial mass function (IMF), through favouring the formation of low-mass stars by fostering cloud fragmentation in low-metallicity environments and inhibiting the formation of massive stars (Omukai et al. 2005). Finally dust also affects the detectability of galaxies, because it absorbs the UV starlight and reradiates it at longer wavelengths. For all these reasons, investigating dust properties and dust evolution is a powerful diagnostic to achieve a more complete view of galaxy evolution throughout cosmic time.

With the launch of ESA’s Herschel¹ Space Observatory (Pilbratt et al. 2010), thanks to its improved sensitivity and angular resolution with respect to previous instruments, it has become possible to investigate dust properties in large samples of galaxies (e.g., Dunne et al. 2011; Buat et al. 2012; Magdis et al. 2012; Magnelli et al. 2014; Symeonidis et al. 2013, and many others). Its two imaging instruments, PACS (Poglitsch et al. 2010) and SPIRE (Griffin et al. 2010), accurately sample the far-infrared (FIR) and submillimetre dust peak from 70 to 500 μm. In this work we use the data collected by two extragalactic surveys, PACS Evolutionary Probe (PEP, Lutz et al. 2011) and Herschel Multi-tiered Extra-galactic Survey (HerMES, Oliver et al. 2012), to investigate the evolution of the dust and gas content in galaxies from the local Universe out to z ~ 2.5.

We first study how the dust content scales with the galaxy stellar content and star formation rate (SFR). Dust mass, stellar mass and SFR are essential parameters for understanding the evolution of galaxies. Since dust is formed in the atmosphere of evolved stars and in SN winds, we expect these parameters to be tightly linked with each other. The scaling relations between dust mass, stellar mass and SFR in the local or relatively nearby (z < 0.35) Universe have been investigated by recent studies based on Herschel data, such as Cortese et al. (2012) and Bourne et al. (2012). In this work, we extend the analysis to higher redshifts, and by enlarging the Herschel detected sample by means of a stacking analysis we gain enough statistics to study the correlations between the dust mass and either the stellar mass or the SFR, by keeping the other parameter fixed within reasonably small intervals. For the first time we investigate the dust scaling relations by disentangling the effects of stellar mass and those of the SFR. This resolves degeneracies associated with the so-called star formation main sequence (MS). The latter is a tight correlation observed between the SFR and the stellar mass from the local Universe out to at least z ~ 3, with a roughly 0.3 dex scatter (e.g., Brinchmann et al. 2004; Noeske et al. 2007; Elbaz et al. 2007; Santini et al. 2009; Karim et al. 2011; Rodighiero et al. 2011; Whitaker et al. 2012, and references therein). Galaxies on the MS are thought to form stars through secular processes by gas accretion from the intergalactic medium. Outliers above the MS are defined as starbursts (e.g., Rodighiero et al. 2011). Star formation episodes in these galaxies are violent and rapid, likely driven by mergers (e.g., Elbaz et al. 2011; Wuyts et al. 2011b; Nordon et al. 2012). Despite the much more vigorous star formation activity observed in starbursts, according to recent studies (e.g., Rodighiero et al. 2011; Sargent et al. 2012; Lamastra et al. 2013a), these galaxies play a minor role in the global star formation history of the Universe, accounting for only ~10% of the cosmic SFR density at z ~ 2. Since at any redshift most of the galaxies are located on the MS, most studies cannot investigate the dependence of physical quantities (e.g., dust content) on stellar mass and SFR independently, since these two quantities are degenerate along the MS. To disentangle the intrinsic dependence on each of these quantities large samples of objects are required to properly investigate the dependence on SFR at any fixed M_star and, vice versa, the dependence on M_star at a fixed SFR.

Knowledge of the dust content can be further exploited to obtain information on the gas content, if the dust-to-gas ratio is known. In the past, most studies on the gas content in high-z galaxies have been based on CO observations (e.g. Tacconi et al. 2010, 2013; Daddi et al. 2010; Genzel et al. 2010). These studies have allowed the investigation of the relation between the molecular gas mass and the SFR, i.e., the Schmidt-Kennicutt relation (S-K, Schmidt 1959; Kennicutt 1998), at different cosmic epochs. However, these observations are time consuming and affected by uncertainties associated with the CO-to-H₂ conversion factor, which is poorly constrained for starburst or metal-poor galaxies (see Bolatto et al. 2013, for a review).

An alternative method to derive the gas content is to exploit the dust masses inferred from FIR-submm measurements and convert them into gas masses by assuming a dust-to-gas ratio (e.g., Eales et al. 2010; Leroy et al. 2011; Magdis et al. 2011; Scoville 2012). We adopt this approach in the second part of this work. We convert the dust mass into gas mass by assuming that the dust-to-gas ratio scales linearly with the gas metallicity and that dust properties are similar to those in the local Universe, where the method is calibrated. We estimate the gas metallicity from our data by exploiting the fundamental metallicity relation (FMR) fitted by Mannucci et al. (2010) on local galaxies and shown to hold out to z ~ 2.5. According to the FMR, the gas metallicity only depends on the SFR and the stellar mass, and does not evolve with redshift (see also Lara-López et al. 2010). With these assumptions, which will be discussed in the text, we study the relation between the SFR and the gas mass and investigate the evolution of the gas fraction out to z ~ 2.5 independently of CO measurements. We note, however, that the two methods for measuring the gas mass (the “dust-method” and CO observations) are cross-calibrated with each other.

A similar approach was adopted by Magdis et al. (2012) by using Herschel data from the GOODS-Herschel survey. We improve over their work by also using the data in the COSMOS field that, thanks to the large number of objects, allows us to greatly expand the stacking technique to a range of galaxy physical parameters not explored by Magdis et al. (2012), and to significantly shrink the uncertainties. Moreover, while Magdis et al. (2012) bin the data in terms of their distance from the MS at any redshift, we bin our data in stellar mass, SFR and redshift, to avoid the inclusion of any a-priori relation between stellar mass and SFR and to study the existing trend as a function of physical parameters.

The paper is organized as follows. After presenting the data set (Sect. 2) and the method used to compute SFR, stellar, dust and gas masses, and gas metallicities (Sect. 3), we present the dust scaling relations in Sect. 4, and the study of the evolution of the gas content in Sect. 5. Finally, we summarize our results in Sect. 6.

In the following, we adopt the Λ-CDM concordance cosmological model (H₀ = 70 km s^-1 Mpc^-1, Ω_M = 0.3 and Ω_Λ = 0.7) and a Salpeter IMF.

2. The data set

For this work we take advantage of the wide photometric coverage available in three extragalactic fields: the two deep GOODS fields (GOODS-S and GOODS-N, ~17′ × 11′ each) and the much larger but shallower COSMOS field (~85′ × 85′). Dealing with these fields together represents an excellent combination of having good statistics on both bright and faint sources from low to high redshift.

Most important for the aim of this work, i.e., essential to derive dust masses, are the FIR observations carried out by Herschel with the shorter wavelength (70, 100 and 160 μm) PACS camera and the longer wavelength (250, 350, 500 μm) SPIRE camera. As anticipated in Sect. 1, we use the data collected by the two extragalactic surveys PEP and HerMES. Catalogue extraction on Herschel maps is based on a PSF fitting analysis that makes use of prior knowledge of MIPS 24 μm positions and fluxes. PACS catalogues are described in Lutz et al. (2011, and references therein) and Berta et al. (2011), while SPIRE catalogues are presented in Roseboom et al. (2010) and are updated following Roseboom et al. (2012). The 3σ limits² at 100, 160, 250, 350 and 500 μm are 1.2, 2.4, 7.8, 9.5, 12.1 mJy in GOODS-S, 3.0, 5.7, 9.2, 12.0, 12.1 mJy in GOODS-N and 5.0, 10.2, 8.1, 10.7, 15.4 mJy in COSMOS, respectively. The only field which was observed at 70 μm is GOODS-S. After testing that the use of 70 μm photometry does not introduce any significant difference in the dust mass estimates, we ignored this band for consistency with the other fields.

In order to infer redshifts and other properties needed for this study, we complement Herschel observations with public multi-wavelength photometric catalogues. For GOODS-S we use the updated GOODS-MUSIC catalogue (Santini et al. 2009; Grazian et al. 2006). For GOODS-N we use the catalogue compiled by the PEP Team and described in Berta et al. (2010) and Berta et al. (2011)³. For the COSMOS field we use the multi-wavelength catalogue presented in Ilbert et al. (2009) and McCracken et al. (2010)⁴. COSMOS data reduction is described in Capak et al. (2007), although the new catalogue uses better algorithms for source detection and photometry measurements. This catalogue is supplemented with IRAC photometry from Sanders et al. (2007) and Ilbert et al. (2009) and 24 μm photometry from Le Floc’h et al. (2009).

All the catalogues are supplemented with either spectroscopic or photometric redshifts. Spectroscopic redshifts are available for ~30%, ~27% and ~3% of the final sample, respectively, in GOODS-S, GOODS-N and COSMOS. For the remaining sources, we adopt the photometric redshift estimates publicly released with the two GOODS catalogues and those computed by the authors for COSMOS and presented in Berta et al. (2011). The latter were computed for all sources rather than for the I-selected subsample released by Ilbert et al. (2009), and show similar quality for the objects in common. Photometric redshifts in GOODS-S are estimated by fitting the multi-wavelength photometry to the PEGASE 2.0 templates (Fioc & Rocca-Volmerange 1997), as presented in Grazian et al. 2006 and updated as in Santini et al. 2009. For GOODS-N and COSMOS, the EAZY code (Brammer et al. 2008) was adopted, as discussed in Berta et al. (2011). We refer to the papers cited above, as well as to Santini et al. (2012b) for more detailed information about spectroscopic and photometric redshifts and their accuracy.

2.1. Sample selection

In order to achieve a reliable estimate of the main physical parameters required for this analysis, we need to apply some selections to the galaxy sample in the three fields.

We firstly require the signal-to-noise ratio in K band to be larger than 10. This selection ensures clean photometry and reliable stellar mass estimates for all sources.

Secondly, in order to estimate the SFR from an IR tracer, independent of uncertain corrections for dust extinction, we require a 24 μm detection for all galaxies (see Sect. 3.2). This is the tightest selection criterion and limits the final sample to galaxies with relatively high star formation (32–52% of the sample, depending on the field). However, although it reduces the dynamical range probed, a SFR cut is not an issue for most of this study, since we analyse trends as a function of SFR or at fixed SFR. In the latter case, the use of narrow SFR intervals prevents strong incompleteness effects within each individual bin.

Finally, we remove all known AGNs from the catalogues (~2.5% of the total final sample), by considering X-ray detected sources (the AGN sample of Santini et al. 2012b), highly obscured AGNs detected through their mid-IR (MIR) excess (following Fiore et al. 2008), and IRAC selected AGNs (Donley et al. 2012). Indeed, besides the cold dust heated by star formation regions, these sources host a warm dust component, which is heated by nuclear accretion processes and which might bias the dust mass estimates. The dust content in AGNs will be investigated in a dedicated forthcoming paper (Vito et al., in prep.).

3. Parameters determination

We describe in this section how the basic ingredients of our analysis, i.e., stellar masses (M_star), SFR, dust masses (M_dust), gas masses (M_gas) and gas metallicities, are obtained.

3.1. Stellar masses

Stellar masses are estimated by fitting observed near-UV to near-IR photometry with a library of stellar synthetic templates (e.g. Fontana et al. 2006). We adopt the same procedure described in Santini et al. (2009): we perform a χ² minimization of Bruzual & Charlot (2003) synthetic models, parameterizing the star formation histories as exponentially declining laws of timescale τ and assuming a Salpeter⁵ IMF. Age, gas metallicity, τ and reddening are set as free parameters, and we use a Calzetti et al. (2000) or SMC extinction curve (whichever provides the best fit). We refer to Santini et al. (2009) and references therein for more details on the stellar template library. In the fitting procedure, each band is weighted with the inverse of the photometric uncertainty. Since Bruzual & Charlot (2003) models do not include emission from dust reprocessing, we fit the observed flux densities out to 5.5 μm rest-frame. The redshift is fixed to the photometric or spectroscopic one, where available.

To ensure reliable stellar mass estimates, in the following we remove all sources with a reduced χ² larger than 10 (~4–13% of the final sample, depending on the field).

Our sample spans a large redshift interval, hence the range of rest-frame wavelengths used to measure stellar masses is not the same for all sources. More specifically, high-z galaxies lack constraints at the longest rest-frame wavelengths. However, Fontana et al. (2006) have shown that the lack of IRAC bands when estimating the stellar mass from multi-wavelength fitting, while producing some scatter, does not introduce any systematics (see also Mitchell et al. 2013). In any case, the rest-frame K band, essential for a reliable stellar mass estimate, is sampled even at the highest redshifts probed by our analysis.

3.2. SFR

Star formation rates are estimated from the total IR luminosity integrated between 8 and 1000 μm (L_IR) and taking into account the contribution from unobscured SF. We use the calibrations adopted by Santini et al. (2009, see references therein): $\begin{array}{ccc}\mathit{SFR}\mathrm{[}{\mathit{M}}_{\mathrm{\odot }}\mathit{/}\mathrm{yr}\mathrm{]}& \mathrm{=}& \mathrm{1.8}\mathrm{\times }{\mathrm{10}}^{-10}\mathrm{\times }{\mathit{L}}_{\mathrm{bol}}\mathrm{\left[}{\mathit{L}}_{\mathrm{\odot }}\mathrm{\right]}\mathrm{;}\\ {\mathit{L}}_{\mathrm{bol}}& \mathrm{=}& \mathrm{2.2}\mathrm{\times }{\mathit{L}}_{\mathrm{UV}}\mathrm{+}{\mathit{L}}_{\mathrm{IR}}\mathit{.}\end{array}$(1)Here L_UV = 1.5 × νL_ν(2700 Å) is the rest-frame UV luminosity derived from the SED fitting and uncorrected for extinction.

Since Herschel detections are only available for ~11–25% (depending on the field) of the sample⁶, in order to have a consistent SFR estimate for a larger number of sources, we estimate L_IR from the 24 μm MIPS band (reaching 3σ flux limits of 20 and 60 μJy in the GOODS fields and in COSMOS, respectively). Most importantly, this approach also avoids any degeneracy with the dust mass estimates, derived from Herschel data. We fit 24 μm flux densities to the MS IR template derived by Elbaz et al. (2011) on the basis of Herschel observations. This template, thanks to an updated treatment of the MIR-to-FIR emission, overcomes previous issues related with the 24 μm overestimate of L_IR and provides a reliable estimate of the SFR for all galaxies (see Fig. 23 of Elbaz et al. 2011). As a further confirmation, in Appendix B we compare the 24 μm-based SFRs with those derived by fitting the full FIR photometry and find very good agreement. This test proves that the adoption of 24 μm-based SFR does not introduce relevant biases in the analysis. Most importantly, it provides a SFR estimate that is independent of the dust and gas mass measurement and therefore allows us to confidently investigate correlations among these quantities.

3.3. Stacking procedure

Dust masses are computed by means of Herschel observations. Only a small fraction of the sources are individually detected by Herschel, and only less than 10%⁷ fulfill the requirements of good FIR sampling adopted for the dust mass estimate (see Sect. 3.5). Therefore, a stacking procedure to estimate the average flux of a group of sources is needed to perform an analysis, which is unbiased towards the brightest IR galaxies. We describe here how average fluxes for subsamples of sources are estimated. In the next section we explain how such subsamples are compiled.

The stacking procedure adopted in this work is similar to that described by Santini et al. (2012b) and also used in Rosario et al. (2012) and Shao et al. (2010). First of all, in each Herschel band we restrict to the area where the coverage (i.e., integration time) is larger than half its value at the centre of the image. This removes the image boundaries where stacking may be less reliable due to the larger noise level. For each z–M_star–SFR bin containing at least 10 sources and for each Herschel band, we stack⁸ on the residual image (i.e., map from which all 3σ detected sources have been subtracted) at the positions of undetected sources (by “undetected” we mean below 3σ confidence level). Each stamp is weighted with the inverse of the square of the error map. The photometry on the stacked PACS images is measured by fitting the PSF, while for SPIRE images we read the value of the central pixel (SPIRE maps are calibrated in Jy/beam), which was suggested by Béthermin et al. (2012) to be more reliable in the case of clustered sources. Uncertainties in the stacked flux densities are computed by means of a bootstrap procedure. The final average flux density $\overline{\mathit{S}}$ is obtained by combining the stacked flux (S_stacked) with the individually detected fluxes (S_i) in the same bin: $\overline{\mathit{S}}\mathrm{=}\frac{{\mathit{S}}_{\mathrm{stacked}}\mathrm{\times }{\mathit{N}}_{\mathrm{stacked}}\mathrm{+}{\sum }_{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\det }}{\mathit{S}}_{\mathit{i}}}{{\mathit{N}}_{\mathrm{tot}}}\mathit{,}$(2)where N_stacked, N_det and N_tot are the number of undetected, detected and total sources, respectively, in the bin.

The stacking procedure implicitly assumes that sources in the image are not clustered. However, in the realistic case sources can be clustered with other sources either included or not included in the stacking sample. This effect may result in an overestimation of the flux in blended sources (see, e.g., Béthermin et al. 2012; or Magnelli et al. 2014). Given the lack of information on sources below the noise level, it is not straightforward to correct for this effect. However, if we are able to recognize its occurrence, we can ignore the bins where the stacking is affected by confusion. For this purpose, an ad hoc simulation has been put into place by the PEP Team. We briefly recall the basic steps of the simulation, and refer the reader to Magnelli et al. (2014) for a more detailed description. Synthetic SPIRE fluxes were estimated through the MS template of Elbaz et al. (2011), given the observed redshifts and SFRs, and simulated catalogues and maps were produced. Whenever we stack on a group of sources on real SPIRE maps, we also stack at the same positions on the simulated maps and obtain a simulated average flux density (${\overline{\mathit{S}}}_{\mathrm{sim}}$). We compare ${\overline{\mathit{S}}}_{\mathrm{sim}}$ with the mean value (${\overline{\mathit{S}}}_{\mathrm{input}}$) of the same flux densities contained in the simulated catalogue (previously used to create the simulated maps). Following Magnelli et al. (2014), if $\mathrm{|}{\overline{\mathit{S}}}_{\mathrm{input}}\mathrm{-}{\overline{\mathit{S}}}_{\mathrm{sim}}\mathrm{|}\mathit{/}{\overline{\mathit{S}}}_{\mathrm{input}}\mathit{>}\mathrm{0.5}$ we reject the corresponding bin⁹. The largest blending effects are seen at low flux densities and in the 500 μm band, as expected. The criterion above implies rejection of ~10% of the stacked fluxes at 250 μm, ~16% at 350 μm and ~33% at 500 μm. We also run our analysis by including these bins, to check that their rejection does not bias our results.

3.4. The z–M_star–SFR grid and combination across fields

The basis of our stacking analysis is to infer an average dust mass for sources showing similar properties. To this aim, we divide the redshift-stellar mass-SFR parameter space into small bins, and run the stacking procedure on all galaxies belonging to each bin. The ranges covered by our grid are 0.05–2.5 in redshift, 9.75–12 in log M_star [M_⊙] and − 0.75–3 in log SFR [M_⊙/yr]. The boundaries of the bins, listed in Tables A.1 together with the abundance of sources per bin, are chosen to provide a fine sampling of the M_star–SFR parameter space and at the same time to have good statistics in each bin. We adopt bins of 0.25 dex in M_star and 0.2 dex in SFR at intermediate M_star and SFR values, where we have the best statistics, and slightly larger bins at the boundaries. This choice strongly limits the level of incompleteness within each individual bin. Incompleteness issues will simply result into bins not populated and therefore missing from our grid (e.g., at low M_star and SFR as redshift increases).

To combine the different fields, we stack on them simultaneously by weighting each stamp with the relative weight map. The total number of sources in each bin and the contribution of each field are reported in Tables A.1. Since the statistics are strongly dominated by the COSMOS field, we do not expect intrinsic differences among the fields to significantly affect our results.

For each bin of the grid we compute the average redshift, M_star and SFR of the galaxies belonging to it, and associate these values to the bin. The standard deviations of the distribution of these parameters within the bin provide the error bars associated with the average values.

3.5. Dust masses

For a population of dust grains at a given temperature and with a given emissivity, the dust mass can be inferred from their global thermal infrared greybody spectrum and, in particular, by its normalization and associated temperature. More generally, the dust thermal emission in galaxies is composed by multiple thermal components. In order to account for this, we use, as a description of the dust emission, the spectral energy distribution (SED) templates of Draine & Li (2007). In doing so, we implicitly assume that the dust properties and emissivities of our sources are similar to those of local galaxies, on which the templates were tested (Draine et al. 2007). Such assumption is supported by the lack of evolution in the extinction curves, at least out to z ~ 4 (Gallerani et al. 2010). It is also supported by the gas metallicity range probed by our sample (≥8.58, see Sect. 3.6) and by the recent results of Rémy-Ruyer et al. (2013), claiming that the gas metallicity does not have strong effects on the dust emissivity index. Moreover, our sample is mostly made of MS galaxies. The Draine & Li (2007) model is also based on the assumption that dust is optically thin, plausibly applicable to our sample, which does not include very extreme sources such as local ULIRGs or high-z sources forming a few thousands of solar masses per year. However, as a sanity check, we also have used the GRASIL model (Silva et al. 1998) which includes extreme optically thick young starburst components, and the final results are unaffected (see below). Finally, Galliano et al. (2011), by studying the Large Magellanic Cloud, found that dust masses may be systematically understimated by ≃50% when computed from unresolved fluxes. The authors ascribe this effect to possible vealing of the cold dust component by the emission of the warmer regions. However, this effect would only introduce an offset without modifying the main results of this analysis.

According to the Draine & Li (2007) model, the interstellar dust is represented as a mixture of amorphous silicate and graphite grains, with size distribution modelled by Weingartner & Draine (2001) and updated as in Draine & Li (2007), mimicking different extinction curves. A fraction q_PAH of the total dust mass is contributed by PAH particles (with <1000 C atoms). Although they only provide a minor contribution to the total dust mass, their abundance has an important effect in shaping the galaxy SED at short wavelengths. The majority (a fraction equal to 1 − γ) of dust grains are located in the diffuse ISM and heated by a diffuse radiation field contributed by many stars. This results in a single radiation intensity U = U_min, where U is a dimensionless factor normalized to the local ISM. The rest of the grains are localized in photodissociation regions close to bright stars, and exposed to multiple and more intense starlight intensities (U_min < U < U_max) distributed as a power law (∝U^− α).

Following the prescriptions of Draine et al. (2007), we build a library of MW-like models with PAH abundances q_PAH in the range 0.47–4.58%, 0.0 < γ < 0.3, α = 2, U_max = 10⁶ and U_min varying between 0.7 and 25. This latter prescription (instead of using U_min ≥ 0.1) prevents the risk of fitting erroneous large dust masses in the absence of rest-frame submillimetre data to constrain the amount of cool dust.

Dust masses are derived by fitting and normalizing the stacked 100-to-500 μm Herschel photometry to this template library. The redshift is fixed to the mean redshift in the bin. The template showing the minimum χ² is chosen, and the normalization of the fit provides a measure of the dust mass.

In the fitting procedure, we require the stacked fluxes to have at least 3σ significance. In order to have a good sampling of the spectrum, especially on the Rayleigh-Jeans side, most sensitive to the dust mass, we only consider bins with significant flux in at least 3 bands, of which at least one is longward of rest-frame 160 μm (Draine et al. 2007). This enables to account for changes in the dust temperature and makes us confident of the resulting M_dust. 26% of the total number of bins are rejected because of these selections. We visually inspect every single bin to ensure the quality of the fits, and conservatively reject 5 of them (~4%), where the stacked fluxes were not satisfactorily reproduced by the best-fit template. An example of our fitting output can be seen in Fig. 1. The best-fits for all bins in the final sample can be seen in Appendix D.

MIR fluxes are not included in the fit so that the dust mass and SFR estimates are totally independent. As a consistency check, we also computed M_dust by including 24 μm flux densities. The resulting dust masses are in very good agreement with our reference estimates: their mean (median) ratio ($\log \mathrm{(}{\mathit{M}}_{\mathrm{dust}}^{\mathrm{24}\mathit{\mu }\mathrm{m}}\mathit{/}{\mathit{M}}_{\mathrm{dust}}^{\mathrm{no}\mathrm{24}\mathit{\mu }\mathrm{m}}\mathrm{)}$) is − 0.001 (0.008), with a scatter of 0.07, and the average error bar (see below) is only ~10% lower than without including the 24 μm band. This ratio shows no trends with either stellar mass, SFR or redshift, except a slightly larger scatter at low-z (here rest-frame wavelengths below 100 μm are not sampled in the absence of 24 μm data).

Errors on M_dust are estimated by allowing the stacked photometry to vary within its uncertainty and the redshift to move around the mean value within its standard deviation in the bin. The uncertainty is given by the minimum and maximum M_dust allowed by templates whose probability according to a χ² test is larger than 32%. All data points whose associated error on M_dust is larger than 1 dex (further ~5% of the available bins) are ignored throughout the analysis, being unable to contribute in understanding the existing trends and only making the plots more crowded without adding information. After all selections, we end up with 122 data points sampling the z–M_star–SFR grid (see Fig. 4). Dust masses measured in each bin of our grid are listed in Tables A.1.

Our dust masses are in very good agreement with those computed by Magnelli et al. (2014) with the same recipe.

In addition to using the Draine & Li (2007) templates, we also fit our data with a library extracted from the chemospectrophotometric model GRASIL (Silva et al. 1998), tested to reproduce the small galaxy sample of Santini et al. (2010), and with a simple modified blackbody, which assumes a single-temperature dust distribution. For consistency with the Draine & Li (2007) model, we build a modified blackbody with emissivity index β = 2 and absorption cross section per unit dust mass at 240 μm of 5.17 cm²/g (Li & Draine 2001; Draine & Lee 1984). We find that the simplified assumption of single-temperature leads to dust masses which are lower by a factor of ~1.5 compared to those obtained with the more realistic assumption of a multi-temperature grain distribution (in agreement with previous studies, e.g., Santini et al. 2010; Magnelli et al. 2012a,b; Dale et al. 2012; Magdis et al. 2012). Indeed, the attempt of reproducing the Wien side and at the same time the Rayleigh-Jeans side of the modified blackbody spectrum has the effect of overestimating the dust temperature and hence underestimating the dust mass, inversely proportional to the blackbody intensity. However, the recent work of Bianchi (2013) ascribes such disagreement to possible inconsistencies in the treatment of dust emission properties between the two approaches. The GRASIL library fits dust masses larger than the Draine & Li (2007) templates by a factor of 1.5 on average. A direct comparison between the parameters assumed by the two models is not possible, since GRASIL computes dust emission by considering the physical properties of each single grain, instead of assuming an average emissivity. One reason for the discrepancy could be that the optically thin assumption of Draine & Li (2007) is not always verified (even if true, this would not affect our results, which would be simply offset). The GRASIL library adopted, however, has not been tested to work in the absence of submillimeter data. Both the fit with GRASIL and with the modified blackbody provide χ² values that are a factor of 1.5–2 larger than the Draine & Li (2007) library. For these reasons we decided to use the dust masses obtained from the Draine & Li (2007) templates. We will expand the GRASIL library by enlarging the parameter space to better reproduce our galaxies in a future analysis. In any case, we note that the effect of choosing one dust model or the other only produces an offset, leaving the main trends outlined below almost unchanged.

3.6. Gas metallicities and gas masses

It is possible to take a further step forward with respect to observables directly measurable from our data and compute gas masses by converting dust masses through the dust-to-gas ratio (e.g., Eales et al. 2010). In order to do that, we need to make some assumptions.

We first assume that the gas metallicity is described by the FMR of Mannucci et al. (2010). The FMR is a 3D relation between gas metallicity¹⁰, stellar mass and SFR, with a very small scatter (0.05 dex). We assume that it does not evolve from the local Universe to z ~ 2.5, as confirmed by a number of recent works (e.g. Mannucci et al. 2010; Cresci et al. 2012; Nakajima et al. 2012; Henry et al. 2013a,b; Belli et al. 2013). More recently, Bothwell et al. (2013) have shown that the FMR is likely a by-product of a more fundamental relation, between H i gas mass, stellar mass and metallicity (H i-FMR). However, it is beyond the scope of the paper to discuss the origin of this relation. Given the average M_star and SFR in each bin of our grid, following Mannucci et al. (2010), we compute the gas metallicity from the linear combination μ_0.32 = log M_star − 0.32log SFR, after converting to a Chabrier IMF (as they adopt) both stellar masses ($\log {\mathit{M}}_{\mathrm{star}}^{\mathrm{Cha}}\mathrm{=}\log {\mathit{M}}_{\mathrm{star}}^{\mathrm{Sal}}\mathrm{-}\mathrm{0.24}$, Santini et al. 2012a) and SFR (log SFR^Cha = log SFR^Sal − 0.15, Davé 2008), using their Eqs. (4) and (5) and the extrapolation for low μ_0.32 values published in Mannucci et al. (2011). The inferred gas metallicities are in the range 8.58–9.07, with a scatter of 0.14 dex around the mean value of ~8.9.

We note that the FMR has not been tested over the entire SFR range studied in this work on large galaxy samples, so the extrapolation to SFRs larger than ~100  M_⊙/yr might in principle result in gas metallicity estimates that are incorrect. Moreover, the detailed shape of the FMR is matter of debate (e.g. Yates et al. 2012; Andrews & Martini 2013). For these reasons, we also tested the robustness of our results by adopting the redshift-dependent mass-metallicity relations published by Maiolino et al. (2008) and verified that all our results are independent of the specific description of the gas metallicity.

As suggested by previous studies, focused either on local (e.g., Draine et al. 2007; Leroy et al. 2011; Smith et al. 2012; Corbelli et al. 2012; Sandstrom et al. 2013), z < 0.5 (e.g., James et al. 2002) or high-z galaxies (e.g., Zafar & Watson 2013; Chen et al. 2013, Cresci et al., in prep.), we consider that a fixed fraction of metals are incorporated in dust. Within the metallicity range probed by our sample, this is true within 0.3 dex at most. Following the parameterization provided by Draine et al. (2007), we assume that the dust-to-gas ratio (δ_DGR) scales linearly with the oxygen abundance through the constant factor k_DGR: $\begin{array}{ccc}{\mathit{\delta }}_{\mathrm{DGR}}& \mathrm{=}& {\mathit{k}}_{\mathrm{DGR}}\mathrm{\times }\mathrm{\left(}\mathrm{O}\mathit{/}\mathrm{H}\mathrm{\right)}\mathrm{=}\mathrm{0.01}\mathrm{\times }\mathrm{\left(}\mathrm{O}\mathit{/}\mathrm{H}\mathrm{\right)}\mathit{/}\mathrm{(}\mathrm{O}\mathit{/}\mathrm{H}{\mathrm{)}}_{\mathrm{MW}}\\ & \mathrm{=}& \end{array}$(3)where Z = 12 + log (O/H) is the gas metallicity and Z_⊙ = 8.69 is the Solar value (Allende Prieto et al. 2001; Asplund et al. 2009). We find almost identical results from our analysis if we apply the linear relation between log δ_DGR and gas metallicity inferred by Leroy et al. (2011).

The universality of the depletion factor of metals into dust is outlined by the recent work of Zafar & Watson (2013). According to their analysis, the dust-to-metal ratio can be considered universal, independent of either column density, galaxy type or age, redshift and metallicity. However, De Cia et al. (2013) claim that the dust-to-metal ratio is significantly reduced with decreasing gas metallicity at Z < 0.1  Z_⊙ and low column densities. Yet, this should not be a concern for our analysis, since our sample does not include such low-metallicity galaxies. In a more recent paper, Chen et al. (2013) combine constraints on the dust-to-gas ratio of lensed galaxies, GRBs and quasar absorption systems, and find support for a simple, linear universal relation between dust-to-gas ratio and metallicity.

The total gas mass (atomic + molecular, M_gas hereafter) can be computed as ${\mathit{M}}_{\mathrm{gas}}\mathrm{=}{\mathit{M}}_{\mathrm{dust}}\mathit{/}{\mathit{\delta }}_{\mathrm{DGR}}\mathit{.}$(4)We can finally compute the gas fraction (f_gas hereafter) as ${\mathit{f}}_{\mathrm{gas}}\mathrm{=}{\mathit{M}}_{\mathrm{gas}}\mathit{/}\mathrm{(}{\mathit{M}}_{\mathrm{gas}}\mathrm{+}{\mathit{M}}_{\mathrm{star}}\mathrm{)}\mathit{.}$(5)The dust content, typically negligible with respect to the gas and stellar mass components (M_dust ≲ 0.01 M_star, see below), is ignored in the computation of f_gas.

4. Dust scaling relations

In this section we investigate the correlations between M_star, SFR and M_dust, and their evolution with redshift.

4.1. Dust content vs. SFR

Figure 2 shows the relation between the SFR and the dust content for galaxies of different M_star at different redshifts. A correlation between the dust content and the star formation activity is evident at all M_star and at all redshifts, although with some scatter, while no clear effect is observed with the stellar mass, with bins of different M_star sometimes overlapping (see also next section).

Before discussing the interpretation of this correlation, we stress here that, not only M_dust and the SFR are estimated from different observed flux densities (Herschel and 24 μm bands, respectively) to avoid any possible degeneracy and with intrinsically independent methods, but also they are not expected to be correlated by definition. The SFR (although in our case measured from 24 μm observations) is in principle linked to the integrated IR luminosity, i.e., it is linked to the normalization of the FIR spectrum. The dust mass comes from a combination of the template normalization and temperature(s), which determines the shape; since the template library that we have used contains multiple heating source components, the dust mass is not trivially proportional to the SFR, though related to it through the dust temperature. To verify that any observed correlation is physical and not an obvious outcome of the relation between correlated variables, we run a simulation that is described in Appendix C, showing that, by starting from a completely random and uncorrelated distribution of dust masses and SFRs, our method does not introduce any artificial correlation.

The correlation observed in Fig. 2 primarily tells us that the dust temperature plays a secondary role. The SFR–M_dust correlation is clearly a consequence of the S-K law, linking the SFR to the gas content. Indeed, as shown in Sect. 3.6, the dust mass is related to the gas mass by means of the dust-to-gas ratio. In other words, we expect the dust mass to be roughly proportional to the gas mass, with the gas metallicity introducing minor effects through the dust-to-gas ratio. Before converting dust masses into gas masses by adopting the appropriate dust-to-gas ratio in the next section, in order to represent the S-K relation on a SFR vs. M_dust plot, for the moment we assume a constant dust-to-gas ratio for all galaxies. By using Eq. (4), the S-K law (in its integrated¹¹ version inferred by Daddi et al. 2010 for local spirals and z ~ 2 BzK galaxies, Daddi et al. 2004) can be written in terms of SFR as a function of M_dust as $\log \mathit{SFR}\mathrm{[}{\mathit{M}}_{\mathrm{\odot }}\mathit{/}\mathrm{yr}\mathrm{]}\mathrm{=}\mathrm{1.31}\mathrm{\times }\log \left(\frac{{\mathit{M}}_{\mathrm{dust}}\mathrm{\left[}{\mathit{M}}_{\mathrm{\odot }}\mathrm{\right]}}{{\mathit{\delta }}_{\mathrm{DGR}\hspace{0.17em}\mathrm{\odot }}}\right)\mathrm{+}\mathrm{7.80}\mathit{,}$(6)where the last term includes the factor (1.8 × 10^-10) used to convert the total infrared luminosity (the original quantity in the expression given in Daddi et al.) into SFR, as well as the offset of 0.15 needed to convert from a Chabrier to a Salpeter IMF (see Sect. 3.6), and δ_{DGR    ⊙} is the dust-to-gas ratio computed from Eq. (3) by assuming a constant Solar metallicity. The dashed line in Fig. 2 shows the inferred S-K relation on the SFR–M_dust diagram.

Our observational points follow reasonably well the trend expected from the S-K law, with some scatter and a systematic trend (flatter slope) at high-z. We will discuss this in Sect. 5.1, where we also account for the variation of the metallicity (hence the variation of the dust-to-gas ratio as a function of metallicity).

4.2. Dust vs. stellar mass content

We plot in Fig. 3 the dust mass as a function of the stellar mass in bins of redshift. When galaxies are separated according to their SFR (coded with different colours), the correlation found by previous authors (e.g., at low redshift by Bourne et al. 2012) becomes much flatter and sometimes even disappears, hinting that this correlation is at least partly an indirect effect driven by other phenomena. More specifically, the M_dust–M_star correlation is partly a consequence of the M_dust–SFR correlation, reported in the previous section, combined with the MS, i.e., the relation between SFR and M_star. When all SFR are combined together, the low mass bins are dominated by low SFR (as a consequence of the MS), which are associated with low M_dust (because of the SFR–M_dust relation). On the other hand, high mass bins are dominated by high SFR and therefore associated with high M_dust. This results into an apparent M_star–M_dust correlation. To better visualize this effect in studies that combine together all galaxies (i.e., without binning in a grid of SFR and M_dust), in Fig. 3 we have marked with black circles the bins closest to the MS (and in every case within 0.3 dex from it). These are the bins where the bulk of the star-forming galaxy population is concentrated, and, as expected, they show a steeper M_star–M_dust trend compared to bins of constant SFR.

The dashed line in Fig. 3 represents the expected maximum amount of metals (${\mathit{M}}_{\mathrm{Z}}\mathrm{=}{\mathit{y}}_{\mathrm{Z}}\mathrm{\times }{\mathit{M}}_{\mathrm{star}}^{\mathrm{tot}}$, where y_Z ≃ 0.014 and ${\mathit{M}}_{\mathrm{star}}^{\mathrm{tot}}$ is the total stellar mass formed, including the final products of stellar evolution¹²) produced by stars and supernovae explosions, associated with the star formation required to account for the observed M_star. This is also the maximum amount of dust that can be associated with a given M_star in a “closed box” scenario and assuming a condensation efficiency in the ejecta close to 100%. More realistically, of these metals only about 30–50% (Draine et al. 2007, grey dotted line in Fig. 3) are expected to be depleted into dust grains. These lines give the maximum amount of dust expected as a function of stellar mass if the galaxy behaves as a “closed box”, and metals are condensed in dust grains with reasonable/high efficiency. Most of the galaxies, in particular the high mass systems, lie below the “closed box” lines. This finding qualitatively agrees with the expectations of theoretical models for the evolution of the dust content: rather flat M_dust–M_star trends, i.e., decreasing dust-to-stellar mass ratios as the gas is consumed and transformed into stars (see, e.g., Eales & Edmunds 1996; Calura et al. 2008; Dunne et al. 2011). Alternatively, this result might indicate that most of the dust in these systems is lost. In support of this scenario, independently of the dust information, it has been acknowledged that massive galaxies have a deficit of metals, by a factor of a few, relative to what must have been produced in the same galaxies (Zahid et al. 2012), which is ascribed to winds that have expelled metal-rich gas out of these massive galaxies. On the contrary, hints can be seen for low M_star galaxies (log M_star [M_⊙] ≲ 9.75) to show a high dust mass, close to the maximum “closed box” limit. Recent studies based on SPIRE data in the local and low-z (z < 0.5) Universe support this evidence: large dust-to-stellar mass ratios were reported by Smith et al. (2012), while anti-correlations between the dust-to-stellar mass ratio and stellar mass were observed by Cortese et al. (2012) and Bourne et al. (2012). Due to the necessity of a careful check of optical counterpart associations to IR galaxies with low M_star, we do not extend this work to such low stellar masses. The dust content in low M_star galaxies will be investigated by means of a dedicated analysis in a forthcoming paper.

4.3. Summary view

To give a global view of these correlations, we show in Fig. 4 the SFR–M_star plane at different redshifts, where each bin is colour coded according to the associated dust mass. We also show MS relations from the literature (from Peng et al. 2010 at z ~ 0; and from Santini et al. 2009 at high-z). This representation gives a quick overview on the scaling relations existing between M_star, SFR and M_dust: a weak and sometimes absent trend of M_dust with M_star and a clear correlation between M_dust and SFR.

It is also worth noting that we observe no evidence for evolution of M_dust across the different redshift ranges at a given M_starand SFR; the main difference between the various redshift panels in Fig. 4 is simply that they are populated differently. To make this more clear, Fig. 5 shows M_dust as a function of M_star, in bins of SFR, where the colour coding identifies different redshift bins (note that, as a by-product, Fig. 5 provides further evidence of weak/absent dependence of M_dust on M_star at a fixed SFR). At a given M_starand SFR, there is no clear evidence for evolution of M_dust with redshift within uncertainties. We note, however, that we cannot firmly exclude a decrease in M_dust by a factor of 2 from low- to high-z, though this trend is in a few cases reversed. However, observational uncertainties on our data do not allow us to claim any redshift evolution.

It is certainly true that, on average, the overall amount of dust in galaxies at high redshift is higher, as a consequence of the overall higher ISM content in the bulk of high-z galaxies (see Sect. 5.5). As a matter of fact, the normalization of the MS, representing the locus where the bulk of the population of star-forming galaxies lies, does increase with redshift (e.g., Santini et al. 2009; Rodighiero et al. 2010; Karim et al. 2011) and, as a consequence, the dominant galaxy population moves towards larger SFR, hence being characterized by larger dust masses (Fig. 2). However, our results indicate that galaxies with the same properties (same SFR and same M_star) do not show any significant difference in terms of dust content across the cosmic epochs, at least out to z ~ 2.5. In other words, dust mass in galaxies is entirely determined by the SFR and, to a lesser extent, by M_star, and it is independent of redshift within uncertainties. Put simply, different cosmic epochs are populated by galaxies with different typical SFR and M_star values, and hence are characterized by different dust masses.

At fixed SFR, a non evolving M_dust translates into a non evolving dust temperature (T_dust). This does not contradict the results of Magnelli et al. (2014), presenting only a very smooth negative evolution in the normalization of the T_dust-specific SFR (SSFR = SFR/M_star) relation. They also find a stronger positive evolution in the normalization of the relation between T_dust and the distance from the MS. However, as discussed above, the normalization of the MS itself increases with redshift, hence different SFR–M_star combinations are probed at different epochs.

Given the lack of any significant redshift evolution in the dust mass at fixed M_starand SFR, it is meaningful to represent all redshift bins on the same SFR–M_star plane (upper right panel of Fig. 4) to provide an overview of the dust content over a wider range of M_star and SFR. Here the dust mass is computed by averaging the values at different redshifts. This further confirms the trends already outlined (M_dust depends strongly on the SFR and weakly on M_star), over a wider dynamic range.

5. The evolution of the gas content in galaxies

We investigate here the relation between the gas content and the SFR, as well as the evolution of the gas fraction, with the aim of understanding the processes driving the conversion of gas into stars in galaxies throughout the cosmic epochs. We recall that gas masses are inferred from dust mass measurements by assuming that the dust-to-gas ratio scales with the gas metallicity, and by computing the latter by means of the FMR of Mannucci et al. (2010) (see Eqs. (3), (4)). We verified that all the results presented below are almost unchanged if the redshift-dependent mass-metallicity relation of Maiolino et al. (2008) is used instead of the FMR.

5.1. The star formation law

We plot in Fig. 6 the values of SFR as a function of gas mass. The colour code identifies bins of different redshift. For the sake of clarity, the data points at the different redshifts are also plotted on separate panels on the right side. This figure is analogous to Fig. 2, except that M_gas, plotted here instead of M_dust, takes into account the dependence of the gas metallicity with stellar mass and SFR (see Sect. 3.6). This, however, introduces only a minor effect (the gas metallicity changes less than a factor of 2–3, while the dust mass spans 2–3 orders of magnitude). The relation shown in Fig. 6 can be referred to as the integrated S-K law, meaning that gas masses and SFRs are investigated values rather than their surface densities, as in the original S-K law, where the SFR surface density is related to the gas surface density by a power law relation. We fit the data points with the relation $\log \mathit{SFR}\mathrm{=}\mathit{a}\mathrm{(}\log {\mathit{M}}_{\mathrm{gas}}\mathrm{-}\mathrm{10}\mathrm{)}\mathrm{+}\mathit{b}\mathit{.}$(7)A standard χ² fit cannot be performed on our data given the asymmetric error bars. Therefore, all over our work, we apply a maximum likelihood analysis by assuming rescaled log-normal shapes for the probability distribution functions of the variables with the largest uncertainties (log M_gas in this case) and by ignoring the uncertainties on the other variables. By fitting the total sample we obtain $\mathit{a}\mathrm{=}\mathrm{1.5}{\mathrm{0}}_{-0.10}^{\mathrm{+}\mathrm{0.12}}$ and $\mathit{b}\mathrm{=}\mathrm{1.8}{\mathrm{2}}_{-0.20}^{\mathrm{+}\mathrm{0.21}}$, where a bootstrap is performed to compute the parameter 1σ errors. The best-fit relation is represented by the black solid line in the left panel of Fig. 6. However, due to inhomogeneous sampling in SFR at different redshifts, the fit might suffer from biases in case there is an evolution in the slope or normalization of the relation. To investigate such effects, we also separately fit the points in each individual redshift bin (coloured solid lines in the right panels of Fig. 6). The inferred slopes monotonically decrease with redshift from $\mathrm{1.4}{\mathrm{5}}_{-0.41}^{\mathrm{+}\mathrm{0.37}}$ in the local Universe to $\mathrm{0.7}{\mathrm{6}}_{-0.13}^{\mathrm{+}\mathrm{0.11}}$ at z ~ 2, while the normalizations increase from $\mathrm{1.5}{\mathrm{5}}_{-0.47}^{\mathrm{+}\mathrm{0.43}}$ to $\mathrm{2.1}{\mathrm{0}}_{-0.52}^{\mathrm{+}\mathrm{0.48}}$. The best-fit parameters are given in the bottom right corner of each panel of Fig. 6.

By following the theoretical model of Davé et al. (2011, 2012) and the observational results of Tacconi et al. (2013), we also attempt to fit our data points with a relation that has a single redshift-independent slope and normalization slowly evolving with redshift, i.e., yielding a cosmological scaling of the depletion time (=M_gas/SFR): $\log \mathit{SFR}\mathrm{=}\mathit{m}\mathrm{\left(}\log {\mathrm{M}}_{\mathrm{gas}}\mathrm{-}\mathrm{10}\mathrm{\right)}\mathrm{+}\mathit{n}\log \mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}\mathrm{+}\mathit{q}\mathit{.}$(8)The best-fit parameters are $\mathit{m}\mathrm{=}\mathrm{1.0}{\mathrm{1}}_{-0.17}^{\mathrm{+}\mathrm{0.14}}$, $\mathit{n}\mathrm{=}\mathrm{1.4}{\mathrm{0}}_{-0.74}^{\mathrm{+}\mathrm{0.85}}$ and $\mathit{q}\mathrm{=}\mathrm{1.2}{\mathrm{8}}_{-0.17}^{\mathrm{+}\mathrm{0.14}}$. The dashed-triple dotted lines in the right panels of Fig. 6 show the inferred relation at the median redshift in each bin. However, this function provides a worse fit to the data in terms of probability of the solution as computed from the likelihood, with respect to Eq. (7).

In both cases, the evolution of the relation with redshift may be partly caused by mixing different stellar masses, whose contribution strongly depends on the SFR and redshift because of the evolution of the MS relation.

5.1.1. Comparison with previous works

The inferred relations agree, on average, well with those fitted by previous work based on CO measurements for normal star-forming galaxies (Daddi et al. 2010, lower dashed grey line in Fig. 6, see also Eq. (6); and Genzel et al. 2010¹³, lower dotted grey line), although we fit a steeper slope on all data points. The values of the best-fit slopes are independent of the galaxy population (i.e., consistent fits are found if starburst galaxies are removed, see below), and of the recipe adopted for the gas metallicity (i.e., consistent results are obtained if we assume no dependence on the SFR and redshift evolution of the mass-metallicity relation). Anyhow, the broad agreement with previous studies for the majority of galaxies (see below) and the small dispersion (the average absolute residual is ~0.15 dex in terms of log M_gas) shown by our data points are impressive, especially given the completely different and independent approaches used to derive the star formation law. This confirms the reliability of our approach of deriving gas mass estimates from dust mass measurements.

We remind the reader that the dust method is supposed to trace both the molecular and atomic gas (the dust-to-gas conversion factor adopted refers to the total gas mass). Bigiel et al. (2008) have measured steeper slopes for the star formation laws in local galaxies when both the molecular and atomic gas components are considered. This may explain our steeper slopes compared to previous CO-based studies (e.g., Genzel et al. 2010; Tacconi et al. 2013). However, the fair agreement with the Daddi et al. (2010) relation (inferred from CO observations, a proxy for molecular hydrogen) is suggesting that, if the latter is correct, the bulk of the gas in these galaxies is in the molecular phase, which is reasonable given that most of these galaxies are vigorously forming stars and will have high pressure ISM conditions (see also Leroy et al. 2009; Magdis et al. 2012). The steeper slopes found at low redshift may be determined by a larger atomic-to-molecular gas ratio at low than at high-z (see below). Another possibility to explain this is the trend for the MS template of Elbaz et al. (2011) to slightly underpredict the SFR in the absence of Herschel data for bright galaxies at high-z (SFR > 100  M_⊙/yr, see Fig. B.1 and Berta et al. 2013); by moving the data points with the largest SFR towards lower SFR values, this effect might be responsible for the shallower slope measured at high redshift. However, as it can be seen in Fig. B.1, this effect is not larger than 0.1–0.2 dex, and is therefore unlikely to affect our other results (on the other side, the fitted slope of the S-K law may be sensitive to small offsets in the SFR). As a matter of fact, the results presented in this paper are very similar if other IR templates (e.g., Dale & Helou 2002) are used to measure the SFR from 24 μm fluxes or from all Herschel bands. Finally, steep slopes for the global star formation law may be explained by the results of Saintonge et al. (2013), who claim that the gas-to-dust ratio may be 1.7 times larger at z > 2 than observed locally. This, however, would only marginally affect our highest redshift bins, whose mean redshift value is around 2.

Note that the fact that the slope of the global S-K relation, as well as those at z < 0.6, are steeper than unity implies that galaxies with high SFRs have higher star formation efficiency (defined as SFE = SFR/M_gas, equal to the inverse of the depletion time), even if they are regular, MS galaxies. The magenta dash-dotted lines in Fig. 6 trace the loci with SFE = 1 Gyr^-1 (lower line) and SFE = 10 Gyr^-1 (upper line). As a consequence of the super-linear slope of the S-K relation, moderate star-forming galaxies (SFR ~ 1  M_⊙/yr) have a SFE approaching 1 Gyr^-1, while strongly star-forming galaxies (SFR ~ several 100 M_⊙/yr) have a SFE approaching 10 Gyr^-1, implying gas depletion timescale of a few times 100 Myr. However, the SFE is more properly defined as the ratio of the SFR over the molecular gas content. Therefore, another possibility to interpret our result is that the SFR/M_molgas stays the same, and the atomic gas content decreases in strongly star-forming galaxies, or, in other words, the latter have a larger molecular to atomic fraction. This would be confirmed by the results of Bauermeister et al. (2010), who observe little evolution in the cosmic H i density, while the molecular component is expected to positively evolve out to the peak of cosmic star formation (z ~ 2–3, Obreschkow & Rawlings 2009; Lagos et al. 2011; Popping et al. 2013).

5.1.2. The star formation law for starburst galaxies

Symbols marked with a black box in Fig. 6 correspond to bins which lie in the starburst region of the SFR vs. M_star diagram according to the selection of Rodighiero et al. (2011). They select starburst galaxies as sources deviating from a Gaussian logarithmic distribution of the SSFR, having SSFR four times higher than the peak of the distribution (associated to MS galaxies). Given the average scatter of 0.3 dex of the MS (Noeske et al. 2007), these galaxies are located >2σ above the MS (see Fig. 4). The effectiveness of this SSFR criterion in selecting starburst galaxies is confirmed by semi-analytical models where starburst events are triggered by galaxy interactions during their merging histories (Lamastra et al. 2013a). Galaxies from our sample located in the starburst regions do seem to follow the same star formation law as all other galaxies. We note that the selection of starburst galaxies above is based on the knowledge of the MS from the literature, rather than computed directly on the present sample. However, this does not affect our conclusions. Indeed, the observed correlation between SFR and M_gas is tight enough that, even in case of small variations in the location of the MS, sources selected as starburst would still follow the same relation (for example, results are unchanged if the MS from Whitaker et al. 2012 is used, despite its shallower slope). Unless indicative of a larger fraction of atomic gas in starbursts, this result is in contrast with what suggested by previous studies, mostly based on CO emission (e.g., Daddi et al. 2010; Genzel et al. 2010; Saintonge et al. 2012; Magdis et al. 2012; Sargent et al. 2013). The latter studies find a normalization of the star formation law ~10 times higher for starburst galaxies, implying a larger SFE. In any case, since the slope of the relation that we infer is larger than unity (except at z > 0.6), our result does not imply a low efficiency in converting gas into stars for galaxies located in the starburst region (see next section): starburst galaxies do have, on average, larger star formation efficiency (i.e., shorter depletion times) than the bulk of star-forming galaxies (typically at lower SFR).

We note that our work does not sample the most extreme objects lying at the bright tail of the SFR distribution (all but one of the bins selected as “starbursts” are located between 2σ and 3σ above the MS). Physical properties of very extreme sources, such as local ULIRGs or high-z SMGs, are not always compliant with local-based expectations (see, e.g., Santini et al. 2010) and need to be treated with ad-hoc techniques (for example, Magdis et al. 2012 claim the need of submm data to reliably estimate dust masses of SMGs). Moreover, larger statistics is needed. We will therefore study such extreme sources in a future work.

5.2. The evolution of the star formation efficiency

The slope of the integrated S-K relation inferred from our data is generally steeper than unity (except possibly at high redshift). As a consequence, the SFE for high redshift galaxies, which are also on average more star-forming, is higher than for local galaxies, or, equivalently, the depletion time is shorter (we here assume negligible atomic fraction for all galaxies, but see comment in Sect. 5.1.1). This is illustrated in Fig. 7, where the SFE is plotted as a function of redshift, and where an increase in the SFE with redshift is indeed observed, although with large scatter. Due to degeneracy between SFR evolution and redshift it is not clear whether the increase in the SFE with redshift truly reflects a cosmic evolution of the SFE, i.e., galaxies of a given SFR convert their gas into stars more efficiently at high-z, or it is simply a by-product of the slope of the S-K relation convolved with the higher SFR characterizing high-z galaxies (higher normalization of the MS). In Fig. 7 galaxies with different SFRs are plotted with different colours, in an attempt to break the degeneracy between redshift and SFR. Galaxies with similar SFR show no clear internal evolution with redshift. However, due to observational biases (difficulties in observing faint sources at high-z as well as paucity of rare bright sources in small volumes at low-z) the redshift spanned by each of these sets of points is very narrow, and the dispersion very high, hence we cannot rule out a real, intrinsic evolution of the SFE in galaxies (at a given SFR). We will investigate this issue in more detail in another paper (Santini et al., in prep.).

In any case, regardless of whether the evolution of the SFE is an intrinsic redshift evolution or driven by the slope of the S-K relation and the evolution of the SFR, the net result is that the bulk of the galaxy population (i.e., galaxies on the MS) at high redshift (z ~ 2) do form stars with a SFE higher by a factor of ~5 than the bulk of the population of local star-forming galaxies. This evolution is roughly consistent with the evolution of the dust mass-weighted luminosity (L_IR /M_dust, proportional to the SFE except for a metallicity correction) found by Magdis et al. (2012) (a factor of ~4 from z ~ 0 to z ~ 2) and only slightly steeper than the evolution of the depletion time (a factor of ~3 in the same redshift range) observed by Tacconi et al. (2013), likely due to the steeper S-K law inferred by us compared to their work.

5.3. The evolution of the gas fraction

Figure 8 shows the gas fraction as a function of the stellar mass, colour coded according to the redshift, in panels of different SFR. The gas fraction decreases with the stellar mass, as expected by the gas conversion into stars in a closed-box model, and increases with the SFR, as a consequence of the S-K relation (see also the results of Magdis et al. 2012; and those of the PHIBSS survey presented in Tacconi et al. 2013). Most interesting is the lack of evolution of the gas fraction with redshift, once galaxies are separated according to their M_star and SFR values. Given the assumptions made to compute the gas mass, hence gas fractions, this finding is the result of the lack of (or marginal) evolution of the dust content in bins of fixed M_starand SFR (see Fig. 5), combined with a minor contribution from the gas metallicity evolution with M_star and SFR (the FMR, Mannucci et al. 2010).

From the lack of redshift evolution of the gas fraction at fixed SFR and M_star, it follows that galaxies within a given population (identified by a combinations of SFR and M_star), convert gas at the same rate regardless of redshift, i.e., the physics of galaxy formation is independent of redshift, at least out to the epochs probed by our work. This is essentially a consequence of the unimodal inferred S-K relation, but Fig. 8 shows the result more neatly by also slicing the relation through the dependence on stellar mass, which is the third fundamental parameter. We note that this does not contradict the evolution of the SFE observed in Fig. 7, where different stellar masses and SFR are mixed together and where selection effects cause the different SFR bins to be populated differently at different redshifts (hence the average at each redshift is certainly biased).

In summary, our result implies that, at fixed stellar mass, the SFR is uniquely driven by the gas fraction via the star formation law. In other words, if two among SFR, M_star and M_gas are known, the third property is completely determined and does not depend on redshift. This provides a powerful tool to overcome the observational difficulties related with the measurement of gas or dust masses and analyse the gas content for much larger samples of galaxies.

5.4. The fundamental f_gas–M_star–SFR relation

Given the lack of evolution with redshift observed for the gas fraction once galaxies with the same M_starand SFR are considered, we can combine all redshift bins together to increase the statistics and infer more clearly the trend of f_gas as a function of M_star in different SFR intervals. Figure 9 shows the resulting global dependence of the gas fraction (given by the colour coding) on the SFR–M_star plane. In each SFR interval, we fit the data points with a linear relation in the logarithmic space: $\log {\mathit{f}}_{\mathrm{gas}}\mathrm{=}\mathit{\alpha }\mathrm{+}\mathit{\beta }\mathrm{(}\log {\mathit{M}}_{\mathrm{star}}\mathrm{-}\mathrm{11}\mathrm{)}\mathit{.}$(9)We shift the stellar masses, placing them across zero, in order to de-correlate the slope and offset parameters in the linear fit result. The best-fit parameters are given in Table 1, and the best-fit curves are shown by the dashed grey lines in Fig. 8 and also by the solid coloured lines in Fig. 10, which provides a direct comparison at different SFRs. We note that the functional form adopted above does not necessarily have physical meaning: it is a purely phenomenological representation of the data to better visualize the observed trends and to interpolate the three physical quantities for later use of this 3D relation.

The three-dimensional f_gas–M_star–SFR relation shown in Fig. 10 is a fundamental relation that does not evolve with redshift, at least out to z ~ 2.5. Galaxies move over this surface during their evolution.

Figure 11 shows a 3D representation of such a relation. Further investigation of this 3D relation and its physical interpretation goes beyond the scope of this paper, and will be discussed in a future work, as well as the relation between the independent quantities M_gas, M_star and SFR. Here we only emphasize that the redshift evolution of the S-K law seems to disappear once sources are divided in bins of M_star. Indeed, the redshift evolution of the SFE illustrated in Fig. 7 is most likely a consequence of the fact that high-z bins are mostly populated by galaxies with high SFR, which are characterized by high SFE, as a consequence of the super-linear slope of the S-K relation.

We note that the fundamental relation presented here is indeed a physical result, rather than just a way of looking at the redshift evolution through the evolution of another parameter (e.g., SFR). In other words, the inclusion of the SFR or stellar mass as parameters is not masking a true underlying redshift evolution. As a matter of fact, no similar relation is obtained if redshift is replaced to either SFR or M_star.

5.5. The evolution of the gas fraction among main sequence galaxies

The finding that the fundamental f_gas–M_star–SFR relation does not evolve with redshift does not contradict the claimed redshift evolution of the gas fraction in galaxies (e.g., Daddi et al. 2010; Tacconi et al. 2010, 2013; Magdis et al. 2012). Indeed, as already mentioned, galaxies do not uniformly populate this 3D surface. As they evolve, the bulk of star-forming galaxies populate different regions of this surface, as a consequence of gas accretion, gas consumption by star formation and gas ejection. The projection of such a distribution onto the M_star–SFR plane yields the MS and its evolution with redshift.

As suggested by various models, the evolution of galaxies is likely driven by the evolution of their gas content. The evolution of the MS is likely a by-product of the gas content through the S-K relation, or more generally through the fundamental f_gas–M_star–SFR relation illustrated above. While the evolution of the MS has been constrained by several observations, its driving process, which is the evolution of the gas content, is still loosely constrained. We can however exploit the observed evolution of the MS to infer the evolution of the gas fraction of the population of galaxies dominating star formation at any epoch, by exploiting the fundamental f_gas–M_star–SFR relation.

We take advantage of the mathematical representation of the gas fraction as a function of M_star at given SFR shown in Fig. 10, and we linearly interpolate these relations onto a finer SFR grid. We then adopt the MS relations reported in Fig. 4 and linearly interpolate them onto a fine redshift grid. At a given M_star and redshift, we use the MS relation to compute the expected SFR, according to which we select the appropriate f_gas parameterization.

The resulting evolution of f_gas with stellar mass at different redshifts (colour coded) is shown in the left panel of Fig. 12. The orthogonal plot, i.e., the redshift evolution of f_gas for different stellar masses (colour coded), is shown in the right panel of Fig. 12. These plots illustrate how the “bulk” of the star-forming galaxy population at various epochs populates the 3D fundamental f_gas–M_star–SFR relation as a function of redshift. Essentially, for a given stellar mass, the average gas content of star-forming galaxies increases steadily with redshift, at least out to z ~ 2.5. The increase rate is steeper for low mass galaxies with respect to massive galaxies. Galaxies with log M_star [M_⊙] ~ 10.6 reach f_gas ~ 0.25 around the peak of cosmic star formation at z ~ 2.5, while massive galaxies, with log M_star [M_⊙] ~ 12 reach a gas fraction of only 0.18 at the same cosmic epoch. This behaviour is consistent with a downsizing scenario (Cowie et al. 1996; Fontanot et al. 2009), where massive galaxies have already consumed most of their gas at high redshift, while less massive galaxies have a larger fraction of gas (though more complex scenarios resulting from the interplay of inflows, outflows and star formation are not excluded). Further, in massive galaxies the gas fraction decreases more steeply moving towards lower redshift (with respect to low mass galaxies) and their gas evolution flattens to low values at z ≲ 1.3. Instead, low mass galaxies show a shallower and more regular decrease of the gas content, moving towards lower redshifts. Both trends are further indications of downsizing.

The f_gas values are somewhat lower by a factor of ~1.5–2 on average (after accounting for the IMF conversion) than inferred by the high-z CO survey of Tacconi et al. (2013). A similar or even larger mismatch with CO-based results was found by Conselice et al. (2013), who compute gas fractions from SFR and galaxy sizes by inverting the S-K law. We ascribe the discrepancy to the combination of the various uncertainties associated with CO studies and with our method. In addition, the underestimate by ≃50% of the dust mass of unresolved sources found by Galliano et al. (2011) may also explain the lower values found by us. The gas fractions derived by us are also lower by a factor of ~2 than those published by Magdis et al. (2012), who adopt a similar method. This might be caused by cosmic variance effects: based on the two GOODS fields only, the analysis of Magdis et al. (2012) may be affected by statistical uncertainty. The inclusion of COSMOS data provides much improved statistics that is crucial in stacking analyses. Indeed, the stacking result is closely related to the number of stacked sources. Even if COSMOS is shallower than the deep GOODS fields, SPIRE observations, on which dust masses mostly rely, are confusion limited. Therefore, the statistics is strongly dominated by COSMOS. To verify whether cosmic variance effects could be responsible for such disagreement, we repeated our analysis by only including the two GOODS fields. Given the limited statistics, we end up with only 10 data points. We compared these with our gas fractions and found that in 30% of the cases the former are indeed larger by a factor of 2–2.5, while the rest of the points are consistent within their error bars. Finally, we note that the disagreement with previous works is reduced when the GRASIL model is adopted instead of Draine & Li (2007).

5.6. Comparison with theoretical predictions

The evolution of the gas fraction is a powerful observable to test the various physical processes at play in galaxies and implemented by theoretical models, such as star formation, gas cooling and feedback. Here we compare our findings for the evolution of the gas fraction with the expectations of the semi-analytical model of galaxy formation developed by Menci et al. (2008, and references therein). This connects, within a cosmological framework, the baryonic processes (gas cooling, star formation, supernova feedback) to the merging histories of the dark matter haloes, computed by means of a Monte Carlo simulation. Gas is converted into stars through two main channels: a steady (or quiescent) accretion mode, in which the cold gas in the galaxy disk is converted into stars on long timescales (~1 Gyr), and an interaction-driven mode, where gas destabilized during major and minor mergers and fly-by events is converted into stars on shorter timescales (~10⁷ yr; see Lamastra et al. 2013a,b, for a more detailed description). AGN activity triggered by the same galaxy interactions and the related feedback processes are also included.

The predicted gas fraction as a function of stellar mass and redshift is shown in Fig. 13. On the same figure we report the extrapolations for MS galaxies based on our observations already shown in Fig. 12. As discussed above, MS galaxies represent the bulk of the galaxy population and can be directly compared to the darkest contours, enclosing the region occupied by most of the galaxies.

Observations are generally well reproduced by the theoretical model, although with some systematic deviations. The trends with both stellar mass and redshift are recovered, as well as the downsizing expectations: a strong evolution can be noticed in low mass galaxies (M_star ≲ 10¹¹  M_⊙/yr), which are gas-rich out to z ~ 1 (bottom right panel), while progressively more massive galaxies have already consumed their gas at this epoch (upper right panel). While a very good agreement is recovered for all stellar masses at high redshift (z ~ 2, upper-left panel), the predicted evolution of the gas fraction is more regular than observed at intermediate redshifts, with a gas fraction in log M_star [M_⊙] ≲ 11.5 galaxies of ~0.2 at z ≳ 0.6, around twice the observed value (central left panels). The overall systematic gas richness of model galaxies compared to the observations relates to the long-standing problem of theoretical models in reproducing the galaxy stellar mass functions at high redshift. Indeed, the number of massive galaxies is underpreticted by the models (e.g. Fontanot et al. 2009; Santini et al. 2012a), consistently with the inefficiency of the gas conversion and mass buildup processes in the distant Universe. Once gas consumption has started, it is not effectively suppressed at late stages. Indeed, the model predicts a fraction of very massive (log M_star [M_⊙] ≳ 11.5) galaxies which are still gas-rich at z < 1, at variance with what observed (lower- and central-left panels and top-right one). Although it can be partly ascribed to fluctuations in the f_gas distribution generated by the low number statistics of such high M_star galaxies, this behaviour is a manifestation of a known problem common to all theoretical models, in which the suppression of the star formation activity is still inefficient, despite the feedback processes at work. This is related to the difficulties in reproducing the fraction of red passive galaxies (Fontana et al. 2009).

For all these reasons, the comparison of observed and predicted gas fraction is of major importance to constrain the physical processes implemented in models of galaxy formation and evolution. A more detailed and complete comparison with theoretical expectations will be tackled in a future work.

6. Summary

We have used Herschel data from both PACS and SPIRE imaging cameras to estimate the dust mass of a large sample of galaxies extracted from the GOODS-S, GOODS-N and COSMOS fields. To explore a wide range of galaxy properties, including low mass and moderate star-forming galaxies, we have performed a stacking analysis on a grid of redshifts, stellar masses and SFR, and considered average values. With these outputs we have studied the scaling relations in place between the dust content of galaxies and their stellar mass and SFR at different redshifts, from the local Universe out to z = 2.5. Our main results are the following.

• No clear evolution of the dust mass with redshift is observed at a given SFR and stellar mass. Although there is a global redshift evolution of the dust content in galaxies, as a consequence of the increased ISM content at high-z, our findings indicate that galaxies with the same properties (same SFR and same M_star) do not show any significant difference in terms of dust content across the cosmic epochs, at least out to z ~ 2.5. In other words dust mass in galaxies is mostly determined by SFR and M_star and is independent of redshift.

• The dust content is tightly correlated with the star formation activity of the galaxy. This correlation is in place at all values of M_star probed and at least out to z ~ 2.5. Under the assumption that the dust content is proportional to the gas content (with a factor scaling with the gas metallicity), the observed correlation is a natural consequence of the Schmidt-Kennicutt (S-K) law.

• The correlation between the dust and stellar mass observed by previous studies (which averaged together all SFR) becomes much flatter or even disappears when taken at a fixed SFR. The M_dust − M_star relation is at least partly a result of the M_dust–SFR correlation combined with the main sequence (MS) of star-forming galaxies.

We have then taken one step further and computed gas metallicities from the stellar mass and the SFR according to the fundamental metallicity relation (FMR) fitted by Mannucci et al. (2010), and estimated gas masses by assuming that the dust-to-gas ratio linearly scales with the gas metallicity. We note that all our results are robust against the specific parameterization chosen to describe the gas metallicity (e.g., FMR against redshift-dependent mass-metallicity relation). This method provides a complementary approach to investigate the galaxy gas content independently of CO observations. Under our assumptions we find the following.

• We fit a power law relation between the SFR and the gas mass, in good agreement with that previously obtained by Daddi et al. (2010), and also broadly consistent with the results of Genzel et al. (2010). This agreement is remarkable, given the completely different approach between our study and the two works above based on CO measurements. We find that all galaxies follow the same star formation law (integrated S-K law), with no evidence of starbursts lying on an offset relation, though our sample lacks the most extremely starbursting sources (such as local ULIRGs and their analogues at high-z). The slope of this relation is on average steeper than unity, implying that strongly star-forming galaxies have higher star formation efficiency (SFE, i.e., the inverse of the depletion time), or shorter depletion time. We also find a mild, but significant evolution of the S-K law with redshift.

• We observe an evolution of the SFE with redshift, by about a factor of 5 from z ~ 0 to z ~ 2.5. This applies to the bulk of the galaxy population dominating star formation at each epoch. However, it is not clear whether such evolution is an intrinsic redshift evolution or is simply a consequence of sampling more star-forming galaxies at high redshift combined with the slope of the integrated S-K relation being on average steeper than unity.

• The measured gas fraction decreases with stellar mass and increases with SFR, as expected. However, when considering bins of given stellar mass and SFR, the gas fraction does not show any redshift evolution, at least out to z ~ 2.5. This primarily results from the non-evolution of the dust mass (within uncertainties), with gas metallicity effects only providing a second-order contribution. The 3D relation between f_gas, M_star and SFR is a fundamental relation that holds at any redshift. It provides a powerful tool to overcome the observational difficulties related with the measurement of gas or dust masses and to analyse the gas content for much larger samples of galaxies. Galaxies populate such a 3D fundamental f_gas–M_star–SFR relation in a different way throughout the cosmic epochs. The projection of galaxies on the 3D fundamental relation onto the M_star–SFR plane gives the MS and its evolution with redshift.

• We “de-project” the MS galaxies, at various cosmic epochs, onto the 3D fundamental f_gas–M_star–SFR relation, to infer the evolution of the gas fraction of “typical” star-forming galaxies as a function of redshift. A clear redshift evolution from z ~ 0 to z ~ 2.5 in the gas fraction is observed for MS galaxies. The evolution of the gas content in massive (M_star ≳ 10¹¹  M_⊙) galaxies is steep between z ~ 2.5 and z ~ 1.3 and flattens to low f_gas values at lower redshifts. Low mass M_star ≲ 10¹¹  M_⊙ galaxies show a less steep and more regular decrease of the gas fraction from z ~ 2.5 to z ~ 0. These trends are in agreement with the downsizing scenario for galaxy evolution, and they are on average well reproduced by the theoretical expectations of the semi-analytical model of Menci et al. (2008), despite a systematic larger gas richness compared to our data.

Online material

Appendix A: Statistics on the z–M_star–SFR grid

We report in Tables A.1 the number of sources in each z–M_star–SFR bin and the associated average dust mass.

Appendix B: Reliability of the SFR estimates

To verify the reliability of the 24 μm-based SFR tracer (SFR₂₄), we compare it with the SFR measured by fitting the average FIR stacked flux densities (SFR_FIR). We fit these flux densities with the Dale & Helou (2002) IR template library to infer an estimate of the total IR luminosity, and account for the unobscured SFR by taking into account the average rest-frame UV luminosity in the bin uncorrected for extinction (see Sect. 3.2 for more details). The same results are found if the MS template of Elbaz et al. (2011) is used instead of the Dale & Helou (2002) library. In the left panel of Fig. B.1 we compare the average SFR₂₄ in each bin of the grid with the average SFR_FIR. The two SFR measurements nicely agree with each other with the only noticeable exception of the lowest redshift bin, making us confident of the method adopted. The small offset observed at low-z does not significantly affect our results.

We also repeated the same test by making use of a “ladder of SFR indicators” (SFR_ladder, e.g. Wuyts et al. 2011a) as input for the grid production. Such a “ladder” approach combines different SFR estimates by using the best available choice for each galaxy. More specifically, a Herschel-based SFR is used for galaxies detected by PACS or SPIRE, the 24 μm-based tracer

is adopted for galaxies undetected by Herschel but detected by MIPS, and the output of the optical-UV SED fitting described in Sect. 3.1 is used for galaxies undetected at IR wavelengths. Most importantly, this approach has the advantage of increasing the number of galaxies for which a SFR estimate is available and enlarging the SFR dynamical range. However, as evident from the right panel of Fig. B.1, the scatter with respect to SFR_FIR is larger than in the previous case. Moreover, the correlation between SFR_ladder and SFR_FIR flattens at low SFR and z > 0.2, exactly at the SFR regime where in principle the “ladder” approach provides an improvement over the 24 μm-based SFR. One possibility to explain the flat behaviour at low SFR (below a redshift-dependent threshold) shown in the right panel of Fig. B.1 is to ascribe it to failures in the associations of optical counterparts for the extremely faint IR galaxies or blending issues mostly affecting the faintest galaxies during the stacking procedure. Alternatively, dust heating by old stellar population might also be responsible for the enhanced IR flux at low SFR. However, investigating the reasons of such disagreement is beyond the scope of the present work. Based on the tests performed, we decide to use SFR₂₄ as a SFR tracer, at the expenses of reducing the SFR dynamical range, in order not to run the risk to introduce systematics in the analysis.

Appendix C: Simulation to test against possible degeneracies in the SFR–M_dust correlation

We run a simulation to verify that the trend observed between the SFR and the dust mass (Fig. 2) is real and not a trivial outcome of the fact that both physical variables are related to the FIR peak of the galaxy SED. Indeed, while the SFR is simply proportional to the integrated light in the dust emission peak, the dust mass depends not only on the normalization of the spectrum but also on the temperature of the grains, which determines its shape.

The aim of the simulation is to verify whether an initially scattered and uncorrelated distribution of mock SFR and M_dust gives rise to a correlation when the two quantities are recomputed according to our methods. To do that, we consider a set of GRASIL templates, each associated to a dust mass (M_{dust  INPUT}) and to a SFR (SFR_INPUT). The latter is computed by integrating the template from 8 and 1000 μm and multiplying by the calibration factor 1.8 × 10^-10 (see Sect. 3.2). In order to sample a wide region of the SFR–M_dust parameter space, we multiply each SFR–M_dust pair by a normalization factor. This corresponds to multiplying the associated SED, since both parameters scale with the SED normalization. We consider a range in SFR and M_dust which mimic that observed in one of the redshift interval most populated by our data, i.e., 0.6 < z < 1. The resulting SFR–M_dust distribution is shown in the left panel of Fig. C.1.

We redshift each mock galaxy to a random redshift within the 0.6-1 interval, and interpolate the associated spectrum with the MIPS 24 μm and Herschel 100–500 μm filters. To mimic the real case, we perturb such mock flux densities by adding a noise randomly extracted from the observed noise distribution in each band. We then measure the SFR (SFR_OUTPUT) and M_dust (M_{dust  OUTPUT}) for each mock galaxy exactly as we have done for the real data. Consistently with what described in Sect. 3.5, we reject sources not compliant with our requirements to ensure reliable dust mass estimates. The resulting measurements are shown in the right panel of Fig. C.1 and show no evidence for any correlation between SFR and M_dust. The absence of any correlation induced by our measures is statistically confirmed by the values of the Pearson coefficients on the input (0.11) and output (0.19) data point distributions. This simulation illustrates that the correlations between M_dust and SFR is not an artefact of the method, but is physical (i.e., the result of the S-K law).

Appendix D: Fits of the FIR SEDs

In Fig. D.1 we report the best fits of Herschel stacked flux densities with Draine & Li (2007) templates computed to estimate the dust mass. The secondary bump around 50 μm which can be seen in few of the best-fit SEDs is due to a warm dust component. This feature gradually disappears when the maximum radiation intensity (U_max) in the Draine & Li (2007) model is set to lower values. However, fixing U_max to a value lower than 10⁶ has the overall effect of making each template slightly warmer and lowering the normalization of the Rayleigh-Jeans side of the spectrum. This has the effect of increasing the inferred dust masses by a factor of 1.5–2, due to a larger normalization of the SED for a given set of observed flux densities. We decide to follow the prescription given by Draine et al. (2007) and fixing U_max to 10⁶. However, we note that an offset would not change our main results.

---

Acknowledgments

P.S. thanks N. Scoville for interesting and useful discussions and A. Marconi and G. Risaliti for helping with fitting routines and statistical issues. This work was supported by grant ASI I/005/11/0. PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/ OAA/OAP/OAT, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain). SPIRE has been developed by a consortium of institutes led by Cardiff University (UK) and including University of Lethbridge (Canada), NAOC (China), CEA, LAM (France), IFSI, University of Padua (Italy), IAC (Spain), Stockholm Observatory (Sweden), Imperial College London, RAL, UCL-MSSL, UKATC, University of Sussex (UK), Caltech, JPL, NHSC, University of Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA).

References

All Tables

All Figures