© ESO, 2018

1. Introduction

The majority of star-forming (SF) galaxies follow a relatively tight relation between stellar mass (M^∗) and star formation rate (SFR; Noeske et al. 2007; Brinchmann et al. 2004), the so-called star-forming main sequence (MS), up to high redshifts (Whitaker et al. 2012; Steinhardt et al. 2014; Tasca et al. 2015). This relation is key to understanding the evolution of galaxies and the interplay between secular effects, like gas depletion and gas infall, and more stochastic effects, such as mergers (Peng et al. 2014; Tacchella et al. 2016). For this reason, in the last decade, several works have been devoted to studying this relation and its evolution with redshift (Oliver et al. 2010; Karim et al. 2011; Rodighiero et al. 2011; Whitaker et al. 2014; Johnston et al. 2015; Shivaei et al. 2015; Kurczynski et al. 2016).

In spite of this wide range of studies, a general consensus on the MS slope and its evolution with redshift has not been reached. This is because different sample selections, SFR derivations, and fitting techniques influence the derived MS making the comparison between different works difficult. For example, as explained by Speagle et al. (2014), the MS derived using criteria that select blue, highly star-forming and not-dusty galaxies will be steeper than the MS derived with other selection mechanisms. Ideally, having no pre-selection of the star-forming population, by using colours or a specific SFR (sSFR) cut, would seem the right approach to avoid influencing the derived MS slope. However, a sample with no pre-selection of SF galaxies could be contaminated by quenched galaxies, particularly at high stellar masses, producing a flatter MS than that derived from samples containing only SF galaxies.

To have a MS that does not require pre-selection and is not contaminated by quenched galaxies, Renzini & Peng (2015) defined the MS as the ridge line in the SFR-stellar mass-number density three-dimensional (3D) distribution surface. With this definition they could clearly distinguish between SF galaxies in the MS and quenched galaxies below the MS at redshift z ~ 0.05, without influencing the slope of the derived MS.

As discussed exhaustively by Calzetti (2013), different SFR indicators are used across the wavelength spectrum, including X-ray, ultra-violet (UV), optical, infrared (IR), radio, and emission lines, based on direct stellar light, dust-processed light or ionised gas emission. SFR derived from the UV light has been used many times because it is dominated by the light of very young stars, and therefore is optimal to trace recent SF. Moreover, it is possible to derive the SFR in a wide range of redshifts observing the rest-frame UV redshifted to optical wavelengths. However, the presence of dust obscuring the UV light makes the derived SFR uncertain.

After the launches of the Spitzer Space Telescope and Herschel Space Telescope, observations at IR wavelengths started to be used more commonly to derive SFR, particularly for dusty galaxies. However, it is not possible to obtain IR data for large galaxy samples, and IR observations are less sensitive and have a smaller angular resolution than UV-optical observations. Moreover, several observations in the IR are necessary to reconstruct the IR spectral energy distribution (SED) at 3–1000 μm and therefore L_IR, which is then linked to the SFR. A series of relations between monochromatic IR luminosities and total IR luminosities have been derived during recent years (e.g. Chary & Elbaz 2001; Bavouzet et al. 2008; Rieke et al. 2009). However, these relations were obtained from different samples, at different redshift ranges, and using different monochromatic IR luminosities, and therefore the derived MS could be influenced by the relation used to convert from monochromatic IR luminosities and L_IR.

The aim of this work is to analyse the SFR–M^∗ plane from z = 0.5 to z = 3, by comparing the traditional MS derivation method, after pre-selecting the SF galaxies with sSFR cut, with the ridge line definition from Renzini & Peng (2015), obtained by fitting three Gaussian components to the sSFR distribution. In addition, we also analysed the effect of different SFR_IR derivation methods on the bright-end of the MS.

The structure of the paper is as follows. In Sect. 2 we present our sample selection, while in Sect. 3 we describe the derivation of photometric redshifts and stellar masses by SED fitting. In Sect. 4 we report the different methods used to derive SFR from the UV and from the IR. In Sect. 5 we analyse the sources on the SFR–M^∗ plane, while in Sect. 6 we summarise our main findings and conclusions. Throughout this paper, we adopt an H₀ = 70 km s^-1 Mpc^-1, Ω_M = 0.27, Ω_Λ = 0.73 cosmology and a Salpeter (1955) initial mass function over (0.1–100) M_⊙.

2. Sample selection

2.1. Data

The GOODS-S field (Giavalisco et al. 2004) has been targeted with deep observations at different wavelengths over more than a decade. It is centred at α(J2000) = 3^h32^m20^s and δ(J2000) = −27°48′′20^s and covers an area of approximately 10′ × 16′. In this study, we collect datasets in this field from the UV to the mid-IR, from both ground-based and space telescopes.

