Rest-frame FUV-to-FIR view of Lyman break galaxies at $z\sim3$: Templates and dust attenuation
J. Álvarez-Márquez D. Burgarella V. Buat O. Ilbert P. G. Pérez-González
astro-ph.GA
 


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

   






^1
Centro de Astrobiología (CSIC-INTA), Carretera de Ajalvir, 28850 Torrejón de Ardoz, Madrid, Spain

javier.alvarez@cab.inta-csic.es

^2
Aix-Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR7326, 13388, France
Descargar



   



 
  

  
This work explores, from a statistical point of view, the rest-frame far-ultraviolet (FUV) to far-infrared (FIR) emission of a population of Lyman-break galaxies (LBGs) at z∼3 that cannot be individually detected from current FIR observations.
  
We performed a stacking analysis over a sample of ∼17000 LBGs at redshift 2.5<z<3.5 in the COSMOS field. The sample is binned as a function of UV luminosity (L_FUV), UV continuum slope (β_UV), and stellar mass (M_*), and then stacked at optical (BVriz bands), near-infrared (YJHKs 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 obtained 30 rest-frame FUV-to-FIR spectral energy distributions (SEDs) of LBGs at z∼3, and analyzed these with the CIGALE SED-fitting analysis code. We were able to derive fully consistent physical parameters, that is, M_*, β_UV, L_FUV, L_IR, A_FUV, star formation rate, and the slope of the dust attenuation law; we built a semiempirical library of 30 rest-frame FUV-to-FIR stacked LBG SEDs as functions of L_FUV, β_UV, and M_*.
  
We used the so-called IR-excess (IRX ≡ L_IR / L_FUV) to investigate the dust attenuation as a function of β_UV and M_*. Our LBGs, averaged as a function of β_UV, follow the well-known IRX-β_UV calibration of local starburst galaxies. Stacks as a function of M_* follow the IRX-M_* relationship presented in the literature at high M_* (log(M_* [M_⊙])> 10). However, a large dispersion is shown in the IRX-β_UV and IRX-M_* planes, in which the β_UV and M_* are combined to average the sample. Additionally, the SED-fitting analysis results provide a diversity of dust attenuation curve along the LBG sample, and their slopes are well correlated with M_*. Steeper dust attenuation curves than Calzetti's are favored in low stellar mass LBGs (log(M_* [M_⊙]) < 10.25), while grayer dust attenuation curves are favored in high stellar mass LBGs (log(M_* [M_⊙]) > 10.25). We also demonstrate that the slope of the dust attenuation curves is one of the main drivers that shapes the IRX-β_UV plane.



   

   Rest-frame far-ultraviolet to far-infrared view of Lyman break galaxies at z∼3: Templates and Dust attenuation
    J. Álvarez-Márquez1,2
D. Burgarella2
V. Buat2
O. Ilbert2
P. G. Pérez-González1

    Received ; accepted 
==============================================================================================================





§ INTRODUCTION


About 20 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 <cit.>. 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, ). These galaxies have been a key population to investigate the mass assembly of galaxies during the first gigayears of the universe <cit.>.

The rest-frame ultraviolet (UV) light, mostly emitted by young and massive stars, has been commonly used as a star formation rate (SFR) tracer. However, the interstellar dust scatters or absorbs the UV light, and 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 general picture of star formation activity across cosmic time peaks around z∼2-4 and significantly drops at z>4 (see  for a review). The obscured star formation dominates the total SFR density (SFRD) over the redshift range 0<z<3, and corresponds to half of the SFRD at z=3.6 <cit.>. Then, it is necessary to combine the UV and IR emission to determine the complete energy budget of star formation. However, because of 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 dust attenuation (e.g. ). 

Only a few LBGs have been directly detected at z∼3 in the mid-infrared (MIR) and FIR <cit.> thanks to Spitzer and Herschel. This sample is biased against submillimeter bright galaxies and is not representative of the LBG population in terms of stellar mass (M_*), dust attenuation, and SFR <cit.>. The recently gain in sensitivity with the Atacama Large Millimeter Array (ALMA) at submillimeter and millimeter wavelengths provides some insight into the obscured SFR at high-z <cit.>. Even so, it is still very difficult to obtain large samples to carry out statistical and representative analyses. Therefore, stacking analysis techniques <cit.> have been applied to relatively large samples of LBGs to derive their FIR/submillimeter emission <cit.>.

