ALMA Lensing Cluster Survey: ALMA-Herschel Joint Study of Lensed Dusty Star-forming Galaxies across z ≃ 0.5 – 6

All of the sources in this work are selected with the ALCS, which will be detailed in S. Fujimoto et al. (2022, in preparation). ALMA Band 6 observations for the 33 clusters were conducted through Program 2018.1.00035.L (ALCS; PI: Kohno), and we also combined archival data from Programs 2013.1.00999.S and 2015.1.01425.S (ALMA Frontier Fields; PI: Bauer; González-López et al. 2017). The list of the ALCS clusters with their coordinates, short names (e.g., M0553 for MACS J0553.4–3342), and HST program names is presented in Table 1. The observations were obtained at a central wavelength of 1.15 mm with a 15 GHz total bandwidth (i.e., two tunings of dual polarization: 250.0–257.5 GHz and 265.0–272.5 GHz). The use of two tunings instead of one allows us to search for line-emitting galaxies over a larger volume, which is another important science goal of ALCS. This, for example, led to the serendipitous discovery of a [C ii] emitter at z = 6.072 (Fujimoto et al. 2021). All the ALMA data were reduced with casa (McMullin et al. 2007), with different pipelines versions for observations obtained in different cycles (e.g., v5.4.0 for 26 clusters observed in Cycle 6 and v5.6.1 for the remaining clusters in Cycle 7). Natural-weighting continuum imaging was performed at both the native (FWHM ∼ 1'') and uv-tapered (∼2'') resolutions with the casa tclean algorithm.

Through a peak pixel identification routine of SExtractor (Bertin & Arnouts 1996) with the ALMA maps at both the native and 2''-tapered resolutions (before primary-beam correction), we securely detected 141 sources that are at either (i) signal-to-noise ratio S/N_nat ≥ 5 in the native-resolution maps or (ii) S/N_tap ≥ 4.5 in the 2''-tapered maps, over an area of ∼134 arcmin² with primary-beam response greater than 0.3 (S. Fujimoto et al. 2022, in preparation). Based on the number count of negative peaks, the number of spurious sources above these S/N cuts is expected to be around 1. We further refer to these 141 secure ALCS sources as the main sample.

A total of 258 sources were tentatively detected at (i) S/N_nat = 4 − 5 in the native-resolution maps and (ii) S/N_tap < 4.5 in the 2''-tapered maps, down to a minimal flux density of ∼0.2 mJy at 1.15 mm. Based on HST and Spitzer/IRAC images, we identified 39 of these sources with NIR/MIR counterparts within a separation of 1''. Given the high source densities in cluster fields (∼0.06 arcsec⁻² in the HST/F160W band; Sun et al. 2021a), we expect 7 ± 3 pairs of random associations among these tentative ALCS sources and cross-matched NIR/MIR counterparts. Assuming that the majority of these HST/Spitzer-matched sources are real, we note that the number found is consistent with the number difference between positive and negative peaks in the ALMA maps in this S/N range (37 in total). We further refer to these 39 tentative ALCS sources as the secondary sample, but we warn that 18% ± 8% of these sources are likely to be false matched.

We note that the continuum S/N of ALCS sources could be boosted by serendipitous emission-line detections, for example, CO (5–4) line at z = 1.11 − 1.17 and 1.24 − 1.31 (e.g., M0553-ID133/190/249 at z = 1.142; Ebeling et al. 2017; Sun et al. 2021b). All the line emitters will be reported by a future paper of the collaboration. However, because of the large bandwidth (15 GHz), the CO line contamination to continuum S/N and flux density is limited to ≲1%–10%. In addition, ALCS can only sample faint high-J CO lines (upper J number at ≥ 7) for sources at z ≳ 2. According to the CO spectral line energy distribution (SLED) of high-redshift SMGs reported in the literature (e.g., Greve et al. 2014; Béthermin et al. 2016; Yang et al. 2017; Birkin et al. 2021), even the high-J CO SLED is as flat as those reported for local active-galactic-nucleus-host (AGN-host) galaxies (e.g., Rosenberg et al. 2015), and the CO contamination will be ∼1%–10% at most. Only one [C ii] emitter was found among all continuum sources (Fujimoto et al. 2021; Laporte et al. 2021), and the continuum flux density of this source was measured on a line-subtracted continuum image.

The Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) on Herschel enabled simultaneous observations at 160 μm (red channel) with the long-wavelength camera and at either 70 or 100 μm (blue/green channel) with the short-wavelength camera. Eighteen out of the 33 clusters were imaged with PACS at both 100 and 160 μm, and two clusters were also observed with PACS at 70 μm (M1149 and AS1063). The analysis of PACS 70 μm data in these cluster fields has been presented by R16, and the only two matched sources are AS1063-ID17 (z = 1.44, f₇₀ = 7.3 ± 0.9 mJy) and AS1063-ID147 (z = 0.610, f₇₀ = 28.8 ± 2.3 mJy; analyzed in detail by Walth et al. 2019).

A total of 16 out of the 18 clusters were observed by PACS as part of the HLS (Egami et al. 2010; Sun et al. 2021b), which combines an Open-Time Key Program (KPOT; program ID: KPOT_eegami_1; nine clusters) and an Open-Time Cycle 2 (OT2; program ID: OT2_eegami_5; seven clusters) Program (both PI: Egami). The remaining two clusters, namely, A370 and RX J1347–1145, were observed as part of the PACS Evolutionary Probe (PEP; program ID: KPGT_dlutz_1, PI: Lutz; Lutz et al. 2011). All of the PACS 100 and 160 μm observations consist of two orthogonal scan maps, each comprising 18–22 repetitions of 13 parallel 4' scan legs. The summary of the PACS observations, including the observation IDs and total scan time for each cluster, is presented in Table 1.

We followed the same data reduction procedure as detailed in R16 for the HFF clusters. The PACS images were generated with UniMap (Piazzo et al. 2015) with a pixel scale of 10 at 100 μm and 20 at 160 μm. The final PACS image products have a typical field of view (FOV) with a radius of ∼4', covering the full ALMA footprints of the 18 clusters. The typical beam sizes are 74 and 114 at 100 and 160 μm, and the depths of the PACS data at the cluster center are presented in Table 2.

The Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al. 2010) on Herschel worked simultaneously at 250, 350, and 500 μm. All 33 clusters were imaged with SPIRE in two observing modes with different depths. The 18 clusters also observed with PACS were scanned with SPIRE in the "deep" mode down to confusion-limited depths (RMS_map ∼ 6 mJy beam⁻¹ at 250 μm; measured as the noise of sky background after sources being filtered out), and the remaining 15 clusters were observed in the "snapshot" mode with a shorter scan duration and thus at shallower depths (RMS_map ∼ 11 mJy beam⁻¹; as visualized in Figure 1).

Zoom In Zoom Out Reset image size

