Free Access
Issue
A&A
Volume 643, November 2020
Article Number A4
Number of page(s) 13
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202038163
Published online 27 October 2020

© ESO 2020

1. Introduction

Over the past decades, extragalactic surveys have provided large observational data of galaxies, covering wide wavelength ranges from the rest-frame ultraviolet (UV) to the rest-frame far-infrared (FIR). These systematic, panchromatic observations have enabled the connection of the build-up of galaxies from the very early stages of their formation at high-redshift to the matured population in the local Universe. In particular, the cosmic star formation rate density (SFRD) is found to increase at a rapid rate until z ∼ 2 − 3, followed by a smooth and slow decline by an order of magnitude until the present day Universe (e.g., Wilkins et al. 2008; Madau & Dickinson 2014; Bouwens et al. 2015; Oesch et al. 2018). While great progress has been made to flesh out this overall picture, the measurement of the total SFRD at z ≳ 3 is still uncertain due to uncertain dust correction factors.

The star formation rates (SFRs) of z ≳ 3 galaxies are typically estimated using the UV emission from massive stars that have a short (∼100 Myr) lifetime. Several studies provide empirical relations to estimate the SFR directly from the UV continuum or from emission lines by atoms ionised by the UV emission (e.g., Kennicutt 1998; Kennicutt & Evans 2012; Madau & Dickinson 2014; Wilkins et al. 2019). As the UV emission is highly sensitive to dust attenuation, it needs to be corrected to estimate the total intrinsic star formation activity (Calzetti et al. 2000; Salim & Narayanan 2020). Thus, an understanding of the dust attenuation properties as a function of redshift and other galaxy properties is one of the most important ingredients to obtain accurate SFRs, and thus to gain an accurate census of galaxy build up across cosmic history.

Absorbed UV photons heat the dust grains, which in turn re-emit the energy as thermal emission at far-infrared (FIR) wavelengths. To correct for the absorbed UV emission, several empirical relations have been established between the dust attenuation and its FIR re-emission. Of particular importance is the relation between the infrared excess (IRX = LIR/LUV) and the UV spectral slope (β: fλ ∝ λβ). This was calibrated using local starburst galaxies (e.g., Meurer et al. 1999, hereafter M99) and is routinely used in the literature to estimate the dust attenuation from high-redshift galaxies, based on the measured UV spectral colors alone.

The relation between IRX and stellar mass (M) is another tool that relates the dust attenuation and stellar masses of galaxies. The stellar mass of galaxies reflects their past star formation activity, which in turn is responsible for producing dust particles. Therefore, the stellar mass is expected to correlate with the dust content of the interstellar medium (ISM), as shown by observations (e.g., Heinis et al. 2014; Santini et al. 2014; Pannella et al. 2015; Bouwens et al. 2016; Whitaker et al. 2017; Fudamoto et al. 2020; Álvarez-Márquez et al. 2016, 2019), and suggested by simulations (e.g., Graziani et al. 2020).

Both the IRX–β and the IRX–M relations are well studied from z ∼ 0 to z ≲ 4, over which most authors find no significant evolution for ensemble averages (e.g., Heinis et al. 2014; Pannella et al. 2015; Bouwens et al. 2016; Koprowski et al. 2018, 2020; Álvarez-Márquez et al. 2016, 2019; Fudamoto et al. 2017, 2020). However, at z ≳ 5, these relations are still uncertain as a statistically significant sample is yet to be obtained. While a few individual galaxies have been detected with luminous dust continuum emission even out to the highest redshifts (e.g., Watson et al. 2015; Laporte et al. 2017; Bowler et al. 2018; Hashimoto et al. 2019; Tamura et al. 2019), the general trend from population averages points to rapid changes of the relation to lower IRX values at z >  5 (e.g., Capak et al. 2015; Bouwens et al. 2016; Smit et al. 2018). In particular, using a sample of only 10 main-sequence galaxies, Capak et al. (2015) and Barisic et al. (2017) report a redshift evolution of the relation at z >  5, which indicates that the IRX of z >  5 galaxies is ∼1 dex lower than expected from the lower-redshift IRX–β and IRX–M relations. If correct, it implies that the use of the “classical” correction would overestimate the true IR luminosity, and hence also the total SFR. At the same time, the rapid decrease of the IR emission implies that at z >  4 the fraction of star formation ongoing in obscured environments (i.e., the obscured fraction of star formation) becomes smaller relative to un-obscured star formation. In previous studies, the fraction of obscured star formation as a function of stellar mass was found to be un-changed over the redshift range between 0 <  z <  2.5 (Whitaker et al. 2017), whereas the fainter IR emission found by Capak et al. (2015) suggests a decrease of the obscured fraction in z >  5 star forming galaxies (see also, Wilkins et al. 2018, for theoretical predictions at z >  6). However, the existing studies at z >  4 − 5 are based only on a handful (∼10) of sources, and the results need to be confirmed with a larger sample size. This is the goal of the present paper.

To examine the dust attenuation properties of z ≳ 5 using a large sample of normal star forming galaxies, we have studied the IRX–β, the IRX–M relation, and the obscured fraction of star forming galaxies using data from ALMA Large Program to INvestigate CII at Early Times (ALPINE, Le Fèvre et al. 2020; Béthermin et al. 2020; Faisst et al. 2020). ALPINE is a ∼70 h program to observe the [CII]158 μm emission lines and dust continuum of 118 normal star forming galaxies at z ∼ 4.4 − 5.8.

The survey provides the largest sample of normal star forming galaxies at 4 <  z <  6, which we here use to study the dust attenuation from a comparison of the IR and UV emission, and to examine how it relates with observable and derived galaxy properties.

This paper is organised as follows: in Sect. 2 we describe our sample and observations, and in Sect. 3 we presents the methods of our analyses. Section 4 shows the results on the IRX–β/M relations and the obscured fraction of star formation obtained from our sample. Finally, we conclude our study in Sect. 5. Throughout this paper, we assume a cosmology with (Ωm, Ωλ, h) = (0.3, 0.7, 0.7), and the Chabrier (Chabrier 2003) initial mass function (IMF) where applicable.

2. Observations

In this section, we briefly discuss the ALPINE galaxy sample and observations that are relevant to this study. We refer to Le Fèvre et al. (2020), Béthermin et al. (2020), and Faisst et al. (2020) for a complete description of the survey objectives, the ALMA data processing, and the multiwavelength ancillary observations, respectively.

2.1. Sample and ancillary data

ALPINE observed 118 star forming galaxies at z ∼ 4.4 − 5.9 (see Le Fèvre et al. 2020). The targets have secure spectroscopic redshifts from two observation campaigns in the COSMOS (105) and the GOODS-South (13) fields using rest-frame UV emission and/or absorption lines (Le Fèvre et al. 2015; Hasinger et al. 2018).

Using the vast amount of ancillary data available for these fields, we performed SED fitting to estimate SFRs, stellar masses, UV luminosities, UV spectral slopes, and other quantities, as described in detail by Faisst et al. (2020). We adopt these SED-based quantities from Faisst et al. (2020) for this paper.

Two galaxies in the sample are confirmed AGN based on an X-ray detection (DEIMOS_COSMOS_845652; Faisst et al. 2020) and from deep optical spectroscopy (CANDELS_GOODS_14; Grazian et al. 2020). As AGNs could outshine the rest-UV emission of stars, we removed these two sources from further analysis. Additionally, a stack of all the Chandra X-ray images of the remaining galaxies showed no detection, confirming that our sample is not significantly AGN-contaminated, on average. While heavily dust-obscured AGNs are not excluded by these rest-UV and X-ray criteria above, they should have little impact on the FIR emission probed by our ALMA data which is dominated by relatively cold dust.

By selection, the galaxies from the ALPINE survey have SFRs and stellar masses, which are consistent with the main-sequence of star formation at z ∼ 5 (see Fig. 1; Steinhardt et al. 2014; Schreiber et al. 2015; Faisst et al. 2020). The ranges of stellar mass and total SFR spanned by our targets are log (M/M)∼8.5 − 11.5 and log (SFR/M yr−1) ∼ 0.5 − 2.5, respectively.

thumbnail Fig. 1.

Stellar mass and SFR diagram of our galaxies estimated using the SED fitting code LePhare (Ilbert et al. 2006; Arnouts et al. 1999). Top and right panels: detection rates (fd) of continuum emission as functions of SFR and stellar mass. Filled and open points in the middle panel represent the continuum detected (> 3.5σ) and the non-detected galaxies (blue: z ∼ 4.5, red: z ∼ 5.5 galaxies), respectively. Solid and dashed lines show two different estimates of the z ∼ 5 main-sequence of star forming galaxies (Lee et al. 2012; Schreiber et al. 2015). Typical uncertainty of SFR and stellar mass are shown in the bottom right inset. The continuum detected galaxies mostly show stellar mass above ∼1010M, and SFR above ∼30 M yr−1.

