# ALMA Lensing Cluster Survey: Deep 1.2 mm Number Counts and Infrared Luminosity Functions at $z \simeq 1-8$

Seiji Fujimoto, Kotaro Kohno, Masami Ouchi, Masamune Oguri, Vasily Kokorev, Gabriel Brammer, Fengwu Sun, Jorge Gonzalez-Lopez, Franz E Bauer, Gabriel B Caminha, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Seiji Fujimoto, Kotaro Kohno, Masami Ouchi, Masamune Oguri, Vasily Kokorev, et al. ALMA Lensing Cluster Survey: Deep 1.2 mm Number Counts and Infrared Luminosity Functions at $z \simeq 1-8$. The Astrophysical Journal Supplement, 2024, 275 (2), pp.36. 10.3847/1538-4365/ad5ae2. hal-04042755

HAL Id: hal-04042755

https://hal.science/hal-04042755v1

Submitted on 17 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# ALMA Lensing Cluster Survey: Deep 1.2 mm Number Counts and Infrared Luminosity Functions at $z \simeq 1-8$

Seiji Fujimoto<sup>1,2,3,47</sup> (D), Kotaro Kohno<sup>4,5</sup> (D), Masami Ouchi<sup>6,7,8</sup> (D), Masamune Oguri<sup>9,10</sup> (D), Vasily Kokorev<sup>2,3,11</sup> (D), Gabriel Brammer<sup>2,3</sup> (D), Fengwu Sun<sup>12</sup> (D), Jorge González-López<sup>13,14</sup> (D), Franz E. Bauer<sup>15,16,17</sup> (D), Gabriel B. Caminha<sup>11</sup> (D), Bunyo Hatsukade<sup>4,6,18</sup> (D), Johan Richard<sup>19</sup> (D), Ian Smail<sup>20</sup> (D), Akiyoshi Tsujita<sup>4</sup> (D), Yoshihiro Ueda<sup>21</sup> (D), Ryosuke Uematsu<sup>21</sup> (D), Adi Zitrin<sup>22</sup> (D), Dan Coe<sup>23,24,25</sup> (D), Jean-Paul Kneib<sup>26</sup> (D), Marc Postman<sup>23</sup> (D), Keiichi Umetsu<sup>27</sup> (D), Claudia del P. Lagos<sup>2,28,29</sup> (D), Gergö Popping<sup>30</sup> (D), Yiping Ao<sup>31</sup> (D), Larry Bradley<sup>23</sup> (D), Karina Caputi<sup>2,11</sup> (D), Miroslava Dessauges-Zavadsky<sup>32</sup> (D), Eiichi Egami<sup>12</sup> (D), Daniel Espada<sup>33,34</sup> (D), R. J. Ivison<sup>35</sup> (D), Mathilde Jauzac<sup>20,36</sup> (D), Kirsten K. Knudsen<sup>37</sup> (D), Anton M. Koekemoer<sup>23</sup> (D), Georgios E. Magdis<sup>2,3,38</sup> (D), Guillaume Mahler<sup>20,36</sup> (D), A. M. Muñoz Arancibia<sup>17,39</sup> (D), Timothy Rawle<sup>40</sup> (D), Kazuhiro Shimasaku<sup>5,41</sup> (D), Sune Toft<sup>2,3</sup> (D), Hideki Umehata<sup>42,43,44</sup> (D), Francesco Valentino<sup>2,3,30</sup> (D), Tao Wang<sup>45,46</sup> (D), and Wei-Hao Wang<sup>27</sup> (D)

<sup>1</sup>Department of Astronomy, The University of Texas at Austin, Austin, TX, USA; fujimoto@utexas.edu

<sup>2</sup>Cosmic Dawn Center (DAWN), Denmark

<sup>3</sup>Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK2200 Copenhagen N, Denmark

<sup>4</sup>Institute of Astronomy, Graduate School of Science, The University of Tokyo, 2-21-1 Osawa, Mitaka, Tokyo 181-0015, Japan

<sup>5</sup>Research Center for the Early Universe, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

<sup>6</sup>National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

<sup>7</sup>Institute for Cosmic Ray Research, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8582, Japan

<sup>8</sup>Kavli Institute for the Physics and Mathematics of the Universe (WPI), University of Tokyo, Kashiwa, Chiba 277-8583, Japan

<sup>9</sup>Center for Frontier Science, Chiba University, 1-33 Yayoi-cho, Inage-ku, Chiba 263-8522, Japan

<sup>10</sup>Department of Physics, Graduate School of Science, Chiba University, 1-33 Yayoi-Cho, Inage-Ku, Chiba 263-8522, Japan

<sup>11</sup>Kapteyn Astronomical Institute, University of Groningen, Postbus 800, 9700 AV Groningen, The Netherlands

<sup>12</sup>Steward Observatory, University of Arizona, 933 N. Cherry Ave., Tucson, AZ 85721, USA

<sup>13</sup>Núcleo de Astronomía de la Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército Libertador 441, Santiago, Chile

<sup>14</sup>Las Campanas Observatory, Carnegie Institution of Washington, Casilla 601, La Serena, Chile

<sup>15</sup>Instituto de Astrofísica, Facultad de Física, Pontificia Universidad Católica de Chile, Campus San Joaquín, Av. Vicuña Mackenna 4860, Macul Santiago, 7820436, Chile

<sup>16</sup>Centro de Astroingeniería, Facultad de Física, Pontificia Universidad Católica de Chile, Campus San Joaquín, Av. Vicuña Mackenna 4860, Macul Santiago, 7820436, Chile

<sup>17</sup>Millennium Institute of Astrophysics, Nuncio Monseñor Sótero Sanz 100, Of 104, Providencia, Santiago, Chile

<sup>18</sup>Graduate Institute for Advanced Studies, SOKENDAI, Osawa, Mitaka, Tokyo 181-8588, Japan

<sup>19</sup>Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230, Saint-Genis-Laval, France

<sup>20</sup>Centre for Extragalactic Astronomy, Department of Physics, Durham University, Durham DH1 3LE, UK

<sup>21</sup>Department of Astronomy, Kyoto University, Kyoto 606-8502, Japan

<sup>22</sup>Physics Department, Ben-Gurion University of the Negev, P.O. Box 653, Be'er-sheva 8410501, Israel

<sup>23</sup>Space Telescope Science Institute (STScI), 3700 San Martin Dr., Baltimore, MD 21218, USA

<sup>24</sup>Association of Universities for Research in Astronomy (AURA), Inc. for the European Space Agency (ESA), USA

<sup>25</sup>Center for Astrophysical Sciences, Department of Physics and Astronomy, The Johns Hopkins University, 3400 N Charles St., Baltimore, MD 21218, USA

<sup>26</sup>Laboratoire d'Astrophysique, Ecole Polytechnique Fédérale de Lausanne, Observatoire de Sauverny, CH-1290 Versoix, Switzerland

<sup>27</sup>Institute of Astronomy and Astrophysics, Academia Sinica (ASIAA), AS/NTU Astronomy-Mathematics Building, No. 1, Sec. 4, Roosevelt Rd., Taipei 10617, Taiwan

<sup>28</sup>International Centre for Radio Astronomy Research (ICRAR), M468, University of Western Australia, 35 Stirling Hwy, Crawley, WA 6009, Australia

<sup>29</sup>ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia

<sup>30</sup>European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748, Garching, Germany

<sup>31</sup>Purple Mountain Observatory and Key Laboratory for Radio Astronomy, Chinese Academy of Sciences, Nanjing, People's Republic of China

<sup>32</sup>Département d'Astronomie, Université de Genève, Chemin Pegasi 51, 1290 Versoix, Switzerland

<sup>33</sup>Departamento de Física Teórica y del Cosmos, Campus de Fuentenueva, Edificio Mecenas, Universidad de Granada, E-18071, Granada, Spain

<sup>34</sup>Instituto Carlos I de Física Teórica y Computacional, Facultad de Ciencias, E-18071, Granada, Spain

<sup>35</sup>European Southern Observatory, Karl-Schwarzschild-Strasse 2, D-85748 Garching, Germany

<sup>36</sup>Institute for Computational Cosmology, Durham University, South Road, Durham DH1 3LE, UK

<sup>37</sup>Department of Space, Earth and Environment, Chalmers University of Technology, SE-412 96 Gothenburg, Sweden

<sup>38</sup>DTU-Space, Technical University of Denmark, Elektrovej 327, 2800, Kgs. Lyngby, Denmark

<sup>39</sup>Center for Mathematical Modeling (CMM), Universidad de Chile, Beauchef 851, Santiago, Chile

<sup>40</sup>European Space Agency (ESA), ESA Office, Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

<sup>41</sup>Department of Astronomy, School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

<sup>42</sup>Institute for Advanced Research, Nagoya University, Furocho, Chikusa, Nagoya 464-8602, Japan

<sup>43</sup>Department of Physics, Graduate School of Science, Nagoya University, Furocho, Chikusa, Nagoya 464-8602, Japan

<sup>44</sup>Cahill Center for Astronomy and Astrophysics, California Institute of Technology, 1200 E California Blvd, MC 249-17, Pasadena, CA 91125, USA

<sup>45</sup>School of Astronomy and Space Science, Nanjing University, Nanjing 210093, People's Republic of China

<sup>46</sup>Key Laboratory of Modern Astronomy and Astrophysics (Nanjing University), Ministry of Education, Nanjing 210093, People's Republic of China

*Received 2023 March 3; revised 2024 June 5; accepted 2024 June 20; published 2024 November 27*

<figure><img src="image_4.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further

distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="47"><sup>47</sup>Hubble Fellow.</footnote>
Abstract

We present a statistical study of 180 dust continuum sources identified in 33 massive cluster fields by the Atacama Large Millimeter/submillimeter Array Lensing Cluster Survey (ALCS) over a total of 133 arcmin² area, homogeneously observed at 1.2 mm. ALCS enables us to detect extremely faint millimeter sources by lensing magnification, including near-infrared (NIR) dark objects showing no counterparts in existing Hubble Space Telescope and Spitzer images. The dust continuum sources belong to a blind sample (N = 141) with signal-to-noise ratio (S/N) ≳ 5.0 (a purity of >0.99) or a secondary sample (N = 39) with S/N = 4.0–5.0 screened by priors. With the blind sample, we securely derive 1.2 mm number counts down to ~7 µJy, and find that the total integrated 1.2 mm flux is $20.7_{-6.5}^{+8.5}$ Jy deg$^{-2}$, resolving ≈80% of the cosmic infrared background light. The resolved fraction varies by a factor of 0.6–1.1 due to the completeness correction depending on the spatial size of the millimeter emission. We also derive infrared (IR) luminosity functions (LFs) at $z = 0.6–7.5$ with the $1/V_{max}$ method, finding the redshift evolution of IR LFs characterized by positive luminosity and negative density evolution. The total (= UV + IR) cosmic star formation rate density (SFRD) at $z > 4$ is estimated to be $161_{-21}^{+25}$% of the Madau and Dickinson measurements mostly based on rest-frame UV surveys. Although our general understanding of the cosmic SFRD is unlikely to change beyond a factor of 2, these results add to the weight of evidence for an additional (≈60%) SFRD component contributed by the faint millimeter population, including NIR-dark objects.

*Unified Astronomy Thesaurus concepts: Galaxy formation (595); Galaxy evolution (594); Starburst galaxies (1570); Millimeter astronomy (1061); Cosmic background radiation (317); Luminosity function (942); Interferometry (808)*

1. Introduction

Extragalactic background light is electromagnetic radiation that fills the Universe, arising from several distinct physical processes and energetically dominated by cosmic microwave background (CMB), cosmic optical background (COB), and cosmic infrared background (CIB) radiation. The CMB and COB are explained by the leftover radiation from the Big Bang and the stellar continuum of galaxies, respectively (e.g., Totani et al. 2001; Planck Collaboration et al. 2014b; cf. Lauer et al. 2022). On the other hand, the origin of the CIB has not yet been fully accounted for as yet, despite its importance implied by the fact that the total energy of the CIB has been known to be comparable to that of the COB since its initial discovery with the Cosmic Background Explorer (COBE) satellite (Puget et al. 1996; Fixsen et al. 1998; Hauser et al. 1998; Hauser & Dwek 2001; Dole et al. 2006). Therefore, concluding whether the CIB is explained by known physical mechanisms and objects in the Universe is an essential open question in modern astrophysics.

Some fraction of the stellar continuum in galaxies is absorbed by dust, reemitted as thermal infrared (IR) emission, and thought to contribute to the CIB. Thus, to understand the origin of the CIB, the first approach is detecting individual IR-emitting galaxies and evaluating their contribution to the CIB. In these studies, one of the most advantageous wavelength regimes is the submillimeter/millimeter, owing to the negative k-correction. Rare IR-bright, high-redshift dusty star-forming galaxies were first identified a few decades ago at submillimeter/millimeter wavelengths, referred to as submillimeter galaxies (SMGs) due to their brightness at the submillimeter/millimeter wavelengths ($S_{1mm} \gtrsim 1$ mJy; e.g., Hughes et al. 1998; Blain et al. 2002; Casey et al. 2014). However, bright SMGs only contribute ~10%–20% to the CIB, based on previous blank field observations with single-dish telescopes (e.g., Perera et al. 2008; Hatsukade et al. 2011; Scott et al. 2012), indicating that the bulk of the CIB is produced by populations distinct from the SMGs. While several previous studies using the stacking analyses with Herschel and SCUBA2 data have pushed the detection limit beyond the SMGs and hinted that faint dusty galaxies may dominantly account for the

CIB at submillimeter wavelengths (e.g., Béthermin et al. 2012; Viero et al. 2013; Zavala et al. 2017), direct individual detection is essential to avoid potential systematics inherent in the stacking approach. Moreover, recent studies suggest the presence of remarkably extended diffuse dust emission beyond individual galaxy scales (e.g., Xiao et al. 2023a), which is challenging to distinguish given the poor spatial resolution of these instruments (~10"–30"). Consequently, deep submillimeter/millimeter observations with high angular resolution are crucial for directly verifying the origin of the CIB, whether it arises from compact emission of individual galaxies or diffuse emission around galaxies, for example.

Atacama Large Millimeter/submillimeter Array (ALMA) observations allow us to explore a faint submillimeter/millimeter regime ($S_{1mm} < 1$ mJy) without uncertainties from source confusion and blending, owing to ALMA's high sensitivity and angular resolution. Faint submillimeter/millimeter sources have been investigated in a large variety of ALMA observations: via dedicated blind surveys of blank (e.g., Aravena et al. 2016; Wang et al. 2016; Dunlop et al. 2017; Umehata et al. 2017; Franco et al. 2018; Hatsukade et al. 2018; González-López et al. 2020; Zavala et al. 2021; Adscheid et al. 2024) and lensing cluster fields (e.g., Fujimoto et al. 2016; Muñoz Arancibia et al. 2019, 2023); and as serendipitous sources in single pointing observations (e.g., Hatsukade et al. 2013; Ono et al. 2014; Carniani et al. 2015; Fujimoto et al. 2016; Zavala et al. 2018; Béthermin et al. 2020) and calibrator fields (e.g., Oteo et al. 2016; Klitsch et al. 2020). These studies show that faint submillimeter/millimeter sources newly identified with ALMA strongly outnumber the SMGs. For instance, the deepest observations, aided by strong gravitational lensing, report that ~70%–100% of CIB is resolved down to ~0.01 mJy at ≈1 mm (e.g., Fujimoto et al. 2016; Muñoz Arancibia et al. 2019). On the other hand, the deepest ALMA blank field surveys so far in Hubble Ultra Deep Field (HUDF) estimate a 1.2 mm total flux density integrated down to zero flux by extrapolation to be $6.8 \pm 0.4$ Jy deg$^{-2}$ (González-López et al. 2020), which corresponds to ~40% of the CIB at 1.2 mm. Because the effective survey areas in these previous deepest ALMA studies are less than a few square arcminutes at the faintest 1 mm flux density regime of 0.01–0.1 mJy, the difference might result from
cosmic variance (e.g., Trenti & Stiavelli 2008), and concluding whether the CIB is fully resolved or not has still been challenging even with ALMA. Naturally, this also advocates the difficulty in constraining the IR luminosity function (LF) and subsequently dust-obscured cosmic star formation rate density (SFRD), especially at high redshifts ($z \gtrsim 4$; e.g., Novak et al. 2017; Liu et al. 2018; Gruppioni et al. 2020).

To enlarge the survey area at the faintest regime and overcome cosmic variance, making full use of gravitational lensing is one of the most effective approaches. Assuming a magnification factor of $\mu$ and a given survey area in the observed frame of A, the effective survey area and the data depth after the lens correction decrease by a factor of $\mu$ from the observed frame. We thus need $\mu$ times observation time to recover the given survey area of $A (=A/\mu \times \mu)$. On the other hand, if we achieve the $\mu$ times deep data at the area of A without the lensing support, we generally need $\mu^2$ times integration time for the observations. Therefore, the gravitational lensing support saves in observation time by a factor of $\sim\mu$ to achieve observations at a given survey area and data depth. Following the successful outcome from the SCUBA-2 lensing cluster survey (e.g., Cowie et al. 2022), a complete ALMA survey toward massive galaxy clusters will allow us to detect large numbers of faint submillimeter/millimeter sources even out to the epoch of reionization (Watson et al. 2015; Laporte et al. 2017; Inami et al. 2022) in the most effective way and offers a unique opportunity to assess the origin of the CIB.

In this paper, we present statistics of 180 faint 1.2 mm sources identified from the large imaging and spectral campaign of 33 massive galaxy clusters in the ALMA Lensing Cluster Survey (ALCS) to assess the origins of the CIB. Making full use of the rich ancillary data sets, including Hubble Space Telescope (HST) and Spitzer IRAC images and photometry, we also study the IR LF at $z=1-8$ and the cosmic SFRD based on the largest statistics so far for the faint submillimeter/millimeter population with ALMA. The structure of this paper is as follows. In Section 2, we overview the survey of ALCS and the data sets. Section 3 outlines the methods of the source extraction, flux measurements, redshift estimates, magnification estimates, corrections for the flux measurements and completeness through Monte Carlo (MC) simulations, and survey area estimates. We characterize our ALCS sources in Section 4. In Section 5, we present the 1.2 mm number counts and the contribution to the CIB. The redshift evolutions of IR LFs and the cosmic SFRD are discussed in Section 6. A summary is presented in Section 7. Throughout this paper, we assume a flat universe with $\Omega_m = 0.3$, $\Omega_\Lambda = 0.7$, and $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$. We adopt a search radius of $1''0$ in crossmatching catalogs unless otherwise specified.48 We take the CMB effect in submillimeter/millimeter observations into account (e.g., da Cunha et al. 2013; Pallottini et al. 2015; Zhang et al. 2016; Lagache et al. 2018), following the recipe presented by da Cunha et al. (2013).

# 2. Observations and Data Processing

## 2.1. Survey Design

ALCS is a cycle-6 ALMA large program (project ID: 2018.1.00035.L; PI: K. Kohno) to map high-magnification

regions in 33 massive galaxy clusters at 1.2 mm in Band 6. The sample is selected from the best-studied clusters drawn from HST treasury programs, i.e., Hubble Frontier Fields (HFFs; Lotz et al. 2017), the Cluster Lensing and Supernova Survey with Hubble (CLASH; Postman et al. 2012), and the Reionization Lensing Cluster Survey (RELICS; Coe et al. 2019). Observations were carried out between 2018 December and 2019 December in compact array configurations of C43-1 and C43-2, where 26 and seven clusters were observed in Cycles 6 and 7, respectively. We adopt two frequency setups with the sky frequencies of 259.4 and 263.2 GHz, accomplishing a 15 GHz wide spectral scan in ranges of 250.0–257.5 GHz and 265.0–272.5 GHz to enlarge the survey volume for line-emitting galaxies. For the five HFF clusters visible with ALMA, ALCS only performed observations in one of the two frequency setups, because observations for the other frequency setup had been taken in previous ALMA HFF surveys (2013.1.00999.S and 2015.1.01425.S; e.g., González-López et al. 2017b). The observations were carried out in mosaic mode, covering highly magnified regions spanning $\sim 1–9 \text{ arcmin}^2$ sky areas per cluster. The survey description is also presented in Kohno et al. (2023).

## 2.2. Data Reduction, Calibration, and Imaging

The ALMA data were reduced and calibrated with the Common Astronomy Software Applications (CASA) package versions 5.4.0 and 5.6.1 (McMullin et al. 2007) for the data taken in Cycles 6 and 7, respectively, with the pipeline script in the standard manner. For the HFF data taken in the previous surveys, we used CASA versions 4.2.2–4.5.3 in order to use the pipeline scripts from the previous cycles. We created several data products (measurement sets, maps, and cubes) for every ALCS cluster. We imaged the calibrated visibilities with natural weighting, a pixel scale of $0.′′16$, and a primary beam limit of 0.2 with the CASA task TCLEAN. For continuum maps, the TCLEAN routines were executed down to the $2\sigma$ level with a maximum iteration number of 100,000 in the automask mode.49 For cubes, we adopted common spectral channel bins of 30 and $60 \text{ km s}^{-1}$ and created the cubes without the CLEAN iteration. We did not identify bright emission per channel in 31 cluster cubes, where we use the cubes in the following analysis. In the other two clusters of MACSJ0553 and MACS1931, we identified significant signals in each channel due to bright line emitters serendipitously detected in the cubes (Section 3.4). We thus performed the CLEAN algorithm for these two cluster cubes. When we found systematic stripe patterns in the products by visual inspection, we applied additional flaggings50 and/or performed TCLEAN using manual mode masking.51 The natural-weighted maps achieved full width at half-maximum (FWHM) sizes of the synthesized beam between $0.′′94 \times 0.′′65–1.′′57 \times 1.′′01$, with $1\sigma$ sensitivities of $46.9–91.6 \text{ µJy beam}^{-1}$ for the continuum and $848.1–1706.4 \text{ µJy beam}^{-1}$ for the $60 \text{ km s}^{-1}$ width channel cube. We show a reduced ALMA map in Figure 1. We

<footnote n="48">48 While Gómez-Guijarro et al. (2022) suggest that $\le 0.′′4$ is the robust counterpart search radius, we adopt the relatively large search radius as a default because the spatial offsets among the different wavelengths can be enhanced in highly magnified sources.</footnote>

<footnote n="49">49 We determined subparameters of tclean in the automask mode based on recommendations from the ALMA automasking guide: https://casaguides.nrao.edu/index.php/Automasking_Guide.</footnote>

<footnote n="50">50 Field ID = 18 in uid_A002_Xd98580_X44f4.ms was flagged before imaging, because the field had been observed for a calibration purpose, but wrongly named as a target field.</footnote>

<footnote n="51">51 We invoked manual mode in TCLEAN for the cubes of MACSJ0553 and MACS1931 by setting rectangle masks ($\sim 10′′ \times 10′′$) around the bright line emitters.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 1. Footprints of ALMA Band 6 (left), HST/F160W (middle), and IRAC/ch1 (right) in one of the ALCS fields. The dashed and solid cyan lines show the relative sensitivity response to the deepest 30% and 50% of the mosaic, respectively. The white lines denote the $\mu = 200$ magnification curve at $z = 2$, estimated from our fiducial lens model of this cluster. The labeled solid and dashed green circles represent the ALMA continuum catalog sources with $S/N_{nat} \ge 5$ and $4-5$, respectively.</figcaption></figure>

summarize the data properties of the continuum map and the cube in Table 1. We also produced lower resolution maps and cubes by applying a uv-taper parameter ($2'' \times 2''$) to recover spatially extended, low-surface brightness emission associated with large, resolved galaxies or due to gravitational lensing effects. We refer to our ALMA maps (cubes) without and with the uv-taper as natural and tapered maps (cubes), respectively. The ALCS products, including the uv-tapered maps and cubes, are publicly available via the dedicated page<fnref n="52" /> and ALMA science data archive.

## 3. Data Analysis

### 3.1. Source Extraction

We conduct a blind 2D source extraction with SExTRACTOR version 2.5.0 (Bertin & Arnouts 1996). We use the natural and tapered maps before primary beam correction for the source extraction and apply the correction to the flux measurements in the following analysis. We extract sources with positive peak counts above $2.5\sigma$ confidence, referred to as a peak catalog. Because the peak catalog should include a large number of spurious sources, after confirming that the noise remains Gaussian and centered at zero flux, we also conduct a negative peak analysis (e.g., Hatsukade et al. 2013; Ono et al. 2014; Carniani et al. 2015; Fujimoto et al. 2016) to evaluate the expected number of real sources. We create inverted maps by multiplying $-1$ to both the natural and tapered maps, run SExTRACTOR again, and extract sources with negative peak counts above $2.5\sigma$ level. In both original and inverted maps, we select sources that are identified in high-sensitivity regions, where the relative sensitivity to the deepest part of the mosaic is greater than 30%.

In Figure 2, we show histograms of the positive and negative sources in our ALMA maps from 33 ALCS fields as a function of the signal-to-noise ratio (S/N) at the peak count. Note that

