© ESO, 2015

1. Introduction

Numerous observational results show a tight relationship between the stellar mass (M⋆) and the star formation rate (SFR) of star-forming galaxies (e.g. Noeske et al. 2007a; Elbaz et al. 2007; 2011; Daddi et al. 2007; Peng et al. 2010; Karim et al. 2011). The star-forming galaxies are distributed in the M⋆-SFR plane along what is commonly called the “star-forming main sequence”. If we do not consider quiescent galaxies, the existence of such a M⋆-SFR relation implies that the galaxies that are currently the most star-forming were also the most star-forming in their past history. Star-forming galaxies are scattered around this relation as expected from the stochasticity in their individual star formation histories (SFHs; e.g. Hopkins et al. 2014; Domínguez et al. 2014) and from the variety of possible SFHs. Extreme events like mergers could decouple the instantaneous SFR from the past star formation history and create outliers to the M⋆-SFR relation, which is one definition of starbursts (e.g. Rodighiero et al. 2011).

While the shape and the scatter of the M⋆-SFR relation already provide deep insights into the galaxy assembly process, its evolution along cosmic time is also of a great interest. Noeske et al. (2007a), Daddi et al. (2007) and Elbaz et al. (2007) find that the M⋆-SFR relation scales with cosmic time, such that the SFR increases with redshift at a given stellar mass. This evolution is also seen as an increase in the specific SFR (hereafter sSFR = SFR/M_⋆) with redshift at a given stellar mass. There is a growing consensus that the sSFR evolution is deeply linked to the hierarchical growth of dark matter structures (e.g. Bouché et al. 2010; Lilly et al. 2013). Assuming that galaxies are fed by fresh gas at a constant fraction of the averaged cosmological accretion rate, the sSFR should evolve as the specific dark matter increase rate (hereafter sMIR_DM) defined as Ṁ_H/M_H with M_H the mass of the dark matter halos (Lilly et al. 2013). However, the galaxies could be fed more efficiently in fresh gas at high redshift than in the local Universe since cold accretion occurs mainly at z> 2 (e.g. Dekel et al. 2009). Therefore, having an accurate characterization of the sSFR evolution with redshift is crucial in order to link the galaxy stellar mass assembly with the growth of the dark matter structures.

Below z< 1−1.5, the sSFR is relatively well measured using robust infrared data (Noeske et al. 2007a; Elbaz et al. 2007, 2011) and radio data (Karim et al. 2011). The sSFR increases steadily from the present day to z ~ 2 (e.g. Daddi et al. 2007; Karim et al. 2011) and its evolution is usually parametrized as sSFR ∝ (1 + z)^b. The values of b in the literature cover the full range between 2.5 and 5 (e.g. Speagle et al. 2014). Most studies assume a linear relation between log (sSFR) and log (M⋆) and characterize the slope and the scatter of this relation. Depending on the survey characteristics and on the SFR tracer, the value of the slope varies significantly in the literature. Noeske et al. (2007a) find a slope of −0.33^{± 0.08}, while several other studies obtain a value close to −0.1 (Elbaz et al. 2007; Daddi et al. 2007; Pannella et al. 2009; Peng et al. 2010). In the compilation of Speagle et al. (2014), the slope varies between −0.65 and −0.05 and depends on the SFR tracer. Some studies show that the slope could depend on the stellar mass and they even show a probable break in the M⋆-sSFR relation (e.g. Noeske et al. 2007a; Bauer et al. 2013; Lee et al. 2015; Whitaker et al. 2014). The parametrization assuming a linear relation between log (sSFR) and log (M⋆) is likely not valid over the full mass range. Finally, the scatter of the M⋆-sSFR relation is also debated in the literature. The scatter is ranging from 0.15 dex (Salmi et al. 2012) to 0.5 dex (Salim et al. 2007) and does not depend on the mass (e.g. Speagle et al. 2014; Lee et al. 2015). While studied in great detail, no consensus has been reached on the evolution, the scatter and the slope of the M⋆-sSFR relation.

Most of the analyses of the M⋆-sSFR relation are based on scatter diagrams (i.e. displays of the location of the individual sources in the M⋆-sSFR plane). However, this method does not provide any quantitative information on how galaxies are distributed around the median sSFR, and does not account for galaxies that could be under-sampled or missed by selection effects. In order to overcome this limitation, one should split the M⋆-sSFR plane in several mass bins and characterize the sSFR distribution in each bin correcting for selection effects. Then, accurate and robust information can be extracted from the analysis of the sSFR distribution.

The sSFR distribution has already been investigated in a few studies. Guo et al. (2013) produced the sSFR distributions per stellar mass bin in a sample similar to ours. Their study is limited to 0.6 <z< 0.8 while we want to explore a large redshift range in order to analyse the sSFR evolution. The work of Rodighiero et al. (2011) is also limited to one redshift slice at z ~ 2. Moreover, Rodighiero et al. (2011) need to rely on the UV light to trace the sSFR for the bulk of the star-forming population. Unfortunately, converting the UV light into SFR introduces uncertainties since it requires an estimate of the UV light absorbed by dust (e.g. Heinis et al. 2013; Rodighiero et al. 2014). By fitting a log-normal function over the sSFR distribution established by Rodighiero et al. (2011), Sargent et al. (2012) find . Moreover, Sargent et al. (2012) include in their fit a population of starbursts, i.e., galaxies having a sSFR higher than expected from the main sequence position. They estimate that 4% of the galaxies could be considered starbursts at z ~ 2. Finally, Kajisawa et al. (2010) measured the sSFR distribution per mass bin at 0.5 <z< 3 but with a sample limited in size.

Here, we estimate the sSFR functions, i.e., the number density of the galaxy per comoving volume and per sSFR bin. We measure the sSFR functions in four stellar mass bins from z = 0.2 to z = 1.4. In order to overcome the limitations of previous studies, we follow the following principles. First, our results rely on one robust SFR tracer, the 24 μm IR data obtained with the Multiband Imaging Photometer (MIPS) camera on board the Spitzer satellite. By limiting the analysis at z< 1.4, the galaxy L_IR can be derived with an accuracy better than 0.15 dex using the MIPS 24 μm data (Elbaz et al. 2010). The advantage of using a 24 μm selected sample is that we reach a lower SFR limit in comparison to a sample selected in one Herschel band. Since we apply one single cut in flux, we can easily correct for selection effect. Second, we limit our analysis to galaxy samples which are complete in stellar mass. These criteria allows us to consider only the SFR limit without having to consider an additional mass limit. Third, we combine the COSMOS (Scoville et al. 2007) and GOODS (Giavalisco et al. 2004) surveys. The large COSMOS area of 1.5 deg² allows us to get rare and massive star-forming sources, while the deep GOODS data allow us to study the shape of the relation at low sSFR and low mass. Therefore, we have a broader view of the main sequence and we deal with selection effects. Finally, we parametrize the shape of the sSFR function to fit the data. We try several options for the parametrization. Based on these fits, we can derive accurate measurements of the median sSFR, or of the width of the sSFR function, which modify some previous findings on the star-forming main sequence.

The paper is organized as follows. The data are introduced in Sect. 2. Since the “main sequence” refers only to star-forming galaxies, we need to carefully select this population, as described in Sect. 3. The method used to estimate the sSFR functions and the associated uncertainties is explained in Sect. 4. We discuss the evolution of the sSFR functions in Sect. 5. We compare our reference sSFR functions with the ones obtained using optical SFR tracers in Sect. 6 and with predictions of a semi-analytical model in Sect. 7. Finally, we discuss our results in Sect. 8 and conclude in Sect. 9.

Throughout this paper, we use the standard cosmology (Ω_m = 0.3, Ω_Λ = 0.7 with H₀ = 70 km s^-1 Mpc^-1). Magnitudes are given in the AB system (Oke 1974). The stellar masses (M⋆) are given in units of solar masses (M_⊙) for a Chabrier (2003) initial mass function (IMF). The sSFR is given in Gyr^-1.

2. The galaxy stellar mass and SFR samples

Our analysis combines the data from the GOODS and the COSMOS surveys and our measurements are based on MIPS selected samples at F_{24 μm}> 20 μJy and F_{24 μm}> 60 μJy, respectively.

In the COSMOS field, we use the i⁺-selected catalogue (limiting magnitude of 26.2 mag at 5σ) created by Capak et al. (2007). We use an updated version of the photometric catalogue including the UltraVISTA DR1 data release (McCracken et al. 2012) and new SPLASH IRAC data at 3.6 and 4.5 μm (Capak et al., in prep.). The photometric redshifts are estimated using 30 bands, as described in Ilbert et al. (2013). Their accuracy is similar to Ilbert et al. (2013) in the redshift range considered in this paper (0.2 <z< 1.4). By comparing these photometric redshifts with 10 800 spectroscopic redshifts from the zCOSMOS bright survey (Lilly et al. 2007), we find a precision of σ_{Δz/ (1 + z)} = 0.008 at i⁺< 22.5 and z< 1.4. Using the spectroscopic samples from Comparat et al. (2015), Capak et al. (in prep.) and the VIMOS Ultra-Deep Survey (Le Fèvre et al. 2014), we find σ_{Δz/ (1 + z)} = 0.03 at i⁺< 24.

