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

# The Evolution of the IR Luminosity Function and Dust-obscured Star Formation over the Past 13 Billion Years

J. A. Zavala¹, C. M. Casey¹, S. M. Manning¹, M. Aravena², M. Bethermin³, K. I. Caputi⁴,⁵, D. L. Clements⁶, E. da Cunha⁷, P. Drew¹, S. L. Finkelstein¹, S. Fujimoto⁵,⁸, C. Hayward⁹, J. Hodge¹⁰, J. S. Kartaltepe¹¹, K. Knudsen¹², A. M. Koekemoer¹³, A. S. Long¹⁴, G. E. Magdis⁵,⁸,¹⁵,¹⁶, A. W. S. Man¹⁷, G. Popping¹⁸, D. Sanders¹⁹, N. Scoville²⁰, K. Sheth²¹, J. Staguhn²²,²³, S. Toft⁵,⁸, E. Treister²⁴, J. D. Vieira²⁵, and M. S. Yun²⁶

¹ The University of Texas at Austin, 2515 Speedway Boulevard, Stop C1400, Austin, TX 78712, USA; jzavala@utexas.edu

² Núcleo de Astronomía, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército 441, Santiago, Chile

³ Aix Marseille Univ, CNRS, LAM, Laboratoire d'Astrophysique de Marseille, Marseille, France

⁴ Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, 9700AV Groningen, The Netherlands

⁵ Cosmic Dawn Center (DAWN), Denmark

⁶ Imperial College London, Blackett Laboratory, Prince Consort Road, London, SW7 2AZ, UK

⁷ International Centre for Radio Astronomy Research, University of Western Australia, 35 Stirling Highway, Crawley, WA 6009, Australia

⁸ Niels Bohr Institute, University of Copenhagen, Lyngbyvej 2, DK-2100 Copenhagen, Denmark

⁹ Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA

¹⁰ Leiden Observatory, Niels Bohrweg 2, 2333 CA Leiden, The Netherlands

¹¹ School of Physics and Astronomy, Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, NY 14623, USA

¹² Chalmers University of Technology, Department of Space, Earth and Environment, Onsala Space Observatory, SE-439 92 Onsala, Sweden

¹³ Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

¹⁴ Department of Physics and Astronomy, University of California, Irvine, CA 92697, USA

¹⁵ DTU-Space, Technical University of Denmark, Elektrovej 327, DK-2800 Kgs. Lyngby, Denmark

¹⁶ Institute for Astronomy, Astrophysics, Space Applications and Remote Sensing, National Observatory of Athens, GR-15236 Athens, Greece

¹⁷ Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON M5S 3H4, Canada

¹⁸ European Southern Observatory, Karl-Schwarzschild-Strasse 2, D-85748, Garching, Germany

¹⁹ Institute for Astronomy, 2680 Woodlawn Drive, University of Hawaii, Honolulu, HI 96822, USA

²⁰ California Institute of Technology, MC 249-17, 1200 East California Boulevard, Pasadena, CA 91125, USA

²¹ NASA Headquarters, 300 E. Street SW, Washington, DC 20546, USA

²² The Henry A. Rowland Department of Physics and Astronomy, Johns Hopkins University, 3400 North Charles Street, Baltimore, MD 21218, USA

²³ Observational Cosmology Lab, Code 665, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

²⁴ Instituto de Astrofísica, Facultad de Física, Pontificia Universidad Católica de Chile, Casilla 306, Santiago 22, Chile

²⁵ Department of Astronomy, University of Illinois, 1002 West Green Street, Urbana, IL 61801, USA

²⁶ Department of Astronomy, University of Massachusetts, Amherst, MA 01003, USA

Received 2020 October 15; revised 2021 January 11; accepted 2021 January 11; published 2021 March 15

## Abstract

We present the first results from the Mapping Obscuration to Reionization with ALMA (MORA) survey, the largest Atacama Large Millimeter/submillimeter Array (ALMA) blank-field contiguous survey to date (184 arcmin²) and the only at 2 mm to search for dusty star-forming galaxies (DSFGs). We use the 13 sources detected above 5σ to estimate the first ALMA galaxy number counts at this wavelength. These number counts are then combined with the state-of-the-art galaxy number counts at 1.2 and 3 mm and with a backward evolution model to place constraints on the evolution of the IR luminosity function and dust-obscured star formation in the past 13 billion years. Our results suggest a steep redshift evolution on the space density of DSFGs and confirm the flattening of the IR luminosity function at faint luminosities, with a slope of $\alpha_{LF} = -0.42^{+0.02}_{-0.04}$. We conclude that the dust-obscured component, which peaks at $z \approx 2-2.5$, has dominated the cosmic history of star formation for the past $\sim 12$ billion years, back to $z \sim 4$. At $z=5$, the dust-obscured star formation is estimated to be $\sim 35\%$ of the total star formation rate density and decreases to $25\%-20\%$ at $z=6-7$, implying a minor contribution of dust-obscured star formation in the first billion years of the universe. With the dust-obscured star formation history constrained up to the end of the epoch of reionization, our results provide a benchmark to test galaxy formation models, to study the galaxy mass assembly history, and to understand the dust and metal enrichment of the universe at early times.

*Unified Astronomy Thesaurus concepts: Galaxy evolution (594); Galaxies (573); Dust continuum emission (412); Millimeter astronomy (1061); Submillimeter astronomy (1647); Star formation (1569); High-redshift galaxies (734); Surveys (1671); Galaxy counts (588); Luminosity function (942)*

## 1. Introduction

Mapping the cosmic history of star formation is of fundamental importance in our understanding of galaxy formation and evolution not only because it contains the galaxy mass assembly history but also because it represents the footprint of the metal enrichment of the universe.

A complete, unbiased determination of the star formation rate density (SFRD) requires a multiwavelength approach to directly probe the stellar emission of newborn stars, as well as the starlight that has been absorbed (then reemitted) in dust-obscured regions. The former is efficiently measured using rest-frame ultraviolet (UV) observations and the latter using
far-infrared (FIR) and (sub)millimeter surveys that trace the dust reprocessed emission from young stars (see Madau & Dickinson 2014 for a review).

The current UV census of star formation reaches out to $z \sim 11$, close to the formation epoch of the first galaxies (Oesch et al. 2018). Nevertheless, despite some individual detections of dust in galaxies up to $z \sim 8$ (Watson et al. 2015; Laporte et al. 2017), studies of the global dust-obscured SFRD from FIR/submillimeter surveys are very limited at $z > 3$.

The highest-redshift estimations of the SFRD from FIR to millimeter surveys, at $z = 4.5–5.5$, are uncertain even at the $\approx 1.0–1.5$ dex level (Rowan-Robinson et al. 2016; Williams et al. 2019; Loiacono et al. 2020) and differ by up to two orders of magnitude (see contrasting results by Rowan-Robinson et al. 2016 and Michalowski et al. 2017). This is mostly due to the low number statistics since most of the samples used to derive these measurements range from only one (Williams et al. 2019; Loiacono et al. 2020) to a handful of objects (Dudzevičiūtė et al. 2020; Gruppioni et al. 2020). This implies that the total amount of dust-enshrouded star formation in the earliest epochs of the universe has remained unknown, and therefore our current picture of the history of cosmic star formation remains incomplete.

This study aims at estimating the history of dust-obscured star formation back to $z \sim 7$ by combining a backward galaxy evolution model of the dusty star-forming galaxy (DSFG) population (Casey et al. 2018a, 2018b) with the state-of-the-art Atacama Large Millimeter/submillimeter Array (ALMA) surveys.

We infer constraints on the prevalence and characteristics of these galaxies through measurements of galaxy number counts at different wavelengths. Despite its relative conceptual simplicity, number counts have proven to be a very powerful tool to test and constrain galaxy formation models (e.g., Baugh et al. 2005; Hayward et al. 2013) and consequently represent an essential measurement in any long-wavelength survey.

The galaxy number counts used in this work include the first arcsecond-resolution interferometric number counts at 2 mm achieved with ALMA and the most recent estimations at 1.2 and 3 mm. The former is from our new 2 mm Mapping Obscuration to Reionization with ALMA (MORA) survey described in Section 2, the first ALMA large map at this wavelength with a total area of $184 \ \text{arcmin}^2$ (an order of magnitude larger than previous contiguous interferometric blind surveys). At 3 mm, the measurements come from the $\sim 200 \ \text{arcmin}^2$ ALMA archival program reported in Zavala et al. (2018), which we revise in this work after finding three false detections in their sample (see Section 3). These number counts are also complemented with those from the $4.2 \ \text{arcmin}^2$ deep ALMA large program ASPECS survey (at 1.2 and 3 mm; González-López et al. 2019, 2020), in which 40 continuum sources were detected. All these surveys represent the state-of-the-art observations at long wavelengths, comprising not only the deepest measurements at those wave bands (beyond those assisted by gravitational amplification) but also the only interferometric surveys at 2 and 3 mm, which are crucial to constraining the infrared (IR) luminosity function (IRLF) at high redshift (see discussion by Bethermin et al. 2015; Casey et al. 2018a, 2018b; Zavala et al. 2018; see also Wilner & Wright 1997, for a pioneering $\sim 3$ mm survey).

The constraints provided by these surveys and by all the FIR and submillimeter data aggregated over the past two decades

allow us to infer the IR luminosity function and dust-obscured star formation in the past 13 billion years of the universe, pushing the current redshift limits from $z \sim 4–5$ to $z \sim 7$. Our constraints are discussed in Section 4 and summarized in Section 5.

Throughout this work we assume $H_0 = 67.7 \ \text{km s}^{-1} \ \text{Mpc}^{-1}$ and $\Omega_\lambda = 0.69$ (Planck Collaboration et al. 2016) and a Chabrier initial mass function (Chabrier 2003) for SFR estimations.

## 2. Mapping Obscuration to Reionization with ALMA (MORA) Survey

