© ESO 2020

1. Introduction

The detailed determination of the conditions that led to the formation of the first galaxies in the early Universe and their subsequent evolution remains a key issue in modern astronomy (e.g., Stark 2016). A truly broadband multi-wavelength perspective is likely required to robustly account for a galaxy’s growth and energy production. However, obtaining such multi-wavelength properties can be challenging due to their faint fluxes and the large distances involved.

A good example of this is the assessment of star formation rates (SFRs) in galaxies, where we must account for extinction by gas and dust in order to extract the intrinsic amount of the ultraviolet (UV) light emitted by the underlying stellar population. Deep near-infrared (NIR), optical, and UV surveys now routinely allow us to estimate unobscured SFRs down to a few M_⊙ yr⁻¹ in galaxies out to z ∼ 6–10 (e.g., Bouwens et al. 2015; McLeod et al. 2016; Santini et al. 2017; Oesch et al. 2018). A straightforward way to measure the extinction from these sources is to estimate the steepness of their UV spectra (e.g., Bouwens et al. 2012, 2014), generally characterized by fitting a power law (f_λ ∼ λ^β) to two or more rest-frame UV bands. A synthetic stellar population with solar metallicity and an age of ≳100 Myr is expected to have intrinsic β values in the range of ∼ − 2.0 to −2.2. Redward (higher β) deviations from this are thought to relate to the amount of dust extinction (reddening) and scattering that light from massive stars suffers after its emission. Blueward (lower β) deviations likely imply a very young or metal-deficient stellar population (e.g., Heap 2012; Stark 2016).

Detailed spectroscopic observations are generally required to break degeneracies between extinction, stellar age, and metallicity (e.g., Stark et al. 2013), all of which ultimately contribute to the observed UV stellar slope β. However, for fainter or more distant galaxies, this remains quite challenging (e.g., Laporte et al. 2017a; Bowler et al. 2017; Hoag et al. 2018; Hashimoto et al. 2018). Such degeneracies become particularly problematic at high redshifts, where the likelihood of young, metal-poor stellar populations and, hence, the uncertainties, are largest (e.g., Anders & Fritze-v. Alvensleben 2003; Schaerer & de Barros 2009; Eldridge et al. 2017).

A second approach for assessing extinction/absorption, as well as to examine the potential for highly or entirely obscured regions of star formation, is to measure the IR luminosity. Until recently, such observations were strongly limited in sensitivity and resolution (spatial and spectral), effectively only probing down to SFRs of ∼10–100 M_⊙ yr⁻¹ at z ∼ 1–2 (e.g., Magnelli et al. 2013). The advent of the Atacama Large Millimetre Array (ALMA), with its large collecting area and high spatial resolution capabilities, now provides the opportunity to narrow considerably the SFR gap between the UV and optical, and FIR and mm bands for galaxies across a large redshift range and, hence, make a fairer comparison between the obscured and visible light being generated.

Numerous observational studies of z ≳ 1 star-forming galaxies have been made over the years, comparing the two approaches above to well-known correlations for local galaxies (e.g., Meurer et al. 1999, hereafter M99; Reddy et al. 2006; Bouwens et al. 2009, 2016, hereafter B16; Boquien et al. 2012; Capak et al. 2015; Álvarez-Márquez et al. 2016; McLure et al. 2018; Koprowski et al. 2018). Many observers have focused on the relationships between the so-called “infrared excess” (IRX ≡ L_IR/L_UV) and UV-continuum slope (β) or stellar mass (M_⋆); such relations are often invoked to make dust attenuation corrections out to high redshifts. Most critically, while such correlations appear to be confirmed out to z ∼ 1–2, based on a variety of multi-wavelength data (e.g., Reddy et al. 2006, 2008, 2010; Daddi et al. 2007a,b; Pannella et al. 2009), it remains unclear how applicable they are at earlier times (e.g., B16).

The goal of our work here is to characterize the IR emission (individually and, given the low number of expected detections, as stacked-averages) for robust samples of Lyman-Break galaxy (LBG) candidates at z = 2–8 found in the Frontier Fields (FFs) survey¹. The FFs were initiated as Hubble (HST) and Spitzer space telescope director’s discretionary campaigns to peer as deeply as possible into the distant universe, leveraging the power of gravitational lensing from six massive high-magnification clusters of galaxies to probe to extremely faint emission levels in the most highly magnified regions (Coe et al. 2015; Lotz et al. 2017).

These fields have since been observed across the electromagnetic spectrum with, for example, Chandra, VLT/MUSE, JVLA and, of course, ALMA. We aim here to assess the IR and UV emission, stellar masses, and star formation properties of these LBG candidates, and to investigate how they compare to z ∼ 0 objects and correlations.

This paper is organized as follows. In Sect. 2, we describe the ALMA FFs observations, the LBG candidates, and their derived properties. In Sect. 3, we explain the selection criteria we applied to our candidates and the stacking procedures we utilized (ALMA image stacking and IRX stacking). In Sect. 4, we present the individual properties that we obtain for our sample, as well as the stacked values for luminosities and IRXs. Section 5 provides a comparison of our results with previously published works, as well as results not covered fully in preceding sections. Finally, we summarize our work and present our conclusions in Sect. 6. Throughout this work, we assume a cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, and Ω_Λ = 0.7.

2. Data and derived quantities

2.1. ALMA data

The inner ∼2′ × 2′ regions of the FFs, centered on the massive clusters to benefit most strongly from the boost from gravitational lensing, were observed in band 6 by ALMA through two projects, 2013.1.00999.S (PI Bauer; cycle 2) and 2015.1.01425.S (PI Bauer; cycle 3). Only five FFs clusters were completely observed by ALMA and, thus, used here. These include, from cycle 2, Abell 2744, MACSJ0416.1−2403, and MACSJ1149.5+2223 observed in 2014 and 2015 (hereafter A2744, MACSJ0416 and MACSJ1149, respectively) and, from cycle 3, Abell 370 and Abell S1063 –also designed as RXJ2248−4431– observed in 2016 (hereafter A370 and AS1063, respectively). As stated in González-López et al. (2017a), MACSJ0717.5+3745 was only partially observed (just 1 out of 9 planned executions) and, given its substantially worse sensitivity and calibration, is not useful for this work.

The mosaic data were reduced and calibrated using the Common Astronomy Software Applications (CASA v4.2.2; McMullin et al. 2007)²; details can be found in González-López et al. (2017b). Automatic reduction with the CASA-generated pipelines for A2744 and MACSJ1149 presented problems and, hence, manual and ad-hoc pipelines were used to reduce the data. For MACSJ0416, A370, and AS1063, the CASA-generated pipelines worked smoothly and were used. Observations from ALMA are characterized as visibilities (u–v plane), which must be Fourier-transformed to obtain image files (image-plane). Each visibility corresponds to an antenna pair or baseline. The visibilities (or baselines) can be weighted to produce different synthetic beamsizes and shapes. To assess the results, we applied two nominal weighting schemes, natural and taper, to the imaged (or CLEANed) datasets using CASA³. For this work, we adopted a taper parameter of $t=1\stackrel{″}{.}5$. Employing both weighting schemes offers more flexibility (and sensitivity) when searching for point-like and extended detections.

Our reductions achieved natural-weight rms⁴ errors of 55, 61, 67, 59 and 71 μJy beam⁻¹ for FFs A2744, A370, AS1063, MACSJ0416 and MACSJ1149, respectively. The resulting maps have relatively uniform rms properties over the central regions due to Nyquist sampling, but exhibit strong attenuation at the edges from the primary beam (PB) pattern. For the purposes of this work, we limited our analysis to regions of each mosaic with a PB-correction factor pbcor >  0.5, designated hereafter as the field of view (FoV) of each observation; regions with pbcor <  0.5 have substantially elevated rms values that are not very constraining. Notably, portions of the MACSJ0416 and MACSJ1149 mosaics exhibit rms variations by as much as ∼15–20% (for details, see Sect. 2.4 and Fig. 4 of González-López et al. 2017b). These variations were captured in the pbcor values used to weight individual sources in our stacking procedure (see Sect. 3.2).

Some basic properties of each dataset, including central position, are listed on Table 1. For reference, the ALMA maps of the FFs are all sufficiently deep to detect exceptional LBGs like Abell 1689-zD1, which has a band 6 flux of 0.56 ± 0.1 mJy (Knudsen et al. 2017), with S/N ∼ 8–10.

2.2. LBG candidates

Deep HST images are available in seven broadband filters as part of the FFs campaign (Lotz et al. 2017): Advanced Camera for Surveys (ACS) filters F435W, F606W, F814W (with $0\stackrel{″}{.}4$ aperture 5-σ depths of 28.8, 28.8 and 29.1 ABmag, respectively); Wide Field Camera 3 (WFC3) IR filters F105W, F125W, F140W, F160W (with $0\stackrel{″}{.}4$ aperture 5-σ depths of 28.9, 28.6, 28.6 and 28.7 ABmag, respectively). Two additional deep images were obtained with WFC3 UVIS filters F275W and F336W (with $0\stackrel{″}{.}4$ aperture 5-σ depths of ≈27.5–28.0 ABmag, depending on the cluster) as part of a supporting UV campaign (PI: Siana; Alavi et al. 2016).

Bouwens et al. (in prep.; hereafter B19) use these images to identify large samples of z ∼ 2, 3, 4, 5, 6, 7, 8, and 9 star-forming galaxies through the LBG selection technique in the FFs. Light from the foreground cluster galaxies and the intracluster medium was removed using GALFIT (Peng et al. 2002) and fitting the background light via median filtering routines, respectively, as described in B19. Source catalogs were then produced using SExtractor (Bertin & Arnouts 1996) by detecting sources in the coadded images of the four WFC3/IR filters. Colors were measured in small scalable apertures using a Kron (1980) factor of 1.2. The small scalable aperture magnitudes were then corrected to total ones based on (1) the relative extra flux seen in larger versus small scalable apertures (Kron factor of 2.5 vs. Kron factor of 1.2) and (2) the point-source encircled energy estimated to lie outside the larger scalable apertures. The correction to the total magnitude was performed based on the detection image constructed by coadding all four WFC3/IR bands. See Bouwens et al. (2015) for more details on the applied photometric procedure. Finally, B19 applied several color and signal-to-noise ratio (S/N) criteria to select LBG candidates in crude redshift bins as well as remove obvious point-like (“stellar”) contamination.

For our purposes, we did not use the B19 z ∼ 4 LBG sample, due to the lack of photometric coverage around ∼5500 Å (e.g., F555W) coupled with the potential for strong contamination by foreground galaxies in four of the five clusters considered.

B19 produced a final list of 3050 LBG candidates based on the HST cluster and parallel observations of the six FFs across all their drop-out bands, with 3029 candidates selected in the bands we use for our study (z ∼ 2, 3, 5, 6, 7, and 8). From this parent sample, we investigate the properties of the 1582 candidates located within the FoVs of five ALMA-observed FFs. Thus, all of our final results are drawn from this subset. We expect the spatial distribution of our LBG candidates to be roughly uniform over the source plane of the selected FFs. This will translate to fewer sources in highly magnified regions (near critical lines on the magnification maps) in the image plane, as we are sampling smaller intrinsic space densities. However, in a critical sense, the magnification means we probe further down the luminosity function in these regions. Thus, we expect the targets to span an interesting range in properties (e.g., magnification, SFR, M_⋆, redshift, etc.). This helps to build a statistically diverse set of LBG properties to study. Distributions of their attributes can be seen in Sect. 2.4 and later sections.

2.3. ALMA stacking considerations

We used STACKER (Lindroos et al. 2015) to perform the stacking of our candidates in the ALMA images (see Sect. 3.2). This program takes, as input, the lists of target positions (RA, Dec) and weights (for the actual stacking process). Weights are drawn from the CASA clean PB-correction map, which corresponds to the sky sensitivity over the field. This initial definition of the weight can be modified by further criteria (see Sect. 3.2). For this work, two schemes were used to weight the stacked signal according to the observed properties of the LBG candidates.

