The Star Formation Main Sequence in the Hubble Space Telescope Frontier Fields

A key step to understanding galaxy evolution is estimating galaxy redshifts and physical properties, in particular the star formation rate (SFR) and the stellar mass (M). Star formation rate and stellar mass directly probe the process of gas conversion into stars and subsequent stellar mass build-up. The analysis of large, statistically significant galaxy samples has allowed the establishment of the existence of a well-defined relation between the SFR and the stellar mass, called main sequence (MS hereafter), a subject of study of a large number of papers in the literature (Brinchmann et al. 2004; Daddi et al. 2007; Elbaz et al. 2007; Noeske et al. 2007; Santini et al. 2009; Peng et al. 2010; see Speagle et al. 2014 for a compilation of results and more references).

The MS has traditionally been parameterized as a power law of the form $\mathrm{log}\,\mathrm{SFR}=\alpha \mathrm{log}\,M+\beta$. In recent years, however, many studies have reported evidence that the MS flattens at high masses ($M\sim {10}^{10\mbox{--}11}\,{M}_{\odot }$) (e.g., Magnelli et al. 2014; Whitaker et al. 2014; Schreiber et al. 2015), with the turnover mass increasing with redshift (Lee et al. 2015; Tasca et al. 2015; Tomczak et al. 2016) and at the same time becoming less evident (Whitaker et al. 2014; Schreiber et al. 2015). Recently, the analysis of Dunlop et al. (2017) (see also Koprowski et al. 2016) ascribes the observed flattening at $z\gtrsim 2$ to the inability of correctly recovering the SFR in deeply obscured high-mass galaxies when rest-frame optical/UV tracers are adopted (although some of the works above are based on FIR estimators, this is qualitatively consistent with the weakening of the mass turnover with redshift).

The existence of the MS is universally recognized and there is general consensus of an increasing normalization with redshift, associated with a higher rate of gas accretion in the early universe. However, the details and in particular the MS slope vary from one study to the other, ranging from ∼0.6 to 1, most strongly depending on the sample selection and SFR tracer adopted (Santini et al. 2009; Rodighiero et al. 2014; Speagle et al. 2014). Nevertheless, the majority of recent results point toward a roughly linear and unevolving slope (Whitaker et al. 2014; Schreiber et al. 2015; Tasca et al. 2015; Tomczak et al. 2016), at least out to $z\sim 4$.

Physically, the existence itself and the tightness of the MS relation (scatter of ∼0.25–0.4 dex, Rodighiero et al. 2011; Speagle et al. 2014; Schreiber et al. 2015) suggests similarity in the gas accretion histories of galaxies. The flattening at high masses has been interpreted as due to either the contribution of the bulge to the stellar mass (while the SFR comes primarily from the disk, Schreiber et al. 2015) or the onset of quenching processes (Tasca et al. 2015). According to some work (e.g., Renzini & Peng 2015; Whitaker et al. 2015), it seems to arise from an incorrect separation of passive and quenched galaxies from star-forming ones.

While the majority of galaxies occupy the locus of the MS, outliers are also observed with intense levels of SFR given their stellar mass. These two populations have been associated with different growth mechanisms (Daddi et al. 2010; Genzel et al. 2010; Elbaz et al. 2011): MS galaxies are thought to grow on long timescales as a consequence of smooth gas accretion from the Intergalactic Medium (IGM), while MS outliers, also called starbursts, seem to be triggered by mergers and form stars with high efficiency, although this view is currently being debated (e.g., Kennicutt & Evans 2012; Narayanan et al. 2012; Santini et al. 2014; Mancuso et al. 2016). The latter, being very rare, despite their high level of SFR, seem to contribute modestly to the cosmic star formation history (SFH; Rodighiero et al. 2011; Sargent et al. 2012; Lamastra et al. 2013).

