GOODS-HERSCHEL: STAR FORMATION, DUST ATTENUATION, AND THE FIR–RADIO CORRELATION ON THE MAIN SEQUENCE OF STAR-FORMING GALAXIES UP TO z ≃ 4^*

The galaxy sample we use in this work is drawn from a K_S-band selected multi-wavelength catalog in the GOODS-North field, spanning 19 passbands from GALEX near-ultraviolet (NUV) to IRAC 8 μm (namely, GALEX NUV, KPNO Mosaic U, Subaru Suprime-Cam B-V-IA624-R-I-z-Y, Canada–France–Hawaii Telescope (CFHT) WIRCam J-H-Ks, Subaru MOIRCS J-H-Ks and Spitzer IRAC 3.6-4.5-5.8 and 8 μm). The GALEX data are part of the GALEX GR6 data release.²⁶ The optical, NIR, and IRAC data are described in Capak et al. (2004), Ouchi et al. (2009), Wang et al. (2010), Kajisawa et al. (2011), L. Lin et al. (2012, 2015, in preparation), M. Dickinson et al. (2015, in preparation), and M. Onodera et al. (2015, in preparation).

We first register all images to a common pixel grid and then use SExtractor (Bertin & Arnouts 1996) in dual image mode to measure photometry. The K_S-band CFHT WIRCAM image has been adopted as the primary detection image because it represents the best compromise, among all available bands, between the need for a robust tracer of galaxy stellar mass and an angular resolution $(\simeq 0\buildrel{\prime\prime}\over{.} 8)$ matching that of most of the other bands, which simplifies catalog assembly and photometry measurements. The whole catalog contains 56,144 objects over the CFHT WIRCAM K_Simage field of approximately 900 arcmin² and down to an AB magnitude of 24.5 (i.e., the image 5σ limiting magnitude, see Wang et al. 2010 for details).

The images used here have very different resolutions. Rather than convolving all images to the worse resolution, which would result in a significant loss of information, we account for this in the estimate of aperture colors by applying point-spread function (PSF)-matching corrections based on the growth-curve of point-like sources. To limit uncertainties in such corrections we use $2\prime\prime$ diameter apertures to sample the galaxy SEDs.

Band-merging with IRAC photometry was done with a cross-match of the NUV–optical–NIR catalog to the K_S–IRAC catalog released by Wang et al. (2010). The latter study used a "real-cleaning" procedure to perform PSF modeling in the IRAC images by assuming as priors the positions of sources detected in the same K_SCFHT WIRCAM image we have used. We refer the reader to Wang et al. (2010) for a more detailed description and assessment of the extracted IRAC photometry. Here we want to stress that because we use the same K_Simage that was publicly released by Wang and collaborators, our cross-matching procedure becomes a very robust way of associating IRAC fluxes to our K_S-based multi-wavelength catalog.

The whole catalog contains 56,144 objects over the CFHT WIRCAM K_Simage field of approximately 900 arcmin² and down to an AB magnitude of 24.5 (i.e., the image 5σ limiting magnitude, see Wang et al. 2010 for details).

We classified 2,072 objects as stars, and as such exclude them from the sample according to their SExtractor stellarity index, at bright (K_S < 20) mag, and to their position in the BzK diagram, as in Daddi et al. (2004).

The morphological selection of stars is also effective in removing those Type 1 (i.e., optically unobscured) active galactic nuclei (AGNs) whose emission from the active nucleus dominates over the host galaxy in almost all bands. In any case, these objects are extremely difficult to deal with in terms of both the photometric redshift derivation and stellar mass estimate (e.g., Salvato et al. 2009; Cisternas et al. 2011; Bongiorno et al. 2012), so it is preferable to remove them from our sample. However, we do not attempt to remove from the star-forming galaxy sample those X-ray sources detected in the 2-Ms Chandra data (Alexander et al. 2003; Bauer et al. 2004) whose optical-to-NIR emission is dominated by the host galaxy. A number of recent studies (see, e.g., Mullaney et al. 2012; Santini et al. 2012; Juneau et al. 2013; Rosario et al. 2013) have shown that these AGNs are mostly Type II obscured AGN and LINERS, which have similar SFR to normal (i.e., non-AGN) star-forming and passive galaxies, respectively. Thus their inclusion has no net effect on our results. In Section 3 we discuss our approach to minimize the contribution of the AGNs present in our sample to the IR budget.

Photometric redshifts were estimated by running EAZY²⁷ (Brammer et al. 2008) on the multi-wavelength catalog. We used EAZY in its standard setup and modeled the observed galaxy SED with a linear combination of the seven standard EAZY galaxy templates (Whitaker et al. 2011) to maximize the likelihood as a function of the galaxy redshift. Following a common procedure for producing accurate photometric redshifts (see, e.g., Capak et al. 2007; Gabasch et al. 2008), we allowed for a photometric offset in each measured band. For each band in our catalog, we computed a photometric zero-point shift by iteratively running the code on a subsample of high fidelity spectroscopic redshifts. For each iteration and for each band, we computed the median of the ratio between the measured flux density in the catalog and the carefully computed model value (i.e., the integral of the best-fit solution template through the theoretical response curve), and applied these median offsets before starting the next iteration. Computing median offsets, instead of mean ones, is an effective way to filter out catastrophic events, such as badly measured SEDs or wrong spectroscopic redshifts determinations. The procedure described here usually converges to a robust solution within a limited number of iterations.

The photometric offsets between observed and modeled data may be attributed to many different effects, including actual zero-point errors, differences of the actual system response curves with respect to adopted ones, and inaccurate PSF-matching corrections, as well as uncertainties and limitations in the template library adopted.

In the following sections, we use a multi-wavelength catalog that has been corrected for these systematic offsets. In all bands these amount to less than 20% of the measured flux so that, while substantially improving the photometric redshift accuracy, applying these flux corrections does not have any major effect on the derived galaxy properties, nor in general on the main results of this study. We will discuss in the relevant sections when notable differences arise from the use of the corrected versus uncorrected photometry.

When comparing to the spectroscopic sample of Barger et al. (2008) and D. Stern et al. (2015, in preparation), we reach a relative (i.e., ${\rm{\Delta }}z=({z}_{\mathrm{phot}}-{z}_{\mathrm{spec}})/(1+{z}_{\mathrm{spec}})$) accuracy of 3% (see Figure 1), with less than 3% catastrophic outliers (i.e., objects with Δ z > 0.2).

Zoom In Zoom Out Reset image size

