J Álvarez-Márquez

D Burgarella

V Buat

O Ilbert

P G Pérez-González

Rest-frame FUV-to-FIR view of Lyman break galaxies at z ∼ 3: dust attenuation

Keywords: Galaxies, starburst -Ultraviolet, galaxies -Infrared, galaxies -Galaxies, high-redshift -Cosmology, early universe

Introduction

About twenty years ago, the Hubble Space Telescope (HST) allied to the 10-m class ground-based telescopes opened up a window on the first 2 Gyr of cosmic times (Madau et al. 1996). Lyman-break galaxies (LBGs) represent the largest sample of star-forming galaxies known at high redshift (z > 2.5) owing to the efficiency of their selection technique in deep broadband imaging survey (Lyman-break or dropout technique, Steidel et al. 1996). They have been a key population to investigate the mass assembly of galaxies during the first Gyrs of the universe (Shapley et al. 2001;Somerville et al. 2001;Madau et al. 1996;Steidel et al. 1996;Giavalisco 2002;Blaizot et al. 2004;Shapley et al. 2005;Baugh 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;Bouwens et al. 2015;Oesch et al. 2015;Roberts-Borsani et al. 2016;Stefanon et al. 2017;Oesch et al. 2018).

The rest-frame ultraviolet (UV) light, mostly emitted by young and massive stars, has been commonly used as a starformation rate (SFR) tracer. However, the interstellar dust scatters and/or absorbs the UV light, hence only a fraction of the energy output from star-formation is observed in the UV. The rest is re-emitted in the infrared (IR) by the heated dust. The gen-eral picture of star-formation activity across cosmic time peaks around z ∼ 2 -4, and significantly drops at z > 4 (see Madau & Dickinson 2014 for a review). The obscured star-formation dominates the total star-formation rate density (SFRD) over the redshift range 0 < z < 3, and corresponds to half of the SFRD at z = 3.6 [START_REF] Burgarella | [END_REF]. Then, it is necessary to combine the UV and IR emission to determine the complete energy budget of star-formation. However, due to the limited sensitivity of far-infrared (FIR) and radio observations, most of the current information at redshifts z > 3 -4 are obtained from UV observations of LBG samples that need to be corrected for the dust attenuation (Madau & Dickinson 2014, more references).

Only a few LBGs have been directly detected at z ∼ 3 in the mid-infrared (MIR) and FIR (Magdis et al. 2012;Casey et al. 2012;Oteo et al. 2013a), thanks to Spitzer and Herschel. This sample is biased against submillimeter (submm) bright galaxies, and are not representative of the LBG population in terms of stellar mass (M * ), dust attenuation, and SFR (Burgarella et al. 2011;Oteo et al. 2013a). The recently gain in sensitivity with the Atacama Large Millimeter Array (ALMA) at submm and millimeter wavelengths provides some insight on the obscured SFR at high-z (Capak et al. 2015;Bouwens et al. 2016;Fudamoto et al. 2017a;McLure et al. 2018). Even so, it is still quite difficult to obtain large samples to carry out statistical and representative analyses. Therefore, stacking analysis techniques (Dole et al. 2006) have been applied to relatively large samples of LBGs to derive their FIR/submm emission (Magdis et al. 2010b,a;Rigopoulou et al. 2010;Coppin et al. 2015;Álvarez-Márquez et al. 2016;Koprowski et al. 2018).

Due to the lack of FIR/submm information of individual LBGs, empirical recipes are used to correct the UV emission for dust attenuation. The most commonly adopted is the relation between the UV continuum slope (β UV , Calzetti et al. 1994) and the so-called IR-excess (IRX ≡ L IR /L FUV ) calibrated on local starburst galaxies by Meurer et al. (M99 hereafter;1999). Despite its general use to infer the dust attenuation at high-z, there are several complications that make the methodology quite uncertainty. The β UV is known to be sensitive to the intrinsic UV spectra of galaxies -which depends on the metallicity (Z), age of the stellar population, and star formation history (SFH) -and the shape of the dust attenuation curve. Stacking analysis of LBG and star-forming galaxies at redshifts of 1.5 < z < 5 have showed different behaviors, where their mean results follow the M99 relation (Magdis et al. 2010a;Reddy et al. 2012;Koprowski et al. 2018;McLure et al. 2018), lie above (Coppin et al. 2015;Bourne et al. 2017), andbelow (Álvarez-Márquez et al. 2016;Bouwens et al. 2016;Reddy et al. 2018). These deviations have been showed to be driven for the stellar masses (Álvarez-Márquez et al. 2016;Bourne et al. 2017), the shape of the dust attenuation curve (Salmon et al. 2016;Lo Faro et al. 2017), and the sample selection (Buat et al. 2015). Direct detections of individual LBGs and IR-selected sources suggest a large scatter in the IRX-β UV plane (e.g. Oteo et al. 2013a;Casey et al. 2014;Fudamoto et al. 2017b). In addition, the stellar mass has been found to correlate with the dust attenuation, and its relation seems not to evolve with redshift (Reddy et al. 2010;Buat et al. 2012;Heinis et al. 2014;Whitaker et al. 2014;Pannella et al. 2015;Álvarez-Márquez et al. 2016;Bouwens et al. 2016;Whitaker et al. 2017), at least in the intermediate range 1.0 < z < 4.0.

With the objective to deeply investigate the dust attenuation in LBGs at high-z, a sample of 22,000 LBGs in the redshift range 2.5 < z < 3.5 were selected in the COSMOS field (Scoville et al. 2007) and presented in a first publication (Álvarez-Márquez et al. 2016, AM16 hereafter). A statistically-controlled stacking analysis from the FIR to mm wavelengths (100µm to 1.1mm) was applied to derive their full IR spectral energy distribution (SEDs). Thanks to the large LBG sample, the stacking analysis was performed in different sub-samples as a function of several parameters: L FUV , β UV , and M * . We investigated the dust properties for each sub-sample of LBGs, and studied their evolution in the IRX-β UV , IRX-M * , and IRX-L FUV planes.

In this paper, we extend the initial work by applying a new stacking analysis from optical to millimeter wavelengths. It allow us to obtain the full rest-frame far-ultraviolet (FUV) to FIR SEDs of LBGs at z ∼ 3, and derive the main physical parameters by mean of SED-fitting analysis techniques. This allows to study the amount of dust attenuation and the shape of the dust attenuation curve that best represents the LBG population at z ∼ 3. The outline of this paper is as follow. In Section 2, we describe the COSMOS field data and photometry used in our analysis. In Section 3, the definition of the LBG sample, and its associated photometric redshift, M * , L FUV and β UV are presented. In Section 4, we define the method to stack the LBG sample from Optical to FIR observation, and obtain the stacked rest-frame FUV-to-FIR LBGs SEDs at z ∼ 3. In Section 5, a SED-fitting analysis with CIGALE SED-fitting analysis code is performed on the stacked LBGs SED at z ∼ 3 to derive their physical properties and build a library of semi-empirical templates of LBGs at z ∼ 3. In Section 6, we discuss the dust attenuation results obtained for the LBG sample in terms of the IRX, and which dust attenuation curve is more likely to be representative for LBGs at z ∼ 3. Finally, Section 8 presents the summary and conclusions.

Throughout this paper we use a standard cosmology with Ω m = 0.3, Ω Λ = 0.7, and Hubble constant H 0 = 70 km s -1 Mpc -1 , and the AB magnitude system. When comparing our data to others studies, we assume no conversion is needed for stellar mass estimates between Kroupa (2001) and Chabrier (2003) IMFs. To convert from Salpeter (1955) to Chabrier (2003) IMFs, we divide M * Salpeter by 1.74 (Ilbert et al. 2010).

Data

We use the available Optical to FIR imaging from COSMOS field, and the Optical/NIR multi-color catalogue (Capak et al. 2007, version 2.0). The specific data sets are presented in the following Section.

UV/Optical/NIR data

We use optical broad-band imaging (B J , V J , r + , i + , z ++ ) from the COSMOS-20 survey (Capak et al. 2007;Taniguchi et al. 2007Taniguchi et al. , 2015)), observed at the prime-focus camera (Suprime-Cam) on the 8.2m Subaru Telescope. In combination with the YJHK s broad-band imaging released in the UltraVISTA DR2 1 (Mc-Cracken et al. 2012), and observed with the VIRCAM instrument on the VISTA telescope. The original images have been homogenized to the same Point Spread Function (PSF) with a Full Width at Half Maximum (FWHM) equal to 0.8 by Laigle et al. (2016, private communication).

We additionally use the Optical/NIR broad-band photometry (Subaru: B J , V J , g + , r + , i + , z + , and VISTA:Y,J, H and K s ) from the COSMOS multi-color catalogue (Capak et al. 2007, version 2.0). The photometry was performed using SExtractor in dualimage mode over an aperture of 3 centered at the position of the i + band detection. The COSMOS multi-color catalogue was updated by Laigle et al. (2016), where they provide a new NIR selection catalogue on the ultra VISTA-DR2 observations. However, we still use the i-band selection catalogue (Capak et al. 2007, version 2.0) because the LBG sample, used on this analysis (see Sect. 3), has been asked to be detected at V J and i + bands (AM16).

Mid and Far Infrared imaging

We use the Infrared Array Camera (IRAC: 3.6, 4.5, 5.8, and 8.0 µm) and the Multiband Imaging Photometer for SIRTF (MIPS: 24 µm) observations from the SCOSMOS survey (Sanders et al. 2007) as part of the Spitzer Cycle 2 and 3 Legacy Programs. The IRAC and MIPS observations present a 5σ flux sensitivity of 0.9, 1.7, 11.3, 14.9, and 80 µJy, respectively.