The MORA survey (ALMA project code: 2018.1.00231.S, PI: C. Casey) was originally designed to cover a contiguous area of $230 \ \text{arcmin}^2$ in two different tunings in the 2 mm band (centered at 147.3 and 139 GHz, respectively). The project was only partially observed, covering $184 \ \text{arcmin}^2$ in two separate mosaics of 156 and $28 \ \text{arcmin}^2$, respectively. The largest mosaic, which covers positions 10–20 from the original design, so-called P10–P20, is centered at $\alpha \approx 10^h 00^m 17^s$, $\delta \approx +02^\circ 22' 30''$ and covers an area of $6'6 \times 23'6$. The smaller one (P03) is centered at $\alpha \approx 10^h 00^m 44^s$, $\delta \approx +02^\circ 22' 30''$ and has $1'2 \times 23'6$ dimensions (see Figure 1). The deepest portion of the mosaics has an rms of $\sigma_{2 \ \text{mm}} = 60 \ \mu\text{Jy}$ beam, with $101 \ \text{arcmin}^2$ covered at or below the proposed map depth of $90 \ \mu\text{Jy}$ beam (see Figures 1 and 2), with a typical beam size of $\theta_{\text{FWHM}} \approx 1'^{\prime}_{\prime} 8 \times 1'^{\prime}_{\prime} 4$. This project represents the largest ALMA blank-field survey to date and the only one at 2 mm to search for DSFGs.

The signal-to-noise ratio (S/N) maps are well modeled by Gaussian statistics within $-4.5 \le S/N \le 4.5$ (see Figure 3). Above $4.5\sigma$ the excess of positive pixels comes from the detection of astronomical sources, as expected. Interestingly, there is also an excess of negative pixels with $S/N < -5$. All these pixels (identified in orange in the figure) come from the same region and are associated with a single noise peak at $-6\sigma$ (plus the effect of beam smearing; see Figure 1). Nevertheless, assuming Gaussian statistics, the probability of finding a $-6\sigma$ noise peak in our mapped area is $\lesssim 0.5\%$. An alternative explanation might be the Sunyaev–Zel'dovich (SZ) effect (Sunyaev & Zeldovich 1970) since its peak decrement is expected close to the observed frequency of our maps ($\sim 145$ GHz). A thorough discussion of the nature of this negative detection will be presented in a future work, after carrying out the required follow-up observations.

C. M. Casey et al. (2021, in preparation) will present further details of the survey along with a description of the source catalog and other physical characteristics of the 2 mm selected galaxies. Briefly, we first create a noise map using two different approaches. In the first one, we measure the standard deviation of all the pixels with primary beam responses between 0.95 and 1, which give us the minimum noise in our mosaic (or maximum depth); then, we multiply the inverse of the primary beam response map by this minimum noise vale, whose result is adopted as the noise map. In the second approach, we use a 2D boxcar-like function across the primary-beam-corrected flux mosaic (after applying a sigma-clipping procedure to remove the bright sources) to estimate the noise in each pixel, measured as the standard deviation of the pixels within the squared region of the boxcar-like function. Both procedures give us consistent results within $\sim 1\%$. Once we have a noise map, we simply divide the primary-beam-corrected flux map by this noise mosaic in order to obtain an S/N map (which is
<figure><img src="image_2.png" /><figcaption>**Figure 1.** MORA 2 mm S/N maps and the 13 sources detected above $5\sigma$ (orange squares). The noise variations across the mosaics are illustrated by the contours, which range from 60 to 240 µJy in steps of 30 µJy. The typical size of the primary beam response at the frequency of our observations ($\theta_{FWHM} \approx 43''$) is represented by the orange circle (note that the synthesized beam size of $\theta_{FWHM} \sim 1''5$ is significantly smaller). The $-6\sigma$ peak is also identified with the pale green diamond.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Cumulative distribution of $1\sigma$ survey depth and over what solid angle that rms depth is achieved. The P03 mosaic alone (covering $28 \text{ arcmin}^2$) is represented in pale green, while the P10–P20 mosaic (covering $156 \text{ arcmin}^2$) is shown in blue. The proposed map depth (90 µJy) is illustrated by the orange dashed line. With a total area of $184 \text{ arcmin}^2$ (orange line), the MORA program represents the largest ALMA blank-field survey to date and the only one at 2 mm to search for DSFGs.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Histogram of pixel values in the S/N maps. The distribution of pixel values within $-4.5 \lesssim \text{S/N} \lesssim 4.5$ is very well fitted by a Gaussian function (pale green solid line), which confirms the Gaussian properties of the map's noise. Above $\sim 4.5\sigma$ the distribution starts to diverge owing to the presence of positive sources. Our $5\sigma$ threshold to identify robust sources is indicated by the dashed orange line and arrow. All the pixels highlighted in orange, which have values below $-5\sigma$, lie within a region spanning approximately a beam size. This negative signal could be associated with a single noise peak at $-6\sigma$ or with a detection of the SZ effect, as discussed in Section 2.</figcaption></figure>

shown in Figure 1). This S/N map is then used to find source candidates by searching for pixels above an S/N threshold. Here we use the sources detected at $>5\sigma$ (where the false detection rate is expected to be $\lesssim 8\%$; see Section 2.1) to calculate the first ALMA number counts at 2 mm, as described below.

## 2.1. The 2 mm Number Counts

The cumulative number counts—i.e., the number of galaxies above a certain flux density threshold, S, per unit area—can be directly estimated by counting the number of detected galaxies
as a function of flux density and making the required corrections for contamination, completeness, and flux boosting.

Previous blind ALMA surveys have estimated a minimum, close-to-zero contamination fraction due to false detections (noise peaks) for S/N thresholds in the range of ≈4–5 (e.g., Fujimoto et al. 2016; Oteo et al. 2016). Nevertheless, despite adopting a conservative threshold of 5σ in our analysis (see Figure 3), a thorough characterization of the contamination rate is required given the large area covered by the MORA survey and the large number of independent beams in the maps.

The false detection rate of the MORA survey is characterized as follows: first, we create noise maps with the same dimensions as the original maps under the assumption that the noise in the ALMA observations is well represented by Gaussian statistics (as has been shown before for other ALMA data, e.g., Franco et al. 2018, and demonstrated in Figure 3). Second, the maps are convolved with a beam representative of the average synthesized beam in our observations. Then, to characterize the contamination rate as a function of S/N, the map is scaled such that the final standard deviation of the whole distribution of pixels is equal to 1. This is equivalent to an S/N map with no astronomical sources. Finally, we simply quantify the number of peaks above the adopted detection threshold using the same algorithm used to find the real sources. This process is repeated 100 times in order to measure the expected number of false sources and its corresponding uncertainty. At our adopted threshold of 5σ, the expected number of false detections in the MORA catalog is 0.8 ± 0.2. This is in good agreement with the findings in C. M. Casey et al. (2021, in preparation), who found near-IR counterparts for all the MORA detections with the exception of the lowest-significance source in the catalog.

To take into account this false detection rate in the number counts calculation, a statistical approach is adopted in which the contamination fraction is distributed among the 13 sources according to their S/Ns with a proportionality given by a normal probability distribution (i.e., estimating the probability that a normal random variable is greater than the sources’ S/Ns). As expected, this procedure gives a higher probability of false detection to the sources with the lowest S/Ns and an (almost) negligible probability to those detected at ≳ 5.5σ.

The survey completeness is calculated in a similar fashion in which artificial sources are inserted in the real maps (after removing the bright detections), followed by computing the ratio between the recovered sources and inserted sources. The artificial sources are inserted, one at a time, at random positions around the whole mosaic, and then they are recovered with the same source extraction procedure used to build the real source catalog. The survey completeness can be estimated as a function of input flux density, recovered flux density, detected S/N, local noise rms, or a combination of these parameters. A source is considered recovered if it is detected within a synthesized beam from the input position (~1′′5). In each flux density bin of 25 µJy, we repeat this process 10,000 times in order to sample the noise variation across the maps.

The same set of simulations are used to examine the flux-boosting effects (meaning sources’ flux densities systematically biased upward by noise and the presence of fainter undetected sources). This is done by estimating the average ratio between the sources’ intrinsic flux density and the recovered flux density in the simulations. Again, this value depends on the flux density and on the local noise of each source (or similarly, S/N), and

therefore several flux-boosting factors can be derived as a function of these parameters. In average, we find flux-boosting factors of around 5% for those sources detected just above our S/N threshold. These factors, which decrease with increasing S/N, are, however, much lower than the uncertainties on the measured flux densities (~20% for a 5σ detection).

The number counts can then be directly estimated by counting the number of detected galaxies as a function of flux density and making the appropriate corrections for contamination, completeness, and flux boosting. We estimate the contribution of a source with a deboosted flux density, Sᵢ, and measured S/N, σᵢ, to the cumulative number counts to be

$$ \eta_i(S_i, \sigma_i) = \frac{1 - f_{\text{cont}}(\sigma_i)}{\zeta(S_i)A_{\text{eff}}} \quad (1) $$

where fcont(σᵢ) is the estimated fraction of contamination at the measured S/N (σᵢ) of the source, ζ is the corresponding completeness for a deboosted flux density Sᵢ, and Aeff is the total area of the MORA maps used for source extraction (i.e., 184 arcmin²; see Figure 2). Finally, the cumulative number counts, N(>S), are estimated by the sum over all sources with a flux density higher than S, i.e., N(>S) = Σᵢηᵢ(Sᵢ, σᵢ).

In order to take into account the uncertainties associated with the correction factors (completeness and false detection) and flux densities in the estimation of the number counts, we perform a Monte Carlo simulation. In each realization, the adopted measured flux density for each source, Sᵢ, is extracted from a Gaussian distribution with a standard deviation equal to the measured error and centered at the observed value. Then, this new flux density is used to estimate a new S/N, σᵢ, by dividing by the local noise. These S/Ns are then used to update the contamination fraction and completeness accordingly, whose new values are drawn from Gaussian distributions. This procedure is repeated 100 times, with the mean values representing the final number counts and the 16th and 84th percentiles their associated confidence interval. Finally, given the relatively small number of sources in our catalog, Poisson uncertainties are added in quadrature according to Gehrels (1986). The MORA 2 mm number counts derived in this work are reported in Table 1 and are shown in Figure 4.

Given that the 2 mm flux density of one of our sources, MORA-10, was found to be contaminated by synchrotron emission (at a ~40% level; C. M. Casey et al. 2021, in preparation), we repeat all the processes described above after removing this source, aiming at providing the 2 mm number counts of star-forming galaxies only (note that the source would have fallen below our detection threshold without the synchrotron emission). This is particularly important since the model used in Section 4.1 does not take into account nonthermal emission; nevertheless, we highlight that the difference between the two estimations is not significant (see Table 1).