One important issue to consider is that we used information from HST and ALMA. It is possible that potential mm and submm emission in the ALMA maps may arise from a somewhat different position than the optical one, given the large span in observed wavelengths and distinct emission and extinction mechanisms at work (e.g., Goldader et al. 2002). In particular, the more dust-rich regions that could give rise to submm continuum emission would tend to attenuate embedded stars, while nearby stars in less dust-rich regions might contribute more to the observed near-IR light.

We argue, however, that such offsets were unlikely to affect our final results (i.e., the stacked flux). For one, the angular sizes of the LBG candidates are generally similar or smaller than the beam sizes of our ALMA observations (∼1″). Secondly, spatial offsets between securely detected bright mm/submm sources (e.g., submillimetre galaxies, or SMGs) and optical/NIR counterparts are generally small (e.g., $0\stackrel{″}{.}17\pm 0\stackrel{″}{.}02$González-López et al. 2017b). The offsets in SMGs, where extinction is so high that the optical emission is not detected, probably represent one extreme, while offsets in less extreme UV-selected LBGs should relatively minimal.

For the reasons above and to simplify calculations, no correction was performed regarding relative positional offsets. That being said, to obtain actually stacked fluxes from the ALMA maps in Sect. 2.10, we will consider values at larger distances due to the influence of the synthesized beam.

2.4. Photometric redshifts

As a cross-check on our LBG candidate selection, we used the photometry for each LBG candidate to obtain a photometric redshift estimate. For this purpose, we used the C++ version⁵ of the code FAST (Fitting and Assessment of Synthetic Templates; Kriek et al. 2009) with a bin size of Δz_ph = 0.001 and 500 Monte Carlo simulations per source to derive confidence levels.

The distribution of the photometric redshifts from our candidates calculated with FAST++ are shown in Fig. 1, color-coded by the drop-out band used to detect them. We find that the sub-samples do not overlap strongly and show roughly flat distributions. Only three sources, all z ∼ 5 dropouts, exhibit strong deviations between their drop-out selection band and FAST++ estimate, with z_ph ∼ 1.5; these three sources were excluded. We did not consider here the error distribution provided by FAST++, which would extend the drop-out distributions shown in Fig. 1 by ∼25%.

Fig. 1. Photometric redshift (zph) values in our sample (see Sect. 3.1) calculated from FAST++. Colors represent each drop-out band.

We also assessed the dropout candidates by comparing them to published spectroscopic redshifts. Unfortunately, only a few fields have extensive published redshift catalogs and most of our candidates are either too faint or were not targetted at appropriate wavelengths to confirm their redshifts in such surveys. Nonetheless, we compared our dropout catalogs with the VLT/MUSE redshift catalogs of Mahler et al. (2018), Lagattuta et al. (2019), Karman et al. (2017), Caminha et al. (2017), Grillo et al. (2016), and Treu et al. (2016) resulting in matches for ∼10% of our candidates (per cluster) within a $0\stackrel{″}{.}5$ circle. Among the 238 matches, only 14 candidates have strong differences between their photometric and spectroscopic redshifts; all 14 were removed.

Given the gap in the photometric redshift distribution of the LBG candidates around z ∼ 4 shown in Fig. 1, for convenience we separated our candidates into two main sub-samples: high (z_ph ≥ 4.0) and low (z_ph <  4.0) redshift. Additionally, we further subdivided the high-redshift sample into two parts: 4.0 ≥ z_ph ≤ 7.0 and z_ph >  7.0. These divisions are used for the rest of the work.

2.5. Magnification factors

Magnification factors were obtained following the procedure from Coe et al. (2015), coded as a public Python script⁶. This code obtains the values from the lensing shear (γ) and mass surface density (κ) maps that are part of the lens models products, to calculate the magnification map for each redshift. Based on FFs mass model comparisons (e.g., Meneghetti et al. 2017; Remolina González et al. 2018), we adopted the CATS (Clusters As TelescopeS) team models for our work (v4; Jauzac et al. 2014; Richard et al. 2014), as their methodology is well-documented and they appear to be among the most reliable mass models and magnifications maps of the publicly available models. With the CATS models and the photometric redshifts of the candidates as input, magnification factors (μ) were obtained from the expression:

$$\begin{array}{c}\hfill \frac{1}{\mu }=|{(1-\kappa )}^2-{\gamma }^2|.\end{array}$$(1)

To assess uncertainties associated with the magnification factors, we calculated both statistical errors using the limits of the 1-σ confidence levels of the photometric redshifts and systematic uncertainties based on the standard deviation of the magnifications of each source using four different version v4 FFs models: CATS, GLAFIC (Oguri 2010; Kawamata et al. 2016), Diego (Diego et al. 2005, 2007), and Williams (Liesenborgs et al. 2006; Jauzac et al. 2014). These uncertainties are presented in Table C.3. We note that the dispersion can be large and asymmetric since some models are not as robust as others; for this reason, we chose to incorporate a systematic error coupled with the CATS team model, rather than find a representative μ value from all the models.

Some targets can lie in positions very close to the critical curves for a given lens model and redshift, leading to extreme magnifications (μ ≳ 1000); see Fig. 2. Given the photometric redshift and lens model uncertainties, as well as the observed compact sizes of most candidates, extreme magnifications should be far less probable than moderate ones. Thus, to avoid possible spurious results when using these targets in calculations (e.g., when stacking with magnification factors as weights; see Sect. 3.2), we capped the magnification factors at μ = 10, even when models predicted larger values. This choice was driven by the fact that, after accounting for both the statistical and systematic uncertainties, > 60% of our μ >  10 candidates are compatible with magnifications of μ ≤ 10 at 1-σ confidence and > 88% at 2-σ confidence level. Coupled with the small probability that candidates can have μ >  10, we consider lower magnifications to be far more likely. As a result of this imposed ceiling, the magnification values in our full sample range from μ = 1.23 to μ = 10, with a manifest over-population at μ = 10 (due to our cap) for all three redshift bins, as seen in Fig. 3. For a comparison, we present a histogram of the unmodified μ values in Fig. 2. Finally, we note that this magnification cap should have no strong effect on our results, as higher magnifications result in no change in IRX or β values, and will only lower stellar masses, pushing candidates into a regime where we expect few detections (see Sect. 2.13); thus, some care should be taken in evaluating detections at lower stellar masses due to highly magnified sources.

Fig. 2. Distribution of magnification factors (μ) for the three photometric redshift bin samples.

Fig. 3. Capped distribution of magnification factors (μ) for the three photometric redshift bin samples. Compare to the full distributions shown in Fig. 2.

2.6. UV-continuum slope

The observed UV-continuum slope, β, is often used to assess the amount of extinction/absorption that a particular stellar population suffers, under the assumption that a nominal intrinsic UV slope is typically β₀ ≈ −2.0 to −2.2, for constant star formation, with high values indicating higher attenuation. For each candidate in our sample, nine-band HST photometry was used to obtain the observed values of β. Several methods have been developed to calculate β from different photometric bands (for a review, see Sect. 2 of Rogers et al. 2013 and Sect. 2.7 of McLure et al. 2018). Here, we adopted a simplistic approach using the bands (and the flux in them) that fall in the expected UV-continuum spectral region (e.g., ≲3000 Å) assuming the previously derived redshift. A power law (F_λ ∝ λ^β) was fit to the rest-frame UV photometry using the Python implementation of the Affine Invariant Markov chain Monte Carlo Ensemble sampler (emcee;Foreman-Mackey et al. 2013). In particular, we adopted the functional form chosen by Castellano et al. (2012), which is:

$$\begin{array}{c}\hfill m_i=-2.5\times (\beta +2.0)\times log({\lambda }_i)+c,\end{array}$$(2)

where m_i is the AB magnitude in the ith band (Oke & Gunn 1983) at an effective wavelength λ_i and c is the intercept. As priors for the model fitting, we used the outputs of a simple maximum likelihood estimator with Eq. (2). For each LBG candidate, 2500 iterations were performed per each one of the 100 “random-walkers” which were set for this procedure. From them, we obtained the most probable β values and the limits of their 1-σ credible intervals.

A comparison of the UV slopes (β) and magnification-corrected magnitudes, as well as their overall distributions, is presented in Fig. 4. The three broad divisions in photometric redshift do not show any particular trend between β and redshift. A finer binning of the targets according to photometric redshift is shown in Fig. 5, where it can be seen that the UV-slopes of our LBG candidates are generally consistent with being more or less constant between z ∼ 1–8, within the large dispersion. Previous works such as Bouwens et al. (2012, 2014); Finkelstein et al. (2012) have reported mild evolution in β for z_ph ≳ 4 LBG candidates due to a possible increase in dust extinction with time. This weak evolution lies within the dispersion of our sample and, thus, we can neither confirm nor reject it.

Fig. 4. Comparison of β and magnification-corrected magnitudes for selected sources. Top and right panels: histograms of UV slope and magnitude distributions. Colors represent photometric redshift subsamples, as described in the legend and in Sect. 3.

Fig. 5. Distribution of β values according to photometric redshift. Heights of boxes represent the 25% and 75% quartiles of the data. Horizontal lines inside the box indicate the median value for each redshift bin. Vertical error bars span the central 2-σ of the data. Numbers above the median in each box state the number of LBG candidates assigned to each bin. Even though there is not a zph ∼ 4 band from drop-out selection, there are candidates in that bin.

For the purposes of this work, we defined the UV flux or luminosity (F_UV, L_UV) to be that measured at 1600 Å (following, among others, Madau & Dickinson 2014, who suggested that UV wavelengths between 1400 Å and 1700 Å provided a reasonable estimate). In our case, we used the photometric band which lies closest to that rest-frame wavelength.

The UV slope can also be related to the dust attenuation factor, A_λ, as in B16 and Calzetti et al. (2000). For this work, we favored the relation found by Calzetti et al. (2000):

$$\begin{array}{c}\hfill A_{1600}=2.31\times \beta +4.85,\end{array}$$(3)

which was similarly assessed at 1600 Å for low-redshift galaxies.

2.7. Stellar masses

Stellar masses were estimated using FAST++, which fits stellar population synthesis templates to photometric data. The input values were the magnitudes from our LBG catalogs as well as the photometric redshifts (also determined with FAST++). For this work, we assumed Bruzual & Charlot (2003) stellar spectral energy distributions (SEDs) with a Chabrier Stellar IMF (Chabrier 2003). We assumed an approximately constant SFR in modeling the star formation history, effectively realized by setting log₁₀(τ/yr) = 11 with an exponentially declining star formation history (SFR ∝ exp(−t/τ)) and a metallicity of 0.2 Z/Z_⊙. Finally, a Calzetti et al. (2000) dust attenuation law with a range of 0.0 ≤ A_V ≤ 1.0 was adopted. The code outputs, apart from other relevant properties, a stellar mass estimate for each target. The above parameter choices have a sizeable impact on inferred quantities such as the stellar population age (> 0.3–0.5 dex) but do not strongly impact the inferred stellar masses (≳0.2 dex).

To obtain the magnification-corrected stellar masses, the values given by FAST++ were divided by the magnification factors.

For the rest of this work, we refer to the magnification-corrected stellar mass simply as stellar mass. The distribution of best-fit values for our three photometric redshift bins can be seen in Fig. 6. Factoring in the 1-σ confidence intervals on the stellar mass (see Fig. 12), the full range spans ∼ 10^5.6 M_⊙ to ∼ 10^10.2 M_⊙.

Fig. 6. Stellar masses in our sample. Sample has been divided according to the photometric redshift bins defined in Sect. 3. Vertical dark line represents the approximate completeness limit from M18 (See Sect. 5.1.5).

2.8. (Specific) Star formation rates

One of the by-products of FAST++ is a SFR estimation. As with stellar mass, SFR values were corrected by the magnification factor (μ).

For the rest of this work, we refer to the magnification-corrected SFR as SFR. From the SFRs and stellar masses, specific SFRs, or sSFRs, are obtained as:

$$\begin{array}{c}\hfill {log}_{10}(\mathrm{sSFR}/{\mathrm{yr}}^{-1})={log}_{10}\left(\frac{\mathrm{SFR}/M_{\odot }{\mathrm{yr}}^{-1}}{M_{\star }/M_{\odot }}\right).\end{array}$$(4)

2.9. ALMA primary-beam corrections