As part of the SED fitting procedure in Faisst et al. (2020), we calculated monochromatic UV luminosities, LUV, without dust attenuation corrections at rest frame 1600 Å, and fitted a power-law function (fλ ∝ λβ) to measure the UV continuum slopes β using the wavelength range 1300Å−2300Å, consistent with previous studies (e.g., Bouwens et al. 2016; McLure et al. 2018; Fudamoto et al. 2020). Uncertainties were estimated from Monte Carlo simulations, by perturbing the fluxes of each filter assuming a Gaussian error distribution in the photometric measurement. We used median values as our best fits, and 16th and 84th percentiles as lower and upper uncertainties of the LUV and β measurements.

Based on the LUV measurements, we calculated UV-based SFRs without dust attenuation corrections by employing the equation from Madau & Dickinson (2014) converted to a Chabrier IMF as follows

(1)

The above conversion provides consistent SFRs with other studies (e.g., Kennicutt 1998), within our LUV measurement errors.

2.2. ALMA observations and data reduction

The details of the ALMA observations and the data reduction are presented in Béthermin et al. (2020). In short, ALMA observed our sample between 07 May, 2018 (Cycle 5) and 10 January, 2019 (Cycle 6) using antenna configurations C43-1 and C43-2. The lower side bands were used for both continuum measurements and the [CII] 158 μm emission line expected from the redshifts, and the upper side bands covered the dust continuum only. The integration times ranged from 15 to 45 min, with an average of 22 min.

After a basic calibration using the pipeline of the Common Astronomy Software Applications package (CASA; McMullin et al. 2007), and additional flagging for bad antennae, continuum maps were produced using the line free channels. In particular, we excluded channels within the ±3σ width of the detected [CII] emission lines (i.e., excluding above ∼1% of the maximum amplitude assuming a Gaussian profile). When the [CII] emissions have complex morphology (e.g., mergers or potential outflow; Jones et al. 2020; Ginolfi et al. 2020), the ±3σ widths may not exclude all components of the [CII] lines. In these cases, we further extended the [CII] line masks by ∼0.1 − 0.2 GHz to prevent the [CII] line from contaminating our continuum maps. Using these [CII] masks, we imaged continuum maps using the CASA task cmttTCLEAN with the natural weighting scheme to maximise sensitivity.

The resulting median (minimum−maximum) point-source sensitivities and resolutions are 41 μJy beam−1 (16.8−72.1 μJy beam−1) and 1.1″ (0.9 − 1.6″), respectively.

2.3. FIR continuum measurements and detection fractions

As described in Béthermin et al. (2020), the detection threshold used for the FIR continuum is 3.5σ within 2″ diameters of the UV counterpart position. Signal-to-noise ratios were determined using peak pixels and background rms. With this conservative approach, the fidelity of our detections is > 95% (Béthermin et al. 2020). In total, 23 of our galaxies have continuum detections at this level. We measured flux densities using 2D Gaussian fitting with our customised routine. Flux measurement uncertainties were estimated using a method provided by Condon (1997), which properly accounts for correlated noise as present in interferometric data.

For nondetections, we used upper limits for our analyses. We determined 3σ upper limits by searching the maximum flux within 2″ diameter of the UV counterpart position, and by adding three times the background rms to that maximum value. These upper limits are more conservative than in most previous works. In particular, these values account for potentially weak continuum signals that are just below the detection threshold (see Béthermin et al. 2020, for a discussion).

Figure 1 summarizes the continuum detections as a function of stellar mass and SFR. Clearly, the detection fractions are strongly mass and SFR dependent. Most of the FIR detected galaxies are massive (M ≳ 1010M), and highly star forming (SFRSED ≳ 30 M yr−1) galaxies. Nevertheless, some galaxies at the massive and the most star forming end did not show significant continuum detections, even though the sensitivity does not largely change. This suggests that the SFRSED and/or stellar mass alone are not ideal indicators of LIR.

3. Analysis

3.1. LIR estimation from ALMA observations

The total infrared luminosities LIR over the wavelength range of λrest = 8 − 1000 μm were estimated using the conversion factor from 158 μm to LIR presented in Béthermin et al. (2020). In particular, Béthermin et al. (2020) constructed a mean stacked continuum FIR SED for ALPINE galaxy analogs in terms of redshift, stellar mass and SFR, using the full multiwavelength dataset available in the COSMOS field (Davidzon et al. 2017). This includes deep Herschel (Lutz et al. 2011; Oliver et al. 2012), AzTEC/ASTE (Aretxaga et al. 2011), and SCUBA2 data (Casey et al. 2013). The mean stacked SED was then fit with a FIR template using the Béthermin et al. (2017) model, which uses an updated version of Magdis et al. (2012) templates based on the evolution of the average dust SEDs. The template has a UV interstellar radiation field ⟨U⟩ = 50, as parameterised in Béthermin et al. (2015). In particular, the template matches a modified black body with a relatively high peak dust temperature (Td ∼ 41 K) with an additional mid-infrared component that reproduces the Herschel fluxes.

This relatively high dust temperature is consistent with the extrapolation of lower redshift trends (e.g., Td >  40 K for z ≳ 3), as reported in several studies (Béthermin et al. 2015; Álvarez-Márquez et al. 2016; Ferrara et al. 2017; Faisst et al. 2017; Schreiber et al. 2018; Fudamoto et al. 2020). We note, however, that we do not have constraints on the FIR SEDs of individual galaxies. If extreme variations of dust temperatures are present within our sample, this will thus introduce a scatter on the derived ALMA continuum to LIR conversion factor. Since we are interested in the population average LIR values, however, we can ignore this scatter and simply adopt the best-fit template to the mean stacked FIR SED as derived in Béthermin et al. (2020).

To derive the LIR values in practice, our template was normalised to the continuum fluxes or 3σ upper limits measured at rest-frame 158 μm, and then integrated over the wavelength range 8 − 1000 μm. In this way, the measured monochromatic luminosity at λrest  =  158 μm was converted to LIR by ν158 μmLν158 μm/LIR  =  0.13 (Béthermin et al. 2020). The uncertainties on LIR were derived directly from the flux measurement uncertainties. Based on this LIR and the attenuation uncorrected LUV derived from the optical photometry, we computed the infrared excess as IRX = LIR/LUV, with the proper uncertainties.

Finally, we converted the LIR to obscured SFRs by employing the equation from Madau & Dickinson (2014) converted to the Chabrier IMF as follows

(2)

which is consistent with the values derived in Kennicutt (1998).

3.2. Stacking analysis

To obtain average properties of our galaxies, we performed a stacking analysis of all ALMA continuum images, including both individual detections and nondetections. We used image based stacking as described in Khusanova et al. (2020), which we briefly summarise here.

We performed stacking using the λrest = 158 μm continuum images centered on the UV counterpart positions. To create stacked images free from projected bright sources, we masked all serendipitous sources detected above 5σ significance based on the serendipitous detection catalog from Béthermin et al. (2020). We used weighted median stacking. That is, after aligning all images to the UV-based phase center, the stacked images were constructed by comparing the distributions of intensities for each pixel and taking a weighted median as follows:

(3)

where is the weighted median flux, fν158 μm, i is the rest-frame 158 μm continuum flux of the ith galaxy image, is the median UV luminosity of all galaxies in the stacking bin, and LUV, i is the UV luminosity of the ith galaxy. The inverse LUV weighting scheme was chosen to provide us with an accurate measure of the IRX (LIR/LUV). For testing purposes, we also computed unweighted median stacks, which did not differ significantly from our weighted medians.

If the UV and IR components within a given galaxy are significantly offset spatially, it could be possible that we lose part of the IR fluxes in our stacking procedure, as we use the UV-based phase center as the stacking reference. To test this, we performed stacks using the individual FIR continuum detections, and found that both the UV-centered stacks and the FIR-centered stacks resulted in identical stacked fluxes. This is consistent with the finding in Fujimoto et al. (2020), who concluded that the UV and FIR continuum positions of our sample are not significantly offset within the relatively large beam size (∼1″) of our observations.

Flux densities were measured from the stacks using aperture photometry (r = 1.5″), and the measurement uncertainties were estimated by boot strapping. As for our individual continuum detections (cf. Sect. 2.3), our detection thresholds for stacks are 3.5σ using the peak pixel flux density, where σ is the pixel background rms. When no continuum is detected in the stacks at this level, we determined conservative 3σ upper limits by searching the maximum pixel value within 2″ of the expected position, and by adding three times the background rms to the local maximum, as done for the individual sources in the ALPINE catalog (Béthermin et al. 2020).

We performed stacking in two different redshift bins of z <  5 and z >  5. In each redshift bin, we further split the sample in bins based both on UV slope and on stellar mass. For the β, LUV, and M values of each bin, we used the median values of the appropriate sample. The results of our stacking analysis are summarised in Table 1, while the stacked images are shown in Figs. 2 and 3.

