© ESO, 2016

1. Introduction

Understanding the formation mechanisms and evolution of galaxies with cosmic time has been, and is still, one of the major goals of both observational and theoretical astronomy. High redshift galaxies play a key role in developing our knowledge about the evolution of galaxies and intergalactic medium (IGM) in the distant universe. Since these galaxies are being studied within a few Gyr after the Big Bang, they provide a unique probe of the physics of one of the first generations of stars. One of the most well-known and efficient techniques for detecting galaxies at redshifts z> 2.5 is the Lyman-break technique, also known as the dropout technique (Steidel et al. 1996). This technique selects galaxies over a certain redshift interval by combining the photometry in three broad bands to detect the Lyman break, which is the drop in the ultraviolet (UV) flux observed in correspondence with the absorption of the intergalactic and interstellar medium at λ< 912Å (the Lyman limit, the wavelength below which the ground state of neutral hydrogen is ionized).

The Lyman-break galaxies (LBGs) represent the largest sample of star-forming galaxies known at high redshifts (z> 2.5) owing to the efficiency of their selection technique. They form a key population used to investigate the mass assembly of galaxies during the first Gyr of the Universe. Since the mid-1990s, their optica/near-infrared (NIR) rest-frame spectrum, combining both photometry and spectroscopy, have been studied by many authors and have provided good measurements of their stellar mass, 10⁹⁻¹⁰M_⊙, and star-formation rate, 10−100 M_⊙ yr^-1 (Shapley et al. 2001; Madau et al. 1996; Steidel et al. 1996; Giavalisco 2002; Blaizot et al. 2004; Shapley et al. 2005; Verma et al. 2007; Magdis et al. 2008; Stark et al. 2009; Chapman & Casey 2009; Lo Faro et al. 2009; Magdis et al. 2010c; Pentericci et al. 2010; Oteo et al. 2013b; Bian et al. 2013). Giving their masses and typical star formation rate (SFR) in the context of the ΛCDM model, the LBGs are believed to be the building blocks from which, by merging processes, a fraction of massive local galaxies (L>L^∗) have formed (Somerville et al. 2001; Baugh et al. 2005). The LBGs with the highest SFRs (>100 M_⊙ yr^-1) are thought to be progenitors of the present-day elliptical galaxies and the passive red galaxies at z ~ 2−3 (Verma et al. 2007; Stark et al. 2009).

The most commonly adopted LBG SFR tracer is from the rest-frame UV, optical, and NIR at high redshift where most of the energy is emitted by young stars (ages around 10 to 100 Myr). However, interstellar dust scatters and/or absorbs the light emitted by young stars, hence only a fraction of the energy output from star formation is observed in the UV and the rest is re-emitted over the full IR range 8−1000 μm. Burgarella et al. (2013) showed that, even at z = 3.6, about half of the star formation still resides in the far-infrared (far-IR), so it is necessary to combine both these tracers to determine the complete energy budget of star formation. However, statistical and representative information on the dust emission of these LBGs are still missing owing to their faintness in the far-IR and submillimeter (submm) wavelengths, which is very likely related to their low dust content (Giavalisco 2002).

Only a few LBGs selected using the dropout technique have been directly detected at z ~ 3 in the mid-infrared (mid-IR) (Magdis et al. 2010c), far-IR (Oteo et al. 2013a; Magdis et al. 2012; Casey et al. 2012), and submm (Capak et al. 2015; Chapman & Casey 2009; Chapman et al. 2000). This could be related to the fact that the most dust obscured objects may be located out of the color-color selection and also because of their intrinsic faintness in the IR. This subsample of LBGs that was detected in the far-IR and submm is extremely biased and not representative of the LBG population in terms of stellar mass, dust attenuation, and SFR. They are more likely the equivalent of submm bright galaxies, as suggested in Burgarella et al. (2011). Previous works have managed to detect them (Magdis et al. 2010b,a; Rigopoulou et al. 2010; Davies et al. 2013) using stacking analysis in the far-IR (using 922, 68, and 48 LBGs, respectively). The stacking analysis combines the signal of multiple sources that have been previously selected in other wavelength observations (Dole et al. 2006; Marsden et al. 2009; Béthermin et al. 2012; Heinis et al. 2013). However, these samples contain a small number of objects, which make it difficult to obtain a representative sampling of the LBGs characteristics. The sample used by Magdis et al. (2010b), Magdis et al. (2010a), and Rigopoulou et al. (2010) contains IRAC and MIPS detected LBGs. They are clearly biased towards the IR-bright, massive, and/or dusty LBGs. The most recent results from Coppin et al. (2015) presented the stacking of LBGs at z ~ 3, 4, and 5 in 850 μm images (with SCUBA-2), where they enlarged the sample of LBGs studied up to 4201 at z ~ 3.

The determination of the dust attenuation and star formation rate using the UV+IR emission is very challenging owing to the lack of deep IR data. Most of the studies of LBGs at high z use empirical recipes to correct the UV emission for dust attenuation. The most well-known is the relation between β_UV and the IR-to-UV luminosity ratio by Meurer et al. (M99; hereafter 1999), which is calibrated on local starburst galaxies. The β_UV and the L_FUV can be easily derived using the rest-frame UV colors when optical/NIR photometry data are available at high z and then the corrected UV luminosity can be computed using the M99 relation. This method has uncertainties, for example, the β_UV value is sensitive to the intrinsic UV spectral slope of galaxies – which depends on the metallicity, age of the stellar population, and star formation history – and the relation between the dust attenuation and β_UV depends on the dust properties and geometry (Calzetti 2001). The M99 relation is derived from local starburst galaxies and might not be valid for more “normal” star-forming galaxies (Buat et al. 2005). Indeed, various recent studies of local star-forming galaxies and high redshift UV-selected galaxies have found that the relation between dust attenuation and β_UV does not follow the M99 relation (Casey et al. 2014; Capak et al. 2015). We might wonder, however, if LBGs at low and high redshifts do.

It is also necessary to investigate the link between dust attenuation and other galaxy properties in addition to β_UV in order to be able to correct for dust attenuation in a statistical analysis of data collected in surveys. On the one hand, we know that the UV luminosity is not closely correlated with dust attenuation (Heinis et al. 2013, 2014). On the other hand, the stellar mass shows a good correlation with dust attenuation and does not appear to evolve with redshift (Heinis et al. 2014; Pannella et al. 2009; Ibar et al. 2013).

The present project makes use of the COSMOS field (Scoville et al. 2007), which covers 1.4 × 1.4 sq. deg. to select a large LBG sample (about 22 000 objects) and stack them in the far-IR and submm wavelength. The size of our sample allows us to divide them into several bins of L_FUV, β_UV, and M_∗ to present the observed statistics of the dust properties of LBGs at z ~ 3, as witnessed by the Herschel Space Observatory from Herschel-PACS+SPIRE data¹ and AzTEC from the Atacama Submillimeter Telescope Experiment (ASTE).

This paper is organized as follows. In Sect. 2, we present how we select the LBGs and gather the data to further analyze our sample. In Sect. 3, we detail the methodology used to stack the LBGs in the far-IR and submm, including incompleteness and clustering corrections applied to obtain valid results. In Sect. 4, we present the dust properties of our LBGs and discuss them in the context of their formation and evolution. Finally, we present our conclusions in Sect. 5. Throughout this paper we use a standard cosmology with Ω_m = 0.3, Ω_Λ = 0.7, H₀ = 70 Km s^-1, and the AB magnitude system. We employ the Chabrier (2003) initial mass function (IMF). When comparing our data to other studies, we assume no conversion is needed for SFR and stellar mass estimates between Kroupa (2001) and Chabrier (2003) IMFs. When converting from Salpeter (1955) IMF to Chabrier (2003) IMF, we divide M_{∗ Salpeter} by 1.74 (Ilbert et al. 2010), and SFR_Salpeter by 1.58 (Salim et al. 2007).

2. Data

2.1. UV/Optical/NIR data

We used optical imaging from the COSMOS field (Capak et al. 2007; Taniguchi et al. 2007), more specifically, the V_J and i⁺ bands. The images in the V_J and i⁺ have obtained from the SUBARU telescope using the Suprime-Cam instrument. They cover the entire COSMOS field reaching a 5σ depth of 26.5 and 26.1, respectively, for a 3′′ aperture.

In addition, we used the multi-color catalogue (Capak et al. 2007, version 2.0). The fluxes were measured in different bands from data taken at the Subaru (broad bands: B_J, V_J, g⁺, r⁺, i⁺, z⁺; intermediate bands: IA427, IA464, IA484, IA505, IA527, IA574, IA624, IA679, IA709, IA738, IA767, IA827; and narrow bands: NB711, NB816), CFHT (u^∗, i^∗, H and K_s bands), UKIRT (J-band), UltraVISTA (Y,J, H and K_s), Spitzer (3.6 - 8 μm), and GALEX (1500−3000 Å). The photometry was performed using SExtractor in dual-image mode where the source detection was run on the deepest image, i⁺. For the UV-NIR data, the point spread function (PSF) varies in the range between 0.5′′ and 1.5′′ from the u^∗ to the K images. In order to obtain accurate colors, all the images were degraded to the same PSF of 1.5′′ following the method described in Capak et al. (2007). The final photometric catalogue containes PSF-matched photometry for all the bands, measured over an aperture of 3′′ diameter at the position of the i⁺ band detection. It also provides the aperture correction to calculate the total flux for each object in the field.

2.2. Far-infrared data

We used observations of the COSMOS field from the ESA Herschel Space Observatory (Pilbratt et al. 2010). The Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) Evolutionary Probe survey (PEP; Lutz et al. 2011) mapped the COSMOS field at 100 and 160 μm with a point-source sensitivities of 1.5 mJy and 3.3 mJy and a PSF full width half maximum (FWHM) of 6.8′′ and 11′′. Also, the Spectral and Photometric Imaging Receiver (SPIRE, Griffin et al. 2010) as a part of the Herschel Multi-Tiered Extragalactic Survey (HerMES, Oliver et al. 2012) has observed the COSMOS field at 250, 350, and 500 μm. These maps have been downloaded from HeDaM². For the SPIRE maps the PSF FWHM is 18.2′′, 24.9′′, and 36.3′′; the 1σ instrumental noise is 1.6, 1.3, and 1.9 mJy beam^-1; and the 1σ confusion noise is 5.8, 6.3, and 6.8 mJy beam^-1 (Nguyen et al. 2010) at 250, 350, and 500 μm, respectively.

2.3. AzTEC data

The AzTEC observations in the COSMOS field have been described in Aretxaga et al. (2011). Data reduction has been performed using the standard AzTEC pipeline (see Scott et al. 2008; Downes et al. 2012). The observations have been taken using the Atacama Submillimeter Telescope Experiment (ASTE; Ezawa et al. 2004) to a depth ~1.3 mJy beam^-1. The FWHM of these observations was set to 33′′ by fitting the post-filtered PSFs of each set of observations with a Gaussian. The AzTEC observation covers only 0.72 deg² of the COSMOS field.