Herschel Space Observatory (Pilbratt et al. 2010) mapped the COSMOS field at 100, and 160 with the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010), and 250, 350, and 500 µm Photometric Imaging Receiver (SPIRE;Griffin et al. 2010) as part of the Evolutionary Probe survey (PEP; Lutz et al. 2011) and Herschel Multi-Tiered Extragalactic Survey (HerMES, Oliver et al. 2012). The PACS images ( 100Álvarez-Márquez, Burgarella, Buat et al.: FUV-to-FIR view of LBGs at z ∼ 3 and 160 µm) present point-source sensitivities of 1.5 mJy and 3.3 mJy and a PSF FWHM of 6.8 and 11 . The SPIRE maps (250, 350, and 500 µm) have PSF FWHM of 18.2 , 24.9 , and 36.3 , 1σ instrumental noise of 1.6, 1.3, and 1.9 mJy beam -1 , and 1σ confusion noise of 5.8, 6.3, and 6.8 mJy beam -1 (Nguyen et al. 2010). These maps have been downloaded from HeDaM2

LBG sample

The LBG sample for this analysis is a sub-sample of the one selected and characterized by AM16. The original LBG sample was selected in the COSMOS field by means of the classical U-dropout technique (Steidel et al. 1996), and using the broadband filters u * , V J and i + . It contains ∼22000 LBGs in the redshift range 2.5 < z < 3.5, and UV luminosities log(L FUV [L ]) > 10.2.

One of the objectives in this analysis is to perform a consistent rest-frame FUV to FIR stacking analysis. However, the available UltraVISTA survey (Y JHK s bands) covers a reduced area of the COSMOS field: 1.5 deg 2 . To be fully consistent and stack the same LBGs at any wavelength, we restricted the original LBGs sample to the limited area given by the UltraV-ISTA survey. Therefore, our final LBGs sample contains ∼17000 LBGs, which are all LBG from the original sample enclosed in the UltraVISTA area.

Photometric redshift and stellar mass

We use the photometric redshifts (photo-z) and stellar masses computed by Ilbert et al. (2009, version 2.0) for i-band detected sources in the COSMOS field. The photo-z in the range 1.5 < z < 4 were tested against the zCOSMOS faint sample and faint DEIMOS spectra, showing an accuracy of ∼3% (Ilbert et al. 2009, version 2.0). The stellar masses were derived by applying SED-fitting techniques to the available optical and NIR photometry. They assume Bruzual & Charlot (2003) single stellar population (SSP), exponentially declining star formation history (SFH), and Chabrier (2003) initial mass function (IMF). Ilbert et al. (2010) suggest that the i-band photometric redshift catalog is 90% complete at 5µJy and 50% complete at 1µJy in the IRAC 3.6µm images.

Definition of L FUV and β UV for each LBG

The UV continuum slope, β UV , is defined following the methodology presented by Finkelstein et al. (2012). They performed a rest-frame UV to optical SED-fitting analysis, and obtained the β UV by fitting a power-law to the derived best-fit synthetic spectral model. In our case, we limit our LBG SEDs to the available optical-to-NIR broad-band photometry (optical: B J , V J , g + , r + , i + , z + and NIR: Y,J, H and K s ) that lie at the rest-frame wavelength range 1000 < λ rest-f rame [Å] < 3500. It allows to perform an SED-fitting analysis over the same rest-frame spectral range for each LBGs independently of its redshift.

We use CIGALE (Code Investigating GALaxy Emission, Burgarella et al. 2005;Noll et al. 2009;Boquien et al. 2019, see Sect. 5) to perform the rest-frame UV-to-Optical SED-fitting analysis and derive the best-fit model spectrum for each LBG of the sample. We assume Bruzual & Charlot (2003) synthetic stellar population libraries and Chabrier (2003) IMF. We vary the metallicity (0.02 < Z [Z ]< 1.00 ), age of the stellar population (1 Myr < t < t H ), dust extinction (0 < A V < 2 mag, using the Calzetti et al. (2000) dust attenuation law) and star formation history (SFH ∝ e -t/τ , with τ = 0.1, 10, 100, 10 3 , 10 4 , 10 5 , -300, -10 3 , -10 4 Myr).

β UV is calculated by power-law fit ( f λ ∝ λ β UV ) to the bestfit synthetic stellar population model that has been derived by CIGALE. The power-law fit is performed at rest-frame wavelength range 1250 < λ rest-f rame [Å] < 2000 following the spectral windows defined by Calzetti et al. (1994). We excluded the range 2000 < λ [Å]< 2600 for two reasons; i) to omit the dust feature at 2175 Å and ii) to have a homogeneous rest-frame wavelength range independently of the redshift of the galaxy. L FUV is computed in a squared bandpass defined as the GALEX FUV-band (λ e f f = 152.8 nm and effective bandwidth = 11.4 nm). We check the uncertainties in the β UV and L FUV estimations using the mock analysis tool from CIGALE (see Sect. 5.2 for details of the mock analysis tool). We obtain a dispersion in the β UV of ∆β UV ∼ 0.2, and a really robust L FUV measure with ∆ log(L FUV [L ]) ∼ 0.03. Additionally, we obtain a systematic offset between the β UV calculated by power-law fit directly to the photometry (AM16), and the one derived from power-law fit to the best-fit model obtained by SED-fitting analysis of the photometry (see Appendix A for a detailed discussion).

Stacking analysis

The stacking analysis is a technique to determine the mean flux density of a population of sources that are individually too dim to be detected in a confusion-limited maps (e.g., Dole et al. 2006;Marsden et al. 2009;Béthermin et al. 2012;Heinis et al. 2013). Our LBG sample is complete at V J and i + bands, where the sample has been selected by mean of the U-dropout technique. However, the detection of individual LBGs at other optical/NIR bands are dropped up to the 70% of completeness, and it gets worse at longer wavelength as far as only few of them could be detected at FIR observations (SPIRE 250µm). Then, we decided to average the signals of multiple LBGs from Optical-to-FIR observations to obtain fully consistent mean rest-frame FUV-to-FIR SEDs of LBGs at z ∼ 3. This allows to statistical detect the LBG population at redder wavelength at expenses of averaging over their individuals properties.

Previous works that computed average rest-frame FUV-to-FIR LBG SEDs (e.g. Magdis et al. 2010a), performed a simple average to the Optical/NIR photometry. Our LBG sample is 70-80% complete at NIR, and 90-100% at Optical bands. Then, we decided to perform a stacking analysis in the Optical/NIR to eliminate a possible bias to the brightest LBGs (if only detected ones are averaged).

We split our LBG sample as a function of L FUV , β UV , and M * in five different ways. The first three are done as a function of the single parameters L FUV , β UV , and M * . For the two others, we split the sample as a function of the combination of both, β UV and M * , in the (β UV , M * ) plane. Table 1 shows the sample binning. It includes the interval used to define each bin, the number of LBGs that are stacked, and the mean values of L FUV , β UV , M * , and redshift that define each bin. The uncertainties associated to the average values are defined by the standard deviation of the mean. The number of bins and the intervals are optimized to obtain a good signal-to-noise ratio (SNR) on the final FIR (SPIRE bands) stacked images.

General method

The general method used to perform the optical-to-FIR stacking analysis is based on the IAS library (Bavouzet 2008 and Béther- 2008;Béthermin et al. 2010;Heinis et al. 2013;AM16). The three points below are of major importance:

1. Correction for clustering of the input catalogue. LBGs are clustered between them and other star-forming galaxies at high redshift. Depending on the PSF size of the available observations, the PSF beam may be contaminated by different nearby sources included in the same sample. Therefore, we could account for the flux of the same source more than once during the stacking analysis. This effect would lead to an overestimation of the mean flux due to the clustered nature of the sources. It is corrected by taking into account the angular correlation function of the input catalogue. The complete formalism is detailed in AM16. 2. Correction 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). If the stacked sample presents a relatively high percentage of incompleteness, part of the faint LBG population (mostly located in the crowded areas of the image 3 http://www.ias.u-psud.fr/irgalaxies/downloads.php or close to bright objects) are missed during the source extraction process at optical images. When we move to longer wavelength observation (e.g. FIR), the optical crowded areas and the bright objects contribute to enlarge the local background level. This means that if we perform the stacking analysis on an incomplete population, we would stack mostly the areas where the background presents lower levels and lose the contribution of the areas with higher ones. It produces a negative contribution in the central part of the stacked image related to the global background level. We corrected this effect using extensive simulations to characterize the detection process and sample selection. We performed these simulations by injecting mock sources in the original images (V J and i + bands), and keeping track of the recovered sources. The full formalism used to correct this effect is deeply explained in AM16.

3. Correction for bright field sources in crowded field observations. A deep observation, not dominated by the confusion noise, generates a crowded field image. The presence of a large amount of bright neighbors may generate an additional noise level that will perturb the final stacked image. We remove the contribution of field sources in our stacking procedure by performing a weighted stacking analysis to mask them. We proceed as follow:

-We generated a segmentation map, using SExtractor (Bertin & Arnouts 1996), on the image where the stacking analysis will be performed. A source extraction was applied to detect all the objects with flux larger than 3σ Back , excluding the LBGs from our sample.

Álvarez-Márquez, Burgarella, Buat et al.: FUV-to-FIR view of LBGs at z ∼ 3

-A weight map is created by assigning a value equal to 1 for background pixels, and equal to 0 for detected source pixels.

-We performed the stacking analysis by using a weighted mean, i I i W i / i W i , for each individual pixel of the cutout image. Where i corresponds to each object inside the sample, I i one of the pixels from the cutout image extracted from the map, and W i one of the pixels from the cutout image extracted from the weighted map. The left panel shows the stacked image using the weighted stacking analysis and the right one the general stacking method. The bottom panel shows the radial profiles estimated from the upper images. The blue line corresponds to the stacking analysis performed using the general method, and the red line correspond to the stacking analysis performed using the weighted method. The lighter blue and red points correspond to the value of each pixel from the stacked images.

Figure 1 shows the improvement that we get in the final stacked image by using a weighted stacking analysis for the 24µm MIPS observations. This example corresponds to the final stacked image, with and without applying the weighted procedure, for one of the bins of the stacking analysis as a function of L FUV . The new methodology provides a considerable improvement in the background noise, which is essential to detect and recover the emission of the stacked population.

Each of these corrections impacts at various levels on the final stacked image, which mainly depends on the image quality (PSF, confusion noise, depth), observed wavelength, and sample selection method. In Sections 4.2, 4.3, and 4.4, we will discuss which of the previous corrections are taken into account during the stacking analysis for each of the data set.

FIR analysis

The following procedure is applied to perform the stacking analysis on PACS (100 and 160 µm) and SPIRE (250, 350 and 500 µm) images. It was previously used and deeply explained in AM16. We use the calibrated PACS and SPIRE maps without cleaning the images from any detected source, and we do not apply the correction for bright field sources in crowded field observations to mask the bright neighbors during the stacking procedure. If a large number of sources were stacked in confusion limit maps, as in our FIR images, the bright neighbors will tend to average homogeneously to a constant background level. We demonstrate this by performing a random stacking analysis at the same number of positions than objects in each of our sample bins.

PACS and SPIRE observations present large PSF sizes, with FWHM PSFs that vary from 6 to 36.3 . Our LBG sample is clustered, and could be defined by an auto-correlation function, ω[θ, φ] ∝ θ -γ , with γ = 0.63 (AM16). The combination of large PSF sizes and the clustered nature of our LBG sample could produce a contamination of our LBG PSF beams due to a nearby LBGs. If we stack them, we will account for the flux of the same source more than once during the stacking analysis. Therefore, we applied the correction for clustering of the input catalogue to deconvolve the emission of our LBG population and the clustering contribution. Additionally, we corrected our stacking analysis for the incompleteness of the input catalogue in the dense regions following the analysis done in AM16.

Mid-Infrared analysis

The following procedure will be applied to perform the stacking analysis for IRAC (3.6, 4.5, 5.8, 8 µm) and MIPS (24 µm) observations. The IRAC images present a PSF FWHM equal to 1.66 , 1.72 , 1.88 , and 1.98 with a pixel size of 0.6 , and the MIPS observation has a PSF FWHM equal to 6 and a pixel size of 1.2 . The improvement of the PSF FWHM size and sensitivity from FIR to Mid-IR observations reduces the confusion noise and improves the quality of the Mid-IR images. It generates a crowded image with a large number of field sources. The presence of the large amount of field sources in the maps will generate an additional noise level and clustering contribution that will perturb our stacked images. To remove this effect from our stacked sources, we applied the weighted stacking analysis previously presented in the Correction for bright field sources in crowded field observations.

The reduction of the PSF FWHM sizes from FIR to MIR maps lowers and makes negligible the clustering effect on the stacked LBGs. Therefore, here, we did not take into account any corrections related to the clustering nature from our input LBGs population in the stacked LBGs. However, the bias induced in our stacked image due to the incompleteness of our input catalogue is still present. We realized that this effect decreases when the stacking procedure is performed closer to the sample selection wavelengths (V J and i + bands). However, it is still important to correct it in the MIR wavelength.

We perform an aperture photometry to calculate the flux density of the stacked LBGs in IRAC and MIPS observations. We make a circular aperture with radius equal to 3.6 and 7 , and a background annulus from 3.6 to 8.4 and 20 to 32 . We then apply an aperture correction of 1.125,1.120,1.135,1.221,and 1.61 at IRAC (3.6,4.5,5.8,8 µm) and MIPS (24 µm) 4 . Bootstrap resampling is used to obtain the mean values and errors. We repeat the above procedure over 3000 random bootstrap samples, and the 1σ of the distribution of the derived fluxes is adopted as the uncertainty for our results. Appendix C presents the derived fluxes and uncertainties of our stacking analysis at FIR and MIR. It includes the stacking analysis performed for AzTEC (1.1mm), and previously presented in AM16.

Optical/NIR analysis

The following procedure is applied to perform the stacking analysis in the optical (B J , V J , r + , i + , z ++ bands) and NIR (Y, J, H, K s ) observations. The optical and NIR images have an homogenized PSF with a FWHM equal to 0.8 , and a large depth improvement with respect to the MIR observations. The optical/NIR observations contain large number of bright field sources that are removed by performing a weighted stacking analysis. We did not apply any corrections related to the clustering nature of our sample and the incompleteness of our input catalogue. These two effects are negligible in our stacked optical/NIR images.

We computed the photometry on the stacked object by fitting a Moffat function. The error is calculated by combining the contribution of different effects; bootstrap, error in the fitting procedure, and Poisson noise. Appendix B presents the derived fluxes and uncertainties of our stacking analysis in Optical and NIR wavelengths. We verify the methodology by comparing the derived stacked fluxes with an average observed photometry of individuals LBGs at V J and i + bands, where the sample is complete. In both cases, we obtain results within the derived uncertainties.

SED-fitting analysis

CIGALE (Burgarella et al. 2005;Noll et al. 2009;Boquien et al. 2019) 5 is a Python-based modular code dedicated to fit SEDs from the UV to the FIR and to create galaxy emission models in the same wavelength range. CIGALE is modular and it allows the user to create models by selecting independent modules for, e.g. the star formation history, the stellar models, the shape of the dust attenuation curve, the nebular emission, the contribution of an active galactic nuclei and the dust emission templates or models.

Initial parameters

Deriving physical properties of galaxies by fitting synthetic SEDs to a multi-band photometry requires making prior assumptions on the initial parameters to build the synthetic SED library. Assumptions such as SFH, IMF, metallicity, and dust attenuation curve, affect the derived galaxy properties, like the SFR, stellar mass, age of the stellar population, amount of dust attenuation, 4 The photometry configuration and aperture correction have been obtained from the instrument handbook of IRAC (http://irsa.ipac.caltech.edu/data/SPITZER/docs/ irac/iracinstrumenthandbook/27/), and MIPS (http: //irsa.ipac.caltech.edu/data/SPITZER/docs/mips/ mipsinstrumenthandbook/50/) 5 Code Investigating GALaxy Emission (CIGALE): https:// cigale.lam.fr/ etc (e.g. Papovich et al. 2001;Maraston et al. 2010;Pforr et al. 2012;Schaerer et al. 2013;Buat et al. 2014;Lo Faro et al. 2017). Table 2 summarizes the initial parameters used to perform the SED-fitting analysis for our stacked LBGs SEDs at z ∼ 3. We adopt the stellar population models from Bruzual & Charlot (2003) and an IMF from Chabrier (2001). Delayed SFHs (SFR ∝ t/τ 2 e -t/τ ), with varying e-folding times are assumed to model the SFH of our LBGs at z ∼ 3 (Lee et al. 2010). We perform a test to check the influence of the chosen SFH on the final derived physical parameters, where we run CIGALE with different SFHs recipes (exponential declining, exponential rising, constant, and delayed) using the same initial parameters from Table 2. The results suggest that our SED-fiting analysis cannot distinguish which of the tested SFHs recipes are more likely to reproduce the LBG population. Nevertheless, our SED-fitting analysis can derive fully consistent physical parameters independently of the SFH chosen, with the exception of the age of the stellar population.

The metallicity is fixed to 0.2 Z (Castellano et al. 2014). The choice to fix the metallicity is motivated by the fact that: i) the SED-fitting technique applied to broad-band photometry can hardly constrain the metallicity due to the age-dust-metallicity degeneracy (López Fernández et al. 2016), and ii) the metallicity of our stacked SEDs are averaged due to the combination of a large number of galaxies. However, there is a possibility to have an evolution of metallicity in our LBG sample due to the well known mass-metallycity relation (e.g. Mannucci et al. 2009). Then, we perform an individual SED-fitting analysis for the stacking analysis as a function of M * with the metallicity as a free parameter. The derived physical parameters are in agreement, within the uncertainties, with the ones from the fixed metallicity. Therefore, we concluded that our choice of fixed metallicity does not have a large impact on the final derived physical parameters.

Dust attenuation is treated using the prescription given by Noll et al. (2009) Eq. 1. They use the Calzetti dust attenuation law (k(λ), Calzetti et al. 2000) as baseline, and give the possibility to: i) include a UV bump at 2175 Å modeled by a Lorentzianlike profile, and ii) modify the slope of the dust attenuation curve by a power law.

A(λ) = A V 4.05 k(λ) + D λ o ,γ,E bump λ λ V δ (1)
where A(λ) is the modified dust attenuation curve presented by Noll et al. (2009), A V is the dust attenuation in the V-band, D λ o ,γ,E bump is the UV bump profile, δ is the power-law slope with respect to the Calzetti one, and λ V = 5500Å. The shape of the dust attenuation curve depends on the star-dust geometry, grain size distribution, etc. (Witt & Gordon 2000a). Complex dust geometries can produce a grayer or shallower dust attenuation curve than the Calzetti one (Buat et al. 2011a;Kriek & Conroy 2013;Zeimann et al. 2015;Salmon et al. 2016;Lo Faro et al. 2017). Our SED fitting takes δ as a free parameter to mimic the variation on the shape of the dust attenuation curve. We did not include any absorption feature produced by dust at 2175 Å (UV bump), because our stacked SEDs do not have enough spectral resolution.

The IR emission is modeled using Draine & Li (2007) dust models. Their validity to reproduce the FIR emission of Herschel detected high-z main sequence galaxies has been confirmed by Magdis et al. (2012). The dust models require the fine tuning of several free parameters; q PAH , U min , U max , α, γ and M dust . However, Draine & Li (2007) showed that fixed values of α = 2 and U max = 10 6 can reproduce the SEDs of galaxies with a wide range of properties. We fix the q PAH value to 3.9%, which is the one that best fit our stacked LBG SEDs in average. This value is higher that the one proposed by Magdis et al. (2012), who found an interval of q PAH values from 1.12% to 3.19%. We consider U min varying between 0.7 to 50. Draine & Li (2007) suggest using U min > 0.7 because lower values will correspond to dust temperatures bellow ∼ 15K, which we cannot constrain by our FIR photometry alone. The fraction of dust enclosed in the PDRs was fixed by (Magdis et al. 2012, γ=0.02), and Draine & Li (2007) suggested that γ value should not be higher than 0.04 (4%). We detect some variation at rest-frame wavelength 20-50µm for our set of stacked SEDs, which suggest different γ values. Therefore, we use a γ range from 0.01 to 0.04. The adopted initial configuration (Table 2) gives a total of 609840 synthetic galaxy models/SEDs. This set of models is the one used to perform the SED-fitting and Bayesian analysis over our thirty rest-frame FUV-to-FIR stacked LBGs SEDs at z ∼ 3.