The stellar masses are estimated using “Le Phare” (Arnouts et al. 2002; Ilbert et al. 2006). We define the stellar mass as the total mass in stars at the considered age (without the mass returned to the interstellar medium by evolved stars). We derive the galaxy stellar masses using a library of synthetic spectra generated using the stellar population synthesis (SPS) model of Bruzual & Charlot (2003). In addition to the library used in Ilbert et al. (2010) assuming exponentially declining SFH, we add two other star formation histories based on delayed SFH (τ^-2te^{− t/τ}) having a maximum SFR peak after 1 and 3 Gyr. For all these templates, two metallicities (solar and half-solar) are considered. Emission lines are added following Ilbert et al. (2009). We include two attenuation curves: the starburst curve of Calzetti et al. (2000) and a curve with a slope λ^0.9 (Appendix A of Arnouts et al. 2013). E(B−V) is allowed to take values as high as 0.7. We assign the mass using the median of the marginalized probability distribution function (PDF). As shown in Mitchell et al. (2013), this procedure allows us to reduce some discontinuities in the mass estimate. The 1σ uncertainties derived from the PDF increase from 0.035 dex at 0.2 <z< 0.4 to 0.055 dex at 1.2 <z< 1.4 for the MIPS selected sample considered in this paper. We also match our own mass estimates with the two independent measurements of the masses published in Brammer et al. (2011) and Muzzin et al. (2013). The three mass catalogues are established for the same sources, but use a different photometry, different photo-z codes and different assumptions to construct the SED templates. Based on this comparison¹, we conclude that the mass uncertainties increase from 0.05 dex at 0.2 <z< 0.4 to 0.07 dex at 1.2 <z< 1.4. Systematic uncertainties on the stellar masses (e.g. due to the IMF or SPS choices) are not included here.

The main part of our analysis is based on a MIPS 24 μm selected catalogue. The deep MIPS S-COSMOS data were taken during Spitzer Cycle 3 and cover the full COSMOS 2-deg² (Sanders et al. 2007). The 24 μm sources are detected with SExtractor (Bertin & Arnouts 1996) and their fluxes measured with a PSF fitting technique (Le Floc’h et al. 2009). Le Floc’h et al. (2009) identified the optical/near-infrared (NIR) counterparts of the 24 μm detection. When possible, we also use the Herschel data observed as part of the PEP survey at 100 and 160 μm (Lutz et al. 2011) and Hermes survey at 250, 350 and 500 μm (Oliver et al. 2012). The Herschel fluxes are extracted using the 24 μm catalogue as prior which makes the cross-identification with the optical sample straightforward.

We estimate the SFR of our sample following exactly the same method as Arnouts et al. (2013). The total SFR is obtained by summing the contribution of the IR and UV light using Eq. (1) of Arnouts et al. (2013), i.e. SFR [M_⊙ yr¹] = 8.6 × 10^-11 (L_IR + 2.3L_NUV). The infrared luminosities are extrapolated from the 24 μm fluxes using the Dale & Helou (2002) library following Le Floc’h et al. (2009). With this method, a given IR luminosity is associated to one template. In order to quantify the uncertainties generated by this extrapolation, we derive using a minimum of three bands (the 24 μm, one band from PACS and one from SPIRE) and we allow any template to be fitted. As shown in the top panel of Fig. 1, we find no systematic offset between and . The dispersion between both measurements increases from 0.03 dex at z< 0.6 to 0.12 at z> 1. Therefore, our extrapolation from the 24 μm flux assuming one SED for a given L_IR does not introduce significant uncertainties or biases. We also compare and ( is computed without using the 24 μm data) in order to use two independent estimates of the L_IR. Based on the scatter of , we expect an uncertainty on the L_IR of 0.06 dex, 0.09 dex, and 0.13 dex at 0.2 <z< 0.6, 0.6 <z< 1.0, and 1 <z< 1.4, respectively. We observe a systematic offset of 0.1 dex between and , showing that one of the two estimates could be biased. As shown in Fig. 1, this offset is present for several sets of templates available in the literature. Such an offset could be partially explained by the combined uncertainties in the absolute calibration of MIPS and/or Herschel data². Still, a systematic shift independent of the redshift does not affect our conclusions. Hereafter, we adopt the Dale & Helou (2002) templates and we use to get the total infrared luminosity.

Since AGN could contaminate the 24 μm emission and bias the stellar mass estimate, we remove the bright X-ray sources detected in XMM (Brusa et al. 2007). We keep the sources identified as IRAC power-laws (Donley et al. 2012), but we also perform the full analysis removing the IRAC power-laws without a noticeable change in our results.

In the GOODS field, we use the FIREWORKS data published by Wuyts et al. (2008). This catalogue reaches K< 24.3 at 5σ over 138 arcmin². We compute the photometric redshifts using Le Phare and the same method as the COSMOS field. We obtain photometric redshifts comparable to the ones of Wuyts et al. (2008) with a precision at σ_{Δz/ (1 + z)} = 0.03 at i⁺< 24. The comparison is based on several spec-z samples compiled by Wuyts et al. (2008), including VVDS data from Le Fèvre et al. (2004) and K20 data from Mignoli et al. (2005). We apply exactly the same method as in COSMOS to derive the stellar masses. This catalogue also includes MIPS data. Following Wuyts et al. (2008), we apply a selection at F_{24 μm}> 20 μJy in this catalogue. We also add the GOODS-Herschel data at 100 and 160 μm (Elbaz et al. 2011). The SFR is estimated following exactly the same method as for the COSMOS field.

3. Selecting the star-forming galaxies

In order to study the evolution of the main sequence, we need to identify star-forming and quiescent galaxies. The presence of a bimodal distribution in color (e.g., Bell et al. 2004; Faber et al. 2007; Franzetti et al. 2007; Smolčić et al. 2008; Fritz et al. 2014) or in the M⋆-SFR plane (e.g., Peng et al. 2010) shows that galaxies are transitioning rapidly from a star-forming main sequence to a red clump. Therefore, a quantitative criterion can be established to select the star-forming galaxies.

Williams et al. (2009) show that the combination of two rest-frame colors (M_U−M_V, M_V−M_J) is sufficient to separate quiescent and star-forming galaxies without mixing galaxies that are red because of dust extinction and the ones with a quenched star formation. We use a modified version of this selection criterion by combining the two rest-frame colors M_NUV−M_R and M_R−M_K following Arnouts et al. (2013). The absolute magnitudes are derived using the method described in Appendix B of Ilbert et al. (2005): in order to minimize the uncertainty induced by the k-correction term, the rest-frame luminosity at a given wavelength λ is derived from the apparent magnitude observed at λ(1 + z). Figure 2 shows the galaxy distribution in the NUV-R-K plane within the COSMOS field. The red clump is clearly isolated from the star-forming sequence by a lower density region in which galaxies transit rapidly. We establish a limit to separate the quiescent and the star-forming galaxies within this lower density region in the NUV-R-K plane. This limit changes with cosmic time because of the evolution of the stellar populations. In order to apply a single criterion at all redshifts to select the star-forming galaxies, we add a time dependent correction C to our selection criterion, with C = −0.17 [t_H(z)−t_H(z = 2)] if z< 2 and t_H the age of the Universe at a given redshift in Gyr. The galaxies with (M_NUV−M_R) + C< 2.6 and (M_NUV−M_R) + C< 2(M_R−M_K) + 1.7 are considered to be star-forming. The considered limit is indicated with the brown dashed lines in Fig. 2. We note that the time correction C is established empirically to produce the cleanest separation between the red and blue regions³. The same criterion is applied to the GOODS sample, providing an equally good separation between the star-forming sequence and the red clump.

We find that 2% of the MIPS sources fall within the quiescent region at 0.2 <z< 1.4 (above the brown dashed lines). This 24 μm flux could be explained by the contribution of an AGN. Moreover, some post-starburst galaxies with a quenched star formation could still be seen in IR (e.g. Hayward et al. 2014). Therefore, we do not include the sources falling in the quiescent locus in our analysis. However, we have checked that including this population would not affect our conclusions.

We emphasize the importance of using a two-color criterion to study the M⋆-SFR relation for the most massive galaxies. We would remove a significant fraction of massive dust-extinguished star-forming galaxies from the main sequence by using a single-color criterion. For instance, we would loose 20% of the galaxies more massive than M⋆ > 10^10.5 M_⊙ with a selection (M_NUV−M_R) > 3.5.

Finally, our color-color selection corresponds approximately to a cut in log (sSFR_SED) at −2 (e.g. Ilbert et al. 2010) with sSFR_SED estimated using the template fitting procedure (blue and red contours in Fig. 2). We check with our dataset that 3–4% of all our star-forming sources at M⋆ > 10^9.5 M_⊙ (not MIPS selected) have log (sSFR_SED) < −2 while 3–4% of the galaxies that we do not classified as star-forming have sSFR_SED> −2. The majority of these sources are located ± 1 dex around log (sSFR) = −2. Therefore, our classification in colors is very similar to a classification in sSFR_SED.