The latest deep multi-band surveys allowed us to observe fainter and fainter galaxies and probe the MS to stellar masses as faint as $\sim {10}^{9}\,{M}_{\odot }$ at $z\gt 3$ (see, e.g., the results of the HST CANDELS survey, Salmon et al. 2016). A different, complementary approach involves exploiting gravitational lensing, i.e., using galaxy clusters as natural telescopes able to amplify the light emitted by background sources. The HST Frontier Fields program (Lotz et al. 2017) combined the capabilities of the two photometric cameras ACS and WFC3 onboard HST with the power of gravitational lensing to realize the deepest images ever produced in six different pointings, each targeting one $z\sim 0.3\mbox{--}0.5$ cluster together with their background galaxies. In parallel, HST produced six images of close-by fields, referred to as parallel fields.

In this paper, we take advantage of the HST Frontier Fields program to investigate the MS relation at $1.3\leqslant z\lt 6$ down to very low masses ($\mathrm{log}\,M/{M}_{\odot }\sim 7.5$ at $z\lesssim 4$ and $\mathrm{log}\,M/{M}_{\odot }\sim 8$ at higher redshifts). The paper is organized as follows. Section 2 describes the data set and the method applied to estimate stellar masses and SFRs. Sections 3–5 present our results, respectively, on the MS relation, its scatter, and the evolution of the SFRs per unit mass. Finally, Section 6 summarizes the results. In the following, we adopt the Λ-CDM concordance cosmological model (H₀ = 70 km s⁻¹ Mpc⁻¹, ${{\rm{\Omega }}}_{M}=0.3$ and ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.7$) and a Salpeter (1955) IMF. All magnitudes are in the AB system.

We use the multiwavelength catalogs of the first four HST Frontier Fields clusters (field 1: Abell2744; field 2: MACS0416; field 3: MACS0717; and field 4: MACS1149) developed within the ASTRODEEP project,⁹ publicly released by our team¹⁰ and also available through an interactive CDS interface.¹¹

The four fields are presented in Merlin et al. (2016, fields 1 and 2) and Di Criscienzo et al. (2017, fields 3 and 4). The 10 bands photometric catalogs include the B435, V606, and I814 ACS/HST bands; the Y105, J125, JH140, and H160 images from WFC3/HST; the K_s band from Hawk-I/VLT (fields 1 and 2) and MOSFIRE/Keck (fields 3 and 4); and IRAC 3.6 and 4.5 μm data. The source detection was carried out with SExtractor (Bertin & Arnouts 1996) on the H band image, with the addition of faint sources detected on a weighted average of the processed Y, J, JH, and H images. Photometry in the other bands has been derived with the template-fitting code T-PHOT (Merlin et al. 2015), which uses galaxy shapes in the detection band as prior information. The method adopted to assemble the multiwavelength catalogs is described in details in Merlin et al. (2016) and Di Criscienzo et al. (2017). For our analysis, we applied a magnitude cut of H < 27.5. This limit corresponds to a detection completeness, assessed through simulations taking into account the variation of the intracluster light across the image (Merlin et al. 2016), of 90%–95% for point-like sources and 50%–80% for extended disks with 02 half-light radius, depending on the field.

Photometric redshifts, when spectroscopy is unavailable (i.e., for 94% of the parent sample), have been obtained as the median value of six different techniques, as detailed in Castellano et al. (2016, fields 1 and 2) and Di Criscienzo et al. (2017, fields 3 and 4). This approach reduces the scatter and the fraction of outliers as well as minimizes systematics possibly associated to each individual method (Dahlen et al. 2013).