thumbnail Fig. 2.

6″ × 6″ cutouts of β binned stacks of ALMA continuum images used to derive the stacked infrared luminosities. Upper and lower panels: stacks of galaxies at z ∼ 4.5 and z ∼ 5.5, respectively. Black solid contours show 2, 3, 4, 5σ and white dashed contours show −3, −2σ, if present. The ranges of β used in each stacks are shown in each cutouts. While in the z ∼ 4.5 bins all stacks have clear detections at > 3.5σ, in the z ∼ 5.5 bins the stack in the bluest bin (lower left panel) remains undetected and we only report a conservative 3σ upper limit (see text for details).

thumbnail Fig. 3.

Same as Fig. 2, but for M binned stacks. The ranges of M used in each stacks are shown in each cutouts (in log solar masses). While we obtained strong detections from most of the stacked images, in the lower mass bin of z ∼ 5.5 galaxies, only a ∼2σ peak is present at the source position. This is lower than our detection threshold for individual images, therefore we provide a conservative 3σ upper limit for this bin (see text).

Table 1.

Results of the stacking analysis.

For the IRX–β analysis the bins were chosen to split the sample roughly equally in the two different redshift bins. Specifically, at z <  5, we used β = [ − 2.65, −1.75], [ − 1.75, −1.25], [ − 1.25, −0.5], while for z >  5 galaxies we used β = [ − 2.65, −2.0], [ − 2.0, −1.7], [ − 1.7, −1.0] (Fig. 2).

From the stacks at 4 <  z <  5, we detected significant continuum emission from all β bins, while from the galaxies at 5 <  z <  6 we detected all but the bluest bin (i.e., β = −2.65 to −2.0). We note that the stellar mass range used in this β binning is comparable with the mass range used in several previous studies, such that our results can be directly compared (Álvarez-Márquez et al. 2019; Fudamoto et al. 2017, 2020).

Stacks based on M bins were used for the IRX–M analysis binned by log(M/M) = [9, 10] and log(M/M) = [10, 12] for both z <  5 and z >  5 galaxies. In total, 111 galaxies are used, while a few very low mass galaxies (7) are not included in these stacks.

From the stacks at 4 <  z <  5, we detected 158 μm continuum from all M bins. In the range 5 <  z <  6, the stack in the high mass bin has a strong detection, while the lower mass bin only shows a tentative signal at ∼3σ, which we report as a nondetection and provide our conservative 3σ upper limit (Fig. 3 and Table 1).

4. Results

4.1. IRX–β relation

We study the IRX–β relations by separating our sample in two different redshift intervals, at 4 <  z <  5 and 5 <  z <  6. We then compare our results with existing studies at z <  4 to investigate the evolution of dust attenuation properties over a wide redshift range between z ∼ 0 and z ∼ 6.

In order to connect the IRX–β diagram to the attenuation properties for an ensemble of galaxies, one has to make an assumption about their intrinsic, dust-free UV continuum slope β0, which depends on stellar population properties (e.g., metallicity, age). Once this is fixed, the slope of the attenuation curve directly determines the position of galaxies in the IRX–β space (see e.g., Salim & Narayanan 2020). In the following, we parameterize the slope of the attenuation curve as the change in FUV attenuation for a given change in UV continuum slope, dAFUV/dβ, following previous analyses (e.g., Meurer et al. 1999; Bouwens et al. 2016). The IRX–β relation can then be written as

(4)

where the prefactor 1.75 comes from the bolometric correction of the UV luminosity (BCUV). Several studies derived BCUV  ∼  1.7 (e.g., Meurer et al. 1999; Bouwens et al. 2016; Koprowski et al. 2018), and here we used BCUV  =  1.75 to make our assumption on the IRX–β relation consistent with the recent ALMA based study (Bouwens et al. 2016). We refer to “Meurer-like” and “SMC-like” attenuation based on the UV reddening slopes of dAFUV/dβ  =  1.99 (Meurer et al. 1999) and dAFUV/dβ  =  1.1, respectively (see also Reddy et al. 2018). Note that we are thus implicitly treating the extinction curve as measured toward stars in the Small Magellanic Cloud (SMC) as an attenuation curve.

4.1.1. The observed IRX–β relation

The measured IRX–β diagram of our z ∼ 4 − 6 galaxy sample is shown in Fig. 4. From the overall distribution of all the individual detections, 3σ upper limits, and stacks, we find that our galaxies generally do not follow the Meurer-like IRX–β relations from local starbursts, or from z ∼ 3 massive (M ≳ 1010M) star forming galaxies (e.g., Bouwens et al. 2019; Álvarez-Márquez et al. 2019; Fudamoto et al. 2017, 2020). The local IRX–β relation updated with photometry using larger apertures (i.e., Overzier et al. 2011; Takeuchi et al. 2012) are more consistent with our individual detections. However, these IRX–β relations require a relatively red or large β0 (e.g., β0 = −1.96, Overzier et al. 2011), which is not consistent with many of our β measurements including nondetections and the theoretically predicted evolution of β0 of metal poor high-redshift galaxies (e.g., Wilkins et al. 2011; Alavi et al. 2014; Reddy et al. 2018).

thumbnail Fig. 4.

Left: IRX–β diagram of ALPINE galaxies with previously determined relations at two different redshifts. Blue and red points show individual FIR continuum detections at 4 <  z <  5 and at 5 <  z <  6, respectively. Open downward triangles show 3σ upper limits of individual IR nondetections. Stacks are shown as blue (at z ∼ 4.5) and as red (at z ∼ 5.5) rectangles. The nondetection of stack is indicated by a downward arrow. In the bottom right of the left panel, a horizontal bar shows the median uncertainty of individual β measurements (Δβ = ±0.18). Solid and dotted lines show the IRX–β relation of local starbursts from Meurer et al. (1999) and that assuming SMC dust attenuation (e.g., Prevot et al. 1984), respectively. We also plot the updated local relation from Overzier et al. (2011) (dotted line). Assuming bluer intrinsic β, Reddy et al. (2018) proposed an IRX–β relation at z ∼ 2 which has a SMC dust extinction curve and bluer β (dot-dashed line). Right: fitting results to the stacks at each redshift assuming an intrinsic β0 of −2.62 are shown with red and blue lines with 1σ uncertainty (orange and blue bands). Small points and downward triangles show representative stacking detections and upper limits of z ∼ 3 galaxies, respectively (Bourne et al. 2017; Koprowski et al. 2018; Álvarez-Márquez et al. 2019; Fudamoto et al. 2020). Both at z ∼ 4.5 and z ∼ 5.5, our stacks show lower IRX than z ∼ 3 results at fixed β. The IRX–β relation from Meurer et al. (1999) is not consistent with our results. Within uncertainties, at z ∼ 4.5 and at z ∼ 5.5, our sample prefers the IRX–β relation from Reddy et al. (2018). Our stacks and fitting results indicate that galaxies follow an IRX–β relation most similar to an SMC type of attenuation.

The lower IRX values of our sample suggest an evolution of the average dust attenuation properties at z >  3, becoming more similar to an SMC-like dust attenuation. However, we note the very large dispersion present in IRX values at fixed β for individual detections, which reaches more than 1 dex relative to the stacked values in some cases. Such scatter to higher IRX has been seen before, and can be explained with geometric effects (e.g., Howell et al. 2010; Casey et al. 2014; Popping et al. 2017; Narayanan et al. 2018; Faisst et al. 2020). Despite this scatter within the population, it is clear that the average IRX values at fixed β, as measured through our stacks, are evolving compared to lower redshift derivations.

In detail, at z ∼ 4.5, more than half of the galaxies which are individually detected (10 out of 15) are located below the M99 relation. Furthermore, all the stacks with β >  −1.75 lie ≳0.5 dex below the M99 relation. Based on these results, we conclude that the IRX–β relation of our sample at z ∼ 4.5 does not follow a Meurer-like IRX–β relation, but requires a steeper attenuation curve.

At z ∼ 5.5, we detected only a small fraction of galaxies in the dust continuum (8 out of 51), but the majority of these individually detected sources lie somewhat below the M99 relation. Although the stack of the redder galaxies (β = [ − 1.75, −1.0]) is perfectly consistent with an SMC-like IRX–β relation, the stack of bluer galaxies (β = [ − 2.0, −1.75]) lies slightly above the SMC-like relation, while the bluest bin only resulted in an upper limit on the IRX. Nevertheless, we conclude that an SMC-like IRX–β relation is more consistent with our z ∼ 5.5 galaxy sample than the Meurer-like relation, as indicated also by other studies (Capak et al. 2015; Barisic et al. 2017). We quantify this in detail in the next section.

4.1.2. Fitting the IRX–β relation