4. Measurement and fit of the sSFR functions

In this section, we describe the method used to derive the sSFR functions per stellar mass bin. We discuss the possible selection effects in the mass-sSFR plan shown in Fig. 3 and correct for them when necessary (mainly the 24 μm flux limit since we are complete in stellar mass). We assume two possible profiles for the sSFR functions (a log-normal function and a double exponential function). Figure 4 shows the sSFR functions, and the best-fit parameters are given in Tables 1 and 2.

4.1. The mass-sSFR scatter diagram

Figure 3 shows the distribution of the sSFR as a function of the stellar mass for star-forming galaxies in the COSMOS field (red crosses) and in the GOODS field (blue triangles). Since GOODS covers a small volume with a deep NIR coverage, this sample includes preferentially low-mass galaxies at z< 1, while COSMOS which covers an area × 30 larger includes rare and massive sources. This difference explains why the GOODS and the COSMOS samples cover a different mass range in Fig. 3. Still, the sSFR values of the COSMOS survey are larger than the values found in the GOODS field for masses M⋆ < 10¹⁰ M_⊙. This effect is explained by the × 3 difference in sensitivity between the two MIPS surveys. While the COSMOS survey includes mostly starbursting sources at low masses, the GOODS survey is able to reach the bulk of the star-forming population. We will discuss in more detail this selection effect in Sect. 7 using a semi-analytical model.

The green solid line corresponds to the relation(1)established using the mass dependency of the sSFR of −0.21 provided by Rodighiero et al. (2011), as well as the normalization of the main sequence at z ~ 2 from the same analysis. We assume an evolution in (1 + z)^2.8 from Sargent et al. (2012). The position of our GOODS data agrees with this relation. However, such parametrization is not suitable for the COSMOS field. It demonstrates the need for a statistical study taking into account selection effects.

4.2. Methodology to estimate the non-parametric sSFR functions

In order to fully characterize the evolution and the shape of the main sequence, we measure the sSFR function, i.e., the number density in a comoving volume (in Mpc^-3) and per logarithmic bin of sSFR (in dex^-1). We derive the sSFR function per stellar mass and redshift bin. We divide the star-forming sample into 6 redshift bins with Δz = 0.2 and four stellar mass bins log (M⋆) = 9.5−10 dex, 10−10.5, 10.5−11, 11−11.5.

We note that the sSFR functions are measured per stellar mass bin of 0.5 dex. Therefore, one has to multiply their normalization by 2 in order to express the density per logarithmic bin of sSFR and per logarithmic bin of M⋆ simultaneously, i.e., in Mpc^-3 dex^-2 (a bivariable galaxy mass and sSFR function).

In order to take into account the flux limit at 24 μm (F_{24 μm}> 20 μJy in GOODS and F_{24 μm}> 60 μJy in COSMOS), we adopt standard estimators as the 1/V_max (Schmidt 1968), the SWML (Efstathiou 1988) and the C⁺ (Lynden-Bell 1971). These estimators are included in the tool ALF used to compute the sSFR function, as described in Appendix B of Ilbert et al. (2005).

Because of the depth of the COSMOS optical and NIR images, we do not need to consider any incompleteness in stellar mass. Indeed, only 4% and 0.5% of the galaxies are fainter than i> 25.5 and m(3.6) > 24 (this magnitude limits are 0.5–1 mag brighter than the magnitude limit of our survey) in the most incomplete bin M⋆ < 10¹⁰ M_⊙ and z> 1.2. Only 2% of the star-forming galaxies would require a 1/V_max correction in this bin for our considered limits in NIR. Therefore, the samples considered in this analysis are complete in mass. Since the GOODS data are deeper than the COSMOS data in optical, the GOODS sample is also complete in stellar mass at M⋆ > 10^9.5 M_⊙ and z< 1.4.

We define sSFR limits, denoted sSFR_complete, above which we can safely correct for selection effects. As shown in Ilbert et al. (2004), if a galaxy population is not observable anymore below a given sSFR, denoted sSFR_complete, the standard estimators cannot correct for this missing population. Moreover, the various estimators are biased differently below sSFR_complete. We adopt the following definition: sSFR_complete is the sSFR for which 90% of the galaxies have their sSFR_limit<sSFR_complete, with sSFR_limit being the lowest sSFR observable for each galaxy given the 24 μm flux limit. Following this procedure, not more than 10% of the galaxies could be missed at sSFR>sSFR_complete. We also restrict our analysis to the sSFR range where the 3 non-parametric estimators produce consistent results. We will use the 1/V_max estimator in this paper, but the results would be the same at sSFR>sSFR_complete using the other estimators. One advantage of the 1/V_max estimator is that it produces a lower limit in density at sSFR<sSFR_complete (Ilbert et al. 2004). Therefore, we conserve this important information and use the 1/V_max estimator as lower limits when possible (the lower limits are shown with arrows in Fig. 4).

Star-forming sources with an extreme dust attenuation could be missing in the parent optical photometric catalogue. Indeed, Le Floc’h et al. (2009) do not find any optical counterpart for 10% of the MIPS sources. Given our survey limits, we establish that these missing sources are redder than m_R−m₂₄> 7.6 (or fake MIPS detections). Dey et al. (2008) and Riguccini et al. (2011) show that >85% galaxies selected with such color criteria are located at z> 1.4, i.e., above the maximum redshift considered in our study.

Finally, the GOODS field (138 arcmin²) covers a smaller area than the COSMOS field (1.5 deg²), the uncertainties due to the cosmic variance reach σ = 0.21−0.37 (for log (M⋆) = 9.75−11.25) at z ~ 0.9 in the GOODS field and σ = 0.09−0.15 in the COSMOS field (Moster et al. 2011). Therefore, we use COSMOS as the anchor for the normalization of the sSFR functions. We adjust the normalization of the sSFR functions in the GOODS field by applying the following factors to the normalization: 1.243, 0.7653, 0.7605, 1.305, 0.9264 and 1.186 in the redshift bins 0.2–0.4, 0.4–0.6, 0.6–0.8, 0.8–1, 1–1.2 and 1.2–1.4, respectively. These factors are derived from the ratio between the redshift distributions of the two fields for a same magnitude selection limit.

4.3. Parametric fit of the sSFR functions

We fit simultaneously the 1/V_max data of the COSMOS and GOODS fields. We consider two possible profiles: a log-normal function and a double exponential function. The log-normal function is parametrized as: (2)with Φ^∗ the normalization factor, sSFR^∗ the characteristic sSFR, and σ the standard deviation. We also consider a double-exponential profile (e.g., Saunders et al. 1990; Le Floc’h et al. 2005), (3)with α the faint-end slope. While the double-exponential profile is not commonly used to describe the sSFR distribution, it allows for a significant density of star-forming galaxies with a low sSFR.

In order to take into account the uncertainties on the sSFR, we convolve these profiles with a Gaussian function having a standard deviation σ = 0.06, 0.07, 0.08, 0.11, 0.14 and 0.17 dex at z = 0.2−0.4, z = 0.4−0.6, z = 0.6−0.8, z = 0.8−1.0, z = 1.0−1.2, and z = 1.2−1.4, respectively. These values are obtained by summing in quadrature the statistical uncertainties expected for the SFR and the stellar masses, as estimated in Sect. 2 (systematic uncertainties are not included). We note that φ_c is the convolved profile.

We add a starburst component to the sSFR function to produce a better fit of the 1/V_max data at high sSFR. Since the contribution of the starbursts cannot be constrained in each individual bin of redshift and stellar mass, we set their contribution following Sargent et al. (2012). We assume that the starbursts are distributed with a log-normal distribution having σ = 0.25 and centered on the median sSFR shifted by + 0.6 dex. We set the starburst contribution to be 3% of the total density of star-forming galaxies.

We fit the 1/V_max data measured in the stellar mass bin [M⋆^min;M_⋆^max] by minimizing the χ² value defined as: (4)with φ_c being the function that we fit, the density at sSFR_i estimated with the 1/V_max estimator, and its associated Poisson errors; N corresponds to the number of bins in sSFR. Equation (4) contains three components:

• the first term of the equation corresponds to the standardχ² minimization method;

• the second term accounts for the lower limits obtained with the 1/V_max estimator below the completeness limit. The negative error indicates that the density is a lower limit;

• the third term includes an additional constraint using the galaxy stellar mass function of the star-forming galaxies, denoted GSMF_SF. Indeed, the sSFR function integrated over the full sSFR range should match the GSMF_SF integrated over the considered mass bin [log (M⋆^min),log (M_⋆^max)] with ΔM = log (M⋆^max)−log (M_⋆^min). Since we slightly modified the method used Ilbert et al. (2013) to compute the stellar masses, and since we are not using the same parent photometric catalogue, we recompute the GSMF_SF in each mass/redshift bin with this sample for the sake of consistency.

The best-fit parameters are given in Tables 1 and 2 for the double-exponential and the log-normal functions, respectively. We also add the values of the median and average sSFR in these tables. We caution that the median and average sSFR are different, even for a log-normal function.

We find that adding the starburst component is necessary in the mass bin 9.5 < log (M⋆) < 10 to reproduce the high sSFR end, while we could ignore it above log (M⋆) > 10.5 (black solid lines and blue dotted lines in Fig. 4).