To obtain images from mosaics of interferometric data, each element of the observation has to be corrected by a combination of the sensitivity of every pointing in the observation and the change of sensitivity across the mosaic. These two elements constitute the PB corrections of the observed maps.

With interferometric data, the deconvolution from the u–v visibility plane to the image plane includes a division (deconvolution) by these primary-beam correction factors. One of the ALMA pipeline data products is a normalized map of sensitivities, which incorporates the primary-beam correction factors, ranging from no (0) to full (1) sensitivity (see, for instance, Thompson et al. 2017; Wilson et al. 2013).

2.10. ALMA peak fluxes

For simplicity, we adopted peak flux measurements, $F{}_{\mathrm{ALMA},\mathrm{peak}}^{\mathrm{indiv},\mathrm{obs}}$, since integrated fluxes require an assumption about the flux distribution shape. To assess these peak fluxes within the ALMA maps, we searched for the pixel with the maximum value within a $0\stackrel{″}{.}5\times 0\stackrel{″}{.}5$ box (i.e., comparable to one synthesized beam) centered at the position of each LBG candidate. This procedure attempts to account for the influence of the synthesized beam, as well as possible extended emission, in the ALMA maps. We corrected this flux for the PB attenuation (i.e., accounting only for the properties of the observations) as follows:

$$\begin{array}{c}\hfill F_{\mathrm{ALMA},\mathrm{peak},\mathrm{pbcor}}^{\mathrm{indiv},\mathrm{obs}}=\frac{F_{\mathrm{ALMA},\mathrm{peak}}^{\mathrm{indiv},\mathrm{obs}}}{pbco{r}_{\mathrm{ALMA}}^{\mathrm{indiv}}}·\end{array}$$(5)

Likewise, we related the rms error at the position of an individual source to the field rms (rms${}_{\mathrm{ALMA}}^{\mathrm{cluster}}$) listed in Sect. 2 for each studied cluster, as

$$\begin{array}{c}\hfill {\mathrm{rms}}_{\mathrm{ALMA},\mathrm{pbcor}}^{\mathrm{indiv}}=\frac{{\mathrm{rms}}_{\mathrm{ALMA}}^{\mathrm{cluster}}}{pbco{r}_{\mathrm{ALMA}}^{\mathrm{indiv}}}·\end{array}$$(6)

The bulk of our candidates have ALMA fluxes comparable to the rms values of their respective maps, but a few are associated with brighter peak fluxes. For this reason, we want to define clearly which targets are detected and for which we only have upper limits. As a first conservative approach, we searched for LBG candidates with S/N above 5.0 in each image, which roughly corresponds to the blind detection limit for the ALMA-FF maps (González-López et al. 2017b). This high S/N limit arises in the context of having large maps with ≈1.7 × 10⁷ pixels yet only a handful of highly secure detections per field. The map noise is approximately Gaussian (González-López et al. 2017b), meaning that there should be roughly 45896, 1077, and 9 pixels above 3, 4, and 5 times the rms, respectively, in each map. Excess numbers of pixels above these expectations imply real sources. We defined here the S/N as:

$$\begin{array}{c}\hfill S/N=\frac{F_{\mathrm{ALMA},\mathrm{peak}}^{\mathrm{indiv},\mathrm{obs}}}{{\mathrm{rms}}_{\mathrm{ALMA}}^{\mathrm{cluster}}}=\frac{F_{\mathrm{ALMA},\mathrm{peak},\mathrm{pbcor}}^{\mathrm{indiv},\mathrm{obs}}}{{\mathrm{rms}}_{\mathrm{ALMA},\mathrm{pbcor}}^{\mathrm{indiv}}}·\end{array}$$(7)

None of our targets fulfills this first condition, with a maximum value of S/N = 4.21 for a candidate in AS1063.

The blind detection limit, however, is with respect to a search of all positions on the map. Nevertheless, since we know the positions of the 1582 LBG candidates and they comprise only a small fraction of the overall map area (≈1.1 × 10⁵ pixels)⁷, a more realistic estimate of the detection significance is to evaluate the False Detection Rate (FDR or p_FDR, Benjamini & Hochberg 1995; Benjamini & Yekutieli 2001) for each ALMA map. As described in Miller et al. (2001) and Hopkins et al. (2002), the FDR is different from other thresholding methods in that it constrains the fraction of false detections compared with the total number of detections rather than the fraction of pixels falsely detected over the total number of pixels. Given its definition, the FDR does not depend on the distribution of sources and, thus, we are not forced to assume a specific behavior for them.

To this end, following the procedures outlined in Muñoz Arancibia et al. (2018), we generated 1000 simulated maps for each ALMA field with a normal distribution in units of signal-to-noise. From these we extracted the same number of simulated peak fluxes per cluster as we did for the LBG candidates, again choosing the highest peak flux within a square of $0\stackrel{″}{.}5$ on a side. We defined p_FDR(S/N) to be the fraction of simulated maps of a specific cluster where at least one sampled pixel was found above a given S/N. Figure 7 shows the FDRs for our five ALMA maps.

Fig. 7. False detection rate, pFDR, for the five ALMA maps. Vertical dashed lines denote the highest detected S/N among the LBG candidates in each ALMA map. The vertical solid line denotes our adopted S/N cutoff of 4.1, which equates to a FDR around 15% among the cluster fields.

Based on the FDRs, we find that sources with S/N ≳ 4.1 have a relatively low (≲15%) chance of being false. For simplicity and uniformity, we considered all LBG candidates above this limit to be detected, while the LBG candidates below this were treated as upper limits. We calculated individual detected peak fluxes following Eq. (5), while n-σ upper limits were calculated as

$$\begin{array}{c}\hfill F_{\mathrm{ALMA},\mathrm{peak},\mathrm{pbcor}}^{\mathrm{indiv},\mathrm{obs}\mathrm{n}\text{-}\sigma \lim }=\frac{(F_{\mathrm{ALMA},\mathrm{peak}}^{\mathrm{indiv},\mathrm{obs}}>0)+n\times {\mathrm{rms}}_{\mathrm{ALMA}}^{\mathrm{cluster}}}{pbco{r}_{\mathrm{ALMA}}^{\mathrm{indiv}}},\end{array}$$(8)

where the > 0 expression indicates the fact that the observed peak flux from the ALMA map is only used if it is greater than zero. This implies that no single candidate will have a 1-σ upper limit lower than the noise level of the map to which it belongs. The incorporation of local map noise, in addition to the average rms, yields a more conservative upper limit.

2.11. UV luminosities

As in Sect. 2.6, we defined the UV flux and luminosity as that at 1600 Å, ensuring that appropriate rest-frame and magnification corrections are applied for the best-fit photometric redshift.

2.12. IR luminosities

With only one point of constraint from ALMA, the procedure to estimate the IR luminosity is more model-dependent than for the UV bands. For this, we fit a graybody spectrum (e.g., Casey 2012; Schaerer et al. 2013) to the ALMA photometric data. We adopted evolving values for dust temperature following the redshift-dependent formulation from Schreiber et al. (2018). We caution, however, that this relation is only fit up to z = 4 and, thus, extrapolations may be problematic. The distribution of dust temperatures of our candidates ranges from ∼30 K to ∼65 K, which is in line with the z ≥ 5 simulations from Ma et al. (2019) and the 2 ≤ z ≤ 4 simulations from Liang et al. (2019). We also considered typical fixed values of 2.0 for the mid-IR power-law slope, and 1.6 for the emissivity (e.g., best-fit values for the GOALS survey; Casey 2012). For simplicity, we adopted the same shape for every LBG candidate. The best-fitted rest-frame SED is integrated between 8 μm and 1000 μm to yield the rest-frame IR luminosity. In practical terms, we defined a scale factor $f_{\mathrm{IR}}^{\mathrm{ALMA},\mathrm{peak},\mathrm{pbcor},\mu \mathrm{cor}}$ to convert observed ALMA peak flux to the magnification-corrected, rest-frame IR luminosity as

$$\begin{array}{c}\hfill f_{\mathrm{IR}}^{\mathrm{ALMA},\mathrm{peak},\mathrm{pbcor},\mu \mathrm{cor}}=\left(\frac{F_{\mathrm{ALMA},\mathrm{peak},\mathrm{pbcor}}/\mu }{F_{\mathrm{SED}[1.14\mathrm{mm}/(1+z)]}}\right)·\end{array}$$(9)

We chose this method over FAST++ or magphys (Multi-wavelength Analysis of Galaxy Physical Properties; da Cunha et al. 2008) SED fitting to obtain IR luminosity estimates due to the fewer number of free parameters (e.g., dust temperature, SED templates), which made for a more straightforward implementation and interpretation. In general, the luminosities derived from the best-fit modified blackbody to the ALMA data are factors of 10–100 higher than rest-frame UV/optical based estimates from FAST++ or magphys. Our estimates are presumably more robust for the few detections, while the upper limits should be considered as very conservative.

To test this method, we calculated IR luminosities for the sources reported by Aravena et al. (2016) using their ALMA (Band 6) flux measurements. Our results lie within ∼0.5 dex of theirs, which were obtained with magphys. These results demonstrate that we can obtain relatively reliable IR luminosities from the graybody spectrum.

In addition to the aforementioned corrections for redshift and magnification, the IR luminosities (or fluxes) have an additional dependence on the redshift of the candidate due to the impact of the CMB temperature on the dust properties. Following the procedure of da Cunha et al. (2013), the derived IR luminosities were divided by the factor

$$\begin{array}{c}\hfill g_{\nu }^{\mathrm{CMB}}=[1-\frac{B_{\nu }(T_{\mathrm{CMB}}(z))}{B_{\nu }(T_{\mathrm{dust}},z)}],\end{array}$$(10)

where B_ν(T_dust) and B_ν(T_CMB, z) correspond to the source and CMB blackbody contributions at the observed frequency and redshift of the source, respectively.

Errors were propagated according to Eq. (6) (in units of luminosity), applying the same corrections (e.g., redshift, magnification, and CMB temperature). We calculated IR luminosity upper limits as the peak value at source location plus n-sigma:

$$\begin{array}{c}\hfill L_{\mathrm{IR},\mathrm{peak},\mathrm{pbcor}}^{\mathrm{indiv},\mathrm{obs}n\text{-}\sigma \lim }=L_{\mathrm{IR},\mu \mathrm{cor}}^{\mathrm{peak},\mathrm{pbcor}}+n\times {\mathrm{rms}}_{\mathrm{IR},\mathrm{up}\lim }^{\mathrm{pbcor},\mu \mathrm{cor}}[L_{\odot }],\end{array}$$(11)

and generally adopted 1-σ as the credible interval used. This upper limit formalism is also adopted for other quantities throughout this work (e.g., IRX).

2.13. IRX relations

Sensitive millimetre facilities such as Herschel and ALMA have only become available in the last decade. Prior to these, it was generally difficult to measure IR luminosities for distant galaxies, and indirect methods were employed to understand and predict the IR emission. Principal among these is the so-called IR excess ratio (IRX), which is loosely defined as the ratio between the IR and UV luminosities (or fluxes) of a source (in this case, a galaxy). One of the most utilized definitions was developed by B16, which relates the UV and IR fluxes as:

$$\begin{array}{c}\hfill \mathrm{IRX}=\frac{F_{\mathrm{IR}}}{F_{\mathrm{UV}}}\end{array}$$(12)

where F_IR is the rest-frame 8–1000 μm IR flux and F_UV is the rest-frame 1600 Å UV flux, both of them corrected for magnification factors. This can be trivially extended for rest-frame luminosities instead of fluxes. These relations were developed using local galaxy data, but have been tested on a variety of distant (mostly massive) galaxy samples.

Similar to the IRX–β relations, there have been a large number of studies arguing that the total stellar mass of a galaxy is strongly related to the degree of dust extinction and, hence, IRX. We highlight four recent published correlations between IRX and stellar mass by Heinis et al. (2014, hereafter H14), Fudamoto et al. (2017, hereafter F17), McLure et al. (2018, hereafter M18), and B16.

Finally, B16 also derived a “consensus” IRX–M_⋆ relation from a variety of previous studies in the redshift range z ∼ 0 to z ∼ 3 (e.g, Pannella et al. 2009; Reddy et al. 2010).