2.4. Photometric redshift and stellar masses

We used the photometric redshifts (photo-z) and stellar masses computed for the COSMOS field by Ilbert et al. (2009, version 2.0) for i-band detected sources. The photometric redshifts in the range 1.5 <z< 4 were tested against the zCOSMOS faint sample and faint DEIMOS spectra, showing that the accuracy of the photo-z is around 3% (Ilbert et al. 2009, version 2.0). The stellar masses were derived from spectral energy distribution (SED) fitting to the available optical and near-infrared photometry, assuming Bruzual & Charlot (2003) single stellar population templates, an exponentially declining star formation history, and the Chabrier (2003) IMF. The X-ray detected active galactic nuclei (AGN) were removed from the catalogue.

This catalogue provides us with three different calculations of the photo-z. In this paper, we make use of the photo-z derived from the median of the probability distribution function (PDF-z) or the one that minimizes the χ² (Chi-z); for more information see Ilbert et al. (2009, version 2.0) and see Sect. 3.1 here to follow the discussion of the photo-z selection for LBGs sample.

3. Sample

This section presents the characterization and selection of the LBG sample. First, we present the selection of the LBG sample in the color-color diagram (Sect. 3.1). Then, we explain how the calculation of the L_FUV and β_UV parameters are computed in the paper (Sect. 3.2). Next, we build a mock catalogue that will be used to define color-color selection criteria for our LBGs and to characterize the completeness as a function of different parameters (Sect. 3.3). Finally, we calculate the UV luminosity function (LF) from our LBGs (Sect.3.4).

3.1. LBG selection

Our LBGs at z ~ 3 are selected by means of the classical U-dropout technique (Steidel et al. 1996) using band filters u^∗, V_J, and i⁺. The LBGs at z ~ 3 must present lower fluxes in the bluest band (U-band), where the Lyman break is located since the selection requires detections in V_J and i⁺ bands.

We select objects brighter than the magnitudes V_J = 26.5 and i⁺ = 26.1 according to a 5σ depth in AB mag (Capak et al. 2007). We also require our objects to have photometric redshifts included in the catalogue of Ilbert et al. (2009, version 2.0), which means that the objects are out of the masked areas and are classified as galaxies. The flux in the u-band should be low because it is affected by the Lyman break; therefore, part of the LBGs are not detected in this band. So, we assign a magnitude equal to 1σ (28.7 mag) to these objects before applying the color selection. By using the mock catalogue (see Sect. 3.3) and photo-z, we derive the following color selection for LBGs at 2.5 <z< 3.5:

$\begin{array}{ccc}& & \mathit{u}\mathrm{-}{\mathit{V}}_{\mathrm{J}}\mathit{>}\mathrm{1}\\ & & {\mathit{V}}_{\mathrm{J}}\mathrm{-}{\mathit{i}}^{\mathrm{+}}\mathit{<}\mathrm{0.8}\\ & & \mathit{u}\mathrm{-}{\mathit{V}}_{\mathrm{J}}\mathit{>}\mathrm{3.2}\mathrm{(}{\mathit{V}}_{\mathrm{J}}\mathrm{-}{\mathit{i}}^{\mathrm{+}}\mathrm{)}\mathit{.}\end{array}$(1)Figure 1 shows the color-color diagram from the objects in the COSMOS field, with spectroscopic redshifts in the range 2.5 ≤ z ≤ 3.5 and the LBG selection defined. It also shows the redshift evolution of the expected colors (0.2 ≤ z ≤ 3.5) for three fiducial star-forming galaxies assuming E(B − V) = 0, 0.2, 0.4 from our mock catalogue.

Fig. 1 Selection of the LBGs in the color-color diagram. The small gray dots are the sample of objects detected in the COSMOS field, which have magnitudes (VJ< 26.5 and i+< 26.1), and photo-z calculated from Ilbert et al. (2009, version 2.0). The pink dots are our sample selected within the limits of Eq. (1) (black lines), with photometric redshift between 2.5 and 3.5 and log(LFUV[L⊙]) > 10.2. The yellow dots are the stars identified by Ilbert et al. (2009, version 2.0). The diamonds correspond to the objects with spectroscopic redshift identify in the range 2.5 <zspec< 3.5, color-coded red to blue as a function of their zspec. The black lines correspond to the expected colors in redshift evolution (small dots, 0.2 ≤ z ≤ 2.5; large dots, 2.5 ≤ z ≤ 3.5) for three fiducial star-forming galaxies (see Sect. 3.3.1) with E(B − V) = 0,0.2, and 0.4. The panel inserted in the top right corner shows the probability of finding an object in the redshift range 2.5 ≤ z ≤ 3.5 for our selection. The blue and red lines are the probability of finding an object with a good photo-z or spectrocopic redshift (available within the COSMOS collaboration), respectively. The back line is the probability of finding an object with a spectroscopic redshift after selecting the objects with photo-z in the range 2.5 ≤ z ≤ 3.5 as we did in our sample selection (which corresponds to our redshift probability as a function of redshift).

To reduce the effect of the incompleteness in the stacking process (see Sect. 4.1.1), we keep the objects with log  L_FUV[L_⊙]) ≥ 10.2. This luminosity corresponds to around 75% of the completeness for the objects detected in V_J and i⁺ bands (See Fig. 3.3.3).

In order to clean the sample of lower z interlopers, we select those galaxies whose PDF-z are within 2.5 ≤ z ≤ 3.5 (20, 819 galaxies). We also keep in our sample the objects for which the PDF-z are not in the adopted redshift range, but whose value of Chi-z and the error of the PDF-z agree with our redshift range (930 galaxies). We find that the fraction of sources that we recover using our color selection with respect to the photo-z is around 80 %. We also tested our sample against the COSMOS spectroscopic master catalogue (available within the COSMOS collaboration), which contain around 500 objects within 2.5 ≤ z_spec ≤ 3.5 and 97% of the sources present a magnitude lower than 25 in the i-band. We recovered 82% of the sources spectroscopically confirmed to be lying in the redshift range 2.5 ≤ z ≤ 3.5, after selecting objects with photo-z in the redshift range 2.5 ≤ z ≤ 3.5 (as we do for our LBG selection). The inset in the Fig. 1 shows the probability as a function of redshift for our selection criteria. This probability is calculated in relation to the total catalogue (number of objects inside our selection criterion for each redshift bin, divided by the number of objects in the total catalogue for the same redshift bin).

The final sample of LBGs contains ~22 000 LBGs with a z_mean = 3.02 ± 0.25. We note that the AzTEC observations cover 0.72 deg²; therefore, the sample used to stack in the AzTEC data is reduced to a subsample of approximately 7, 700 LBGs. Because the source selection criteria is exactly the same inside and outside the covered area, this should not introduce any bias.

The main goal of this work is to characterize the dust properties of LBGs at redshift z ~ 3 as a function of different physical parameters. The relatively large number of LBGs included in the sample allows us to study the LBGs in several bins as a function of their L_FUV, β_UV, and M_∗ to better investigate their variety. We split the sample as a function of L_FUV, β_UV, and M_∗, defining size bins of 0.3 dex, 0.4 dex, and 0.25 dex, respectively (see Table 1 for the number of the objects and interval values for each bin selected).

3.2. Far-UV luminosity and slope of the UV continuum

The values of β_UV and L_FUV throughout this paper are computed using the broad bands (B_J, V_J, g⁺, r⁺, i⁺), intermediate bands (IA464, IA484, IA505, IA527, IA574, IA624, IA679, IA709, IA738, IA767, IA827), and narrow bands (NB711, NB816) from the Capak et al. (2007, version 2.0) catalogue. We consider the UV rest-frame wavelength range, 1250 Å < λ < 2000 Å (Calzetti et al. 1994), to calculate β_UV and L_FUV. We exclude the range 2000 Å < λ < 2600 Å for two reasons, to omit the relevant dust feature at 2175 Å and to have a homogeneous rest-frame wavelength range independently of the redshift of the galaxy. We impose that the bands used in the analysis must be detected with a signal-to-noise (S/N) > 3σ. The values of β_UV and L_FUV are obtained by fitting the photometry to a simple power-law SED, f_λ ∝ λ^β_UV. The L_FUV is calculated at 1600 Å.

If we consider the rest-frame wavelength range used in this work, there are at least 11 bands available to calculate β_UV and L_FUV. The large number of bands reduces the error in the determination of β_UV and L_FUV. The uncertainty in the photometric redshift should influence our β_UV and L_FUV calculations. We consider that the UV spectrum follows a simple power-law. If the photo-z are perturbed according to their uncertainties, the slope of the linear fit has to remain the same. For our LBG sample (see Sect. 3.1) we find ⟨ σ_βUV ⟩ = 0.3 ± 0.1. On the other hand, the uncertainty in the photometric redshift produces an influence on the L_FUV value. However, this influence is not taken into account because it is smaller than the uncertainty in our L_FUV calculation (relative error equal to 20 ± 10%).

We would like to emphasize that our analysis is based on a statistical study of LBGs where we split the sample into different bins of β_UV, L_FUV, and M_∗. The errors showed above for the β_UV and L_FUV calculations are smaller than the bin sizes used in the binning of the sample. Therefore, they do not have a strong effect on the final results of our stacking analysis.

3.3. Mock catalogue

A mock catalogue is used to characterize the completeness of the LBG sample in different parts of this work. First, the photometric catalogue is cut at fixed L_FUV to reduce the incompleteness of the final LBGs sample (see Sect. 3.1). Then, we characterize the completeness of the LBG sample as a function of L_FUV and redshift to compute the UV LF (see Sect. 3.4). Finally, we estimate the incompleteness correction for each bin in the stacking as a function of L_FUV, β_UV, and M_∗ (see Sect. 4.1.1). Additionally, the mock catalogue is used to define the color selection criteria for our LBGs at z ~ 3 (see Sect. 3.1).

3.3.1. Model galaxies

We set up our fiducial galaxy SED model using CIGALE (Burgarella et al. 2005, in prep.; Noll et al. 2009; Ciesla et al. 2015; Boquien et al., in prep.)³. We simulate star-forming galaxies by considering a constant star formation history (SFH), an age of 100 Myr (as suggested by van der Burg et al. 2010), a subsolar metallicity (0.2 Z_⊙ as in Castellano et al. 2014), Bruzual & Charlot (2003) stellar population libraries, emission lines, and a Chabrier (2003) IMF. Computed template SEDs are then reddened by means of the standard Calzetti et al. (2000) attenuation law for starburst with a Gaussian distribution in E(B − V) following ⟨ E(B − V) ⟩ = 0.2 ± 0.1 and E(B − V) ≥ 0. The chosen E(B − V) distribution matches the β_UV distribution from our sample and the ones from our fiducial templates. We note that the assumed distribution for the E(B − V) is also in agreement with previous studies (e.g, Papovich et al. 2001; Reddy et al. 2008; Castellano et al. 2014), which suggests a mean between 0.1 and 0.2.