Despite the combination of GOODS and COSMOS data, we are not able to directly constrain the full shape of the sSFR function. In most of the redshift and mass bins, the sSFR function is incomplete below the peak in sSFR. The lower limits obtained with the 1/V_max estimator at low sSFR (the arrows in Fig. 4) indicate the minimum possible densities, which is important in order to discriminate between a log-normal and a double-exponential profile. In most of the mass and redshift bins, we do not sample sufficiently low sSFR to see any advantage of using either one or the other parametrization. The fit with a double-exponential function is more suitable than the log-normal function at 0.2 <z< 0.6 and 10 < log (M⋆) < 11 (see Fig. 4). In these bins, the lower limits favor a double-exponential parametrization. Moreover, the position of the best-fit sSFR peak is in better agreement with the non-parametric data using a double-exponential profile.

Adding the GSMF_SF information into Eq. (4) brings an important constraint on the sSFR distribution, not visible by a simple examination of the fit. For instance, a higher value of the slope α of the sSFR function in the case of a double-exponential fit will overproduce the density of star-forming galaxies. Even with this constraint, the uncertainties on α remain large and α varies between −1 and 0.5 when left free. Therefore, we arbitrarily set its value at α = −0.5 which is suitable in all the mass/redshift bins.

5. Evolution of the sSFR functions

In this section, we analyze the evolution of the median sSFR derived from the sSFR functions obtained in Sect. 4.3. We obtain a mass-dependent increase in the sSFR as a function of redshift and a decrease in log (sSFR) as −0.17 M⋆. We also combine all our sSFR functions correcting for time evolution, showing that the width of the sSFR distribution increases with the stellar mass.

5.1. Evolution of the median sSFR

Figure 5 shows the evolution of the median sSFR. The median sSFR is obtained from the best-fit functions (see Sect. 4.3) to avoid selection biases. We observe a clear increase in the sSFR as a function of redshift (left panel) and a decrease with M⋆ (right panel).

We adopt the following parametrization of the sSFR evolution as a function of redshift and mass, (5)with a the normalization, β the dependency on the mass and b the dependency on redshift. The best-fit parameters are given in Table 3. Assuming that the redshift evolution of the sSFR does not depend on the mass, we find β = −0.172 ± 0.007 and b = 3.14 ± 0.07. The result is shown with solid lines in Fig. 5. Then, we relax the assumption that the parameter b is independent of the stellar mass and we fit each stellar mass bin independently. We obtain b = 2.88 ± 0.12, b = 3.31 ± 0.10, b = 3.52 ± 0.15 and b = 3.78 ± 0.60 in the stellar mass bins log (M⋆) = 9.5−10 dex, 10−10.5, 10.5−11, and 11−11.5, respectively. The result is shown with dashed lines in Fig. 5 (left). It suggests that the evolution is faster for the massive galaxies, which is in agreement with a downsizing pattern (Cowie et al. 1996). These values are obtained assuming a double-exponential profile, but the results are similar if we consider a log-normal profile.

The parameter b is directly comparable with several values from the literature using the same parametrization in ∝(1 + z)^b. In Karim et al. (2011), b varies between 3.42 and 3.62 at 10.2 < log (M⋆) < 11.1 which is consistent with our results in the same mass range. Elbaz et al. (2011) find an evolution with b = 2.8, based on deep GOODS data. Therefore, a dependency of b on the stellar mass could explain the differences between the various values found in the literature.

With our parametrization, log (sSFR) is proportional to M⋆. The parameter β that we obtain cannot be directly compared with values from the literature. In most of the studies, a linear dependency with log (M⋆) is assumed (dotted lines in Fig. 5, right panel) while we assume a linear dependency with M⋆ (solid lines). Our parametrization in M⋆ produces a more rapid decrease in the sSFR at M⋆ > 10^10.5 M_⊙ and less evolution at lower mass than a parametrization in log (M⋆). As shown in Fig. 5, a parametrization in log (M⋆) is not suitable for the considered mass range. It could explain why the slope values in the literature depends on the considered mass range when log (sSFR) ∝ log (M_⋆) (Lee et al. 2015; Whitaker et al. 2014). We emphasize that the choice of having log (sSFR) proportional to M⋆ is physically motivated since “τ models” converge toward such parametrization at high masses (e.g. Noeske et al. 2007b; Bauer et al. 2013).

Finally, we would like to note that including the GSMF_SF within the χ² expression (Eq. (4)) brings a decisive constraint over the median sSFR. We illustrate this point in Fig. 6, showing the evolution of the median sSFR with redshift in the mass bin 9.5 < log (M⋆) < 10. In this mass bin, the non-parametric sSFR function is primarily composed of lower limits below the sSFR peak. By using only the information contained in the non-parametric data to perform the fit, we would not constrain the median sSFR, as shown by the green error bars. Adding the GSMF_SF as an additional observable within the χ² expression breaks the degeneracy between the best-fit parameters and allows an accurate estimate of the median sSFR (red and blue error bars).

5.2. Broadening of the sSFR function

Figure 7 shows the evolution of σ as a function of redshift, in the case of a log-normal fit. We find that σ increases with mass. We also find that σ is consistent with being constant with redshift at z< 1.4.

Since log (sSFR) is not linearly proportional to log (M⋆), it could create an artificial broadening of the sSFR function, especially if a large redshift range is considered. In order to test this effect, we compute the sSFR functions in smaller mass bins of Δ(log M⋆) = 0.2 rather than Δ(log M⋆) = 0.5. We still find that σ is close to 0.3 at 10.0 < log (M⋆) < 10.2 and 0.45 at 11.0 < log (M) < 11.2.

We use the σ value from the log-normal function since this value can be directly compared with other values from the literature. Our value of σ is higher than previous studies which converge to an rms of 0.2 dex (e.g. Peng et al. 2010; Sargent et al. 2012; Salmi et al. 2013; Speagle et al. 2014). Noeske et al. (2007a) find an rms of 0.35 dex, before deconvolution, which is close to our value for the intermediate mass range. Almost no study finds a scatter of >0.4 dex as we get for M⋆ > 10¹¹ M_⊙ galaxies, except Salim et al. (2007).

An attractive interpretation is that the different mass ranges covered by each survey could explain the various rms measured in the literature. However, Salim et al. (2007) and Whitaker et al. (2012) find that the scatter of the main sequence decreases with M⋆, which is at odds with our result. Moreover, Lee et al. (2015) find an rms of 0.35 dex almost constant with the M⋆ using similar data.

If we measure the rms of our own M⋆−sSFR scatter diagram (Fig. 3), we obtain an rms below 0.25 dex at log (M⋆) = 11−11.5, much lower than the σ measured using the sSFR function. One interpretation is that the dynamical sSFR range covered by the data is not sufficiently large to correctly estimate the r.m.s. from a scatter diagram. We demonstrate this effect with a simulated catalogue in Sect. 7. The advantage of using the sSFR distribution is to extrapolate the shape of the function over the full sSFR range, even if the data span a restricted sSFR range.

5.3. Shape of the combined sSFR functions

We correct for the redshift evolution of the sSFR and we combine all the measurements at z = 0. Figure 8 shows the combined 1/V_max data per stellar mass bin. The shape of the sSFR distribution appears invariant with time: the dispersion between the data points is around 0.1 dex as shown in the insets of Fig. 8. We observe that the double-exponential and the log-normal profiles provide a good fit of the combined data (solid red line and blue dashed line, respectively). Even so, the χ² values are smaller for double-exponential fit at 10 < log (M⋆) < 11.

As shown in Fig. 7, the broadening of the sSFR function is also visible in the combined sSFR functions. If we let the fraction of starbursts⁴ vary while we fit the combined data, we obtain that the fraction of starbursts is consistent with 0 at log (M⋆) > 10.5, but the uncertainties are consistent with a contribution of 1%. At lower mass, we find a fraction of starbursts of 2 ± 1% and 4 ± 1% assuming a double-exponential profile at log (M⋆ /M_⊙) = 10−10.5 and 9−9.5. The associated uncertainties are underestimated since all the parameters describing the shape of the starburst contribution are set, except the normalization. There is a hint that the fraction of starburst increases at low masses. However, we would need a survey covering a larger volume to cover the high sSFR range of the distribution, since massive starbursting galaxies are rare galaxies.

6. The sSFR function using other SFR tracers

When using the SFR_UV+IR tracer, the density of star-forming galaxies below the sSFR peak relies on the extrapolation of the best-fit profile (double-exponential or log-normal). The UV and optical SFR tracers allow us to cover the full sSFR range and could bring some information at low sSFR. In this section, we derive the sSFR function using SFR tracers based on the stellar emissivity only (without using the IR data), as shown in Fig. 9.