In particular, in the optical and near-IR, the GOODS-S field has been observed with the Hubble Space Telescope (HST), the Very Large Telescope (VLT) and Spitzer as part of the CANDELS survey (Grogin et al. 2011; Koekemoer et al. 2011). In this work, we use the CANDELS GOODS-S multi-wavelength photometric catalogue (Guo et al. 2013) that includes 17 bands: U Blanco/CTIO, U VLT/VIMOS, F435W, F606W, F775W, F814W and F850LP HST/ACS, F098M, F105W, F125W and F160W HST/WFC3, K_s VLT/ISAAC, K_s VLT/HAWK-I (Fontana et al. 2014) and 3.6, 4.5, 5.8 and 8.0 μm Spitzer/IRAC (Ashby et al. 2015). The sources of this catalogue are detected in the F160W band and we correct all fluxes for galactic extinction.

The HAWK-I K band data available in the Guo et al. (2013) catalogue do not cover the full GOODS-S field, therefore we complement these observations with other K VLT/HAWK-I observations (P.I. N. Padilla) that are shallower than the observations in the Guo et al. (2013) catalogue in the same band but cover the full field. To extract photometry in this band, we run SExtractor (Bertin & Arnouts 1996) using a 1.5σ threshold over four contiguous pixels and an aperture of 3′′. A low detection threshold is used to find a counterpart for all secure detections at 1.6 μm. We correct the 3′′ fluxes to total fluxes and ensure no offset with respect to Guo et al. (2013).

The GOODS-S field has been observed also at mid- and far-IR wavelengths. In particular, we include in our data set observations taken with the multiband imaging photometer (MIPS; Rieke et al. 2004) on Spitzer at 24 μm and observations done by the photodetector array camera and spectrometer (PACS; Poglitsch et al. 2010) on Herschel (Pilbratt et al. 2010) at 70 μm, 100 μm and 160 μm as part of the PACS (Poglitsch et al. 2010) evolutionary probe (PEP, Lutz et al. 2011) and GOODS-Herschel (Elbaz et al. 2011) programs.

2.2. Sample selection and counterpart identifications

Our main catalogue consists of 34 930 galaxies that are part of the CANDELS catalogue and are detected in the F160W band, with the addition of HAWK-I K band observations that cover the full GOODS-S sample and observations in the mid- and far-IR with Spitzer and Herschel. The CANDELS catalogue can be considered to correspond to a mass-limited sample, because galaxies are selected in F160W and the mass-to-light ratio is quite stable at these wavelengths. On the other hand, the subsample with mid- and far-IR observations is biased towards highly star-forming and dusty galaxies.

First, we add the HAWK-I/K band observations to the CANDELS catalogue by cross-matching the catalogue with the wide HAWK-I/K band observations with a 1′′ search radius. When more than one counterpart is present in the K-band (~ 7% of the catalogue), mainly because of the low threshold used for the source extraction in this band, we consider the brightest K band source as a counterpart, that is usually also the closest one.

For the Spitzer/MIPS 24 μm catalogue we consider only sources brighter than 80 μJy, for which the catalogue completeness is ~ 80% and the spurious sources are less than 10% (e.g. Papovich et al. 2004). In order to identify the optical counterpart of each 24 μm-detected source, we cross-match the MIPS catalogue and the CANDELS+HAWK-I/K catalogue using a 2′′ searching radius respect to the HAWK-I/K source positions. Among the initial 24 μm-detected galaxies, 733 (98%) have an optical-IR counterpart. The remaining 8 objects, detected at 24 μm, do not have any CANDELS counterpart within 2′′. However, visual analysis of these sources shows that they could be associated with extended sources or multiple objects around the MIPS detection. We do not consider 24 μm fluxes when multiple optical counterparts are associated with the same 24 μm source. To conclude, 644 sources have photometry also at 24 μm, corresponding to the 88% of the total MIPS catalogue and ~ 2% of the total CANDELS catalogue.

At far-IR wavelengths, Magnelli et al. (2013) derived two different catalogues from PEP observations: in one of them, sources are extracted blindly in the three bands and then they are crossmatched with a 24 μm-detected source catalogue; in the other one, sources are extracted using the 24 μm source positions as prior. We cross-match our catalogue with both PEP catalogues with a radius of 1′′ by comparing the 24 μm positions. Around 90% (585) of objects detected at 24 μm are also detected in at least one of the three PEP bands, in particular ~ 50% of the 24 μm sources are detected at 70 μm, 86% at 100 μm and 78% at 160 μm.

3. Photometric redshifts and stellar masses