Spectral synthesis modeling and external multiwavelength information also indicate that the rest-frame UV wavelength range can be described for z ~ 2−3 galaxies by starbursts phases and constant star formation (e.g., Shapley et al. 2005). For instance, we are aware that some degeneracy exists, for example with respect to the age of the stellar populations and the amount of dust attenuation, that made no unique assumptions about modelling the colors of LBGs at z ~ 3. However, the chosen parameters provide representative ultraviolet colors of galaxies at z ~ 3, like LBGs (see van der Burg et al. 2010; Reddy et al. 2008).

3.3.2. Building the mock catalogue

We take advantage of the Monte Carlo approach to determine the transformation between the intrinsic properties of galaxies (UV luminosity, redshift, and reddening) and their observed rest-frame UV colors. We simulate our fiducial templates by assuming the same redshift distribution as that for the observed sample (i.e., 2.5 <z< 3.5). Then 200 000 fiducial galaxy SEDs are created using a flat distribution in redshift and the E(B − V) distribution aforementioned. CIGALE provides us with SEDs and intrinsic properties normalized to 1 M_⊙. To create the mock catalogue, we compute the L_FUV and the UV slope for the fiducial SEDs normalized to 1 M_⊙ using the same procedure exploited for our sample and we rescale them to a flat distribution in L_FUV (9 < log (L_FUV [ L_⊙ ] ) < 12). The Galactic extinction is added to the mock catalogue using the Cardelli et al. (1989) extinction curve. We add a stellar mass value for each object in the mock catalogue using the rest-frame UV light to mass ratio, keeping in mind the limitations (see next paragraph). The stellar mass is calculated by rescaling the value of the stellar mass, given by CIGALE for the models normalized to 1 M_⊙, by the same factor used to obtain the L_FUV. We also add a Gaussian scatter with a sigma equal to 0.5 dex to the stellar mass calculation.

The emission of our galaxies (modeled by assuming both a single constant SFH and one stellar population with a given distribution in E(B − V)) allows us to predict the UV rest-frame colors of z ~ 3 LBGs (e.g., Shapley et al. 2005), but not their stellar masses. Sawicki (2012) found for BX galaxies at z ~ 2.3 a stellar mass – L_FUV relation with a scatter of about 0.5 dex. Also, Hathi et al. (2013) found for LBGs at z ~ 1−5 a relation with a scatter of about 0.3 dex. For this, we chose to include a 0.5 dex dispersion in the calculation of the stellar mass to break the degeneracy in the method. We note that we only wish to obtain a realistic stellar mass scale to see the impact on the stacking analysis (see Sect. 4.1.1).

It is well known that the K-band rest-frame magnitude provides a good correlation with the stellar mass (Kauffmann & Charlot 1998). When we created fiducial galaxy SEDs with CIGALE, we also included the 8 μm IRAC band (K-band rest-frame at z ~ 3). We calculated the stellar mass using the mass-to-light ratio given by Magdis et al. (2010c), including the uncertainty in their relation. The stellar mass distributions of the mock catalogue estimated from the two methods are found to be consistent. This was done only to confirm that the first method gives a realistic stellar mass scale.

3.3.3. Completeness

To quantify the completeness as a function of different parameters (L_FUV, β_UV, M_∗, and redshift), which we will use in our sample selection, UV LF calculation, and stacking analysis, we insert the objects from the mock catalogue into i⁺ and V_J images. These two images are the ones where we perform the selection of our LBG sample (see Sect. 3.1). LBGs have typical half-light radii of r ~ 0.1′′−0.3′′ (Giavalisco 2002) and thus are unresolved by our observations. They can be treated as point sources. Futhermore, we assume a Gaussian profile in agreement with the PSF of the images to inject the objects. We inject 10 000 simulated objects, each time in random position on the image. We then attempt to recover these simulated galaxies using the same procedure as in Capak et al. (2007). We repeat it for the 200 000 objects of the mock catalogue. Figure 2 shows the completeness as a function of the L_FUV for objects detected in i⁺ and V_J images.

3.4. Luminosity function

We derive here the UV luminosity function (LF) to test whether our LBG sample is so far representative of the LBG population obtained in previous works. The UV LF is calculated by using the V_max method (Schmidt 1968) for our total LBG sample. The effective volumes (V_max) of our survey are given by $\begin{array}{ccc}& & {\mathit{\phi }}^{\mathit{i}}\mathrm{\left(}{\mathit{M}}_{\mathrm{UV}}\mathrm{\right)}\mathrm{d}\mathit{M}\mathrm{=}\sum _{\mathit{i}}^{{\mathit{N}}_{{\mathit{M}}_{\mathrm{UV}}^{\mathit{i}}\mathit{<}{\mathit{M}}_{\mathrm{UV}}\mathit{<}{\mathit{M}}_{\mathrm{UV}}^{\mathit{i}\mathrm{+}\mathrm{1}}}}\frac{\mathit{p}\mathrm{\left(}{\mathit{M}}_{\mathrm{UV}}^{\mathit{i}}\mathit{,z}\mathrm{\right)}}{{\mathit{V}}_{\mathrm{tot}}^{\mathit{i}}}\\ & & {\mathit{V}}_{\mathrm{tot}}\mathrm{=}\mathit{A}\mathrm{\int }\frac{\mathrm{d}{\mathit{V}}_{\mathrm{c}}}{\mathrm{d}\mathit{z}}\mathrm{d}\mathit{z,}\end{array}$(2)where ${\mathit{N}}_{{\mathit{M}}_{\mathrm{UV}}^{\mathit{i}}\mathit{<}{\mathit{M}}_{\mathrm{UV}}\mathit{<}{\mathit{M}}_{\mathrm{UV}}^{\mathit{i}\mathrm{+}\mathrm{1}}}$ is the number of objects in the bin of M_UV, A is the field area in deg², and $\frac{\mathrm{d}{\mathit{V}}_{\mathrm{c}}}{\mathrm{d}\mathit{z}}$ is the comoving volume per deg². The $\mathit{p}\mathrm{\left(}{\mathit{M}}_{\mathrm{UV}}^{\mathit{i}}\mathit{,z}\mathrm{\right)}$ is a function calculated using the mock catalogue and the completeness (see Sect. 3.3). This corresponds to the number of sources recovered with an observed magnitude in the interval [M_UV; M_UV+ΔM], which are selected as dropouts and divided by the number of injected sources with an intrinsic magnitude in the same interval [M_UV; M_UV+ΔM] and a redshift in the interval [z; z+Δz]. The magnitude is measured at rest-frame wavelength 1600 Å. The resulting luminosity function is binned to Δmag = 0.25.

The uncertainties in the luminosity function are derived combining the Poisson noise term, cosmic variance, and the determination of p(m,z). The cosmic variance in the COSMOS field using the mean mass of our sample and a redshift bin size, Δz = 1, amounts to 5% of LF (Moster et al. 2011). In the calculation of p(m,z) we estimate an uncertainty of around 5%, which comes from the Poisson noise in the determination of the completeness in M_FUV and z. We add the three uncertainties to compute the error on the LF.

Fig. 2 Completeness as a function of LFUV. We use the simulation from Sect. 3.3.1 to compute the completeness of our catalogue as a function of LFUV for objects detected in the VJ and i+ bands in the SUBARU observations. The error bars are calculated from the Poisson noise term.

We fit a Schechter function (Schechter 1976) to the binned data points, $\mathit{\phi }\mathrm{\left(}\mathit{M}\mathrm{\right)}\mathrm{d}\mathit{M}\mathrm{=}\mathrm{0.4}\ln \mathrm{\left(}\mathrm{10}\mathrm{\right)}{\mathit{\phi }}^{\mathrm{\ast }}{\mathrm{10}}^{\mathrm{0.4}\mathrm{(}\mathit{\alpha }\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}{\mathit{M}}^{\mathrm{\ast }}\mathrm{-}\mathit{M}\mathrm{)}}\exp \mathrm{(}\mathrm{-}{\mathrm{10}}^{\mathrm{0.4}\mathrm{(}{\mathit{M}}^{\mathrm{\ast }}\mathrm{-}\mathit{M}\mathrm{)}}\mathrm{)}$(3)with M^∗ being the characteristic magnitude, α the faint-end slope, and φ^∗ the overall normalization. We find φ^∗ = (2.65 ± 0.41) × 10^-3, M^∗ = −20.94 ± 0.08, and α = −1.84 ± 0.09. As can be seen in Fig. 3, our results agree within the 1σ region of the UV LF at 2.7 <z< 3.4 by Reddy & Steidel (2009) and van der Burg et al. (2010), but we find the UV LF systematically higher than these previous results, and that this excess increases toward the faint-end. In Sect. 3.1, we find that the fraction of sources recovered, combining both the color and the photo-z selections, with respect to the spectroscopic sample is around 82% (which means that we are including around 20% of sources that have a photo-z within the redshift range, but the spectroscopic redshifts are out of the target redshift range). The uncertainty in the photo-z calculation also increases for fainter objects. This is one of the possible explanations for the overdensity of objects observed in our LF. However, we should take into account that some differences can be also found due to differences in the selection of the sample and purity corrections. We do not use the same filters and color-selection to define our LBG sample as in the previous works (UBr bands), which can generate differences in the distribution of redshift between the different LBG samples. The previous LFs are computed in a shorter redshift range (2.7 <z< 3.4). We note that we clean the sample by using a photo-z selection (see Sect. 3.1). Nevertheless, the other mentioned LFs use simulations to obtain the fraction of stars and low redshift interlopers that contaminate the sample. The cosmic variance between the different fields can also produce differences in the LFs. Despite all of this, our M^∗ and α estimations of the LF are found to be within 1σ, and φ^∗ within 1.5σ with respect to Reddy & Steidel (2009) estimates. Therefore, we conclude that our LBG sample is roughly representative of the LBG population.

Fig. 3 UV LF for our LBGs sample in the redshift bin 2.5 <z< 3.5. Our data points and best-fit Schechter functions are shown in red. We compare them with the 1σ region for the UV LF of LBGs calculated by Reddy & Steidel (yellow, 2009) and van der Burg et al. (blue, 2010) in the redshift bin 2.7 <z< 3.4.

4. Stacking measurements

Stacking is a technique that combines the signal from multiple sources selected from observations at other wavelengths (e.g., Dole et al. 2006; Marsden et al. 2009; Béthermin et al. 2012; Heinis et al. 2013). With this stacking technique we are able to obtain a statistically significant measure of physical parameters in faint galaxies at high redshift, which are otherwise lost beneath the noise levels, at the expense of averaging over individual properties. We simultaneously stack 30 × 30 pixel cut-outs using the IAS library (Bavouzet 2008; and Béthermin et al. 2010)⁴ in PACS (100 and 160 μm images), SPIRE (250, 350, and 500 μm images), and AzTEC (1.1 mm image).

4.1. Stacking corrections