Owing to the lack of FIR/submillimeter 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, ) and the so-called IR excess (IRX ≡ L_IR / L_FUV) calibrated on local starburst galaxies by  (M99 hereafter; 1999). Despite its general use to infer the dust attenuation at high-z, there are several complications that make this method uncertain. 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 analyses of LBG and star-forming galaxies at redshifts of 1.5 < z < 5 have shown different behaviors; some of these results follow the M99 relation <cit.> and others lie above <cit.> or below this relation <cit.>. These deviations have been shown to be driven by the stellar masses <cit.>, the shape of the dust attenuation curve <cit.>, and the sample selection <cit.>. Direct detections of individual LBGs and IR-selected sources have suggested a large scatter in the IRX-β_UV plane (e.g., ). In addition, the stellar mass has been found to correlate with the dust attenuation, and its relation does not seem to evolve with redshift <cit.> at least in the intermediate range 1.0 < z < 4.0.

With the objective to investigate deeply 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 <cit.> and presented in a first publication (, AM16 hereafter). A statistically controlled stacking analysis from the FIR to millimeter 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 subsamples as a function of several parameters: L_FUV, β_UV, and M_*. We investigated the dust properties for each subsample 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. This approach allows 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 us 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 follows. In Section <ref>, we describe the COSMOS field data and photometry used in our analysis. In Section <ref>, the definition of the LBG sample and its associated photometric redshift, M_*, L_FUV, and β_UV are presented. In Section <ref>, we define the method to stack the LBG sample from optical to FIR observations and obtain the stacked rest-frame FUV-to-FIR LBGs SEDs at z∼3. In Section <ref>, 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 semiempirical templates of LBGs at z∼3. In Section <ref>, 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 <ref> 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 <cit.> and <cit.> IMFs. To convert from <cit.> to <cit.> IMFs, we divide M_* _Salpeter by 1.74 <cit.>.





§ DATA


We used the available optical to FIR imaging from the COSMOS field, and the optical/NIR multicolor catalog (, version 2.0). The specific data sets are presented in the following Section.



 §.§ Optical and near-infrared data


We used optical broadband imaging (B_J, V_J, r^+, i^+, z^++) from the COSMOS-20 survey <cit.>, observed at the prime-focus camera (Suprime-Cam) on the 8.2 m Subaru Telescope. In combination with the YJHK_s broadband imaging released in the UltraVISTA DR2[<www.eso.org/sci/observing/phase3/data_release/uvista_dr2.pdf>] <cit.> and observed with the VISTA InfraRed CAMera (VIRCAM) instrument on the Visible and Infrared Survey Telescope for Astronomy (VISTA). 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  (2016, private communication). 