Our catalogue is different, in the optical and near-IR, from the Guo et al. (2013) catalogue for the presence of the HAWKI-K band observations which covers the full field. For this reason we derive both the photometric redshift and stellar mass using our own catalogue. In particular, to determine simultaneously the photometric redshift and the stellar mass of each galaxy, we run the code LePhare (Arnouts et al. 1999; Ilbert et al. 2006) to obtain the best SED templates based on a χ² fitting procedure. We use Bruzual & Charlot (2003) templates considering solar metallicity, a range of exponentially declining star formation histories and delayed star formation histories, both with different characteristic time scales τ from 0.1 to 10 Gyr, and ages ranging from 0.05 to 13.5 Gyr. We allow for redshift between 0 and 6, we apply the Calzetti et al. reddening law (Calzetti et al. 2000) for the galactic extinction with A_V values from 0 to 4 and we allow LePhare to include nebular emission lines. We fit models by using all the optical bands and the two shortest-wavelength IRAC bands, because the two longest wavelengths have a lower signal-to-noise ratio (S/N) than the shortest wavelength bands and they could be contaminated by polycyclic aromatic hydrocarbon (PAH) emissions at low redshift.

We compare our photometric redshifts with spectroscopic redshifts available in the literature in the GOODS-S field. We use secure spectroscopic redshifts (z_spec) present in the ESO public compilation in the Chandra deep field south (CDFS; Cristiani et al. 2000; Croom et al. 2001; Bunker et al. 2003; Dickinson & GOODS Team 2004; Stanway et al. 2004a,b; Strolger et al. 2004; Szokoly et al. 2004; van der Wel et al. 2004; Doherty et al. 2005; Le Fèvre et al. 2005; Mignoli et al. 2005; Ravikumar et al. 2007; Popesso et al. 2009; Balestra et al. 2010; Silverman et al. 2010; Kurk et al. 2013; Vanzella et al. 2008) and spectroscopic redshifts from Morris et al. (2015), for a total of 2635 matches in the CANDELS catalogue. The deviation of our photometric redshifts from the spectroscopic ones is δz = | z_phot−z_spec | /(1 + z_spec) = 0.039 ± 0.032 for the full sample, excluding outliers with δz> 0.15 that correspond to ~ 11% of the sample (Fig. 1). These outliers are distributed in a wide range of stellar masses and SFRs and do not occupy specific regions of the SFR–M^∗ plane, therefore they do not have significant influence on the results of this paper (Appendix A).

We divide our sample into three redshift bins, in order to study the position of galaxies on the SFR–M^∗ plane: , and . In these three redshift bins there are 6881, 9578 and 8004, respectively, making a total of 24 463 galaxies.

In order to calculate the stellar mass completeness of our sample, we follow the empirical approach presented in Schreiber et al. (2015). In particular, the CANDELS sample is selected at 1.6 μm, where the light-to-mass ratio is quite stable. Indeed, a tight relation is present between the stellar mass and the observed luminosity at 1.6 μm that can be described as a power-law plus a Gaussian scatter in the log space. We derive this relation and the related scatter directly from the data in different redshift bins. From the 3σ detection limit in the F160W band it is possible to derive the corresponding 3σ limit at each redshifts for the L_{1.6 μm/(1 + z)}. Then, we randomly extract points inside the Gaussian scatter of this relation in a wide range of stellar masses and luminosities and we derive the percentages of galaxies above the 3σ detection limit in each 0.1 dex stellar mass bin. The 90% completeness limits are at log₁₀(M^∗/M₀) = 7.35, 7.95 and 8.25 at , and , respectively.

4. Star formation rates

For the full CANDELS sample, we do not use the SFR derived from the SED fitting, which highly depend on the assumed star formation history (Pforr et al. 2012), but we derive the SFR from the rest-frame UV continuum flux corrected for dust extinction using the Calzetti et al. (2000) reddening law and the A_v value derived from the SED fitting. Depending on the redshift, we use fluxes in the U-band, F435W, F606W or F775W band. These UV fluxes are then converted to SFR using the relation derived by Kennicutt (1998): (1)For the subsample of galaxies detected at 24 μm, we also derive the SFR from the dust-processed stellar light in the IR with the addition of the UV-based SFR explained previously, but not corrected for dust extinction. In particular, the total IR luminosity from 3 to 1000 μm is related to the SFR, as it was found by Kennicutt (1998): (2)However, we do not observe directly the total IR luminosity and several methods are present in the literature to estimate this total luminosity from a monochromatic luminosity. In this work, we consider five different methods to derive the total IR luminosity:

• the relation by Chary & Elbaz (2001) betweenIR and monochromatic luminosities, using the observed24 μm flux;

• the relation by Bavouzet et al. (2008) between IR and monochromatic luminosities, using the observed flux at 24 μm;

• the relation by Bavouzet et al. (2008) between IR and monochromatic luminosities, using the observed flux at 70 μm;

• the relation by Rieke et al. (2009) between IR and observed flux at 24 μm;

• a direct fit to the available IR observations, from 24 μm to 160 μm, by using IR templates presented in Rieke et al. (2009).

When IR observations are not available, we consider SFR derived from the UV and corrected for dust-extinction.