The various IRX–M_⋆ relations have relatively similar slopes and exhibit a typical dispersion of up to ∼1 dex, excluding the strong deviation of H14 above 10¹⁰ M_⊙. As such, they provide a potentially useful means of predicting dust attenuation as a function of stellar mass.

3. Methods

3.1. Target final sample

With all of the derived quantities in hand (Sect. 2), we now address the selection of the LBG candidate sample, in order to improve the reliability and trustworthiness of the estimated physical properties and stacking results.

We began by discarding a handful (7) of LBG candidates with UV-slopes β <  −4.0 or β ≥ 1.5 (see Fig. 4). These extremely low or high values arise at faint magnitudes, have large error bars, and are physically implausible. This is qualitatively comparable to a (UV) color selection.

Before stacking, we also excluded 408 LBG candidates in close proximity but unrelated to any ≥4-σ detected sources in the ALMA maps in order to avoid contamination in the stacked signal. We conservatively adopted a circular exclusion region equal to five times the major axis of the natural-weighted synthesized beam for each map (i.e., $3\stackrel{″}{.}2$–$7\stackrel{″}{.}6$). We additionally removed all LBG candidates with primary-beam correction factors lower than 0.5 (see Sect. 2.9), as the edges of the ALMA maps have considerably higher noise and other observational artifacts that can adversely affect the sensitivity of the stacking.

Based on the FDR assessment in Sect. 2.10, we also identified two LBG candidates associated with ALMA detections at S/N ≳ 4.1, adopting a matching radius of 0.5 times the major axis of the natural-weighted synthesized beam for each map (i.e., $0\stackrel{″}{.}3$–$0\stackrel{″}{.}8$). These sources, along with their key attributes are listed in Table 4 and they were not included in the main stacking and have been treated separately. For comparison, the typical positional uncertainties between ALMA and HST sources are $\lesssim 0\stackrel{″}{.}1$ (e.g., ≲10% of the beam size in González-López et al. 2017b).

Finally, we considered whether LBG candidates were multiply imaged. We did not want to double-count the same source, as this could have potentially distorted our stacking results. Thus, we removed all multiple images. To determine whether a candidate was multiply imaged, we matched the positions of our LBG candidates against the multiple-image catalogs from the CATS team (v4; see Sect. 2.5), which comprise a compilation of secure multiple images found via HST or ground-based spectroscopic confirmation (e.g., Smith et al. 2009; Merten et al. 2011; Zitrin et al. 2011, 2013; Jauzac et al. 2014; Richard et al. 2014; Kawamata et al. 2016, 2018; Caminha et al. 2017; Lagattuta et al. 2017; Mahler et al. 2018). In total, we removed 53 LBG candidates with positions conservatively lying within $0\stackrel{″}{.}5$ radius of a known multiple images (23 lie within $0\stackrel{″}{.}25$).

We summarize our selection criteria in Table 2, which resulted in a sample of 1580 undetected LBG candidates to stack: 383 from A2744; 369 from MACSJ0416; 315 from MACSJ1149; 121 from A370; and 391 from AS1063. For some specific results below, to avoid problems related to combining values spanning several orders of magnitude (e.g., the weights from Sect. 3.2), we restricted the sample even further; for instance, when considering stacking in bins of M_⋆, we discarded a handful of very low-mass LBGs and only considered 1569 candidates.

3.2. Stacking

To perform the stacking process for our ALMA data, we used the STACKER code developed by Lindroos et al. (2015). It can stack interferometric data in both the u–v (visibilities) and image domains. For the image domain, the code uses median or mean stacking with weights. These weights can be fixed a priori or obtained from the PB-correction data present in ALMA datasets. The product of this stacking process is an ALMA image file. In the u–v domain, the stack aligns the phases and then adds up the weighted visibilities.

We adopted four different weighting schemes for the stacking code and further analysis: no or equal weights for all sources; PB correction pbccor-weighting; (magnification-corrected) UV flux F_UV and pbcor weighting; and magnification μ and pbcor-weighting. For the equal weight scenario, the weight factor ($W_k^{\mathrm{no}}$) is simply a constant of unity for all k sources.

For the pbcor-weighting scenario, the sensitivity maps were used, with the weight factor given by:

$$\begin{array}{c}\hfill W_k^{\mathit{pbcor}}={\left(pbco{r}_{\mathrm{ALMA}}^{\mathrm{indiv}}\right)}^2.\end{array}$$(13)

This scheme simply counteracts the effects of the primary-beam correction on the determination of ALMA peak fluxes and, hence, enhances the contributions from the sources with the lowest rms values.

For the UV-flux F_UV weighting scenario, the factor has the form:

$$\begin{array}{c}\hfill W_k^{\mathrm{UV}}={\left(pbco{r}_{\mathrm{ALMA}}^{\mathrm{indiv}}\right)}^2\times F_{\mathrm{UV}}^2.\end{array}$$(14)

This scheme should enhance the contribution from sources that show a higher ultraviolet flux and, by extension, higher star formation activity (and possibly stellar masses due to the star-formation main sequence), in addition to the pbcor correction. We caution that this scheme could bias the stacking results toward sources that are less obscured and are more likely to lie closer to the B16 IRX–β relation.

Likewise, for the magnification μ-weighting scenario, the weight factor is:

$$\begin{array}{c}\hfill W_k^{\mu }={\left(pbco{r}_{\mathrm{ALMA}}^{\mathrm{indiv}}\right)}^2\times {\mu }^2.\end{array}$$(15)

This weight configuration takes advantage of the magnification power of the galaxy clusters, which can amplify the influence of faint or less obscured sources in the final results, in addition to the pbcor correction.

We expected some S/N variations among the different weighted stacks since they include different contributions of ALMA flux into the final results. The adopted weighting schemes might have inadvertently downweighted contributions from LBG candidates with higher individual S/N values. For instance, by favoring properties that are not directly expressed in the ALMA data, we may have been selecting against the most dust enshrouded candidates.

This stacking produces, ultimately, an image file. In this image, the stacked flux from the candidates is present in the central pixel if the objects are point-like. If highly extended or offset sources are part of stacked targets, other considerations must be taken into account; for instance, if extended, we would want to adopt an appropriate beam shape, or if offset, we would want to calculate the center of each target from the ALMA observation itself, rather than adopting the HST catalog position. As stated in Sect. 2.3, we did not expect UV and IR offsets to be a preponderant issue here and, thus, calculated the stacking results adopting the individual UV (HST) positions of the LBG candidates.

After STACKER was run for each data configuration, every stacked image was inspected to determine if a detection has been achieved. We calculated the detection levels for each stacked image using the procedure described by González-López et al. (2017b), in which peaks (sources) with S/N > 5-σ are iteratively discarded until we arrive at a stable rms noise value.

On the other hand, to obtain stacked values of IRX, a different method must be employed in which the stacking of ALMA observations is not directly utilized.

Following previous discussions from Bourne et al. (2017) and Koprowski et al. (2018), and taking into account the weights we are using, the appropriate method to determine stacked IRX values is

$$\begin{array}{c}\hfill \overline{\mathrm{IRX}}=\overline{\left(\frac{L_{\mathrm{IR}}}{L_{\mathrm{UV}}}\right)}\end{array}$$(16)

for each subsample in bins of redshift, stellar mass and UV-slope.

We adopted this indicator since it is non-trivial to know, a priori, how the UV and IR luminosities are related. Thus, we stacked the individual IRX values and not the separate luminosities. The calculated IR luminosities were provided using the procedure described in Sect. 2.12.

In the case of upper limits, we stacked, separately, the peak IRX values and their 3-σ error values. Then, we combined them to obtain the final stacked upper limits. That is:

$$\begin{array}{c}\hfill \overline{{\mathrm{IRX}}^{n\text{-}\sigma \lim }}=\overline{\left(\frac{L_{\mathrm{IR},\mu \mathrm{cor}}^{\mathrm{peak},\mathrm{pbcor}}}{L_{\mathrm{UV}}}\right)}+n\times \mathrm{\Delta }\overline{\left(\frac{L_{\mathrm{IR},\mu \mathrm{cor}}^{\mathrm{peak},\mathrm{pbcor}}}{L_{\mathrm{UV}}}\right)}·\end{array}$$(17)

Finally, to investigate the relation between IRX and other parameters, the target stacking was binned as a function of three different quantities; UV-slope, stellar mass, and redshift.

With UV-slope, targets were stacked in five bins and, for stellar mass, in nine bins. Candidates with stellar masses less than 10^6.0 M_⊙ were excluded from stacking calculations because of their very low expected luminosities and low numbers. For redshift, three sub-samples were utilized. These divisions were adopted considering the apparent distribution of redshift values shown in Fig. 1.

The choice of bin widths was made as a compromise between having sufficient numbers of sources to reap the benefits of stacking and using equal-width bins in parameter space to facilitate interpretation. For the latter reason, we did not attempt to have a similar number of elements per bin. The bins are presented in Table 3, while the number of sources per bin are presented in Col. 3 of Tables B.1–B.4. We can see that the uncertainties for the β and M_⋆ (Tables C.3 and C.4) are small enough to not pose major problems to the binning of the sources.

3.3. Considerations on stacking weighting

Stacking of the ALMA data and IRX values can potentially constrain the average properties of a sample well below the formal detection limits for individual sources. The obtained values, however, should be regarded with some reservations. For one, the average properties can be skewed by a few outliers, since we are not individually detecting objects. Secondly, we employed μ and F_UV weighting schemes (see Sect. 3.2) with the aim to improve our sensitivity. The downside of weighting, however, is that our stacked result can be biased toward the candidates with the highest weights.

As an example, consider the case of IRX stacking with F_UV weighting. We can expect that stacking results will be skewed toward candidates with higher UV luminosities and, hence, lower stacked IRX values, which is not, necessarily, an expression of the behavior of most LBG candidates. Thus, any stacked IRX value has to be considered as a manifestation of the influence of the candidates with the highest weights and not as a true expression of the overall trend from the full studied sample.

4. Results

We describe below the main results obtained both for the individually detected sources reported in Sect. 2 and from the stacking of the ALMA and IRX values of our sample.

4.1. Individual results

Based on the individual luminosities obtained using the graybody SED and our HST photometry, we derive IRX values (or upper limits) and compare them with previously calculated properties for each LBG candidate. We focus our comparisons on the UV slopes and stellar masses of the candidates. Some key properties for our ALMA detections are listed in Table 4. A broader set of properties for all our LBG candidates are listed in the tables of Appendix C.

4.1.1. ALMA peak fluxes

The mean and peak S/N distributions for the 1582 LBG candidates are shown in Fig. 8. As already mentioned in Sect. 2.10, all our targets exhibit S/N values lower than | ± 5.0|. The mean S/N distributions for each redshift bin are centered around ∼0 as expected, while the peak S/N distributions are centered around ∼1 as a result of selecting the peak pixel which arises within half a beamwidth; this conservatively biases the maximum flux associated with a candidate to higher values. Both distributions appear roughly Gaussian.

Fig. 8. Mean (left) and peak (right) signal-to-noise ratios (S/N) for our candidates in the ALMA maps. The LBG candidates are separated into three photometric redshift sub-samples, represented by distinct colors. The mean value is centered around S/N ∼ 0 (vertical dark line) and is roughly Gaussian. The peak values are centered around S/N ∼ 1, rather than S/N ∼ 0 (vertical dark line), due to the selection of the peak pixel which arises within half a beamwidth; this conservatively biases the maximum flux associated with a candidate to higher values.

From our sample, we find two (2) candidates with ${S/N}_{\mathrm{peak}}^{\mathrm{indiv}}>4.1$ (see Table 4). Based on the results from Sect. 2.10, we expect ≈0.3 candidates to be false positives at this S/N (ρ_FDR = 0.15) and, thus, consider the two detections to be real.

4.1.2. UV and IR luminosities

Following the steps described in Sects. 2.11 and 2.12, we utilized HST photometry to calculate UV luminosities for each LBG candidate and a graybody SED to calculate the IR luminosities, re-scaled by the individual ALMA peak fluxes. The vast majority of the latter are upper limits. The distributions of the individual UV and IR luminosities (3-σ upper limits) are shown in Figs. 9 and 10, respectively.