Observed physical parameters, such as stellar masses, have been derived with a SED fitting approach. We used Bruzual & Charlot (2003) stellar population models and included nebular continuum and line emission following Schaerer & de Barros (2009) as detailed in Castellano et al. (2014). We adopted delayed-τ SFHs of the form $\psi (t)\propto {t}^{2}/{\tau }^{3}\cdot \exp (-t/\tau )$. These are rising-declining laws peaking at $t=2\tau$. We set $\tau \geqslant 0.6$ Gyr, forcing high-redshift galaxies ($z\gtrsim 4$) to experience the increasing phase. This choice is motivated by the results of Salmon et al. (2016), who measured on their high-z galaxy sample a SFH increasing with time as a power law with index 1.4 or maybe higher, due to possible incompleteness at low stellar mass. However, because stellar masses are only mildly sensitive to the choice of the SFH (Santini et al. 2015, although there may be exceptions for peculiar galaxy populations, see also, e.g., Michałowski et al. 2012, 2014), this choice does not significantly affect our results. Stellar metallicities range from $0.02\,{Z}_{\odot }$ to solar, the dust reddening $E(B-V)$ is comprised between 0 and 1.1 and the extinction law can be either Calzetti et al. (2000) or SMC (Prevot et al. 1984). 1σ uncertainties on the stellar masses were computed by accounting for all the solutions within ${\chi }_{\min }+1$.

K and IRAC fluxes affected by heavy blending issues (Merlin et al. 2015) have not been considered in the SED fitting, neither for photo-z (Castellano et al. 2016) nor for physical parameters. We checked that their removal has no systematic effect on the inferred stellar masses. We excluded from our analysis sources at $z\geqslant 4$ whose K and IRAC fluxes have been ignored as their SEDs result highly unconstrained and, as a consequence, their inferred properties are highly unreliable. They amount to ∼7% of the H-selected sample.

SFRs were estimated from observed UV rest-frame photometry using the same technique as Castellano et al. (2012). Briefly, the UV slope β was obtained by fitting the observed photometric points adopting a power-law approximation for the 1280–2600 Å spectral range, and the Meurer et al. (1999) relation was assumed to infer the extinction used to correct the 1600 Å luminosity. Finally, the dust-corrected luminosity was converted into a SFR estimate with the Kennicutt & Evans (2012) factor. This method avoids issues associated with SED parameter degeneracy, particularly serious for the SFR estimate. However, it limits the analysis to $z\geqslant 1.3$, where the available photometry samples the UV rest-frame spectral region. After visual inspection of sources with extreme values of the UV slope, we removed sources with $\beta \gt 1$ or $\beta \leqslant -3.5$, mostly caused by noisy photometry (∼8% of the sample in the redshift range analyzed). 1σ errors on the SFRs were computed by considering the uncertainty in fitting the UV slope due to photometric errors as well as a scatter of 0.55 dex around the Meurer et al. (1999) relation (see also Fudamoto et al. 2017), with the two contributions summed in quadrature. We note that recent ALMA results have found consistency with the Meurer et al. (1999) law at z ∼ 3, while suggesting a possible evolution toward lower attenuations at z larger than 4 or 5 (Capak et al. 2015; Bouwens et al. 2016; Fudamoto et al. 2017). The nature of the attenuation law at high redshift is however still debated in the recent literature, and consistency with a Calzetti/Meurer law is found by some other studies (e.g., Scoville et al. 2015; Cullen et al. 2017).

Intrinsic physical parameters were inferred after correcting for the magnification due to gravitational lensing. The procedure adopted to estimate the magnification factor from the shear and mass surface density maps provided by different teams¹² is described in detail in Castellano et al. (2016) (see Meneghetti et al. 2016 and Priewe et al. 2017 for references and for a comparison between different models). We adopted the median magnification, thus excluding possible outlier values associated with a particular model. For fields 1 and 2, at variance with the released catalogs, we only adopted the most updated maps available to date (i.e., v3). This results into the adoption of 6, 8, 7, and 7 models, respectively, for the four cluster fields. The magnification factor for the entire sample has a median value of 2, but can reach values as high as ∼60, as shown in Castellano et al. (2016).