this differs from the histogram of all pixels in the map, as an island of the emission is counted as a single source. In these histograms, the excess of the positive to the negative sources indicates the expected number of the real sources, known as the purity (e.g., González-López et al. 2020; Gómez-Guijarro et al. 2022), defined as

$$p = \frac{N_{\text{pos}} - N_{\text{neg}}}{N_{\text{pos}}}, \qquad (1)$$

where $N_{\text{pos}}$ and $N_{\text{neg}}$ represent the number of positive and negative sources at a given S/N, respectively. We find that $p$ shows $> 0.99$ at S/N $\ge 5.0$ and $\ge 4.5$ in the natural (S/N$_{nat}$) and tapered maps (S/N$_{tap}$), respectively. The difference is explained by less areas relative to the beam size in the tapered maps, where the significant count caused by the noise fluctuation should be reduced. We also find that $p$ maintains an excess above zero down to S/N$_{nat} = 4.0$, indicating that a certain number of sources are likely real down to S/N$_{nat} = 4.0$. Based on these procedures, we obtain a source candidate catalog consisting of 399 objects with S/N$_{nat} \ge 4.0$ from the peak catalog. In this candidate catalog, the number of the real sources down to S/N$_{nat} = 4.0$ is expected to be $177_{-14}^{+12}$ (where the uncertainty is calculated from Poisson statistics) by summing the excess of the positive to the negative sources down to the S/N limit.

### 3.2. Catalog

We analyze the source candidate catalog consisting of the 399 objects to construct a reliable source catalog. Because the positive and negative source histograms of Figure 2 show a rapid decrease of purity below S/N$_{nat} = 5.0$ and S/N$_{tap} = 4.5$, we make two catalogs from blind and prior-based approaches. First, in the blind approach, we adopt thresholds of S/N$_{nat} \ge 5.0$ or S/N$_{tap} \ge 4.5$ in the source candidate catalog, yielding 141 ALMA sources that achieve a purity of $p > 0.99$. We refer to these 141 sources as the blind catalog. Second, in the prior-

<footnote n="52">⁵² http://www.ioa.s.u-tokyo.ac.jp/ALCS/</footnote>
<table><caption>**Table 1**
ALCS Data Properties for 33 Lensing Clusters</caption><thead><tr><th>Cluster Name<br>(1)</th><th>Area<br>(arcmin<sup>2</sup>)<br>(2)</th><th>Beam<br>(arcsec)<br>(3)</th><th>σ (cont)<br>(μJy beam<sup>-1</sup>)<br>(4)</th><th>σ (cube, tune1)<br>(μJy beam<sup>-1</sup>)<br>(5)</th><th>σ (cube, tune2)<br>(μJy beam<sup>-1</sup>)<br>(6)</th><th>N<br>(7)</th></tr></thead><tbody><tr><td colspan="7" style="text-align: center; font-weight: bold;">HFF</td></tr><tr><td>A370</td><td>5.44</td><td>1.12 × 0.89</td><td>47</td><td>673</td><td>848</td><td>7</td></tr><tr><td>A2744</td><td>5.13</td><td>0.94 × 0.65</td><td>51</td><td>839</td><td>643</td><td>11</td></tr><tr><td>A1063</td><td>5.01</td><td>1.15 × 0.96</td><td>53</td><td>765</td><td>974</td><td>4</td></tr><tr><td>M0416.1</td><td>5.24</td><td>1.48 × 0.85</td><td>55</td><td>879</td><td>1158</td><td>7</td></tr><tr><td>M1149.5</td><td>5.52</td><td>1.24 × 1.1</td><td>64</td><td>778</td><td>975</td><td>4</td></tr><tr><td colspan="7" style="text-align: center; font-weight: bold;">CLASH</td></tr><tr><td>A209</td><td>1.20</td><td>1.31 × 1.07</td><td>63</td><td>946</td><td>1193</td><td>1</td></tr><tr><td>A383</td><td>1.39</td><td>1.13 × 0.87</td><td>61</td><td>989</td><td>1245</td><td>3</td></tr><tr><td>M0329</td><td>3.03</td><td>1.57 × 1.01</td><td>71</td><td>1083</td><td>1354</td><td>1</td></tr><tr><td>M0429</td><td>1.27</td><td>1.58 × 0.96</td><td>92</td><td>1397</td><td>1706</td><td>3</td></tr><tr><td>M1115</td><td>1.67</td><td>1.16 × 1.06</td><td>63</td><td>925</td><td>1166</td><td>5</td></tr><tr><td>M1206</td><td>2.87</td><td>1.13 × 0.98</td><td>53</td><td>789</td><td>1008</td><td>8</td></tr><tr><td>M1311</td><td>1.38</td><td>1.10 × 0.98</td><td>62</td><td>906</td><td>1159</td><td>2</td></tr><tr><td>M1423</td><td>1.88</td><td>1.34 × 1.07</td><td>65</td><td>951</td><td>1266</td><td>4</td></tr><tr><td>M1931</td><td>2.64</td><td>1.34 × 1.03</td><td>56</td><td>817</td><td>1023</td><td>6</td></tr><tr><td>M2129</td><td>2.54</td><td>1.23 × 0.95</td><td>47</td><td>773</td><td>940</td><td>3</td></tr><tr><td>R2129</td><td>1.01</td><td>1.26 × 1.04</td><td>40</td><td>720</td><td>876</td><td>2</td></tr><tr><td>R1347</td><td>3.57</td><td>1.15 × 1.01</td><td>53</td><td>795</td><td>1037</td><td>6</td></tr><tr><td colspan="7" style="text-align: center; font-weight: bold;">RELICS</td></tr><tr><td>A3192</td><td>5.17</td><td>1.44 × 0.96</td><td>73</td><td>1125</td><td>1472</td><td>7</td></tr><tr><td>A2163</td><td>2.14</td><td>1.23 × 1.07</td><td>50</td><td>751</td><td>953</td><td>1</td></tr><tr><td>A2537</td><td>3.03</td><td>1.34 × 1.09</td><td>69</td><td>971</td><td>1262</td><td>5</td></tr><tr><td>A295</td><td>4.63</td><td>0.90 × 0.81</td><td>74</td><td>1219</td><td>1481</td><td>2</td></tr><tr><td>AC0102</td><td>6.02</td><td>1.17 × 0.91</td><td>72</td><td>1048</td><td>1336</td><td>14</td></tr><tr><td>M0035</td><td>3.38</td><td>1.42 × 1.03</td><td>52</td><td>846</td><td>1016</td><td>4</td></tr><tr><td>M0159</td><td>3.25</td><td>1.38 × 1.08</td><td>63</td><td>942</td><td>1192</td><td>4</td></tr><tr><td>M0257</td><td>2.54</td><td>1.46 × 0.94</td><td>83</td><td>1238</td><td>1587</td><td>1</td></tr><tr><td>M0417</td><td>6.54</td><td>1.47 × 0.95</td><td>84</td><td>1248</td><td>1566</td><td>8</td></tr><tr><td>M0553</td><td>9.28</td><td>1.21 × 0.92</td><td>65</td><td>1057</td><td>1275</td><td>14</td></tr><tr><td>P171</td><td>5.40</td><td>1.32 × 0.99</td><td>73</td><td>1231</td><td>1595</td><td>4</td></tr><tr><td>RJ2211</td><td>6.83</td><td>1.33 × 1.07</td><td>78</td><td>1078</td><td>1389</td><td>4</td></tr><tr><td>R0032</td><td>8.17</td><td>1.17 × 1.13</td><td>71</td><td>1030</td><td>1306</td><td>20</td></tr><tr><td>R0600</td><td>7.91</td><td>1.22 × 0.95</td><td>57</td><td>822</td><td>1043</td><td>5</td></tr><tr><td>R0949</td><td>3.62</td><td>1.22 × 1.18</td><td>62</td><td>903</td><td>1144</td><td>7</td></tr><tr><td>SM0723</td><td>2.52</td><td>1.01 × 0.77</td><td>66</td><td>969</td><td>1242</td><td>4</td></tr><tr><td>AVERAGE</td><td>3.98</td><td>1.26 × 0.98</td><td>63</td><td>955</td><td>1195</td><td>5.5</td></tr><tr><td>TOTAL</td><td>131.22</td><td>...</td><td>...</td><td>...</td><td>...</td><td>180</td></tr></tbody></table>

**Note.** Column (1): Names of the 33 ALCS target clusters. We shorten the prefixes of Abell(S), MACS(J), PCLKG, RX(C)J, and SMACSJ to A, M, P, R, and SM, respectively. The numbers of declinations in the names are omitted. Column (2): sky area observed in ALCS within the relative sensitivity above 30% to the deepest part of the mosaic. Column (3): FWHM of the synthesized beam in the natural-weighted map. Column (4): data depth evaluated by the standard deviation of the pixel with the CASA task IMSTAT with the by weight algorithm. Column (5): data depth evaluated in the same manner as column (4), but for the 60 km s⁻¹ width cube obtained in the frequency setup with the central sky frequency at 259.4 GHz. Column (6): same as column (5), but for the cube with the central sky frequency at 263.2 GHz. Column (7): number of continuum sources identified.

based approach, we perform source extraction in Spitzer/IRAC channel 2 maps and retrieve real sources from the source candidates with S/N<sub>nat</sub> down to 4.0 by crossmatching them with the IRAC sources. This is because, both in previous ALMA studies and in the ALCS, it is found that the vast majority (≈93%) of real ALMA sources (e.g., from the blind sample with *p* = 1.0) have IRAC channel (ch) 2 counterparts. We use secure IRAC sources with S/N<sub>IRAC</sub> ≥5, which corresponds to a typical detection limit of ~23 mag among the IRAC maps taken in the ALCS clusters. The details of the IRAC data reduction and related source extraction are described in Kokorev et al. (2022) and Sun et al. (2022),

respectively. We associate 40 ALMA sources with S/N<sub>nat</sub> = 4–5 with IRAC counterparts. Among these, the probability of the chance projection is estimated to be 4.6% (0.05%) at the offset of 1′′0 (0′′1) (Downes et al. 1986) based on the typical number density in the IRAC ch2 map of 0.015 arcsec⁻² (Sun et al. 2021). More precisely, we calculate the probability for each 40 ALMA source according to its spatial offset from the associated IRAC source and obtain the expected number of chance projections to be ~1 among the 40 ALMA sources. In fact, from visual inspection, we find that one source is likely caused by the chance projection, and we do not include the source in the following analysis. We refer to the remaining 39
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Differential number of positive (red) and negative (blue) sources as a function of peak S/N, summed over the 33 ALCS fields. Based on the excess of the positive to the negative sources, a total of $177_{-14}^{+12}$ sources are expected to be real down to S/N = 4.0 in the natural map (black solid line). The purity *p* exceeds 0.99 at S/N = 5.0 and 4.5 in the natural and tapered maps, respectively (black dashed line). We identify 141 sources at *p* > 0.99, referred to as the blind catalog.</figcaption></figure>

sources as the prior sample. From both approaches and the visual inspection, a total of 180 (=141 + 39) ALMA sources are identified as reliable sources, which is in excellent agreement with the expected number of the real sources of $177_{-14}^{+12}$ estimated from the excess of the positive to negative sources down to our S/N limit (Section 3.1). We list these 180 ALMA sources in Table 2, referred to as ALCS continuum source catalog.

### 3.3. Flux Measurement

We estimate the flux density for our 180 ALCS sources in several different ways. The peak count in the natural maps

($S_{\text{nat}}$) is obtained from the source extraction procedure (Section 3.1). When the source is unresolved, the peak count (Jy beam$^{-1}$) should equal the total flux density. However, as ALCS sources may be spatially resolved, we further measure the flux density in the following three methods. In the first method, we measure the peak count in the tapered map ($S_{\text{tap}}$). Here, we regard a peak above $3\sigma$ level in the uv-tapered map identified within a radius of $1''0$ from the peak in the natural map as the corresponding peak count of the source in the tapered map. If we do not identify any peaks above $3\sigma$ level within the $1''0$ radius, we use the count in the tapered map at the peak pixel position in the natural map. In the second method, we measure an enclosed flux density with an aperture radius of $2''0$ with the natural map ($S_{\text{ap}}$). We analyze the enclosed flux density growth curve as a function of the aperture radius for our ALCS sources to confirm that the $2''0$ radius sufficiently exceeds the peak of the growth curve for the majority of our ALCS sources. In the third method, we measure a spatially integrated flux density with the 2D elliptical Gaussian fitting with the CASA task IMFIT ($S_{\text{imf}}$). In the IMFIT fitting, we use the peak count, position, and beam size in the natural map for the initial values, but all parameters are ultimately fitted as free parameters. We limit the fitting area to $5'' \times 5''$. Table 2 summarizes all flux density measurements described above for our ALCS sources.

In Figure 3, we compare the flux density measurements. To understand secure trends, we only show the ALCS sources with $S/N_{\text{nat}} \ge 5$. We find that the peak count measurement in the natural map is generally lower than the $2''0$ aperture measurement. This trend implies that most of the ALCS sources are spatially resolved. Given that the typical beam size of $\sim 1''0$ in our natural maps is sufficiently large for the emission from high-redshift galaxies (FWHM $\ge 0''2-0''3$; see, e.g., Ikarashi et al. 2015; Simpson et al. 2015; Hodge et al. 2016; González-López et al. 2017a; Fujimoto et al. 2017, 2018; Gómez-Guijarro et al. 2022), this trend may at least be partially related to elongation from gravitational lensing. We also find no significant difference among the latter three measurements. To obtain secure results, we adopt $S_{\text{imf}}$ for the source with $S/N_{\text{nat}} \ge 10$, but with no fitting errors in IMFIT (e.g., fitting not converged, significant residual). We use $S_{\text{ap}}$ for the remaining sources, while we adopt $S_{\text{tap}}$ or $S_{\text{nat}}$ in several cases if the aperture encloses the emission from nearby sources or $S_{\text{ap}}$ shows negative due to the increased noise fluctuation in the aperture measurement. We choose $S_{\text{tap}}$ or $S_{\text{nat}}$, which provides a higher S/N than the other.

### 3.4. Serendipitous Line Detection

The 15 GHz frequency coverage of the ALCS data sets sometimes detects emission lines from the sources serendipitously, which could affect the continuum flux measurements. To investigate this possibility, we extract the 15 GHz spectra with $1''0$ diameter aperture for our ALCS sources in the two ALCS cubes with velocity channel widths of 30 and 60 km s$^{-1}$, respectively. We systematically calculate integrated fluxes with integration velocity ranges of 120–1200 km s$^{-1}$ for all channels, where the contribution of the continuum is subtracted by determining it with the median of the spectrum. We identify 13 line emitters whose integrated fluxes exceed the $5\sigma$ level in both velocity-channel spectra. Among these is the strong [C II] line emitter at $z = 6.0719$ reported in recent ALCS studies (Fujimoto et al. 2021; Laporte et al. 2021). We evaluate the equivalent width (EW) of these emission lines and find they contribute to
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Comparison between our various photometric methods described in Section 3.3 for S/N ≥ 5 sources. The red dots are our measurements, and the dashed line is the one-to-one relation. From left to right, we show the comparison between the 2''/0 radius aperture flux ($S_{ap}$) and the peak flux in the natural-weighted map ($S_{nat}$), the peak flux in the uv-tapered map ($S_{tap}$), and the 2D fit flux with IMFIT ($S_{imf}$).</figcaption></figure>

<table><caption>**Table 2**
Example of ALCS Continuum Sources</caption><thead><tr><th>ID<br>(1)</th><th>R.A.<br>(deg)<br>(2)</th><th>Decl.<br>(deg)<br>(3)</th><th>S/N<sub>nat</sub> (S/N<sub>tap</sub>)<br>(4)</th><th>PB<br>S<sub>nat</sub><br>(mJy beam<sup>-1</sup>)<br>(5)</th><th>S<sub>tap</sub><br>(mJy beam<sup>-1</sup>)<br>(6)</th><th>S<sub>ap</sub><br>(mJy)<br>(7)</th><th>S<sub>imf</sub><br>(mJy)<br>(8)</th><th>Flag<br>(9)</th></tr><tr><td colspan="9">Blind Catalog (N = 141)</td></tr></thead><tbody><tr><td>AC0102-C11</td><td>15.7477863</td><td>-49.2904001</td><td>5.4 (3.1)</td><td>0.96</td><td>0.39 ± 0.07</td><td>0.34 ± 0.11</td><td>0.28 ± 0.12</td><td>0.35 ± 0.1</td><td>0</td></tr><tr><td>AC0102-C22</td><td>15.7626821</td><td>-49.2864566</td><td>7.9 (7.4)</td><td>0.79</td><td>0.71 ± 0.09</td><td>0.97 ± 0.13</td><td>1.23 ± 0.14</td><td>0.96 ± 0.19</td><td>2</td></tr><tr><td>AC0102-C50</td><td>15.7475095</td><td>-49.2821233</td><td>5.8 (5.6)</td><td>0.99</td><td>0.41 ± 0.07</td><td>0.58 ± 0.1</td><td>0.66 ± 0.11</td><td>0.67 ± 0.16</td><td>2</td></tr><tr><td>AC0102-C52</td><td>15.7561196</td><td>-49.2823856</td><td>5.8 (3.2)</td><td>0.99</td><td>0.41 ± 0.07</td><td>0.33 ± 0.11</td><td>0.04 ± 0.11</td><td>0.4 ± 0.12</td><td>0</td></tr><tr><td>AC0102-C118</td><td>15.7128798</td><td>-49.2607919</td><td>34.3 (26.5)</td><td>0.92</td><td>2.63 ± 0.08</td><td>2.99 ± 0.11</td><td>3.29 ± 0.12</td><td>2.93 ± 0.14</td><td>2</td></tr><tr><td>AC0102-C160</td><td>15.7508527</td><td>-49.2677044</td><td>4.7 (4.9)</td><td>0.98</td><td>0.34 ± 0.07</td><td>0.52 ± 0.11</td><td>0.72 ± 0.12</td><td>0.87 ± 0.23</td><td>2</td></tr><tr><td>AC0102-C215</td><td>15.7288441</td><td>-49.2540865</td><td>42.4 (32.7)</td><td>0.99</td><td>3.0 ± 0.07</td><td>3.4 ± 0.1</td><td>3.83 ± 0.11</td><td>3.6 ± 0.18</td><td>2</td></tr><tr><td>AC0102-C223</td><td>15.7051647</td><td>-49.2524655</td><td>7.9 (6.7)</td><td>0.96</td><td>0.58 ± 0.07</td><td>0.72 ± 0.11</td><td>1.01 ± 0.12</td><td>0.78 ± 0.26</td><td>2</td></tr><tr><td>AC0102-C224</td><td>15.7320137</td><td>-49.2525191</td><td>91.8 (78.4)</td><td>0.98</td><td>6.57 ± 0.07</td><td>8.24 ± 0.11</td><td>9.5 ± 0.12</td><td>8.99 ± 0.26</td><td>2</td></tr><tr><td>AC0102-C241</td><td>15.742206</td><td>-49.2489458</td><td>5.4 (4.5)</td><td>0.45</td><td>0.85 ± 0.16</td><td>1.05 ± 0.23</td><td>1.08 ± 0.25</td><td>0.53 ± 0.16</td><td>1</td></tr><tr><td>AC0102-C251</td><td>15.7294397</td><td>-49.2386615</td><td>5.1 (4.4)</td><td>0.63</td><td>0.56 ± 0.11</td><td>0.72 ± 0.16</td><td>0.97 ± 0.18</td><td>0.75 ± 0.24</td><td>2</td></tr><tr><td>AC0102-C276</td><td>15.7053793</td><td>-49.2439317</td><td>9.9 (8.6)</td><td>0.81</td><td>0.86 ± 0.09</td><td>1.09 ± 0.13</td><td>1.41 ± 0.14</td><td>1.16 ± 0.18</td><td>2</td></tr><tr><td>AC0102-C294</td><td>15.7056994</td><td>-49.2514717</td><td>14.6 (19.1)</td><td>0.97</td><td>1.06 ± 0.07</td><td>2.03 ± 0.11</td><td>3.1 ± 0.12</td><td>3.16 ± 0.27</td><td>2</td></tr><tr><td>M0417-C46</td><td>64.3888354</td><td>-11.9175841</td><td>48.0 (36.7)</td><td>0.97</td><td>4.04 ± 0.08</td><td>4.26 ± 0.12</td><td>4.29 ± 0.12</td><td>4.3 ± 0.13</td><td>1</td></tr><tr><td>M0417-C49</td><td>64.4172807</td><td>-11.9168995</td><td>36.1 (28.6)</td><td>0.71</td><td>4.18 ± 0.12</td><td>4.57 ± 0.16</td><td>4.61 ± 0.16</td><td>3.32 ± 0.17</td><td>2</td></tr><tr><td>M0417-C58</td><td>64.3985559</td><td>-11.9147976</td><td>30.4 (24.3)</td><td>0.95</td><td>2.63 ± 0.09</td><td>2.9 ± 0.12</td><td>3.18 ± 0.12</td><td>2.89 ± 0.14</td><td>2</td></tr><tr><td>M0417-C121</td><td>64.4042424</td><td>-11.9055111</td><td>20.7 (16.3)</td><td>0.89</td><td>1.91 ± 0.09</td><td>2.07 ± 0.13</td><td>2.14 ± 0.13</td><td>1.93 ± 0.14</td><td>2</td></tr></tbody></table>

**Note.** The full list is presented in Table 12. Column (1): source ID, including the prefix of the cluster name + “-C.” Column (2): source coordinate of the continuum peak in the natural-weighted map. Column (3): signal-to-noise ratio of the peak pixel in the natural-weighted map. Column (4): primary beam sensitivity in the mosaic map. Column (5): source flux density at 1.2 mm measured by the peak pixel count in the natural-weighted map. Column (6): source flux density at 1.2 mm measured by the peak pixel count in the uv-tapered (2'' × 2'') map. Column (7): source flux density at 1.2 mm measured by a 2''/0 radius aperture in the natural-weighted map. Column (8): source flux density at 1.2 mm measured by 2D elliptical Gaussian fitting with CASA IMFIT with the natural maps. Column (9): flag of the IMFIT output (0, point source; 1, marginally resolved, i.e., resolved in only one direction; 2, fully resolved; -1, fitting error; -2, near the edge of the map).

the 15 GHz width continuum flux measurement at the ~1%–12% level. We summarize the line properties and their contributions to the continuum flux density in Appendix A, Table 9. We correct the continuum flux measurements for the 15 ALCS sources with detected line emission by subtracting these line contributions in the following analyses.

In Figure 4, we show three of these line emitters detected in the M0553 field as an example. The left panel presents the 15 GHz spectra, where the three lines are consistently detected at ~269 GHz. The right panel shows an HST color map overlaid with the line intensity in red contours. The counterparts of these three-line emitters are triplet multiple images spectroscopically

confirmed at $z=1.14$ (Ebeling et al. 2017), such that the millimeter line corresponds to CO(5–4). The line properties and spectra of the remaining line candidates are summarized in Appendix A. In addition to lines associated with ALCS continuum source positions, we also identify dozens of “blind” line candidates using a blind line search algorithm (e.g., Fujimoto et al. 2021), which will be presented in a separate paper.

### 3.5. Source Redshift

We estimate the source redshift for our ALCS sources. We adopt five categories for our redshift estimates and describe the details of those categories in the following subsections.
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Example of the line emitters serendipitously detected among the ALCS continuum sources. We evaluate the EW and subtract the line contributions to the continuum flux measurement in our analysis (Section 3.4). Left: ALCS 15 GHz wide spectrum in ALMA Band 6. Lines are clearly detected at ~269 GHz among three-line emitters identified in the M0553 cluster. Right: HST color image for M0553-C133/190/249, using F814W (blue), F125W (green), and F160W (red) filters. Red contours show the line intensity drawn at 6σ, 8σ, and 10σ. The counterparts of the line emitters are a triplet multiple-imaged system spectroscopically confirmed at z = 1.14 (Ebeling et al. 2017), concluding that the millimeter line is CO(5-4).</figcaption></figure>

### 3.5.1. Spectroscopic Redshift

The first category (1) is spectroscopic redshift ($z_{\text{spec}}$). We crossmatch the ALCS sources with spectroscopic catalogs in the literature (Soucail et al. 1999; Sand et al. 2004; Ebeling et al. 2009; Gómez et al. 2012; Wold et al. 2012; Biviano et al. 2013; Schmidt et al. 2014; Treu et al. 2015; Wang et al. 2015; Dessauges-Zavadsky et al. 2017; Ebeling et al. 2017; Cerny et al. 2018; Caminha et al. 2019; Caputi et al. 2021; Fujimoto et al. 2021; Richard et al. 2021; Kokorev et al. 2022; Muñoz Arancibia et al. 2023), identifying 51 ALCS sources in this category. We also newly determine $z_{\text{spec}}$ for eight sources from the serendipitous line detections in 13 ALCS sources (Section 3.4). In addition, single or multiple lines have been successfully detected in line scan follow-up observations, and we further obtain $z_{\text{spec}}$ for another eight sources. For all sources with only single-line detections, we adopt the most plausible line associated with their $z_{\text{phot}}$ distribution function. In total, 67 (= 51 + 8 + 8) sources have been regarded as $z_{\text{spec}}$ and categorized as this category. In Appendix A, we individually describe the process of the redshift determination for these line-detected sources.