4.1. K-corrections

Every method explained in the previous paragraph converts a monochromatic luminosity at a reference wavelength into total IR luminosity. Because the IR data are available in specific bands, we need to apply a k-correction to obtain the luminosity at the reference wavelength of each method. We start from models of pure starburst galaxies (Rieke et al. 2009) that have IR luminosities from 10^9.75 to 10¹³L_⊙. We convolve each model with the 24 μm, 70 μm, 100 μm and 170 μm filters. We then calculate the k-correction by dividing the observed flux by the flux at the reference wavelength. We use this ratio to calculate the total IR luminosity and we compare it with the one associated to the used template. We iterate the process until the output IR luminosity corresponds to the IR luminosity of the used model.

4.2. IR luminosity derivation

We derive total IR luminosities by using four different relations between monochromatic and total IR luminosities and, independently, by fitting SED models to the total fluxes at 24 μm, 70 μm, 100 μm and 160 μm.

• 1.

  Chary & Elbaz (2001, here after C&E) studieddifferent samples of nearby galaxies and for galaxies withL_IR> 10¹⁰L_⊙ they derived three relations between total IR luminosity and monochromatic luminosity at 6.7 μm, 12 μm and 15 μm: (3)Close to each relation we show the redshift range in which we implement it, in order to apply the smallest k-correction possible from the flux at 24 μm and the flux at the wavelength considered in each relation.

• 2.

  In Bavouzet et al. (2008), the total IR luminosity is related to the monochromatic luminosity at 8, 24 and 70 μm by: (4)In this case we use observations at 24 μm (hereafter Ba24) and 70 μm (hereafter Ba70), when available. In order to apply the smallest k-correction possible, in Ba24, we use the first relation at z> 0.5 and the second one for closer objects. Similarly, in Ba70, we use the third relation at , the second one at and the first one otherwise.

• 3.

  In Rieke et al. (2009; hereafter Ri09), SFR is directly related to 24 μm observed flux by: (5)where D_L is the luminosity distance and values for A and B are listed in the same paper for different redshifts below 3. We do not apply this relation at z> 3. We derive IR luminosity for Ri09 by reversing Eq. (2).

Finally, when at least two observations are present in the IR, we use again the templates from Rieke et al. (2009) to directly fit observations at 24, 70, 100 and 160 μm using a χ² fitting procedure to estimate the best solutions among all templates. However, among different SED templates presented in the literature there is a large scatter on the PAH features strength. In order for our results to not be too dependent on the SED template used, we multiply by two the error bars associated to 24 μm fluxes, which contained PAH features between z ~ 0.8 and z ~ 2.8. We also derive L_IR by SED fitting without including observations at 24 μm, but the difference is less than 20% for 99% of the sample with observations in at least two bands.

In order to estimate the errors associated with the IR luminosity for each monochromatic relation, we consider both the propagation of the flux and photometric redshift errors inside each relation and the intrinsic scatter of each relation. The mean errors are 0.35 dex for C&E, 0.32 dex for Ba24, 0.43 dex for Ba70 and 0.21 dex for Ri09. The mean error of Ba70 is higher than the mean error of Ba24, because S/N associated with the flux at 70 μm are generally smaller than the ones at 24 μm.

On the other hand, in order to estimate IR luminosity errors for the SED fitting technique, we analyse the probability distribution, P(L_IR) ∝ e^{− χ2/ 2}, and we derive the full width half maximum associate to this distribution for each galaxy. The mean IR luminosity error associated with the SED fitting method is 0.09 dex.

In Fig. 2, we show the relation between IR luminosity and photometric redshift for each conversion formulae used to derive the bolometric infrared luminosity L_IR. Both Ri09 and C&E methods produce L_IR> 10¹⁴L_⊙, that are not present in other methods. The IR luminosities derived with Ba24 and Ba70 have the shallowest increment with redshift among the methods tested here and, in particular, IR luminosities remain below 10¹³L_⊙ for the majority of cases. IR luminosities derived fitting the IR photometry are always L< 10¹³L_⊙, because this is the maximum luminosity present in the used templates. Differences among the derived IR luminosities are due to small variations in the K-corrections, but mainly to differences among the used conversion formulae because of the different galaxy samples used to derive them.

In Fig. 3, we show the comparison between UV, dust-corrected SFR and SFR_IR + SFR_UV, where the IR SFR is derived by fitting the IR SED the UV SFR is not corrected for dust extinction, for the galaxies which have mid- and far-IR observations. The two SFRs are, overall, consistent, in particular at . However, the SFR_UV tends to be underestimated at high SFR and redshifts, particularly at . This comparison is different when comparing the SFR_UV with SFR_IR derived with other methods, but these quantities remain in overall agreement.

5. Analysis of sources on the SFR–M^∗ plane

