# Mapping Obscuration to Reionization with ALMA (MORA): 2 mm Efficiently Selects the Highest-redshift Obscured Galaxies

Caitlin M. Casey, Jorge A. Zavala, Sinclaire M. Manning, Manuel Aravena, Matthieu Béthermin, Karina I. Caputi, Jaclyn B. Champagne, David L. Clements, Patrick Drew, Steven L. Finkelstein, et al.

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

► To cite this version:

Caitlin M. Casey, Jorge A. Zavala, Sinclaire M. Manning, Manuel Aravena, Matthieu Béthermin, et al. Mapping Obscuration to Reionization with ALMA (MORA): 2 mm Efficiently Selects the Highest-redshift Obscured Galaxies. The Astrophysical Journal, 2021, 923, 32 pp. 10.3847/1538-4357/ac2eb4. hal-03585959

HAL Id: hal-03585959

https://hal.science/hal-03585959

Submitted on 24 Feb 2022

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>
<figure><img src="image_3.png" /></figure>

# Mapping Obscuration to Reionization with ALMA (MORA): 2 mm Efficiently Selects the Highest-redshift Obscured Galaxies

Caitlin M. Casey¹, Jorge A. Zavala¹,², Sinclaire M. Manning¹,³,<fnref n="23" />, Manuel Aravena⁴, Matthieu Béthermin⁵, Karina I. Caputi⁶,⁷, Jaclyn B. Champagne¹, David L. Clements⁸, Patrick Drew¹, Steven L. Finkelstein¹, Seiji Fujimoto⁷,⁹, Christopher C. Hayward¹⁰, Anton M. Dekel¹¹, Vasily Kokorev⁷,⁹, Claudia del P. Lagos⁷,¹²,¹³, Arianna S. Long¹⁴, Georgios E. Magdis⁷,⁹,¹⁵, Allison W. S. Man¹⁶, Ikki Mitsuhashi²,¹⁷, Gergö Popping¹⁸, Justin Spilker¹,<fnref n="23" />, Johannes Staguhn¹⁹,²⁰, Margherita Talia²¹, Sune Toft⁷,⁹, Ezequiel Treister²², John R. Weaver⁷,⁹, and Min Yun³

¹ Department of Astronomy, The University of Texas at Austin, 2515 Speedway Blvd Stop C1400, Austin, TX 78712, USA; cmcasey@utexas.edu

² National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

³ Department of Astronomy, University of Massachusetts Amherst, 710 N. Pleasant St, Amherst, MA 01003, USA

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

⁵ Aix-Marseille Université, CNRS, CNES, LAM, Marseille, France

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

⁷ Cosmic Dawn Center (DAWN)

⁸ Imperial College London, Blackett Laboratory, Prince Consort Rd, London, SW7 2AZ, UK

⁹ Niels Bohr Institute, University of Copenhagen, Lyngbyvej 2, DK-2200 Copenhagen, Denmark

¹⁰ Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Ave, New York, NY 10010, USA

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

¹² International Centre for Radio Astronomy Research (ICRAR), M468, University of Western Australia, 35 Stirling Hwy, Crawley, WA 6009, Australia

¹³ ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D)

¹⁴ 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

¹⁶ Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Rd, Vancouver, BC V6T 1Z1, Canada

¹⁷ Department of Astronomy, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan

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

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

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

²¹ Dipartimento di Fisica e Astronomia, Università di Bologna, Via Gobetti 93/2, I-40129, Bologna, Italy

²² Instituto de Astrofísica and Centro de Astroingeniería, Facultad de Física, Pontificia Universidad Católica de Chile, Casila 306, Santiago 22, Chile
Received 2021 July 28; revised 2021 October 1; accepted 2021 October 9; published 2021 December 23

## Abstract

We present the characteristics of 2 mm selected sources from the largest Atacama Large Millimeter/submillimeter Array (ALMA) blank-field contiguous survey conducted to date, the Mapping Obscuration to Reionization with ALMA (MORA) survey covering 184 arcmin² at 2 mm. Twelve of 13 detections above 5σ are attributed to emission from galaxies, 11 of which are dominated by cold dust emission. These sources have a median redshift of $\langle z_{2 \text{ mm}} \rangle = 3.6_{-0.3}^{+0.4}$ primarily based on optical/near-infrared photometric redshifts with some spectroscopic redshifts, with 77% ± 11% of sources at $z > 3$ and 38% ± 12% of sources at $z > 4$. This implies that 2 mm selection is an efficient method for identifying the highest-redshift dusty star-forming galaxies (DSFGs). Lower-redshift DSFGs ($z < 3$) are far more numerous than those at $z > 3$ yet are likely to drop out at 2 mm. MORA shows that DSFGs with star formation rates in excess of 300 $M_{\odot} \text{yr}^{-1}$ and a relative rarity of $\sim 10^{-5} \text{ Mpc}^{-3}$ contribute $\sim 30\%$ to the integrated star formation rate density at $3 < z < 6$. The volume density of 2 mm selected DSFGs is consistent with predictions from some cosmological simulations and is similar to the volume density of their hypothesized descendants: massive, quiescent galaxies at $z > 2$. Analysis of MORA sources' spectral energy distributions hint at steeper empirically measured dust emissivity indices than reported in typical literature studies, with $\langle \beta \rangle = 2.2_{-0.4}^{+0.5}$. The MORA survey represents an important step in taking census of obscured star formation in the universe's first few billion years, but larger area 2 mm surveys are needed to more fully characterize this rare population and push to the detection of the universe's first dusty galaxies.

*Unified Astronomy Thesaurus concepts: Millimeter astronomy (1061); Submillimeter astronomy (1647); Dust continuum emission (412); High-redshift galaxies (734); Active galaxies (17); Infrared galaxies (790); Starburst galaxies (1570)*

## 1. Introduction

Half of all extragalactic radiation is absorbed by dust and re-emitted at long wavelengths (e.g., Fixsen et al. 1998). Decades of progress, both technological and observational, have taught us that the obscured emission emanates from very different galaxies than unobscured light; the former is largely from massive, star-forming galaxies while the latter is from lower-mass galaxies (Whitaker

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

<footnote n="23">²³ NHFP Hubble Fellow.</footnote>

<footnote>Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.</footnote>
et al. 2017). So while the need to take census of star formation has been a key focus of extragalactic astrophysics for some time (e.g., Madau & Dickinson 2014), it has been clear that the very deep surveys of cosmic star formation conducted in the rest-frame ultraviolet and optical may not adequately capture the full picture. Surveys of submillimeter-luminous dusty star-forming galaxies (DSFGs; e.g., Smail et al. 1997; Blain et al. 2002; Casey et al. 2014)—star-forming galaxies with star formation rates (SFRs) $\gtrsim 100 M_{\odot}$ yr$^{-1}$ whose stellar emission is over $\gtrsim 95\%$ obscured by dust—have been the primary method of unveiling the universe’s obscured contribution to the star formation rate density (SFRD). This approach is in direct contrast with the strategy of measuring the total SFRD by taking a census of UV-selected galaxies and correcting estimates for dust attenuation (e.g., as in Bouwens et al. 2020).

A key limitation in all surveys of distant, dust obscured galaxies is the difficulty in identifying their redshifts. Unlike Lyman Break Galaxies (LBGs), whose selection method directly indicates their redshifts, DSFGs’ spectral shape in the (sub)millimeter regime is highly degenerate with redshift solutions spanning $1 \lesssim z \lesssim 12$. Similar efforts to characterize obscured emission indirectly via synchrotron radio emission are faced with similar challenges, though primarily limited to $1 \lesssim z \lesssim 4$ (Novak et al. 2017). Following up obscured sources in the optical or near-infrared for characteristic emission lines present in star-forming galaxies is difficult due to significant obscuration (Chapman et al. 2003, 2005; Swinbank et al. 2004). Pursuing millimeter spectroscopy has been prohibitive for large samples until recently due to technological and sensitivity limitations in available instrumentation.<fnref n="24" /> In addition, the large beam sizes of many single-dish (sub)millimeter facilities can further obfuscate redshift identification via uncertain astrometry and source confusion, requiring another intermediate stage of interferometric follow-up to constrain positions (e.g., Hodge et al. 2013; Karim et al. 2013).

The complexities in identifying DSFGs’ redshifts has led to great difficulty in measuring the volume density of highly obscured galaxies beyond $z \gtrsim 3-4$. A small error in redshift measurement for a small fraction of a uniformly selected DSFG sample (typically selected at $\lambda \le 1$ mm) can result in very different inferred volume densities for the population at these redshifts (e.g., see the wide variety of high-$z$ volume density measurements in Rowan-Robinson et al. 2016; Koprowski et al. 2017; Dudzevičiūtė et al. 2020; Gruppioni et al. 2020; Khusanova et al. 2021; Loiacono et al. 2021). This is primarily because the peak in the redshift distribution for 850 $\mu$m selected DSFGs (as well as 1 mm selected DSFGs) is between $2 < z < 3$, and sources at higher redshifts are quite rare per unit solid angle on the sky in comparison. This has been demonstrated throughout the literature via the measurement of the redshift distribution of DSFGs selected at $\sim 850 \mu$m-$1.2$ mm (e.g., Smolčić et al. 2012; Hatsukade et al. 2013; Brisbin et al. 2017). Recent work by Dudzevičiūtė et al. (2020) points out that only 6% of 850 $\mu$m selected galaxies sit at $z > 4$. Thus, the accurate identification of such systems—often having degenerate submillimeter colors with sources at lower redshifts—is effectively equivalent to searching for a needle in a haystack.

Some efforts have focused on selecting the highest-redshift DSFGs using submillimeter colors, identifying characteristically

“red” spectral energy distributions (SEDs) across Herschel bands (e.g., Dowell et al. 2014; Ivison et al. 2016; Bakx et al. 2018; Donevski et al. 2018; Duivenvoorden et al. 2018; Yan et al. 2020). However, Herschel bands do not benefit from the negative K-correction because they probe emission near the peak of the dust SED rather than the Rayleigh-Jeans tail (Blain et al. 2002; Casey et al. 2014); thus Herschel data sets tend to have reduced sensitivity to un lensed DSFGs beyond $z \sim 2-3$. Furthermore, this technique may select against high-redshift DSFGs with atypically warm SEDs, due to the degeneracy between dust temperature (or SED peak wavelength, $\lambda_{\text{peak}}$) and redshift. There have also been some attempts to constrain the $z \gtrsim 3$ DSFG population through measurements of anisotropies in the Cosmic Infrared Background (CIB; Maniyar et al. 2018, 2021), though Herschel is largely insensitive to the high-redshift tail.

This paper presents data from a new large ALMA mosaic conducted at 2 mm (band 4), whose aim is to efficiently select DSFGs at $z \gtrsim 3-4$ and measure their volume density at these epochs with more precision than has been done before. This program is called the Mapping Obscuration to Reionization with ALMA (MORA) Survey, for its focus on this especially early epoch of DSFG formation.

The MORA survey is based on the hypothesis that 2 mm dust continuum is an efficient “filter” for $z \gtrsim 3-4$ systems. Surveys conducted at 850 $\mu$m-$1.2$ mm benefit from the very negative K-correction, meaning that their expected flux density is not redshift dependent for a given fixed IR luminosity (i.e., at fixed $L_{\text{IR}}$, $S_{850} \approx C$, where C is a constant without redshift dependence). Moreover, surveys at 2 mm have an even more extreme negative K-correction, such that 2 mm flux density *increases* with redshift (i.e., at fixed $L_{\text{IR}}$, $S_{2\text{mm}} \propto (1+z)^n$, where $n \approx 0.5-1.0$). As a result of this extreme negative K-correction, $z \gtrsim 3$ galaxies of matched luminosity should appear *brighter* at 2 mm than those at $z \sim 1-3$. This contrasts with nearly every other wave band in which galaxies are observed, from the X-ray through the radio, where more distant objects are expected to be fainter than objects closer to us. If a blank-field 2 mm survey depth is adjusted appropriately, it should be sensitive to detecting DSFGs at $z \gtrsim 3$ while the much more common $z \sim 1-3$ DSFGs should be undetected (see detailed modeling work in Casey et al. 2018a, 2018b; Zavala et al. 2018a). This is the premise for the design of the MORA survey.

A parallel work to this paper is presented by Zavala et al. (2021), hereafter referred to as Z21, who conduct a number counts analysis of the MORA Survey 2 mm mosaics and study implications for the integrated cosmic SFRD at $z \gtrsim 3$.

This paper presents the characteristics of the individual sources detected in the MORA Survey and what can be inferred about their redshifts, masses, SEDs, and descendants. Section 2 presents the details of the MORA Survey design and data acquisition. Section 3 then presents the identification of the robust sample of 2 mm selected galaxies, and describes details of each source from what is known in the literature; this section also presents a discussion of sources with low signal-to-noise ratio (S/N) detections below the formal $5\sigma$ detection threshold. Section 4 describes models of the 2 mm universe, including semi-empirical models based in cosmological simulations and empirical models. Section 5 presents the main calculations and results of our manuscript, including analysis of the 2 mm population redshift distribution, SEDs, and other unresolved physical characteristics. Section 6 presents a discussion, including results of the contribution of the MORA sample to the cosmic SFRD, an extrapolation of what the MORA galaxy sample will evolve to

<footnote n="24">²⁴ And recently, large samples are still quite challenging to confirm via millimeter spectral scans, as it requires anywhere from 30 minutes to 3 hr of integration time per source, even for luminous (unlensed) DSFGs (see Vieira et al. 2013).</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 1.** The positions of the MORA mosaics shown in teal relative to the entire COSMOS field. The moon is provided for scale in the lower left. Shown are the HST/ACS F814W imaging (Koekemoer et al. 2007; Scoville et al. 2007), the Chandra-COSMOS deep and medium deep survey areas in red (Civano et al. 2016), the CANDELS-COSMOS area in yellow (Grogin et al. 2011; Koekemoer et al. 2011), the deep 450 and 850 µm SCUBA-2 pointings in orange from the STUDIES program (Wang et al. 2017) and S2CLS (Roseboom et al. 2013; Geach et al. 2017), and the MORA mosaics in teal. We also show the proposed MORA area (dotted teal square) and the forthcoming PRIMER and COSMOS-Web deep NIRCam imaging programs (dotted purple lines). Other multiwavelength data are not shown in this figure for clarity, though many cover the full COSMOS field (e.g., deep 3 GHz radio continuum from Smolčić et al. 2017). Most spectroscopy in the field is concentrated within the central $1\deg^2$ around the Chandra-Deep and COSMOS-Web footprints.</figcaption></figure>

become, discussion of the measured emissivity spectral index, and the potential impact of cosmic variance on our measurements. Section 7 then presents our conclusions. Throughout we assume a standard $\Lambda$CDM cosmology adopting the Planck-measured parameters, with $H_0 = 67.4 \text{ km s}^{-1} \text{ Mpc}^{-1}$ (Planck Collaboration et al. 2020), and where SFRs are mentioned, we assume a Kroupa initial mass function (Kroupa & Weidner 2003) and scaling relations drawn from Kennicutt & Evans (2012).

## 2. Data and Observations

MORA Survey observations were originally designed to cover $230 \text{ arcmin}^2$ in two tunings, both in ALMA band 4 at 2 mm in the center of the Cosmic Evolution Survey (COSMOS) field (Capak et al. 2007; Koekemoer et al. 2007; Scoville et al. 2007). Figure 1 shows the context of the proposed and observed MORA mosaics in the larger COSMOS field relative to other key data sets in the field. The COSMOS field was chosen as the location of the mosaic due to its rich multiwavelength data (discussed more in Section 3.1) and, specifically, the CANDELS portion of the COSMOS field (Grogin et al. 2011; Koekemoer et al. 2011) was chosen for its even deeper near-infrared imaging with Hubble Space Telescope (HST) WFC3. The near-infrared depth will be significantly enhanced with the addition of James Webb Space Telescope (JWST) data from the PRIMER and COSMOS-Web surveys.

The target continuum rms for the entire MORA mosaic was $90 \mu\text{Jy beam}^{-1}$ at $1\sigma$. This exact tuning configuration in band 4 was chosen based on the program's secondary goal: to search the

known $z \sim 2.5$ protocluster structure, “Hyperion”, which spatially overlaps with this map (Casey et al. 2015; Chiang et al. 2015; Diener et al. 2015; Cucciati et al. 2018), for a blind search of molecular and neutral gas emitters. The first tuning was centered on a local oscillator (LO) frequency of 147.28 GHz (referred to as “Tune147”) and covered the frequency ranges 139.5–143.2 GHz and 151.4–155.2 GHz. This tuning is sensitive to the detection of $C I(1-0)$ at $2.44 < z < 2.52$. The second tuning is centered on an LO frequency of 139.03 GHz and covered the frequency ranges 131.2–134.9 GHz and 143.2–146.9 GHz (referred to as “Tune139”). It is tuned to enable the detection of $CO(4-3)$ at $2.42 < z < 2.51$. This resulted in 21 scheduling blocks (SBs) of 149 pointings each for each tuning, resulting in 42 total SBs. Each SB was spatially distributed as 2 columns of pointings with fixed R.A. and 74–75 pointings in decl. spanning a decl. range $+02:11:04$ to $+02:33:56$. Each SB’s fixed R.A. position is referred to as a position and a number in this text, e.g., “PXX” where XX ranges from 03–20. As proposed, the mosaic would have covered a total of 3129 pointings at two frequency settings each. The individual pointings of the mosaic were spaced by $19''/3$, which is 0.47 times the primary beam FWHM at the highest frequency of data acquisition, 155.2 GHz. This spacing leads to a slightly more compact mosaic than the default Nyquist spacing for mosaics; the same spacing was used for both tunings to make data processing more straightforward. This consequently resulted in a greater depth of observations than proposed (as Nyquist sampling was used to derive the on-source time).

Observations with the Atacama Large Millimeter and submillimeter Array (ALMA) were carried out under program 2018.1.00231.S from 2019 March 27 through 2019 April 3 in the C43-3 configuration for a total of 14.6 hr including overheads and calibrations. Data were acquired under an average precipitable water vapor of PWV = 5 mm with conditions ranging from 2 mm < PWV < 6.5 mm.

The program was observed in part only: 14 of the 42 SBs were executed, 11 of which were at the higher frequency tuning, Tune147, and three of which were at the lower frequency tuning, Tune139. One of the higher frequency SBs, Tune147 for position P16, did not pass QA0 due to poor weather conditions, but was processed after the fact manually as “semi-pass” data and folded into the final mosaic after flagging problematic antennae. Table 1 lists the observational conditions and data characteristics of each observed SB and the final mosaics.

There are two final mosaics produced from these data that are spatially distinct: one elongated mosaic represents observations taken in the “P03” position (two pointings wide) with both tunings over a total area of $28 \text{ arcmin}^2$. P03 is too far spatially offset from the rest of the data to be joined in one mosaic. The other mosaic represents all other data from the spatially adjacent positions “P10–P20” over a total area of $156 \text{ arcmin}^2$. We refer to these as the “P03” and the “P10–P20” mosaics, respectively. The reason there are two mosaics rather than one is because the program was only partially completed and not all data were taken.

We imaged these data using natural weighting (Briggs weighting with robust = 2) to optimize source signal-to-noise ratios. The synthesized beam across all observations was broadly consistent with the beam size of the final mosaics: $1''91 \times 1''41$ for P03 and $1''83 \times 1''43$ for P10–P20. This beam size is larger than the characteristic scale of dust in high-redshift galaxies ($\sim 0''6$, e.g., Hodge et al. 2016) but smaller than the minimum anticipated scale of source confusion at
<table><caption>**Table 1**
MORA Survey Observed Scheduling Blocks, Weather Conditions, and Noise Characteristics</caption><thead><tr><th>Position</th><th>Tuning (GHz)</th><th>R.A.</th><th>PWV (mm)</th><th>On-source Time (min)</th><th>rms (μJy beam<sup>-1</sup>)</th><th>Synth. Beam Size</th></tr></thead><tbody><tr><td>P03</td><td>139</td><td>10:00:43.83</td><td>5.87</td><td>44.00</td><td>62.8</td><td>1"85 × 1"46 79°</td></tr><tr><td>P03</td><td>147</td><td>10:00:42.83</td><td>4.92</td><td>48.82</td><td>62.8</td><td>2"11 × 1"32 72°</td></tr><tr><td>P10</td><td>139</td><td>10:00:29.19</td><td>6.19</td><td>43.98</td><td>98.3</td><td>1"77 × 1"39 68°</td></tr><tr><td>P11</td><td>147</td><td>10:00:25.96</td><td>5.05</td><td>48.75</td><td>83.5</td><td>1"51 × 1"45 30°</td></tr><tr><td>P12</td><td>147</td><td>10:00:23.73</td><td>4.91</td><td>48.78</td><td>88.0</td><td>1"74 × 1"40 54°</td></tr><tr><td>P13</td><td>147</td><td>10:00:21.50</td><td>5.68</td><td>49.75</td><td>87.7</td><td>1"56 × 1"36 86°</td></tr><tr><td>P14</td><td>147</td><td>10:00:19.26</td><td>5.99</td><td>49.77</td><td>90.9</td><td>1"68 × 1"33 67°</td></tr><tr><td>P15</td><td>147</td><td>10:00:17.03</td><td>2.25</td><td>48.78</td><td>68.5</td><td>2"27 × 1"44 72°</td></tr><tr><td>P16</td><td>147</td><td>10:00:14.80</td><td>4.72</td><td>48.77</td><td>112.1</td><td>2"05 × 1"30 60°</td></tr><tr><td>P17</td><td>147</td><td>10:00:12.57</td><td>3.66</td><td>48.77</td><td>69.2</td><td>1"94 × 1"44 68°</td></tr><tr><td>P18</td><td>139</td><td>10:00:10.34</td><td>5.98</td><td>43.98</td><td>60.4</td><td>2"03 × 1"42 68°</td></tr><tr><td>P18</td><td>147</td><td>10:00:10.34</td><td>2.98</td><td>48.75</td><td>60.4</td><td>2"44 × 1"65 61°</td></tr><tr><td>P20</td><td>147</td><td>10:00:05.87</td><td>5.31</td><td>48.83</td><td>91.1</td><td>1"82 × 1"36 77°</td></tr><tr><td>P03 MOSAIC</td><td>...</td><td>...</td><td>...</td><td>92.82</td><td>62.8</td><td>1"91 × 1"41 83°</td></tr><tr><td>P10-P20 MOSAIC</td><td>...</td><td>...</td><td>...</td><td>528.91</td><td>varies</td><td>1"83 × 1"43 88°</td></tr></tbody></table>

**Note.** The range of declinations of pointing centers in each scheduling block is uniform across all SBs and is +02:11:04.16 to +02:33:55.82. The positions ("PXX") correspond to distinct RAs of the mosaic. Two positions were observed at both frequencies: P03 and P18, while all other positions were only observed with one of the two tunings. Only P10–P20 are spatially adjacent such that they have been stitched together in one contiguous mosaic, and the P03 data is a separate mosaic of its own. This table states the final continuum rms achieved in each position of the final mosaic product, yet the synthesized beam is measured for a representative subsample of pointings in each individual data set. The last two rows state the resulting synthesized beams and total on-source time spent for the two end-product mosaics: the P03 mosaic and the P10–P20 mosaic.