Lyman-break galaxies are clustered with other LBGs and other populations of star-forming galaxies at high redshift (Hickox et al. 2012), which causes a nonhomogeneous background in the stacked image (Bavouzet 2008; Béthermin et al. 2010; Heinis et al. 2013). To get valid and reliable results, we have to correct our measurements for two effects, incompleteness and clustering of the sample.

4.1.1. Correcting for incompleteness of the input catalogue in the dense regions

A bias is produced in the stacked image when the population of sources is not complete (Dole et al. 2006; Bavouzet 2008; Béthermin et al. 2010; Heinis et al. 2013; Viero et al. 2013). During the detection process in the optical images, we miss some of the faint objects located in the dense areas or close to bright objects. Therefore, we lose the contribution of the high background areas when we stack in the far-IR images, causing a negative flux contribution near to the stacked object in relation to the global background.

Fig. 4 Radial profile of the bias correction maps used for each bin of the stacking as a function of the LFUV (top) and M∗ (middle) at the SPIRE 250 μm. The bottom panel shows the bias correction used for each bin of the stacking as a function of the βUV because we did not find any difference in the bias correction as a function of the βUV.

We evaluate the contribution of this bias in the stacking using the method from Heinis et al. (2013). They showed that a bias effect of this nature increases inversely proportional to L_FUV. They selected the sample from the u-band, which corresponds to the UV (λ_UV = 1600 Å) in the rest-frame wavelength at z ~ 1.5. Our selection sample is different; the UV rest-frame wavelength is located in the r band at z ~ 3, and we select the objects in the V_J and i⁺ bands. Therefore, our stacking bias effects are related to the detection of the objects in the V_J and i⁺ images. To quantify the impact of this bias, we use the mock catalogue created to estimate the completeness (see Sect. 3.3.1). We split the recovered simulated sources in the V_J and i⁺ images from the mock catalogue in the same way as our LBGs sample. To build the bias correction map in each bin, in the far-IR and submm images we stack the random positions where we previously recovered the simulated sources in the V_J and i⁺ images from the mock catalogue. In order to correct our stacking measurement for this effect, we subtract the bias maps from the maps obtained by stacking at the position of the true LBGs.

Figure 4 shows the radial profile of the bias correction maps as a function of the L_FUV and M_∗. In agreement with Heinis et al. (2013), we find that the amplitude of this effect increases for objects with fainter L_FUV. This result is related to the well-known effect according to which the detection efficiency is lower in the denser of the UV images used for the color-color selection, and in particular for faint objects that are close to the brighter ones. This effect is also dependent on M_∗, but it presents a lower difference in the amplitude between the different bins of M_∗ than for the L_FUV. However, we did not detect any dependence in the UV-slope. The last two effects are related to the UV-luminosity dispersion in each bin when we stack as a function of M_∗ and β_UV. Considering that we did not find any evolution as a function of the UV slope, a bias correction made using the full recovered mock catalogue is used to correct each bin of the stacking as a function of the β_UV values (see bottom panel in Fig. 4).

4.1.2. Correcting for clustering of the input catalogue

The large PACS, SPIRE, and AzTEC beams can be contaminated by neighboring sources owing to the large PSF. This effect would lead to an overestimation of the mean flux of our sample due to the clustered nature of the sources. We use the formalism developed by Bavouzet (2008) and Béthermin et al. (2010) to address this. They assumed that the excess probability of finding another galaxy from a sample compared to a randomly distributed population is proportional to the angular auto-correlation function. The two-dimensional profile of the resulting stacking can then be written as $\mathit{I}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\mathrm{=}\mathit{\alpha }\mathrm{\times }\mathrm{PSF}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\mathrm{+}\mathit{\beta }\mathrm{\times }\mathit{w}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\mathrm{\ast }\mathrm{PSF}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\mathrm{+}\mathit{\gamma }\mathit{.}$(4)Here I(θ,φ) is the stacked map after bias correction, α is the average flux of the stacked population, PSF(θ,φ) is the point spread function (PSF) at the stacked wavelength, w(θ,φ) is the angular auto-correlation function of the input catalogue, β is a parameter related to the density of the input population, and γ is the constant related to the background in the stacked image.

We measure w(θ,φ) using the Landy & Szalay (1993) estimator in the position of our LBGs sample. We fit it with a power law w(θ,φ) ∝ θ^{− δ}, finding that the correlation function is well modeled with δ = 0.63. We check the slope of the auto-correlation function in each bin as a function of different parameters and we do not find any significant change (Heinis et al. 2013). Hence, we consider only the best fit to the auto-correlation function of the full sample. The auto-correlation function diverges when θ = 0, therefore, we assume that it takes effect out of 3 arcsec in relation to the center of the stacked image.

We solve the matrix system of Eq. (4) to obtain the flux density of our stacked objects. Bootstrap resampling is used to obtain the mean values and errors. We repeat the above procedure on 3000 random bootstrap samples, and the 1σ of the distribution of the derived fluxes is adopted as the uncertainty of our results. Figure 5 presents an example in profile for the solution of Eq. (4) in a specific L_FUV bin for the 250 μm image and its uncertainty. We also illustrate the contribution of the incompleteness and the clustering correction from the input catalogue.

Fig. 5 Example for the solution of Eq. (4) as a radial profile where we present the contribution of the corrections for incompleteness and clustering of the stacked galaxies. We show here the results for the third bin of the stacking as a function of LFUV in the 250 μm band (10.80 < log(LFUV[L⊙]) > 11.1, LBG-L3; see Table 1). We obtain the parameters of the fit: α = 0.99 ± 0.26, β = 0.46 ± 0.16, and γ = −0.12 ± 0.08. The black dashed line is the profile of the stacked LBGs population inside of the LFUV bin. The red line corresponds to the profile of the stacked LBGs population inside of the LFUV bin corrected for the incompleteness of the input catalogue, I (θ,φ). The green dashed line is the contribution for the clustering of the input catalogue, β × w(θ, φ) ∗ PSF(θ, φ). The blue line corresponds to the sum of the real emission of the object, α × PSF(θ , φ); the contribution for the clustering of the input catalogue, β × w(θ, φ) ∗ PSF(θ, φ); and the background, γ. The blue region shows the uncertainty of the results using bootstrap resampling. This plot is only a 1D illustration for the solution of Eq. (4), but we have solved it in 2D to obtain the stacked images for this work.

4.2. Stacking results

Figure 6 shows the stacked images for a specific mass bin in the PACS, SPIRE, and AzTEC bands; the stacked images are corrected for the bias maps. It indicates a significant detection of the emission for all the PACS, SPIRE, and AzTEC bands of the LBG population.

Table 1 shows the results of the stacking as a function of L_FUV, β_UV, and M_∗. We present the mean values of the parameters characterizing each population bin: L_FUV, β_UV, M_∗, photo-z, number of stacked LBGs, and their respective uncertainties (estimated by the standard deviation). The fluxes from the stacking measurement in each band and their errors are also listed. We also give the S/N calculated using the flux of our stacked LBGs and the σ⁵ of the stacked images, which are presented for each bin and band bracket.

Our selection was carried out as a function of L_FUV after applying the color selection criteria. For the stacking as a function of L_FUV, we find statistical detections for most of the stacked bands and bins by making use of the total sample. For the stacking as a function of β_UV and M_∗, we obtain no detection for the bins with β_UV< −1.9 and log(M_∗[L_⊙ ] ) < 9.75. We would need a larger number of stacked LBGs in each bin to reduce the contribution of the background and obtain a statistical detection.

Fig. 6 Example of the stacked images in the PACS, SPIRE, and AzTEC bands. Here we present the bin LBG-M3 (see Table 1) of the stacking as a function of stellar mass. The pixel sizes correspond to 2′′, 3′′, 6′′, 8.3′′, 12′′, and 3′′.

The AzTEC observations cover around three times less area than the entire COSMOS field. We stack only sources in the covered region to compute the AzTEC (1.1 mm) mean flux densities. We obtain a significant detection for the bins with, log(M_∗[L_⊙ ] ) > 10 in the stacking as a function of stellar mass. As we noted in Sect. 3.1, the source selection criteria are exactly the same inside and outside the covered area. We confirmed that this does not introduce any bias in stacking analysis by also stacking the AzTEC subsample in the PACS and SPIRE bands. The result of this stacking shows that the IR luminosity giving by the SED-fitting (Sect. 5.1) are within the uncertainties in both cases.

Fig. 7 Average IR luminosity obtained for our stacked LBGs as a function of the LFUV (left), βUV (middle), and M∗ (right).

5. Results and discussion

5.1. IR luminosity

We use the code CIGALE to fit the average stacked SED obtained for our LBGs to calculate the IR luminosity, which is estimated by integrating the best-fit Dale et al. (2014) templates over the range 8 <λ< 1000 μm. X-ray AGN are removed from our LBG sample since the mid-infrared emission of these objects could be strongly affected. X-ray obscured AGN might still be present in the sample. In this work, we do not have enough mid-IR data to constrain the AGN contribution for our stacked LBGs. However, Béthermin et al. (2015) have suggested that the possible presence of AGN emission appears to have a limited impact on the stacking analysis of star-forming galaxies. For this reason, we use the Dale et al. (2014) templates with an AGN contribution fixed to 0% and α_Dale as a free parameter to compute the IR luminosity. Where α is the exponent of the power law of the dust mass distribution with the radiation field intensity dU/ dM ∝ U^{− α}. Varying α provides a wide range of dust temperatures. The uncertainties are estimated by running Monte Carlo simulations where each stacked LBG is fit 3000 times. For each realization the flux in each band and the mean redshift of each stacked LBG are both randomly perturbed from the actual values by drawing from a Gaussian distribution. The width is given by the flux uncertainty in the specific band in the first case and by the standard deviation of the redshift distribution calculated using the photo-z for the objects in each bin in the second case. We take the standard deviation of the IR luminosities and α_Dale obtained in each iteration as the error. Fig. 8 shows IR SEDs and best models for our LBGs stacked as a function of their L_FUV, their β_UV, and their M_∗ values.

The average IR luminosities obtained for our stacked LBGs cover the range between 3 × 10¹⁰ to 3.23 × 10¹²L_⊙ (see Table 2) and thus have luminosities in the range of luminous infrared galaxies (LIRGs, 10¹¹<L_IR/L_⊙< 10¹²). However, there are some stacked LBGs with average IR emission in the ultra-luminous infrared galaxy (ULIRG, 10¹²<L_IR/L_⊙< 10¹³) range. Figure 7 shows the results of the average IR luminosities as a function of the L_FUV, β_UV, and M_∗ values. The three relations follow a power law where the IR luminosity increases as a function of the three parameters.

The shape of the average IR SEDs and the best-fit models (Fig. 8) suggest that there is an evolution of the dust temperature as a function of L_FUV, β_UV, and M_∗. Since α_Dale is inversely proportional to the dust temperature (Dale & Helou 2002; Chapman et al. 2003), the value of α_dale obtained reflect an evolution of the dust temperature (see Table 2). On the one hand, for the stacking as a function of L_FUV and M_∗, α_Dale increases; therefore, the average dust temperature for our LBGs decreases with these parameters. On the other hand, the average dust temperature shows an increasing trend with β_UV.