### 3.5.2. Photometric Redshift

Because of the difficulty of the heterogeneous availability of optical-near-IR (NIR) data among our ALCS sources (e.g., sources fall outside HST footprints), we divide the photometric redshift ($z_{\text{phot}}$) further into four categories: (2, 3, 4, and 5). As a summary, there are

(2) seven sources with robust multiple images that refine $z_{\text{photo}}$ from their spatial configurations with the lens model;

(3) 80 sources with $z_{\text{photo}}$ from Kokorev et al. (2022), constrained by the optical-NIR spectral energy distribution (SED) analysis for the sources detected in both HST and IRAC;

(4) 18 and four sources that are not detected in HST and fall outside the HST footprints, respectively, and whose $z_{\text{phot}}$ are constrained by an SED template fitting;

(5) four sources that have no HST nor IRAC, so $z_{\text{photo}} = 3.0 \pm 2.0$ is assumed.

For (2), we crossmatch the ALCS sources with catalogs of the multiple images used in constructing the lens models in the literature. We identify seven ALCS sources in this category. Here, we adopt the same redshift uncertainty of $\Delta z = 0.5$ used to constrain the lens models.

For (3), we crossmatch them with the HST+IRAC catalog of the ALCS 33 fields (Kokorev et al. 2022). We identify 80 ALCS sources crossmatched with the HST+IRAC catalog and adopt the $z_{\text{phot}}$ estimate obtained from the SED fitting using the EAZY code (Brammer et al. 2008). This catalog uses a template set consisting of 12 templates derived from the Flexible Stellar Population Synthesis library (Conroy et al. 2009; Conroy & Gunn 2010), which enables to reproduce a much larger library spanning a range of dust attenuation, ages, mass-to-light ratios, and realistic star formation histories (e.g., bursty, slowly rising, slowly falling). The fitting details are presented in Kokorev et al. (2022). Because the mosaics for all available HST filters are combined and used for the detection image in this HST +IRAC catalog, the resulting crossmatches indicate that these 80 ALCS sources have counterparts in one or more HST band(s).

For (4), of the remaining sources, 18 lack counterparts in HST/WFC3 bands to faint limits, and are considered so-called NIR-dark (or optical-infrared dark) objects (e.g., Fujimoto et al. 2016; Franco et al. 2018; Wang et al. 2019b; Casey et al. 2019; Williams et al. 2019; Yamaguchi et al. 2019; Romano et al. 2020; Fudamoto et al. 2021; Talia et al. 2021; Fujimoto et al. 2022; Gómez-Guijarro et al. 2022; Manning et al. 2022; Barrufet et al. 2023b; Xiao et al. 2023b; Giulietti et al. 2023), while the sky positions for the remaining eight fall outside the HST/WFC3 footprints. For these 26 sources, we extract IRAC photometry using a 3''0 diameter aperture forced at the ALMA continuum positions. To mitigate contamination from nearby sources, we use the residual IRAC maps, whereby sources detected in the HST maps are modeled with the IRAC point-spread function and subtracted (Kokorev et al. 2022). Given the limited optical-NIR information, we also include Herschel photometry obtained in Sun et al. (2022) and evaluate $P(z)$ for these sources by calculating $\chi^2$ with a composite SED model obtained from 707 high-redshift dusty star-forming galaxies (Dudzevičiūtė et al. 2020). Based on these procedures, the redshift estimates fall within a range of 0.45-4.25 for the remaining ALCS sources lacking HST counterparts/data. In Appendix B, we show an example of $P(z)$ estimate with the composite SED model for these unique sources, based on upper limits along with our ALMA detection.

For (5), among the 26 sources in the category of (4), four blind-sample sources remain undetected (< 2σ-3σ) in both Herschel photometry and the forced-aperture photometry in the residual IRAC maps, despite their secure ALMA detections. These sources are likely associated with either a very faint dusty system at $z \sim 1-2$ (e.g., Aravena et al. 2020), a significantly dust-attenuated system at $z \sim 3-5$ (e.g., Umehata et al. 2020; Smail et al. 2021), or very high-redshift dusty system at $z > 7$ (Fudamoto et al. 2021; Fujimoto et al. 2022). Based on these possibilities, we assume source redshifts of $z = 3.0 \pm 2.0$ for these four sources in the following analyses. We list our redshift estimate of our ALCS sources in Table 3.
<table><caption>**Table 3**
Example of ALCS Continuum Source Properties</caption><thead><tr><th>ID</th><th>z<sub>spec</sub></th><th>z<sub>phot</sub></th><th>S<sub>obs</sub><br>(mJy)</th><th>L<sub>IR</sub><br>(L<sub>☉</sub>)</th><th>μ</th><th>z type</th><th>References</th><th>Note</th></tr><tr><td>(1)</td><td>(2)</td><td>(3)</td><td>(4)</td><td>(5)</td><td>(6)</td><td>(7)</td><td>(8)</td><td></td></tr><tr><td colspan="9">Primary Catalog (N = 125<sup>a</sup>)</td></tr></thead><tbody><tr><td>AC0102-C11</td><td>...</td><td>(3.0 ± 2.0)</td><td>0.30 ± 0.10</td><td>12.00 ± 0.19</td><td>1.6<sup>+0.7</sup><sub>-0.6</sub></td><td>5</td><td>...</td><td>No counterparts in HST, IRAC, Herschel</td></tr><tr><td>AC0102-C22</td><td>...</td><td>3.45<sup>+0.17</sup><sub>-0.17</sub></td><td>1.11 ± 0.13</td><td>12.53 ± 0.19</td><td>2.1<sup>+0.8</sup><sub>-0.8</sub></td><td>4</td><td>Appendix B</td><td>No counterparts in HST</td></tr><tr><td>AC0102-C50</td><td>...</td><td>2.38<sup>+0.10</sup><sub>-0.11</sub></td><td>0.60 ± 0.10</td><td>12.85 ± 0.07</td><td>2.4<sup>+1.0</sup><sub>-1.0</sub></td><td>3</td><td>K22</td><td>...</td></tr><tr><td>AC0102-C52</td><td>...</td><td>2.02<sup>+0.35</sup><sub>-0.16</sub></td><td>0.30 ± 0.09</td><td>12.36 ± 0.11</td><td>2.3<sup>+1.0</sup><sub>-0.9</sub></td><td>3</td><td>K22</td><td>...</td></tr><tr><td>AC0102-C160</td><td>...</td><td>1.99<sup>+0.18</sup><sub>-0.14</sub></td><td>0.65 ± 0.10</td><td>12.61 ± 0.13</td><td>2.3<sup>+0.9</sup><sub>-0.9</sub></td><td>3</td><td>K22</td><td>...</td></tr><tr><td>AC0102-C224</td><td>4.320</td><td>...</td><td>3.19 ± 0.15</td><td>13.07 ± 0.14</td><td>5.1<sup>+3.6</sup><sub>-3.6</sub></td><td>1</td><td>C21</td><td>Multiple images (AC0102-C118/215/224)</td></tr><tr><td>AC0102-C241</td><td>...</td><td>2.82<sup>+0.05</sup><sub>-0.67</sub></td><td>0.97 ± 0.23</td><td>12.70 ± 0.14</td><td>2.0<sup>+0.8</sup><sub>-0.8</sub></td><td>3</td><td>K22</td><td>...</td></tr><tr><td>AC0102-C251</td><td>...</td><td>(3.0 ± 2.0)</td><td>0.90 ± 0.16</td><td>11.84 ± 0.51</td><td>1.0<sup>+0.4</sup><sub>-0.0</sub></td><td>5</td><td>...</td><td>No counterparts in HST, IRAC, Herschel</td></tr><tr><td>AC0102-C276</td><td>...</td><td>3.46<sup>+0.30</sup><sub>-0.43</sub></td><td>1.41 ± 0.14</td><td>12.95 ± 0.10</td><td>5.2<sup>+2.7</sup><sub>-2.7</sub></td><td>3</td><td>K22</td><td>...</td></tr><tr><td>AC0102-C294</td><td>...</td><td>4.0<sup>+0.5</sup><sub>-0.5</sub></td><td>2.12 ± 0.28</td><td>12.54 ± 0.13</td><td>3.7<sup>+2.1</sup><sub>-2.1</sub></td><td>2</td><td>Ce18</td><td>Multiple images (AC0102-C223/294)</td></tr><tr><td>M0417-C49</td><td>...</td><td>3.16<sup>+0.10</sup><sub>-0.10</sub></td><td>4.69 ± 0.24</td><td>13.19 ± 0.25</td><td>2.0<sup>+0.8</sup><sub>-0.8</sub></td><td>3</td><td>K22</td><td>...</td></tr><tr><td>M0417-C121</td><td>3.652</td><td>...</td><td>3.05 ± 0.15</td><td>13.09 ± 0.14</td><td>4.0<sup>+2.3</sup><sub>-2.3</sub></td><td>1</td><td>Appendix A</td><td>Multiple images (M0417-C46/58/121)</td></tr><tr><td>M0417-C204</td><td>...</td><td>1.80<sup>+0.22</sup><sub>-0.22</sub></td><td>0.94 ± 0.13</td><td>12.57 ± 0.16</td><td>2.8<sup>+1.1</sup><sub>-1.2</sub></td><td>4</td><td>Appendix B</td><td>Outside of HST/WFC3</td></tr><tr><td>M0417-C218</td><td>...</td><td>2.00<sup>+0.07</sup><sub>-0.07</sub></td><td>1.42 ± 0.18</td><td>12.76 ± 0.29</td><td>4.8<sup>+1.9</sup><sub>-1.9</sub></td><td>4</td><td>Appendix B</td><td>No counterparts in HST</td></tr></tbody></table>

**Notes.** The full list is presented in Table 13. Column (1): ALCS continuum source ID. Column (2): Spectroscopic redshift (Section 3.5). Column (3): Photometric redshift (Section 3.5). Column (4): Observed flux density, including the corrections for the primary beam, Eddington bias (Section 3.7), and the line contamination (Section 3.4), but without the lens correction. The errors in the flux measurement and corrections are propagated. For the multiple-image system, we show an observed flux density calculated by multiplying the average values of the intrinsic flux densities and magnifications. Column (5): Infrared luminosity ($L_{IR}$) estimated by our SED fitting (see Section 6.1). Column (6): Magnification factor based on our fiducial lens model, where the error is evaluated from the propagation of the redshift uncertainty and the systematic uncertainty from different lens models (Section 3.6). Column (7): Redshift category. 1, $z_{spec}$; 2, $z_{phot}$ + lens model constraint from the multiple image positions in the literature; 3, $z_{phot}$ with EAZY; 4, $z_{phot}$ based on the composite SED from 707 dusty galaxies in AS2UDS for the sources without HST counterparts/data (Appendix B); 5, $z = 3.0 \pm 2.0$ assumed for four sources with no counterparts in all bands other than ALMA. Column (8): Reference: C21 (Caputi et al. 2021), Ce18 (Cerny et al. 2018), K22 (Kokorev et al. 2022), and Appendix A.

### 3.6. Lens Model and Magnification Correction

We construct lens models for the ALCS 33 lensing clusters. We select cluster member galaxies as well as multiple images behind clusters based on the photometric redshift, colors, and morphology of galaxies in the HST images and apply them to independent algorithms including GLAFIC (Oguri 2010; Kawamata et al. 2016, 2018; Okabe et al. 2020), LENSTOOL (Jullo et al. 2007; Richard et al. 2014; Caminha et al. 2016, 2017a, 2017b, 2019; Niemiec et al. 2020), Light-Traces-Mass (LTM; Zitrin et al. 2013, 2015), and Pseudo-Isothermal Elliptical Mass Distribution plus elliptical Navarro-Frenk-White (dPIEeNFW; Zitrin et al. 2013, 2015). In Appendix C, we summarize the number of multiple images, the accuracy of the model predictions of the multiple images, and lens models available in each of the ALCS 33 lensing clusters. In this paper, we adopt the lens model of GLAFIC as a fiducial model for our analyses, although we also use other available models to evaluate uncertainties in magnification factors. The lens model of GLAFIC is constructed in the same manner as Kawamata et al. (2016), wherein interested readers can find more specific lens modeling procedures using GLAFIC.

We apply lensing corrections to all ALCS sources located behind their respective clusters, defined as sources with photometric (spectroscopic) redshift estimates exceeding the cluster redshift by more than 0.2 (0.1). Among the ALCS sample, six groups of the sources have been spectroscopically confirmed as multiple images. Based on source positions, lens model predictions, and similar SED properties, we also classify ACT0102-C223/294 and R0032-C127/131/198 as multiple images at $z = 4.0 \pm 0.5$ and $z_{spec} = 2.391$, respectively. After

removing the cluster member galaxies and correcting the counts for the multiple images, we identify 146 out of 180 ALCS sources in the blind+prior sample as the unique sources behind the clusters. We then evaluate the magnification factor based on the peak pixel position of the ALMA continuum emission and the redshift estimate. We then calculate the intrinsic flux density by dividing the observed flux density by this magnification factor. For the multiply imaged systems, we adopt the average value of the intrinsic flux densities among the multiple images. Table 3 includes these 146 ALCS sources and their intrinsic flux density after the lens correction. The median magnification is 2.70 for the unique ALCS sources behind the clusters.

To evaluate the systematic uncertainty of the lens model, we calculate $\Delta\mu/\mu$, where $\mu$ is the magnification factor from the fiducial model, and $\Delta\mu = |\mu_{other} - \mu|$ is the difference between magnification factors of the fiducial and other available lens models ($\mu_{other}$) at random positions in the cluster field. Because the accuracy of the lens model generally depends on the richness of the multiwavelength data in the cluster field for identifying the multiple images and the cluster member galaxies, we separately evaluate the $\Delta\mu/\mu$ values for HFF, CLASH, and RELICS clusters. In Appendix D, we compare the magnification factors between our fiducial and other lens models at random pixel positions of the maps. We find a trend of increasing $\Delta\mu/\mu$ toward high $\mu$ from ~20%–40% at $\mu = 3$ to ~50%–60% at $\mu = 10$, where the trends in CLASH and RELICS fields are almost comparable. Among the CLASH and RELICS clusters, the total number of multiple images is generally the same, while the fraction of spectroscopic redshifts

<footnote><sup>a</sup> This is after correcting the lensing effect for the multiply imaged sources, resulting in a slightly lower number than that in Table 2.</footnote>
among the multiple images in CLASH is higher (see Table 11). The small differences between the CLASH and RELICS trends are consistent with the fact that the magnification uncertainty is more driven by the total number of multiple images rather than the spectroscopic redshift fraction (Johnson & Sharon 2016). We also confirm that the increasing trend of the magnification uncertainty is generally consistent with the previous HFF lens model comparison results (Meneghetti et al. 2017), where $\Delta\mu/\mu$ is estimated to be ~10% at $\mu = 3$ and ~30% at $\mu = 10$. Based on our comparison results, we adopt 20% (40%), 50%, 60%, and 80% for the systematic uncertainty in $\mu$ among the different lens models at $\mu \le 5$, $5 < \mu \le 10$, $10 < \mu \le 30$, and $\mu > 30$ in HFF (CLASH and RELICS).

We caution that choosing the median magnification from multiple lens models dismisses highly magnified sources. This is because even slight differences in the critical curve positions among different lens models will smear out the highly magnified regions when the median is sampled. In Appendix D, we also show the survey area by producing the median magnification map in AS1063 and confirm that the effective survey area is underestimated at $\mu \gtrsim 10$ than the areas of each model individually. To retain the full range of the lensing magnification, we thus use the fiducial lens model and include the systematic uncertainties from the different models in the following analysis. We further discuss potential systematics in our results in Section 5.6.

### 3.7. Simulations for Flux Measurements and Completeness

We perform MC simulations to investigate potential systematics in our flux density measurements due to Eddington bias or other unknown factors. First, we produce a pure noise mosaic map to inject artificial sources. In ALMA data cubes, multiplying every other channel by a factor of $-1$ will make the continuum emission disappear in the collapsed map. We thus collapse the $60 \text{ km s}^{-1}$ width cube after multiplying by $(-1)^t$ for every $t$th channel and regard this as the noise map. Second, we create 1000 artificial sources with uniform distributions in the total flux density and the source size for each ALCS cluster field. For the total flux density, we adopt a range of 2.5–80 times larger than $\sigma_{\text{cont}}$, where $\sigma_{\text{cont}}$ is the continuum data depth in each field. For the source size, potential correlations between the IR-emitting region size and the IR luminosity have been reported. However, it is still debated whether the correlation is positive, negative, or null (e.g., González-López et al. 2017a; Fujimoto et al. 2017, 2018; Tadaki et al. 2018; Smail et al. 2021). Given the lack of a definitive conclusion, we fix the intrinsic source size at a typical value in the literature of $0.1''/15$ with a circularly symmetric Gaussian morphology. Third, we inject the artificial sources at random positions into the noise map and calculate the expected source distortion in the observed frame based on our fiducial lens model (Section 3.6), assuming a source redshift of $z = 2$. Fourth, we perform mock observations for the artificial sources in the observed frame with CASA SIMOBSERVE. To obtain the same beam profile according to each ALCS cluster, we set the same configuration and sky coordinate as the observations for each ALCS cluster in SIMOBSERVE and produce the uv-visibility of the artificial source. We then create CLEANed maps of the artificial sources with the same CLEAN parameters as our ALCS maps and inject them at random positions in the noise map. We run these procedures for both the natural and tapered maps.

<figure><img src="image_9.png" /><figcaption>**Figure 5.** MC simulation results for input and output flux densities as a function of output $S/N_{\text{nat}}$ to correct the effects of the flux boosting and Eddington bias. We compare the four types of flux measurements ($S_{\text{nat}}$, $S_{\text{tap}}$, $S_{\text{ap}}$, and $S_{\text{imf}}$) as well as the photometry finally used in our analysis ($S_{\text{use}}$) (Section 3.3). We find that $S_{\text{use}}$ deviates from the input value by ~10% below $S/N_{\text{nat}} = 8$, where we apply the correction of 10% for our flux density measurement.</figcaption></figure>

In Figure 5, we show the comparison between output and input flux densities as a function of $S/N_{\text{nat}}$ from the MC simulations. For the output, we show the four different photometries of $S_{\text{nat}}$, $S_{\text{tap}}$, $S_{\text{ap}}$, and $S_{\text{imf}}$, in addition to the final photometry employed in those combinations, $S_{\text{use}}$ (Section 3.3). The error bars represent the 16–84th percentile in the MC simulations. In general, $S_{\text{nat}}$ and $S_{\text{tap}}$ underestimate the input value, while $S_{\text{ap}}$ and $S_{\text{imf}}$ are consistent with the input value within the $1\sigma$ range. We find that $S_{\text{use}}$ is in excellent agreement with the input value down to $S/N = 8$. Although $S_{\text{use}}$ still agrees with the input value within the $1\sigma$ range below $S/N = 8$, $S_{\text{use}}$ overestimates the input value by ~10% probably due to the Eddington bias. We thus apply a correction to our flux density measurement for the sources with $S/N < 8$ by 10%.

We also evaluate the completeness of our source detection, which is affected by the flux density and the spatial size of the sources in the observed frame. The background color in Figure 6 shows the completeness estimated with the output from the above MC simulations. We regard the sources recovered in the MC simulations if the artificial sources show a positive peak count within $1.0''$ (≈beam size) from the injected position with $S/N_{\text{nat}} \ge 5.0$ or $S/N_{\text{tap}} \ge 4.5$ for the blind sample and $S/N_{\text{nat}} \ge 4.0$ for the blind+prior sample. For comparison, Figure 6 also shows our ALCS sources in the blind sample. For the sources whose structure is well resolved with $S/N_{\text{nat}} \ge 10$ (red filled circles), the IMFIT results are used for the source size estimate, but for the other sources (white circles), we use the expected source sizes via the lensing distortion according to the source position, our fiducial lens model, the source redshift, and the intrinsic source size assumption of $0.1''/15$. We find that the majority of our ALCS sources fall in the parameter space with the completeness of >50%, but that a few sources show very low completeness of <10%. Given the potential
<figure><img src="image_10.png" /><figcaption>**Figure 6.** MC simulation results for completeness. The red circles and arrows present the spatial size in half width at half maximum (HWHM) and upper limit evaluated from CASA IMFIT, respectively, for the ALCS sources whose structure is well resolved with $S/N_{\text{nat}} \ge 10$. The white circles show the remaining ALCS sources whose spatial sizes are estimated from the lensing distortion based on the source redshifts and our fiducial lens models by assuming the intrinsic effective radius of $0''_{15}$. As an example, the background color scale denotes the completeness of our sample selection for the blind sample ($S/N_{\text{nat}} \ge 5.0$ or $S/N_{\text{tap}} \ge 4.5$; Section 3.1) based on the ALCS map taken in M2129 whose data depth and synthesized beam size are the typical value among our 33 ALCS lensing clusters (Table 1).</figcaption></figure>

uncertainty in the completeness estimate, and in order not to make our results significantly affected by the few sources with very low completeness, we perform the following analyses by using 10% as the bottom value for the completeness. Note that we confirm the predicted sizes in the above method generally being consistent with the observed sizes by IMFIT within ~40% by comparing these estimates among the spatially resolved high-S/N sources shown in red-filled circles.

### 3.8. Survey Area

We estimate the survey areas of each ALCS field by counting the areas with the primary beam (PB) sensitivity in the mosaic maps of ≥30%, where we perform the source extraction (Section 3.1). Figure 7 presents the total survey areas of the 33 ALCS lensing cluster fields before and after the lens correction based on our fiducial lens models for $z=2$ and $z=6$. Since the sensitivity and the magnification are not spatially uniform across the ALMA maps, the survey area varies according to the intrinsic flux density. Although the magnification at a given sky position should change as a function of redshift, the total survey areas show a negligible change between redshifts $z=2$ and $z=6$. In Appendix D, we also evaluate the systematic uncertainty of the survey area from the choice of the lens model, which is also confirmed to be negligible. With our fiducial model, we create the magnification maps for each cluster from $z=0$ to $z=6$, with a step of 0.5, and calculated the effective survey area at each redshift. In the following analysis, we use the closest redshift calculation result in each

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Survey areas of ALCS ($S/N_{\text{nat}} \ge 4.0$, PB ≥30%) and other large ALMA surveys in the literature (González-López et al. 2020; Casey et al. 2021; Gómez-Guijarro et al. 2022; Muñoz Arancibia et al. 2023). The light red line denotes the survey area before the lensing correction. The red dashed and solid lines represent the effect survey area after the lensing correction, assuming the redshift at $z=2.0$ and $z=6.0$, respectively, with our fiducial lens models. ALCS explores the unique parameter space toward faint and wide regimes, compared to previous ALMA surveys, demonstrating the power of the lensing support. In addition to the little difference in the results between $z=2.0$ and $z=6.0$, we confirm that the difference in the survey areas with different lens models is negligible in Appendix D. We scale the flux densities using the methods described in Section 5.2.</figcaption></figure>

ALCS source, while we adopt the $z=0$ result for the sources in the foreground of the clusters.

For comparison, we also present survey areas of other ALMA blind surveys in previous studies (González-López et al. 2020; Casey et al. 2021; Gómez-Guijarro et al. 2022; Muñoz Arancibia et al. 2023). Critically, we find that the ALCS survey area after the lens correction (red line) explores the widest and deepest parameter spaces among the ALMA blind surveys so far performed (color lines). For example, the ALCS survey area decreases at $z=6$ down to $\sim 20 \text{ arcmin}^2$ and $\sim 4 \text{ arcmin}^2$ at detection limits of 0.1 and 0.04 mJy, which are still larger than those in ASPECS (gray line) at the same detection limits by factors of $\sim 4.8$ and 2.0, respectively. This demonstrates the power of gravitational lensing, which allows us to identify unique objects that have been missed in previous ALMA surveys (e.g., Caputi et al. 2021; Fujimoto et al. 2021; Laporte et al. 2021).

## 4. Properties of ALCS Sources

Based on the analyses presented in Section 3, we characterize our ALCS sources by comparing them with other recent ALMA sources identified in similarly large surveys (total observation times $\ge 40$ hr). Figure 8 displays the intrinsic flux density ($S_{\text{int}}$) and redshift distributions of our ALCS sources (red circles). To obtain unbiased results, we only show the 146 ALCS sources after removing the cluster member galaxies and correcting for duplicated counts among the
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Intrinsic ALMA flux density ($S_{\text{int}}$) and redshift distributions of the ALCS sources after the lens correction. The red-filled circles indicate the 146 ALCS sources after removing the cluster member galaxies and adopting the average intrinsic flux density for multiply imaged systems. The spectroscopically confirmed sources are marked with red open circles. The color crosses denote ALMA sources identified in other recent ALMA surveys—707 sources from Band 7 (870 µm) snapshot follow-up observations for SCUBA2 sources in UDS (AS2UDS; Dudzevičiūtė et al. 2020), 34 sources from Band 6 deep blind observations in HUDF (ASPECS; González-López et al. 2020), 88 sources from Band 6 wide blind observations in GOODS-S (A2GS; Gómez-Guijarro et al. 2022), and 12 sources from Band 4 (2 mm) wide blind observations in COSMOS (MORA; Casey et al. 2021). For the AS2UDS and MORA sources, we scale their flux density by assuming a modified blackbody based on their best-fit $T_d$ and $\beta_d$ values presented in the literature. The top and right panels present the probability distribution for the redshift and intrinsic flux density, where the peaks are normalized. The solid lines are the median value for each sample. The dashed and dotted lines show the median values for the spectroscopic and photometric samples in the ALCS sources, respectively, which show good agreement with each other and indicate that their potential difference is negligible owing to the lensing support.</figcaption></figure>