At this wavelength, the only determinations of the number counts reported in the literature beyond our MORA survey measurements come from the GISMO/IRAM surveys reported in Staguhn et al. (2014) and Magnelli et al. (2019), covering 31 and 250 arcmin², respectively. The latter are in very good agreement with our estimations, while the former lie above by a factor of ≈2.0–2.5. This overestimation is thought to be caused by the uncertainties in the flux deboosting factors in the confusion-limited map of Staguhn et al. (2014), as discussed in Magnelli et al. (2019), although cosmic variance might also be
<table><caption>**Table 1**
MORA Survey 2 mm Number Counts</caption><thead><tr><th colspan="4">Raw Number Counts</th><th colspan="4">Star-forming Galaxy Number Counts</th></tr><tr><th>S<sub>2 mm</sub><br>(mJy)</th><th>N(&gt; S<sub>ν</sub>)<br>(deg<sup>-2</sup>)</th><th>δN<sup>-</sup><br>(deg<sup>-2</sup>)</th><th>δN<sup>+</sup><br>(deg<sup>-2</sup>)</th><th>S<sub>2 mm</sub><br>(mJy)</th><th>N(&gt; S<sub>ν</sub>)<br>(deg<sup>-2</sup>)</th><th>δN<sup>-</sup><br>(deg<sup>-2</sup>)</th><th>δN<sup>+</sup><br>(deg<sup>-2</sup>)</th></tr></thead><tbody><tr><td>0.33</td><td>440</td><td>120</td><td>60</td><td>0.33</td><td>380</td><td>110</td><td>150</td></tr><tr><td>0.36</td><td>370</td><td>110</td><td>140</td><td>0.36</td><td>320</td><td>100</td><td>130</td></tr><tr><td>0.39</td><td>310</td><td>90</td><td>130</td><td>0.39</td><td>270</td><td>80</td><td>110</td></tr><tr><td>0.41</td><td>260</td><td>80</td><td>110</td><td>0.48</td><td>220</td><td>70</td><td>100</td></tr><tr><td>0.48</td><td>220</td><td>70</td><td>100</td><td>0.52</td><td>185</td><td>64</td><td>91</td></tr><tr><td>0.52</td><td>185</td><td>64</td><td>91</td><td>0.54</td><td>157</td><td>59</td><td>85</td></tr><tr><td>0.54</td><td>157</td><td>59</td><td>85</td><td>0.57</td><td>131</td><td>53</td><td>79</td></tr><tr><td>0.57</td><td>131</td><td>53</td><td>79</td><td>0.61</td><td>106</td><td>47</td><td>72</td></tr><tr><td>0.61</td><td>106</td><td>47</td><td>72</td><td>0.68</td><td>83</td><td>41</td><td>65</td></tr><tr><td>0.68</td><td>83</td><td>41</td><td>65</td><td>0.78</td><td>60</td><td>34</td><td>59</td></tr><tr><td>0.78</td><td>60</td><td>34</td><td>59</td><td>0.88</td><td>41</td><td>27</td><td>53</td></tr><tr><td>0.88</td><td>41</td><td>27</td><td>53</td><td>1.03</td><td>21</td><td>17</td><td>47</td></tr><tr><td>1.03</td><td>21</td><td>17</td><td>47</td><td></td><td></td><td></td><td></td></tr></tbody></table>

Note. The star-forming galaxy number counts are estimated after removing source MORA-10, whose flux density falls below our detection threshold after removing the contribution from synchrotron emission.

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

**Figure 4.** Cumulative number counts derived from our main sample (S/N > 5σ) are represented by the blue filled circles. Previous estimations of the number counts using the GISMO camera on the IRAM telescope are represented by the pale green and orange squares (Magnelli et al. 2019) and the yellow stars (Staguhn et al. 2014). Additionally, the predictions from the SHARK galaxy evolution model (Lagos et al. 2020) are illustrated by the gray line.

important given the relatively small mapped area. The predictions from the SHARK semianalytic model of galaxy formation (Lagos et al. 2020) are also plotted in the figure, which are indeed in very good agreement with our estimations.

### 3. Revising the 3 mm Number Counts

Zavala et al. (2018) compiled archival ALMA band 3 observations toward three extragalactic legacy fields: COS-MOS, CDF-S, and the UDS, resulting in a total of 135 individual maps adding up to a total area of ≈200 arcmin². After masking out the original targets of these observations and any galaxies that were potentially physically associated with them, they derived the first ALMA galaxy number counts at 3 mm using a total of 13 sources detected above 5σ.

Our recent ALMA 2 mm + 3 mm follow-up observations on this sample (ALMA projects: 2018.1.00478.S and 2019.1.00838. S; PI: J. Zavala) show that three of these sources (ALMA-3mm.14, ALMA-3mm.15, and ALMA-3mm.16) are not recovered, meaning that they are likely spurious detections. Two of these sources were indeed noted to have a spectral index that might be inconsistent with thermal emission in Zavala et al. (2018); nevertheless, without further data at the time, the low S/Ns (5.0–5.2) prevented a firm conclusion, and therefore they were included in the previous number counts estimation.

Thus, we have revised the false detection rate to be higher than the value reported in the original work of Zavala et al. (2018), for which only one source was expected to be false. This highlights the complexity of the interferometric data, particularly when using observations with different beam sizes, integration times, and array configurations, as is typical for data sets derived from archival projects.

Here we revised the 3 mm number counts following exactly the same procedure as described in Zavala et al. (2018) but removing the three spurious sources. Additionally, we updated the flux density of the source ALMA-3mm.03 (aka ASPECS-3mm.1) since the reported value in Zavala et al. (2018) seems to be contaminated by a bright CO emission line (González-López et al. 2019); we thus adopt the flux density reported in the ASPECS catalog.

The revised 3 mm number counts are reported in Table 2 and shown in Figure 5. The new values are significantly different at the bright end (S<sub>2 mm</sub> > 0.2 mJy) since two of the false detections were the brightest galaxies in the sample.<fnref n="27" /> At fainter flux densities, the updated number counts are a factor of ~1.5× lower and, thus, in better agreement with the ASPECS results (González-López et al. 2019), although still a factor of ~2× higher. Nevertheless, the two number counts are consistent with each other within the uncertainties. The difference is thus not statistically significant.

<footnote n="27">²⁷ Given the inhomogeneous observations in the archival data, a low S/N does not necessarily imply a low flux density since each map has a different noise rms depth. The two brightest false detections were indeed detected in areas with large rms values.</footnote>
<table><caption>**Table 2**
Revised ALMA Archival 3 mm Number Counts</caption><thead><tr><th>S<sub>3</sub> mm<br>(μJy)</th><th>N(> S<sub>3</sub>)<br>(deg<sup>-2</sup>)</th><th>δN<sup>-</sup><br>(deg<sup>-2</sup>)</th><th>δN<sup>+</sup><br>(deg<sup>-2</sup>)</th></tr></thead><tbody><tr><td>33</td><td>4440</td><td>1200</td><td>1990</td></tr><tr><td>52</td><td>2140</td><td>640</td><td>1000</td></tr><tr><td>63</td><td>1420</td><td>470</td><td>710</td></tr><tr><td>79</td><td>990</td><td>360</td><td>540</td></tr><tr><td>98</td><td>750</td><td>290</td><td>450</td></tr><tr><td>109</td><td>550</td><td>230</td><td>370</td></tr><tr><td>117</td><td>410</td><td>190</td><td>320</td></tr><tr><td>127</td><td>290</td><td>150</td><td>280</td></tr><tr><td>137</td><td>170</td><td>110</td><td>230</td></tr><tr><td>167</td><td>74</td><td>55</td><td>170</td></tr></tbody></table>

## 4. Constraining the IR Luminosity Function and Dust-obscured Star Formation Rate Density

In this section we use a backward evolution model in combination with the state-of-the-art FIR/millimeter surveys to constrain the IR luminosity function. The model draws spectral energy distributions (SEDs) from the known breadth and characteristics of dusty galaxies and then works backward to discriminate between different luminosity function scenarios using the galaxy number counts as constraints. Finally, the dust-obscured SFRD is estimated by integrating the best-fit IR luminosity function (see Zavala et al. 2018 for a similar analysis).

### 4.1. Model Description

The adopted backward evolution model is described in detail in Casey et al. (2018a, 2018b), where the specifics of the model and the assumed values for each parameter can be found. A summary of the salient characteristics and assumptions follows.

The model combines a parameterized evolving galaxy IR luminosity function with the thermal SED properties of galaxies’ dust emission to make predictions for galaxy (sub) millimeter surveys.

Galaxies’ SEDs are described by a modified blackbody function of the form $S_ν \propto (1 - e^{-(ν/v_0)^{\beta_{em}}})B_ν(T_{dust})$, with an extra mid-IR power law to account for hotter dust emission at $T \gg T_{dust}$ (e.g., Casey 2012). The SED library follows the well-known luminosity-temperature relationship in terms of $L_{IR} - \lambda_{peak}$ and takes into account the observed dispersion in this relation. We highlight that using the peak wavelength rather than dust temperature minimizes the impact of the assumed effective dust opacity in the model ($\tau = 1$ at $\lambda_{rest} = 100$ µm), which could change the derived dust temperature but has a minor impact on the measured $\lambda_{peak}$. Note that although this relation is not assumed to evolve with redshift (which is supported by the current data at least up to $z \sim 4.5$; Casey et al. 2018b; Dudzevičiūtė et al. 2020), if one selects galaxies on the main sequence, a redshift evolution of $\lambda_{peak}$ (proportionally to dust temperature) naturally arises. This evolution, which is driven by the underlying evolution of the main sequence of star-forming galaxies, is in line with recent results from the literature (e.g., Bethermin et al. 2015; Schreiber et al. 2015; Magdis et al. 2017).

As is common practice in the literature (e.g., Sanders et al. 2003; Magnelli et al. 2011, 2013; Lim et al. 2020), the assumed IR luminosity function in the model is described by a double

power law of the form

$$ \Phi(L, z) = \begin{cases} \Phi_{*}(z) \left( \frac{L}{L_{*}(z)} \right)^{\alpha_{LF}}, & \text{if } L < L_{*}, \\ \Phi_{*}(z) \left( \frac{L}{L_{*}(z)} \right)^{\beta_{LF}}, & \text{if } L \ge L_{*}, \end{cases} \quad (2) $$

where $\alpha_{LF}$ and $\beta_{LF}$, represent the slopes at faint and bright luminosities, respectively, while $L_*$ and $\Phi_*$ are the characteristic galaxy luminosity and characteristic number density, two fundamental parameters that are allowed to evolve with redshift.

Following previous works from the literature, the characteristic number density is assumed to evolve as

$$ \Phi_{*} \propto \begin{cases} (1+z)^{\Psi_{1}}, & \text{if } z < z_{\text{turn}}, \\ (1+z)^{\Psi_{2}}, & \text{if } z \ge z_{\text{turn}}. \end{cases} \quad (3) $$