Stellar masses were determined with FAST²⁸ (Kriek et al. 2009) on the U to 4.5 μm PSF-matched aperture photometry, using Bruzual & Charlot (2003) delayed exponentially declining SFHs ($\psi (t)\propto \displaystyle \frac{t}{{\tau }^{2}}\mathrm{exp}(-t/\tau )$) with 0.01 < τ < 10 Gyr, solar metallicities (${Z}_{\odot }$ = 0.02), Salpeter IMF, and the Calzetti et al. (2000) reddening law with A_V up to 4 mag.

The choice of a fixed solar metallicity has essentially no impact on the derived stellar masses, because the age–metallicity degeneracy leaves mass-to-light ratios basically unchanged. Moreover, solar metallicities are still a fair assumption for our sample at high redshift given that, due to the increasing mass completeness with redshift, it contains at all cosmic epochs the most massive and most metal-rich systems.

Derived (aperture) masses were extrapolated to "total" masses using the ratio between the total (FLUX_AUTO) and aperture flux in the Ks-band detection image. While this approach corrects for the bulk of the flux loss, it is based on one band only and thus neglects any color gradient within the galaxy.

With respect to the choice of the SFH, we note that it has been shown that other forms of SFHs are more appropriate for star-forming galaxies at high redshift (Pannella et al. 2009a; Maraston et al. 2010; Papovich et al. 2011). For what directly concerns this work, rising or constant (possibly truncated) SFHs would change the stellar masses of our star-forming galaxy sample by an amount that is well within the median estimated accuracy of about 0.2 dex (see, e.g., Papovich et al. 2006, 2011; Buat et al. 2014; de Barros et al. 2014; and the results presented in Appendix C).

As a further accuracy test, we compared stellar mass estimates derived with FAST with the SED-fitting code described in Drory et al. (2004, 2009). The results from runs with the different codes and different adopted SFHs show no systematic differences in the stellar mass estimates, except for the highest stellar masses (${M}_{*}\gtrsim 2\times {10}^{11}\;{M}_{\odot }$), where the two-component fitting from Drory et al. (2004) predicts masses that are larger by 0.2 dex, with a global scatter of about 0.2 dex. We provide more details on the estimated stellar mass errors and on the comparison between the two codes in Appendix C.

The main aim of this study is to characterize the relation between star formation, UV dust attenuation, and the stellar mass content of galaxies over cosmic time. To do so, it is important not to mix galaxies that might still have some residual ongoing star formation (but are substantially quenched), with galaxies that are genuinely star-forming and reddened by dust attenuation. A selection based on direct SFR indicators, such as UV or line emission, inevitably mixes the two populations. Ideally, one would like to remove galaxies with extremely low specific star-formation rates (SSFRs) from the sample, but this is not easily achieved for individual galaxies because it requires accurate dust attenuation corrections or, alternatively, FIR/radio derived star-formation rate estimates. Following established procedures, we excluded quiescent galaxies from the sample by using the U − V versus V − J rest-frame color plot (Wuyts et al. 2007; Williams et al. 2009). This approach builds on two main ingredients: (1) the intrinsic rest-frame U − V color is a good proxy for the SSFR of a galaxy (e.g., Salim et al. 2005; Pannella et al. 2009b); (2) the color–color plot allows us to break the age-dust degeneracy and hence efficiently split the sample into two classes of galaxy, which are either essentially quiescent or still actively forming stars. We selected star-forming galaxies at all redshifts as

$$\begin{eqnarray*}&&(U-V)\lt 1.3,(V-J)\gt 1.6,(U-V)\lt (V-J)\\ &&\quad \times 0.88+0.59.\end{eqnarray*}$$

These UVJ selection limits were originally defined by Williams et al. (2009) in order to maximize the difference in SSFRs between regions. However, the rest-frame color distributions might be slightly different than Williams et al. (2009) due to photometric coverage, band selection, and the specific redshift where the analysis is performed (the rest-frame color of quiescent galaxies is becoming redder with decreasing redshift). As a result, different studies have used slightly different dividing lines, sometimes changing with redshift (e.g., Cardamone et al. 2010; Brammer et al. 2011; Whitaker et al. 2011; Muzzin et al. 2013; Strazzullo et al. 2013; Viero et al. 2013). The differences between these studies are always lower than 0.2 mag in color and often comparable with the same accuracy to which the rest-frame colors are known. In order to be conservative and minimize the contamination of quiescent galaxies in our sample, we chose to use the same selection region at all redshifts. We also checked whether slightly changing our assumptions on the selection limits (i.e., applying ±0.1 mag shifts on the rest-frame colors) would quantitatively change our results.

We checked our selection against the BzK selection of star-forming galaxies at z ≃ 2 (Daddi et al. 2004) (see Figure 2), and the results are in good agreement over the common redshift range. Finally, we tested whether all the results of this paper would remain the same if instead of using the UVJ selection to separate the active from passive/quenched galaxies we used a cut in SSFR, namely selecting as star-forming galaxies all the objects with SSFR ≥ 10⁻¹¹ yr⁻¹ as derived from the FAST SED-fitting output.

Zoom In Zoom Out Reset image size

Hereafter, we will concentrate on the UVJ-selected subsample of 14,483 star-forming galaxies that fall within the GOODS-Herschel area, the deepest FIR imaging available in the northern sky (Elbaz et al. 2011).

Following Brammer et al. (2011), rest-frame magnitudes are computed with EAZY for all objects in the catalog from the best-fit SED model integrated through the theoretical filters. The β slope of the UV continuum was calculated as

$$\begin{eqnarray}&&\beta =\displaystyle \frac{\ \mathrm{log}({f}_{{\lambda }_{1}}/{f}_{{\lambda }_{2}})}{\ \mathrm{log}({\lambda }_{1}/{\lambda }_{2})}=\displaystyle \frac{0.4({m}_{2}-{m}_{1})}{\ \mathrm{log}({\lambda }_{1}/{\lambda }_{2})}-2,\end{eqnarray} \tag{ 1 }$$

where m₁ and m₂ are the magnitudes at wavelengths λ₁ and λ₂, respectively (see, e.g., Overzier et al. 2011; Nordon et al. 2013). In this study, we adopt the GALEX FUV (λ_c ≃ 1530 Å) and NUV (λ_c ≃ 2315 Å) response curves to compute the photometric measure of β, but we also checked that the results obtained would remain largely unchanged by using different rest-frame bands sampling the slope, or by obtaining a direct fit of the form ${f}_{\lambda }\propto {\lambda }^{\beta }$ to the best-fit SED in the rest-frame range 1250 and 2600 Å. In Figure 3 we show the results of this comparison using the best-fit SEDs of FAST. The GALEX-derived estimate of β provides an accurate description of the UV slope, with a median offset of 0.1 and a scatter of around 0.2. This result is in good agreement with the finding reported in Finkelstein et al. (2012).