Among the total 18 clusters in the "deep" mode, 16 of them were observed as part of the HLS. Observations of the nine clusters through KPOT_eegami_1 consisted of 20 repetitions in the large scan map mode, each with two 4' scans and cross-scans (total scan time as t_scan ∼ 1.7 hr per cluster). The other seven clusters observed through Open-Time Cycle 1/2 Programs (OT1_eegami_4 and OT2_eegami_5; both PI: Egami) were imaged through 11-repetition small scan maps (1 in OT1 and 10 in OT2), and each repetition consisted of one scan and one cross-scan of 4' length (t_scan ∼ 0.4 hr per cluster). The remaining two clusters, A370 and R1347, were observed as part of the Herschel Multi-tiered Extragalactic Survey (HerMES; program ID: KPGT_soliver_1; PI: Oliver; Oliver et al. 2012). Both of these clusters were observed with eight small scan maps (six repetitions per each covering the cluster core) and three large scan maps (one repetition per each with 38' length covering a wider area), and the total scan time is 3.5 hr per cluster. The final SPIRE map products of these two clusters have a wider FOV, but the central depths are comparable to those of the HLS data.

All of the 15 clusters in the "snapshot" mode were observed as part of the HLS through single-repetition small scan maps, and each repetition consisted of one scan and one cross-scan of a 4' length (${t}_{\mathrm{scan}} \tilde 3\,\mathrm{minutes}$ per cluster). Fourteen of them were observed by the OT1 program OT1_eegami_4, and the remaining one was observed by the OT2 program OT2_eegami_6 (PI: Egami).

Table 1 summarizes the IDs and total scan times of all the SPIRE observations. All of the SPIRE data were processed by the standard reduction pipeline in HIPE v12.2 (Ott 2010), which is also detailed in R16. The output pixel sizes of the final image products are 6'', 9'', and 12'' at 250, 350, and 500 μm. The typical radii of the SPIRE FOVs are 7' for the 15 clusters observed in the "snapshot" mode, 8' for the seven clusters observed in OT1/2, 11' for the nine clusters observed in KPOT, and 33' for A370 and R1347. The full survey area of the ALCS was covered by these SPIRE images. The typical beam sizes are 18'', 24'', and 35'' in these three bands, and the depths of the SPIRE data at the cluster center are presented in Table 2.

For the purpose of enhancing the astrometric accuracy of Herschel data, we include the Spitzer/IRAC data of these 33 cluster fields in our analysis obtained from the NASA/IPAC Infrared Science Archive (IRSA ³⁷ ). We also include the HST data of all the cluster fields but only for comparing the positions of dust continuum sources with the stellar components. We defer the study of optical/NIR counterparts and panchromatic SED modeling of ALCS sources to another paper from the collaboration.

To supply accurate redshifts for FIR SED modeling (Section 4), we cross-matched the ALCS source sample with the spectroscopic redshift (z_spec) catalogs made available by the CLASH-VLT spectroscopic survey (Biviano et al. 2013), Grism Lens-Amplified Survey from Space (GLASS; Schmidt et al. 2014; Treu et al. 2015; Wang et al. 2015) and recent Very Large Telescope/MUSE surveys of massive cluster fields by Caminha et al. (2019) and Richard et al. (2021). A maximum separation of 15 is allowed for cross-matching, which is comparable to the FWHM of the IRAC point-spread function (PSF). We also include redshifts for a few sources reported by various studies in the literature (e.g., M0553 triply lensed system at z = 1.14; Ebeling et al. 2017) or private communication (e.g., M0417-ID46/58/121, an HST H-faint triply lensed system at z = 3.65; K. Kohno et al. 2022, in preparation). In addition, we also include two ALMA-HFF sources reported by Laporte et al. (2017), with their z_spec values derived from the GLASS detection of the 4000 Å break, and a triply lensed ALCS source system that belongs to a MUSE-confirmed galaxy group at z = 4.32 (Caputi et al. 2021). Spectroscopic redshifts are available for 60 ALCS sources in both the main and secondary samples.

We also utilized the HST photometric redshift (z_phot) catalogs of optical/NIR sources tabulated by CLASH (Molino et al. 2017), HFF (Shipley et al. 2018), and RELICS (Coe et al. 2019) groups. Sources are cross-matched by their coordinates, and a maximum separation of 15 is allowed. Fujimoto et al. (2016) reported a median offset of 025 between the HST and ALMA centroids of ALMA sources, and such an observed offset could be larger in cluster fields because of the lensing magnification. We also apply visual inspections of the HST F814W, F105W, and F160W images to remove any conspicuous mismatch. We identified cataloged HST z_phot measurements for 125 ALCS sources in both the main and secondary samples, including 49 sources with additional spectroscopic redshifts.

3.1.1. Image Alignment

3.1.2. Input ALCS Source Catalog

3.1.3. Background Subtraction

3.1.4. Neighborhood Examination

Herschel source extraction was performed with an iterative PSF photometry approach in an increasing order of wavelength using photutils (Bradley et al. 2019). We adopted the PSF models of PACS and SPIRE from the Herschel Science Archive. ³⁸ We also applied the spacecraft orientation angle to calculate the realistic PSF for the data taken in each cluster field.

3.2.1. Initial Guess of Flux Densities

3.2.2. Iterations of PSF Photometry

The PSF photometry was performed in three rounds of iterations in decreasing order of the significance of detection. Because most of the ALCS sources are compact in spatial extent (FWHM ≲ 1''), we assumed a point-like profile for all of the sources to be extracted. In each iteration, we only kept the results of those sources with positive extracted flux densities. The uncertainties of extracted flux densities were estimated from the covariance matrix of least-squares fitting.

In the first iteration, we tentatively extract sources that were (i) ALCS sources at S/N_ALMA ≥ 10 or (ii) ALMA-faint Herschel sources described above. These two types of objects should be the brightest sources seen in a given Herschel map. Therefore, an accurate flux density modeling of these sources will provide helpful guesses for the final combined source models in the whole data frame. We applied the rotated PSFs and modeled the flux densities at given source positions. The DAOGROUP algorithm (Stetson 1987) was adopted to group sources within a separation of one beam FWHM. The best-fit model was then stored and used as prior information for the next iteration.

In the second iteration, we extracted sources that were (i) ALCS sources in the main sample or (ii) ALMA-faint Herschel sources. These two types of objects should be secure sources, and thus their flux densities should be positive in any Herschel band. The flux priors were given by the first iteration or aperture photometry if the sources were not modeled previously. In SPIRE bands, the initial flux guesses of sources at S/N_ALMA < 10 were assumed as RMS_map. With a similar PSF photometry routine, we modeled and subtracted sources with positive best-fit flux densities. In this intermediate step, most of the real sources in the current Herschel band have been extracted. The best-fit model and residual maps were recorded for the next iteration.

In the last iteration, we extracted tentative ALCS sources in the secondary sample. Source extraction was performed on the residual map of the second iteration, and the initial guesses of the flux densities were set to zero. Only sources modeled with positive flux densities were kept. The best-fit model and residual were then recorded as part of the final products.

Figures 2 and 3 display the Herschel scientific and residual images (i.e., before and after the source extraction) of two cluster fields observed in both the "deep" and "snapshot" modes, namely, A2744 (Figure 2; "deep") and M0417 (Figure 3; "snapshot"). Secure and tentative ALCS sources and ALMA-faint Herschel sources extracted at S/N > 2 are shown as open magenta, cyan, and green circles, respectively. No significant residual in SPIRE bands can be found within the ALCS footprint, shown as the region enclosed by the red solid line, although weak ring-shape residuals can be identified in PACS bands for a few very bright sources (f ≳ 50 mJy). This could be caused by the invalidity of point-source assumption or mismatch of PSF models, but our examination in Appendix A suggests no gain or loss of PACS flux densities through this PSF photometry routine.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a consequence of the limited depths of the Herschel data when compared with the deep ALCS data, a significant number of sources were not successfully extracted in the Herschel bands (see statistics in Section 3.5). Therefore, we only provide 3σ upper limits of their flux densities. As pointed out by R16 and other works, because the source positions are known from the ALMA data, the actual Herschel flux limit of nondetections is lower than the nominal confusion noise limit (e.g., Nguyen et al. 2010).

Based on the flux densities and their uncertainties modeled with PSF photometry, we calculated the median flux density uncertainty of the extracted secure sources in each band and each cluster field (main ALCS sample, Herschel S/N > 2). These uncertainties were obtained based on positional priors using the covariance matrix of least-squares PSF modeling, and thus we define them as RMS_prior. Table 2 presents all the measured RMS_prior, along with the RMS_map, which is directly measured from the 2D uncertainty map within the ALCS footprint. We also compare the RMS_prior with RMS_map measured for all our Herschel data in Figure 4, and we find a median ratio of RMS_prior to RMS_map of 0.43 ± 0.01 in the Herschel/SPIRE bands. A similar value can also be found for the PACS bands. This means that with prior knowledge of source positions, the actual 3σ limit of nondetection is around 1.3 times the local background RMS_map. Such an upper limit is adopted for all the Herschel nondetections in this work.

Zoom In Zoom Out Reset image size

The median 3σ depths derived for the 18 clusters observed in the "deep" mode are 7.5, 7.6, and 8.2 mJy at 250, 350, and 500 μm, respectively. These are consistent with the 3σ rms of deblended flux densities using the cross-identification procedure reported by Roseboom et al. (2010, for HerMES fields) and only slightly lower than the reported value in Swinbank et al. (2014) at 500 μm, which also included ALMA positional priors for deblending.

Several ALCS sources show secure close companions (i.e., angular separation less than 6'', which is one-third of the beam FWHM at 250 μm) at S/N ≥ 5 in the ALMA maps. Five of these seven systems have already been confirmed as lensed arcs or multiply lensed systems with HST or ALMA data. Due to the coarse resolution of the Herschel data, especially those of SPIRE, the flux density ratios among these source groups might be incorrectly modeled in Section 3.2.

In this step, we redistributed the Herschel flux densities of these source groups according to their ALMA flux density ratios. If a source was resolved on the native-resolution ALMA map with a major-axis FWHM less than 3'' (morphological parameters modeled with casa imfit; S. Fujimoto et al. 2022, in preparation), we adopted the ALMA flux density measured with a circular aperture of r = 2''. For sources with larger FWHMs, we adopted the best-fit ALMA flux densities derived from surface brightness profile modeling (assuming 2D Gaussian profile with imfit). For unresolved sources, we used the peak flux densities per beam measured on the 2''-tapered maps. We note that the redistribution of Herschel flux densities assumes a fixed FIR SED shape among blended sources in each group. Only one source from each blended group is considered for the discussions of dust temperature in Section 6.4.

We redistributed the Herschel fluxes for all the ALCS sources in the main sample and within a separation of 6''. This includes (i) ACT0102-ID215/224 (lensed galaxy pair known as "la Flaca"; Lindner et al. 2015; Wu et al. 2018; Caputi et al. 2021), (ii) ACT0102-ID223/294 (lensed galaxy pair; Wu et al. 2018), (iii) M0417-ID221/223 (two ALMA sources at S/N ≳ 6 with a separation of 23; K. Kohno et al. 2022, in preparation), (iv) M0553-ID133/190/249 (triply lensed arc; Ebeling et al. 2017; Sun et al. 2021b), (v) M1206-ID55/60/61 (known as the "Cosmic Snake"; Ebeling et al. 2009; Cava et al. 2018; Dessauges-Zavadsky et al. 2019), (vi) R0032-ID53/55/57/58 (lensed arc; Dessauges-Zavadsky et al. 2017), and (vii) R1347-ID145/148 (IR-bright lensed arc in R1347). In Figure 5 we display the HST, ALMA, and SPIRE images of all these source groups. Note that we only redistributed the 500 μm flux densities of ACT0102-ID215/224 because of a moderate separation (93).

Zoom In Zoom Out Reset image size

By comparing our PSF photometric results with those derived with aperture photometry, PSF photometry using MIR priors (R16) and different software (xid+; Hurley et al. 2017) as presented in Appendix A, we confirmed the quality of our Herschel flux density measurements. In Table C1 we present the Herschel photometric catalog of 180 ALCS sources in both the main and secondary samples (Section 2.1). The definition of these samples is also illustrated in Figure 6.

Zoom In Zoom Out Reset image size

Figure 7 shows the Herschel detection rates as functions of ALMA S/N cut from 100 to 500 μm. The rates of 100 and 250 μm detection (>2σ) are almost constant at S/N_ALMA ≥ 5. However, the detection rates at long wavelength (350 and 500 μm) are clearly correlated with the significance of ALMA sources. This is likely the consequence of (i) a larger beam size and stronger blending effect toward longer wavelength and (ii) a decreasing fraction of high-redshift sources (z > 3) toward lower 1.15 mm flux density (e.g., Béthermin et al. 2015; Casey et al. 2018; Popping et al. 2020). We also note one caveat that certain extragalactic ALMA surveys of rest-frame UV/optical-selected galaxies may have a selection bias against highly dust-obscured sources at z ≳ 3. For such surveys, the most accurate measurement of the redshift distribution can be obtained after the sample is spectroscopically complete (see Reuter et al. 2020; Chen et al. 2022).

Zoom In Zoom Out Reset image size

Among the total of 141 secure ALCS sources at S/N_nat ≥ 5 or S/N_tap ≥ 4.5 (58 of which fall in the PACS coverage), we successfully extracted 40, 44, 94, 77, and 55 sources at 100, 160, 250, 350, and 500 μm at above a 2σ significance, respectively. The SPIRE detection rate in the "deep"-mode clusters is higher than that in the "snapshot"-mode clusters by ∼10%. A total of 113 (99) ALCS sources were detected at >2σ ( >3σ) in one Herschel band at least (including eight brightest cluster galaxies (BCGs)), and 91 sources were detected at >2σ in two Herschel bands at least.

Only 28 secure ALCS sources (20% of the total 141 sources) cannot be extracted at S/N >2 in any Herschel band. The 16th, 50th, and 84th percentiles of 1.15 mm flux densities of these sources (0.45, 0.92, and 1.41 mJy) are smaller than those of Herschel-detected sources (0.66, 1.22, and 3.00 mJy). Except for M1206-ID58 as a BCG at z = 0.440, we refer to the remaining 27 sources as Herschel-faint galaxies in later analysis (Section 4.3; also called Herschel-dropout galaxies in Boone et al. 2013).

Among the 39 tentative ALCS sources at S/N_nat = 4 − 5 and S/N_tap < 4.5, 22 of them can be extracted above a 2σ significance in at least one Herschel band, including one BCG and one mismatched source. The remaining 17 sources are undetected in any Herschel band, including one BCG (R0032-ID162). These sources are excluded for further analysis because of a higher false-ID rate.

We further justify such a 2σ detection threshold by calculating the joint probability of spurious sources through a χ² statistic of the detection significance in all available Herschel bands. For all ALCS sources extracted at S/N >2 in any Herschel band, only five sources (4% ± 2%) exhibit p-values of spurious detection at above 0.01, including three secure ALCS sources (ACT0102-ID118, P171-ID69, P171-ID177) and two tentative sources (AS295-ID269, M1115-ID33). The largest p-value = 0.03 is seen for source P171-ID69 with S/N_ALMA = 19.7. Therefore, we conclude that with the single-band 2σ detection threshold the number of spurious Herschel sources will be on the order of unity.

We perform FIR SED modeling of our sample with the best available redshifts (z_best) using magphys (da Cunha et al. 2008, 2015). Here the z_best is either spectroscopic redshift (z_spec), published HST-derived z_phot (Section 2.5), or FIR z_phot (priority from high to low). Redshift uncertainty is propagated into the uncertainties of derived physical properties through a Monte Carlo sampling of z_phot likelihood when z_spec is not available. In order to derive and validate FIR z_phot for sources without z_spec or HST z_phot, we also perform simultaneous FIR SED fitting and photometric redshift estimates of our sample using magphys+photoz (Battisti et al. 2019), the photoz extension of magphys.

magphys assumes a Chabrier (2003) initial mass function (IMF), a continuous delayed exponential star formation history with random starburst, and energy balance between dust absorption in the UV and reemission in the IR. At FIR wavelengths, the dust model assumed by magphys includes a warm (30–80 K) and a cold (20–40 K) component, and the prior distribution of luminosity-weighted dust temperature (T_dust) peaks around 37 K with a 1σ dispersion of ∼20%. Such a dust temperature is comparable to those of widely adopted spectral templates at around median L_IR ∼10¹² L_⊙, including Chary & Elbaz (2001) and Magdis et al. (2012, which is based on the model of Draine & Li 2007). For a full description of the models and parameters assumed by magphys, see da Cunha et al. (2008, 2015) and Battisti et al. (2019). Here we only include five or three bands of Herschel data and ALMA 1.15 mm flux densities for our SED modeling. Further optical/NIR counterpart matching, uniform HST and Spitzer photometry, and panchromatic SED fitting will be presented by another paper of the ALCS collaboration, and certain conclusions on T_dust and redshift distribution depending on the FIR z_phot values may be further revised.

Here we highlight several caveats of our SED modeling obtained with magphys. First of all, the accuracy of FIR z_phot is subject to a well-known degeneracy between redshift and dust temperature, typically showing an error around Δz ∼ 0.2(1 + z). In addition, magphys+photoz assumes a prior redshift distribution peaking at z ∼ 1.7, and in practice we find that such a prior will lead to an artificial shift of z_phot estimate toward such a redshift. To address this issue, we adopt a uniform redshift prior instead. Furthermore, the nonthermal emission of BCGs seen at 1.15 mm cannot be properly modeled, and thus their boosted ALMA flux densities (e.g., M1931-ID41; Fogarty et al. 2019) will lead to a wrong estimate of IR luminosity and SFR. Therefore, we do not perform SED modeling for all the known BCGs. We also note that only three lensed ALCS sources are detected in X-ray among 31 cluster fields with publicly available Chandra data (A370-ID110, M0416-ID117, M0329-ID11; these sources will be discussed by Uematsu et al. from the collaboration). We also estimate the upper limit of the AGN contribution to the derived IR luminosities. We assume the SKIRTOR model (Stalevski et al. 2012, 2016) for a type 2 AGN SED with an inclination angle of 70°. To estimate the upper limit of X-ray luminosity, we used a simple absorbed power-law model with a photon index of 1.9 and an intrinsic absorption of $\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{2})=23$. The median X-ray luminosity is L_X < 2 × 10⁴³ erg s⁻¹ for X-ray-undetected sources, corresponding to an IR luminosity of ≲6 × 10⁹ L_⊙. Therefore, the AGN contamination will not be a concern for the majority (≳95%) of ALCS sources, but we also note that in the case of a Compton-thick AGN ($\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{2})\gt 24$), the upper limit on L_X can be larger by more than an order of magnitude. Finally, magphys can only provide weak constraints on the physical properties of Herschel-faint galaxies individually, which are specifically discussed in Section 4.3.