We additionally used optical/NIR broadband photometry (Subaru: B_J, V_J, g^+, r^+, i^+, z^+, and VISTA: Y,J, H and K_s) from the COSMOS multicolor catalog (, version 2.0). The photometry was performed using SExtractor in dual-image mode over an aperture of 3^'' centered at the position of the i^+ band detection. The COSMOS multicolor catalog was updated by <cit.>, in which they provide a new NIR selection catalog on the ultra VISTA-DR2 observations. However, we still used the i-band selection catalog (, version 2.0) because the LBG sample, used in this analysis (see Sect. <ref>), has been proposed 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 <cit.> 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 <cit.> mapped the COSMOS field at 100 and 160 with the Photodetector Array Camera and Spectrometer (PACS;  ), and 250, 350, and 500 μm with the Photometric Imaging Receiver (SPIRE; ) as part of the Evolutionary Probe survey (PEP; ) and Herschel Multi-Tiered Extragalactic Survey (HerMES, ). The PACS images (100 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 <cit.>. These maps were downloaded from HeDaM[Herschel Database in Marseille: <http://hedam.lam.fr/HerMES/>]





§ LBG SAMPLE


The LBG sample for this analysis is a subsample of that selected and characterized by AM16. The original LBG sample was selected in the COSMOS field by means of the classical U-dropout technique <cit.> using the broadband filters u^*, V_J, and i^+. This original sample 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 (YJHK_s bands) covers a reduced area of the COSMOS field of 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 UltraVISTA 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 used the photometric redshifts (photo-z) and stellar masses computed by  (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% ( 2009, version 2.0). The stellar masses were derived by applying SED-fitting techniques to the available optical and NIR photometry, and assuming; <cit.> single stellar population (SSP), exponentially declining SFH, and <cit.> initial mass function (IMF). <cit.> suggested 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 method presented by <cit.>. 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 limited our LBG SEDs to the available optical-to-NIR broadband 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-frame [Å] < 3500. This allows us to perform an SED-fitting analysis over the same rest-frame spectral range for each LBGs independent of its redshift. 

We used CIGALE (Code Investigating GALaxy Emission; , see Sect. <ref>) 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 <cit.> synthetic stellar population libraries and <cit.> IMF. We varied 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 <cit.> dust attenuation law), and SFH (SFH ∝ e^-t/τ, where τ = 0.1, 10, 100, 10^3, 10^4, 10^5, -300, -10^3, -10^4 Myr).

The value β_UV is calculated by power-law fit (f_λ∝λ^β_UV) to the best-fit synthetic stellar population model that has been derived by CIGALE. The power-law fit is performed at rest-frame wavelength range 1250 < λ_rest-frame [Å]  < 2000 following the spectral windows defined by <cit.>. 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 independent of the redshift of the galaxy. The value L_FUV is computed in a squared bandpass defined as the GALEX FUV band (λ_eff = 152.8 nm and effective bandwidth = 11.4 nm). We checked the uncertainties in the β_UV and L_FUV estimations using the mock analysis tool from CIGALE (see Sect. <ref> for details of the mock analysis tool). We obtained a dispersion in the β_UV of Δβ_UV∼0.2 and a very robust L_FUV measure with Δlog(L_FUV[L_⊙])∼0.03. Additionally, we obtained a systematic offset between the β_UV calculated by power-law fit directly to the photometry (AM16) and that derived from power-law fit to the best-fit model obtained by SED-fitting analysis of the photometry (see Appendix <ref> for a detailed discussion).  





§ STACKING ANALYSIS


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., ). Our LBG sample is complete at V_J and i^+ bands, where the sample has been selected by means of the U-dropout technique. However, the detection of individual LBGs at other optical/NIR bands are dropped up to 70% of completeness, and it gets worse at longer wavelengths as far as only few of the LBGs could be detected at FIR observations (SPIRE 250 μm). We then 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 us to detect statistically the LBG population at redder wavelength at the expense of averaging over their individuals properties.

Previous works that computed average rest-frame FUV-to-FIR LBG SEDs (e.g., ) performed a simple average to the optical/NIR photometry. Our LBG sample is 70-80% complete at NIR and 90-100% complete 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 LBGs 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  β_UV and M_* in the (β_UV, M_*) plane. Table <ref> 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 with 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 (S/N) 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 )[<http://www.ias.u-psud.fr/irgalaxies/downloads.php>]. Given a specific catalog and field image, the library extracts square cutout images centered on each source and stores these images in a cube. Then, it generates a stacked image by averaging them. Additionally, it rotates each image by +π/2 with respect to the previous image 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 nonhomogeneous background on the stacked image (Bavouzet 2008; ; ; AM16). The following three points are of major importance: 



  * Correction for clustering of the input catalog. The 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 catalog. The complete formalism is detailed in AM16.

  * Correction for incompleteness of the input catalog in the dense regions. A bias is produced in the stacked image when the population of sources is not complete (; Bavouzet 2008; ; ; ). 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 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 mostly stack the areas in which the background presents lower levels and lose the contribution of the areas with higher levels. 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.


  * 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 number of bright neighbors may generate an additional noise level that perturbs the final stacked image. We removed the contribution of field sources in our stacking procedure by performing a weighted stacking analysis to mask them. We proceeded as follows:



  * We generated a segmentation map, using SExtractor <cit.>, on the image where the stacking analysis is performed. A source extraction was applied to detect all the objects with flux larger than 3σ_Back, excluding the LBGs from our sample. 

  * We created a weight map 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, ∑_iI_i W_i / ∑_iW_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. 



  

Figure <ref> shows the improvement that we get in the final stacked image 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. This new method 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 the final stacked image, which mainly depends on the image quality (PSF, confusion noise, and depth), observed wavelength, and sample selection method. In Sections <ref>, <ref>, and <ref>, we discuss which of the previous corrections are taken into account during the stacking analysis for each of the data sets.



 §.§ Far-infrared 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 explained well in AM16. We used the calibrated PACS and SPIRE maps without cleaning the images from any detected source, and we did 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 tend to average homogeneously to a constant background level. We demonstrated this by performing a random stacking analysis at the same number of positions than objects in each of our sample bins.  

The 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 autocorrelation 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 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 catalog to deconvolve the emission of our LBG population and the clustering contribution. Additionally, we corrected our stacking analysis for the incompleteness of the input catalog in the dense regions following the analysis done in AM16.



 §.§ Mid-Infrared analysis


The following procedure is 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 MIR observations reduces the confusion noise and improves the quality of the MIR images. This generates a crowded image with a large number of field sources. The presence of the large amount of field sources in the maps generates an additional noise level and clustering contribution that perturbs 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 catalog 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 for this effect in the MIR wavelength.

We performed an aperture photometry to calculate the flux density of the stacked LBGs in the IRAC and MIPS observations. We made 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 applied 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)[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/>)]. We used bootstrap resampling to obtain the mean values and errors. We repeated the above procedure over 3000 random bootstrap samples and adopted 1σ of the distribution of the derived fluxes as the uncertainty for our results. Appendix <ref> presents the derived fluxes and uncertainties of our stacking analysis at FIR and MIR. This appendix includes the stacking analysis performed for AzTEC (1.1mm), and previously presented in AM16.



 §.§ Optical/near-infrared 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 a 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 catalog. 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 <ref> presents the derived fluxes and uncertainties of our stacking analysis in optical and NIR wavelengths. We verified the method 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