Zoom In Zoom Out Reset image size

An important issue, which is often not given careful consideration in high redshift studies dealing with star-forming galaxies, is the assessment of the sample's mass completeness. This is because the amount of dust attenuation plays an important role in the selection function: low mass galaxies (and the definition of low depends of course on the selection band and detection limit of the specific data set) can only be detected if they suffer from low dust attenuation. In other words, as redshift increases, magnitude-limited samples become more and more biased against objects with lower masses and/or high M_*/L ratios, such as galaxies that are not currently forming stars or are highly dust obscured. Because we are interested in dissecting dust attenuation properties, we want to make sure that this issue has no impact on our results. Following Rodighiero et al. (2010), we estimated our completeness mass by using a stellar population with a 1 Gyr old constant SFH and highly attenuated ($E[B-V]$ = 0.6). As a function of redshift, there is a one-to-one correspondence between the observed Ks magnitude and the stellar mass built up by this model. We hence define as the completeness mass at a given redshift the mass corresponding to the completeness magnitude (≃23.8 AB mag) of our Ks selected catalog. All galaxies more massive than our completeness mass will have a brighter flux at Ks and hence be detected in our survey; the same is true for galaxies suffering less dust attenuation. Obviously, this definition of mass completeness is dependent on the specific stellar population parameters adopted, and can only apply to galaxies that can be likely represented by such an SFH (i.e., star-forming galaxies). Figure A.1 of Rodighiero et al. (2010) illustrates how the actual mass completeness varies by adopting different ages and attenuations for the stellar population models. However we note that our modeling is very conservative and the quoted mass completenesses (see Table 1) should be regarded as safe and robust, although theoretical, values. As previously stated, our sample, which is composed of a priori selected star-forming galaxies, is also cut in SSFR (at approximately 10⁻¹¹ yr⁻¹, see discussion above), meaning that, even if we are complete down to a certain mass and able to robustly trace the SFR-M* made by the bulk of the star-forming galaxy population, we would still be able to detect galaxies with extremely low SSFR if they suffer low dust attenuation. This can be translated in an SFR limit at the completeness mass that we also quote in Table 1. We want to stress here that, as our primary selection is a mass selection, the secondary limits in SFR have no influence on our results because they apply to galaxies well below the MS locus.

Table 1.  Stellar Mass and Star Formation Rate Completeness vs. Redshift

z	0.7	1	1.3	1.7	2.3	3.3
log M*/${M}_{\odot }$	9.0	9.5	10.0	10.3	10.5	10.7
SFR (${M}_{\odot }$ yr−1)	1	2	10	20	30	90

Download table as:  ASCIITypeset image

The main goal of this work is to define, in a consistent way, the locus of the MS of star-forming galaxies over a significant range of cosmic time. Many studies have been published on this topic in recent years (e.g., Rodighiero et al. 2010; Karim et al. 2011). We use data probing the full FIR SED (using all Herschel bands), so as not to rely on the ${L}_{\mathrm{IR}}$ derived from the bolometric correction of a single band. In Figure 4 we show the MS of star-forming galaxies (i.e., the correlation between the galaxy stellar mass and SSFR) in different redshift bins. Only results from mass-complete samples are shown. We fit a linear trend in the first three redshift bins and find a slope consistent with −0.2 ± 0.08, corresponding to a slope of 0.8 ± 0.08 for the proper logM_*–logSFR correlation, up to redshift z ≃ 1.5. This value of the slope is consistent with a number of previous studies (e.g., Santini et al. 2009; Rodighiero et al. 2011, 2014; Lin et al. 2012; Heinis et al. 2014; Renzini & Peng 2015). We then assume the same slope for the higher redshift bins, and fit only the normalization of the MS up to z ≃ 4.

Zoom In Zoom Out Reset image size

In Figure 5 we show the evolution of the MS normalization for a stellar mass of log(M_*/${M}_{\odot }$) = 10.5 as obtained from this work, together with some previously published results (Brinchmann et al. 2004; Daddi et al. 2007, 2009; Elbaz et al. 2007; Noeske et al. 2007; Pannella et al. 2009a; Magdis et al. 2010) and the analytical descriptions of the MS evolution from Pannella et al. (2009a) and Elbaz et al. (2011).²⁹

Zoom In Zoom Out Reset image size

Our median SSFR values are slightly lower than previous studies at z ≥ 1.5. In principle this could be explained by the fact that we are dealing with a mass-complete sample, while previous studies were biased toward brigther (${K}_{{\rm{S}}}\lt 23$ in Pannella et al. 2009a) or star-formation selected (24 μm detected sources in Elbaz et al. 2011) samples that almost by construction are bound to overshoot the median value of the whole star-forming population. At the same time we cannot exclude the possibility that at least part of the offset could be due to some misclassified galaxies that entered the UVJ-selected star-forming sample.

Finally, it is worth noticing how the SSFR tends to keep increasing up to the highest redshift without any evidence of the possible flattening at z ≥ 2 obtained in previous studies (e.g., Daddi et al. 2009; Stark et al. 2009; González et al. 2010), which were based on SFR estimates derived from dust-attenuation-corrected UV luminosities. We will discuss this in more detail in the following sections, but based on our results this might be related to a change in the dust and/or stellar population properties of star-forming galaxies with cosmic time.

Taking advantage of the deep 1.4 GHz map available in the GOODS-N field (Morrison et al. 2010), we perform a stacking analysis of the radio continuum in the same bins of mass and redshift as done in the FIR Herschel bands in order to study the evolution of the radio–FIR correlation for a uniformly selected, mass-complete sample of star-forming galaxies over cosmic time. Sargent et al. (2010) discussed the impact of sample selection, specifically the FIR versus radio selection, in explaining the controversial results obtained in literature studies on the evolution of the FIR–radio correlation. Based on Spitzer 24 μm-derived SFR estimates, they showed how the correlation was not evolving up to redshift 1 and possibly even at higher redshift, although this was plagued by the uncertain extrapolation from the measured 24 μm to total IR luminosity at z ≥ 1.5. Here we have the advantage of being able to stack mass-complete samples.

