Free Access
Issue
A&A
Volume 573, January 2015
Article Number A113
Number of page(s) 17
Section Cosmology (including clusters of galaxies)
DOI https://doi.org/10.1051/0004-6361/201425031
Published online 07 January 2015

© ESO, 2015

1. Introduction

Galaxy properties evolve rapidly across cosmic time. In particular, various studies have shown that the mean star formation rate (SFR) at fixed stellar mass increases by a factor of about 20 between z = 0 and z = 2 (e.g., Noeske et al. 2007; Elbaz et al. 2007, 2011; Daddi et al. 2007; Pannella et al. 2009, 2014; Magdis et al. 2010; Karim et al. 2011; Rodighiero et al. 2011; Whitaker et al. 2012; Heinis et al. 2014). This very high SFR can be explained by either larger reservoirs of molecular gas or a higher star formation efficiency (SFE). Large gas reservoirs have been found in massive galaxies at high redshift (e.g., Daddi et al. 2008, 2010a; Tacconi et al. 2010, 2013; Aravena et al. 2013), which could imply high SFRs with SFE similar to that of normal star-forming galaxies in the local Universe. On the other hand, follow-up of bright submillimeter galaxies (SMGs) revealed that their very intense SFR (~1000 M/yr) is also driven by a SFE boosted by a factor of 10 with respect to normal star-forming galaxies in the local Universe (e.g., Greve et al. 2005; Frayer et al. 2008; Daddi et al. 2009a,b), likely induced by a major merger. This difference can be understood if we consider that galaxies are driven by two types of star formation activity: smooth processes fed by large reservoirs of gas in normal star-forming galaxies and boosted star-formation in gas rich mergers (Daddi et al. 2010b; Genzel et al. 2010).

Using models based on the existence of this main-sequence of star-forming galaxies, i.e., a tight correlation between SFR and stellar mass, and outliers of this sequence with boosted sSFRs (SFR/M) called starbursts hereafter, Sargent et al. (2012) showed that the galaxies with the highest SFR mainly correspond to starbursts, while the bulk of the star formation budget (~85%) is hosted in normal star-forming galaxies. This approach allows us to better understand the heterogeneous characteristic of distant objects concerning their gas fraction and their SFE (Sargent et al. 2014). The quick rise of the sSFR would thus be explained by larger gas reservoirs in main-sequence galaxies. However, the most extreme SFRs observed in high-redshift starbursts would be caused by a SFE boosted induced by major mergers.

At high redshift, the gas mass is difficult to estimate. Two main methods are used. The first is based on the measurement of the intensity of rotational transitions (generally with Jupper< 3) of 12CO and an assumed CO-to-H2 conversion factor (Daddi et al. 2008; Tacconi et al. 2010, 2013; Saintonge et al. 2013). The main limitation of this method is the uncertainty on this conversion factor, which is expected to be different from the local calibrations in high-redshift galaxies with strongly sub-solar metallicities (Bothwell et al. 2010; Engel et al. 2010; Genzel et al. 2012, 2014; Tan et al. 2013). The second method is based on the estimate of the dust mass, which is then converted into gas mass using the locally-calibrated relation between the gas-to-dust ratio and the gas metallicity (e.g., Muñoz-Mateos et al. 2009; Leroy et al. 2011; Rémy-Ruyer et al. 2014). The main weakness of this method is the need of an accurate estimate of the gas metallicity and the possible evolution in normalization and scatter of the relation between gas-to-dust ratio and gas metallicity. This method was applied on individual galaxies at high redshift by Magdis et al. (2011, 2012a) and Scoville et al. (2014), but also on mean spectral energy distributions (SEDs) measured through a stacking analysis (Magdis et al. 2012a; Santini et al. 2014). This method has not been applied at redshifts higher than ~2. The aim of this paper is to extend the studies of dust emission and gas fractions derived from dust masses to z ~ 4 and analyze possible differences in trends as redshift increases.

In this paper, we combine the information provided by the Herschel data and a mass-selected sample of galaxies built from the UltraVISTA data (Ilbert et al. 2013) in COSMOS to study the mean dust emission of galaxies up to z = 4 (Sect. 2). We measure the mean SED of galaxies on the main sequence and strong starbursts using a stacking analysis. We then deduce the mean intensity of the radiation field and the mean dust mass in these objects using the Draine & Li (2007) model (Sect. 3). We discuss the observed evolution of these quantities in Sect. 4 and the consequences on the nature of star formation processes at high redshift in Sect. 5. Throughout this paper, we adopt a ΛCDM cosmology with Ωm = 0.3, ΩΛ = 0.7, H0 = 70 km s-1Mpc-1 and a Chabrier (2003) initial mass function (IMF).

2. Data

2.1. Stellar mass and photometric redshift catalog using UltraVISTA data

Deep Y, J, H, and Ks data (mAB,5σ ~ 25 for the Y band and 24 for the others) were produced by the UltraVISTA survey (McCracken et al. 2012). The photometric redshift and the stellar mass of the detected galaxies were estimated using Le PHARE (Arnouts et al. 1999; Ilbert et al. 2006) as described in Ilbert et al. (2013). The precision of the photometric redshifts at 1.5 <z< 4 is σΔz/(1 + z) = 0.03. According to Ilbert et al. (2013), this catalog is complete down to 1010.26M at z< 4. X-ray detected active galactic nuclei (AGNs) are also removed from our sample of star-forming galaxies, since the mid-infrared emission of these objects could be strongly affected the AGN. Luminous X-ray obscured AGNs might still be present in the sample. However, their possible presence appear to have limited impact on our work as no mid-infrared excess is observed in the average SEDs measured by stacking (see Figs. 4 and 5 and Sect. 4).

As this paper studies star-forming galaxies, we focused only on star-forming galaxies selected following the method of Ilbert et al. (2010) based on the rest-frame NUVr+ versus r+J and similar to the UVJ criterion of (Williams et al. 2009). The flux densities in each rest-frame band are extrapolated from the closest observer-frame band to minimize potential biases induced by the choice of template library. At z> 1.5, 4060% of the objects classified as passive by this color criterion have a sSFR>10-11 yr-1 according to the SED fitting of the optical/near-IR data (Ilbert et al. 2013, their Fig. 3). However, the sSFRs obtained by SED fitting are highly uncertain, because of the degeneracies with the dust attenuation. These peculiar objects are at least 10 times less numerous than the color-selected star-forming sample in all redshift bins. Including them or not in the sample has a negligible impact (~0.25σ) on the mean SEDs measured by stacking (see Sect. 3). We thus based our study only on the color-selected population for simplicity.

2.2. Spitzer/MIPS data

The COSMOS field (2 deg2) was observed by Spitzer at 24 μm with the multiband imaging photometer (MIPS). A map and a catalog combined with the optical and near-IR data was produced from these observations (Le Floc’h et al. 2009). The 1σ point source sensitivity is ~15 μJy and the full width at half maximum (FWHM) of the point spread function (PSF) is ~6′′.

2.3. Herschel/PACS data

The PACS (photodetecting array camera and spectrometer, Poglitsch et al. 2010) evolutionary probe survey (PEP, Lutz et al. 2011) mapped the COSMOS field with the Herschel1 space observatory (Pilbratt et al. 2010) at 100 and 160 μm with a point-source sensitivity of 1.5 mJy and 3.3 mJy and a PSF FWHM of 7.7′′ and 12′′, respectively. Sources and fluxes of the PEP catalog were extracted using the position of 24 μm sources as a prior. This catalog is used only to select strong starbursts up to z ~ 3. The 24 μm prior should not induce any incompleteness of the strong-starburst sample, since their minimum expected 24 μm flux is at least 2 times larger than the detection limit at this wavelength2.

2.4. Herschel/SPIRE data

We also used Herschel data at 250 μm, 350 μm, and 500 μm taken by the spectral and photometric imaging receiver (SPIRE, Griffin et al. 2010) as part of the Herschel multitiered extragalactic survey (HerMES, Oliver et al. 2012). The FWHM of the PSF is 18.2′′, 24.9′′, and 36.3′′, the 1σ instrumental noise is 1.6, 1.3, and 1.9 mJy/beam, and the 1σ confusion noise is 5.8, 6.3, and 6.8 mJy/beam (Nguyen et al. 2010) at 250 μm, 350 μm, and 500 μm, respectively. In this paper, we used the sources catalog extracted using as a prior the positions, the fluxes, the redshifts, and mean colors measured by stacking of 24 μm sources, as described in Béthermin et al. (2012b).