Table C2 presents a summary of the best-fit galaxy properties of 125 ALCS sources detected at S/N ≥ 4, including 47 sources that are spectroscopically confirmed and an additional 42 sources with cataloged HST z_phot. This sample, further referred to as the ALCS-Herschel joint sample, includes 105 secure (i.e., the main sample) and 20 tentative ALCS sources (the secondary sample) detected above 2σ in one Herschel band at least, except for 11 BCGs and one special source (A2537-ID06) due to the poorness of SED fitting. A2537-ID06 is only detected with ALMA at S/N = 4.2 and offset from a passive cluster dwarf galaxy by 09, and therefore it is likely a false detection with Herschel/SPIRE fluxes coming from an ALMA-faint Herschel source. The definition of this sample is also visualized as the red filled regions in Figure 6. The postage stamp images and best-fit FIR SEDs of these 125 sources are shown in Appendix B.

We further model the dust temperature of sources in Table C2 with the Herschel and ALMA data. We fit the dust continuum emissions of all sources with a modified blackbody (MBB) using the best available redshifts. The dust absorption coefficient is assumed to be κ = 0.40 × (ν/250)^β in units of cm² g⁻¹, where ν is the frequency in GHz in the rest frame. We assume a fixed dust emissivity of β = 1.8, which was widely adopted in previous studies (e.g., Díaz-Santos et al. 2017; Dudzevičiūtė et al. 2020; Sun et al. 2021b) and supported by a recent 2 mm study of SMGs at z ≃ 1 − 3 (da Cunha et al. 2021). Following previous work including Greve et al. (2012) and Sun et al. (2021b), we only fit the SED over a rest-frame wavelength of 50 μm to avoid an optically thick regime and eliminate any possible contribution of warm dust components at shorter wavelength. We note that luminous SMGs (L_IR ≃ 10^12.5 − 10¹³ L_⊙) are found to be optically thick at λ_thick ∼ 100 μm (e.g., Spilker et al. 2016; Simpson et al. 2017; Dudzevičiūtė et al. 2020), but we argue that less luminous ALCS sources have lower dust mass densities (e.g., Sun et al. 2021b) and therefore smaller optical depths at the same wavelength in the rest frame. The best-fit dust temperatures and masses (M_dust) are also presented in Table C2. The uncertainty of redshift is propagated into that of the dust temperature if the z_spec is unknown. The dust masses derived from this fitting procedure are also consistent with those from magphys.