Some recent Herschel-FIR based studies have claimed a redshift evolution of the correlation (e.g., Ivison et al. 2010; Casey et al. 2012; Magnelli et al. 2012, but see also, e.g., Strazzullo et al. 2010; Bourne et al. 2011; Barger et al. 2012, 2014; Del Moro et al. 2013 for different conclusions), with increasing radio luminosities for the same IR luminosity at higher redshift.

On the other hand, simple theoretical arguments predict the opposite trend (e.g., Condon 1992; Carilli et al. 2008), because of the increased inverse Compton cooling of the relativistic electron population due to the scattering off of the increasing cosmic microwave background at high redshift, ${U}_{\mathrm{CMB}}\propto {(1+z)}^{4}$.

We derive 1.4 GHz monochromatic luminosities from the measured radio flux densities assuming a synchrotron spectral index, α, of –0.8 and the median redshift of each stacked subsample, as

$$\begin{eqnarray}{L}_{1.4} & = & 1.19\times {10}^{14}{\mathrm{DM}}^{2}\times {\text{}}{S}_{1.4}\\ & & \times {(1+{z}_{\mathrm{med}})}^{-(\alpha +1)}\ ({\rm{W}}{\mathrm{Hz}}^{-1}),\end{eqnarray} \tag{ 4 }$$

where DM is the distance in Mpc at the median redshift z_med of the stacked subsample, and S_1.4 is the measured flux density in the stack in units of μJy. We notice that our assumption of a nonevolving spectral index is justified by most observational studies dealing with statistical samples of star-forming galaxies (see e.g., Ibar et al. 2009; Bourne et al. 2011). We adopt the Yun et al. (2001) relation to convert such luminosities in star-formation rates:

$$\begin{eqnarray}&&{\text{}}{\mathrm{SFR}}_{1.4}=5.9\times {10}^{-22}\;{\text{}}{\text{}}{L}_{1.4}\ ({M}_{\odot }\ {\mathrm{yr}}^{-1}).\end{eqnarray} \tag{ 5 }$$

In order to study the evolution of the correlation, we compute, following Helou et al. (1985), the luminosity ratio:

$$\begin{eqnarray}&&{q}_{\mathrm{IR}}=\mathrm{log}({L}_{\mathrm{IR}}/(3.75\times {10}^{12}\;{\rm{W}}))-\mathrm{log}({L}_{1.4}/{\rm{W}}{\mathrm{Hz}}^{-1})\end{eqnarray} \tag{ 6 }$$

for all redshift and mass bins where we have performed our stacking analysis. We show the results in Figure 6 together with the ±1σ confidence locus of the local universe q_IR value (e.g., Yun et al. 2001; Bell 2003). We find striking agreement with the local correlation at all redshifts and masses explored.

Zoom In Zoom Out Reset image size

We stress that our result is the first so far obtained for a mass-complete uniformly selected sample of star-forming galaxies and hence is expected to be significantly more robust against selection bias, as compared to all previous analyses.

This finding has some non-trivial implications. First, it shows that a FIR–radio correlation is well defined up to z ≃ 4 with basically the same q_IR value measured in the local universe, and thus radio continuum data are an ideal tool for estimating dust-unbiased SFRs at high redshift.

Second, since radio luminosity is only contributed by young stellar populations, our result suggests that the IR emission of star-forming galaxies in the redshift range 0.5 < z < 4 can only have, in a statistical sense, a marginal contribution from old stellar populations, and likely smaller than the 20% estimated in the local universe (see, e.g., Law et al. 2011).

In order to study the redshift evolution of the UV dust attenuation, we define as an effective measurement of dust attenuation the quantity:

$$\begin{eqnarray}&&{A}_{\mathrm{UV}}=2.5\times \mathrm{log}({\mathrm{SFR}}_{\mathrm{IR}}/{\mathrm{SFR}}_{\mathrm{UV}}^{\mathrm{obs}}+1),\end{eqnarray} \tag{ 7 }$$

(i.e., the magnitude correction to the measured UV emission). This means that, by construction, we obtain

$$\begin{eqnarray*}&&{{SFR}}_{\mathrm{UV}}^{\mathrm{corr}}={{SFR}}_{\mathrm{total}}={{SFR}}_{\mathrm{IR}}+{{SFR}}_{\mathrm{UV}}^{\mathrm{obs}}.\end{eqnarray*}$$

For each subsample stacked in each redshift and mass bin, we estimate the median FUV emission and UV spectral slope (see Equation (1)).

Assuming the Calzetti et al. (2000) empirical correlation between UV dust attenuation and slope β:

$$\begin{eqnarray}&&{A}_{\mathrm{UV}}^{\beta }=4.85+2.31\times \beta ,\end{eqnarray} \tag{ 8 }$$

we then define the quantity

$$\begin{eqnarray}&&{\mathrm{SFR}}_{{\mathrm{UV}}_{}^{\beta }}={\mathrm{SFR}}_{\mathrm{total}}^{\beta }={\mathrm{SFR}}_{\mathrm{UV}}^{\mathrm{obs}}\times {10}^{{A}_{\mathrm{UV}}^{\beta }/2.5}\end{eqnarray} \tag{ 9 }$$

which provides an estimate of the total SFR by correcting the observed UV luminosity assuming the measured β slope and the Calzetti et al. (2000) law, which is a standard prescription adopted for high redshift studies.

In the following we discuss how the measured UV dust attenuation, A_UV, depends on galaxy properties, such as the stellar mass and the UV spectral slope over the redshift range 0.5–4.

4.3.1.  ${A}_{{UV}}$ Versus M_*

In the top-left panel of Figure 7 we show the correlation between dust attenuation and galaxy stellar mass. The correlation evolves only marginally, less than 0.3 mag in A_UV, over the redshift range explored:

$$\begin{eqnarray}&&{A}_{\mathrm{UV}}=1.6\times \mathrm{log}{M}_{*}-13.5\;\mathrm{mag}\;\;\mathrm{for}\;{\text{}}z=1.2-4.\end{eqnarray} \tag{ 10 }$$

Zoom In Zoom Out Reset image size