2.5. LABOCA data

The COSMOS field was mapped at 870 μm by the large APEX bolometer Camera (LABOCA) mounted on the Atacama Pathfinder Experiment (APEX) telescope3 (PI: Frank Bertoldi, Navarrete et al., in prep.). We retrieved the raw data from the ESO Science Archive facility and reduced them with the publicly available CRUSH (version 2.12–2) pipeline (Kovács 2006, 2008). We used the algorithm settings optimized for deep field observations4. The output of CRUSH includes an intensity map and a noise map. The mapped area extends over approximately 1.4 square degrees with a non-uniform noise that increases toward the edges of the field. In this work we use the inner ~0.7 deg2 of the map where a fairly uniform sensitivity of ~4.3 mJy/beam is reached (Pannella et al. in prep.) with a smoothed beam size of ~27.6′′. Contrary to SPIRE data, which are confusion limited, LABOCA data are noise limited and the maps are thus beam-smoothed to minimize their RMS.

2.6. AzTEC data

An area of 0.72 deg2 was scanned by the AzTEC bolometer camera mounted on the Atacama submillimeter telescope experiment (ASTE). The sensitivity in the center of the field is 1.23 mJy RMS and the PSF FWHM after beam-smoothing is 34′′ (Aretxaga et al. 2011).

3. Methods

3.1. Sample selection

thumbnail Fig. 1

Stellar mass distribution of our samples of star-forming galaxies in the various redshift bins we used. Only galaxies more massive than our cut of 3 × 1010M are represented. The first bin contain fewer objects than the second one because our cut fall at the middle of the first one. The arrows indicate the mean stellar mass in each redshift bin.

Open with DEXTER

