# Characterization of Two 2 mm detected Optically Obscured Dusty Star-forming Galaxies

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

► To cite this version:

Sinclaire M. Manning, Caitlin M. Casey, Jorge A. Zavala, Georgios E. Magdis, Patrick M. Drew, et al.. Characterization of Two 2 mm detected Optically Obscured Dusty Star-forming Galaxies. The Astrophysical Journal, 2022, 925, 10.3847/1538-4357/ac366a. insu-03667484

HAL Id: insu-03667484

https://insu.hal.science/insu-03667484

Submitted on 13 May 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_1.png" /></figure>

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

# Characterization of Two 2 mm detected Optically Obscured Dusty Star-forming Galaxies

Sinclaire M. Manning¹,²,<fnref n="24" />, Caitlin M. Casey¹, Jorge A. Zavala¹,³, Georgios E. Magdis⁴,⁵,⁶, Patrick M. Drew¹, Jaclyn B. Champagne¹, Manuel Aravena⁷, Matthieu Béthermin⁸, David L. Clements⁹, Steven L. Finkelstein¹, Seiji Fujimoto⁴,⁵, Christopher C. Hayward¹⁰, Jacqueline A. Hodge¹¹, Olivier Ilbert⁸, Jeyhan S. Kartaltepe¹², Kirsten K. Knudsen¹³, Anton M. Koekemoer¹⁴, Allison W. S. Man¹⁵,¹⁶, David B. Sanders¹⁷, Kartik Sheth¹⁸, Justin S. Spilker¹, Johannes Staguhn¹⁹,²⁰, Margherita Talia²¹,²², Ezequiel Treister²³, and Min S. Yun²

¹Department of Astronomy, The University of Texas at Austin, 2515 Speedway Boulevard Stop C1400, Austin, TX 78712, USA; smanning.astro@gmail.com

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

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

⁴Cosmic Dawn Center (DAWN), Jagvej 128, DK-2200 Copenhagen N, Denmark

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

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

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

⁸Aix Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, F-13388, Marseille, France

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

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

¹¹Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands

¹²School of Physics and Astronomy, Rochester Institute of Technology, Rochester, NY 14623, USA

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

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

¹⁵Dunlap Institute for Astronomy & Astrophysics, 50 St. George Street, Toronto, ON M5S 3H4, Canada

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

¹⁷Institute for Astronomy, 2680 Woodlawn Drive, University of Hawai'i, Honolulu, HI 96822, USA

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

¹⁹Department of Physics and Astronomy, Johns Hopkins University, 3400 North Charles Street, Baltimore, MD, 21218, USA

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

²¹Department of Physics and Astronomy “Augusto Righi” (DIFA), University of Bologna, Via Gobetti 93/2, I-40129, Bologna, Italy

²²INAF—Osservatorio di Astrofisica e Scienza dello Spazio, Via Gobetti 93/3, I-40129, Bologna, Italy

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

Received 2021 June 30; revised 2021 October 22; accepted 2021 November 2; published 2022 January 21

## Abstract

The 2 mm Mapping Obscuration to Reionization with ALMA (MORA) Survey was designed to detect high-redshift ($z \gtrsim 4$), massive, dusty star-forming galaxies (DSFGs). Here we present two likely high-redshift sources, identified in the survey, whose physical characteristics are consistent with a class of optical/near-infrared (OIR)-invisible DSFGs found elsewhere in the literature. We first perform a rigorous analysis of all available photometric data to fit spectral energy distributions and estimate redshifts before deriving physical properties based on our findings. Our results suggest the two galaxies, called MORA-5 and MORA-9, represent two extremes of the “OIR-dark” class of DSFGs. MORA-5 ($z_{phot} = 4.3^{+1.5}_{-1.3}$) is a significantly more active starburst with a star formation rate (SFR) of $830^{+340}_{-190} M_{\odot} \text{ yr}^{-1}$ compared to MORA-9 ($z_{phot} = 4.3^{+1.3}_{-1.0}$), whose SFR is a modest $200^{+250}_{-60} M_{\odot} \text{ yr}^{-1}$. Based on the stellar masses ($M_{\star} \approx 10^{10-11} M_{\odot}$), space density ($n \sim (5 \pm 2) \times 10^{-6} \text{ Mpc}^{-3}$, which incorporates two other spectroscopically confirmed OIR-dark DSFGs in the MORA sample at $z=4.6$ and $z=5.9$), and gas depletion timescales (<1 Gyr) of these sources, we find evidence supporting the theory that OIR-dark DSFGs are the progenitors of recently discovered $3 < z < 4$ massive quiescent galaxies.

Unified Astronomy Thesaurus concepts: Starburst galaxies (1570); High-redshift galaxies (734); Dust continuum emission (412)

## 1. Introduction

Our current understanding of star formation and galaxy evolution within the first two billion years after the Big Bang is severely limited by a lack of infrared (IR) constraints and associated sample incompleteness at $z > 4$ (Gruppioni et al. 2013; Casey et al. 2014; Madau & Dickinson 2014;

Casey et al. 2018b). Very recent studies (Zavala et al. 2021; Casey et al. 2021) have worked to extend IR measurements up to $z \sim 7$, but uncertainties are still dominant. The census of cosmic star formation out to the highest redshifts is also biased toward unobscured star formation tracers, relying on rest-frame ultraviolet (UV) continuum measurements to seek out Lyman-break galaxies (LBGs). While LBG-based studies have yielded the most robust results out to $z \sim 6$ and beyond and have contributed heavily to analyses of the UV luminosity function (Reddy & Steidel 2009; Finkelstein et al. 2015; Ono et al. 2018), it is well understood that they miss both heavily and moderately dust-obscured star formation by design (Magnelli et al. 2012; Gruppioni et al. 2013). Therefore, the

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

<footnote n="24">²⁴NASA Hubble Fellow.</footnote>
true abundance of dusty star-forming galaxies (DSFGs) in the
early universe remains unclear. Knowing the prevalence of
dust-obscured star formation is particularly important at $z > 4$
(within the first 1.5 billion years of the universe), when
cosmic time becomes a significant constraint on the physical
processes which produce the dust, metals, and stars observed
in galaxies.

Much debate exists over the ubiquity of DSFGs at high
redshifts and their contribution to the total cosmic star
formation within the first billion years after the Big Bang.
The detection of these dust-obscured systems require observa-
tions at far-IR and (sub)millimeter wavelengths which trace the
dust reprocessed UV/optical emission from young stars: an
historically tricky task given the small area of deep pencil-
beam observations and conversely the shallow depths reached
by wider surveys. While some works propose DSFGs as the
dominant source of star formation at $z > 4$ (Rowan-Robinson
et al. 2016; Gruppioni et al. 2020), others find the contribution
from DSFGs to be insignificant (Koprowski et al. 2017; Smith
et al. 2017). In an attempt to reconcile these two disparate
theories, we look to identify individual DSFGs and constrain
their source density at $z > 4$.

Models (Casey et al. 2018a, 2018b) and observations
(Staguhn et al. 2014) show that 2 mm observations offer an
effective strategy to detect high-redshift dusty galaxies at $z \gtrsim 4$
while simultaneously filtering out lower-redshift DSFGs at
$z \lesssim 2.5$. The Mapping Obscuration to Reionization ALMA
(MORA) Survey is the largest (184 arcmin²) 2 mm mosaic
mapped with the Atacama Large Millimeter/submillimeter
Array (ALMA), whose aim is to efficiently select $z \gtrsim 3-4$
DSFGs. For a full description of the MORA Survey, see Zavala
et al. (2021) and Casey et al. (2021). Twelve total sources were
robustly detected ($>5\sigma$) in the MORA survey. Two of these
sources were then identified as being of particular interest given
their robust mm detections, but lack of optical/near-IR (OIR)
counterparts. Systems with these characteristics have the
potential to be the highest-redshift galaxies in the sample or
exceedingly obscured (or both). Additionally, recent works
have unearthed a population of massive galaxies undetected in
deep (5$\sigma$ point-source limit H-band $>27$ mag) *Hubble Space
Telescope (HST)* observations (dubbed *HST*-dark or rest-frame
H/K-band dropouts). As such, we focus on these two sources,
MORA-5 and MORA-9, as they appear to be representative of
this unique class of heavily obscured galaxies. “OIR-dark” is
used as an identifier throughout this paper in reference to
DSFGs with nondetections at UV, optical, and near-IR
wavelengths shorter than 2.2 $\mu$m.

In this paper we present the physical characteristics of two
OIR-dark sources found in the 2 mm MORA survey and
discuss their potential to belong to the parent population of
high-$z$ massive quiescent galaxies. Section 2 presents the
ALMA observations from MORA and archival searches as well
as ancillary data from the Cosmic Evolution Survey Field
(COSMOS) survey (Scoville et al. 2007). Section 3 describes
the various spectral energy distribution (SED) fitting techni-
ques used to fit the photometric data. Section 4 lists the
physical properties derived from SED fitting. Section 5
discusses how these galaxies fit into the broader DSFG
population and compare to known OIR-dark sources in the
literature. Finally, we assess the potential evolutionary pathway
of OIR-dark DSFGs as the progenitors of massive quiescent
galaxies at early times. Throughout this paper we assume a

Planck cosmology with $H_0 = 67.8$ km s⁻¹ Mpc⁻¹ and
$\Omega_m = 0.308$ (Planck Collaboration et al. 2016) and a Kroupa
initial mass function (Kroupa & Weidner 2003) when referring
to SFRs.

## 2. Data & Observations