The stellar mass content of a star-forming galaxy turns out to be a robust proxy of the dust attenuation affecting its UV emission (see, e.g., Pannella et al. 2009a, 2013; Schaerer & de Barros 2010; Buat et al. 2012; Kashino et al. 2013; Heinis et al. 2014; Oteo et al. 2014). Previous results focused on a relatively small redshift range, and sometimes on a UV-selected sample, which, by construction, tends to be biased against the most massive and most obscured star-forming galaxies. Here, we test this correlation over a wide redshift range and in a mass-complete way. Testing the real dispersion of a correlation obtained through stacking is not trivial and often impractical. To at least obtain an estimate of the correlation scatter, we used all Herschel detected sources in the field (Elbaz et al. 2011) up to z ≃ 1.3, where the GOODS-Herschel data are deep enough to allow a good statistical description of the parent sample, and estimate a dispersion of ≃1 mag in ${A}_{\mathrm{UV}}$.

The mild evolution of this correlation lends support to a number of earlier studies claiming that the same amount of star formation suffered less dust extinction at high redshift compared to the local universe (e.g., Buat et al. 2007; Reddy et al. 2012). Because of the evolution of the MS with redshift, the same amount of star formation is hosted in galaxies that are less and less massive as redshift increases, and hence they will suffer a correspondingly lower UV dust extinction.

On the other hand, dust attenuation is known to correlate with both the ISM metallicity and the mode and geometry of star formation in the host galaxy (i.e., normal versus starbursting galaxies; e.g., Heckman et al. 1998; Cortese et al. 2006). The ISM metallicity is expected to decrease with increasing redshift at a fixed mass, according to the well-studied evolution of the mass–metallicity relation (e.g., Tremonti et al. 2004; Savaglio et al. 2005; Erb et al. 2006; Zahid et al. 2013a; Steidel et al. 2014; Wuyts et al. 2014). By assuming that the ISM conditions of star-forming galaxies do not change with redshift, one would see an effective decrease with redshift of UV attenuation at a fixed stellar mass; if anything we measure a weak trend going in the opposite direction. This suggests that the ISM conditions of MS galaxies are evolving, becoming somehow more extreme with redshift, which might be understood as a combination of both geometrical effect (star-forming galaxies are becoming smaller with redshift) and physical effects (gas fractions and SSFRs are also increasing), leading to an overall enhanced volume density of the UV radiation field, as also suggested by their higher dust temperature (e.g., Hwang et al. 2010; Magdis et al. 2012; Symeonidis et al. 2013). We will come back to this topic in Section 5.

In the top-right panel of Figure 7 we plot the median UV luminosities for the same stacking samples as in the top-left panel. We show that at all redshifts explored, there is basically no correlation between stellar mass—a robust proxy of the ongoing star formation and UV dust attenuation—and the emerging UV photons. The amount of UV radiation that escapes dust absorption rises systematically with redshift, mimicking the rise of the general SFR level. This is consistent with a number of previous studies (e.g., Pannella et al. 2009a; Buat et al. 2012; Heinis et al. 2013). Here we explore a wider redshift range and find basically the same result: the emerging UV light is not (or only marginally in the lowest redshift bin) correlated to the actual SFR present in a galaxy, because the tight correlation between stellar mass (a proxy for the SFR; i.e., the intrinsic UV luminosity) and the dust attenuation conspires so that the emerging, average UV light is basically the same for very different SFR levels.

In the bottom-left panel of Figure 7 we plot the median values of the galaxy UV spectral slope, β, per bin of stellar mass at different redshifts. The two quantities, β and M_*, clearly correlate at all redshifts, but while the attenuation is fairly constant, or slightly increasing with redshift, β values are becoming systematically bluer with redshift. The fact that high mass galaxies at high redshift have a similar dust attenuation compared to similar mass galaxies at lower redshift, but have bluer UV slopes, has important implications for UV-derived SFRs in the high redshift universe. Using standard recipes to correct the observed UV emission by means of the β slope would systematically underpredict galaxy SFRs. This is better quantified in the bottom-right panel of Figure 7 where we compare the MS evolution plot derived from Herschel data (black squares, as in Figure 4) to the one derived from UV data that have been corrected by dust attenuation according to Equation (9). The comparison shows three main features mainly dependent on redshift: (1) at redshift lower than ≃1, the ${\mathrm{SFR}}_{{\mathrm{UV}}_{}^{\beta }}$ tends to overpredict the real SFR, possibly because old stellar populations present in star-forming galaxies contribute to create redder UV slopes, which is in agreement with previous results (see, e.g., Kong et al. 2004 and Overzier et al. 2011); (2) in the redshift range 1.5–2.5 the ${\mathrm{SFR}}_{{\mathrm{UV}}_{}^{\beta }}$ provides a nice match to the real SFR, which is in agreement with the results of Daddi et al. (2007), Pannella et al. (2009a), and more recently Rodighiero et al. (2014); (3) at redshift z ≥ 3, ${\mathrm{SFR}}_{{\mathrm{UV}}_{}^{\beta }}$ are systematically lower than the real SFR and underpredict the SFR by a factor of around two. Intriguingly enough, this underprediction of SFR conspires to create a plateau of the UV-derived SSFR values at z ≥ 2, which was indeed observed in previous studies (see, e.g., Daddi et al. 2009; Stark et al. 2009; González et al. 2010). Moreover, the underestimate of SFR has an impact on the estimates of the global star-formation rate density (SFRD) at z ≥ 3, which were obtained by applying standard recipes as in Equation (9), which should be revised upward by a factor of around two. A similar conclusion was already presented in de Barros et al. (2014) and Castellano et al. (2014) for a sample of low mass Lyman-break selected galaxies. We extend their result to the high mass end of the star-forming galaxy population at high redshift, putting forward a more compelling case for an upward revision of the SFRD estimates obtained so far at redshift z ≥ 3. We will further discuss this point in the next section.

4.3.2. A_UV Versus β Slope

In this section we explore in more detail the correlation between UV dust attenuation and UV spectral slope by repeating the stacking procedure in bins of β slope. This correlation is possibly the most widely used method in the literature to derive UV dust attenuation and hence SFR from UV emission. Daddi et al. (2004) calibrated it for a sample of star-forming galaxies at z ≃ 2 using the observed B − z color, which is a proxy for the UV spectral slope, and this calibration has been shown to behave well in a number of studies (e.g., Daddi et al. 2007; Pannella et al. 2009a; Rodighiero et al. 2014). Previously, Meurer et al. (1999) had shown that local starburst galaxies and z ≃ 3 Lyman Break Galaxies (LBGs) were following the same tight correlation between dust attenuation and UV spectral slope. After Meurer et al. (1999), numerous investigations have tried to understand whether the so-called Meurer law would be followed by galaxies at all redshifts. Different studies have reported contradicting results on this topic. One thing that became clear within a few years after the Meurer et al. (1999) study was that local normal star-forming galaxies do not follow the Meurer relation, but instead have "redder" spectra for the amount of dust attenuation they suffer (e.g., Kong et al. 2004; Buat et al. 2005; Overzier et al. 2011; Boquien et al. 2012 and Grasha et al. 2013).

