Star Formation in Galaxies at z ∼ 4–5 from the SMUVS Survey: A Clear Starburst/Main-sequence Bimodality for Hα Emitters on the SFR–M* Plane

Studying star formation in galaxies at high redshifts is crucial to understanding the early stages of galaxy evolution. Over the past 10 yr, a picture has emerged indicating that the global cosmic star formation rate (SFR) density increased after the big bang, until reaching a peak about 10 billion years ago, and then declined until today (e.g., Hopkins & Beacom 2006; Behroozi et al. 2013; Madau & Dickinson 2014). That peak is likely the product of two effects: mainly the net increase in the number density of galaxies that make the bulk of star formation in the first few billion years and possibly, but this is less clear, the fact that the SFR of individual galaxies may increase over that time. Indeed, our current knowledge of how star formation and stellar mass buildup proceeded over the first few billion years is still very sparse.

In most observed star-forming galaxies up to at least $z\sim 3$, the instantaneous SFR appears to correlate (within some scatter) with the stellar mass. This correlation on the SFR–M* plane is the so-called "galaxy star formation main sequence" (e.g., Noeske et al. 2007; Rodighiero et al. 2010; Elbaz et al. 2011). For these galaxies the specific SFRs ($\mathrm{sSFR}\equiv \,\mathrm{SFR}/{M}^{* }$) are roughly constant, which implies the existence of a scaling relation between gas consumption and galaxy growth (see, e.g., Popping et al. 2015; Lagos et al. 2016). In addition, there exists a minority of star-forming galaxies that are characterized by significantly higher sSFRs and are usually called "starbursts." The starburst phase is presumably a temporary state in which the galaxy is taken out of the main sequence, due to some kind of perturbation that temporarily enhances the star formation. Starburst galaxies are very rare in the local universe and have been found to constitute a small fraction of the dusty star-forming galaxies observed at redshifts $z\sim 0.3$ to $z\sim 2$ (Rodighiero et al. 2011; Sargent et al. 2012). This inferred weak evolution in the starburst fraction with redshift is based on the analysis of the most luminous dusty galaxies, as only these galaxies could be detected by the last generation of infrared (IR) telescopes, particularly the Herschel Space Telescope. Therefore, these results mainly concern massive galaxies, and it is not obvious whether they can be extrapolated to all star-forming galaxies.

A main limiting factor to understanding galaxy evolution in the high-z universe has been the lack of deep galaxy surveys over significantly large areas of the sky. Such surveys could provide sufficient statistics and dynamic range to investigate how star formation, as well as different physical properties, varies among different stellar-mass galaxies. For star formation, in particular, another limitation is that the only tracer currently available for large galaxy samples is the rest-UV flux shifted into observed optical wavelengths at high z. UV fluxes are extremely sensitive to dust extinction, and thus the derived SFRs can be very uncertain. Far-IR and (sub)millimeter photometry is the ideal complement to UV photometry for the purpose of obtaining total SFR estimates, but single-dish far-IR telescopes are insufficiently sensitive to systematically study representative galaxy samples in the high-z universe. As an alternative, Balmer lines, especially ${\rm{H}}\alpha \,\lambda 6563$, which is much less affected by dust than the UV spectral continuum, are very suitable SFR tracers widely used for galaxies up to z ∼ 2–3. However, Hα is shifted into the mid-IR regime ($\lambda \gtrsim 3\,\mu {\rm{m}}$) at $z\gt 3$, making its detection prohibitive for current spectrographs.

When emission lines have sufficiently large equivalent widths (EWs), however, their presence can be inferred even from broadband photometry. At redshifts $3.9\lesssim z\lesssim 4.9$, the Hα line is encompassed by the $3.6\,\mu {\rm{m}}$ filter passband at the Infrared Array Camera (IRAC; Fazio et al. 2004) on board the Spitzer Space Telescope (Werner et al. 2004). This fact has been exploited by different authors to investigate the presence of intense Hα emitters at these redshifts. Shim et al. (2011) analyzed a sample of 74 galaxies at similar redshifts and concluded that $\sim 70 \%$ of these galaxies have an excess flux at $3.6\,\mu {\rm{m}}$ with respect to the stellar continuum. More recently, Stark et al. (2013) and Smit et al. (2016) conducted similar analyses, based on galaxy samples with spectroscopic and photometric redshifts, and confirmed the presence of intense Hα emitters at $z\sim 4\mbox{--}5$. The percentage of galaxies displaying an "Hα excess" and the derived rest EWs varied according to the different selection effects, but in general they found median values of Hα EW $\approx 300\mbox{--}400\,\mathring{\rm A}$. These EWs are on average much larger than those observed in the local universe and are broadly consistent with an increase with redshift, as determined by Fumagalli et al. (2012), Faisst et al. (2016), and Mármol-Queraltó et al. (2016). All these studies have been very valuable in raising awareness on the increasing importance of nebular emission up to at least $z\sim 5$. However, none of them have analyzed a sufficiently representative galaxy sample at $z\sim 4\mbox{--}5$, which would allow for a more complete investigation of the implications of these results within the context of galaxy evolution.

The Spitzer Matching Survey of the UltraVISTA ultra-deep Stripes (SMUVS; M. Ashby et al. 2017, in preparation) is a Spitzer Exploration Science Program that has collected ultra-deep IRAC 3.6 and $4.5\,\mu {\rm{m}}$ images over a significant part of the COSMOS field (Scoville et al. 2007), making it the largest quasi-contiguous Spitzer field to analyze the high-z universe. The unprecedented level of statistics of SMUVS allows for a detailed study of galaxy properties and evolution over more than three decades in stellar mass at high z. In this paper we analyze the large SMUVS galaxy sample containing almost 6000 sources at $3.9\leqslant z\leqslant 4.9$ to investigate the presence of prominent Hα emitters as a function of stellar mass, along with the implications for the SFR versus M* relation at those high redshifts.

This paper is organized as follows. In Section 2 we describe the utilized data sets, source catalog construction, photometric redshifts, and stellar mass derivations. In Section 3 we explain our selection criteria for "Hα excess" galaxies among all SMUVS galaxies at $3.9\leqslant z\leqslant 4.9$. In Section 4 we present our results on the derived nebular line EW, SFRs and resulting SFR–M* and sSFR–M* relations, and the inferred cosmic SFR density (SFRD) at $z\sim 4\mbox{--}5$. We also discuss a rare $z\sim 4\mbox{--}5$ overdensity in the SMUVS field. Finally, in Section 5 we summarize our findings and discuss our conclusions. Throughout this paper we adopt a cosmology with ${H}_{0}=70\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$, ${{\rm{\Omega }}}_{M}=0.3$, and ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.7$. All magnitudes and fluxes are total, with magnitudes referring to the AB system (Oke & Gunn 1983). Stellar masses correspond to a Chabrier (2003) initial mass function (IMF), except where explicitly stated otherwise (see Section 4.3).

The SMUVS program (PI Caputi; M. Ashby et al. 2017, in preparation) has collected ultra-deep Spitzer 3.6 and 4.5 μm data over the region of the COSMOS field (Scoville et al. 2007) overlapping the three UltraVISTA ultra-deep stripes (McCracken et al. 2012) with deepest optical coverage from the Subaru telescope (Taniguchi et al. 2007). The SMUVS mosaics utilized here correspond to the almost final survey depth, which reaches an average integration time of $\sim 25\,{\rm{h}}$r/pointing (including previously available IRAC data in COSMOS; Sanders et al. 2007; Ashby et al. 2013; Steinhardt et al. 2014; Ashby et al. 2015). The considered UltraVISTA data correspond to the third data release (DR3), which in the ultra-deep stripes reaches an average depth of ${K}_{s}=24.9\pm 0.1$ and $H=25.1\pm 0.1$ (2'' diameter; $5\sigma$).⁵