A Python-based modular code, CIGALE <cit.>[Code Investigating
GALaxy Emission (CIGALE): <https://cigale.lam.fr/>] is  dedicated to fitting SEDs from the UV to the FIR and to creating galaxy emission models in the same wavelength range. The CIGALE code is  modular and allows the user to create models by selecting independent modules for parameters such as the SFH, stellar models, shape of the dust attenuation curve, nebular emission, contribution of an active galactic nuclei, and dust emission templates or models.



 §.§ Initial parameters 


Deriving physical properties of galaxies by fitting synthetic SEDs to a multiband 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 (e.g., ). Table <ref> summarizes the initial parameters used to perform the SED-fitting analysis for our stacked LBGs SEDs at z∼3.



We adopted the stellar population models from <cit.> and an IMF from <cit.>. Delayed SFHs (SFR ∝ t/τ^2 e^-t/τ), with varying e-folding times are assumed to model the SFH of our LBGs at z∼3 <cit.>. We performed a test to check the influence of the chosen SFH on the final derived physical parameters, where we ran CIGALE with different SFHs recipes (exponential declining, exponential rising, constant, and delayed) using the same initial parameters from Table <ref>. The results suggest that our SED-fitting 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 independent of the SFH chosen with the exception of the age of the stellar population. 

The metallicity was fixed to 0.2 Z_⊙ <cit.>. We fixed the metallicity because i) the SED-fitting technique applied to broadband photometry can hardly constrain the metallicity due to the age-dust-metallicity degeneracy <cit.>, 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-metallicity relation (e.g., ). Then, we performed 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 those 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 <cit.> in Eq. <ref>. These authors used the Calzetti dust attenuation law (k(λ), ) as a baseline. This provides the possibility to i) include a UV bump at 2175 Å modeled by a Lorentzian-like 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)^δ
    ,
 

where A(λ) is the modified dust attenuation curve presented by <cit.>, 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's law, and λ_V = 5500Å. The shape of the dust attenuation curve depends on the stardust geometry, grain size distribution, etc. <cit.>. Complex dust geometries can produce a grayer or shallower dust attenuation curve than Calzetti's law <cit.>. Our SED fitting takes δ as a free parameter to mimic the variation on the shape of the dust attenuation curve. The absorption feature produced by dust at 2175 Å, UV bump, is set to zero because the stacked SEDs do not have enough spectral resolution to constrain it. 

The IR emission is modeled using <cit.> dust models. The validity of these models to reproduce the FIR emission of Herschel detected high-z main sequence galaxies was confirmed by <cit.>. The dust models require the fine-tuning of several free parameters: the mass fraction of PAH (polycyclic aromatic hydrocarbons, q_PAH), minimum and maximum radiation field (U_min, U_max), power slope dU/dM∝ U^-α (α), dust fraction in PDRs (γ), and dust mass (M_dust). However, <cit.> showed that fixed values of α = 2 and U_max = 10^6 can reproduce the SEDs of galaxies with a wide range of properties. The stacked SEDs only have the information of MIPS 24 μm (rest-frame 6μm) in the PAH part of the spectrum, which is not enough to constrain the q_PAH. After some initial tests with the q_PAH as a free parameter, we fixed q_PAH to 3.9%. This value best fit our SEDs by eye even though it is higher than that proposed by <cit.>, who found an interval from 1.12% to 3.19%. We consider U_min to vary between 0.7 to 50. <cit.> suggested using U_min > 0.7 because lower values correspond to dust temperatures below ∼ 15K, which we cannot constrain by our FIR photometry alone. The fraction of dust enclosed in the PDRs was fixed by <cit.> to γ=0.02, and <cit.> suggested that the γ value should not be higher than 0.04 (4%). We detected some variations at rest-frame wavelength 20-50 μm for our set of stacked SEDs, which suggest different γ values. Therefore, we used a γ range from 0.01 to 0.04. 

The adopted initial configuration (Table <ref>) gives a total of 609840 synthetic galaxy models/SEDs. This set of models is used to perform the SED fitting and Bayesian analysis over our 30 rest-frame FUV-to-FIR stacked LBGs SEDs at z∼3.
  


 §.§ Physical parameters for LBGs at z∼3
 