We divide our sample into three redshift bins, , and , in order to study the SFR–M^∗ plane and the MS slope at different redshifts, and we use two different fitting techniques to derive it.

Firstly, we remove quenched galaxies from our sample, selecting only galaxies with log₁₀(sSFR/ yr^-1) > −9.8 and we derive the MS by fitting a linear regression using all non-quenched galaxies in each redshift bin. This method is similar to what is usually found in the literature and it allows us to make a direct comparison with previous work.

Secondly, we follow Renzini & Peng (2015) who analysed the 3D SFR–M^∗ relation (SFR vs. stellar mass vs. number of galaxies inside each SFR–M^∗ bin) and identified two peaks, one for quenched galaxies and one for MS galaxies, and defined the MS as the ridge line of the MS distribution. This definition has the advantage of not requiring a pre-selection to derive the MS position and, at the same time, is not influenced by the presence of quenched galaxies when deriving the MS. Following their approach, we divide galaxies into stellar mass bins of 0.25 dex and we fitted the sSFR distribution inside each stellar mass bin with three different Gaussian distributions, corresponding to the quenched galaxies (QG), the MS galaxies, and the starburst galaxies (SB). In Renzini & Peng (2015) the third component (starburst galaxies) was not evident because this population is rare at the low redshifts that they considered (0.02 <z< 0.085). We use the sSFR instead of the SFR, because the three components are more clearly separated in the sSFR distribution than in the SFR one. Other studies have analysed the sSFR distribution at different stellar masses (e.g. Ilbert et al. 2015), but they did not model this distribution to separate among SF, quenched, and starburst galaxies, and derive the MS. Finally, we fit an orthogonal distance regression (ODR) to the derived peak positions of the three Gaussian components in each stellar mass bin, in order to take into account the stellar mass bin size as well as the errors associated with each Gaussian peak position. The sequences derived by using this method are sensitive to the adopted stellar mass bin, however, the results are all consistent.

5.1. Single fitting to SFR–M^∗ plane

Here we discuss the SFR–M^∗ plane derived using a single linear regression fitting to all SF galaxies, as is commonly found in the literature. In Fig. 4 we show the linear regression for all SF galaxies, which are selected by log₁₀(sSFR/ yr^-1) > −9.8. The slopes that we derive are 0.88 ± 0.01, 0.84 ± 0.01 and 0.79 ± 0.01, at , and , respectively.

Figure 5 shows the comparison between our derived slopes and some values present in the literature. The scatter among different methods is very large and our values are within this range. In particular, our values are above the values derived by Speagle et al. (2014), which combine results from numerous previous works, and are more similar to recent values derived by Kurczynski et al. (2016). However, in Speagle et al. (2014) the MS is analysed for log₁₀(M^∗/M₀) > 9.5, therefore we re-derive the MS in a similar stellar mass range to compare with our work. The derived slopes are 0.71 ± 0.01, 0.67 ± 0.01 and 0.74 ± 0.01 at , and 2, respectively, which are consistent with the median inter-publication scatter around the relation found by Speagle et al. (2014). This analysis makes evident to which extent the derived MS slope is influenced by the considered stellar mass range.

Using IR-based SFR for the MS, we find that the slopes remain almost unchanged. This is because the galaxies at low-intermediate stellar masses are dominating the fits with this method because they are the most numerous, and for these galaxies the SFR is derived from the UV, because there are no observations in the IR.

5.2. Multi-Gaussian decomposition of the SFR–M^∗ plane

Here we discuss the SFR–M^∗ plane by separating SB, SF and QG galaxies fitting the sSFR distribution with three Gaussian components in different stellar mass and redshift bins.

In Fig. 6 we show the best Gaussian least squares fits of the sSFR distribution for our sample, in the three considered redshift bins and for stellar masses between log₁₀(M^∗/M₀) = 7.5 and 11.25. We analyse only stellar mass bins with more than 50 galaxies. The three Gaussian components used in the fit, one for each galaxy mode (QG, MS and SB), have different relative positions and relative heights, changing both with stellar masses and redshift. In particular, at high stellar masses the three components are more separated one from the other than at low stellar masses.

The SB component in the fit is necessary to describe the observed sSFR distribution in some stellar mass and redshift ranges where SB are not negligible, as analysed in Appendix B by comparing the akaike information criterion (AIC; Akaike 1973) for a fit with and without the SB component.

5.2.1. The three modes in the SFR–M^∗ plane