We now investigate how the median UV spectral slope relates to dust attenuation. We plot at all redshifts the predicted correlations between dust attenuation and β slope according to Meurer et al. (1999), Calzetti et al. (2000), Daddi et al. (2004), and Overzier et al. (2011) in the top-left panel of Figure 8. The first three determinations overlap to a large extent at blue UV slopes, whereas they spread out at red slopes. The determination of Overzier stands out, clearly predicting redder slopes at a fixed attenuation. However, we remind the reader that the sample used in Overzier et al. (2011) was strictly based on blue star-forming objects with no data points having a UV slope redder than beta −0.5. For this reason, the derived relation (R. Overzier 2015, private communication) should not be used for very attenuated galaxies.

Zoom In Zoom Out Reset image size

When binned in β values, the median dust attenuation of the MS galaxy population follows the prediction of the Calzetti law already in our lowest redshift bin with a possible hint of a dust attenuation overprediction only in the reddest slope bins and up to redshift z ≃ 1. This latter feature might be due to old stellar populations of massive galaxies already contributing partially to the measured red slope. The departure from the Meurer relation of normal star-forming galaxies in the local universe has been thoroughly investigated in the last years and it is now well accepted to explain it as an age effect of the stellar populations hosted in local galaxies (see e.g., Kong et al. 2004; Boquien et al. 2012 and Grasha et al. 2013). In the lowest redshift bin we show, for illustrative purpose, the relation obtained by Boquien et al. (2012) for a sample of seven face-on and spatially resolved local star-forming galaxies drawn from the Herschel Reference Survey.

Already at z ≃ 1.3 and up to z ≤ 3 the correlation matches well the observed data points, agreeing with the results of Daddi et al. (2007), who found the UV slope derived SFR to be in good agreement with those derived from 24 μm and radio continuum data. In the last redshift bin, between redshift three and four, we see that UV dust attenuation tends to overshoot the prediction of the Calzetti law, so that galaxies tend to have bluer UV slope (an offset of about −0.5) compared to that predicted by the Calzetti law at a fixed amount of dust attenuation. These results are overall in very good agreement with the ones discussed in the previous section.

There are several possible explanations for such a finding. At face value, our result seems to support the suggestion of Maiolino et al. (2004) that high redshift galaxies have a grayer attenuation law than that predicted by the Calzetti law. More recently, Castellano et al. (2012, 2014), Alavi et al. (2014), and de Barros et al. (2014) have suggested that the bluer slopes can be explained by galaxy stellar populations with young ages and/or low metallicities. On the other hand, Buat et al. (2012) and Kriek & Conroy (2013) have suggested that the presence of a dust bump, the feature at 2175Å, in the spectra of high redshift galaxies may result in underestimated values for the UV slope. Finally, blue spectral slopes at high redshift have been reported previously by a number of authors as the possible outcome of a measurement bias when estimating slopes directly from the observed photometry rather than, as we actually do in this study, using the best-fit SED (see, e.g., Buat et al. 2011; Bouwens et al. 2012; Finkelstein et al. 2012).

We remind the reader that our β estimate is directly proportional to the difference between two rest-frame magnitudes, and as such, one might wonder about the possible impact of the shifts applied to the photometric catalog on our conclusions. The only sensible difference we find by using the uncorrected photometry is in the first redshift bin. The uncorrected photometry would produce redder slopes (${\rm{\Delta }}\beta \;\sim \;$ 0.3) that would possibly move our dust attenuation further away from the Meurer relation and possibly closer to the correlation found in the local universe (see e.g., Boquien et al. 2012). All the other redshift bins remain unchanged with possibly a weak trend, but well within the estimated β uncertainty, to obtain slightly bluer values compared to the ones we used. If anything, this would reinforce our conclusions.

In the top-right panel of Figure 8 we show the evolution of the correlation between β slope and total star-formation rate at different redshifts. In all panels the dotted line is the best fit to the correlation in the lowest redshift bin, and is kept constant at all redshifts to aid comparison. As redshift increases, the amount of star-formation rate at fixed UV slope increases steadily. This can be interpreted as further evidence that at a fixed star-formation rate there is less and less UV dust attenuation with increasing redshift.

In the bottom panel of Figure 8 we plot the median measured UV absolute magnitudes from stacks in bins of β. When binning in β values, again basically in UV dust attenuation, we find an anti-correlation between the measured UV luminosities and UV slopes: galaxies with redder UV slopes (i.e., that are more dust attenuated) tend to have fainter UV-measured luminosities, and galaxies that are less attenuated have brighter luminosities. Pannella et al. (2009a) found a similar result at z ≃ 1.7; here we can extend this at all probed redshifts (see also Heinis et al. 2013, 2014 for similar conclusions).

Overall, the anti-correlation is consistent with an identical slope at all redshifts, but with an increasing normalization, with the UV-observed luminosity increasing by about 2 mag (more than a factor of six in luminosity) over the redshift range probed.

Reddy et al. (2009) found that low luminosity objects are less attenuated, statistically speaking, than more UV luminous ones. This is contrary to what we found here. The difference between our result and the study of Reddy et al. (2009) might easily be explained by the different selection bands used to build the galaxy samples. We use a K band selection, which is very close to a stellar mass selection at these redshifts and is therefore sampling the mass function in a fairly complete way at the high mass end. Massive galaxies are intrinsically UV-bright but also extremely dust attenuated, so it is reasonable to expect these objects to be very faint in the UV and even fainter than lower mass galaxies. On the other hand, the LBG-like selection of the Reddy et al. (2009) sample, which is mostly an observed frame UV selection, will sample the mass function in a sparse way, giving more weight toward the low-to-intermediate range of the stellar mass function, thereby missing almost completely the high mass end as well as very obscured galaxies, which is a very well-known drawback of UV-selected samples. These UV-selected low mass galaxies preferentially have low star-formation rates, low dust attenuation, and faint UV output, and they outnumber by several orders of magnitude the high mass galaxies.