To study the MS relation and inferring its slope, it is important to consider a sample that is complete above a given stellar mass. Our sample is selected in the H band, but there is a dispersion in the relation between H band magnitude and stellar mass. Therefore, we assessed the mass completeness limit of our sample by following the procedure outlined in Fontana et al. (2004). Briefly, considering a passively evolving dust-free model, we compute the minimum mass given the magnitude limit of the data and the M/L distribution distribution allowed by the adopted stellar library. Then we use the observed M/L distribution close to the magnitude limit at different redshifts to account for the galaxies that are not observed. We set the mass limit as the mass above which 90% of objects are observed. This redshift-dependent limit is shown as a red vertical line in Figure 1. However, a fraction of objects intrinsically below the mass limit are actually observed above this threshold, thanks to gravitational lensing, boosting the flux from intrinsically fainter, less massive objects. The magnification only depends on the relative position of the lensed and lensing galaxies and not on their physical properties. In particular, it does not depend on their SFR, which will be randomly distributed, in a statistical sense. Based on this argument, we included in our analysis galaxies with intrinsic stellar mass below the mass limit, pushed above the threshold by gravitational lensing magnification. Figure 1 shows how gravitational lensing allows us to extend the analysis up to a factor of ∼10 below the strict observational mass limit (hatched region of the final mass distribution represented by the blue histogram).

Zoom In Zoom Out Reset image size

After excluding 17 objects at $z\gt 3$ which have been a posteriori visually inspected due to their suspect very high SFR with respect to the MS, and which turned out to have bad fits due to very noisy photometry, the final sample includes 1711 sources in the redshift range $1.3\leqslant z\lt 6$.

3.1. Observed Main Sequence

We show in Figure 2 the relation between SFR and stellar mass, i.e., the MS, in different redshift bins, from z = 1.3 to z = 6. By means of a 2σ-clipping¹³ analysis, we compute average values in bins of stellar mass (reported in Table 1). We recover the mass turnover at high masses, at least out to z = 3, while above this value the turnover moves to higher masses where we lack the required statistics due to the small sky area covered. At the opposite tail, we observe an extension of the MS relation down to stellar masses only rarely probed in previous studies.

Zoom In Zoom Out Reset image size

Table 1.  SFR—Stellar Mass Relation Average Values for the Cluster Frontier Fields

$1.3\leqslant z\lt 2$			$2\leqslant z\lt 3$			$3\leqslant z\lt 4$			$4\leqslant z\lt 5$			$5\leqslant z\lt 6$
(1)	(2)	(3)	(1)	(2)	(3)	(1)	(2)	(3)	(1)	(2)	(3)	(1)	(2)	(3)
7.85	−0.84	0.54	7.84	−0.81	0.36	7.74	−0.33	0.62	7.90	0.04	0.33	⋯	⋯	⋯
8.25	−0.33	0.53	8.29	−0.25	0.45	8.28	0.07	0.39	8.36	0.21	0.41	8.31	0.63	0.45
8.72	0.02	0.41	8.73	0.08	0.39	8.78	0.38	0.31	8.78	0.51	0.36	8.78	1.03	0.63
9.23	0.44	0.47	9.21	0.61	0.37	9.23	0.79	0.30	9.25	0.83	0.37	9.22	1.13	0.32
9.71	1.08	0.32	9.70	1.19	0.19	9.68	1.14	0.21	9.72	0.96	0.24	9.76	1.54	0.86
10.16	1.40	0.31	10.17	1.44	0.76	10.13	1.22	0.54	10.24	1.66	0.74	10.20	2.76	0.91
10.68	1.00	0.69	10.76	1.68	0.38	10.68	2.28	0.24	10.77	1.59	0.68	⋯	⋯	⋯

Note. (1): $\langle \mathrm{log}\,M/{M}_{\odot }\rangle$, (2): $\langle \mathrm{log}\,\mathrm{SFR}/({M}_{\odot }\,{\mathrm{yr}}^{-1})\rangle$, and (3): ${\sigma }_{\mathrm{log}\mathrm{SFR}}$. Average values and associated scatter have been obtained with a 2σ clipping procedure.

Download table as:  ASCIITypeset image

In the bottom right panel of Figure 2, we report the average SFR values in bins of stellar mass at different redshifts. We observe an overall mild increase in the normalization with redshift and a mild flattening of the slope (at least out to $z\sim 5$, although we note that the highest redshift bin is noisier as it includes fainter sources and is more affected by poor number statistics).