For each galaxy mode, we derive the Gaussian parameters from the non-linear least squares fit. We calculate the errors associated to each Gaussian parameter by adding in quadrature the errors estimated from the fit and the errors derived by performing a bootstrap analysis. Then, we use the peak positions and its error to fit an ODR and derive the sequence for SB, SF and quenched galaxies, as shown in Fig. 7. The SB sequence has a slope close to unity up to log₁₀(M^∗/M₀) = 11. The MS has a slope between 0.93 ± 0.03 and 0.83 ± 0.05, decreasing with the redshift, but a flattening is present at high stellar masses, in particular at z< 2. This flattening has been observed previously in other works with different methods (Whitaker et al. 2014; Tomczak et al. 2016) and, because this flattening is present in our data after selecting only SF galaxies, it is not due to contamination by quenched galaxies. However, the inclusion of IR observations influences this flattening, as is explained in detail in Sect. 5.2.3. When fitting only data points with log₁₀(M^∗/M₀) < 10.25, below the high mass bending, the slope is close to unity at all redshifts and with a mild evolution. On the other hand, the QG sequence is less defined than the other two, with the high-mass part (log₁₀(M^∗/M₀) > 10) well separated from the MS and the low-mass part very close to the MS, resulting in a slope between 0.54 and 0.56. A similar behaviour has been observed with low-z galaxies by Renzini & Peng (2015).

5.2.2. The evolution of the starburst and quenched fractions

In Fig. 8 we show the sSFR probability density distribution in the three analysed redshift bins, above each stellar mass completeness. In this plot it is already evident that the number of starburst and quenched galaxies varies with redshift. This evolution could be partially due to selection effects for the QG, because they could be below the detection limit. On the other hand, there is a clear increment in the number of SB between z ~ 0.75 and z ~ 1.5.

To further investigate this evolution, we study the fraction of the three different galaxy modes at different stellar masses and redshifts. In particular, in Fig. 9 we show the ratio between the number of starburst (quenched) galaxies with respect to the total number galaxies for different stellar masses and redshift bins. Galaxies in our sample are classified depending on the three Gaussian fits of the sSFR distribution in each stellar mass bin. The number of galaxies in each mode is derived by integrating the associated Gaussian component in each stellar mass and redshift bin. The errors associated to each fraction are derived propagating the Gaussian parameter errors, which are obtained adding in quadrature the errors derived from the fit and the ones derived from the bootstrap analysis.

In general, MS galaxies are the dominant mode, however, at log₁₀(M^∗/M₀) > 10.5 and z< 2, quenched galaxies are >20%. At log₁₀(M^∗/M₀) > 10, the number of QG decreases with increasing redshift, while it is almost constant at low-intermediate stellar masses (≲20%). Above this stellar mass (log₁₀(M^∗/M₀) ~ 10), QG are also more separated from the MS, while below they are very close to the MS and are not a clear separate peak. This could be a stellar mass of transition between two different dominant quenching mechanisms. A dip is present in the fraction of QG at log₁₀(M^∗/M₀) ~ 10, however, large error bars are presented at low stellar masses due to the proximity of QG and the MS. Therefore, it is not possible to understand if this dip is due to uncertainties or it is an underpopulated region between two different quenching modes. At z> 1, the fraction of QG decreases at log₁₀(M^∗/M₀) < 9, however, this fraction is underestimated in this redshift and stellar mass range because QG are below the 90% stellar mass completeness.

On the other hand, the fraction of SB galaxies increases with the redshift, particularly at log₁₀(M^∗/M₀) < 9. Moreover, this fraction increases by ≳2 between redshift ~0.75 and ~1.5. The fraction of SB at the stellar mass range of log₁₀(M^∗/M₀) > 10.5 could be contaminated by active galactic nuclei (AGN). However, these SB galaxies, when detected at 24 μm, do not show extreme IR fluxes and, moreover, only one galaxy in each stellar mass bin is classified as AGN at X-rays (Luo et al. 2017). Overall, looking at the total mass range above the 90% completeness (log₁₀(M^∗/M₀) = 8.25–11.25), SB galaxies are 5 ± 1% of the total number of galaxies at , 12 ± 1% at , and 16 ± 1% at .

For comparison with other works, Rodighiero et al. (2011) defined a starburst as a galaxy with sSFR 2.5σ above the MS and found that they correspond the 2–3% of galaxies with log₁₀(M^∗/M₀) = 10–11.5 at . In our work the SB fraction is generally higher in the same stellar mass range, probably because of the different definitions of starburst. Indeed, applying a cutoff at a defined sSFR does not take into account wings in the SB distribution as well as variation in the slopes or in the separation between the MS and the SB sequence with stellar masses. In agreement with our work, Sargent et al. (2012), by using the same data of Rodighiero et al. (2011) but a technique similar to the one use in this work, found a higher fraction of SB galaxies than Rodighiero et al. (2011).

Caputi et al. (2017) found that ~ 15% of all galaxies with log₁₀(M^∗/M₀) > 9.2 at z = 4−5 are SB. Although at these higher redshifts, SB galaxies are characterised by much higher sSFR. This fraction is slightly higher than the fraction of SB galaxies we found at z = 2–3 (11 ± 1%) in the same stellar-mass range, showing that the evolution of the fraction of SB in this stellar-mass range seems to continue up to z = 5.