Significant observational efforts have provided good constraints on the evolution of $\Phi_*$ at low redshifts, revealing a flat trend with $\Psi_1$ close to zero and a turnover redshift of $z_{turn} \approx 2$, in line with the peak of the cosmic SFRD (CSFRD; Magnelli et al. 2011; Gruppioni et al. 2013; Lim et al. 2020). Although it is clear that at higher redshifts the evolution is significantly steeper, the current estimations are highly uncertain and limited to $z \lesssim 4$, with values ranging from $\Psi_2 \sim -6$ (or lower; Koprowski et al. 2017) to $-2$ (or even higher; Rowan-Robinson et al. 2016), with a few measurements in between. Indeed, Zavala et al. (2018) found an intermediate value of $\Psi_2 \approx -4.2$ using the initial 3 mm sample.

In this work, we aim to constrain two of the most important unknowns of the galaxy luminosity function: the faint-end slope, $\alpha_{LF}$, and the redshift evolution of the characteristic number density at high redshift driven by $\Psi_2$ (since $\Phi_* \propto (1+z)^{\Psi_2}$ for $z \gtrsim 2$). We also explore different values of the dust emissivity index, $\beta_{em}$, in order to minimize any possible bias and to test possible degeneracies between the different parameters.

For the rest of the model parameters, we adopt the same values as in Casey et al. (2018b) and Zavala et al. (2018). We refer the reader to Casey et al. (2018b, and particularly to its Appendix A.1) for a detail discussion about these choices and their impact on the model. In short, the bright-end slope of the luminosity function is fixed to be $\beta_{LF} = -3.0$. This value is relatively well constrained (at least up to $z \sim 4$) with a measured dispersion of $1\sigma \approx 0.15$. The evolution of the characteristic luminosity, $L_*$, follows $L_* \propto (1+z)^{\gamma_1}$ for $z < z_{turn}$ and $L_* \propto (1+z)^{\gamma_2}$ for $z \ge z_{turn}$. While the evolution at lower redshifts ($z < z_{turn}$) is well characterized (we adopt $\gamma_1 = 2.8$), the evolution of $L_*$ beyond $z_{turn}$ is more uncertain. The model adopts $\gamma_2 = 1.0$, implying that $L_*$ continues to evolve upward toward higher redshifts, in line with expectations for hierarchical structure formation and cosmic downsizing. Additionally, adopting $\gamma_2 = 1.0$ makes the $L_*$ evolution consistent with the $L_*$ of the quasar luminosity function (Hopkins et al. 2007). Adopting a reversal evolution ($\gamma_2 \le 0$) underpredicts the number counts at wavelengths longer than $\sim 850$ µm, regardless of the adopted number density (for any $\Psi_2 < -1.5$), while a more rapid evolution with $\gamma_2 \ge 1.5$ would imply very bright values for $L_*$ at $z \sim 5$, which might even exceed the luminosity of the brightest DSGFs known to date. Moreover, in order to match the contribution from LIRGs and
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Cumulative number counts—galaxies above a given flux density per unit area—at 1.2 mm (left panel), 2 mm (middle panel), and 3 mm (right panel). The data points used for the fitting analysis are illustrated as blue filled circles, while other measurements from the literature are plotted as green squares. The best-fit number counts from the model, represented by the gray lines, nicely reproduce the number counts at the three different wavelengths simultaneously, as well as number counts at shorter wavelengths (see Appendix A), spanning more than a decade in wavelength and orders of magnitude in flux density. The best-fit number counts from the model are also broken down into two ranges of luminosity and three redshift bins. Galaxies with $L_{IR} > 10^{12} L_{⊙}$ are shown by solid lines, while those with $L_{IR} < 10^{12} L_{⊙}$ are shown by dotted lines. Similarly, the three different redshift ranges are illustrated by different colors: blue for $z < 3$, orange for $3 < z < 6$, and gold for $z > 6$.</figcaption></figure>

ULIRGs to the cosmic SFRD at $z < 2.5$ to those reported in the literature, a modest dependence of $z_{\text{turn}}$ on $\Psi_2$ of the form $z_{\text{turn}} = 1.6-0.09\Psi_2$ has to be adopted (see Casey et al. 2018b for details).

This model has some caveats that the reader should keep in mind. The exact values of $L_0$ and $\Phi_0$ are correlated with $\gamma_1$ and $\Psi_1$, respectively.<fnref n="28" /> The adopted values ($L_0 = 10^{11.1} L_\odot$, $\Phi_o = 10^{-3.5} \text{Mpc}^{-3} \text{dex}^{-1}$, $\gamma_1 = 2.8$, and $\Psi_1 = 0$) were chosen to simultaneously reproduce (i) the IRLF at $z \lesssim 2$, (ii) the reported values of $L_*$ and $\Phi_*$ from the literature, and (iii) the relative contribution from LIRGs and ULIRGs to the total SFRD within $0 < z < 2$. This combination of values is, however, not necessarily unique. Another caveat that the reader should keep in mind is the assumption of a fixed emissivity spectral index at all redshifts, and although several values were explored in the analysis (from $\beta = 1.6$ to 2.6; see Figure 6), each realization adopts a single value across all redshifts for all sources. This is not necessarily in disagreement with previous literature, given the current limited data. Although a spread on $\beta$ has been found across different galaxies, these variations have not yet been shown to correlate with redshift or other galaxy characteristics that may impact the model. Indeed, there are clear examples of $\beta \sim 1.5$ to $\sim 2.0$ in both local galaxies (e.g., Remy-Ruyer et al. 2013) and high-redshift systems (e.g., Jin et al. 2019). The same is true for the adopted slopes of the luminosity function at both the faint and bright ends ($\alpha_{LF}$ and $\beta_{LF}$, respectively). Although a redshift evolution is possible given that we know that similar parameters do evolve for the UV luminosity function (e.g., Finkelstein et al. 2015), the lack of samples of galaxies with spectroscopic redshifts covering a wide range of IR luminosities prevents its confirmation (or rejection). Indeed, nonevolving

slopes have been adopted in several works in the literature (e.g., Magnelli et al. 2011, 2013; Gruppioni et al. 2013; Lim et al. 2020). Therefore, in the absence of more information, we choose to fix these quantities. The impact of these assumptions on our results is further discussed below.

## 4.2. Fitting Methodology and Data Constraints

To find the best-fit model parameters, first, the different evolutionary scenarios of the IRLF are combined with galaxies' SEDs to create mock observations that resemble the real surveys in terms of wavelength, noise depth, and angular resolution, with the CMB heating effects taken into account (following da Cunha et al. 2013). Second, the respective number counts are estimated in a similar fashion to that in the real observational works. These simulated number counts are then compared to the measured number counts in a joint analysis that combines multiple observations at different wavelengths simultaneously. Finally, the best-fit model parameters that better reproduce all the observed number counts are derived.

The cumulative number counts are preferred over the differential number counts since the latter suffer from larger uncertainties due to the relatively small number of sources in each flux density bin. However, we test that adopting the differential number counts instead does not change significantly our results.

As discussed in Casey et al. (2018a, 2018b), observations at long wavelengths ($\lambda \gtrsim 1$ mm) are ideal to distinguish between the different scenarios of the model since the majority of constraining data sets at shorter wavelengths only inform about the evolution of the IRLF at $z \lesssim 3.0$, where the model parameters are already relatively well understood. Therefore, here we only use number counts at 1.2, 2, and 3 mm as constraints, which also reduces significantly the computational cost of the fitting analysis, which would be prohibitive

<footnote n="28">²⁸ We refer the reader to Equations (8) and (9) from Casey et al. (2018b) for the exact parameterization of the evolution of $L_*$ and $\Phi_*$ and the role of $L_0$ and $\Phi_0$.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Constraints on the IRLF model parameters $\alpha_{LF}$ (the faint-end slope of the luminosity function), $\Psi_2$ (which governs the evolution of $\Phi_*$ at high redshift given $\Phi_* \propto (1+z)^{\Psi_2}$), and $\beta_{em}$ (the dust emissivity index). The confidence regions at 68%, 95%, and 99.7% for each parameter derived using the maximum likelihood approach are represented by the blue contours (from dark blue to light blue). Additionally, the best-fit values extracted from 100 realizations of the Nelder-Mead multidimensional minimization algorithm (amoeba) are illustrated with the green filled circles. On the top of each column, we show the respective 1D marginalized probability distribution derived from the maximum likelihood approach (solid blue line) and the histogram of best-fit values found by the Nelder-Mead method (green dashed histograms). The best-fit values derived from both methods and their corresponding $1\sigma$ uncertainties (68% C.I.) are shown on the top of each row. Consistent results are derived from the two different approaches.</figcaption></figure>

otherwise. We adopt the MORA 2 mm number counts described above (see Section 2.1), the updated 3 mm number counts reported in this work (see Section 3), and the 1.2 and 3 mm number counts from the ASPECS survey (González-López et al. 2019, 2020). Altogether, these surveys add a total of ~400 arcmin² surveyed area and around 60 detected galaxies. Note that active galactic nucleus contribution or nonthermal emission (e.g., synchrotron emission) is expected to be negligible at these wave bands and at the flux density ranges explored here<fnref n="29" /> (see discussion in Casey et al. 2018a; Zavala et al. 2018).

Once the fiducial number counts have been defined, to derive the best-fit model parameters that better reproduce these number counts, we use two different methods: a maximum likelihood approach and a multidimensional minimization algorithm.

The maximum likelihood approach is done in a similar fashion to that reported by Zavala et al. (2018). We create a three-dimensional grid, representing the parameter space formed by $\alpha_{LF}$, $\Psi_2$, and $\beta_{em}$, with values of $\alpha_{LF} = [-0.6, -0.2]$ in steps of 0.05, $\Psi_2 = [-8.8, -2.2]$ in steps of 0.4, and $\beta_{em} = [1.5, 2.6]$ in steps of 0.05. For each combination of the parameter's subset (corresponding to different points in the grid), the likelihood function of the measured cumulative number counts relative to the simulated number counts is estimated. Finally, the best-fit model is assumed to be the one with the highest probability, and the confidence intervals at 68%, 95%, and 99.7% are obtained by integrating the normalized likelihood distribution.

The Nelder-Mead optimization method (aka the Amoeba method; Nelder & Mead 1965; Press et al. 1992) relies on a downhill simplex algorithm to perform a multidimensional minimization of a given parameter, in this case, the square differences ($\chi^2$) between the real and the modeled cumulative number counts. In each realization, we randomly vary the