We fit our entire sample (not just the average values) with the linear relation

$$\begin{eqnarray}&&\mathrm{log}\,\mathrm{SFR}=\alpha \mathrm{log}(M/{M}_{9.7})+\beta \end{eqnarray} \tag{ 1 }$$

where ${M}_{9.7}={10}^{9.7}\,{M}_{\odot }$, avoiding the high-mass region where the linearity breaks up. To determine this threshold, we used the redshift-dependent best-fit turnover mass obtained by Tomczak et al. (2016) and its value at z = 4 for higher redshift. We adopt a 2σ-clipping procedure to remove the outliers: we iteratively reject all points more than 2σ apart from the best-fit, where σ is calculated as the standard deviation of the residuals. The best-fits are shown as dashed curves in Figure 2. The evolution of the MS slope and normalization is shown in Figure 3 (black solid boxes). To verify the robustness of our results against details associated with the magnification correction, in the left panel of Figure 3 we also show the best-fit MS parameters obtained with each individual lensing model instead of the median magnification. The global picture is unaffected, and the scatter among the various results can be regarded as a measure of the observational uncertainty on the MS parameterization to be associated with lensing modeling. To further confirm that our results are not driven by uncertainties in the magnification corrections, we also compare the MS parameters with those obtained from the parallel fields, where the magnification is very modest.

Zoom In Zoom Out Reset image size

Overall, we apparently observe a decline in the slope from ∼1 at $z\lt 3$ to ∼0.6 at $z\sim 6$, associated with an increase in the normalization across the entire redshift range studied.

3.2. Correction of the Eddington Bias

Since measurement errors in the physical parameters can be high, especially for high-redshift faint sources, we tried to take into account their influence on the inferred MS relation. In other words, we correct here the fitted MS for the Eddington bias, that is, we try to recover the intrinsic, "true" MS, that, once convolved with measurement errors, gives rise to the observed relation.

We run a Monte Carlo simulation comprising the following steps:

• 1.  
  Starting from an arbitrary MS, we populate the relation with a realistic log N–log S distribution and a given Gaussian scatter. To this aim, we start from the initial observed sample without applying any magnitude or mass cut. To increase the statistics, we replicated the observed distribution 100 times. Given the observed redshift and stellar mass, we assign to each source the SFR predicted by the MS, starting from the observed relation and the observed scatter (the latter is plotted as a black line in Figure 4). For stellar masses lower than $\mathrm{log}\,M/{M}_{\odot }=8$ and larger than $\mathrm{log}\,M/{M}_{\odot }=10.3$, where our measurement of the scatter is affected by poor number statistics, we adopt a standard value of 0.3 dex.
• 2.  
  We randomly perturb each source in stellar mass and SFR adding Gaussian noise based on the 1σ uncertainties on these parameters associated with each object.
• 3.  
  We fit the resulting sample using exactly the same method adopted on real data, i.e., applying the magnitude and mass cuts and adopting a 2σ-clipping procedure to infer the MS.
• 4.  
  We iteratively modify the input MS, both in slope and normalization, until the best-fit simulated relation matches the observed MS within a tolerance of 0.02 in both parameters. Similarly, we modify the input scatter until the simulated scatter is consistent the observed one. At $\mathrm{log}\,M/{M}_{\odot }\lt 8$ and $\mathrm{log}\,M/{M}_{\odot }\gt 10.3$ we assume an input scatter equal to the scatter in the $8\lt \mathrm{log}\,M/{M}_{\odot }\lt 8.8$ and $9.6\,\lt \mathrm{log}\,M/{M}_{\odot }\lt 10.3$ mass bins, respectively, as our data prevent us to reliably evaluate the scatter outside this mass range.
• 5.  
  We iterate the procedure 100 times. The resulting parameters and associated uncertainties of the "true" MS as well as its scatter in bins of stellar mass are obtained by computing the average and standard deviation of the best-fit values and of the standard deviation of log(SFR) over the 100 iterations.

Zoom In Zoom Out Reset image size