5.2.3. Impact of using different IR conversion formulae to derive the bolometric IR luminosity

In Fig. 10, we show the SFR–M^∗ plane, with SFR_IR derived using the five different L_IR conversion formulae analysed in this paper. When IR observations are available we use the SFR_IR + SFR_UV, with UV SFR not corrected for dust extinction. On the other hand, when IR observations are not available, we use the dust-corrected, UV SFR. The MS is derived in each case by fitting an ODR to the Gaussian mean in the sSFR distribution in each 0.25 dex stellar mass bin. The best fit values derived using this method are listed in Table 1.

Observations in the IR are available only for galaxies with log₁₀(M^∗/M₀) ≳ 10, therefore in this work we can only analyse the impact of different L_IR conversion formulae on the high-mass end of the MS and, in particular, on its flattening. The MS slopes change by less than 0.11 with all methods and at all redshifts. In particular, at , the Ba24 method shows the flattest slope, 0.87 ± 0.04, while the SED fit method has the steepest slope, 0.94 ± 0.03. At , adding the IR observations decreases the slope by 0.01–0.02 dex for all methods, excluding the Ba24 and SED fit methods. At , the differences are slightly higher than in the lower redshift bin between −0.01 and +0.10 dex. However, taking into account the errors associated to each slope, all slopes derived by using only UV data or including IR observations are consistent with each other.

When including IR observations, the flattening of the MS changes depending on the L_IR conversion formulae used and, in some cases, the MS is consistent with no flattening, even at the lowest redshift bin (i.e. Ri09). Therefore, a careful inclusion of IR observations is necessary to properly quantify the flattening of the MS at high stellar masses.

6. Summary and conclusions

In this work we studied the SFR–M^∗ plane for a sample of 24 463 galaxies at with photometry from the UV through the near-IR. In particular, we firstly derived the MS with a linear fitting and, secondly, we studied quenched, star-forming, and starburst galaxies, separating them by fitting their sSFR with three Gaussian components. In addition, we analysed the impact of different bolometric IR luminosities conversion formulae on the MS of SF galaxies for a subsample with IR observations.

In our first method, we derived the MS of star-forming galaxies applying a method commonly used in the literature. Specifically, we made a linear regression for all SF galaxies, defined as all galaxies with log₁₀(sSFR/ yr^-1) > −9.8. The MS slope derived using this method depends on the considered stellar mass range, because of contamination from starbursts and quenched galaxies, whose fraction depends on the analysed stellar masses, and a possible flattening of the MS at high stellar masses. Therefore, it is important to carefully distinguish between quenched, star-forming, and starburst galaxies, in particular when considering a narrow stellar mass range, and, moreover, it is important to consider the same stellar mass range when comparing different works. Indeed, our results are consistent with the values found in the literature once we considered the same stellar mass range.

As a second step, we decomposed the sSFR distribution (in each stellar mass bin). Three components are evident (quenched, SF, and SB galaxies) with different importance depending on the stellar mass and the redshift considered. In particular, the SB population increases in number with increasing redshift, particularly for galaxies with log₁₀(M^∗/M₀) < 9 between z ~ 0.75, where starburst galaxies are less than 20%, and z ~ 1.5, where they are between 20 and 30%, with respect to the total number of galaxies in the same stellar mass range. The SB fraction increases also with decreasing stellar mass. Overall, SB galaxies are ~ 5% at z = 0.5–1, ~ 12% at z = 1–2 and ~ 16% at z = 2–3, with respect to the total number of galaxies with log₁₀(M^∗/M₀) = 8.25−11.25 of those redshifts. On the other hand, the quenched galaxies increase in number with stellar masses and they are a substantial number at log₁₀(M^∗/M₀) > 10.5 at z < 2 (20–60%). The fraction of quenched galaxies also decreases with redshift. The slope of the MS changes from 0.93 ± 0.03 at 0.5 ≤ z < 1 to 0.83 ± 0.05 at 2 < z < 3, decreasing with the redshift. However, a slight bending is present at high stellar masses (log₁₀(M^∗/M₀) > 10.25) and the MS at lower masses is consistent with a slope of ~1 and a mild evolution with redshift. A tight sequence of SB galaxies is also present, with a slope close to unity. On the other hand, quenched galaxies show a less clear sequence, they seem separated from the MS at log₁₀(M^∗/M₀) > 10 and closer to the MS at lower stellar masses.

In addition, we analysed the impact of using different IR luminosity derivation methods on the high-mass end of the MS. In particular, we considered two methods to convert 24 μm luminosity to L_IR (C&E and Ba24), one to convert 70 μm luminosity to L_IR (Ba70), one to convert flux at 24 μm to SFR (Ri09), and a fifth one that derives the L_IR by fitting 24–160 μm observations with empirical IR SED templates.