We ran CIGALE for the sample of 30 rest-frame FUV-to-FIR stacked LBGs SEDs at z∼3. We used the adopted initial configuration listed in Table <ref>. The CIGALE code estimates and derives the physical parameters, such as stellar mass, dust mass, age of stellar populations, SFR, rest-frame IR and FUV luminosities, FUV attenuation, and β_UV, as well as the input parameters from Table <ref>. Additionally, CIGALE provides the synthetic model that best fits the observational SEDs by χ^2 minimization.

The CIGALE code uses a Bayesian analysis to derive the physical parameters. This code builds the probability distribution function (PDF) for each parameter, and takes its mean and standard deviation <cit.>. We used a mock analysis to check the robustness, accuracy, and parameter degeneracy of our estimates. The mock analysis consists in generating a mock catalog of artificial SEDs, similar to the input LBGs, for which the physical parameters have been previously derived, and adding a Gaussian noise where the dispersion is taken from the observed uncertainty. Then, we compared these values with the CIGALE estimates (more details: ). This procedure allows us to check our ability in estimating and constraining the individual output parameters given the information provided by the detailed shape of the broadband SED. The results are summarized in Figure <ref>. We call the parameter value provided by the best-fit model “true” and used this value to compute the mock SEDs. This “true” parameter value, on the x-axis, is compared with the PDF value, on the y-axis, which is computed by SED-fitting analysis over the mock SEDs catalog. From top left to bottom right, Figure <ref> 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 independent 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 <ref>, 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 <ref> lists the input and additional physical parameters derived by SED-fitting analysis in the 30 stacked LBGs SEDs. The β_UV and L_UV are computed in a separate FUV-optical SED-fitting analysis; the input parameters and procedure are already explained in Section <ref>. 



 §.§ Templates of LBGs at z∼3