<footnote n="29">²⁹ From all the sources used in this work, we found a nonnegligible contamination from synchrotron emission only in one source, namely, MORA-10. This source has been removed from our analysis (see Section 2.1).</footnote>
initial starting point of the search, with each search limited to 100 evaluations. During this analysis, we introduce a fourth parameter, $z_{cutoff}$, following Zavala et al. (2018), which represents the redshift above which no more dusty galaxies exist (with a range of $z_{cutoff} = 5.5–9.0$). Nevertheless, we find that $z_{cutoff}$ is not well constrained by the data and has a minor impact on the analysis.

As shown in Figure 5, the best-fit number counts recovered from the model are in good agreement with the measured values, reproducing simultaneously not only the data used in the fitting but also the rest of the measurements reported in the literature, as discussed below.

At 1.2 mm, beyond the ASPECS number counts, the model predictions nicely reproduce the number counts from the 1.6 deg$^2$ AzTEC surveys reported in Scott et al. (2012) and those from the ALMA follow-up survey of SCUBA-2 galaxies detected over ~1 deg$^2$ reported by Stach et al. (2018),<fnref n="30" /> which probe the brightest flux densities. At fainter flux densities, the model also reproduces the more recent ALMA results from Dunlop et al. (2017) and Umehata et al. (2017). The number counts from the GOODS-ALMA survey (Franco et al. 2018; covering the flux density range of $S_{1.2\,mm} \approx 0.5–2\,mJy$) are, however, lower than our estimations, and those reported by Fujimoto et al. (2016) (which probe the faintest flux densities in the figure) lie above our estimations. Note, however, that the measurements reported by Fujimoto et al. (2016) come from gravitationally lensed fields, and therefore these discrepancies might be caused by the uncertainties in the magnification factors and small survey area (see discussion in González-López et al. 2020). The discrepancy with the Franco et al. (2018) results is likely related to the relatively low completeness associated with their high angular resolution observations ($\theta_{FWHM} = 0.2''$29 and $\theta_{FWHM} = 0.2''$60 after filtering the data).

At 2 mm, the only determinations of the number counts reported in the literature beyond our MORA survey measurements come from the GISMO/IRAM surveys reported in Staguhn et al. (2014) and Magnelli et al. (2019), covering 31 and 250 arcmin$^2$, respectively. The latter are in very good agreement with the model predictions, while the former lie above by a factor of ≈2.0–2.5. This difference is thought to be caused by the uncertainties in the flux deboosting factors in the confusion-limited map of Staguhn et al. (2014), as discussed in Magnelli et al. (2019), plus the possible effects of cosmic variance (see Section 2.1).

At 3 mm, the only measurement plotted in Figure 5 that was not used in our analysis is the brightest bin of the ASPECS survey, which represents only an upper limit. This value is, however, in good agreement with the model predictions.

The success of the model can also be illustrated by its predicting power at other wavelengths. The model does reproduce the number counts at shorter wavelengths, from 70 to 850 µm spanning over two decades of observations over hundreds of square degrees (see Appendix A).

### 4.3. The IRLF and Best-fit Parameters

As is clearly shown in Figure 6, both methods provide consistent results and good constraints on the model parameters

governing the IRLF while successfully reproducing the observed number counts (see Figures 5 and 9).

The faint-end slope of the IR luminosity function is found to be flat, with $\alpha_{LF} = -0.42^{+0.02}_{-0.04}$. This is in line with recent studies based on the deepest single-dish observations and interferometric surveys (although those were limited to lower redshifts), with reported values of $\alpha_{LF} = -0.4$ (Koprowski et al. 2017) and $-0.5 \pm 0.7$ (Lim et al. 2020), and with the observed flattening in the 1.1 mm luminosity function (Popping et al. 2020).

Although the faint-end slope is fixed for all redshifts in our model, we highlight that an evolution to steeper values with increasing redshift is inconsistent with our data and with the very low number of $z>4$ galaxies detected in the deepest ALMA surveys (Dunlop et al. 2017; Hatsukade et al. 2018; Aravena et al. 2020). Indeed, most of the sources in the ASPECS survey lie at $z<3$ (Aravena et al. 2020), in good agreement with our model predictions (see Figure 5). This IR flat slope contrasts with the steep faint-end slope of the UV luminosity function at high redshift, which ranges from -1.6 at $z \sim 4$ to -2.0 at $z \sim 7$ (e.g., McLure et al. 2013; Finkelstein et al. 2015). As a consequence, deep pencil-beam observations at (sub)millimeter wavelengths would not significantly increase the number of detected sources, as is commonly the case in the rest-frame UV/optical observations (e.g., Ferguson et al. 2000).

Regarding the evolution of the characteristic number density at high redshift, the best-fit $\Psi_2$ value of $-6.5^{+0.8}_{-1.8}$ implies a steep redshift evolution of the IR number density ($\Phi_* \propto (1+z)^{-6.5^{+0.8}_{-1.8}}$), which disfavors the dust-rich hypothetical model discussed in Casey et al. (2018b), for which a value of $\Phi_* \propto (1+z)^{-2.5}$ was adopted. This indicates that DSFGs are indeed rare at early epochs. This steep evolution is also in line with the rapid drop-off of the quasar luminosity function at high redshift (Hopkins et al. 2007) and similar to the number density evolution of UV-bright ($M_{UV} = -21$) galaxies ($\propto (1+z)^{-5.9}$; Finkelstein et al. 2015), which might suggest that these galaxies occupy similar dark matter halos.

There is, though, a caveat related to the effective dust optical opacity of our SEDs that the reader should keep in mind. Although, as mentioned before, the SEDs are parameterized in terms of $\lambda_{peak}$ rather than dust temperature, incorporating the CMB effects on the heating and detectability of the sources requires an estimation of the dust temperature (da Cunha et al. 2013). Assuming an optically thin opacity form would result in lower dust temperatures for the same $\lambda_{peak}$ than those derived from the optically thick model, decreasing the contrast between the dust emission and the CMB. This implies that, in order to match the measured number counts, a higher number of galaxies would be required compared to the optically thick model, i.e., a higher $\Psi_2$ value (see discussion in Zavala et al. 2018). Nevertheless, this effect is only important if the number counts are dominated by very high redshift ($z > 6$) sources, which our results suggest is likely not the case.

Finally, Figure 6 shows that the dust emissivity spectral index is well constrained to be $\beta_{em} = 1.8 \pm 0.1$, in very good agreement with the values reported in the literature (e.g., Dunne et al. 2011; Galliano et al. 2018; although we note that measurements of this parameter at high redshifts are scarce). This figure also reveals, unsurprisingly, a mild correlation between $\beta_{em}$ and $\Psi_2$ since higher values of $\beta_{em}$ imply lower

<footnote n="30">³⁰ The results from AzTEC observations and other ALMA surveys at 1.1 mm were scaled by a factor of 0.8, while those from the ALMA/SCUBA-2 survey at 850 µm were scaled by a factor of 0.4.</footnote>
flux densities at long wavelengths (for a given IR luminosity, redshift, and dust temperature).

We highlight that despite the large range of values explored for the three different parameters and despite the caveats described above, α$_{LF}$ and β$_{em}$ are well constrained, and while the uncertainties on Ψ$_2$ seem large, our constraints rule out values grater than Ψ$_2$ ~ 4, which would imply a significantly larger number of DSFGs at high redshifts.

### 4.4. The History of Dust-obscured Star Formation

The SFRD as a function of redshift can now be calculated by integrating the best-fit IR galaxy luminosity function$^{31}$ and propagating the associated uncertainties.

An important source of uncertainty is the field-to-field variation due to the large-scale structure of the universe, which is known as cosmic variance. To infer its impact in our estimations, we adopt the results from a model for the dust continuum number counts of galaxies (Popping et al. 2020) that builds on the UNIVERSEMACHINE model (Behroozi et al. 2019), from which 1000 different light cones of 400 arcmin$^2$ (approximately the combined area the MORA, ASPECS, and the 3 mm archival surveys) are used to measure the variance in the number of detected sources. The cosmic variance is then defined following Moster et al. (2011):

$$ \sigma^2 \equiv \frac{\langle N^2 \rangle - \langle N \rangle^2}{\langle N \rangle^2} - \frac{1}{\langle N \rangle}, \quad (4) $$

where $\langle N^2 \rangle$ and $\langle N \rangle$ are the variance and mean number of sources in the light cones, respectively. At the most critical redshifts of our work ($z > 3$), the cosmic variance is estimated to be ≈35%. At $z \approx 7.0$, the cosmic variance increases significantly and our data sets suffer a loss of constraining power. Furthermore, an increasing fraction of the star formation activity is expected to be in highly clustered structures (i.e., galaxy protoclusters; Chiang et al. 2017) that might have been missed in our surveys. Therefore, the model predictions presented in this work are limited to $z \lesssim 7.0$.

The inferred dust-obscured star formation history and its associated uncertainty, including that from cosmic variance, are presented in Figure 7. As can be seen, it is dominated by galaxies with $L_{IR} \approx 10^{11}-10^{13}$, a population that our surveys are particularly sensitive to (see Figure 10 in Appendix B). Figure 7 also includes other determinations of the CSFRD from the literature, which are split into two groups: FIR/submillimeter and UV/optical-based measurements. The first set includes the works by Magnelli et al. (2011, 2013, 2019), Gruppioni et al. (2013), Casey et al. (2012), Swinbank et al. (2014), Bourne et al. (2017), Koprowski et al. (2017), Liu et al. (2018), Williams et al. (2019), Wang et al. (2019b), Lim et al. (2020), and Dudzevičiūtė et al. (2020), and the second set the results reported in Wyder et al. (2005), Schiminovich et al. (2005), Robotham & Driver (2011), Cucciati et al. (2012), Dahlen et al. (2007), Reddy & Steidel (2009), Bouwens et al. (2012), Schenker et al. (2013), and Finkelstein et al. (2015). We remind the reader that our model is not a fit to those data points but is shown as a mean for comparison.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Inferred dust-obscured star formation history—illustrated by the orange shaded region in the bottom panel. For comparison, we plot independent measurements from the literature based on IR/submillimeter and UV surveys (orange circles and blue squares, respectively) and the average unobscured star formation derived from rest-frame UV optical surveys (i.e., not corrected for dust attenuation; blue shaded region; Finkelstein et al. 2015). The total inferred SFRD derived in this work is shown in gray. The uncertainties in our estimation include those from the best-fit parameters and cosmic variance. The middle panel represents the fraction of obscured star formation, SF$_{obs}$/ (SF$_{obs}$ + SF$_{unobs}$), and its associated uncertainty (lighter shaded area). The contribution of dust-obscured galaxies, which dominates the cosmic star formation history through the past ~12 Gyr, rapidly decreases beyond its maximum, reaching values that are comparable to the unobscured star formation traced by the rest-frame UV/optical surveys by $z \approx 4-5$. The top panel represents the contribution from galaxies with different luminosity ranges to the dust-obscured SFRD, being dominated by ULIRGs (ultraluminous IR galaxies; $10^{12} L_{\odot} < L_{IR} < 10^{13} L_{\odot}$) and LIRGs ($10^{11} L_{\odot} < L_{IR} < 10^{12} L_{\odot}$).</figcaption></figure>