Given the stacked IRX values as a function of β derived above, we can quantify the slope of the attenuation curve of our high-redshift galaxies by fitting the IRX–β relation. As shown in Eq. (4), the IRX–β relation depends both on the shape of the dust curve, through dAFUV/dβ, and on the intrinsic UV continuum slope β0. While, the classical value β0 = −2.23 derived in M99 for local galaxies is heavily used in the literature also for higher-redshift sources, several authors have pointed out that the physical conditions in early galaxies are likely very different, which affect the intrinsic slope, including younger stellar populations, lower metallicities, or potential changes to the initial mass function (e.g., Wilkins et al. 2011; Alavi et al. 2014; Castellano et al. 2014). In particular, Reddy et al. (2018) derived a more appropriate value of β0 = −2.62 for early galaxies, based on the binary population and spectral synthesis model (BPASS; Eldridge & Stanway 2012; Stanway et al. 2016) with a stellar metallicity of Z = 0.14 Z.

Indeed, this bluer β0 is clearly more consistent with the β distribution of our sample, which includes several galaxies that are significantly bluer than the nominal dust-free value of M99. In the following, we thus use β0 = −2.62 as our baseline value for our fits of the IRX–β relation, but we discuss how our results change by using the classical value of β0 = −2.23 derived in M99.

Using the IRX–β relation (Eq. (4)) with fixed β0 = −2.62, we fit the stacked IRX values as a function of β, which results in dAFUV/dβ = 0.71 ± 0.15 at z ∼ 4.5 and dAFUV/dβ = 0.48 ± 0.13 at z ∼ 5.5 (see Fig. 4 and Table 2). This is significantly steeper than a Meurer-like attenuation of dAFUV/dβ = 1.99, meaning that less UV attenuation is required to result in a given reddening of the UV slope. Therefore, in both redshift bins, our data reject a Meurer-like attenuation curve at > 5σ. The best-fit values even lie below the SMC extinction curve of dAFUV/dβ = 1.1, as can be appreciated from the comparison of our best-fit curves in Fig. 4 with the line derived in Reddy et al. (2018) using the same β0.

Table 2.

Fitting results of the IRX–β relation.

Our current data are unfortunately not good enough to constrain β0 of our sample, as such constraints require extremely deep IR observations of blue galaxies with accurate β measurements. However, the exact β0 value does not affect our conclusion of the SMC-like dust properties at z >  4. In particular, even if we change β0 to the classical value from M99 of β0 = −2.23, we still derive a best-fit dAFUV/dβ = 1.01 ± 0.20 for 4 <  z <  5 galaxies and dAFUV/dβ = 0.74 ± 0.24 values for 5 <  z <  6 galaxies (see Table 2). These values also exclude a Meurer-like attenuation at ∼4σ and at ∼5σ, respectively. Thus, our observations indicate that the attenuation curve of z >  4 UV-selected main-sequence galaxies are SMC-like, irrespective of the intrinsic UV slopes. In the appendix, we summarise the predicted IR luminosities based on our best fit IRX–β relations in Table A.1, together with the measured LIR for detections.

As shown in Fig. 4, previous studies found that the IRX–β relations of UV selected star forming galaxies at z ∼ 3 − 4 are consistent with the M99 relation using stacks of Herschel images (Heinis et al. 2014; Álvarez-Márquez et al. 2016, 2019; Koprowski et al. 2018). Also from individual detections and stacking analyses of ALMA observations, studies showed that massive UV selected, star forming galaxies at z ∼ 3 − 4 are consistent with the M99 relation, while lower mass galaxies show a relation close to the SMC-like IRX–β (Bouwens et al. 2016; Fudamoto et al. 2020). These previous z ∼ 3 − 4 observations are based on UV selected star forming galaxies (such as Lyman Break galaxies) similar to our analysis here, and they included both detections and nondetections in their stacks. Hence, they should be directly comparable to our work. Nevertheless, at the higher redshift probed by our sample here, our results indicate that SMC type relations are more applicable even for massive galaxies. Even within our sample itself, we find tentative evidence for a redshift evolution of the attenuation properties, given the different values derived for dAFUV/dβ at the 2σ level between z ∼ 4.5 and z ∼ 5.5. This is consistent with the previous analysis at z ∼ 5.5 (Capak et al. 2015; Barisic et al. 2017). Overall, these results indicate an evolution of the average dust properties (such as grain size distribution and/or composition) of UV-selected main-sequence galaxies at z ≳ 3.

4.1.3. Systematic trends with morpho-kinematics

In addition to the overall relation, we investigated if there is any correlations between a galaxy’s location in the IRX–β diagram and its morphology and/or kinematic conditions. Based on the detected [CII] emission lines kinematics and morphology, our sample was classified as rotators (9 galaxies), mergers (31 galaxies), extended dispersion dominated (15 galaxies), and compact dispersion dominated (8 galaxies) galaxies (Le Fèvre et al. 2020). We did not find any clear trends that systematically correlate with the locations of galaxies in the IRX–β diagram. Nevertheless, we noted that none of the compact, dispersion dominated galaxies (≲10% of our whole sample) are individually detected in the continuum. However, due to the small sample statistics, it is not yet possible to conclude, if this morpho-kinematic class of galaxies is systematically faint in the FIR.

4.2. IRX–M relation

In this section, we study the IRX–M relation of our galaxies by splitting the sample in the two redshift bins 4 <  z <  5 and 5 <  z <  6 again, and we compare our results with previously determined IRX–M relations at z ≲ 4.

The observed IRX–M diagram is shown in Fig. 5. As is evident, we find a strong redshift dependence of the relation at z >  4, despite the presence of a very large dispersion within the population. The individual detections and the stacked IRX values of our sample are generally lower (on average by ∼0.2 dex, but reaching up to ∼1 dex) than most of the previously determined IRX–M relations at lower redshifts z ∼ 1.5 − 4.0. In particular, in many cases, the 3σ upper limits of individual nondetections at the massive end of our sample are not consistent with previous relations. While the steep IRX–M relation presented in Fudamoto et al. (2020) using z ∼ 3 galaxies is still consistent with our z ∼ 4.5 sample, this is no longer the case at z ∼ 5.5 as the average IRX further decreases by ∼0.38 dex from redshift z ∼ 4.5 to z ∼ 5.5 (Fig. 5).

thumbnail Fig. 5.

IRX–M diagram of our galaxies compared to previously determined relations at different redshifts. Blue and red points show individual FIR continuum detections at 4 <  z <  5 and at 5 <  z <  6, respectively. Open downward triangles show 3σ upper limits of individual IR nondetections. Stacks are shown rectangles (blue at z ∼ 4.5, red for z ∼ 5.5). The nondetection of the stack is indicated by a downward arrow. The inset shows the median uncertainty of individual M estimations (log (ΔM/M) = ± 0.34). Lines show IRX–M relations derived in previous studies at z ∼ 2 − 4 (Álvarez-Márquez et al. 2016; Bouwens et al. 2016; Koprowski et al. 2018; Fudamoto et al. 2020). The gray band shows the 1σ uncertainty of the Fudamoto et al. (2020) relation. With the exception of the high mass bin at z ∼ 5.5, the population average stacked IRX values are consistent with this relation. However, the intrinsic dispersion from galaxy to galaxy at fixed stellar mass is clearly extremely large.

To compare our results in detail with the IRX–M relations previously determined at lower redshift, we use relations that are derived from UV selected star forming galaxies at 1.5 <  z <  4. The IR observations of these studies are based on Herschel (Heinis et al. 2014; Álvarez-Márquez et al. 2016, 2019; Koprowski et al. 2018), or ALMA (Fudamoto et al. 2020). Heinis et al. (2014) found no evolution of the IRX–M relation up to z ∼ 4. However, they do point out that the IRX values decrease at fixed stellar mass as a function of UV luminosity in a way that UV luminous galaxies show systematically lower IRX values than UV fainter ones. To present a fair comparison, we use the IRX–M relation from Heinis et al. (2014) obtained by stacks of the highest LUV bin (log LUV/L = 10.64 − 10.94), which is close to the typical LUV of our sample (median UV luminosities are log LUV/L = 10.91 for z ∼ 4.5 and log LUV/L = 10.86 at z ∼ 5.5). Thus, the sample selection bias of our rest-frame UV luminous galaxies has little, if any, effect on the comparison.

The previously derived IRX–M relations up to z ∼ 4 are mostly consistent with each other, and thus several studies agree on the nonevolution of the IRX–M relation over a wide redshift range. However, in a recent study, Fudamoto et al. (2020) found a significantly steeper IRX–M relation using an unbiased sample of star forming galaxies at z ∼ 3 by exploring the entire public ALMA archive data in COSMOS (A3COSMOS, Liu et al. 2019). The steeper slope in Fudamoto et al. (2020) could reflect a difference in the FIR SED used in previous studies and/or a potential measurement bias in previous stacking analyses, which relied on low-resolution data (e.g., Herschel), however the exact reason is not yet clear.