The "true" MS are plotted as solid lines in Figure 2. The best-fit parameters are reported in Table 2, where the uncertainties have been obtained by adding in quadrature the uncertainties on the "true" MS with those on the observed MS. The evolution of the "true" MS is shown by the hatched region in Figure 3. Once corrected for the Eddington bias, the slope of the MS is consistent with unity across the 1.3–6 redshift range, while the normalization increases with redshift.

Table 2.  Main Sequence Best-fit Parameters (Corrected for the Eddington Bias)

Redshift	α	β
$1.3\leqslant z\lt 2$	1.04 ± 0.03	1.01 ± 0.04
$2\leqslant z\lt 3$	1.16 ± 0.03	1.22 ± 0.03
$3\leqslant z\lt 4$	1.02 ± 0.04	1.37 ± 0.03
$4\leqslant z\lt 5$	0.94 ± 0.06	1.37 ± 0.05
$5\leqslant z\lt 6$	0.92 ± 0.15	1.99 ± 0.13

Note. Data have been fitted to the functional shape $\mathrm{log}\,\mathrm{SFR}=\alpha \mathrm{log}(M/{M}_{9.7})\,+\beta$, where ${M}_{9.7}={10}^{9.7}\,{M}_{\odot }$, with a 2σ clipping procedure.

Download table as:  ASCIITypeset image

The Eddington bias has a negligible effect at $z\lesssim 3$, but it is responsible for a flattening of the relation at higher redshift. It is therefore crucial to take into account the effect of measure uncertainties to reliably estimate the slope of the MS.

3.3. Comparison with Previous Works

The existence of a relation on the SFR–stellar mass diagram has been interpreted as evidence that galaxy growth is likely regulated by cold gas accretion from the IGM on long timescales (e.g., Dekel et al. 2009). In this scenario, the tightness of the MS is related to the level of similarity of the gas accretion histories (see e.g., Shimakawa et al. 2017). It is interesting to investigate whether the scatter around the MS depends on the stellar mass.

The observed scatter around the MS is computed as the standard deviation of log(SFR) after a 2σ-clipping procedure to remove the MS outliers, and it is shown by the black lines in Figure 4 as a function of the stellar mass (we do not show the highest redshift bin due to the higher noise level and poorer statistics, which make it difficult to accurately evaluate the scatter). If we focus on the common mass range ($9\lesssim \mathrm{log}\,M/{M}_{\odot }\lesssim 10$), the observed scatter is comparable with that measured by previous works (e.g., Salmon et al. 2016; Schreiber et al. 2015).

The intrinsic scatter around the MS is smaller than the observed one, as the true scatter is convolved with the MS evolution within each bin as well as with the scatter arising from uncertainties in deriving the redshift and the physical parameters (Speagle et al. 2014). To account for the effect of observational uncertainties, we consider the scatter assumed as input to the simulation described in Section 3.2 (blue stars and lines in Figure 4), such that the simulated scatter (red circles) matches the observed one. This can be regarded as the intrinsic scatter, except for the fact that it is convolved with the evolution of the MS within each redshift bin. In two of the redshift-mass bins the simulated scatter turns out to be significantly larger than the observed one even by assuming no intrinsic scatter in the simulation. This suggests that the large observational uncertainties completely dominate the scatter of the data points around the MS and/or the simplified assumptions in our simulation (e.g., the assumption of Gaussian noise) hamper the evaluation of the intrinsic scatter in these bins.

Overall, however, we observe an indication, although not very robust, that the intrinsic scatter decreases with the stellar mass, although this trend seems to be reversed (or is impossible to evaluate, as discussed above) at $3\lt z\lt 4$. This tentative trend is in contrast with recent studies reporting no dependency on the stellar mass (Whitaker et al. 2012; Salmon et al. 2016; Schreiber et al. 2015; Kurczynski et al. 2016), although not extending to masses lower than ${10}^{9}\,{M}_{\odot }$ (except for the latter study). It instead agrees with the results of Salim et al. (2007), who finds that in the local universe the scatter declines by 0.11 dex per order of magnitude in stellar mass from 10⁸ to ${10}^{10.5}\,{M}_{\odot }$.