The dust-obscured component traced by the FIR to millimeter surveys has dominated the cosmic history of star formation for the past ~12 billion years, with a peak era in the range of $z = 2-2.5$ (~10–11 Gyr ago) and contributing around ~80% of total SFRD (see middle panel in Figure 7). Beyond this peak redshift, the dust-obscured star formation rapidly decreases, following the strong evolution of the number density of the IR luminosity function (see Section 4.3), with values that are comparable to the unobscured star formation traced by the rest-frame UV/optical surveys at $z \sim 4$.

At higher redshifts, the dust reprocessed SFRD becomes less dominant than the unobscured star formation. This is because of the combination of the flat faint-end slope of the IR luminosity function and the steep redshift evolution of its number density as compared to the UV luminosity function. At $z = 5$, the dust-obscured star formation represents $35\%^{+10\%}_{-25\%}$ of the total SFRD and decreases to $25\%^{+15\%}_{-20\%}$ at $z = 6$. Given that massive, IR-bright galaxies dominate the obscured component (as shown in the top panel of Figure 7), the observed decline of the dust-obscured star formation likely reflects the dearth of massive galaxies at high redshifts.

<footnote>$^{31}$ The IR luminosity function is integrated over the interval log ($L/L_{\odot}$): [9, 13.8]. A change of these limits has a minor impact on our results since the majority of the contribution arises from galaxies with luminosities in the range of $L_{IR} = 10^{11}-10^{13} L_{\odot}$ (see Figure 7).</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Total (obscured + unobscured) SFRD derived in this work (gray region) in comparison to the predictions from galaxy evolution models. The solid blue line represents the predictions from the cosmological hydrodynamical IllustrisTNG simulations (Pillepich et al. 2018), the dashed champagne line those from the SHARK semianalytic model (Lagos et al. 2018), and the dashed-dotted green line the results from the SIDES simulations (Béthermin et al. 2017). Additionally, we include the density evolution of luminous quasars from Wang et al. (2019a) (scaled for better visualization), which strongly resemble the shape of the SFRD from bright DSFGs with $L_{IR} > 10^{12} L_{\odot}$ illustrated by the orange dashed region.</figcaption></figure>

### 4.4.1. Comparison to Other Measurements and Model Predictions

The SFRD described above is in line with previous results from (sub)millimeter surveys, although most of them were limited to $z \lesssim 4-5$. For example, Dunlop et al. (2017) reported a transition from unobscured-dominated to obscured-dominated star formation at $z \approx 4$, in very good agreement with our results. Nevertheless, the small area of their survey ($\sim 4.5$ arcmin$^2$) prevented them from deriving conclusions beyond this redshift. The results from the larger ALMA surveys presented by Hatsukae et al. (2018) and Franco et al. (2018), covering 26 and 69 arcmin$^2$, respectively, also indicate a minor contribution from DSFGs in the $z \approx 4-5$ range (see also Yamaguchi et al. 2019). More recently, Dudzevičiūtė et al. (2020) used ALMA observations to investigate the properties of $\sim 700$ DSFGs detected over $\sim 1$ deg$^2$ and inferred an SFRD in the range of $\rho_{\text{obs}} \approx 3 \times 10^{-3}$ at $z \approx 4-5$ for galaxies with $S_{850 \ \mu m} > 1$ mJy (see also Koprowski et al. 2017), in very good agreement with our results (particularly if we look at the contribution from $L_{IR} > 10^{12} L_{\odot}$ galaxies; see Figure 8).

The “wedding cake” structure of all these surveys ensures that the contributions from the more abundant faint galaxies and the rare bright sources are accounted for (with the possible exception of the most extreme galaxies with $\log(L_{IR}/L_{\odot}) \gtrsim 13.5$; see Appendix B), suggesting that we are not missing any significant population of galaxies. This is also supported by the fact that the number counts predicted by the model are in good agreement with the deepest ALMA observations achieved to date and with the large-area single-dish telescope surveys (see Figures 5 and 9).

Other studies based on stacking analysis on (sub)millimeter maps and using samples of UV/optically selected galaxies have also concluded that the fraction of star formation that is obscured by dust decreases at high redshifts (e.g., Capak et al. 2015; Bouwens et al. 2016; Fudamoto et al. 2020). Indeed, Bouwens et al. (2020) complemented previous results from the literature with their dust-corrected SFRs to estimate the obscured and unobscured components of the cosmic history of star formation and found that the CSFRD transitions from

being primarily unobscured to obscured at $z \sim 5$, in relatively good agreement with our estimations$^{32}$ from direct millimeter-selected samples.

Our results are also compatible with estimations from other tracers such as radio observations. For example, using the VLA-COSMOS 3 GHz radio survey, Novak et al. (2017) inferred lower limits for the SFRD up to $z \sim 5$ by integrating the radio luminosity function after converting the radio luminosities to SFRs. These VLA-COSMOS 3 GHz radio constraints are consistent with our measurements. The authors also provided a completeness-corrected estimation of this quantity by extrapolating the luminosity function to account for the faintest star-forming galaxies. Their estimates show good consistency with our total SFRD, with the possible exception of their last bin at $z \sim 5$. Nevertheless, the large extrapolations involved in this process plus the systematic uncertainties that go into calculating SFRD from radio data introduce very large uncertainties in these measurements. The most recent estimates of SFRD from the VLA-COSMOS 3 GHz survey presented by Leslie et al. (2020) are in better agreement with our results.

There are, however, a few other studies that have proposed a different picture, with the obscured component significantly dominating the CSFRD back to $z \sim 5$. Rowan-Robinson et al. (2016) estimated the obscured SFRD using a sample of Herschel-selected galaxies over $\sim 20$ deg$^2$, finding values far higher than the UV estimations. Their constraints, however, cannot rule out the possibility that the UV-based measurements are dominant given the large uncertainties on their estimations, as mentioned by the authors. Additionally, it is possible that their calculations might be contaminated by the effect of gravitationally lensing or by overestimated Herschel flux densities (in view of their extreme SFRs that extend to $20,000 M_{\odot} \text{yr}^{-1}$).

More recently, Gruppioni et al. (2020) derive the dust-obscured SFRD using the serendipitously detected sources in the ALPINE survey, finding values in excess of those derived from UV/optical surveys even at $z \sim 5$. While the use of highly confused Herschel observations might also overestimate the derived luminosities and SFRs, we think that the discrepancy is mainly due to the possible clustering of serendipitous detections around the original targets. The clustering is expected since the observations’ original targets are massive galaxies ($\log(M/M_{\star}) \gtrsim 10.5$) at $z \approx 4-6$. In this case, their measurements would be more representative of an overdense region in the large-scale structure of the universe. We note, however, that a significant fraction of their estimates comes from the extrapolation of the IR luminosity function since the reported total SFRD is a factor of $\sim 5$ greater than the SFRD estimated when using only the detected sources. Therefore, it is also possible that the assumptions on the extrapolation of the IRLF are responsible for part of the observed discrepancy.

In Figure 8 we compare the SFRD derived in this work with results from galaxy evolution models, including the predictions from the cosmological hydrodynamical IllustrisTNG simulations (Pillepich et al. 2018), the SHARK semianalytic model

<footnote>$^{32}$ Note, however, that in Bouwens et al. (2020) most of the dust-obscured star formation at $z > 4$ is produced by faint (UV-selected) galaxies with $L_{IR} < 10^{12} L_{\odot}$. Nevertheless, those estimations were done assuming SED templates whose dust temperatures increase with redshift. If there is no significant evolution on the dust temperature (e.g., Dudzevičiūtė et al. 2020), then the SFRs derived for these galaxies would be lower by a factor of $\sim 2.5$, as discussed by the authors. This would decrease the contribution of faint galaxies and bring our results into better agreement.</footnote>
(Lagos et al. 2018), and the results from the SIDES simulations (Béthermin et al. 2017). Generally speaking, and taking into account the associated uncertainties in these values, our results are in good agreement with the aforementioned studies, pointing toward a convergence between galaxy evolution models and observations.

Finally, we include the density evolution of luminous ($M_{1450} < -26$) quasars (Wang et al. 2019a) in Figure 8 (scaled for visualization). The shape of the space density of bright quasars strongly resembles that from the bright DSFGs with $L_{IR} > 10^{12} L_{\odot}$ (shaded orange region in the figure), suggesting a connection between these two populations and therefore between the onset of star formation and the growth of their massive black holes (e.g., Wall et al. 2005).

## 5. Discussion and Conclusions

Exploiting the FIR and (sub)millimeter data aggregated over the past two decades, and particularly the state-of-the-art ALMA blind surveys at 1.2, 2, and 3 mm, we have constrained the evolution of the IR luminosity function and dust-obscured star formation in the past 13 billion years, back to $z \sim 7$. This is achieved by combining a model of the DSFG population with those long-wavelength observations, inferring constraints on the prevalence and characteristics of these galaxies through measurements of galaxy number counts at different wavelengths.

By using a library of SEDs and assuming different evolutionary scenarios for the IR luminosity function, the model makes predictions for galaxy (sub)millimeter surveys and then works backward to discriminate between the different luminosity function scenarios through the galaxy number counts. Finally, the dust-obscured SFRD is estimated by integrating the IR luminosity function (see also Casey et al. 2018a, 2018b; Zavala et al. 2018).

The model's constraints used in this work include our 2 mm galaxy number counts derived as part of the 2 mm MORA survey, the largest ALMA survey to date sensitive enough to detect DSFGs and dust-obscured star formation up to the epoch of reionization, and the only one carried out at a 2 mm wavelength (see also C. M. Casey et al. 2021, in preparation and S. Manning et al. 2021, in preparation). Additionally, we use the 3 mm number counts reported in Section 3, which are based on our ALMA follow-up observations on the sample of 3 mm selected galaxies reported by Zavala et al. (2018). Finally, we also include the number counts from the ASPECS survey (González-López et al. 2019, 2020), a deep ALMA large program at 1.2 and 3 mm. Altogether, these surveys add a total of ~400 arcmin² of deep observations with arcsecond resolution, representing the state-of-the-art blind ALMA surveys to date.