Fig. 9. UV luminosities in our sample. The LBG candidates are separated into three photometric redshift sub-samples, represented by distinct colors.

Fig. 10. IR luminosities (3-σ upper limits from our ALMA maps) in our sample. The LBG candidates are separated into three photometric redshift sub-samples, represented by distinct colors.

The magnification-corrected observed UV luminosity 3-σ upper limits of the LBG candidates span a range from ∼10^7.8–10^10.8 L_⊙, effectively probing apparent SFRs between ∼0.02 and 20 M_⊙ yr⁻¹ (e.g., Calzetti 2013). We see a peak at around ∼10⁹ L_⊙ for the two lower redshift bins (z <  4 and 4 ≤ z <  7), while we see a relatively flat distribution between ∼10^8.5 and 10^10.5 L_⊙ for the higher redshift bin. In general, the UV luminosities probed here are lower than the values presented in other works (e.g., Narayanan et al. 2018; Reddy et al. 2018).

The magnification-corrected IR luminosity limits of the LBG candidates exhibit a somewhat different behavior from the UV luminosities. Due to the nature of the K-correction on the long-wavelength side of the graybody SED, high redshifts probe somewhat lower IR luminosities. Specifically, we find that the z <  4, 4 ≤ z <  7, and z ≥ 7 bins are centered around values of ∼10^11.9, ∼10^11.7, and ∼10^11.7 L_⊙, respectively. Given our imposed maximum magnification of 10, coupled with the relatively uniform rms limits, we see that each photometric redshift subsample spans roughly 1.5 dex in luminosity (without accounting for outliers). Thus, all our redshift bins probe IR luminosity limits of ∼10^11.1–10^12.5 L_⊙, or equivalently 20–400 M_⊙ yr⁻¹ (e.g., Hao et al. 2011; Calzetti 2013).

Comparing the UV and IR luminosity limits, it is clear that the UV data generally probes to much lower effective SFRs. Thus, our current individual ALMA constraints are only able to rule out the possibility of rather extreme obscured star formation events associated with any of the LBG candidates.

The two detected LBGs have UV and IR luminosities in the range L_UV ∼ 10^8.7–10^8.8 L_⊙ and L_IR ∼ 10^11.3–10^11.9 L_⊙, respectively. Relating these in terms of SFRs, the detected LBGs have ∼2–3 dex more obscured than unobscured star formation present.

4.1.3. IRX–β relation

With the UV and IR luminosities in hand, we can compare IRX limits to the UV-slope β, as shown in Fig. 11. We color-code the LBG candidates as functions of redshift, magnification, sSFR, M_⋆, and L_UV, as well as show the local IRX–β relations presented in Sect. 2.13.

Fig. 11. Comparison of infrared excess (IRX) 3-σ upper limits and UV-slopes (β) for our LBG candidates. Downward arrows have 1-σ length. From top to bottom panels: colors represent: photometric redshift (zph), magnification factor (μ), star formation rate (SFR), stellar mass (M⋆) and UV Luminosity (LUV). Local IRX–β relations presented in Sect. 2.13 are shown for reference. Blue crosses represent the two detections.

The main trends we see in the IRX–β diagram are with the FAST-derived quantities SFR and M_⋆ (third and fourth panels), where stronger upper limits tend to lie to the lower right, closer to the local relations (and weaker limits tend to lie further away from local relations). This is due in part to observation bias, coupled with the M_⋆–SFR main sequence relation. We detect LBG candidates spanning ∼3 dex in m_UV or L_UV (bottom panel), while our IR limits only span 1 dex. Thus, the highest M_⋆–SFR sources have the lowest IRX limits, and the lowest ones have the highest IRX limits. This trend extends into the z_ph and μ panels with lower redshift and higher μ sources (i.e., lower L_UV candidates) having higher IRX limits. There appears to be a mild intrinsic trend between higher (redder) β values and higher M_⋆ (see Sect. 5.1.5 for further details).

While the vast majority of limits lie above the local relations, we find 3 LBG candidates located completely below at least one relation. Given the dispersion in these local relations, however, all we can say is that our individual limits remain consistent with the relations.

4.1.4. IRX–M_⋆ relation

We can also compare the IRX limits and stellar masses M_⋆ of our LBG candidates, depicted in Fig. 12. Again, we color-code the LBG candidates as functions of redshift, magnification, sSFR, β, and L_UV, and show several IRX–M_⋆ relations from Sect. 2.13.

Fig. 12. Comparison of infrared excess (IRX) 3-σ upper limits and stellar masses M⋆ for our LBG candidates. Downward arrows have 1-σ length. From top to bottom panels: colors represent: photometric redshift (zph), magnification factor (μ), star formation rate (SFR), UV slope (β) and UV Luminosity (LUV). Local IRX–M⋆ relations presented in Sect. 2.13 are shown for reference. Blue crosses represent the two detections.

We see a number of trends in the IRX–M_⋆ diagram as functions of μ (second panel), sSFR (third panel), β (fourth panel), and L_UV (fifth panel). Unsurprisingly, higher magnifications allow us to probe lower stellar masses. M_⋆ is related to sSFR and L_UV following the star-formation main sequence. Here we now see more clearly a M_⋆ and β trend, such that more massive systems (which have built up more metals and dust) tend to show higher extinction.

In this case, unlike the IRX–β trends, all of our 3-σ upper limits lie completely above the relations. Factoring in the dispersion in these relations, our individual limits remain consistent with the relations. The massive and luminous LBG candidates that lie closest to the relations all have high (z ≳ 5) photometric redshifts and low magnifications and, hence, comprise the rare, bright end of the high-z population.

4.2. Stacking results

To gain further insights into the LBG population, we used STACKER to perform u–v stacking on all five ALMA cluster datasets. Some tests were applied to the stacking method and their details are presented in Appendix A.

Importantly, our tests demonstrate the capabilities of STACKER in substantially reducing the noise levels compared to the nominal natural-weight CLEANing rms (e.g., from 55 μJy to 90 μJy to stacked rms errors as low as 2 μJy, which is close to the theoretical limit). Comparable results are achieved with image stacking, and give us confidence in the LBG stacking results presented below.

From here, we turned to stacking the undetected LBG candidates in the three broad photometric redshift bins as functions of UV-slope binning and stellar mass binning. The u–v stacking results are presented in Tables B.1–B.3, respectively. Stacked image stamps for two example bins are presented in Fig. 13 (4.0 ≤ z <  7.0 and −2.0 ≤ β <  −1.0) and Fig. 14 (4.0 ≤ z <  7.0 and 9.0 ≤ log(M_⋆/M_⊙) <  9.5). With the large number of undetected LBG candidates in some bins, we achieve stacked rms values as low as ≈5 μJy. This highlights the power of stacking to reduce the errors and increase the signal strength (S/N) accordingly by $\sim \sqrt{N}$.

Fig. 13. Example u–v stacked image stamps for 227 undetected LBG candidates in the range of −2.0 ≤ β <  −1.0 and 4.0 ≤ z <  7.0. Details same as Fig. A.1. Color scale spans −125 μJy to +125 μJy range.