We estimate the SFR from the SED fitting procedure (the same method as for the stellar mass, see Sect. 2). When comparing SFR_SED and SFR_UV+IR, we find some bias reaching 0.25 dex and a scatter between 0.25 and 0.35 dex (the scatter increases both with the mass and the redshift). The comparison between our reference sSFR_UV+IR functions (black solid lines) and the sSFR_SED functions (red solid lines) shows that the sSFR_SED functions are much flatter⁵. We convolve the reference sSFR_UV+IR functions by a Gaussian function with σ = 0.3 dex (black dotted lines) to mimic the expected uncertainties. After this convolution, the agreement between the convolved sSFR_UV+IR function (dotted lines) and the sSFR_SED function is better. The sSFR_SED function would favor a fit with a double-exponential profile. At high mass M⋆ > 10^10.5 M_⊙, the density of low sSFR_SED galaxies even exceeds what we expect from the double-exponential extrapolation.

Arnouts et al. (2013) have developed a new method for estimating the SFR from optical data. This SFR – denoted SFR_NRK – is estimated by correcting the UV intrinsic luminosity ℒ_UV by the infrared excess IRX = ℒ_IR/ ℒ_UV, directly estimated from the position of the galaxy into the NUV-R-K plane. We use the parametrization of the IRX from Arnouts et al. (2013), slightly modified by Le Floc’h et al. (in prep.)⁶. The dispersion between SFR_NRK and SFR_IR+UV is around 0.15 dex (could reach 0.2 dex). Figure 9 shows the excellent agreement between the sSFR_NRK functions (blue open stars) and the sSFR_IR+UV functions (black solid lines) at M⋆ < 10¹¹ M_⊙. The positions of the peak of the sSFR functions are similar (within 0.2 dex). Using the sSFR_NRK tracer, we find a better agreement with the log-normal profile since the sSFR_NRK density falls sharply in the lowest sSFR bin. Still, the double-exponential profile is more appropriate in some bins (for instance, in the redshift bin 0.2 <z< 0.4). Given a possible bias in sSFR_NRK at low sSFR (Arnouts et al. 2013) leading to an underestimation of the sSFR_NRK, we cannot conclude that this sharp cut at low sSFR_NRK is real.

To summarize, given the large uncertainties affecting the UV and optical SFR tracers, it is still challenging to constrain the low sSFR end.

Finally, we note that our reference SFR tracer combining UV and IR could be overestimated, since dust could be heated by the old stellar populations (Utomo et al. 2014). This bias would mainly affect galaxies with sSFR [Gyr^-1] < −1. Correcting for such a bias would modify the shape of our sSFR function: the slope of the low sSFR side of the sSFR function obtained with the double-exponential profile would be even flatter than the observed one.

7. Comparison with a semi-analytical model

We now compare our results with the predictions of a semi-analytical model. The mock catalogue is based on ΛCDM simulations from Wang et al. (2008) and the galaxy properties were generated using the galaxy formation model, as detailed in De Lucia & Blaizot (2007) and Wang & White (2008). The light cone survey covers an area of 1.4 × 1.4 deg² similar to COSMOS. The redshift, the SFR and the stellar mass are available for all galaxies in the simulation, as well as the observed magnitudes expected for these galaxies. We select the star-forming galaxies using a criterion similar to our selection in the NUV-R-K plane.

We first test if we can reproduce the same selection effects as discussed in Sect. 4.1. We apply a K-band selection in the simulation similar to the ones applied in the data (K< 24 and K< 24.3 for the COSMOS and GOODS surveys, respectively). The selection at 24 μm creates an observational limit in the redshift-SFR plane. We apply the same SFR limits in the simulation as the ones established for the COSMOS and the GOODS surveys. Finally, we select galaxies over an area of 1.5 deg² for COSMOS and 138 arcmin² for GOODS. The blue and red points in Fig. 11 show the distribution of the simulated sources in the M⋆−sSFR plane for the GOODS-like and COSMOS-like surveys, respectively. We reproduce exactly the same selection effect as the ones discussed in Sect. 4.1. The predicted COSMOS-like and GOODS-like surveys do not cover the same area in the M⋆−sSFR plane. Even with the GOODS-like survey, the MIPS data are not sufficiently deep to provide a representative sample of low-mass galaxies in terms of sSFR.

We also test that the width of the main sequence is not correctly measured using simply the rms of the sample without taking into account selection effects. For instance, the rms of the sSFR without any selection (orange points) is 0.38 dex at 0.8 <z< 1 and 9.5 < log (M⋆) < 10, but we only measure an rms of 0.18 dex and 0.23 dex in the COSMOS-like and GOODS-like survey, respectively. It illustrates the necessity of taking into account selection effects in SFR limited surveys, as discussed in Rodighiero et al. (2014) and Kelson (2014). In particular, any study looking at the evolution of the sSFR with the mass would be biased.

Finally, we directly compute the predicted sSFR functions from the simulated catalogue. A comparison with the sSFR functions predicted by the models (green lines) and the observed ones is shown in Fig. 10. A qualitative comparison shows that the predicted shape of the sSFR functions is similar to the observed one. A parametrization with a double-exponential profile is perfectly suitable for the simulation. In specific redshift and mass bins, the agreement with the data is remarkable (e.g. 0.4 <z< 0.6 and log (M⋆) < 10.5). The slope of the predicted sSFR function is in excellent agreement with the double-exponential profile. The predicted density of low sSFR star-forming galaxies exceeds the density extrapolated from the log-normal profile. Therefore, an extrapolation with the double-exponential profile is more natural than a log-normal profile on the theoretical point of view.

The agreement between the predicted and observed sSFR functions breaks down for galaxies more massive than log (M⋆) > 10.5, but also at z> 1. As a global trend, the galaxies with the largest sSFR are missed in the simulation (e.g. 0.8 <z< 1 and log (M⋆) > 10.5). At z> 1, the predicted distribution is shifted at lower sSFR in comparison to the data. We will discuss in Sect. 8.1 the predicted evolution of the median sSFR compared to the observed one.

Finally, we note that we do not use the most recent SAMs. We keep the same SAM as in our previous works in COSMOS. For instance, we use the same SAM to compare predicted and observed GSMF_SF as in Ilbert et al. (2013). However, more detailed comparisons with recent numerical simulations will be necessary in the future.

8. Discussion

In this section, we discuss our main results: 1) the mass-dependent increase in the sSFR with redshift; 2) the decrease in log (sSFR) as −0.17 M⋆; and 3) the broadening of the sSFR function with mass. We discuss here the numerous complex processes that can reduce the sSFR as the stellar mass increases, from the hot halo quenching mode to secular evolution of galaxy disks. The diversity of these processes could explain the broadening of the sSFR functions with mass, and their complexity could reduce the ability of the SAM to reproduce the sSFR evolution for the most massive galaxies.

8.1. Increase in sSFR with redshift: link with the cosmological accretion rate

We compare here the evolution of the median sSFR with the specific mass increase rate sMIR_DM (Ṁ_H/M_H following Lilly et al. 2013) and with the predictions of semi-analytical models.

Assuming the gas inflow rate is driven by the cosmological accretion rate of the dark matter structures, we expect that the sSFR follows the evolution of the sMIR_DM (in the following, we implicitly divide the sMIR_DM by 1−R with R the return fraction, as discussed in Appendix A). In simple models in which galaxies reach a quasi-steady state (Bouché et al. 2010; Lilly et al. 2013), the evolution of the sSFR is coupled with the evolution of the sMIR_DM. Based on N-body simulations and extended Press-Schechter formalism, Neistein & Dekel (2008) show that sMIR_DM evolves as ∝0.047(M_H/ 10¹² M_⊙)^0.15 × (1 + z + 0.1(1 + z)^-1.25)^2.5, which could explain why the sSFR increases with redshift. The green shaded areas in Figs. 12 and 13 show the evolution of the sMIR_DM, after having determined the value of M_H using the stellar-to-halo mass ratio from Coupon et al. (2015).

We first discuss the sample of low-mass galaxies at 9.5 < log (M⋆) < 10. In Fig. 12, we show the evolution of sMIR_DM and we add the sSFR evolution predicted by the SAM from Weinmann et al. (2011; red solid line) and Wang et al. (2008; brown dashed line). In this mass range, the evolution of the sSFR predicted by the SAM follows closely the evolution of sMIR_DM (Weinmann et al. 2011). This statement is also true even with the latest results from the hydrodynamical Illustris simulation (Sparre et al. 2015). We also show in Fig. 12 the data compilations from Weinmann et al. (2011) and from Behroozi et al. (2013) from various measurements available in the literature (gray and magenta shaded areas). As discussed by Weinmann et al. (2011), the observed sSFR from the literature are well above the predictions of the SAM at z< 1.5. We add in Fig. 12 our own measurements of the median sSFR. Our values are located in the lowest part of the Weinmann et al. (2011) and from Behroozi et al. (2013) compilations. Therefore, we find a much better agreement between the observed and theoretical evolution of the sSFR, as expected if the gas feeding is directly driven by the cosmological accretion rate. There are several reasons for the difference with previous results: 1) we take into account selection effects that lead to a lower median sSFR value than the ones obtained directly from a SFR limited survey; 2) the previous compilations do not differentiate between median and average sSFR which could modify the sSFR values by 0.2 dex; or 3) a systematic uncertainty of −0.1 dex could affect our SFR measurements as discussed in Sect. 2. Error bars in Fig. 12 include these systematic uncertainties, as well as a possible ±0.1 dex systematic uncertainty on the stellar mass.