We used the results from the panchromatic stacking and SED-fitting analysis to build a library of empirical rest-frame FUV-to-FIR templates of LBGs at z∼3.[Templates  and  SEDs  are  only  available  at  the  CDS via  anonymous  ftp  to <cdsarc.u-strasbg.fr> (130.79.128.5)  or  via  <http://cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/630/A153>] The library contains 30 empirical SEDs and their best-fit synthetic model spectrum. The empirical SEDs correspond to the rest-frame FUV-to-FIR stacked LBGs SEDs derived by stacking analysis in Section <ref>. These SEDs result from combining ∼17000 individual LBGs in different bins as a function of L_FUV, M_*, β_UV, and a combination of 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 S/N, which is not comparable with any individual detected LBG at z∼3. The set of synthetic model spectrum associated with each stacked LBG SED is obtained with CIGALE.

Thanks to the binning configuration used to split our large LBG sample, we 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 < SFR  [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 versatile in the sense that our templates fit a wide range of physical properties for star-forming galaxies at z∼3.

Figure <ref> shows the four SEDs and templates derived in the stacking analysis as a function of L_FUV. If a sample binning is carried out as a function of L_FUV, the subsample becomes 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 tend to average to a similar mean SED with a different normalization factor, i.e., L_FUV.  Figure <ref> shows 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 shifted 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 <cit.>. 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.

     

      
  
The top left panel of Figure <ref> illustrates the six SEDs and templates derived from the stacking analysis as a function of M_*. Their rest-frame FUV-to-NIR SEDs present an increase of the rest-frame NIR flux, which is well known to correlate with M_* <cit.>, and a reddening of the UV/optical part of the spectrum. No significant variations are found in the shape of the FIR emission, but the IR luminosity increases proportionally with M_*. The top right panel of Figure <ref> shows the four SEDs and templates derived from the stacking analysis 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, which give an IR peak strongly shifted to longer wavelengths. Blue galaxies are very faint on FIR, and we suggest that our FIR stacking analysis on LBG-β1 is very 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 <ref> illustrates the six SEDs and templates derived from the LBG-Mβ1 stacking analysis. As a reminder, we split the sample into two large bins of M_*, and split each of these bins into three bins of β_UV. Their rest-frame FUV-to-NIR SEDs seem to subdivide 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 are found, except for an increase of the IR luminosity with M_* and β_UV. The bottom right panel of Figure <ref> 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 into β_UV, red and blue bins, gives access to two templates for each M_*, with similar NIR shapes, but in a 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 means of the IRX. At that time, we had access to the full FIR/submillimeter part of the spectrum to compute the average L_IR. However, the β_UV and M_* values associated with 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 us to derive the  β_UV and M_*  directly from the stacked LBG SEDs. We then reviewed the information in the IRX-β_UV and IRX-M_* planes using the results obtained from the rest-frame FUV-to-FIR SED-fitting analysis on the stacked LBGs SEDs.



 §.§ Dust attenuation as a function of β_UV


Figure <ref> 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-known IRX-β_UV relation calibrated from local starburst galaxies (M99) and the relation recomputed by <cit.> 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 <cit.> and 2<z<3 <cit.>, star-forming galaxies at redshifts, 2<z<3 <cit.> and 1.5<z<2.5 <cit.>, and a mass selected sample at redshift z∼3 <cit.>. 

    

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  an IRX value that is ∼2-3 times higher for a given β_UV than the FIR values. The difference resides on the method used to derive the β_UV in each of the analyses. In Appendix <ref>, we demonstrate that β_UV-power (AM16) are biased to redder values in comparison with β_UV-SED (this work, see Section <ref>). If the IRX-β_UV results from the FIR analysis are corrected by the relation between the β_UV-power and β_UV-SED, obtained in Appendix <ref>, they are in agreement with the IRX-β_UV relation obtained from the FUV-to-FIR analysis. <cit.> have also confirmed that discrepancies in the IRX-β_UV relation are due to bias in the method used to determine the β_UV. These authors 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 <ref>. 

Our FUV-to-FIR stacking and SED-fitting analysis results as a function of β_UV follow the IRX-β_UV local relation (M99) within the uncertainties. These results are in agreement with previous stacking analyses done on star-forming galaxies at high-z. <cit.> found that the dust-corrected UV-SFR derived from the M99 relation presents a good match with the FIR and radio SFR estimators for stacked IRAC spectropically confirmed LBGs at z∼3. <cit.>, by stacking a sample of ∼4000 LBGs at redshifts, 3<z<5, have concluded that LBGs are consistent with the M99 IRX-β_UV relation and do not present a redshift evolution. Similar results were obtained by <cit.>, who studied a sample of star-forming galaxies at redshifts, 2<z<3, from a deep ALMA 1.3 mm continuum data. Additional stacking analysis of spectroscopically confirmed z∼2 LBGs in FIR <cit.>, z∼4 LBGs in the radio continuum (1.4GHz; ), and direct detections of LBGs in PACS at lower redshift (z∼1; ) also agree with the M99 relation. <cit.> have 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). These authors have obtained similar results as for local starburst galaxies in M99, although their results deviate to higher IRX values at high stellar masses, which is similar to the analysis presented in Sect. <ref>. However, <cit.> have used ∼3500 star-forming galaxies at redshifts, 1.5<z<2.5, and have concluded that the M99 IRX-β_UV local relation overpredicts the IRX by a factor of ∼3 at a given β_UV. <cit.> and <cit.> used ALMA observations of LBGs at z∼2-3, and also showed that the M99 relation overpredicts the IRX. The ALMA observations contain small samples that could make it difficult to constrain the IRX-β_UV relation statistically. In particular, <cit.> stacked a LBG sample at z=2-3 in bins of M_* and β_UV. Their results suggested an IRX-β_UV values below the M99 relation, but the high M_* bin (log(M_* [M_⊙]) > 9.75) is in agreement within the uncertainties.



 §.§ Dust attenuation as a function of M_*


Figure <ref> 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 <cit.> for star-forming galaxies at redshifts 2<z<3 and defined from the combination of different works done in the literature <cit.>. This relation includes that 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 <cit.> and a complete sample of star-forming galaxies up to z∼4 <cit.>. We also show the IRX-M_* results from LBGs at redshifts, z∼3 <cit.> and 2<z<3 <cit.>, and star-forming galaxies at redshifts, 2<z<3 <cit.> and 1.5<z<2.5 <cit.>. 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 <cit.>, the incompleteness of our LBG sample at z∼3 in terms of M_* is 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 <ref>). <cit.> showed that the i-band photometric redshift catalog (, version 2.0), used to build our LBG sample in AM16 is 90% complete 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 are some indications in the literature that the origin of the departure might not be due to M_* bias. For instance, <cit.> showed 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, however. It could also be noticed in <cit.> 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, <cit.> 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 <ref> 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 <ref>, 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 <cit.>. If this population is removed to perform the stacking analysis and derive the mean IRX-M_* relation (light red filled square in Figure <ref>), the results are closer to the IRX-M_* relations from <cit.> and <cit.>, and the star-forming galaxies at 2<z<3 stacked in ALMA observation <cit.>.

   
   


 §.§ 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 is split on the (M_*, β_UV) plane in two different ways. For LBG-Mβ1, the LBG sample is divided into two bins of M_*, and each of these is split as a function of β_UV. And for LBG-Mβ2, it is divided into two bins of β_UV and each of these is split as a function of M_*. Figure <ref> 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 <cit.> in the IRX-β_UV plane. These authors 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 <cit.>, 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. In particular, 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 <ref>.