Physical parameters for LBGs at z ∼ 3

We run CIGALE for the sample of thirty rest-frame FUV-to-FIR stacked LBGs SEDs at z ∼ 3. We use the adopted initial configuration listed in Table 2. CIGALE estimates and derives the physical parameters, such as the stellar mass, dust mass, age of stellar populations, SFR, rest-frame IR and FUV luminosities, FUV at-tenuation, and β UV , as well as the input parameters from Table 2. Additionally, CIGALE provides the synthetic model that best fits the observational SEDs by χ 2 minimization.

CIGALE uses a Bayesian analysis to derive the physical parameters. It builds the probability distribution function (PDF) for each parameter, and takes its mean and standard deviation (Noll et al. 2009;Boquien et al. 2019). We use a mock analysis to check the robustness, accuracy, and parameter degeneracy of our estimates. The mock analysis consists in generating a mock catalogue of artificial SEDs, similar to the input LBGs, for which the physical parameters have been previously derived, addd a Gaussian noise where the dispersion is taken from the observed uncertainty. Then, we compare with CIGALE estimates (more details: Buat et al. 2012;Ciesla et al. 2015;Lo Faro et al. 2017;Boquien et al. 2019). This procedure allows to check our ability in estimating/constraining the individual output parameters given the information provided by the detailed shape of the broad band SED. The results are summarized in Figure 2. We call true the parameter value provided by the best-fit model and used to compute the mock SEDs. This true parameter value, on the x-axis, is compared with the PDF one, on the y-axis, which is computed by SED-fitting analysis over the mock SEDs catalogue. From top-left to bottom-right, Figure 2 shows the results for the age of the stellar population, e-folding timescale of delay-τ SFH, SFR, FUV dust attenuation, E(B-V), change of the slope of the dust attenuation law with respect to Calzetti -δ, stellar mass, total IR luminosity, and U min . On the one hand, the results emphasize that we are not able to determine the shape of the SFH and age of the stellar population by SED-fitting analysis. The estimated PDFs of the e-folding timescale, τ, appear to average out to a constant value independently of the true value used to compute the mock SEDs. The age of the stellar population is slightly constrained for the particular case of delay-τ SFH model, but it is biased against the chosen SFH. On the other hand, the results suggest that our procedure well constrain the physical parameters: stellar mass, IR luminosity, FUV dust attenuation, SFR, and change of the slope of the dust attenuation law with respect to Calzetti -δ, with the chosen delay-τ SFH. In Section 5.1, we also conclude that these physical parameters are not affected by the chosen SFH. Therefore, despite of the degeneracy for the SFH and ages of the stellar population, our SED-fitting analysis provides consistent physical parameters.

Appendix D lists the derived input and additional physical parameters by SED-fitting analysis for our thirty stacked LBGs. The β UV and L UV are computed in a separate FUV-Optical SEDfitting analysis, following the procedure already explained in Section 3.2.

Templates of LBGs at z ∼ 3

We use the results from the panchromatic stacking and SEDfitting analyses to build a library of empirical rest-frame FUVto-FIR templates of LBGs at z ∼ 3. The library contains thirty empirical SEDs and their best-fit synthetic model spectrum. On the one hand, the empirical SEDs correspond to the rest-frame FUV-to-FIR stacked LBGs SEDs derived by stacking analysis in Section 4. They are the result from combining ∼17000 individual LBGs in different bins as a function of L FUV , M * , β UV , and a combination of both M * and β UV . More specifically, the number of individuals LBGs used to obtain each stacked LBG SED varies from about 9000 to 100. This statistics strengthens our library in the sense that our SEDs have an exceptional SNR, which is not comparable with any individual detected LBG at z ∼ 3. The set of synthetic model spectrum associated to each stacked LBG SED are obtained with CIGALE.

Thanks to the binning configuration used to split our large LBG sample, we have derived a set of templates with a large variety of physical properties. Our library contains templates within an interval of stellar mass, 9.2 < log(M * [M ]) < 11.4, SFR, 20 < S FR [M yr -1 ] < 300, β UV , -1.8 < β UV < -0.2, FUV dust attenuation, 1.5 < A FUV [mag] < 4.8, IR luminosity, 11.2 < log(L IR [L ]) < 12.7, and FUV luminosity, 10.4 < log(L FUV [L ]) < 11.2. This variety makes our library quite versatile, in the sense that our templates fit a wide range of physical properties for star-forming galaxies at z ∼ 3.

Figure 3 shows the four SEDs and templates derived in the stacking analysis as a function of L FUV . If a sample binning is done as a function of L FUV , the sub-sample will be a mix of red, blue, high, and low mass galaxies. For example, a massive red galaxy with a large dust attenuation may have the same L FUV than a low-mass blue galaxy. Therefore, if a stacking analysis is applied to a mix population of LBGs, the final results will tend to average to a similar mean SED with a different normalization factor, the L FUV . This is what we can see in Figure 3, the rest-frame FUV-to-NIR part of the SEDs presents a similar behavior except for the lowest L FUV bin, which has a slightly bluer β UV . We can also see this from the derived physical properties on the SED-fitting analysis, where we only find differences in the bolometric luminosity, M * and SFR, which is associated with an increase of the normalization factor. However, The FIR peak has progressively shifts to longer wavelengths for higher L FUV , which suggests a lower mean dust temperature. If we consider that the U min values are correlated with the location of the FIR peak (Draine & Li 2007;Magdis et al. 2012). This behavior is reinforced by the variation of the U min values (38.5 < U min < 19.2), derived from our SED-fitting analysis, for the lowest to the highest L FUV bin. as a function of β UV . Their rest-frame FUV-to-NIR SEDs have higher rest-frame NIR flux for redder β UV , and the slope of the UV spectrum follows the mean β UV of the stacked population. The FIR emissions present similar shapes with an increase of L IR for redder galaxies, except for the bluest β UV that give an IR peak strongly shifted to longer wavelength. Blue galaxies are quite faint on FIR, and we suggest that our FIR stacking analysis on LBG-β1 is quite uncertainty. Besides, if this IR behavior is real we should see it on LBG-Mβ1 stacking analysis, and this is not happening.

The bottom-left panel of Figures 4 illustrates the six SEDs and templates derived from the LBG-Mβ1 stacking analysis. As a reminder, we split the sample in two large bins of M * , and each of them in three bins of β UV . Their rest-frame FUV-to-NIR SEDs seem to sub-divide into two different sets of templates, high and low M * , with variations in the UV part of the SED. No variations in the shape of the FIR part is found, except for an increase of the IR luminosity with M * and β UV . The bottom-right panel of Figure 4 shows the ten SEDs and templates derived from the LBG-Mβ2 stacking analysis. In this case, we split the sample in two large bins of β UV , and each of them in five bins of M * . Their rest-frame FUV-to-NIR SEDs present similar shapes in comparison with the stacking analysis as a function of M * : an increase of the rest-frame NIR flux and a reddening of UV/Optical part of the spectrum for massive LBGs. However, the split in β UV , red and blue bins, gives access to two templates for each M * , with similar NIR shape, but strongly different UV part of the spectrum. The FIR shape does not present a large variety.

Dust attenuation

We analyzed in AM16 how the dust attenuation varies as a function of β UV and M * for our LBG sample by mean of the IRX. At that time, we had access to the full FIR/sub-mm part of the spectrum to compute the average L IR . However, the β UV and M * values, associated to each stacked LBG SEDs, were derived by averaging over the individual values of each LBG inside the stacked bin. The new rest-frame FUV-to-FIR stacking analysis allows to derive the β UV and M * directly from the stacked LBG SEDs. We then review the information in the IRX-β UV and IRX-M * planes using the results obtained from the rest-frame FUVto-FIR SED-fitting analysis on the stacked LBGs SEDs.

Dust attenuation as a function of β UV

Figure 5 shows the IRX-β UV plane, where we plot the results from our rest-frame FUV-to-FIR stacking analysis as a function of β UV (LBG-β). This is compared with the well-know IRX-β UV relation calibrated from local starburst galaxies (M99), and the one re-computed by Takeuchi et al. (2012) using the same sample and new photometric data from AKARI and GALEX. We also include our previous IRX-β UV relation derived from a FIR stacking analysis over the same LBG sample (AM16). In addition, we plot the IRX-β UV results from LBGs at redshifts, z ∼ 3 (Koprowski et al. 2018) and 2 < z < 3 (Bouwens et al. 2016), star-forming galaxies at redshifts, 2 < z < 3 (McLure et al. 2018) and 1.5 < z < 2.5 (Reddy et al. 2018), and a mass selected sample at redshift z ∼ 3 (Bourne et al. 2017). Takeuchi et al. (2012, dashed line), and the one derived from a FIR stacking analysis over the same LBG sample (AM16, dotted blue line). In addition, we plot the IRX-β UV results from LBGs at redshifts, z ∼ 3 (Koprowski et al. 2018, orange dots) and 2 < z < 3 (Bouwens et al. 2016, gray dots), star-forming galaxies at redshifts, 2 < z < 3 (McLure et al. 2018, green dots) and 1.5 < z < 2.5 (Reddy et al. 2018, black dots), and a mass selected sample at redshift z ∼ 3 (Bourne et al. 2017, blue dots). The downward triangles represent upper limits.