While the sSFR evolution matches the SAM predictions and follows the sMIR_DM evolution for our low-mass sample, this agreement breaks down at higher masses. Figure 13 shows the evolution of the median sSFR predicted by the Wang et al. (2008) model as well as the evolution of the sMIR_DM in several stellar mass bins. We first note that the evolution of sMIR_DM no longer corresponds to the evolution of the sSFR in the SAM. Indeed, AGN feedback is included in the SAM in order to quench the star formation in massive halos (e.g. Croton et al. 2006; Cattaneo et al. 2006; Wang et al. 2008). While these recipes are sufficient to recover a broad agreement with the observed sSFR, we obtain that the median sSFR evolves faster in our data than in the SAM of Wang et al. (2008). In the data, b varies from 2.9 to 3.8 from low-mass to high-mass galaxies. We find the reverse trend in the simulation. The simulation predicts that b decreases with mass: b = 2.3, 2.1, 1.9, 1.5 at log (M⋆) = 9.5−10 dex, 10−10.5, 10.5−11 and 11−11.5, respectively. We also observe that the width of the sSFR function is smaller in the model than in the data for the massive galaxies. The simulated scatter of the sSFR distribution is 0.22 dex at M⋆ < 10¹¹ M_⊙ but reaches 0.16 dex for the most massive galaxies. Therefore, the trend is the reverse of the observed one.

In Ilbert et al. (2013), the low-mass end of star-forming GSMF_SF is accurately reproduced by the SAM model of Wang et al. (2008) while the model under-predicts the density of massive star-forming galaxies (their Fig. 14). Here, we also show that the evolution of the sSFR with redshift is in agreement with the evolution predicted by the SAMs for low-mass galaxies, but complex physical processes that could affect the SFH in massive galaxies, such as quenching or secular evolution, need to be modeled more accurately. In particular, galaxies with the highest sSFR are missing in the simulation at z ~ 1 as shown in Sect. 7.

8.2. Gradual decline of the sSFR with the mass: quenching processes and/or lower efficiency of the star formation

One of our main results is that the full sSFR distribution is shifted toward lower sSFR as the mass increases, with log (sSFR) ∝ −0.17 M_⋆. We discuss here possible mechanisms that could create such a dependency on the stellar mass.

8.2.1. Quenching processes

A first hypothesis is that all massive galaxies are on their way to being quenched and we observe galaxies transitioning toward an even lower sSFR.

In some scenarios, such as Hopkins et al. (2008), a major merger could trigger a burst of star formation and then quench a galaxy in less than 0.3 Gyr. This quenching process cannot be ongoing for all massive star forming galaxies simultaneously: we would observe the density of massive star-forming galaxies dropping rapidly with time, while the high-mass end of the GSMF_SF does not evolve significantly at z< 1 (e.g. Arnouts et al. 2007; Ilbert et al. 2010; Boissier et al. 2010). With such a short quenching timescale, star-forming galaxies would be removed almost instantaneously from our considered sample.

Galaxy could also be quenched by an exhaustion of the cold gas supply as the DM halos grow. For instance, cold accretion across filaments is suppressed in massive halos at z< 2 (Dekel et al. 2009) which reduces the supply of cold gas and then in the star formation. Hydrodynamical simulations predict the formation of a virial shock in dark matter halos with M_H> 10¹² M_⊙. These massive halos can be maintained “hot” with the radio AGN feedback mode or extreme star formation feedback (e.g., Croton et al. 2006; Cattaneo et al. 2006; Wang et al. 2008). According to Gabor & Davé (priv. comm.), mass quenching and environment quenching would be the consequence of the same process: the starving of the galaxies falling in a halo more massive than 10¹² M_⊙. The simple model of Noeske et al. (2007b) reflects the SFH in such gas exhaustion case. In Noeske et al. (2007b), the decreases in the sSFR with mass is reproduced by assuming exponentially declining SFH with τ having an inverse dependency with mass – τ ∝ 1/ℳ. The “stage” model of Noeske et al. (2007b) reproduces well the turn-over at high mass that we observe. For a galaxy as massive as log (M⋆) = 11.3 at z = 0.5, this model associates an exponentially declining SFH with a τ value as large as 4 Gyr. Therefore, the bending of the sSFR with mass could be explained by gas exhaustion over long timescales >3–4 Gyr (τ value decreases with stellar mass). Such timescales are longer than the one usually adopted in simulations to quench star formation in hot halos (typically 1.2 ± 0.5 Gyr for Gabor & Davé, priv. comm.). Therefore, if the quenching in the hot halo mode explains the bending of the sSFR at high mass, it should occur on a longer timescale than usually assumed in simulation. However, if the quenching in the hot halo mode occurs with a timescale <1 Gyr, this process does not explain the bending of the sSFR with mass (although this process could be perfectly relevant to explain the formation of the quenched galaxies over time).

While these quenching processes are probably crucial to generating the quiescent population, they do not likely explain the evolution of the sSFR with mass in the star-forming population since they act on too short a timescale.

8.2.2. Declining efficiency of the star formation toward massive systems – impact of the bulge

A second possibility is that the efficiency in forming new stellar populations declines slowly as the stellar mass increases, without necessarily quenching. Kassin et al. (2012) show that the massive galaxies are on average more kinetically settled at 0.2 <z< 1.2. They find a possible threshold around 10^10.4 M_⊙ to move from a disordered to settled disk. If we speculate that random motion in the gas is going in the direction of a higher star formation efficiency, it could explain a decrease in the sSFR above 10^10.4 M_⊙. Sheth et al. (2008) show that most massive spiral galaxies have a higher fraction of bars, associated with a bulge and having redder colors. They suggest that massive systems are more mature.

Abramson et al. (2014) show in the SDSS that the decrease in the sSFR with stellar mass is explained by the increase in the bulge-mass fraction with stellar mass. The bulge is less efficient in forming stars, which explains a decrease in sSFR with mass. At 0.5 <z< 2.5, Lang et al. (2014) show that the mass fraction within the bulge increases from 30% in 10¹⁰ M_⊙ star-forming galaxies to 50% in 10¹¹ M_⊙ star-forming galaxies. Surprisingly, these ratios remain consistent between z ~ 1 and z ~ 2. Therefore, the mass contribution of the bulge to the total mass increases with M⋆ at all redshifts.

The presence of the bulge could lower the star formation efficiency. In the local Universe, Saintonge et al. (2012) show that the depletion timescale of molecular gas is longer when the galaxy is bulge-dominated, pointing to a lower star formation efficiency when a bulge is present. Using hydrodynamical simulation, Martig et al. (2009) show that a bulge stabilizes the disk against fragmentation and this process suppresses the formation of massive star-forming clumps in the inner part of the galaxy. Genzel et al. (2014) show that the Toomre parameter Q increases at the galaxy center for a sample of z ~ 2 massive disk galaxies, shutting off the gravitational instability and reducing the star formation efficiency in the inner part of the disk. Finally, Förster-Schreiber et al. (2014) show that the presence of AGN-driven massive outflows in the nuclear region that are visible only for their most massive disk galaxies at z ~ 2 (M⋆ > 10¹¹ M_⊙). Such outflows could clear the inner region from the gas and suppress the star formation in the bulge.

We note that the decline (or even the shut down) of the star formation in the inner region of the galaxy does not imply a quenching of the star formation in the entire galaxy. We take as an example the case of the Milky Way (MW). Snaith et al. (2014) and van Dokkum et al. (2013) analyze the SFH of the MW. For a lookback time of 6 Gyr, which corresponds to z ~ 0.7, these studies expect a SFR below 3 M_⊙/yr and log (M⋆) = 10.6. From Fig. 4, the MW falls in the declining part of the sSFR function with a log (sSFR [Gyr^-1] ) = −1.1 dex. Six Gyr later, the MW is not yet a quiescent galaxy and still forming 1 M_⊙/yr (e.g. van Dokkum et al. 2013). There is no reason why the MW should quench on a timescale of a few Gyr. Therefore, these massive galaxies with a low sSFR are not necessarily quenching but could simply be quietly forming stars along cosmic time, as in the MW. A significant density of low sSFR star-forming galaxies is expected in the SAM (see Sect. 7), in agreement with a double-exponential profile for the sSFR function. Unfortunately, the small dynamical sSFR range covered by our SFR_UV+IR tracer in COSMOS and GOODS (see Sect. 5.3), as well as the large uncertainties within the SFR_SED and SFR_NRK tracers (see Sect. 6) do not allow us to definitively come to a conclusion about the density of the low sSFR galaxies not yet quenched.

Finally, we note that the bulge formation could be done through two channels, through secular evolution and by major and minor mergers. In the former case, a bulge could form along time under the action of a bar (e.g. Perez et al. 2013), or through gravitational disk instabilities with large star forming clumps moving inward (e.g. Immeli et al. 2004; Bournaud et al. 2007; Genzel et al. 2008; Bournaud et al. 2011; Perez et al. 2013). Our analysis provides useful information on the SFH of the galaxy which evolves secularly. However, if the bulge originates from galaxy mergers (e.g. Kauffmann & Haehnelt 2000; Martig et al. 2009), the stellar mass has not been formed in situ which makes the ratio SFR/ (M_bulge + M_disk) difficult to interpret in term of SFH (Abramson et al. 2014). In general, mergers would bring stellar mass created ex-situ, leading to a stellar mass growth. If the SFR is not triggered to a higher value during the merger, it could lead to a growth in mass and then a decrease in the sSFR (Peng et al. 2014).