The dependence of M_* in the IRX-β_UV plane and β_UV in the IRX-M_* plane makes us think about the utility and efficiency of providing a new empirical IRX relation combining the β_UV and M_*. It may reduce the uncertainty to derive the dust attenuation for galaxies for which only UV/optical/NIR observations are 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 <ref> 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 star-forming galaxies at high-z (e.g., ). 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 Å <cit.>. However, the feature is weaker in the Large Magellanic Cloud (LMC; ) and absent in the Small Magellanic Cloud (SMC; ). This diversity is also reflected by systematic changes in the FUV slope, which gets steeper from MW to SMC. <cit.> deduced a mean dust attenuation curve from a sample of 39 local UV-bright starburst galaxies. Their measurements are characterized by a grayer FUV slope than both MW and LMC extinction curves and an absence of the 2175 Å absorption feature. At higher redshifts, a variety can be also seen. <cit.> found for a UV-selected galaxy at z>1 a clear present of 2175 Å bump and evidence for a steeper rise of the dust attenuation curve in comparison to the Calzetti law. Other studies support this model, demonstrating that the 2175 Å bump is typically visible in normal star-forming galaxies (e.g., ). Other authors have also suggested that dust attenuation curves vary significantly at any redshift (e.g., ). 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 within a galaxy, the line-of-sight galaxy orientation, and the galaxy type <cit.>.


The results from the SED-fitting analysis over our stacked LBGs SEDs suggest 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 <ref>), which produce dust attenuation curves steeper and grayer than Calzetti law (δ=0), respectively. Similar results were reported by <cit.>, who found a range of δ's from -0.5 to 0.2 in star-forming galaxies at z∼1.5-3. Figure <ref> 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's law, and the large stellar mass population of LBGs (log(M_* [M_⊙]) > 10.25) favors grayer Calzetti's law. A similar trend with the stellar mass is seen in the local Universe by <cit.>; these authors have suggested that the apparent dependence with the stellar mass is due to a relation between the stellar mass and optical opacity. They conclude that opaque galaxies have shallower dust attenuation laws. <cit.> studied a sample of star-forming 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 at lower masses. However, <cit.> performed an analysis of star-forming galaxies selected by their rest-frame optical 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 Calzetti's, which differs from our findings.

  

Figure <ref> illustrates the derived dust attenuation curves in the stacked LBGs SEDs as a function of M_* (LBG-M). This figure shows how dust attenuation curves get steeper when the stellar mass decreases, up to the lowest stellar mass bin that presents a slope closer to the SMC dust attenuation curve <cit.>.  Other mean dust attenuation laws derived from UV-selected galaxies at z > 1 <cit.> and star-forming galaxies at z∼2 to 6 in the COSMOS field <cit.> are also compared; they are in agreement with the derived dust attenuation curves for our stacked LBGs SEDs.
  
We note in our stacked LBGs SEDs that some of the z-band (λ_rest∼ 2200Å) data points have lower observed flux than its consecutive data points (i-band and Y-band). This might be a consequence of a larger dust attenuation owing to the coincidence with the UV bump at 2175Å. As we commented in Section <ref>, the spectral resolution of the FUV-to-FIR stacked LBGs SEDs is not sufficient to constrain the amplitude and shape of the UV bump, and its amplitude has been set to zero in our SED-fitting analysis. However, the presence or not of a UV bump in the dust attenuation law might modify the derived δ values. We performe an additional SED-fitting analysis with the amplitude of the UV bump as a free parameter to check the validity of the derived δ values. We obtain a mean differences, ⟨δ_D_λ_o, γ, E_bump-free - δ_D_λ_o, γ, E_bump=0⟩ = -0.013 ± 0.008, which are within the derived δ's uncertainties of the initial SED-fitting analysis.



 §.§ IRX-β_UV plane and the slope of the dust attenuation curve


We previously reported a large dispersion in the IRX-β_UV plane associated with a M_* variations. We also obtain a diversity of dust attenuation curves along our LBG sample, which correlates with the M_*. We check in this section if the dispersion on the IRX-β_UV plane is a consequence of a variety in M_* or in the shape of the dust attenuation curve. A few recent studies showed that the diversity in the shape of the dust attenuation curves has a strong impact on shaping the IRX-β_UV plane (e.g., ; ). 

   
  
The observed UV part of the spectrum, where the β_UV is calculated, depends on different physical properties of a galaxy (e.g., SFH, metallicity,  and dust attenuation law).  If we focus on the dust attenuation law and fix the others properties, when the UV shape of the dust attenuation curve changes the observed UV part of the spectrum would  consequently be modified along with its associated β_UV. To test this, 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 <ref> illustrates the IRX-β_UV plane color coded as a function of δ. The continuous lines represent the simulations performed with CIGALE for the different δ values. Dust attenuation curves steeper than the Calzetti’s law, with δ < 0, are located below the M99 relation, and dust attenuation curves grayer than Calzetti’s law, with δ > 0, are located above. We note that the IRX-β_UV relations simulated for each δ are not unique as we show in this work. This is a consequence of fixing the metallicites and SFH. These relations should be broader if we consider different SFHs and metallicities, as in <cit.> and <cit.>. We overplot the results from the SED-fitting and stacking analysis as functions of β_UV and M_* and the combination of both parameters (LBG-β, LBG-M, LBG-Mβ1, and LBG-Mβ2). Taking 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 dust attenuation curves are one of the main drivers that shapes the IRX-β_UV plane, following previous studies in other galaxy population (e.g., ; ; ).  