Hα line emission is a well-known and robust SFR tracer (e.g., Kennicutt 1998). Lying at 6563 Å, Hα suffers much less dust extinction compared to the UV rest-frame of a galaxy spectrum, and hence the uncertainty on SFR due to the poorly known dust attenuation can be an order of magnitude less severe as compared to UV-based estimates. On the other hand, Hα is mainly produced by extremely massive (≥10 ${M}_{\odot }$), short lived (≃10 Myr) stars that are deeply embedded in the giant molecular cloud (GMC) H ii regions, whereas stars producing the UV rest-frame stellar continuum are less massive (<10 ${M}_{\odot }$), shine over timescales 10 times longer, and have time to migrate out of the dense H ii regions. The net outcome of this process is that Hα emission suffers from an extra attenuation that has, as we just described, an origin driven by the distribution and density of H ii regions within the galaxy itself.

The extra amount of attenuation suffered by nebular emission has been quantified in the local universe by Calzetti et al. (2000) to be a factor of about 1.7. This can be derived by accounting for the fact that the stellar continuum suffers roughly half of the reddening suffered by the ionized gas, encoded into E(B-V)_star = 0.44 E(B-V)_gas, and for the fact that Calzetti et al. (2000) used two different extinction curves for the nebular (Fitzpatrick 1999) and continuum (Calzetti et al. 2000) emission,³⁰ which have similar shapes but different normalizations (i.e., R_V = 3.1 and 4.0, respectively) and thus for example:

$$\begin{eqnarray}&&\displaystyle \frac{{{A}_{{\rm{V}}}}_{\mathrm{gas}}}{{{A}_{{\rm{V}}}}_{\mathrm{star}}}=\displaystyle \frac{{{R}_{{\rm{V}}}}_{\mathrm{gas}}}{{{R}_{{\rm{V}}}}_{\mathrm{star}}}\displaystyle \frac{E{[B-V]}_{\mathrm{gas}}}{E{[B-V]}_{\mathrm{star}}}\sim 1.7.\end{eqnarray} \tag{ 11 }$$

Here we assume that this ratio is constant in wavelength. Although this is obviously not strictly true, because it assumes that the two mentioned extinction curves have exactly the same shape and only differ in their normalizations, it turns out not to be a bad approximation at the wavelengths of interest for this study, namely at the FUV (1500 Å) and Hα (6563 Å) wavelengths.

At higher redshift, things are much less constrained. Erb et al. (2006) and Reddy et al. (2010) claimed a vanishing diference (i.e., a factor close to 1) at $z\;\simeq \;2$, but this has not been confirmed by more recent studies (e.g., Förster Schreiber et al. 2009; Mancini et al. 2011; Wuyts et al. 2013).

The main problem for high redshift studies is the lack of a robust SFR indicator to directly compare with the Hα and UV stellar continuum measurements. Most of the quoted studies rely on indirect arguments by estimating the SFR from SED-fitting or from the β-corrected UV luminosity. This approach is, by construction, extremely uncertain and model dependent, and the end product of the comparison is closer to a correction factor, which produces in the end a consistent result, rather than a real physical measure.

In this study we are able to look at the redshift evolution of the extra attenuation affecting H ii regions in a way that is relatively model-independent. The idea is to compare the A_UV–M_* relation we found in Section 4.3.1 to the one that similarly relates ${A}_{{\rm{H}}\alpha }$ to M_* as a function of redshift. The latter relation was first established in the local universe by Garn & Best (2010), and more recently confirmed by Zahid et al. (2013b), by using the observed Balmer decrement (i.e., ${{\rm{H}}}_{\alpha }$/H_β) measured in the SDSS spectra, as previously done in Brinchmann et al. (2004).

In the last few years, several different studies have looked at the redshift evolution, at least up to z ≃ 1.5, of the ${A}_{{\rm{H}}\alpha }-\mathrm{log}$M_* relation (e.g., Sobral et al. 2012; Domínguez et al. 2013; Ibar et al. 2013; Kashino et al. 2013; Momcheva et al. 2013; Price et al. 2014). Compared to the SDSS studies mentioned, these are affected by additional uncertainties and systematics mainly due to the fact that the Hα line shifts into the observed NIR for $z\;\geqslant \;0.5$, which this often prevents a robust measurement of the Balmer decrement. Despite the diverse caveats, most of the results obtained so far, up to z ≃ 1.5, suggest that there is no redshift evolution of the ${A}_{{{\rm{H}}}_{\alpha }}-\mathrm{log}$M_* correlation (see, e.g., Figure 6 in Price et al. 2014, for a compilation of available results).

In order to investigate the nebular exinction, we start from the ${A}_{{{\rm{H}}}_{\alpha }}-\mathrm{log}$M_* correlation defined in Garn & Best (2010):

$$\begin{eqnarray}&&{A}_{{{\rm{H}}}_{\alpha }}=0.91+0.77\cdot X+0.11\cdot {X}^{2}-0.09\cdot {X}^{3}\end{eqnarray} \tag{ 12 }$$

where the quantity X = (log M_* − 10 + 0.23) is the stellar mass term in units of 10¹⁰ ${M}_{\odot }$, which has been rescaled from the Chabrier (2003) IMF used in Garn & Best (2010) to the Salpeter (1955) IMF we are using in this study. We then assume a Calzetti et al. (2000) attenuation law to derive the line reddening $E{[B-V]}_{\mathrm{gas}}$ as

$$\begin{eqnarray}&&E{[B-V]}_{\mathrm{gas}}={A}_{{\rm{H}}\alpha }/3.33\end{eqnarray} \tag{ 13 }$$

and the dust attenuation for nebular emission in the rest-frame UV as

$$\begin{eqnarray}&&{A}_{\mathrm{UV}-\mathrm{gas}}=10.4\;E{[B-V]}_{\mathrm{gas}}.\end{eqnarray} \tag{ 14 }$$