Overall, only one bin among all of the stacks achieves a S/N high enough to be considered a detection (227 F_UV-weighted sources in the range 4.0 ≤ z <  7.0 and −2.0 ≤ β <  −1.0 with ${S/N}_{\mathrm{peak}}^{\mathrm{stack}}=4.24$ for the natural-weight CLEANing (Fig. 13). We treat this result with caution, however, since is not replicated in any other weighting schemes and CLEANing configurations for the same targets, which yield ${S/N}_{\mathrm{peak}}^{\mathrm{stack}}=-1.14$ to 3.04. As seen in Tables B.1–B.3, there are only a few bins with even marginally significant signal (i.e., S/N ≥ 3.00), the highest being ${S/N}_{\mathrm{peak}}^{\mathrm{stack}}=3.92$ for 19 F_UV-weighted candidates with stellar masses and reshifts in the ranges 9.0 ≤ log(M_⋆/M_⊙) <  9.5 and 4.0 ≤ z <  7.0 (Fig. 14). In general, the equal and pbcor weighting schemes achieve lower rms values in each bin, but the S/N values are modestly higher in some bins with F_UV and μ-weighting, mirroring the results from stacking all sources combined. For instance, when using F_UV (μ) weighting, we find that ∼4% (21%) of the binning configurations with more than one candidate deliver better S/N values than the pbcor or equal weighting cases.

Fig. 14. Example u–v stacked image stamps for 19 undetected LBG candidates in the range of 9.0 ≤ log(M⋆/M⊙) <  9.5 and 4.0 ≤ z <  7.0. Details same as Fig. A.1. Color scale spans −125 μJy to +125 μJy range.

Given past efforts (e.g., B16), it is somewhat surprising that we do not find significant stacked signal from LBG candidates with stellar masses in excess of or close to 10¹⁰ M_⊙ (Table B.3). In part, this is a consequence of the small number of sources in the highest mass bin (only three candidates, one per each redshift bin). Moreover, for most of the configurations in this range, the stacks show relatively high noise levels (rms ≳ 115 μJy), which arise because the targets are rare and generally lie close to the border of the ALMA maps and, hence, have higher noise due to beam attenuation. For this reason, the stacking results for these bins provide only relatively weak constraints.

4.2.1. Stacked IRX–β relation

We next consider the stacking constraints on the IRX–β relation, which are presented in Fig. 15. We apply all four weighting schemes and list the full results in Table B.1. Here we split the sources into several β bins for three distinct photometric redshift ranges. For completeness, we plot the ALMA detected LBG candidates alongside the stacking results. We omit β bins which contain no LBG candidates or resulted in a negative IRX stacked value.

Fig. 15. Stacked observed infrared-excess (IRX) 3-σ limits as a function of UV-slope (β). For each β bin, the weighted median IRX upper limit is shown (orange triangles for equal weighting, black diamonds for pbcor-weighting, purple squares for magnification-weighting, and light rose circles for UV flux weighting; see Sect. 3.2). Results are separated into three photometric redshift bins all of which are upper limits. For comparison, we also show the local IRX–β relations (B16, SMC, Smit et al. 2016; Takeuchi et al. 2012; Casey et al. 2014) and the locations of three well-known local star-forming galaxies (M 82, NGC 7552 and NGC 7714). The downward arrows have 1-σ length. Horizontal errorbars indicate the 16th and 84th percentiles of the distribution of LBG candidates for each UV-slope bin. Yellow crosses show the two individual ALMA detections.

In all three redshift bins, we see that the F_UV-weighting generally produces much lower average IRX constraints than the other weighting schemes. This is perhaps no surprise, given the previously mentioned correlation between stellar mass and L_UV (or equivalently F_UV over limited redshift ranges) in Sect. 4.1.4. Indeed, the most massive and UV-luminous LBG candidates in Fig. 11 were the ones with limits closest to the local relations. In contrast, the pbcor and equal weighting schemes generate substantially weaker IRX constraints because they average in more of the high individual IRX constraints, which result from low L_UV detections and high L_IR limits. Finally, as can be seen from the μ color-coded panel of Fig. 11, the individual high-magnification LBG candidates generally have some of the highest IRX limits, which translates into high IRX limits for the μ-weighted stacked bins too. We additionally see that the lowest β bins (β <  −2) have systematically higher IRX limits, mirroring the trend seen in the individual limits of Fig. 11.

For comparison, we show the two ALMA-detected LBG candidates and three local star-forming galaxies: M 82 (Sheth et al. 2010; Förster Schreiber et al. 2003; Dale et al. 2007; Greco et al. 2012), NGC 7552 (Sheth et al. 2010; Dale et al. 2007; Wood et al. 2015) and NGC 7714 (Sheth et al. 2010; González Delgado et al. 1999; Brandl et al. 2004; Brown et al. 2014). These local galaxies have range of M_⋆ ∼ 10^10.3–10^10.7 M_⊙ and SFR ∼ 1–10 M_⊙ yr⁻¹, with M 82 being perhaps the most reasonable counterpart to the more massive LBG candidates. The detected LBG candidates generally have lower UV-slopes (less extinction), much lower stellar masses, and higher or comparable IRX values to the local objects. The limits for the F_UV-weighted limits are systematically lower than the two ALMA detections and, in the higher redshift bins, show similar or lower IRX values than the local objects despite having similar stellar masses.

The limits at z <  4 lie well above the local IRX–β relations, demonstrating that at least 1 dex deeper IR observations are needed to start placing meaningful constraints on even the most luminous z ∼ 2–4 LBGs, and 2–3 dex more for the bulk of the population. At higher redshifts, the results appear more auspicious, as the limits on the most UV-luminous objects are approaching those of the local relations. Unfortunately, the low numbers of sources in these high-redshift bins mean the results are subject to small number statistical uncertainties and, thus, we can only say that they remain marginally consistent with the local IRX–β relations at the depths we probe.

4.2.2. Stacked IRX–M_⋆ relation

Finally, we consider the IRX–M_⋆ relation, the results of which are presented in Fig. 16. Again, we apply all four weighting schemes and list the full results in Tables B.2 and B.3. Here we split the sources into several M_⋆ bins spanning three photometric redshift ranges. We omit stellar mass bins that contain no LBG candidates or resulted in a negative IRX stacked value. For comparison, we plot the ALMA-detected LBG candidates and local galaxies alongside the stacking results.

Fig. 16. Stacked observed infrared-excess (IRX) 3-σ limits as a function of stellar mass (M⋆). For each M⋆ bin, the weighted median IRX upper limit is shown (same colors and markers as in Fig. 15). Results are separated into three photometric redshift bins, all of which are upper limits. We also show a number of previously reported IRX–M⋆ relations (Consensus, F17, M18, B16 and H14) and the locations of three well-known local star-forming galaxies (M 82, NGC 7552 and NGC 7714). Downward arrows have 1-σ length. Horizontal errorbars indicate the 16th and 84th percentiles of the distribution of LBG candidates for each M⋆ bin. Yellow crosses show the two individual ALMA detections.

As with the IRX–β results, we find that the F_UV-weighting produces lower median IRX constraints compared to the other weighting schemes, although the distinction between these is, in general, less pronounced than with β. In most bins, the limits lie 1–2 dex above the consensus relations, although we do see a trend wherein the highest mass bins (e.g., ≳10⁹ M_⊙) have lower IRX limits than the lower M_⋆ bins and fall near or below the F17, M18 and consensus relations. It is important to note here, however that these limits comprise only 1–2 of the most extreme individual LBG candidates and, hence, cannot be considered representative of the full sample.

Fig. 17. UV-slope (β) vs. stellar mass (M⋆) for our selected LBG candidates. Colors represent our three photometric redshift bins (both point and polynomial fitting lines). Blue dashed line represents fit from M18 for their sample with stellar masses log(M⋆/M⊙) ≥ 8.5 and green solid line shows our third-order polynomial fit (Eq. (18)). We also include the first-order polynomial from Eq. (19) in gray.

The ALMA-detected LBG candidates have comparable IRX and stellar mass values to the stacked limits (regardless of weighting scheme) and are presumably the extreme tail of the distribution.

5. Discussion

5.1. Individual constraints

We compare the properties from our sample (UV-slope, stellar mass, UV magnitude) with the distributions from B16, where 330 LBGs were studied and six 2-σ tentative ALMA detections were obtained. In particular, B16 present histograms for these properties as a function of drop-out bins in their Fig. 2, while we present our sample distributions in Figs. 4 and 6.

The β values in both samples show a peak near β ∼ −2.2 and similar distribution shapes, with the bulk of sources located in the range −3.5 ≲ β ≲ 0.5. The stellar mass distributions both peak at around log(M_⋆/M_⊙) = 8. However, our sample effectively probes one dex lower in mass due to the magnifying power of the galaxy clusters. Finally, the (magnification-corrected) apparent UV magnitude distributions share similar maximum (m_UV ∼ 24) and peak values (28 ≲ m_UV ≲ 29.5) but our candidates probe two magnitudes deeper (m_UV ∼ 32) than in B16, again due to the galaxy cluster lensing. Notably, without the capped magnification factors, our distributions would extend to even smaller values.

5.1.1. ALMA expectations

Our ALMA observations (both detections and upper limits) can be contrasted with previous studies of mm and submm emission from LBG candidates over comparable redshift ranges. Several works have stacked multi-band IR photometry to generate empirical SEDs or fit against templates to derive average physical properties (e.g., Elbaz et al. 2011; Magdis et al. 2011a,b; Oteo et al. 2013a,b; Coppin et al. 2015; Schreiber et al. 2017; Faisst et al. 2017; Bowler et al. 2018). These works studied star-forming galaxies or LBGs with redshifts ranging from z ∼ 1 to z ∼ 7 and derived IR luminosities ranging from log(L_IR/L_⊙) = 9.9 (in Magdis et al. 2011a) to log(L_IR/L_⊙) = 12.5 (in Oteo et al. 2013b). The majority of these works focus on extreme or rare LBG candidates with masses above (log M_⋆/M_⊙ ≳ 10.0) and SFRs ≳ 100. In these cases, there is essentially no overlap in M_⋆ or SFR distributions compared to our sample, making comparisons and interpretations of detections or limits impossible. Two exceptions that share some overlap are Bowler et al. (2018) and Coppin et al. (2015), which we discuss below.

Bowler et al. (2018) reported on ALMA band 6 observations of 6 z ∼ 7 LBGs with log M_⋆/M_⊙ = 9.0–9.6 and log L_UV/L_⊙ = 11.3–11.6, selected from 1.0 deg² of UltraVISTA imaging in the COSMOS field; one confirmed LBG is marginally detected in ALMA Band 6 with 168 ± 56 μJy while the stacked limit for the remaining five candidates is 100 ± 50 μJy. We only have 1 LBG candidate which overlaps with this stellar mass at z >  7, but has a much fainter UV luminosity; our limits are consistent with the ones reported by Bowler et al. (2018).

Coppin et al. (2015) stacked 850 μm SCUBA data for 5138 z ∼ 3 − 5 LBG candidates with (log M_⋆/M_⊙ = 9.0–11.0) and log L_UV/L_⊙ = 10.0–11.6, selected from 0.62 deg² in the UKIDSS-UDS field; they reported 850 μm detections of 250 ± 29, 411 ± 64, 875 ± 229 μJy for their z ∼ 3, 4, and 5 samples, respectively (adopting an emissivity of 1.6, these equate to 1.1 mm detections of ≈165, 270, and 580 μJy). Our stacked limits for 72 (z <  4) and 25 (4 <  z <  7) similarly massive LBG candidates (log M_⋆/M_⊙ = 9.0–10.0; Table B.3) lie between 24 ± 9 and 62 ± 23 for z <  4 and 22 ± 19 to 27 ± 30 μJy for 4 <  z <  7. Even our two tentative detections at ≈285 μJy (Table 4) lie a factor of 2 below the stacked average at z ∼ 5, and have masses that are more than two orders of magnitude lower than those from Coppin et al. (2015). This strong discrepancy, by factors > 3−20, implies that a few strong sources (e.g., obscured AGN or dusty star-forming galaxies), are very likely biasing their results. Coppin et al. (2015) do note that the fitted SED models show lower radio fluxes than the measured values, possibly implying an AGN contribution in part of their sample.

Another method to test our IR expectations involves deriving IR luminosities with the use of the UV luminosities (see Sect. 2.11) and their SFR values (Sect. 2.8) separately. Several calibrations have been developed to extract a star formation rate estimate, with the one from Kennicutt (1998) and the reformulation from Bell & Kennicutt (2001) being among the most utilized. The difference between the corrected and uncorrected UV SFRs corresponds to the IR SFR values (Daddi et al. 2007a). From the derivations of Oteo et al. (2013c) and Elbaz et al. (2018), we derived the corrected and uncorrected UV SFRs using the Calzetti et al. (2000) law (see Sect. 3).

Thus, we obtained an estimate for the 3-σ IR luminosity upper limits and compared them with our values derived from ALMA observations. We find that UV-based L_IR estimates spread along a wider range (∼10^8.5–10^13.2 L_⊙) than the ALMA ones (∼10^11.1–10^13.4 L_⊙). For the vast majority of LBG candidates, the UV-based estimates lie well below the ALMA-based ones, essentially in lines with their locations on the IRX diagrams. For a handful of a few high-luminosity and high-redshift LBG candidates at high redshift, the constraints are similar; these are the same sources that lie near the local IRX relations. In summary, we find good consistency between this method and our IRX analysis.

5.1.2. Detected LBGs

Overall, we detect only two sources with S/N ≳ 4.1, both from the AS1063 field. The latter fact may imply that these detections are not representative of the LBG population as a whole. One source is highly magnified (μ = 39.0), while the other is not (μ = 2.9). The best-fit UV properties of the detected sources appear normal for LBGs, but quite atypical of known ALMA sources, with high redshifts (z ∼ 5.4–5.5), low stellar masses (≈6.6–6.9 log M_⋆/M_⊙), low UV fluxes (≈8.7–8.8 log L_UV/L_⊙), and modest UV slopes (β ≈ −1.2 to −1.6) yet relatively high IR luminosities (≈11.3–11.9 log(L_IR/L_⊙)). These properties imply that both sources may be spurious, or that the ALMA detections are coming from highly obscured regions within these galaxies that are not accounted for by the HST data. Despite the low number of individual detections, stacking provides better insights into our LBG candidates and their properties.

With so few detected sources, stacking does not offer any further insight. The non-detections in all but one of the ALMA-FFs, despite the potential for strong lensing, indicate that the intrinsic IR emission from all LBGs is faint. Our results are consistent with many past works (B16; G17; M18; H14; Casey et al. 2014; Takeuchi et al. 2012; Álvarez-Márquez et al. 2016; Barisic et al. 2017; Bowler et al. 2018; Bourne et al. 2017).

Regardless of their scarcity, we can compare the detections to the stacked IRX upper limits depicted in Figs. 15 and 16. The IR luminosities of the detections are higher than an important fraction of our sample but seem to lie within the sample distribution, consistent with expectations.

5.1.3. IRX–β upper limits and previous works

Figure 11 shows the distribution of our sample over the IRX–β plane color-coded by five different quantities. Along with this distribution, we include four different relations from past works: M99 IRX–β relation by Meurer et al. (1999, also introduced in Sect. 2.13), Small Magellanic Cloud (SMC) IRX–β relation (Smit et al. 2016), Takeuchi et al. (2012) and Casey et al. (2014) relations.

We can see that the vast majority of our candidates have limits well above the relations already mentioned. Only three (3) points are located below the M99 IRX relation. Thus, most of our upper limits remain compatible with all of the studied IRX–β relations. It is pertinent to mention that the candidates situated between the M99 and SMC relations have high photometric redshifts; they lie in the interval z_ph = 5.69–8.12 with a mean value of z_ph = 6.77. The UV luminosities from these candidates all skew toward the high end of the distribution, while their IR luminosities skew toward the low end. This combination leads to lower IRX values, pushing them below several published IRX–β relations.

Given the known dispersion in the local relations, we would have naively expected to detect at least a few among the handful of sources that lie between the M99 and SMC relations. The fact that we have no detections at least hints at the possibility that high-redshift relations may have systematically lower IRX values.

Targets that lie below β = −2.23, which represents the intrinsic, non-dust-obscured, UV-slope value from Meurer et al. (1999), cannot be compared directly with the mentioned relations as they do not cover the same region of the parameter space. Extremely low β sources must be compared with previous relations, as mentioned in Sect. 2.13, developed specifically with such galaxies in mind. Leaving these issues aside, it is clear that with current instrumentation, we are unable to probe to sufficiently low IR luminosities yet to study the behavior of typical LBGs with extremely low β values and, thus, extend known relations. Furthermore, having obtained only upper limits makes it impossible to search for a meaningful correlation between IRX and either β or M_⋆.

Most of the recent works mentioned here agree on the fact that star-forming galaxies up to z ∼ 3 − 4 follow the B16 IRX–M_⋆ relation more closely than an SMC-like curve. Our IRX–M_⋆ upper limits (both individual and F_UV-weighted stacks) suggest that the most luminous and highest redshift sources (β ≳ −1.0 and z ≳ 4) are pushing below the B16 relation and may be more compatible with an SMC-like curve (as suggested by, e.g., Koprowski et al. 2018, F17).

Regarding a possible evolution of the IRX relations with redshift as mentioned by, for example, F17, we are unable to establish this given that we only find IRX upper limits in the stacked emission. We can only observe (uppermost panels in Figs. 11 and 12) that, roughly, upper limits with higher redshifts tend to exhibit lower IRX ratios, and the lack of any individual detections hints at some evolution, but the statistics are currently too limited to say more.

Finally, we compare our upper limits in the IRX–β space with the recent values presented by Salim & Boquien (2019), who examined more than 20 000 low-redshift galaxies (z <  0.3) from GALEX-SDSS-WISE Legacy Catalog 2 (GSWLC-D2). Notably, they find that a majority of their sources in the range −2.0 ≲ β ≲ −0.5 lie below the B16 relation, while nearly all β ≳ −0.5 sources lie well below it, indicating a less abrupt slope comparable to the other relations presented in Sect. 2.13. We see that our best 10 upper limits are currently consistent with the IRXs of these possible low redshift analogs, although given that our constraints are only limits, there remains the potential for some mild tension with even the SMC-like relation, which appears to act as a lower bound on the local Salim & Boquien (2019) sample.

5.1.4. IRX–M_⋆ upper limits and previous works

Individual IRX and stellar mass values are shown in Fig. 12, along with five relations from previous works (see Sect. 2.13 for details; F17; M18; B16; H14). These relations are extrapolated down to stellar masses of log(M_⋆/M_⊙) = 6.0 where applicable, to match our lowest stacking bins. For the same reasons exposed in Sect. 5.1.3, we do not fit any relations to our upper limits.

All our upper limits lie above the curves. The LBGs that fall very close to the higher IRX–M_⋆ relations correspond to the candidates with the highest stellar masses. Similar to the observed behavior with the UV-slope, lower stellar mass candidates do not probe to as low IRX values as their higher mass counterparts.

While our candidates share similar distributions in several important properties with B16, the resulting IRX limits for the majority of candidates remain ≈0.3 dex higher due to our shallower ALMA maps (factor of ≈4–5 higher rms) and the low average magnifications (⟨μ⟩≈4.3) of our candidates. For the rare high-magnification targets, predominantly low-mass LBG candidates (as seen in panel 2 of Fig. 12), our maps provide modestly deeper IRX constraints.

With respect to the lower part of the IRX distribution, our z <  4.0 sample lies ∼1.0 dex above the respective bins from B16.

Comparing our results with those of F17, we find that our IRX upper limits, for comparable β and M_⋆ ranges, are ∼0.5 dex higher than their sources. Since both ALMA observations reach similar noise levels, the discrepancy must lie in the fact that their sample is composed of sources up to much higher stellar masses (M_⋆ ∼ 10^10.7 M_⊙) and much higher L_UV values (i.e., their most stringent constraints arise from the sample of Capak et al. 2015, who targeted the brightest LBGs over the much larger 2 deg² COSMOS field).

5.1.5. M_⋆–β correlation

As mentioned in Sect. 4.1.4, we see a fairly clear trend between M_⋆ and β. To place this in better context, we plot in Fig. 17 the relation between UV-slope and stellar mass directly for our sample. Under the assumption that all star-forming galaxies have similar intrinsic UV slopes, M18 used the values of β as a proxy for the UV attenuation (A₁₆₀₀). They fit a polynomial to a mass-complete sample of star-forming galaxies selected from the Hubble Ultra Deep Field (HUDF; Beckwith et al. 2006; Dunlop et al. 2017) and obtained the third-order relation plotted in the dashed blue line in Fig. 17 (for stellar masses in the range 8.5 <  log(M_⋆/M_⊙) <  11.5, with a 1.1 mm rms of 35 μJy beam⁻¹).

We performed the same experiment using our full LBG sample down to a mass of log(M_⋆/M_⊙) = 6.0. We caution that this limit is likely substantially below the nominal mass completeness threshold in the Hubble Frontier Fields (HFFs), which should be similar to that of the HUDF at z ∼ 3 (e.g., ≳10^8.5 M_⊙; see Fig. 6). Due to the lensing amplification, we do expect to find at least some representative sources among the lower mass LBGs in our sample, but we could have strong selection effects that bias the resulting fitted relations at low stellar masses.

Nonetheless, applying a third-order polynomial fit to our M_⋆ and β values, we find

$$\begin{array}{c}\hfill \beta =-0.99+0.62X+0.13X^2+0.02X^3,\end{array}$$(18)

in which X = log(M_⋆/10¹⁰ M_⊙). Similar exercises were performed binning the sample in our adopted redshift bins. The z <  4.0 trend is nearly identical to the full sample, due in large part to the fact that such low-redshift sources account for the majority of our sample (1068 candidates). However, the trends found for the higher redshift bins remain consistent within the dispersion. Within the mass-complete range of 8.5 ≲ log(M_⋆/M_⊙) ≲ 10, our fits appear to be consistent with that of M18, particularly in the low-redshift bin (within 0.25 dex), which is most comparable to the range they studied.

We note that fitting a simpler first-degree polynomial (power law) to our full sample (M_⋆ − β space) yields the following line:

$$\begin{array}{c}\hfill \beta =-1.08+0.40X,\end{array}$$(19)

with X = log(M_⋆/10¹⁰ M_⊙). This line, shown in gray in Fig. 17, is nearly identical to the third-order polynominal fit above, demonstrating that there is no unexpected oscillatory behavior in the former curve and, thus, corroborating the trend seen in M18.

Pushing below stellar masses of ∼10^8.5 M_⊙, we observe a smooth trend toward lower (bluer) β values, consistent with expectations from increasingly metal-poor stellar populations. Some caution must be exercised nonetheless since both properties, stellar mass and UV-slope, have been derived from the same data (HST photometry) and they are not, consequently, completely independent.

5.2. The role of dust temperature

To understand the effects of our evolving dust temperature prescription on our results, we compared against a model assuming that all LBG candidates have a constant temperature of T_dust = 35 K (following, for example, Kovács et al. 2006; Coppin et al. 2008). Such low dust temperatures predominantly affect the properties of high-redshift (z >  4) candidates, resulting in IR luminosity and IRX values drops of ∼0.5 dex, such that more candidate upper limits (both individually and stacked) are pushed below the B16 IRX–β relation (with some even approaching the SMC relation) and a few candidates fall below the IRX–M_⋆ consensus relation.

5.3. ALMA-FF LBG sample overview

To place our LBG sample in context with other samples, we estimated in very crude terms the source density of LBGs per angular area in the source plane. We shifted to the source plane because the high magnification by massive clusters strongly affects the number density of observed background sources. This is non-trivial, however, since the exact magnification depends on the redshifts of the sources.

We estimated the source-plane area as follows. Following the methodology presented in Sect. 2.5, we calculated the source-plane beam area for each candidate LBG as the area of the synthesized beam for the observed ALMA map centered on the LBG position, divided by the adopted magnification factor. We then summed the individual LBG source and image-plane areas and divided the totals to obtain a “demagnification” factor, which in our case is 0.35. Finally, we multiplied the total (image-plane) area of the ALMA FFs observations by this ratio to estimate crudely the total source-plane area covered by our ALMA observations. Each FF cluster was observed over a $\approx 2\stackrel{\prime }{.}1\times 2\stackrel{\prime }{.}2$ ALMA mosaic (González-López et al. 2017b), summing up to an image-plane area of ∼23 arcmin² used in this work. Applying the factor of 0.35, we should have, on average, an effective source-plane region of $\approx 1\stackrel{\prime }{.}24\times 1\stackrel{\prime }{.}3$ per cluster, and a total source-plane area of ∼8 arcmin² over the five FFs clusters.

With this, we can obtain an estimate for the intrinsic density of LBGs (regardless of their redshift) per unit area. A simple ratio of our 1582 studied sources over the effective area covered by ALMA gives a value of ≈200 LBGs per arcmin². To establish a reference with other ALMA observations, B16 studied 330 LBGs over a 1 arcmin² region of HUDF. And F17 examined 67 star-forming galaxies in an area of 39″ × 39″, which corresponds to a density of ≈160 sources per arcmin². Our survey appears to be intermediate between these two.

6. Summary and conclusions

In this paper, we utilize ALMA 1.1 mm mosaic images for five of the six FFs clusters, with rms values between ≈55 and 71μJy, to place constraints on the IR excesses of 1582 UV-selected LBGs as functions of their UV-slopes (β), stellar masses (M_⋆), sSFRs and photometric redshifts. After correcting for magnification, the source plane area of the five clusters is ∼8 arcmin², probing LBG candidates with rest-frame UV magnitudes ranging from ∼23 to 32 ABmag. We summarize our results as follows:

1. The rms levels in the ALMA maps, coupled with the likely faint intrinsic fluxes of the LBG candidates, result in very few outright detections. With a detection threshold of 4.1-σ (equivalent to a 15% false detection rate), only two LBG candidates are considered detected (both located in AS1063). The rest are treated as upper limits. Comparing our 1580 IRX 3-σ upper limits with previous IRX relations (IRX–β and IRX–M_⋆), the vast majority lie above the local and Consensus relations; only 3 LBGs are constrained to lie below the M99 IRX–β relation and none of them below any of the discussed IRX–M_⋆ relations. Our lowest IRX limits appear consistent with the known dispersion around these relations.

2. We divided the 1580 undetected LBG candidates into bins of stellar mass (6.0 ≤ log(M_⋆/M_⊙) ≤ 11), UV-slope (−4.0 ≤ β ≤ 1.5) and photometric redshift (< 4.0, 4.0 − 7.0, > 7.0), and stacked their ALMA data using the STACKER software. We implemented four weighting schemes for the uv-stacking: equal weighting; pbcor-weighting only; pbcor and UV-flux F_UV-weighting; and pbcor and magnification μ-weighting. With these configurations, we stacked the ALMA observations and computed the stacked IRX values, obtaining upper limits in all but one bin (i.e., S/N ≲ 3.5). This single detected bin yields a S/N ≈ 4.24, although this was only obtained for one weighting scheme and CLEANing configuration; as such we do not consider it to be a robust result. The bulk of our stacked IRX values remain above those cited in most of the literature for IRX–β and IRX–M_⋆ relations, although a few bins (β >  −1.0 with z_ph ≥ 7.0 and M_⋆ >  10^9.5 M_⊙ with z_ph ≥ 7.0) push below the B16 IRX–β and Consensus relations. Since these limits only represent a small portion of the overall LBG population and there is substantial known dispersion in the relations themselves, we can only say at present that these high-redshift LBGs appear consistent with the relations. That being said, the lack of any detections in the vicinity of these relations hints a possible evolution of these relations, which could be linked to an evolution in metallicity.

3. We also investigated the correlation between β and stellar mass for our candidates. Despite the significant dispersion, as well as growing incompleteness below stellar masses of ∼10^8.5 M_⊙, we observe a clear and smooth trend that extends to lower masses and bluer (lower) β values, consistent with expectations from previous works regarding metallicity-driven evolution.

To improve upon our results would require either the stacking of substantially larger LBG samples or ALMA observations with at least a factor of several and ideally > 1 dex lower average rms. Since covering larger areas or reducing noise levels to this extent would be costly in terms of ALMA observation time (e.g., Franco et al. 2018 achieve ≈180 μJy rms in the ∼69 arcmin² GOODS-S mosaic at 1.1 mm in 18.5 h; Gonzalez-Lopez et al., in prep., achieve ≈14 μJy rms in a ∼5 arcmin² HUDF mosaic at 1.1 mm in 85 h), a better strategy may be to target the most highly magnified LBG candidates with deep, single-pointing observations, as opposed to the mosaic observations used for this work. In this way, significantly fewer pointings would be required but for longer durations.

Acknowledgments

We acknowledge support from CONICYT grants Basal-CATA AFB-170002 (RC, FEB, ET, RD, LG), FONDECYT Regular 1141218 (RC, FEB) 1160999 (ET) and 1190818 (ET), Programa de Astronomia FONDO ALMA 2016 31160033 (LG), and the Ministry of Economy, Development, and Tourism’s Millennium Science Initiative through grant IC120009, awarded to The Millennium Institute of Astrophysics, MAS (FEB). This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00999.S, ADS/JAO.ALMA#2015.1.01425.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. Part of our work is based on observations obtained with the NASA/ESA Hubble Space Telescope, retrieved from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute (STScI). STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555. This work utilizes gravitational lensing models produced by PIs Bradač, Natarajan & Kneib (CATS), Merten & Zitrin, Sharon, Williams, Keeton, Bernstein and Diego, and the GLAFIC group. This lens modeling was partially funded by the HST Frontier Fields program conducted by STScI. STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555. The lens models were obtained from the Mikulski Archive for Space Telescopes (MAST). This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013; Price-Whelan et al. 2018), matplotlib, a Python library for publication quality graphics (Hunter 2007), APLpy, an open-source plotting package for Python (Robitaille & Bressert 2012), and NASA’s Astrophysics Data System.

References

Appendix A: u–v stacking tests

As a cross-check to confirm that STACKER behaves as expected, we first performed image and u–v stacking for the twelve (12) robustly detected dusty star-forming galaxies from A2744, MACSJ0416, and MACSJ1149, which are described in González-López et al. (2017b) and in Laporte et al. (2017b); the individual detections range from 5.1-σ to 25.9-σ. These stacking results are presented in Figs. A.1 (u–v stacking) and A.2 (image stacking).

Fig. A.1. u–v stacked image stamps for the 12 detected sources from González-López et al. (2017b). Panels denote specific weighting configurations (from left to right: equal, pbcor, FUV, and μ) and CASACLEANing procedures (upper: Natural; lower: Taper). Color scale spans −455 μJy to +455 μJy range and contours are drawn for every 5 × rms level. White and black ellipses represent the synthesized beam size, while white bars in the right corner denote 5″ scale. The number of sources used for the stacked bin is denoted at bottom, as well as the resultant rms. (a) Weight: equal, (b) weight: pbcor, (c) weight: FUV, (d) weight: μ.

Fig. A.2. Image stacked stamps for the 12 detected sources from González-López et al. (2017b). Details same as Fig. A.1.

The u–v stacking (Fig. A.1), adopting natural-weight CLEANing, yields stacked detections of 39.1-σ, 37.3-σ, 29.9-σ, and 7.1-σ, for equal, pbcor, F_UV, and μ-weighting, respectively. Unsurprisingly the unweighted and pbcor-weighted stacks achieve much better stacked S/N and rms in this scenario because all of the sources have high S/N and lie within the high pbcor regions. The strong variation between the μ and F_UV weighted stacks arises because of the dominating presence of the UV-bright, low-magnification source M1149-ID01.

As a second test, we stacked all 1569 sources (upper limits with stellar masses higher than 10⁶ M_⋆) into a single image, adopting the four weighting schemes (equal, pbcor, F_UV, and μ), obtaining rms levels down to ∼2 μJy. S/N values are S/N_no = 1.5, S/N_pbcor = 1.5, S/N_FUV = 3.1, and S/N_μ = 2.0 for equal, pbcor, F_UV, and μ-weighting, respectively. In this case, the μ-weighting achieves only slightly higher rms values than the equal and pbcor-weighting, despite weighting far fewer sources. On the other hand, the F_UV-weighting achieves a much worse rms yet a relatively high S/N since it is optimized for candidates with higher star formation rates. Stamps of the stacked images are shown in Fig. A.3.

Fig. A.3. u–v stacked stamps for the 1569 upper limits with stellar masses in excess of 106 M⋆. Details same as Fig. A.1. Color scale spans −30 μJy to +30 μJy range.

Appendix B: Stacking results

We include here the results of the stacking process for the LBG candidates as functions of UV-slope β and stellar mass M_⋆. For each configuration we list four results pairing the different stacking weights and CLEANing methods that were used.

B.1. UV slope binning

LBG candidates were divided into five UV-slope (β) bins as described in Sect. 3. Targets with stellar masses below 10⁶ M_⊙ were excluded to avoid poorly constrained sources with low-quality photometry and, hence, maintain consistency with stellar mass stacking. We list results in Table B.1 for both magnification (μ) and UV flux (F_UV) stacking weights, as well as natural and taper CLEANing methods (see Sect. 3).

B.2. Stellar mass binning

The values obtained after stacking each different configuration are shown in Tables B.2 and B.3.

B.2.1. Candidates with log (M_⋆ / M_⊙) < 6.0

Among the 1580 non-detected LBG candidates in our initial sample, 11 have stellar masses below 10⁶ M_⊙, which were not considered for the results presented in Sect. 4.2. Regarding photometric redshifts, eight (8) low-mass candidates lie at z_ph <  4.0, three (3) between 4 ≤ z_ph <  7.0 and none at z_ph ≥ 7. With so few candidates, and considering their relatively low UV luminosities, we do not expect to detect anything then via stacking. Nonetheless, for completeness, we report results of their uv-stacking in Table B.4. The highest absolute S/N obtained is 1.18 for low-redshift candidates with taper CLEANing, confirming our expectation.

Appendix C: Individual properties

The individual properties of the first 10 of the 1582 selected LBG candidates are presented in Tables C.1 (HST photometry), C.2 (ALMA properties), C.3 (derived properties), C.4 (FAST++ derived properties) and C.5 (observed luminosities and flux densities). Properties for the remaining targets are available at the CDS.

All Tables

All Figures

	Fig. 1. Photometric redshift (zph) values in our sample (see Sect. 3.1) calculated from FAST++. Colors represent each drop-out band.
In the text

	Fig. 2. Distribution of magnification factors (μ) for the three photometric redshift bin samples.
In the text

	Fig. 3. Capped distribution of magnification factors (μ) for the three photometric redshift bin samples. Compare to the full distributions shown in Fig. 2.
In the text

	Fig. 4. Comparison of β and magnification-corrected magnitudes for selected sources. Top and right panels: histograms of UV slope and magnitude distributions. Colors represent photometric redshift subsamples, as described in the legend and in Sect. 3.
In the text

Fig. 5. Distribution of β values according to photometric redshift. Heights of boxes represent the 25% and 75% quartiles of the data. Horizontal lines inside the box indicate the median value for each redshift bin. Vertical error bars span the central 2-σ of the data. Numbers above the median in each box state the number of LBG candidates assigned to each bin. Even though there is not a zph ∼ 4 band from drop-out selection, there are candidates in that bin.

In the text

	Fig. 6. Stellar masses in our sample. Sample has been divided according to the photometric redshift bins defined in Sect. 3. Vertical dark line represents the approximate completeness limit from M18 (See Sect. 5.1.5).
In the text

	Fig. 7. False detection rate, pFDR, for the five ALMA maps. Vertical dashed lines denote the highest detected S/N among the LBG candidates in each ALMA map. The vertical solid line denotes our adopted S/N cutoff of 4.1, which equates to a FDR around 15% among the cluster fields.
In the text

Fig. 8. Mean (left) and peak (right) signal-to-noise ratios (S/N) for our candidates in the ALMA maps. The LBG candidates are separated into three photometric redshift sub-samples, represented by distinct colors. The mean value is centered around S/N ∼ 0 (vertical dark line) and is roughly Gaussian. The peak values are centered around S/N ∼ 1, rather than S/N ∼ 0 (vertical dark line), due to the selection of the peak pixel which arises within half a beamwidth; this conservatively biases the maximum flux associated with a candidate to higher values.

In the text

	Fig. 9. UV luminosities in our sample. The LBG candidates are separated into three photometric redshift sub-samples, represented by distinct colors.
In the text

	Fig. 10. IR luminosities (3-σ upper limits from our ALMA maps) in our sample. The LBG candidates are separated into three photometric redshift sub-samples, represented by distinct colors.
In the text

Fig. 11. Comparison of infrared excess (IRX) 3-σ upper limits and UV-slopes (β) for our LBG candidates. Downward arrows have 1-σ length. From top to bottom panels: colors represent: photometric redshift (zph), magnification factor (μ), star formation rate (SFR), stellar mass (M⋆) and UV Luminosity (LUV). Local IRX–β relations presented in Sect. 2.13 are shown for reference. Blue crosses represent the two detections.

In the text

Fig. 12. Comparison of infrared excess (IRX) 3-σ upper limits and stellar masses M⋆ for our LBG candidates. Downward arrows have 1-σ length. From top to bottom panels: colors represent: photometric redshift (zph), magnification factor (μ), star formation rate (SFR), UV slope (β) and UV Luminosity (LUV). Local IRX–M⋆ relations presented in Sect. 2.13 are shown for reference. Blue crosses represent the two detections.

In the text

	Fig. 13. Example u–v stacked image stamps for 227 undetected LBG candidates in the range of −2.0 ≤ β <  −1.0 and 4.0 ≤ z <  7.0. Details same as Fig. A.1. Color scale spans −125 μJy to +125 μJy range.
In the text

	Fig. 14. Example u–v stacked image stamps for 19 undetected LBG candidates in the range of 9.0 ≤ log(M⋆/M⊙) <  9.5 and 4.0 ≤ z <  7.0. Details same as Fig. A.1. Color scale spans −125 μJy to +125 μJy range.
In the text

Fig. 15. Stacked observed infrared-excess (IRX) 3-σ limits as a function of UV-slope (β). For each β bin, the weighted median IRX upper limit is shown (orange triangles for equal weighting, black diamonds for pbcor-weighting, purple squares for magnification-weighting, and light rose circles for UV flux weighting; see Sect. 3.2). Results are separated into three photometric redshift bins all of which are upper limits. For comparison, we also show the local IRX–β relations (B16, SMC, Smit et al. 2016; Takeuchi et al. 2012; Casey et al. 2014) and the locations of three well-known local star-forming galaxies (M 82, NGC 7552 and NGC 7714). The downward arrows have 1-σ length. Horizontal errorbars indicate the 16th and 84th percentiles of the distribution of LBG candidates for each UV-slope bin. Yellow crosses show the two individual ALMA detections.

In the text

Fig. 16. Stacked observed infrared-excess (IRX) 3-σ limits as a function of stellar mass (M⋆). For each M⋆ bin, the weighted median IRX upper limit is shown (same colors and markers as in Fig. 15). Results are separated into three photometric redshift bins, all of which are upper limits. We also show a number of previously reported IRX–M⋆ relations (Consensus, F17, M18, B16 and H14) and the locations of three well-known local star-forming galaxies (M 82, NGC 7552 and NGC 7714). Downward arrows have 1-σ length. Horizontal errorbars indicate the 16th and 84th percentiles of the distribution of LBG candidates for each M⋆ bin. Yellow crosses show the two individual ALMA detections.

In the text

Fig. 17. UV-slope (β) vs. stellar mass (M⋆) for our selected LBG candidates. Colors represent our three photometric redshift bins (both point and polynomial fitting lines). Blue dashed line represents fit from M18 for their sample with stellar masses log(M⋆/M⊙) ≥ 8.5 and green solid line shows our third-order polynomial fit (Eq. (18)). We also include the first-order polynomial from Eq. (19) in gray.

In the text

Fig. A.1. u–v stacked image stamps for the 12 detected sources from González-López et al. (2017b). Panels denote specific weighting configurations (from left to right: equal, pbcor, FUV, and μ) and CASACLEANing procedures (upper: Natural; lower: Taper). Color scale spans −455 μJy to +455 μJy range and contours are drawn for every 5 × rms level. White and black ellipses represent the synthesized beam size, while white bars in the right corner denote 5″ scale. The number of sources used for the stacked bin is denoted at bottom, as well as the resultant rms. (a) Weight: equal, (b) weight: pbcor, (c) weight: FUV, (d) weight: μ.

In the text

	Fig. A.2. Image stacked stamps for the 12 detected sources from González-López et al. (2017b). Details same as Fig. A.1.
In the text

	Fig. A.3. u–v stacked stamps for the 1569 upper limits with stellar masses in excess of 106 M⋆. Details same as Fig. A.1. Color scale spans −30 μJy to +30 μJy range.
In the text