We also note that many works use the rest-frame wavelength of FIR SED peak (λ_peak) to quantify the luminosity-weighted dust temperature (e.g., Casey et al. 2018; Reuter et al. 2020; Burnham et al. 2021). This is because λ_peak is less dependent on dust opacity assumption compared with T_dust. Under the dust absorption coefficient assumption that we adopt, the conversion between λ_peak and T_dust is ${\lambda }_{\mathrm{peak}}=3\times {10}^{3}{T}_{\mathrm{dust}}^{-1}$ μmK.

To validate the FIR photometric redshifts derived with magphys+photoz, we first compare the FIR z_phot values of 47 spectroscopically confirmed sources with their z_spec values in the left panel of Figure 8. The median redshift discrepancy is (z_phot − z_spec)/(1 + z_spec) = −0.01 ± 0.04, indicating an excellent agreement. We apply a Kolmogorov–Smirnov (K-S) test on the discrepancy between z_phot and z_spec divided by the uncertainty of FIR z_phot, in comparison with the standard normal distribution. We conclude that the standard deviation of (z_phot − z_spec)/(1 + z_spec) is well predicted by the uncertainty of FIR z_phot (p-value = 0.59).

Zoom In Zoom Out Reset image size

The comparison of FIR and HST z_phot values is shown in the right panel of Figure 8. We find a general agreement of z_phot estimate between HST and FIR SED modeling, albeit with a significant dispersion. We note that the best-fit linear slope of z_phot,FIR(z_phot,HST) is only 0.31 ± 0.08. However, this is mainly contributed by roughly four tentative sources in the secondary sample with z_phot,HST ∼ 4 but z_phot,FIR ∼ 2, which may be subject to a high false-ID rate (∼18%). The discrepancy between the two z_phot values, defined as (z_phot,HST −z_phot,FIR)/[1 + (z_phot,HST + z_phot,FIR)/2], is observed to be ${0.01}_{-0.23}^{+0.27}$ (16th, 50th, and 84th percentile) with a typical uncertainty of 0.21. We perform a K-S test on the photometric redshift discrepancies divided by their uncertainties of the cross-matched sources. The null hypothesis that the FIR-HST z_phot discrepancies relative to their uncertainties are drawn from a standard normal distribution cannot be ruled out (p-value = 0.33), reinforcing the agreement between these two photometric redshift estimates.

Among the 27 Herschel-faint galaxies reported in Section 3.5, two of the sources have been spectroscopically confirmed. They are R0600-ID164 as a [C ii]-emitting lensed arc at z = 6.072 (Fujimoto et al. 2021; Laporte et al. 2021) and R0032-ID32 as the faintest component of a multiply lensed arc at z = 3.631 (Dessauges-Zavadsky et al. 2017). In addition to this, eight sources exhibit HST counterparts with tabulated photometric redshifts (median z_phot = 2.0 ± 1.0; Molino et al. 2017; Coe et al. 2019). Postage stamp images of these 10 sources are shown in Appendix B (Figure B3). FIR SEDs of these galaxies are also modeled with magphys, and we note that the adopted IR spectral templates are essentially MBB spectra at around T_dust = 35 ± 6 K without exceeding the Herschel nondetection limits except for SM0723-ID93 (S/N_ALMA = 4.6 in 2''-tapered map), which is likely a random association between a cluster dwarf galaxy and spurious ALMA source.

The remaining 17 sources do not have cross-matched HST z_phot because they are intrinsically faint shortward of 1.6 μm and/or out of HST/WFC3-IR coverage (Figure B4 in Appendix B). Such an NIR-dark (also often called "H-dropout/faint") nature suggests that they are likely dust-obscured SFGs at z ∼ 4 that have raised general interest in recent studies (e.g., Simpson et al. 2014; Franco et al. 2018; Alcalde Pampliega et al. 2019; Wang et al. 2019; Williams et al. 2019; Yamaguchi et al. 2019; Dudzevičiūtė et al. 2020; Gómez-Guijarro et al. 2022; Smail et al. 2021; Sun et al. 2021a). For each individual source, the nondetections in the HST and Herschel bands prevent us from deriving useful constraints of their redshifts and physical properties.

To address this issue, we stack the Herschel residual images of Herschel-faint galaxies without spectroscopic or HST photometric redshifts. We note that four out of five sources at S/N_ALMA < 5.5 do not show any counterpart in HST or Spitzer bands (ACT0102-ID11, M2129-ID24, ACT0102-ID251, and R1347-ID51). Therefore, these sources could be spurious detections, or highly obscured galaxies at very high redshift (i.e., similar to R0600-ID67 and R0949-ID19, the brightest Herschel-faint galaxies in Figure B4 that do not show any HST or Spitzer counterpart, and also the ALMA-only [C ii] emitters at z >6 reported recently by Fudamoto et al. 2021). Our stacking analysis suggests that including these sources will lead to a lower S/N in SPIRE 350 and 500 μm bands, and therefore we only present the stacked SEDs of 12 sources at S/N_ALMA >5.5.

We first normalize the Herschel/SPIRE residual and uncertainty images of all sources by their ALMA flux densities. Here the residual images are the scientific images with all the other Herschel sources subtracted assuming point-source models as described in Section 3.2. PACS images are not stacked because of the unavailability for most sources. We stack all the images in each SPIRE band using an inverse variance weighting method. The stacked SPIRE images are presented in Figure 9. We measure the flux densities of stacked sources using an aperture of r_aper = 18'' with appropriate aperture correction factors. The sky background is subtracted using the median of sigma-clipped local annulus, and photometric uncertainty is estimated from the rms of that.

Zoom In Zoom Out Reset image size