multiply imaged systems (Section 3.6). We also present the ALMA sources identified in other recent ALMA surveys (Dudzevičiūtė et al. 2020; González-López et al. 2020; Casey et al. 2021; Gómez-Guijarro et al. 2022). The normalized probability distributions of $z$ and $S_{\text{int}}$ for each sample are drawn in the top and right panels, with the median value indicated by the solid line.

Compared to the other ALMA samples highlighted here, we find that the ALCS sources are the most widely distributed in both parameter spaces, spanning $S_{\text{int}} \simeq 0.007-3.8$ mJy and $z \simeq 0-6$. The median $S_{\text{int}}$ and $z$ values are estimated to be $S_{\text{int}} = 0.35$ mJy and $z = 2.03$, respectively. These median $S_{\text{int}}$ and $z$ values of our ALCS sources fall between those of the ASPECS and A2GS samples, but our ALCS sources explore fainter $S_{\text{int}}$ and higher $z$ parameter spaces than those two samples. The ALCS survey provides the largest sample number ($N=146$) among the aforementioned blind surveys ($N=12-88$). Our results demonstrate the capability of the lensing boost to increase sensitivity even at high redshifts, as well as the efficiency of the wide-lensing survey scheme (Section 1). The median values of $S_{\text{int}}$ and $z$ are $S_{\text{int}} = 0.32$ and $0.37$ mJy and $z = 2.05$ and $1.96$, for the photometric and spectroscopic ALCS samples, respectively. The small differences between these two samples suggest a general agreement

of their physical properties, strengthening our interpretations for the photometric sample in the statistical sense. Overall, our ALCS sources represent an optimal statistical sample to study the faint millimeter population in wide flux and redshift ranges.

Comparing the median values among the different samples, we find a general positive correlation between the source flux and redshift. A similar positive correlation is reported in the bright SMG population with the submillimeter flux density (e.g., Stach et al. 2018; Simpson et al. 2020), indicating that this correlation also exists among fainter millimeter sources than SMGs down to $S_{1.2\text{mm}} \sim 0.1-0.01$ mJy. These correlations are likely in line with the increasing luminosity evolution of the IR LFs toward high redshifts reported in previous studies (e.g., Gruppioni et al. 2013; Koprowski et al. 2017). We further discuss the redshift evolution of IR LFs in Section 6.1.

# 5. Number Counts and CIB

## 5.1. Number Counts at 1.2 mm

We derive number counts at 1.2 mm based on the most common method: directly counting the sources, correcting their purity and completeness, and obtaining the number density per survey area (e.g., Hatsukade et al. 2013; Ono et al. 2014; Carniani et al. 2015; Aravena et al. 2016; Fujimoto et al. 2016;
Dunlop et al. 2017; Umehata et al. 2017; Franco et al. 2018; Hatsukade et al. 2018; Muñoz Arancibia et al. 2018; González-López et al. 2020; Gómez-Guijarro et al. 2022). Note that the lensing effect is corrected in the measurements of the flux density and the survey area. A total of 104 and 146 sources are accounted for in our number counts' analysis based on the blind and blind+prior samples, after removing the sources that are regarded as the cluster member galaxies and correcting the counts from the multiple images (Section 3.6).

A contribution to the number counts from an identified source, $\xi$, is given by

$$ \xi(S) = \frac{p(S)}{C(S, D)A_{\text{eff}}(S)}, \qquad (2) $$

where S, C, D, A<sub>eff</sub>, and p(S) are the intrinsic flux density, completeness, intrinsic source size, survey area, and purity (Equation (1)), respectively. Then, a sum of the contributions for each flux bin is computed by

$$ n(S) = \frac{\Sigma\xi(S)}{\Delta S}, \qquad (3) $$

where $\Delta S$ is the bin width with the unit of mJy to normalize the difference due to the bin size. To evaluate the uncertainty of n (S), we include Poisson statistical errors of the source counts per bin and the uncertainty of the intrinsic flux density estimate. For the Poisson error, we use the values presented in Gehrels (1986) that are applicable to a small number of statistics. For the intrinsic flux density uncertainty, we take the following contributions into account: the random noise, the error of the absolute accuracy of the ALMA Band 6 flux calibration, and errors of the lens correction that are contributed by uncertainties from the lens model and the source redshift. Given the negligible impact of the Poisson statistical errors on the purity estimate at our high purity cut (Section 3.1), the purity uncertainty is not included in our procedure.

We perform MC simulations to include all the uncertainties described above and derive realistic number counts. We make a mock catalog of the ALCS sources whose flux densities follow Gaussian probability distributions. The standard deviations of the Gaussians are given by the combination of the random noise and the uncertainties of the absolute flux accuracy and lens correction. We adopt the measurement uncertainty for the random noise and 10% for the absolute flux accuracy.<sup>53</sup> We evaluate the $1\sigma$ uncertainty of the lens correction by propagating the magnification error from the source redshift uncertainty (Section 3.5) and the systematic lens model uncertainty among different models (Section 3.6). We produce 1000 mock catalogs and derive the number counts for each catalog in the same manner. We then evaluate the average and the 16–84th percentile of the number counts per bin, where the uncertainty of the intrinsic flux density estimate is fully taken into account. Propagating the Poisson error per bin based on the median in the MC iterations, we finally obtain the differential number counts and the associated $1\sigma$ uncertainties.

In Figure 9, we show our differential number-count estimates based on the blind sample and the blind+prior sample with the red open circles and pentagons, respectively. The estimate after applying the MC simulations based on the blind sample is also

shown in the red-filled circles. We list all the estimates in Table 4. Figure 9 shows that our study successfully covers 2.5 dex in flux density and explores the 1.2 mm number counts down to ~7 µJy, owing to the wide-area mapping of ALCS toward 33 massive lensing clusters. We find that our number-count estimates based on the blind and blind+prior samples are consistent within the $1\sigma$ uncertainties in the wide flux range.

To characterize the shape of our number counts, we fit Schechter (Schechter 1976) and double power-law (DPL) functions to our differential and cumulative number counts. The Schechter function form is given by

$$ \phi(S_{\nu}) = \phi_{0} \left(\frac{S_{\nu}}{S_{0}}\right)^{\alpha} \exp\left(-\frac{S_{\nu}}{S_{0}}\right), \qquad (4) $$

where $\phi_0$, $S_0$, and $\alpha$ are the normalization, characteristic flux density, and faint-end slope power-law index, respectively, and the DPL function form is given by

$$ \phi(S_{\nu}) = \phi_{0} \left[ \left( \frac{S_{\nu}}{S_{0}} \right)^{\alpha} + \left( \frac{S_{\nu}}{S_{0}} \right)^{\beta} \right]^{-1}, \qquad (5) $$

where the definition of $\phi_0$, $S_0$, and $\alpha$ are the same as those of Equation (4), and $\beta$ is the bright-end slope. We search the best-fit model with the MCMC method using EMCEE (Foreman-Mackey et al. 2013). Given the statistically sufficient number in our ALCS sample, we only use our ALCS sample for the fitting to remove systematic uncertainties due to different measurement approaches among the studies.

We show the best-fit Schechter and DPL functions with the solid red and dashed magenta curves in Figure 9, and summarize the best-fit parameters with the $\chi^2$/degrees of freedom (dof) values in Table 5. We also perform the Schechter and DPL function fitting in the 1000 realizations of the MC simulations, which yield the red and magenta shaded regions representing the 16–84th percentile of the 1000 best-fit Schechter and DPL functions. Figure 9 and the $\chi^2$/dof values indicate that the number counts are well represented by Schechter and DPL. In the following comparison and analyses, we adopt the best-fit Schechter function for the entire 33 ALCS clusters without any weights among the clusters as our fiducial measurement. We further discuss the impact of the faint-end slope estimate with different assumptions, subsamples, and weights according to the different qualities of the lens model in Sections 5.5 and 5.6.

## 5.2. Comparison with Previous ALMA Observations

In Figure 9, we also show previous ALMA measurements of submillimeter–millimeter number counts in the literature to compare with our estimate (Hatsukade et al. 2013; Ono et al. 2014; Carniani et al. 2015; Aravena et al. 2016; Fujimoto et al. 2016; Oteo et al. 2016; Dunlop et al. 2017; Umehata et al. 2017; Franco et al. 2018; Hatsukade et al. 2018; Muñoz Arancibia et al. 2018; Zavala et al. 2018; González-López et al. 2020; Gómez-Guijarro et al. 2022). For measurements observed at different wavelengths from 1.2 mm, we scale the flux density by assuming a typical far-IR (FIR) SED shape based on a single modified blackbody with a dust temperature $T_{dust} = 35$ K, spectral index $\beta_d = 1.8$, and $z=2$. For several measurements at 870 µm in Band 7 (Karim et al. 2013; Stach et al. 2018; Béthermin et al. 2020; Simpson et al. 2020), we use the predicted 870 and 1150 µm fluxes from the MAGPHYS

<footnote n="53">53 https://almascience.eso.org/documents-and-tools/cycle6/alma-proposers-guide</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Differential number counts at 1.2 mm. The red open circles and pentagons are the number counts derived from our ALCS blind and blind+prior samples without the MC simulations, respectively, where the error bars only take Poisson uncertainties into account. The filled red circles represent the number counts corrected by the MC simulations that implement all relevant uncertainties, such as redshift, magnification, and flux density measurements, where the error bars indicate the 16–84th percentile in the 1000 MC realizations. The solid red line and the shaded region denote our fiducial Schechter function estimate and the associated $1\sigma$ error, while the dashed magenta line presents our fiducial DPL function estimate. In the fitting, we use our number-count estimates alone to avoid potential systematic uncertainty in different measurements from different surveys. For comparison, we also show previous ALMA submillimeter/millimeter number counts in the literature (Karim et al. 2013; Hatsukade et al. 2014; Ono et al. 2014; Carniani et al. 2015; Fujimoto et al. 2016; Umehata et al. 2017; Franco et al. 2018; Hatsukade et al. 2018; Muñoz Arancibia et al. 2018; Stach et al. 2018; Béthermin et al. 2020; González-López et al. 2020; Gómez-Guijarro et al. 2022; Muñoz Arancibia et al. 2023). The black dashed line is the best-fit Schechter function using the results obtained in GOODS-S (Gómez-Guijarro et al. 2022) and HUDF (ASPECS; González-López et al. 2020). The black solid line is taken from González-López et al. (2020) obtained from the *P(D)* analysis that indicates the presence of the flattening below ~0.1 mJy; the black shaded region denotes the associated $1\sigma$ error. When comparing previous measurements of the number counts at different wavelengths to our 1.2 mm band, we scale the flux densities using the methods described in Section 5.2.</figcaption></figure>

modeling of the 707–870 µm selected SMGs in the AS2UDS survey from Dudzevičiūtė et al. (2020) to derive a relation between $S_{870\mu m}/S_{1150\mu m}$ and $S_{870\mu m}$ of the form $S_{870\mu m}/S_{1150\mu m} = 2.3 - 0.04 \times S_{870\mu m}$. This flux ratio variation is primarily driven by the evolution in the median redshift of the SMG population with submillimeter flux density (e.g., Stach et al. 2018; Simpson et al. 2020).

Although some uncertainty may remain in the flux conversion to 1.15 mm, we find that our measurements are consistent with previous studies within the errors in general. At $S_{1.2\text{mm}} > 1.0$ mJy, our measurements are in excellent agreement with the previous results from the follow-up ALMA observations for bright SMGs at 870 µm (Karim et al. 2013; Stach et al. 2018; Simpson et al. 2020). For intermediate flux densities of $S_{1.2\text{mm}} \simeq 0.1–1.0$ mJy, the ALCS number counts

fall between the scatter of previous studies and are mostly consistent with the previous measurements within the $1\sigma$ uncertainties. At $S_{1.2\text{mm}} < 0.1$ mJy, our measurements remain consistent with the results obtained in the HUDF (ASPECS; González-López et al. 2020) down to $\simeq 0.04$ mJy and with Fujimoto et al. (2016) to $\simeq 0.02$ mJy, within their $1\sigma$ uncertainties.

To compare the shape of the number counts, we also fit the DPL function in the same manner as results observed in the latest ALMA blind surveys toward GOODS-S (A2GS; Gómez-Guijarro et al. 2022) and HUDF (ASPECS; e.g., González-López et al. 2020) whose area is a part of GOODS-S. For a fair comparison, we fit the results that are also obtained by the direct counts, instead of the *P(D)* analysis presented in ASPECS. The black dashed line shows the best-fit DPL
<table><caption>**Table 4**
Differential Number Counts at 1.2 mm</caption><thead><tr><th>S<br>(mJy)<br>(1)</th><th>log(dN/dS)<br>(deg<sup>-2</sup> mJy<sup>-1</sup>)<br>(2)</th><th>⟨N⟩<br>(3)</th></tr></thead><tbody><tr><td>0.007–0.012</td><td>7.76<sup>+0.24</sup><sub>-99.99</sub></td><td>0.8</td></tr><tr><td>0.012–0.020</td><td>7.51<sup>+0.15</sup><sub>-0.25</sub></td><td>2.7</td></tr><tr><td>0.020–0.029</td><td>6.89<sup>+0.25</sup><sub>-1.04</sub></td><td>1.8</td></tr><tr><td>0.029–0.043</td><td>6.29<sup>+0.29</sup><sub>-1.03</sub></td><td>1.7</td></tr><tr><td>0.043–0.063</td><td>6.00<sup>+0.19</sup><sub>-0.26</sub></td><td>2.8</td></tr><tr><td>0.063–0.093</td><td>5.47<sup>+0.24</sup><sub>-0.95</sub></td><td>1.9</td></tr><tr><td>0.093–0.136</td><td>5.16<sup>+0.20</sup><sub>-0.47</sub></td><td>3.3</td></tr><tr><td>0.14–0.20</td><td>5.06<sup>+0.13</sup><sub>-0.19</sub></td><td>9.7</td></tr><tr><td>0.20–0.27</td><td>4.76<sup>+0.15</sup><sub>-0.20</sub></td><td>11.5</td></tr><tr><td>0.27–0.36</td><td>4.54<sup>+0.13</sup><sub>-0.19</sub></td><td>15.0</td></tr><tr><td>0.36–0.47</td><td>4.27<sup>+0.12</sup><sub>-0.18</sub></td><td>14.6</td></tr><tr><td>0.47–0.63</td><td>3.99<sup>+0.11</sup><sub>-0.16</sub></td><td>13.6</td></tr><tr><td>0.63–0.84</td><td>3.73<sup>+0.09</sup><sub>-0.12</sub></td><td>12.0</td></tr><tr><td>0.84–1.12</td><td>3.38<sup>+0.11</sup><sub>-0.12</sub></td><td>7.7</td></tr><tr><td>1.12–1.50</td><td>3.03<sup>+0.16</sup><sub>-0.20</sub></td><td>4.8</td></tr><tr><td>1.50–2.00</td><td>2.79<sup>+0.13</sup><sub>-0.26</sub></td><td>3.7</td></tr><tr><td>2.00–3.16</td><td>2.16<sup>+0.17</sup><sub>-0.30</sub></td><td>2.0</td></tr><tr><td>3.16–5.00</td><td>1.36<sup>+0.30</sup><sub>-1.36</sub></td><td>0.5</td></tr></tbody></table>

**Note.** Column (1): flux range used in the differential number counts. Column (2): Differential number counts (dN/dS) in the logarithm scale. The 1σ uncertainties are estimated from the combination of the number-count Poisson statistical errors and the uncertainty of the intrinsic flux density estimate associated with the measurement, redshift, and magnification errors (see the text). Column (3): average source number in each flux bin among the 1000 MC realizations.

<table><caption>**Table 5**
Best-fit Parameters for Our 1.2 mm Number Counts</caption><thead><tr><th></th><th>α</th><th>β</th><th>log(S<sub>s</sub>)<br>(mJy)</th><th>log(φ<sub>s</sub>)<br>(deg<sup>-2</sup>)</th><th>χ<sup>2</sup>/dof</th></tr></thead><tbody><tr><td>Schechter</td><td>-2.05<sup>+0.12</sup><sub>-0.10</sub></td><td>...</td><td>0.60<sup>+0.38</sup><sub>-0.33</sub></td><td>2.85<sup>+0.50</sup><sub>-0.52</sub></td><td>4.9/14</td></tr><tr><td>DPL</td><td>2.12<sup>+0.15</sup><sub>-0.13</sub></td><td>3.81<sup>+2.60</sup><sub>-1.65</sub></td><td>2.97<sup>+0.89</sup><sub>-0.59</sub></td><td>0.44<sup>+0.48</sup><sub>-0.59</sub></td><td>4.2/13</td></tr></tbody></table>

function obtained from the previous GOODS+HUDF results. Over the entire flux range probed, the best-fit shape of our number counts is higher than that in the GOODS+HUDF. Around the characteristic luminosity of ~1 mJy, the number counts in A2GS are a factor of ≈2 lower than our estimate. Popping et al. (2020) evaluate the effects of cosmic variance on the number-count measurements for the same wavelength and area as the ASPECS 1.2 mm image (González-López et al. 2020) to have a 2σ scatter of 1.5. However, the calculation in Popping et al. (2020) uses the survey area of 4.2 arcmin², where the PB sensitivity decreases down to 0.1. When the same calculation is performed based on the survey area of 1.8 (2.4) arcmin² with a PB sensitivity of ≥0.9 (≥0.7), the 2σ scatter increases to factors of ~5 (~2). Notably, previous UVLF and active galactic nucleus (AGN) studies also suggest that the galaxies are underdense in HUDF (e.g., Cowie et al. 2002; Moretti et al. 2003; Bauer et al. 2004; Oesch et al. 2007). Thus, we attribute the differences to cosmic variance.

Over the entire flux regime, in general, the larger survey area and the 33 independent lines of sight in our ALCS survey help to mitigate the cosmic variance compared to the previous

surveys, likely making our number counts fall between the scatter among previous studies. Moreover, our conservative purity cut of 0.99, compared to typical purity cuts of ~0.3–0.5 in previous studies (e.g., Hatsukade et al. 2013; Ono et al. 2014; Carniani et al. 2015; Fujimoto et al. 2016; Oteo et al. 2016; González-López et al. 2020), makes our results relatively immune to purity uncertainties. Additionally, our intrinsic source size assumption, instead of a point source, allows us to perform realistic completeness corrections. These aspects also likely contribute to converging within the scatter of the previous studies. We further discuss the impact of the different source size assumptions on the completeness correction in Section 5.5.

González-López et al. (2020) report that the faint-end slope has a flattening shape below $S_{1.2mm} \approx 0.1$ mJy from the triple power-law (TPL) function fit based on the $P(D)$ analysis, which we plot in Figure 9. This TPL from the $P(D)$ analysis predicts fewer number counts than the DPL from the direct counts' analysis by almost 2 dex, even with the same data, suggesting that it remains challenging to conclude the existence of any flattening at the faint end of the 1.2 mm number counts and that some caution should be exercised in the choice of the methodology. Given the average sensitivity of our ALCS maps (~60 µJy) and the secure detection cut of $S/N_{nat} = 5.0$, we require magnifications in the range of $\mu \sim 3–30$ to detect sources as faint as 0.01–0.1 mJy in our survey. This requirement still falls within the regime of lower systematic uncertainty in lens corrections, as reported in previous UVLF studies (e.g., Bouwens et al. 2017). Although the redshift uncertainty in photometric sources enhances the magnification uncertainty, we account for both uncertainties through MC simulations in our estimates.

We note that there are seven sources with secure redshift estimates ($z_{spec}$ or $z_{phot}$ + constraints from lens models due to their multiple images) below 0.06 mJy that mitigate the magnification uncertainty. In Figure 10, we summarize these seven sources labeled with their redshift and intrinsic source flux. All seven sources show distorted morphology in the HST maps whose shear orientation agrees with the predicted gravitational lensing effects, supporting the high-magnification estimates. In the prior sample, we also identify two sources (A383-C50 and A2537-C24) whose redshifts are similarly and securely estimated, falling in the same faintest regime with high-magnification estimates ($\mu > 10$) supported by strongly distorted morphology in the HST maps (see Figure 27). For example, R0600-C164 is spectroscopically confirmed at $z = 6.072$, and its local and global scale magnification factors of $163^{+27}_{-13}$ and $29^{+4}_{-7}$ are securely ensured by flux ratios and positions of its multiple images with three independent models (Fujimoto et al. 2021; Laporte et al. 2021). The intrinsic flux density is estimated to be ~0.01 mJy, such that, even if we only count R0600-C164 in the flux density bin at 0.01 mJy, the number-count estimate and its 2σ lower limit from the single-sided Poisson uncertainty Gehrels (1986) are estimated to be ~3 × 10⁷ mJy⁻¹ deg⁻² and 1.5 × 10⁶ mJy⁻¹ deg⁻², respectively; this source by itself rules out the predicted flattening shape at $S_{1.2mm} \sim 0.01$ mJy by > 0.5 dex. Naturally, when considering all seven robustly highly magnified sources summarized in Figure 10, the gap from flattening shape becomes even more significant.

Overall, even accounting for possible lensing magnification uncertainties, our results disfavor the scenario whereby the
<figure><img src="image_14.png" /><figcaption>**Figure 10.** HST/F160W 30" × 30" cutout images for seven ALCS sources falling in the intrinsically faintest millimeter regime of 0.01–0.06 mJy in the blind sample (S/N<sub>nat</sub> > 5.0 or S/N<sub>tot</sub> > 4.5), whose redshifts are spectroscopically confirmed or securely constrained by the lens models from its multiple image positions and their photometric redshifts. These secure redshift estimates mitigate the magnification uncertainty and do not support any strong flattening at ~0.1 mJy (see Section 5.2). The green contour indicates the ALMA Band 6 continuum at 2σ, 3σ, 4σ, and 5σ. The white curves represent the contours for the lensing magnification of µ = 200.</figcaption></figure>

1.2 mm number counts start flattening at $S_{1.2mm} \sim 0.1$ mJy. In González-López et al. (2020), the *P(D)* analysis is performed with the dirty image. Given the sparse density of dusty galaxies and the superb resolution of ALMA, the pixel counts below the direct detection limits might be more affected by the side lobes remained in the dirty map and the noise fluctuation rather than the weak signals from faint sources, which could be one of the causes in the different faint-end slope derived from the direct counts and the *P(D)* analysis. The point source assumption for injected sources in the *P(D)* analysis may also be another possible reason, making the completeness overestimated and thus favoring a shallow faint-end slope.

### 5.3. Comparison with Simulations

In Figure 11, we also display theoretical predictions from cosmological simulations that do not include the tweak in the initial mass function (IMF). We show semianalytical simulations from SHARK (Lagos et al. 2018, 2020) and SIDES (Béthermin et al. 2017; Bethermin et al. 2022), and a hydrodynamical simulation from EAGLE (Schaye et al. 2015). We also show a data-driven model from UniverseMachine (Behroozi et al. 2019), using empirical scaling relations to connect the star formation rate (SFR) and stellar mass of galaxies to their dust continuum emission (Popping et al. 2020). We note that SHARK and EAGLE are galaxy formation physics simulations. Hence, the IR number counts and LFs are predictions of these models. In contrast, SIDES and UniverseMachine, being empirical models, are tuned to fit some of the observations in the literature. Cowley et al. (2019) reported that EAGLE simulations underestimate the abundance of bright dusty galaxies, which we confirm at $S_{1.2mm} \gtrsim 0.2$ mJy. In Appendix E, we describe how we calculate the 1.2 mm flux density from the outputs of the EAGLE simulation.

<figure><img src="image_15.png" /><figcaption>**Figure 11.** Comparison of differential number counts at 1.2 mm with the simulation predictions. The red-filled circles, red line, red shaded region, and black dashed line are assigned as in Figure 9. The other colored lines show several model predictions for number counts (Lagos et al. 2020; Popping et al. 2020; Trayford et al. 2020; Bethermin et al. 2022). These model predictions are generally consistent with the individual data points of ALCS within $1\sigma-2\sigma$.</figcaption></figure>