The IRX-β UV results from the FUV-to-FIR (this work) and only FIR (AM16) stacking analyses disagree, even if both are computed with the same LBG sample. The FUV-to-FIR results gives ∼2-3 time higher IRX value for a given β UV than the FIR ones. The difference resides on the methodology to derive the β UV used in each of the analysis. In Appendix A, we demonstrate that β UV-power (AM16) are biased to redder values in comparison with β UV-SED (this work, see Section 3.2). If the IRXβ UV results from the FIR analysis are corrected by the relation between the β UV-power and β UV-SED , obtained in Appendix A, they are in agreement with the IRX-β UV relation obtained from the FUV-to-FIR analysis. Koprowski et al. (2018) also confirmed that discrepancies in the IRX-β UV relation are due to bias in the methodology used to determine the β UV . They suggested that these inconsistencies are driven by scatter in measured values of β UV from limited photometry that serves to artificially flatten the IRX-β UV relation. They obtained that the scatter is significantly reduced by determining β UV from SED-fitting analysis, which is the same conclusion presented in Appendix A.

Our FUV-to-FIR stacking and SED-fitting analysis results as a function of β UV follows the IRX-β UV local relation (M99) within the uncertainties. They are in agreement with previ-ous stacking analyses done on star-forming galaxies at high-z. Magdis et al. (2010a) found that the dust corrected UV-SFR derived from M99 relation presents a good match with the FIR and radio SFR estimators for stacked IRAC spectropically-confirmed LBGs at z ∼ 3. Koprowski et al. (2018) by stacking a sample of ∼4000 LBGs at redshifts, 3 < z < 5, conclude that LBGs are consistent with the M99 IRX-β UV relation, and they do not present a redshift evolution. Similar results were obtained by McLure et al. (2018), who studied a sample of star-forming galaxies at redshifts, 2 < z < 3, from a deep ALMA 1.3mm continuum data. Additional stacking analysis of spectroscopicallyconfirmed z ∼ 2 LBGs in FIR (Reddy et al. 2012), z ∼ 4 LBGs in the radio continuum (1.4GHz, To et al. 2014), and direct detections of LBGs in PACS at lower redshift (z ∼ 1, Oteo et al. 2013b;Burgarella et al. 2011) also lie on the M99 relation. Bourne et al. (2017) used a deconfusion methodology to reach below the confusion limit maps (SCUBA2 and Herschel) for massive galaxies at redshifts, 0.5 < z < 6, detected in the CANDELS fields (AEGIS, COSMOS, and UDS). They obtained similar results than for local starburst galaxies in M99, although it deviates to higher IRX values at high stellar masses, which is similar to the analysis presented in Sect. 6.3. However, Reddy et al. (2018) used ∼ 3500 star-forming galaxies at redshifts, 1.5 < z < 2.5, and conclude that the M99 IRX-β UV local relation overpredicts the IRX at a given β UV by a factor of ∼ 3. They used the β UV-power definition, which may be the cause of the differences as we already commented in Appendix A. Additional ALMA observations suggested that galaxies at z ∼ 2-3 present lower IRX values than the M99 relation (Fujimoto et al. 2016;Bouwens et al. 2016).

Dust attenuation as a function of M *

Figure 6 shows the IRX-M * plane, where the results from our rest-frame FUV-to-FIR stacking analysis as a function of M * (LBG-M * ) are plotted. These are compared with the relationship presented by Bouwens et al. (2016) for star-forming galaxies at redshifts 2 < z < 3, and defined from the combination of different works done in the literature (Reddy et al. 2010;Whitaker et al. 2014;Álvarez-Márquez et al. 2016). This relation includes the one previously derived from a FIR stacking analysis over our LBG sample (AM16). We overplot the relations obtained from a stacking analysis of UV-selected galaxies at redshift z ∼ 3 (Heinis et al. 2014), and a complete sample of star-forming galaxies up to z ∼ 4 (Pannella et al. 2015). We also show the IRX-M * results from LBGs at redshifts, z ∼ 3 (Koprowski et al. 2018) and 2 < z < 3 (Bouwens et al. 2016), and star-forming galaxies at redshifts, 2 < z < 3 (McLure et al. 2018) and 1.5 < z < 2.5 (Reddy et al. 2018). The IRX results from the FUV-to-FIR stacking analysis show a good agreement at the highest stellar mass bins (log(M * [M ]) > 10). However, the IRX results disagree at low stellar masses bins (log(M * [M ]) < 10) presenting larger IRX values for a given M * .

We suggest that an incompleteness of our LBG sample in terms of M * might be the origin of the IRX-M * discrepancies. If M * is proportional to the rest-frame NIR emission (Kauffmann & Charlot 1998), the incompleteness of our LBG sample at ∼ 3 in terms of M * will be proportional to the detectability of our LBGs at IRAC-3.6µm observations. The mean fluxes of the LBG population with stellar masses, log(M * [M ]) < 10, are between 1.2 to 3.6 µJy (see Appendix C). Ilbert et al. (2010) showed that the i-band photometric redshift catalog (Ilbert et al. 2009, version 2.0), used to buid our LBG sample in AM16, is 90% complete The blue line presents the consensus relationship for star-forming galaxies at redshifts 2 < z < 3 (Bouwens et al. 2016), and defined by a combination of different works done in the literature (Reddy et al. 2010;Whitaker et al. 2014;Álvarez-Márquez et al. 2016). This relation includes the one previously derived from a FIR stacking analysis over our LBG sample (AM16, dotted blue line). We overplot the relations obtained from a stacking analysis of UV-selected galaxies at redshift z ∼ 3 (Heinis et al. 2014, dashed line), and a complete sample of star-forming galaxies up to z ∼ 4 (Pannella et al. 2015, dot-dashed line). We also show the IRX-M * results from LBGs at redshifts, z ∼ 3 (Koprowski et al. 2018, orange dot) and 2 < z < 3 (Bouwens2016, gray dot), and star-forming galaxies at redshifts, 2 < z < 3 (McLure et al. 2018, green dots) and 1.5 < z < 2.5 (Reddy et al. 2018, black dots). The downward triangle represent upper limits. at 5µJy and 50% complete at 1µJy in the IRAC 3.6µm observations. Then, the low stellar mass bins (log(M * [M ]) < 10) of LBG-M stacking analysis are between 80% to 50% complete in term of M * . However, there is some indication in the literature that the origin of the departure might not be due to M * bias. For instance, Whitaker et al. (2017) show that the region corresponding to galaxies at log(M * [M ]) < 10 and high IRX values is populated by galaxies. This effect does not seem to be at z = 0, though. It could also be noticed in Heinis et al. (2014) that there is a trend for galaxies with low-L UV to have larger IRX. We can wonder whether this effect does not hide the fact that these low-L UV only trace dust-free stars in galaxies that might be otherwise dusty (meaning that most of the stars do not contribute to the UV luminosity). Finally, Spitler et al. (2014) studied a sample of 57 galaxies with log(M * [M ]) < 10.6, and they find that massive star-forming galaxies can be found at all M UV . So, we propose that the trend we observe at log(M * [M ]) < 10 in Figure 6 might be real and not due to an observational bias. The impact of this hypothesis will be further studied in Bogdanoska et al. (in prep.). In Section 6.3, when LBGs are stacked as a function of M * and β UV , a population of LBGs appear to locate at the same region of the IRX-M * diagram than Whitaker et al. (2017). If this population is removed to perform the stacking analysis and derive the mean IRX-M * relation (light red filled square in Figure 7), the results locate closer to the IRX-M * relations from Pannella et al. (2015) and Heinis et al. (2014), and the star-forming galaxies at 2 < z < 3 stacked in ALMA observation (McLure et al. 2018).

Dispersion on the IRX-M * and IRX-β UV planes

We use our stacking analyses (LBG-Mβ) to investigate the dispersion on the IRX-β UV and IRX-M * planes, and the validity of the IRX-β UV and IRX-M * mean relations. On the LBG-Mβ stacking analyses, the LBG sample has been split on the (M * , β UV ) plane in two different ways. For LBG-Mβ1, the LBG sample is divided in two bins of M * , and each of them is split as a function of β UV . And for LBG-Mβ2, it is divided in two bins of β UV and each of them is split as a function of M * . Figure 7 shows the results from the LBG-Mβ stacking analyses in the IRX-β UV and IRX-M * planes. The results suggest large dispersion on both planes, as a consequence of the M * evolution in the IRX-β UV plane, and the β UV in the IRX-M * plane. On the one hand, high stellar mass LBGs present higher IRX values than the mean IRX-β UV relation. Similar results were reported by Bourne et al. (2017) in the IRX-β UV plane. They showed an evolution of the IRX-β UV relation from low to high stellar masses using a sample of massive galaxies at redshifts 0.5 < z < 6. Some studies of IR selected galaxy samples have found IRX-β UV values located above M99 relation (Casey et al. 2014;Oteo et al. 2013a;Buat et al. 2015), which suggest that the criterion to select the sample has a strong impact on the mean dust attenuation of a galaxy population. On the other hand, redder LBGs tend to have higher IRX values than the mean IRX-M * relation. Specially, a population of LBGs with β UV > -1 and stellar masses (log(M * [M ]) < 10) shows large IRX values that flatten the mean IRX-M * presented in Figure 6.

The dependence of M * in the IRX-β UV plane, and β UV in the IRX-M * plane, makes us think about the utility and efficiency to provide a new empirical IRX relation combining the β UV and M * . It may reduce the uncertainty to derive the dust attenuation for galaxies where UV/Optical/NIR observations are only available. We use the results from the SED-fitting and stacking analysis as a function of β UV and M * , and the combination of the two (LBG-β, LBG-M, LBG-Mβ1, and LBG-Mβ2) to attempt performing a plane fit in the 3D IRX-β UV -M * diagram. We obtain a mean relation equal to log(IRX) = (0.51 ± 0.06)β UV + (0.37 ± 0.08) log(M * [M ]) -1.89 ± 0.40. Figure 8 shows the dispersion between the IRX of each stacked LBGs obtained from SED-fitting analysis, and the IRX derived from the mean IRXβ UV -M * relation.

Shape of the dust attenuation curve for LBGs at z ∼ 3