At z ∼ 4.5, in our ALPINE sample, we find that several individual detections are still consistent with previous relations within 3σ uncertainty (e.g., Heinis et al. 2014; Bouwens et al. 2016; Fudamoto et al. 2020). However, for individual nondetections at M >  1010M where the sensitivity of our observations provide strong constraints, all of our upper limits lie below the previous relations except for 2 galaxies whose upper limits are still consistent with the Fudamoto et al. (2020) relation. While the z ∼ 4.5 stacks show significant detections (upper panels of Fig. 3), and are still consistent with the steep IRX–M relation from Fudamoto et al. (2020), the z ∼ 4.5 stacks show much lower IRX than the previously determined relations at z ∼ 3 including the most UV luminous bin from Heinis et al. (2014). Comparing our stacks with the IRX–M relation of Bouwens et al. (2016), the IRX of z ∼ 4.5 galaxies, on average, are ∼0.63 dex lower at a fixed stellar mass.

Based on these considerations, we conclude the IRX–M relation of our sample at z ∼ 4.5 is consistent with that of Fudamoto et al. (2020), and deviate from the other previously found relations. This suggests that the IRX–M relation of main-sequence galaxies at 4 <  z <  5 either rapidly evolves and shows 0.6 dex lower IRX at a fixed stellar mass than the z ≲ 4 relations, or is consistent with the steep IRX–M relation of Fudamoto et al. (2020) found at z ∼ 3.

At z ∼ 5.5 this situation changes rapidly as almost all the individual detections lie below the IRX–M relations at z ≲ 4, and, at the massive end of our z ∼ 5.5 sample (i.e., M >  1010M), all the individual 3σ upper limits are below the z ≲ 4 IRX–M relations, except for one upper limit that is still consistent with Fudamoto et al. (2020). Stacking results emphasise the discrepancies between our z ∼ 5.5 sample and the z ≲ 4 IRX–M relations. While the stack of the lower mass bin (M = 109 − 1010M) does not show a significant detection (lower left panel of Fig. 3), its 3σ upper limit lies ∼0.4 dex below the Heinis et al. (2014) and the Bouwens et al. (2016) IRX–M relations. The stack of higher mass bin (M = 1010 − 1011M) shows a detection (lower right panel of Fig. 3), however its IRX is ≳1 dex below all the previously estimated IRX–M relations, suggesting that previously known IRX–M relations could over-predict the IR luminosities of z >  5 star forming galaxies by ∼1 dex.

The overall comparisons above demonstrate that the IRX–M relations from our observations start to become inconsistent with the previously determined IRX–M relations from z <  4 (except with steeper relation as in Fudamoto et al. 2020), and is no longer consistent with all the previously derived IRX–M relations by z ∼ 5.5. The rapid decrease of the IRX from z ≲ 4 to z >  4.5 in massive galaxies suggests a rapid evolution of dust attenuation properties of star formation in the high-redshift Universe, consistent with our conclusions from the IRX–β diagram.

4.3. The obscured fraction of star formation

Having estimated both the UV-based SFR and the infrared-based SFRs for our galaxies, we can estimate the obscured fraction of star formation occurring at z >  4. This obscured fraction, fobs, represents the fraction of star formation activity observable from IR continuum emission (i.e., dust obscured star formation activity) relative to its total amount, and is defined as fobs = SFRIR/SFRtot, where SFRIR is the SFR observed from IR (Eq. (2)) and SFRtot is the total star formation rate obtained by adding UV and IR based SFR estimates (i.e., Eqs. (1) and (2)).

Using Spitzer MIPS 24 μm observations of a mass complete sample at log M*/M ≳ 9, Whitaker et al. (2017) show that fobs is highly mass dependent, with more than 80% of star formation being obscured in massive galaxies with log M*/M ≳ 10. Remarkably, Whitaker et al. (2017) find that this fobs − M relation remains constant over the full redshift range z ∼ 2.5 to z ∼ 0. Since fobs is directly related to the IRX, the nonevolution of the fobs − M relation can be considered a product of the nonevolution of the IRX–M relation observed at z ∼ 1.5 − 3 (e.g., Heinis et al. 2014). In the same way, our finding of an evolving IRX–M relation at z >  4.5 (Sect. 4.2) thus implies an evolution of the obscured fraction of star formation at z >  4.

Figure 6 presents the fobs − M diagram of our galaxy sample. While individual detections and upper limits generally lie between fobs = 60 − 90%, given our ALMA sensitivity limits, the population average stacked values are significantly lower for 3 out of our 4 stacks. In particular, at M >  1010M, our stacks show lower values than the fobs − M relation of Whitaker et al. (2017). While z = 0 − 2.5 galaxies reach fobs ≳ 80% at these masses, we find at z ∼ 4.5 and only at z ∼ 5.5. Remarkably, in this high-mass bin at z ∼ 5.5, even all the individual detections and upper limits lie significantly below the lower redshift fobs − M relation.

thumbnail Fig. 6.

Obscured fraction of star formation as a function of stellar mass (fobs − M relation) of our UV selected sample. Blue and red points show individual FIR continuum detections at 4 <  z <  5 and at 5 <  z <  6, respectively. Triangles show 3σ upper limits for IR nondetections. Stacks are shown by blue (at z ∼ 4.5) and by red (at z ∼ 5.5) rectangles. The nondetection of the stack is indicated by a downward arrow. Gray points and lines show the observed relation at redshifts between z = 0 and z ∼ 2.5 (Whitaker et al. 2017). The solid line shows the constant fobs − M relation of z ∼ 2.5 to z ∼ 0 using a template from Dale & Helou (2002), and the dashed line shows the same using Béthermin et al. (2015) or Magdis et al. (2012) templates. At M <  1010M, our z ∼ 4.5 stacks are potentially consistent with the fobs at z ∼ 2.5 to z ∼ 0 using Béthermin et al. (2015) or Magdis et al. (2012) templates. However, at M >  1010M, our stacks show decreasing fobs from z ∼ 2.5 to z ∼ 5.5, suggesting a rapid evolution of dust obscured star formation activity in main-sequence galaxies at z >  4. However, we caution that UV-selected samples at z >  4 may be incomplete at very high masses (M >  1010.5M) given the absence of deep rest-frame optical imaging, as as been shown recently by the detection of a significant population of UV-undetected, massive galaxies (e.g., Wang et al. 2019; Alcalde Pampliega et al. 2019).

At M <  1010M, Whitaker et al. (2017) discussed that the estimated LIR systematically changes depending on the FIR SED templates used. Using their default template from Dale & Helou (2002), they find a significantly shallower trend with mass compared to a template set based on Béthermin et al. (2015) and Magdis et al. (2012), which incorporate an evolution of dust temperature as a function of redshift. Figure 6 shows both these trends. While the mass trend changes, Whitaker et al. (2017) show that the nonevolution of the fobs − M relation up to z ∼ 2.5 is conserved, independent of the template choice.

Our stacked data points at M <  1010M result in a mean of at z ∼ 4.5 and a 3σ upper limit of fobs <  0.43 at z ∼ 5.5. As seen in Fig. 6, these numbers are consistent with the lower redshift fobs − M relation of Whitaker et al. (2017) using the Béthermin et al. (2015) and Magdis et al. (2012) templates. However, the individual detection and the stacked data points at M >  1010M show rapid decrease of the obscured fractions at z >  4.5 for both template sets. This may thus indicate a different redshift evolution for low and high mass galaxies. It is also clear that the obscured fraction of star formation varies significantly from galaxy to galaxy, as is evident from the many individual continuum detections of low-mass galaxies at z ∼ 4.5, which imply fobs >  0.7, while the stacked median value is only .

One likely caveat of our analysis is that our sample is not perfectly mass complete. While ALPINE galaxies are selected to lie on the main-sequence, they were all required to have spectroscopic redshift measurements, which in most cases are based on rest-frame UV spectra. This could potentially bias our sample to somewhat more UV-luminous, less obscured systems. However, Faisst et al. (2020) show that the ALPINE sample only exhibits a weak bias toward bluer UV continuum slopes compared to a mass-selected parent sample with photometric redshifts zphot = 4 − 6. Additionally, the UV selection should affect the z ∼ 4.5 or z ∼ 5.5 galaxies in a similar way. Simple tests comparing the ALPINE sample with a mass-matched COSMOS parent sample indeed do not find systematic differences in IRX measurements as a function of UV slope.

Nevertheless, it is clear that extremely obscured, dusty sources, which still lie on the main sequence, would be missing from our sample. The existence of such a population of very massive, dusty galaxies at z >  3, which can remain undetected at rest-frame UV wavelengths, has recently been suggested in the literature based on ALMA or Spitzer detections (e.g., Franco et al. 2018; Wang et al. 2019; Williams et al. 2019; Casey et al. 2019; Gruppioni et al. 2020). The exact contribution to the total SFRD of such sources is still uncertain, given the currently small sample sizes and the difficulty of measuring their exact redshifts and stellar masses. However, it is clear that such UV-faint galaxies have the potential to dominate the SFRD at the massive end of the z >  3 galaxy population, at M ≳ 1010.5M. Solving this question will likely have to await the advent of the James Webb Space Telescope, which will provide much deeper rest-frame optical observations.