In total, 13 robust sources ($>5\sigma$) were detected in the 2 mm
mosaic. We direct the reader to the accompanying MORA
papers for a more thorough overview of the survey and
complete discussions on 2 mm number counts (Zavala et al.
2021) and source characterization (Casey et al. 2021). What
follows is an explanation of the ALMA data reduction and
imaging, including available archival data, as well as the
existing ancillary data for MORA-5 and MORA-9, which are
the focus of this paper. All available photometry is compiled in
Table 1.

### 2.1. ALMA Data

The 2 mm data were obtained during Cycle 6 under program
2018.1.00231.S (PI: Casey), MORA-5 and MORA-9 both sit in
the contiguous 156 arcmin² mosaic, where the rms reached
89 µJy and 74 µJy at their respective source positions. We select
all pointings contributing more than 5% to the total sensitivity at
the given positions (21 pointings total for both sources) before
reducing and imaging the data using the Common Astronomy
Software Application<fnref n="25" /> (CASA), version 5.6.1, following the
standard reduction pipeline scripts. We experimented with
different manually defined clean boxes during the cleaning
process, of order the size of the emission. We adopted a 1" aperture centered on the source position as our clean box,
encompassing the peak of emission without including any
noise. This is a slightly different reduction independent from
the large mosaic presented in Zavala et al. (2021), customized
for the individual sources given our ability to test multiple
setups, but the results are consistent with the full mosaic.

Band 4 (2 mm) observations covered frequencies
139.4–143.2 GHz and 151.5–155.2 GHz with a restoring beam
of 1.167 ′′ × 1.141. The continuum rms over the 7.5 GHz bandwidth
reached 85.1 µJy/beam. Various visibility weights were explored
for imaging as we experimented with both natural and
Briggs weighting (robust = 0.0, 0.5, 2.0), the latter in an attempt
to detect any spatially resolved components, though in the end we
determined both sources to be unresolved. No spectral lines
were detected in the covered bandwidth after close inspection.
Based on the continuum imaging, the aggregate band 4 (central
frequency of 147.3 GHz) flux densities of MORA-5 and MORA-
9 are 610 ± 90 µJy and 380 ± 70 µJy, respectively.

Band 6 (1.3 mm) and band 7 (870 $\mu$m) continuum data
exists from other ALMA programs 2016.1.00279.S (PI: Oteo)
and 2016.1.00463.S (PI: Matsuda) for MORA-5. These data were
downloaded from the archive and reimaged. Band 7 observations
covered frequencies 335.6–339.6 GHz and 347.6–351.4 GHz
with a restoring beam of 0.94 ′′ × 0.84 ′′. Band 6 observations
covered frequencies 223–227 GHz and 239–243 GHz with a
restoring beam of 0.90 ′′ × 0.76 ′′. The emission in both bands is
unresolved and the reported band 6 and band 7 flux densities for
MORA-5 are 2.26 ± 0.17 mJy and 6.49 ± 0.22 mJy, respectively.

Comparatively, MORA-9 lacks coverage at other frequencies,
so we secured follow-up Atacama Compact Array (ACA)

<footnote n="25">²⁵ http://casa.nrao.edu</footnote>
<table><caption>**Table 1**
Photometry of MORA-5 and MORA-9</caption><thead><tr><th>Band</th><th>Wavelength</th><th>Units</th><th>MORA-5</th><th>MORA-9</th><th>Reference</th></tr></thead><tbody><tr><td>HST-F606W</td><td>606 nm</td><td>nJy</td><td>(-11 ± 24)</td><td>(50 ± 24)</td><td>Koekemoer et al. (2011)</td></tr><tr><td>HST-F814W</td><td>814 nm</td><td>nJy</td><td>(8 ± 32)</td><td>(-21 ± 32)</td><td>Koekemoer et al. (2011)</td></tr><tr><td>HST-F125W</td><td>1.25 µm</td><td>nJy</td><td>(-12 ± 44)</td><td>(-28 ± 43)</td><td>Koekemoer et al. (2011)</td></tr><tr><td>HST-F140W</td><td>1.40 µm</td><td>nJy</td><td>(145 ± 83)</td><td>—</td><td>Brammer et al. (2012)</td></tr><tr><td>HST-F160W</td><td>1.60 µm</td><td>nJy</td><td>(40 ± 42)</td><td>(40 ± 43)</td><td>Koekemoer et al. (2011)</td></tr><tr><td>UVISTA-Ks</td><td>2.2 µm</td><td>nJy</td><td>(59 ± 40)</td><td>269 ± 40</td><td>McCracken et al. (2012)</td></tr><tr><td>IRAC-CH1</td><td>3.6 µm</td><td>nJy</td><td>680 ± 140</td><td>760 ± 150</td><td>Ashby et al. (2015)</td></tr><tr><td>IRAC-CH2</td><td>4.5 µm</td><td>nJy</td><td>790 ± 160</td><td>880 ± 170</td><td>Ashby et al. (2015)</td></tr><tr><td>MIPS24</td><td>24 µm</td><td>µJy</td><td>(1 ± 20)</td><td>(-2 ± 10)</td><td>Le Floc'h et al. (2009)</td></tr><tr><td>PACS</td><td>100 µm</td><td>mJy</td><td>(-0.2 ± 1.7)</td><td>(-0.5 ± 1.7)</td><td>Lutz et al. (2011)</td></tr><tr><td>PACS</td><td>160 µm</td><td>mJy</td><td>(-0.1 ± 3.7)</td><td>(-0.3 ± 3.7)</td><td>Lutz et al. (2011)</td></tr><tr><td>SPIRE</td><td>250 µm</td><td>mJy</td><td>(3.8 ± 5.8)</td><td>(0 ± 5.8)</td><td>Oliver et al. (2012)</td></tr><tr><td>SPIRE</td><td>350 µm</td><td>mJy</td><td>(7.3 ± 6.3)</td><td>(0 ± 6.3)</td><td>Oliver et al. (2012)</td></tr><tr><td>SCUBA-2</td><td>450 µm</td><td>mJy</td><td>(8.53 ± 4.13)</td><td>—</td><td>Casey et al. (2013)</td></tr><tr><td>SPIRE</td><td>500 µm</td><td>mJy</td><td>(7.1 ± 6.8)</td><td>(0 ± 6.8)</td><td>Oliver et al. (2012)</td></tr><tr><td>ALMA-B8</td><td>740 µm</td><td>mJy</td><td>—</td><td>3.62 ± 0.66</td><td>THIS WORK</td></tr><tr><td>SCUBA-2</td><td>850 µm</td><td>mJy</td><td>6.80 ± 0.53</td><td>2.61 ± 0.59</td><td>Simpson et al. (2019)</td></tr><tr><td>ALMA-B7</td><td>868 µm</td><td>mJy</td><td>6.49 ± 0.22</td><td>—</td><td>THIS WORK</td></tr><tr><td>ALMA-B7</td><td>874 µm</td><td>mJy</td><td>—</td><td>2.54 ± 0.37</td><td>THIS WORK</td></tr><tr><td>AzTEC</td><td>1100 µm</td><td>mJy</td><td>3.7 ± 0.9</td><td>—</td><td>Aretxaga et al. (2011)</td></tr><tr><td>ALMA-B6</td><td>1287 µm</td><td>mJy</td><td>2.26 ± 0.17</td><td>1.24 ± 0.25</td><td>THIS WORK</td></tr><tr><td>ALMA-B4</td><td>2036 µm</td><td>µJy</td><td>(610 ± 90)</td><td>380 ± 70</td><td>THIS WORK</td></tr><tr><td>VLA-3 GHz</td><td>10 cm</td><td>µJy</td><td>(10.1 ± 3.4)</td><td>(4.3 ± 2.4)</td><td>Smolčić et al. (2017a)</td></tr><tr><td>VLA-1.4 GHz</td><td>21.4 cm</td><td>µJy</td><td>(20.0 ± 12.5)</td><td>(-4.2 ± 10.8)</td><td>Schinnerer et al. (2007)</td></tr></tbody></table>

Note. Available photometric measurements for MORA-5 and MORA-9. Nondetections and measurements with <3σ significance are denoted by parentheses. MORA-9 sits just outside the coverage area of the HST-F140W observations. Neither source is detected in the intermediate band filters available in COSMOS. ALMA band 7 and band 6 were downloaded from the archive and reimaged for MORA-5. Band 6, band 7, and band 8 data for MORA-9 were obtained with follow-up ACA observations.

observations under program 2019.2.00143 (PI: Casey) to obtain band 6, band 7, and band 8 observations for improved dust continuum constraints. The band 6 observations resulted in a detection just under 5σ with flux density of 1.24 ± 0.25 mJy. The aggregate band 7 continuum has a flux density of 2.54 ± 0.37 mJy. This is consistent with SCUBA-2 data (see Section 2.2), but is more precise both astrometrically and in the measurement of its flux density. The band 8 continuum flux density for MORA-9 is 3.62 ± 0.66 mJy. Given the lower resolution of the ACA data (FWHM = 5.5 × 3.3) we adopt a 2" aperture centered on the source position as our clean box. Figures 1 and 2 show the ALMA continuum maps for MORA-5 and MORA-9. Natural weighting (robust = 2) is used for all MORA-5 imaging as well as the 740 µm (because of its higher rms) and 2 mm continuum map of MORA-9. Briggs weighting (robust = 0) is utilized for the 870 µm map of MORA-9 to maximize spatial resolution since the ACA data mitigates any loss in flux expected by not using natural weighting.

## 2.2. Ancillary Archival COSMOS Data sets

Both sources lie in the Cosmic Evolution Survey Field (COSMOS; Capak et al. 2007; Scoville et al. 2007) within the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011; Koekemoer et al. 2011) area. We explore a range of apertures centered on the ALMA 2 mm positions and extract photometry from the available OIR data. Apertures of 1.5" diameter are placed down on the HST-F606W, F814W, F125W, F140W, and F160W images resulting in no detections.