The Calzetti law is frequently assumed when studying starforming galaxies at high-z (e.g. Bouwens et al. 2009;Schaerer et al. 2013;Finkelstein et al. 2015). From observations in the local universe, it is well known that dust attenuation curves are not universal across galaxies. Most sight-lines in the Milky Way (MW) have a strong extinction bump near 2175 Å (Stecher 1965;Cardelli et al. 1989). However, the feature is weaker in the Large Magellanic Cloud (LMC; Koornneef & Code 1981) and absent in the Small Magellanic Cloud (SMC; Prevot et al. 1984). This diversity is also reflected by systematic changes in the FUV slope, it gets steeper from MW to SMC. Calzetti et al. (1994) deduced a mean dust attenuation curve from a sample of 39 local UV-bright starburst galaxies. Their measurements are character- ized by a grayer FUV slope than both MW and LMC extinction curves, and absence of the 2175 Å absorption feature. At higher redshifts, a variety can be also seen. Buat et al. (2011bBuat et al. ( , 2012) ) found for a UV-selected galaxies at z > 1 a clear present of 2175 Å bump, and evidence for a steeper rise of the dust attenuation curve in comparison to Calzetti law. Other studies support this model, demonstrating that the 2175 Å bump is typically visible in "normal" star-forming galaxies (e.g., Burgarella et al. 2005;Noll et al. 2007;Conroy 2010). Other authors also suggest that dust attenuation curves vary significantly at any redshift (e.g. Kriek & Conroy 2013;Zeimann et al. 2015;Salmon et al. 2016;Lo Faro et al. 2017). This diversity in the shape of the dust attenuation curves could be produced by the differences in the dust grain properties, the relative geometrical distribution of stars and dust withing a galaxy, the line-of-sight galaxy orientation, and the galaxy type (Witt & Gordon 2000b;Pierini et al. 2004;Tuffs et al. 2004;Chevallard et al. 2013).

The results from the SED-fitting analysis over our stacked LBGs SEDs suggests a large variation in the shape of the dust attenuation curve across the LBG population at z ∼ 3. We obtain a range of δ's from -0.3 to 0.3 (see Appendix D), that produce dust attenuation curves steeper and grayer than Calzetti law (δ = 0), respectively. Similar results have been reported by Salmon et al. (2016), who found a range of δ's from -0.5 to 0.2 for star-forming galaxies at z ∼ 1.5 -3. Figure 9 illustrates the dependence of δ as a function of M * for the stacking analysis as a function of β UV and M * , and the combination of these both parameters (LBG-β, LBG-M, LBG-Mβ1, and LBG-Mβ2). We obtain that δ is well correlated with the M * . The low stellar mass population of LBGs (log(M * [M ]) < 10.25) favors steeper dust attenuation curves than the Calzetti one, and the large stellar mass population of LBGs (log(M * [M ]) > 10.25) favors grayer ones. Cullen et al. (2018) studied a sample of starforming galaxies at redshifts, 3 < z < 4, and stellar mass range, 8.2 < log(M * [M ]) < 10.6. Their mean FUV-to-Optical dust attenuation curve is in agreement with Calzetti's, and gets steeper for lower masses. However, Zeimann et al. (2015) performed an analysis of star-forming galaxies selected by their rest-frame op- tical emission lines at redshifts, 1.9 < z < 2.35, and mass range, 7.2 < log(M * [M ]) < 10.2. They found grayer dust attenuation curves than the Calzetti one, being in opposition to us. They use a low mass galaxies selected by detection of optical emission lines, and their results might be biased against high optical emitters galaxies that could present different behavior than our LBG population.

Figure 10 illustrates the derived dust attenuation curves from the SED-fitting analysis over the stacked LBGs SEDs obtained from the stacking analysis as a function of M * . It shows how dust attenuation curves get steeper when the stellar mass decreases. LBGs with log(M * [M ] ∼ 10 or lower present a dust attenuation curve similar to the SMC (Pei 1992). Other mean dust attenuations laws derived from UV-selected galaxies at z > 1 Buat et al. (2011b), and star-forming galaxies at z ∼ 2 to 6 in the COSMOS field Scoville et al. (2015) are also compared, and they are in agreement with the dust attenuation area defined by our SEDfitting results.

We note that some of the z-band (λ rest ∼ 2200Å) stacked data points have lower observed flux than its consecutive ones (i-band and Y-band) in our stacked LBGs SEDs. It could be a consequence of larger dust attenuation due to the coincidence with the bump at 2175Å that appear in some dust attenuation laws (e.g. MW). However, we do not have enough spectral resolution to constrain its amplitude and shape for our FUV-to-FIR stacked LBGs.

IRX-β UV plane and the slope of the dust attenuation curve

We previously reported a large dispersion in the IRX-β UV plane associated to a M * variations. We also obtain a diversity of dust attenuation curves along our LBG sample, which correlates with the M * . We check here, if the dispersion on the IRX-β UV plane is a consequence of a variety in M * or in the shape of the dust (Pei 1992) and two additional attenuation curve derived from; UV-selected galaxies at z > 1 (Buat et al. 2011b) and star-forming galaxies at z ∼ 2 to 6 (Scoville et al. 2015). attenuation curve. A few recent studies showed that the diversity in the shape of the dust attenuation curves has a strong impact to shape the IRX-β UV plane (e.g. Salmon et al. 2016;Lo Faro et al. 2017).

The observed UV part of the spectrum, where the β UV is calculated, depends on different physical properties of a galaxy (SFH, metallicity, dust attenuation law, ...). If we focus on the dust attenuation law and fix the others ones, when the UV shape of the dust attenuation curve changes the observed UV part of the spectrum will be consequently modified, and its associated β UV . To test it, we use CIGALE in its SED-modeling mode to compute synthetic SEDs with a fixed configuration (Delayed-SFH with τ=100Myr and age = 300Myr, metallicity = 0.02Z , and fixed IR model), and different dust attenuation laws (-0.4 < δ < 0.4) and reddening (0 < E(B-V) < 1.5).

Figure 11 illustrates the IRX-β UV plane color coded as a function of δ. The continuous lines represent the simulations performed with CIGALE for the different δ values. We can see that dust attenuation curves steeper than the Calzetti one (δ < 0) tend to locate bellow the M99 relation, and grayer ones (δ > 0) above. Note that the IRX-β UV relations simulated for each δ are not unique as we show here. This is a consequence of fixing the metallicites and SFH. These relations should be broader if we consider different SFHs and metallicities, as in Salmon et al. (2016). We over plot the results from the SED-fiting and stacking analysis as a function of β UV and M * , and the combination of both parameters (LBG-β, LBG-M, LBG-Mβ1, and LBG-Mβ2). Taken into account the uncertainties in the SFH and metallicity, our stacked LBGs are in agreement with our simulations. Our stacked LBGs with associated δ ∼ 0 follow M99, and when δ increases the objects present bluer β UV and higher IRX values than M99 relation. We conclude that one of the main drivers to shape the IRX-β UV plane is the dust attenuation curves, as previ- ous studies in other galaxy population (e.g. Salmon et al. 2016;Lo Faro et al. 2017).

Summary and conclusions

We investigate the full rest-frame FUV-to-FIR emission of LBGs at z ∼ 3 by stacking analysis at the Optical (BVriz bands), NIR (Y JHK s bands), IRAC (3.6, 4.5, 5.6 and 8.0 µm), MIPS (24µm), PACS (100 and 160 µm), SPIRE (250, 350, and 500 µm), and AzTEC (1.1mm) observations. We use a sub-sample of ∼17000 LBGs from the previously ones selected and characterized in AM16. We split our LBG sample as a function of the single parameters L FUV , β UV , and M * , and the combination of both β UV and M * in the (β UV , M * ) plane. This allows us to build thirty restframe FUV-to-FIR LBGs SEDs at z ∼ 3, and we investigates the evolution of their physical properties as a function of the binning parameters.

We use CIGALE, a physically-oriented spectral synthesis and SED-fitting code, to analyze our rest-frame FUV-to-FIR stacked LBGs SEDs. CIGALE provides us with: i) the synthetic model spectra that better fit our stacked LBGs SEDs by χ 2 mization, and ii) the mean physical parameters that characterize each of the stacked LBG SEDs by applying a Bayesian analysis. After performing exhaustive checks on the validity of the derived physical parameters, we conclude that the SED-fitting analysis can derive fully consistent physical parameters (M * , L IR , A FUV , SFR, and change of the slope of the dust attenuation law with respect to Calzetti -δ).

We use the stacked LBGs SEDs and their associated modelled spectra to build a library of thirty SEDs and templates of LBGs at z ∼ 3. Thanks to the binning configuration used to perform the stacking analysis, we derive a set of templates with a large variety of physical properties. The library contains templates within an interval of stellar mass, 9.2 < log(M * [M ]) < 11.4, SFR, 20 < S FR [M yr -1 ] < 300, β UV , -1.8 < β UV < -0.2, FUV dust attenuation, 1.5 < A FUV [mag] < 4.8, IR luminosity, 11.2 < log(L IR [L ]) < 12.7, and FUV luminosity, 10.4 < log(L FUV [L ]) < 11.2. This diversity makes our library quite versatile, in the sense that our templates fit the physical properties of a large population of star-forming galaxies at z ∼ 3.

We use the mean physical parameters derived from CIGALE, for each of our rest-frame FUV-to-FIR stacked LBGs SEDs, to investigate their dust attenuation properties. We conclude the following:

-Our LBG sample follows the well known IRX-β UV calibration from local starburst galaxies (M99), when LBGs are stacked as a function of β UV (LBG-β), being in agreement with recent works done on star-forming galaxies at high-z (Koprowski et al. 2018;McLure et al. 2018;Bourne et al. 2017). -Combining simulation and the stacked results from our LBG sample, we corroborate that the slope of the dust attenuation curve is one of the main drivers to shape the IRX-β UV plane.