A complete description of our SMUVS source multiwavelength catalog construction and spectral energy distribution (SED) fitting is provided in S. Deshmukh et al. (2017, in preparation). Here we only summarize our main steps. First, we extracted sources on UltraVISTA HKs average stack mosaics of the three relevant ultra-deep stripes, using the software SExtractor (Bertin & Arnouts 1996) with a detection threshold of 1.5σ over 5 contiguous pixels. The positions of these sources were considered as priors to perform iterative PSF-fitting photometric measurements on the Spitzer SMUVS 3.6 and 4.5 μm mosaics, using the DAOPHOT package on IRAF with empirical point-spread functions (PSFs) obtained from stars in the field (in each stripe separately). The PSF-fitting algorithm converged for $\sim 70 \%$ of the sources. For the remaining ones, we measured directly IRAC aperture fluxes in $2.4\,\mathrm{arcsec}$ diameter circular apertures at the position of the UltraVISTA sources and corrected them to total fluxes multiplying by a factor of 2.13, which was determined from the curves of growth of stars in the field. Overall, we found that $\sim 95 \% \mbox{--}96 \%$ of the UltraVISTA ultra-deep sources are detected in at least one IRAC band (and 93%–94% in both bands). The comparison of the resulting IRAC number counts with those obtained in the deeper S-CANDELS survey (Ashby et al. 2015) indicates that our resulting SMUVS catalog is 80% (50%) complete at [3.6] and [4.5] = 25.5 (26.0) total magnitudes.

For all these sources, we measured 2 arcsec diameter circular photometry on 26 broad, intermediate, and narrow bands from the U through the K_s bands, using SExtractor on dual-image mode with the UltraVISTA HKs stacks as detection images, and applied corresponding point-source aperture corrections in each band. After cleaning for galactic stars based on a (J–[3.6]) versus (B–J) color–color diagram (e.g., Caputi et al. 2011) and masking regions of contaminated light around the brightest sources, our final catalog contains ∼291,300 UltraVISTA sources with at least one IRAC-band detection over a net area of 0.66 deg². This is our SMUVS parent catalog with 28-band photometry (U through K_s + IRAC) for SED-fitting analysis.

The PSF-fitting technique assumes that all sources are point-like: on the IRAC images, this is indeed a reasonable assumption for virtually all sources with $[3.6]\gt 21$ mag (see Figure 25 in Ashby et al. 2013). Besides, all the multiwavelength photometry that we considered for SED fitting has been measured on circular apertures (and corrected to total), so the IRAC photometry based on PSF fitting is consistent with this procedure.

We performed the source SED fitting using the ${\chi }^{2}$-minimization code LePhare (Arnouts et al. 1999; Ilbert et al. 2006) over our 2 arcsec based total flux source catalog. We used Bruzual & Charlot (2003) templates corresponding to a single stellar population and different exponentially declining star formation histories (SFHs), all with solar metallicity (${Z}_{\odot }$), allowing for the addition of emission lines. In Appendix A, we discuss the impact of considering two possible metallicities (${Z}_{\odot }$ and $0.2\,{{\rm{Z}}}_{\odot }$) for the SED fitting. To account for internal dust extinction, we used the Calzetti et al. (2000) reddening law. Throughout this paper we will consider that the extinction derived in this way from the SED fitting of each galaxy is the same for the spectral continuum and lines at the same wavelengths, which should be a reasonable assumption for the majority of our galaxies (e.g., Reddy et al. 2010; Shivaei et al. 2015). In Appendix B, we explore the effects of assuming a different dust extinction law that is directly dependent on the UV slope of each galaxy.

As in Caputi et al. (2015), in the case of nondetections we adopted $3\sigma$ flux upper limits in the broad bands and ignored narrow and intermediate bands. For the minority ($\lt 14 \%$) of sources corresponding to [3.6] or [4.5] $\gt \,22$ galaxies with a [3.6] or [4.5] $\lt \,23$ neighbor within a 3 arcsec radius, we performed the SED fitting only on 26 bands (ignoring the IRAC bands).

We obtained zero-point corrections to the photometry from a first LePhare run and then performed a second, definitive run taking into account these corrections. From LePhare's output we obtained photometric redshifts and stellar mass estimates for $\gt 99.9 \%$ of our sources. The COSMOS field has a large amount of spectroscopic data that are very useful to assess the quality of the obtained photometric redshifts (e.g., Lilly et al. 2007): the resulting dispersion of the $| {z}_{\mathrm{phot}}-{z}_{\mathrm{spec}}| /(1+{z}_{\mathrm{spec}})$ distribution in our sample, based on ∼14,000 galaxies with reliable spectroscopic redshifts, is $\sigma =0.026$ with 5.5% outliers (S. Deshmukh et al. 2017, in preparation). This general photometric redshift quality is similar to that restricted to the redshift range analyzed here, i.e., $3.9\leqslant z\leqslant 4.9$. At these high redshifts we have 55 sources with spectroscopic redshifts and obtain $\sigma =0.023$ with 7.3% outliers.

From the obtained SMUVS photometric redshift catalog, we excluded $\sim 1 \%$ of sources because their best-fit ${z}_{\mathrm{phot}}$ were incompatible with their detection at short wavelengths—see criteria in Caputi et al. (2015). Our final SMUVS catalog with photometric redshift determinations contains ∼288,000 galaxies, including 5925 at $3.9\leqslant z\leqslant 4.9$, which is the relevant redshift range in this work. As part of LePhare's output, we also have stellar mass estimates for all these sources.

At redshifts $3.9\leqslant z\leqslant 4.8$ the entire (Hα λ6563 + [N ii] λλ6548, 6583 + [S ii] λλ6716, 6730) emission line complex is encompassed within the IRAC 3.6 μm filter and the (Hα λ6563 + [N ii] λλ6548, 6583) complex is present up to $z\approx 4.9$. At these same redshifts, no prominent emission line is present in the IRAC $4.5\,\mu {\rm{m}}$ filter. At $z\gtrsim 4.4$, the [O ii] $\lambda \lambda 3727$ enters the K_s band, although its EW is typically significantly smaller than the Hα EW (e.g., Mouhcine et al. 2005). The (${\rm{H}}\beta$ $\lambda 4861$ + [O iii] $\lambda \lambda 5007$) complex lies in the gap between the K_s and $3.6\,\mu {\rm{m}}$ bands.

Only emission lines or line complexes with sufficiently large EW can produce a significant flux excess in a broadband. For the SED model templates considered in this work, we determined that producing a $3.6\,\mu {\rm{m}}$ magnitude brightening of at least 0.1 mag with respect to the continuum at $3.9\lesssim z\lesssim 4.9$ requires an (Hα + [N ii] + [S ii]) rest-frame $\mathrm{EW}\gtrsim 150\,\mathring{\rm A}$.

With all these considerations in mind, in order to determine which of the 5925 galaxies at $3.9\leqslant z\leqslant 4.9$ have a significant $3.6\,\mu {\rm{m}}$ photometric excess with respect to the continuum, we remodeled their SEDs running LePhare again without emission lines, excluding the $3.6\,\mu {\rm{m}}$ band, and adopting the fixed ${z}_{\mathrm{phot}}$ that were previously determined in the original run (with all bands and emission lines). We selected the "Hα excess" sources by imposing the following:

• (1)  
  ${[3.6]}^{\mathrm{obs}}-{[3.6]}^{\mathrm{mod}}\lt -0.1$, where ${[3.6]}^{\mathrm{obs}}$ is the observed $3.6\,\mu {\rm{m}}$ magnitude and ${[3.6]}^{\mathrm{mod}}$ is the filter-convolved model magnitude obtained from the new LePhare SED fitting without emission lines;
• (2)  
  ${[4.5]}^{\mathrm{obs}}-{[4.5]}^{\mathrm{mod}}\gt -0.1$, so there is no excess light at $4.5\,\mu {\rm{m}};$
• (3)  
  ${K}_{s}^{\mathrm{obs}}\mbox{--}{K}_{s}^{\mathrm{mod}}\gt -0.1$, so there is no excess light in the K_s band either;
• (4)  
  if the source has any IRAC magnitude $\gt 22$, then it should not have any neighbor with an IRAC magnitude $\lt 23$ within a 3 arcsec radius.

These conditions altogether make for a conservative approach that guarantees that the continuum around the (Hα + [N ii] + [S ii]) complex is well modeled and we do not overestimate the $3.6\,\mu {\rm{m}}$-band flux excess. Although we did not set upper limits for the ${K}_{s}^{\mathrm{obs}}\mbox{--}{K}_{s}^{\mathrm{mod}}$ and ${[4.5]}^{\mathrm{obs}}-{[4.5]}^{\mathrm{mod}}$ differences, we note that their median values are ∼0.12 and $\sim -0.005$ mag, respectively, i.e., comparable to or lower than the typical photometric error bars, i.e., the observed fluxes are consistent with the continuum level.

We found that 1904 galaxies satisfy conditions (1) to (4), i.e., about one-third of our SMUVS galaxy sample at $3.9\leqslant z\leqslant 4.9$. In the rest of this paper, we refer to these galaxies as the "Hα excess" sample.

For both the "Hα excess" and "non-Hα excess" galaxy samples analyzed below, the integral of the total normalized P(z) (i.e., sum of P(z) for all the galaxies) at $3.9\leqslant z\leqslant 4.9$ is $\sim 95 \%$, indicating that the vast majority of our galaxies are in the correct redshift range to study the presence of Hα, and no conclusion in this paper is affected by the possible presence of a minor fraction of redshift interlopers.

4.1. (Hα + [N ii] + [S ii]) Emission versus Stellar Mass in Galaxies at $z=3.9\mbox{--}4.9$

We converted the $3.6\,\mu {\rm{m}}$ flux excess of each galaxy with "Hα excess" into the corresponding (Hα + [N ii] + [S ii]) rest-frame EW. For galaxies at $z\gt 4.8$, the excess corresponds to (Hα + [N ii]) only. We constructed a grid of $3.6\,\mu {\rm{m}}$ photometric excesses corresponding to different rest $\mathrm{EW}$ for each galaxy SED template (characterized by an SFH and age) at different redshifts. To do this, we modeled each line complex as a single Gaussian, with rest-frame widths of 40 and $180\,\mathring{\rm A}$ for the (Hα + [N ii]) and (Hα + [N ii] + [S ii]), respectively. Note that the exact values of these widths are irrelevant, as they cancel out when recovering the corresponding line EW and fluxes, as discussed by Shim et al. (2011). We convolved the galaxy SED templates, each with an added line complex of different rest EW, with the IRAC $3.6\,\mu {\rm{m}}$ filter transmission curve. In each case, we measured the resulting IRAC $3.6\,\mu {\rm{m}}$ magnitude and worked out the magnitude difference with respect to the original SED template with no emission lines. Finally, we considered each real galaxy photometric redshift and best SED model parameters to infer the corresponding line complex $\mathrm{EW}$ from the observed $3.6\,\mu {\rm{m}}$ flux excess, by interpolating the values in the model grid.

Figure 1 (left) shows the derived (Hα + [N ii] + [S ii]) ((Hα + [N ii]) at $z\gt 4.8$) rest EW versus stellar mass for our "Hα excess" galaxies. The stellar masses shown in this and all the remaining plots in this paper are those obtained from the original LePhare SED fitting (i.e., those including the $3.6\,\mu {\rm{m}}$ band and line emission). We see that these rest $\mathrm{EW}$ span values between $\sim 150\,\mathring{\rm A}$ and several thousand $\mathring{\rm A}$. The median values vary with stellar mass: they are $\langle \mathrm{EW}\rangle \approx 1000\,\mathring{\rm A}$ for galaxies with ${M}^{* }=(1\mbox{--}4)\times {10}^{9}\,{M}_{\odot }$, but only $\langle \mathrm{EW}\rangle \approx 400\,\mathring{\rm A}$ for galaxies with ${M}^{* }=(1.5\mbox{--}4)\times {10}^{10}\,{M}_{\odot }$. These values are broadly consistent with those obtained by Smit et al. (2016).

Zoom In Zoom Out Reset image size

The incidence of "Hα excess" sources also changes with stellar mass. As we can see from Figure 1 (right), the fraction of sources with "Hα excess" reaches 37%–40% at ${M}^{* }={10}^{9}\mbox{--}{10}^{10}\,{M}_{\odot }$ but decreases to $\approx 18 \%$ at ${M}^{* }\sim 5\,\times {10}^{10}\,{M}_{\odot }$. As we discuss in the next section, the fact that the percentage of galaxies with large Hα EW is higher among intermediate-mass galaxies indicates that the ongoing star formation (with respect to past start formation) is more important in them than in more massive galaxies at $3.9\leqslant z\leqslant 4.9$.

Interestingly, the "Hα excess" incidence rises again at ${M}^{* }\gtrsim 6\times {10}^{10}\,{M}_{\odot }$, albeit with moderate median EW. This reversing trend is puzzling and may suggest that star formation does not subside in the most massive galaxies before it does in the more typical massive ones. As we discuss below, an alternative and more likely explanation is that the $3.6\,\mu {\rm{m}}$ flux excess in the most massive galaxies is linked to the presence of active galactic nuclei (AGNs).

4.2. The SFR–M* and sSFR–M* Planes

4.2.1. SFR from Hα Emission

From the (Hα + [N ii] + [S ii]) EW we can derive an Hα-based SFR for each "Hα excess" galaxy. To do this, we obtained emission-line fluxes assuming that the modeled SED continuum is constant over the entire line complex wavelength range. As explained in Section 2, we also assumed that the extinction is the same for the continuum and the lines at the same wavelengths, which should be reasonable for the bulk of our galaxies (e.g., Reddy et al. 2010; Shivaei et al. 2015). This extinction has been obtained from the SED fitting of each galaxy. To obtain the net Hα contribution to these fluxes, we considered that f(Hα) = 0.63f(Hα + [N ii] + [S ii]) and f(Hα) = 0.81f(Hα + [N ii]), which is valid for solar metallicities (Anders & Fritze-van Alvensleben 2003). Then we converted the clean, derived ${\rm{H}}\alpha$ luminosities into SFRs using the Kennicutt (1998) relation:

$$\begin{eqnarray}&&\mathrm{SFR}({\rm{H}}\alpha )\,({M}_{\odot }\,{\mathrm{yr}}^{-1})=7.936\times {10}^{-42}\times {L}_{{\rm{H}}\alpha }\,(\mathrm{erg}\,{{\rm{s}}}^{-1}),\end{eqnarray} \tag{ 1 }$$

where the resulting SFR corresponds to a Salpeter (1955) IMF over $(0.1\mbox{--}100)\,{M}_{\odot }$, so we divided them by a factor of 1.69 to rescale them to a Chabrier IMF.