8.3. Broadening of the M_⋆ – SFR relation: star formation stochasticity and diversity in SFHs

As discussed in Sect. 5.2 and as shown in Fig. 7, we find that the intrinsic scatter of the main sequence increases with mass. In particular, the standard deviation found for the most massive galaxies 11 < log (M⋆) < 11.5 reaches σ ~ 0.45, which is well above the values commonly found in the literature. In Appendix B, we show that the intrinsic sSFR evolution and the criterion used to select star-forming galaxies do not artificially create a broadening of the sSFR function with the mass.

The intrinsic scatter of the M⋆-SFR relation indicates how tightly the instantaneous star formation is determined by the past star formation history of the galaxies. Numerous processes could scatter the relation: the intrinsic scatter of the sMIR, galaxy mergers, the variety of the possible SFHs, or the variation of the star formation efficiency within the galaxy itself.

The dynamics of the gas and star content within a galaxy could create SFR variations over million-year timescale. These variations create a natural scatter around the M⋆-SFR relation (Domínguez et al. 2014). Hopkins et al. (2014) analyze the variability of eight galaxies using hydrodynamical simulations and show that the variability could easily reach 50% for M⋆^∗ galaxies when considering a timescale of 20 millions years. In particular, SN feedback has an important impact on this rapid variation. These stochastic fluctuations result from variations in the star formation efficiency over short timescales, generated mainly by the local impact of SN feedback. Based on hydrodynamical simulations, Domínguez et al. (2014) show that these fluctuations generate an intrinsic scatter in the M⋆-SFR relation reaching 0.5 dex for the dwarf galaxies at M⋆ ~ 10⁷ M_⊙ which decreases at 0.2 dex for intermediate mass galaxies at M⋆ ~ 10¹⁰ M_⊙. Since our lowest mass range is 9.5 < log (M⋆) < 10, we cannot detect such a decrease of the scatter with the mass. Still, the intrinsic scatter that we measure for our less massive galaxies 10^9.5<M⋆ /M_⊙< 10¹⁰ could be explained by the stochasticity of the star formation.

The intrinsic scatter in the M⋆-SFR relation induced by the star formation stochasticity in individual galaxies decreases with the stellar mass (Domínguez et al. 2014). Therefore, this process does not explain the increase in σ that we find at M⋆ > 10¹⁰ M_⊙. We also do not expect the scatter of the sMIR to increase with the halo mass. Indeed, we do not detect an increase in the sSFR scatter in the SAM (see Sect. 8.1). Another possibility is that the diversity in the possible SFH increases with the mass. As discussed in Sect. 8.2, numerous processes could affect the SFH and tend to reduce the sSFR as the mass increases. Given the variety of these processes, and their possible dependency on the halo mass and on the galaxy morphology (e.g. the growth of a bulge), the impact on the SFH could vary significantly. Therefore, the same processes could be simultaneously responsible for the increasing diversity in the SFHs (i.e., the scatter of the relation) and for the decrease in the sSFR with the mass.

Kelson (2014) defines a statistical framework using the central limit theorem to predict the width of the M⋆-SFR relation. In this paper, the Hurst parameter H (Kelson 2014, and references therein) determines the behavior of the stochastic fluctuations in the SFR. For H = 0, there is no stochastic fluctuation. For a value of H = 0.5, there is no covariance between the stochastic changes in SFR, i.e., the expectation of the SFR at t_{i + 1} has the same value as the SFR at t_i. Using the central limit theorem, Kelson (2014) predicts a width of 0.3 dex for the log (sSFR) distribution when H = 0.5. For a value of H = 1, i.e., the stochastic changes in SFR are strongly correlated with previous values (if the SFR decreases at a given timestep, the SFR is more likely to decrease again the following timestep). For a value of H = 1, Kelson (2014) find that the width of the distribution reaches 0.43 dex. If H changes with the mass, it could explain why we observe a variation of σ from 0.3 to 0.45 dex between our two extreme considered bins. It implies that stochastic changes in SFR are more correlated as the mass increases, which could be seen as a larger diversity of SFH as the mass increases. Kelson (2014) shows that the difference between the median and the averaged sSFR could be used to establish the value of H. We find that the differences between the median and the average sSFR are around 0.1–0.15 dex in the lowest mass bin while it increases at 0.2–0.25 dex in the highest mass bin (see Tables 1 and 2). Therefore, we measure that H increases with stellar mass.

8.4. Interpreting the evolution of the SFR function and infrared luminosity functions

In this section, we convert the sSFR functions into SFR functions and we discuss the evolution of the IR luminosity function based on our results. This approach is complementary to Sargent et al. (2012) and Bernhard et al. (2014), who combined the GSMF_SF and an universal sSFR-distribution based on Rodighiero (2011) to interpret the evolution of the SFR function.

The SFR function in a given mass bin is easily obtained by simply adding the median of log (M⋆) of the considered bin⁷ to the sSFR function. We sum the SFR functions computed in several stellar mass bins to obtain the total SFR function, as shown with a thin black line in Fig. 14.

Still, our data are limited to M⋆ > 10^9.5 M_⊙ and we need to account for the contribution of the low-mass population when we derive the global SFR function (thick black line in Fig. 14). Therefore, we assume that:

• the density of star-forming galaxies (in log) evolvesproportionally to −0.4log (M⋆), as derived from the GSMF_SF of star-forming galaxies (e.g. Peng et al. 2010; Baldry et al. 2012; Ilbert et al. 2013; Tomczak et al. 2014)⁸;

• the shape of the sSFR function at log (M⋆) = 9.5−10 is conserved at lower mass. The width of the main sequence found in our lowest mass bin is similar to the one found in the deepest surveys (e.g. Whitaker et al. 2012);

• our parametrization of the median sSFR evolution holds at M⋆ < 10^9.5 M_⊙. With our parametrization, log (sSFR) does not depend significantly on M⋆ at M⋆ < 10^9.5 M_⊙, in agreement with Lee et al. (2015) and Whitaker et al. (2014).

We reconstruct the total SFR functions, as shown with the thick solid lines in Fig. 14.

In Fig. 14, we compare our total SFR function with direct measurements from the literature. We convert the IR luminosity functions into SFR functions following Kennicutt (1998). We find an excellent agreement between our SFR functions and the ones derived directly from the IR luminosity functions.

In the remaining part of this section, we use our results on the sSFR to interpret the behavior of the IR luminosity function discussed in the literature.

The total SFR function could be characterized as the combination of two power laws (e.g. Magnelli et al. 2009), with a change in slope at SFR_knee. Assuming that the slope and the characteristic M⋆^∗ of the star-forming GSMF_SF do not evolve with time which is reasonable at z< 1.4 (e.g. Arnouts et al. 2007; Ilbert et al. 2010; Boissier et al. 2010; Peng et al. 2010; Moustakas et al. 2013), the SFR_knee should evolve as the sSFR. Using the sSFR evolution in 3.2(1 + z) found in Sect. 5.1, we indeed reproduce the position of the knee in Fig. 14 (vertical red thick line at top of each panel). By fitting the evolution of the knee of IR LFs, Le Floc’h et al. (2005), Magnelli et al. (2011, 2013) find a consistent evolution of , , and , respectively.

Figure 15 shows the contribution of the galaxies of a given stellar mass range to the full star-forming population above a given SFR threshold. If we use SFR_knee as the SFR threshold, we obtain that the contribution of a given stellar mass range remains stationary over the full redshift range (solid lines). Figure 15 also shows that galaxies at M⋆ > 10¹⁰ M_⊙ dominate the SFR function above SFR_knee, while galaxies with M⋆ < 10¹⁰ M_⊙ contribute to less than 5%.

Several studies tried to reconcile the fact that the GSMF_SF is accuretely represented by a Schechter function while the IR luminosity function is better represented by a double-exponential (e.g. Sargent et al. 2012). We propose here the following interpretation: the SFR function could be seen as the GSMF_SF convolved with the sSFR function. The density of star-forming galaxies drops above M⋆^∗ in the GSMF_SF, and the contribution of M⋆ > 10¹¹ M_⊙ galaxies stays below 20%. The high star-forming end is dominated by galaxies around M⋆^∗. The SFR of the galaxies around M⋆^∗ will be spread following their distribution in sSFR. Therefore, the shape of the SFR function above SFR^∗ is driven by the width of the sSFR function of M⋆^∗ galaxies.

Finally, the faint-end slope of the SFR functions should also be a power-law with the same slope as the star-forming GSMF_SF if the sSFR does not depend on the mass at M⋆ < 10^9.5 M_⊙ (the term depending on the mass in Eq. (5) becomes negligible). Gruppioni et al. (2013) find a slope of −0.2 and Magnelli et al. (2009) a slope of −0.6 for the IR luminosity functions, while we would expect a value around −0.4 from the GSMF_SF (e.g., Peng et al. 2010; Baldry et al. 2012; Ilbert et al. 2013; Tomczak et al. 2014).