Fig. 8 For the first time, we are able to obtain well-sampled stacked SED of the LBGs for each bin in LFUV (top), M∗ (middle) and the βUV (bottom). Data points are the average flux densities calculated from the stacking procedure in each band. The 3σ upper limits (from bootstrap resampling) are shown by an arrow. We overplot the best-fit SED from the Dale et al. (2014) templates obtained using CIGALE.

5.2. Dust attenuation

It has been shown (e.g., Buat & Xu 1996; Gordon et al. 2000) that the IR-to-UV luminosity ratio (L_IR/L_FUV ≡ IRX) is a robust tracer of the dust attenuation in star-forming galaxies. We define the FUV dust attenuation, A_FUV, as a function of the IRX according to the prescription by Meurer et al. (1999), ${\mathit{A}}_{\mathrm{FUV}}\mathrm{=}\mathrm{2.5}\log \left(\frac{\mathrm{B}{\mathrm{C}}_{\mathrm{dust}}}{\mathrm{B}{\mathrm{C}}_{\mathrm{FUV}\mathit{,}\mathrm{\ast }}}\mathit{IRX}\mathrm{+}\mathrm{1}\right)\mathit{,}$(5)where A_FUV is the dust attenuation in the FUV, BC_{FUV, ∗} is the bolometric correction to the total light emitted by stars (BC_{FUV, ∗} = 1.68, M99), and BC_dust is the bolometric correction to the total light emitted by dust (BC_dust = 1). Despite the exact value of the bolometric correction for our specific sample of galaxies, we have chosen them to compare our best-fit A_FUV to other works that use similar notation.

5.2.1. IRX-β_UV relation

The β_UV slope has been found to be a good tracer of the UV dust attenuation (e.g., Meurer et al. 1999). It provides an estimate of the dust attenuation from the rest-frame UV without requiring far-IR data or spectral line diagnostics. Calibrations have been derived from spectro-photometric samples of starburst galaxies at low redshift (Meurer et al. 1999). Additional works (Overzier et al. 2011; Takeuchi et al. 2012) have provided new estimates of the original M99 relation using Galaxy Evolution Exploter (GALEX) data to compute the UV luminosity and solve the problem of the small aperture/field of view from the International Utraviolet Explorer (IUE). Takeuchi et al. (2012) also included a new L_IR estimation using AKARI data. For the same sample they found lower IRX and redder β_UV values than Meurer et al. (1999). A more recent work by Casey et al. (2014) proposed a new IRX-β_UV relation based on 1236 nearby galaxies observed by GALEX, for which they derive a much redder relation than the one presented by (Meurer et al. 1999) and roughly consistent with the aperture corrected M99 measured by Takeuchi et al. (2012). Our intention is to investigate whether the LBG population at z ~ 3 follows, on average, the same IRX-β_UV relation by stacking them at far-IR.

The M99 relation has been widely used to estimate dust-corrected SFRs at high redshift as SFR estimates based on IR, X-ray, or radio data are only available for the brightest objects (Bouwens et al. 2009). However, galaxies of different types are located differently in the IRX-β_UV plane. Young, metal-poor galaxies like the SMC and LMC are redder and less dusty and lie below the M99 relation. Instead, dusty star-forming galaxies (DSFG) related with the IR-bright galaxy population lie above the M99 relation (Casey et al. 2014; Oteo et al. 2013a). The dust geometry and star formation properties are thus seen to play a key role in the dispersion of the IRX-β_UV.

Figure 9 shows the IRX-β_UV relation for our stacked LBGs as a function of β_UV (x-axis) and dust attenuation (right y-axis). The calibration from local starburst galaxies (M99), the aperture photometry corrected the M99 by Overzier et al. (2011, O11; hereafter we use the IRX_M99,total derived in this paper), Takeuchi et al. (2012, T12), and the new local calibration by Casey et al. (2014) are also shown. In addition, we plot the LBGs detected in PACS at z ~ 3 (Oteo et al. 2013a) and the results of UV-selected galaxies at z ~ 1.5 by Heinis et al. 2013). Our data points are, within the errors, in excellent agreement with the correction of the M99 relation by (Takeuchi et al. 2012). The comparison of our stacked LBGs with the stacking results of a large sample (42, 184) of UV-selected galaxies at z ~ 1.5 by Heinis et al. (2013) confirms that we do not see a large evolution in the IRX-β_UV relation for average populations of galaxies selected from UV colors (LBGs and UV-selected), and from redshift 1.5 to 3. However, the IRX-β_UV relation derived from local galaxies by Casey et al. (2014) presents a different slope, but most of our data points show agreement within the uncertainty.

Fig. 9 IRX-βUV diagram. The right-hand axis shows the equivalent attenuation in the FUV band, in magnitudes, computed using Eq. (5). Our data points and the best fit to Eq. (5) are shown by red squares and the red solid line, respectively. The gray circles are the results of the stacking as a function of M∗, the size of the uncertainty in the x-axis is shown in the legend. Lines show various IRX-βUV relations: the local calibration of M99 (triple-dot-dashed line), the aperture correction of the M99 relation by T12 (dotted line), O11 (solid line), and the new local calibration by Casey et al. (2015; dashed line). We show measurements at z ~ 3 from direct LBGs detection in PACS by Oteo et al. (2013; blue open square), and the results of UV-selected galaxies at z~1.5 (Heinis et al. 2013). The results from this work are in good agreement with Takeuchi et al. (2012). Oteo’s detected LBGs seem to be extracted from biased IR-bright LBGs, not representative of the average population. It is interesting to note that AFUV continuously increases with increasing stellar masses.

The dust attenuation in Eq. (5) can be written as a function of β_UV: A_FUV = C₀ + C₁β_UV (Meurer et al. 1999; Overzier et al. 2011). We obtain C₀ = 3.15 ± 0.12 and C₁ = 1.47 ± 0.14 as best-fitting parameters for LBGs at z ~ 3. This equation implies that the UV slope of the dust-free objects is β_dust−free = −2.2 ± 0.3 and is in agreement, within the errors, with the value expected from stellar population models and a constant SFH mode (Leitherer & Heckman 1995). We note that the β_UV range study in this work does not allow the bluer part of the IRX-β_UV relation to be constrained; this is very important in order to determine the UV slope of a dust-free population.

Previous works have found the local starburst relation (M99) to hold for LBGs at various redshifts. Magdis et al. (2010a) find that the dust corrected UV-SFR derived from M99 relation presents a good match with the far-IR and radio SFR estimators by stacking an IRAC spectroscopically confirmed LBG sample at z ~ 3 in MIPS (24 μm), AzTEC (1.1 mm), and radio (1.4 GHz). This stacking analysis of spectroscopically confirmed z ~ 2 LBGs in far-IR (Reddy et al. 2012), z ~ 4 LBGs at radio continuum (1.4 GHz, To et al. 2014), and direct detection of LBGs in PACS at lower redshift (z ~ 1, Oteo et al. 2013b) also lie on the M99 relation. However, if we compare the IRX-β_UV relation obtained here with the original M99 relation, the previous works present bluer colors and/or higher dust attenuation than our mean LBG population. Buat et al. (2015) showed that the selection sample has an influence on the mean dust attenuation. They obtained differences of up to 2 mag between the UV-selected and IR-selected sample. In the particular case of Magdis et al. (2010a), their sample is IRAC-selected LBGs which is probably the origin of the higher dust attenuation. In Sect. 5.2.3, we will investigate in greater detail the effect on the dust attenuation due to the definition of our sample.

The LBGs detected in PACS at z ~ 3 (Oteo et al. 2013a) are found to be outliers of the IRX-β_UV relation. These galaxies are mainly ULIRG and hyper-luminous infrared galaxies (HLIRG, 10¹³ < L_IR/L_⊙ < 10¹⁴), similar to the DSFG and thus not representative of the average LBG population. Casey et al. (2014) showed that there is a deviation of this relation towards bluer colors for galaxies with L_IR> 10^11−11.5 due to the presence of recent and fast episodes of star formation that produce more prominent population of young O-B stars (contributing to the rest-frame far-UV emission) than galaxies of more modest SFRs. If we analyze the stacking results as a function of stellar mass in the IRX-β_UV plane, we find that they present bluer colors than our IRX-β_UV relation for the bins with log (M_∗ [ L_⊙ ] ) > 10.25 and L_IR> 7 × 10¹¹L_⊙, showing the same behavior as the DSFG by Casey et al. (2014). The stellar mass content in the galaxy could be the main driver of the dispersion of the IRX-β_UV plane to bluer colors or could simply be a consequence of the fact that more massive LBGs have associated larger infrared luminosity at high redshift.

5.2.2. IRX-M_∗ relation

The stellar mass has been shown to correlate with the dust attenuation in LBGs (Reddy et al. 2010), UV-selected galaxies (Buat et al. 2012; Heinis et al. 2014), and a mass-complete sample of galaxies (Pannella et al. 2015). The IRX-M_∗ relation presents no significant evolution with redshift, Pannella et al. (2015) found less than 0.3 mag of difference in the A_FUV from redshift ~0.7 to ~3.3, and Heinis et al. (2014) also showed consistent results from z ~ 1.5 to 4. In Sect. 5.2.1, we show that β_UV presents a correlation with the dust attenuation. The average β_UV for a population of galaxies and M_∗ clearly correlate at all redshifts (also for our LBG sample, see Table 1), but while the attenuation is fairly constant, or slightly increasing with redshift, β_UV becomes systematically bluer. The fact that high mass galaxies at high redshift have a similar dust attenuation but a bluer β_UV than similar mass galaxies at lower redshift has important implications for UV-derived SFRs in the high redshift Universe (Pannella et al. 2015). This leads to an inconsistency between dust attenuation measurements obtained using β_UV or M_∗. Figure 10 shows our measurements for the IRX as a function of M_∗. We also plot the relation from UV-selected galaxies by Heinis et al. (2014), the relation from a complete sample of star-forming galaxies by Pannella et al. (2015), and the LBGs detected in PACS at z ~ 3 (Oteo et al. 2013a).

Fig. 10 IRX versus M∗. The right-hand axis shows the equivalent attenuation in the FUV band in magnitudes, obtained using Eq. (5). Our data points and the best fit to Eq. (6) are shown as red squares and a red solid line, respectively. We show previous measurements at z ~ 3 from the stacking of UV selected galaxies by Heinis et al. (2014; dashed line), star-forming galaxies valid for redshift range 0.5−4 by Pannella et al. (2015; dotted line) and direct LBGs detections in PACS by Oteo et al. (2013; blue open square).

We assume here the following relation between the IRX and stellar mass, $\log \mathrm{\left(}\mathit{IRX}\mathrm{\right)}\mathrm{=}\mathit{\alpha }\log \left(\frac{{\mathit{M}}_{\mathrm{\ast }}}{{\mathrm{10}}^{\mathrm{10.35}}}\right)\mathrm{+}\mathit{IR}{\mathit{X}}_{\mathrm{0}}\mathit{,}$(6)and we obtain as best-fitting parameters, α = 0.84 ± 0.11 and IRX₀ = 1.17 ± 0.05.