Apart from the underestimate of the EAGLE simulation, these simulations show overall good agreement within $1\sigma-2\sigma$ errors with the individual data points of our 1.2 mm number counts. On the other hand, below ~0.06 mJy, we find that our best-fit Schechter function (red line) exceeds all these simulations, except for EAGLE (green line). These offsets might indicate missing physical mechanisms in the current simulations to reproduce a steeper faint-end slope in the
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Integration of the 1.2 mm number counts as a function of integration limit $S_{\text{ALMA, limit}}$. The red (black) solid and dashed lines indicate the integration of our fiducial Schechter function estimate with the ALCS (GOODS+HUDF) number counts down to the nominal detection limit of $S_{\text{limit}} = 7~\mu\text{Jy}$ (40 $\mu$Jy) and its extrapolation beyond the detection limit, respectively. The magenta (black) horizontal dotted line denotes the CIB estimate in Odegard et al. (2019) by interpolating COBE/FIRAS measurements at 217 and 353 GHz according to the central frequency of the ALCS (GOODS+HUDF) survey after subtracting the Galactic foreground from a linear combination of Galactic HI and $H\alpha$. The shaded region shows the possible range of the interpolation of the CIB by taking the estimates from the lower and upper limits of the observed frequency range in the ALMA surveys.</figcaption></figure>

1.2 mm number counts, such as a different dependence between the dust mass and metallicity at the low-mass regime. Another possibility is that the slope at the faintest end does not have the flattening at ~0.1 mJy (Section 5.2), but might be slightly shallower, and these offsets might be dismissed, as the uncertainties in the redshift estimates and lens models are reduced in future observations. We further discuss potential additional uncertainties in the faint-end slope measurement in Sections 5.5 and 5.6.

Note that these simulations also exceed the best-fit DPL obtained with the data points in GOODS+HUDF (solid black line) down to ~0.2 mJy, supporting the argument that the GOODS+HUDF region is underdense. We also note that the flattening below ~0.1 mJy is not reproduced by any of these current simulations.

## 5.4. Resolving CIB

With our fiducial Schechter parameter set (Section 5.1), we calculate the integrated flux densities, $\int_{S_{\text{limit}}}^{\infty} \phi(S) dS$, down to the flux limit of $S_{\text{limit}}$. We adopt $S_{\text{limit}} = 7~\mu\text{Jy}$ from the faintest bin of our number-count data (Figure 9). We calculate the integrated flux density to be $20.7_{-6.5}^{+8.5}~\text{Jy deg}^{-2}$ at 1.2 mm. Figure 12 shows the integrated flux density as a function of $S_{\text{limit}}$. With the integrated value, we also evaluate the fraction contributed to the CIB measurement at 1.2 mm by interpolating the COBE/FIRAS measurements at 217 and 353 GHz where the Galactic foreground emission is subtracted with a linear combination of Galactic HI and $H\alpha$ (Odegard et al. 2019). We find that the individual sources, down to $S_{\text{limit}}$ of the ALCS sample, contribute to $76.5_{-24.0}^{+31.4}\%$ of the CIB at 1.2 mm. In

Section 5.5, we discuss the impact of the different source size assumptions on this resolved fraction of the CIB.

As discussed in González-López et al. (2020), we caution that the uncertainty of the CIB itself is large due to the model dependence of the CMB and Galactic dust emission, which are subtracted from the observations to estimate the CIB. Still, the same interpolation of the COBE/FIRAS measurements and the same $S_{\text{limit}}$ suggest that ~40% of the CIB is resolved with the best-fit function for the GOODS+HUDF data points. This relatively small resolved fraction supports the interpretation that the GOODS+HUDF regions are underdense.

The uncertainty of the CIB measurement notwithstanding, our best-fit shape of the number counts ($\alpha \sim -2.0$) indicates that integration exceeds the upper limit of the CIB from 0.2 $\mu$Jy, instead of reaching down to 0 Jy. This implies that the faint-end slope of the 1.2 mm number counts may flatten or turnover within ~2 dex of our detection limit (e.g., Fujimoto et al. 2016; González-López et al. 2020). If any, the flattening or turnover flux density ($S_{\text{turn}}$) is expected to be below ~20 $\mu$Jy (and above 0.2 $\mu$Jy), and we explore the physical origins of the potential flattening/turnover as follows. Based on the typical modified blackbody assumption at the median redshift of our ALCS sources ($z \sim 2$, Section 4), the 0.2–20 $\mu$Jy flux corresponds to a dust mass $M_{\text{dust}}$ of $\sim 10^{6.5-9.5}$, $M_{\odot}$. Assuming a typical dust-to-stellar mass ratio of 0.001 (e.g., Santini et al. 2014) and the $M_{\text{star}}$–metallicity relation (e.g., Mannucci et al. 2010; Sanders et al. 2016; Iyer et al. 2018), this $M_{\text{dust}}$ range corresponds to 12+log(O/H) ~7.5–8.5. Interestingly, the gas-to-dust mass ratio (GDR) is known to increase at low metallicities of 12+log(O/H) $\lesssim 8.0$ (e.g., Asano et al. 2013; Rémy-Ruyer et al. 2014), suggesting that the observed $S_{\text{turn}}$ of 0.2–20 $\mu$Jy may be related to the break in the GDR–metallicity relation. On the other hand, an integration with a shallower faint-end slope of $\alpha \gtrsim -1.8$ does not diverge and thus does not require flattening nor turnover. Our fiducial Schechter estimate of $\alpha = -2.05_{-0.10}^{+0.12}$ is still consistent with $\alpha = -1.8$ within ~2$\sigma$. If we adopt a different size assumption for the dust continuum in the completeness correction (Section 5.5) and add more weight to the cluster fields with higher quality lens models in the number-count calculations (Section 5.6), the $\alpha$ measurement becomes close to $\alpha \sim -1.8$. Therefore, our results remain consistent with the possibility of $\alpha \sim -1.8$ and no turnover/flattening.

## 5.5. Impact from the Dust Continuum Size Distribution

In Section 5.1, we derive the number counts based on the completeness estimate with the assumption that the intrinsic source size is constant (Section 3.7). However, recent ALMA studies suggest that there exists a negative (e.g., González-López et al. 2017a; Tadaki et al. 2020; Smail et al. 2021) or positive (Fujimoto et al. 2017, 2018) correlation between the dust continuum size and IR luminosity. Interestingly, recent deep ALMA follow-up observations for a strongly ($\mu \sim 9$) lensed galaxy at $z = 7.13$ also report hints of a very extended dust continuum structure beyond the stellar continuum observed in the rest-frame UV (Akins et al. 2022). The different assumptions of the intrinsic source size distribution affect the completeness correction, especially for the strongly lensed galaxies (e.g., Kawamata et al. 2018), which may have an impact on the shape of the 1.2 mm number counts and the contribution to the CIB. Therefore, we also derive the 1.2 mm number counts with the following two different assumptions
<figure><img src="image_17.png" /><figcaption>**Figure 13.** Impact of the completeness correction from different intrinsic size distributions. Left: Faint end of the number counts. The red open circles and solid line are the same assignments as Figure 9 based on the constant size assumption for the dust continuum emission. The blue and green circles and solid lines show the number counts and the best-fit Schechter functions obtained by rerunning the same MC simulations for the corrections of the completeness and flux measurements with the assumptions of the negative (e.g., González-López et al. 2017a; Smail et al. 2021) and positive (Fujimoto et al. 2017, 2018) correlations between $L_{IR}$ and the dust continuum size, respectively. Right: integration of the 1.2 mm number counts down to our detection limit according to the best-fit Schechter functions with the different size assumptions for the dust continuum emission, suggesting that the resolved fraction of the CIB at 1.2 mm could change by factors of ~0.6–1.1 due to the different size assumptions.</figcaption></figure>

for the intrinsic source size distribution based on the literature: (1) negative correlation of $r_{e, \text{dust}} \propto L_{\text{IR}}^{-0.53}$, and (2) positive correlation of $r_{e, \text{dust}} \propto L_{\text{IR}}^{0.28}$, where $L_{\text{IR}}$ is the IR luminosity, and $r_{e, \text{dust}}$ is the effective radius of the dust continuum emission. We carry out the same MC simulations for correcting the flux measurements, completeness (Section 3.7), and associated uncertainties in the derivation of the number counts (Section 5.1) with these different size assumptions. To convert the 1.2 mm flux density to $L_{\text{IR}}$, we assume a single modified blackbody with a dust temperature $T_{\text{dust}} = 35$ K (e.g., Coppin et al. 2008) and $\beta_{\text{d}} = 1.8$ (e.g., Planck Collaboration et al. 2011).

In Figure 13, we show the 1.2 mm number counts with the different size assumptions. The red, green, and blue circles with the error bars and curves represent the average and 16–84th percentile from the MC simulation and the best-fit Schechter functions for the results that are all obtained in the same manner as Section 5.1, but with the size assumptions of constant, negative, and positive correlations, respectively. We find that the different size assumptions can change the faint-end slope. This is because, in the positive (negative) correlation, the fainter source is more compact (extended), which requires the completeness correction less (more) than that in the constant size assumption case. We also find that the different size assumption provides negligible effects on the other parameters of the DPL function.

With the best-fit Schechter parameters, we also integrate the 1.2 mm number counts and evaluate the contributions to the CIB with the different size assumptions. We find that the resolved fraction reaches ~90% and ~50% with the negative and positive correlation cases, respectively. Compared to our fiducial estimate with the constant size assumption, we find that the resolved fraction of the CIB can be changed by factors of ~0.6–1.1 by the different size assumption. To evaluate the faint end of the 1.2 mm number counts with less than the above precision, our results suggest the importance of constraining the

source size distribution for the dusty galaxies, where little has still been known at $L_{IR} \lesssim 10^{12} L_{\odot}$ (e.g., González-López et al. 2017a; Fujimoto et al. 2017, 2018; Tadaki et al. 2020; Smail et al. 2021; Gómez-Guijarro et al. 2022).

## 5.6. Potential Caveats

Although we evaluate the magnification uncertainty separately in the HFF, CLASH, RELICS clusters (Section 3.6) and include these magnification uncertainties via the MC method in the number counts' analysis (Section 5.1), we caution that there may remain further systematic uncertainties related to the different qualities of the lens models among the clusters. In general, as more multiple images are identified in deeper data, the lens models have more granularity and small-scale structure, which in turn increases the estimated source plane area as the high magnification becomes larger, compared to a less detailed model based on shallower data and fewer multiple images of the same cluster. In fact, Jauzac et al. (2015) compare the surface areas in the source plane above a given threshold magnification in one of the HFF clusters of A2744 and find that the lens models constructed before the HFF underestimate the surface area compared to the lens models constructed with the HFF data by factors of ~1.4–1.8 at $\mu \gtrsim 3$ (see Figure 4 in Jauzac et al. 2015). This indicates that the survey areas at high magnification in the CLASH and RELICS clusters (≈shallower HST data than HFF) are likely underestimated. On the other hand, the increased area at high magnifications may lead to high-magnification estimates in some galaxies currently with moderate magnification estimates. Therefore, the number density of intrinsically faint sources might not be changed much, and the impact on the measurement of the faint-end slope is uncertain.

To examine the potential impact from the different qualities of the lens model among the clusters, Figure 14 shows the 1.2 mm number counts derived in the same manner as
<figure><img src="image_18.png" /><figcaption>**Figure 14.** Comparison of differential number counts at 1.2 mm among different subsamples of ALCS clusters. The red circles and line are the same assignments with Figure 9. The open blue, orange, and green circles and lines indicate the number counts and their best-fit Schechter functions derived only with the HFF, CLASH, and RELICS clusters, respectively. The red and gray shaded regions represent the $1\sigma$ areas evaluated with the Bootstrap tests by sampling 33 clusters from all ALCS clusters (including the duplication) and only five clusters from the CLASH and RELICS clusters, respectively. The blue circles and upper limits denote the previous ALMA Band 6 survey results in HFF (Muñoz Arancibia et al. 2023), where the data are shallower than ALCS.</figcaption></figure>

Section 5.1, but separately with the HFF, CLASH, and RELICS clusters. The dashed line indicates the best-fit Schechter function for each subsample, and its faint-end slope is summarized in Table 6. The faint-end slopes are steeper in CLASH and RELICS ($\alpha_{CLASH} = -2.24^{+0.24}_{-0.20}$, $\alpha_{RELICS} = -1.91^{+0.17}_{-0.14}$) than that in HFF ($\alpha_{HFF} = -1.71^{+0.34}_{-0.29}$), while the large uncertainty in the HFF results still makes all these results consistent within $1\sigma$ due to its small statistics from the small survey area. This indicates that the underestimate of the high-magnification area in less detailed lens models might be related to the steeper faint-end slope in CLASH and RELICS than in HFF, while the effect from the cosmic variance might still be a more dominant factor. In fact, the faintest data point in HFF deviates from the best-fit Schechter function with all 33 ALCS clusters (red line) by a factor of ~6. The sum of the survey areas of the CLASH and RELICS clusters contributes to ≈80% of the total 33 clusters (Appendix D), and the practical factor of the underestimate is ~1.1–1.5. This is much smaller than the above factor to fill the gap, suggesting that the steeper faint-end slope results in CLASH, RELICS, and the entire 33 clusters than in HFF are unlikely solely due to the underestimate of the high-magnification area.

To evaluate the cosmic variance effect in the five HFF clusters, we also perform the Bootstrap test by randomly sampling five clusters from CLASH and RELICS. In Figure 14, the gray shaded region shows the $1\sigma$ region from the 1000 realizations of the Bootstrap test. For comparison, we also show the red shaded region, corresponding to the $1\sigma$ region obtained from another Bootstrap test by sampling 33 clusters among all ALCS clusters, including the duplication. We find that the best-fit Schechter function in HFF is included in the $1\sigma$ region from the five random CLASH+RELICS clusters. We thus conclude that we cannot rule out the possibility that the

<table><caption>**Table 6**
Faint-end Slope $\alpha$ in Various Data and Assumptions</caption><thead><tr><th>Data Set<br>(1)</th><th>FIR Size<br>(2)</th><th>$\alpha$<br>(3)</th></tr></thead><tbody><tr><td>ALL ($N=33$)</td><td>Const.</td><td>$-2.05^{+0.12}_{-0.10}$</td></tr><tr><td>ALL ($N=33$, with weight)</td><td>Const.</td><td>$-1.96^{+0.14}_{-0.13}$</td></tr><tr><td>HFF only ($N=5$)</td><td>Const.</td><td>$-1.71^{+0.34}_{-0.29}$</td></tr><tr><td>CLASH only ($N=12$)</td><td>Const.</td><td>$-2.24^{+0.24}_{-0.20}$</td></tr><tr><td>RELICS only ($N=16$)</td><td>Const.</td><td>$-1.91^{+0.17}_{-0.14}$</td></tr><tr><td>Bootstrap ($N=5$)</td><td>Const.</td><td>$-2.02^{+0.16}_{-0.15}$</td></tr><tr><td>Bootstrap ($N=33$)</td><td>Const.</td><td>$-2.05^{+0.13}_{-0.13}$</td></tr><tr><td>$\mu \lesssim 10$ only</td><td>Const.</td><td>$-1.94^{+0.18}_{-0.15}$</td></tr><tr><td>$\mu \lesssim 3$ only</td><td>Const.</td><td>$-2.00^{+0.34}_{-0.23}$</td></tr><tr><td>ALL ($N=33$)</td><td>$\propto L_{IR}^{a}$</td><td>$-1.86^{+0.16}_{-0.14}$</td></tr><tr><td>ALL ($N=33$)</td><td>$\propto L_{IR}^{b}$</td><td>$-2.12^{+0.13}_{-0.11}$</td></tr></tbody></table>

**Note.** Column (1): various data sets for the $\alpha$ measurement. Column (2): FIR size assumption in the completeness correction estimate. Column (3): derived faint-end slope of $\alpha$ in the Schechter function.

<sup>a</sup> Assuming positive and negative correlations with $L_{IR}$ (see Section 5.5).

different faint-end slope results among different clusters are still dominated by the cosmic variance.

As an alternative approach to take the different qualities of the lens models among the clusters, we also add weights to the number counts according to the quality of the lens model. Based on the number of the multiple images ($N_{img}$) identified in each cluster field (Table 11), we classify the five HFF clusters and the R1347 and M1206 as a group of the good model with $N_{img} > 100$. Following the results of the relative magnification uncertainty among the lens models as a function of the total number of the multiple presented in Johnson & Sharon (2016), we add 4 times more weights on the number counts of the good model clusters than the others. This yields $-1.96^{+0.14}_{-0.13}$, which is consistent with other results within $1\sigma$.

To reduce any systematic uncertainty from lensing effects, we also estimate the faint-end slope only with the data points at $>0.03$ mJy and $>0.1$ mJy from all 33 clusters, which are equal to using only the sources with $\mu \lesssim 10$ and $\mu \lesssim 3$, given the typical $5\sigma$ detection limit of 0.3 mJy ($=5 \times 0.06$ mJy; see Table 1). We obtain the best-fit $\alpha$ of $-1.94^{+0.18}_{-0.15}$ and $-2.00^{+0.34}_{-0.23}$ for the sources only with $\mu \lesssim 10$ and $\mu \lesssim 3$, respectively. Both $\alpha$ measurements are consistent with our fiducial estimate ($-2.05^{+0.12}_{-0.10}$) down to ~0.01 mJy within $1\sigma$.

Table 6 summarizes the $\alpha$ measurement for each subsample above. We also list the $\alpha$ measurements by adopting the different FIR size assumptions with all 33 ALCS clusters (Section 5.5). We find that all these $\alpha$ measurements are consistent within the ~$1\sigma-2\sigma$ errors, regardless of the subsample and the FIR size assumptions. These various test results suggest the $\alpha$ measurement with a conservative uncertainty would be $\approx -2.0^{+0.2}_{-0.3}$.

# 6. Luminosity Function and Cosmic SFR Density

## 6.1. IR Luminosity Function at $z \simeq 1-8$

While the IR LFs have been extensively studied with Herschel (e.g., Gruppioni et al. 2013; Magnelli et al. 2013; Wang et al. 2019a), the results are inevitably affected by the large beam size, which causes the source blending and the confusion limit in the source detection. The IR LFs have also
been studied with recent ALMA observations (e.g., Koprowski et al. 2017; Hatsukade et al. 2018; Gruppioni et al. 2020; Zavala et al. 2021), while only a limited luminosity and redshift range have been constrained by the small statistics due to the small field of view (FOV) of ALMA.

Taking advantage of the large blind sample of our ALCS sources (Section 4), we also analyze the IR LFs based on their secure redshift estimates that benefit from the homogeneous HST and IRAC data sets. Because we confirm the consistency between the blind and blind+prior samples in the number counts (Section 5.1), we use the blind+prior sample to obtain the statistically reliable results in the following analyses.

To calculate $L_{IR}$ for our ALCS sources, we use a panchromatic SED fitting tool STARDUST (Kokorev et al. 2021) based on the IRAC (Kokorev et al. 2022), Herschel (Sun et al. 2022), and ALMA photometry (Section 3.3), including the upper limits ($2\sigma$). The STARDUST code utilizes Draine & Li (2007; hereafter DL07) templates for the IR-millimeter SED with the dust emission model, including the additional updates from Draine et al. (2014; see also Aniano et al. 2020). These models describe the contribution from warm dust and polycyclic aromatic hydrocarbon (PAH) features in the photodissociation regions, together with cold dust in the diffuse part of the interstellar matter. To ensure our fits are physically meaningful, we restrict our distribution of the fitting parameters as follows. The minimum radiation field intensity ($U_{min}$) acts as a proxy for $T_{dust}$ in the DL07 set. We restrict $U_{min}$ to lie in the $3.0 < U_{min} < 50.0$ range, which corresponds to a range of $20\,K \lesssim T_{dust} \lesssim 50\,K$. We additionally restrict the fraction of the total dust mass locked in PAHs ($q_{PAH}$) between 0% and 10%. We have fixed $U_{max} = 10^6$ and $\alpha = 2$, following both Magdis et al. (2012) and Kokorev et al. (2021). Among our 180 ALCS sources, 125 sources are detected at least with one band in Herschel/SPIRE or PACS (Sun et al. 2021), providing moderately secure FIR SED models based on the ALMA 1.2 mm detection and the $2\sigma$ upper limits in the nondetection bands. Although caution may arise, we also infer $L_{IR}$ for the remaining sources without any Herschel detection based on the $2\sigma$ upper limits and the secure ALMA 1.2 mm detection.

Based on the associated X-ray emission, three ALCS sources are reported to be lensed AGNs (A370-C110, M0416-C117, and M0329-C11; Uematsu et al. 2023), and we use the $L_{IR}$ estimates after subtracting the AGN component evaluated in Uematsu et al. (2023) for these three ALCS sources. For the other ALCS sources, we do not take the AGN component into account in the IR-millimeter SED fitting with STARDUST, while, after the fitting, we systematically subtract the typical fraction of 20%, which is estimated from deep stacking and dedicated SED analysis for faint millimeter sources (Dunlop et al. 2017). Our $L_{IR}$ measurements are summarized in Table 12. We confirm that these measurements show no systematic offsets from the estimates presented in Sun et al. (2022) that use another fitting code of MAGPHYS (da Cunha et al. 2015) within $1\sigma$ uncertainties. The consistency in the FIR SED fitting outputs between MAGPHYS and CIGALE (Boquien et al. 2019) has also been confirmed in our ALCS sources (Uematsu et al. 2024).

We adopt the $V_{max}$ method (Schmidt 1968) to evaluate the IR LFs. This method uses the maximum observable volume of each source. The LF gives the number of the ALMA sources in

a comoving volume per logarithm of luminosity, given by

$$ \Phi(L, z) = \frac{1}{\Delta L} \sum \frac{1}{C V_{\text{max}, i}}, \qquad (6) $$

where $V_{max,i}$ is the maximum observable volume of the $i$th source, $C$ is the completeness, and $\Delta L$ is the width of the luminosity bin in the logarithm scale. We calculate $V_{max}$ as the integration of comoving volume spherical shells, given by

$$ V_{\text{max},i} = \int_{z_{\text{min}}}^{z_{\text{max}}} \frac{\Omega(S_i(z)) dV}{4\pi dz} dz, \qquad (7) $$

where $S_i(z)$ is the flux density of $i$th source observed at $z$, $\Omega$ is the effective survey area that can detect the flux density of the $i$th source, and $z_{min}$ and $z_{max}$ are maximum and minimum redshifts of a redshift bin. To study the faint ALMA sources that are gravitationally lensed behind the massive lensing clusters ($z \lesssim 0.5$), we focus on the redshift range beyond $z = 0.6$ and derive the IR LFs with the redshift bins of $0.6 \le z < 1.0$ ($N=13$), $1.0 \le z < 2.0$ ($N=46$), $2.0 \le z < 3.0$ ($N=39$), $3.0 \le z < 4.0$ ($N=30$), $4.0 \le z < 6.0$ ($N=12$), and $6.0 \le z < 7.5$ ($N=1$). To include all relevant uncertainties, such as the source redshift and the lensing magnification, we similarly perform the MC simulations as Section 5.1. We make a mock catalog of the ALCS sources whose flux densities and redshifts are randomly redistributed based on the Gaussian probability distribution and recalculate the magnification according to the redshift. We divide each luminosity and redshift bin and repeat this procedure 1000 times. In Figure 15, we show our IRLFs results at $z = 0.6-7.5$. The red circles and error bars denote the average and 16–84 percentile from the MC realization.

For comparison, Figure 15 also presents previous IR LF measurements in the literature from Herschel (Gruppioni et al. 2013; Magnelli et al. 2013; Wang et al. 2019a) and ALMA observations (Koprowski et al. 2017; Hatsukade et al. 2018; Gruppioni et al. 2020; Fudamoto et al. 2021; Barrufet et al. 2023a). We adopt the central values of the redshift ranges in the previous measurements to assign them to our redshift bins. When there are multiple measurements in a redshift bin from one literature source, we show the previous measurement whose central value of the redshift range is the closest to our redshift bin.