§ 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 (YJHKs 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 subsample of ∼17000 LBGs from those previously  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 30 rest-frame FUV-to-FIR LBGs SEDs at z∼3 and investigate 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. The CIGALE code provides us with  the synthetic model spectra that better fit our stacked LBGs SEDs by χ^2 minimization and 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, i.e., 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 modeled spectra to build a library of 30 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 < SFR  [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 very 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 so-called IRX and the mean physical parameters derived for each of the rest-frame FUV-to-FIR stacked SEDs to investigate the dust attenuation properties of LBGs at z∼3. We conclude as follow:

 

  * 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-β), which is in agreement with recent works done on star-forming galaxies at high-z <cit.>.


  * 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, ). 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), the sample shows a large dispersion on the IRX-M_* and IRX-β_UV planes. We find that the evolution to higher M_* or redder β_UV makes them differ from the main IRX-M_* and IRX-β_UV relations. New empirical IRX relation combining  β_UV and M_* is suggested.


Finally, we investigate which dust attenuation curve is more likely to reproduce our LBG sample. We use <cit.> dust attenuation law recipe and derive a variation of δ (change in the slope of the dust attenuation law with respect to Calzetti's law) from -0.3 to 0.3, which is well correlated with M_*. Steeper dust attenuation curves than Calzetti's are favored in low stellar mass LBGs (log(M_* [M_⊙]) < 10.25), while grayer dust attenuation curves than Calzetti's are favored in high stellar mass LBGs (log(M_* [M_⊙]) > 10.25). Combining the simulation and 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.



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/>.).




aa


 



§ COMPARISON OF METHODS TO DERIVE THE Β_UV ON OUR LBG SAMPLE


We find a large discrepancy between the β_UV calculations from Section <ref> and those presented by AM16 on the same LBG sample. Both methodologies use the same wavelength interval to perform the power-law fit. However, in Section <ref>, we use the best-fit synthetic spectral model obtained on the rest-frame UV-to-optical SED-fitting analysis (β_UV - SED), and the rest-frame UV photometry  in AM16 (β_UV - Power). Figure <ref> shows, in red, the comparison between β_UV - Power (y-axis) and β_UV - SED (x-axis) for each LBG of the sample. This figure 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. 

<cit.> used 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 presents essentially no bias at all redshifts, while the results from power-law method are biased toward redder values. <cit.> 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 SED-fitting analysis over the full LBG sample to derive their best-fit spectrum and SED (photometric data points of the best-fit model). On the one hand, we use the best-fit spectrum to calculate the β_UV - SED associated with the SED-fitting method. On the other hand, we use the best-fit SED to obtain the β_UV - Power associated with 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. Figure <ref> shows the β_UV comparison, in green, when the both β_UV 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 <ref> (Δβ_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 recompute the β_UV - Power. The scatter has been defined as equal to the observational photometric error for each band and LBG from our sample. Figure <ref> illustrates the comparison between the β_UV - SED and the β_UV - Power derived from the best-fit model SED + scatter in blue. That shows a similar trend than that obtained when we compared the both β_UV in our LBG sample. This suggests that the differences between the two methods are basically due to the scatter in the photometric data points as a consequence of the photometric uncertainty.  <cit.> also found similar discrepancies and suggested that these inconsistencies are driven by scatter in measured values of β_UV from limited photometry that serves to flatten IRX-β_UV relation artificially. These authors obtained that the scatter is significantly reduced by determining β_UV from SED-fitting analysis. 




§ FLUXES AND UNCERTAINTIES AT OPTICAL AND NIR WAVELENGTHS FOR THE 30 STACKED LBGS SEDS AT Z 3




 





§ FLUXES AND UNCERTAINTIES AT FIR AND MIR WAVELENGTH FOR THE 30 STACKED LBGS SEDS AT Z 3. IT ALSO INCLUDES THE STACKING ANALYSIS PERFORMED FOR AZTEC (1.1MM), AND PREVIOUSLY PRESENTED IN AM16.




 




§ PHYSICAL PARAMETERS DERIVED BY SED-FITTING ANALYSIS FOR THE 30 STACKED LBGS SEDS AT Z 3





 