Figure 2 shows the resulting SFR and sSFR versus stellar mass diagrams, which span more than three decades in stellar mass. In these plots we include the SFR and sSFR based on the derived ${\rm{H}}\alpha$ luminosities for the "Hα excess" galaxies and upper limits for all other galaxies at $3.9\leqslant z\leqslant 4.9$. These upper limits are based on the minimum line complex EW that we can detect from flux excess in the $3.6\,\mu {\rm{m}}$ band. The vertical dashed lines indicate the 50% stellar-mass completeness level of our sample, which is ${\mathrm{log}}_{10}{M}^{* }\approx 9.2$ at $3.9\leqslant z\leqslant 4.9$. Note that the 80% stellar-mass completeness level is ${\mathrm{log}}_{10}{M}^{* }\approx 9.6$, but no main conclusion in this paper would change if we restricted our analysis to this higher stellar mass limit.

Zoom In Zoom Out Reset image size

The asterisks on the SFR–M* plot indicate galaxies that have a $24\,\mu {\rm{m}}$ counterpart within a 2 arcsec radius in the Mid-Infrared Photometer for Spitzer (Rieke et al. 2004) S-COSMOS catalog (Sanders et al. 2007). These sources correspond mostly to galaxies with ${\mathrm{log}}_{10}({M}^{* })\gt 10.8$. Given the limited depth of the $24\,\mu {\rm{m}}$ catalog, the detection of sources at $3.9\leqslant z\leqslant 4.9$ suggests that they are likely AGNs. Unfortunately, an SED power-law analysis in the IRAC bands is not useful to study this issue further because of a k-correction effect: the maximum contribution of an AGN mid-IR power law is beyond the IRAC bands at such high redshifts, so only the hottest dust AGNs could be manifested as IR power-law sources in IRAC (Caputi 2013). But one of the $24\,\mu {\rm{m}}$ detected most massive galaxies in our sample is indeed an X-ray source (Civano et al. 2016) spectroscopically confirmed to be at z = 4.596, which indicates its AGN nature. Therefore, as a matter of precaution, we flagged the $24\,\mu {\rm{m}}$ detected sources in our sample and excluded all galaxies with ${\mathrm{log}}_{10}({M}^{* })\gt 10.8$ from further analysis. This high fraction of $24\,\mu {\rm{m}}$ detections among the most massive "Hα excess" galaxies supports our AGN hypothesis also based on the reversing trend discussed in Section 4.1, which shows that the $3.6\,\mu {\rm{m}}$ incidence becomes higher at the highest-mass end, after a continuous decrease from intermediate to high stellar masses.

4.2.2. The Bimodal Distribution of Hα Emitters on the SFR–M* Plane

Figure 2 shows that the "Hα excess" galaxies form two distinct clouds on the SFR–M* and sSFR–M* planes. We identify these clouds as the so-called main sequence of star formation (e.g., Noeske et al. 2007; Rodighiero et al. 2010) and a starburst cloud, which appears to be more prominent than what is observed at any lower redshift. The two-cloud separation becomes completely evident when we plot the fraction of "Hα excess" galaxies versus sSFR (Figure 3). From this figure we can empirically derive that main-sequence galaxies are those with ${\mathrm{log}}_{10}\mathrm{sSFR}\lesssim -8.05$, while starburst galaxies have ${\mathrm{log}}_{10}\mathrm{sSFR}\gtrsim -7.60$.

Zoom In Zoom Out Reset image size

In addition, there is an sSFR local minimum at $-8.05\lesssim \mathrm{sSFR}\lesssim -7.60$ that here we denominate the "star formation valley," which likely contains galaxies in transition between the two star formation modes. This is analogous to the so-called "green valley," which in the literature refers to galaxies in transition from the main-sequence to the passive regime (e.g., Balogh et al. 2011; Zehavi et al. 20011; Renzini & Peng 2015).

On the sSFR–M* plot in Figure 2 we indicate upper limits to the median sSFR per stellar mass bin, obtained from all the galaxies at $3.9\leqslant z\leqslant 4.9$. For the "Hα excess" galaxies we adopted the derived sSFR, while for the other galaxies we considered sSFR upper limits derived from the line complex EW (and corresponding Hα flux) upper limits. In that plot we clearly see that the upper limits to the median sSFR decrease with increasing stellar mass in the stellar mass range $9.2\leqslant {\mathrm{log}}_{10}({M}^{* })\leqslant 10.8$. This trend is consistent with what is observed at lower z (e.g., Karim et al. 2011).

We fitted the SFR–M* relation in the two clouds formed by the "Hα excess" galaxies in Figure 2 (left) with simple linear regressions and obtained

$$\begin{eqnarray}&&{\mathrm{log}}_{10}(\mathrm{SFR})=(0.69\pm 0.01)\times {\mathrm{log}}_{10}({M}^{* })-(5.44\pm 0.13)\end{eqnarray} \tag{ 2 }$$

for the main sequence and

$$\begin{eqnarray}&&{\mathrm{log}}_{10}(\mathrm{SFR})=(0.89\pm 0.02)\times {{log}}_{10}({M}^{* })-({6.18}_{-0.15}^{+0.16})\end{eqnarray} \tag{ 3 }$$

for the starburst cloud. In these relations, we have obtained the error bars on the slopes and intercepts through Markov Chain Monte Carlo fittings, assuming a crude 30% average error on both the SFR and stellar mass.

The "Hα excess" galaxies identified as starbursts here constitute all the starbursts in our entire SMUVS sample at $3.9\leqslant z\leqslant 4.9$, as virtually no SFR upper limit (corresponding to the non-"Hα excess" galaxies) lies on the starburst cloud in Figure 2. Instead, the main sequence defined by the "Hα excess" galaxies and fitted by Equation (2) is not the complete main sequence at those redshifts. Therefore, we cannot use Equation (2) along with analogous relations from the literature at lower redshifts to draw conclusions regarding the redshift evolution of the galaxy main sequence. Nonetheless, to put our results in context, in Figure 4 we show our separate main-sequence and starburst cloud fittings at $z=3.9\mbox{--}4.9$, along with the median main-sequence determination by Speagle et al. (2014) at similar redshifts, which is given by ${\mathrm{log}}_{10}(\mathrm{SFR})\approx 0.80\times {\mathrm{log}}_{10}({M}^{* })-6.36$ at $z\sim 4.3$. For reference, we also show some recent main-sequence determinations at $2\lt z\lt 3$ (Speagle et al. 2014; Whitaker et al. 2014; Bisigello et al. 2017). Among these literature results, only the Bisigello et al. (2017) curve has been obtained after performing a starburst/main-sequence separation as we do here. The others fit all star-forming galaxies onto a single sequence.

Zoom In Zoom Out Reset image size

The analysis of the different curves in Figure 4 tells us the following: first, the Speagle et al. (2014) main-sequence curve at $z\sim 4.3$ is above our own, even when ours is an upper limit to the main sequence that is biased for fitting only the "Hα excess" galaxies. This is simply because Speagle et al. (2014) do not separate main-sequence and starburst galaxies as we do, but rather fit a single relation for all star-forming galaxies. Second, our "Hα excess" galaxy main sequence broadly coincides with the main sequences by Speagle et al. (2014) and Whitaker et al. (2014) at $z\sim 2$. We caution the reader against a wrong interpretation: this coincidence is just fortuitous, as our own main sequence is an upper limit and the Speagle et al. (2014) and Whitaker et al. (2014) main sequences include all star-forming galaxies, so any direct comparison could be misleading. Our main-sequence determination can be more fairly compared with that of Bisigello et al. (2017) at $z\sim 2.5$, which has been obtained after a starburst segregation. Taking into account that our own main sequence at $z=3.9\mbox{--}4.9$ is an upper limit to the true main sequence at these redshifts, we can conclude that there is probably little or no evolution of this sequence between $z\sim 2.5$ and $z=3.9\mbox{--}4.9$.