Fig. 1 :

 1 Fig. 1: Upper maps are the final MIPS stacked images for the first bin as a function of the L FUV (10.2 < log(L FUV [L ]) > 10.5).The left panel shows the stacked image using the weighted stacking analysis and the right one the general stacking method. The bottom panel shows the radial profiles estimated from the upper images. The blue line corresponds to the stacking analysis performed using the general method, and the red line correspond to the stacking analysis performed using the weighted method. The lighter blue and red points correspond to the value of each pixel from the stacked images.

  Fig. 2: CIGALE mock analysis results. On the x-axis the true parameter values are provided by the best-fit model + noise consitite the mock SED catalogue. These true values are compared with the PDF one estimated by CIGALE, on the y-axis. The 1-1 correlation line is shown as a black long dashed line in each panel. From top-left to bottom-right the results for the age of the stellar population, e-folding timescale of delay-τ SFH, SFR, FUV dust attenuation, E(B-V), slope of the dust attenuation curve -δ, stellar mass, total IR luminosity, and U min . The regression lines for each assumed configuration are also plotted as red solid lines.

Fig. 3 :

 3 Fig. 3: Stacked LBGs SEDs and the best-fit models for the stacking as a function of L UV .

β

  Fig.5: IRX-β UV diagram. Lines show different IRX-β UV relations: the well-known local starburst galaxies calibration (M99, triple-dot-dashed line), the aperture correction of M99 relationTakeuchi et al. (2012, dashed line), and the one derived from a FIR stacking analysis over the same LBG sample (AM16, dotted blue line). In addition, we plot the IRX-β UV results from LBGs at redshifts, z ∼ 3(Koprowski et al. 2018, orange dots) and 2 < z < 3(Bouwens et al. 2016, gray dots), star-forming galaxies at redshifts, 2 < z < 3 (McLure et al. 2018, green dots) and 1.5 < z < 2.5(Reddy et al. 2018, black dots), and a mass selected sample at redshift z ∼ 3(Bourne et al. 2017, blue dots). The downward triangles represent upper limits.

Fig. 6 :

 6 Fig.6: IRX-M * diagram. Our rest-frame FUV-to-FIR stacking analysis as a function of M * (LBG-M) results and the power law fit are shown by red squares and solid line. The blue line presents the consensus relationship for star-forming galaxies at redshifts 2 < z < 3(Bouwens et al. 2016), and defined by a combination of different works done in the literature(Reddy et al. 2010;Whitaker et al. 2014;Álvarez-Márquez et al. 2016). This relation includes the one previously derived from a FIR stacking analysis over our LBG sample (AM16, dotted blue line). We overplot the relations obtained from a stacking analysis of UV-selected galaxies at redshift z ∼ 3(Heinis et al. 2014, dashed line), and a complete sample of star-forming galaxies up to z ∼ 4(Pannella et al. 2015, dot-dashed line). We also show the IRX-M * results from LBGs at redshifts, z ∼ 3(Koprowski et al. 2018, orange dot) and 2 < z < 3 (Bouwens2016, gray dot), and star-forming galaxies at redshifts, 2 < z < 3 (McLure et al. 2018, green dots) and 1.5 < z < 2.5(Reddy et al. 2018, black dots). The downward triangle represent upper limits.

Fig. 7 :

 7 Fig. 7: Results for LBG-Mβ1 (green) and LBG-Mβ2 (red) in the IRX-β UV and IRX-M * plane. Left panel: IRX-β UV diagram, for reference see Fig. 5. The LBG-Mβ2 results are shown in red, the tonalities represent the increase of the M * from 9.75 (light red) to 11.50 (dark red). The LBG-Mβ1 are showed in green, the tonalities represent the two different bins in M * , 9.75 < log(M * [M ]) < 10.65 (light green filled square) and 10.65 < log(M * [M ]) < 11.5 (dark green filled circles). The filled squares represent the bluer β UV bins and the filled circles the redder β UV bins. Right panel: IRX-M * diagram, for reference see Fig. 6. The LBG-Mβ2 results areshown in red, the tonalities represent the two different bins in β UV , -1.7 < β UV < -0.5 (light red filled square) and -0.5 < β UV < 0.5 (dark red filled circles). LBG-Mβ1 are shown in green, the tonalities represent the increase of the β UV from -1.7 (light green) to 0.5 (dark green). The filled squares represent the low M * bins and the filled circles the high M * bins.

Fig. 8 :

 8 Fig. 8: Empirical IRX calibration by combination the β UV and M * . The plane fit in the 3D IRX-β UV -M * diagram provides us with the equation log(IRX) = (0.51 ± 0.06)β UV + (0.37 ± 0.08) log(M * [M ]) -1.89 ± 0.40 that can be used to correct the dust attenuation when β UV and M * are known. Y-axis is the IRX obtained from SED-fitting analysis over the Stacked LBGs SEDs. X-axis corresponds to the IRX calculated from the plane fit equation. Black points are the results from the SED-fitting and stacking analysis as a function of β UV and M * , and the combination of the two (LBG-β, LBG-M, LBG-Mβ1, and LBG-Mβ2). The red line is the 1-to-1 relation.

  Fig. 9: Change of the slope of the dust attenuation law with respect to Calzetti (δ) as a function of stellar mass. The different colors shows each stacking analysis as a function of β UV and M * , and the combination of both parameters (LBG-β, LBG-M, LBG-Mβ1, and LBG-Mβ2). It suggests a correlation between the stellar mass of a galaxy and the shape of the dust attenuation curve, being steeper for LBGs with log(M * [M ]) < 10.25 and grayer LBGs with log(M * [M ]) > 10.25 than calzetti law.

Fig. 10 :

 10 Fig.10: Dust attenuation curves derived from the SED-fitting analysis over the stacked LBGs SEDs obtained from the stacking analysis as a function of M * . They are color coded as a function of M * , light blue represents larger stellar masses being grayer than Calzetti law (red line), and dark blue correspond to lower stellar masses with are steeper than Calzetti law. As a comparison, we also show the SMC dust extinction curve(Pei 1992) and two additional attenuation curve derived from; UV-selected galaxies at z > 1(Buat et al. 2011b) and star-forming galaxies at z ∼ 2 to 6(Scoville et al. 2015).

δ

  Fig. 11: IRX-β UV diagram seeing form a δ point of view. Continuous lines show the simulation performed by CIGALE with a fixed SFH and metallicity for a different the dust attenuation curves (see text). The square points represent the results from the SED-fiting analysis and stacking analysis as a function of β UV and M * , and the combination of the two (LBG-β, LBG-M, LBG-Mβ1, and LBG-Mβ2). Both are color-coded as a function δ. The dished line represents the classical M99 relation.

-

  Our LBG sample, stacked as a function of M * (LBG-M), is in agreement with most of the IRX-M * relations presented in the literature at high stellar masses (log(M * [M ]) > 10, Pannella et al. 2015; Heinis et al. 2014; Bouwens et al. 2016; Álvarez-Márquez et al. 2016). However, it provides higher IRX values for a given M * at low stellar masses (log(M * [M ]) < 10). -If our LBG sample is stacked as a function of a combination of both M * and β UV (LBG-Mβ1 and LBG-Mβ2) shows a large dispersion on the IRX-M * and IRX-β UV planes. We find that the evolution to higher M * or redder β UV makes differ them from the main IRX-M * and IRX-β UV relations. New empirical IRX relation combining both, β UV and M * , is suggested. -We investigate which dust attenuation curve is more likely to reproduce our LBG sample. We use Noll et al. (2009) dust attenuation law recipe, and derive a variation of δ (change in the slope of the dust attenuation law with respect to Calzetti's) from -0.3 to 0.3, which is well correlated with M * . Stepper slopes than Calzetti's are favoured in low stellar mass LBGs (log(M * [M ]) < 10.25), while grayer ones are favoured in high stellar mass LBGs (log(M * [M ]) > 10.25).

Table 1 :

 1 Sample binning

	ID	Range	log(L FUV [L ]) log(M * [M ])	β UV	z	N gal
		Stacking as a function of L FUV (LBG-L)
	LBG-L1	10.2 -10.5	10.35±0.09	9.67±0.47	-1.51±0.45 3.00±0.23 9145
	LBG-L2	10.5 -10.8	10.63±0.09	9.78±0.44	-1.56±0.41 3.08±0.25 5385
	LBG-L3	10.8 -11.1	10.91±0.08	9.94±0.37	-1.54±0.36 3.11±0.25 1510
	LBG-L4	11.1 -11.4	11.20±0.08	10.14±0.38	-1.54±0.32 3.15±0.24 199
		Stacking as a function of β UV (LBG-β)
	LBG-β1	-1.7 --1.1	10.49±0.24	9.75±0.44	-1.42±0.17 2.99±0.25 8531
	LBG-β2	-1.1 --0.7	10.41±0.22	9.88±0.62	-0.94±0.11 2.94±0.24 2445
	LBG-β3	-0.7 --0.3	10.30±0.20	10.00±0.66	-0.54±0.10 2.91±0.21 593
	LBG-β4	-0.3 -0.3	10.26±0.21	10.17±0.72	-0.10±0.21 2.89±0.21 114
		Stacking as a function of stellar mass (LBG-M)
	LBG-M1	9.75-10.00	10.51±0.24	9.86±0.08	-1.45±0.38 3.05±0.25 4517
	LBG-M2	10.00-10.25	10.55±0.26	10.10±0.08	-1.29±0.41 3.02±0.26 2257
	LBG-M3	10.25-10.50	10.58±0.29	10.35±0.08	-1.14±0.43 2.98±0.26 1041
	LBG-M4	10.50-10.75	10.55±0.27	10.60±0.08	-1.08±0.45 2.95±0.26 372
	LBG-M5	10.75-11.00	10.52±0.31	10.86±0.08	-0.95±0.43 2.94±0.25 160
	LBG-M6	11.00-11.50	10.47±0.28	11.16±0.12	-0.77±0.51 2.96±0.24	55
		Stacking as a function of (β UV , M * ) -LBG-Mβ1
	LBG-Mβ1-M 0 β 0	(-1.7 --1, 9.75-10.65)	10.55±0.25	10.00±0.20	-1.35±0.20 3.02±0.25 4799
	LBG-Mβ1-M 0 β 1	(-1 --0.5, 9.75-10.65)	10.41±0.22	10.12±0.23	-0.82±0.14 2.93±0.24 1343
	LBG-Mβ1-M 0 β 2	(-0.5 -0.3, 9.75-10.65)	10.26±0.19	10.18±0.23	-0.32±0.16 2.87±0.21 208
	LBG-Mβ1-M 1 β 0	(-1.7 --1, 10.65-11.50)	10.62±0.31	10.83±0.16	-1.28±0.19 2.97±0.26 146
	LBG-Mβ1-M 1 β 1	(-1 --0.5, 10.65-11.50)	10.44±0.20	10.84±0.16	-0.78±0.14 2.89±0.23 122
	LBG-Mβ1-M 1 β 2 (-0.5 -0.3, 10.65-11.50)	10.35±0.20	10.95±0.20	-0.41±0.05 2.90±0.22	40
		Stacking as a function of (β UV , M * ) -LBG-Mβ2
	LBG-Mβ2-M 0 β 0 (-1.7 --0.8, 9.75-10.125)	10.51±0.23	9.91±0.11	-1.31±0.23 3.01±0.25 4056
	LBG-Mβ2-M 0 β 1 (-0.8 -0.3, 9.75-10.125)	10.27±0.16	9.95±0.11	-0.60±0.19 2.94±0.22 338
	LBG-Mβ2-M 1 β 0 (-1.7 --0.8, 10.125-10.5)	10.61±0.28	10.27 ±0.10	-1.22±0.23 2.99±0.26 1373
	LBG-Mβ2-M 1 β 1 (-0.8 -0.3, 10.125-10.5)	10.37±0.21	10.29±0.10	-0.56±0.21 2.88±0.21 337
	LBG-Mβ2-M 2 β 0	(-1.7 --0.8, 10.5-10.75)	10.59±0.27	10.61±0.07	-1.19±0.25 2.96±0.26 260
	LBG-Mβ2-M 2 β 1	(-0.8 -0.3, 10.5-10.75)	10.40±0.22	10.60±0.08	-0.52±0.23 2.87±0.25	86
	LBG-Mβ2-M 3 β 0	(-1.7 --0.8, 10.75-11)	10.59±0.31	10.85±0.07	-1.18±0.23 2.94±0.25 102
	LBG-Mβ2-M 3 β 1	(-0.8 -0.3, 10.75-11)	10.35±0.18	10.86±0.08	-0.49±0.27 2.91±0.23	53
	LBG-Mβ2-M 4 β 0	(-1.7 --0.8, 11-11.5)	10.58±0.29	11.14±0.13	-1.15±0.24 3.04±0.26	27
	LBG-Mβ2-M 4 β 1	(-0.8 -0.3, 11-11.5)	10.40±0.24	11.18±0.12	-0.38±0.26 2.88±0.18	26

min et al. 2010)