To obtain statistically stable results in the highest-redshift bin of $z = 6.0-7.5$, we also calculate the volume density of several IR-detected sources serendipitously identified in blind surveys in the literature. Watson et al. (2015) report the faint dust emission from a Lyman-break galaxy (LBG) at $z = 7.13$ (Wong et al. 2022) in an ALMA blind continuum survey in A1689. Fujimoto et al. (2016) estimate the effective survey area of this ALMA survey in A1689 in the same manner as Section 3.8, in which we use the survey volume calculation. Endsley et al. (2022) identify 41 luminous LBGs at $z \simeq 6.6-6.9$ by combining the narrowband and broadband filters in $1.5 \, \text{deg}^2$ COSMOS field, where one of the LBGs turns out to be obscured AGN population and with $L_{IR}$ of $9.0 \times 10^{12} L_{\odot}$ after subtracting the AGN contribution and spectroscopically confirmed at $z = 6.853$ (Endsley et al. 2023). We compute the survey volume with the area of $1.5 \, \text{deg}^2$ and the redshift range of $z = 6.6-6.9$. Fujimoto et al. (2022) report the discovery of a red quasar at $z = 7.19$, GNz7q, enshrouded by a dusty starburst host in the systematic analysis of all publicly
<figure><img src="image_19.png" /><figcaption>**Figure 15.** ALCS IR LF measurements at $z = 0.6–7.5$. The red-filled circles represent our ALCS sources corrected from the MC simulations that implement all relevant uncertainties, such as redshift, magnification, and flux density measurements, where the error bars indicate the 16–84th percentile in the 1000 MC realizations (see Section 5.1). Previous IR LF measurements with Herschel (Gruppioni et al. 2013; Magnelli et al. 2013; Wang et al. 2019a) and ALMA(+SCUBA2) are presented in gray and colored symbols, respectively. The red line presents the best-fit DPL functions, where we include our measurements and other ALMA blind survey results (Koprowski et al. 2017; Hatsukade et al. 2018) in the fitting. The other color lines show the predictions from simulations (Lagos et al. 2020; Trayford et al. 2020; Bethermin et al. 2022; Vijayan et al. 2022) and an empirically calibrated model (Zavala et al. 2021).</figcaption></figure>
<table><caption>**Table 7**
ALCS IR LFs</caption><thead><tr><th rowspan="2">log(L<sub>IR</sub>)[L<sub>☉</sub>]</th><th colspan="6">log(Φ)[Mpc<sup>-1</sup> dex<sup>-1</sup>]</th></tr><tr><th>0.6 ≤ z < 1.0</th><th>1.0 ≤ z < 2.0</th><th>2.0 ≤ z < 3.0</th><th>3.0 ≤ z < 4.0</th><th>4.0 ≤ z < 6.0</th><th>6.0 ≤ z < 7.5</th></tr></thead><tbody><tr><td>9.9–10.4</td><td>...</td><td>-0.85<sup>+0.52</sup><sub>-0.76</sub></td><td>...</td><td>...</td><td>...</td><td>-2.49<sup>+0.52</sup><sub>-0.76</sub></td></tr><tr><td>10.4–10.9</td><td>-3.72<sup>+2.35</sup><sub>-99.9</sub></td><td>-1.93<sup>+0.3</sup><sub>-0.34</sub></td><td>...</td><td>-1.94<sup>+0.3</sup><sub>-0.34</sub></td><td>-3.24<sup>+0.52</sup><sub>-0.76</sub></td><td>...</td></tr><tr><td>10.9–11.4</td><td>-2.13<sup>+0.4</sup><sub>-1.24</sub></td><td>-2.54<sup>+0.2</sup><sub>-0.22</sub></td><td>-2.86<sup>+0.3</sup><sub>-0.34</sub></td><td>-2.96<sup>+0.25</sup><sub>-0.28</sub></td><td>...</td><td>...</td></tr><tr><td>11.4–11.9</td><td>-3.38<sup>+0.54</sup><sub>-0.25</sub></td><td>-2.95<sup>+0.13</sup><sub>-0.15</sub></td><td>-3.05<sup>+0.1</sup><sub>-0.12</sub></td><td>-3.98<sup>+0.26</sup><sub>-0.28</sub></td><td>-3.78<sup>+0.22</sup><sub>-0.25</sub></td><td>...</td></tr><tr><td>11.9–12.2</td><td>-3.59<sup>+0.57</sup><sub>-0.28</sub></td><td>-3.16<sup>+0.1</sup><sub>-0.17</sub></td><td>-3.59<sup>+0.12</sup><sub>-0.17</sub></td><td>-3.68<sup>+0.2</sup><sub>-0.22</sub></td><td>-4.3<sup>+0.3</sup><sub>-0.34</sub></td><td>...</td></tr><tr><td>12.2–12.6</td><td>-4.28<sup>+0.52</sup><sub>-99.9</sub></td><td>-3.62<sup>+0.2</sup><sub>-0.22</sub></td><td>-3.91<sup>+0.22</sup><sub>-0.25</sub></td><td>-3.93<sup>+0.25</sup><sub>-0.28</sub></td><td>-4.51<sup>+0.37</sup><sub>-0.45</sub></td><td>...</td></tr><tr><td>12.6–12.9</td><td>...</td><td>-4.59<sup>+0.52</sup><sub>-99.9</sub></td><td>-4.59<sup>+0.52</sup><sub>-99.9</sub></td><td>-4.27<sup>+0.37</sup><sub>-0.45</sub></td><td>...</td><td>...</td></tr><tr><td>12.9–13.2</td><td>...</td><td>-4.5<sup>+0.52</sup><sub>-0.76</sub></td><td>-4.52<sup>+0.52</sup><sub>-99.9</sub></td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

available HST archive (~3 deg²). After subtracting the AGN contribution, the host galaxy is estimated to have $L_{IR}$ of $1.1 \times 10^{13}$. Given that GNZ7q is identified from the unambiguous Lyman-break feature in the HST bands in the first place, we adopt the redshift range of $\Delta z = 1.0$, which is a typical range of the dropout technique to identify $z \sim 7$ galaxies, and survey area of 3 deg² to compute its survey volume. The errors of these estimates are obtained from the Poisson confidence limit of 84.13% in Gehrels (1986). Because the sources in Endsley et al. (2022) and Fujimoto et al. (2022) are identified initially in the optical-NIR surveys, we place the lower limits for these estimates. Chen et al. (2022) perform the ALMA line-scan observations for the brightest submillimeter sources at 850 µm($S_{850\mu m} > 5$ mJy) and identify no sources at $z > 6$, despite the redshift range covered in the line scan including $z = 6.0–7.2$. We thus place the 84.13% upper limit from Poisson statistic (Gehrels 1986), where $L_{IR}$ is calculated by assuming the single modified blackbody at $z = 6.5$ with $S_{850\mu m} = 5$ mJy with the typical $T_{dust}$ of 35 K (e.g., Coppin et al. 2008) and $\beta_d = 1.8$ (e.g., Planck Collaboration et al. 2011).

In Figure 15, we find that our IR LF estimates explore the faint end by ~1–2 dex in most redshift bins, compared to the previous studies. We also find that our IR LF estimates are generally consistent with the previous studies within the $1\sigma$ uncertainties.

To evaluate the shape of IR LFs, we perform the DPL fit to our IR LF estimates. For the LF form, we define the DPL by

$$ \Phi(L) = \Phi_{*} \left[ \left( \frac{L}{L_{*}} \right)^{\alpha} + \left( \frac{L}{L_{*}} \right)^{\beta} \right]^{-1}, \quad (8) $$

where $L_*$ is the characteristic IR luminosity, and $\Phi_*$ is the volume density at $L_*$. In the fitting, we make the faint- and bright-end slopes $\alpha$ and $\beta$ free parameters in the redshift bins at $z = 0.6–4.0$. On the other hand, we fix both of them at $z \ge 4$ by using the best-fit results at $z = 1.0–2.0$, where the IR LF shape is the most securely determined with the largest sample size. We include the data points in the fitting that are obtained from other ALMA blind surveys in the literature (Koprowski et al. 2017; Hatsukade et al. 2018). In the highest-redshift bin at $z = 6.0–7.5$, we also use measurements from Fudamoto et al. (2021) and our calculation from Watson et al. (2015) for the fitting, given their serendipitous discovery nature, while we do not include the measurements from Barrufet et al. (2023a) due to its targeted observation aspect. In addition, because the bright end $z = 6.0–7.5$ is also constrained both from the lower

<table><caption>**Table 8**
Best-fit DPL Parameters of ALCS IR LFs</caption><thead><tr><th>z Range<br>(1)</th><th>α<br>(2)</th><th>β<br>(3)</th><th>log(Φ<sub>*</sub>)<br>(4)</th><th>log(L<sub>*</sub>)<br>(5)</th></tr></thead><tbody><tr><td>0.6–1.0</td><td>0.94 ± 0.39</td><td>(3.72)</td><td>-3.65 ± 0.12</td><td>12.15 ± 0.06</td></tr><tr><td>1.0–2.0</td><td>0.94 ± 0.43</td><td>3.72 ± 1.10</td><td>-3.38 ± 0.62</td><td>12.21 ± 0.27</td></tr><tr><td>2.0–3.0</td><td>0.93 ± 0.26</td><td>(3.72)</td><td>-3.86 ± 0.3</td><td>12.52 ± 0.12</td></tr><tr><td>3.0–4.0</td><td>1.04 ± 0.20</td><td>(3.72)</td><td>-4.49 ± 0.36</td><td>12.73 ± 0.15</td></tr><tr><td>4.0–6.0</td><td>(0.94)</td><td>(3.72)</td><td>-4.77 ± 0.19</td><td>12.72 ± 0.09</td></tr><tr><td>6.0–7.5</td><td>(0.94)</td><td>(3.72)</td><td>-5.26 ± 0.43</td><td>12.71 ± 0.26</td></tr></tbody></table>

**Note.** Column (1): redshift bin range of ALCS IR LFs. Column (2): Faint-end slope in DPL function. It is fixed in the fitting in the redshift bins of $z = 4.0–6.0$ and $z = 6.0–7.5$. Column (3): Bright-end slope in DPL function. It is fixed in the fitting in the redshift bins, except for $z = 1.0–2.0$. Column (4): characteristic volume density in DPL function. Column (5): characteristic IR luminosity in DPL function.

and upper limits (Endsley et al. 2022; Chen et al. 2022; Fujimoto et al. 2022), we also calculate the central value from these upper and lower limits and use it for the fitting, where we regard the difference between the average and the upper/lower limits as the errors. Our IR LF measurements are summarized in Table 7.

In Figure 15, we present the best-fit DPL functions with the red solid curves that are extrapolated with the dashed curves below the faintest luminosity bin. The best-fit parameters in each redshift bin are summarized in Table 8. We find that ours and previous IR LF measurements are well represented by the DPL functions. However, we find that our best-fit values of the faint-end slope, $\alpha \simeq 0.9–1.0$, are higher than the previous reports of $\alpha \simeq 0.4–0.6$ (Gruppioni et al. 2013; Magnelli et al. 2013; Koprowski et al. 2017; Casey et al. 2018a; Zavala et al. 2021), while our results are close to the reports in the UVLF studies of $\simeq 1.5–2.0$ at similar redshifts (e.g., Bouwens et al. 2022; Finkelstein & Bagley 2022). A wide dynamic range is essential to determine the IR LF shapes securely. Compared to our measurements, the previous IR LF measurements generally cover only the part of the IR luminosity range close to $L_*$, which could be the reason for the difference. We note for the previous Herschel studies that the faint-end slope is generally fixed with a canonical value estimated at their lowest redshift bin at $z \lesssim 0.3$ (Gruppioni et al. 2013; Magnelli et al. 2013). The UV LF studies suggest that the faint-end slope becomes steeper toward high redshifts (e.g., Bouwens et al. 2022; Finkelstein & Bagley 2022). Besides, the physical process of triggering the dusty star-forming activity may be much different between the
local and high-redshift universe due to the rapid decrease of the gas fraction toward lower redshifts (e.g., Casey et al. 2014; Tacconi et al. 2020). The superb sensitivity and spatial resolution of ALMA would also contribute to overcoming the confusion limit and reduce the impact of the incompleteness of identifying the faint sources. We thus think that the difference in the faint-end slope estimate would be explained by the combination of these reasons described above.

For comparison, Figure 15 also shows predictions from recent simulations (Lagos et al. 2020; Trayford et al. 2020; Bethermin et al. 2022; Vijayan et al. 2022) and an empirically calibrated model (Zavala et al. 2021) in colored lines. Note that the turnover at faint luminosities in EAGLE (green) and FLARES (brown) simulations is caused by the incompleteness due to the mass resolution limit at ~$10^9 M_\odot$ in their calculations. We find that our measurements and the best-fit DPL functions are generally consistent with these recent simulations within the $1\sigma-2\sigma$ errors, except for the model of Zavala et al. (2021). Since the Zavala et al. (2021) model is calibrated to the faint-end 1.2 mm number counts only from the ASPECS results in HUDF (González-López et al. 2020), its deviation from ours and other measurements and simulation predictions is likely explained by the underdense of HUDF (Section 5.2). We also find that our best-fit DPL functions start to deviate from the predictions of the simulations at the faint end of $L_{IR} < 10^{11} L_\odot$ in the redshift bins of $z=0.6-4.0$. These trends are consistent with the results in the 1.2 mm number counts (Figure 11). These results might indicate that our current understanding of the dependence between dust mass and metallicity among low-mass galaxies might be different. Still, since the different size distribution may raise the same impact on the completeness correction in the faint-end regime as the 1.2 mm number counts (Section 5.5), the positive correlation between the IR luminosity and size may mitigate the deviation by a factor of ~0.6.

## 6.2. Redshift Evolution of IR LF

In Figure 16, we compare the best-fit DPL functions from $z=0.6$ to $7.5$. We find that $\Phi_*$ and $L_*$ decrease and increase toward high redshifts, respectively. These trends are thought to represent the apparent “down-sizing” of the millimeter source population (e.g., Koprowski et al. 2017), which is consistent with the evolution in median redshift with the 1 mm flux density (Section 4). In Figure 17, we also summarize the redshift evolution of $\Phi_*$ and $L_*$ to quantitatively analyze these trends. For comparison, we also show the previous measurements with Herschel and SCUBA2 in the gray symbols (Gruppioni et al. 2013; Magnelli et al. 2013; Wang et al. 2019a; Lim et al. 2020) and with ALMA+SCUBA2 in the magenta symbol (Koprowski et al. 2017). The solid magenta lines denote the best-fit empirical functions in Koprowski et al. (2017), where the dashed line presents their extrapolations. The blue lines also show the best-fit empirical function obtained in the latest UV LF study (Bouwens et al. 2022).

We confirm that our measurements are consistent with the general trends of the evolution of IR LFs characterized by positive $L_*$ and negative $\Phi_*$ toward high redshifts reported in previous studies. We find that our measurements are the most closely consistent with the measurements obtained with ALMA +SCUBA2 (Koprowski et al. 2017), while we also find that our measurements are lower (higher) than the previous studies with single-dish telescopes (Gruppioni et al. 2013;

<figure><img src="image_20.png" /><figcaption>**Figure 16.** Evolution of IR LFs based on our best-fit DPL functions presented in Figure 15. We find the overall trends characterized by the positive luminosity evolution coupled with the negative density evolution.</figcaption></figure>

Magnelli et al. 2013; Wang et al. 2019a; Lim et al. 2020). As discussed in Section 6.1, the difference is likely explained by the combination of the partial coverage of $L_{IR}$ and the coarse spatial resolutions in the previous single-dish surveys.

Our measurements suggest a flattening in the correlation between the redshift and $L_*$ at $z \sim 4$. A similar flattening has been reported in the UV LF studies at $z \sim 2.5$ (e.g., Bouwens et al. 2022), which is thought to be explained by the maximum value of UV luminosity due to the increasing importance of dust extinction in the highest stellar mass and SFR sources. These results might indicate that there is also a maximum value in the IR luminosity, especially at such high redshifts due to the relatively short cosmic timescale for the required dust mass formation and evolution. The maximum IR luminosity might also be contributed by the decreasing trend of the dust continuum size toward high redshifts (e.g., Fujimoto et al. 2017; Gómez-Guijarro et al. 2022), where the Eddington limit would regulate the high-surface density of IR luminosity.

## 6.3. Cosmic SFRD

By integrating our IR LF measurements, we calculate the volume average of the integrated IR luminosity and thus investigate the obscured side of the cosmic SFR density (SFRD) $\psi_{SFRD(\text{IR})}$. We integrate the best-fit DPL functions down to our detection limit of $L_{IR} = 10^{10} L_\odot$. For conversion from $L_{IR}$ to SFR, we use $\text{SFR} = KL_{IR}$, where $K = 4.5 \times 10^{-44} M_\odot \text{yr}^{-1} \text{erg}^{-1} \text{s}$, which is the same calibration as Madau & Dickinson (2014), but we multiply a factor of 0.63 to apply the correction from the Salpeter IMF to the Chabrier IMF.

The left panel of Figure 18 shows the redshift evolution of $\rho_{SFR(\text{IR})}$ from our IR LF measurements. For comparison, we also present two possible scenarios of “Dust Rich” and “Dust Poor” in the universe presented in Casey et al. (2018b) and the model prediction (Lagos et al. 2020). We also show previous measurements presented in Zavala et al. (2021) that perform a backward evolution model to constrain the evolution of the IR LFs to reproduce the observed number counts at 1.2 and 3 mm.
<figure><img src="image_21.png" /><figcaption>**Figure 17.** Evolution of $\Phi_*$ (left) and $L_*$ (right). The red circles are our IR LF measurements with the best-fit DPL functions (Table 8). We note in our measurements that all parameters in the DPL function are free in the fitting for the $z = 1.0-2.0$ bin, while the bright-end (and faint-end) slope is fixed for the other bins. For comparison, we also show the previous measurements with Herschel (Gruppioni et al. 2013; Magnelli et al. 2013; Wang et al. 2019a), SCUBA2 (Lim et al. 2020), and SCUAB2+ALMA (Koprowski et al. 2017). The magenta solid and dashed lines are the best-fit measurements presented in Koprowski et al. (2017) and their extrapolations, respectively.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 18.** Left: Cosmic IR SFRD measurements from ALCS IR LFs at $z = 0.6-7.5$. The red circles are our measurements with a conversion of SFR [$M_{\odot}\,\text{yr}^{-1}$] = $4.5 \times 10^{-44} L_{\text{IR}}$ [erg s$^{-1}$], which is the same calibration as those from Madau & Dickinson (2014), but we multiply a factor of 0.63 to apply the correction from the Salpeter IMF to the Chabrier IMF. In both panels, the filled and open red circles indicate the integration range of IR LFs down to our detection limit of $L_{\text{IR}} = 10^{10} L_{\odot}$ and the potential turnover or flattening IR luminosity of $L_{\text{IR}} = 10^8 L_{\odot}$ (see Section 6.3). The black solid and dashed lines are two empirical models of “Dust Rich” and “Dust Poor” scenarios presented in Casey et al. (2018a), respectively. The gray shaded region denotes previous measurements of Zavala et al. (2021) that perform a backward evolution model to constrain the evolution of the IR LFs to reproduce the observed number counts at 1.2 and 3 mm. The gray line is the prediction from the SHARK simulation (Lagos et al. 2018, 2020), where the integrated $L_{\text{IR}}$ range is matched with our measurements. Right: Cosmic total (= IR + UV) SFRD measurements. The black line shows the previous best estimate of the empirical model (Equation (9)) fitting to the SFRD$_{\text{total}}$ and SFRD$_{\text{corr}}$ measurements presented in Madau & Dickinson (2014). The blue and brown shaded regions present the $1\sigma$ range of the previous SFRD$_{\text{UV}}$ estimates without and with the dust correction, respectively, drawn by fitting the empirical model (Equation (9)) to the measurements presented in Bouwens et al. (2020). The red filled circles represent our SFRD$_{\text{total}}$ measurements estimated from our SFRD$_{\text{IR}}$ and the previous SFRD$_{\text{UV}}$ estimate without the dust correction (Bouwens et al. 2020). The red line and shaded region are the best-fit empirical model (Equation (9)) for our SFRD$_{\text{total}}$ measurements and its $1\sigma$ range. For presentation purposes, we slightly shift the open red circles along the abscissa.</figcaption></figure>

We find that our measurements are overall consistent with the previous measurements and the prediction of the simulation within the errors. At $z > 2$, our measurements fall between the two possible scenarios in Casey et al. (2018b). Although our measurements are generally higher than those of Zavala et al.

(2021), this is probably because the fitting at the faint regime ($S_{1.2\text{mm}} < 0.1$ mJy) in their backward model relies on the ASPECS results that may be obtained from the underdense region (Section 5.1). Our measurements are well consistent with the prediction from the simulation of Lagos et al. (2020).
In the right panel of Figure 18, we also show the evolution of the total (= IR + UV) SFRD $\psi_{SFRD}$ (IR+UV) by combining our measurements and the unobscured SFRD estimates without the dust extinction correction presented in Bouwens et al. (2020). The red line and shaded region denote the best-fit evolution and its $1\sigma$ uncertainty range by using an empirical function form in Madau & Dickinson (2014) given by

$$ \psi(z) = A \frac{(1+z)^a}{1 + [(1+z)/B]^b}, \qquad (9) $$

where the evolution is proportional to $(1+z)^a$ at low redshift, and $(1+z)^{a-b}$ at high redshift, with A and B as coefficients. Because our measurements do not have the constraint at $z < 0.6$, we calculate $\psi(z=0)$ with the best-fit parameters presented in Madau & Dickinson (2014) and include it in our fitting by assuming its error bar scale similar to our measurements. We obtain the best-fit parameters of $a = 2.8 \pm 0.5$, $b = 6.6 \pm 0.5$, $A = 0.010 \pm 0.003$, and $B = 3.3 \pm 0.3$ for our SFRD measurements.

Overall, we find that our measurements are generally consistent with the previous total SFRD measurements within the errors. At $z > 2$, although we find that the central values in our measurements are higher than the previous measurements in both Madau & Dickinson (2014) and Bouwens et al. (2020), the optical–NIR-dark objects that have been recently reported even up to $z \simeq 7$ (e.g., Fujimoto et al. 2016; Franco et al. 2018; Wang et al. 2019b; Casey et al. 2019; Williams et al. 2019; Yamaguchi et al. 2019; Romano et al. 2020; Fudamoto et al. 2021; Talia et al. 2021; Fujimoto et al. 2022; Gómez-Guijarro et al. 2022; Manning et al. 2022; Barrufet et al. 2023b; Xiao et al. 2023b; Giulietti et al. 2023; Rodighiero et al. 2023) are not included in the optical–NIR surveys used in these previous measurements. Our blind survey for the millimeter sources includes the optical–NIR-dark objects and thus provides a complete estimate of the total SFRD measurement, which is likely the cause of the slight excess in our estimate. Similarly, at $z \lesssim 2$, we also find that the central values in our measurements are higher than the previous studies. Although Madau & Dickinson (2014) include those measurements obtained from the IR surveys that unlikely miss the optical–NIR-dark objects, those measurements are obtained with a flat faint-end slope assumption in the IR LFs that has not been well constrained in the previous Herschel studies at $z > 0.3$. Therefore, the potential difference at $z \lesssim 2$ is attributed to the faint-end slope in the IR LFs.

Note that these cosmic SFRD estimates can vary with the integration range. In Section 5.4, we find that the faint-end slope in the 1.2 mm number counts may have a turnover or flattening ($S_{turn}$) above $\sim 0.2$ µJy not to make the integrated 1.2 mm flux density exceed the upper limit of the CIB. Given that $S_{turn} = 0.2$ µJy is $\sim 2$ dex deeper than our detection limit in the direct counts of $S_{1.2mm} \sim 7$ µJy, we assume the turnover or flattening is also taking place in the IR LFs at $\sim 2$ dex deeper than the detection limit of $L_{IR} = 10^{10} L_{\odot}$. We thus regard $L_{IR} = 10^8 L_{\odot}$ as the turnover or flattening IR luminosity $L_{IR, turn}$ and compute the integration of IR LFs down to $L_{IR, turn}$, below which the integrated values are thought to be negligible. In Figure 18, the open circles present our SFRD measurements with $L_{IR, turn}$. We find that the integrated values are not increased more than a factor of $\sim 2$ in the entire redshift range, and at $z > 4$, the increase is almost negligible. At $z > 4$, the total cosmic SFRD measurements with $L_{IR, turn}$ are estimated to

be $161_{-21}^{+25}\%$ of the previous measurement of Madau & Dickinson (2014). Our results indicate that the general understanding of the cosmic SFRD is unchanged by a factor of $\sim 2$, while there may be an additional ($\sim 60\%$) SFRD component contributed by the faint millimeter population, including NIR-dark objects.

In addition to the NIR-dark objects, recent ALMA observations have unveiled the presence of highly obscured star formation even in low-mass UV-selected galaxies at $z = 4-7$ (e.g., Akins et al. 2022; Mizener et al. 2024). These results suggest that the previous CSFRD estimates, generally relying on UV-selected sources with dust correction based on the UV and continuum slope relation (e.g., Bouwens et al. 2014), may be underestimated. Although the additional SFRD component could potentially introduce further tension with the cosmic stellar mass density (SMD) measurements at $0 \le z \le 3$ (e.g., Madau & Dickinson 2014; Leja et al. 2015), these discrepancies may be explained by various causes such as an overestimate in SFR, underestimate in stellar mass, uncertainty in the faint end of UVLF and/or SMF, unknown IMF shapes at high $z$, and the mass-loss factor of stellar mass within galaxy (see, e.g., Section 5.3 in Madau & Dickinson 2014). Therefore, although actual contributions of NIR-dark galaxies to the CSFRD need to be further quantified by systematic ALMA follow-up observations for more statistics and/or more accurate measure of $L_{IR}$ via multiband observations, our findings of a potential $\sim 60\%$ excess in CSFRD measurements may be consistent with this recent discovery of NIR-dark objects and do not violate the cosmic SMD. However, we acknowledge that our results could be affected by the small number statistics, which is propagated to the uncertainty in our final CSFRD estimates. Therefore, we emphasize that our SFRD measurements are still consistent with previous measurements within the $1\sigma-2\sigma$ errors.

## 7. Summary