In this paper, we base our analysis on mass-selected samples of star-forming galaxies (see Sect. 2.1). We chose the same stellar mass cut of 3 × 1010M at all redshifts to be complete up to z ~ 4. We could have used a lower mass cut at lower redshifts, but we chose this single cut for all redshifts to be able to interpret the observed evolution of the various physical parameters of the galaxies in our sample in an easier way. This cut is slightly higher than the 90% completeness limit at z ~ 4 cited in Ilbert et al. (2013, 1.8 × 1010M and implies an high completeness of our sample, which limits potential biases induced by the input catalog on the results of our stacking analysis (e.g., Heinis et al. 2013). The exact choice of our stellar mass cut has negligible impact on the mean SEDs measured by stacking: we tested a mass cut of 2 × 1010M and 5 × 1010M and found that, after renormalization at the same LIR, the SEDs are similar ( and 0.79, respectively). These results agree with Magdis et al. (2012a), who did not find any evidence of a dependence of the main-sequence SED on stellar mass at fixed redshift. The mass distribution of star-forming galaxies does not vary significantly with redshift, except in normalization (Ilbert et al. 2013 and Fig. 1). The average stellar mass at all redshifts is between 1010.75M and 1010.80M (Fig. 1 and Table 2).

Star-forming galaxies whose stellar mass is larger than our cut do not correspond to the same populations at z = 4 and z = 0. The massive objects at z = 4 are formed in dense environments, corresponding to the progenitors of today’s clusters and massive groups (e.g., Conroy et al. 2009; Moster et al. 2010; Behroozi et al. 2013; Béthermin et al. 2013, 2014). Most of these objects are in general quenched between z = 4 and z = 0 (e.g., Peng et al. 2010). In contrast, our mass cut at z = 0 corresponds to Milky-Way-like galaxies. At all redshift, this cut is just below the mass corresponding to the maximal efficiency of star formation inside halos (defined here as the ratio between stellar mass and halo mass, Moster et al. 2010, 2013; Behroozi et al. 2010; Béthermin et al. 2012b; Wang et al. 2013).

Our stellar mass cut is slightly below the knee of the mass function of star-forming galaxies (Ilbert et al. 2013). The population we selected thus hosts the majority (>50%) of the stellar mass in star-forming galaxies. Since there is a correlation between stellar mass and SFR, we are thus probing the population responsible for a large fraction the star formation (4065% depending on the redshift according to the Béthermin et al. 2012b model; see also Karim et al. 2011). Our approach is thus different from Santini et al. (2014) who explore in detail how the SEDs evolve at z< 2.5 in the SFR-M plane using a combination of UV-derived and 24 μm-derived SFRs. We aim to push our analysis to higher redshifts and we thus use this more simple and redshift-invariant selection to allow an easier interpretation and to limit potential selection biases. In addition to this mass selection, we divide our sample by intervals of redshift. The choice of their size is a compromise between large intervals to have a good signal-to-noise ratio at each wavelength and small intervals to limit the broadening of the SEDs because of redshift evolution within the broad redshift bin.

We also removed strong starbursts from our sample (sSFR> 10sSFRMS) and studied them separately. These objects are selected using the photometric catalogs described in Sect. 2. For the sources which are detected at 5σ at least in two Herschel bands, we fitted the SEDs with the template library of Magdis et al. (2012a) allowing the mean intensity of the radiation field U to vary by ±0.6 dex (3σ of the scatter used in the Béthermin et al. 2012a model). These criteria of two detections at different wavelengths and the high reliability of the detections prevent biasing of the starbursts towards positive fluctuations of the noise in the maps and limit the flux boosting effect. We then estimated the SFR from the infrared luminosity, LIR, using the Kennicutt (1998) relation. We performed a first analysis using the same evolution of the main-sequence (sSFRMS versus z) as in Béthermin et al. (2012a) to select sSFR> 10sSFRMS objects. We then fit the measured evolution of the main-sequence found by a first stacking analysis (see Sects. 3.2 and 3.3) to prepare the final sample for our analysis. We could have chosen a lower sSFR cut corresponding to 4 times the value at the center of the main-sequence as in Rodighiero et al. (2011), but the sample would be incomplete at z> 1 because of the flux limit of the infrared catalogs.

thumbnail Fig. 2

The thick red solid line represents the luminosity limit corresponding to a criterion of a 5σ detection in at least two Herschel bands. The other solid lines are the limits for a detection at only one given wavelength (purple for 100 μm, blue for 160 μm, turquoise for 250 μm, green for 350 μm, orange for 500 μm). The dashed, dot-dash, and three-dot-dash lines indicate the infrared luminosity of a galaxie of 3 × 1010M (our mass cut) at the center of the main sequence, a factor of 4 above it, and a factor of 10 above it, respectively.

Open with DEXTER

Figure 2 shows the luminosity limit corresponding to a detection at 5σ at two wavelengths or more. This was computed using both the starburst and the main-sequence templates of the Magdis et al. (2012a) SED library. This library contains different templates for main-sequence and starburst galaxies. The main-sequence template evolves with redshift, but not the starburst one. The lines correspond to the highest luminosity limit found using these two templates for each wavelength, which is the most pessimistic case. We also computed the infrared luminosity associated with a galaxy of 3 × 1010M, i.e., our mass limit, on the main sequence (dashed line), a factor of 4 above it (dot-dash line), and a factor of 10 above it (three-dot-dash line). All the M> 3 × 1010M strong starbursts (sSFR> 10sSFRMS) should thus be detected in two or more Herschel bands below z = 4. There is only one starburst detected in the 3 <z< 4 bin. We thus do not analyze this bin, because of its lack of statistical significance. The other bins contain 3, 6, 6, and 8 strong starbursts, respectively, by increasing redshift.

The sample of main-sequence galaxies is contaminated by the starbursts which have sSFR< 10sSFRMS. We expect that this contamination is negligible, since the contribution of all starbursts to the infrared luminosity density is lower than 15% (Rodighiero et al. 2011; Sargent et al. 2012). To check this hypothesis, we statistically corrected for the contribution of the remaining starbursts with sSFR< 10sSFRMS based on the Béthermin et al. (2012b) counts model. We assumed both the SED library used for the model and the average SED of strong starbursts found in this study. We found that this statistical subtraction only affected our measurements at most at the 0.2σ level. Consequently, we have neglected this contamination in the rest of our study.

3.2. Stacking analysis

We use a similar stacking approach as in Magdis et al. (2012a) to measure the mean SEDs of our sub-samples of star-forming galaxies from the mid-infrared to the millimeter domain. Different methods are used at the various wavelength to optimally extract the information depending if the data are confusion or noise limited. At 24 μm, 100 μm, and 160 μm, we produced stacked images using the IAS stacking library (Bavouzet 2008; Béthermin et al. 2010a). The flux is then measured using aperture photometry with the same parameters and aperture corrections as Béthermin et al. (2010a) at 24 μm. At 100 μm and 160 μm, we used a PSF fitting technique. A correction of 10% is applied to take into account the effect of the filtering of the data on the photometric measurements of faint, non-masked sources (Popesso et al. 2012). At 250 μm, 350 μm, and 500 μm, the photometric uncertainties are not dominated by instrumental noise but by the confusion noise caused by neighboring sources (Dole et al. 2003; Nguyen et al. 2010). We thus measured the mean flux of the sources computing the mean flux in the pixels centered on a stacked source following Béthermin et al. (2012b). This method minimizes the uncertainties and a potential contamination caused by the clustering of galaxies (Béthermin et al. 2010b). Finally, we used the same method, but on the beam-convolved map, for LABOCA and AzTEC data as they are noise limited and lower uncertainties are obtained after this beam smoothing. LABOCA and AzTEC maps do not cover the whole area. We thus only stack sources in the covered region to compute the mean flux densities of our various sub-samples. The source selection criteria being exactly the same inside and outside the covered area, this should not introduce any bias.

These stacking methods can be biased if the stacked sources are strongly clustered or very faint. This bias is caused by the greater probability of finding a source close to another one in the stacked sample compared to a random position. This effect has been discussed in detail by several authors (e.g., Bavouzet 2008; Béthermin et al. 2010b, 2012b; Kurczynski & Gawiser 2010; Bourne et al. 2012; Viero et al. 2013). In Magdis et al. (2012a), the authors estimated that this bias is lower than the 1σ statistical uncertainties and was not corrected. The number of sources to stack in COSMOS compared to the GOODS fields used by Magdis et al. (2012a) is much larger and hence the signal-to-noise ratio is much better. The bias caused by clustering is thus non-negligible in COSMOS. Because of the complex edge effects caused by the absence of data around bright stars, the methods using the position of the sources to deblend the contamination caused by the clustering cannot be applied (Kurczynski & Gawiser 2010; Viero et al. 2013). Consequently, we developed a method based on realistic simulations of the Spitzer, Herschel, LABOCA, and AzTEC maps to correct this effect, which induces biases up to 50% at 500 μm around z ~ 2. The technical details and discussion of these corrections are presented in Appendix A.

The uncertainties on the fluxes are measured using a bootstrap technique (Jauzac et al. 2011). This method takes into account both the errors coming from the instrumental noise, the confusion, and the sample variance of the galaxy population (Béthermin et al. 2012b). These uncertainties are combined quadratically with those associated with the calibration and the clustering correction.

3.3. Mean physical properties from SED fitting

We interpreted our measurements of the mean SEDs using the Draine & Li (2007) model as in Magdis et al. (2012a). This model, developed initially to study the interstellar medium in the Milky Way and in nearby galaxies, takes into account the heterogeneity of the intensity of the radiation field. The redshift slices we used have a non-negligible width. To account for this, we convolve the model by the redshift distribution of the galaxies before fitting the data. The majority of the redshifts in our sample are photometric. We thus sum the probability distribution function (PDF) of the redshifts of all the sources in a sub-sample to estimate its intrinsic redshift distribution. The uncertainties on the physical parameters are estimated using the same Monte Carlo method as in Magdis et al. (2012a). The uncertainties on each parameter takes into account the potential degeneracies with the others, i.e., they are the marginalized uncertainties on each individual parameters. Our good sampling of the dust SEDs (8 photometric points between 24 μm and 1.1 mm including at least six detections) allows us to break the degeneracy between the dust temperature and the dust mass which is found if only (sub-)mm datapoints are used.

Instead of using the three parameters describing the distribution of the intensity of the radiation field U of the Draine & Li (2007) model (the minimal radiation field Umin, the maximal one Umax, and the slope of the assumed power-law distribution between these two values α), we considered only the mean intensity of the radiation field U for simplicity. The other parameters derived from the fit and used in this paper are the bolometric infrared luminosity integrated between 8 and 1000 μm (LIR) and the dust mass (Md). The SFR is derived from LIR using the Kennicutt (1998) conversion factor ( after conversion from Salpeter to Chabrier IMF), since the dust-obscured star formation vastly dominates the unobscured component at z< 4 given the mass-scale considered (Heinis et al. 2013, 2014; Pannella et al. 2014). The sSFR is computed using the later SFR and the mean stellar mass extracted from the Ilbert et al. (2013) catalog. The uncertainties on the derived physical parameters presented in the various figures and tables of this paper are the uncertainties on the average values. The dispersion of physical properties inside a population is difficult to measure by stacking and we did not try to compute it in this paper (see Sect. 5).

The residuals of these fits are presented in Appendix B. Tables 1 and 2 summarize the average photometric measurements and the recovered physical parameters, respectively.

Table 1

Summary of our flux density measurements by stacking.

Table 2

Summary of the average physical parameters of our samples.

4. Results

4.1. Evolution of the mean SED of star-forming galaxies

thumbnail Fig. 3

Mean flux density as a function of wavelength (observed wavelength in the top panels and rest-frame wavelength in the bottom panels) at various redshifts (see color coding). The left panels show the mean SEDs of the main-sequence sample and the right panels those of the strong starbursts.

Open with DEXTER

thumbnail Fig. 4

Rest-frame mean spectral energy distribution of our selection of massive, star-forming galaxies at various redshift measured by stacking analysis. The data points are fitted using the Draine & Li (2007) model. This model is convolved with the redshift distribution of the sources before being compared to the data. The black and blue lines represent the intrinsic and convolved SEDs, respectively. The bottom right corner summarizes the redshift evolution seen in our data.

Open with DEXTER

thumbnail Fig. 5

Rest-frame mean spectral energy distribution of our selection of strong starbursts at various redshift measured by stacking analysis. The data points are fitted using the Draine & Li (2007) model. This model is convolved with the redshift distribution of the sources before being compared to the data. The black and red lines show the intrinsic and convolved SEDs, respectively.

Open with DEXTER

Figure 3 summarizes the results of our stacking analysis. For the main-sequence sample, the flux density varies rapidly with redshift in the PACS 100 μm band, while it is almost constant in the SPIRE 500 μm band. The peak of the flux density distribution in the rest frame moves from ~120 μm to 70 μm between z = 0 and z = 4. This shift with redshift was already observed at z ≲ 2 for mass-selected stacked samples (Magdis et al. 2012a) or a Herschel-detected sample (Lee et al. 2013; Symeonidis et al. 2013). We found no evidence of an evolution of the position of this peak (~70 μm) for the sample of strong starbursts.

Figures 4 and 5 show the mean intrinsic luminosity (in νLν units, the peak of the SEDs is thus shifted toward shorter wavelengths compared with Lν units) of our samples of massive star-forming galaxies (since this sample is dominated by main-sequence galaxies, hereafter we call it main-sequence sample) and the fit by the Draine & Li (2007) model. We also observe a strong evolution of the position of the peak of the thermal emission of dust in main-sequence galaxies from ~80 μm at z ~ 0.4 to ~30 μm at z ~ 3.75 in νLν units. The SEDs of strong starbursts have a much more modest evolution (from 50 μm at to 30 μm). The mean luminosity of the galaxies also increases very rapidly with redshift for both main-sequence and strong starburst galaxies.

At z> 2, we find that the peak of the dust emission tends to be broader than at lower redshift. The broadening of the mean SEDs induced by the size of the redshift bins has a major impact only on the mid-infrared, where the polycyclic aromatic hydrocarbon (PAH) features are washed out (see black and blue lines in Figs. 4 and 5), and cannot fully explain why the far-IR peak is broader at higher redshifts. The Draine & Li (2007) model reproduces this broadening by means of a higher γ coefficient, i.e., a stronger contribution of regions with a strong heating of the dust. This is consistent with the presence of giant star-forming clumps in high-redshift galaxies (e.g., Bournaud et al. 2007; Genzel et al. 2006). The best-fit models at high z presents two breaks around 30 μm and 150 μm, which could be artefacts caused by the sharp cuts of the U distribution at its extremal values in the Draine & Li (2007) model.

4.2. Evolution of the specific star formation rate

From the fit of the SEDs, we can easily derive the evolution of the mean specific SFR of our mass-selected sample with redshift. The results are presented in Fig. 6. The strong starbursts lie about a factor of 10 above the main-sequence, demonstrating that this population is dominated by objects just above our cut of 10 sSFRMS. Our results can be fitted by an evolution in redshift as (0.061±0.006 Gyr-1) × (1 + z)2.82 ± 0.12 at z< 2 and as (1 + z)2.2 ± 0.3 at z> 2. We compared our results with the compilation of measurements of Sargent et al. (2014) at M = 5 × 1010M. At z< 1.5, our results agree well with the previous measurements. Between z = 1.5 and z = 3.5, our new measurements follow the lower envelop of the previous measurements. This mild disagreement could have several causes.

First of all, the clustering effect was not taken into account by the previous analyses based on stacking. This effect is stronger at high redshift, because the bias5 of both infrared and mass-selected galaxies increases with redshift (e.g., Béthermin et al. 2013). In addition, the SEDs peak at a longer wavelength, where the bias is stronger owing to beam size (see Sect. A.1). The tension with the results based on UV-detected galaxies could be explained by a slight incompleteness of the UV-detected samples at low sSFR or a small overestimate of the dust corrections. There could also be effects caused by the different techniques and assumptions used to determine the stellar masses in the various fields (star formation histories, PSF-homogenized photometry or not, presence of nebular emission in the highest redshift bins, template libraries, etc.). Finally, this difference could also be an effect of the variance. These discrepancies on the estimates of sSFRs will be discussed in detail in Schreiber et al. (2015).

thumbnail Fig. 6

Evolution of the mean sSFR in main-sequence galaxies (blue triangles) and strong starbursts (red squares). The gray diamonds are a compilation of measurements at the same mass performed by Sargent et al. (2014). The blue line is the best fit to our data.

Open with DEXTER

4.3. Evolution of the mean intensity of the radiation field

thumbnail Fig. 7

Evolution of the mean intensity of the radiation field U in main-sequence galaxies (blue triangles) and strong starbursts (red squares). The black diamonds are the measurements presented in Magdis et al. (2012a) based on a similar analysis but in the GOODS fields. The orange asterisk is the mean value found for the local ULIRG sample of da Cunha et al. (2008; see also Magdis et al. 2012a). The black circle is the average value in HRS galaxies (Ciesla et al. 2014). The solid and dashed lines represent the evolutionary trends expected for a broken and universal FMR, respectively (see Sect. 4.3). The blue dotted line is the best fit of the evolution of the main-sequence galaxies ((3.0 ± 1.1) × (1 + z)1.8 ± 0.4) and the red dotter line the best fit of the strong starburst data by a constant (31 ± 3).

Open with DEXTER

The evolution of the mean intensity of the radiation field has different trends in main-sequence galaxies than in strong starbursts (see Fig. 7). This quantity is strongly correlated to the temperature of the dust. We found a rising U with increasing redshift in main-sequence galaxies up to z = 4 ((3.0 ± 1.1) × (1 + z)1.8 ± 0.4), confirming and extending the finding of Magdis et al. (2012a) at higher redshift. Other studies (e.g., Magnelli et al. 2013; Genzel et al. 2014) found an increase of the dust temperature with redshift in mass-selected samples.

The evolution of U we found can be understood from a few simple assumptions on the evolution of the gas metallicity and the (SFE) of galaxies. As shown by Magdis et al. (2012a), U is proportional to LIR/Mdust. We can also assume that (1)where the left-side of the proportionality is the well-established Kennicutt (1998) relation. The right-side of the proportionality is the integrated version of the Schmidt-Kennicutt relation which links the SFR to the mass of molecular gas in a galaxy (Mmol). Sargent et al. (2014) found a best-fit value for s of 0.83 compiling a large set of public data about low- and high-redshift main-sequence galaxies. The molecular gas mass can also be connected to the gas metallicity Z and the dust mass (e.g., Leroy et al. 2011; Magdis et al. 2012a), (2)where Z(M,SFR) is the gas metallicity which can be connected to M and SFR through the fundamental metallicity relation (FMR, Mannucci et al. 2010). There is recent evidence showing that this relation breaks down at high redshifts. For instance, Troncoso et al. (2014) measured a ~0.5 dex lower normalization at z ~ 3.4 compared to the functional form of the FMR at low redshift. Amorín et al. (2014) found the same offset in a lensed galaxy at z = 3.417. At z ~ 2.3, Steidel et al. (2014, see also Cullen et al. 2014 found an offset of 0.34–0.38 dex in the mass-metallicity relation and only half of this difference can be explained by the increase of SFR at fixed stellar mass using the FMR. Finally, a break in the metallicity relation is also observed in low mass (log(M/M) ~ 8.5) damped Lyman α absorbers around z = 2.6 (Møller et al. 2013). In our computations, we consider two different relations: a universal FMR where metallicity depends only on M and SFR, and a FMR relation with a correction of 0.30 × (1.7 − z) dex at z> 1.7 (hereafter broken FMR), which agrees with the measurements cited previously. Combining these expressions, we can obtain the following evolution: (3)We computed the expected evolution of U using the fit to the evolution of sSFR presented in Sect. 4.2 and assuming the mean stellar mass of our sample is 6 × 1010M, the average mass of the main-sequence sample6. We used the value of Magdis et al. (2012a) at z = 0 to normalize our model. The results are presented in Fig. 7 for a universal and a broken FMR. The broken FMR is compatible with all of our data points at 1σ. The universal FMR implies a significant underestimation of U at high redshifts (3 and 2σ in the two highest redshift bins).

We checked that the dust heating by the cosmic microwave background (CMB) is not responsible for the quick rise the quick rise in main-sequence galaxies. The CMB temperature at z = 4 is 13.5 K. The dust temperature that our high-redshift galaxies would have for a virtually z = 0 CMB temperature, , is estimated following da Cunha et al. (2013)(4)where is the temperature of the CMB at z = 0 and is the measured dust temperature at high redshift. This temperature is estimated fitting a gray body with an emissivity of β = 1.8 to our photometric measurements at λrest> 50μm. The CMB has a relative impact which is lower than 2 × 10-4 at all redshifts and thus this effect is negligible. These values are small compared to da Cunha et al. (2013), who assumed a dust temperature of 18 K. The warmer dust temperatures we measured suggests that the CMB should be less problematic than anticipated.

Concerning the evolution of U in strong starbursts, we found no evidence of evolution ((1 + z)− 0.1 ± 1.0) and our results can be fitted by a constant U of 31 ± 3. Our value of U at 0.5 <z< 3 is similar to the measurements on a sample of local ULIRGs (da Cunha et al. 2008). This suggests that high-redshift strong starbursts are a more extended version of the nuclei of local ULIRGs, as also suggested by the semi-analytical model of Lagos et al. (2012). At z ~ 2.5, the main-sequence galaxies and the strong starbursts have similar U values. However, we do not interpret the origins of these high values of U in the same way (see Sects. 4.4, 4.5, and 5). At z> 2.5, we cannot constrain with our analysis if U in strong starbursts rises as in main-sequence galaxies or stays constant.

4.4. Evolution of the ratio between dust and stellar mass

We also studied the evolution of the mean ratio between the dust and the stellar mass in the main-sequence galaxies and the strong starbursts. The results are presented Fig. 8. In main-sequence galaxies, this dust-to-stellar-mass ratio rises up to z ~ 1 and flattens above this redshift. Strong starbursts typically have 5 times higher ratio. Our measurements are compatible within 2σ with the slowly rising trend of (1 + z)0.05 found by Tan et al. (2014) for a compilation of individual starbursts. However, our data favors a steeper slope.

thumbnail Fig. 8

Mean ratio between dust and stellar mass as a function of redshift in main-sequence galaxies (blue triangles) and strong starbursts (red squares). The orange asterisk is the mean value found for the local ULIRG sample of da Cunha et al. (2008; see Magdis et al. 2012a). The black circle is the average value in HRS galaxies (Ciesla et al. 2014). The solid and dashed lines represent the evolutionary trends expected for a broken and universal FMR, respectively (see Sect. 4.3). The red dot-dashed line is the best-fit of the evolution found for a sample of individually-detected starbursts of Tan et al. (2014). The predictions of the models of Lagos et al. (2012) and Lacey et al. (in prep.) after applying the same mass cut and sSFR selection are overplotted with a three-dot-dash line and a long-dash line, respectively, with the same color code as the symbols.

Open with DEXTER

We modeled the evolution of this ratio in main-sequence galaxies using the same simple considerations as in Sect. 4.3. The evolution of the mean dust-to-stellar-mass ratio can be written as (5)One can see that Mdust/M is the result of a competition between the rising SFR with increasing redshift and the decreasing gas metallicity. The results are compatible with the broken FMR at 1σ. The relation obtained with the universal FMR rises too rapidly at high redshift.

We also compared our results with predictions of two semi-analytical models. The Lagos et al. (2012) and Lacey et al. (in prep.) models are based on GALFORM. The main difference between these two models is that (Lagos et al. 2012) adopt a universal IMF (a Galactic-like IMF; Kennicutt 1983), while Lacey et al. (in prep.) adopt a non-universal IMF. In the latter star formation taking place in galaxy disks has a Galactic-like IMF, while starbursts have a more top-heavy IMF. This is done to reproduce the number counts of SMGs found by surveys.

We select galaxies in the models in the same way we do in the observations based on their stellar mass and distance from the main sequence. An important consideration is that to derive stellar masses in the observations we fix the IMF to a Chabrier IMF, which is different to the IMFs adopted in both models. In order to correct for this we multiply stellar masses in the Lagos et al. (2012) model by 1.1 to go from a Kennicutt IMF to a Chabrier IMF. However, this is non-trivial for the Lacey et al. (in prep.) model, since it adopts two different IMFs. In order to account for this we correct the fraction of the stellar mass that was formed in the disk by the same factor of 1.1, and divide the fraction of stellar mass that was formed during starbursts by 2. The latter factor is taken as an approximation to go from their adopted top-heavy IMF to a Chabrier IMF, but this conversion is not necessarily unique, and it depends on the dust extinction and stellar age (see Mitchell et al. 2013 for details). In this paper we make a unique correction, but warn the reader that a more accurate approach would be to perform SED fitting to the predicted SEDs of galaxies and calculating the stellar mass in the same way we would do for observations.

Compared to the observations of main-sequence galaxies, the Lagos et al. (2012) model reproduces observations well in the redshift range 1 <z< 3, while at z< 1 and z> 3 it overpredicts the dust-to-stellar mass ratio. There are different ways to explain the high dust-to-stellar mass ratios: high gas metallicities, high gas masses or stellar masses being too low for the dust masses. In the case of the Lagos et al. (2012) model the high dust-to-stellar mass ratios are most likely coming from massive galaxies being too gas rich since their metallicities are close to solar, which is what we observe in local galaxies of the same stellar mass range. The Lacey et al. (in prep.) model predicts dust-to-stellar mass ratios that are twice too high compared to the observations in the whole redshift range. In this case this is because the gas metallicities of MS galaxies in the Lacey model are predicted to be supersolar on average (close to twice the solar metallicity, 12 + log(O/H)~9.0), resulting in dust masses that are higher than observed.

In the case of starbursts, the high values inferred for the dust-to-stellar mass ratio in the observations are difficult to interpret. The Lagos et al. (2012) model underpredicts this quantity by a factor of ~5 and the Lacey et al. (in prep.) model by a factor of ~2. At first the ratio of 1.52% inferred in the observations seems unphysical. However, since the gas fraction (defined here as Mmol/ (Mmol + M)) in these high-redshift starbursts is around 50% (see Sect. 4.5, but also, e.g., Riechers et al. 2013 and Fu et al. 2013), the high values observed for the dust-to-stellar mass ratio can be reached if the gas-to-dust ratio is 5067. Values similar to the latter are observed in metal-rich galaxies (12 + log(O/H)~9, e.g., Rémy-Ruyer et al. 2014). This high metal enrichment in strong starbursts compared to main-sequence galaxies could be explained by several mechanisms:

  • the transformation of gas into stars is quicker and the metals are notdiluted by the accretion of pristine gas;

  • a fraction of the external layers of low-metallicity gas far from the regions of star formation could be ejected by the strong outflows caused by these extreme starbursts;

  • a top-heavy IMF could produce quickly lots of metals through massive stars without increasing too rapidly the total stellar mass because of mass losses.

This high ratio in strong starbursts is discussed in details in Tan et al. (2014).

When it comes to the comparison with the models, one can understand the lower dust-to-stellar mass ratios predicted by the model as resulting from the predicted gas metallicities. Lagos et al. (2012) predict that the average gas metallicity in strong starbursts is close to 0.4 solar metallicities (12 + log(O/H)~8.3), which is about 4 times lower than we can infer from a gas-to-dust mass ratio of 50 (see previous paragraph). While the Lacey et al. (in prep.) model predicts gas metallicities for starbursts that are on average close to solar metallicity (12 + log(O/H)~8.7), 2 times too low for the inferred metallicity of the strong starbursts we observe. We note that both models predict main sequence galaxies having higher metallicities than bright starbursts of the same stellar masses. This seems to contradict the observations and may be at the heart of why the models struggle to get the dust-to-stellar mass ratios of both the main sequence and starburst populations at the same time.

4.5. Evolution of the fraction of molecular gas

Finally, we deduced the mean mass of molecular gas from the dust mass using the same method following Magdis et al. (2011) and Magdis et al. (2012a). They assumed that the gas-to-dust ratio depends only on gas metallicity and used the local relation of Leroy et al. (2011)7: (6)Given the relatively high stellar mass of our samples, and the rising gas masses and ISM pressures to high redshifts (Obreschkow & Rawlings 2009), we expect the contribution of atomic hydrogen to the total gas mass to be negligible and we neglect it in the rest of the paper, considering total gas mass or molecular gas mass to be equivalent. For main-sequence galaxies, the gas metallicity is estimated using the FMR as explained in Sect. 4.3. We converted the values provided by the FMR from the KD02 to the PP04 metallicity scale using the prescriptions of Kewley & Ellison (2008) before using it in Eq. (6).

The gas metallicity in strong starbursts cannot be estimated using the FMR. Indeed, this relation predicts that, at fixed stellar mass, objects forming more stars are less metallic. This effect is expected in gas regulated systems, because a higher accretion of pristine gas involves a stronger SFR, but also a dilution of metals (e.g., Lilly et al. 2013). This phenomenon is not expected to happen in starbursts, since their high SFRs are not caused by an excess of accretion, but more likely by a major merger. These high-redshift starbursts are probably progenitors of current, massive, elliptical galaxies (e.g., Toft et al. 2014). We thus assumed that their gas metallicity is similar and used a value of 12 + log(O/H) = 9.1±0.2 (see a detailed discussion in Magdis et al. 2011 and Magdis et al. 2012a).

We then derived the molecular gas fraction in main-sequence galaxies, defined in this paper as Mmol/ (M + Mmol). The results are presented in Fig. 9. We found a quick rise up to z ~ 2. At higher redshifts, the recovered trend depends on the assumptions on the gas metallicity. The rise of the gas fraction in main-sequence galaxies continues at higher redshift if we assume the broken FMR favored by the recent studies, but flattens with a universal FMR. If the broken FMR scenario is confirmed, there could thus be no flattening or reversal of the molecular gas fraction at z> 2 contrary to what is claimed in Magdis et al. (2012b), Saintonge et al. (2013), and Tan et al. (2013). Our estimations agree with the previous estimates of Magdis et al. (2012a) at z = 1, but are 1σ lower at z = 2, because the bias introduced by clustering was corrected in our study. Our results also agree at 1σ with the analysis of Santini et al. (2014) at the same stellar mass up to z = 2.5 after converting the stellar mass from a Salpeter to a Chabrier IMF convention. However, our estimates are systematically higher than theirs and agree better with the CO data. Our measurements also agree with the compilation of CO measurements of Saintonge et al. (2013) and the two galaxies studied at z ~ 3 by Magdis et al. (2012b). These measurements are dependent on the assumed αCO conversion factor, and on the completeness corrections. The good agreement with this independent method is thus an interesting clue to the reliability of these two techniques.

Strong starbursts have molecular gas fractions 1σ higher than main-sequence galaxies, but follows the same trend. Sargent et al. (2014) predicted that starbursts on average should have a deficit of gas compared to the main sequence (but that gas fraction are expected to rise continuously as the sSFR-excess with respect to the MS increases). Here we selected only the most extreme starbursts with an excess of sSFR of a factor of 10 instead of the average value of ~4. These extreme starbursts may only be possible by the mergers of two gas-rich galaxies galaxies already above the main-sequence before the merger. This could explain this small positive offset compared to the main-sequence sample.

thumbnail Fig. 9

Evolution of the mean molecular gas fraction in massive galaxies (>3 × 1010M). The starbursts are represented by red squares and the main-sequence galaxies by blue triangles or light blue diamonds depending on wether the gas fraction is estimated using a broken or an universal FMR, respectively. These results are compared with previous estimate using dust masses of Magdis et al. (2012a, black plus) and Santini et al. (2014, gray area), using CO for two z> 3 galaxies (Magdis et al. 2012b, black crosses), and the compilation of CO measurements of Saintonge et al. (2013, black asterisks). The predictions of the models of Lagos et al. (2012) and Lacey et al. (in prep.) for the same mass cut are overplotted with a three-dot-dash line and a long-dash line, respectively.

Open with DEXTER

We also compared our results with the models of Lagos et al. (2012) and Lacey et al. (in prep.) presented in Sect. 4.4. Both models agree well with our measurements of the gas fraction for starburst galaxies at all redshifts and main-sequence galaxies at 1.5 <z< 3. Both the Lagos et al. (2012) and Lacey et al. (in prep.) models overpredict the molecular gas fraction at z< 0.5 at a 12σ level. At reshifts z> 3, the Lacey et al. (in prep.) model agrees better with the universal FMR scenario at z> 3, while the Lagos et al. (2012) model is more compatible with the broken FMR. The fact that both models predict molecular gas fractions that in overall agree with the observations supports our interpretation in Sect. 4.4, which points to the model of metal enrichment as the source of discrepancy in the dust-to-stellar mass ratios.

4.6. Evolution of the depletion time

We estimated the mean depletion time of the molecular gas, defined in our analysis as the ratio between the mass of molecular gas and the SFR. Figure 10 shows our results. The depletion time in strong starbursts does not evolve with redshift and is compatible with 100 Myr, the typical timescale of the strong boost of star formation induced by major mergers (e.g., Di Matteo et al. 2008). This timescale is longer in main-sequence galaxies and slightly (1σ) evolves with redshift at z< 1. It decreases from 1.3 Gyr at z ~ 0.375 to ~500 Myr around z ~ 1.5 and is stable at higher redshift in the case of a broken FMR (but continues to decrease with redshift for a universal FMR). This timescale is similar to the maximum duration high-redshift massive galaxies can stay on the main-sequence before reaching the quenching mass around 1011M (Heinis et al. 2014). The mass of molecular gas and stars contained in these high-redshift objects is already sufficient to reach this quenching mass without any additional accretion of gas.

thumbnail Fig. 10

Evolution of the mean molecular gas depletion time. The symbols are the same as in Fig. 9.

Open with DEXTER

5. Discussion

5.1. What is the main driver of the strong evolution of the specific star formation rate?

thumbnail Fig. 11

Relation between the mean SFR rate and the mean molecular gas mass in our galaxy samples, i.e., integrated Kennicutt-Schmidt law. The solid line and the dashed line are the center of the relation fitted by Sargent et al. (2014) on a compilation of data for main-sequence galaxies and starbursts, respectively. The dotted lines represent the 1σ uncertainties on these relations. The triangles and diamonds represent the average position of massive, main-sequence galaxies in this diagram assuming a broken FMR and an universal FMR, respectively. The squares indicates the average position of strong starbursts.

Open with DEXTER

We checked the average position of our selection of massive galaxies in the integrated Kennicutt-Schmidt diagram (SFR versus mass of molecular gas) to gain insight on their mode of star formation. In this diagram, normal star-forming galaxies and starbursts follow two distinct sequences. For comparison, we used the fit of a recent data compilation performed by Sargent et al. (2014). The results are presented in Fig. 11.

The average position of our sample of strong starbursts is in the 1σ confidence region of Sargent et al. (2014) for starbursts. They are systematically below the central relation, but the uncertainty is dominated by the systematic uncertainties on their gas metallicity. In addition, Sargent et al. (2014) suggested that the SFEs of starbursts follow a continuum of values depending on their boost of sSFR. Our objects are thus not expected to be exactly on the central relation. The interpretation of the results for main-sequence galaxies is dependent on the hypothesis on the gas metallicity. In the scenario of a broken FMR favored by recent observations, the average position of main-sequence galaxies at all redshifts falls on the relation of normal star-forming galaxies. This suggests that the star formation is dominated by galaxies forming their stars through a normal mode at all redshifts below z = 4. In the case of a non-evolving FMR, the massive high-redshift galaxies do not stay on the normal star-forming sequence and have higher SFEs.

If the scenario of a broken FMR favored by the most recent observations is consolidated, the strong star-formation observed in massive high-redshift galaxies would thus be caused by huge gas reservoirs probably fed by an intense cosmological accretion. This strong accretion of primordial gas dilute the metals produced by the massive stars (e.g., Bouché et al. 2010; Lilly et al. 2013). Consequently, the gas-to-dust ratio is much lower at high redshift than at low redshift. Since the SFE is only slowly evolving (), the number of UV photons absorbed per mass of dust is thus higher and the dust temperature is warmer as observed in our analysis (see Sect. 4.3). This scenario provides thus a consistent interpretation of evolution of both the sSFR and the dust temperature of massive galaxies with redshift.

5.2. Limitations of our approach

Our analysis provided suggestive results. However, it relies on several hypotheses, which cannot be extensively tested yet. In this section, we discuss the potential limitation of our analysis.

The evolution of the metallicity relations at z> 2.5 was measured only by a few pioneering works, which found that the normalization of the FMR evolves at z> 2.5. We used a simple renormalization depending on redshift to take this evolution into account. The redshift sampling of these studies is relatively coarse and we used a simple linear evolution with redshift. Future studies based on larger samples will allow a finer sampling of the evolution of the gas metallicity in massive galaxies at high redshift. However, the current results are very encouraging. The current assumption of a broken FMR allows us to recover naturally both the evolution of the U parameter and the integrated Schmidt-Kennicutt relation at high redshift.

The gas metallicity of strong starbursts was more problematic to set. We can reasonably guess it assuming they are progenitors of the most massive galaxies. However, direct measurements of their gas metallicity are difficult to perform using optical/near-IR spectroscopy because of their strong dust attenuation. The millimeter spectroscopy of fine-structure lines with ALMA will be certainly an interesting way to determine the distribution of gas metallicity of strong starbursts in the future (e.g., Nagao et al. 2011).

The validity of the calibration of the gas-to-dust ratio versus gas metallicity relation in most extreme environment is also uncertain and difficult to test with the current data sets. Saintonge et al. (2013) found an offset of a factor 1.7 for a population of lensed galaxies and discussed the possible origins of the tension between the gas content estimated from CO and from dust. However, we found no offset with the integrated Kennicutt-Schmidt relation in our analysis and a good agreement with the compilation of CO measurements of gas fractions. The lensed galaxies of Saintonge et al. (2013) could be a peculiar population because they are UV-selected and then biased toward dust-poor systems. They could also be affected differential magnification effects or Herschel-selection biases. The hypotheses performed to estimate the gas metallicity are also different between their and our analysis (standard mass-metallicity relation versus broken FMR).

Finally, the stacking analysis only provides an average measurement of a full population. Thus it is difficult to estimate the heterogeneity of the stacked populations. Bootstrap techniques can be applied to estimate the scatter on the flux density at a given wavelength (Béthermin et al. 2012b). However, because of the correlation between U and Md, this technique cannot be applied to measure the scatter on each of these parameters.

6. Conclusion

We used a stacking analysis to measure the evolution of the average mid-infrared to millimeter emission of massive star-forming galaxies up to z = 4. We then derived the evolution of the mean physical parameters of these objects. Our main findings are the following.

  • The mean intensity of the radiation fieldU in main-sequence galaxies, which is strongly correlated with their dust temperature, rises rapidly with redshift: U ⟩ = (3.0 ± 1.1) × (1 + z)1.8 ± 0.4. This evolution can be interpreted considering the decrease in the gas metallicity of galaxies at constant stellar mass with increasing redshift. We found no evidence for an evolution of U in strong starbursts up to z = 3.

  • The mean ratio between the dust mass and the stellar mass in main-sequence galaxies rises between z = 0 and z = 1 and exhibit a plateau at higher redshift. The strong starbursts have a higher ratio by a factor of 5.

  • The average fraction of molecular gas (Mmol/ (M + Mmol)) rises rapidly with redshift and reaches ~60% at z = 4. A similar evolution is found in strong starbursts, but with slightly higher values. These results agree with the pilot CO surveys performed at high redshift.

  • We compare with two state-of-the-art semi-analytic models that adopt either a universal IMF or a top-heavy IMF in starbursts and find that the models predict molecular gas fractions that agree well with the observations but the predicted dust-to-stellar mass ratios are either too high or too low. We interpret this as being due to the way metal enrichment is dealt with in the simulations. We suggest different mechanisms that can help overcome this issue. For instance, outflows affecting more metal depleted gas that is in the outer parts of galaxies.

  • The average position of the massive main-sequence galaxies in the integrated Kennicutt-Schmidt diagram corresponds to the sequence of normal star-forming galaxies. This suggests that the bulk of the star-formation up to z ~ 4 is dominated by the normal mode of star-formation and that the extreme SFR observed are caused by huge gas reservoirs probably induced by the very intense cosmological accretion. The strong starbursts follow another sequence with a 5–10 times higher star-formation efficiency.


1

Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.

2

The minimum expected flux for our mass-selected sample of strong starbursts is computed using the three-dot-dash curve in Fig. 2 and the Magdis et al. (2012a) starburst template.

3

APEX project IDs: 080.A–3056(A), 082.A–0815(A) and 086.A–0749(A).

4

More details on the CRUSH settings can be found at: http://www.submm.caltech.edu/~sharc/crush/v2/README

5

The bias b is defined by wgal = b2wDM, where wgal and wDM are the projected two-point correlation function of galaxies and dark matter, respectively. The higher the bias is, the stronger is the clustering density of galaxies compared to dark matter.

6

We could have used the mean stellar masses in each redshift bin provided in Table 2. However, assuming a single value of the stellar mass at all redshift has a negligible impact on the results and the tracks are smoother.

7

Converted to PP04 convention.

Acknowledgments

We thank the anonymous referee for providing constructive comments. We acknowledge Morgane Cousin, Nick Lee, Nick Scoville, and Christian Maier for their interesting discussions/suggestions, Laure Ciesla for providing an electronic table of the physical properties of the HRS sample, and Amélie Saintonge for providing her compilation of data. We gratefully acknowledge the contributions of the entire COSMOS collaboration consisting of more than 100 scientists. The HST COSMOS program was supported through NASA grant HST-GO-09822. More information on the COSMOS survey is available at http://www.astro.caltech.edu/cosmos. ased on data obtained from the ESO Science Archive Facility. Based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under ESO programme ID 179.A-2005 and on data products produced by TERAPIX and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium. M.B., E.D., and M.S. acknowledge the support of the ERC-StG UPGAL 240039 and ANR-08-JCJC-0008 grants. A.K. acknowledges support by the Collaborative Research Council 956, sub-project A1, funded by the Deutsche Forschungsgemeinschaft (DFG).

References

Online material

Appendix A: Estimation and correction on the bias caused by the galaxy clustering on the stacking results

As explained in Sect. 3.2, the standard stacking technique can be strongly affected by the bias caused by the clustering of the galaxies. We use two independent methods to estimate and correct it.

Appendix A.1: Estimation of the bias using a simulation based on the real catalog

We performed an estimate of the bias induced by the clustering using a realistic simulation of the COSMOS field based on the positions and stellar masses of the real sources. The flux of each source in this simulation is estimated using the ratio between the mean far-IR/(sub-)mm fluxes and the stellar mass found by a first stacking analysis. The galaxies classified as passive are not taken into account in this simulation. This technique assumes implicitly a flat sSFR-M relation, since we use a constant SFR/M ratio versus stellar mass at fixed redshift. However, we checked that using a more standard sSFR relation (e.g., Rodighiero et al. 2011) has a negligible impact on the results. We applied no scatter around this relation in our simulation for simplicity. As mean stacking is a linear operation, the presence or not of a scatter has no impact on the results (Béthermin et al. 2012b).

A simulated map is thus produced using all the star-forming galaxies of the Ilbert et al. (2013) catalog. In order to avoid edge effects (absence of sources and thus a lower background caused by the faint unresolved sources in the region covered by the optical/near-IR data), we fill the uncovered regions drawing with replacement sources from the UltraVISTA field and putting them at a random position. The number of drawn sources is chosen to have exactly the same number density inside and outside the UltraVISTA field.

Finally, we measured the mean fluxes of the M> 3 × 1010M sources by stacking in the simulated maps, using exactly the same photometric method as for the real data. We finally computed the relative bias between the recovered flux and the input flux (Sout/Sin − 1). The results are shown Fig. A.1 (blue triangles). The uncertainties are computed a bootstrap method. As expected, the bias increases with the size of the beam. We can see a rise of the bias with redshift up to z ~ 2. This trend can be understood considering the rise of the clustering of the galaxy responsible for the cosmic infrared background (Planck Collaboration XXX 2014) and a rather stable number density of emitters especially below z = 1 (Béthermin et al. 2011; Magnelli et al. 2013; Gruppioni et al. 2013). At higher redshift, we found a slow decrease. This trend is probably driven by the decrease in the infrared luminosity density at high redshift (Planck Collaboration XXX 2014; Burgarella et al. 2013) combined with the decrease in the number density of infrared emitters (Gruppioni et al. 2013).

thumbnail Fig. A.1

Relative bias induced by the clustering as a function of redshift at the various wavelengths we used in our analysis. The FWHM of the beam is provided in brackets. The blue triangles are the estimations from the simulation (Sect. A.1) and the red diamonds are provided by the fit of the clustering component in map space (Sect. A.2). These numbers are only valid for a complete sample of M> 3 × 1010M galaxies.

Open with DEXTER

Appendix A.2: Estimation of the bias fitting the clustering contribution in the stacked images

The method presented in the previous section only takes into account the contamination of the stacks by known sources. However, faint galaxy populations could have a non negligible contribution, despite their total contribution to the infrared luminosity and their clustering are expected to be small. We thus used a second method to estimate the bias caused by the clustering which takes into account a potential contamination by these low-mass galaxies. This method is based on a simultaneous fit in the stacked images of three components: a point source at the center of the image, a clustering contamination, and a background. This technique was already successfully used by several previous works based on Herschel and Planck data (Béthermin et al. 2012b; Heinis et al. 2013, 2014; Welikala et al. 2014).

In presence of clustering, the outcome of a stacking is not only a PSF with the mean flux of the population and a constant background (corresponding to the surface brightness of all galaxy populations i.e., the cosmic infrared background). There is in addition a signal coming from the greater probability of finding another neighboring infrared galaxy compared to the field because of galaxy clustering. The signal in the stacked image can thus be modeled by Bavouzet (2008) and Béthermin et al. (2010b)(A.1)where m is the stacked image, PSF the point spread function, and w the auto-correlation function. The symbol represents the convolution. α, β, and γ are free parameters corresponding to the intensity of the mean flux of the population, the clustering signal, and the background, respectively. This method works only if the PSF is well-known, the extension of the sources is negligible compared to the PSF, and the effects of the filtering are small at the scale of the stacked image. Consequently, we applied this method only to the SPIRE data for which these hypotheses are the most solid. The uncertainties on the clustering bias (β/α for the photometry we chose to use for SPIRE data) are estimated fitting the model described previously on a set of stacked images produced from 1000 bootstrap samples. The results are shown in Fig. A.1 (red diamonds).

Appendix A.3: Corrections of the measurements

In Fig. A.1, we can see that the two methods provide globally consistent estimates. This confirms that the low-mass galaxies not included in the UltraVISTA catalog have a minor impact. We found few outliers for which the two methods disagree. In particular, in the 1.5 <z< 1.75 bin, the estimation from the simulation is higher than the trend of the redshift evolution at all wavelengths, and the results from the profile fitting are lower. This could be caused, as instance, by a structures in the field or a systematic effect in the photometric redshift. Because of these few catastrophic outliers, we chose to use a correction computed

from a fit of the redshift evolution of the bias instead of an individual estimate in each redshift slice.

The evolution of the bias with redshift is fitted independently at each wavelength. We chose to use a simple, second-order, polynomial model (az2 + bz + c). We used only the results from the simulation to have a consistent treatment of the various wavelengths. The scatter of the residuals is larger than the residuals, probably because bootstrap does not take into account the variance coming from the large-scale structures. We thus used the scatter of the residuals to obtain a conservative estimate of the uncertainties on the bias. In Fig. A.1, the best fit is represented by a solid line and the 1σ confidence region by a dashed line.

In a few case, the bias at z> 3 can converge to unphysical negative values. We then apply no corrections, but combine the typical uncertainty on the bias to the error bars. A special treatment is also applied to the samples of strong starbursts. Their flux is typically 10 times brighter in infrared by construction (their sSFR is 10 times larger than the main sequence). In contrast, the clustering signal is not expected to be significantly stronger, because the clustering of massive starbursts and main-sequence galaxies is relatively similar (Béthermin et al. 2014). We thus divide the bias found for the full population of galaxy by a factor of 10 to estimate the one of the starbursts for simplicity.

Appendix A.4: Testing another method

We also tried to apply the simstack algorithm (Viero et al. 2013) to our data. This algorithm is adapted from Kurczynski & Gawiser (2010) and uses the position of the known sources to deblend their contamination. Contrary to Kurczynski & Gawiser (2010), simstack can consider a large set of distinct galaxy populations. The mean flux of the each population is used to estimate how sources contaminate their neighbors. All populations are treated simultaneously. This is the equivalent of PSF-fitting codes but applied to a full population instead of each source individually. Unfortunately, this method is not totally unbiased in our case. We found biases up to 15% running simstack on the simulation presented in Sect. A.1, probably because the catalog of mass-selected sources is not available around bright sources. At the edge of the optical/near-IR-covered region, the flux coming from the sources outside the covered area is not corrected, when the flux from all neighbors is taken into account at the middle of zone where the mass catalog is extracted. Indeed, the algorithm works correctly if we put on the simulation only sources present in the input catalog.

Appendix B: Fit residuals

Figures B.1 and B.2 shows the residuals of the fits of our mean SEDs derived by stacking. We did not find any systematic trend, except a 2σ underestimation of the millimeter data in main-sequence galaxies at z> 3.

thumbnail Fig. B.1

Residuals of our fit of mean SEDs of main-sequence galaxies by the Draine & Li (2007) model.

Open with DEXTER

thumbnail Fig. B.2

Residuals of our fit of mean SEDs of strong starbursts by the Draine & Li (2007) model.

Open with DEXTER

All Tables

Table 1

Summary of our flux density measurements by stacking.

Table 2

Summary of the average physical parameters of our samples.

All Figures

thumbnail Fig. 1

Stellar mass distribution of our samples of star-forming galaxies in the various redshift bins we used. Only galaxies more massive than our cut of 3 × 1010M are represented. The first bin contain fewer objects than the second one because our cut fall at the middle of the first one. The arrows indicate the mean stellar mass in each redshift bin.

Open with DEXTER
In the text
thumbnail Fig. 2

The thick red solid line represents the luminosity limit corresponding to a criterion of a 5σ detection in at least two Herschel bands. The other solid lines are the limits for a detection at only one given wavelength (purple for 100 μm, blue for 160 μm, turquoise for 250 μm, green for 350 μm, orange for 500 μm). The dashed, dot-dash, and three-dot-dash lines indicate the infrared luminosity of a galaxie of 3 × 1010M (our mass cut) at the center of the main sequence, a factor of 4 above it, and a factor of 10 above it, respectively.

Open with DEXTER
In the text
thumbnail Fig. 3

Mean flux density as a function of wavelength (observed wavelength in the top panels and rest-frame wavelength in the bottom panels) at various redshifts (see color coding). The left panels show the mean SEDs of the main-sequence sample and the right panels those of the strong starbursts.

Open with DEXTER
In the text
thumbnail Fig. 4

Rest-frame mean spectral energy distribution of our selection of massive, star-forming galaxies at various redshift measured by stacking analysis. The data points are fitted using the Draine & Li (2007) model. This model is convolved with the redshift distribution of the sources before being compared to the data. The black and blue lines represent the intrinsic and convolved SEDs, respectively. The bottom right corner summarizes the redshift evolution seen in our data.

Open with DEXTER
In the text
thumbnail Fig. 5

Rest-frame mean spectral energy distribution of our selection of strong starbursts at various redshift measured by stacking analysis. The data points are fitted using the Draine & Li (2007) model. This model is convolved with the redshift distribution of the sources before being compared to the data. The black and red lines show the intrinsic and convolved SEDs, respectively.

Open with DEXTER
In the text
thumbnail Fig. 6

Evolution of the mean sSFR in main-sequence galaxies (blue triangles) and strong starbursts (red squares). The gray diamonds are a compilation of measurements at the same mass performed by Sargent et al. (2014). The blue line is the best fit to our data.

Open with DEXTER
In the text
thumbnail Fig. 7

Evolution of the mean intensity of the radiation field U in main-sequence galaxies (blue triangles) and strong starbursts (red squares). The black diamonds are the measurements presented in Magdis et al. (2012a) based on a similar analysis but in the GOODS fields. The orange asterisk is the mean value found for the local ULIRG sample of da Cunha et al. (2008; see also Magdis et al. 2012a). The black circle is the average value in HRS galaxies (Ciesla et al. 2014). The solid and dashed lines represent the evolutionary trends expected for a broken and universal FMR, respectively (see Sect. 4.3). The blue dotted line is the best fit of the evolution of the main-sequence galaxies ((3.0 ± 1.1) × (1 + z)1.8 ± 0.4) and the red dotter line the best fit of the strong starburst data by a constant (31 ± 3).

Open with DEXTER
In the text
thumbnail Fig. 8

Mean ratio between dust and stellar mass as a function of redshift in main-sequence galaxies (blue triangles) and strong starbursts (red squares). The orange asterisk is the mean value found for the local ULIRG sample of da Cunha et al. (2008; see Magdis et al. 2012a). The black circle is the average value in HRS galaxies (Ciesla et al. 2014). The solid and dashed lines represent the evolutionary trends expected for a broken and universal FMR, respectively (see Sect. 4.3). The red dot-dashed line is the best-fit of the evolution found for a sample of individually-detected starbursts of Tan et al. (2014). The predictions of the models of Lagos et al. (2012) and Lacey et al. (in prep.) after applying the same mass cut and sSFR selection are overplotted with a three-dot-dash line and a long-dash line, respectively, with the same color code as the symbols.

Open with DEXTER
In the text
thumbnail Fig. 9

Evolution of the mean molecular gas fraction in massive galaxies (>3 × 1010M). The starbursts are represented by red squares and the main-sequence galaxies by blue triangles or light blue diamonds depending on wether the gas fraction is estimated using a broken or an universal FMR, respectively. These results are compared with previous estimate using dust masses of Magdis et al. (2012a, black plus) and Santini et al. (2014, gray area), using CO for two z> 3 galaxies (Magdis et al. 2012b, black crosses), and the compilation of CO measurements of Saintonge et al. (2013, black asterisks). The predictions of the models of Lagos et al. (2012) and Lacey et al. (in prep.) for the same mass cut are overplotted with a three-dot-dash line and a long-dash line, respectively.

Open with DEXTER
In the text
thumbnail Fig. 10

Evolution of the mean molecular gas depletion time. The symbols are the same as in Fig. 9.

Open with DEXTER
In the text
thumbnail Fig. 11

Relation between the mean SFR rate and the mean molecular gas mass in our galaxy samples, i.e., integrated Kennicutt-Schmidt law. The solid line and the dashed line are the center of the relation fitted by Sargent et al. (2014) on a compilation of data for main-sequence galaxies and starbursts, respectively. The dotted lines represent the 1σ uncertainties on these relations. The triangles and diamonds represent the average position of massive, main-sequence galaxies in this diagram assuming a broken FMR and an universal FMR, respectively. The squares indicates the average position of strong starbursts.

Open with DEXTER
In the text
thumbnail Fig. A.1

Relative bias induced by the clustering as a function of redshift at the various wavelengths we used in our analysis. The FWHM of the beam is provided in brackets. The blue triangles are the estimations from the simulation (Sect. A.1) and the red diamonds are provided by the fit of the clustering component in map space (Sect. A.2). These numbers are only valid for a complete sample of M> 3 × 1010M galaxies.

Open with DEXTER
In the text
thumbnail Fig. B.1

Residuals of our fit of mean SEDs of main-sequence galaxies by the Draine & Li (2007) model.

Open with DEXTER
In the text
thumbnail Fig. B.2

Residuals of our fit of mean SEDs of strong starbursts by the Draine & Li (2007) model.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.