3

. Given a specific catalogue and field image, the library extracts square cutout images centered on each source, and store them in a cube. Then, it generates a stacked image by averaging them. Additionally, it rotates each image by +π/2 with respect to the previous one, to cancel out the large-scale background gradients. To get valid and reliable results, the general method has to be corrected for different effects that generate a non homogeneous background on the stacked image (Bavouzet

Table 2 :

 2 Initial parameters used in the SED-fitting analysis

	Parameter	Symbol Range
	SFH (Delay-τ)
	age e-folding timescales	age τ	10 to 2000 Myr 20 to 1500 Myr
	SSP (Bruzual & Charlot 2003)
	Initial mass function	IMF	Chabrier (2001)
	Metallicity	Z	Fixed to 0.2Z
	Dust attenuation (Noll et al. 2009)
	Color excess Slope of the power law	E(B-V) 0.025 to 1.5 δ -0.5 to 0.5
	UV bump at 2175Å		Not included
	Dust emission (Draine & Li 2007)
	Mass fraction of PAH Minimum radiation field Power slope dU/dM ∝ U -α Dust fraction in PDRs	p PAH U min α γ	fixed to 3.9 0.7 to 50 fixed to 2.0 0.01 to 0.04

www.eso.org/sci/observing/phase3/data_release/ uvista_dr2.pdf

Herschel Database in Marseille: http://hedam.lam.fr/HerMES/

Article number, page 5 of 19 A&A proofs: manuscript no. output

Article number, page 8 of 19 Álvarez-Márquez, Burgarella, Buat et al.: FUV-to-FIR view of LBGs at z ∼ 3

Article number, page 16 of 19 Álvarez-Márquez, Burgarella, Buat et al.: FUV-to-FIR view of LBGs at z ∼ 3

Acknowledgements. JAM acknowledges support from the Spanish Ministry of Science, Innovation and Universities through project number ESP2017-83197. 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/.).

Appendix A: Comparison of methodologies to derive the β UV on our LBG sample

We find a large discrepancy between the β UV calculations from Section 3.2, and the ones presented by AM16 on the same LBG sample. Both methodologies use the same wavelength interval to perform the power-law fit. However, in Section 3.2, we use the best-fit synthetic spectral model obtained on the restframe UV-to-Optical SED-fitting analysis (β UV-SED ), and the rest-frame UV photometry in AM16 (β UV-Power ). Figure A.1 shows, in red color, the comparison between β UV-Power (y-axis) and β UV-SED (x-axis) for each LBG of the sample. It illustrates that the β UV-Power are redder than the β UV-SED in the range of β UV used in that work (1.9<β UV <0.3), and that their differences increase at redder values. Finkelstein et al. (2012) use a sample of z = 4 to 8 galaxies selected in CANDELS fields to compare the β UV-SED and β UV-Power . They found that the results from the SED-fitting method present essentially no bias at all redshifts, while the ones from power-law method are biased toward redder values. Reddy et al. (2015) used a sample of z ∼ 2 star-forming galaxies to compare the two different β UV calculations, concluding that the two measurements are highly correlated. However, Table 2 of their paper shows that the mean values of β UV-Power are 0.35 redder than the mean values of β UV-SED for the same population of galaxies. If the SED-fitting analysis provides a good model fit to the photometric data points, why would the results from the two methods be different? we perform a simulation to try to understand the origin of that differences. We perform the SEDfitting analysis over the full LBG sample to derive their bestfit spectrum and SED (photometric data points of the best-fit model). On the one hand, we use the best-fit spectrum to cal-culate the β UV-SED associated to the SED-fitting method. On the other hand, we use the best-fit SED to obtain the β UV-Power associated to the power-law method. In that case, we are working with a simulated spectrum and its associated SED, then, if both methods are consistent we should obtain similar β UV 's. Figure A.1 shows the β UV comparison, in green, when the both β UV 's have been calculated from the same best-fit model, using the associated spectrum and SED. We can see that both methods are consistent as the differences, even if systematic, are within the uncertainty showed in Section 3.2 (∆β UV ∼ 0.2).

We perform an additional test, where we include uncertainty at each photometric data point of the SED. We add a Gaussian scatter to the derived best-fit SED and re-compute the β UV-Power . The scatter has been defined equal to the observational photometric error for each band and LBG from our sample. 2018) also found similar discrepancies, and suggested that these inconsistencies are driven by scatter in measured values of β UV from limited photometry that serves to artificially flatten IRX-β UV relation. They obtained that the scatter is significantly reduced by determining β UV from SED-fitting analysis. .291 LBG-Mβ1-M 0 β 0 0.213±0.012 0.408±0.022 0.480±0.026 0.534±0.030 0.524±0.029 0.576±0.033 0.675±0.038 1.077±0.061 1.399±0.078 LBG-Mβ1-M 0 β 1 0.140±0.008 0.303±0.017 0.366±0.021 0.409±0.025 0.465±0.027 0.597±0.038 0.754±0.047 1.248±0.076 1.609±0.097 LBG-Mβ1-M 0 β 2 0.096±0.009 0.223±0.018 0.278±0.022 0.368±0.033 0.384±0.029 0.477±0.052 0.654±0.068 1.137±0.109 1.464±0.143 LBG-Mβ1-M 1 β 0 0.331±0.032 0.598±0.050 0.729±0.066 0.860±0.077 0.786±0.069 0.812±0.092 1.237±0.115 2.654±0.244 4.012±0.311 LBG-Mβ1-M 1 β 1 0.187±0.022 0.391±0.037 0.463±0.045 0.548±0.050 0.572±0.047 0.836±0.089 1.118±0.107 2.594±0.223 3.761±0.293 LBG-Mβ1-M 1 β 2 0.098±0.014 0.271±0.039 0.364±0.047 0.446±0.054 0.512±0.057 0.764±0.119 1.188±0.154 2.708±0.302 4.560±0.495 LBG-Mβ2-M 0 β 0 0.187±0.010 0.361±0.020 0.422±0.023 0.465±0.026 0.465±0.026 0.508±0.029 0.586±0.034 0.900±0.051 1.143±0.064 LBG-Mβ2-M 0 β 1 0.080±0.006 0.205±0.014 0.237±0.015 0.317±0.020 0.321±0.022 0.450±0.044 0.535±0.055 0.753±0.069 0.997±0.085 LBG-Mβ2-M 1 β 0 0.254±0.015 0.490±0.028 0.587±0.034 0.653±0.040 0.665±0.039 0.774±0.047 0.957±0.058 1.620±0.097 2.169±0.127 LBG-Mβ2-M 1 β 1 0.126±0.010 0.279±0.019 0.344±0.024 0.418±0.028 0.472±0.031 0.620±0.055 0.782±0.060 1.391±0.101 1.743±0.127 LBG-Mβ2-M 2 β 0 0.275±0.024 0.520±0.038 0.623±0.045 0.712±0.052 0.691±0.051 0.835±0.071 1.108±0.085 2.175±0.156 2.908±0.197 LBG-Mβ2-M 2 β 1 0.133±0.015 0.299±0.030 0.392±0.037 0.504±0.057 0.536±0.044 0.780±0.090 1.104±0.101 2.227±0.172 3.205±0.252 LBG-Mβ2-M 3 β 0 0.325±0.036 0.607±0.055 0.708±0.066 0.844±0.077 0.800±0.074 0.896±0. 112 Appendix D: Physical parameters derived from the Bayesian analysis in the SED-fitting analysis

Appendix B: Optical to NIR stacking results

Stacking results at