4.2.3. The Origin of the sSFR Bimodality

A priori it may seem surprising that the bimodality displayed by the Hα emitters on the SFR–M* and sSFR–M* planes is not observed in the rest EW versus stellar mass diagram (Figure 1, left). Fumagalli et al. (2012) derived that the sSFR of an Hα emitter is proportional to the ratio between the line rest EW and the R-band mass-to-light ratio, i.e., $\mathrm{sSFR}\propto \mathrm{EW}/({M}^{* }/{L}_{R})$, where ${L}_{R}$ is the R-band luminosity. Therefore, one could reasonably expect that the distribution of ${M}^{* }/{L}_{R}$ ratios is different for main-sequence and starburst Hα emitters and that this difference is to some extent responsible for the bimodal sSFR behavior.

Figure 5 shows the R-band mass-to-light ${M}^{* }/{L}_{R}$ ratio versus sSFR for our Hα emitters with stellar mass $9.2\leqslant {\mathrm{log}}_{10}({M}^{* })\leqslant 10.8$. The R-band luminosities ${L}_{R}$ considered here correspond to continuum luminosities and have been obtained from the best-fit templates in LEPHARE's run excluding the $3.6\,\mu {\rm{m}}$ photometry. Although there is no bimodality in the ${M}^{* }/{L}_{R}$ ratios, a decreasing trend of these values with increasing sSFR is evident, with starburst galaxies having a tight distribution of small ${M}^{* }/{L}_{R}$ values. For the starbursts, we find a median ${M}^{* }/{L}_{R}=0.12$, with 16th–84th percentiles of 0.09–0.18. Instead, main-sequence Hα emitters typically have larger ${M}^{* }/{L}_{R}$, following a wide distribution: the median is ${M}^{* }/{L}_{R}=0.31$, and the 16th–84th percentiles are 0.17–0.45.

Zoom In Zoom Out Reset image size

The small ${M}^{* }/{L}_{R}$ ratios for starbursts are mainly the consequence of their systematic young ages, which typically are a few × 10⁷ yr, according to their best SED-fitting templates. However, starbursts are not the only galaxies in our sample that are characterized by such young ages: about 20% of the "non-Hα excess" galaxies with stellar mass $9.2\leqslant {\mathrm{log}}_{10}({M}^{* })\leqslant 10.8$ at $3.9\leqslant z\leqslant 4.9$ have equally young best-fit ages in their SED fitting. Besides, starbursts are also characterized by having the largest extinction values among the "Hα excess" galaxies at fixed stellar mass, but there is no bimodality in the color excess $E(B-V)$ distribution, as can be seen in Figure 6. All these facts indicate that the sSFR bimodality for the Hα emitters does not arise as a consequence of a single galaxy property, but is rather the consequence of the existence of a galaxy population (i.e., the starbursts) with a particular combination of properties: large Hα EW, young ages, and mostly high dust extinctions.

Zoom In Zoom Out Reset image size

Although dust extinction and age can be degenerate in the SED fitting, this does not appear to be a significant problem for our starburst galaxies. If we analyze the probability density distribution as a function of age, marginalized over all other variables, we find that ∼90% of the starbursts have 84th percentiles at ages <10⁸ yr, i.e., they are truly galaxies dominated by young stellar populations, in contrast to the main-sequence Hα emitters, whose ages are mostly >10⁸ yr.

4.2.4. The Importance of Starburst Galaxies at z ∼ 4–5

The existence of different regions on the SFR–M* plane, corresponding to different modes of star formation, has been analyzed in the literature at different redshifts, from the local universe (Renzini & Peng 2015) to $z\sim 3$ (e.g., Santini et al. 2009; Kajisawa et al. 2010; Wuyts et al. 2011; Ilbert et al. 2015). Here we show, for the first time, the existence of a prominent starburst sequence along with the main sequence on the SFR–M* plane at $3.9\leqslant z\leqslant 4.9$.

Starbursts constitute a small fraction (15%) of all the SMUVS galaxies at $3.9\leqslant z\leqslant 4.9$, but the percentage that we find here is significantly higher than the percentages found at $z\sim 2\mbox{--}3$ (e.g., Rodighiero et al. 2011; Sargent et al. 2012; Schreiber et al. 2015). This is the case in spite of defining starbursts with an sSFR cut that is quite higher than the value ${\mathrm{log}}_{10}(\mathrm{sSFR})\gtrsim -8.1$ adopted by Rodighiero et al. (2011) at $z\sim 2$. At these lower redshifts, the percentage of galaxies with ${\mathrm{log}}_{10}({M}^{* })\gt 9$ and ${\mathrm{log}}_{10}(\mathrm{sSFR})\gt -7.60$ is negligible.

Starbursts and main-sequence galaxies are similarly important in number among our "Hα excess" galaxies with $9.2\leqslant {\mathrm{log}}_{10}({M}^{* })\leqslant 10.8$ at $3.9\leqslant z\leqslant 4.9$. However, since starbursts make only $\sim 15 \%$ of all the SMUVS galaxies in these stellar mass and redshift bins, the median sSFR upper limits shown in Figure 2 (right) lie all on the star formation main sequence.

Figure 7 shows how the fraction of starburst galaxies at $3.9\leqslant z\leqslant 4.9$ varies with stellar mass. These fractions range from ∼0.25 at ${\mathrm{log}}_{10}({M}^{* })\sim 9.3$ to $\lt 0.05$ at ${\mathrm{log}}_{10}({M}^{* })\gt 10.2$, showing that the starburst phenomenon is much more common among intermediate-mass than massive galaxies at $z\sim 4\mbox{--}5$, consistently to what is predicted by theoretical galaxy models (e.g., Cowley et al. 2017) and observed at lower redshifts (e.g., Bisigello et al. 2017).

Zoom In Zoom Out Reset image size

It is interesting to do a back-of-the-envelope calculation of how much stellar mass these starbursts could accumulate at these redshifts. A ${\mathrm{log}}_{10}(\mathrm{sSFR})=-7.2$ value implies a stellar-mass doubling time of $1/\mathrm{sSFR}\approx 1.6\times {10}^{7}\,\mathrm{yr}=0.016\,\mathrm{Gyr}$. The elapsed time between z = 4.9 and z = 3.9 is 0.379 Gyr. Besides, from Figure 7 we see that the fraction of galaxies with ${\mathrm{log}}_{10}({M}^{* })\sim 9.5$ that are starbursts is 0.20. So if we assume that all galaxies with similar stellar mass would spend a similar amount of time in the starburst phase, then the ${\mathrm{log}}_{10}({M}^{* })\sim 9.5$ galaxies would spend one-fifth of the above-mentioned elapsed time as starbursts, i.e., $\sim 0.076\,\mathrm{Gyr}$, which is almost five times the stellar-mass doubling time. Therefore, assuming a recycled fraction of 50%, we infer that these galaxies could grow their stellar mass by a factor of ∼2.5 at those redshifts. This means that a typical ${M}^{* }\sim 3\times {10}^{9}\,{M}_{\odot }$ galaxy at $3.9\leqslant z\leqslant 4.9$ would become a galaxy with ${M}^{* }\sim 8\times {10}^{9}\,{M}_{\odot }$ at $z\lt 3.9$. This crude estimate shows that the large sSFR values derived here are consistent with a rapid, but plausible, growth of intermediate-mass galaxies happening $\sim 1.5\,\mathrm{Gyr}$ after the big bang.

4.3. The Inferred Cosmic SFRD at z = 3.9–4.9