The stacked source can be detected at ∼3σ in SPIRE 350 and 500 μm band while remaining undetected at 250 μm (<2σ). The stacked FIR SED is shown below the stacked SPIRE images in Figure 9. We also show the median FIR SEDs of ALCS sources that are detected with Herschel to visualize the clear difference in the continuum strength at below 500 μm. Here all the photometric data are normalized to the median ALMA flux density (1.09 mJy). With magphys+photoz, we derive a median z_phot of 4.2 ± 1.2 (uncertainties denote the 16th and 84th percentiles of the likelihood distribution), IR luminosity of 10^12.1±0.3 μ⁻¹ L_⊙, and SFR of ${100}_{-50}^{+100}$ μ⁻¹ M_⊙ yr⁻¹ before lensing magnification correction. The derived redshift is consistent with those of HST H-faint galaxies in previous studies (e.g., Simpson et al. 2014; Wang et al. 2019; Alcalde Pampliega et al. 2019; Dudzevičiūtė et al. 2020; Sun et al. 2021a).

Using the composite SED templates of AS2UDS SMGs at z >3 (Dudzevičiūtė et al. 2020) and ALESS SMGs at z >3.5 (da Cunha et al. 2015), we derive best-fit z_phot of 3.8 ± 0.4 and 5.2 ± 0.6, respectively. Given the large scattering of T_dust for sources at given luminosity and redshift (e.g., Schreiber et al. 2018; Dudzevičiūtė et al. 2020), the redshift uncertainty can be significantly underestimated with single template matching techniques. In addition to this, the median IR luminosity of SMGs in the ALESS z >3.5 sample is ∼10 times higher than that of stacked Herschel-faint sources and therefore likely exhibits a higher T_dust. This likely leads to an overestimated z_phot through template matching for the ALESS z >3.5 sample.

We calculate the lensing magnification factor (μ) based on the best available redshifts. We adopt two sets of parametric lens models: the so-called Zitrin-NFW lens models (Zitrin et al. 2013, 2015) for all the HFF and CLASH clusters, and GLAFIC lens models (Oguri 2010; Okabe et al. 2020) for the RELICS clusters. The lensing magnification is derived using the maps ³⁹ of projected cluster mass surface density (κ) and weak-lensing shear (γ) at the centroid of ALCS source as $\mu ={[{(1-\kappa \cdot {D}_{{ls}}/{D}_{s})}^{2}-{(\gamma \cdot {D}_{{ls}}/{D}_{s})}^{2}]}^{-1}$, where D_ls is the angular diameter distance between the lens and the source and D_s is the angular diameter distance to the source. We assume no magnification (μ = 1) for sources within or in front of the cluster fields (z_s < z_cl + 0.1), following R16. For 16 Herschel-faint galaxies without cataloged redshifts (Section 4.3), we calculate their magnification at z_s = 4.2 uniformly.

The distribution of magnification factors of the ALCS-Herschel joint sample and Herschel-faint sample is shown in Figure 10. The 16th, 50th, and 84th percentiles of the distribution of μ are 1.8, 2.6, and 5.2, and seven sources exhibit a strong magnification with μ >10. The median lensing magnification suggests that the ALCS has reached a great depth that typical surveys in blank fields would require a ∼7 × longer observing time to achieve.

Zoom In Zoom Out Reset image size

Among the 125 sources in our ALCS-Herschel joint sample, six groups of sources have been spectroscopically confirmed as multiply imaged systems. These include (i) ACT0102-ID118/215/224 (Caputi et al. 2021), (ii) M0417-ID46/58/121 (K. Kohno et al. 2022, in preparation), (iii) M0553-ID133/190/249 (Ebeling et al. 2017), (iv) M1206-ID27/55/60/61 (Ebeling et al. 2009), (v) R0032-ID53/55/57/58 (Dessauges-Zavadsky et al. 2017), and (vi) R1347-ID145/148 (Richard et al. 2021). Additionally, M1931-ID47/55/61, R0032-ID208/281/304, and R0032-ID127/131/198 are also multiply lensed candidates yet to be spectroscopically confirmed, including HST H-faint sources M1931-ID47 and R0032-ID208/281/304, which are part of the Herschel-faint sample. This reduces the number of independent sources in Table C2 to 109. Multiply imaged sources of the same system are shown separately in diagrams for individual sources (e.g., Figure 8) but only counted once in all statistics of physical properties (e.g., redshifts, SFRs, and dust temperatures) in Section 6.

The detailed characterization of the lensing magnification uncertainty (σ_μ ) will be presented by a forthcoming source count paper of the ALCS collaboration (S. Fujimoto et al. 2022, in preparation), and here we only present a qualitative estimate with a few quick methods. First of all, we estimate the σ_μ based on the uncertainty of photometric redshift for sources without spectroscopic confirmation. The typical lensing magnification error propagated from z_phot uncertainty is σ_μ /μ ∼ 5%. In addition to this, we also compare the derived magnification factors using ALMA and HST source centroids, and the typical difference is found to be less than 1%. In order to quantify the σ_μ caused by the extended source profile, we also measure the average magnification factor within a radius of 06 from the ALMA source centroid. Such an effect is negligible (σ_μ /μ ≲ 2%) in most cases except for M0159-ID24 and R0032-ID57 (Figure 5) because of the galaxy−galaxy lensing effect, which provides a stronger magnification gradient over a smaller angular scale.

As an alternative method, we also analyze the uncertainty map of magnification presented by Zitrin et al. (2013, 2015) in the CLASH cluster fields, obtained through Markov Chain Monte Carlo (MCMC) fitting routines when the cluster mass models were constructed. We find that the σ_μ /μ is around ∼5% at μ = 3 assuming a fixed sky position and redshift. However, this uncertainty could be as large as ∼10% at μ = 10 and ∼50% at μ = 100, indicating that sources with larger magnifications are subject to a larger relative uncertainty, and therefore their intrinsic physical properties (e.g., L_IR and SFR) are more uncertain.

Finally, based on the standard deviation of magnifications predicted by different lens models of the same clusters produced by different methods and groups (i.e., GLAFIC, CATS, Zitrin-NFW, and Zitrin Light-Traces-Mass (LTM) for Frontier Field clusters; Oguri 2010; Richard et al. 2014; Zitrin et al. 2015; Kawamata et al. 2016, 2018), we find a typical magnification uncertainty of σ_μ /μ ∼ 20%. This is comparable to the σ_μ reported in R16 (σ_μ = 0.5, ∼20% of the median magnification for sources in this work).

We further study the FIR colors of Herschel-detected ALCS sources. Figure 11 displays four color–color diagrams of our sample, covering the full wavelength range from 100 μm to 1.15 mm. Based on the average IR spectral templates in Rieke et al. (2009), we superimpose the redshift evolution tracks of a typical SFG (L_IR = 10^10.25 L_⊙, corresponding to an SFR_IR of ∼2 M_⊙ yr⁻¹ assuming the conversion factor in Kennicutt & Evans 2012) and a ULIRG (L_IR = 10¹² L_⊙, corresponding to an SFR_IR of ≳100 M_⊙ yr⁻¹) on all of the four color–color diagrams. The dust temperatures of the SFG and ULIRG templates are around 20 and 40 K. We also compare the distribution of our ALMA-selected sample with the R16 sample in two of the color–color diagrams where the 1.15 mm flux density is not invoked.

Zoom In Zoom Out Reset image size

6.1.1. Selection Bias

6.1.2. Comparison with Rawle et al. (2016)

The top left panel of Figure 12 displays the redshift distribution of sources in our ALCS-Herschel joint sample. The 16th, 50th, and 84th percentiles of redshifts are 1.11, 1.90, and 2.95 (1.05, 1.81, and 2.96) for the main (full) joint sample, respectively, and no obvious difference in redshift distribution can be found among sources in clusters observed with "deep" and "snapshot" mode. The uncertainty of photometric redshifts is accounted for by the derivation of percentiles through a Monte Carlo sampling of the probability distributions. Among the joint sample, spectroscopically confirmed sources are generally at slightly lower redshifts. The 16th, 50th, and 84th percentiles of z_spec values are 1.06, 1.55, and 2.90 (1.03, 1.49, and 2.70) for 26 (31) independent sources in the main (full) sample, respectively. However, we note that this ALCS-Herschel joint sample does not include 27 Herschel-faint galaxies that likely reside at higher redshifts (Section 4.3).

Zoom In Zoom Out Reset image size

To estimate the median redshift of the full ALCS sample, we must include the Herschel-faint sources. Based on the probability distribution of 17 optical/NIR-dark Herschel-faint galaxies (intrinsically 15 sources after removing multiply lensed images) derived with magphys+photoz (z_phot =4.2 ± 1.2), the 16th, 50th, and 84th percentiles of redshifts are 1.15, 2.08, and 3.59, respectively, for secure ALCS sources in the main sample (middle left panel of Figure 12; note that cluster member galaxies are not included).