In general, the slope of the MS changes by <0.11 when including the SFR_IR, with all different L_IR derivation methods, both with a simple linear regression and decomposing the sSFR distribution with the Gaussian components. However, a proper analysis and inclusion of the SFR_IR is essential when analysing the MS slope at high stellar mass and its possible flattening, which is indeed visible with UV-based SFR. Indeed, both the curvature and the starting stellar mass of the bending depend on the method used to derive the L_IR and, remarkably, with some of the methods analysed in this work the MS is consistent with no flattening at high stellar masses even at the lowest redshift bin.

Overall, using different fitting techniques and different L_IR derivation methods, we derived slopes between 0.87 and 0.93, at , between 0.84 and 0.87, at , and between 0.79 and 0.93 at . A dependence of the MS slope on the stellar mass is evident, therefore a flatter (steeper) MS slope could be derived limiting the study to high (low) stellar masses. Depending on the analysed redshift and stellar mass range, it could be necessary to distinguish between quenched, SF, and starburst galaxies in order to derive the MS without contaminants.

Acknowledgments

This paper is based on observations taken by the CANDELS Multi-Cycle Treasury Program with NASA/ESA HST, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. L.B. and K.I.C. acknowledge the support of The Nederlandse Onderzoekschool voor de Astronomie (NOVA). K.I.C. also acknowledges funding from the European Research Council through the award of the Consolidator Grant ID 681627-BUILDUP. We thank N. D. Padilla for useful discussions.

References

Appendix A: Analysis of redshift outliers

When comparing the photometric redshifts of the main sample used in this work with some spectroscopic redshifts presented in the literature, there are ~11% of outliers with δz = | z_phot−z_spec | /(1 + z_spec) > 0.15 (see Sect. 3). In this appendix we analyse the distribution of these outliers in the SFR–M^∗ plane and their effect on the general conclusions of the paper.

Figure A.1 shows the position of all galaxies with spectroscopic redshift and of all redshift outliers in the SFR–M^∗ plane. Outliers are 8 ± 1% at redshift 0.5 ≤ z < 1 and 7 ± 1% at redshift 1 ≤ z < 2 and 2 ≤ z< 3, with no evident evolution in the analysed redshift range. From Fig. A.1 is evident that the redshift outliers do not systematically occupy a specific locus on the plane. Indeed, at 0.5 ≤ z < 1 outliers are distributed among starburst, MS galaxies and quenched galaxies. Only at 1 ≤ z < 2, it is evident that there are almost no outliers among the quenched galaxies, therefore the fraction of quenched galaxies at high stellar masses could be even higher than what has been derived in the paper. At the highest redshift bin the statistic starts to be poor and outliers are distributed among both SB and MS galaxies.

Overall, outliers do not occupy specific regions of the SFR–M^∗ plane, therefore they do not qualitatively change the general conclusions of the paper. Only the fraction of QG with log₁₀(M^∗/M₀) > 10 at 1 ≤ z < 2 may be underestimated. However, a more quantitative analysis of the effect of these outliers is not possible because galaxies with spectroscopic redshift are not representative of the full sample, they have generally log₁₀(M^∗/M_⊙) ≳ 9; they are a minor fraction, only ~8% of the CANDELS sample at 0.5 <z ≤ 3; and they do not have an homogeneous selection criteria.

Appendix B: Analysis of the need for a SB component

In this paper, we fit the sSFR distribution in each stellar mass and redshift bin using three Gaussian components. In this appendix we repeat the fit using two Gaussian components, one for the MS galaxies and another for the quenched galaxies, in order to see whether the SB component is really necessary to describe the data. To compare the two fits we derive the akaike information criterion (AIC; Akaike 1973) defined as: (B.1)where ℒ is the likelihood function and k is the number of parameters in the model. This criterion takes into account the goodness of the fit and the complexity of the model at the same time. Among different models, the one that best represents the data with the minimum amount of parameters corresponds to the model with the smallest AIC value. We calculate the AIC value for both the three-component (Fig. 6) and the two-component fits (Fig. B.1) and we compare the AIC values in each stellar mass and redshift bin.

Unsurprisingly, the three-component fit is favoured by the data in the lowest-stellar-mass bins that also correspond to the stellar-mass bins where SB are more numerous. The stellar-mass bin limit below which the third component is necessary to fit the sSFR distribution increases with redshift from log(M^∗/M₀) < 8.25 at 0.5 ≤ z< 1 to log(M^∗/M₀) < 9 at 2 ≤ z< 3, following the SB fraction evolution with redshift. At high stellar masses galaxies are mainly in the MS or quenched, therefore a third component is not necessary.

As a summary, at low stellar masses and at high redshifts, SB galaxies are not negligible and a third component is necessary to properly describe the sSFR distribution.

All Tables

All Figures