We compare our results with previous published laws. We find that our IRX-M_∗ relation presents a steeper slope than those of Heinis et al. (2014) and Pannella et al. (2015). However, the Pannella et al. (2015) relation presents a good agreement with most of our data points in the range from 10 < log (M_∗ [ M_⊙ ] ) < 11.25. The difference could stem from the incompleteness in the low stellar mass bin that could present some variation in the dust attenuation and SFR with respect to a complete sample of galaxies.

In this case, the IRX-M_∗ relation tends to be aligned with the most massive (log (M_∗ [ M_⊙ ] ) > 11) LBGs detected in PACS by Oteo et al. (2013a). However, the lower stellar mass objects (log (M_∗ [ M_⊙ ] ) < 11) are located above the IRX-M_∗ relation showing high dust attenuation and/or the possible dispersion in the IRX-M_∗ plane owing to the different nature of the objects.

As already mentioned in Sect. 5.2.1 and in the first part of this section, the results obtained as a function of β_UV and M_∗ present different behaviors in the IRX-β_UV and IRX-M_∗ planes. The more massive LBGs in the stacking as a function of M_∗ suggest higher dust attenuation and bluer β_UV than the stacking as a function of β_UV.

5.2.3. Comparison between IRX-M_∗ and IRX-β_UV relations

Fig. 11 Results for stack 1 (green) and stack 2 (red) in the IRX-βUV and IRX-M∗ plane. Left panel: IRX-βUV diagram, for reference see Fig. 9. The stack 1 results are shown as diamonds, the tonalities of green represent the increase of the M∗ from 9.75 (light green) to 11.50 (dark green). The filled red circles are the results from stack 2, the tonalities of red represent the two different bins in M∗, 9.75 < log (M∗ [ M⊙ ] ) < 10.5 (light red) and 10.5 < log (M∗ [ M⊙ ] ) < 11.5 (dark red). Right panel: IRX-M∗ diagram, for reference see Fig. 10. The stack 1 results are shown as diamonds, the tonalities of green represent the two different bins in βUV, −1.7 <βUV< −0.5 (light green) and −0.5 <βUV< 0.5 (dark green). The filled red circles are the results from stack 2, the tonalities of red represent the increase of the βUV from −1.7 (light red) to 0.5 (dark red).

In Sect. 5.2.1 and 5.2.2, we present the mean IRX-M_∗ and IRX-β_UV relations for our LBG population. These previous stacking analyses do not provide any information to quantify the possible dispersion in the IRX relations. We also note that the two relations are not fully consistent; the more massive LBGs do not follow the IRX-β_UV relation presenting higher dust attenuation and bluer colors. We investigate here the dispersion and nature of the differences between the IRX-β_UV and IRX-M_∗ relations. We perform two new stacking analyses by building subsamples in the (M_∗, β_UV) plane in different ways. On the one hand, we split our sample into six bins of M_∗ and each of these bins is further separated into two bins of β_UV (stack 1)⁶. On the other hand, we split our sample into four bins of β_UV and each of these bins is further separated into two bins of M_∗ (stack 2)⁷. In these new stacking analyses, we limit the sample in the intervals as a function of M_∗ (9.75 < log (M_∗ [ M_⊙ ] ) < 11.5) and β_UV (−1.7 <β_UV< 0.5) obtaining around 9000 LBGs. We stack following the same procedure as in Sect. 4. We define this binning with the objective that all the stacked LBGs should be detected in SPIRE, because these bands are the main bands to compute the IR luminosity. We also stack in PACS and find a detection in part of them, but we did not use AzTEC in this analysis because we reach no detection. Figure 11 presents the results for both stacks in the IRX-β_UV and IRX-M_∗ planes.

In the IRX-β_UV plane, stack 1 presents a trend from the mean IRX-β_UV relation to higher IRX values and/or bluer colors. The lowest stellar mass bin has the same IRX as the mean IRX-β_UV relation, and the highest stellar mass bin tends to lie near the LBGs detected in PACS by Oteo et al. (2013a). This trend is well defined in the bluer bins, but the redder bins present a large dispersion owing to the low number of objects. Stack 2 shows the same behavior as stack 1: when the stellar mass increases, the IRX presents large values. These results suggest that the mean IRX-β_UV relation is not a well-defined relation valid for any object, but rather presents a large dispersion, which could be due to the differences in the stellar mass content in the galaxy for a given β_UV. Casey et al. (2014) proposed that the origin of the effect might be sought in the IR luminosity, but we show here that the M_∗ can also be the main driver. We should consider that the IR luminosity increases with M_∗; therefore, the IR luminosity and M_∗ correlate for our LBGs at z ~ 3 (see Fig. 7). So, both the IR luminosity and the stellar mass have an influence in the dispersion of the IRX-β_UV diagram.