We can use our derived SFR to obtain an estimate of the cosmic SFRD at $3.9\leqslant z\leqslant 4.9$. Considering only the "Hα excess" galaxies gives us a lower limit to the SFRD. Instead, if we added the SFR contributions of all the other SMUVS galaxies at these redshifts, for which we only have SFR upper limits, we could estimate an upper limit to the SFRD. However, we recall that SMUVS significantly loses completeness at stellar masses ${\mathrm{log}}_{10}({M}^{* })\lt 9.2$, so even if we apply completeness corrections, we can only extrapolate the missing galaxy SFRs based on the galaxies that we detect. Therefore, estimating a more secure upper limit requires accounting somehow for the galaxies that we do not see. This can be done, at least in a crude manner, taking into account the faint-end slope of the galaxy stellar mass function (GSMF) at those redshifts and assuming that the SFR versus stellar mass trends shown in Figure 2 can be extrapolated to lower stellar masses.

Figure 8 shows the redshift evolution of the SFRD in the so-called Lilly–Madau diagram (Lilly et al. 1996; Madau et al. 1996). In this plot we have included our derived lower and upper limits to the SFRD at a median redshift $\langle z\rangle =4.29$, as well as a compilation of values recently reported in the literature, which are based on different galaxy surveys and individual galaxy SFRs from a variety of SFR tracers: UV fluxes (Bouwens et al. 2015), spectral line emission from narrowband surveys (Sobral et al. 2014), far-IR fluxes (Gruppioni et al. 2013), and a combination of UV and IR fluxes (Kajisawa et al. 2010; Burgarella et al. 2013; Dunlop et al. 2017). Multiple other previous works have presented SFRD estimates at high z, especially based on UV fluxes. However, UV-based SFRs are much more affected by dust extinction corrections than those based on any other tracer, making the resulting SFRD values particularly uncertain at least to $z\sim 6$ (see, e.g., Castellano et al. 2014).

Zoom In Zoom Out Reset image size

In the derivation of our own SFRD, we have explicitly excluded all galaxies with ${\mathrm{log}}_{10}({M}^{* })\gt 10.8$ to minimize any plausible AGN contamination (see discussion in Section 4.2). By considering only the SFR of "Hα excess" galaxies, we obtained a robust SFRD lower limit of $0.066\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{Mpc}}^{-3}$. To obtain an SFRD upper limit, we estimated separately the contributions of all galaxies with $9.2\leqslant {\mathrm{log}}_{10}{M}^{* }\leqslant 10.8$ and $8\leqslant {\mathrm{log}}_{10}{M}^{* }\lt 9.2$ and added them up, as follows. For the $9.2\leqslant {\mathrm{log}}_{10}{M}^{* }\leqslant 10.8$ objects we directly adopted our derived SFR ("Hα excess" galaxies) or SFR upper limits (for the other galaxies). To account for the contribution of lower-mass galaxies, we considered that the extrapolation of the GSMF at z = 4–5 (Caputi et al. 2015; Grazian et al. 2015) indicates that we should expect ∼5 times more galaxies with $8\leqslant {\mathrm{log}}_{10}{M}^{* }\lt 9.2$ than with ${\mathrm{log}}_{10}{M}^{* }\geqslant 9.2$. And according to Figure 2, the $8\leqslant {\mathrm{log}}_{10}{{\rm{M}}}^{* }\lt 9.2$ objects should have on average SFR that are about 1 dex lower than the more massive galaxies. Thus, the lower-mass galaxies should roughly add $\sim 50 \%$ to the SFRD value calculated based only on the $9.2\leqslant {\mathrm{log}}_{10}{M}^{* }\leqslant 10.8$ sources. With all these considerations, we obtain an SFRD upper limit of $0.106\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{Mpc}}^{-3}$. All the SFRD quoted here and in Figure 8 correspond to a Salpeter (1955) IMF over stellar masses $(0.1\mbox{--}100)\,{M}_{\odot }$ to facilitate the comparison with the literature, in which this IMF is mostly used. Our resulting range of SFRD at $\langle z\rangle =4.29$ is in good agreement with the other recent observational determinations in the literature.

The different curves in Figure 8 show best fits to compilations of observational data (Behroozi et al. 2013; Madau & Dickinson 2014) and theoretical predictions (Dayal et al. 2014; Khaire & Srianand 2015). We see that the majority of the most recent SFRD estimates at high z, including our own, are at least 0.10–0.15 dex above the best fits provided by Madau & Dickinson (2014) and Behroozi et al. (2013). This is probably because most of the observational works considered for these fits calculated SFRs based on UV fluxes, which are very sensitive to dust extinction corrections, as mentioned above. A revision of these corrections (Bouwens et al. 2015) or, even better, a consideration of UV+IR data (Kajisawa et al. 2010; Burgarella et al. 2013) or SFRs based on line emission like those we obtain here should provide better constraints to the cosmic SFRDs at high z at least to $z\sim 6$. Dust extinction corrections are expected to be less critical at earlier cosmic times, but this will have to be confirmed in the future, with systematic IR galaxy surveys that can trace Balmer line emission in the early universe, which will happen with the James Webb Space Telescope (JWST). The theoretical predictions appear to be in good agreement with our SFRD lower limit at $\langle z\rangle =4.29$, as well as some of the other observational data points.

The contribution of starburst galaxies to the total SFRD budget up to $z\sim 3$ has been discussed in the literature (e.g., Rodighiero et al. 2011; Sargent et al. 2012; Schreiber et al. 2015). The derived contributions vary according to the selection effects and the starburst definition adopted by different groups, but in all cases it has been found that starbursts make $\lesssim 15 \%$ of the total SFRD. As we discussed before, all these studies have been based on the analysis of ${M}^{* }\gtrsim {10}^{10}\,{M}_{\odot }$ galaxies. Here, we obtain that the "Hα excess" galaxies defined as starbursts by ${\mathrm{log}}_{10}(\mathrm{sSFR})\gtrsim -7.60$ account for $\gt 50 \%$ of our upper limit to the SFRD (and 84% of our lower limit based only on the "Hα excess" galaxies) at $3.9\leqslant z\leqslant 4.9$. These percentages are substantially higher than all previous determinations. The key reason for this higher percentage is that here we are including galaxies down to lower stellar masses, which are much more numerous than massive galaxies and have a higher fraction of starbursts among them.

4.4. A >50σ Overdensity of 3.9 ≤ z ≤ 4.9 Galaxies

Within our analysis of SMUVS galaxies at $3.9\leqslant z\leqslant 4.9$, we find that six of these sources lie in a very small region of only $\sim 0.20\times 0.20$ arcmin², with a median redshift $z\approx 4.2$. We have a total of 5925 SMUVS galaxies at those redshifts over a total area of 0.66 deg², which makes for an average surface density of ∼2.5 arcmin⁻². Therefore, the six galaxies concentrated in the tiny $\sim 0.20\times 0.20$ arcmin² area constitute a remarkable $\gt 50\sigma$ overdensity of sources at $3.9\leqslant z\leqslant 4.9$.

Figure 9 shows the UltraVISTA HKs stack and SMUVS $3.6\,\mu {\rm{m}}$ images of this galaxy overdensity. The properties of the six individual galaxies, as well as those of another galaxy at a similar (albeit less secure) redshift, which lies only ∼0.3 arcmin apart, are given in Table 1.

Zoom In Zoom Out Reset image size

Table 1.  Properties of Individual Sources in the $z\approx 4.2$ Galaxy Overdensity