In this paper, we identify 180 ALMA continuum sources with a flux density of $\sim 0.007-3.8$ mJy in the 33 lensing clusters homogeneously observed in the large ALMA project of ALCS. The ALCS explores the faintest layer of millimeter wavelengths most comprehensively and efficiently with gravitational lensing support, resulting in a large, blind sample of faint millimeter sources out to $z \sim 6$. This allows for the identification of unique sources, such as a strongly lensed and multiply imaged ($\mu \sim 30$ and $\sim 160$) galaxy at $z = 6.07$, as well as NIR-dark populations with no counterparts in deep HST and even Spitzer images. Using this large sample, we derive the 1.2 mm number counts and discuss the contributions of these faint sources to CIB light, including potential uncertainties in the faint-end shape. We also derive the IR LFs and the total (obscured and unobscured) cosmic SFRD at $z = 0.6-7.5$, in conjunction with the recent serendipitous detection of millimeter sources at $z > 7$ in the literature. The main findings of this paper are summarized as follows:

1. The target clusters were selected from the best-studied 33 massive galaxy clusters drawn from HST treasury programs of HFF, CLASH, and RELICS, where the deep HST and Spitzer maps are available. The observations were carried out with two frequency setups in ALMA Band 6, spanning a total of 15 GHz width with the center at 261.25 GHz ($\sim 1.15$ mm) to enlarge the blind
identifications of the line emitters, in addition to the continuum sources. The total survey area above a relative sensitivity of 30% in the mosaic map covers 133 arcmin² in the observed frame. After the lens correction at z = 6, this is decreased down to 20 and 4 arcmin² at detection limits of 0.1 and 0.04 mJy, which are larger than those in ASPECS at the same detection limits by factors of ~4.8 and ~2.0, respectively.

2. We produce natural-weighted and uv-tapered maps ($2''0 \times 2''0$) for each cluster not to miss strongly elongated millimeter sources due to the lensing effect. Based on the positive and negative source histogram as a function of peak S/N in the maps, we adopt thresholds of S/N = 5.0 and 4.5 for the natural-weight and uv-tapered maps, respectively, to maintain a high purity >0.99. This yields 141 ALMA sources, referred to as the blind sample. From crossmatching with counterparts in the Spitzer/IRAC maps, we also identify 39 additional ALMA sources with S/N down to 4.0 in the natural-weighted maps. Among the blind sample, we find 18 ALMA sources with no HST counterparts, some showing no counterparts even in Spitzer/IRAC ch1 and ch2 maps.

3. After the lens correction and removing the dust emission from the cluster member galaxies, we obtain 146 ALCS sources, which span in the intrinsic flux density range of ≈0.007–3.8 mJy and $z \simeq 0–6$ with the median values of 0.35 mJy and $z = 2.03$ and are the most widely distributed in both parameter spaces among the sources identified in recent ALMA surveys. These median values are consistent with the general trend of the increasing median flux density at submillimeter/millimeter wavelengths in the increasing median redshift. The median magnification factor in our ALCS sources is estimated to be 2.7.

4. With careful corrections for the flux measurements and completeness via dedicated MC simulations based on a constant FIR size assumption, we derive the 1.2 mm number counts by using the blind sample, exploring its faint end down to ~7 µJy. We find that the number counts are well represented by both Schechter and DPL functions and generally consistent with previous measurements in a wide flux range over a ~3 dex scale. We estimate a faint-end slope of $\alpha = -2.05_{-0.10}^{+0.12}$ for our fiducial model. These results suggest little flattening in the counts below ~0.1 mJy as previously claimed by earlier ALMA studies.

5. To assess the sensitivity of the completeness correction on the intrinsic source size, we derive 1.2 mm number counts by assuming a positive or negative correlation of FIR size with IR luminosity $L_{IR}$ as reported in the literature and evaluate the impact of this on the form of the faint number counts. We find that the faint-end slope varies depending on the FIR size assumptions, potentially changing the CIB resolved fraction by factors of ~0.6–1.1.

6. We also test the potential impact on the derived number counts of the varying precision of the lens model among the clusters by separately estimating the 1.2 mm counts for the HFF, CLASH, and RELICS clusters. The highest quality lens models, those for HFF, yield the shallowest faint-end slope, $\alpha = -1.71_{-0.29}^{+0.34}$ (but with a larger cosmic variance uncertainty). This might indicate that the high-magnification areas in the CLASH and RELICS clusters

are underestimated due to their less detailed mass models with fewer multiple images, which decreases the granularity and small-scale structures included in the mass models. Combining all the clusters weighted by the quality of the lens models yields a slightly shallower value of $\alpha = -1.96_{-0.13}^{+0.14}$. We also note the estimate of the faint-end slope excluding the highly magnified sources ($\mu \gtrsim 3–10$) also yields $\alpha \sim -2.0$, similar to our fiducial model. Hence, despite the potential remaining uncertainties, all faint-end slope estimates agree within ~1$\sigma$–2$\sigma$, suggesting a conservative $\alpha$ range of $-2.0_{-0.3}^{+0.2}$, regardless of the clusters used or the different FIR size assumptions.

7. We integrate our fiducial Schechter function fit (which assumes a constant FIR size and no weights among clusters) to the 1.2 mm counts and obtain an integrated intensity of $20.7_{-6.5}^{+8.5}$ Jy deg$^{-2}$ down to a flux limit of 7 µJy. This integrated flux density corresponds to ~80% of the CIB measured by COBE/FIRAS assuming the model for Galactic foreground emission of Odegard et al. (2019). Extrapolating our fiducial Schechter function to $\lesssim 7$ µJy will exceed the CIB measurements unless there is a turnover or flattening in the counts at fainter flux densities or a shallower slope of $\alpha \gtrsim -1.8$ whose integration down to zero does not diverge. However, we note that there is also significant systematic uncertainty in the CIB constraint due to the subtraction of the foreground emission.

8. We calculate $L_{IR}$ of our ALCS sources by modeling the IR–millimeter SED with the Spitzer, Herschel, and ALMA photometry and also derive the IR LFs out to $z = 7.5$, combining our ALCS sources with the faint millimeter sources serendipitously detected in recent ALMA studies at $z > 7$. We obtain the faint-end slope of ≈0.9–1.0 in the DPL function, which is higher than previous IR LF measurements of ≈0.4–0.6 at $z \lesssim 0.3$, but close to reports in the UVLF studies of ≈1.5–2.0. The difference may be explained by the relatively large luminosity range (~2–3 dex) covered by our IR LFs owing to the lensing support, compared to the previous studies, and/or by the different triggering mechanisms of dusty starbursts in the local universe due to the significant decrease of the gas fraction. Our ALCS IR LF results are generally consistent with the predictions from simulations within ~1$\sigma$–2$\sigma$.

9. We integrate the IR LFs and derive the redshift evolution of the obscured and total (= obscured + unobscured) SFRD. We find that our total cosmic SFRD measurement at $z = 4–8$ is estimated to be $161_{-21}^{+25}$% of the previous measurements of Madau & Dickinson (2014), where the measurements at $z > 4$ relied only on optical–NIR surveys. This indicates that our current understanding of the cosmic SFRD might be underestimated by a factor of ≈1.6 due to missing the SFRD component contributed by both the very faint millimeter populations and NIR–dark galaxies.

Acknowledgments

We thank the anonymous referee for the thorough review, insightful comments, and suggestions. We thank Carlos Gómez-Guijarro and Aswin Vijayan for sharing their measurements and
the compilation of the literature. We are also grateful to Caitlin Casey for useful discussions on the evolution of the IR LFs, Steven Finkelstein, and Peter Behroozi for the useful comments on the cosmic SFR and SMD evolutions, and Soh Ikarashi for the helpful comments on the impact of the size assumptions to the faint end of the 1.2 mm number counts. This paper makes use of the ALMA data: ADS/JAO; ALMA Nos. 2018.1.00035. L, 2019.1.00237.S, 2019.1.00949.S, 2021.1.01246.S, and 2021.1.00407.S. ALMA is a partnership of the ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by the ESO, AUI/NRAO, and NAOJ. This work is based on observations and archival data made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA, along with archival data from the NASA/ESA Hubble Space Telescope. This research also made use of the NASA/IPAC Infrared Science Archive (IRSA), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This work was supported in part by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan, and JSPS KAKENHI grant Nos. 24H00004, 22H01260, 20H05856, 20H00181, 17H06130, 22K21349, 20H00180, 21H04467, 20H01953, 22KK0231, and 23K20240. The Cosmic Dawn Center is funded by the Danish National Research Foundation under grant No. 140. This work was supported by the joint research program of the Institute for Cosmic Ray Research (ICRR), University of Tokyo. This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 847523 "INTERACTIONS" and from NASA through the NASA Hubble Fellowship grant No. HST-HF2-51505.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for

Research in Astronomy, Incorporated, under NASA contract NAS5-26555. S.F. acknowledges support from the European Research Council (ERC) Consolidator Grant funding scheme (project ConTExt, grant No. 648179). K. Kohno acknowledges support from the NAOJ ALMA Scientific Research grant No. 2017-06B and JSPS KAKENHI grant Nos. JP17H06130, JP22H04939, and JP23K20035. G.E.M. acknowledges the Villum Fonden research grants 13160 and 37440. I.R.S. acknowledges support from dSTFC (ST/T000244/1). F.E.B. acknowledges support from ANID-Chile BASAL CATA FB210003, FONDECYT Regular 1200495 and 1190818, and Millennium Science Initiative Program-ICN12_009. D.E. acknowledges support from a Beatriz Galindo senior fellowship (BG20/00224) from the Spanish Ministry of Science and Innovation, projects PID2020-114414GB-100 and PID2020-113689GB-100 financed by MCIN/AEI/10.13039/501100011033, project P20-00334 financed by the Junta de Andalucía, and project A-FQM-510-UGR20 of the FEDER/Junta de Andalucía-Consejería de Transformación Económica, Industria, Consejería Universidades. K. Kirsten acknowledges support from the Knut and Alice Wallenberg Foundation (2017.0292 and 2019.0443) and the Swedish Research Council (2015-05580).

**Software:** CASA (v5.4.0, v5.6.1; CASA Team et al. 2022), SOURCEXTRACTOR (v2.5.0; Bertin & Arnouts 1996), EAZY (Brammer et al. 2008), GLAFIC (Oguri 2010), LENSTOOL (Jullo et al. 2007), LTM (Zitrin et al. 2015).

## Appendix A
Serendipitous and Follow-up Line Detection

We extract the ALMA Band 6 spectra with a 1''0 radius aperture for the 180 ALCS sources with the 15 GHz ALCS data cube and identify the serendipitous line emission from 13 ALCS continuum sources (Section 3.4). Three lines are the CO emission line from the brightest cluster galaxies (BCGs) in the clusters, and the remaining 10 lines are interpreted as [C II] and CO emission lines from the sources behind the clusters. We

<table><caption>**Table 9**
Serendipitous Line Detection from ALCS Continuum Sources in ALMA Band 6</caption><thead><tr><th>ID<br>(1)</th><th>S/N<br>(2)</th><th>ν<br>(GHz)<br>(3)</th><th>FWHM<br>(km s<sup>-1</sup>)<br>(4)</th><th>I<sub>line</sub><br>(Jy km s<sup>-1</sup>)<br>(5)</th><th>f<sub>cont</sub><br>(%)<br>(6)</th><th>Line<br>(7)</th><th>z<br>(8)</th><th>Note<br>(9)</th></tr></thead><tbody><tr><td>M1931-C41</td><td>46.0</td><td>255.67 ± 0.00</td><td>224 ± 11</td><td>2.05 ± 0.14</td><td>1.7</td><td>CO(3–2)</td><td>0.352</td><td>BCG</td></tr><tr><td>A370-C31</td><td>18.9</td><td>251.75 ± 0.00</td><td>168 ± 10</td><td>0.63 ± 0.05</td><td>12.0</td><td>CO(3–2)</td><td>0.375</td><td>BCG</td></tr><tr><td>M0035-C94</td><td>9.5</td><td>254.40 ± 0.02</td><td>332 ± 50</td><td>0.48 ± 0.1</td><td>10.9</td><td>CO(3–2)</td><td>0.359</td><td>BCG</td></tr><tr><td>R0600-C164</td><td>9.1</td><td>268.68 ± 0.01</td><td>162 ± 30</td><td>0.32 ± 0.08</td><td>4.9</td><td>[C II]</td><td>6.072</td><td>HST-detected (F21, L21)</td></tr><tr><td>M0553-C133</td><td>15.2</td><td>269.07 ± 0.02</td><td>557 ± 62</td><td>1.53 ± 0.22</td><td>3.2</td><td>CO(5–4)</td><td>1.142</td><td>HST-detected, multiple images (E17)</td></tr><tr><td>M0553-C190</td><td>32.8</td><td>269.08 ± 0.02</td><td>562 ± 55</td><td>3.36 ± 0.43</td><td>4.1</td><td>CO(5–4)</td><td>1.142</td><td>HST-detected, multiple images (E17)</td></tr><tr><td>M0553-C249</td><td>24.2</td><td>269.06 ± 0.02</td><td>592 ± 62</td><td>2.14 ± 0.29</td><td>4.0</td><td>CO(5–4)</td><td>1.142</td><td>HST-detected, multiple images (E17)</td></tr><tr><td>M0416-C120</td><td>5.4</td><td>267.30 ± 0.02</td><td>157 ± 55</td><td>0.25 ± 0.12</td><td>6.7</td><td>CO(8–7)</td><td>2.448</td><td>HST-detected, z<sub>phot</sub> + single line</td></tr><tr><td>R0600-C13</td><td>7.8</td><td>253.76 ± 0.02</td><td>325 ± 45</td><td>0.49 ± 0.09</td><td>4.5</td><td>CO(5–4)</td><td>1.271</td><td>HST-detected, z<sub>phot</sub> + single line</td></tr><tr><td>M0553-C303</td><td>6.2</td><td>250.92 ± 0.02</td><td>260 ± 56</td><td>0.34 ± 0.1</td><td>3.4</td><td>CO(5–4)</td><td>0.836</td><td>HST-detected, z<sub>phot</sub> + single line</td></tr><tr><td>R0032-C127</td><td>9.2</td><td>271.74 ± 0.04</td><td>674 ± 119</td><td>0.83 ± 0.19</td><td>1.6</td><td>CO(8–7)</td><td>2.391</td><td>HST-detected, two lines,<sup>a</sup> multiple images</td></tr><tr><td>R0032-C131</td><td>6.3</td><td>272.04 ± 0.02</td><td>210 ± 54</td><td>0.28 ± 0.1</td><td>0.9</td><td>CO(8–7)</td><td>2.391</td><td>HST-detected, two lines,<sup>a</sup> multiple images</td></tr><tr><td>R0032-C198</td><td>5.0</td><td>271.93 ± 0.08</td><td>551 ± 219</td><td>0.33 ± 0.17</td><td>2.1</td><td>CO(8–7)</td><td>2.391</td><td>HST-detected, two lines,<sup>a</sup> multiple images</td></tr></tbody></table>

**Notes.** Column (1): ALCS continuum source ID. Column (2): S/N of the 1''0 diameter aperture photometry in the velocity-integrated map for the line emission. Column (3): central frequency of the line emission. Column (4): FWHM of the line width. Column (5): velocity-integrated line flux. Column (6): contribution of the line flux to the 1.2 mm continuum flux density estimate, which is corrected in Section 3.3. Column (7): line classification. Column (8): spectroscopic redshift from the line. Column (9): references: E17 (Ebeling et al. 2017), F21 (Fujimoto et al. 2021), and L21 (Laporte et al. 2021). The values in columns (3)–(5) are measured from a single Gaussian fitting in the spectrum.

<footnote><sup>a</sup> The second line is securely detected in follow-up Band 3 line-scan observations (Table 10) at a consistent frequency from these three sources.</footnote>
<figure><img src="image_23.png" /><figcaption>**Figure 19.** ALCS Band 6 continuum-subtracted spectra for the 10 serendipitous line detection in the ALCS continuum sources. The S/N, central frequency, and line width are summarized in Table 9.</figcaption></figure>

<table><caption>**Table 10**
Line Detection in Follow-up Line Scan in ALMA Band 3, 4 for NIR-dark Galaxies in ALCS</caption><thead><tr><th>ID<br>(1)</th><th>z<br>(2)</th><th>Band<br>(3)</th><th>Line<br>(4)</th><th>Note<br>(5)</th><th>References<br>(6)</th></tr></thead><tbody><tr><td>M0553-C67</td><td>4.932</td><td>3</td><td>CO(5–4)</td><td>No counterparts in HST/WFC3, z<sub>phot</sub> + single line</td><td>App. A</td></tr><tr><td>R0032-C127</td><td>2.391</td><td>3, 4</td><td>CO(3–2), CO(4–3), [C I](1–0)</td><td>Two lines, multiple images</td><td>T23</td></tr><tr><td>R0032-C131</td><td>2.391</td><td>3, 4</td><td>CO(3–2), CO(4–3), [C I](1–0)</td><td>Two lines, multiple images</td><td>T23</td></tr><tr><td>R0032-C198</td><td>2.391</td><td>3, 4</td><td>CO(3–2), CO(4–3), [C I](1–0)</td><td>Two lines, multiple images</td><td>T23</td></tr><tr><td>R0032-C208</td><td>2.985</td><td>3, 4</td><td>CO(3–2), CO(5–4)</td><td>No counterparts in HST/WFC3, two lines, multiple images</td><td>T23</td></tr><tr><td>R0032-C281</td><td>2.985</td><td>3, 4</td><td>CO(3–2), CO(5–4)</td><td>No counterparts in HST/WFC3, two lines, multiple images</td><td>T23</td></tr><tr><td>R0032-C304</td><td>2.985</td><td>3, 4</td><td>CO(3–2), CO(5–4)</td><td>No counterparts in HST/WFC3, two lines, multiple images</td><td>T23</td></tr><tr><td>M0417-C46</td><td>3.652</td><td>3, 4</td><td>CO(4–3), CO(6–5)</td><td>No counterparts in HST/WFC3, two lines, multiple images</td><td>T23</td></tr><tr><td>M0417-C58</td><td>3.652</td><td>3, 4</td><td>CO(4–3), CO(6–5)</td><td>No counterparts in HST/WFC3, two lines, multiple images</td><td>T23</td></tr><tr><td>M0417-C121</td><td>3.652</td><td>3, 4</td><td>CO(4–3), CO(6–5)</td><td>No counterparts in HST/WFC3, two lines, multiple images</td><td>T23</td></tr><tr><td>R0600-C67</td><td>(4.80)</td><td>3, 4</td><td>CO(7–6)?, CO(6–5)?</td><td>No counterparts in HST, IRAC, two tentative line + z<sub>phot</sub><sup>a</sup></td><td>B23</td></tr><tr><td>R0949-C10</td><td>5.753</td><td>3, 4</td><td>CO(6–5), CO(5–4)</td><td>No counterparts in HST, IRAC, two lines</td><td>B23</td></tr><tr><td>R0949-C19</td><td>(5.53)</td><td>3, 4</td><td>CO(5–4)?</td><td>No counterparts in HST, IRAC, single tentative line + z<sub>phot</sub><sup>a</sup></td><td>B23</td></tr></tbody></table>

**Notes.** Column (1): ALCS continuum source ID. Column (2): redshift determined by the emission lines. Column (3): ALMA band used in the follow-up line-scan observations. Column (4): lines detected in the follow-up observations. Column (5): notes. Column (6): references: T23 (Tsujita et al. 2024) and B23 (F. Bauer et al. 2024, in preparation).

<sup>a</sup> The upper limit of the flux ratio at 4.5 µm(IRAC ch2) and 1.2 mm (ALMA) suggests the source redshift at $z > 4.5–5.0$ based on the SED template of the high-z dusty galaxy. We regard the line feature observed with the highest S/N (~5) in the Band 3 + 4 spectrum as the CO line, which satisfies the lowest-z solution at $z > 4.5$ (Appendix A).

summarize these 13 source properties in Table 9. The ALMA Band 6 line spectra for three sources are presented in Figure 4, and Figure 19 presents the remaining 10 sources.

Several ALMA CO line scan follow-up programs have been ongoing, and we additionally detect emission lines from six ALCS sources, which is summarized in Table 10. Figure 20 shows an ALMA Band 3 line spectrum for an NIR-dark lensed object of M0553-C67 that has been observed in a follow-up line scan program. We provide further details for each line emission as follows.

**M1931-C41 (serendipitous).** The emission line is serendipitously identified in the ALCS data cube at 255.69 GHz with S/N = 46.0. The source position corresponds to one of the BCGs in M1931, and the emission line is determined as CO(3–2) at $z = 0.352$, which is in excellent agreement with the cluster redshift of M1931.

**A370-C31 (serendipitous).** The emission line is serendipitously identified in the ALCS data cube at 251.77 GHz with S/N = 18.9 in the velocity-integrated line flux density over 210 km s<sup>-1</sup>. The source position corresponds to one of the BCGs in A370, and the emission line is determined as CO(3–2) at $z = 0.375$, which is in excellent agreement with the cluster redshift of A370.

**M0035-C94 (serendipitous).** The emission line is serendipitously identified in the ALCS data cube at 254.43 GHz with S/N = 9.5 in the velocity-integrated line flux density over 420 km s<sup>-1</sup>. The source position corresponds to one of the BCGs in M0035, and the emission line is determined as CO(3–2) at $z = 0.359$, which is in excellent agreement with the cluster redshift of M0035.

**R0600-C13 (serendipitous).** The emission line is serendipitously identified in the ALCS data cube at 253.77 GHz with
<figure><img src="image_24.png" /><figcaption>**Figure 20.** Zoom-in spectra of the 27 GHz wide (~85–112 GHz) follow-up ALMA Band 3 line cube (left) and the ALCS 15 GHz wide Band 6 line cube (right) for a NIR-dark lensed galaxy of M0553-C67. The Band 3 follow-up line-scan observations yield a robust single line detection (S/N = 12.1) at ~97.1 GHz. Based on the photometric redshift of $z_{phot} \gtrsim 4$ and the nondetection of neighboring lines, we regard the line as CO(5–4) at z = 4.93 in this paper (see Appendix A). The blue and green dashed vertical lines indicate the line frequency of [C I](1-0) at z = 3.74 and [C II]158µm at z = 6.12 in case the 97.1 GHz line corresponds to CO(4–3) and CO(6–5), respectively.</figcaption></figure>

S/N = 7.8 in the velocity-integrated line flux density over 420 km s⁻¹. The source falls outside of the HST/WFC3 images, but the template fitting with another multiwavelength suggests $z_{phot} = 1.00_{-0.05}^{+0.05}$ (Appendix B). An independent SED analysis based on the FIR photometry with MAGPHYS suggests $z_{phot} = 1.23_{-0.45}^{+0.40}$. We thus regard this line as CO(5–4) at z = 1.271, while this interpretation might be updated with a second-line detection in future follow-up observations.

**R0600-C67 (follow-up).** Although the ALMA continuum is very brightly and securely detected (~2.3 mJy, S/N = 62), M0600-C67 shows no counterparts in HST/WFC3 and even IRAC maps (F160W $\gtrsim$27 mag, IRAC $\gtrsim$25 mag, $3\sigma$; Coe et al. 2019). Because of its uniquely NIR-faint nature even in the IRAC bands, despite the brightness in the ALMA Band 6, a follow-up ALMA Band 3 + 4 line-scan observation was carried out (No. 2021.1.00407.S, PI: F. Bauer). Only marginal line features are detected (S/N = 4.6–5.0), where no convincing redshift solutions are found among the combinations of these marginal line features. We thus regard the marginal line feature in the Band 3 + 4 spectrum with the highest S/N of 5.0 detected at 138.979 GHz as CO(7-6) at z = 4.80 and use this redshift interpretation in our analysis throughout the paper. However, future follow-up observations might update this interpretation with a secure second-line detection. Further details of the follow-up line-scan observation results will be presented in F. Bauer et al. (2024, in preparation).

**R0949-C10 and C19 (follow-up).** Similar to R0600-C13, R0949-C10 and R0949-C19 are characterized as very bright in ALMA Band 6 (~5.5 mJy and ~1.6 mJy, S/N = 62 and 13), but relatively very faint in NIR with no counterparts in HST/WFC3 and even IRAC maps (F160W $\gtrsim$27 mag, IRAC $\gtrsim$25 mag, $3\sigma$; Coe et al. 2019). The upper limits of the flux ratio at 4.5 µm (IRAC ch2) and 1.2 mm (ALMA) indicate their source redshifts at z > 5.5 and z > 4.5 in R0949-C10 and R0949-C19, respectively, based on the composite SED of the high-z dusty galaxies (Dudzevičiūtė et al. 2020). A follow-up ALMA Band 3 + 4 line-scan observation was conducted (No. 2021.1.00407.S, PI: F. Bauer). The multiple lines are successfully detected in R0949-C10, determining the source redshift of z = 5.753, which meets the above lower limit of the

redshift. On the other hand, only marginal line features are detected in R0949-C19 (S/N = 4.3–4.7), where no convincing redshift solutions are found among the combinations of these marginal line features. We thus regard the marginal line feature in the Band 3 + 4 spectrum with the highest S/N of 4.7 detected at 88.266 GHz as CO(5–4) at z = 5.53 and use this redshift interpretation for R0949-C19 in our analysis throughout the paper. However, future follow-up observations might update this interpretation with a secure second-line detection. Further details of the follow-up line-scan observation results will be presented in F. Bauer et al. (2024, in preparation).