The median redshift of secure ALCS sources is higher than that of the main ASPECS sample (z_med = 1.80 ± 0.15) in the Hubble Ultra Deep Field (HUDF; Aravena et al. 2020) but smaller than those of slightly shallower surveys, including AS2UDS (z_med = 2.61 ± 0.08 down to ∼1 mJy at 850 μm; Dudzevičiūtė et al. 2020), ASAGAO (z_med = 2.38 ± 0.14 down to ∼0.2 mJy at 1.2 mm; note that NIR-dark ALMA sources at z ∼ 4 are not included; Yamaguchi et al. 2019, 2020), and GOODS-ALMA (z_med = 2.56 ± 0.13 down to ∼0.5 mJy at 1.1 mm in their main sample; Franco et al. 2018; Gómez-Guijarro et al. 2022). The median redshift of ALCS sources is also smaller than those of millimeter-selected strongly lensed SMGs with wider but shallower surveys (e.g., z_med = 2.9 ± 0.1 with Planck's dusty GEMS sample, Cañameras et al. 2015; z_med = 3.9 ± 0.2 with the SPT sample, Reuter et al. 2020), as well as that of unlensed sources selected with the MORA survey at 2 mm (z_med = 3.6 ± 0.3 down to ∼ 0.3 mJy; Casey et al. 2021).

The median redshift of millimeter and submillimeter sources as a function of flux density limit is a key test for the evolution model of dust-obscured star formation history of the universe (see a review by Hodge & da Cunha 2020). The right panel of Figure 12 shows the best available redshifts versus intrinsic 1.15 mm flux densities of all ALCS sources in the main sample (excluding cluster member galaxies; filled squares) and Herschel-detected ALCS sources in the secondary sample (open downward-pointing triangles). Assuming the probability distribution of redshifts for optical/NIR-dark Herschel-faint galaxies derived with magphys+photoz, we are able to derive the median redshift of secure ALCS sources above the given 1.15 mm flux density cuts (i.e., z_med( >f₁₁₅₀); magenta circles). The uncertainty of photometric redshift is propagated into the uncertainty of z_med( >f) through a Monte Carlo sampling of z_phot likelihood.

The median redshifts of secure ALCS sources decrease from z_med = 2.40 ± 0.29 at f₁₁₅₀ >1 mJy to z_med = 2.04 ± 0.12 at f₁₁₅₀ >0.1 mJy, where our survey is ∼50% complete given the depth and median lensing magnification. We also perform a K-S test for the redshift distributions of ALCS and ASPECS sources above 0.1 mJy, and no obvious difference can be found (p-value = 0.64). For the spectroscopically confirmed sample, the decrease of median redshift is not conspicuous because of a smaller sample size (z_med = 1.71 ± 0.87 at >1 mJy to z_med = 1.60 ± 0.21 at >0.1 mJy). We note that a higher-redshift assumption (e.g., z_med ∼ 6) of optical/NIR-dark Herschel-faint galaxies will not change our measurements of z_med( >f) despite a larger standard error. We also compute the median redshifts of ALCS sources in logarithmic flux density bins from 0.1 to 2 mJy (bin size is 0.1 dex), and the null hypothesis that there is no monotonic relation between z_med and f₁₁₅₀ can be ruled out (p-value = 0.03, computed from Spearman's rank correlation coefficient ρ = 0.56).

Further compared with previous 1.1–1.2 mm surveys of SMGs (Michałowski et al. 2012; Yun et al. 2012; Brisbin et al. 2017; Dunlop et al. 2017; Miettinen et al. 2017; Aravena et al. 2020; Gómez-Guijarro et al. 2022), the derived z_med( >f₁₁₅₀) function suggests an increasing fraction of z ≃ 1 − 2 galaxies at f₁₁₅₀ < 1 mJy (also shown with ASPECS and the semiempirical model presented in Popping et al. 2020). A linear fitting to the ${z}_{\mathrm{med}}[\gt \mathrm{log}({f}_{1150})]$ measurements suggests that the positive correlation is significant ( >5σ). We note that z_med( >f₁₁₅₀) below 0.5 mJy was poorly probed with previous surveys because of either limited volume (ASPECS-like) or relatively shallower depth (e.g., GOODS-ALMA and single-dish surveys). Similar z_med( >f) trends were also reported at 870 μm (Ivison et al. 2007; Stach et al. 2019; Simpson et al. 2020; Birkin et al. 2021; Chen et al. 2022).

The overall z_med( >f₁₁₅₀) function obtained with the main ALCS sample is lower than that predicted by the shark semianalytic model presented in Lagos et al. (2020; the offset is Δz ∼ − 0.4) but very close to that empirically modeled by Casey et al. (2018; Δz ∼ 0.1) assuming a negligible contribution (<10%) of SMGs to the cosmic SFR density at z >4, also known as the "dust-poor" scenario. Nevertheless, the observed z_med( >f₁₁₅₀) of our sample is much lower than that of the "dust-rich" scenario modeled by Casey et al. (2018, Δz ∼ −0.7), in which SMGs dominate the star formation ( ∼90%) at z >4. This may suggest that the majority of cosmic star formation at z >4 will be hosted in an unobscured environment, similar to the conclusions drawn in Dudzevičiūtė et al. (2020), Bouwens et al. (2020), Casey et al. (2021), and Zavala et al. (2021). However, we note that the decreasing obscured fraction of cosmic SFR density at z >4 cannot be considered as the unique cause of positive z_med( >f₁₁₅₀) relation. Other scenarios, e.g., a steep faint-end slope of dust mass function at z = 1 ∼ 2, can lead to a similar z_med( >f₁₁₅₀) function to that observed. Detailed characterization of IR luminosity function and cosmic SFR density evolution will be presented in a future work.

We study the distribution of the intrinsic SFRs (total SFR derived by magphys and corrected for lensing magnification) for ALCS sources in Figure 13. For all ALCS sources in the main sample (excluding cluster member galaxies and including Herschel-faint galaxies), the 16th, 50th, and 84th percentiles of the distribution are 40, 94, and 178 M_⊙ yr⁻¹, respectively, slightly larger than those of the GOALS sample (16th, 50th, and 84th percentiles of SFR as 25, 45, and 175 M_⊙ yr⁻¹, respectively; Armus et al. 2009; Howell et al. 2010), which is mostly composed of LIRGs in the local universe. Herschel-detected sources in the secondary sample exhibit LIRG-like SFRs with 16th, 50th, and 84th percentiles of the distribution of 10, 25, and 97 M_⊙ yr⁻¹, respectively.

Zoom In Zoom Out Reset image size

We also compare the intrinsic SFRs with those of ALMA sources reported by AS2UDS (707 SMGs at z_med = 2.6; Stach et al. 2019; Dudzevičiūtė et al. 2020), GOODS-ALMA (35 sources at z_med = 2.7 studied in Franco et al. 2020; note that the SFRs are recomputed assuming a Chabrier 2003 IMF), and ASPECS survey (32 sources in the main sample at z_med = 1.8; Aravena et al. 2020; González-López et al. 2020). Although the SFR distribution is similar in the image plane (lensing uncorrected), the median source-plane SFR (intrinsic) of ALCS sources is lower than that of conventional SMGs (i.e., f >1 mJy at 850 μm) in the AS2UDS sample by a factor of ∼3 (median SFR = 236 ± 8 M_⊙ yr⁻¹; although 92 sources in the AS2UDS sample are at <100 M_⊙ yr⁻¹ with a median redshift of 2.0 ± 0.1). Below 30 M_⊙ yr⁻¹, which is the median SFR of 32 sources in the main ASPECS sample down to a fidelity of 50%, our sample contains more sources than ASPECS (N = 16), including 11 Herschel-detected and ∼10 Herschel-faint sources at S/N ≥ 5, as well as 12 tentative sources in the secondary sample. Because of the same software (magphys) and SFR tracer (FIR) being used, these comparisons are fair. The median ratio between SFR and L_IR is found to be ${10}^{-10.1}\,{M}_{\odot }\ {\mathrm{yr}}^{-1}\ {L}_{\odot }^{-1}$ for both the main samples of ALCS and ASPECS, which is only offset from the conversion factor in Kennicutt & Evans (2012) by 0.1 dex (assuming Chabrier IMF). We also note that the ASPECS sample presented by Aravena et al. (2020) did not enforce any Herschel detection, and no SPIRE flux density information was given because of blending issues, making it difficult to directly and reliably constrain the dust temperature from the FIR SED. This highlights the uniqueness of our ALCS-Herschel joint sample as the probe of (sub-)LIRG population at z ≲ 2.

Figure 14 shows the distribution of intrinsic IR luminosity and SFR as functions of redshift, highlighting the less vigorous star formation among the ALCS sources compared with unlensed SMGs (Dudzevičiūtė et al. 2020). We calculate the nominal detection limits of infrared galaxies in the ALMA (>0.3 mJy at 1.15 mm) and Herschel ( >2σ in the "deep" mode; ≳5 mJy) bands assuming a ULIRG template at L_IR = 10¹² L_⊙ (Rieke et al. 2009). We find that the intrinsic luminosities of 29% (37%) of sources in the ALCS-Herschel joint sample are lower than the nominal detection limit of Herschel (ALMA). This demonstrates the benefit of a sky survey in lensing cluster fields, especially for Herschel because the lensing magnification will allow us to extract and study sources below the nominal confusion noise limit.

Zoom In Zoom Out Reset image size

Sources whose observed flux densities are below either of the ALMA or Herschel detection thresholds are excluded from the ALCS-Herschel joint sample (Table C2). These sources are referred to as ALMA-faint Herschel sources (Section 3.2) or Herschel-faint galaxies (Section 4.3; also included in Figure 14) in our study. Based on the left panel of Figure 14, we find out that the ALMA-faint Herschel sources are likely distributed at z ≃ 0 − 2, consistent with the typical redshifts of Herschel sources detected in the Frontier Fields (R16). ALCS sources without secure Herschel detection primarily reside at z ≳ 3, consistent with the redshifts of Herschel-faint galaxies presented in Section 4.3.

We also find that ∼77% of ALCS sources (excluding cluster member galaxies) host a lower SFR than the star formation "main sequence" (MS) at a fixed stellar mass of M_star = 10¹¹ M_⊙ (Speagle et al. 2014). Assuming that the majority of ALCS sources are on the MS (e.g., Aravena et al. 2020; Sun et al. 2021b), this suggests a lower intrinsic stellar mass than unlensed SMG samples in the literature (e.g., median M_star = 10^11.1 M_⊙ in Dudzevičiūtė et al. 2020). The median stellar mass is likely around 10^10.6±0.2 M_⊙ based on a comparison with the MS across z ≃ 1 − 4 (Speagle et al. 2014), which will be further analyzed and constrained by a future work from the collaboration. The 1σ scattering of measured SFRs around the SFR_MS(M_star = 10^10.6 M_⊙) is 0.4 dex.

We study the relation between the dust temperature and intrinsic IR luminosity of sources in the ALCS-Herschel joint sample as shown in the left panel of Figure 15. Sources with redistributed Herschel fluxes (Section 3.3) are considered only once among those in each blended group. We also compare our sample with a variety of galaxies, including nearby galaxies (KINGFISH sample; Skibba et al. 2011; Hunt et al. 2015), local (U)LIRGs (GOALS sample; Díaz-Santos et al. 2017), low-redshift LIRGs (HerMES/PEP sample at z ≲ 1; Symeonidis et al. 2013), 850 μm selected SMGs at z_med = 2.6 (AS2UDS sample; Dudzevičiūtė et al. 2020), 450 μm selected SMGs at z_med = 1.8 (STUDIES sample; Lim et al. 2020), and galaxy-lensed SMGs at z_med = 3.9 (SPT sample; Spilker et al. 2016; Strandet et al. 2016; Reuter et al. 2020).

Zoom In Zoom Out Reset image size

The median dust temperature is 32.0 ± 0.5 K for the ALCS-Herschel joint sample (with or without FIR-z_phot galaxies included). Spectroscopically confirmed sources have a slightly higher median dust temperature (33.6 ± 0.9 K) than that of sources without z_spec values (31.9 ± 0.5 K). We identify a weak positive L_IR–T_dust relation as ${T}_{\mathrm{dust}}\,=(2.2\pm 1.7)\mathrm{log}({L}_{\mathrm{IR}}/{10}^{12})+(32.6\pm 0.7)$, where the units of L_IR and T_dust are L_⊙ and K, respectively. Such a weak L_IR–T_dust relation (Spearman's ρ = 0.22, p-value = 0.03; uncertainties of T_dust are considered) is not consistent with the strong relations drawn from most of the high-redshift reference samples (see also Burnham et al. 2021). It is, however, consistent with some cosmological simulation results, including Liang et al. (2019).

We argue that a weak observed L_IR–T_dust relation is a consequence of an inhomogeneous L_IR(T_dust) selection threshold caused by the lensing effect. The L_IR–T_dust relation is a joint effect of both physics (the Stefan–Boltzmann law) and selection. If lensing magnification is not applied, the detection threshold of L_IR at each given redshift and millimeter flux density (e.g., z = 2 and f₁₁₅₀ = 0.3 mJy) will be exactly a monotonic function of T_dust (e.g., Lim et al. 2020). Combined with the selection bias toward high-L_IR and high-T_dust sources at higher redshifts in the Herschel bands, a positive L_IR–T_dust relation could be identified. In the image plane, we do find a positive μ L_IR–T_dust correlation at 4σ significance. However, with the lensing magnification, we are able to detect sources with lower L_IR at given T_dust inhomogeneously, leading to a weaker L_IR–T_dust relation in the source plane. This is also seen with the strongly lensed SPT sources (Spilker et al. 2016; Reuter et al. 2020), where the significance of L_IR–T_dust relation is also around 2σ, despite that SPT sources are far more luminous. We also note that a strong L_IR–T_dust relation can be identified for the joint sample of ALCS and SPT. This is because the sources of two surveys are selected in distinct ranges of IR luminosities, and the increase of T_dust over a wider range of L_IR becomes significant enough (T ∝ L^0.25 R^−0.5 according to the Stefan–Boltzmann law).

We find that at an intrinsic IR luminosity between 10¹¹ and 10¹² L_⊙ (i.e., LIRGs; $z={1.5}_{-0.4}^{+0.5}$ for 16th, 50th, and 84th percentiles of the redshift distribution) the dust temperatures of ALCS sources resemble those of local analogs (Díaz-Santos et al. 2017) and low-redshift LIRGs selected by Herschel/SPIRE in cosmological deep fields (Symeonidis et al. 2013). This indicates no or weak evolution of the dust temperature of LIRG-like galaxies from z ∼ 2 to the local universe, consistent with the conclusion made based on lensed HLS sources at z_med = 1.9 on the Σ_IR–T_dust plane (Sun et al. 2021b).

At L_IR ≳ 10¹² L_⊙, ALCS sources ($z={2.3}_{-0.9}^{+0.8}$ for 16th, 50th, and 84th percentiles of the redshift distribution) show lower dust temperatures than those of local ULIRGs, resembling 850 μm selected SMGs at z_med = 2.6 (AS2UDS sample; Dudzevičiūtė et al. 2020) except for one warm outlying system (R0032 lensed arc at z = 3.63, T_dust ∼ 50 K; Figure 5 and Dessauges-Zavadsky et al. 2017). As a result, the difference in T_dust between LIRGs and ULIRGs in the ALCS-Herschel joint sample at z ≳ 1 is not significant, which is reflected by the weak slope of the L_IR–T_dust relation. Previous works have reported a lower dust temperature or larger optical depth in z ≃ 1 − 3 SMGs compared to their local analogs (e.g., Symeonidis et al. 2009, 2013; Hwang et al. 2010; Swinbank et al. 2014). As also pointed out by a few of these studies, the high IR luminosity with relatively low dust temperature is likely caused by a larger size of the star-forming region in SMGs (R_e = 1 ∼ 2 kpc, e.g., Ikarashi et al. 2015; Simpson et al. 2015; Hodge et al. 2016; Rujopakarn et al. 2016; Fujimoto et al. 2017; Elbaz et al. 2018; Gullberg et al. 2019; Lang et al. 2019; Sun et al. 2021b; Gómez-Guijarro et al. 2022), in contrast to the compact size often seen in local ULIRGs (e.g., ∼0.1 kpc in Arp 220; Soifer et al. 2000; Barcos-Muñoz et al. 2017; Sakamoto et al. 2017).

We further study the redshift evolution of dust temperature in the IR luminosity bins of 10^11.5 − 10¹² L_⊙ (LIRGs) and 10¹²–10^12.5 L_⊙ (ULIRGs). Because 84% of ALCS sources have lensing magnification factors greater than 1.8, our joint sample is ≳80% complete for galaxies in the LIRG bin at z ≲ 2.5 and T_dust ≲ 35 K and galaxies in the ULIRG bin at z ≲ 4 and T_dust ≲ 45 K. The right panel of Figure 15 shows the dust temperatures versus redshifts of galaxies in the two L_IR bins, with the incomplete region on the T_dust–z plane visualized with gray shading. Note that we only include galaxies with uncertainties of ΔT/T less than 20%, which only consist of sources with z_spec or accurate HST z_phot.

To test the existence of any monotonic redshift evolution of T_dust, we compute the Spearman's rank correlation coefficient and p-value for sources within the redshift ranges of high completeness. A weak positive T_dust(z) relation can be tentatively drawn for sources in the L_IR = 10^11.5–10¹² L_⊙ bin (p-value = 0.06), with sources at z ≲ 2 showing lower T_dust than that of local LIRGs (Díaz-Santos et al. 2017). We argue that the selection bias against low-redshift sources with warm dust temperatures (T_dust ≳ 35 K; discussed in Section 6.1) is likely the cause of this weak T_dust(z) relation.

We note that a similar trend is also seen with Symeonidis et al. (2013) for galaxies with L_IR = 10^11.6–10^11.8 L_⊙ at z < 1, and a positive T_dust(z) trend has been suggested by certain stacking analyses (e.g., Magnelli et al. 2014; Schreiber et al. 2018; Simpson et al. 2019) and simulations (e.g., Liang et al. 2019; also Lagos et al. 2020, but a weak T_dust(z) evolution). However, the samples used by most of the stacking analyses are stellar mass selected, which are different in total SFR from L_IR- or M_dust-selected samples obtained with submillimeter/millimeter surveys. This means that the evolution of T_dust(z) at fixed stellar mass can be a combined effect of the weak evolution of the L_IR–T_dust relation and strong evolution of star-forming MS (i.e., increasing SFR/M_star and thus L_IR/M_star toward higher redshifts; see also a recent study by Drew & Casey 2022). The positive T_dust(z) trend is not seen in 450/850 μm selected SMGs in the AS2UDS and STUDIES samples (Dudzevičiūtė et al. 2020, 2021; Lim et al. 2020).

For ULIRGs in the bin of L_IR = 10¹²–10^12.5 L_⊙, no redshift evolution of T_dust can be identified (p-value = 0.71). A similar conclusion can also be drawn with the STUDIES and AS2UDS SMGs with relatively accurate T_dust measurements (uncertainty <15%; typically requires ≥1-band SPIRE detection) in this L_IR bin.

Assuming a canonical gas-to-dust mass ratio (GDR) of δ_GDR = 100, we can estimate a median gas depletion timescale for sources in the ALCS-Herschel joint sample as ${t}_{\mathrm{dep}}={\delta }_{\mathrm{GDR}}{M}_{\mathrm{dust}}/\mathrm{SFR}={190}_{-95}^{+266}$ Myr (error bar denotes 1σ dispersion). Sources with higher T_dust generally show shorter t_dep (equivalently higher star-forming efficiency, SFE =SFR/M_gas). For sources only with FIR z_phot, their dust temperatures are close to the peak of prior T_dust distribution as assumed by magphys, and therefore the t_dep is nearly constant. Despite such a degeneracy, the derived median t_dep does not change significantly if we exclude sources without z_spec (166 ± 25 Myr) or HST z_phot (219 ± 30 Myr) from the statistics. We further note that t_dep is a conserved quantity independent of the lensing magnification effect.

The derived t_dep of our sample is broadly consistent with those derived with the AS2UDS sample (${295}_{-185}^{+307}$ Myr; Dudzevičiūtė et al. 2020), ASPECS main sample (${299}_{-167}^{+431}$ Myr; Aravena et al. 2020) and HLS bright lensed source sample (${226}_{-73}^{+196}$ Myr; Sun et al. 2021b) using the same method (i.e., δ_GDR = 100 with M_dust derived from magphys SED modeling). Although the requirement of Herschel detections will result in a selection bias toward sources with higher T_dust and shorter t_dep, the use of Herschel data provides characterization of obscured SFR and t_dep in higher precision.

Figure 16 shows the gas depletion timescale as a function of redshift for sources in the ALCS-Herschel joint sample. We also compare ALCS sources with the ASPECS main sample (Aravena et al. 2020), AS2UDS (Dudzevičiūtė et al. 2020), GOODS-ALMA (Franco et al. 2020), and PHIBSS (Tacconi et al. 2018) across comparable redshift ranges. The best-fit models of t_dep(z) at M_star = 10^10.6 M_⊙ and MS SFR (Speagle et al. 2014), based on the prescriptions of Scoville et al. (2017), Tacconi et al. (2018), and Liu et al. (2019), are also plotted for comparison. Sources in our sample exhibit a shorter depletion timescale than the model predictions. This is caused by the IR-selection nature (i.e., favoring sources above the star-forming MS with higher L_IR and T_dust) and potentially low-T_dust assumptions for single-band millimeter continuum sources in previous studies (e.g., 25 K in Scoville et al. 2017). Additionally, δ_GDR is observed to be a function of metallicity (e.g., Leroy et al. 2011; δ_GDR ∼ 200 at half-solar metallicity), which could introduce further uncertainty to M_gas and t_dep.

Zoom In Zoom Out Reset image size

Although a declining trend of t_dep toward higher redshift can be identified, we argue that this is a selection effect of sources with higher L_IR and T_dust toward higher redshift. If we restrict the sample to an intrinsic IR luminosity of 10¹²–10^12.5 L_⊙ at z ≃ 1 − 4, where the ALMA and Herschel/SPIRE detections are ≳80% complete (Section 6.4 and the lower right panel of Figure 15), we find no significant redshift dependence of t_dep (the p-value of Spearman's rank correlation is 0.21). This is consistent with the conclusion of Dudzevičiūtė et al. (2020) drawn upon the L_IR-limited sample with both ALMA and SPIRE detections.

R0600-ID164 is a [C ii]-emitting lensed arc at z = 6.072 discovered blindly by the ALCS (named RXC J0600–2007 z6.1/z6.2 in Fujimoto et al. 2021). Laporte et al. (2021) reported the ALMA dust continuum detection of this source (S/N = 4.84 on the 2''-tapered map). Although this source is not detected in any Herschel band, the upper limits of flux densities at 250–500 μm can be used to constrain its dust temperature.

Figure 17 displays the observed FIR SED of R0600-ID164 before lensing magnification correction (red open circles; upper limits are at 3σ). We compare the SED with MBB spectral templates at intrinsic T_dust = 50 − 125 K with the CMB effects considered following the prescription of da Cunha et al. (2013). Based on the dust mass, dust continuum size reported in Laporte et al. (2021), and dust absorption coefficient adopted in Section 4.1, we find a low dust mass surface density of (6 ± 3) × 10⁶ M_⊙ kpc⁻². Therefore, the optical depth of dust continuum is on the order of unity at λ_thick ≲ 30 μm in the rest frame. The dust emissivities are assumed as β = 1.6 and 2.0. We also compare R0600-ID164 with galaxies at z ≃ 6 − 7 with normal dust temperatures (40–50 K; shown as green circles) confirmed with ALMA continuum detections in two bands at least, including J1211–0118, J0217–0208 (z = 6.03 and 6.20, respectively; Harikane et al. 2020), A1689-zD1 (z =7.13; Watson et al. 2015; Knudsen et al. 2017; Inoue et al. 2020; Bakx et al. 2021), and B14–65666 (z = 7.15; Hashimoto et al. 2019; Sugahara et al. 2021).

Zoom In Zoom Out Reset image size

Bakx et al. (2020) showed that MACS0416 Y1, a lensed galaxy at z = 8.31 (black open squares in Figure 17; also Tamura et al. 2019), exhibits abnormally warm dust temperature (T_dust >80 K, 90% confidence) and/or steep dust emissivity index (β >2). With the three-band SPIRE flux density upper limits, we calculate the χ² of nondetections assuming a T_dust = 80 K, β = 2.0, and λ_thick = 30 μm dust continuum model. The derived χ² is 6.6 (reduced ${\chi }_{\nu }^{2}=2.2$), indicating that such a model can be ruled out at >90% confidence level. Although assuming lower dust emissivity and longer λ_thick will lead to a smaller χ², we can draw the conclusion that the T_dust and β of R0600-ID164 are not as extreme as those observed for MACS0416 Y1. However, we also note that the differential lensing effect seen around the caustic line (μ = 30 − 160; Fujimoto et al. 2021) may introduce further uncertainty on the dust temperature constraint.