An increase in the scatter around the MS at low stellar masses as well as low redshift is a natural outcome of theoretical models of galaxy formation, and arises from the hierarchical structure formation and/or from the highly efficient stellar feedback processes in small halos (e.g., Merlin et al. 2012; Lamastra et al. 2013; Hopkins et al. 2014). To visualize this effect, we plot in Figure 4 the predictions of the semi-analytical model of Menci et al. (2014), which, not surprisingly, forecasts a smaller scatter than the observations as it is not affected by measurement errors. In a hierarchical scenario, high-mass galaxies have a large number of progenitors already collapsed at high-z, as they formed in high density regions of the primordial dark matter density field. The SFHs of these progenitor galaxies are peaked at high redshift since the gas cooling and the star formation efficiency are extremely efficient at early epochs. In addition, due to the higher compactness and lower virial temperatures of high-z dark matter halos, supernova feedback is less efficient at heating/expelling the gas from the galaxy. Conversely, the smaller number of progenitors of low-mass galaxies, and their different collapse time allowing also more prolonged star formation activity, lead to a larger variety of their SFHs.

We investigate here the evolution of the average, or typical, specific SFR, i.e., the SFR per unit stellar mass (sSFR = SFR/M), which captures information on the mass build-up process across cosmic time. Indeed, theoretical models in which galaxy growth is dominated by cold accretion predict that the sSFR should increase with redshift as ${(1+z)}^{2.25}$ (e.g., Dekel et al. 2009; Davé et al. 2011a).

As the sSFR is not constant with the stellar mass (i.e., the slope of the MS is not unitary at all redshifts), we derived the average sSFR at constant stellar mass. We considered the $\mathrm{log}\,M/{M}_{\odot }\,=9.5\mbox{--}10$ mass bin (average value $\lt \mathrm{log}\,M/{M}_{\odot }\gt \sim 9.7$). Despite adopting a fixed stellar mass at all redshifts implies considering a different galaxy population at different times, this allows an easy comparison with previous results computed at the reference mass of $\mathrm{log}\,M/{M}_{\odot }=9.7$, and allows us to extend the analysis of the evolution of the sSFR outside the redshift range probed by the present work (i.e., at $z\lt 1.3$ and $z\gt 6$). We computed the observed average values by means of a 2σ-clipping procedure directly on the galaxies in this mass bin, without using the best-fit MS relation (black solid squares). Given that the MS slope at high-z is sometimes found to be (slightly) sub-linear, the sSFR is a decreasing function of the stellar mass. Not to run the risk of underestimating the average sSFR, it is therefore essential to ensure completeness at the chosen mass bin. The chosen mass bin is well above the completeness level in our data. Finally, to take into account the effect of observational uncertainties, we make use of the Eddington-corrected MS estimated from the simulation (Section 3.2) and compute the typical sSFR at $\mathrm{log}\,M/{M}_{\odot }=9.7$ (black open squares). We adopt the same procedure for the data from the literature when the average value is not available at the relevant mass.

Our results are in line with the compilation of data from previous works (see Figure 5). Until a few years ago, observations appeared to show that the sSFR increased steeply from $z\sim 0$ to $z\sim 2$ (e.g., Noeske et al. 2007; Santini et al. 2009; Karim et al. 2011) and flattened at higher redshift (e.g., Stark et al. 2009; González et al. 2010; Labbé et al. 2010; McLure et al. 2011; Reddy et al. 2012), in sharp contrast with theoretical predictions (see e.g., Weinmann et al. 2011). The improvement in the dust-correction estimates and the inclusion of nebular lines in SED modeling have led to higher sSFR estimates at high redshift, alleviating the tension with theoretical models. A moderate increase in the sSFR compared to its value at $z\sim 2$ is now found in many studies (Stark et al. 2013; González et al. 2014; Tasca et al. 2015; Koprowski et al. 2016; Mármol-Queraltó et al. 2016), especially at $z\gtrsim 5$. Nevertheless, the evolution seems overall milder than that expected in an accretion-dominated scenario, although some results are getting closer (de Barros et al. 2014; Faisst et al. 2016). We note, however, that given the scatter in the MS, a flat trend cannot be ruled out.