While neither source was initially found to have counterparts in the original Ks-selected catalog (Muzzin et al. 2013) within 1", a Ks-band detection for MORA-9 is now reported in the updated COSMOS2020 catalog (Weaver et al. 2022) after the inclusion of the fourth data release (DR4) of the UltraVISTA survey (McCracken et al. 2012) which folds in deeper (up to 1 mag) Ks-band imaging. MORA-5, however, still lacks a robust detection down to the 25.2 mag (3σ) limit with a signal-to-noise ratio (SNR) ~1.5.

Both sources do have detections in the deep S-CANDELS Spitzer IRAC data (Ashby et al. 2015) at 3.6 µm and 4.5 µm. Neither MORA-5 nor MORA-9 are detected in Spitzer 24 µm imaging (Le Floc'h et al. 2009), Herschel Photodetector Array Camera and Spectrometer (PACS) 100 µm and 160 µm imaging (Lutz et al. 2011), or Spectral and Photometric Imaging Receiver (SPIRE) 250 µm, 350 µm, and 500 µm imaging (Oliver et al. 2012).

The sources were also both detected with SCUBA-2 as reported in Geach et al. (2017) and then again in the deeper, expanded S2COSMOS survey presented in Simpson et al. (2019). The ratio of 850 µm flux density to 2 mm flux density for both MORA-5 and MORA-9 is what initially indicated these sources were likely to be high redshift at z ≥ 3 ($S_{850\mu m}/S_{2\text{mm}}$=1.1 ± 2.5 and 6.9 ± 2.8, respectively). Specifically, the ratio for MORA-9 is significantly lower than the expected value for emission on the Rayleigh-Jeans tail of the cold dust blackbody. This suggests that the SCUBA-2 measurement at 850 µm is near the peak of the dust SED rather than lower down on the Rayleigh-Jeans tail.

Lastly, neither source was found to have a detection significant enough to be included in the 1.4 GHz Karl G. Jansky Very
<figure><img src="image_4.png" /><figcaption>**Figure 1.** 10" × 10" cutouts of 870 µm, 1.3 mm, and 2 mm continuum imaging for MORA-5 using natural weighting (robust = 2). Contours in the 870 µm image show 5, 10, 15, 20, and 25σ, in contrast to the 1.3 mm and 2 mm images which depict 3, 5, 7, 9, and 11σ.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** 20" × 20" and 10" × 10" cutouts of 740 µm (ACA), 870 µm (ACA) and 2 mm continuum imaging for MORA-9. Contours represent 3, 5σ in all three images. The star denotes the 2 mm ALMA position of MORA-9 and the white dashed box depicts the 10" × 10" cutout to emphasize the varying scales between the images. Briggs weighting (robust = 0) is used to create the 870 µm image to maximize the spatial resolution of the ACA data while natural weighting (robust = 2) is used for the 740 µm and 2 mm images with lower signal-to-noise.</figcaption></figure>

Large Array (VLA)-COSMOS catalog ($3\sigma \sim 30$ µJy beam$^{-1}$; Schinnerer et al. 2007) or the 3 GHz VLA-COSMOS catalog ($5\sigma \sim 12$ µJy beam$^{-1}$; Smolčić et al. 2017a). After manually extracting fluxes from the radio imaging at the 2 mm source positions, we find marginal detections of $S_{3GHz} = 10.1 \pm 2.4$ µJy for MORA-5 and $S_{3GHz} = 4.3 \pm 3.4$ µJy for MORA-9.

### 3. SED Fitting Techniques

We explore the use of several SED fitting techniques in relation to the available multiwavelength data on MORA-5 and MORA-9 to derive redshift constraints and physical properties. We elect to input the photometry exactly as it appears in Table 1 into our SED fitting, and therefore our photometric redshift estimates, preserving all formal non-detections rather than substituting them with upper limits. We also note here for clarity that the final redshift estimates adopted for MORA-5 and MORA-9 are a combination of the photometric redshift probability distribution functions (PDFs) of each technique and are referred to as $z_{PDF}$ throughout the paper. Due to the highly obscured nature of these galaxies, it would be possible to have an extremely obscured, optically thick central region. Implementing energy balance techniques in SED fitting works best when data sampling the full SED, both the OIR and long-wavelength regimes, is available,

and the galaxy in question does not have much optically thick dust. For these reasons, we fit the SEDs of the OIR and far-IR (FIR)/millimeter regimes separately, but also fit the full SED for a fair comparison of the differing techniques. We recognize the challenge of fitting SEDs without well-constrained redshifts, and as we are lacking spectroscopic redshifts for MORA-5 and MORA-9, we approach the process iteratively. Here we describe the tools used for photometric redshift and SED fitting; Section 4.1 describes our results with respect to redshift constraints.

#### 3.1. FIR/millimeter SED and MMPZ

We fit the FIR/millimeter-obscured SED with a piecewise function comprised of a single modified blackbody and mid-IR power law using a Metropolis Hastings Markov Chain Monte Carlo (MCMC). For a complete description of the FIR/millimeter SED fitting technique, called MCIRSED, we encourage the reader to see the forthcoming publication by P. Drew et al. (submitted), but briefly explain the model here. At short wavelengths, the FIR/millimeter SED is dominated by a power law with mid-IR slope, $\alpha$, representing warmer dust emission from star-forming regions and/or active galactic nuclei (AGN). At longer wavelengths, a modified Planck function dominates the SED and represents the cold dust. An
<figure><img src="image_6.png" /><figcaption>**Figure 3.** 8" × 8" multiwavelength cutouts of MORA-5 (top) and MORA-9 (bottom) from HST ACS/F606W, HST WFC3/F125W, HST WFC3/F160W, UltraVISTA Ks-band, Spitzer 3.6 µm and 4.5 µm, ALMA 2 mm, and VLA-3 GHz. Orange contours in each frame show the 5σ contours (3 and 5σ for MORA-9) from the 2 mm imaging.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** $L_{IR}-\lambda_{peak}$ tracks traced by SEDs fit to the available photometry with associated uncertainty from $z=2-8$. We show results using $\beta = 1.8$ (teal) and 2.2 (purple) to assess if the value of $\beta$ affects where these tracks fall in the $L_{IR}-\lambda_{peak}$ parameter space. Inset plots show the corresponding photometric redshift distributions from the results of MMPZ. The MMPZ redshift is dictated by the intersection of the SED fit tracks and the intrinsic $L_{IR}-\lambda_{peak}$ relation in galaxies (black line with 1,2σ scatter in gray).</figcaption></figure>

MCMC is then employed to sample the posteriors of the fit parameters and provide confidence intervals. We fix the mid-IR power-law slope, $\alpha = 4$, and emissivity spectral index, $\beta = 1.8$, following Casey et al. (2019). We also fix the wavelength where opacity equals unity, $\lambda_0$, to 200 µm as is commonly adopted in the literature. Finally, we leave the dust temperature ($T_{dust}$), wavelength corresponding to the peak in the FIR/millimeter SED ($\lambda_{peak}$), and IR luminosity ($L_{IR(8-1000\mu m)}$) as free parameters.

Casey et al. (2018b) (see Figure 3 therein) finds an empirical relationship between galaxies’ intrinsic $L_{IR(8-1000\mu m)}$ and their observed $\lambda_{peak}$ given the typical SEDs of cold dust emission in DSFGs. Comparing these two properties provides us with an approximation for the range of redshift solutions which are feasible based on existing observations. This empirical relationship is the basis for the MMPZ fitting technique (Casey 2020) which we use to obtain FIR photometric redshift estimates. The results of this method are illustrated in Figure 4. The $L_{IR}-\lambda_{peak}$ trend is shown by the black line with its associated 1 and 2σ scatter from measured data in gray. Across a broad span of redshifts, the teal and purple tracks trace the range of possible SEDs constrained by the measured photometry given $\beta$ values of 1.8 and 2.2 (fixed only to $\beta = 1.8$ for simplicity in Casey 2020). As we can

see, the choice of $\beta$ has a negligible effect on the redshift solution, so we adopt the $\beta = 1.8$ value for consistency. The associated errors reflect the number of photometric data points available to sample the SEDs, resulting in larger errors for MORA-9.

## 3.2. EAZY OIR SED

We fit the sparse OIR data with the photometric redshift fitting software EAZY (v1.3) (Brammer et al. 2008) using the provided set of templates. The eight EAZY templates result from a linear combination of ~500 Bruzual & Charlot (2003) synthetic galaxy photometry models and includes a dusty starburst model, an important addition for this work. Due to the limited OIR data, particularly for MORA-5, our results show very broad photometric redshift probability distributions, as expected. The EAZY SEDs are shown as blue lines in Figure 5 and their associated PDFs correspond to the blue distributions in Figure 6.

## 3.3. MAGPHYS Energy Balance SED

Lastly, the full SED is fit using MAGPHYS and the updated MAGPHYS+photo-z code (da Cunha et al. 2008, 2015; Battisti et al. 2019). The original MAGPHYS incorporates an
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Composite SED (black) comprised of two components: OIR stellar emission (blue) from EAZY and thermal dust emission from MMPZ (orange). Black arrows depict $3\sigma$ upper limits for illustrative purposes, while the exact photometric values listed in Table 1 are used for SED creation. Light and dark gray lines are MAGPHYS SEDs, having provided the original MAGPHYS with MMPZ $z_{\text{mode}}$ (light gray) and $z_{\text{PDF}}$ (dark gray) as input redshifts. Note that the offset between MAGPHYS SEDs is due to the conversion from flux density to $\nu L_{\nu}$ being dependent on luminosity distance (i.e., redshift). The OIR SED, FIR/millimeter SED, data points, and upper limits have been converted given the luminosity distance at $z_{\text{PDF}}$. Bottom: the purple line shows the MAGPHYS SED of MORA-9 based on the redshift reported in the COSMOS2020 catalog.</figcaption></figure>

energy balance technique taking into account emission from stellar populations as well as absorption and emission by dust in galaxies. We opt to use MAGPHYS+photo-z results rather than the classic MAGPHYS given its inclusion of IR, (sub) millimeter, and radio data, plus new star formation history and dust temperature priors which are more appropriate for high-redshift DSFGs. Specifically, we adopt stellar masses and rest-frame V-band dust attenuation ($A_{\nu}$) values from this updated version. However, we obtain the SEDs shown in Figure 5 by running the original MAGPHYS, setting the redshift to that determined from MMPZ as well as the final $z_{\text{PDF}}$. The light gray lines are the MAGPHYS SEDs given $z_{\text{mode}}$ from MMPZ while the dark gray shows the MAGPHYS SED assuming the $z_{\text{PDF}}$. For MORA-9, we also run MAGPHYS given the COSMOS2020 redshift, and this is shown as the purple line in the bottom panel. We note that the SEDs produced by MAGPHYS+photo-z are also consistent with the original MAGPHYS SEDs when $z_{\text{PDF}}$ is provided as the redshift. Thus, we do not show the MAGPHYS+photo-z SEDs in Figure 5.

# 4. Results

We rely on photometry to estimate redshifts and several physical properties for MORA-5 and MORA-9. Here we discuss these physical properties, based on the $z_{\text{PDF}}$ derived from the various SED fitting techniques and computed in Section 4.1, and list them in Table 2.

## 4.1. Photometric Redshifts

Photometric redshift PDFs from all SED fitting techniques utilized are displayed in Figure 6. The redshift determined from each one and associated $1\sigma$ errors are depicted by the triangles at the top of the plot. We sum all available PDFs (three for MORA-5 and four for MORA-9 after the inclusion of the COSMOS2020 distribution) and adopt the value from marginalizing over the PDF as our final photometric redshift, which we refer to as $z_{\text{PDF}}$ (black distribution). Combining photometric redshift PDFs from different codes in this way has proven to be advantageous (Dahlen et al. 2013). For MORA-5, $z_{\text{PDF}} = 4.3_{-1.3}^{+1.5}$ and for MORA-9, $z_{\text{PDF}} = 4.3_{-1.0}^{+1.3}$. The available
<figure><img src="image_9.png" /><figcaption>**Figure 6.** PDFs of the three methods used to determine photometric redshifts for MORA-5 (left) and MORA-9 (right): OIR photometry with EAZY (blue), energy balance technique with MAGPHYS+photo-z (gray), and FIR/millimeter photometry with MMPZ (orange). For MORA-9 (right), we also show the PDF result from the LePhare template fitting code (Ilbert et al. 2006), following the same method as Laigle et al. (2016) as reported in the COSMOS2020 catalog (magenta). $1\sigma$ errors are shown at the top with triangles denoting the photometric redshifts. For our EAZY and MAGPHYS+photo-z (shortened in the legend to MAGPHYS for space) results, $z_{50\text{th}}$ signifies the redshifts determined from marginalizing over their respective distributions. We report the final redshift estimates ($z_{\text{PDF}}$) for these sources as the 50th percentile value of the combined distribution of all three (four for MORA-9) PDFs, illustrated by the black lines.</figcaption></figure>

<table><caption>**Table 2**
Derived Properties of MORA-5 and MORA-9</caption><thead><tr><th>Property</th><th>Units</th><th>MORA-5</th><th>MORA-9</th></tr></thead><tbody><tr><td>RA</td><td>—</td><td>10:00:24.157</td><td>10:00:17.298</td></tr><tr><td>DEC</td><td>—</td><td>+02:20:05.39</td><td>+02:27:15.80</td></tr><tr><td>z<sub>PDF</sub></td><td>—</td><td>4.3<sup>+1.5</sup><sub>-1.3</sub></td><td>4.3<sup>+1.3</sup><sub>-1.0</sub></td></tr><tr><td>L<sub>IR(8–1000µm)</sub></td><td>L<sub>☉</sub></td><td>(5.6<sup>+2.3</sup><sub>-1.3</sub>) × 10<sup>12</sup></td><td>(1.4<sup>+1.7</sup><sub>-0.4</sub>) × 10<sup>12</sup></td></tr><tr><td>SFR</td><td>M<sub>☉</sub> yr<sup>-1</sup></td><td>830<sup>+340</sup><sub>-190</sub></td><td>200<sup>+250</sup><sub>-60</sub></td></tr><tr><td>λ<sub>peak (rest)</sub></td><td>µm</td><td>83<sup>+9</sup><sub>-8</sub></td><td>106<sup>+19</sup><sub>-20</sub></td></tr><tr><td>T<sub>dust</sub></td><td>K</td><td>60<sup>+8</sup><sub>-7</sub></td><td>43<sup>+14</sup><sub>-9</sub></td></tr><tr><td>M<sub>☉</sub></td><td>M<sub>☉</sub></td><td>(1.5<sup>+1.0</sup><sub>-0.7</sub>) × 10<sup>11</sup></td><td>(4.1<sup>+1.8</sup><sub>-1.4</sub>) × 10<sup>10</sup></td></tr><tr><td>M<sub>dust</sub></td><td>M<sub>☉</sub></td><td>(3.6 ± 0.5) × 10<sup>9</sup></td><td>(2.2 ± 0.4) × 10<sup>9</sup></td></tr><tr><td>M<sub>gas(2 mm)</sub></td><td>M<sub>☉</sub></td><td>(2.6 ± 0.4) × 10<sup>11</sup></td><td>(1.6 ± 0.3) × 10<sup>11</sup></td></tr><tr><td>A<sub>V</sub></td><td>—</td><td>4.3<sup>+0.9</sup><sub>-0.7</sub></td><td>2.9<sup>+0.4</sup><sub>-0.5</sub></td></tr><tr><td>q<sub>IR</sub></td><td>—</td><td>> 1.9</td><td>> 1.6</td></tr></tbody></table>

**Note.** Derived properties based on the separately fit OIR and FIR/millimeter SEDs. The $z_{\text{PDF}}$ reported comes from combining the photometric redshift PDFs of all SED fitting techniques discussed in Section 3. $L_{\text{IR(8–1000µm)}}$ is the derived IR luminosity integrated from 8–1000 µm. SFRs are determined directly from our $L_{\text{IR}}$ results using the associated conversion factors reported in Kennicutt & Evans (2012). $\lambda_{\text{peak}}$ is the rest-frame wavelength where the FIR/millimeter dust SED peaks. Stellar masses and $A_V$ are taken from the MAGPHYS+photo-z energy balance SED results. Dust and gas masses are derived directly from the 2 mm dust continuum measurements. $q_{\text{IR}}$ is computed directly from the FRC and is reported as a lower limit for MORA-5 given its marginal detection at 3 GHz.

ALMA observations only offer a few, narrow frequency ranges for molecular lines to be detected, and a detection at these frequencies would be indicative of a $z < 4$ redshift solution, providing some evidence for our adopted estimates. Next, we describe the PDFs and associated photometric redshifts which contribute to the $z_{\text{PDF}}$ values.

From MMPZ (orange distribution), $z_{\text{mode}}$ is the mode of the PDF derived from the intersection of the galaxy's redshift track with the $L_{\text{IR}}-\lambda_{\text{peak}}$ trend from the literature as shown in Figure 4. For MORA-5, $z_{\text{mode}} = 3.4_{-0.9}^{+1.1}$ and for MORA-9, $z_{\text{mode}} = 5.2_{-1.0}^{+0.9}$. This is the highest redshift estimate for MORA-9. Low-redshift

solutions for this source would imply implausibly cold dust for $z \lesssim 3$ ($\lambda_{\text{peak}} > 163$ µm or $T_{\text{dust}} < 18$ K). This would be an extremely cold system considering the luminosity of MORA-9 ($L_{\text{IR}} \sim 10^{12} L_{\odot}$, which is largely insensitive to redshift) and as such we find a higher-redshift solution to be more likely given its FIR-mm colors.

Results from EAZY are illustrated in blue in both panels of Figure 6. We show the redshifts determined from marginalizing over the full EAZY PDF ($z_{50\text{th}} = 4.6 \pm 1.9$ for MORA-5 and $z_{50\text{th}} = 3.4_{-0.7}^{+0.8}$ for MORA-9). The broad distribution of the EAZY PDF is expected given the limited OIR photometric data, particularly for MORA-5, which lacks a Ks-band detection and is only constrained by two IRAC detections.

The photometric redshift distributions from MAGPHYS+photo-z are displayed in gray. Similar to our EAZY results, we show the 50th percentile redshift values for each source. For MORA-5 the redshift estimate is $z_{50\text{th}} = 4.8_{-1.1}^{+1.0}$, while for MORA-9 it is $z_{50\text{th}} = 4.2 \pm 0.6$.

As noted in Section 2, MORA-9 was found to have a Ks-band counterpart in the most recent UltraVISTA data and consequently is included in the COSMOS2020 catalog with a corresponding redshift. COSMOS2020 includes photometric redshift estimates using both EAZY and the same method as cited in Laigle et al. (2016), which utilizes the template-fitting code LePhare (Arnouts et al. 2002; Ilbert et al. 2006). We refer the reader to Weaver et al. (2022) for a full explanation, but include these results for comparison and fold in the COSMOS2020 PDF into the final $z_{\text{PDF}}$ estimate for MORA-9. The redshift reported for MORA-9 in the Classic LePhare COSMOS2020 catalog is $z_{50\text{th}} = 4.57_{-0.89}^{+0.87}$.

## 4.2. IR Luminosity and Star Formation Rate

Total IR luminosities ($L_{\text{IR}}$) are determined by integrating the MMPZ FIR SEDs from 8–1000 µm. In Table 2 we list the $L_{\text{IR}}$ assuming the combined $z_{\text{PDF}}$ values. With these values we then calculate SFRs utilizing calibrators from the literature (Kennicutt & Evans 2012 with references to Murphy et al. 2011; Hao et al. 2011 therein) and obtain SFR estimates of $830_{-190}^{+34} M_{\odot} \text{ yr}^{-1}$ for MORA-5 and $200_{-60}^{+250} M_{\odot} \text{ yr}^{-1}$ for
MORA-9. If we were to assume the MMPZ-determined redshifts and corresponding $L_{IR}$ instead, the SFRs shift to $680_{-120}^{+140} M_{\odot}$ yr$^{-1}$ and $340_{-100}^{+350} M_{\odot}$ yr$^{-1}$ respectively, remaining consistent to the previous SFRs given the large uncertainties. Since the IR luminosities are redshift dependent, the SFRs increase if we assume higher-redshift solutions. For MORA-5, MMPZ produces a lower redshift than the result of the combined PDFs, whereas MORA-9 is the inverse case and MMPZ produces a higher redshift. Had we used the least-constrained photometric redshifts produced by EAZY, the opposite would occur: the higher-redshift estimate of MORA-5 would push it to an extreme IR luminosity ($\sim 10^{13} L_{\odot}$) and subsequent SFR, while the lower-redshift estimate of MORA-9 would suggest the system more closely resembles a normal, luminous infrared galaxy (LIRG; Sanders & Mirabel 1996) found in the local universe with $L_{IR} \sim 5 \times 10^{11} L_{\odot}$.

### 4.3. Stellar Mass and Attenuation

We report the stellar masses ($M_*$) and absolute magnitudes of extinction ($A_V$) values inferred from MAGPHYS+photo-z. As shown in Battisti et al. (2019), $M_*$ can be underestimated by a factor of two for DSFGs, and most divergent at higher stellar masses ($\gtrsim 10^{10} M_{\odot}$) when using codes which rely only on rest-frame UV-NIR data such as EAZY and the original MAGPHYS. Given this, we report the stellar masses for MORA-5 and MORA-9 to be $(1.5_{-0.7}^{+1.0}) \times 10^{11} M_{\odot}$ and $(4.1_{-1.4}^{+1.8}) \times 10^{10} M_{\odot}$ while $A_V$ is $4.3_{-0.7}^{+0.9}$ and $2.9_{-0.5}^{+0.4}$, respectively. We note that both estimates are only marginally dependent on redshift within each galaxy's PDF. Any further improvements in our stellar mass estimates will require better constraints on both redshifts and the near/mid-IR SED from, for example, future JWST observations.

### 4.4. Rest-frame Peak Wavelength and Dust Temperature

As described in Section 3.1, $\lambda_{peak}$ is the wavelength at which the best-fit model peaks after fitting the FIR/millimeter data with a mid-IR power law plus a modified blackbody. The relationship between the measurable quantity, $\lambda_{peak}$, and dust temperature, $T_{dust}$, is dependent on the dust opacity model assumed (see Figure 20 of Casey et al. 2014), and here we assume $\tau = 1$ at $\lambda_{rest} = 200 \mu m$. The dust temperature is set as a free parameter in the model, and we extract our adopted values from the best-fit FIR/millimeter SED. Assuming $z_{PDF}$, $\lambda_{peak}$ and $T_{dust}$ are $83_{-8}^{+9} \mu m$ and $60_{-7}^{+8} K$ for MORA-5 and $106_{-20}^{+19} \mu m$ and $43_{-9}^{+14} K$ for MORA-9.

### 4.5. Dust Mass

We directly infer dust masses from the associated 2 mm dust continuum detections on the Rayleigh-Jeans tail of blackbody emission (at $\lambda_{rest} \gtrsim 300 \mu m$ where dust emission is likely to be optically thin) following the framework of Scoville et al. (2016). The dust mass is proportional to the mass-weighted dust temperature (set to 25 K, representative of the bulk of dust dispersed throughout the entire galaxy) and observed flux density at $\nu_{obs}$. Given the potentially high-redshift solutions, we assume cosmic microwave background (CMB) heating of the dust in these sources is nonnegligible, so we employ Equation (1) from Casey et al. (2019) for our dust mass calculation as it incorporates a correction factor for suppressed flux density against the CMB background. We estimate dust masses of $(3.6 \pm 0.5) \times 10^9 M_{\odot}$ for MORA-5 and $(2.2 \pm 0.4) \times 10^9 M_{\odot}$

for MORA-9 with the associated errors propagated from the measurement error of the 2 mm fluxes. The dust-to-stellar mass ratios ($M_{dust}/M_*$) for our sources are $0.024 \pm 0.018$ and $0.054 \pm 0.031$, respectively. When compared to a recent compilation of $M_{dust}/M_*$ for DSFGs (Donevski et al. 2020), we find that both MORA-5 and MORA-9 lie above the determined scaling relation at their respective redshifts, which is identified to rise up to $z \sim 2$ before a mild decline/flattening at $z \gtrsim 2$ is observed. Contrary to low-redshift galaxies, we find our higher $M_{dust}/M_*$ values to be unsurprising considering the high dust content is assumed to be primarily responsible for their defining features as OIR-dark and IR-luminous sources.

### 4.6. Gas Mass

Functionally, our gas mass estimates are determined in the same manner as dust mass, as they are both inferred from the 2 mm data with gas mass being scaled via a dust-to-gas ratio. We adopt a CO-to-H$_2$ conversion factor of $\alpha_{CO} = 6.5 M_{\odot}$ (K km s$^{-1}$ pc$^{-2}$)$^{-1}$, accounting for the mass of both H$_2$ and He gas. This follows the methodology of Scoville et al. (2016) and Equation (3) of Casey et al. (2019) and accounts for CMB heating. This method makes use of the empirically calibrated conversion factor from 850 $\mu m$ luminosity to interstellar medium (ISM) mass, $\alpha_{850} = (6.7 \pm 1.7) \times 10^{19} \text{ erg s}^{-1} \text{ Hz}^{-1} M_{\odot}^{-1}$, as determined in Scoville et al. (2016), and intrinsic to this calculation is the assumed CO-to-H$_2$ conversion factor. We adopt the same mass-weighted dust temperature of 25 K as we did in Section 4.5. Ultimately, we find $M_{gas} = (2.6 \pm 0.4) \times 10^{11} M_{\odot}$ for MORA-5 and $M_{gas} = (1.6 \pm 0.3) \times 10^{11} M_{\odot}$ for MORA-9.

## 5. Discussion

### 5.1. OIR-dark Galaxies: A Subsample of the SMG/DSFG Population

We suspect OIR-dark DSFGs make up part of the high-z tail of the population of dust-obscured submillimeter galaxies (SMGs; Hughes et al. 1998; Blain et al. 2002; Casey et al. 2014). Canonically selected at $850 \mu m$ with $S_{850} \gtrsim 2$ mJy, SMGs are the most well-studied DSFGs, with over two decades of work dedicated to understanding their redshift distributions (Chapman et al. 2005; Zavala et al. 2014; Brisbin et al. 2017), number counts (Scott et al. 2002), spatial density (Aravena et al. 2010; Smolčić et al. 2017b), and intrinsic physical properties (Miettinen et al. 2017; Dudzevičiūtė et al. 2020). SMGs on average have reported stellar and gas masses $\sim 10^{11} M_{\odot}$, SFRs in the range of 100–1000 $M_{\odot}$ yr$^{-1}$, and have often been cited as the progenitors of massive quiescent galaxies, e.g., ellipticals in the centers of local galaxy groups and clusters (Hainline et al. 2011). Similarly, $z \gtrsim 3$ SMGs have been suggested to be the parent source for intermediate-redshift ($z \sim 2$), compact quiescent galaxies (Toft et al. 2014). Moving forward, we group SMGs under the more generalized term of DSFGs for simplicity, which encompasses all galaxies at high-z originally selected at FIR or submillimeter wavelengths. The OIR-dark galaxies discussed in this work are part of this population, having also been first detected at (sub)mm wavelengths and exhibiting similarly broad ranges of stellar masses and star formation activity.

The majority of known DSFGs at $z > 4$ are almost all extreme star-forming systems (forming upwards of 1000 $M_{\odot}$ yr$^{-1}$ and with stellar masses exceeding $10^{11} M_{\odot}$) as they
are predominantly detected by Herschel. As with early DSFG observations, initial OIR-dark DSFG discoveries were of the most extreme, rare starbursts (Walter et al. 2012), and not necessarily termed OIR-dark. Gas-/dust-rich, optically undetected DSFGs have also been discovered as companion systems to (sub)millimeter bright, high-redshift quasars (Fogasy et al. 2017, 2020). Now, with the advent of longer wavelength interferometric surveys covering wider areas, independent of single-dish studies, more moderate systems are being readily discovered and are believed to make up the bulk of the population of these highly obscured galaxies (Fujimoto et al. 2016; Oteo et al. 2016; Franco et al. 2018; González-López et al. 2020; Zavala et al. 2021). Interestingly, Valentino et al. (2020) suggest galaxies with lower (sub)millimeter fluxes ($S_{850} < 3.5$ mJy) better reproduce the $M_*$ and SFR distributions expected for the progenitors of high-$z$ quiescent galaxies, further supporting our interest in uncovering such systems. The extremely dust-obscured sources (undetected in the deepest $Ks$-band surveys) are proposed to make up 20% of all DSFGs according to Dudzevičiūtė et al. (2020), which is in line with our results from the MORA survey as 23% of the detected DSFGs are OIR-dark.

Wang et al. (2019) reported 39 OIR-dark DSFGs (selected as H-band dropouts) at $z_{phot} > 3$. They were detected via ALMA at observed-frame 870 µm and determined to be DSFGs via their robust millimeter detections. The median stellar mass of this population is $M_* \sim 10^{10.6} M_{\odot}$ with a characteristic IR luminosity of $L_{IR(8-1000\mu m)} = (2.2 \pm 0.3) \times 10^{12} L_{\odot}$ and SFRs $\sim 200 M_{\odot}$ yr$^{-1}$ as determined via stacking. Similarly, Williams et al. (2019) reported serendipitous detections of two previously unknown $z > 3$ sources found in a small 8 arcmin$^2$ 3 mm survey in the COSMOS field, one of which being an “ALMA-only” source (named 3MM-1) with $z_{phot} = 5.5_{-1.1}^{+1.2}$ and OIR counterparts only detected ($\sim 3\sigma$) in deep, stacked IRAC 3.6 µm and 4.5 µm observations (see also Zavala 2021, which reports a tentative CO(6-5) line indicative of a $z = 5.857$ redshift solution). Both sources have SFRs akin to “normal” main sequence (MS) galaxies (a few $100 M_{\odot}$ yr$^{-1}$) and are in line with the H-band dropouts on the MS. OIR-dark DSFGs have also been identified starting from deep radio surveys. For example, Talia et al. (2021) present 197 VLA sources ($S_{3GHz} > 12.65$ µJy beam$^{-1}$) which do not have a COSMOS2015 counterpart, though we note MORA-5 and MORA-9 would not be included given this selection criteria.

Based on our current observations and SFR calculations, MORA-9 appears to belong to the moderately star-forming subsample of OIR-dark DSFGs, while MORA-5 displays physical characteristics consistent with other rare, extreme OIR-dark starbursts. Figure 7 shows SFR versus $M_*$ of our results of MORA-5 and MORA-9 (blue and teal diamonds) along with MORA-4, a.k.a. MAMBO-9, a spectroscopically confirmed OIR-dark DSFG (see Casey et al. 2019). The H-band dropouts (open circles), 3 mm detected OIR-dark DSFGs (filled squares), and $z > 3$ DSFGs (orange triangles) from the ALESS survey (da Cunha et al. 2015) are shown for comparison. Stellar masses of the ALESS sample are reduced by 0.3 dex in Wang et al. (2019) (depicted by the orange arrow) to account for different assumptions in their respective mass estimations. Since stellar masses for the MORA sources are estimated in the same manner as the ALESS sample (through MAGPHYS), we do not apply this offset in our figure. To compare apples to apples, the H-band dropout sample should

<figure><img src="image_10.png" /><figcaption>**Figure 7.** SFR vs. $M_*$ for MORA-5 (teal diamond) and MORA-9 (blue diamond) given their $z_{PDF}$ photometric redshifts. MORA-4 (a.k.a. MAMBO-9) is shown for comparison (small maroon diamond) as one of the other OIR-dark DSFGs detected in the MORA survey. MORA-3 (a.k.a. AzTEC-2) is not depicted here due to its lack of stellar mass estimate, a result of its IRAC photometry severely blending with foreground sources. We plot our results on top of H-band dropouts from Wang et al. (2019) (open circles) and 3 mm detected OIR-dark DSFGs (Williams et al. 2019) (filled squares) at similar redshifts from the literature. We also show $z > 3$ ALESS DSFGs (da Cunha et al. 2015) (orange triangles) and the MS of star formation at $z = 3-5$ Schreiber et al. (2015) (gray region). The length of the orange arrow depicts the 0.3 dex offset in $M_*$ applied to the ALESS sample in Wang et al. (2019) due to the different methods used to calculate stellar masses.</figcaption></figure>

have a factor of $\sim 2$ higher mass than shown. These results are all plotted on top of the MS of galaxies at $z = 3-5$ (filled gray region) as parameterized in Schreiber et al. (2015).

## 5.2. The Dynamic Range of OIR-dark DSFGs

This work, along with those mentioned in Section 5.1, illustrates the heterogeneous nature of OIR-dark DSFGs. The MORA survey, with its relatively large coverage area ($184 \text{ arcmin}^{-2}$ in total) for an interferometric survey, has shown to be effective at detecting both rarer, bright OIR-dark DSFGs with high submillimeter fluxes as well as fainter sources—a unique aspect of the survey design. We defined OIR-dark to simply include any galaxy (likely a DSFG) lacking detections at wavelengths shortward of 2.2 µm, and it is clear this encompasses a broad range of star formation and stellar mass.

While we lack some clarity on why some DSFGs are OIR-dark and others are not, it seems likely that a combination of factors are at play. Perhaps the most evident component contributing to the OIR-dark nature of DSFGs is the high redshifts of these sources compared to other DSFGs with comparable (sub)millimeter fluxes across $S_{870} = 3-10$ mJy. Smail et al. (2021) find that K-faint (K-band magnitude $>25.3$) DSFGs from the AS2UDS survey (Dudzevičiūtė et al. 2020) have a significantly higher median photometric redshift compared to their brighter K-detected sample: $z = 3.44 \pm 0.06$ versus $z = 2.36 \pm 0.11$, respectively. This work also shows a higher median dust attenuation ($A_V = 5.2 \pm 0.4$) for the K-faint population compared to their redshift-matched control sample ($A_V = 2.9 \pm 0.1$), which is attributed to the smaller dust continuum sizes observed (i.e., higher specific SFR, $\Sigma$SFR).
Finally, Smail et al. (2021) point to the mixture and relative scales of the obscured versus unobscured components, suggesting OIR-dark DSFGs may exhibit an absence of spatially extended, less obscured star formation. The lower $A_V$ value for MORA-9 and its associated K-band detection suggests it has a stellar component which is geometrically distinct from the dust emission, whereas those same components for MORA-5 may be coupled.

### 5.3. Do OIR-dark DSFGs Host Active Galactic Nuclei?

We searched for X-ray emission in the Chandra COSMOS Legacy data to infer the presence of AGN in MORA-5 and MORA-9, but neither source is detected. However, given the suspected high redshifts of these sources and the depth of the X-ray observations, the current data cannot rule out the presence of an AGN.

The IR-to-1.4 GHz radio luminosity ratio, $q_{IR}$, quantifies the far-IR/radio correlation (FRC). Various works have looked to determine if a redshift evolution of $q_{IR}$ exists as well as how it may differ in star-forming galaxies (SFGs) compared to AGN (Magnelli et al. 2012; Delhaize et al. 2017). Assuming a synchrotron slope of $\alpha = -0.8$ (Condon 1992), the marginal (3$\sigma$) 3 GHz detection of MORA-5 implies a lower limit of $q_{IR} > 1.9$, which is in agreement with what is predicted by the evolution of $q_{IR}$ for a $z=4.3$ SFG ($q_{IR} = 2.1$; Delhaize et al. 2017). Examining the radio emission of MORA-9 produces similar results and we find a lower limit of $q_{IR} > 1.6$.

When examined separately, AGN are found to have lower $q_{IR}$ values overall, and at this redshift a galaxy with a luminous AGN could be expected to have $q_{IR} < 1.5$ given the evolution of the FRC found by Delhaize et al. (2017). Even with a marginal radio detection in MORA-5, the uncertainty in the source's redshift makes it impossible to conclude whether or not an AGN is present; for example, it may sit at lower redshift with radio emission from star formation or at higher redshift with a radio-loud AGN. Confirming the redshifts of the two sources via an emission line search will be critical in determining whether or not they host luminous AGN.

### 5.4. Potential Evolutionary Tracks

The discovery of OIR-dark DSFGs at $z>3$ has now been evolutionarily linked to the recent detections of high-redshift (up to $z \sim 4$), massive ($M_* \sim 10^{11} M_\odot$) quiescent galaxies (Spitler et al. 2014; Glazebrook et al. 2017; Schreiber et al. 2018; Valentino et al. 2020; Stevans et al. 2021). The prevalence of OIR-dark galaxies in the early universe has large implications for the cosmic SFR density, stellar mass function (SMF), and theoretical models which have broadly excluded this population of galaxies due to their poor observational constraints (Dudzevičiūtė et al. 2020). These discoveries also drastically compress the formation and quenching timescales of quiescent galaxies to just a few billion years and raise questions regarding the galaxies they evolved from. While a handful of DSFGs have been discovered out to $z \sim 7$ (Marrone et al. 2018), they represent rare systems and are not thought to make up the bulk of the population of DSFGs at $z > 3$. Ultimately, the number density of massive, high-redshift galaxies required to reproduce the observed population of quiescent galaxies at $z=2-4$ is not supported by the summation of $z > 4$ massive LBGs and known extreme DSFGs alone (Straatman et al. 2014). To assess the viability of

OIR-dark DSFGs as the progenitors of high-$z$ massive quiescent galaxies, we follow arguments based on number density and stellar mass evolution.

#### 5.4.1. Space Density of OIR-dark DSFGs

In total, four OIR-dark DSFGs were detected in the MORA survey. Along with MORA-5 and MORA-9, two additional sources which already have spectroscopic confirmations, MORA-3 (a.k.a. AzTEC-2; Jiménez-Andrade et al. 2020) and MORA-4 (a.k.a. MAMBO-9; Casey et al. 2019), are folded into our source density calculations. We rely on these spectroscopic redshifts (MORA-3: $z=4.63$ and MORA-4: $z=5.85$) to aid in our estimate of the upper bound of the selection volume, ultimately adopting $z=6$. Nominally, the maximum redshift of a 2 mm detected DSFG could be much higher given the negative K-correction, even after accounting for CMB heating. Realistically, though, we know sources are extremely rare at $z>6$ based on current surveys and our understanding of the IR luminosity function (Zavala et al. 2021). We then set the lower bound to $z=4$, given our photometric redshifts for MORA-5 and MORA-9. Examining the contiguous 156 arcmin$^2$ area where all four sources were detected across $4 < z < 6$ results in a number density of $0.03 \text{ arcmin}^{-2}$ and comoving volume density of $n \sim (5 \pm 2) \times 10^{-6} \text{ Mpc}^{-3}$. The uncertainty on the volume density is generated by running MC trials to create the distribution of expected sources given both cosmic variance (35% between $3 < z < 6.5$ as determined in Zavala et al. 2021) and Poisson noise. These distributions are then summed and the resulting $1\sigma$ errors are reported. The area of the survey is sufficiently large, despite the small sample size, for this first-order estimate given that cosmic variance only begins to dominate beyond $z > 6.5$ (Zavala et al. 2021). A larger survey area, however, would of course improve our determination of number density.

This result is $\sim 4 \times$ lower than the total H-band dropout space density reported in Wang et al. (2019) ($n \sim 2 \times 10^{-5} \text{ Mpc}^{-3}$), which is said to be in agreement with that of massive quiescent galaxies found at $z \sim 3$ (Straatman et al. 2014). We note that our sample is ALMA selected, whereas $\sim 37\%$ of the H-band dropouts are ALMA undetected. Limiting the calculation to H-band dropouts with ALMA detections lowers the volume density of the Wang et al. (2019) sample by a factor of 2. Furthermore, if we look only at the $20 \ 4 < z < 6$ sources in the sample, the density drops down to $n \sim 6 \times 10^{-6} \text{ Mpc}^{-3}$ and falls in line with what we observe in MORA across this redshift range.

Whether or not a source is observed to be OIR-dark, all DSFGs are expected to quench and evolve into passive elliptical galaxies. Indeed, being OIR-dark is likely very dependent on redshift and geometry as discussed in Section 5.2. With this in mind, we also calculate the volume density of all 2 mm detected DSFGs in the MORA sample from $3 < z < 6$ and find $n \sim 1 \times 10^{-5} \text{ Mpc}^{-3}$. While it is reassuring to see agreement across different samples, the measurements of spatial density presented in the literature as well as this work both suffer from uncertainties in the comoving volume due to the reliance on photometric redshifts, highlighting the need for future spectroscopic follow-up of high-$z$ DSFGs.
### 5.4.2. Stellar Mass Growth

Both MORA-5, $M_\star = (1.5_{-0.7}^{+1.0}) \times 10^{11} M_\odot$, and MORA-9, $M_\star = (4.1_{-1.4}^{+1.8}) \times 10^{10} M_\odot$, have stellar masses already consistent with $z \sim 3$ massive quiescent galaxies (Schreiber et al. 2018), albeit on the lower mass end for MORA-9. High-z DSFGs are commonly thought to undergo a short ($\sim$50–150 Myr), bursty period of star formation mediated by strong galactic feedback (Mihos & Hernquist 1996; Cox et al. 2008; Swinbank et al. 2014; Toft et al. 2014; Aravena et al. 2016; Spilker et al. 2020), potentially triggered by an initial gas-rich major merger or interaction. Conversely, consistent gas infall may lead to significantly longer duty cycles of up to a gigayear (Hayward et al. 2013; Narayanan et al. 2015). Based on the SFRs and gas masses estimated in Section 4, MORA-5 and MORA-9 have star-forming gas depletion timescales of $\tau_{depl} \sim 310$ Myr and $\tau_{depl} \sim 800$ Myr, respectively. This is significantly longer than DSFGs at $z \sim 2–3$. Assuming these gas depletion times, the stellar mass of MORA-5 increases to $\sim 4 \times 10^{11} M_\odot$ by $z=3.6$ and the stellar mass of MORA-9 increases to $\sim 2 \times 10^{11} M_\odot$ by $z=2.9$. These evolved masses are still in line with the mass range of the quiescent sample, however it should be noted that we cannot exclude the possibility of accretion refilling the gas reservoirs and extending the depletion timescales. Valentino et al. (2020) (Figure 6 therein) illustrate the wide range in predicted number densities of quiescent galaxies with $\log(M_\star/M_\odot) > 10.6$ from both observations and theory across $3 \lesssim z \lesssim 4$. This mass threshold is well aligned with the entire MORA sample for which we expect to produce galaxies with stellar masses in excess of $\log(M_\star/M_\odot) > 10.6$. (Casey et al. 2021, Figure 13). For the COSMOS field specifically, the SMF determined in Davidzon et al. (2017) produces a number density of $n \sim 2 \times 10^{-6}$ Mpc$^{-3}$, which then drops precipitously as the mass threshold increases. The lack of consensus in the number density of $z > 3$ quiescent galaxies speaks to the need for future observations with the JWST to detect and constrain this galaxy population. Luckily, this is a science goal for several of the approved Cycle 1 proposals and, as such, we eagerly await those results.

If we consider the current stellar mass estimates alone in this closed-box scenario, the MORA OIR-dark DSFGs are viable progenitors to high-z quiescent galaxies. For now, our conjecture stops there until more accurate stellar masses and greater understanding of potential feedback and quenching mechanisms for these two OIR-dark DSFGs can be obtained. Future observations with JWST Mid-Infrared Instrument (MIRI) imaging would provide constraints on the rest-frame 1.6 µm stellar bump, thus elucidating the stellar masses of these heavily obscured systems.

## 6. Conclusions

We present photometric redshifts and physical characterization of MORA-5 and MORA-9, two OIR-dark DSFGs detected in the 2 mm MORA survey. These sources, of interest for their potential to be the highest-redshift objects in MORA, lacked secure detections in any photometric bands short of 2.2 µm, consistent with other OIR-dark systems reported in recent literature.

The photometric redshifts reported in this work, $z_{PDF} = 4.3_{-1.3}^{+1.5}$ for MORA-5 and $z_{PDF} = 4.3_{-1.0}^{+1.3}$ for MORA-9, are the result of combining the redshift distributions determined via

OIR, FIR/millimeter, and energy balance SEDs. It is this $z_{PDF}$ value which informs the physical properties derived for these two systems with the exception of stellar mass and $A_V$. These two properties are instead determined by MAGPHYS +photo-z, given its inclusion of IR, (sub)millimeter, and radio data, which is especially beneficial for high-z, dust-obscured galaxies.

Based on our current observations, MORA-5 is the more active of the two systems. Its high stellar mass, $1.5_{-0.7}^{+1.0} \times 10^{11} M_\odot$, and star formation rate, SFR $\approx 830 M_\odot$ yr$^{-1}$, suggests it is part of a class of rarer, more extreme OIR-dark galaxies. MORA-9 is modest in comparison with $(4.1_{-1.4}^{+1.8}) \times 10^{10} M_\odot$ and SFR $\approx 200 M_\odot$ yr$^{-1}$. From our 2 mm dust continuum observations we determine MORA-5 and MORA-9 to have gas masses of $M_{gas} = (2.6 \pm 0.4) \times 10^{11} M_\odot$ and $M_{gas} = (1.6 \pm 0.3) \times 10^{11} M_\odot$ and gas depletion timescales of $\tau_{depl} \sim 310$ Myr and $\tau_{depl} \sim 800$ Myr, respectively. Evolving the stellar masses according to our estimated gas depletion timescales results in both systems remaining consistent with the masses of known $z \sim 3$ quiescent galaxies.

It is promising to see convergence toward $z > 3$ solutions in all the redshift PDFs given the decision to observe at 2 mm to effectively weed out $z \lesssim 2.5$ interlopers thanks to the very negative K-correction at millimeter wavelengths. However, we still highlight the need for spectroscopic follow-up of OIR-dark DSFGs to confirm redshifts and allow for accurate calculations of their spatial density and contribution to the cosmic SFR density as a function of redshift. Secure redshifts are also essential for placing constraints on dust production mechanisms and buildup within the first billion years after the Big Bang. Redshift confirmation via spectral scan to search for CO transitions with ALMA could be achieved in less than eight hours for both sources combined, taking into account the required rms sensitivity determined by the $L_{IR}-L_{CO}$ relation (Bothwell et al. 2013; Greve et al. 2014).

Given the uncertainties in redshift (and, thus, SFR) and stellar masses of OIR-dark DSFGs, it is difficult to compare them directly with high-z quiescent galaxies. Based on our current understanding of the space density of both populations and estimates of stellar mass growth, our work is in agreement with previous studies and supports the notion that OIR-dark DSFGs are the progenitors of $z=2-4$, massive, quiescent galaxies. Resolved (sub)millimeter observations providing effective radius measurements of OIR-dark DSFGs are necessary for determining an evolutionary link to compact quiescent galaxies in terms of their comparative spatial extents. They would also allow for measurements of the average dust column density (alluding to the ISM conditions in OIR-dark DSFGs) as well as star formation surface density. Future observations, either with JWST NIRSpec campaigns or emission line surveys (e.g., with ALMA, NOEMA, LMT), would establish accurate redshifts, while JWST NIRCam and MIRI imaging would constrain the stellar emission and offer insights into the morphologies of these systems. Finally, we will look to expand with wide-field 2 mm imaging as an effective strategy to detect high-z, OIR-dark DSFGs.

We thank the anonymous reviewer for their thoughtful suggestions which helped improve this manuscript. S.M.M. thanks Isabella Cortzen, Seiji Fujimoto, and Francesco Valentino for their hospitality and helpful discussions while visiting the Cosmic DAWN Center. S.M.M. also thanks her
parents who provided a home, companionship, and a desk to write this paper during the lockdown brought on by the COVID-19 pandemic. Finally, S.M.M. thanks the National Science Foundation for support through the Graduate Research Fellowship under grant No. DGE-1610403, the NSF GROW travel grant, Georgios Magdis for being her site host at DAWN, and the NRAO SOS grant. Support for this work was also provided by NASA through the NASA Hubble Fellowship grant No. \#HST-HF2-51484 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. C.M.C. thanks the National Science Foundation for support through grant Nos. AST-1814034 and AST-2009577 and thanks RCSA for a Cottrell Scholar Award, sponsored by IF/THEN, an initiative of the Lyda Hill Philanthropies. G.E.M. acknowledges the Villum Fonden research grant No. 37440, “The Hidden Cosmos and the Cosmic Dawn” Center of Excellence funded by the Danish National Research Foundation under then grant No. 140. M.A. acknowledges support from FONDECYT grant No. 1211951, “CONICYT + PCI + INSTITUTO MAX PLANCK DE ASTRONOMIA MPG190030” and “CONICYT+PCI+REDES 190194”. K. K. acknowledges support from the Knut and Alice Wallenberg Foundation. M.T. acknowledges the support from grant PRIN MIUR 2017. This paper makes use of the following ALMA data: ADS/JAO.ALMA\#2019.2.00143.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), da Cunha, E., Walter, F., Smail, I. R., et al. 2015, ApJ, 806, 110 Dahlen, T., Mobasher, B., Faber, S. M., et al. 2013, ApJ, 775, 93 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 Donevski, D., Lapi, A., Małek, K., et al. 2020, A&A, 644, A144 Dudzevičiūtė, U., Smail, I., Swinbank, A. M., et al. 2020, MNRAS, 494, 3828 Finkelstein, S. L., Ryan, R. E. J., Papovich, C., et al. 2015, ApJ, 810, 71 Fogasy, J., Knudsen, K. K., Drouart, G., Lagos, C. D. P., & Fan, L. 2020, MNRAS, 493, 3744 Fogasy, J., Knudsen, K. K., Lagos, C. D. P., Drouart, G., & Gonzalez-Perez, V. 2017, A&A, 597, A123 Franco, M., Elbaz, D., Béthermin, M., et al. 2018, A&A, 620, A152 Fujimoto, S., Ouchi, M., Ono, Y., et al. 2016, ApJS, 222, 1 Geach, J. E., Dunlop, J. S., Halpern, M., et al. 2017, MNRAS, 465, 1789 Glazebrook, K., Schreiber, C., Labbé, I., et al. 2017, Natur, 544, 71 González-López, J., Novak, M., Decarli, R., et al. 2020, ApJ, 897, 91 Greve, T. R., Leonidaki, I., Xilouris, E. M., et al. 2014, ApJ, 794, 142 Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35 Gruppioni, C., Béthermin, M., Loiacono, F., et al. 2020, A&A, 643, A8 Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 432, 23 Hainline, L. J., Blain, A. W., Smail, I., et al. 2011, ApJ, 740, 96 Hao, C.-N., Kennicutt, R. C., Johnson, B. D., et al. 2011, ApJ, 741, 124 Hayward, C. C., Narayanan, D., Kereš, D., et al. 2013, MNRAS, 428, 2529 Hughes, D. H., Serjeant, S., Dunlop, J., et al. 1998, Natur, 394, 241 Ilbert, O., Arnouts, S., McCracken, H. J., et al. 2006, A&A, 457, 841 Jiménez-Andrade, E. F., Zavala, J. A., Magnelli, B., et al. 2020, ApJ, 890, 171 Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531 Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36 Koprowski, M. P., Dunlop, J. S., Michałowski, M. J., et al. 2017, MNRAS, 471, 4155 Kroupa, P., & Weidner, C. 2003, ApJ, 598, 1076 Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24 Le Floc'h, E., Aussel, H., Ilbert, O., et al. 2009, ApJ, 703, 222 Lutz, D., Poglitsch, A., Altieri, B., et al. 2011, A&A, 532, A90 Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 Magnelli, B., Lutz, D., Santini, P., et al. 2012, A&A, 539, A155 Marrone, D. P., Spilker, J. S., Hayward, C. C., et al. 2018, Natur, 553, 51 McCracken, H. J., Milvang-Jensen, B., Dunlop, J., et al. 2012, A&A, 544, A156 McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser., 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

Software: CASA (v5.6.1; McMullin et al. 2007), EAZY (v1.3; Brammer et al. 2008), MAGPHYS (da Cunha et al. 2008), MAGPHYS+photoz (Battisti et al. 2019), LePhare (Arnouts et al. 2002; Ilbert et al. 2006).

## ORCID iDs

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

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

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

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

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

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

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

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

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

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

Olivier Ilbert https://orcid.org/0000-0002-7303-4397

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

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

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

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

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

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

Justin S. 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

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

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

## References

Aravena, M., Bertoldi, F., Carilli, C., et al. 2010, ApJL, 708, L36

Aravena, M., Spilker, J. S., Bethermin, M., et al. 2016, MNRAS, 457, 4406

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

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

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

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

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

Bothwell, M. S., Smail, I., Chapman, S. C., et al. 2013, MNRAS, 429, 3047

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

Brammer, G. B., van Dokkum, P. G., Franx, M., et al. 2012, ApJS, 200, 13

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

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., Hodge, J., Zavala, J. A., et al. 2018a, ApJ, 862, 78

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

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

Casey, C. M., Zavala, J. A., Manning, S. M., et al. 2021, ApJ, 923, 215

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

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

Condon, J. J. 1992, ARA&A, 30, 575

Cox, T. J., Jonsson, P., Somerville, R. S., Primack, J. R., & Dekel, A. 2008, MNRAS, 384, 386

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

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

Dahlen, T., Mobasher, B., Faber, S. M., et al. 2013, ApJ, 775, 93

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

Donevski, D., Lapi, A., Małek, K., et al. 2020, A&A, 644, A144

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

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

Fogasy, J., Knudsen, K. K., Drouart, G., Lagos, C. D. P., & Fan, L. 2020, MNRAS, 493, 3744

Fogasy, J., Knudsen, K. K., Lagos, C. D. P., Drouart, G., & Gonzalez-Perez, V. 2017, A&A, 597, A123

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

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

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

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

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

Greve, T. R., Leonidaki, I., Xilouris, E. M., et al. 2014, ApJ, 794, 142

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

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

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

Hainline, L. J., Blain, A. W., Smail, I., et al. 2011, ApJ, 740, 96

Hao, C.-N., Kennicutt, R. C., Johnson, B. D., et al. 2011, ApJ, 741, 124

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

Hughes, D. H., Serjeant, S., Dunlop, J., et al. 1998, Natur, 394, 241

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

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

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

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

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

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

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

Le Floc'h, E., Aussel, H., Ilbert, O., et al. 2009, ApJ, 703, 222

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

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

Magnelli, B., Lutz, D., Santini, P., et al. 2012, A&A, 539, A155

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

McCracken, H. J., Milvang-Jensen, B., Dunlop, J., et al. 2012, A&A, 544, A156

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser., 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127
Miettinen, O., Delvecchio, I., Smolčić, V., et al. 2017, A&A, 606, A17

Mihos, J. C., & Hernquist, L. 1996, ApJ, 464, 641

Murphy, E. J., Condon, J. J., Schinnerer, E., et al. 2011, ApJ, 737, 67

Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJS, 206, 8

Narayanan, D., Turk, M., Feldmann, R., et al. 2015, Natur, 525, 496

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

Ono, Y., Ouchi, M., Harikane, Y., et al. 2018, PASJ, 70, S10

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

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

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

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

Sanders, D. B., & Mirabel, I. F. 1996, ARA&A, 34, 749

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

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

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

Scott, S. E., Fox, M. J., Dunlop, J. S., et al. 2002, MNRAS, 331, 817

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

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

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

Smail, I., Dudzevičiūtė, U., Stach, S. M., et al. 2021, MNRAS, 502, 3426

Smith, D. J. B., Hayward, C. C., Jarvis, M. J., & Simpson, C. 2017, MNRAS, 471, 2453

Smolčić, V., Miettinen, O., Tomičić, N., et al. 2017b, A&A, 597, A4

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

Spilker, J. S., Aravena, M., Phadke, K. A., et al. 2020, ApJ, 905, 86

Spitler, L. R., Straatman, C. M. S., Labbé, I., et al. 2014, ApJL, 787, L36

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

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

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

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

Talia, M., Cimatti, A., Giulietti, M., et al. 2021, ApJ, 909, 23

Toft, S., Smolčić, V., Magnelli, B., et al. 2014, ApJ, 782, 68

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

Walter, F., Decarli, R., Carilli, C., et al. 2012, Natur, 486, 233

Wang, T., Schreiber, C., Elbaz, D., et al. 2019, Natur, 572, 211

Weaver, J. R., Kauffmann, O. B., Ilbert, O., et al. 2022, ApJS, 258, 11

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

Zavala, J. A. 2021, RNAAS, 5, 15

Zavala, J. A., Aretxaga, I., & Hughes, D. H. 2014, MNRAS, 443, 2384

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