Keeping this potential sample bias at very high masses in mind, we can nevertheless conclude that our UV selected, main-sequence galaxies at M >  1010M show a much lower obscured fraction than galaxies at z ≲ 3. This implies a very rapid build-up of dust in such massive galaxies in the early universe, as the obscured fraction increases from ∼45% to ≳80% between z ∼ 5.5 to z ∼ 2.5. In contrast, at lower masses, our stacked measurements are consistent with the constant lower redshift fobs − M relation of Whitaker et al. (2017) with Béthermin et al. (2015) SED templates, which imply fobs ≲ 45% at M <  1010M.

5. Summary and conclusions

We have examined the IRX–β relation, the IRX–M relation, and the obscured fraction of star formation as a function of stellar mass (fobs − M relation) of UV-selected main-sequence star forming galaxies at z ∼ 4.4 − 5.8 using λrest = 158 μm continuum observations of 118 galaxies from the ALMA large program, ALPINE. The sample has secure spectroscopic redshifts, but is nevertheless representative of star forming galaxies, that is to say the distribution of SFRs and M are consistent with normal main-sequence galaxies in the observed redshift range: log (M/M)∼8.5 − 11.5 and log (SFR/M yr−1) ∼ 0.5 − 2.5 (see Faisst et al. 2020). From individual FIR measurements and stacks of both detections and nondetections we found:

  • (i)

    The IRX–β relation of our sample (Fig. 4) is generally located below the Meurer relation (Meurer et al. 1999) with average IRX values > 0.5 dex lower for galaxies redder than βUV >  −1.5, as hinted at by earlier studies (e.g., Capak et al. 2015; Barisic et al. 2017). Individual measurements of the UV spectral slope β suggest that the intrinsic (or dust free) UV spectral slope β0 is smaller than the locally estimated value of β0 = −2.23 (Meurer et al. 1999), and more consistent with a bluer β0, as expected for younger, more metal poor galaxies (e.g., β0 = −2.62, Reddy et al. 2018). We fit for the slope of the UV attenuation curve, finding it to be significantly steeper than that of local star forming galaxies, and similar to, or even steeper than that of the SMC (treating the SMC extinction curve as an attenuation curve; see Sect. 4.1).

  • (ii)

    Most of the IRX–M relations previously reported in the literature at z ≲ 3 are not consistent with the population average relation of our galaxies at z = 4.4 − 5.8 (Fig. 5). At z ∼ 4.5, our sample is still consistent with the steep IRX–M relation found by Fudamoto et al. (2020) for z ∼ 3 galaxies, even though individually detected galaxies show a very large scatter around the mean relation. However, at z ∼ 5.5, all the previously found IRX–M relations over-predict the IR luminosities by ∼1 dex, in particular for galaxies at the high-mass end of our sample, log (M/M) > 10.0.

  • (iii)

    The fraction of dust-obscured star formation among our UV-selected sample shows a decrease from fobs ∼ 65% at z ∼ 4.5 to only ∼45% at z ∼ 5.5 in massive (M ∼ (1 − 3)×1010M) galaxies with a potentially large scatter. This average obscured fraction is significantly lower than the ≳80% found for lower redshift galaxies, and provides direct evidence that the obscured fraction of star forming galaxies rapidly evolves in the early universe.

Taken together, our results indicate that the dust attenuation properties evolve rapidly between z ∼ 6, at the end of cosmic reionization, to z ∼ 2 − 3, around the peak of cosmic star formation. It is possible that we are seeing a transition to supernovae (SNe) driven dust production at z ≳ 3, which is predicted theoretically given the limited time available for dust production in lower mass stars (e.g., Todini & Ferrara 2001; Nozawa et al. 2003; Schneider et al. 2004). Under certain conditions, such SNe dust might also be consistent with the steeper dust curve we find for our z ∼ 4 − 6 galaxy sample compared to the M99-like attenuation inferred at z <  3 (e.g., Maiolino et al. 2004; Hirashita et al. 2005; Stratta et al. 2007; Gallerani et al. 2010). We refer to a future paper to analyze this possibility in detail.

Importantly, our results indicate that the previous dust attenuation corrections calibrated at z <  4 will over-predict the LIR of higher-redshift UV selected, star forming galaxies by up to ∼1 dex, especially for red and relatively massive galaxies. Future multiwavelength observations of UV-red galaxies (i.e., with β >  −1.5) will be crucial to further constrain the dust attenuation properties of galaxies at 5 <  z <  6 and to obtain a complete census of the total SFR density in the early universe. This is one of the key questions that can be addressed in the near future by exploiting the synergy between the FIR observations from ALMA and the rest-frame optical data coming from the James Webb Space Telescope.

Acknowledgments

This paper is dedicated to the memory of Olivier Le Fèvre, PI of the ALPINE survey. We thank the anonymous referee for insightful comments that improved the quality of this manuscript. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.00428.L. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan) and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. Based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under ESO programme ID 179.A-2005 and on data products produced by TERAPIX and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium. This work was supported by the Swiss National Science Foundation through the SNSF Professorship grant 157567 “Galaxy Build-up at Cosmic Dawn”. G.C.J. acknowledges ERC Advanced Grant 695671 “QUENCH” and support by the Science and Technology Facilities Council (STFC). D.R. acknowledges support from the National Science Foundation under grant number AST-1614213 and from the Alexander von Humboldt Foundation through a Humboldt Research Fellowship for Experienced Researchers. A.C., F.P. and M.T. acknowledge the grant MIUR PRIN2017. LV acknowledges funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Grant agreement No. 746119. S.F. is supported by and the Cosmic Dawn Center of Excellence funded by the Danish National Research Foundation under then grant No. 140. R.A. acknowledges support from FONDECYT Regular Grant 1202007.