In Figure 9 we overplot the rescaled Garn & Best (2010) relation divided by different scaling factors (1.7, 1.3 and 1) to the measured ${A}_{\mathrm{UV}}$–log M_* relation at different redshifts. The factor 1.7 corresponds to the value found by Calzetti et al. (2000) in the local universe. We find that this factor is clearly smaller at higher redshift, being 1.3 at z ≃ 1, and becomes close to unity at higher z (i.e., there is a vanishing difference between the reddening of the H ii regions and that affecting the young massive, but slightly older stars responsible for the rest-frame continuum UV emission). This can be understood through the fact that, on the one hand, galaxies are becoming slightly smaller in physical size (e.g., Ferguson et al. 2004; Elbaz et al. 2011; M. Pannella et al. 2015, in preparation), but at the same time their SFR density is dramatically increasing from the local universe out to z ≃ 4, by more than a factor of 40 according to the evolution of the normalization of the log M_*–logSFR correlation. For this reason, the volume density of H ii regions is becoming so high that it almost entirely fills up the available galaxy surface. This suggests that the ISM of star-forming galaxies was more and more dense and opaque to UV radiation with increasing redshift, eventually resembling, at z ≥ 1.5, the physical properties of local H ii regions, as has already been suggested by Reddy et al. (2010). More recently Wild et al. (2011), by analysing a sample of 15,000 galaxies drawn from the Sloan Digital Sky Survey, have found a correlation between the galaxy SSFR and the ratio between line and continuum reddening, such that the ratio becomes smaller with increasing SSFR (i.e., as galaxies move toward the starburst regime locus). The correlations found by Wild et al. (2011) nicely explain our findings and also support our conclusions.

Zoom In Zoom Out Reset image size

On the other hand, Kreckel et al. (2013) have shown, for a sample of eight nearby spatially resolved star-forming galaxies, that the ratio between line and continuum reddening increases with increasing SFR density, which seems to be in contradiction with the Wild et al. (2011) result and hence our conclusions. Very recently, Boquien et al. (2015) obtained results on M33 that are in nice agreement with those presented in our study, and suggest possible causes explaining the controversial results of Kreckel et al. (2013).

A more physical way to look at ISM conditions in star-forming galaxies and their possible evolution with redshift is to consider the correlation between dust attenuation and ISM metal content for star-forming galaxies. In the local universe, the correlation was first explored by Heckman et al. (1998) for a sample of starburst galaxies, and found to yield a tight correlation for star-forming galaxies, with more dusty galaxies being more metal-rich. The correlation is expected, as metals and more complex molecules are needed to build up the actual dust that absorbs the UV radiation.

In a more recent study, Cortese et al. (2006) looked at the properties of normal star-forming galaxies in the local universe, which is a sample of objects much closer to what in this paper we are calling MS galaxies. They found that these objects were also following a tight correlation between dust attenuation and ISM metallicity, but the correlation was offset with respect to the one of starburst galaxies: at a fixed metallicity, starburst galaxies are more obscured compared to normal star-forming galaxies. The offset in the correlation can be interpreted in terms of the different ISM conditions found in normal star-forming and starburst galaxies. The latter host star-formation rate densities that are much higher than the former—both because star formation happens on smaller scales and because of the higher overall levels of star formation—and they can enrich most of their ISM much earlier and more quickly than normal smoothly star-forming MS galaxies.

Assembling mass-complete samples of galaxies at high redshift with measured metallicities is practically impossible with the presently available facilities. Most of the existing data in the high redshift universe faces the unavoidable bias of being limited by line sensitivity, and hence selected against the more dusty, more massive, but also more metal-rich systems. As already discussed, dust attenuation has an impact on the source detection itself, which becomes more and more severe with increasing redshift. This bias can in principle become extreme when it comes to metallicity measurements; where easiest, and possibly only, the objects for which it is possible to measure metallicity will be the ones with the lowest metallicities. This effect, which seems to nicely fit and explain some recent observations of metallicity and gas masses at high redshift (see, e.g., Ouchi et al. 2013; Tan et al. 2013; Troncoso et al. 2014), should be carefully taken into account when dealing with high redshift sources and deriving more general conclusions on the global galaxy population.

Here we estimate the median ISM metallicity for our samples of star-forming galaxies in bins of stellar mass and redshift, by applying the fundamental metallicity relation, FMR, found by Mannucci et al. (2010). The relation is defined as a tri-dimensional surface linking the SFR, metallicity, and stellar mass of star-forming galaxies. It was first defined in a sample of local SDSS galaxies and then found to be essentially identical at higher redshift at least out to z ≃ 2. The relation implies that galaxies obey an anti-correlation between SFR and metallicity: at a fixed stellar mass, the drop in metal content of star-forming galaxies at higher redshift (e.g., Tremonti et al. 2004; Savaglio et al. 2005; Erb et al. 2006; Zahid et al. 2013a) is compensated by the increase of SFR.

As already stated, the FMR has been verified to remain unchanged up to z ≃ 2 (but see also, e.g., Cullen et al. 2014; Steidel et al. 2014; Wuyts et al. 2014), whereas a sudden change has been claimed recently by the same authors in Troncoso et al. (2014) for a sample of z ≃ 3 galaxies. Of course, we have no way to test their claim here because of how the mass completeness of their sample could possibly affect their conclusions according to that discussed in the previous paragraph. Very recently, Maier et al. (2014) have shown that, when using the formalism put forward in Lilly et al. (2013) or equivalently the Equation (2) of Mannucci et al. (2010), as opposed to the Equation (5) of the same paper, to describe the FMR, measured high redshift metallicities are overall well described by the FMR predictions. Here we assume that the FMR is still valid in our last redshift bin at z ≥ 3.

We show in Figure 10 the evolution of dust attenuation (here shown as the classical quantity log(${L}_{\mathrm{IR}}$/L_UV) instead of A_UV as in the previous figures, to be consistent with the studies of Heckman et al. 1998 and Cortese et al. 2006) versus the ISM metallicity at different redshifts. We see a clear evolution in redshift, with MS galaxies moving progressively from the locus occupied locally by normal star-forming galaxies (solid line), to overshooting the one occupied locally by starburst galaxies (long-dashed line). We stress here that our result does not depend dramatically on the use of the FMR, but that a qualitatively similar outcome would have been obtained by simply using the evolution of the mass–metallicity relation. In other words, the shift in redshift we are seeing is mostly driven by the fact that galaxies of a fixed stellar mass are becoming less and less metal-rich with increasing time while keeping the same level of UV dust attenuation.

Zoom In Zoom Out Reset image size

Galaxies are becoming overall more compact in sizes at higher redshift and at the same time they are also significantly increasing their SFRs, and this trend makes them more and more similar to local starbursts in terms of ISM conditions. This similarity is also confirmed by the reported evolution of the IR SED shape of MS galaxies up to redshift $z\;\simeq \;2$ by Magdis et al. (2012) and more recently by Magnelli et al. (2014). It is extended at higher redshift by C. Schreiber et al. (2015, in preparation), with the median dust temperature becoming higher with increasing redshift and indeed more similar to the dust temperatures of highly star-forming galaxies seen in the local universe.