9. Conclusions

We characterize the shape and the evolution of the star-forming main sequence by measuring the sSFR function, i.e., the number density in a comoving volume (in Mpc^-3) and per logarithmic bin of sSFR (in dex^-1) of star-forming galaxies. We combine the data from the GOODS and the COSMOS surveys and we derive the sSFR functions at 0.2 <z< 1.4 in four stellar mass bins between 10^9.5 M_⊙<M_⋆< 10^11.5 M_⊙. We show that the GOODS and the COSMOS surveys do not cover the same area in the M⋆-sSFR plane, which demonstrates the importance of taking into account selection effects in the study of the main sequence.

We base our analysis on a MIPS 24 μm selected catalogue, adding Herschel data when possible. We estimate the SFR by summing the contribution of the IR and UV light. While our conclusions are based on the sSFR_UV+IR functions, we also measure the sSFR functions using optically based tracers of the SFR.

We estimate the sSFR functions of star-forming galaxies using several non-parametric estimators. We select the star-forming population using the presence of a bimodal distribution in the M_NUV−M_R/M_R−M_K plane and we check that our conclusions are not too sensitive to the exact position of the selection criterion. We fit the non-parametric sSFR functions by considering two possible profiles: a log-normal function and a double-exponential function. We add a starburst component to the sSFR function and we also add an additional constraint in the fitting procedure using the GSMF_SF.

Based on our sSFR functions, we derive the evolution of the median and average sSFR. We obtain a clear increase in the sSFR as a function of redshift as ∝(1 + z)^b. Assuming that the sSFR evolution does not depend on the mass, we find b = 3.18 ± 0.06. If we allow b to depend on the mass, we obtain that the evolution is faster for massive galaxies: b varies from b = 2.88 ± 0.12 at M⋆ = 10^9.5 M_⊙ to b = 3.78 ± 0.60 at M⋆ = 10^11.5 M_⊙. Our observed evolution of the sSFR is consistent with the evolution of the sMIR_DM for M⋆ < 10¹⁰ M_⊙ galaxies, but deviates from it at higher masses.

We also compare our results with the predictions of a semi-analytical model from Wang et al. (2008). While the predicted sSFR functions could be parametrized by a double-exponential profile and matches our results at M⋆ < 10^10.5 M_⊙ reasonably well, we observe that the agreement breaks down for massive galaxies at high sSFR. The description of the recipes impacting the SFH of massive galaxies should probably be improved in this SAM.

We note that even at z< 1, it is challenging to constrain the full shape of the sSFR functions. Dust-free tracers of the SFR do not reach a sufficiently deep SFR limit to sample well below the peak in sSFR, while tracers based on the optical are prone to large biases because of uncertain dust corrections. Still, we combine all non-parametric estimates of the sSFR functions at z = 0. We find that the shape of the sSFR distribution seems invariant with time at z< 1.4 but depends on the mass. We observe a broadening of the main sequence with M⋆. Assuming a log-normal distribution, we find that σ does not vary with redshift at z< 1.4, and increases from 0.28 to 0.46 dex between log (M⋆) = 9.5−10 and log (M⋆) = 11−11.5 dex. While the stochasticity of the star formation in individual galaxies could explain the width of the sSFR function at low mass, it cannot explain an increase in this width with M⋆. A possibility is that the SFHs become more diverse as the mass increases, as a result of the numerous processes that reduce the star formation in massive galaxies.

We also show that the evolution of the median sSFR in a logarithmic scale decreases as −0.17M⋆. We note that the commonly adopted linear relation between log (sSFR) and log (M⋆) is not suitable for our data. Such a dependency with M⋆ at high mass could be reproduced by assuming exponentially declining SFH with τ having an inverse dependency with mass τ ∝ 1/ℳ. Several processes could reduce the sSFR as the stellar mass increases. Accretion of cold gas can be suppressed in hot gas halos (M_H> 10¹² M_⊙) leading to gas exhaustion in the central galaxies. This process should occur on longer timescales (>3–4 Gyr) than usually assumed to explain our observed trend. Another possibility is that the efficiency of the star formation is declining toward massive sources: disks are more settled and stabilized against fragmentation as the mass increases. The presence of a bulge could be crucial in reducing the star formation efficiency. Finally, a combined analysis of the sSFR functions and of the quiescent GSMF could constrain the relative importance between secular and quenching processes by setting the quenching timescale with the GSMF evolution.

---

Acknowledgments

We are grateful to the referee for the careful reading of the manuscript and his/her useful suggestions. We thank Samuel Boissier, Véronique Buat, Andrea Cattaneo, Jared Gabor, Sylvain De La Torre and Mark Sargent for useful discussions. We thank Carlota Gruppioni for providing her 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. This paper is based on 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. LMD acknowledges support from the Lyon Institute of Origins under grant ANR-10-LABX-66. A.K. acknowledges support by the Collaborative Research Council 956, sub-project A1, funded by the Deutsche Forschungsgemeinschaft (DFG).

References

Online material

Appendix A: Comparison between sSFR and sMIR_DM: impact of the mass loss

If the feeding in gas is driven by the growth of the dark matter halos and if the fraction of the infalling gas converted into stars stays constant, we expect that . However, only a fraction of the mass created stays trapped in an old stellar population. Therefore, is different from the sSFR defined as with R the return fraction depending on the age of the stellar populations (Renzini & Buzzoni 1986).

Assuming a constant R value, we expect sSFR′ = sMIR_DM/ (1−R). For the Chabrier (2003) IMF, the maximum value of R is 0.5. Therefore, we use this value to define the upper bound of the green area in Figs. 12 and 13.

However, R depends on time. The value of the stellar mass lost by a given galaxy depends on its SFH. In order to determine the mass lost along the galaxy history, for each galaxy in our sample we measure the difference between the total mass obtained by integrating the SFH (without taking into account the stellar mass loss) and the stellar mass. Then, we measure the median of these differences as a function of redshift. Figure A.1 shows the median stellar mass effectively lost as a function of redshift R_med(z). The lower bound of the green shaded area in Figs. 12 and 13 corresponds to sMIR_DM/ (1−R_med(z)) measured for low-mass galaxies. We check that the averaged stellar mass effectively lost in the higher mass bin would be within our two boundaries.

We note that the evolution of is flatter than sMIR_DM/ 2 (i.e., taking R as a constant), since R_med(z) is smaller at high redshift than at low redshift (stellar populations are younger).

Appendix B: Additional selection effects linked to the broadening of the M_⋆ – SFR relation

In this appendix, we investigate wether the criterion used to select star-forming galaxies could artificially create a broadening of the sSFR function with the mass.

We first check that the intrinsic evolution of the sSFR does not enlarge our estimate of σ significantly (e.g., Speagle et al. 2014). Using a simple model with sSFR ∝ (1 + z)^3.8 (our extreme value of b), we find that the broadening cannot be overestimated by more than 0.02 dex in our redshift bins owing to the sSFR intrinsic evolution.

Since the rest-frame colors are closely correlated with the sSFR, one could artificially modify the shape of the sSFR function depending on the rest-frame color cut used to separate quiescent and star-forming galaxies. In particular, it is unclear if the population that is transitioning from the star-forming to the quiescent population should be included in the analysis, and how it affects the sSFR function. As shown in Fig. B.1, the most massive galaxies tend to lie much closer to the transitioning area than the other star-forming galaxies. This result is not surprising since Peng et al. (2010) show that the probability of a star-forming galaxy being quenched is proportional to its SFR and mass (their Eq. (17)). In order to investigate the impact of the selection criterion, we move down and up the selection criterion by 0.3 mag which are the extreme values we could adopt (dotted lines in Fig. B.1). We find that this change has no impact on the sSFR functions of galaxies less massive than M⋆ < 10¹¹ M_⊙ and the value of σ remains above 0.4 dex at M⋆ > 10¹¹ M_⊙. While the measurements are slightly modified for the most massive galaxies, the overall shape of the sSFR functions is not affected. Therefore, our conclusions are not sensitive to the adopted limit used to select star-forming galaxies.

We also investigate if the uncertainties associated with the stellar mass could artificially create such broadening of the sSFR function. Indeed, uncertainties in the stellar mass could move the galaxies from one mass bin to another. Since the sSFR depends on the mass, it could artificially broaden the sSFR distribution. The galaxies at M⋆ < 10¹¹ M_⊙ could contaminate the most massive galaxies and artificially add galaxies with a larger sSFR since the sSFR decreases with the stellar mass. In order to test this effect, we use the semi-analytical model described in Sect. 7. We add random errors using a Gaussian distribution having a standard deviation of 0.1 dex (already larger than our expected uncertainties) to the predicted SFR and to the predicted mass. We find that the sSFR functions predicted by the model are almost unmodified. If we adopt uncertainties of 0.2 dex in mass and in SFR, we get unrealistic predictions at all masses. We note that systematic uncertainties are not discussed here since they shift the full distribution.

All Tables

All Figures