Zoom In Zoom Out Reset image size

We observe an evolution by a factor of ∼2 from $z\sim 5.5$ to $z\sim 1.6$. Although a steeper trend is found when correcting for the Eddington bias, the large scatter associated to the highest redshift bin prevent us from drawing firm conclusions in this regard. The observed decrease is lower than that predicted by theoretical models in the same redshift interval (a factor of ∼8), implying that some sort of star formation suppression at high redshift is required to match the observations (feedback effects can modify the trend, see e.g., Davé et al. 2011b), that the effect of metallicity should be taken into account (Krumholz & Dekel 2012) and that other processes besides pure accretion are at play in assembling galaxies, such as mergers (e.g., Tasca et al. 2015; Faisst et al. 2016). The evolution of the sSFR is a critical observable to understand all these processes regulating star formation and galaxy growth. We note that a tension between observations and theoretical models is also observed in terms of absolute values of the sSFR (see the predictions of Menci et al. 2014 semi-analytical model). This effect is closely connected with the well-known overall lower normalization of the predicted MS compared to the observed one (e.g., Santini et al. 2009; Lamastra et al. 2013).

We studied the main sequence relation of star-forming galaxies, i.e., the relation between galaxy SFR and stellar mass, in the first four HST Frontier Fields. Thanks to gravitational lensing amplification of faint sources from massive foreground galaxy clusters, we could extend the analysis to masses lower than has usually been possible with the deepest data available before the Frontier Fields program and a factor of ∼10 lower than most studies. We investigated the redshift range $1.3\leqslant z\lt 6$ and probed stellar masses down to $\mathrm{log}\,M/{M}_{\odot }\sim 7.5$ at $z\lesssim 4$ and $\mathrm{log}\,M/{M}_{\odot }\sim 8$ at higher redshifts. We find that the MS relation extends to such low masses. At the opposite side, we recover the mass turnover at high masses at $z\lt 3$, i.e., where we have enough statistics at $\mathrm{log}\,M/{M}_{\odot }\gtrsim 10.5$. We run an accurate Monte Carlo simulation to take into account the effect of measure uncertainties on stellar masses and SFRs and correct for the Eddington bias. Such step is crucial to reliably measure the slope of the MS, which tends to be flattened by observational errors at $z\gtrsim 3$. We find increasing normalizations with redshift and an approximately linear unevolving slope.

The combination of deep HST observations and gravitational lensing allows us to observe tentative evidence that the scatter around the MS relation increases at low stellar masses, although we cannot make a strong claim. If verified, this result implies a higher level of uniformity in the SFHs at high masses, in agreement with observations in the local universe and with the predictions of theoretical models of galaxy evolution.

We confirm the previously observed mild increase in the average sSFR from $z\sim 1.6$ to $z\sim 5.5$ by a factor of ∼2, in tension with the steeper trend predicted by accretion-driven models (a factor of ∼8 in the same redshift interval). Being able to reproduce the evolution of the sSFR is crucial to understand the processes regulating the star formation and galaxy growth at different cosmic epochs and assessing the relative importance of gas accretion versus merging as well as the role of feedback mechanisms.

James Webb Space Telescope observations in the next future will allow to extend the analysis of the MS to even lower masses. Discerning between a MS that keeps extending to faint sources and a broken relation will provide hints on physical mechanisms such as feedback and reionization, or cold versus warm dark matter scenarios.

We thank the anonymous referee for the detailed review and useful suggestions. The research leading to these results has received funding from the European Union Seventh Framework Programme ASTRODEEP (FP7/2007-2013) under grant agreement n° 312725. MJM acknowledges the support of the National Science Centre, Poland through the POLONEZ grant 2015/19/P/ST9/04010. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 665778.