ID	R.A. (J2000)	Decl. (J2000)	${z}_{\mathrm{phot}}$	${\mathrm{log}}_{10}\,{M}^{* }\,({M}_{\odot })$	SFR $({M}_{\odot }\,{\mathrm{yr}}^{-1})$	ID (L+16)	${z}_{\mathrm{phot}}$ (L+16)
SM2_103948	10:00:37.43	+02:25:40.25	4.83 ± 0.10	${9.71}_{-0.08}^{+0.08}$	36.9 ± 7.2	⋯	⋯
SM2_103968	10:00:37.55	+02:25:40.97	4.57 ± 0.10	${10.22}_{-0.50}^{+0.08}$	8.6 ± 8.5	⋯	⋯
SM2_104007	10:00:37.73	+02:25:40.37	4.07 ± 0.10	${9.80}_{-0.09}^{+0.08}$	53.3 ± 13.9	743559	4.05
SM2_104067	10:00:37.91	+02:25:44.92	4.22 ± 0.10	${9.06}_{-0.08}^{+0.08}$	<9.8	⋯	⋯
SM2_104073	10:00:37.81	+02:25:44.37	3.93 ± 0.10	${9.02}_{-0.08}^{+0.08}$	83.5 ± 9.1	744162	4.10
SM2_104100	10:00:37.62	+02:25:46.00	3.97 ± 0.10	${8.77}_{-0.09}^{+0.08}$	28.2 ± 3.2	744344	4.13
SM2_104510	10:00:37.10	+02:25:57.71	${4.26}_{-3.81}^{+0.15}$	${7.79}_{-0.10}^{+0.11}$	$\lt 0.5$	746450	0.36

Note. Source coordinates correspond to the HK_s mosaic. The SFR values are based on the inferred Hα rest EWs and fluxes derived from the $3.6\,\mu {\rm{m}}$ photometry. The last two columns contain ID and photometric redshifts from the COSMOS 2015 catalog by Laigle et al. (2016; L+16). The last row corresponds to a galaxy outside the main overdensity region, but only $\sim 0.3\,\mathrm{arcmin}$ apart from it.

Download table as:  ASCIITypeset image

In the main overdensity region, three out of six galaxies have redshifts $3.93\leqslant z\leqslant 4.07$. Within the error bars, this suggests that these galaxies are likely part of the same bound structure. These three galaxies are also present in the COSMOS 2015 catalog by Laigle et al. (2016), whose photometry and redshifts have been independently determined based on a previous release of the UltraVISTA data (DR2) and shallower Spitzer data (see Table 1). Their photometric redshifts are consistent with ours within the error bars. These independent results support our claim that a $z\sim 4$ overdensity is present in that region of the sky. The other three galaxies that we find in the overdensity region are at higher redshifts $4.22\leqslant z\leqslant 4.83$, so spectroscopic data are needed to determine whether they belong to the same structure or not.

Among the six galaxies in the overdensity, five are classified as "Hα excess" sources, although for one of them this classification is marginal. The other four have SFRs corresponding to several tens ${M}_{\odot }\,{\mathrm{yr}}^{-1}$. Note, however, that only one of them (ID SM2_104073) is a starburst galaxy, while the others are main-sequence galaxies. All these galaxies have intermediate stellar masses, with only one of them having ${\mathrm{log}}_{10}{M}^{* }\gt 10$.

In this paper we have studied a sample of 5925 galaxies at $3.9\leqslant z\leqslant 4.9$ from the SMUVS survey over 0.66 deg² of the COSMOS field. We have analyzed the presence of flux excess in the IRAC $3.6\,\mu {\rm{m}}$ band to identify "Hα excess" galaxies at these redshifts. From the inferred Hα EW we have derived SFR and sSFR for these galaxies and obtain upper limits for all the other galaxies at these redshifts.

The incidence of "Hα excess" galaxies in our $3.9\leqslant z\leqslant 4.9$ sample decreases from stellar masses ${\mathrm{log}}_{10}({M}^{* })\sim 9.2$ through ${\mathrm{log}}_{10}({M}^{* })\sim 10.8$, indicating that intense star formation is relatively more important among intermediate-mass than massive galaxies at these redshifts. At higher stellar masses, surprisingly, we see a reversing trend that seems difficult to explain by simply invoking enhanced star formation for the most massive systems. As a large fraction of the most massive galaxies with $3.6\,\mu {\rm{m}}$ flux excess at $3.9\leqslant z\leqslant 4.9$ are $24\,\mu {\rm{m}}$ detected, we suggest that AGN activity could at least be partly responsible for this flux excess. This statement is supported by a spectroscopically confirmed AGN among the most massive "Hα excess" galaxies.

The most important result in this paper is that we found that the "Hα excess" galaxies form two distinct sequences on the SFR–M* plane, which we recognize as the so-called star formation main sequence and a starburst sequence. Although the exact location of galaxies on this plane depends on the SED-fitting assumptions, we have shown that the sSFR bimodality is independent of the assumed metallicities and dust reddening law, under reasonable assumptions. Starburst galaxies are characterized by a unique combination of large Hα EW, young ages, and mostly high dust extinctions, which result in very high sSFR.

This suggests the existence of two distinct star formation modes, one corresponding to secular baryonic accretion within dark matter halos (giving rise to the main sequence), and another one much more effective for rapid galaxy growth, possibly linked to galaxy mergers or external perturbations that produce the starburst phase. A fundamental intrinsic assumption in all these results is that the same conversion from gas into stars given by Kennicutt's law (Equation (1)) is valid for all star-forming galaxies at high z. Demonstrating this general validity is extremely difficult, but recent studies of the link between gas content and star formation in high-z galaxies suggest that a universal conversion may indeed hold (e.g., Scoville et al. 2016).

The star formation main sequence at different redshifts has been analyzed from multiple observational data sets (e.g., Brinchmann et al. 2004; Noeske et al. 2007; Rodighiero et al. 2010; Elbaz et al. 2011; Speagle et al. 2014; Whitaker et al. 2014; Tasca et al. 2015). Theoretical studies have attempted to explain the main sequence as a mere consequence of gas accretion within dark matter halos (e.g., Somerville et al. 2008; Lu et al. 2014; Cousin et al. 2015), but have difficulty in reproducing the evolution of its normalization with cosmic time (Dutton et al. 2010; Genel et al. 2014; Furlong et al. 2015; Sparre et al. 2015)—see Mitchell et al. (2014) for a thorough discussion of this problem. Other groups have investigated the role of galaxy–black hole coevolution in shaping the star formation main sequence (e.g., Mancuso et al. 2016; Kaviraj et al. 2017). In any case, a complete theoretical explanation of this relation between stellar mass and instantaneous SFR is still missing.

In contrast to the widely discussed main sequence, most observational studies have failed to recognize the presence of a separate starburst cloud on the SFR–M* plane (its existence has been suggested by Cassarà et al. [2016], but these authors referred to it as a "second main sequence"). The failure to recognize the starburst cloud is mainly due to the relatively small galaxy samples analyzed in most of the literature, with starburst galaxies typically identified and studied only among massive galaxies (Rodighiero et al. 2011; Sargent et al. 2012). Here we show, for the first time, that starburst galaxies make a clearly distinct sequence on the SFR–M* plane and constitute a significant percentage ($\sim 15 \%$) of all galaxies with $9.2\leqslant {\mathrm{log}}_{10}{M}^{* }\leqslant 10.8$ at $3.9\leqslant z\leqslant 4.9$. Only a slightly smaller percentage is obtained at $z\sim 2\mbox{--}3$ when a broad dynamic range in stellar mass is considered, albeit with significantly lower sSFR (Bisigello et al. 2017). We find that the fraction of starburst galaxies has a strong dependence on stellar mass, varying from ∼0.25 at ${\mathrm{log}}_{10}({M}^{* })\sim 9.3$ to $\lt 0.05$ at ${\mathrm{log}}_{10}({M}^{* })\gt 10.2$. This is also expected from galaxy formation models (e.g., Vogelsberger et al. 2015).