All these data provide good constraints on the model parameters and thus on the IR luminosity function and its evolution with redshift. Based on our best-fit model, which simultaneously reproduces the FIR and (sub)millimeter data from single-dish telescopes and interferometric surveys, we constrained the faint-end slope of the IR luminosity function to be flat, with a value of $\alpha_{LF} = -0.42_{-0.04}^{+0.02}$ (see also Koprowski et al. 2017; Lim et al. 2020). This implies that deep pencil-beam observations at (sub)millimeter wavelengths would not significantly increase the number of detected sources, in line with previous results from the literature (e.g., Popping et al. 2020). The characteristic number density of the luminosity function, $\Phi_{\star}$, decreases as $\Phi_{\star} \propto (1+z)^{-6.5_{-1.8}^{+0.8}}$ at $z \gtrsim 2$, in a

similar fashion to the quasar luminosity function and the density evolution of UV-bright galaxies (Hopkins et al. 2007; Finkelstein et al. 2015), which might suggest that all these galaxies occupy similar dark matter halos.

Our constraints on the dust-obscured star formation indicate that the cosmic history of star formation had a peak at $z \approx 2-2.5$ and has been dominated by the dust-obscured component during the past 12 billion years, back to $z \sim 4$, when the unobscured and obscured contributions were comparable. Beyond this epoch, the dust reprocessed SFRD was less dominant than the visible star formation, contributing around $35_{-25}^{+10}$% at $z=5$ and $25_{-20}^{+15}$% at $z=6$. This suggests that the bulk (>$80$%) of the star formation activity in the first billion years of the universe was not dust enshrouded. Given the massive nature of DSFGs, this drop-off of the obscured component is in line with the decreasing number of high-mass galaxies with increasing redshift (see also Dunlop et al. 2017; Bouwens et al. 2020).

Our picture of the history of the cosmic star formation is consistent with previous results from long-wavelength (FIR to millimeter) surveys (e.g., Bourne et al. 2017; Dunlop et al. 2017; Hatsukade et al. 2018; Dudzevičiūtė et al. 2020), although most of those were limited to $z \lesssim 4$. Hence, our results represent significant progress on our understanding of the prevalence of DSFGs during the first 1.5 billion years of the universe and complement the significant efforts carried out using UV/optically selected galaxies (e.g., Bouwens et al. 2020). The inferred SFRD is also in broad agreement with the most recent predictions from galaxy evolution models (like IllustrisTNG and SHARK; Pillepich et al. 2018; Lagos et al. 2020), which point toward a convergence between models and observations.

With estimations for both the obscured and unobscured components, the shape of the cosmic history of star formation is now constrained out to the end of the epoch of reionization. This measurement, which preserves the galaxy mass assembly history of the universe, provides a benchmark against which to compare galaxy formation models and simulations, as well as a step forward in our understanding of the dust and metal enrichment of the early universe.

We thank the anonymous referee for a careful review of our manuscript and his/her valuable suggestions.

This paper makes use of the following ALMA data: ADS/JAO. ALMA#2018.1.00231.S, ADS/JAO.ALMA#2018.1.00478.S, and ADS/JAO.ALMA#2019.1.00838.S. ALMA is a partnership of 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 ESO, AUI/NRAO, and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Dunlap Institute is funded through an endowment established by the David Dunlap family and the University of Toronto. J.A.Z. and C.M.C. thank the University of Texas at Austin College of Natural Sciences for support. C.M.C. also thanks the National Science Foundation for support through grants AST-1714528 and AST-1814034 and support from the Research Corporation for Science Advancement from a 2019 Cottrell Scholar Award sponsored by IF/THEN, an initiative of Lyda Hill Philanthropies. M.A. has been supported by the grant "CONICYT + PCI + INSTITUTO MAX PLANCK DE
ASTRONOMIA MPG190030" and "CONICYT+PCI+REDES 190194." K.I.C. acknowledges funding from the European Research Council through the award of the Consolidator Grant ID 681627-BUILDUP. J.H. acknowledges support of the VIDI research program with project No. 639.042.611, which is (partly) financed by the Netherlands Organisation for Scientific Research (NWO). K.K. acknowledges support from the Knut and Alice Wallenberg Foundation. G.E.M. acknowledges the Villum Fonden research grant 13160 "Gas to stars, stars to dust: tracing star formation across cosmic time" and the Cosmic Dawn Center of Excellence funded by the Danish National Research Foundation under then grant No. 140. E.T. acknowledges support from FONDECYT Regular 1190818, CONICYT PIA ACT172033, and Basal-CATA AFB170002 grants.

Finally, our deepest gratitude to all the people who continue
to provide essential services during the COVID-19 pandemic.

Data and Materials Availability

This paper makes use of the following ALMA data:
ADS/JAO. ALMA#2018.1.00231.S, ADS/JAO.ALMA
#2018.1.00478.S, and ADS/JAO.ALMA#2019.1.00838.S,
archived at https://almascience.nrao.edu/alma-data/archive.

**Appendix A**

Number Counts Model Predictions at Shorter Wavelengths

As discussed in Casey et al. (2018a, 2018b), the number
counts at $\lambda \lesssim 850$ µm mostly inform about the evolution of the

IRLF at $z \lesssim 3.0$, where the model parameters are already relatively well understood. Given that our primary goal is to constrain the IRLF at earlier epochs, only the number counts at 1.2, 2, and 3 mm were used in our analysis. Additionally, incorporating all the available number counts in the fitting procedure, particularly those from single-dish telescope surveys, would increase the computational cost significantly since the simulated area would increase by two orders of magnitude. This would make our fitting analysis prohibitive.

Nevertheless, as mentioned in the main text, our best-fit
model nicely reproduces the number counts at shorter
wavelengths, from $\lambda = 70$ to 850 µm, covering a large dynamic
range of flux density. In Figure 9, three different examples are
shown: the number counts at 250, 450/500, and 850 µm. The
remaining number counts can be found in Casey et al. (2018b).
Since most of the galaxies detected at these short wavelengths
are relatively bright galaxies at $z \lesssim 3$, the best-fit parameters
studied in this work (which govern the faint and high-redshift
population) do not significantly change the model predictions
at those wavelengths compared to Casey et al. (2018b).

**Appendix B**

Luminosity Limits of the Adopted Surveys

The luminosity sensitivity limits of the surveys used in this
work to constrain the IRLF and dust-obscured SFRD are
illustrated in Figure 10.

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

**Figure 9.** Differential number counts at 250 µm (left panel), 450 µm/500 µm (middle panel), and 850 µm (right panel). The predictions from the best-fit model are represented by gray lines, while data points from the literature are illustrated as blue circles. Note that in order to minimize the computational cost, only predictions from a model with $\alpha_{LF} = -0.42$, $\Psi_2 = -6.5$, and $\beta_{em} = 1.8$ are included, i.e., without taking into account the uncertainties in the best-fit model parameters. Generally speaking, the model reproduces all the FIR and submillimeter data aggregated over the past two decades from both single-dish telescopes and interferometric surveys spanning more than a decade in wavelength and orders of magnitude in flux density. At both 250 µm and 500 µm, the data come from BLAST (Patanchon et al. 2009; Béthermin et al. 2010) and Herschel (Clements et al. 2010; Oliver et al. 2010; Béthermin et al. 2012). The 450 µm data points (which have not been scaled relative to the 500 µm flux density) come from Casey et al. (2013), Chen et al. (2013), Geach et al. (2013), Hsu et al. (2016), Wang et al. (2017), and Zavala et al. (2017). Finally, the 850 µm panel includes the works of Chapman et al. (2002), Webb et al. (2003), Coppin et al. (2006), Scott et al. (2006), Beelen et al. (2008), Knudsen et al. (2008), Weiß et al. (2009), Casey et al. (2013), Chen et al. (2013), Karim et al. (2013), Hsu et al. (2016), and Zavala et al. (2017).
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Luminosity sensitivity limits of the three different surveys adopted in this work. Given that we use directly the galaxy number counts to constrain our model, these luminosity limits were calculated adopting the minimum flux density bin of each of the galaxy number counts (see Figure 5), respectively. The deepest number counts at 1.1 mm are sensitive to galaxies down to $L_{IR} \sim 2 \times 10^{10} L_{\odot}$ up to $z \sim 8$, while the larger-area surveys at 2 and 3 mm are sensitive to brighter and rarer galaxies with $L_{IR} > 10^{11} - 10^{12}$. Note that, although the cosmic SFRD presented in Figure 7 was calculated by integrating the best-fit IR luminosity function over the interval $\log(L/L_{\odot}) = [9, 13.8]$ (i.e., extrapolating beyond the surveys’ luminosity limits), the majority of the contribution arises from galaxies with luminosities in the range of $L_{IR} = 10^{11} - 10^{13} L_{\odot}$, well within the luminosity limits of the adopted data.</figcaption></figure>

## ORCID iDs

J. A. Zavala https://orcid.org/0000-0002-7051-1100

C. M. Casey https://orcid.org/0000-0002-0930-6466

S. M. Manning https://orcid.org/0000-0003-0415-0121

M. Aravena https://orcid.org/0000-0002-6290-3198

M. Bethermin https://orcid.org/0000-0002-3915-2015

K. I. Caputi https://orcid.org/0000-0001-8183-1460

D. L. Clements https://orcid.org/0000-0002-9548-5033

E. da Cunha https://orcid.org/0000-0001-9759-4797

P. Drew https://orcid.org/0000-0003-3627-7485

S. L. Finkelstein https://orcid.org/0000-0001-8519-1130

S. Fujimoto https://orcid.org/0000-0001-7201-5066

C. Hayward https://orcid.org/0000-0003-4073-3236

J. Hodge https://orcid.org/0000-0001-6586-8845

J. S. Kartaltepe https://orcid.org/0000-0001-9187-3605

K. Knudsen https://orcid.org/0000-0002-7821-8873

A. M. Koekemoer https://orcid.org/0000-0002-6610-2048

A. S. Long https://orcid.org/0000-0002-7530-8857

G. E. Magdis https://orcid.org/0000-0002-4872-2294

A. W. S. Man https://orcid.org/0000-0003-2475-124X

G. Popping https://orcid.org/0000-0003-1151-4659

D. Sanders https://orcid.org/0000-0002-1233-9998

N. Scoville https://orcid.org/0000-0002-0438-3323

K. Sheth https://orcid.org/0000-0002-5496-4118