2 mm (>20"; Staguhn et al. 2014).<fnref n="25" /> The synthesized beam is ideally matched to our science goals allowing for the detection of unresolved point sources; therefore, no tapering or alternative data weighting was needed.

Several channels covering a 140 MHz wide frequency range centered on an atmospheric absorption feature at 142.2 GHz were flagged for removal in the Tune147 data sets, while no channels were flagged in the Tune139 data sets. We determined that the channel flagging in the Tune147 data sets improved the rms depth of the map by 3% on average (see also Zavala et al. 2021).

To analyze the computational time required to produce the full mosaic, we tested time binning the data by 5, 10, and 30 s. Time binning did substantially speed up the process of combining the visibilities with tclean and the resulting mosaic images were consistent with one another. For our final analysis, we use the 10 s time averaged maps.

Our final mosaics<fnref n="26" /> cover 184 arcmin² of the proposed 230 arcmin². Most of the area is in the P10–P20 mosaic (156 arcmin²) while the remaining 28 arcmin² is in the P03 mosaic. Of the full area, 101 arcmin² is covered at or below the proposed depth of 90 µJy beam⁻¹ (with the deepest part of the map reaching 60 µJy beam⁻¹). See Z21 for more complete information on the noise characteristics of the maps.

Figure 2 shows signal-to-noise ratio maps and rms maps of both mosaics. For context, we have overlaid (in contours) the S2COSMOS SCUBA-2 850µm signal-to-noise ratio map from Simpson et al. (2019).

### 3. The 2 mm Selected Sample

Extensive tests of large extragalactic mosaics from ALMA (Umehata et al. 2015; Aravena et al. 2016; Dunlop et al. 2017; Hatsukade et al. 2016; Walter et al. 2016; Franco et al. 2018, 2020; González-López et al. 2019, 2020) show that ALMA deep fields exhibit Gaussian noise in the absence of bright (S/N > 10) sources. This is demonstrated for the MORA mosaics in our accompanying paper, Z21, which analyzed the number counts and noise characteristics of this data set in depth. In Casey et al. (2018b) we simulated completeness and contamination rates for mock ALMA data sets using the assumption that the noise is Gaussian. As argued in Section 3.1 of Casey et al. (2018b), the measured contamination rates and completeness of simulated ALMA sources do not depend on the wavelength or underlying number density of sources in the map because they are not confusion limited.

Z21 tested that confusion is, indeed, not a concern for the MORA mosaics by masking significantly detected sources (of which there are few across the large mosaic) and injecting fake sources through Monte Carlo trials throughout the rest of the map, measuring both completeness and contamination rates as a function of signal-to-noise ratio and 2 mm flux density. Z21 completed the same procedure for fake maps that have Gaussian noise and the same heterogeneous rms characteristics of our data and found identical rates of contamination and source completeness. Furthermore, we find that sources are recovered at the same rates as simulated in Figure 6 of Casey et al. (2018b), despite the differences in simulated wavelength and the different beam size of observations.

Note that the synthesized beam size of our observations, ~1"9 × 1"4, exceeds the expected FWHM of obscured emission for galaxies at $z > 1$ (~1–5 kpc FWHM, corresponding to scales <0.6, e.g., Simpson et al. 2015; Hodge et al. 2016; Fujimoto et al. 2017); therefore we do not expect any sources to be resolved out beyond the scale of one synthesized beam. Thus the treatment of these sources as point sources rather than sources with extended

<footnote n="25">²⁵ While source confusion is often discussed in the context of single-dish (sub) millimeter maps, it should be noted that the density of 2 mm sources on the sky is much lower than at 850 µm or shorter wavelengths, and therefore confusion would only set in for very low resolution (~1' beam), deep (<1 mJy rms) 2 mm maps.</footnote>

<footnote n="26">²⁶ Mosaics and Measurement Sets available to download at www.as.utexas.edu/čmcasey/downloads/mora.html.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Signal-to-noise ratio (S/N; left) and rms (right) maps of the two MORA mosaics: P03 on the left, and P10–20 on the right. Sources detected at >5σ significance are enclosed in boxes and numbered, corresponding to the sources listed in Table 3 in decreasing order of S/N. The SCUBA-2 850 µm signal-to-noise ratio map from Simpson et al. (2019) is shown in orange contours, denoting 3.5σ, 5σ, and 6.5σ significance. Sources detected in the 2 mm map between 4 < σ < 5 are shown as small diamonds. The rms maps are overlaid with contours beginning at 60 µJy beam⁻¹ and increasing in 30 µJy beam⁻¹ steps; the distribution of rms depths across the full map is shown in Figure 2 of the accompanying MORA Survey paper, Z21. Note that irregularities in the rms map and the discontinuity between the two mosaics are due to the lack of completion of the MORA program.</figcaption></figure>

emission is appropriate and flux densities are measured at the sources' peak. This contrasts with the recent GOODS-ALMA survey covering 69 arcmin² at 1.1 mm, whose synthesized beam size of 0.′′6 is similar to the expected size of sources; this led to incompleteness with regard to spatially extended sources in that work and the need to taper observations to recover extended emission (Franco et al. 2018, 2020).

To further diagnose source contamination, we search the MORA maps for significant negative peaks and find 106 sources below −4σ significance (i.e., strong negative peaks), 15 lower than −4.5σ, and 2 lower than −5σ significance. With ~3 × 10⁵ beams in the

MORA maps, we expect $48_{-5}^{+6}$, $7 \pm 3$, and $1 \pm 1$ negative sources to arise at these respective significances (>4, >4.5, and >5σ). The uncertainties on the expected number of negative noise peaks are determined by generating several fake noise maps with the same noise characteristics as the existing data. The number of negative sources found in the map skews somewhat higher than expected, despite the consistencies of the maps' noise characteristics with modeled Gaussian noise. Z21 demonstrates that the noise characteristics of the MORA mosaics are indeed Gaussian. We do not suspect that the atypical number of negative sources is from the sidelobes of nearby bright sources. Instead the excess of noise
peaks could be due to slight imperfections in modeling the noise. Z21 highlights that the most significant negative detection in the map is found at −6σ, which only has a ≲0.5% chance of being generated in a map of this size due to noise. Whether or not this negative detection is of genuine astronomical origin (i.e., potentially a decrement in the Cosmic Microwave Background (CMB) caused by inverse Compton scattering) is briefly discussed in Z21, but requires further observational follow-up to refute or confirm.

Our tests are consistent with our previous findings in Casey et al. (2018b): sources identified at >5σ significance have little to no contamination with false noise peaks; only 1 ± 1 false source is expected across both MORA mosaics above this threshold from Gaussian noise. Contamination rises sharply at 4.5 < σ < 5 significance, consisting of ~40% false noise peaks, and noise peaks come to dominate the sample detected between 4 < σ < 4.5 significance. In this paper, we present the robust sample of >5σ detected sources and then proceed to analyze the marginal 4 < σ < 5 sample in conjunction with prior identification at other complementary wavelengths. Thirteen >5σ sources are identified, one of which is thought to be false, translating to a >5σ purity of 12/13 ≈ 92% (a figure that would be higher if real sources were more common).

### 3.1. COSMOS Ancillary Data

We make extensive use of the rich ancillary data available in the COSMOS field including the two most recent generations of the optical/near-infrared photometric catalog, COSMOS2015 (Laigle et al. 2016) and COSMOS2020 (Weaver et al. 2021). There are over 30 bands of optical/near-infrared (OIR) imaging that make up these photometric databases, from deep broadband coverage (with depths of 26–28 magnitudes for 5σ point sources) to intermediate and narrowband imaging campaigns (with depths of 25–26 magnitudes).

COSMOS also contains a wealth of multiwavelength data, from the X-ray (Civano et al. 2012, 2016) through the radio (Schinnerer et al. 2007; Smolčić et al. 2017). We include details of detections at these other wave bands where relevant. Particularly important for this work are the other millimeter and submillimeter data sets in the field, including SCUBA-2 at 850 µm (Casey et al. 2013; Geach et al. 2017; Simpson et al. 2019), 450 µm (Casey et al. 2013; Roseboom et al. 2013)<fnref n="27" />, Herschel PACS and SPIRE at 100–500 µm (Lutz et al. 2011; Oliver et al. 2012), AzTEC at 1.1 mm (Scott et al. 2008; Aretxaga et al. 2011), and Spitzer at 24 µm (Le Floc'h et al. 2005). In addition, a variety of ALMA archival data sets have built up in the field—see Liu et al. (2018) for details on the A3COSMOS project—at a range in observed frequencies, with the most common tunings in band 6 (1.2 mm) and band 7 (870 µm). We make use of these data to refine the spectral energy distribution fits for galaxies detected in the MORA maps.

### 3.2. Photometric Redshift Fitting

Photometric redshifts are fit to this existing photometry using the LEPHARE<fnref n="28" /> package and Bruzual & Charlot (2003) templates for star-forming and quiescent galaxies, as in Ilbert et al. (2013). Where OIR photometric redshifts exist for MORA-detected sources, we provide the most recent estimates

from COSMOS2020 when available, or in some instances, we use photometric redshifts from elsewhere in the literature. Several photometric redshift catalogs are presented in the COSMOS2020 compilation; we make use of the CLASSIC SExtractor<fnref n="29" /> photometry and LEPHARE photometric redshifts, as we find they minimize χ² for the redshift estimates. Note that the redshifts overall are indistinguishable from one another for this sample and would not change our results. See Weaver et al. (2021) for more detail on the differences between these catalogs.

For each source in our sample, we also employ the MMPZ technique<fnref n="30" /> for far-infrared (FIR) and millimeter photometric redshift fitting described in Casey (2020) as an independent check on redshift constraints at other wavelengths. The MMPZ technique uses the aggregate FIR/millimeter flux density measurements for a source to derive a probability density distribution in redshift. Rather than basing the redshift fit on a single long-wavelength template, MMPZ presumes the galaxy is most likely to lie on the correlation between galaxies' IR luminosities and their rest-frame peak wavelengths, i.e., in the $L_{IR}-\lambda_{peak}$ plane (as shown in Lee et al. 2013; Strandet et al. 2016; Casey et al. 2018a), and the algorithm determines the redshift range over which the galaxy's SED is likely to be most consistent with that empirical relationship.

### 3.3. The 5σ Subsample

Thirteen sources are identified in our maps above 5σ significance. They are numbered in order of decreasing S/N at 2 mm; their physical characteristics are given in Table 2; and their basic detection characteristics are listed in Table 3. Twelve of the 13 have been detected at other wavelengths and reported in the literature in various forms; in particular, those 12 sources are also detected with SCUBA-2 at 850 µm (Simpson et al. 2019) and with Spitzer IRAC at 3.6 µm. Two of the sources in the sample (sources 5 and 9) have particularly uncertain redshifts and are undetected in deep near-infrared H-band imaging from the CANDELS-COSMOS survey (Koekemoer et al. 2011); these two sources and their properties are described in greater detail in the accompanying paper, Manning et al. (2021).

Here, we provide a brief summary of each of the >5σ sources and the extent to which they have already been characterized in the literature. Table 3 includes additional photometric measurements of each source's integrated flux density at all available wavelengths. Figure 3 shows multiwavelength cutouts of all 13 sources from the optical through the radio.

#### 3.3.1. MORA-0 (850.04 or MAMBO-1)

The highest signal-to-noise ratio source in the MORA maps is detected at nearly 8σ, and has been previously identified as a DSFG through detection at 1.2 mm (MAMBO-1 in Bertoldi et al. 2007), at 1.1 mm (AzTECC7 in Aretxaga et al. 2011), and at 850 µm (850.04 in Casey et al. 2013). The source is only marginally detected with Herschel SPIRE and SCUBA-2 at 450 µm (Oliver et al. 2012; Casey et al. 2013). To date, it does not have a reliable spectroscopic confirmation despite being targeted repeatedly in near-infrared campaigns; the H-band MOSFIRE spectrum presented in Casey et al. (2017) is spatially offset from the ALMA source by 1". The closer (and much fainter) OIR counterpart, only 0.3 offset from the ALMA centroid, has an

<footnote n="27">²⁷ Though deeper 450 µm data exist in a subset of the MORA field from Wang et al. (2017), these data are not publicly accessible and are therefore not included in our analysis.</footnote>

<footnote n="28">²⁸ Software: LEPHARE (Arnouts et al. 2002; Ilbert et al. 2006).</footnote>

<footnote n="29">²⁹ SExtractor (Bertin & Arnouts 1996).</footnote>

<footnote n="30">³⁰ Software: MMPZ (Casey 2020).</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Multiwavelength cutouts (8'' × 8'', where north is up and east is left) of the 13 sources identified at >5σ significance at 2 mm. Cutouts are from HST ACS/F814W (Koekemoer et al. 2007; Scoville et al. 2007), HST WFC3/F125W and WFC3/F160W (Grogin et al. 2011; Koekemoer et al. 2011), Ultra-VISTA Ks band (Laigle et al. 2016), Spitzer 3.6 µm and 4.5 µm (Ashby et al. 2015), 2 mm (this work), and VLA 3 GHz (Smolčić et al. 2017). The sources that sit outside of the CANDELS WFC3 coverage area have J- and H-band cutouts from Ultra-VISTA shown instead. The 5σ contour from the 2 mm imaging is overlaid in each panel for reference. Some sources that sit adjacent to the 2 mm detected galaxies are labeled with redshifts with more details given in the text.</figcaption></figure>

OIR-based photometric redshift from the Laigle et al. (2016) COSMOS2015 catalog of $z_{phot} = 3.31_{-0.81}^{+0.76}$; the source is missing from the COSMOS2020 photometric catalogs for no obvious reason other than its marginal detection near the detection limit of the near-infrared catalogs.

MORA-0 has an FIR photometric redshift from Brisbin et al. (2017) of $z_{FIR} = 4.4 \pm 1.0$. Using the MMPZ FIR/millimeter photometric redshift technique, we derive a millimeter-based redshift of $z_{mm} = 3.4_{-0.6}^{+0.6}$; both FIR/millimeter photometric redshifts are higher than, though consistent with, the OIR photometric redshift. The source is detected in both the 1.4 GHz and 3.0 GHz radio maps. There are several ALMA programs that have obtained continuum data on MORA-0 in band 6 (1287 and 1250 µm) and in band 7 (870 µm). While two radio galaxies were originally thought to be associated with this source (Bertoldi et al. 2007), these are now thought to sit at

different redshifts: one with $z_{spec} = 1.436$ (reported in Casey et al. 2017) and MORA-0 with a higher photometric redshift. The source at $z = 1.436$ is not detected in any of the ALMA data, while MORA-0 is detected in all of bands 4, 6, and 7 where data exist. No millimeter spectroscopy exists for MORA-0.

The lack of detection in the COSMOS2020 photometric redshift catalog, yet detection in COSMOS2015, casts some doubt on the quality of OIR constraints; however, the broad consistency of the COSMOS2015 constraint with the MMPZ redshift is reassuring. We adopt the COSMOS2015 photometric redshift throughout the rest of the text.

### 3.3.2. MORA-1 (AzTEC-5)

The second most significant source, MORA-1, has been studied under many names, the most widely used of which is
<figure><img src="image_8.png" /><figcaption>Figure 3. (Continued.)</figcaption></figure>

AzTEC-5 for its initial detection in Scott et al. (2008) at 1.1 mm. Magnelli et al. (2019) includes a nice discussion of its known characteristics and redshift constraints, which we summarize here for completeness. AzTEC-5 is Herschel SPIRE, SCUBA-2 (450 and 850 µm) and ALMA detected, with continuum measurements in bands 6 and 7.

The literature alludes to a spectroscopic identification for AzTEC-5 of $z = 3.791$ based on $\text{Ly}\alpha$ emission in a DEIMOS spectrum (Capak et al. 2010; Smolčić et al. 2012); however, that solution was revealed to be inaccurate by subsequent ALMA follow-up that failed to detect emission lines to corroborate the redshift.

Gómez-Guijarro et al. (2018) identify four components of AzTEC-5 in deep near-infrared imaging, which they dubbed AzTEC5-1, 5–2, 5–3, and 5–4. MORA-1 is coincident with the position of AzTEC5-1, which is the only component lacking direct redshift constraints. Gómez-Guijarro et al. (2018) do not fit a photometric redshift for AzTEC5-1; it is absent from the COSMOS2020 photometric catalogs. The photometric redshifts for the other components are $z = 3.63_{-0.15}^{+0.14}$ for AzTEC5-2, $z = 4.02 \pm 0.08$ for AzTEC5-3, and $z = 3.66_{-0.43}^{+0.40}$ for AzTEC5-4. The redshift constraints on existing components is shown in the third panel (second row) of Figure 3. In addition to AzTEC5-1, AzTEC5-2 also appears to have associated 870 µm emission,

separated from AzTEC5-1 by $0.7''$; ATEC5-2 is not detected at $>5\sigma$ significance in our mosaics. Our millimeter-derived photometric redshift for the aggregate photometry of all components of AzTEC-5 gives $z_{\text{mm}} = 2.6_{-0.8}^{+2.2}$, while AzTEC5-1 alone has $z_{\text{mm}} = 4.8_{-2.1}^{+3.7}$. Due to the significant uncertainties on both, they are consistent with the photometric redshift constraints of the other three components from Gómez-Guijarro et al. (2018).

Having detected AzTEC-5 at 2 mm with GISMO, albeit with far worse spatial resolution than the MORA map, Magnelli et al. (2019) argue that the redshift of the entire system is likely $z \sim 3.6$, with substantial obscuration in AzTEC5-1 making it difficult to spectroscopically confirm. Indeed, the close spatial separation $<1''$ between the two sources (AzTEC5-1 and AzTEC5-2) would support this claim, as the likelihood of identifying two 870 µm sources with $>3$ mJy at different redshifts separated by $<1''$ is exceedingly low (0.02% based on 850 µm number counts; Geach et al. 2017; Simpson et al. 2019). While spectroscopic confirmation is needed for MORA-1$^{31}$, we determine that it is best, for the purposes of this work, to adopt the median redshift

<footnote>$^{31}$ We note that alternate band 4 ALMA spectral line observations of MORA-1 exist under the source name “GalD” in program 2018.1.01824.S; however, the tuning would only capture CO(6–5) from redshifts $z = 3.65–3.77$ and $z = 4.06–4.20$. No CO line is detected in the data set to a depth of $1\sigma \sim 0.4$ mJy beam$^{-1}$ in 100 km s$^{-1}$ channels.</footnote>
of the three other components of the AzTEC-5 system given in Gómez-Guijarro et al. (2018), which is $z_{phot} = 3.78^{+0.5}_{-0.32}$.

### 3.3.3. MORA-2 (COS850.0020)

The third source in the MORA sample, MORA-2, is detected at a signal-to-noise ratio of 7.6 at 2 mm. It has no match in the earliest SCUBA-2 surveys (Roseboom et al. 2012; Casey et al. 2013) but is detected in the wider-field surveys of COSMOS from Geach et al. (2017) and Simpson et al. (2019). Both ALMA bands 6 and 7 data exist for MORA-2, which was selected both as a submillimeter source and a dropout source in the ZFOURGE survey (Spitler et al. 2012). Its OIR photometry gives a photometric redshift from COSMOS2020 of $z_{phot} = 3.36^{+0.60}_{-0.28}$. There are no spectroscopic constraints for MORA-2; our millimeter-derived photometric redshift for the source is $z_{mm} = 2.8^{+0.6}_{-0.6}$, which is consistent with the OIR photometric redshift.

### 3.3.4. MORA-3 (AzTEC-2 or 450.03)

The fourth source, MORA-3, is more widely known as AzTEC-2 (Younger et al. 2007) detected with the SMA, and at 450 µm and 850 µm by SCUBA-2 (where it is known as 450.03; Casey et al. 2013). It was also previously detected at 2 mm as GISMO-C2 (Magnelli et al. 2019). It is detected by Herschel SPIRE and has ALMA continuum data in bands 6 and 7. The interferometric data reveal two millimeter sources, the primary (coincident with the 2 mm emission) is AzTEC2-A and the secondary fainter source (detected at $2.6\sigma$ significance in the 2 mm map) is AzTEC2-B. While a spectroscopic identification based on a possible OIR counterpart existed at $z = 1.123$ (Smolčić et al. 2012; Casey et al. 2017) 1" to the south of the primary source A, that redshift has since been shown to be associated with a foreground galaxy.

The spectroscopic redshift of AzTEC-2 is now confirmed through the detection of [C II] and CO(5-4) at $z = 4.63$ by Jiménez-Andrade et al. (2020); this was independently confirmed in Simpson et al. (2020). AzTEC2-A has a spectroscopic redshift of $z = 4.625$, while AzTEC2-B is at $z = 4.633$. Our millimeter-derived photometric redshift for this source is consistent with its spectroscopic identification, $z_{mm} = 3.3^{+1.0}_{-0.8}$. Both components of AzTEC-2 are undetected in existing HST deep imaging, and the Spitzer imaging is highly confused with two foreground galaxies: the $z = 1.123$ galaxy 1" to the south of AzTEC2-A and an elliptical galaxy at $z = 0.34$ at 1" south of the $z = 1.123$ system. Both components of AzTEC-2 are detected at 3 GHz in radio continuum. See Jiménez-Andrade et al. (2020) for a more thorough discussion of this source.

Given the proximity of the foreground elliptical galaxy at $z = 0.34$, Jiménez-Andrade et al. (2020) estimate that the luminosity of AzTEC2-A, or MORA-3, is gravitationally lensed by a magnification factor of $\mu = 1.5$. We scale physical quantities proportional to luminosity by this magnification factor for MORA-3 for the rest of this paper. Note that this source sits in a large-scale overdensity at $z = 4.6$ identified and described in Mitsuhashi et al. (2021), further corroborating the conjecture that high star formation rate galaxies at $z > 4$ are highly clustered and good signposts for the most massive overdense structures in the universe (Casey 2016; Chiang et al. 2017).

### 3.3.5. MORA-4 (MAMBO-9)

The fifth source of the sample, MORA-4, is known as MAMBO-9. It was originally detected at 1.2 mm by the MAMBO instrument (Bertoldi et al. 2007), and subsequently detected at 1.1 mm (as AzTEC/C148; Aretxaga et al. 2011) and by SCUBA-2 at 850 µm (as 850.43 and COS.0059 in Casey et al. 2013 and Geach et al. 2017, respectively). Several teams identified the source as potentially high-z based on being undetected in the Herschel SPIRE bands; Jin et al. (2019) initially report a spectroscopic identification of $z = 5.850$ based on a low signal-to-noise ratio 3 mm spectral scan, confirmed by detection of $^{12}$CO($J=6 \to 5$) and p-H$_2$O($2_{1,1} \to 2_{0,2}$) in Casey et al. (2019). MAMBO-9 is composed of two galaxies separated by 6 kpc (1") and both confirmed at $z = 5.850$. Our millimeter-derived photometric redshift for MAMBO-9 is $z_{mm} = 5.1^{+0.6}_{-0.8}$, only in slight tension with the measured spectroscopic redshift. MAMBO-9 is the most distant un lensed DSFG found to date, and we refer the reader to Casey et al. (2019) for a more thorough characterization of the MAMBO-9 system.

### 3.3.6. MORA-5 (850.13 or AzTEC C114)

The sixth source, MORA-5, has been identified at both 850 µm and 1.1 mm (Aretxaga et al. 2011 and Casey et al. 2013, named 850.13 and AzTEC C114, respectively). There is no spectroscopic redshift for MORA-5, and there is no OIR-based photometric redshift from either the COSMOS2015 or COSMOS2020 catalogs due to a lack of a counterpart in the near-infrared. Brisbin et al. (2017) offer an FIR photometric redshift of $z_{FIR} = 5.3 \pm 3.2$, while noting that a radio-FIR-based photometric redshift is consistently lower, based on the source's detection at both 1.4 GHz and 3 GHz (respective radio-FIR photometric redshifts of $z_{1.4GHz} = 2.9^{+4.2}_{-0.4}$ and $z_{3.0GHz} = 1.9^{+1.1}_{-0.3}$). Our millimeter-derived photometric redshift is $z_{mm} = 3.4^{+1.1}_{-0.9}$. This source is undetected at all wavelengths shortward of 3.6 µm, including deep CANDELS H-band and Ultra-VISTA K-band imaging. Our accompanying paper, Manning et al. (2021), presents a more detailed analysis of this source and calculates a hybrid photometric redshift estimate for the source of $z_{phot} = 4.3^{+1.5}_{-1.3}$ combining the MMPZ redshift with direct extraction and refitting of OIR constraints using both EAZY$^{32}$ and MAGPHYS$^{33}$ approaches (see Manning et al. 2021, for more details). We adopt the Manning et al. hybrid photometric redshift for the rest of this paper.

### 3.3.7. MORA-6 (450.00)

The seventh source, MORA-6, was detected as the brightest 450 µm source, named 450.00, in the 394 arcmin$^2$ map in Casey et al. (2013). This source is also detected at 850 µm with SCUBA-2 and has ALMA continuum follow-up in both bands 6 and 7. It lacks a spectroscopic redshift but does have a fairly well-constrained OIR-based photometric redshift of $z_{phot} = 3.34^{+0.13}_{-0.12}$ from COSMOS2020. Our millimeter-derived photometric redshift is $z_{mm} = 2.5^{+1.6}_{-0.7}$, consistent with the OIR photometric redshift.

<footnote>$^{32}$ Software: EAZY (Brammer et al. 2008).</footnote>

<footnote>$^{33}$ MAGPHYS (da Cunha et al. 2008).</footnote>
### 3.3.8. MORA-7 (850.53)

The eighth source, MORA-7, has been previously identified at 850 µm in both Casey et al. (2013) as 850.53 (lacking a corresponding detection at 450 µm) and Geach et al. (2017) as COS850.0035. This source lacks a spectroscopic redshift, but is detected with Spitzer and Ultra-VISTA, rendering an OIR photometric redshift estimate of $z_{phot} = 2.85_{-0.33}^{+0.24}$. The source has dust continuum observations from ALMA in both bands 6 and 7. Our millimeter photometric redshift is $z_{mm} = 2.3_{-0.8}^{+2.7}$, consistent with the OIR photometric redshift.

### 3.3.9. MORA-8 (850.09)

The ninth source of our sample, MORA-8, has been previously identified at both 450 and 850 µm in Casey et al. (2013), wherein it was referred to as 850.09, and in Geach et al. (2017), where it was named COS850.0016. The source lacks spectroscopic confirmation, but the OIR photometric redshift from COSMOS2020 is fairly well constrained as $z_{phot} = 2.29_{-0.08}^{+0.12}$. This source also has ALMA band 7 continuum data, from which our millimeter-derived redshift is $z_{mm} = 3.0_{-1.0}^{+2.2}$. Both ALMA 2 mm and 870 µm sources are well aligned with the OIR counterpart.

### 3.3.10. MORA-9

The tenth source of the sample, MORA-9, has been detected at 850 µm in Simpson et al. (2019) at a $4.4\sigma$ significance; it had not been detected in the earlier compilation of Geach et al. (2017). Aside from its detection at 2 mm and 850 µm, the source is also detected at 3.6 µm and 4.5 µm from Spitzer and at low significance in Ultra-VISTA K-band imaging. This source is undetected in all other data sets, but its K-band counterpart renders it present in the COSMOS2020 photometric catalogs with an OIR photometric constraint of $z_{OIR} = 4.57_{-1.22}^{+1.33}$. It is one of only two $>5\sigma$ sources not already surveyed by ALMA at other wavelengths (the other being MORA-12). The ratio of 850 µm to 2 mm flux density is highly suggestive of a high-z solution; we derive a millimeter photometric redshift of $z_{mm} = 5.5_{-0.8}^{+0.6}$ for this source. Manning et al. (2021) describe this source's characteristics and derive a hybrid photometric redshift for MORA-9 of $z_{phot} = 4.3_{-1.0}^{+1.3}$, which includes this millimeter photometric redshift along with OIR constraints using photometric redshift fitting techniques EAZY and MAGPHYS. We adopt the Manning et al. hybrid photometric redshift for the rest of this paper.

### 3.3.11. MORA-10 (450.09)

The eleventh source of our sample, MORA-10, is well characterized as a known DSFG detected at both 450 and 850 µm from SCUBA-2 (named 450.09; Casey et al. 2013), with a near-infrared spectroscopic redshift of $z = 2.472$ as reported in Casey et al. (2015). The source is one of many DSFGs in the COSMOS field that sits in a protocluster environment at $z \sim 2.5$, the same structure that motivated the dual spectral tunings for the MORA program (Casey et al. 2015; Chiang et al. 2015; Diener et al. 2015; Cucciati et al. 2018). ALMA data exist for 450.09 in both bands 7 and 3, where the band 3 data confirm its spectroscopic redshift via detection of CO(3–2). The millimeter photometric redshift for the source is $z_{mm} = 4.5_{-1.0}^{+1.4}$, which is $2\sigma$ discrepant with the measured spectroscopic redshift (this discrepancy originates from the galaxy appearing to be a bit colder than the average SED).

MORA-10 is unique among the MORA Survey sample for being particularly luminous in its radio continuum ($S_{1.4 GHz} = 5.7$ mJy and $S_{3.0 GHz} = 3.2$ mJy) with a rest-frame radio luminosity of $L_{178 MHz} = 1.9 \times 10^{27}$ W Hz$^{-1}$, nearly bright enough to fit the local Fanaroff–Riley class II radio-loud active galactic nucleus (AGN; Fanaroff & Riley 1974) definition at its redshift. Fitting the existing radio continuum measurements to a power law, we derive a synchrotron slope of $\alpha = 0.85$; extended to the observed 2 mm band data, we estimate a total synchrotron contribution of $S_{2 mm,synch} = 159 \pm 60$ µJy toward the total observed $S_{2 mm} = 405 \pm 78$ µJy. This implies that $\approx 39 \pm 17\%$ of the total measured 2 mm flux density is due to synchrotron processes. While this clearly does not dominate the total flux density, an absence of synchrotron emission in this source would have rendered it below the $5\sigma$ detection limit of our sample. Even if the synchrotron slope varied somewhat, the source would have not been detected at a high significance from its dust emission alone were its synchrotron component subtracted from the total 2 mm flux density. Because the primary goal of this work is to identify thermal dust emission at 2 mm, we will exclude MORA-10 from analysis of population statistics, such as its contribution to the 2 mm selected galaxy redshift distribution and SFRD. For MORA-10's far-infrared/millimeter SED fit, we adjust the 2 mm flux density to only account for the estimated dust continuum component, removing the synchrotron component, as we are only fitting the dust SEDs in this paper.

### 3.3.12. MORA-11 (AzTECC76 or COS.0194)

The penultimate source of the sample, detected at $5.2\sigma$ significance, is MORA-11. It is detected at 850 µm in both Geach et al. (2017), where it was named COS.0194, and in Simpson et al. (2019). It is also detected at 1.1 mm from Aretxaga et al. (2011) and in the Herschel SPIRE bands at 250–500 µm. Broadly detected in the near-infrared through radio, MORA-11 has a reported medium-band survey redshift of $z = 3.17 \pm 0.12$ from the NEWFIRM Medium Band Survey (Whitaker et al. 2011); the source has additional Keck-NIRSPEC K-band spectroscopic observations from Marsan et al. (2017), though no emission lines were detected. This redshift is in agreement with its OIR photometric redshift from COSMOS2020 of $z_{phot} = 3.00_{-0.13}^{+0.10}$. We measure a millimeter photometric redshift of $z_{mm} = 4.3_{-1.2}^{+4.9}$ for this source. We adopt the medium-band photometric redshift for MORA-11 for the rest of our analysis in lieu of the COSMOS2020 photometric redshift due to the improved precision offered by the spectrophotometric analysis of Marsan et al. (2017). Note that spectral analysis of our own MORA mosaic has tentatively identified an emission line at 140.85 GHz, which could be CO(5–4) at $z = 3.091$, consistent with the adopted redshift of $z = 3.17 \pm 0.12$; further analysis of this line identification is in progress (Mitsuhashi et al. 2021, in preparation).

### 3.3.13. MORA-12

The last source in our sample, MORA-12, is detected at $5.02\sigma$ significance; unlike the rest of the $>5\sigma$ sample, it has no corresponding detection at 850 µm. While there is a 850 µm source (with strong 24 µm emission, named 850.77 in Casey et al. 2013) $9''$6 away from the ALMA 2 mm position, we think it is unlikely that the two are associated. MORA-12 is also not
detected in the Herschel COSMOS maps, as well as the
GISMO map from Magnelli et al. (2019). This source is the
only source in the >5σ sample to lack detection in Spitzer
IRAC at 3.6 µm or 4.5 µm. There is no additional ALMA data
at any other wavelength at the position of MORA-12. With
only one detection at one wavelength, we are unable to derive a
millimeter photometric redshift for this source. Because it has
not been detected at any other wavelength and it sits on the
boundary of the P10–P20 mosaic, and given our expected
contamination rate of 1 ± 1 false source detected above the 5σ
threshold, we conclude it is most likely that MORA-12 is not
real and that it is a positive noise peak.

### 3.3.14. Summary Characteristics of >5σ Sample

Of the 13 sources detected at >5σ, we determine that 12 of
them are real 2 mm detected galaxies, while the last and least
significant (MORA-12) is likely a positive noise peak. Of the
12 real sources, all are detected with both Spitzer IRAC and
SCUBA-2 at 850 µm. Despite all being previously detected,
not all sources had been identified as high-z candidates. One of
the 12 sources (MORA-10) is thought to have a substantial
contribution (39% ± 17%) from synchrotron radio emission to
its 2 mm flux density, rendering the sample of purely dust-
selected 2 mm sources with only 11 galaxies.

Three of the 12 sources are spectroscopically confirmed at
z = 2.472 (MORA-10, known as 450.09, with the synchrotron
component), z = 4.625 (MORA-3, known as AzTEC-2), and
z = 5.850 (MORA-4, known as MAMBO-9); of the remaining
nine sources, eight have some form of OIR-based photometric
redshift (Laigle et al. 2016; Marsan et al. 2017; Gómez-Guijarro
et al. 2018; Weaver et al. 2021). The last source (MORA-5) lacks
an OIR counterpart and any redshift constraint; MORA-5 along
with MORA-9, whose OIR phot-z is highly uncertain, have
hybrid photometric redshift fits provided in our accompanying
manuscript (Manning et al. 2021).

Four of the 11 dust-selected sources (36%) are “OIR-dark,”
meaning they lack near-infrared H-band counterparts in deep
WFC3 imaging, which in this case is CANDELS-COSMOS
data reaching a 5σ point-source depth of 27.15 (Koekemoer
et al. 2011). These sources are MORA-3 (known as AzTEC-2,
spectroscopically confirmed at z = 4.63), MORA-4 (known as
MAMBO-9, spectroscopically confirmed at z = 5.85), MORA-
5, and MORA-9. As these represent the potential highest-
redshift subset of the 2 mm selected sample, a more thorough
analysis of them is given in the accompanying paper by
Manning et al. (2021). It appears that MORA-3 (AzTEC-2) is
the only galaxy in the sample that is gravitationally lensed
(μ = 1.5; Jiménez-Andrade et al. 2020).

The majority of the MORA 2 mm selected galaxy sample is
consistent with relatively little AGN activity, or AGN activity
that does not dominate the galaxies’ bolometric luminosities.
We investigate the sample’s AGN content by analyzing X-ray
imaging, radio emission, and mid-infrared emission. None of
the >5σ sample is detected in the deep COSMOS Chandra or
XMM data, though the sensitivity of such X-ray surveys is
rather shallow at z > 3. The galaxies’ radio luminosities
measured at 3 GHz (Smolčić et al. 2017), and in the somewhat
shallower 1.4 GHz data, are in line with expectation for
synchrotron emission generated via star formation processes
instead of AGNs (e.g., Yun et al. 2001; Ivison et al. 2010;
Delhaize et al. 2017). The one source that proves an exception
to this is MORA-10, which appears to be radio loud and whose

2 mm emission is partially dominated by such nonthermal
emission mechanisms. In the mid-infrared, five of 11³⁴ (5/
11 = 45%) are 24 µm detected. Aside from MORA-10, which
clearly has an AGN, three of the other four 24 µm detected
galaxies are the lowest-redshift sources in the sample (MORA-
6, −7, and −8), which, at those redshifts, could be due to either
AGN or star formation via polycyclic aromatic hydrocarbon
(PAH) emission. The last 24 µm detection, MORA-1, is a
blend of several sources at z ~ 3.6–4.0, and the centroid of the
emission is not precisely constrained to the 2 mm source in
question. While AGNs overall seem nondominant in this
sample of 2 mm selected DSFGs, it should be noted that recent
modeling work has shown that AGNs can contribute
substantially to the heating of host galaxy scale dust, even in
IR luminous galaxies without prominent AGNs (McKinney
et al. 2021); such effects are difficult to account for without
constraining observations; thus we do not account for it directly
in this work.

### 3.4. Marginal Sources with <5σ Significance

Below the 5σ threshold, contamination from positive noise
peaks becomes a significant concern. There are 87 sources
identified in the MORA P10–20 mosaic and 11 sources found
in the P03 mosaic at 4 < σ < 5 significance, or 98 in total.
Their positions and characteristics are given in Table 4. We
determine which of these sources are most likely to be real by
testing for coincident identification at other wavelengths, for
example, in the COSMOS2015 catalog, the COSMOS2020
catalog, the 3 GHz radio continuum survey (Smolčić et al.
2017), or at 850 µm from SCUBA-2 (Simpson et al. 2019).

Given the high source density of galaxies in the OIR catalogs,
there is an 8% chance of a random point in our map aligning with
a COSMOS2015 counterpart within 1" (which is a conservative
maximum physical scale on which we see spatial offsets of
obscured and unobscured emission in galaxies, e.g., Biggs &
Ivison 2008, with more characteristic scales of 0"′1–0"′3, as
in Cochrane et al. 2019). The probability of a chance alignment with
a source in COSMOS2020 is slightly higher within 1", ~13%,
given its increased depth at near-infrared wavelengths. Out of a
sample of 98 marginal sources, this would suggest a total of 8 ± 3
sources matched at random to COSMOS2015 and 13 ± 4 sources
matched at random to COSMOS2020. We find a total of 12
matches to COSMOS2015 and 17 matches to COSMOS2020
within the MORA 4 < σ < 5 2 mm sample (11/12 overlap
between COSMOS2015 and COSMOS2020), suggesting that
there are very few real sources (<1/3) in this marginal sample and
none that can be directly identified reliably.

We also tested the correlation between the marginal sample
with 850 µm selected DSFGs observed by SCUBA-2 (Simpson
et al. 2019); we find that there should be a 2.9% chance of a
random alignment between a >3.5σ SCUBA-2 source (within
the SCUBA-2 15" FWHM beam) and a marginal source in our
2 mm catalog. The rate of false positives is lower for 850 µm
counterparts than for OIR counterparts because the sky density
is much lower for 850 µm sources. We find 12 sources that are
spatially coincident with a 850 µm SCUBA-2 source (a total of
12% of the marginal sample), well in excess of the anticipated
~3%. Contrary to our findings using OIR counterparts, this
demonstrates a true excess of real sources within the marginal

<footnote n="34">³⁴ MORA-3, known as AzTEC-2, is too severely blended with foreground galaxies to discern whether or not it is 24 µm luminous.</footnote>
sample. Of these 12 sources, three have secure IRAC 3.6 µm counterparts matched to the positions of the 2 mm emission, corroborating their identification as real sources.

We further investigate radio continuum counterparts for the marginal sample using the 3 GHz radio continuum map in COSMOS presented in Smolčić et al. (2017). Some of the sources are also covered by the deeper COSMOS-XS survey (Algera et al. 2020; van der Vlugt et al. 2021). We find that there is a 0.6% chance of a random alignment between a marginal MORA source and a 3 GHz $>5\sigma$ detected source within 1" of one another. The rate of false positives is lowest for radio counterparts due to their increased rarity on the sky, in addition to their precisely measured positions. We find that 6 MORA 2 mm sources are spatially coincident with a 3 GHz radio continuum source (a total of 6% of the marginal sample), a factor of 10 higher than the anticipated ~0.6%. The strength of this excess is such that the 3 GHz detected subset can be reliably identified as having real 2 mm emission. Nevertheless, the sample is somewhat limited in size to analyze in further detail.

We include an abbreviated table of marginally detected $4 < \sigma < 5$ 2 mm sources in Table 4 for reference and note which sources have counterparts at which wavelengths. However, due to high contamination rates and a similarly high incompleteness rate in this flux density regime, we do not analyze the sample further. Later in Section 5.5 we analyze the 2 mm emission properties of any DSFGs that have been independently observed by ALMA (not a part of MORA) in the field, some of which overlap with this marginal sample.

The lack of reliability of the marginal sample is further verified by analyzing the number of negative peaks in the mosaics between 4–5$\sigma$, of which there are a total of 106. If there were a significant population of real sources embedded in the marginally identified sample, the positive peaks (98) would likely outnumber the negative peaks (106).

## 4. Models of the 2 mm Universe

We make use of several cosmological semi-empirical and empirical models of the 2 mm luminous universe to draw comparisons with the MORA data set. A brief description of each model data set follows.

First, we compare to the Simulated Infrared Dusty Extra-galactic Sky (SIDES) model (Bethermin et al. 2017), which builds galaxies’ SEDs from their stellar masses and star formation rates (assuming a bimodal population of main-sequence galaxies and starbursts), and which updates the 2SFM (2 Star Formation Modes) galaxy evolution model (Béthermin et al. 2012; Sargent et al. 2012) to analyze the impact of clustering on IR map analysis. The 2SFM model posits that there is a bimodal population of star-forming galaxies: those that are on the main sequence and starbursts that have elevated specific star formation rates; the 2SFM model uses this framework to model all galaxies. We make use of the full SIDES 2 deg$^2$ light cone in our analysis, both to sample cosmic variance and to understand possible trends for 2 mm selected galaxies on angular scales larger than the MORA survey.

Second, we compare our results to the SHARK semi-analytic model of galaxy formation (Lagos et al. 2018). By using the SED code PROSPECT (Robotham & Bellstedt 2020)$^{35}$ and the radiative transfer analysis of the EAGLE hydrodynamical simulations of Trayford et al. (2020), the

SHARK model was successfully able to predict the ultraviolet to far-infrared emission of galaxies over a wide range assuming a universal IMF (Lagos et al. 2019).$^{36}$ Lagos et al. (2020) presents detailed predictions for the (sub)millimeter galaxy population, including 2 mm number counts and redshift distributions. We make use of the full 108 deg$^2$ SHARK light cone for our comparisons.

Third, we compare our work to the semi-empirical model for dust continuum emission published by Popping et al. (2020); that work primarily focused on comparison of 850 µm and 1.1 mm number counts and redshift distributions with the ASPECS survey results (González-López et al. 2019, 2020; Aravena et al. 2020). Using the same methodology of Popping et al. (2020), several light cones from the UNIVERSEMACHINE framework (Behroozi et al. 2019), which is grounded in the Bolshoi–Planck dark matter simulation (Klypin et al. 2016; Rodríguez-Puebla et al. 2016), are stitched together to form a 7.7 deg$^2$ light cone. Dark matter halos are populated with galaxies calibrated to observationally constrained relations (in stellar mass and SFR distributions), and dust continuum characteristics are applied using the relations derived in Hayward et al. (2011, 2013) using the SUNRISE (Jonsson 2006) dust radiative transfer code as a function of SFR and dust mass.

Lastly, we compare to the empirical model predictions for the submillimeter sky presented by Casey et al. (2018a, 2018b) and expanded on in Zavala et al. (2018a) and, finally, in Z1. The key difference between the Casey et al. (2018a) models and the cosmological semi-empirical models is the built-in flexibility to test different hypothetical evolving infrared luminosity functions against data. Not tethering this model directly to any cosmological simulation renders it a tool to interpret data that may be discrepant with such simulations. Casey et al. (2018a) use it to present the hypothetical (sub) millimeter sky in two diametrically opposed universes: the “dust poor” universe (model A) in which there is steep evolution from $z \sim 7$ to $z \sim 2$ in the characteristic number density of the IRLF ($\Phi_*$), and the “dust rich” universe (model B), in which the evolution of the characteristic number density is much more shallow. The “dust poor” universe effectively translates to a very minor contribution of intense, dusty starbursts to cosmic star formation beyond $z \gtrsim 5$ (<10%), while they would dominate cosmic star formation at similar epochs in the “dust rich” universe (>90%). These models are built to capture two extremes. Zavala et al. (2018a) use 3 mm number counts from the ALMA archive to argue for a solution between these two extremes. Our accompanying paper, Zavala et al. (2021), presents an update to Zavala et al. (2018a) using MORA 2 mm number counts and updated 3 mm archival counts, as well as deep 1 mm number counts from the ASPECS survey (González-López et al. 2019, 2020; Aravena et al. 2020). The expanded data set used to constrain the model in Zavala et al. (2021) relative to Zavala et al. (2018a) has resulted in a change of the predicted $z \gtrsim 2$ evolution of $\Phi_*$, the number density of DSFGs at early times, from $\Phi_* \propto (1+z)^{-4.2}$ to $\Phi_* \propto (1+z)^{-6.5}$ (see Zavala et al. 2021, for more details).

<footnote>$^{35}$ Software: PROSPECT (Robotham & Bellstedt 2020).</footnote>

<footnote>$^{36}$ See also Lovell et al. (2020) and Hayward et al. (2021) for similar analyses using the SIMBA and ILLUSTRIS-TNG hydrodynamical simulations from Davé et al. (2019) and the radiative transfer code POWDERDAY from Narayanan et al. (2021).</footnote>
<figure><img src="image_9.png" /></figure>

**Figure 4.** Comparison of the redshift probability density distributions derived from: the FIR/millimeter using the MMPZ technique (teal distributions), OIR photometric redshifts (blue distributions), and spectroscopic constraints (black lines). The hashed regions of those same colors indicate the 68th percentile inner confidence interval. Sources are ordered by redshift, from highest (top) to lowest (bottom). The two unconfirmed OIR-dark sources with the most uncertain redshift constraints are MORA-5 and MORA-9, described in more detail in our accompanying manuscript, Manning et al. (2021). The redshifts for these two galaxies are derived jointly using OIR and millimeter constraints, here called hybrid.

# 5. Results

## 5.1. Redshift Constraints and Distribution

Redshift constraints for our sample are heterogeneous, ranging from firm spectroscopic confirmations to limited far-infrared/millimeter photometric constraints. We exclude MORA-10 (which has $z_{spec} = 2.472$) from analysis of the sample's redshift distribution due to the suspected contribution of synchrotron emission to its 2 mm flux density, without which it would not have been detected above $5\sigma$ in our data.

From the most to least reliable, two (2/11 = 18%) have spectroscopic redshifts, seven (7/11 ≈ 64%) have OIR photometric redshifts, and two (2/11 ≈ 18%) have hybrid FIR/millimeter and OIR photometric redshift constraints. It is somewhat interesting to note that the two spectroscopic redshifts are the highest-redshift sources in the sample. Figure 4 shows the redshift constraints for each source (including MORA-10) in order of increasing redshift from bottom to top and how consistent existing constraints are with FIR/millimeter-derived redshifts from the MMPZ tool described in Casey (2020). The millimeter photometric redshifts serve as a sanity check on the tighter constraints given by other methods.

Figure 5 shows the cumulative redshift distribution for the entire sample in two panels. Given the heterogeneous constraints, each source's probability density distribution in $z$ is coadded and shown as a cumulative distribution function (CDF) to make clear the relative fraction of the sample below or above a certain redshift threshold. Accounting for the uncertainties in the individual redshift constraints through Monte Carlo draws from the CDF, we measure the median redshift of the sample as $\langle z_{2\,\text{mm}} \rangle = 3.6_{-0.3}^{+0.4}$. The variance in the redshift CDF is shown in gray and encompasses a 68% ($\pm 1\sigma$) minimum confidence interval. We find that 77% $\pm 11\%$ of the distribution lies at $z > 3$ and 38% $\pm 12\%$ lies at $z > 4$.

The median redshift of 2 mm selected galaxies has been measured twice before in the literature, using the GISMO instrument on the single-dish IRAM 30 m telescope. Staguhn et al. (2014) measure a median redshift of $\langle z_{2\,\text{mm}} \rangle = 2.9 \pm 0.9$ for sources in the GISMO Deep Field in GOODS-N, while Magnelli et al. (2019) measure a median redshift of $z \sim 4$, though only for five sources and four sources, respectively, with $S/N > 4$. Both are in agreement with our measured median redshift of $\langle z_{2\,\text{mm}} \rangle = 3.6_{-0.3}^{+0.4}$.

Figure 5 also shows the predicted cumulative redshift distributions for the MORA data set from several models in the literature, described in Section 4. For all models, we have generated the cumulative redshift distribution by sampling simulations over the $184 \, \text{arcmin}^2$ area of the MORA Survey. Some of the models have been simulated over larger areas (e.g., SHARK, SIDES, and UNIVERSEMACHINE) while the Casey et al. and Zavala et al. empirical models have multiple realizations the same size as the MORA survey. Sources in each model data set then have an rms noise assigned to them following the heterogeneous distribution of rms noise in the MORA maps to best mimic the data. We retain simulated sources that were detected at or above $5\sigma$ significance. Sources in these simulations are assumed to be point sources, as the probability of them being spatially resolved on scales larger than $1'' \approx 10 \, \text{kpc}$ is unlikely at these redshifts. This procedure is repeated 1000 times to constrain the uncertainties of the total number of selected sources and their redshift distribution.

The predictions of the Casey et al. (2018a) models A and B, as well as the Zavala et al. (2018a) model, are shown in the left panel of Figure 5. All are an overestimate of the number of sources in the true data, though the redshift distribution for the model A from Casey et al. (2018a) does agree with our data within uncertainties. The right panel of Figure 5 shows four models, one of which is the updated empirical model from our accompanying paper, Z21. Z21 uses number counts from the MORA Survey, alongside 1.2 mm number counts from ASPECS and 3 mm number counts, to derive constraints on the high-$z$ IRLF, even in the absence of direct redshift measurements of individual galaxies. Thus it is important to emphasize that the predicted redshift
<figure><img src="image_10.png" /><figcaption>**Figure 5.** The measured cumulative redshift distribution for the 11 dust continuum dominated 2 mm galaxies detected in the MORA maps compared to literature models. Each galaxy's redshift probability density distribution is coadded here to reflect the uncertainty in the total CDF; we measure the uncertainty on the median redshift estimate, $\langle z_{2\,mm} \rangle = 3.6_{-0.3}^{+0.4}$, through Monte Carlo draws from the CDF. Literature model redshift CDFs are measured by drawing galaxies from the model with matched flux densities to our real sources (within uncertainties) and inferring the redshift distributions of those sources. The left panel compares to three hypothetical universe predictions using the Casey et al. (2018a) model framework: the Casey et al. “dust poor” model A universe (blue), the “dust rich” model B universe (orange), and the Zavala et al. (2018a) 3 mm number counts based model (green). The right panel compares with three other models: the SIDES model (purple; Bethermin et al. 2017), the SHARK model (red; Lagos et al. 2020), and the most recent 1 mm/2 mm/3 mm updated number counts based model from Z21. Note that the Z21 model is largely based on our MORA maps, but it does not use any source redshift information as input, and thus its redshift distribution predictions are derived independent of these data.</figcaption></figure>

distribution from Z21 constitutes an independent prediction, despite the use of MORA number counts to generate the model constraints.

Also shown on the right side of Figure 5 are the results from the three semi-empirical models grounded in cosmological simulations: SIDES (Bethermin et al. 2017), SHARK (Lagos et al. 2020), and UNIVERSEMACHINE (Popping et al. 2020). All models, with exception of the Popping et al. (2020) model, which is only in slight tension, accurately predict the number of sources to be found in MORA within uncertainties due to cosmic variance. Similarly, their redshift distributions are also in broad agreement. Nevertheless, within the uncertainties of the redshift distributions, there is a systematic offset where most models predict median redshifts of $\langle z_{2\,mm} \rangle = 2.2–3.2$ versus the observed $\langle z_{2\,mm} \rangle = 3.6$. The SHARK model comes closest to the measured median redshift, though we note the SIDES model has a more prominent high-redshift tail to its distribution, similar to the MORA distribution. Only further data can reduce the uncertainties on the redshift distribution and discriminate between these models.

In summary, comparing the measured redshift distribution and number counts of our MORA results with models suggests that, indeed, the prevalence of DSFGs beyond $z \gtrsim 5$ is inherently low. This is well aligned with predictions from cosmological simulations, which fundamentally limit the buildup of massive star-forming galaxies in the first ~2 Gyr of the universe’s history directly from the volume density of massive halos at that epoch.

Despite this verification of models, there are still subtleties to these measurements worth exploring that could help refine early universe DSFG volume densities further. For example, Figure 6 shows the relationship between 2 mm flux density and redshift for the MORA sample and model predictions. The Popping et al. (2020) model predicts a low average redshift for a 2 mm luminous sample that does not vary as a function of flux density. While SIDES and SHARK accurately predict the redshift distribution and source density in the MORA map, they

offer different predictions of which sources are most likely to sit at the highest redshifts; in other words, while SIDES predicts sources to sit at higher redshifts with brighter 2 mm flux density, SHARK predicts that the average redshift for bright ($S_{2\,mm} \gtrsim 0.5$ mJy) 2 mm sources is relatively low compared to faint ($S_{2\,mm} \lesssim 0.5$ mJy) sources. Our data—though limited severely by the small sample size—suggest that brighter sources tend to sit at higher redshifts (see also Koprowski et al. 2017). This is also in line with the predictions of both the Casey et al. (2018a) Model A and Z21 adjusted number counts based MORA model, both of which show the highest average redshift for $S_{2\,mm} \gtrsim 1$ mJy sources.

The origins of these second-order discrepancies between model predictions and data will inevitably require larger samples of 2 mm galaxies identified over wider areas. Nevertheless, we discuss possible origins of such discrepancies (and other potential degeneracies in our conclusions) later in Section 6.5; possible causes include evolution in the faint-end slope of the IRLF, evolution of galaxies’ dust emissivity spectral indices, or evolution of galaxies’ bulk luminosity-weighted dust temperatures.

## 5.2. Direct SED Fits Using Redshift Priors

We fit the FIR through millimeter spectral energy distributions of the 2 mm detected sample using a modified blackbody plus a mid-infrared power law; this procedure is a modified version of the fitting technique described in Casey (2012) and will be described in full in a forthcoming paper (Drew et al., in preparation). The difference with the Casey (2012) analytical approximation is that the blackbody and mid-infrared power law are added together as a piecewise function (where the mid-infrared power law is joined at the point where the slope of the blackbody is equal to the power-law index $\alpha_{MIR}$). The
<figure><img src="image_11.png" /><figcaption>**Figure 6.** Top: 2 mm flux density against best redshift constraints for the MORA sample; sources are labeled by their IDs. Bottom: the cumulative median redshift of the MORA sample above a given 2 mm flux density (thick black line; dashed lines enclose the inner 68th percentile). The median redshifts of five models are overplotted for comparison: the SIDES model from Betherrin et al. (2017) in purple (dashed line), the SHARK model from Lagos et al. (2020) in red (dashed line), the UNIVERSEMACHINE model from Popping et al. (2020) in yellow (dotted-dashed line), the dust poor model from Casey et al. (2018a) in light blue (solid), and the Z21 model in teal (solid). The dearth of sources >1 mJy in this data set limits our ability to distinguish between models, with exception of the Popping et al. (2020) model that underestimates the redshifts of 2 mm selected galaxies of all flux densities.</figcaption></figure>

functional form of the fit used is:

$$S_{\nu}(\lambda) = \begin{cases} N_{\text{pl}} \lambda^{\alpha} & : \lambda \mid \frac{\partial \log S_{\nu}}{\partial \log \lambda} > \alpha \\ N_{\text{bb}} \frac{(1 - e^{-(\lambda_{0}/\lambda)^{\beta}}) \lambda^{-3}}{e^{hc/\lambda kT} - 1} & : \lambda \mid \frac{\partial \log S_{\nu}}{\partial \log \lambda} \le \alpha \end{cases}. \quad (1)$$

Here, $\lambda$ is the rest-frame wavelength, $T$ the modeled dust temperature, $\lambda_0$ the wavelength where the dust opacity is $\tau = 1$, $\alpha$ is the slope of the mid-infrared power law, and $\beta$ is the observed integrated emissivity spectral index. $N_{\text{pl}}$ and $N_{\text{bb}}$ are normalization constants whose values are tied to one another such that the SED is contiguous at the point of transition, i.e., where

<figure><img src="image_12.png" /><figcaption>**Figure 7.** An illustrative example of the two-dimensional joint posterior distributions for the IR SED fit to source MORA-2. The parameters fit are redshift $z$, the mid-infrared power-law slope $\alpha$, the emissivity spectral index $\beta$, the IR luminosity $L_{\text{IR}}$, and the rest-frame peak wavelength $\lambda_{\text{peak}}$. MORA-2 shows strong covariance between $z$, $L_{\text{IR}}$, and $\lambda_{\text{peak}}$. The resulting best-fit SED for MORA-2 is shown in Figure 8.</figcaption></figure>

$\partial \log S_{\nu} / \partial \log \lambda = \alpha$. Both $N_{\text{pl}}$ and $N_{\text{bb}}$ are tied to $L_{\text{IR}}$, defined as the integral of the SED in $S_{\nu}$ in the range 8–1000 µm.

We adopt the general opacity model for all SEDs; for lack of direct constraints on the dust opacity in these systems, we adopt $\lambda_0 = 200$ µm. This is broadly consistent with what is seen in high-$L_{\text{IR}}$ systems where $\lambda_0$ can be measured (e.g., Conley et al. 2011; Spilker et al. 2016); even if this presumption is incorrect in the case of these systems, the exact adopted value does not impact the results of this work. Specifically, it only impacts the functional relationship between $T$ and $\lambda_{\text{peak}}$, where the latter quantity is insensitive to choice of $\lambda_0$; for example, Cortzen et al. (2020) show that the measured dust temperature for GN20 may vary by 20 K depending on the adopted $\lambda_0$, whereas the measured $\lambda_{\text{peak}}$ remains unimpacted. So while it is known there are degeneracies in measured quantities $\beta$, $T$, and $\lambda_0$ that cannot be resolved with the limited data we have on this sample in full, it is worth emphasizing that measurement of $\lambda_{\text{peak}}$ is independent and not sensitive to the choice of opacity model. For the purposes of SED fitting, the photometric points are allocated an additional 10% calibration error, representing the degree to which the absolute value of the flux density across the (sub)millimeter regime is known.

We converge on best-fit SEDs using a Markov Chain Monte Carlo (MCMC) routine. In all cases, we use the best available redshift probability distribution as a prior. Each SED is handled individually, whereby the choice to fix a variable or let it vary depends on the degree of photometric constraints relevant to each variable; for example, $\alpha_{\text{MIR}}$ is fixed to $\alpha_{\text{MIR}} = 4$ if there are no detections at rest-frame mid-infrared wavelengths, and $\beta$ is fixed to $\beta = 1.8$ (Planck Collaboration et al. 2011) if there are no more than 2–3 photometric points on the Rayleigh–Jeans tail of the blackbody emission in the millimeter regime. The choice of fixed $\alpha_{\text{MIR}} = 4$ is motivated by fits to lower-redshift DSFGs with constrained mid-IR emission (e.g., at $z \sim 2$; Pope
<figure><img src="image_13.png" /><figcaption>**Figure 8.** The best-fit dust SEDs for the $>5\sigma$ 2 mm detected sample; SEDs are fit as a modified blackbody joined with a mid-infrared power law. The modified blackbody fixes the opacity such that $\tau = 1$ at $\lambda_0 = 200~\mu$m near the intrinsic peak of the dust SED; this assumption is not meant to be physically interpreted, but is rather fixed for convenience, since our data cannot independently constrain $\lambda_0$. Fixing $\lambda_0$ has no impact on the measured $\lambda_{peak}$. Four of these sources are detected at 100–160 $\mu$m with Herschel PACS, leading to a more prominent mid-infrared component than those that lack rest-frame mid-infrared detections. The uncertainty of the fit is shown via the light blue SED fits, drawn randomly from the successful MCMC trials. Sources without rest-frame mid-infrared data have fixed mid-infrared power-law slopes of $\alpha_{MIR} \equiv 4$, and MORA-9, the only source with limited photometry on the Rayleigh–Jeans tail, has a fixed $\beta \equiv 1.8$. Measured SED characteristics are given in Table 2 and the sources’ photometry is given in Table 3.</figcaption></figure>

et al. 2008), while the choice of $\beta = 1.8$ is motivated by well-constrained Rayleigh–Jeans emission in similar systems (e.g., Scoville et al. 2016). The free parameters for our least-constrained SED (MORA-9) are redshift, $L_{IR}$, and $\lambda_{peak}$ (where both $\alpha_{MIR}$ and $\beta$, the emissivity spectral index, have been fixed). By contrast, the best-fit SEDs constrain up to five free parameters: $z$, $L_{IR}$, $\lambda_{peak}$, $\alpha_{MIR}$, and $\beta$.

A graphical representation of the two-dimensional posterior distributions for each variable is shown in Figure 7 for MORA-2.

Figure 8 shows the resulting SED fits against measured source photometry. We see a variety of SED shapes presented here, from sources with prominent mid-infrared power laws to quite steep Rayleigh–Jeans tails ($\beta > 2$). The uncertainties are

illustrated by the range of accepted MCMC trial SED solutions, shown as light blue curves in Figure 8. Extracted fit characteristics are given in Table 2.

### 5.3. Distribution in $L_{IR}$, $\lambda_{peak}$, $M_{ISM}$, and $\beta$

Figure 9 shows the distribution of the full sample of 12 galaxies in $L_{IR}$, $\lambda_{peak}$, $M_{ISM}$, and $\beta$, demonstrating the relative heterogeneous nature of the sample. Sources are color coded by the quality of their redshift constraints, and contours denote 1–2$\sigma$ confidence intervals in the given parameter space.

The dynamic range in MORA sources’ IR luminosities is about a factor of ~25, from the most extreme, MORA-1, topping $2 \times 10^{13} L_{\odot}$, to the least extreme, MORA-9, ~9 $\times 10^{11} L_{\odot}$ (the
<table><caption>**Table 2**
Measured Physical Characteristics of the 2 mm Detected Sample</caption><thead><tr><th>NAME</th><th>z</th><th>z-TYPE</th><th>L<sub>IR</sub><br>(L<sub>☉</sub>)</th><th>SFR<br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>λ<sub>peak</sub><br>(μm)</th><th>β</th><th>α<sub>MIR</sub></th><th>M<sub>dust</sub><br>(M<sub>☉</sub>)</th><th>M<sub>*</sub><br>(M<sub>☉</sub>)</th></tr></thead><tbody><tr><td>MORA-0</td><td>3.3 ± 0.8</td><td>OIR</td><td>(4.7<sup>+1.6</sup><sub>-1.0</sub>) × 10<sup>12</sup></td><td>690<sup>+230</sup><sub>-150</sub></td><td>133<sup>+15</sup><sub>-12</sub></td><td>2.3<sup>+0.4</sup><sub>-0.1</sub></td><td>3.6<sup>+1.1</sup><sub>-0.7</sub></td><td>(1.4<sup>+0.5</sup><sub>-0.4</sub>) × 10<sup>9</sup></td><td>(2.2<sup>+0.8</sup><sub>-0.8</sub>) × 10<sup>10</sup></td></tr><tr><td>MORA-1</td><td>3.78<sup>+0.27</sup><sub>-0.32</sub></td><td>OIR</td><td>(2.0 ± 0.3) × 10<sup>13</sup></td><td>3000<sup>+500</sup><sub>-400</sub></td><td>68 ± 4</td><td>2.4<sup>+0.2</sup><sub>-0.4</sub></td><td>5.2<sup>+1.2</sup><sub>-1.6</sub></td><td>(3.6<sup>+1.6</sup><sub>-1.3</sub>) × 10<sup>8</sup></td><td>(5.8<sup>+1.5</sup><sub>-1.9</sub>) × 10<sup>10</sup></td></tr><tr><td>MORA-2</td><td>3.36<sup>+0.60</sup><sub>-0.28</sub></td><td>OIR</td><td>(9.7<sup>+1.4</sup><sub>-2.7</sub>) × 10<sup>12</sup></td><td>1400<sup>+700</sup><sub>-400</sub></td><td>84<sup>+11</sup><sub>-12</sub></td><td>2.0 ± 0.3</td><td>5.6<sup>+1.8</sup><sub>-2.9</sub></td><td>(6.0<sup>+1.2</sup><sub>-2.3</sub>) × 10<sup>8</sup></td><td>(1.4<sup>+0.4</sup><sub>-0.3</sub>) × 10<sup>11</sup></td></tr><tr><td>MORA-3</td><td>4.63</td><td>Spec</td><td>(1.17<sup>+0.18</sup><sub>-0.16</sub>) × 10<sup>13</sup></td><td>1730<sup>+270</sup><sub>-240</sub></td><td>85<sup>+4</sup><sub>-5</sub></td><td>3.11<sup>+0.15</sup><sub>-0.26</sub></td><td>4.4<sup>+1.3</sup><sub>-1.1</sub></td><td>(2.1<sup>+0.5</sup><sub>-0.6</sub>) × 10<sup>8</sup></td><td>...</td></tr><tr><td>MORA-4</td><td>5.85</td><td>Spec</td><td>(3.6<sup>+1.1</sup><sub>-0.9</sub>) × 10<sup>12</sup></td><td>530<sup>+160</sup><sub>-130</sub></td><td>96<sup>+12</sup><sub>-11</sub></td><td>2.1<sup>+0.1</sup><sub>-0.3</sub></td><td>≡ 4</td><td>(6.1<sup>+1.7</sup><sub>-0.7</sub>) × 10<sup>8</sup></td><td>(3.2<sup>+1.0</sup><sub>-1.5</sub>) × 10<sup>9</sup></td></tr><tr><td>MORA-5</td><td>4.3<sup>+1.5</sup><sub>-1.3</sub></td><td>Hybrid</td><td>(4.1<sup>+1.8</sup><sub>-1.8</sub>) × 10<sup>12</sup></td><td>610<sup>+100</sup><sub>-270</sub></td><td>99 ± 21</td><td>2.2<sup>+0.3</sup><sub>-0.5</sub></td><td>≡ 4</td><td>(6.2<sup>+3.2</sup><sub>-3.0</sub>) × 10<sup>8</sup></td><td>(1.5<sup>+1.0</sup><sub>-0.7</sub>) × 10<sup>11</sup></td></tr><tr><td>MORA-6</td><td>3.34<sup>+0.13</sup><sub>-0.12</sub></td><td>OIR</td><td>(1.43<sup>+0.21</sup><sub>-0.18</sub>) × 10<sup>13</sup></td><td>2120<sup>+310</sup><sub>-260</sub></td><td>80 ± 4</td><td>2.4<sup>+0.3</sup><sub>-0.2</sub></td><td>7 ± 2</td><td>(5.2<sup>+1.1</sup><sub>-2.0</sub>) × 10<sup>8</sup></td><td>(7.0<sup>+1.7</sup><sub>-1.4</sub>) × 10<sup>10</sup></td></tr><tr><td>MORA-7</td><td>2.85<sup>+0.24</sup><sub>-0.33</sub></td><td>OIR</td><td>(3.5<sup>+1.5</sup><sub>-0.9</sub>) × 10<sup>12</sup></td><td>510<sup>+220</sup><sub>-130</sub></td><td>105 ± 13</td><td>2.2<sup>+0.4</sup><sub>-0.3</sub></td><td>6<sup>+3</sup><sub>-1</sub></td><td>(4.4<sup>+1.6</sup><sub>-1.6</sub>) × 10<sup>8</sup></td><td>(2.5<sup>+0.5</sup><sub>-0.5</sub>) × 10<sup>11</sup></td></tr><tr><td>MORA-8</td><td>2.29<sup>+0.12</sup><sub>-0.08</sub></td><td>OIR</td><td>(2.2 ± 0.5) × 10<sup>12</sup></td><td>330<sup>+100</sup><sub>-70</sub></td><td>131<sup>+12</sup><sub>-10</sub></td><td>1.8<sup>+0.3</sup><sub>-0.2</sub></td><td>6.6<sup>+2.0</sup><sub>-2.2</sub></td><td>(1.0<sup>+0.4</sup><sub>-0.2</sub>) × 10<sup>9</sup></td><td>(1.1<sup>+0.3</sup><sub>-0.3</sub>) × 10<sup>11</sup></td></tr><tr><td>MORA-9</td><td>4.3<sup>+1.3</sup><sub>-1.0</sub></td><td>Hybrid</td><td>(1.2<sup>+1.0</sup><sub>-0.4</sub>) × 10<sup>12</sup></td><td>180<sup>+150</sup><sub>-60</sub></td><td>123<sup>+26</sup><sub>-23</sub></td><td>≡ 1.8</td><td>≡ 4</td><td>(8.2<sup>+5.1</sup><sub>-3.8</sub>) × 10<sup>8</sup></td><td>(4.1<sup>+1.8</sup><sub>-1.4</sub>) × 10<sup>10</sup></td></tr><tr><td>MORA-10</td><td>2.472</td><td>Spec</td><td>(1.1<sup>+0.1</sup><sub>-0.1</sub>) × 10<sup>13</sup></td><td>1570<sup>+160</sup><sub>-130</sub></td><td>102<sup>+10</sup><sub>-13</sub></td><td>1.7<sup>+0.4</sup><sub>-0.2</sub></td><td>1.3<sup>+0.15</sup><sub>-0.14</sub></td><td>(7.1<sup>+2.2</sup><sub>-2.3</sub>) × 10<sup>8</sup></td><td>(7.3<sup>+2.0</sup><sub>-2.0</sub>) × 10<sup>10</sup></td></tr><tr><td>MORA-11</td><td>3.17 ± 0.05</td><td>OIR</td><td>(3.4<sup>+1.1</sup><sub>-1.1</sub>) × 10<sup>12</sup></td><td>500<sup>+150</sup><sub>-170</sub></td><td>99<sup>+13</sup><sub>-14</sub></td><td>2.1<sup>+0.4</sup><sub>-0.4</sub></td><td>≡ 4</td><td>(3.5<sup>+2.9</sup><sub>-2.9</sub>) × 10<sup>8</sup></td><td>(1.7<sup>+0.5</sup><sub>-0.5</sub>) × 10<sup>11</sup></td></tr></tbody></table>

Note. Derived physical properties of the 2 mm detected MORA Sample. Variables that have been fixed are denoted with ≡ notation. Estimates for ISM masses (not given in this table) can be obtained by multiplying the dust masses in this table by a factor of 125 (Rémy-Ruyer et al. 2014). A stellar mass estimate is not available for MORA-3 (AzTEC-2) due to blending of near-infrared imaging with foreground galaxies. All uncertainties in this table indicate the inner 68% minimum credible interval of posterior distributions.

corresponding star formation rates range 140–3100 $M_{\odot}$ yr$^{-1}$). Given the selection of these sources on the Rayleigh–Jeans tail of dust blackbody emission at 2 mm, the dynamic range in ISM masses is much more narrow, mirroring the somewhat narrow dynamic range in 2 mm flux densities (both a factor of ~3). The second panel of Figure 9 shows the measured ISM masses of the sample, scaled up by a factor of 125 from the SED-inferred dust mass (to account for the total mass of gas in the ISM), which is consistent with the expected gas-to-dust ratio for near solar metallicity galaxies (Rémy-Ruyer et al. 2014). Note that a maximum ISM mass, as a function of redshift, is set in our survey by the size of the survey itself given a ΛCDM universe (purple shaded regions in the $M_{\text{ISM}}-z$ panel of Figure 9, at 1–3σ confidence intervals); this peak is determined by the maximum halo mass as a function of redshift and survey volume from Harrison & Hotchkiss (2013), using an ISM-to-halo mass ratio of 1/20 (as is done in Marrone et al. 2018).

While the nominal 2 mm detection limit in $L_{\text{IR}}$ with redshift follows the strong negative K-correction (and thus is more sensitive to sources at $z \gtrsim 3$ than at $z = 1–2$; denoted by gray bands), the detection limit in $M_{\text{ISM}}$ is approximately flat. These detection curves are not absolute limits, as they are sensitive both to galaxies’ luminosity-weighted dust temperature and the emissivity spectral index. For example, a galaxy of ISM mass $\sim 2 \times 10^{10} M_{\odot}$ at $z = 4$ may only be detectable in MORA if its emissivity spectral index is lower than $\beta \lesssim 2.2$, or a $\sim 2 \times 10^{12} L_{\odot}$ system at $z = 3$ may only be detectable in MORA if its rest-frame peak wavelength falls within $95 \lesssim \lambda_{\text{peak}} \lesssim 175 \ \mu\text{m}$.

The galaxies’ distribution in rest-frame peak wavelength, $\lambda_{\text{peak}}$ (a proxy for the luminosity-weighted dust temperature) is shown in the third panel of Figure 9 relative to the aggregate DSFG population fit found in Casey et al. (2018a), with intrinsic scatter of $1–2\sigma$ shown. The MORA galaxies are largely in agreement with the global DSFG trend, with two sources appearing to be somewhat anomalously cold (MORA-0 and MORA-8). For a sample this small, we would only expect at most one source to sit $>2\sigma$ offset from the global trend given that galaxies are distributed in a Gaussian about the mean $L_{\text{IR}}-\lambda_{\text{peak}}$ relation. The cold SED for MORA-0 is somewhat uncertain given the significant uncertainty on the source’s redshift; if this source is

indeed at the higher-redshift end of its redshift probability distribution function (PDF), its SED would be better aligned with the overall trend for DSFGs’ peak wavelengths. Nevertheless, it is worth highlighting that 2 mm selection itself may skew the temperatures of the sample a bit cold; the gray region in the third panel of Figure 9 shows the parameter space beyond reach of the MORA survey at a fixed redshift of $z = 3$. If the survey had a flat selection in $L_{\text{IR}}$ around $\sim 3 \times 10^{12} L_{\odot}$, MORA-9 and MORA-8 may not have made the cut, as they are the least luminous sources in the sample; excluding them, only one source, MORA-0, is anomalously cold, in line with expectation for a sample this size.

The average value of the emissivity spectral index of MORA galaxies in the sample is $\langle\beta\rangle = 2.2_{-0.4}^{+0.5}$, which is skewed toward higher values than are often typically assumed for galaxies in the absence of direct measurements. The fact that the MORA sample is 2 mm selected would potentially impact the average measured $\beta$. However, one would assume it would do so in the opposite sense, by more efficiently identifying sources with shallower emissivity spectral indices (or $\beta < 1.8$) due to higher relative 2 mm flux densities for a given $L_{\text{IR}}$. While nominally substantial heating from the CMB at high $z$ ($z \gtrsim 5$) could effectively steepen the Rayleigh–Jeans slope (and artificially increase $\beta$ as discussed in Jin et al. 2019), here we account for CMB heating directly and find its effect negligible. There is no evidence in the MORA sample of evolution in $\beta$; because the sample size is relatively small, the mean could also not be precisely constrained. Our interpretation of $\beta$, relative to other literature samples and other DSFGs falling within the MORA footprint, is offered in Section 6.3.

## 5.4. Stellar Masses and the SFR–$M_{\star}$ Relation

Stellar masses are derived using the suite of OIR COSMOS photometric constraints available for each source (Weaver et al. 2021). While Weaver et al. use LePhare (Arnouts et al. 2002; Ilbert et al. 2006) and a range of 19 stellar population templates from Bruzual & Charlot (2003), we only make use of their posterior redshift probability density distributions for the subset of MORA sources detected in the COSMOS2020 catalog. We choose to remodel the galaxies’ stellar populations using a wider range of templates inclusive of extreme starbursts. This refitting is carried out
<figure><img src="image_14.png" /><figcaption>**Figure 9.** The distribution of the 2 mm detected sample in derived physical parameter space. Sources are color coded by the origins of their redshift constraints: gray (spectroscopic), blue (OIR photometric), and teal (hybrid OIR/millimeter photometric). Contours denote sources' 1 and 2σ confidence intervals for the measured parameters, highlighting parameter covariance for the given set of photometry. Regions of the planes inaccessible to the MORA Survey are shaded in gray in the first three panels: the cutoffs depend precisely on parameters of the SED fit. For example, in the $L_{IR}-z$ plane (top left), the gray curves trace out the detection limit of a source detected at $5\sigma$ ($1\sigma = 70~\mu	ext{Jy at 2 mm}$) with $\lambda_{\text{peak}} \equiv 100~\mu\text{m}$ and $\beta = 1.8$, $\beta = 2.2$ and $\beta = 2.6$ (with $\beta = 2.6$ corresponding to the flattest curve). In the $M_{\text{ISM}}-z$ plane (top right) the lower detection threshold corresponds to fits with $\beta = 1.8-2.6$ (where the steepest curve corresponds to $\beta = 2.6$). The light purple shaded regions in the $M_{\text{ISM}}-z$ plane also shows the theoretical maximum $M_{\text{ISM}}$ mass limit as a function of redshift based on the Harrison & Hotchkiss (2013) halo mass survey estimator; the halo mass is scaled down to $M_{\text{ISM}}$ by a factor of 20. The $L_{IR}-\lambda_{\text{peak}}$ plane (bottom left) shows the aggregate, unevolving $L_{IR}-\lambda_{\text{peak}}$ relation from Casey et al. (2018a; gray line with dashed lines enclosing $1\sigma$ and $2\sigma$ scatter) and agreement with the MORA sample SEDs. The gray region marks the detection limit for a $z=3$ system (the curves shift in $L_{IR}$ with $z$ following the gray curves in the $L_{IR}-z$ panel). The distribution of sources' emissivity spectral indexes ($\beta$) is shown against redshift (bottom right) with the canonical value of $\beta = 1.8$ marked with the black dashed line. No evidence for redshift evolution of $\beta$ is seen in this (small) sample; the average value is $\langle\beta\rangle = 2.2_{-0.4}^{+0.5}$.</figcaption></figure>

using the MAGPHYS energy balance code (da Cunha et al. 2008), Bruzual & Charlot (2003) stellar population synthesis templates, and an updated, wider range of star formation histories compatible with DSFGs (da Cunha et al. 2015). Redshift uncertainties are accounted for by iteratively sampling the sources’ redshift PDFs; we find that stellar masses are largely insensitive to redshift uncertainties. Similarly, the stellar masses derived from MAGPHYS show no systematic offset with those reported by Weaver et al., though significant uncertainty on all derived stellar masses remains.

There are a few exceptions where our technique is approached with even more caution: MORA-3, MORA-4, MORA-5, and MORA-9. The first, MORA-3 (AzTEC-2), lacks any OIR counterpart shortward of $2.2~\mu\text{m}$ and at longer

wavelengths is spatially confused with foreground galaxies, thus rendering any stellar mass constraint impossible without future spatially resolved JWST observations. Though the MAGPHYS predicted stellar mass for MORA-4 (MAMBO-9) is reasonable on a standalone basis, Casey et al. (2019) present a detailed argument as to why it is most likely an overestimate given the measured gas mass and implied halo mass. Therein they provide an alternate stellar mass estimate for MAMBO-9 following the methodology outlined in Finkelstein et al. (2015), which uses stellar population modeling plus the contribution from nebular emission. Note that the necessity of this approach is thought to be unique to that source, given that MAMBO-9 is likely the highest-redshift source in the sample and is
<figure><img src="image_15.png" /><figcaption>**Figure 10.** Placement of MORA 2 mm selected sample in the galaxy M<sub>☉</sub>-SFR relation, or galaxy “main sequence.” Colored lines and ±0.3 dex shaded regions represent the z = 2 (purple), z = 3 (blue), z = 4 (green), z = 5 (orange), and z = 6 (red) average relations from Speagle et al. (2014). The dashed navy line represents the lower envelope for starbursts as defined in Caputi et al. (2017, 2021). Galaxies are labeled by their MORA IDs and color coded according to the Δz = 1 interval in which they most likely sit (e.g., MORA-0 with z<sub>OIR</sub> = 3.3 is blue to denote it sits at 3 < z < 4). All sources except MORA-3 (AzTEC-2) have stellar mass estimates, and thus that source's SFR is noted as a horizontal green band with an unconstrained stellar mass; MORA-3 lacks an estimate due to severe blending of IRAC photometry with foreground galaxies. With respect to the main sequence, the MORA sample splits roughly in half: five out of 12 galaxies sit on the galaxy main sequence, while seven out of 12 are significantly elevated, the most extreme of which is MORA-4 (MAMBO-9; Casey et al. 2019), which is (likely) the highest-redshift source of the sample.</figcaption></figure>

intrinsically massive; thus the survey volume itself sets an
upper limit to its mass.

The last two exceptions to the standard MAGPHYS stellar mass fitting approach are sources MORA-5 and MORA-9, which are the two OIR-dark sources explored in more detail in Manning et al. (2021), as they lack good redshift constraints. Their stellar masses are fit using the MAGPHYS+photoz code (Battisti et al. 2019), which is a further update to MAGPHYS allowing for redshift uncertainty; this was deemed necessary in the case of MORA-5 and MORA-9 in particular, given the significant uncertainties on both sources’ photometric redshifts. The resulting best-fit stellar masses are reported in Table 2.

Figure 10 shows the distribution of the MORA sample in the galaxy SFR–M<sub>☉</sub> plane relative to observed trends for the bulk galaxy population at various redshifts (i.e., the galaxy “main sequence”). The trend lines are drawn from Speagle et al. (2014) at fixed redshifts and sources are color coded by redshift, and the dashed navy line demarcates the lower envelope for “starbursts” (specific star formation rate, sSFR > 25 Gyr<sup>-1</sup>; Caputi et al. 2017, 2021). What we observe is a rather heterogeneous sample: 7/12 (58%) MORA-selected sources are significantly elevated above the “main sequence” with specific star formation rates in excess of 10 Gyr<sup>-1</sup> (with 5/12, 42%, above the 25 Gyr<sup>-1</sup> threshold).<sup>37</sup> The remaining

5/12 (42%) are embedded in the galaxy main sequence itself, with 1 Gyr⁻¹ < sSFR < 10 Gyr⁻¹. There does not appear to be a strong correlation between sources' redshifts and whether or not they have elevated sSFRs. The fact that the sample is somewhat heterogeneous in the SFR–M* plane is not surprising given that the 2 mm selection (and submillimeter selection more broadly) roughly corresponds to an SFR cut, modulo variation due to SED dust temperature and emissivity. Later in Section 6.2 we discuss some of the implications of the sample's measured stellar content in relation to the formation of the first massive quiescent galaxies.

With a median stellar mass of $7 \times 10^{10} M_{\odot}$, we note that the median dust-to-stellar mass ratio in this sample is $M_{\text{dust}}/M_{*} = (7_{-3}^{+30}) \times 10^{-3}$. This is nominally a bit higher than expected for DSFGs, though one would expect a 2 mm sample to be slightly biased in that regard, selecting more galaxies rich in dust per unit stellar mass than would selection at shorter wavelengths.

5.5. 2 mm Characteristics of Known (Sub)Millimeter Sources

Here, we analyze the 2 mm characteristics of known submillimeter sources in the MORA footprint, taken from the Simpson et al. (2019) deep SCUBA-2 850 µm map. We select sources from that work detected above a signal-to-noise ratio of 4, roughly corresponding to a flux threshold of $S_{850} \gtrsim 2.3$ mJy (i.e., rms ≈ 0.57 mJy beam⁻¹), minimizing spurious detections at lower significance. Only 55/98 (≈56%) sources have interferometric follow-up from ALMA other than our band 4 data, allowing identification of their precise astrometric positions in the MORA map. The remainder (43/98 ≈ 44%) lack precise position constraints and are not well known to better than the James Clerk Maxwell Telescope (JCMT) beam (14′′8). For those sources, we search for possible corresponding marginally detected 2 mm sources by identifying the highest significance peak within the FWHM of the JCMT beam.

Out of 98 >4σ SCUBA-2 sources overlapping with the main MORA mosaic, 17/98 (=17.3%) have 2 mm peaks detected at >4σ (including the 12 sources in our 5σ sample). This compares to 3.7% of randomly placed JCMT beam “apertures” containing >4σ sources. Most of this excess signal comes from our >5σ sample. However, our data still show an excess between 4 < σ < 5 significance with 5/86 (=5.8%) found in total compared to the random “aperture” rate of 3.7%. We find 56/98 (=57% ± 7%) SCUBA-2 sources have between 3 < σ < 4 significance peaks in the 2 mm map, compared to 48% ± 2% of randomly placed apertures. Of those sources with 3 < σ < 4 2 mm counterparts, the median 850 µm flux density is 3.4 mJy, lower than the median for the 4 < σ < 5 sample with 6.8 mJy. While the excess signal over random apertures is of relatively low significance due to significant potential for contamination, it is still indicative of real 2 mm emission associated with these known 850 µm sources.

Fifty-five of the 98 SCUBA-2 sources sitting in the MORA footprint have other ALMA continuum data covered by the archival A3COSMOS project (described in detail in Liu et al. 2018), and 45/55 have a confirmed ALMA counterpart detected at or above 3σ in their corresponding archival ALMA data (of varying observed frequencies, most at 870 µm or 1.2 mm); seven of the 45 sources have two components or sources (separated by more than 1") detected by ALMA within the JCMT beam (and two out of seven have one of their multiples

<footnote n="37"><sup>37</sup> Note that this likely includes MORA-3, which does not have a stellar mass constraint, but is unlikely to have a stellar mass in excess of ~3 × 10<sup>11</sup> M<sub>☉</sub>, which would be required for it to dip below this threshold. See Jiménez-Andrade et al. (2020) for more details on this source.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 11.** The 870 µm to 2 mm flux density ratio against redshift and in histogram for DSFGs from various samples. Left: the S<sub>870</sub>/S<sub>2000</sub> ratio for DSFGs that have secure redshift constraints, including >5σ MORA sources (dark gray points) and marginal MORA sources (light gray points); we have also overplotted two other samples of well-characterized DSFGs in the literature: the SPT (1.4 mm selected) sample from Reuter et al. (2020), and the ALESS (870 µm selected) sample from da Cunha et al. (2021). Overlaid are tracks of specific SEDs with variable dust temperature (λ<sub>peak</sub> = 95 µm or 125 µm, solid vs. dashed) and variable emissivity spectral index (β = 1.8, 2.2, or 2.6, in increasingly dark shades of orange, respectively). Right: histograms of S<sub>870</sub>/S<sub>2000</sub> for sources without redshifts that have 870 µm ALMA data (dark gray), sources with 870 µm data only from SCUBA-2 (light gray), and the full distributions of the SPT sample (purple dotted–dashed) and ALESS sample (dotted blue). Downward arrows mark the maximum color ratio that could correspond to a given value of β at z = 2.5; higher ratios are not allowed for the given β, as the ratio of flux densities only diminishes with redshift. Our data suggest a range of β values apply to this sample, from some being fully consistent with the often-presumed value of β = 1.8, to some being more consistent with steeper slopes, β ≈ 2.2–2.6.</figcaption></figure>

detected at 2 mm at >5σ). The highest S/N 2 mm detection within
the JCMT beam corresponds to the correct, independently identified
(single) ALMA counterpart in 24 of 45 cases (≈53%).

Out of the 98 sources, 55 also have some form of redshift
constraint (derived independently from the millimeter photo-
metry), ranging from spectroscopic confirmation in the
millimeter or optical/near-infrared, or through an optical/
near-infrared photometric redshift. These 55 sources with
redshift constraints are not the same 55 with independent
ALMA data; only 37 sources have both redshifts and secure
ALMA counterparts.

Figure 11 shows the distribution of millimeter colors for the sample of known SCUBA-2 sources in the field, i.e., S<sub>870</sub>/S<sub>2000</sub>, or the ratio of ALMA measured 870 µm flux density (or SCUBA-2 850 µm) to 2 mm flux density as measured by MORA. Sources are split into those with redshift constraints (left panel) and those without (right panel). We also include several model tracks in color space as a function of dust temperature (or λ<sub>peak</sub>) and emissivity spectral index (β) as well as literature samples of DSFGs with well-constrained Rayleigh-Jeans SEDs, such as the SPT 1.4 mm selected sample from Reuter et al. (2020) and the ALESS 870 µm selected sample from da Cunha et al. (2021).

While naively one might think this 870 µm to 2 mm color
ratio is indicative purely of β, the emissivity spectral index,
given that both 870 µm and 2 mm flux density measurements
are likely to sit on the Rayleigh–Jeans tail of dust blackbody
emission, the color ratio has both a temperature and redshift
dependence that can result in a significant deviation from
expectation given a fixed value of β. For example, the
approximation that the flux density on the Rayleigh–Jeans
tail of a blackbody follows S<sub>ν</sub> ∝ ν<sup>2+β</sup> is often used, while in

practice the SED only asymptotically approaches this at rest-
frame wavelengths longward of rest-frame ~2 mm. Even at
moderate redshifts (z ~ 1), our observed bands probe shorter
rest-frame wavelengths where this approximation no longer
holds. This is why there is a redshift dependence on
S<sub>870</sub>/S<sub>2000</sub> at redshifts well below the point that either
wavelength probes the peak of the SED. The temperature
dependence may also be somewhat intuitive: at colder
intrinsic temperatures, the evolution toward lower
S<sub>870</sub>/S<sub>2000</sub> ratios is steeper at lower redshifts, as the 870
µm band more quickly probes the peak of the SED than it
would for a source with a hotter SED.

We see that both MORA and literature DSFG samples roughly
follow the model trend evolution of S<sub>870</sub>/S<sub>2000</sub> with redshift;
however, at higher redshifts, the >5σ 2 mm sample (points
enclosed in black circles) also suggest consistency with steeper
values of β than the nominal β= 1.8. Figure 11 also marks the
maximum color value (in S<sub>870</sub>/S<sub>2000</sub>) allowed for a given fixed
value of β, which is set by the measured color for z = 2.5 systems
(in other words, at all higher redshifts, the color will be
substantially lower). Though there is a significant degeneracy
between S<sub>870</sub>/S<sub>2000</sub> color and redshift, β and λ<sub>peak</sub>, as indicated by
the tracks on the figure, it appears that many sources skew toward
higher colors than would be allowable for β = 1.8: both from the
sample with redshift constraints and from the distribution of colors
for sources without redshifts. This indicates that at least a
subsample of the population likely skews toward steeper values of
β. This is consistent with our finding for individual sources in the
>5σ 2 mm sample that have ⟨β⟩ = 2.2, described earlier in
Section 5.2. Further discussion of the implications of steeper
values of β are given in Section 6.3.
<figure><img src="image_17.png" /></figure>

**Figure 12.** The contribution of 2 mm sources detected in the MORA Survey to the cosmic star formation rate density (black stars; open stars where incomplete at $z \lesssim 3$). Our measurements are shown relative to the reported literature values from the Madau & Dickinson (2014) review (orange points indicating direct IR measurements, and blue points indicating rest-frame UV/optical measurements, corrected for dust attenuation). The results from our accompanying paper, Z21, for the integrated IR/obscured component is shown as an orange, shaded band out to $z \sim 8$ including the impact of cosmic variance (the black stars do not). The Z21 band is inferred from a joint analysis of 1, 2, and 3 mm number counts and the empirical model described in Casey et al. (2018a). We compare the contribution of this 2 mm selected sample with the measured contribution from 870 µm selected sources from the AS2UDS Survey above 1 mJy (gray circles) in Dudzevičiūtė et al. (2020). The measurements in this paper (black stars) account directly and only for the 11 sources found above $>5\sigma$ in the MORA Survey, and do not represent an extrapolation of a fitted luminosity function; the uncertainties in sources' redshifts and SFRs are accounted for. The dearth of sources at $z \lesssim 3$ is a direct consequence of the 2 mm survey design, meant to efficiently "filter out" the majority of DSFGs at $1 \lesssim z \lesssim 3$ that dominate cosmic star formation at its peak.

## 6. Discussion

Prior to the MORA Survey described herein, constraints on the prevalence of dust obscured star-forming galaxies beyond $z \gtrsim 3-4$ were relatively weak; refining the measurement of DSFGs' number densities and characteristics during this epoch was the survey's primary goal. Here, we discuss the implications of our measurements, from the observed contribution of DSFGs to the cosmic SFRD beyond $z \gtrsim 3$, to the buildup of massive galaxies in the first few Gyr, and to the measured dust emissivity index in early universe DSFGs. We then present the limitations of the survey, particularly given the survey area and small sample size of detected sources, as well as various model degeneracies that could be broken with future millimeter-wavelength observational campaigns.

### 6.1. Contribution to the SFRD

The question of DSFGs' contribution to the cosmic SFRD at early times has been difficult to constrain. However, the MORA Survey's design is meant to effectively filter out lower-redshift DSFGs, enabling direct measurement of the obscured galaxy contribution to the SFRD. For example, the contrast between a "dust poor" universe (e.g., Koprowski et al. 2017; Dudzevičiūtė et al. 2020) and a "dust rich" universe (Rowan-Robinson et al. 2016; Gruppioni et al. 2020) would result in very different manifestations of the MORA data set; the former predicts $\lesssim 15$ sources with a median redshift $z \sim 3.4$, and the latter predicts $\gtrsim 50$ sources with a median redshift $z \sim 5.4$. So what does our data set imply for the SFRD?

We provide direct estimates of the MORA sample's contribution to cosmic star formation in Figure 12. The SFRD contribution for this sample is estimated in the following way. First, we determine the area over which each of the 11 sources would be detectable in the MORA map above a $>5\sigma$ significance. As a reminder, we exclude MORA-10 for its synchrotron component and MORA-12 as a likely false positive. Then we sample both the redshift and SFR probability density distributions for each source through Monte Carlo trials; this method accounts for their relative covariance. Sources' SFRs are converted to an SFRD by dividing by the appropriate survey volume corresponding to the area over which a given source is accessible. A total SFRD, binned by $\Delta z = 1$ intervals, is then measured for each Monte Carlo trial, and the average and 68% minimum confidence interval of all trials is used to infer the MORA source contribution to the SFRD, as shown in Figure 12. Though redshifts in our sample only nominally span $z = 2.2-5.9$, our SFRD estimates span $z = 1.5-6.5$ from the tails in the redshift PDFs for some of the less well-constrained sources.

Overall, we find the prevalence of $z < 2$ sources and $z > 6$ sources in this 2 mm selected sample is quite low compared to the $2 < z < 6$ sample. Here, our $z < 3$ SFRD measurements are lower than that measured from other DSFG samples due to 2 mm selection serving as an effective filter for $1 < z < 3$ DSFGs; thus those points are shown as open stars to contrast with filled stars at $z > 3$, where the 2 mm sensitivity is more complete for $\gtrsim 2 \times 10^{12} L_\odot$ DSFGs, as shown in Figure 9.

The total contribution of MORA 2 mm detected DSFGs at $3 < z < 6$ appears to be roughly $\sim 30\%$ of the total cosmic
SFRD at these epochs, shown as the blue band in Figure 12. Note that what is shown in Figure 12 represents a direct accounting for the >5σ MORA-detected sample only, not an extrapolation from a luminosity function. This is clearly seen in the apparent deficit of the SFRD at $z<3$, where we know the MORA Survey was designed to filter sources out.

Our accompanying paper, Z21, provides a detailed analysis of MORA constraints on the infrared luminosity function (IRLF) and the implied contribution of obscured galaxies to the SFRD out to $z \sim 7$. This is done using a joint analysis of 1 mm, 2 mm, and 3 mm number counts, plus our empirical model (described at length in Casey et al. 2018a, 2018b; Zavala et al. 2018a) to measure the redshift evolution of the IRLF. In particular, the free parameters of the model are the evolution of $\Phi_*$ beyond $z \sim 2$, the faint-end slope $\alpha_{LF}$ (assumed not to evolve with redshift for simplicity), and the average emissivity spectral index of DSFGs $\beta$ (assumed not to evolve for simplicity). The data used to constrain the model are all aggregate (sub)millimeter number counts, with particular emphasis on 2 and 3 mm in their ability to capture the high-$z$ redshift evolution of $\Phi_*$. Z21 finds that $\Phi_* \propto (1+z)^{\psi_2}$, where $\psi_2 = -6.5_{-1.8}^{+0.8}$, and the resulting total contribution of obscured emission to the SFRD is shown by the orange band in Figure 12.

What is of particular note in our SFRD estimate is that our direct accounting of MORA-detected galaxies agrees rather well with the best estimate of the integrated IRLF inferred by number counts. This effectively means that the sources found in our map are the *only* obscured sources to be found at these high redshifts, and there is not likely to be a population of fainter sources lurking just below the detection limit (or at least any population that contributes significantly to cosmic star formation at those epochs). The same cannot be said of rest-frame ultraviolet luminosity functions (UVLF), for which there is often a significant discrepancy between the integrated SFRD contribution between directly detected galaxies and the inferred contribution from an extrapolation down the luminosity function (Finkelstein et al. 2015; Finkelstein 2016). The key difference between the IRLF and UVLF is, of course, their faint-end slopes; while the UVLF has a steep faint-end slope, suggestive of a large population of low-mass, UV-luminous galaxies, the IRLF has a very shallow faint-end slope (see also González-López et al. 2020 and Popping et al. 2020, for measurement of the faint-end slope of the IRLF from the ASPECS survey). This effect also manifests in most of the CIB having been resolved into individual, bright point sources (Béthermin et al. 2012). Our findings—agreement between direct accounting of MORA-detected sources’ contribution to the SFRD and integrating the extrapolated IRLF—are thus consistent with $\gtrsim 10^{12} L_\odot$ sources being the dominant source of obscured emission. This is functionally equivalent to most dust luminosity living in massive galaxies, if indeed high-mass galaxies roughly correspond to high-luminosity galaxies.

Also of note in our SFRD estimates is the relative agreement between an estimate based *solely* on number counts (that from Z21) and our results here, which incorporate sources’ inferred luminosities (thus SFRs) as well as redshift constraints. It is plausible that the additional redshift information introduced here would lead to some discrepancies with the Z21 model, either showing a flatter (or steeper) SFRD contribution with redshift, but the agreement of the Z21 model with these data holds. However, it should not be entirely surprising that the two agree so well across $3 < z < 6$ given the root assumption behind the Z21 model: that there are unlikely to be abrupt kinks or changes in the IRLF with redshift. While we nominally do not find any $z > 6$ DSFGs

in the MORA Survey, we know they exist even if they are rare (Marrone et al. 2018; Zavala et al. 2018b). Measuring their volume density will require larger area 2 mm campaigns, discussed further in Section 6.4.

We also draw comparison with the 870 µm selected DSFGs from the AS2UDS survey in Dudzevičiūtė et al. (2020), who present the characteristics of 870 µm SCUBA-2 selected DSFGs across a 1 deg$^2$ area survey and followed up with ALMA. The contribution of sources above 1 mJy at 870 µm is shown in Figure 12 as gray points. A direct comparison of the 870 µm selection versus the 2 mm selection technique is aligned with expectation: 870 µm efficiently selects DSFGs above $z \gtrsim 1.5$, with a high-redshift tail extending to $z \sim 6$, while 2 mm selection is weighted toward the higher redshifts only ($z > 3$). Both the AS2UDS and MORA results are well aligned with findings from the SHARK model that predicts 2 mm selected galaxies contribute ~15%–28% of the cosmic star formation at $3 < z < 6$ (Lagos et al. 2020).

Lastly, it is worth noting that our measurements are in disagreement with recent constraints on the $z > 4$ IRLF from the ALPINE survey, which infers an integrated SFRD contribution from dust obscured sources a factor of ~10 higher than found in this paper; ALPINE used serendipitous detections in very deep ALMA pointings of $z \sim 4–5$ galaxies (the primary goal being the measurement and characterization of galaxies’ [C II] characteristics) to place constraints on the high-$z$ IRLF (Gruppioni et al. 2020) and cosmic SFRD more broadly (Khusanova et al. 2021; Loiacono et al. 2021). This discrepancy could be due to a bias in the survey area used—whereby sources are physically associated with the ALPINE primary targets, even if that association is not directly known—or alternatively, could be due to increased cosmic variance in the relatively small area of the survey (~25 arcmin$^2$ versus 184 arcmin$^2$ in MORA).

Recent findings from a similar higher-redshift line survey, the REBELS program, led to the discovery of two SFR $\sim 70 M_\odot \text{yr}^{-1}$ systems at $z \sim 7$ (Fudamoto et al. 2021), whose OIR emission is completely obscured. Fudamoto et al. (2021) estimate their volume density and contribution to the SFRD, attempting to correct for the bias introduced by the targeted survey approach. They predict a much higher contribution of lower IR luminosity ($\log L_{IR} \sim 11.5–12.2$) sources to the obscured SFRD, of order the same contribution as we estimate for higher luminosity sources ($\log L_{IR} \gtrsim 12.2$). Indeed, the Z21 model would predict that such lower IR luminous sources ($\log L_{IR} \sim 11.5–12.2$) would have an order of magnitude lower contribution to the SFRD than measured by Fudamoto et al. (2021). The origin of this discrepancy is unclear and cannot be resolved without a more extensive census of the lower luminosity regime. Nevertheless, detection of such “OIR-dark” obscured galaxies within the Epoch of Reionization represents a significant leap, indicating that even modest SFR galaxies may indeed be dust rich at $z > 7$.

Because ALPINE and REBELS are intrinsically targeted surveys and not blank-field work, we omit their estimates from Figure 12 for clarity.

## 6.2. Progenitors of First Massive, Quiescent Galaxies?

An important consequence of our measurement of the prevalence of DSFGs at $z > 3$ is its implications for the formation of the first massive quiescent galaxies, already well established three billion years after the Big Bang (at $z=2$, placing the formation redshift of their stellar mass at $z > 3$; e.g., Kriek et al. 2009; Toft et al. 2014; Schreiber et al. 2018a; Merlin et al. 2019;
Marsan et al. 2020; Santini et al. 2021). Furthermore, the discovery of a substantial population of $M_* \sim 10^{11} M_\odot$ galaxies at $z \sim 4$ (Straatman et al. 2014; Marsan et al. 2020; Sherman et al. 2020; Valentino et al. 2020; Stevans et al. 2021) with a volume density of $(1.8 \pm 0.7) \times 10^{-5} \text{ Mpc}^{-3}$ requires a population of star-forming progenitors at $z \sim 5$ with SFRs $\sim 100 M_\odot \text{ yr}^{-1}$. Though constraints on the quiescent population's volume density are yet uncertain, recent observational works have constrained the number density to a few times $10^{-6}-10^{-5} \text{ Mpc}^{-3}$ (see the recent compilation of quiescent galaxy number densities in Valentino et al. 2020). The population of DSFGs selected at 2 mm (or similarly, 3 mm; see Zavala et al. 2018a; Williams et al. 2019) are such star-forming systems; does their volume density, as a function of redshift, match the quiescent galaxy population? While not all DSFGs will quench in time to transition and become quiescent by $z \sim 3-3.5$, a comparison of their measured volume densities serves as a useful benchmark.

The raw volume density of MORA-detected galaxies (with a rough luminosity cut of $\gtrsim 2 \times 10^{12} L_\odot$ or $300 M_\odot \text{ yr}^{-1}$) at $4 < z < 6$ is $\sim (3.8^{+1.3}_{-0.7}) \times 10^{-6} \text{ Mpc}^{-3}$. This requires a correction for the DSFGs' relatively short duty cycle, $\sim 100$ Myr (Ivison et al. 2011; Bothwell et al. 2013; Swinbank et al. 2014), shorter than the period of the redshift interval over which the volume density is calculated ($\sim 600$ Myr from $4 < z < 6$). The corrected volume density is $\sim (2.2^{+0.8}_{-0.4}) \times 10^{-5} \text{ Mpc}^{-3}$. Note that this roughly agrees with other previous constraints on the $z > 4$ DSFG volume density (e.g., da Cunha et al. 2015; Michałowski et al. 2017; Miettinen et al. 2017), though here the constraints are somewhat more precise. Note that it is, however, unclear if all DSFGs have universally short duty cycles, as the star formation histories of quiescent galaxies suggest some long DSFG lifetimes (e.g., Schreiber et al. 2018a; Forrest et al. 2020).

Within uncertainties, the volume density of MORA 2 mm detected DSFGs is well aligned with constraints on $z \sim 3.5$ quiescent galaxies. Improved statistics from larger field 2 mm surveys will substantially reduce the uncertainty on the DSFG progenitor volume density at these epochs though there remains a substantial relative uncertainty in the quiescent population. This may be due, in large part, to discrepancies between what galaxies qualify as quiescent or star-forming. For example, the density of $z \sim 4$ quiescent systems from Straatman et al. (2014) presumes that all color-selected $z \sim 4$ sources not detected in the FIR with Herschel or Spitzer are, indeed, quiescent. If some additional fraction of their sample are star-forming, the MORA-inferred progenitor population would fall into alignment with observed volume densities. This is likely, as both Herschel and Spitzer have shallow sensitivity at $z = 4$, whereas a longer wavelength selection (starting at 850 µm) can more effectively segregate between quiescent and star-forming systems through FIR/millimeter\ detection, particularly at $z > 2$. Further, well-studied identified quiescent systems (Glazebrook et al. 2017) can reveal low levels of hidden star formation at the sensitivity of deep ALMA observations (Simpson et al. 2017; Schreiber et al. 2018b). Indeed, alternative estimates of the quiescent population number density (Muzzin et al. 2013; Davidzon et al. 2017; Girelli et al. 2019) are quite a bit lower ($\sim$few $\times 10^{-6} \text{ Mpc}^{-3}$), depending on the selection criteria.

Aside from the simple comparison of volume densities between MORA DSFGs and quiescent galaxies at high $z$, we can more closely consider the established stellar masses of MORA-selected galaxies, their potential for future star

formation, and their likely descendant population. In other words, what stellar masses do we expect MORA-descendant galaxies to have after eventually ceasing their star formation? And at what redshift are those descendants fully formed? Figure 13 shows the extrapolated masses and redshifts for two hypothetical quiescent descendant populations from the MORA sample of DSFGs. The first assumes the lion's share of the ISM mass (in the form of H$_2$) is converted into stars at the current (observed) SFR until that gas is fully depleted without replenishment.<fnref n="38" /> On average, the gas depletion time for the MORA sample is $\tau_{\text{depl}} = 110_{-60}^{+30}$ Myr. The second adopts an average starburst timescale of $\sim 100$ Myr (in line with measurement of other DSFG samples, as well as the MORA sample herein, e.g., Bothwell et al. 2013) and presumes star formation will continue at the observed high rate for that fixed period of time before ceasing. The redshifts of the descendant population of quiescent systems is offset from the measured MORA redshifts by the measured or adopted gas depletion time. Note that there is no systematic offset of one population of hypothetical descendants from the other; half of the sample has higher masses and lower redshifts calculated with one method.<fnref n="39" /> Depending on the existing stellar reservoir, we find that anywhere between $\approx 20\%–100\%$ of the eventual stellar mass in MORA descendants will form in the observed star-forming episode. The final stellar masses of these systems is then $\sim 2 \times 10^{11} M_\odot$, with the entire MORA 2 mm sample producing galaxies with stellar masses $> 4 \times 10^{10} M_\odot$; this is well aligned with the measured mass limits and redshift regimes of high-$z$ quiescent galaxy surveys.

### 6.3. Is the Emissivity Spectral Index Steep, and/or Does It Evolve?

The dust emissivity spectral index, or $\beta$, governs the frequency dependence of the emissivity of dust grains per unit mass, such that $\kappa_\nu = \kappa_0(\nu/\nu_0)^\beta$. The mass absorption coefficient, $\kappa_\nu$, traces the chemical and optical properties of dust grains. A lower value of $\beta$ observationally manifests in a shallower Rayleigh–Jeans falloff of the SED in the millimeter, while higher values result in steeper Rayleigh–Jeans falloffs.

Data on high-redshift galaxies—where the FIR/millimeter SED is spatially unresolved and poorly sampled—is of limited use in constraining the intrinsic dust emissivity spectral index, which relates to the fundamental dust composition and likely varies on the scale of molecular clouds. For example, Arendt et al. (2019) measure an average value of $\beta = 2.25$ for the central molecular zone in the Galactic plane of the Milky Way using new 2 mm GISMO observations; this is broadly in line with findings from Planck, which revealed $\beta \approx 1.6$ at high Galactic latitudes but increasingly steep values, $\beta = 1.8–2.0$, toward the inner Galactic plane ($|b| \ge 10^\circ$; Planck Collaboration et al. 2014). Despite our inability to constrain the underlying physical quantity $\beta$, measurement of the slope of the Rayleigh–Jeans falloff, as captured by the galaxy-integrated $\beta$ calculated herein, is useful. Elucidating our view of the

<footnote n="38">³⁸ While the conversion of gas to stars is far from 100% efficient (Evans et al. 2009), the timescale of individual star formation episodes is much shorter than the duration of DSFGs in their elevated SFR phase; thus such substantial conversions are reasonable to presume over long timescales (see also Walter et al. 2020).</footnote>

<footnote n="39">³⁹ The median redshift and mass of the first (fixed gas supply) are $\langle z \rangle = 3.31 \pm 0.02$ and $(1.9 \pm 0.2) \times 10^{11} M_\odot$, while for the second (fixed starburst timescale) they are $\langle z \rangle = 3.24 \pm 0.03$ and $(2.2^{+0.3}_{-0.2}) \times 10^{11} M_\odot$.</footnote>
<figure><img src="image_18.png" /><figcaption>**Figure 13.** Stellar mass with redshift for both the MORA sample of 2 mm selected DSFGs (blue stars) and the extrapolated descendants. Two hypothetical descendants are shown per MORA source: one with a fixed “starburst” duration of 100 Myr where star formation proceeds at the observed rate (light orange squares), and one with a fixed gas supply, such that the current star formation proceeds at the observed rate until the existing gas supply is depleted (dark orange circles). Both descendants for each source are connected with an orange line, and back to the MORA progenitor by a gray triangle. Gray dashed lines show the mass growth of halos of different rarities (spanning $10^{-3}-10^{-6}$ Mpc$^{-3}$; Behroozi et al. 2019) where halo masses have been scaled to stellar masses assuming a stellar mass to halo mass ratio of 2%. Note that the stellar mass of MORA-3 (AZTEC-2) at $z = 4.6$ is unknown and here assumed to be extremely low, for lack of other information. Overall we find that MORA descendant galaxies sit at or above $10^{11} M_{\odot}$ and are well established at those masses before $z = 2$.</figcaption></figure>

average galaxy-integrated $\beta$ is useful for relating galaxies’ dust SEDs to the (sub)millimeter surveys used to identify DSFGs in the first instance, simulating galaxies’ SEDs and SED fitting to high-$z$ unresolved sources.

Our measurement of $\langle\beta\rangle = 2.2^{+0.5}_{-0.4}$ within this sample is consistent with the often assumed fixed value of $\beta \equiv 1.8$, though it is also suggestively a bit steeper (compared to the prevailing theory based on ISM composition that would suggest values of $\beta = 1-2$; Draine 2011). The analysis of SCUBA-2 selected sources sitting inside the MORA mosaic (in Section 5.5) are similarly consistent with a range in $\beta$ that skews a bit higher than $\beta > 1.8$, as shown in Figure 11. How do we interpret these relatively steep values of $\beta$ relative to shallower values in integrated SEDs of well-constrained local universe dust (Dunne & Eales 2001; Clements et al. 2010)? While this could be suggestive of a higher proportion of large silicate grains in higher-redshift dusty galaxies, it could also mark different underlying ISM geometries.

Theoretically, low values (e.g., $\beta \approx 1$) correspond to dust comprised of small amorphous carbons with a mix of underlying dust temperatures from warm to cold; this can be easily understood from the SED that would result from coadding several blackbodies of different temperatures yet similar masses: the net unresolved SED would have a shallower Rayleigh–Jeans tail than a dust distribution of any one temperature. Indeed, coadding spatially distinct SEDs over the scale of an entire galaxy would always result in a *shallowing* of the Rayleigh–Jeans tail and not a steepening. One conjecture in the literature (e.g., Jin et al. 2019) is that the

steepening of the Rayleigh–Jeans tail may be due to CMB heating of the SED at high redshifts, which would cause a greater reduction in the observed flux density at longer wavelengths relative to shorter wavelengths (da Cunha et al. 2013). However, CMB heating is likely only a significant effect for galaxies at $z \gtrsim 5$ at very low dust temperatures ($T < 30$ K) where the CMB temperature is nonnegligible; furthermore, our SEDs do account for this CMB heating (however negligible), and our derived $\beta$ values still hold after accounting for this effect. While steeper values of $\beta$ can also indicate the presence of large crystalline silicate grains that reach overall lower temperatures (Draine & Lee 1984; Agladze et al. 1996; Meny et al. 2007), the degeneracies between spatial distribution of the ISM, dust temperature, opacity, and $\beta$, as well as the relative measurement uncertainties on our measurements from a dearth of data, prevent any direct constraints on dust composition. Still it is an interesting puzzle that our galaxy-integrated $\beta$ seems to skew high relative to that in the prevailing literature (see also Kato et al. 2018), that has, for lack of direct constraints, fixed the value to either $\beta = 1.5$ (Paradis et al. 2009) or $\beta = 1.8$ (Planck Collaboration et al. 2011) from measurements of the Milky Way’s ISM in both atomic and molecular states.

Our measurements of $\beta = 2.2^{+0.5}_{-0.4}$, along with other recent direct measurements of $\beta$ in high-$z$ DSFGs (see da Cunha et al. 2021, who measure $\beta = 1.9 \pm 0.4$ for the ALESS sample via dedicated 2 mm follow-up observations), suggest that a higher $\beta \approx 2$ is more appropriate than $\beta = 1.5$ or $\beta = 1.8$ in instances where a DSFGs’ Rayleigh–Jeans tail is not directly constrained, and a choice of $\beta$ is needed to model sources’ SEDs.

## 6.4. Potential Impact of Cosmic Variance on $z > 3$ DSFG Sample

Though our MORA Survey area is among the largest mosaics stitched together by ALMA to date, a chief concern of our analysis is that cosmic variance may impact our measurement of the ubiquity of DSFGs beyond $z > 3$. In other words, our relatively small area—small in comparison to the relative rarity of DSFGs themselves—may oversample or undersample the average density of the universe at any given epoch. Our parallel paper, Z21, tests the impact of cosmic variance using dust emission models applied to galaxies inside large volume simulations (see Z21, their Section 4.4, as well as Popping et al. 2020). In those simulations, Z21 determine that the MORA survey volume is robust in the measurement of the number counts of 2 mm sources out to $z \sim 7$.

Here, we extend the cosmic variance analysis of Z21 to analyze the cosmic variance specifically for the $z > 4$ high-redshift tail of the 2 mm observed population, and also incorporating predictions from the SHARK semi-empirical model (Lagos et al. 2020). The definition of cosmic variance as given in Moster et al. (2011) is

$$ \sigma_v^2 = \frac{\langle N^2 \rangle - \langle N \rangle^2 - \langle N \rangle}{\langle N \rangle^2}, \qquad (2) $$

where $N$ represents the number of galaxies detectable in a given survey solid angle above a redshift $z$. In other words, we generate Monte Carlo trial mock surveys with increasing area $\Omega$ within the Popping et al. (2020) and Lagos et al. (2020) light cones and assess how many sources in these light cones would be 2 mm detectable (at the rough sensitivity of the MORA survey) above a given redshift. At a fixed survey area, $\Omega$, we
<figure><img src="image_19.png" /><figcaption>**Figure 14.** Cosmic variance, $\sigma_v^2$, as a function of survey solid angle for 2 mm detectable galaxies above $S_{2\,mm} > 0.3$ mJy and a range of redshifts drawn from: the 108 deg$^2$ SHARK light cone (solid lines; Lagos et al. 2020) and the 7.7 deg$^2$ UNIVERSEMACHINE light cone (dotted-dashed lines; Popping et al. 2020). The gray shaded regions denote $\sigma_v^2 > 0.15$ (and $>0.20$), a threshold representing significant (severe) cosmic variance. We show cosmic variance for sources above redshifts $z=4$ (purple lines), $z=4.5$ (green lines), $z=5$ (orange lines), and $z=5.5$ (red lines). Sources at higher redshifts are exceedingly rare ($\sim 1$ deg$^{-2}$) in both SHARK and UNIVERSEMACHINE, such that direct estimates of cosmic variance $\sigma_v^2$ cannot be reliably estimated. The size of the MORA Survey is indicated with the vertical gray dashed line. We find that 2 mm sources are roughly a factor of $\sim 2$ times more rare in the UNIVERSE- MACHINE light cone relative to SHARK, and so the cosmic variance is appreciably higher. This work suggests a significant jump in survey area is needed—exceeding 0.2–0.3 deg$^{-2}$—to take an accurate census of DSFGs beyond $z>5.5$, while surveys of order a few tenths of a square degree are adequate at lower redshifts. It is likely that surveys exceeding 1 deg$^{-2}$ will be needed to sample $z>6$ DSFGs well.</figcaption></figure>

generate a distribution of *N* that satisfies $S_{2\,mm} > 0.3$ mJy and $z > 4$ (or other higher-redshift minima); from that distribution of *N*, we are able to easily compute $\langle N \rangle$ and $\langle N^2 \rangle$ from which we calculate $\sigma_v^2$.

Our results are shown in Figure 14 where surveys with $\sigma_v^2 < 0.15$ are less prone to the effects of cosmic variance. We find that, due to the relative rarity of dusty sources in both SHARK and UNIVERSEMACHINE, the MORA Survey size ($\sim 0.05$ deg$^2$) is prone to cosmic variance above $z > 4$ (according to UNIVERSEMACHINE) and $z > 5$ (according to SHARK). The UNIVERSEMACHINE generated light cone has roughly a factor of 2–3 times fewer high SFR (SFR $\gtrsim 100 M_\odot$ yr$^{-1}$) galaxies than does SHARK at these redshifts and so the estimated cosmic variance is higher. Both models imply that significantly larger areas are needed to survey the $z > 5.5$ universe, of order 0.2–0.3 deg$^2$ or greater to suppress cosmic variance below $\sigma_v^2 < 0.15$. Given the dearth of $z > 6$ DSFGs in both models, it is hard to characterize what survey area would be necessary to mitigate cosmic variance at those epochs, though naturally we may expect it to be $\gtrsim 1$ deg$^2$.

Note too that the uncertainties for the predicted model redshift distributions and number counts, shown in Figure 5, account for the same cosmic variance discussed here by sampling different regions of either light cone over the 184 arcmin$^2$ MORA Survey footprint. Those estimates account for the heterogeneous noise characteristics of the MORA

Survey, while the cosmic variance estimates in Figure 14 assume a uniform flux density threshold of $S_{2\,mm} > 0.3$ mJy.

## 6.5. Potential Degeneracies in Interpretation

This paper has so far argued that 2 mm continuum observations can provide great clarity on the earliest epoch of dust obscured galaxies at $z \gtrsim 4$, but how might our conclusions be impacted by astrophysical unknowns? The MORA Survey represents an important step in constraining the obscured universe in this early epoch, but further refinements will require data to break certain degeneracies, including unknowns about a possible evolving shape in the IRLF, galaxies’ evolving dust SEDs, and dust emissivities.

For example, there is substantial evidence to suggest that the faint-end slope of the UVLF evolves, such that it is much steeper at higher redshifts (Finkelstein 2016). The physical interpretation of this is that there are far more low-mass galaxies at early times, and with time the distribution of galaxies in the luminosity function is distributed across a larger dynamic range. Does the IRLF evolve similarly, differently, in the opposite sense, or not at all?

As presented in Figure 12, our data are consistent with an unevolving faint-end slope of the IRLF, though the constraints do not rule out such an evolution. Being able to measure such an evolution will require more painstakingly deep observations (as can be provided by ALMA with deep 1–2 mm surveys at $< 1$ mJy depth, like ASPECS; González-López et al. 2020) over significantly larger areas ($\gtrsim 0.2$ deg$^2$). Furthermore, source redshifts—likely most efficiently obtained through large-area blind searches for millimeter spectral lines—will be paramount to this measurement. At present, no large field of view millimeter spectrometer exists with the needed sensitivity. Once in hand, the combination of a well-measured UVLF and IRLF, as well as an understanding of the nature of their evolution, can directly inform models of early universe dust formation and attenuation.

The question of evolution in galaxies’ dust SEDs has come up frequently in recent literature (e.g., Liang et al. 2019; Ma et al. 2019) and may have a significant impact on the results discussed herein and in Z21. For example, if $z > 5$ galaxies contain, on average, substantially hotter dust than lower-redshift DSFGs as such works suggest (with a correspondingly lower dust mass per fixed $L_{IR}$), then it is likely that 2 mm dust continuum maps could miss upwards of half of all DSFGs at early epochs. In other words, if the average dust SED is hotter at high redshifts, then those SEDs peak at shorter rest-frame wavelengths and are likely to have lower flux densities on the Rayleigh–Jeans tail (2 mm observed) than would colder DSFGs. This would lead to a smaller fraction of DSFGs being detected at 2 mm, and it would be logical to conclude that our results have underestimated the total volume density of early universe DSFGs (a point raised as a hypothetical in Casey et al. 2018a).

If that is the case, then some tension may exist between the inferred volume density of all DSFGs and the measured volume density of lower-redshift quiescent galaxies (as discussed in Section 6.2). Nevertheless, with exception of very luminous lensed systems (e.g., Reuter et al. 2020) the spectral energy distributions of DSFGs at $z > 4$ are not well constrained as an aggregate population. Measuring them will provide invaluable insight into ISM physics in the first few Gyr as well as further constraints on the volume density measurements of the population themselves.
<table><caption>Table 3
Photometry of the >5σ 2 mm Sample</caption><thead><tr><th>NAME</th><th>POSITION<br>(α<sub>2000</sub>, δ<sub>2000</sub>)</th><th>S/N<sub>2 mm</sub></th><th>S<sub>2 mm</sub><br>(μJy)</th><th>H-band<br>(AB)</th><th>S<sub>3.6</sub><br>(nJy)</th><th>S<sub>24</sub><br>(μJy)</th><th>S<sub>100</sub><br>(mJy)</th><th>S<sub>160</sub><br>(mJy)</th><th>S<sub>250</sub><br>(mJy)</th><th>S<sub>350</sub><br>(mJy)</th><th>S<sub>450</sub><br>(mJy)</th><th>S<sub>500</sub><br>(mJy)</th><th>S<sub>850</sub><br>(mJy)</th><th>S<sub>1.2 mm</sub><br>(mJy)</th><th>S<sub>3 GHz</sub><br>(μJy)</th></tr></thead><tbody><tr><td>MORA-0</td><td>10:00:15.617,<br>+02:15:49.00</td><td>7.90</td><td>818 ± 103</td><td>—</td><td>1714 ± 40</td><td>(21 ± 27)</td><td>(-1.5 ± 1.7)</td><td>(0.6 ± 3.7)</td><td>(8.9 ± 5.8)</td><td>(18.7 ± 6.3)</td><td>(10.7 ± 5.8)</td><td>(20.3 ± 6.8)</td><td>13.56 ± 0.12<sup>†</sup></td><td>4.13 ± 0.18</td><td>19.5 ± 2.6</td></tr><tr><td>MORA-1</td><td>10:00:19.740,<br>+02:32:03.80</td><td>7.74</td><td>703 ± 91</td><td>22.68 ± 0.02</td><td>6958 ± 94</td><td>189 ± 13</td><td>(0.0 ± 1.7)</td><td>31.4 ± 3.7</td><td>37.6 ± 5.8</td><td>48.4 ± 6.3</td><td>29.2 ± 5.2</td><td>35.7 ± 6.8</td><td>8.67 ± 0.06<sup>†</sup></td><td>2.67 ± 0.09</td><td>74.8 ± 0.9<sup>*</sup></td></tr><tr><td>MORA-2</td><td>10:00:10.146,<br>+02:13:35.00</td><td>7.58</td><td>529 ± 70</td><td>24.24 ± 0.04</td><td>5223 ± 100</td><td>(18 ± 25)</td><td>(3.9 ± 1.7)</td><td>13.1 ± 3.8</td><td>15.4 ± 5.8</td><td>31.3 ± 6.3</td><td>(17.2 ± 9.4)</td><td>27.7 ± 6.8</td><td>8.82 ± 0.71</td><td>2.66 ± 0.05</td><td>30.1 ± 2.9</td></tr><tr><td>MORA-3</td><td>10:00:08.037,<br>+02:26:12.20</td><td>7.25</td><td>1220 ± 168</td><td>...</td><td>...</td><td>...</td><td>(0.3 ± 1.7)</td><td>(1.0 ± 3.7)</td><td>(16.8 ± 5.8)</td><td>30.8 ± 6.3</td><td>23.5 ± 4.8</td><td>29.1 ± 6.8</td><td>16.75 ± 0.15<sup>†</sup></td><td>4.62 ± 0.11</td><td>15.0 ± 2.4</td></tr><tr><td>MORA-4</td><td>10:00:26.359,<br>+02:15:28.00</td><td>6.68</td><td>557 ± 83</td><td>—</td><td>87 ± 29</td><td>(10 ± 18)</td><td>(0.5 ± 1.5)</td><td>(-0.6 ± 2.8)</td><td>(2.9 ± 5.8)</td><td>(2.9 ± 6.3)</td><td>(2.3 ± 5.8)</td><td>(4.9 ± 6.8)</td><td>5.908 ± 0.052<sup>†</sup></td><td>2.05 ± 0.11</td><td>(10.6 ± 4.1)</td></tr><tr><td>MORA-5</td><td>10:00:24.157,<br>+02:20:05.40</td><td>6.63</td><td>584 ± 88</td><td>—</td><td>406 ± 24</td><td>(-8 ± 26)</td><td>(-0.2 ± 1.7)</td><td>(-0.1 ± 3.7)</td><td>(3.8 ± 5.8)</td><td>(7.3 ± 6.3)</td><td>(10.7 ± 4.1)</td><td>(7.1 ± 6.8)</td><td>6.80 ± 0.53</td><td>2.30 ± 0.10</td><td>(10.1 ± 3.4)</td></tr><tr><td>MORA-6</td><td>10:00:28.723,<br>+02:32:03.40</td><td>6.37</td><td>615 ± 97</td><td>23.75 ± 0.03</td><td>4332 ± 30</td><td>198 ± 13</td><td>(-0.4 ± 1.7)</td><td>19.3 ± 3.3</td><td>34.7 ± 5.8</td><td>35.0 ± 6.3</td><td>40.6 ± 5.2</td><td>24.8 ± 6.8</td><td>12.57 ± 0.15<sup>†</sup></td><td>3.13 ± 0.10</td><td>45.5 ± 0.8<sup>*</sup></td></tr><tr><td>MORA-7</td><td>10:00:11.574,<br>+02:15:05.20</td><td>6.29</td><td>378 ± 60</td><td>23.52 ± 0.02</td><td>11,956 ± 80</td><td>175 ± 16</td><td>(0.6 ± 1.8)</td><td>(-0.3 ± 3.8)</td><td>(16.1 ± 5.8)</td><td>(17.0 ± 6.3)</td><td>20.7 ± 6.6</td><td>(9.5 ± 6.8)</td><td>5.31 ± 0.16<sup>†</sup></td><td>1.51 ± 0.13</td><td>12.2 ± 2.3</td></tr><tr><td>MORA-8</td><td>10:00:25.292,<br>+02:18:46.20</td><td>5.75</td><td>489 ± 85</td><td>23.21 ± 0.02</td><td>11,156 ± 79</td><td>204 ± 29</td><td>(-1.7 ± 1.9)</td><td>(0.6 ± 3.8)</td><td>(16.2 ± 5.8)</td><td>24.7 ± 6.3</td><td>15.2 ± 4.3</td><td>26.8 ± 6.8</td><td>7.12 ± 0.54</td><td>...</td><td>34.7 ± 2.9</td></tr><tr><td>MORA-9</td><td>10:00:17.298,<br>+02:27:15.80</td><td>5.55</td><td>379 ± 68</td><td>26.75 ± 0.30</td><td>469 ± 33</td><td>(1 ± 13)</td><td>(-0.5 ± 1.7)</td><td>(-0.3 ± 3.7)</td><td>(0.0 ± 5.8)</td><td>(0.0 ± 6.3)</td><td>(1.1 ± 4.1)</td><td>(0.0 ± 6.8)</td><td>2.59 ± 0.37<sup>†</sup></td><td>...</td><td>(4.3 ± 2.4)</td></tr><tr><td>MORA-10</td><td>10:00:16.578,<br>+02:26:38.00</td><td>5.21</td><td>405 ± 78</td><td>22.77 ± 0.01</td><td>10,766 ± 66</td><td>890 ± 17</td><td>5.8 ± 1.7</td><td>18.9 ± 4.0</td><td>(15.6 ± 5.8)</td><td>(17.4 ± 6.3)</td><td>20.4 ± 4.1</td><td>(8.5 ± 6.8)</td><td>6.07 ± 0.07<sup>†</sup></td><td>...</td><td>3212.0 ± 1.6<sup>*</sup></td></tr><tr><td>MORA-11</td><td>10:00:12.922,<br>+02:12:11.40</td><td>5.16</td><td>356 ± 69</td><td>22.64 ± 0.01</td><td>13,322 ± 199</td><td>(22 ± 44)</td><td>(0.5 ± 1.6)</td><td>(0.4 ± 3.1)</td><td>35.7 ± 15.8</td><td>40.0 ± 20.3</td><td>(7.6 ± 14.4)</td><td>41.9 ± 25.9</td><td>5.10 ± 0.78</td><td>1.49 ± 0.07</td><td>12.1 ± 2.4</td></tr><tr><td>MORA-12</td><td>10:00:04.713,<br>+02:29:55.20</td><td>5.03</td><td>950 ± 190</td><td>—</td><td>(1 ± 30)</td><td>(1 ± 12)</td><td>(0.7 ± 1.7)</td><td>(0.4 ± 3.7)</td><td>(13.8 ± 5.8)</td><td>(14.2 ± 6.3)</td><td>(1.2 ± 4.8)</td><td>(5.7 ± 6.8)</td><td>(0.22 ± 0.69)</td><td>...</td><td>(1.0 ± 2.3)</td></tr></tbody></table>

Note. Photometry enclosed in brackets represent <3σ detections. Ellipses (...) indicate an absence of data, while long dashes (—) represent a nondetection in the COSMOS2020 catalog. Note that MORA-3 has highly confused Spitzer imaging, blending MORA-3 with two lower-redshift foreground sources within the Spitzer beam at both 3.6 and 24 µm. A dagger (†) indicates that the flux density in the 850 µm column has been replaced by improved 870 µm photometry from ALMA rather than SCUBA-2. An asterisk (*) indicates that the 3 GHz flux density comes from the deeper COSMOS-XS maps of van der Vlugt et al. (2021) and Algera et al. (2020) rather than the shallower/wider coverage of Smolčić et al. (2017). Note that some of the sample have ALMA band 6 data from archival data sets (recorded here at 1.2 mm).
<table><caption>**Table 4**
Positions and Multiwavelength Counterparts for the Marginal $4 < \sigma < 52$ mm Sample</caption><thead><tr><th>NAME</th><th>POSITION<br>(αJ2000, δJ2000)</th><th>S/N<sub>2</sub> mm</th><th>S<sub>2</sub> mm<br>(μJy)</th><th>H-band<br>(AB)</th><th>S<sub>850</sub><br>(mJy)</th><th>S<sub>3</sub> GHz<br>(μJy)</th></tr></thead><tbody><tr><td>MORA-P03-0</td><td>10:00:42.629, +02:14:39.00</td><td>4.68</td><td>615 ± 131</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-P03-1</td><td>10:00:42.843, +02:14:59.80</td><td>4.51</td><td>493 ± 109</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-P03-2</td><td>10:00:43.643, +02:22:05.20</td><td>4.37</td><td>308 ± 70</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-P03-3</td><td>10:00:42.415, +02:27:18.60</td><td>4.28</td><td>703 ± 164</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-P03-4</td><td>10:00:46.019, +02:28:29.20</td><td>4.24</td><td>506 ± 119</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-P03-5</td><td>10:00:42.629, +02:18:29.20</td><td>4.19</td><td>551 ± 131</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-P03-6</td><td>10:00:42.362, +02:26:56.60</td><td>4.07</td><td>711 ± 174</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-P03-7</td><td>10:00:45.685, +02:24:26.40</td><td>4.07</td><td>378 ± 93</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-P03-8</td><td>10:00:45.818, +02:13:06.20</td><td>4.07</td><td>416 ± 102</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-P03-9</td><td>10:00:44.511, +02:15:49.20</td><td>4.04</td><td>256 ± 63</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-P03-10</td><td>10:00:42.936, +02:11:51.00</td><td>4.01</td><td>408 ± 102</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-13</td><td>10:00:24.398, +02:21:45.00</td><td>4.78</td><td>419 ± 88</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-14</td><td>10:00:26.119, +02:19:07.80</td><td>4.74</td><td>397 ± 84</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-15</td><td>10:00:25.118, +02:15:37.20</td><td>4.68</td><td>398 ± 85</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-16</td><td>10:00:25.505, +02:20:30.20</td><td>4.68</td><td>397 ± 85</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-17</td><td>10:00:27.761, +02:24:43.40</td><td>4.61</td><td>406 ± 88</td><td>26.74 ± 0.24</td><td>...</td><td>...</td></tr><tr><td>MORA-18</td><td>10:00:23.103, +02:12:27.40</td><td>4.59</td><td>401 ± 88</td><td>26.19 ± 0.14</td><td>...</td><td>...</td></tr><tr><td>MORA-19</td><td>10:00:12.000, +02:23:09.80</td><td>4.58</td><td>283 ± 62</td><td>25.00 ± 0.16</td><td>4.00 ± 0.57</td><td>259.0 ± 2.9</td></tr><tr><td>MORA-20</td><td>10:00:24.572, +02:32:03.00</td><td>4.58</td><td>404 ± 88</td><td>...</td><td>...</td><td>30.3 ± 3.9*</td></tr><tr><td>MORA-21</td><td>10:00:06.543, +02:21:31.60</td><td>4.54</td><td>414 ± 91</td><td>26.81 ± 0.50</td><td>...</td><td>...</td></tr><tr><td>MORA-22</td><td>10:00:12.960, +02:32:31.00</td><td>4.54</td><td>304 ± 67</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-23</td><td>10:00:12.575, +02:14:44.20</td><td>4.52</td><td>298 ± 66</td><td>21.30 ± 0.01</td><td>7.19 ± 0.64</td><td>...</td></tr><tr><td>MORA-24</td><td>10:00:07.837, +02:21:05.20</td><td>4.52</td><td>656 ± 145</td><td>24.88 ± 0.08</td><td>...</td><td>...</td></tr><tr><td>MORA-25</td><td>10:00:23.970, +02:17:49.80</td><td>4.52</td><td>398 ± 88</td><td>18.94 ± 0.01</td><td>10.48 ± 0.55</td><td>96.8 ± 2.8</td></tr><tr><td>MORA-26</td><td>10:00:22.716, +02:25:31.80</td><td>4.51</td><td>397 ± 88</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-27</td><td>10:00:22.636, +02:23:04.20</td><td>4.46</td><td>392 ± 88</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-28</td><td>10:00:24.945, +02:24:05.60</td><td>4.46</td><td>385 ± 86</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-29</td><td>10:00:23.650, +02:21:55.40</td><td>4.34</td><td>389 ± 89</td><td>23.21 ± 0.02</td><td>7.87 ± 0.55</td><td>256 ± 2.4</td></tr><tr><td>MORA-30</td><td>10:00:14.563, +02:16:45.00</td><td>4.34</td><td>422 ± 97</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-31</td><td>10:00:24.038, +02:29:48.20</td><td>4.32</td><td>383 ± 89</td><td>26.14 ± 0.16</td><td>6.12 ± 0.62</td><td>...</td></tr><tr><td>MORA-32</td><td>10:00:20.434, +02:32:57.80</td><td>4.29</td><td>397 ± 93</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-33</td><td>10:00:05.274, +02:23:58.40</td><td>4.29</td><td>529 ± 123</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-34</td><td>10:00:27.748, +02:26:45.40</td><td>4.29</td><td>377 ± 88</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-35</td><td>10:00:20.768, +02:13:16.20</td><td>4.28</td><td>382 ± 89</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-36</td><td>10:00:27.280, +02:20:16.80</td><td>4.28</td><td>364 ± 85</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-37</td><td>10:00:15.870, +02:24:46.00</td><td>4.27</td><td>424 ± 99</td><td>...</td><td>5.38 ± 0.58</td><td>...</td></tr><tr><td>MORA-38</td><td>10:00:18.219, +02:18:51.60</td><td>4.27</td><td>311 ± 73</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-39</td><td>10:00:25.438, +02:15:58.80</td><td>4.27</td><td>360 ± 84</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-40</td><td>10:00:05.863, +02:16:08.40</td><td>4.26</td><td>416 ± 98</td><td>...</td><td>...</td><td>6.4 ± 1.0*</td></tr><tr><td>MORA-41</td><td>10:00:05.875, +02:26:16.20</td><td>4.25</td><td>414 ± 97</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-42</td><td>10:00:22.783, +02:16:07.80</td><td>4.25</td><td>373 ± 88</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-43</td><td>10:00:30.736, +02:15:47.00</td><td>4.25</td><td>1140 ± 268</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-44</td><td>10:00:27.560, +02:21:58.20</td><td>4.25</td><td>368 ± 87</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-45</td><td>10:00:15.603, +02:25:50.40</td><td>4.24</td><td>458 ± 108</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-46</td><td>10:00:20.434, +02:28:03.00</td><td>4.23</td><td>389 ± 92</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-47</td><td>10:00:15.964, +02:20:27.80</td><td>4.22</td><td>403 ± 96</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-48</td><td>10:00:22.743, +02:32:41.00</td><td>4.21</td><td>373 ± 89</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-49</td><td>10:00:30.097, +02:31:21.40</td><td>4.20</td><td>646 ± 154</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-50</td><td>10:00:05.394, +02:29:44.20</td><td>4.20</td><td>487 ± 116</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-51</td><td>10:00:25.333, +02:33:32.20</td><td>4.19</td><td>369 ± 88</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-52</td><td>10:00:04.420, +02:28:00.80</td><td>4.19</td><td>1090 ± 261</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-53</td><td>10:00:27.321, +02:31:40.60</td><td>4.19</td><td>358 ± 86</td><td>...</td><td>5.95 ± 0.68</td><td>...</td></tr><tr><td>MORA-54</td><td>10:00:16.751, +02:23:59.40</td><td>4.18</td><td>313 ± 75</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-55</td><td>10:00:10.679, +02:21:41.60</td><td>4.17</td><td>258 ± 62</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-56</td><td>10:00:28.108, +02:27:10.40</td><td>4.15</td><td>377 ± 91</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-57</td><td>10:00:07.184, +02:14:58.20</td><td>4.14</td><td>427 ± 103</td><td>23.69 ± 0.03</td><td>...</td><td>...</td></tr><tr><td>MORA-58</td><td>10:00:06.823, +02:21:18.60</td><td>4.13</td><td>389 ± 94</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-59</td><td>10:00:29.109, +02:29:10.80</td><td>4.13</td><td>427 ± 103</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-60</td><td>10:00:05.637, +02:11:56.40</td><td>4.12</td><td>434 ± 105</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-61</td><td>10:00:20.114, +02:18:01.40</td><td>4.12</td><td>375 ± 91</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-62</td><td>10:00:30.684, +02:31:26.00</td><td>4.12</td><td>1030 ± 251</td><td>23.21 ± 0.01</td><td>...</td><td>...</td></tr><tr><td>MORA-63</td><td>10:00:27.614, +02:22:32.40</td><td>4.11</td><td>358 ± 87</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-64</td><td>10:00:23.170, +02:28:51.40</td><td>4.11</td><td>364 ± 89</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>
<table><caption>**Table 4**
(Continued)</caption><thead><tr><th>NAME</th><th>POSITION<br>(α<sub>J2000</sub>, δ<sub>J2000</sub>)</th><th>S/N<sub>2 mm</sub></th><th>S<sub>2 mm</sub><br>(μJy)</th><th>H-band<br>(AB)</th><th>S<sub>850</sub><br>(mJy)</th><th>S<sub>3 GHz</sub><br>(μJy)</th></tr></thead><tbody><tr><td>MORA-65</td><td>10:00:13.762, +02:14:00.60</td><td>4.10</td><td>320 ± 78</td><td>...</td><td>2.51 ± 0.66</td><td>...</td></tr><tr><td>MORA-66</td><td>10:00:28.214, +02:19:53.20</td><td>4.09</td><td>380 ± 93</td><td>...</td><td>3.10 ± 0.53</td><td>...</td></tr><tr><td>MORA-67</td><td>10:00:28.923, +02:32:18.60</td><td>4.09</td><td>407 ± 99</td><td>23.38 ± 0.02</td><td>...</td><td>...</td></tr><tr><td>MORA-68</td><td>10:00:16.978, +02:28:13.00</td><td>4.09</td><td>290 ± 71</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-69</td><td>10:00:25.332, +02:29:31.00</td><td>4.09</td><td>352 ± 86</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-70</td><td>10:00:27.867, +02:19:51.80</td><td>4.08</td><td>366 ± 90</td><td>...</td><td>3.10 ± 0.53</td><td>...</td></tr><tr><td>MORA-71</td><td>10:00:28.175, +02:31:48.20</td><td>4.08</td><td>372 ± 91</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-72</td><td>10:00:07.396, +02:23:44.60</td><td>4.08</td><td>460 ± 113</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-73</td><td>10:00:12.934, +02:23:45.60</td><td>4.08</td><td>275 ± 68</td><td>...</td><td>2.64 ± 0.57</td><td>93.7 ± 2.6</td></tr><tr><td>MORA-74</td><td>10:00:14.002, +02:23:33.40</td><td>4.08</td><td>336 ± 82</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-75</td><td>10:00:06.104, +02:10:50.40</td><td>4.07</td><td>754 ± 185</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-76</td><td>10:00:28.963, +02:30:13.60</td><td>4.07</td><td>408 ± 100</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-77</td><td>10:00:05.410, +02:12:16.40</td><td>4.06</td><td>470 ± 116</td><td>19.15 ± 0.01</td><td>...</td><td>112.0 ± 2.8</td></tr><tr><td>MORA-78</td><td>10:00:22.997, +02:28:09.40</td><td>4.06</td><td>359 ± 88</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-79</td><td>10:00:09.370, +02:32:14.40</td><td>4.06</td><td>443 ± 109</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-80</td><td>10:00:04.500, +02:28:58.60</td><td>4.06</td><td>963 ± 237</td><td>...</td><td>...</td><td>25.7 ± 3.3*</td></tr><tr><td>MORA-81</td><td>10:00:26.492, +02:15:27.00</td><td>4.05</td><td>338 ± 83</td><td>...</td><td>5.90 ± 0.60</td><td>...</td></tr><tr><td>MORA-82</td><td>10:00:23.170, +02:17:38.60</td><td>4.05</td><td>357 ± 88</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-83</td><td>10:00:22.410, +02:30:26.20</td><td>4.04</td><td>354 ± 88</td><td>22.72 ± 0.01</td><td>3.88 ± 0.64</td><td>42.0 ± 3.5</td></tr><tr><td>MORA-84</td><td>10:00:08.983, +02:32:08.00</td><td>4.04</td><td>617 ± 152</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-85</td><td>10:00:21.902, +02:22:46.00</td><td>4.04</td><td>355 ± 88</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-86</td><td>10:00:05.982, +02:25:12.40</td><td>4.04</td><td>384 ± 95</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-87</td><td>10:00:07.368, +02:33:58.00</td><td>4.03</td><td>655 ± 162</td><td>26.73 ± 0.48</td><td>...</td><td>...</td></tr><tr><td>MORA-88</td><td>10:00:30.122, +02:20:38.60</td><td>4.03</td><td>637 ± 158</td><td>...</td><td>...</td><td>57.1 ± 0.9*</td></tr><tr><td>MORA-89</td><td>10:00:29.629, +02:24:10.20</td><td>4.03</td><td>490 ± 122</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-90</td><td>10:00:20.341, +02:16:49.60</td><td>4.03</td><td>366 ± 91</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-91</td><td>10:00:08.129, +02:33:01.60</td><td>4.03</td><td>728 ± 181</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-92</td><td>10:00:29.309, +02:25:54.20</td><td>4.02</td><td>438 ± 109</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-93</td><td>10:00:26.773, +02:20:44.00</td><td>4.02</td><td>336 ± 84</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-94</td><td>10:00:07.571, +02:16:31.40</td><td>4.01</td><td>492 ± 122</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-95</td><td>10:00:25.745, +02:18:23.00</td><td>4.01</td><td>338 ± 84</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-96</td><td>10:00:14.482, +02:28:48.40</td><td>4.01</td><td>385 ± 96</td><td>20.81 ± 0.01</td><td>...</td><td>...</td></tr><tr><td>MORA-97</td><td>10:00:06.729, +02:25:36.00</td><td>4.01</td><td>373 ± 93</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-98</td><td>10:00:09.412, +02:22:02.20</td><td>4.00</td><td>419 ± 105</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MORA-99</td><td>10:00:19.660, +02:31:04.60</td><td>4.00</td><td>360 ± 90</td><td>23.17 ± 0.02</td><td>...</td><td>...</td></tr></tbody></table>

Note. Marginal 2 mm sources and their multiwavelength counterparts, identified in the COSMOS2020 photometric catalog within 1" of the 2 mm position (with H-band magnitude from Ultra-VISTA; Weaver et al. 2021), at 850 µm with SCUBA-2 within 8" of the 2 mm position (Simpson et al. 2019), or at 3 GHz within 1" of the 2 mm position (Smolčić et al. 2017). Sources marked with asterisks (*) have 3 GHz flux density measurements from the deeper COSMOS-XS survey (Algera et al. 2020; van der Vlugt et al. 2021).

Variation in the dust emissivity index may also impact our volume density constraints in a similar fashion to variation in the aggregate dust temperatures of early DSFGs. With a steeper $\beta$, flux densities on the Rayleigh-Jeans tail drop; for example, the 2 mm flux density of a galaxy at $z=4$ would be three times lower for a very steep $\beta=3.5$ than for a shallower value $\beta=1.5$, presuming a fixed value of $L_{IR}$ and $\lambda_{peak}$. As discussed in Section 6.3, our sample does show moderate evidence for steeper emissivity spectral indices across the sample. If DSFGs on a whole were to exhibit steeper Rayleigh-Jeans tails it may imply that our volume density measurements here, and in Z21, are underestimated.

Lastly, it is crucial to point out that many MORA source redshifts are not yet well constrained. Given our small sample size, even one catastrophically misidentified redshift would shift our integrated SFRD results by ~10%, and thus has the potential to alter our perception of the IRLF's evolution beyond $z>3$.

## 6.6. Future Surveys of $z \gtrsim 4$ Obscured Sources

Future 2 mm surveys of large areas of sky will be crucial to overcome limitations of small sample sizes and cosmic variance, particularly in pursuit of rare, dusty starbursts at $z>4$. In particular, testing the divergent redshift distribution predictions for bright DSFGs ($S_{2 mm} > 1$ mJy; see Figure 6) in the simulations will be an essential next step. While the SHARK and UNIVERSEMACHINE simulations predict that the brightest 2 mm sources have a relatively low median redshift ($\langle z \rangle \sim 2.5$), the SIDES, Casey et al. and Zavala et al. models suggest a higher average redshift ($\langle z \rangle \gtrsim 3.5$). This difference originates with the fundamental problem of how massive, bright galaxies assemble *en masse* so quickly after the Big Bang. Sampling these bright sources in statistically large samples well will truly require a different scale of millimeter continuum survey, of order $\gtrsim 5$ deg$^2$, as such bright high-$z$ systems are predicted to be extremely rare in the models (e.g., SHARK predicts five $z>4$
sources per square degree above 1 mJy at 2 mm, and only one such source above $z > 5$ per square degree). While large-field searches of Herschel SPIRE imaging for 500 µm peaking “red” sources have been fruitful (Ivison et al. 2016; Duivenvoorden et al. 2018; Yan et al. 2020), the population is quite rare (~$10^{-7}$ Mpc$^{-3}$) relative to the MORA-detected DSFGs presented herein, due to the relative shallow sensitivity of Herschel at high redshifts. Instruments like GISMO, which has pioneered much of the 2 mm blank-field mapping work to date, and TolTEC, which in the future has the potential to expand 2 mm sky coverage by orders of magnitude, will play essential roles in pushing this frontier. Similarly, the South Pole Telescope survey has already demonstrated the effectiveness of >1 mm surveys in recovering an intrinsically high-z population (see Vieira et al. 2013; Reuter et al. 2020), and the next generation SPT survey (“SPT3G”) will push the depth such that unlensed sources will be detected at a high significance.

One complementary observing strategy to large surveys would be dedicated 2–3 mm continuum follow-up of DSFGs already identified at shorter wavelengths, allowing for a more efficient redshift selection of the highest-z sources using millimeter colors. Every source in the MORA Survey (save one that we attribute to a false noise peak) is detected above $3.5\sigma$ significance at 850 µm, and many at other (sub)millimeter wavelengths from existing single-dish surveys. As Figure 11 shows, the MORA sources with lower $S_{870}/S_{2000}$ ratios sit at the highest redshifts as would be expected for sources whose 850 µm emission is closer to the peak of dust emission. Though strategic 2 mm or 3 mm follow-up of 850 µm or 1 mm selected sources is more biased than the blank-field survey strategy, it is far more observationally efficient for samples that have already been identified, requiring only a few minutes of on-source time with ALMA. It is an effective, and observationally cheap filter through which needles ($z > 4$ DSFGs) can be spotted in the haystack (dominated by $1 < z < 3$ DSFGs), complementing 2 mm blank-field surveys like MORA.

Lastly, extensions of MORA stand to make significant progress; over the next year, MORA is approved for observations to expand coverage over $0.2\,\text{deg}^2$ (within the COSMOS-Web survey footprint) to the same depth of our current mosaics. These data will likely lead to the discovery of ~20 DSFGs at $z > 4$ and enable more detailed studies of DSFG clustering.

## 7. Conclusions

We have presented the MORA Survey, the first 2 mm extragalactic blank-field map with ALMA covering 184 arcmin$^2$. Our accompanying paper, Zavala et al. (2021), presents the MORA Survey number counts and a number count constrained obscured SFRD measurement out to $z \sim 7$, while this paper focuses on what is known about the individual galaxies that have been detected. We detect 13 sources above >$5\sigma$ significance, one of which may have significant contribution from radio synchrotron emission and another we believe to be a false-positive noise peak, leaving a total of 11 robust detections dominated by thermal dust emission. The redshifts of the sources span $z = 2.2–5.9$ with a median redshift $\langle z \rangle = 3.6_{-0.3}^{+0.4}$. Sources span IR luminosities $10^{12}$–few $\times 10^{13}$ L$_\odot$, ISM masses few $\times 10^{10}$–few $\times 10^{11}$ M$_\odot$, and sit both above and embedded within the galaxy SFR–stellar mass relation.

Our results paint a plausible picture of the buildup of dust obscured galaxies at the earliest epochs ($z > 4$): overall, DSFGs at this epoch are rare and contribute a sizable, though nondominant, fraction toward cosmic star formation. These DSFGs are the most likely progenitors of the universe’s first massive, quiescent galaxies. Our primary conclusions are as follows:

1. ALMA Band 4 (2 mm) is an effective selection wavelength for high-z DSFGs; here, we find $77\% \pm 11\%$ of our 2 mm selected sample lies above $z > 3$ and $38\% \pm 12\%$ above $z > 4$. This contrasts with more traditional methods of DSFG selection at 850 µm or 1 mm, for which the proportion of sources at these redshifts is much smaller (e.g., 6% of 850 µm selected sources sit at $z > 4$; Dudzevičiūtė et al. 2020).

2. “OIR-dark” sources (i.e., sources lacking counterparts in very deep optical and near-infrared data, $<1.6\,\mu$m) make up a sizable fraction of 2 mm selected sources (4/11 ~ 36%). Our accompanying paper, Manning et al. (2021) further explores their physical characteristics.

3. Several semi-empirical and empirical models from the literature accurately reproduce the redshift distribution and volume density of 2 mm detected galaxies. However, several details are not yet resolved by data, namely the epoch when the first DSFGs turn on, and how common such systems are at the highest redshifts, $z \gtrsim 5.5$. Different cosmological models infer that MORA may be impacted by cosmic variance above $z > 5$, with the results highly dependent on dust radiation prescriptions in early universe galaxies. Either way, larger area 2 mm continuum surveys will be essential to gather sufficient statistics to constrain the volume density of DSFGs at the highest-redshift regimes.

4. Dusty galaxies with star formation rates in excess of $300\,\text{M}_\odot\,\text{yr}^{-1}$ contribute ~30% to SFRD between $3 < z < 6$. Furthermore, at these epochs, the IRLF is dominated by these sources, with little contribution from galaxies with lower levels of obscured star formation ($\lesssim 100\,\text{M}_\odot\,\text{yr}^{-1}$). This implies that the vast majority of cosmic dust lives in such ultraluminous galaxies at these epochs, with relatively little to no dust impacting the bolometric output of less luminous, less massive systems.

5. DSFGs at $z > 3$, like those detected in the MORA Survey, are often claimed to be the progenitors of the universe’s earliest quiescent systems; we find that the number density of quiescent galaxies from the literature agrees well with the measured volume density of DSFGs measured at $z > 3$ (~$10^{-5}$ Mpc$^{-3}$). Furthermore, the extrapolated characteristics of MORA DSFGs’ descendants are also well aligned with detected quiescent populations, all well established with stellar masses >$10^{11}\,\text{M}_\odot$ above $z > 2$.

6. Our MORA data hint at a possible higher value for the galaxy-integrated dust emissivity index, $\beta$, than is often assumed in the literature: the sample has a measured $\langle \beta \rangle = 2.2_{-0.4}^{+0.5}$ versus $\beta = 1.8$. We recommend future works modeling SEDs of high-z DSFGs adopt a value of $\beta \approx 2$ as needed in lieu of direct constraints.

The final frontier of this work would lead to the detection of the universe’s first dusty galaxies, the home of the universe’s first substantial dust reservoirs. Do such systems exist beyond $z > 7$?
z > 8? While MORA has taken a tantalizing first step in completing the census of high-z DSFGs, we know higher-redshift DSFGs (specifically, three with $SFR \gtrsim 100 M_{\odot}yr^{-1}$) do exist than are found in this work (e.g., Cooray et al. 2014; Strandet et al. 2017; Marrone et al. 2018; Zavala et al. 2018b). Yet, their rarity—itself not yet well constrained—implies that finding more requires large-area millimeter surveys sufficiently sensitive to detect unlensted galaxies, similar conceptually to the all-sky searches for the highest-redshift quasars found by, e.g., the Sloan Digital Sky Survey (Fan et al. 2001, 2003). And though rare, the earliest dusty starbursts yet to be found serve as unique laboratories for our understanding of the assembly of the first galaxies and, fundamentally, the formation of dust, on which much of what we know about the modern-day universe relies.

We thank the anonymous reviewer for helpful suggestions that greatly improved the manuscript. This paper makes use of the ALMA program ADS/JAO.ALMA #2018.1.00231. 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. C.M. thanks the National Science Foundation for support through grants AST-1714528, AST-1814034, and AST-2009577 and additionally the University of Texas at Austin College of Natural Sciences for support. In addition, C.M.C. acknowledges support from the Research Corporation for Science Advancement from a 2019 Cottrell Scholar Award sponsored by IF/THEN, an initiative of Lyda Hill Philanthropies. The Flatiron Institute is supported by the Simons Foundation. M.T. acknowledges the support from grant PRIN MIUR 2017. E.T. acknowledges support from CATA-Basal AFB-170002, FONDECYT Regular grant 1190818, ANID Anillo ACT172033 and Millennium Nucleus NCN19_058 (TITANs). M.A. acknowledges support from FONDECYT grant 1211951, “CONICYT+PCI+INSTITUTO MAX PLANCK DE ASTRONOMIA MPG 190030” and “CONICYT+PCI+REDES 190194.” The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant No. 140. S.T. and J.W. acknowledge support from the European Research Council (ERC) Consolidator Grant funding scheme (project ConTExt, grant No. 648179). G.E.M. acknowledges the Villum Fonden research grant 13160 Gas to stars, stars to dust: tracing star formation across cosmic time, grant 37440, “The Hidden Cosmos,” and the Cosmic Dawn Center. A.W.S.M. acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC). K.I.C. acknowledges funding from the European Research Council through the award of the Consolidator Grant ID 681627-BUILDUP.

## ORCID iDs

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

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

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

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

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

Jaclyn B. Champagne https://orcid.org/0000-0002-6184-9097

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

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

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

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

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

Vasily Kokorev https://orcid.org/0000-0002-5588-9156

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

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

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

Ikki Mitsuhashi https://orcid.org/0000-0001-7300-9450

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

Justin Spilker https://orcid.org/0000-0003-3256-5615

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

Margherita Talia https://orcid.org/0000-0003-4352-2063

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

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

John R. Weaver https://orcid.org/0000-0003-1614-196X

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

## References

Agladze, N., Sievers, A., Jones, S., Burlitch, J., & Beckwith, S. 1996, ApJ, 462, 1026

Algera, H. S. B., van der Vlugt, D., Hodge, J. A., et al. 2020, ApJ, 903, 139

Aravena, M., Decarli, R., Walter, F., et al. 2016, ApJ, 833, 71

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

Arendt, R. G., Staguhn, J., Dwek, E., et al. 2019, ApJ, 885, 71

Aretxaga, I., Wilson, G., Aguilar, E., et al. 2011, MNRAS, 415, 3831

Arnouts, S., Moscardini, L., Vanzella, E., et al. 2002, MNRAS, 329, 355

Ashby, M., Willner, S., Fazio, G., et al. 2015, ApJS, 218, 33

Bakx, T. J., Eales, S., Negrello, M., et al. 2018, MNRAS, 473, 1751

Battisti, A., da Cunha, E., Grasha, K., et al. 2019, ApJ, 882, 61

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

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Bertoldi, F., Carilli, C., Aravena, M., et al. 2007, ApJS, 172, 132

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

Bethermin, M., Wu, H.-Y., Lagache, E., et al. 2017, A&A, 607, A89

Biggs, A., & Ivison, R. 2008, MNRAS, 385, 893

Blain, A., Smail, I., Ivison, R., Kneib, J., & Frayer, D. 2002, PhR, 369, 111

Bothwell, M., Aguirre, J., Chapman, S., et al. 2013, ApJ, 779, 67

Bouwens, R., González-López, J., Aravena, M., et al. 2020, ApJ, 902, 112

Brammer, G., van Dokkum, P., & Coppi, P. 2008, ApJ, 686, 1503

Brisbin, D., Miettinen, O., Aravena, M., et al. 2017, A&A, 608, A15

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Capak, P., Aussel, H., Ajiki, M., et al. 2007, ApJS, 172, 99

Capak, P. L., Scoville, N., Sanders, D., et al. 2010, AAS Meeting, 215, 410.05

Caputi, K., Deshmukh, S., Ashby, M., et al. 2017, ApJ, 849, 45

Caputi, K., Caminha, G., Fujimoto, S., et al. 2021, ApJ, 908, 146

Casey, C. 2016, ApJ, 824, 36

Casey, C., Hodge, J., Zavala, J., et al. 2018b, ApJ, 862, 78

Casey, C., Narayanan, D., & Cooray, A. 2014, PhR, 541, 45

Casey, C., Cooray, A., Capak, P., et al. 2015, ApJL, 808, L33

Casey, C., Cooray, A., Killi, M., et al. 2017, ApJ, 840, 101

Casey, C., Zavala, J., Spilker, J., et al. 2018a, ApJ, 862, 77

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

Casey, C. M. 2020, ApJ, 900, 68

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

Casey, C. M., Aravena, J. A., Aravena, M., et al. 2019, ApJ, 887, 55

Chapman, S., Blain, A., Smail, I., & Ivison, R. 2005, ApJ, 622, 772

Chapman, S., Helou, G., Lewis, G., & Dale, D. 2003, ApJ, 588, 186

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

Chiang, Y.-K., Overzier, R., Gebhardt, K., et al. 2015, ApJ, 808, 37

Civano, F., Elvis, M., Brusa, M., et al. 2012, ApJS, 201, 30

Civano, F., Marchesi, S., Comastri, A., et al. 2016, ApJ, 819, 62

Clements, D., Dunne, L., & Eales, S. 2010, MNRAS, 403, 274
Cochrane, R., Hayward, C., Angles-Alcazar, D., et al. 2019, MNRAS, 488, 1779

Conley, A., Cooray, A., Viera, J., et al. 2011, ApJL, 732, L35

Cooray, A., Calanog, J., Wardlow, J., et al. 2014, ApJ, 790, 40

Cortzen, I., Magdis, G. E., Valentino, F., et al. 2020, A&A, 634, L14

Cucciati, O., Lemaux, B., Zamorani, G., et al. 2018, A&A, 619, 49

da Cunha, E., Charlot, S., & Elbaz, D. 2008, MNRAS, 388, 1595

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

da Cunha, E., Walter, F., Smail, I., et al. 2015, ApJ, 806, 110

da Cunha, E., Hodge, J., Casey, C., et al. 2021, ApJ, 919, 30

Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827

Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70

Delhaize, J., Smolčić, V., Delvecchio, I., et al. 2017, A&A, 602, A4

Diener, C., Lilly, S., Ledoux, C., et al. 2015, ApJ, 802, 31

Donevski, D., Buat, V., Boone, F., et al. 2018, A&A, 614, A33

Dowell, C., Conley, A., Glenn, J., et al. 2014, ApJ, 780, 75

Draine, B. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton, N.J.: Princeton Univ. Press)

Draine, B., & Lee, H. 1984, ApJ, 285, 89

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

Duivenvoorden, S., Oliver, S., Scudder, J., et al. 2018, MNRAS, 477, 1099

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

Dunne, L., & Eales, S. 2001, MNRAS, 327, 697

Evans, N., II, Dunham, M., Jørgensen, J., et al. 2009, ApJS, 181, 321

Fan, X., Strauss, M., Schneider, D., et al. 2001, AJ, 121, 54

Fan, X., Strauss, M., Schneider, D., et al. 2003, AJ, 125, 1649

Fanaroff, B., & Riley, J. 1974, MNRAS, 167, 31P

Finkelstein, S. 2016, PASA, 33, e037

Finkelstein, S., Ryan, R., Jr., Papovich, C., et al. 2015, ApJ, 810, 71

Fixsen, D., Dwek, E., Mather, J., Bennett, C., & Shafer, R. 1998, ApJ, 508, 123

Forrest, B., Marsan, Z. C., Annunziatella, M., et al. 2020, ApJ, 903, 47

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

Franco, M., Elbaz, D., Zhou, L., et al. 2020, A&A, 643, A30

Fudamoto, Y., Oesch, P. A., Schouws, S., et al. 2021, Natur, 597, 489

Fujimoto, S., Ouchi, M., Shibuya, T., & Nagai, H. 2017, ApJ, 850, 83

Geach, J., Dunlop, J., Halpern, M., et al. 2017, MNRAS, 465, 1789

Girelli, G., Bolzonella, M., & Cimatti, A. 2019, A&A, 632, A80

Glazebrook, K., Schreiber, C., Labbé, I., et al. 2017, Natur, 544, 71

Gómez-Guijarro, C., Toft, S., Karim, A., et al. 2018, ApJ, 856, 121

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

Grogin, N., Kocevski, D., Faber, S., et al. 2011, ApJS, 197, 35

Gruppioni, C., Béthermin, M., Loiacono, F., et al. 2020, A&A, 643, A8

Harrison, I., & Hotchkiss, S. 2013, Journal of Cosmology and Astro-Particle Physics, 7, 022

Hatsukade, B., Ohta, K., Seko, A., Yabe, K., & Akiyama, M. 2013, ApJL, 769, L27

Hatsukade, B., Kohno, K., Umehata, H., et al. 2016, PASJ, 68, 36

Hayward, C. C., Behroozi, P. S., Somerville, R. S., et al. 2013, MNRAS, 434, 2572

Hayward, C. C., Kereš, D., Jonsson, P., et al. 2011, ApJ, 743, 159

Hayward, C. C., Sparre, M., Chapman, S. C., et al. 2021, MNRAS, 502, 2922

Hodge, J., Karim, A., Smail, I., et al. 2013, ApJ, 768, 91

Hodge, J., Swinbank, A., Simpson, J., et al. 2016, ApJ, 833, 103

Ilbert, O., Arnouts, S., McCracken, H. J., et al. 2006, A&A, 457, 841

Ilbert, O., McCracken, H., Le Fèvre, O., et al. 2011, A&A, 556, A55

Ivison, R., Papadopoulos, P., Smail, I., et al. 2011, MNRAS, 412, 1913

Ivison, R., Alexander, D., Biggs, A., et al. 2010, MNRAS, 402, 245

Ivison, R., Lewis, A., Weiss, A., et al. 2016, ApJ, 832, 78

Jiménez-Andrade, V., Zavala, J., Magnelli, B., et al. 2020, ApJ, 890, 171

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

Jonsson, P. 2006, MNRAS, 372, 2

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

Kato, Y., Matsuda, Y., Iono, D., et al. 2018, PASJ, 70, L6

Kennicutt, R., & Evans, N. 2012, ARA&A, 50, 531

Khusanova, Y., Bethermin, M., Le Fèvre, O., et al. 2021, A&A, 649, A152

Klypin, A., Yepes, G., Gottlöber, S., Prada, F., & Hess, S. 2016, MNRAS, 457, 4340

Koekemoer, A., Aussel, H., Calzetti, D., et al. 2007, ApJS, 172, 196

Koekemoer, A., Faber, S., Ferguson, H., et al. 2011, ApJS, 197, 36

Koprowski, M., Dunlop, J., Michałowski, M., et al. 2017, MNRAS, 471, 4155

Kriek, M., van Dokkum, P., Labbé, I., et al. 2009, ApJ, 700, 221

Kroupa, P., & Weidner, C. 2003, ApJ, 598, 1076

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

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

Lagos, C. d. P., Robotham, A. S., Trayford, J. W., et al. 2019, MNRAS, 489, 4196

Laigle, C., McCracken, H., Ilbert, O., et al. 2016, ApJS, 224, 24

Le Floc'h, E., Papovich, C., Dole, H., et al. 2005, ApJ, 632, 169

Lee, N., Sanders, D., Casey, C., et al. 2013, ApJ, 778, 131

Liang, L., Feldmann, R., Kereš, D., et al. 2019, MNRAS, 489, 1397

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

Loiacono, F., Decarli, R., Gruppioni, C., et al. 2021, A&A, 646, A76

Lovell, C. C., Geach, J. E., Davé, R., Narayanan, D., & Li, Q. 2020, arXiv:2006.15156

Lutz, D., Poglitsch, A., Altieri, B., et al. 2011, A&A, 532, A90

Ma, X., Hayward, C., et al. 2019, MNRAS, 487, 1844

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

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

Maniyar, A., Béthermin, M., & Lagache, G. 2018, A&A, 614, A39

Maniyar, A., Béthermin, M., & Lagache, G. 2021, A&A, 645, A40

Manning, S., Casey, C., Zavala, J., et al. 2021, arXiv:2111.02428

Marrone, D., Spilker, J., Hayward, C., et al. 2018, Natur, 553, 51

Marsan, Z., Marchesini, D., Brammer, G., et al. 2017, ApJ, 842, 21

Marsan, Z. C., Muzzin, A., Marchesini, D., et al. 2020, arXiv:2010.04725

McKinney, J., Hayward, C. C., Rosenthal, L. J., et al. 2021, ApJ, 921, 55

Meny, C., Gromov, V., Boudet, N., et al. 2007, A&A, 468, 171

Merlin, E., Fortuni, F., Torelli, M., et al. 2019, MNRAS, 490, 3309

Michałowski, M. J., Dunlop, J., Koprowski, M., et al. 2017, MNRAS, 469, 492

Miettinen, O., Delvecchio, I., Smolčić, V., et al. 2017, A&A, 606, A17

Mitsuhashi, I., Matsuda, Y., Smail, I., et al. 2021, ApJ, 907, 122

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

Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18

Narayanan, D., Turk, M. J., Robitaille, T., et al. 2021, ApJS, 252, 12

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

Oliver, S., Bock, J., Altieri, B., et al. 2012, MNRAS, 424, 1614

Paradis, D., Bernard, J.-P., & Mény, C. 2009, A&A, 506, 745

Planck Collaboration, Abergel, A., Ade, P., et al. 2011, A&A, 536, A25

Planck Collaboration, Abergel, A., Ade, P., et al. 2014, A&A, 571, A11

Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6

Pope, A., Chary, R., Alexander, D., et al. 2008, ApJ, 675, 1171

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

Rémy-Ruyer, A., Madden, S., Galliano, F., et al. 2014, A&A, 563, A31

Reuter, C., Vieira, J., Spilker, J., et al. 2020, arXiv:2006.14060

Robotham, A., Bellstedt, S., Lagos, C. d. P., et al. 2020, MNRAS, 495, 905

Rodríguez-Puebla, A., Behroozi, P., Primack, J., et al. 2016, MNRAS, 462, 893

Roseboom, I., Ivison, R., Greve, T., et al. 2012, MNRAS, 419, 2758

Roseboom, I., Dunlop, J., Cirasuolo, M., et al. 2013, MNRAS, 436, 430

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

Santini, P., Castellano, M., Merlin, E., et al. 2021, A&A, 652, A30

Sargent, M., Béthermin, M., Daddi, E., & Elbaz, D. 2012, ApJL, 747, L31

Schinnerer, E., Smolčić, V., Carilli, C., et al. 2007, ApJS, 172, 46

Schreiber, C., Glazebrook, K., Nanayakkara, T., et al. 2018a, A&A, 618, A85

Schreiber, C., Labbé, I., Glazebrook, K., et al. 2018b, A&A, 611, A22

Scott, K., Austermann, J., Perera, T., et al. 2008, MNRAS, 385, 2225

Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1

Scoville, N., Sheth, K., Aussel, H., et al. 2016, ApJ, 820, 83

Sherman, S., Jogee, S., Florez, J., et al. 2020, MNRAS, 499, 4239

Simpson, J., Smail, I., Swinbank, A., et al. 2015, ApJS, 807, 112

Simpson, J., Smail, I., Swinbank, A., et al. 2019, ApJ, 880, 43

Simpson, J., Smail, I., Dudzviciute, U., et al. 2020, MNRAS, 495, 3409

Simpson, J. M., Smail, I., Wang, W.-H., et al. 2017, ApJL, 844, L10

Smail, I., Ivison, R., & Blain, A. 1997, ApJL, 490, L5

Smolčić, V., Aravena, M., Navarrete, F., et al. 2012, A&A, 548, A4

Smolčić, V., Novak, M., Bondi, M., et al. 2017, A&A, 602, A1

Speagle, J., Steinhardt, C., Capak, P., & Silverman, J. 2014, ApJS, 214, 15

Spilker, J., Marrone, D., Aravena, M., et al. 2016, ApJ, 826, 112

Spitler, L., Labbé, I., Glazebrook, K., et al. 2012, ApJL, 748, L21

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

Stevans, M. L., Finkelstein, S. L., Kawinwanichakij, L., et al. 2021, ApJ, 921, 58

Straatman, C. M., Labbé, I., Spitler, L. R., et al. 2014, ApJL, 783, L14

Strandet, M., Weiss, A., Vieira, J., et al. 2016, ApJ, 822, 80

Strandet, M., Weiss, A., De Breuck, C., et al. 2017, ApJL, 842, L15

Swinbank, A., Smail, I., Chapman, S., et al. 2004, ApJ, 617, 64

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

Toft, S., Smolčić, V., Magnelli, B., et al. 2014, ApJ, 782, 68
Trayford, J. W., Lagos, C. d. P., Robotham, A. S., & Obreschkow, D. 2020, MNRAS, 491, 3937

Umehata, H., Tamura, Y., Kohno, K., et al. 2015, ApJL, 815, L8

Valentino, F., Tanaka, M., Davidzon, I., et al. 2020, ApJ, 889, 93

van der Vlugt, D., Algera, H. S. B., Hodge, J. A., et al. 2021, ApJ, 907, 5

Vieira, J., Marrone, D., Chapman, S., et al. 2013, Natur, 495, 344

Walter, F., Decarli, R., Aravena, M., et al. 2016, ApJ, 833, 67

Walter, F., Carilli, C., Neeleman, M., et al. 2020, ApJ, 902, 111

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

Weaver, J., Kauffmann, O., Ilbert, O., et al. 2021, arXiv:2110.13923

Whitaker, K., Pope, A., Cybulski, R., et al. 2017, ApJ, 850, 208

Whitaker, K. E., Labbé, I., van Dokkum, P. G., et al. 2011 ApJ, 735, 86

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

Yan, H., Ma, Z., Huang, J.-S., & Fan, L. 2020, ApJS, 249, 1

Younger, J., Fazio, G., Huang, J.-S., et al. 2007, ApJ, 671, 1531

Yun, M., Reddy, N., & Condon, J. 2001, ApJ, 554, 803

Zavala, J., Casey, C., & da Cunha, E. 2018a, ApJ, 869, 71

Zavala, J., Montaña, A., Hughes, D., et al. 2018b, NatAs, 2, 56

Zavala, J., Casey, C., Manning, S., et al. 2021, ApJ, 909, 165