The mean IRX-β_UV relation obtained in Sect. 5.2.1 is dominated for low M_∗ population (9.70 < log (M_∗ [ M_⊙ ] ) < 10.12, see Table 1). In the IRX-β_UV plane, the low stellar mass bin (log (M_∗ [ M_⊙ ] ~ 10.05) for stack 2 presents a good correlation with the mean IRX-β_UV relation. However, the high stellar mass bin (log (M_∗ [ M_⊙ ] ~ 10.75) shows a flat behavior with large IRX values. This means that the blue and high stellar mass LBGs lie above the mean IRX-β_UV relation and are located in the same area as the LIRGs and ULIRGs DSFG from Casey et al. (2014).

In the IRX-M_∗ plane, for stack 1, the bluer bins are in agreement with the previous results in the stacking as a function of M_∗. The redder bins are systematically above the mean IRX-M_∗ relation, the lower stellar mass bins show larger differences than the higher stellar mass bins, which are in agreement with the mean IRX-M_∗ relation. For stack 2, the low M_∗ bins present a dust attenuation departure from the IRX-M_∗ relation up to 1.3 mag higher. However, the high M_∗ bins are scattered around the mean IRX-M_∗ relation.

In stack 1, the dust attenuation shows a departure of up to 2.8 mag above the mean IRX-β_UV relation for the same β_UV bin when the log (M_∗ [ M_⊙ ]) increases from 9.75 to 11.5. In stack 2, the dust attenuation also shows a departure of up to 1.3 mag above the mean IRX-M_∗ relation for the same M_∗ bin when β_UV increases from −1.7 to 0.5. Previous works also showed that IR-selected galaxies have a mean dust attenuation that is ~2 magnitudes higher than the UV-selected galaxies at z ~ 0−2 (Buat et al. 2015). We suggest that the criterion in the selection of the sample has a strong impact on the mean dust attenuation of the population.

We would like to note that the IRX-M_∗ plane presents a lower dispersion than the IRX-β_UV for the M_∗ and β_UV intervals investigated here. However, we find that the objects with low stellar mass LBGs (log (M_∗ [ M_⊙ ] ) < 10.5) and red β_UV (β_UV> −0.7), i.e., ~15% of the total sample, present higher dust attenuation than the mean IRX-M_∗, but they are in agreement with the mean IRX-β_UV relation. We suggest that we have to combine both the IRX-β_UV and IRX-M_∗ relations to obtain the best estimation of the dust attenuation from the UV and NIR properties of the galaxies (L_FUV, β_UV, M_∗).

5.2.4. IRX-L_UV relation

We investigate how the stacking as a function of L_UV holds in the IRX-L_UV plane. Previous works show that the IRX remains more or less constant for different bins of L_UV in the average population of UV-selected galaxies. Xu et al. (2007) presented a sample of 600 UV-selected galaxies in GALEX with z ~ 0.6 and the stacking analysis in SWIRE; they obtained a constant IRX (log  (IRX) ~ 0.8) for the different bins in L_UV. More recent works (Heinis et al. 2013, 2014) have studied a large sample of UV-selected galaxies at z ~ 1.5, 3, and 4 in the COSMOS field by stacking them in SPIRE images. They find that the IRX remains more or less constant for the different bins in L_UV at z ~ 1.5, but they obtain a trend to higher dust attenuation for lower L_UV at high redshift. However, the IRX-L_UV plane presents a large dispersion for individual galaxies. Buat et al. (2015) has shown that IR selected have a dust attenuation ~2 mag higher than that the UV-selected galaxies at z ~ 0−2.

Figure 12 shows the IRX−L_FUV relation for our stacking sampling as a function of L_FUV. The IRX is found to be roughly constant over the range of L_FUV that we probe, with a mean of 7.9. This was already expected from Fig. 7, which shows a linear relation between L_FUV and L_IR. Our results are in agreement with the stacking analysis of a sample of UV-selected galaxies at z ~ 3 by Heinis et al. (2014). They suggest that the IRX-L_FUV relation presents a trend, but we cannot independently confirm it because the IRX uncertainty is in agreement with constant behavior and we are missing the faintest L_FUV bin. Coppin et al. (2015) present a stacking analysis of the whole sample (around 4200) of z ~ 3 LBGs in SCUBA2, SPIRE, 24 μm, and radio (1.4 GHz) finding an average IRX around 8. However, the spectroscopic sample of LBGs by Reddy et al. (2012) has lower IRX (7.1 ± 1.1), but consistent within 1σ error. We are also in agreement with the calculation of the dust attenuation of Burgarella et al. (2013), who found an IRX equal to 7.59 ± 6.99 and 5.54 ± 5.59 for galaxies at redshifts ~2.71 and 3.15, respectively.

Fig. 12 IRX versus LFUV.The right-hand axis shows the equivalent attenuation in the FUV band in magnitudes using Eq. (5). Our data points are shown as red squares. We also show previous measurements at z ~ 3 from stacking of UV-selected galaxies in the COSMOS field by Heinis et al. (2014; cyan open squares) and LBGs detected in PACS by Oteo et al. (2013; blue open square). The dashed lines represent the region where LIRG and ULIRG lie.

These results suggest that studies investigating the average population of LBGs can use a constant dust attenuation as a function of L_FUV to correct the observed L_FUV. For example, our mean dust attenuation could be used to derive the contribution of the LBGs at z ~ 3 to the star formation rate density by correcting the observed UV luminosity function (e.g., Madau & Dickinson 2014). However, there is a large dispersion on the IRX-L_FUV plane and, as mentioned above, the different sample selections present different dust attenuation (Buat et al. 2015). If we compare the stacking results as a function of the β_UV and M_∗ where we explore the IRX-L_FUV plane along the IRX axis, they present a variation on the dust attenuation from 1.4 to 3.8 and 0.5 to 2.9 mag, respectively. This gives a good view of dispersion of the IRX-L_FUV plane, where we use different selections from the same sample of galaxies; however, we lose the information of the L_FUV axis. This dispersion occurs because we can find the same L_FUV for a galaxy with large IR emission and dust attenuation (red galaxy) as for a galaxy that has low IR emission and dust attenuation (blue galaxy).

As we already discussed, our results suggest that the dust attenuation is roughly constant with L_FUV or equivalently with M_FUV. However, there is a debate on whether the β_UV evolves with the L_FUV at high redshift (Bouwens et al. 2014, 2012, 2009; Finkelstein et al. 2012). On the one hand, Bouwens et al. (2009, 2014) show a trend between the β_UV as a function of the M_FUV from redshift 2 to 8; they relate this to the metallicity and/or dust attenuation evolution. On the other hand, Finkelstein et al. (2012) found a flat distribution of the β_UV as a function of the M_FUV, but a trend with the stellar mass. For our sample of LBGs at z ~ 3, we show that the β_UV value is roughly constant for each bin in the stacking as a function of the L_FUV, consistent with the Finkelstein et al. (2012) results. Therefore, after seeing the differences between our results and those from Bouwens et al. (2009, 2014) in the β_UV − L_FUV relation, we can conclude that we cannot determine whether the trend in the β_UV − L_FUV relation comes from evolution in dust attenuation and/or a bias in the sample selection.

5.3. Star formation rate (SFR)

Since we have all the information to estimate reliable infrared luminosity for LBGs, we can combine the average estimates of L_IR and the observed uncorrected average UV luminosities to compute the total star formation rates (SFR = SFR_FUV + SFR_IR). We adopt the calibrations of Kennicutt (1998)

$\begin{array}{ccc}& & {\mathit{SFR}}_{\mathrm{IR}}\mathrm{\left(}{\mathit{M}}_{\mathrm{\odot }}\hspace{0.17em}{\mathrm{yr}}^{-1}\mathrm{\right)}\mathrm{=}\mathrm{1.10}\mathrm{\times }{\mathrm{10}}^{-10}{\mathit{L}}_{\mathrm{IR}}\mathrm{\left(}{\mathit{L}}_{\mathrm{\odot }}\mathrm{\right)}\\ & & {\mathit{SFR}}_{\mathrm{FUV}}\mathrm{\left(}{\mathit{M}}_{\mathrm{\odot }}\hspace{0.17em}{\mathrm{yr}}^{-1}\mathrm{\right)}\mathrm{=}\mathrm{1.82}\mathrm{\times }{\mathrm{10}}^{-10}{\mathit{L}}_{\mathrm{FUV}}\mathrm{\left(}{\mathit{L}}_{\mathrm{\odot }}\mathrm{\right)}\end{array}$rescaled from Salpeter (1955) to a Chabrier (2003) IMF.

5.3.1. SFR − stellar mass relation

Many recent studies have found evidence that the SFR in galaxies correlates with stellar mass along a main sequence (MS) relation, which evolves with redshift and represents a steady mode of SF. A linear relation seems to represent well the main sequence at intermediate masses; however, there are studies showing that the slope of this relation varies with the stellar mass or with different selections (e.g., Guzman et al. 1997; Brinchmann & Ellis 2000; Bauer et al. 2005; Bell et al. 2005; Papovich et al. 2006; Reddy et al. 2006; Daddi et al. 2007; Elbaz et al. 2007; Noeske et al. 2007; Pannella et al. 2009; Peng et al. 2010; Rodighiero et al. 2010, 2011; Karim et al. 2011; Whitaker et al. 2012, 2014; Rodighiero et al. 2014; Schreiber et al. 2015; Lee et al. 2015; Ilbert et al. 2015).

Figure 13 shows the average SFR − M_∗ relation computed for our stacked LBGs as a function of M_∗, along with previous stacking studies and observed relations (references in Fig. 13). Our results are in excellent agreement with the stacking at 1.4 GHz by Karim et al. (2011) and the stacking in far-IR by Schreiber et al. (2015) using a complete sample of star-forming galaxies. The far-IR stacking for UV selected galaxies by Heinis et al. (2014) presents higher SFRs than ours. It is interesting to note that we agree for the largest mass bin with the stacking at 850 μm of LBGs by Coppin et al. (2015); however, their lower mass bins show higher SFR. We propose that the differences between our results and theirs come from the absence of correction for clustering in their stacking analysis. This matches their overestimation of the fluxes of their stacked sources where the relative contribution of clustering becomes important in the faint sources. We would like to emphasize here that performing a careful stacking analysis is very important. The difference between the results from Coppin et al. (2015) and ours is an illustration that a statistically controlled stacking analysis is necessary to obtain reliable results.

The SFR-mass relation that we observe is described well by a power law with an average slope of 0.81 ± 0.09 in our mass range. We are in excellent agreement with the last results (0.8 ± 0.08) of Pannella et al. (2015) who use a complete sample of star-forming galaxies. We are also within the uncertainties with previous studies of individual LBGs detected in IRAC (Magdis et al. 2010c) at z ~ 3 and those lower redshift relations (Daddi et al. 2007; Elbaz et al. 2007; Peng et al. 2010), but using a higher normalization factor where they found a slope of ~0.9. We do not clearly detect any flattening of this relation at high M_∗ (Lee et al. 2015; Ilbert et al. 2015). We are also in good agreement with Schreiber et al. (2015) who showed this dependence for all redshifts in their analysis.

Fig. 13 Total SFR versus stellar mass for our stacking sample as a function of stellar mass. The red filled squares show our results. We compare them with previous stacking studies using different sample selections at redshift z ~ 3 for star-forming galaxies (Karim et al. 2011 and Schreiber et al. 2015), UV selected galaxies (Heinis et al. 2014), and LBGs (Coppin et al. 2015); they are all plotted with open symbols. The various lines show previous observed sequences at z = 0 (Peng et al. 2010), z = 1Elbaz et al. (2007), z = 2Daddi et al. (2007), and z = 3 (from Magdis et al. 2010c, based on IRAC-detected LBGs). When necessary, the SFRs and stellar mass values have been converted to a Chabrier (2003) IMF.

5.3.2. Difference between total SFR and corrected SFR by IRX-β_UV relation

We investigate here the differences between the total SFR obtained from UV+IR by stacking a population of galaxies and with SFR corrected by the mean IRX-β_UV relation. We calculate for each individual object the corrected SFR by SFR_corr = SFR_UV × 10^0.4A_FUV, where A_FUV = 3.15 + 1.47 β_UV. Then we split the sample into the same bins as the stacking as a function of M_∗, and we compute the average of the SFR_corr values. Figure 14 shows the comparison between the mean of the two SFR estimators for each bin in stellar mass and the SFR_corr calculate for each LBGs. We obtain that the SFR_corr presents lower SFR values for log (M_∗ [ M_⊙ ] ) > 10 than SFR_total. The differences increase for higher stellar mass. However, we find a good agreement for both SFR estimators in the lowest stellar mass bins. This confirms the conclusion in Sect. 5.2.3 that the SFR corrected from the IRX-β_UV relation underestimates the SFR of the population of galaxies with high stellar mass, as Pannella et al. (2015) found as well. But, it can be a good estimator of the average SFR for a populations of galaxies with low stellar mass.

The IRX relations derived in Sect. 5.2 are obtained from an average population of galaxies. As already shown, these relations present a dispersion that is due to the different nature of the galaxies. Therefore, their use could be useful but limited to works with an average population of galaxies, Applying them to individual galaxies can lead to strong under- or overestimation of the true dust attenuation and therefore SFR. We also show that the selection of the sample might have an impact on the resulting dust attenuation in the IRX planes, which means that depending on the sample selection (M_∗, β_UV, L_UV) and on the science objectives the most valid dust attenuation correction must be chosen.

Fig. 14 Comparison between SFRtotal and SFRcorr as a function of M∗. The red filled squares show the SFRtotal from the stacking results as a function of M∗. The gray point are the individual SFRcorr values for each LBG of our sample. The blue filled squares show the average of the SFRcorr; the error bars represent the error on the mean (σ/$\sqrt{{\mathit{N}}_{\mathrm{objects}}}$).

6. Summary and conclusions

We present in this paper a stacking analysis to study the rest-frame far-IR average properties of LBGs at z ~ 3. We combine the available COSMOS multiwavelength dataset to select a large sample of LBGs (around 22 000) using the dropout technique and photo-z. Thanks to the large number of objects included in our sample, we are able to split it into several bins allowing us to explore the evolution of the dust attenuation and star formation rate as a function of the galaxy properties (L_FUV, β_UV, and M_∗). We perform a stacking analysis in PACS (100 and 160 μm) images from the PACS Evolutionary Probe (PEP) team, SPIRE (250, 350, and 500 μm) images from the Herschel Multi-Tiered Extragalactic Survey (HerMES) programs and the AzTEC (1.1 mm) image from ASTE, obtaining an average flux density for the population in each bin and band. Our main results can be summarized as follows:

• 1.

  We compute the full infrared SEDs and we derive the average IRluminosity for our LBGs as a function of their L_FUV, β_UV, and M_∗. The obtainedIR luminosities cover the range 3 × 10¹⁰ to 3.23 × 10¹²L_⊙. We find that most ofthe stacked LBGs present an average L_IR similar to LIRGs, but themost massive and redder ones are ULIRGs. We find a power-lawcorrelation between the average L_IR and the parameters L_FUV, β_UV, and M_∗ that weuse to split the sample in each stacking analysis.

• 2.

  The average IRX (or dust attenuation) appears to be correlated with the β_UV slope. The relation follows the equation, A_FUV = (3.15 ± 0.12) + (1.47 ± 0.14) β_UV. It lies below the relation derived from local starburst galaxies by M99, but is in agreement with the correction to the M99 relation by T12, suggesting that compact blue starbursts chosen by M99 provide a good estimate of the dust attenuation for our average LBG population.

• 3.

  The average IRX (or dust attenuation) is correlated with the M_∗, following the relation: $\log {}^{\mathrm{(}}\mathit{IR}{\mathit{X}}^{\mathrm{)}}\mathrm{=}\mathrm{(}\mathrm{0.84}\mathrm{\pm }\mathrm{0.11}\mathrm{)}\mathrm{\times }\log {}^{\mathrm{(}}\mathit{M}_{\mathrm{\ast }}\mathit{/}{\mathrm{10}}^{\mathrm{10.35}}{}^{\mathrm{)}}\mathrm{+}\mathrm{1.17}\mathrm{\pm }\mathrm{0.05}$. The IRX-β_UV and IRX-M_∗ relations present an inconsistency; the more massive LBGs do not follow the IRX-β_UV relation and they show high dust attenuation and bluer color.

• 4.

  We perform a stacking analysis where we split our LBG sample in the (M_∗, β_UV) plane in two different ways. The objective of these new stackings is to study the dispersion in the IRX-β_UV and IRX-M_∗ planes and the inconsistency between the mean IRX-β_UV and IRX-M_∗ relations. The IRX-β_UV plane presents a large dispersion caused by differences in M_∗. The dust attenuation shows a departure of up to 2.8 mag above the mean IRX-β_UV relation for the same β_UV bin when the log (M_∗ [ M_⊙ ]) increases from 9.75 to 11.5. The IRX-M_∗ plane is less dispersed, and the dust attenuation also shows a departure of up to 1.3 mag above the mean IRX-M_∗ relation for the same M_∗ bin when β_UV increases from −1.7 to 0.5. However, the low stellar mass LBGs (log (M_∗ [ M_⊙ ] ) < 10.5) and red β_UV (β_UV> −0.7), which is 15% of the total sample, present higher dust attenuation than the mean IRX-M_∗, but they are in agreement with the mean IRX-β_UV relation.

• 5.

  We suggest that we have to combine both the IRX-β_UV and IRX-M_∗ relations to obtain the best estimation of the dust attenuation from the UV and NIR properties of the galaxies (L_FUV, β_UV, M_∗).

• 6.

  The IRX (and dust attenuation) is roughly constant over the L_FUV range for the average population of LBGs with a mean of 7.9 (1.8 mag). The IRX-L_FUV plane presents a large dispersion compared with the results in the stacking as a function of β_UV and M_∗ that can reach up to 2 mag.

• 7.

  The average SFR-M_∗ relation is approximated well by a power law with a slope of 0.81 ± 0.09 in our stellar mass range. We show that our LBG sample is consistent with the main sequence of star formation.

• 8.

  If we compare the total SFR (IR+UV) obtained using stacking analysis and the one calculated by correcting the L_FUV using the IRX-β_UV relation (SFR_corr), we demonstrate that the IRX-β_UV relation underestimates the SFR for high stellar mass LBGs, but it gives a good estimation for lower stellar mass LBGs.

• 9.

  The above relations provide phenomenological recipes to correct the observed L_FUV for dust attenuation, given a stellar mass or a β_UV. These recipes are useful in the absence of observed far-infrared data. However, even if they provide fairly good estimates of the amount of dust attenuation, we note that the results must be seen as statistical and valid for average populations of galaxies. The methods should not be applied to individual galaxies without assuming large uncertainties in the result.

Acknowledgments

PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/ OAA/OAP/OAT, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain). SPIRE has been developed by a consortium of institutes led by Cardiff Univ. (UK) and including: Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK); and Caltech, JPL, NHSC, Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA). The authors acknowledge financial contribution from the contracts PRIN-INAF 1.06.09.05 and ASI-INAF I00507/1 and I005110. The data presented in this paper will be released through the Herschel Database in Marseille (HeDaM; http://hedam.oamp.fr/HerMES). E. Ibar acknowledges funding from CONICYT/FONDECYT postdoctoral project N°:3130504. We gratefully acknowledge the contributions of the entire COSMOS collaboration consisting of more than 100 scientists. This work makes use of TOPCAT (http://www.star.bristol.ac.uk/?mbt/topcat/).

References

All Tables

All Figures

Fig. 1 Selection of the LBGs in the color-color diagram. The small gray dots are the sample of objects detected in the COSMOS field, which have magnitudes (VJ< 26.5 and i+< 26.1), and photo-z calculated from Ilbert et al. (2009, version 2.0). The pink dots are our sample selected within the limits of Eq. (1) (black lines), with photometric redshift between 2.5 and 3.5 and log(LFUV[L⊙]) > 10.2. The yellow dots are the stars identified by Ilbert et al. (2009, version 2.0). The diamonds correspond to the objects with spectroscopic redshift identify in the range 2.5 <zspec< 3.5, color-coded red to blue as a function of their zspec. The black lines correspond to the expected colors in redshift evolution (small dots, 0.2 ≤ z ≤ 2.5; large dots, 2.5 ≤ z ≤ 3.5) for three fiducial star-forming galaxies (see Sect. 3.3.1) with E(B − V) = 0,0.2, and 0.4. The panel inserted in the top right corner shows the probability of finding an object in the redshift range 2.5 ≤ z ≤ 3.5 for our selection. The blue and red lines are the probability of finding an object with a good photo-z or spectrocopic redshift (available within the COSMOS collaboration), respectively. The back line is the probability of finding an object with a spectroscopic redshift after selecting the objects with photo-z in the range 2.5 ≤ z ≤ 3.5 as we did in our sample selection (which corresponds to our redshift probability as a function of redshift).

In the text

	Fig. 2 Completeness as a function of LFUV. We use the simulation from Sect. 3.3.1 to compute the completeness of our catalogue as a function of LFUV for objects detected in the VJ and i+ bands in the SUBARU observations. The error bars are calculated from the Poisson noise term.
In the text

	Fig. 3 UV LF for our LBGs sample in the redshift bin 2.5 <z< 3.5. Our data points and best-fit Schechter functions are shown in red. We compare them with the 1σ region for the UV LF of LBGs calculated by Reddy & Steidel (yellow, 2009) and van der Burg et al. (blue, 2010) in the redshift bin 2.7 <z< 3.4.
In the text

	Fig. 4 Radial profile of the bias correction maps used for each bin of the stacking as a function of the LFUV (top) and M∗ (middle) at the SPIRE 250 μm. The bottom panel shows the bias correction used for each bin of the stacking as a function of the βUV because we did not find any difference in the bias correction as a function of the βUV.
In the text

Fig. 5 Example for the solution of Eq. (4) as a radial profile where we present the contribution of the corrections for incompleteness and clustering of the stacked galaxies. We show here the results for the third bin of the stacking as a function of LFUV in the 250 μm band (10.80 < log(LFUV[L⊙]) > 11.1, LBG-L3; see Table 1). We obtain the parameters of the fit: α = 0.99 ± 0.26, β = 0.46 ± 0.16, and γ = −0.12 ± 0.08. The black dashed line is the profile of the stacked LBGs population inside of the LFUV bin. The red line corresponds to the profile of the stacked LBGs population inside of the LFUV bin corrected for the incompleteness of the input catalogue, I (θ,φ). The green dashed line is the contribution for the clustering of the input catalogue, β × w(θ, φ) ∗ PSF(θ, φ). The blue line corresponds to the sum of the real emission of the object, α × PSF(θ , φ); the contribution for the clustering of the input catalogue, β × w(θ, φ) ∗ PSF(θ, φ); and the background, γ. The blue region shows the uncertainty of the results using bootstrap resampling. This plot is only a 1D illustration for the solution of Eq. (4), but we have solved it in 2D to obtain the stacked images for this work.

In the text

	Fig. 6 Example of the stacked images in the PACS, SPIRE, and AzTEC bands. Here we present the bin LBG-M3 (see Table 1) of the stacking as a function of stellar mass. The pixel sizes correspond to 2′′, 3′′, 6′′, 8.3′′, 12′′, and 3′′.
In the text

	Fig. 7 Average IR luminosity obtained for our stacked LBGs as a function of the LFUV (left), βUV (middle), and M∗ (right).
In the text

Fig. 8 For the first time, we are able to obtain well-sampled stacked SED of the LBGs for each bin in LFUV (top), M∗ (middle) and the βUV (bottom). Data points are the average flux densities calculated from the stacking procedure in each band. The 3σ upper limits (from bootstrap resampling) are shown by an arrow. We overplot the best-fit SED from the Dale et al. (2014) templates obtained using CIGALE.

In the text

Fig. 9 IRX-βUV diagram. The right-hand axis shows the equivalent attenuation in the FUV band, in magnitudes, computed using Eq. (5). Our data points and the best fit to Eq. (5) are shown by red squares and the red solid line, respectively. The gray circles are the results of the stacking as a function of M∗, the size of the uncertainty in the x-axis is shown in the legend. Lines show various IRX-βUV relations: the local calibration of M99 (triple-dot-dashed line), the aperture correction of the M99 relation by T12 (dotted line), O11 (solid line), and the new local calibration by Casey et al. (2015; dashed line). We show measurements at z ~ 3 from direct LBGs detection in PACS by Oteo et al. (2013; blue open square), and the results of UV-selected galaxies at z~1.5 (Heinis et al. 2013). The results from this work are in good agreement with Takeuchi et al. (2012). Oteo’s detected LBGs seem to be extracted from biased IR-bright LBGs, not representative of the average population. It is interesting to note that AFUV continuously increases with increasing stellar masses.

In the text

Fig. 10 IRX versus M∗. The right-hand axis shows the equivalent attenuation in the FUV band in magnitudes, obtained using Eq. (5). Our data points and the best fit to Eq. (6) are shown as red squares and a red solid line, respectively. We show previous measurements at z ~ 3 from the stacking of UV selected galaxies by Heinis et al. (2014; dashed line), star-forming galaxies valid for redshift range 0.5−4 by Pannella et al. (2015; dotted line) and direct LBGs detections in PACS by Oteo et al. (2013; blue open square).

In the text

Fig. 11 Results for stack 1 (green) and stack 2 (red) in the IRX-βUV and IRX-M∗ plane. Left panel: IRX-βUV diagram, for reference see Fig. 9. The stack 1 results are shown as diamonds, the tonalities of green represent the increase of the M∗ from 9.75 (light green) to 11.50 (dark green). The filled red circles are the results from stack 2, the tonalities of red represent the two different bins in M∗, 9.75 < log (M∗ [ M⊙ ] ) < 10.5 (light red) and 10.5 < log (M∗ [ M⊙ ] ) < 11.5 (dark red). Right panel: IRX-M∗ diagram, for reference see Fig. 10. The stack 1 results are shown as diamonds, the tonalities of green represent the two different bins in βUV, −1.7 <βUV< −0.5 (light green) and −0.5 <βUV< 0.5 (dark green). The filled red circles are the results from stack 2, the tonalities of red represent the increase of the βUV from −1.7 (light red) to 0.5 (dark red).

In the text

Fig. 12 IRX versus LFUV.The right-hand axis shows the equivalent attenuation in the FUV band in magnitudes using Eq. (5). Our data points are shown as red squares. We also show previous measurements at z ~ 3 from stacking of UV-selected galaxies in the COSMOS field by Heinis et al. (2014; cyan open squares) and LBGs detected in PACS by Oteo et al. (2013; blue open square). The dashed lines represent the region where LIRG and ULIRG lie.

In the text

Fig. 13 Total SFR versus stellar mass for our stacking sample as a function of stellar mass. The red filled squares show our results. We compare them with previous stacking studies using different sample selections at redshift z ~ 3 for star-forming galaxies (Karim et al. 2011 and Schreiber et al. 2015), UV selected galaxies (Heinis et al. 2014), and LBGs (Coppin et al. 2015); they are all plotted with open symbols. The various lines show previous observed sequences at z = 0 (Peng et al. 2010), z = 1Elbaz et al. (2007), z = 2Daddi et al. (2007), and z = 3 (from Magdis et al. 2010c, based on IRAC-detected LBGs). When necessary, the SFRs and stellar mass values have been converted to a Chabrier (2003) IMF.

In the text

	Fig. 14 Comparison between SFRtotal and SFRcorr as a function of M∗. The red filled squares show the SFRtotal from the stacking results as a function of M∗. The gray point are the individual SFRcorr values for each LBG of our sample. The blue filled squares show the average of the SFRcorr; the error bars represent the error on the mean (σ/$\sqrt{{\mathit{N}}_{\mathrm{objects}}}$).
In the text