J. Staguhn https://orcid.org/0000-0002-8437-0433

S. Toft https://orcid.org/0000-0003-3631-7176

E. Treister https://orcid.org/0000-0001-7568-6412

J. D. Vieira https://orcid.org/0000-0001-7192-3871

M. S. Yun https://orcid.org/0000-0001-7095-7543

## References

Aravena, M., Boogaard, L., Gómez-López, J., et al. 2020, ApJ, 901, 79

Baugh, C. M., Lacey, C. G., Frenk, C. S., et al. 2005, MNRAS, 356, 1191

Beelen, A., Omont, A., Bavouzet, N., et al. 2008, A&A, 485, 645

Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143

Bethermin, M., Daddi, E., Magdis, G., et al. 2015, A&A, 573, A113

Béthermin, M., Dole, H., Cousin, M., & Bavouzet, N. 2010, A&A, 516, A43

Béthermin, M., Le Floc'h, E., Ilbert, O., et al. 2012, A&A, 542, A58

Béthermin, M., Wu, H.-Y., Lagache, G., et al. 2017, A&A, 607, A89

Bourne, N., Dunlop, J. S., Merlin, E., et al. 2017, MNRAS, 467, 1360

Bouwens, R., Gonzalez-Lopez, J., Aravena, M., et al. 2020, ApJ, 902, 112

Bouwens, R. J., Aravena, M., Decarli, R., et al. 2016, ApJ, 833, 72

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2012, ApJL, 752, L5

Capak, P. L., Carilli, C., Jones, G., et al. 2015, Natur, 522, 455

Casey, C. M. 2012, MNRAS, 425, 3094

Casey, C. M., Berta, S., Béthermin, M., et al. 2012, ApJ, 761, 140

Casey, C. M., Chen, C.-C., Cowie, L. L., et al. 2013, MNRAS, 436, 1919

Casey, C. M., Hodge, J., Zavala, J. A., et al. 2018a, ApJ, 862, 78

Casey, C. M., Zavala, J. A., Spilker, J., et al. 2018b, ApJ, 862, 77

Chabrier, G. 2003, PASP, 115, 763

Chapman, S. C., Scott, D., Borys, C., & Fahlman, G. G. 2002, MNRAS, 330, 92

Chen, C.-C., Cowie, L. L., Barger, A. J., et al. 2013, ApJ, 776, 131

Chiang, Y.-K., Overzier, R. A., Gebhardt, K., & Henriques, B. 2017, ApJL, 844, L23

Clements, D. L., Rigby, E., Maddox, S., et al. 2010, A&A, 518, L8

Coppin, K., Chapin, E. L., Mortier, A. M. J., et al. 2006, MNRAS, 372, 1621

Cucciati, O., Tresse, L., Ilbert, O., et al. 2012, A&A, 539, A31

da Cunha, E., Groves, B., Walter, F., et al. 2013, ApJ, 766, 13

Dahlen, T., Mobasher, B., Dickinson, M., et al. 2007, ApJ, 654, 172

Dudzevičiūtė, U., Smail, I., Swinbank, A. M., et al. 2020, MNRAS, 494, 3828

Dunlop, J. S., McLure, R. J., Biggs, A. D., et al. 2017, MNRAS, 466, 861

Dunne, L., Gomez, H. L., da Cunha, E., et al. 2011, MNRAS, 417, 1510

Ferguson, H. C., Dickinson, M., & Williams, R. 2000, ARA&A, 38, 667

Finkelstein, S. L., Ryan, R. E. J., Papovich, C., et al. 2015, ApJ, 810, 71

Franco, M., Elbaz, D., Béthermin, M., et al. 2018, A&A, 620, A152

Fudamoto, Y., Oesch, P. A., Faisst, A., et al. 2020, A&A, 643, A4

Fujimoto, S., Ouchi, M., Ono, Y., et al. 2016, ApJS, 222, 1

Galliano, F., Galametz, M., & Jones, A. P. 2018, ARA&A, 56, 673

Geach, J. E., Chapin, E. L., Coppin, K. E. K., et al. 2013, MNRAS, 432, 53

Gehrels, N. 1986, ApJ, 303, 336

González-López, J., Decarli, R., Pavesi, R., et al. 2019, ApJ, 882, 139

González-López, J., Novak, M., Decarli, R., et al. 2020, ApJ, 897, 91

Gruppioni, C., Bethermin, M., Loiacono, F., et al. 2020, A&A, 643, A8

Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 432, 23

Hatsukade, B., Kohno, K., Yamaguchi, Y., et al. 2018, PASJ, 70, 105

Hayward, C. C., Narayanan, D., Kereš, D., et al. 2013, MNRAS, 428, 2529

Hopkins, P. F., Richards, G. T., & Hernquist, L. 2007, ApJ, 654, 731

Hsu, L.-Y., Cowie, L. L., Chen, C.-C., Barger, A. J., & Wang, W.-H. 2016, ApJ, 829, 25

Jin, S., Daddi, E., Magdis, G. E., et al. 2019, ApJ, 887, 144

Karim, A., Swinbank, A. M., Hodge, J. A., et al. 2013, MNRAS, 432, 2

Knudsen, K. K., van der Werf, P. P., & Kneib, J. P. 2008, MNRAS, 384, 1611

Koprowski, M. P., Dunlop, J. S., Michalowski, M. J., et al. 2017, MNRAS, 471, 4155

Lagos, C. d. P., da Cunha, E., Robotham, A. S. G., et al. 2020, MNRAS, 499, 1948

Lagos, C. d. P., Tobar, R. J., Robotham, A. S. G., et al. 2018, MNRAS, 481, 3573

Laporte, N., Ellis, R. S., Boone, F., et al. 2017, ApJL, 837, L21

Leslie, S., Schinnerer, E., Liu, D., et al. 2020, ApJ, 899, 58

Lim, C.-F., Wang, W.-H., Smail, I., et al. 2020, ApJ, 889, 80

Liu, D., Daddi, E., Dickinson, M., et al. 2018, ApJ, 853, 172

Loiacono, F., Decarli, R., Gruppioni, C., et al. 2020, arXiv:2006.04837

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Magdis, G. E., Rigopoulou, D., Daddi, E., et al. 2017, A&A, 603, A93

Magnelli, B., Elbaz, D., Chary, R. R., et al. 2011, A&A, 528, A35

Magnelli, B., Karim, A., Staguhn, J., et al. 2019, ApJ, 877, 45

Magnelli, B., Popesso, P., Berta, S., et al. 2013, A&A, 553, A132

McLure, R. J., Dunlop, J. S., Bowler, R. A. A., et al. 2013, MNRAS, 432, 2696

Michalowski, M. J., Dunlop, J. S., Koprowski, M. P., et al. 2017, MNRAS, 469, 492

Moster, B. P., Somerville, R. S., Newman, J. A., & Rix, H.-W. 2011, ApJ, 731, 113

Nelder, J. A., & Mead, R. 1965, CompJ, 7, 308

Novak, M., Smolčić, V., Delhaize, J., et al. 2017, A&A, 602, A5

Oesch, P. A., Bouwens, R. J., Illingworth, G. D., Labbé, I., & Stefanon, M. 2018, ApJ, 855, 105
Oliver, S. J., Wang, L., Smith, A. J., et al. 2010, A&A, 518, L21

Oteo, I., Zwaan, M. A., Ivison, R. J., Smail, I., & Biggs, A. D. 2016, ApJ, 822, 36

Patanchon, G., Ade, P. A. R., Bock, J. J., et al. 2009, ApJ, 707, 1750

Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13

Popping, G., Walter, F., Behroozi, P., et al. 2020, ApJ, 891, 135

Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes in C: The Art of Scientific Computing (2nd ed.; Cambridge: Cambridge Univ. Press)

Reddy, N. A., & Steidel, C. C. 2009, ApJ, 692, 778

Remy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2013, A&A, 557, A95

Robotham, A. S. G., & Driver, S. P. 2011, MNRAS, 413, 2570

Rowan-Robinson, M., Oliver, S., Wang, L., et al. 2016, MNRAS, 461, 1100

Sanders, D. B., Mazzarella, J. M., Kim, D. C., Surace, J. A., & Soifer, B. T. 2003, AJ, 126, 1607

Schenker, M. A., Robertson, B. E., Ellis, R. S., et al. 2013, ApJ, 768, 196

Schiminovich, D., Ilbert, O., Arnouts, S., et al. 2005, ApJL, 619, L47

Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74

Scott, K. S., Wilson, G. W., Aretxaga, I., et al. 2012, MNRAS, 423, 575

Scott, S. E., Dunlop, J. S., & Serjeant, S. 2006, MNRAS, 370, 1057

Stach, S. M., Smail, I., Swinbank, A. M., et al. 2018, ApJ, 860, 161

Staguhn, J. G., Kovács, A., Arendt, R. G., et al. 2014, ApJ, 790, 77

Sunyaev, R. A., & Zeldovich, Y. B. 1970, CoASP, 2, 66

Swinbank, A. M., Simpson, J. M., Smail, I., et al. 2014, MNRAS, 438, 1267

Umehata, H., Tamura, Y., Kohno, K., et al. 2017, ApJ, 835, 98

Wall, J. V., Jackson, C. A., Shaver, P. A., Hook, I. M., & Kellermann, K. I. 2005, A&A, 434, 133

Wang, F., Yang, J., Fan, X., et al. 2019a, ApJ, 884, 30

Wang, L., Pearson, W. J., Cowley, W., et al. 2019b, A&A, 624, A98

Wang, W.-H., Lin, W.-C., Lim, C.-F., et al. 2017, ApJ, 850, 37

Watson, D., Christensen, L., Knudsen, K. K., et al. 2015, Natur, 519, 327

Webb, T. M., Eales, S. A., Lilly, S. J., et al. 2003, ApJ, 587, 41

Weiß, A., Kovács, A., Coppin, K., et al. 2009, ApJ, 707, 1201

Williams, C. C., Labbe, I., Spilker, J., et al. 2019, ApJ, 884, 154

Wilner, D. J., & Wright, M. C. H. 1997, ApJL, 488, L67

Wyder, T. K., Treyer, M. A., Milliard, B., et al. 2005, ApJL, 619, L15

Yamaguchi, Y., Kohno, K., Hatsukade, B., et al. 2019, ApJ, 878, 73

Zavala, J. A., Aretxaga, I., Geach, J. E., et al. 2017, MNRAS, 464, 3369

Zavala, J. A., Casey, C. M., da Cunha, E., et al. 2018, ApJ, 869, 71