**M0416-C120 (serendipitous).** The line is serendipitously identified in the ALCS data cube at 267.34 GHz with S/N = 5.3 in the velocity-integrated line flux density over 240 km s⁻¹. The source has a counterpart in the HST/F160W map and IRAC/ch1 and ch2, while it is not detected in any Herschel bands. The SED fitting with EAZY suggests $z_{phot} = 9.62_{-0.06}^{+0.07}$ (Kokorev et al. 2022). Although our ALCS 15 GHz data cube has the potential to detect bright FIR emission lines of [C II] 158µm at z = 5.97–6.17 and 6.38–6.60 as well as [O III] 88 µm at z = 11.45–11.80 and 12.18–12.57, the redshift solution falls outside of all of these redshift ranges. Moreover, the morphology in the HST/F160W map implies a large system (Figure 26), which deviates from typical galaxy sizes at z > 9 (e.g., effective radius of <0.′′1; e.g., Kawamata et al. 2018) even with the lensing effect ($\mu$ = 1.6 at z = 9.6) taken into account. We thus interpret that the Lyman-break and the 4000 Å break features are degenerated in the SED fitting, and $z_{phot} \sim 2.2$, instead of beyond 9. A public HST+IRAC catalog in the same field also shows a consistent estimate of $z_{phot} = 2.34 \pm 0.20$ for this source (Shipley et al. 2018). We thus regard this line as CO(8–7) at 2.448, while this interpretation might be updated with a second-line detection in future follow-up observations.

**M0553-C67 (follow-up).** Although the ALMA continuum is brightly and securely detected (~0.8 mJy, S/N = 10), M0553-C67 shows no HST/WFC3 counterparts, but with a faint counterpart in IRAC ch1 and ch2 (F160W $\gtrsim$27 mag, IRAC ~24 mag; Kokorev et al. 2022). The template fitting with the IRAC and ALMA photometry (Section 3.5) suggests its photometric redshift of $z_{phot} = 4.2 \pm 0.3$, while it could be
higher because the IRAC photometry might be overestimated due to a nearby bright object. A follow-up ~27 GHz wide ALMA Band 3 line-scan observation was carried out (No. 2019.1.00237.S, PI: S. Fujimoto), and we reduced the data in the same manner as the ALCS data cubes (Section 3). A robust single line is detected at 97.140 ± 0.003 GHz with a line width of FWHM = 224 ± 25 km s⁻¹ and S/N = 12.1 in the velocity-integrated map. In Figure 20, we show the robust single-line detection from M0553-C67 in the follow-up Band 3 observations. Given the $z_{phot}$ estimate and nondetection of the secondary lines in the ~27 GHz wide Band 3 data cube, we interpret the line at 97.1 GHz as CO(4–3) at z = 3.74, CO(5–4) at z = 4.93, or CO(6–5) at z = 6.12. In the case of CO(4–3) at z = 3.74, [C I](1-0) should also be detected with S/N > 3 in the same Band 3 data cube, given the typical line ratio of log([C I] (1-0)/CO(4-3)) ranging in [-0.6:0.5] among high-z star-forming galaxies (Valentino et al. 2020). In the case of CO(6–5) at z = 6.12, the [C II] line should be detected with S/N > 5 in the ALCS Band 6 data cube, given the typical ratio of $L_{[C II]}/L_{IR}$. Since we do not identify such secondary line features in the Band 3 and Band 6 data cubes (Figure 20), we regard the single line at 97.1 GHz as CO(5–4) and determine its redshift of z = 4.9323. This interpretation might be updated with a second-line detection in future follow-up observations.

**R0600-C164 (serendipitous).** The line is serendipitously identified in the ALCS data cube at 268.68 GHz with S/N = 7.1 in the velocity-integrated line flux density over 150 km s⁻¹. The source has an HST counterpart with an arc shape elongated over ~6'' scale, and the S/N is increased out to 9.1 with an optimized aperture along with the arc shape (Fujimoto et al. 2021). The arc shows a clear Lyman-break feature at z ~ 6 in the HST bands. The follow-up Gemini/GMOS spectroscopy detects the continuum break and confirms the source redshift at z = 6.19<sup>+0.06</sup><sub>-0.16</sub> (Laporte et al. 2021), and the line is determined as [C II] 158 µm at z = 6.072. The two multiple images of the arc are identified with the line detection at the same frequency and the same SED properties, while the dust continuum is not detected (<4σ) from these multiple images because of different magnifications among multiple images. Three independent lens models consistently suggest that an outskirt region of the galaxy behind the cluster is crossing the caustic line, which is strongly lensed into the arc shape in the image plane with a magnification factor of 163<sup>+27</sup><sub>-13</sub> (29<sup>+4</sup><sub>-7</sub>) in the local (global) scale (Fujimoto et al. 2021). In this paper, we adopt the magnification factor of 29<sup>+4</sup><sub>-7</sub> in our analyses to estimate the intrinsic flux density of the host galaxy on the global scale, given the scopes in the number-count and IRLF measurements.

**M0553-C133/190/249 (serendipitous).** The lines are serendipitously and consistently identified in the ALCS data cube from nearby three ALCS sources at 269.12 GHz with S/N ~ 15–30 in the velocity-integrated line flux density over 660 km s⁻¹ (Figure 4). Their HST counterparts are known to be a multiply imaged system spectroscopically confirmed at z = 1.14 (Ebeling et al. 2017), and thus, these lines are determined as CO(5–4) at z = 1.14. M0553-C113 falls in the FOV of follow-up ALMA Band 3 line-scan observations for M0553-C67 (No. 2019.1.00237.S, PI: S. Fujimoto), and we also detect the CO(2-1) line from M0553-C113 in the Band 3 data cube.

**M0553-C303 (serendipitous).** The line is serendipitously identified in the ALCS data cube at 250.95 GHz with

S/N = 9.5 in the velocity-integrated line flux density over 270 km s⁻¹. The source falls outside of the HST/WFC3 footprint, while it is detected in the IRAC ch1 and ch2 maps with the forced-aperture photometry (Appendix B) and the Herschel/SPIRE maps with the deblended technique (Sun et al. 2022). The template fitting with these IRAC, Herschel, and ALMA photometry suggests $z_{phot} = 0.65_{-0.05}^{+0.05}$ (Appendix B), which is consistent with the independent FIR SED fitting results of 0.98<sup>+0.39</sup><sub>-0.39</sub> (Sun et al. 2022) within the errors. We thus regard this line as CO(4–3) at z = 0.836, while this interpretation might be updated with a second-line detection in future follow-up observations. Note that the cluster redshift of M0553 is 0.430, and M0553-C303 is a lensed galaxy behind the cluster.

**R0032-C127/131/198 (serendipitous).** The lines are serendipitously consistently identified in the ALCS data cube at ~272 GHz from nearby three ALCS sources with S/N ~5.0–9.2 in the velocity-integrated line flux density over 800 km s⁻¹. The counterparts are very faint, but identified in HST/WFC3 at their source positions, and secondary emission lines are securely detected at 101.97 GHz in a follow-up ALMA Band 3 and 4 line scan program in this field (No. 2021.1.01246.S, PI: K. Kohno). Although the central frequencies of the line features observed in the ALCS data cube are slightly different, this is likely explained by either the sensitivity or the differential magnification effect. Together with the similar NIR-dark properties, we interpret these three sources as the multiply imaged system. Based on the frequencies of these two lines, we regard them as CO(8–7) in Band 6 and CO(3–2) in Band 3 at z = 2.391. Our fiducial lens model also confirms the interpretation of the multiple images. Further details of the follow-up line-scan observation results will be presented in Tsujita et al. (2024).

**R0032-C208/281/304 (follow-up).** Although the ALMA continuum is brightly and securely detected (>1.0 mJy, S/N ≥ 10), no HST/WFC3 counterparts are identified at their source positions (F160W ≥ 27 mag, 3σ; Coe et al. 2019). The multiple lines are securely detected in follow-up ALMA Band 3 and 4 line-scan observations (No. 2021.1.01246.S, PI: K. Kohno) exactly at the same frequencies at 86.80 and 144.60 GHz. We conclude that these lines are CO(3–2) and CO(5–4) at z = 2.985 and that these three sources are the multiply imaged system, which is consistent with the predictions from our fiducial lens model. Further details of the follow-up line-scan observation results will be presented in Tsujita et al. (2024).

**M0417-C46/58/121 (follow-up).** Similar to the multiple images of R0032-C208/281/304, no HST/WFC3 counterparts are identified at their source positions (F160W ≥ 27 mag, 3σ; Coe et al. 2019), although the ALMA continuum is brightly and securely detected (>2–4 mJy, S/N > 20). A single line is robustly detected consistently at 99.11 GHz from all three sources in follow-up ALMA Band 3 line-scan observations (No. 2021.1.01246.S, PI: K. Kohno). The secondary line is also securely detected at 148.65 GHz from M0417-C46 in the Band 4 line-scan observations in the same program, while the same observations were not carried out for the other two sources. Our fiducial lens model suggests that these three sources are multiple images, and given that the line at 99.1 GHz is detected exactly at the same frequency among these three sources, we conclude that these are multiple images. Based on the two emission lines, we conclude that their redshift at
$z = 3.652$. Further details of the follow-up line-scan observation results are presented in Tsujita et al. (2024).

## Appendix B

### $z_{phot}$ Estimate by Template Fitting

In Section 3.5, we identify 26 ALCS sources lacking the HST counterparts due to their faintness in the HST bands ($N = 18$) and their sky positions outside of the HST footprints ($N = 8$), respectively. We obtain the forced-aperture photometry in IRAC ch1 and ch2 at the ALCS source position with a diameter of 3.0''0 on the residual IRAC maps after subtracting the models for the HST-detected sources as well as the deblended Herschel photometry if detected (Sun et al. 2022). There remain four ALCS sources that are not detected in any bands in IRAC and Herschel (ACT0102-C11, ACT0102-C251, R1347-C51, and M2129-C24). These four ALCS sources are all included in the primary catalog. With our purity cut for the primary catalog (>99%) based on the positive and negative source distribution (Figure 2), one of these four ALCS sources might be spurious. Still, the majority of them are most likely real.

Given the limited number of photometry data points, we simply perform the $\chi^2$ minimization with a fixed template to evaluate $P(z)$ for these 23 ($=27-4$) ALCS sources. We use the

SED fitting code of STARDUST (Kokorev et al. 2021) by implementing the SED template of Arp220 and a composite SED model from 707 dusty galaxies identified in the AS2UDS survey (e.g., Dudzevičiūtė et al. 2020). We find that the AS2UDS template suggests a lower redshift solution than Arp220 in general. This is because the AS2UDS template is fainter in the rest-frame optical than Arp220 due to a lower typical stellar mass or much higher obscuration (Dudzevičiūtė et al. 2020).

Figure 21 shows an example of the NIR–millimeter photometry and the best-fit SED (left) and $P(z)$ (right) with the Arp220 and AS2UDS templates for M0417-C121. M0417-C121 is robustly detected with ALMA ($S/N_{nat} = 20.7$) without any HST counterparts, and follow-up ALMA Band 3 and 4 line-scan observations determine its redshift at $z_{spec} = 3.652$ with two secure CO line detection (Tsujita et al. 2024). We find that these SED templates suggest the source redshift at $z \sim 3–4$, where the peak of $P(z)$ with the AS2UDS template shows an excellent agreement with $z_{spec}$ of M0417-C121. Given the intrinsically faint (~less massive) and potentially high dust obscuration aspects without any HST counterparts, the SEDs of these ALCS sources may be closer to the AS2UDS template. We thus use $P(z)$ obtained with the AS2UDS template for the 23 ALCS sources.

<figure><img src="image_25.png" /><figcaption>**Figure 21.** Example of our template fitting for the ALCS sources with no HST counterparts. Left: the best-fit SEDs for M0417-C121 based on the templates of Arp220 (black line) and the composite of 707 dusty galaxies in AS2UDS (Dudzevičiūtė et al. 2020). Right: $P(z)$ from the $\chi^2$ minimization. The red dashed line indicates $z_{spec} = 3.652$ of M0417-C121 determined from follow-up ALMA Band 3 and 4 line-scan observations with multiple line detection (Tsujita et al. 2024).</figcaption></figure>
## Appendix C
### Lens Models

We construct lens models using cluster member galaxies as well as multiple images behind clusters based on the photometric redshift, colors, and morphology of galaxies in the HST images with independent algorithms including GLAFIC (Oguri 2010; Kawamata et al. 2016, 2018; Okabe et al. 2020), LENSTOOL (Jullo et al. 2007; Richard et al. 2014; Caminha et al. 2016, 2017a, 2017b, 2019; Niemiec et al. 2020), LTM

(Zitrin et al. 2013, 2015), and dPIEeNFW (Zitrin et al. 2013, 2015). We use the lens model of GLAFIC as a fiducial model for our analyses throughout this paper, while other available models are also used to evaluate uncertainties in magnification factors (Appendix D). Note that we do not use the median magnification estimate among different lens models, because the process of choosing the median estimate systematically reduces the effective survey area at high magnifications of $\mu \gtrsim 10$ (Appendix D). In Table 11, we

<table><caption>**Table 11**
Summary of Our Fiducial and Available Lens Models in the 33 ALCS Clusters</caption><thead><tr><th>Cluster (1)</th><th>z (2)</th><th>N<sub>src,all</sub> (3)</th><th>N<sub>src,spec</sub> (4)</th><th>N<sub>img</sub> (5)</th><th>rms<sub>img</sub> (6)</th><th>References (7)</th><th>Caminha (8)</th><th>CATS (9)</th><th>Zitrin (10)</th></tr></thead><tbody><tr><td colspan="10" style="text-align:center; font-weight:bold;">HFF</td></tr><tr><td>A2744</td><td>0.308</td><td>45</td><td>24</td><td>132</td><td>0.42</td><td>K18</td><td>N</td><td>Y</td><td>Y</td></tr><tr><td>A370</td><td>0.375</td><td>49</td><td>19</td><td>135</td><td>0.50</td><td>K18</td><td>N</td><td>Y</td><td>Y</td></tr><tr><td>M0416</td><td>0.396</td><td>75</td><td>34</td><td>202</td><td>0.44</td><td>K18</td><td>Y</td><td>Y</td><td>Y</td></tr><tr><td>M1149</td><td>0.543</td><td>36</td><td>16</td><td>108</td><td>0.31</td><td>K16</td><td>N</td><td>Y</td><td>Y</td></tr><tr><td>AS1063</td><td>0.348</td><td>53</td><td>19</td><td>141</td><td>0.38</td><td>K18</td><td>Y</td><td>Y</td><td>Y</td></tr><tr><td colspan="10" style="text-align:center; font-weight:bold;">CLASH</td></tr><tr><td>A209</td><td>0.206</td><td>3</td><td>0</td><td>7</td><td>0.52</td><td>Z15</td><td>N</td><td>Y</td><td>Y</td></tr><tr><td>A383</td><td>0.187</td><td>8</td><td>6</td><td>23</td><td>0.41</td><td>Z15</td><td>N</td><td>Y</td><td>Y</td></tr><tr><td>M0329</td><td>0.45</td><td>11</td><td>10</td><td>28</td><td>0.62</td><td>Ca19, R21</td><td>Y</td><td>Y</td><td>Y</td></tr><tr><td>M0429</td><td>0.399</td><td>3</td><td>2</td><td>11</td><td>0.32</td><td>Ca19</td><td>Y</td><td>N</td><td>Y</td></tr><tr><td>M1115</td><td>0.355</td><td>3</td><td>1</td><td>9</td><td>0.44</td><td>Ca19</td><td>Y</td><td>Y</td><td>Y</td></tr><tr><td>M1206</td><td>0.439</td><td>37</td><td>37</td><td>109</td><td>0.47</td><td>Ca17, R21</td><td>Y</td><td>Y</td><td>Y</td></tr><tr><td>M1311</td><td>0.494</td><td>3</td><td>1</td><td>8</td><td>0.57</td><td>Ca19</td><td>Y</td><td>N</td><td>Y</td></tr><tr><td>R1347</td><td>0.451</td><td>40</td><td>38</td><td>129</td><td>0.65</td><td>U18, Ca19, R21</td><td>Y</td><td>Y</td><td>Y</td></tr><tr><td>M1423</td><td>0.545</td><td>3</td><td>2</td><td>12</td><td>0.65</td><td>Z15</td><td>N</td><td>Y</td><td>Y</td></tr><tr><td>M1931</td><td>0.352</td><td>7</td><td>7</td><td>19</td><td>0.39</td><td>Ca19</td><td>Y</td><td>Y</td><td>Y</td></tr><tr><td>M2129</td><td>0.57</td><td>11</td><td>11</td><td>38</td><td>0.73</td><td>Ca19</td><td>Y</td><td>Y</td><td>Y</td></tr><tr><td>R2129</td><td>0.234</td><td>7</td><td>7</td><td>22</td><td>0.43</td><td>Ca19</td><td>Y</td><td>Y</td><td>Y</td></tr><tr><td colspan="10" style="text-align:center; font-weight:bold;">RELICS</td></tr><tr><td>R0032</td><td>0.396</td><td>12</td><td>2</td><td>36</td><td>0.49</td><td>A20</td><td>Y</td><td>Y</td><td>N</td></tr><tr><td>M0035</td><td>0.352</td><td>5</td><td>0</td><td>13</td><td>0.24</td><td>...</td><td>Y</td><td>N</td><td>N</td></tr><tr><td>AC0102</td><td>0.87</td><td>10</td><td>0</td><td>28</td><td>0.52</td><td>Ce18</td><td>Y</td><td>N</td><td>N</td></tr><tr><td>M0159</td><td>0.405</td><td>4</td><td>0</td><td>10</td><td>0.47</td><td>Ci18</td><td>Y</td><td>Y</td><td>N</td></tr><tr><td>A295</td><td>0.3</td><td>6</td><td>4</td><td>18</td><td>0.22</td><td>Ci18</td><td>N</td><td>N</td><td>N</td></tr><tr><td>M0257</td><td>0.505</td><td>4</td><td>0</td><td>12</td><td>0.36</td><td>Z11</td><td>N</td><td>Y</td><td>N</td></tr><tr><td>P171</td><td>0.27</td><td>5</td><td>0</td><td>16</td><td>0.34</td><td>A18</td><td>N</td><td>N</td><td>Y</td></tr><tr><td>A3192</td><td>0.425</td><td>5</td><td>2</td><td>16</td><td>0.51</td><td>H13</td><td>N</td><td>N</td><td>N</td></tr><tr><td>M0417</td><td>0.443</td><td>21</td><td>8</td><td>57</td><td>0.45</td><td>M19 and ALCS<sup>a</sup></td><td>Y</td><td>Y</td><td>N</td></tr><tr><td>M0553</td><td>0.43</td><td>10</td><td>2</td><td>30</td><td>0.73</td><td>E17</td><td>N</td><td>Y</td><td>N</td></tr><tr><td>R0600</td><td>0.46</td><td>8</td><td>5</td><td>26</td><td>0.66</td><td>F21</td><td>Y</td><td>Y</td><td>N</td></tr><tr><td>SM0723</td><td>0.39</td><td>4</td><td>0</td><td>13</td><td>0.27</td><td>...</td><td>N</td><td>Y</td><td>N</td></tr><tr><td>R0949</td><td>0.383</td><td>4</td><td>0</td><td>11</td><td>0.26</td><td>...</td><td>N</td><td>Y</td><td>N</td></tr><tr><td>A2163</td><td>0.203</td><td>4</td><td>3</td><td>15</td><td>0.3</td><td>Ce18, Re20</td><td>Y</td><td>N</td><td>N</td></tr><tr><td>RJ2211</td><td>0.397</td><td>3</td><td>1</td><td>11</td><td>0.14</td><td>Ce18</td><td>N</td><td>N</td><td>N</td></tr><tr><td>A2537</td><td>0.297</td><td>8</td><td>1</td><td>29</td><td>0.44</td><td>Ce18</td><td>N</td><td>Y</td><td>N</td></tr></tbody></table>

**Notes.** Column (1): cluster name. Column (2): cluster redshift. Column (3): number of the multiple system identified in the cluster. Column (4): number of multiple system spectroscopically confirmed identified in the cluster. Column (5): number of multiple images identified in the image plane used for constructing our fiducial mass model with GLAFIC (Oguri 2010). Column (6): rms of positional differences between observed and model-predicted multiple images, in units of arcseconds. Column (7): references for the multiple images: K18 (Kawamata et al. 2018), Z15 (Zitrin et al. 2015), Ca19 (Caminha et al. 2019), R21 (Richard et al. 2021), Ce18 (Cerny et al. 2018), Ci18 (Cibirka et al. 2018), Z11 (Zitrin et al. 2011), M19 (Mahler et al. 2019), E17 (Ebeling et al. 2017), R20 (Rescigno et al. 2020), U18 (Ueda et al. 2018), H13 (Hsu et al. 2013), F21 (Fujimoto et al. 2021). Columns (8)–(10): Availability (Y, yes; N, no) of other mass models so far constructed by G. Caminha with LENSTOOL (Caminha et al. 2016, 2017a, 2017b, 2019), CATS team with LENSTOOL (Richard et al. 2010, 2011; Livermore et al. 2012; Newman et al. 2013; Jauzac et al. 2016; Repp et al. 2016; Dessauges-Zavadsky et al. 2017; Ebeling et al. 2017; Rexroth et al. 2017; Lagattuta et al. 2019; Mahler et al. 2019; Jauzac et al. 2021; Richard et al. 2021), and A. Zitrin with LTM and dPIEeNFW algorithms (Zitrin et al. 2013, 2015) in the ALCS 33 clusters. All model links will be available on the ALCS website (http://www.ioa.s.u-tokyo.ac.jp/ALCS/). Interested readers may find more lens models on HFF (https://archive.stsci.edu/prepds/frontier/), CLASH (https://www.stsci.edu/~postman/CLASH/), and RELICS (https://relics.stsci.edu/) websites.

<footnote><sup>a</sup> The multiple images of the ALCS sources without HST counterparts, M0417-C46/58/121, are also used for refining our fiducial lens model.</footnote>
summarize the number of the multiple images used for constructing the GLAFIC models and the rms of positional differences between the observed and model-predicted multiple images.

## Appendix D
### Lens Models Uncertainty

We compare the magnification factors between our fiducial and other lens models at random pixel positions of the maps and evaluate the systematic uncertainty from the choice of the lens model. To evaluate the systematic uncertainty, we calculate $\Delta\mu/\mu$ ($\equiv|\mu_{other} - \mu|/\mu$). The accuracy of the lens model generally depends on the richness of the multiwave-length data in the cluster field for identifying the multiple images and the cluster member galaxies. Thus, we separately calculate $\Delta\mu/\mu$ for HFF, CLASH, and RELICS. The clusters of A2744, M1115, and R0032 are used for the calculation of $\Delta\mu/\mu$ as the representatives of HFF, CLASH, and RELICS, respectively, and we fix the redshift at $z=2$ for the lens model.

In Figure 22, we show $\Delta\mu/\mu$ at random positions in the cluster fields. The median values as a function of $\mu$ are also

plotted with black squares. We find an increasing trend in $\Delta\mu/\mu$ as a function of $\mu$. We also find that $\Delta\mu/\mu$ is the smallest in HFF at $\mu \le 5$, while there is no large difference in $\Delta\mu/\mu$ among HFF, CLASH, and RELICS at $\mu > 5$. Based on these results, we adopt 20% (40%), 50%, 60%, and 80% for the systematic uncertainty of the different lens models at $\mu \le 5$, $5 < \mu \le 10$, $10 < \mu \le 30$, and $\mu > 30$ in HFF (CLASH and RELICS).

We also evaluate the difference in the effective survey area among different lens models. In Figure 23, we show the survey area after correcting the lensing effect with seven different lens models in AS1063. The redshift is fixed again at $z=2$ for all these lens models. We find that the survey area is unchanged by more than ~10% in the wide magnification ranges. This indicates that the overall shape of the magnification distribution is not significantly changed among the different lens models, although the small variations should remain. On the other hand, however, if we create the median magnification map from these seven different lens models, the survey area is underestimated at high magnifications ($\mu \gtrsim 10$). This is because such highly magnified regions are smeared out even by the small variation of the critical curve positions among the different models.

<figure><img src="image_26.png" /><figcaption>**Figure 22.** Difference of the magnification factor at random positions among different lens models. We evaluate the difference with $\Delta\mu/\mu$ and $\Delta\mu = |\mu - \mu_{other}|$, where $\mu$ and $\mu_{other}$ are the magnification factors with our fiducial and other lens models, respectively. The different color plots indicate different lens models for $z=2$, where there are 8, 3, and 4 publicly available lens models in A2744 (HFF, left), M1115 (CLASH, middle), and R0032 (RELICS, right), respectively, including our fiducial lens model constructed with GLAFIC (Oguri 2010). The black squares denote the median of $\Delta\mu/\mu$ in each $\mu$ bin.</figcaption></figure>