The starburst galaxies in our sample are characterized by ${\mathrm{log}}_{10}(\mathrm{sSFR})\gt -7.60$. At redshifts $z\lesssim 4$, these very high sSFR values are very unusual and have only been found to be common among low stellar mass (${M}^{* }\lesssim {10}^{8}\,{M}_{\odot }$) young galaxies (see, e.g., Figure 11 in Karman et al. 2017; see also Amorín et al. 2017; Vanzella et al. 2017). Here we show that these high sSFRs are relatively common among star-forming, intermediate-mass galaxies at $3.9\leqslant z\leqslant 4.9$.

Another key result of this paper is that starbursts make for at least 50% of the total SFRD budget at $z\sim 4\mbox{--}5$. This percentage is substantially higher than any previous determination found in the literature. This is because previous starburst studies only analyzed galaxies with ${M}^{* }\gtrsim {10}^{10}\,{M}_{\odot }$ (up to $z\sim 3$). The inclusion of lower-mass galaxies is essential to recognize the importance of the starburst phase and how much it contributes to the total cosmic SFRD, particularly at high redshifts.

We have also discovered an unusually high significance galaxy overdensity at these high redshifts. A total of six galaxies in our high-z sample reside in a tiny area of $0.20\times 0.20$ arcmin², five of which correspond to "Hα excess" sources. The spectroscopic confirmation of such an overdensity would unveil one of the most concentrated active sites of star formation known at $z\sim 4\mbox{--}5$.

All the results presented in this paper have exploited the unique combination of area and depth provided by the SMUVS survey at mid-IR wavelengths. This survey has provided us with an unprecedented level of statistics and dynamic range, which are fundamental to revealing unknown aspects of galaxy evolution in the young universe.

Based in part on observations carried out with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. Also based on data products from observations conducted with ESO Telescopes at the Paranal Observatory under ESO program ID 179.A-2005 and on data products produced by TERAPIX and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium. Also based on observations carried out by NASA/ESA Hubble Space Telescope, obtained and archived at the Space Telescope Science Institute, and the Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. This research has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA.

We thank an anonymous referee for a constructive report. K.I.C., S.D., and W.C. acknowledge funding from the European Research Council through the award of the Consolidator Grant ID 681627-BUILDUP.

Facilities: Spitzer - Spitzer Space Telescope satellite, VISTA - , Subaru - Subaru Telescope.

Software: SExtractor, IRAF, LePhare.

In this paper we have modeled the SEDs of our $3.9\leqslant z\leqslant 4.9$ galaxies assuming a solar metallicity. Other authors have preferred to adopt a much lower metallicity ($0.2\,{{\rm{Z}}}_{\odot }$) to study "Hα excess" galaxies at similar redshifts (e.g., Smit et al. 2016). However, in a stellar-mass-selected galaxy sample like SMUVS, the adoption of subsolar metallicities for all galaxies is likely not correct (e.g., Scoville et al. 2016). The most likely scenario is that our galaxies span a range of metallicities between a fraction of the solar and the entire solar values.

To test the impact of this effect on our results, we have rerun the SED fitting of our "Hα excess" galaxies at $3.9\leqslant z\leqslant 4.9$ allowing LePhare to use both templates with $0.2\,{Z}_{\odot }$ and ${Z}_{\odot }$ metallicities, covering the same grid for all the other parameter values (i.e., SFHs, ages, extinctions). By doing this, we found that (i) 96% of all our original "Hα excess" galaxies are confirmed to be in the redshift range $3.9\leqslant z\leqslant 4.9$ and thus confirmed as Hα emitting sources, and (ii) about 34% of the original "Hα excess" galaxies have a best-fit SED with $0.2\,{Z}_{\odot }$, while $\sim 62 \%$ have a best-fit SED with solar metallicity.

Taking into account these results, we recomputed the Hα rest EW and fluxes for the "Hα excess" galaxies that have a best-fit SED with $0.2\,{Z}_{\odot }$ in the new LePhare run. At this low metallicity, the net Hα contribution to the total $3.6\,\mu {\rm{m}}$ flux excess is given by (Anders & Fritze-van Alvensleben 2003) $f({\rm{H}}\alpha )=0.81f({\rm{H}}\alpha +[{\rm{N}}\,{\rm{II}}]+[{\rm{S}}\,{\rm{II}}])$ and $f({\rm{H}}\alpha )=0.92f({\rm{H}}\alpha +[{\rm{N}}\,{\rm{II}}])$, which are valid at $3.9\leqslant z\leqslant 4.8$ and $4.8\lt z\leqslant 4.9$, respectively. Finally, we derived the corresponding clean Hα luminosities and SFR in the same way as explained in Section 4.2.

Figure 10 (left) shows the resulting SFR–M* plane for "Hα excess" galaxies when considering two metallicities for the SED fitting (we only plotted here the 96% of the original "Hα excess" sources that stayed in the $3.9\leqslant z\leqslant 4.9$ redshift range in the new SED-fitting run). Qualitatively, this plot looks similar to the analogous plot in Figure 2, with two distinct regimes (the galaxy star formation main sequence and starburst cloud) clearly visible. Figure 10 (right) shows the corresponding sSFR distribution compared to the original distribution derived adopting only solar metallicity (see Figure 3). This plot confirms the existence of two sSFR regimes, separated by a star formation valley corresponding to an underdensity of sources with $-8.05\lesssim \mathrm{sSFR}\lesssim -7.60$. Therefore, we can conclude that all the main results in this paper are robust against our SED metallicity assumptions.

Zoom In Zoom Out Reset image size

Another parameter that could potentially influence our results is the choice of dust reddening law in the SED fitting and recovery of intrinsic Hα luminosities. In this work we have adopted the Calzetti et al. (2000) reddening law for the SED fitting and applied the derived color excess from the best-fit SED of each "Hα excess" galaxy to determine its intrinsic Hα flux. Here we investigate the impact of a different assumption, namely, that the internal extinction of each galaxy is directly related to its UV spectral slope β, which is defined as the slope in ${f}_{\lambda }(\lambda )\propto {\lambda }^{\beta }$ for rest-frame $1500\lt \lambda \lt 2500\,\mathring{\rm A}$.

We computed the UV slope β of each of our galaxies by fitting this power-law functional form to the photometric fluxes tracing rest $1500\lt \lambda \lt 2500\,\mathring{\rm A}$, which correspond to about 10 photometric bands in the COSMOS field (as has been done by, e.g., Fudamoto et al. 2017). Then we assumed that the internal dust extinction was given by the Meurer et al. (1999) law, ${A}_{1600}=4.43+1.99\beta$, and considered a conversion ${A}_{6563}\approx {A}_{1600}/3.2$. We applied this new extinction correction (which we considered to be the same for the continuum and line) to each galaxy, and we recomputed its Hα rest EW and luminosity and derived the corresponding SFR.

Figure 11 is the analog of Figure 10, but in this case we show the SFR–M* plane and sSFR distribution that result from considering the new dust extinction correction (at fixed solar metallicity, as in Section 4). The star formation main sequence and starburst cloud are clearly distinct in these plots, showing that the sSFR bimodality is still present when we consider a single, β-dependent extinction relation rather than the best-fit extinction obtained from the individual SED fitting of each of our galaxies. We conclude that the observed sSFR bimodality is robust to reasonable changes in the modeling assumptions.

Zoom In Zoom Out Reset image size