References

  1. Alavi, A., Siana, B., Richard, J., et al. 2014, ApJ, 780, 143 [NASA ADS] [CrossRef] [Google Scholar]
  2. Alcalde Pampliega, B., Pérez-González, P. G., Barro, G., et al. 2019, ApJ, 876, 135 [CrossRef] [Google Scholar]
  3. Álvarez-Márquez, J., Burgarella, D., Heinis, S., et al. 2016, A&A, 587, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Álvarez-Márquez, J., Burgarella, D., Buat, V., Ilbert, O., & Pérez-González, P. G. 2019, A&A, 630, A153 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Aretxaga, I., Wilson, G. W., Aguilar, E., et al. 2011, MNRAS, 415, 3831 [Google Scholar]
  6. Arnouts, S., Cristiani, S., Moscardini, L., et al. 1999, MNRAS, 310, 540 [NASA ADS] [CrossRef] [Google Scholar]
  7. Barisic, I., Faisst, A. L., Capak, P. L., et al. 2017, ApJ, 845, 41 [NASA ADS] [CrossRef] [Google Scholar]
  8. Béthermin, M., Daddi, E., Magdis, G., et al. 2015, A&A, 573, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Béthermin, M., Wu, H.-Y., Lagache, G., et al. 2017, A&A, 607, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Béthermin, M., Fudamoto, Y., Ginolfi, M., et al. 2020, A&A, 643, A2 [CrossRef] [EDP Sciences] [Google Scholar]
  11. Bourne, N., Dunlop, J. S., Merlin, E., et al. 2017, MNRAS, 467, 1360 [NASA ADS] [Google Scholar]
  12. Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34 [Google Scholar]
  13. Bouwens, R. J., Aravena, M., Decarli, R., et al. 2016, ApJ, 833, 72 [NASA ADS] [CrossRef] [Google Scholar]
  14. Bouwens, R. J., Stefanon, M., Oesch, P. A., et al. 2019, ApJ, 880, 25 [CrossRef] [Google Scholar]
  15. Bowler, R. A. A., Bourne, N., Dunlop, J. S., McLure, R. J., & McLeod, D. J. 2018, MNRAS, 481, 1631 [NASA ADS] [CrossRef] [Google Scholar]
  16. Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682 [NASA ADS] [CrossRef] [Google Scholar]
  17. Capak, P. L., Carilli, C., Jones, G., et al. 2015, Nature, 522, 455 [NASA ADS] [CrossRef] [Google Scholar]
  18. Casey, C. M., Chen, C.-C., Cowie, L. L., et al. 2013, MNRAS, 436, 1919 [Google Scholar]
  19. Casey, C. M., Scoville, N. Z., Sanders, D. B., et al. 2014, ApJ, 796, 95 [Google Scholar]
  20. Casey, C. M., Zavala, J. A., Aravena, M., et al. 2019, ApJ, 887, 55 [Google Scholar]
  21. Castellano, M., Sommariva, V., Fontana, A., et al. 2014, A&A, 566, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
  23. Condon, J. J. 1997, PASP, 109, 166 [NASA ADS] [CrossRef] [Google Scholar]
  24. Dale, D. A., & Helou, G. 2002, ApJ, 576, 159 [NASA ADS] [CrossRef] [Google Scholar]
  25. Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Eldridge, J. J., & Stanway, E. R. 2012, MNRAS, 419, 479 [NASA ADS] [CrossRef] [Google Scholar]
  27. Faisst, A. L., Capak, P. L., Yan, L., et al. 2017, ApJ, 847, 21 [NASA ADS] [CrossRef] [Google Scholar]
  28. Faisst, A. L., Schaerer, D., Lemaux, B. C., et al. 2020, ApJS, 247, 61 [Google Scholar]
  29. Ferrara, A., Hirashita, H., Ouchi, M., & Fujimoto, S. 2017, MNRAS, 471, 5018 [NASA ADS] [CrossRef] [Google Scholar]
  30. Franco, M., Elbaz, D., Béthermin, M., et al. 2018, A&A, 620, A152 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Fudamoto, Y., Oesch, P. A., Schinnerer, E., et al. 2017, MNRAS, 472, 483 [NASA ADS] [CrossRef] [Google Scholar]
  32. Fudamoto, Y., Oesch, P. A., Magnelli, B., et al. 2020, MNRAS, 491, 4724 [CrossRef] [Google Scholar]
  33. Fujimoto, S., Silverman, J. D., Béthermin, M., et al. 2020, ApJ, 900, 1 [CrossRef] [Google Scholar]
  34. Gallerani, S., Maiolino, R., Juarez, Y., et al. 2010, A&A, 523, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Ginolfi, M., Jones, G. C., Béthermin, M., et al. 2020, A&A, 633, A90 [CrossRef] [EDP Sciences] [Google Scholar]
  36. Grazian, A., Giallongo, E., Fiore, F., et al. 2020, ApJ, 897, 94 [CrossRef] [Google Scholar]
  37. Graziani, L., Schneider, R., Ginolfi, M., et al. 2020, MNRAS, 494, 1071 [CrossRef] [Google Scholar]
  38. Gruppioni, C., Béthermin, M., Loiacono, F., et al. 2020, A&A, 643, A8 [CrossRef] [EDP Sciences] [Google Scholar]
  39. Hashimoto, T., Inoue, A. K., Mawatari, K., et al. 2019, PASJ, 71, 71 [NASA ADS] [CrossRef] [Google Scholar]
  40. Hasinger, G., Capak, P., Salvato, M., et al. 2018, ApJ, 858, 77 [NASA ADS] [CrossRef] [Google Scholar]
  41. Heinis, S., Buat, V., Béthermin, M., et al. 2014, MNRAS, 437, 1268 [NASA ADS] [CrossRef] [Google Scholar]
  42. Hirashita, H., Nozawa, T., Kozasa, T., Ishii, T. T., & Takeuchi, T. T. 2005, MNRAS, 357, 1077 [NASA ADS] [CrossRef] [Google Scholar]
  43. Howell, J. H., Armus, L., Mazzarella, J. M., et al. 2010, ApJ, 715, 572 [Google Scholar]
  44. Ilbert, O., Arnouts, S., McCracken, H. J., et al. 2006, A&A, 457, 841 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Jones, G. C., Béthermin, M., Fudamoto, Y., et al. 2020, MNRAS, 491, L18 [Google Scholar]
  46. Kennicutt, R. C., Jr. 1998, ARA&A, 36, 189 [Google Scholar]
  47. Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531 [NASA ADS] [CrossRef] [Google Scholar]
  48. Khusanova, Y., Bethermin, M., Le Fèvre, O., et al. 2020, A&A, submitted [Google Scholar]
  49. Koprowski, M. P., Coppin, K. E. K., Geach, J. E., et al. 2018, MNRAS, 479, 4355 [NASA ADS] [CrossRef] [Google Scholar]
  50. Koprowski, M. P., Coppin, K. E. K., Geach, J. E., et al. 2020, MNRAS, 492, 4927 [CrossRef] [Google Scholar]
  51. Laporte, N., Ellis, R. S., Boone, F., et al. 2017, ApJ, 837, L21 [NASA ADS] [CrossRef] [Google Scholar]
  52. Lee, K.-S., Ferguson, H. C., Wiklind, T., et al. 2012, ApJ, 752, 66 [Google Scholar]
  53. Le Fèvre, O., Tasca, L. A. M., Cassata, P., et al. 2015, A&A, 576, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Le Fèvre, O., Béthermin, M., Faisst, A., et al. 2020, A&A, 643, A1 [CrossRef] [EDP Sciences] [Google Scholar]
  55. Liu, D., Lang, P., Magnelli, B., et al. 2019, ApJS, 244, 40 [Google Scholar]
  56. Lutz, D., Poglitsch, A., Altieri, B., et al. 2011, A&A, 532, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [NASA ADS] [CrossRef] [Google Scholar]
  58. Magdis, G. E., Daddi, E., Béthermin, M., et al. 2012, ApJ, 760, 6 [NASA ADS] [CrossRef] [Google Scholar]
  59. Maiolino, R., Schneider, R., Oliva, E., et al. 2004, Nature, 431, 533 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  60. McLure, R. J., Dunlop, J. S., Cullen, F., et al. 2018, MNRAS, 476, 3991 [NASA ADS] [CrossRef] [Google Scholar]
  61. McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Data Analysis Software and Systems XVI, eds. R. A. Shaw, F. Hill, & D. J. Bell, ASP Conf. Ser., 376, 127 [Google Scholar]
  62. Meurer, G. R., Heckman, T. M., & Calzetti, D. 1999, ApJ, 521, 64 [NASA ADS] [CrossRef] [Google Scholar]
  63. Narayanan, D., Davé, R., Johnson, B. D., et al. 2018, MNRAS, 474, 1718 [NASA ADS] [CrossRef] [Google Scholar]
  64. Nozawa, T., Kozasa, T., Umeda, H., Maeda, K., & Nomoto, K. 2003, ApJ, 598, 785 [NASA ADS] [CrossRef] [Google Scholar]
  65. Oesch, P. A., Bouwens, R. J., Illingworth, G. D., Labbé, I., & Stefanon, M. 2018, ApJ, 855, 105 [NASA ADS] [CrossRef] [Google Scholar]
  66. Oliver, S. J., Bock, J., Altieri, B., et al. 2012, MNRAS, 424, 1614 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Overzier, R. A., Heckman, T. M., Wang, J., et al. 2011, ApJ, 726, L7 [NASA ADS] [CrossRef] [Google Scholar]
  68. Pannella, M., Elbaz, D., Daddi, E., et al. 2015, ApJ, 807, 141 [NASA ADS] [CrossRef] [Google Scholar]
  69. Popping, G., Puglisi, A., & Norman, C. A. 2017, MNRAS, 472, 2315 [NASA ADS] [CrossRef] [Google Scholar]
  70. Prevot, M. L., Lequeux, J., Maurice, E., Prevot, L., & Rocca-Volmerange, B. 1984, A&A, 132, 389 [Google Scholar]
  71. Reddy, N. A., Oesch, P. A., Bouwens, R. J., et al. 2018, ApJ, 853, 56 [NASA ADS] [CrossRef] [Google Scholar]
  72. Salim, S., & Narayanan, D. 2020, ARA&A, submitted [arXiv:2001.03181] [Google Scholar]
  73. Santini, P., Maiolino, R., Magnelli, B., et al. 2014, A&A, 562, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Schaerer, D., Ginolfi, M., Béthermin, M., et al. 2020, A&A, 643, A3 [CrossRef] [EDP Sciences] [Google Scholar]
  75. Schneider, R., Ferrara, A., & Salvaterra, R. 2004, MNRAS, 351, 1379 [NASA ADS] [CrossRef] [Google Scholar]
  76. Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Schreiber, C., Elbaz, D., Pannella, M., et al. 2018, A&A, 609, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Smit, R., Bouwens, R. J., Carniani, S., et al. 2018, Nature, 553, 178 [NASA ADS] [CrossRef] [Google Scholar]
  79. Stanway, E. R., Eldridge, J. J., & Becker, G. D. 2016, MNRAS, 456, 485 [NASA ADS] [CrossRef] [Google Scholar]
  80. Steinhardt, C. L., Speagle, J. S., Capak, P., et al. 2014, ApJ, 791, L25 [Google Scholar]
  81. Stratta, G., Maiolino, R., Fiore, F., & D’Elia, V. 2007, ApJ, 661, L9 [NASA ADS] [CrossRef] [Google Scholar]
  82. Takeuchi, T. T., Yuan, F.-T., Ikeyama, A., Murata, K. L., & Inoue, A. K. 2012, ApJ, 755, 144 [NASA ADS] [CrossRef] [Google Scholar]
  83. Tamura, Y., Mawatari, K., Hashimoto, T., et al. 2019, ApJ, 874, 27 [NASA ADS] [CrossRef] [Google Scholar]
  84. Todini, P., & Ferrara, A. 2001, MNRAS, 325, 726 [NASA ADS] [CrossRef] [Google Scholar]
  85. Wang, T., Schreiber, C., Elbaz, D., et al. 2019, Nature, 572, 211 [Google Scholar]
  86. Watson, D., Christensen, L., Knudsen, K. K., et al. 2015, Nature, 519, 327 [NASA ADS] [CrossRef] [Google Scholar]
  87. Whitaker, K. E., Pope, A., Cybulski, R., et al. 2017, ApJ, 850, 208 [NASA ADS] [CrossRef] [Google Scholar]
  88. Wilkins, S. M., Trentham, N., & Hopkins, A. M. 2008, MNRAS, 385, 687 [Google Scholar]
  89. Wilkins, S. M., Bunker, A. J., Stanway, E., Lorenzoni, S., & Caruana, J. 2011, MNRAS, 417, 717 [NASA ADS] [CrossRef] [Google Scholar]
  90. Wilkins, S. M., Feng, Y., Di Matteo, T., et al. 2018, MNRAS, 473, 5363 [CrossRef] [Google Scholar]
  91. Wilkins, S. M., Lovell, C. C., & Stanway, E. R. 2019, MNRAS, 490, 5359 [CrossRef] [Google Scholar]
  92. Williams, C. C., Labbe, I., Spilker, J., et al. 2019, ApJ, 884, 154 [CrossRef] [Google Scholar]

Appendix A: Measured and predicted infrared luminosities

Using the best fit IRX–β relations (Eq. (4) and Table 2), we estimated the IR luminosity of the galaxies observed in the ALPINE programme (Table A.1). While the IRX–β assuming β0 = −2.62 is our fiducial relation (Sect. 4.1), we list LIR estimated using the IRX–β relation with β0 = −2.32. When a β is bluer than β0 for each case, the galaxy is consistent to be dust-free, meaning that the estimated LIR is consistent with zero. The LIR estimated with our fiducial IRX–β relation are used to estimate the total SFRs in Schaerer et al. (2020).

Table A.1.

Summary of the total IR luminosity (λ = 8 − 1000 μm) estimations of our sample.

All Tables

Table 1.

Results of the stacking analysis.

Table 2.

Fitting results of the IRX–β relation.

Table A.1.

Summary of the total IR luminosity (λ = 8 − 1000 μm) estimations of our sample.

All Figures

thumbnail Fig. 1.

Stellar mass and SFR diagram of our galaxies estimated using the SED fitting code LePhare (Ilbert et al. 2006; Arnouts et al. 1999). Top and right panels: detection rates (fd) of continuum emission as functions of SFR and stellar mass. Filled and open points in the middle panel represent the continuum detected (> 3.5σ) and the non-detected galaxies (blue: z ∼ 4.5, red: z ∼ 5.5 galaxies), respectively. Solid and dashed lines show two different estimates of the z ∼ 5 main-sequence of star forming galaxies (Lee et al. 2012; Schreiber et al. 2015). Typical uncertainty of SFR and stellar mass are shown in the bottom right inset. The continuum detected galaxies mostly show stellar mass above ∼1010M, and SFR above ∼30 M yr−1.

In the text
thumbnail Fig. 2.

6″ × 6″ cutouts of β binned stacks of ALMA continuum images used to derive the stacked infrared luminosities. Upper and lower panels: stacks of galaxies at z ∼ 4.5 and z ∼ 5.5, respectively. Black solid contours show 2, 3, 4, 5σ and white dashed contours show −3, −2σ, if present. The ranges of β used in each stacks are shown in each cutouts. While in the z ∼ 4.5 bins all stacks have clear detections at > 3.5σ, in the z ∼ 5.5 bins the stack in the bluest bin (lower left panel) remains undetected and we only report a conservative 3σ upper limit (see text for details).

In the text
thumbnail Fig. 3.

Same as Fig. 2, but for M binned stacks. The ranges of M used in each stacks are shown in each cutouts (in log solar masses). While we obtained strong detections from most of the stacked images, in the lower mass bin of z ∼ 5.5 galaxies, only a ∼2σ peak is present at the source position. This is lower than our detection threshold for individual images, therefore we provide a conservative 3σ upper limit for this bin (see text).

In the text
thumbnail Fig. 4.

Left: IRX–β diagram of ALPINE galaxies with previously determined relations at two different redshifts. Blue and red points show individual FIR continuum detections at 4 <  z <  5 and at 5 <  z <  6, respectively. Open downward triangles show 3σ upper limits of individual IR nondetections. Stacks are shown as blue (at z ∼ 4.5) and as red (at z ∼ 5.5) rectangles. The nondetection of stack is indicated by a downward arrow. In the bottom right of the left panel, a horizontal bar shows the median uncertainty of individual β measurements (Δβ = ±0.18). Solid and dotted lines show the IRX–β relation of local starbursts from Meurer et al. (1999) and that assuming SMC dust attenuation (e.g., Prevot et al. 1984), respectively. We also plot the updated local relation from Overzier et al. (2011) (dotted line). Assuming bluer intrinsic β, Reddy et al. (2018) proposed an IRX–β relation at z ∼ 2 which has a SMC dust extinction curve and bluer β (dot-dashed line). Right: fitting results to the stacks at each redshift assuming an intrinsic β0 of −2.62 are shown with red and blue lines with 1σ uncertainty (orange and blue bands). Small points and downward triangles show representative stacking detections and upper limits of z ∼ 3 galaxies, respectively (Bourne et al. 2017; Koprowski et al. 2018; Álvarez-Márquez et al. 2019; Fudamoto et al. 2020). Both at z ∼ 4.5 and z ∼ 5.5, our stacks show lower IRX than z ∼ 3 results at fixed β. The IRX–β relation from Meurer et al. (1999) is not consistent with our results. Within uncertainties, at z ∼ 4.5 and at z ∼ 5.5, our sample prefers the IRX–β relation from Reddy et al. (2018). Our stacks and fitting results indicate that galaxies follow an IRX–β relation most similar to an SMC type of attenuation.

In the text
thumbnail Fig. 5.

IRX–M diagram of our galaxies compared to previously determined relations at different redshifts. Blue and red points show individual FIR continuum detections at 4 <  z <  5 and at 5 <  z <  6, respectively. Open downward triangles show 3σ upper limits of individual IR nondetections. Stacks are shown rectangles (blue at z ∼ 4.5, red for z ∼ 5.5). The nondetection of the stack is indicated by a downward arrow. The inset shows the median uncertainty of individual M estimations (log (ΔM/M) = ± 0.34). Lines show IRX–M relations derived in previous studies at z ∼ 2 − 4 (Álvarez-Márquez et al. 2016; Bouwens et al. 2016; Koprowski et al. 2018; Fudamoto et al. 2020). The gray band shows the 1σ uncertainty of the Fudamoto et al. (2020) relation. With the exception of the high mass bin at z ∼ 5.5, the population average stacked IRX values are consistent with this relation. However, the intrinsic dispersion from galaxy to galaxy at fixed stellar mass is clearly extremely large.

In the text
thumbnail Fig. 6.

Obscured fraction of star formation as a function of stellar mass (fobs − M relation) of our UV selected sample. Blue and red points show individual FIR continuum detections at 4 <  z <  5 and at 5 <  z <  6, respectively. Triangles show 3σ upper limits for IR nondetections. Stacks are shown by blue (at z ∼ 4.5) and by red (at z ∼ 5.5) rectangles. The nondetection of the stack is indicated by a downward arrow. Gray points and lines show the observed relation at redshifts between z = 0 and z ∼ 2.5 (Whitaker et al. 2017). The solid line shows the constant fobs − M relation of z ∼ 2.5 to z ∼ 0 using a template from Dale & Helou (2002), and the dashed line shows the same using Béthermin et al. (2015) or Magdis et al. (2012) templates. At M <  1010M, our z ∼ 4.5 stacks are potentially consistent with the fobs at z ∼ 2.5 to z ∼ 0 using Béthermin et al. (2015) or Magdis et al. (2012) templates. However, at M >  1010M, our stacks show decreasing fobs from z ∼ 2.5 to z ∼ 5.5, suggesting a rapid evolution of dust obscured star formation activity in main-sequence galaxies at z >  4. However, we caution that UV-selected samples at z >  4 may be incomplete at very high masses (M >  1010.5M) given the absence of deep rest-frame optical imaging, as as been shown recently by the detection of a significant population of UV-undetected, massive galaxies (e.g., Wang et al. 2019; Alcalde Pampliega et al. 2019).

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.