# Two Massive, Compact, and Dust-obscured Candidate z 8 Galaxies Discovered by JWST

Hollis B. Akins, Caitlin M. Casey, Natalie Allen, Micaela B. Bagley, Mark Dickinson, Steven L. Finkelstein, Maximilien Franco, Santosh Harish, Pablo Arrabal Haro, Olivier Ilbert, et al.

► To cite this version:

Hollis B. Akins, Caitlin M. Casey, Natalie Allen, Micaela B. Bagley, Mark Dickinson, et al. Two Massive, Compact, and Dust-obscured Candidate z 8 Galaxies Discovered by JWST. The Astrophysical Journal, 2023, 956, 10.3847/1538-4357/acef21. insu-04479051

HAL Id: insu-04479051

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

Submitted on 10 Mar 2024

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
Two Massive, Compact, and Dust-obscured Candidate $z \approx 8$ Galaxies Discovered
by JWST

Hollis B. Akins¹, Caitlin M. Casey¹,², Natalie Allen²,³, Micaela B. Bagley¹, Mark Dickinson⁴, Steven L. Finkelstein¹, Maximilien Franco¹, Santosh Harish⁵, Pablo Arrabal Haro⁴, Olivier Ilbert⁶, Jeyhan S. Kartaltepe⁵, Anton M. Koekemoer⁷, Daizhong Liu⁸, Arianna S. Long¹,⁴,¹⁰, Henry Joy McCracken⁹, Louise Paquereau⁹, Casey Papovich¹⁰,¹¹, Nor Pirzkal¹², Jason Rhodes¹³, Brant E. Robertson¹⁴, Marko Shuntov²,³, Sune Toft²,³, Guang Yang¹⁵,¹⁶, Guillermo Barro¹⁷, Laura Bisigello¹⁸,¹⁹, Véronique Buat⁶, Jaclyn B. Champagne²⁰, Olivia Cooper¹,⁴², Luca Costantin²¹, Alexander de la Vega²², Nicole E. Drakos¹⁴, Andreas Faisst²³, Adriano Fontana²⁴, Seiji Fujimoto¹,⁴¹, Steven Gillman²,²⁵, Carlos Gómez-Guijarro²⁶, Ghassem Gozaliasl²⁷, Nimish P. Hathi⁷, Christopher C. Hayward²⁸, Michaela Hirschmann²⁹,³⁰, Benne W. Holwerda³¹, Shuowen Jin²,²⁵, Dale D. Kocevski³², Vasily Kokorev¹⁵, Erini Lambrides³³,⁴³, Ray A. Lucas⁷, Georgios E. Magdis²,²,²⁵, Benjamin Magnelli²⁶, Jed McKinney¹, Bahram Mobasher²², Pablo G. Pérez-González²¹, R. Michael Rich³⁴, Lise-Marie Seillé⁶, Margherita Talia³⁵,³⁶, C. Megan Urry³⁷, Francesco Valentino²,³⁸, Katherine E. Whitaker²,³⁹, L. Y. Aaron Yung³³,⁴⁴, and Jorge Zavala⁴⁰
the COSMOS-Web and CEERS teams

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

² Cosmic Dawn Center (DAWN), Denmark

³ Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200, Copenhagen N, Denmark

⁴ NSF's National Optical-Infrared Astronomy Research Laboratory, 950 North Cherry Avenue, Tucson, AZ 85719, USA

⁵ Laboratory for Multiwavelength Astrophysics, School of Physics and Astronomy, Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, NY 14623, USA

⁶ Aix Marseille Université, CNRS, CNES, LAM, Marseille, France

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

⁸ Max-Planck-Institut für Extraterrestrische Physik (MPE), Giessenbachstr. 1, D-85748 Garching, Germany

⁹ Institut d'Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Université, 98 bis boulevard Arago, F-75014 Paris, France

¹⁰ Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843-4242, USA

¹¹ George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, Texas A&M University, College Station, TX 77843-4242, USA

¹² ESA/AURA Space Telescope Science Institute, USA

¹³ Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91001, USA

¹⁴ Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA

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

¹⁶ SRON Netherlands Institute for Space Research, Postbus 800, 9700 AV Groningen, The Netherlands
¹⁷ Department of Physics, University of the Pacific, Stockton, CA 90340, USA

¹⁸ INAF-Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, I-35122, Padova, Italy

¹⁹ Dipartimento di Fisica e Astronomia "G.Galilei," Università di Padova, Via Marzolo 8, I-35131 Padova, Italy
²⁰ Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA

²¹ Centro de Astrobiología (CAB/CSIC-INTA), Ctra. de Ajalvir km 4, Torrejón de Ardoz, E-28850 Madrid, Spain

²² Department of Physics and Astronomy, University of California Riverside, 900 University Avenue, Riverside, CA 92521, USA
²³ Caltech/IPAC, 1200 East California Boulevard, Pasadena, CA 91125, USA

²⁴ INAF, Osservatorio Astronomico di Roma, via di Frascati 33, I-00078 Monte Porzio Catone, Italy
²⁵ DTU-Space, Technical University of Denmark, Elektrovej 327, DK-2800 Kgs. Lyngby, Denmark

²⁶ Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, F-91191, Gif-sur-Yvette, France
²⁷ Department of Physics, University of Helsinki, P.O. Box 64, FI-00014 Helsinki, Finland

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

²⁹ Institute for Physics, Laboratory for Galaxy Evolution and Spectral Modelling, Ecole Polytechnique Federale de Lausanne, Observatoire de Sauverny, Chemin Pegasi 51, 1290 Versoix, Switzerland

³⁰ INAF, Osservatorio Astronomico di Trieste, Via Tiepolo 11, I-34131 Trieste, Italy

³¹ Physics & Astronomy Department, University of Louisville, Louisville, KY 40292, USA

³² Department of Physics and Astronomy, Colby College, Waterville, ME 04901, USA

³³ NASA Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771, USA

³⁴ Department of Physics and Astronomy, UCLA, PAB 430 Portola Plaza, Box 951547, Los Angeles, CA 90095-1547, USA

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

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

³⁷ Physics Department and Yale Center for Astronomy & Astrophysics, Yale University, P.O. Box 208120, New Haven, CT 06520, USA

³⁸ European Southern Observatory, Karl-Schwarzschild-Strasse 2, D-85748, Garching bei München, Germany

³⁹ Department of Astronomy, University of Massachusetts, Amherst, MA 01003, USA
⁴⁰ National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
Received 2023 April 23; revised 2023 July 25; accepted 2023 July 31; published 2023 October 6

## Abstract

We present a search for extremely red, dust-obscured, $z > 7$ galaxies with JWST/NIRCam+MIR imaging over the first 20 arcmin² of publicly available Cycle 1 data from the COSMOS-Web, CEERS, and PRIMER surveys. Based on their red color in F277W–F444W (~2.5 mag) and detection in MIRI/F770W (~25 mag), we identify two galaxies, COS-z8M1 and CEERS-z7M1, that have best-fit photometric redshifts of $z = 8.4_{-0.4}^{+0.3}$ and $7.6_{-0.1}^{+0.1}$, respectively. We perform spectral energy distribution fitting with a variety of codes (including BAGPIPES, PROSPECTOR, BEAGLE, and CIGALE) and find a >95% probability that these indeed lie at $z > 7$. Both sources are compact ($R_{\text{eff}} \lesssim 200$ pc) and highly obscured ($A_V \sim 1.5–2.5$) and, at our best-fit redshift estimates, likely have strong [O III]+Hβ emission contributing to their $4.4 \ \mu$m photometry. We estimate stellar masses of $\sim 10^{10} M_\odot$ for both sources; by virtue of detection in MIRI at $7.7 \ \mu$m, these measurements are robust to the inclusion of bright emission lines, for example, from an active galactic nucleus. We identify a marginal (2.9σ) Atacama Large Millimeter/submillimeter Array detection at 2 mm within $0.5''$ of COS-z8M1, which, if real, would suggest a remarkably high IR luminosity of $\sim 10^{12} L_\odot$. These two galaxies, if confirmed at $z \sim 8$, would be extreme in their stellar and dust masses and may be representative of a substantial population of highly dust-obscured galaxies at cosmic dawn.

*Unified Astronomy Thesaurus concepts: Galaxy evolution (594); Galaxy formation (595); High-redshift galaxies (734); Galaxies (573)*

## 1. Introduction

The launch of JWST has immensely widened our view of the $z > 8$ Universe, pushing observations within a mere 600 Myr of the Big Bang. Within the first year of observations, dozens of candidate $z \sim 8–10$ galaxies have been identified based on photometric redshifts (Finkelstein et al. 2022; Naidu et al. 2022a; Donnan et al. 2023), and several have been spectroscopically confirmed with JWST/NIRSpec (e.g., Curtis-Lake et al. 2022; Arrabal Haro et al. 2023; Fujimoto et al. 2023a; Larson et al. 2023; Saxena et al. 2023; Tang et al. 2023). Virtually all of these candidates have been selected via the Lyman break (Steidel et al. 1996), a technique that is most effective for intrinsically blue and UV-luminous sources but fails to capture fainter, reddened objects.

Galaxy selection via the Lyman break in the rest-frame UV poses a particular challenge for measuring accurate stellar masses (e.g., Papovich et al. 2022) and placing the candidates in a cosmological context. The most massive galaxies in the early Universe can provide key constraints on physical models of galaxy formation and evolution (e.g., Boylan-Kolchin 2022; Ferrara et al. 2022a; Menci et al. 2022; Harikane et al. 2023a; Lovell et al. 2023) and will, by nature, be more chemically evolved and may have substantial dust reservoirs (Whitaker et al. 2017). Indeed, a strong candidate population for the most massive galaxies in the early Universe is the dusty star-forming galaxies (DSFGs), which are ubiquitous at cosmic noon ($z \sim 2$), where they dominate the star formation rate (SFR) density of the Universe (Casey et al. 2014a; Madau & Dickinson 2014). While DSFGs have been detected out to $z \sim 6–7$ with obscured SFRs in excess of $1000 M_\odot$ yr$^{-1}$

(Marrone et al. 2018; Zavala et al. 2018; Casey et al. 2019; Endsley et al. 2022b; Fujimoto et al. 2022a), their volume density at this epoch remains largely unconstrained due to the difficulty of constructing a complete sample. In particular, accurate photometric redshift estimates and systematic spectroscopic follow-up for faint DSFGs is challenging due to their faint rest-UV and optical emission (e.g., Talia et al. 2021) and degeneracies with lower-redshift galaxies (which are $\sim 20–100$ times more numerous; e.g., Smolčić et al. 2012).

Various techniques have been developed to constrain the population of dusty galaxies at $z \gtrsim 4–5$, including focusing on gravitationally lensed objects (e.g., Marrone et al. 2018; Zavala et al. 2018; Fujimoto et al. 2021) or observing at longer wavelengths (e.g., 2–3 mm) to efficiently filter out lower-redshift galaxies (e.g., Casey et al. 2021; Zavala et al. 2021; Cooper et al. 2022). However, at the highest redshifts ($z > 7$), these methods have thus far only proven sensitive to the brightest, most extreme ($L_{IR} \gtrsim 10^{12–13} L_\odot$) DSFGs, as confirming lower-luminosity sources at the highest redshifts is like searching for a “needle in a haystack,” with a relatively low yield at $z > 4$ (Casey et al. 2018a, 2018b). At present, most $z > 7$ dust continuum detections come from follow-up of UV-selected Lyman-break galaxies (LBGs; e.g., Watson et al. 2015; Fujimoto et al. 2022a; Bouwens et al. 2022) or serendipitous detections (e.g., Fudamoto et al. 2021) with only one robust dust continuum-detected object at $z > 8$: MACS 0416-Y1 at $z = 8.31$ (Tamura et al. 2019).<fnref n="45" /> It remains unclear whether the perceived rarity of DSFGs in the first $\sim 750$ Myr of the Universe is intrinsic, due to a gradual buildup of dust/stellar mass, or artificial, due to the difficulty of detecting and confirming these objects.

For the first time, the unprecedented sensitivity and continuous infrared wavelength coverage of JWST makes it possible to detect the rest-frame optical emission from extremely obscured galaxies at $z > 8$. While far-infrared (FIR) constraints at this epoch are limited, the near- and mid-

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="41">⁴¹ Hubble Fellow.</footnote>

<footnote n="42">⁴² NSF Graduate Research Fellow.</footnote>

<footnote n="43">⁴³ NPP Fellow.</footnote>

<footnote n="44">⁴⁴ NASA Postdoctoral Fellow.</footnote>

<footnote n="45">⁴⁵ Object A2744-YD4 had previously been reported at $z = 8.38$ (Laporte et al. 2017), but new JWST spectroscopic observations confirm its redshift as $z = 7.88$ (Morishita et al. 2022).</footnote>
infrared (MIR) have the potential to identify the most massive, obscured galaxies based on the reddened stellar continuum, rest-frame Balmer break, and contributions from bright emission lines (Pérez-González et al. 2022; Labbé et al. 2023a; Rodighiero et al. 2023). In particular, the MIR coverage of JWST/MIRI allows constraints on the full rest-frame optical spectral energy distribution (SED) out to $z \sim 9$. A unique advantage of long-wavelength MIRI imaging comes in identifying and characterizing the earliest dust-obscured galaxies, which are most heavily obscured at rest-frame UV wavelengths. These objects will be critical to our understanding of the assembly of massive galaxies (Narayanan et al. 2015; Long et al. 2022), as well as the physical processes responsible for the buildup of large dust reservoirs, whether it be from asymptotic giant branch (AGB) stars, supernovae (SNe), or efficient grain growth in the interstellar medium (ISM; Dwek & Cherchneff 2011; Jones et al. 2013; Michałowski 2015; Leśniewska & Michałowski 2019).

In this paper, we present a search for extremely red, $z \gtrsim 7$ galaxies across the overlapping NIRCam+MIRI coverage in three Cycle 1 Treasury surveys: COSMOS-Web (Casey et al. 2022), Cosmic Evolution Early Release Science (CEERS; Finkelstein et al. 2022), and Public Release Imaging for Extragalactic Research (PRIMER) COSMOS (PI: J. Dunlop; GO No. 1837). We report the detection of two galaxies, in CEERS and COSMOS-Web, with remarkably red colors and photometric redshifts at $z \sim 8$. Both are among the reddest galaxies identified in the entirety of their respective surveys and detected in 7.7 µm MIRI imaging, which provides robust constraints on their stellar mass.

The paper is organized as follows. In Section 2, we describe an overview of the imaging data used and the construction of photometric catalogs. In Section 3, we describe our sample selection and vetting of individual candidates. In Section 4, we present and characterize the two robust candidates, including their photometric redshifts (Section 4.1), sizes (Section 4.2), stellar masses (Section 4.3), and FIR emission (Section 4.4). Finally, in Section 5, we discuss the implications of these discoveries for galaxy formation within a ΛCDM framework. Throughout this work, we assume a Planck cosmology (Planck Collaboration 2020) and a Kroupa (2002) stellar initial mass function (IMF). All quoted magnitudes are in the AB system (Oke 1974).

## 2. Data

We utilize 1–8 µm JWST/NIRCam+MIRI imaging from three publicly available Cycle 1 programs to identify candidate massive, dusty galaxies at $z \gtrsim 8$. Table 1 provides the approximate $5\sigma$ point-source depths and the effective area (i.e., the NIRCam+MIRI overlapping area) for each survey.

### 2.1. JWST Observations and Data Reduction

#### 2.1.1. COSMOS-Web

COSMOS-Web is a large Cycle 1 Treasury program imaging a contiguous 0.54 deg² in the COSMOS field with NIRCam and 0.2 deg² with MIRI in parallel (Casey et al. 2022). As of this writing, six of the 152 visits had been completed during observations executed in 2023 January, constituting a contiguous 77 arcmin² (~4% of the overall area), with 8.7 arcmin² of overlap between the NIRCam and MIRI coverage. The COSMOS-Web imaging includes four NIRCam filters—

<table><caption>Table 1</caption><caption>Combined NIRCam+MIRI Area and Approximate 5σ Depths for the Three Surveys Included in This Work</caption><thead><tr><th>Survey</th><th>Area<br>(arcmin<sup>2</sup>)</th><th>F277W<br>(5σ)</th><th>F444W<br>(5σ)</th><th>F770W<br>(5σ)</th></tr></thead><tbody><tr><td>COSMOS-Web<sup>a</sup></td><td>8.7</td><td>28.3</td><td>28.2</td><td>26.0</td></tr><tr><td>CEERS</td><td>7.8</td><td>29.2</td><td>28.6</td><td>27.1</td></tr><tr><td>PRIMER-COSMOS<sup>a</sup></td><td>4.1</td><td>28.9</td><td>28.7</td><td>26.0</td></tr></tbody></table>

**Notes.** Depths are computed in $0''/3$ ($0''/6$) diameter apertures for NIRCam (MIRI).

<sup>a</sup> The table reflects the data available as of 2023 March. Total NIRCam+MIRI overlap upon survey completion will be 650 arcmin² for COSMOS-Web and 140 arcmin² for PRIMER-COSMOS.

**References.** COSMOS-Web (Casey et al. 2022; GO No. 1727), CEERS (Finkelstein et al. 2022; DD-ERS No. 1345), PRIMER (PI: J. Dunlop; GO No. 1837).

F115W, F150W, F277W, and F444W—and one MIRI filter, F770W, at an approximate $5\sigma$ depth of 26 AB mag.

The full details of the NIRCam and MIRI reduction process will be presented in upcoming papers (M. Franco et al. 2023, in preparation, and S. Harish et al. 2023, in preparation, respectively) but are briefly described here. The raw NIRCam imaging was reduced by the JWST Calibration Pipeline version 1.8.3, with the addition of several custom modifications (as has also been done for other JWST studies; e.g., Finkelstein et al. 2022), including the subtraction of 1/f noise and sky background. We use the Calibration Reference Data System<fnref n="46" /> pmap 0989, which corresponds to the NIRCam instrument mappingzap 0232. The final mosaics are created in stage 3 of the pipeline with a pixel size of $0''/03$ pixel⁻¹. Astrometric calibration is conducted via the JWST TWEAKREG procedure, with a reference catalog based on a Hubble Space Telescope (HST) F814W $0''/03$ pixel⁻¹ mosaic in the COSMOS field with astrometry tied to Gaia-EDR3 (Gaia Collaboration 2018). The median offset in R.A. and decl. between our reference catalog and the NIRCam mosaic is less than 5 mas. The MIRI/F770W observations were reduced using version 1.8.4 of the JWST Calibration Pipeline, along with additional steps for background subtraction that was necessary to mitigate the instrumental effects. The resulting mosaic was resampled onto a common output grid with a pixel scale of $0''/06$ pixel⁻¹ and aligned with ancillary HST/F814W imaging of the region.

#### 2.1.2. CEERS

CEERS is one of 13 early release science surveys designed to obtain and release reduced data in early Cycle 1. CEERS consists of a mosaic of 10 NIRCam and nine MIRI pointings in the CANDELS Extended Groth Strip (EGS) field, alongside spectroscopy with NIRSpec and NIRCam WFSS. Each CEERS/NIRCam pointing includes seven filters: F115W, F150W, F200W, F277W, F356W, F410M, and F444W. The MIRI pointings include a range of filters from F560W to F2100W; here we only use MIRI pointings 3, 6, 7, and 9, as the other pointings either have no NIRCam overlap (1 and 2) or do not include F770W (5 and 8), which is essential for the robust selection of our $z > 7$ dusty galaxy candidates. These four pointings provide 7.8 arcmin² of overlap between MIRI/

<footnote n="46">⁴⁶ jwst-crds.stsci.edu</footnote>
F770W and NIRCam imaging at an approximate $5\sigma$ depth in F770W of ~27 mag. We utilize the NIRCam and MIRI reductions produced by the CEERS team. The NIRCam reduction is described in detail in Bagley et al. (2022), and the MIRI reduction will be described in Yang et al. (2023).

### 2.1.3. PRIMER-COSMOS

The PRIMER survey (PI: J. Dunlop; GO No. 1837) is a large Cycle 1 Treasury Program to image two HST CANDELS Legacy Fields (COSMOS and UDS) with NIRCam+MIRI. PRIMER is conducted with MIRI as the prime instrument and NIRCam in parallel, with observations split between two windows with opposite observational position angles. This configuration maximizes the overlap between the MIRI and NIRCam coverage, though at the time of this writing, only the first epoch had been observed, constituting just 4.1 arcmin² of overlap between NIRCam+MIRI. PRIMER imaging includes eight NIRCam bands (equivalent to CEERS plus F090W) plus two MIRI bands (F770W and F1800W). As part of the reduction of COSMOS-Web data, the COSMOS-Web team has conducted an independent processing of the PRIMER data in the COSMOS field; thus, we include it in the analysis here. We note that while the PRIMER-COSMOS field is contained entirely within the COSMOS-Web footprint, there will not be significant overlap between the two surveys until the completion of COSMOS-Web in 2024 January.

## 2.2. Archival Ground-based and HST Data

In addition to the new JWST data, we utilize the existing multiwavelength imaging in the COSMOS and EGS fields. For COSMOS (encompassing PRIMER-COSMOS), we include the *grizy* imaging from the Subaru/Hyper Suprime-Cam (HSC; Aihara et al. 2019) and *YJHK*, imaging from the UltraVISTA survey (McCracken et al. 2012). We additionally utilize the HST Advanced Camera for Surveys (ACS) F814W imaging covering the entire COSMOS field (Koekemoer et al. 2007) and Spitzer/IRAC imaging from the Cosmic Dawn Survey (Euclid Collaboration et al. 2022). This is the same data as used in Weaver et al. (2022), with the exception of the UltraVISTA data, which are updated to the newest public data release 5.

For CEERS, we include the HST/ACS and WFC3 imaging in F606W, F814W, F125W, and F160W from the CANDELS survey (Grogin et al. 2011; Koekemoer et al. 2011).

## 2.3. Multiwavelength Catalogs

### 2.3.1. COSMOS-Web and PRIMER-COSMOS

For both COSMOS-Web and PRIMER-COSMOS, we conduct source detection, perform model-based photometry, and construct multiband catalogs using *SourceXtractor++* (hereafter SE++; Bertin et al. 2020; Kümmel et al. 2020), an updated version of the popular *SExtractor* package (Bertin & Arnouts 1996). The use of SE++ model-based photometry is motivated by the desire to make full use of the depth and filter coverage of seeing-limited ground-based data in COSMOS, as well as high-resolution near-infrared JWST imaging. To perform source detection, we construct a $\chi^2$ detection image (Szalay et al. 1999) from the four NIRCam bands using the *CHI_MEAN* option in *SWarp* (Bertin et al. 2002). This detection method incorporates information from all four bands to optimally identify faint sources and has proven successful in

past surveys (e.g., Weaver et al. 2022). For each detected source, SE++ then fits a Sérsic model convolved with the filter-specific point-spread function (PSF) in each measurement band, adopting priors on the source centroid based on the $\chi^2$ image. Here we use model PSFs from WebbPSF (Perrin et al. 2012, 2014). The Sérsic model parameters (centroid position, $n$, $R_{\text{eff}}$, $b/a$) are fit jointly between all bands, weighted by their respective signal-to-noise ratios (S/Ns), such that a single model is constructed to best match all available imaging data. The total flux is fit independently in each band, yielding model-based photometric measurements. The details of the SE++ catalogs for COSMOS-Web will be presented in M. Shuntov et al. (2023, in preparation).

Computing photometric uncertainties in model-based photometry is not trivial for sources that are undetected in a given band, where the SE++ model is below the noise. This can occasionally result in significantly underestimated errors in the dropout bands, where the source is not detected; as such, we set a noise floor for each band equal to the rms measured from random-position aperture photometry. In particular, we adopt the measured depths from Weaver et al. (2022) and Casey et al. (2022) using aperture diameters of 0.3" (for ACS/NIRCam), 0.6" (for MIRI), and 2" (for ground-based+Spitzer data). These apertures are generally two to three times the PSF size. We report depths in the relevant JWST filters in Table 1. The SE++ catalogs yield 1789 (585) objects detected with both NIRCam and MIRI for COSMOS-Web (PRIMER-COSMOS).

### 2.3.2. CEERS

For CEERS, we adopt the multiband *SExtractor* catalog previously described in Finkelstein et al. (2022). Source detection for this catalog is done on an inverse variance-weighted sum of the F277W and F356W images. The catalog includes all available HST/ACS and WFC3 data and all JWST/NIRCam bands but not MIRI photometry. We therefore perform independent source detection on the MIRI F770W images using *astropy/photutils*. Before performing detection, we convolve the image with a 5 × 5 pixel Gaussian smoothing kernel with an FWHM of 2 pixels in order to better identify faint objects. We set a detection threshold of 1.1 times the background rms and a minimum area of 4 pixels. In order to be consistent with the NIRCam catalog, we compute photometry on the MIRI F560W and F770W images in small elliptical apertures using the *SourceCatalog* task in *photutils*. We use a Kron factor of 1.1 to restrict the aperture to the central region of each galaxy, maximizing the S/N, and apply a correction based on the median ratio of the flux in these small apertures to equivalent apertures using a Kron factor of 2.5. Similar to the NIRCam photometry (Finkelstein et al. 2022), we find a correction of ~1.5. We then cross-match the catalog of MIRI-detected sources based on the measured centroid positions with the NIRCam catalog. This procedure yields 1190 objects detected with both NIRCam and MIRI.

## 2.4. ePSF Construction

Accurate PSF models are key to characterizing galaxy morphologies and measuring deconvolved sizes. While publicly available JWST PSF models exist (e.g., from WebbPSF; Perrin et al. 2012, 2014), and indeed were used in our catalogs, such models are known to yield somewhat narrower profiles
than are observed for stars (see, e.g., Weaver et al. 2023). This is likely due to the nature of the dithering/mosaicking process and/or scattered-light phenomena, which can broaden the idealized PSF. As such, for the purposes of morphological measurements, we construct empirical PSFs (ePSFs) for each JWST band by stacking bright stars in the field. The ePSF construction for CEERS is described in Finkelstein et al. (2022); we adopt the same PSFs used in that work. For COSMOS-Web, we construct ePSFs from a catalog of bright stars in the COSMOS field (as described in Weaver et al. 2022). We inspect cutouts at the position of each star and remove saturated stars and compact galaxies that were previously misidentified. We then construct ePSFs using the EPSFBuilder module in astropy/photutils (Bradley et al. 2022), which follows the prescriptions of Anderson & King (2000). The ePSFs are indeed slightly broader than the corresponding WebbPSF models but otherwise exhibit similar curve-of-growth shapes.

### 3. Sample Selection and Vetting

We determine initial photometric redshift estimates for the entire NIRCam+MIRI sample across all three imaging surveys using EAzY (Brammer et al. 2008). EAzY computes linear combinations of predefined templates to derive probability distribution functions (PDFs) for the redshift based on the $\chi^2$ of the templates. We fit to all available ground-based, HST, NIRCam, and MIRI photometry based on the multiwavelength catalogs described in Section 2.3. The template set we use includes the standard tweak_fsps_QSF_12_v3 set of 12 FSPS (Conroy et al. 2010) templates, as well as the six templates from Larson et al. (2022). We allow the redshift to vary from zero to 15 with a step size of $\Delta z = 0.01$. Though we refine the photometric redshifts for individual sources of interest later, this first-pass photo-$z$ run allows us to explore the relationship between observed-frame colors and redshift and outline the sample selection of massive, dusty, and red $z > 7$ galaxies in the context of the full catalog.

We identify high-redshift, dusty galaxy candidates based on their colors in $m_{444}-m_{770}$ and $m_{277}-m_{444}$. In particular, we design our selection to specifically target galaxies that are red in $m_{277}-m_{444}$ but not as red in $m_{444}-m_{770}$; this helps mitigate contamination from lower-redshift, extremely obscured sources or a dusty active galactic nucleus (AGN) whose SED continues to rise in the MIR due to hot dust in the AGN torus. Figure 1 shows the color-color diagram for MIRI/F770W-detected galaxies in COSMOS-Web (circles), CEERS (squares), and PRIMER-COSMOS (diamonds). Specifically, we require $S/N_{770} > 5$, $S/N_{444} > 5$, and $S/N_{277} > 2$ to be included in our color-color plot and eventual sample. This ensures robust detection in F770W and F444W and at least a marginal detection in F277W despite the red $m_{277}-m_{444}$ color. We additionally show in Figure 1 two BAGPIPES models (10 Myr old, 20% solar metallicity) reddened by a Calzetti dust law with $A_V = 2.5$ (solid line) and 3.5 (dashed line). We also plot a QSO model by reddening the Sloan Digital Sky Survey (SDSS) QSO composite spectrum (Vanden Berk et al. 2001; Glikman et al. 2006) with a Calzetti law with $A_V = 4$ from $1 < z < 8$.

We outline rough color-color selection criteria to encompass the reddest galaxies in our sample, which all have EAzY photometric redshifts $>7$. Specifically, we use

$$m_{277}-m_{444} > 1.8 \quad \text{and} \quad (1)$$

$$m_{277}-m_{444} > 1.5(m_{444}-m_{770}) - 0.6. \quad (2)$$

While the second criterion does not exclude any galaxies in this sample (i.e., this same selection could be done with just $m_{277}-m_{444}$), we include it to indicate that we require the $m_{444}-m_{770}$ color to be appreciably bluer than $m_{277}-m_{444}$, motivated by the need to reject objects with red MIR SEDs, which are likely obscured AGN. Indeed, the reddened QSO template shown in Figure 1 has $m_{277}-m_{444} > 1.8$ at $z > 2.5$ but is omitted from our selection criteria due to its redder color in $m_{444}-m_{770}$. These color-color criteria are optimized for selecting $z \gtrsim 6-7$ obscured galaxies with young ages and bright emission lines. Future work will refine these selection criteria with larger samples and more thorough modeling.

We find four galaxies satisfying our color-color selection based on the initial photometry: three in COSMOS-Web, one in CEERS, and none in PRIMER-COSMOS. In order to vet these candidates further, we compute photometry in custom circular apertures on each galaxy. This is intended to reject spurious detections or underestimated uncertainties that may be present in the catalog photometry used in Figure 1. In particular, we use the astropy/photutils package to perform aperture photometry in 0.5" diameter circular apertures. All aperture measurements are corrected for the PSF flux falling outside the circular aperture.

With this photometry, we rerun EAzY and inspect the resulting solutions. Two of the four candidates have highly uncertain redshift probability distributions as estimated by EAzY; both are in COSMOS-Web. We show the cutouts and SEDs for these two poorly constrained candidates in the Appendix, but for the remainder of the paper, we focus on the two remaining candidates, which are significantly brighter and, therefore, much better constrained to $z > 7$ with EAzY.

To explore any possible low-redshift redundancy to our EAzY photometric redshift solutions, we explore fits with four other SED fitting codes: PROSPECTOR (Leja et al. 2017; Johnson et al. 2021), BAGPIPES (Carnall et al. 2018), CIGALE (Burgarella et al. 2005; Noll et al. 2009; Boquien et al. 2019), and BEAGLE (Chevallard & Charlot 2016; Gutkin et al. 2016). For BAGPIPES and CIGALE, we adopt a delayed-$\tau$ star formation history (SFH) model, while for BEAGLE, we adopt a constant SFH to be consistent with recent work (Endsley et al. 2022a; Furtak et al. 2022; Whitler et al. 2023). We additionally include a late starburst in all three codes in order to allow for bright emission lines from H II regions. This is implemented as an instantaneous burst with a fixed age of 10 Myr (for BEAGLE) and a flexible age from 10 to 100 Myr (for BAGPIPES and CIGALE). We allow $A_V$ to vary from zero to 4 and adopt a Calzetti et al. (2000) attenuation law for BAGPIPES, a Charlot & Fall (2000) law for CIGALE, and the Chevallard et al. (2013) model for BEAGLE. For PROSPECTOR, we use the PROSPECTOR-$\beta$ model, which includes a nonparametric SFH and informed, joint priors encoding empirical constraints on the redshifts, stellar masses, and SFHs of observed galaxies (see Table 1 of Wang et al. 2023). Finally, we include nebular emission in all of our fits. BEAGLE uses nebular emission templates from Gutkin et al. (2016) that combine the latest Bruzual & Charlot (2003) stellar population models with CLOUDY nebular emission (Ferland et al. 2013). BAGPIPES and PROSPECTOR implement nebular emission via updated CLOUDY models from Byler et al. (2019), while CIGALE uses an updated grid of CLOUDY models as described in Boquien et al. (2019). While each of these fits use slightly different physical assumptions, the
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The $m_{277}-m_{444}$ vs. $m_{444}-m_{770}$ color–color diagram indicating selection of $z \gtrsim 7$ dusty galaxies. Points indicate objects detected at 5σ significance in F444W and F770W and at least 2σ significance in F277W. The points are colored by their best-fit photometric redshift from first-pass EAZY runs, and contours indicate the point density. The dashed lines indicate our proposed color selection criterion, which captures 4 objects, all with photometric redshifts $\gtrsim 7$. We show in gray two model tracks generated from BAGPIPES for young (10 Myr) stellar populations with $A_V = 2.5$ and 3.5 from $z \sim 1$ to 9. We additionally show in red a QSO model with $A_V = 4$, which, at $z > 3$, is redder in $m_{444}-m_{770}$ than the galaxy models. The color–color selection shown here is optimized for selecting $z \gtrsim 6-7$ obscured galaxies with young ages and bright emission lines but rejecting objects with steeply rising MIR SEDs (e.g., obscured AGN).</figcaption></figure>

breadth of approaches here serves as a valuable test of the security of the candidates as genuine $z > 7$ galaxies.

# 4. Results

## 4.1. Two Remarkable Dust-obscured, $z \sim 8$ Candidates

Based on the sample selection and candidate vetting described in Section 3, we identify two robust $z \sim 8$ dusty galaxies with very red colors ($m_{277}-m_{444} \sim 2.5$), which we denote as COS-z8M1 and CEERS-z7M1 for their selection via MIRI detections. Figures 2 and 3 show the best-fit SEDs and redshift probability distributions for the two candidates. We show 1.5 arcsec$^2$ cutouts in all available HST/ACS, JWST/NIRCam, and JWST/MIRI bands and, additionally, a stack of ground-based imaging for the COSMOS source. We plot the SEDs from the BAGPIPES, BEAGLE, CIGALE, and PROSPECTOR fits; for clarity, we show the $P(z)$ from EAZY but not the best-fit SED. We note that the EAZY $\chi^2$ is similar to the other codes, and forcing a $z < 7$ solution with EAZY yields a significantly higher $\chi^2$ than the best-fit redshift (i.e., $\Delta\chi^2 > 10$).

Object COS-z8M1 (Figure 2) is only detected in F277W, F444W, and F770W and marginally in Spitzer/IRAC [3.6]. The extremely red $m_{277}-m_{444}$ color despite a relatively blue $m_{444}-m_{770}$ color drives the redshift solutions to $z > 7$. This SED shape could be due to the redshifted rest-frame 4000 Å break (e.g., Labbé et al. 2023a), which would suggest a maximally old stellar population that formed most of its mass by $z \sim 15$. Alternatively, the SED shape could be due to the contribution of [O III] 5007 Å emission to the F444W flux. Bright emission lines have been shown to play a significant role in elevating broadband photometry in high-redshift galaxies (e.g., Shim et al. 2011; Stark et al. 2013;

Faisst et al. 2016; Fujimoto et al. 2022b; Naidu et al. 2022b; McKinney et al. 2022; Arrabal Haro et al. 2023; Zavala et al. 2023), and longer-wavelength MIRI data are critical to constrain the underlying continuum (Papovich et al. 2022). Indeed, the emission line-dominated scenario is preferred by our SED fitting codes, as shown in Figure 2. At $7 < z < 9$, [O III] 5007 Å emission falls into F444W, while $H\alpha$ falls blueward of F770W, yielding an extreme color differential. A high [O III]+$H\beta$ equivalent width (EW) of $\sim 800$ Å is needed to contribute sufficiently to the F444W flux. This is consistent with recent spectroscopic results for $z > 7$ galaxies, for which the [O III] 5007 Å emission line is known to be particularly bright (Fujimoto et al. 2023a; Saxena et al. 2023; Trump et al. 2023).

Similarly, for CEERS-z7M1, the redshift solution is significantly constrained by the elevated flux in NIRCam/F410M+F444W and MIRI/F560W, likely due to bright [O III] and $H\alpha$ emission, respectively. The source is well detected in all NIRCam long wavelength and MIRI bands (S/N $> 10$), and the contribution from bright emission lines constrains the redshift to $z \sim 7.6$. Importantly, the redshift probability distribution is consistent, albeit broader, if we fit only the COSMOS-Web filter set at the appropriate depth (e.g., only F277W, F444W, and F770W). This is promising for the fidelity of finding similar objects in large surveys like COSMOS-Web.

We do note that CEERS-z7M1 is marginally detected in NIRCam F115W, F150W, and F200W, which may suggest the presence of a less-obscured, UV-luminous component to the galaxy. This is not inconsistent with the SED for COS-z8M1, for which the short-wavelength data are too shallow to constrain any UV emission. Since there is no signal in ACS/F606W or F814W, we consider the $z > 7$ solution robust; indeed, fitting EAZY to only $\lambda < 2 \mu$m data gives $z \sim 6-8$. While none of the SED fitting codes capture this blue component, this is to be expected, as these codes assume a uniform dust screen attenuating all the starlight. In reality, some UV emission may escape unattenuated in the case of patchy dust.

Table 2 gives the physical properties derived from SED fitting for COS-z8M1 and CEERS-z7M1 for each SED fitting code used. The different codes generally agree on a high stellar mass for this epoch ($\log M_*/M_\odot > 9.5$), a steeply rising SFH (with SFR$_{10 Myr}$/SFR$_{100 Myr} \sim 10$), and a high dust attenuation ($A_V \sim 1.5-2.5$ mag) for both candidates. This steeply rising SFH (i.e., young stellar age) is what drives the strong emission lines in these fits. Modeling this effect is only possible due to the inclusion of a recent starburst in the SED fitting parameters, though we note that the nonparametric SFH from PROSPECTOR independently yields a late starburst. For the sake of simplicity, we adopt the photometric redshifts and physical parameters from PROSPECTOR for the remainder of this paper, though we add an additional 0.2 dex uncertainty to the stellar mass to capture the differences between different fits. However, given that the different fits are consistent, despite differences in their assumptions for stellar population synthesis/dust attenuation and different sampling algorithms, these results are not particularly sensitive to the modeling assumptions.

## 4.2. Rest-frame Optical Sizes

We use the 2D image fitting code IMFIT$^{47}$ (Erwin 2015) to characterize the rest-frame optical sizes of COS-z8M1 and CEERS-z7M1. In particular, we fit PSF-convolved models to

<footnote>$^{47}$ https://www.mpe.mpg.de/~erwin/code/imfit/</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Optical through MIR SED of COS-z8M1. The top panels show 1.5 arcsec² cutouts in the available HST/JWST bands plus stacked ground-based imaging from Subaru/HSC and UltraVISTA. The bottom panel shows the measured photometry (or 2σ upper limits) from all available JWST bands (shown in dark red), as well as HST/ACS, Spitzer/IRAC, and ground-based bands (light red). We additionally show maximum a posteriori model SEDs from CIGALE, BEAGLE, BAGPIPES, and PROSPECTOR (for clarity, we do not show the EAzY SED). The inset shows the full redshift probability distributions P(z) from each photo-z code. Robust detection in NIRCam/F444W and MIRI/F770W with a sharp break between F444W and F277W mandates a massive, dusty, $z > 7$ solution.</figcaption></figure>

the NIRCam/F444W images. While F444W has the largest PSF of all of the NIRCam bands, it is also the highest S/N detection for both sources. We use the F444W ePSF used in our photometry as described in Section 2.3.

We find that both sources are well fit by a point-source model. Figure 4 shows the results of the point-source fitting for the two candidates (data, model, and residual). While fitting a Sérsic model yields a marginal improvement in the reduced $\chi^2$ statistic, the resulting Sérsic parameters are poorly constrained. In order to provide a constraint on $R_{\text{eff}}$, we fix the Sérsic index $n=1$, the axis ratio $b/a=1$, and the position angle $\theta=0$ (i.e., an exponential disk profile). We fit a series of models with $R_{\text{eff}} \sim 0''01-0''07$ (50–350 pc at $z=8$). We find that the resulting residuals are significant at the $3\sigma$ level for $R_{\text{eff}} \gtrsim 200$ pc but consistent with the background rms for $R_{\text{eff}} \lesssim 200$ pc, so we adopt this as an upper limit on the true size.

We note that such compactness is relatively common in high-redshift galaxies observed with JWST (Ono et al. 2022; Roberts-Borsani et al. 2022; Robertson et al. 2022; Tacchella et al. 2023). Nonetheless, given that neither source is resolved, we discuss the possibility that the source's emission is dominated by AGN in Section 5.2. However, we note that the derived stellar masses are likely robust to any contamination from strong emission lines from an AGN.

### 4.3. Stellar Masses

We compare the stellar masses derived from SED fitting for these two sources to estimates for similar objects in the literature. Figure 5 shows the stellar mass versus redshift for COS-z8M1, CEERS-z7M1, and numerous dust-obscured, spectroscopically confirmed galaxies in the literature (Marrone et al. 2018; Tamura et al. 2019; Endsley et al. 2022b; Ferrara et al. 2022b). We additionally show the sample of NIRCam-selected, unobscured, spectroscopically confirmed CEERS objects from Fujimoto et al. (2023a). Both COS-z8M1 and CEERS-z7M1 represent the extreme end of the dust-obscured high-$z$ population, with stellar masses $\gtrsim 3$ times higher than other known dust continuum-detected objects at the same redshift. We note that the various SED fitting codes used in this work all yield consistently large stellar masses, even when accounting for extremely high EW emission lines and differing dust attenuation laws (as discussed in Section 4.1).

The detection of these two sources at $z \sim 7-9$ across 20 arcmin² of combined NIRCam+MIRI imaging suggests an approximate volume density of $\sim 2 \times 10^{-5} \text{ Mpc}^{-3}$. We estimate the typical stellar mass for sources of this rarity based on the evolving halo mass function. We compute the halo mass function using the python package `hmf` (Murray et al. 2013). We adopt a Tinker et al. (2008) parameterization, modified with the redshift-dependent parameters from Rodríguez-Puebla et al. (2016) to be
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Same as Figure 2 but for CEERS-z7M1. The NIRCam short wavelength bands reveal a blue UV slope ($\beta \sim -2.5$) at 1–2 µm despite the very red continuum at $>2$ µm. This could be from star formation in an unobscured line of sight or an unobscured AGN component. The additional filters constrain the redshift to precisely $z = 7.6$.</figcaption></figure>

<table><caption>**Table 2**

Physical Properties Derived from SED Fitting for the Two Dust-obscured, $z \sim 8$ Sources</caption><thead><tr><th>ID<br>(1)</th><th>R.A., decl.<br>(J2000)<br>(2)</th><th>Photo-z<br>Code<br>(3)</th><th>Photometric<br>Redshift<br>(4)</th><th>log M<sub>*</sub><br>(M<sub>☉</sub>)<br>(5)</th><th>SFR<sub>10 Myr</sub><br>(M<sub>☉</sub> yr<sup>-1</sup>)<br>(6)</th><th>SFR<sub>100 Myr</sub><br>(M<sub>☉</sub> yr<sup>-1</sup>)<br>(7)</th><th>A<sub>V</sub><br>(mag)<br>(8)</th></tr></thead><tbody><tr><td rowspan="6">COS-z8M1</td><td rowspan="6">09<sup>h</sup>59<sup>m</sup>40<sup>s</sup> 18, +02<sup>d</sup>17<sup>m</sup>32<sup>s</sup> 19</td><td>PROSPECTOR-β</td><td>8.4<sup>+0.3</sup><sub>-0.4</sub></td><td>9.8<sup>+0.1</sup><sub>-0.1</sub></td><td>520<sup>+110</sup><sub>-100</sub></td><td>59<sup>+22</sup><sub>-14</sub></td><td>1.6<sup>+0.3</sup><sub>-0.2</sub></td></tr><tr><td>BAGPIPES</td><td>8.4<sup>+0.4</sup><sub>-0.5</sub></td><td>9.9<sup>+0.2</sup><sub>-0.2</sub></td><td>320<sup>+140</sup><sub>-90</sub></td><td>73<sup>+51</sup><sub>-30</sub></td><td>2.4<sup>+0.2</sup><sub>-0.2</sub></td></tr><tr><td>BEAGLE</td><td>8.2<sup>+0.4</sup><sub>-0.5</sub></td><td>9.5<sup>+0.2</sup><sub>-0.2</sub></td><td>270<sup>+120</sup><sub>-90</sub></td><td>27<sup>+13</sup><sub>-8</sub></td><td>1.9<sup>+0.3</sup><sub>-0.3</sub></td></tr><tr><td>CIGALE</td><td>8.7<sup>+0.2</sup><sub>-0.4</sub></td><td>10.3<sup>+0.4</sup><sub>-0.3</sub></td><td>480<sup>+290</sup><sub>-460</sub></td><td>57<sup>+35</sup><sub>-24</sub></td><td>1.1<sup>+0.3</sup><sub>-0.1</sub></td></tr><tr><td>EAZY</td><td>8.4<sup>+0.3</sup><sub>-0.3</sub></td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CEERS-z7M1</td><td>14<sup>h</sup>19<sup>m</sup>43<sup>s</sup> 08, +52<sup>d</sup>53<sup>m</sup>16<sup>s</sup> 48</td><td>PROSPECTOR-β</td><td>7.63<sup>+0.06</sup><sub>-0.08</sub></td><td>10.3<sup>+0.1</sup><sub>-0.1</sub></td><td>1590<sup>+320</sup><sub>-420</sub></td><td>179<sup>+40</sup><sub>-46</sub></td><td>3.0<sup>+0.2</sup><sub>-0.2</sub></td></tr><tr><td></td><td></td><td>BAGPIPES</td><td>7.62<sup>+0.08</sup><sub>-0.06</sub></td><td>10.2<sup>+0.2</sup><sub>-0.2</sub></td><td>600<sup>+120</sup><sub>-120</sub></td><td>151<sup>+103</sup><sub>-60</sub></td><td>3.5<sup>+0.1</sup><sub>-0.2</sub></td></tr><tr><td></td><td></td><td>BEAGLE</td><td>7.68<sup>+0.04</sup><sub>-0.03</sub></td><td>10.0<sup>+0.1</sup><sub>-0.1</sub></td><td>940<sup>+110</sup><sub>-120</sub></td><td>103<sup>+13</sup><sub>-14</sub></td><td>3.5<sup>+0.1</sup><sub>-0.2</sub></td></tr><tr><td></td><td></td><td>CIGALE</td><td>8.2<sup>+0.1</sup><sub>-0.6</sub></td><td>10.3<sup>+0.2</sup><sub>-0.1</sub></td><td>1160<sup>+270</sup><sub>-210</sub></td><td>122<sup>+28</sup><sub>-21</sub></td><td>1.6<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td></td><td></td><td>EAZY</td><td>8.1<sup>+0.3</sup><sub>-0.5</sub></td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

**Note.** Stellar masses and SFRs from BEAGLE and CIGALE are corrected by a factor of 1.12 to be consistent with our assumption of a Kroupa (2002) IMF.

consistent with a Planck cosmology (see also discussion in Yung et al. 2023). We compute the halo mass $M_{\text{halo}}$ associated with a volume density $n(>M_{\text{halo}}) \approx 2 \times 10^{-5}$ Mpc$^{-3}$. This halo mass is then converted to a stellar mass assuming a cosmic baryon fraction of 15.8% and a baryon conversion efficiency, $\epsilon$. The orange stripe in Figure 5 shows the expected stellar masses for $\epsilon \sim 20\%–30\%$; we do not include uncertainties (e.g., from cosmic variance) in the width of the stripe but rather just use it to illustrate the expected range of masses for different halo growth histories from $6.5 < z < 9$.

### 4.4. FIR Constraints and Dust Mass

Given that both COS-z8M1 and CEERS-z7M1 are very red, and our SED fits imply significant dust obscuration, they should be luminous in the FIR. Both the COSMOS and EGS fields have coverage with Spitzer/MIPS, Herschel/PACS and SPIRE, and JCMT/SCUBA-2, and we investigate these maps for any emission at the position of the NIRCam source. Specifically, we measure flux densities at the expected position of the source and find no significant emission in any of these FIR/submillimeter observations. This is consistent with
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Results of 2D point-source profile fitting to the NIRCam/F444W imaging for both sources. The columns show the data, best-fit model, and residuals for each source. Both sources are well characterized by a point-source model, suggesting extremely compact sizes. We derive an upper limit on the true sizes of $R_{\text{eff}} \lesssim 200$ pc, as discussed in the text.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Inferred stellar mass vs. redshift for COS-z8M1 (red) and CEERS-z7M1 (green), adopting the results from PROSPECTOR. We additionally show measurements for spectroscopically confirmed, dust continuum-detected star-forming galaxies from the literature (Marrone et al. 2018; Tamura et al. 2019; Endsley et al. 2022b; Ferrara et al. 2022b), as well as recent UV-bright galaxies confirmed with JWST/NIRSpec (Fujimoto et al. 2023a). The orange stripe shows the range of stellar masses expected for halos with a constant number density, computed from an evolving halo mass function assuming a baryon conversion efficiency $\epsilon \sim 20\%–30\%$. Both COS-z8M1 and CEERS-z7M1 appear to be among the most massive dust-obscured galaxies at this epoch.</figcaption></figure>

expectations for a ULIRG at $z>7$ given the depth of these data.

However, COS-z8M1 is covered in the Ex-MORA survey (A. Long et al. 2023, in preparation), a blind 2 mm Atacama Large Millimeter/submillimeter Array (ALMA) survey of the COSMOS field designed to identify high-redshift submillimeter galaxies (an extension of the original MORA survey presented in Casey et al. 2021; Zavala et al. 2021).48 While Ex-MORA is still relatively shallow ($5\sigma \sim 1$ mJy), it is deeper than the available SCUBA-2 data and has a smaller beam size ($\theta_{\text{beam}} \sim 1''5$). At 2 mm, such a data set is optimal for

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Contours of the marginal 2 mm emission from Ex-MORA overlaid on the JWST/F444W image. The location of the 2.9$\sigma$ peak is marked with a cross, offset ~0.5" from the NIRCam source but consistent within the positional uncertainty given the large beam size. No other $\pm 2\sigma$ signal falls within 2" of the source. The Ex-MORA flux, if real, suggests a high IR luminosity of $\log L_{\text{IR}}/L_{\odot} \sim 11.9_{-0.3}^{+0.4}$, comparable to local ULIRGs.</figcaption></figure>

identifying DSFGs at higher redshifts. We find a marginal detection near the position of COS-z8M1, with a peak S/N of 2.9$\sigma$ and flux density $S_{2\,\text{mm}} = 0.19 \pm 0.07$ mJy offset by ~0.5" from the NIRCam source. Figure 6 shows the contours of this marginal detection overlaid on the NIRCam/F444W image. While FIR positional offsets of ~0.5"–1.0" are not uncommon in high-redshift DSFGs (e.g., Chapman et al. 2004; Biggs & Ivison 2008; Hodge et al. 2012; Inami et al. 2022), these are typically observed between the rest-frame UV and the FIR (rather than the rest-frame optical and FIR) and in more extended sources. However, we note that the 2 mm centroid position is uncertain to ~0.5" (Condon 1997; Ivison et al. 2007a) due to the low S/N and ~1.5" beam size, making it feasibly associated with the NIRCam source; as such, we do not consider the offset significant.

Though it requires follow-up verification, we examine the implications of this marginal detection for the FIR luminosity $L_{\text{IR}}$. The only constraining power in the FIR SED comes from the marginal 2 mm Ex-MORA detection and the 850 $\mu$m upper limit from SCUBA-2 ($3\sigma = 2.9$ mJy). We fit the FIR data to piecewise functions with an MIR power law and an FIR modified blackbody (as in Casey 2012; Drew & Casey 2022). We use a custom Markov Chain Monte Carlo (MCMC) routine (based on MCIRSED; Drew & Casey 2022) with flat priors on $L_{\text{IR}}$, $T_{\text{dust}}$, and $\beta$. In the absence of significant FIR constraints beyond our one 2 mm data point, we fix $\alpha_{\text{MIR}}$ (the slope of the MIR power law) to 2.3 and $\lambda_0$ (the wavelength at which $\tau = 1$) to 200 $\mu$m (following the recommendations in Drew & Casey 2022). We allow $T_{\text{dust}}$ to vary from 26 K ($\approx T_{\text{CMB}}$ at $z = 8.5$) to 90 K and $\beta$ to vary from 1.5 to 2.4. Due to the very negative $k$-correction in the submillimeter (Casey et al. 2014a), the model SED is insensitive to the precise redshift in the 2 mm regime, so we adopt a fixed redshift of $z = 8.5$. We account for the effects of heating by and decreasing contrast against the

<footnote n="48">48 CEERS-z7M1 is not covered by any archival NOEMA pointing in the EGS field.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Illustration of the various galaxy+AGN possibilities associated with the composite SED shape of CEERS-z7M1, with a blue UV slope but red optical colors. Left: galaxy-only composite model with a massive ($\log M_{\star}/M_{\odot} = 10$) obscured ($A_{\nu} = 4$) galaxy combined with a low-mass ($\log M_{\star}/M_{\odot} = 7$) unobscured galaxy. Physically, this corresponds to a patchy distribution of dust with “holes” allowing some UV light to escape unattenuated. Middle: galaxy+quasar model with a massive obscured galaxy combined with a faint unobscured type I quasar. We adopt the SDSS QSO composite spectrum (Vanden Berk et al. 2001; Glikman et al. 2006). Right: galaxy+quasar model with a low-mass unobscured galaxy combined with a bright but heavily dust-reddened type I quasar.</figcaption></figure>

cosmic microwave background (CMB) following da Cunha et al. (2013). We note that this MCMC fitting is not intended to constrain $T_{\text{dust}}$ or $\beta$ but rather to estimate $L_{\text{IR}}$ marginalized over the uncertainty in the other SED parameters.

This fit yields an IR luminosity of $\log L_{\text{IR}}/L_{\odot} \sim 11.9_{-0.3}^{+0.4}$.
Based on the Murphy et al. (2011) calibration, this corresponds to an obscured SFR of $\sim 110_{-60}^{+160} M_{\odot} \text{ yr}^{-1}$, which is broadly consistent with the SFR$_{100 Myr}$ values in Table 2 (but lower than SFR$_{10 Myr}$, which primarily traces the emission line strengths). This implies a large fraction of obscured star formation ($\sim 99\%$), as the upper limits in the rest-frame UV (uncorrected for dust) imply $\text{SFR}_{\text{UV}} \lesssim 1 M_{\odot} \text{ yr}^{-1}$.

We additionally apply the methodology outlined in Scoville et al. (2016; see also Casey et al. 2019, Section 3.3) to provide an order-of-magnitude estimate for the dust mass in COS-z8M1.<fnref n="49" /> This method depends on the monochromatic flux at some wavelength (here 2 mm), the emissivity index, and the mass-weighted dust temperature, which is not the same as the luminosity-weighted temperature that can be derived in SED fitting. The mass-weighted temperature is typically colder than the SED temperature, as the total mass is dominated by cold dust (see discussion in Liang et al. 2019; Lower et al. 2023). While previous studies at $z < 6$ have adopted a mass-weighted temperature of 25 K (e.g., Scoville et al. 2016), this is colder than the CMB temperature at $z = 8$; the actual dust temperature must be hotter than the CMB, and indeed, simulations predict higher mass-weighted temperatures with increasing redshift (Liang et al. 2019). For temperatures $T \sim 30–40$ K and $\beta \sim 1.5–2.4$, the observed 2 mm flux corresponds to a dust mass of $\log M_{\text{dust}}/M_{\odot} = 8.5_{-0.5}^{+0.5}$.

The dust mass and IR luminosity measurements are consistent with the inferred UV attenuation. In particular, assuming an average attenuation curve with $A_{\text{UV}} = 2.6 A_{\nu}$, we can estimate the infrared excess as $\text{IRX} = 1.68(10^{0.4A_{\text{UV}}} - 1)$ (e.g., Casey et al. 2014b). This implies $\log L_{\text{IR}}/L_{\odot} \sim 12.4$ and a dust mass of $\log M_{\text{dust}}/M_{\odot} \sim 8.7$ assuming $T = 35$ K and $\beta = 1.95$, consistent with the measurement derived from the ALMA 2 mm flux. This dust mass is more than an order of magnitude larger than the only known $z > 8$ galaxy with a dust continuum detection (Tamura et al. 2019). Deeper millimeter observations will be needed to verify this dust mass

measurement and place constraints on the obscured SFR in this object.

## 5. Discussion

Imaging with JWST is already leading to the identification of many interesting high-z candidates, some of which have unexpectedly red rest-frame optical colors. Labbé et al. (2023a) presented a sample of red, $z_{\text{phot}} > 8$ candidates from early CEERS imaging that exhibited “double breaks,” i.e., both the Lyman and Balmer break. They interpreted the red $2–4 \ \mu\text{m}$ color as tracing the 4000 Å break and derived very large masses of $M_{\star} \gtrsim 10^{10} M_{\odot}$, potentially in excess of limits from $\Lambda$CDM cosmology (Boylan-Kolchin 2022; Menci et al. 2022). Recent work has suggested that the stellar masses of high-z candidates may be overestimated due to differences in the IMF (e.g., Steinhardt et al. 2022) or contamination from strong emission lines from AGN (as discussed in Endsley et al. 2022a; Labbé et al. 2023a; Brummel-Smith et al. 2023; Yuan et al. 2023). Indeed, JWST spectroscopy has already revealed a number of reddened AGN at $z > 5$ based on their broadened Balmer emission lines (e.g., Harikane et al. 2023b; Kocevski et al. 2023; Matthee et al. 2023), suggesting an abundant population of these objects in the early Universe. Similar to this work, Furtak et al. (2022) identified an extremely red, compact object at $z \sim 7.7$ in deep JWST/NIRCam imaging as part of the UNCOVER survey (Bezanson et al. 2022). Aided by lensing magnification, they constrained the effective radius to $\lesssim 35$ pc, providing strong evidence for a low-luminosity quasar with strong emission lines driving the red colors.

These results raise the question of whether these red, compact, high-z objects are massive dust-obscured galaxies, reddened quasars, or both. In the absence of spectroscopy, many have explored various SED models to characterize these unique sources. In particular, many of these red, compact objects show a unique SED feature: a red continuum at $\lambda_{\text{rest}} > 3000$ Å but a blue UV slope at shorter wavelengths. This SED shape is similar to CEERS-z7M1 and consistent with COS-z8M1 given the shallower depth of the NIRCam data. It has been interpreted as a composite galaxy+AGN signature, with either the blue or red component comprising the host galaxy or AGN (e.g., Barro et al. 2023; Kocevski et al. 2023), or as a signature of a reddened quasar, with the blue component arising from scattered light from the AGN broad-line region (e.g., Labbé et al. 2023b). Regardless of the origin, this unique

<footnote n="49">⁴⁹ We note that this method assumes that the dust is optically thin, which may not be the case at 2 mm (rest-frame ~200 $\mu$m); however, in this case, the dust mass would be underestimated.</footnote>
SED shape appears ubiquitous across a range of magnitudes and redshifts and has even been observed in faint $z \sim 9-10$ candidates from the NGDEEP survey (Leung et al. 2023). We show in Figure 7 an illustration of the various galaxy/AGN possibilities consistent with this SED shape, in particular, a galaxy–galaxy composite model (left), a blue quasar + red galaxy model (middle), and a blue galaxy + red quasar model (right). In the following subsections, we discuss the likelihood and implications of the two scenarios (dust-obscured galaxy versus quasar) for the candidates presented in this work.

## 5.1. The High-z Extreme of the DSFG Population?

The $z>7$ regime is notoriously difficult for DSFG identification (Casey et al. 2018a, 2018b); therefore, FIR spectroscopic follow-up has been limited to the brightest DSFGs (with $L_{IR} \sim 10^{13} L_{\odot}$; Marrone et al. 2018; Fujimoto et al. 2022a; Endsley et al. 2022b). At the same time, FIR follow-up of UV-luminous galaxies, which are far more numerous and easily characterized than DSFGs, has found significant evidence for obscured star formation at $z \sim 6-8$ (Schouws et al. 2022; Algera et al. 2023; Barrufet et al. 2023). These two samples represent two complementary approaches, from the rest-UV and the rest-FIR, to constraining the dust content of the early Universe. However, constraints on the population in between—obscured in the rest-UV but not so bright as to be readily detectable in the IR—have been limited. The selection method presented in this paper, targeting extremely red objects in JWST/NIRCam+MIRI imaging, represents an alternative approach.

The two candidates presented in this paper appear to represent the high-redshift extreme of the population of moderately obscured DSFGs at $z \sim 2-6$ (Chapman et al. 2005; Aretxaga et al. 2011; Smolčić et al. 2012; McKinney et al. 2023). This would place these sources among the progenitor population of the rare, extreme star-forming factories at $z \sim 5-7$ (e.g., Marrone et al. 2018; Zavala et al. 2018; Casey et al. 2019; Endsley et al. 2022b). Indeed, the estimated volume density for these sources, $n \sim 2 \times 10^{-5} \text{ Mpc}^{-3}$, is comparable to the (albeit poorly constrained) volume density of $z \sim 7$ luminous infrared galaxies ($L_{IR} > 10^{11} L_{\odot}$), for which integrating the IR luminosity function yields $n \sim 0.5-10 \times 10^{-5} \text{ Mpc}^{-3}$ (Zavala et al. 2021; Fujimoto et al. 2023b).

While neither candidate presented in this work is resolved in NIRCam/F444W imaging, such compact galaxies is consistent with expectations for early galaxies. Indeed, to build up extreme stellar masses by $z \sim 8$ requires efficient funneling of gas into cold, dense clouds, and compact starbursts are often observed in bright submillimeter galaxies (Condon et al. 1991; Ma et al. 2016; Jin et al. 2022). Theoretical work has suggested that massive galaxies at ultrahigh redshift may be able to form efficiently via feedback-free starbursts (Dekel et al. 2023), which predict compact, massive objects by $z \sim 8-10$. Indeed, many $z > 9$ galaxies being confirmed by JWST are ultracompact with effective radii of $\lesssim 200$ pc (Ono et al. 2022; Roberts-Borsani et al. 2022; Robertson et al. 2022; Tacchella et al. 2023). The stellar mass densities for the objects in this work ($\sim 10^{11} M_{\odot} \text{ kpc}^{-2}$) are similar to the values measured for local and high-redshift elliptical galaxies and globular clusters, which span a large dynamic range in stellar mass (Hopkins et al. 2010). These objects may therefore represent the

progenitors of massive ellipticals at $z \sim 0$, forming a dense core early on and then growing inside-out from there (see, e.g., Baggen et al. 2023). We note, however, that the F444W photometry is likely dominated by the [O III] emission at $z \sim 7.5-8.5$, which may not be the best tracer of the morphology of the stellar continuum.

The composite blue+red SED of CEERS-z7M1 could be due to a patchy distribution of dust in an overall very dust-obscured galaxy. In fact, many observed submillimeter galaxies at lower redshift ($z \lesssim 5$; where we could efficiently probe the rest-UV continuum pre-JWST) show blue UV slopes despite significant infrared excess (e.g., Casey et al. 2014b). Theoretical work has shown that a patchy dust geometry could allow a faint blue component to shine through “holes” in the ISM dust screen despite most of the stellar light being obscured (Popping et al. 2017; Narayanan et al. 2018). Indeed, we show in the left panel of Figure 7 an illustration of the SED in this “patchy dust” scenario, in which a low-mass unobscured stellar population ($\log M_{\star}/M_{\odot} \sim 7$, $A_{V} = 0$) emits alongside a dust-obscured massive galaxy ($\log M_{\star}/M_{\odot} \sim 10$, $A_{V} \sim 4$), with an areal fraction of holes in the dust screen of $\sim 0.1\%$. This patchy geometry can be captured in part by adopting a very “gray” attenuation curve (e.g., Witt & Gordon 1996; Charlot & Fall 2000). The fact that more and more $z > 5$ galaxies are being identified with this unique SED shape is perhaps to be expected, as the sensitivity of JWST allows us to observe, for the first time, the faint rest-UV emission from high-z dust-obscured galaxies, which are expected to have complex star-dust geometry (Ma et al. 2019).

### 5.1.1. Implications for Dust Production at $z > 8$

The interpretation of these candidates as high-z dust-obscured galaxies would imply an early buildup of dust. In particular, the marginal 2 mm flux from the Ex-MORA survey, if real and associated with COS-z8M1, suggests a high IR luminosity ($\log L_{IR}/L_{\odot} \sim 12$) and dust mass ($\log M_{\text{dust}}/M_{\odot} \sim 8.5$). These estimates are highly uncertain and require direct submillimeter follow-up for confirmation. However, if confirmed, they would have strong implications for the buildup of dust in the early Universe. The implied dust-to-stellar mass ratio of $M_{\text{dust}}/M_{\star} \approx 0.03$ is significantly higher than that of low-redshift DSFGs (e.g., Dunne et al. 2011) but broadly consistent with the observed increase of this ratio at higher redshift (Calura et al. 2014, 2017).

Furthermore, measurements of the dust mass at high redshift can provide stringent constraints on the relative contribution of different dust production mechanisms. There is abundant evidence that AGB stars, which are the dominant producers of dust later in cosmic time, are not able to produce large dust masses in $\lesssim 1$ Gyr (Valiante et al. 2009; Dwek & Cherchneff 2011; Asano et al. 2013); instead, dust production in SNe has been invoked to explain the large dust masses in $z > 7$ galaxies. The recent detection of signatures of carbonaceous dust composition at $z = 6.7$ (Witstok et al. 2023) implies a large dust mass that requires either that significant star formation occurred at $z > 10$ or, more likely, that faster dust production channels dominate.

Even in the absence of FIR constraints, we can derive maximal dust masses based on the dust yield per AGB star or SN. In the dust-obscured galaxy interpretation, the stellar mass estimates for both sources in this paper are constrained thanks to MIRI imaging, making possible estimates of the number of AGB stars/SNe.
Following the method outlined in Michałowski (2015), we estimate $N_{AGB}$ and $N_{SN}$ by integrating the IMF from 3 to 8 and 8 to 40 $M_{\odot}$, respectively. Here we assume a Kroupa (2002) IMF to be consistent with our PROSPECTOR-derived masses. We assume a theoretical maximum dust yield of $1.3 M_{\odot}$ per SN$^{50}$ and $0.04 M_{\odot}$ per AGB star (Michałowski 2015). Based on the derived stellar masses of COS-z8M1 and CEERS-z7M1 (Table 2), we derive maximum dust masses of $\sim 10^8 M_{\odot}$ from SN production and $\sim 10^7 M_{\odot}$ from AGB stars.

Other dust production mechanisms may therefore be needed to explain the high dust masses in these early galaxies, if confirmed. Asano et al. (2013) found that above a certain metallicity threshold ($\sim 0.3 Z_{\odot}$), ISM grain growth can dominate over stellar dust production and form $\sim 10^7 M_{\odot}$ of dust in $\sim 100$ Myr. Even more exotic, the maximal dust yields may be higher in unique cases such as dust produced in supershells (Martínez-González et al. 2021), the wake around Wolf-Rayet stars (Lau et al. 2021, 2022), winds around an AGN accretion disk (Sarangi et al. 2019), or red supergiant winds of high-mass Population III stars (Nozawa et al. 2014). Moreover, a top-heavy IMF (which has been suggested to be common in dust-obscured starbursts and high-$z$ star formation in general; e.g., Zhang et al. 2018) could increase the maximum dust mass by increasing the number of high-mass stars overall, especially for Population III stars. A higher early SN rate may also drive accelerated grain growth in the ISM by the contribution of additional seed metals at $z > 10$. Taken together, these results highlight the importance of deep submillimeter follow-up of these objects and future samples of similarly selected objects to constrain the dust masses and physical processes responsible for the buildup of dust in the early Universe.

### 5.1.2. Implications for Stellar Mass Assembly at $z \sim 8$

The apparent ubiquity of massive ($M_{*} > 10^{10} M_{\odot}$) galaxies identified at $z > 8$ with JWST has produced tensions with ΛCDM (Boylan-Kolchin 2022; Menci et al. 2022). In the dust-obscured galaxy interpretation, the observed photometry implies stellar masses of $\sim 10^{10} M_{\odot}$ for COS-z8M1 and CEERS-z7M1, even after correcting for the contribution from extremely bright emission lines. Compared to the implied halo rarity, this suggests a high efficiency of converting baryons into stars, $\epsilon \sim 25\%$. This is consistent with results of Inayoshi et al. (2020), who derived constraints on the star formation efficiency based on the $z > 10$ candidates identified in early JWST imaging (Naidu et al. 2022a; Finkelstein et al. 2022; Harikane et al. 2022, 2023a; Donnan et al. 2023). Inayoshi et al. noted that, alternatively, such high stellar masses could be explained by a low baryon conversion efficiency in a metal-free stellar population with a top-heavy IMF. While this explanation may prove true for the UV-luminous population, the implied dust obscuration in the candidates presented here suggests a relatively metal-rich stellar population several generations beyond Population III. Therefore, if the stellar masses of these sources prove robust to constraints on any significant AGN contribution, they may suggest highly efficient stellar mass buildup in the early Universe.

To explore the implications of these large stellar masses, we derive the cumulative stellar mass density $\rho_{*}(>M_{*})$ versus $M_{*}$ based on the two candidates presented in this work. We estimate

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Cumulative stellar mass density $\rho_{*}$ vs. $M_{*}$. The red points indicate the cumulative stellar mass density inferred from the detection of COS-z8M1 and CEERS-z7M1. We additionally plot the result for $z \sim 8$ from Labbé et al. (2023a) and integrated Schechter function fits to UV-luminous LBGs from Stefanon et al. (2021). In tabulating these results, we adopt the same redshift bin, $z \sim 7-9$ (see main text), and remove the $z = 5.6$ AGN from the Labbé et al. (2023a) sample. The light blue points show the Labbé et al. (2023a) sample with the most massive candidate removed. The candidates presented in this work are consistent with the Stefanon et al. (2021) $z \sim 7-9$ stellar mass function, though the uncertainties are large given the identification of just two sources.</figcaption></figure>

the effective volume as the differential comoving volume integrated over a redshift bin from $z \sim 7$ to 9 and scaled to the total survey area of 20.6 arcmin$^2$ (Table 1). Error bars include the cosmic variance uncertainty, which we compute from halo number counts in the DREaM simulation (Drakos et al. 2022), as well as the Poisson counting uncertainty. We note that the Poisson uncertainty dominates over cosmic variance ($\sim 70\%–90\%$ versus $\sim 35\%$) given the detection of just two sources.

The resulting estimates are shown in Figure 8 alongside the equivalent measurements from Labbé et al. (2023a); we compute $\rho_{*}$ in the same manner using their full $z \sim 7–9$ sample after removing the $z = 5.6$ AGN identified in Kocevski et al. (2023). The high stellar mass density from Labbé et al. (2023a) is driven by the one $\log M_{*}/M_{\odot} \sim 10.9$ candidate at $z \sim 7.5$; we show in light blue the $\rho_{*}$ estimates with this candidate removed, just to highlight its impact on $\rho_{*}$. The line and shaded region show the stellar mass density corresponding to the Schechter fit stellar mass function from Stefanon et al. (2021), derived from samples of UV-luminous, IRAC-detected galaxies. To produce a curve comparable to our derived $\rho_{*}$, we integrate the $\sim 7, 8, \text{ and } 9$ stellar mass functions and take a weighted mean. We weight each curve by the volume associated with redshift bins from $z \sim 7$ to 7.5, 7.5 to 8.5, and 8.5 to 9, respectively.

We find that the stellar mass density inferred from this work is formally consistent with the $z \sim 7–9$ stellar mass function inferred from the UV-luminous, IRAC-detected population, with a possible slight excess driven by the mass of CEERS-z7M1. The large uncertainties prohibit a robust determination of the stellar mass density, and the relatively small area probed in this work prohibits constraints across a large dynamic range. We note that the Labbé et al. (2023a) sample is also fully consistent with the Stefanon et al. (2021) stellar mass function

<footnote>$^{50}$ This assumes no dust destruction, likely an unphysical assumption. A more realistic yield of $\sim 0.1–0.15 M_{\odot} \text{ SN}^{-1}$ (Leśniewska & Michałowski 2019) would give a lower dust mass by a factor of $\sim 10$.</footnote>
if the $\log M_*/M_\odot \sim 10.9$ candidate were in fact a lower stellar mass (as suggested by Endsley et al. 2022a). Larger samples (e.g., from the remaining COSMOS-Web and PRIMER imaging) are needed to constrain the contribution of obscured galaxies to the cosmic stellar mass density at this epoch.

## 5.2. An Early Population of Obscured AGN/Reddened Quasars?

Given the compact morphology and likely contribution of bright nebular lines, the observed emission in these candidates could be dominated by nuclear activity. This could impact not just the EWs of optical emission lines but also, potentially, continuum emission. Based on size-mass and size-z scaling relationships for known $z \sim 7-10$ star-forming galaxies (Ono et al. 2013; Holwerda et al. 2015), we might expect $R_{\text{eff}} \sim 800$ pc at $M_* \sim 10^{10} M_\odot$, significantly larger than the upper limit inferred from the F444W imaging ($\sim 200$ pc). However, this is by no means conclusive evidence for an AGN, as the dispersion in the size-mass relation is large at this epoch ($\sim 200-300$ pc), and (as discussed in Section 5.1) compact star formation appears common at $z \sim 8$ (e.g., Ono et al. 2022).

The existing multiwavelength data could provide some constraints on the AGN contribution in these objects. Neither source presented in this work is detected in existing VLA radio data at 1.4, 3, and 5 GHz (Willner et al. 2006; Ivison et al. 2007b; Schinnerer et al. 2007; Smolčić et al. 2017) or X-ray imaging from Chandra (Laird et al. 2009; Nandra et al. 2015; Civano et al. 2016; Marchesi et al. 2016). While the radio nondetections are not particularly constraining, as only $\sim 10\%$ of the AGN are expected to be radio-loud (Liu et al. 2021), the X-ray nondetections do provide upper limits on the AGN bolometric luminosity and black hole mass. For COS-z8M1, adopting the 2–10 keV hard band flux limit of the COSMOS-Legacy Survey (Civano et al. 2016),<fnref n="51" /> we find $L_{2-10\,\text{keV}} < 3 \times 10^{44} \text{ erg s}^{-1}$ (assuming a photometric redshift of $z = 8.5$). Similarly, for CEERS-z7M1, adopting the flux limit of the AEGIS-X survey (Nandra et al. 2015), we find $L_{2-10\,\text{keV}} < 4.6 \times 10^{44} \text{ erg s}^{-1}$ (for $z = 7.5$). Applying the X-ray bolometric correction from Shen et al. (2020), we find $L_{\text{bol}} < 6 \times 10^{45} \text{ erg s}^{-1}$ for COS-z8M1 and $< 9.2 \times 10^{45} \text{ erg s}^{-1}$ for CEERS-z7M1. Assuming $\lambda_{\text{Edd}} = 0.1$, this would suggest black hole masses of $M_{\text{BH}} < 5 \times 10^8$ and $< 7 \times 10^8 M_\odot$ for COS-z8M1 and CEERS-z7M1, respectively. These upper limits are on the order of the masses of bright $z \sim 6-7$ quasars (Shen et al. 2019). The X-ray limits, therefore, do not rule out low-mass AGN similar to the recent spectroscopically confirmed red AGN from recent JWST/NIRSpec data (e.g., Harikane et al. 2023b; Kocevski et al. 2023).

While the present data therefore cannot yield a clear answer on the AGN contribution in these sources, here we discuss the different scenarios in which an AGN would impact the observed SED. First, a heavily obscured (i.e., type II) AGN would emit strongly in the rest-frame MIR due to hot torus dust. However, this is likely not significant shortward of rest-frame $\sim 2 \, \mu$m, which is unconstrained by the present depth and filter coverage. Indeed,

the PROSPECTOR-$\beta$ model includes emission from the dusty torus but yields highly unconstrained values of $f_{\text{AGN}}$ and $\tau_{\text{AGN}}$. At this epoch, constraining the rest-frame MIR SED will be difficult and require ultradeep MIR imaging in the reddest wavelengths, which are also the least sensitive.

A second scenario involves strong emission lines from an AGN contaminating the broadband photometry (e.g., Endsley et al. 2022a). While we already account for the potential for strong [O III] and H$\alpha$ emission (albeit from star-forming H II regions with slightly different emission line ratios) in our fits, strong emission lines from an AGN could also contribute significantly. To check this, we run PROSPECTOR, including an empirical, scalable template for emission lines from the AGN narrow-line region (NLR; based on data from Richardson et al. 2014). We first force the redshift to $z < 7$ to examine the likelihood of lower-z AGN interlopers (i.e., Kocevski et al. 2023); the resulting fits favor $z \sim 5$ but require unphysically high emission line EWs to match the observed photometry ([O III]+H$\beta$ EW $\sim 22,000 \AA$) and still achieve poor fits (e.g., $\chi^2 \sim 1.8$ for COS-z8M1 and 5.5 for CEERS-z7M1).

We then adopt a redshift prior based on the PROSPECTOR photometric redshifts reported in Table 2. We find consistent photometric redshifts and slightly ($\sim 0.2 \, \text{dex}$) lower stellar masses than our galaxy-only fits. This difference in stellar mass comes largely from the fact that the current PROSPECTOR AGN NLR models do not include the associated nebular continuum emission, whereas the galaxy-only SED fits do. Regardless, as the MIRI 7.7 $\, \mu$m flux is unaffected by the inclusion of strong AGN emission lines in the model, we consider this a robust constraint on the underlying continuum emission. Stellar mass estimates are therefore robust to this effect; however, an emission line contribution from AGN would impact the current SFRs derived from SED fitting. These AGN emission lines could be coming from an obscured type II AGN or an unobscured type I AGN. In the latter case (i.e., the picture presented in the middle panel of Figure 7), these sources may be in a transition stage between dustenshrouded starbursts and unobscured luminous quasars (e.g., Fu et al. 2017; Fujimoto et al. 2022a).

Finally, given the compact nature of these sources and the relative expected rarity of $M_* \sim 10^{10} M_\odot$ systems, one might suspect that their reddened continuum emission is in fact not stellar in origin but rather dominated by a highly reddened quasar. This interpretation is shown in the right panel of Figure 7 and discussed further in Barro et al. (2023) and Labbé et al. (2023b). In this case, the red continuum would be dominated by thermal emission from a dust-reddened accretion disk, and any residual UV emission could be from the host galaxy (as shown in Figure 7) or due to scattered light from the AGN (as in Labbé et al. 2023b). Unfortunately, constraining the relative contribution from stellar versus quasar continuum will require rest-frame MIR diagnostics or deeper X-ray data, beyond what is feasible with current facilities.

Identification of such luminous reddened quasars in extremely low-mass galaxies would be unexpected; for example, the rightmost SED in Figure 7 would have an implied black hole mass of $M_{\text{BH}} \gtrsim 10^8 M_\odot$ (assuming $\lambda_{\text{Edd}} = 0.1$) with a comparable host galaxy stellar mass of $M_* \sim 10^7 M_\odot$, the ratio of which is well outside of expectation (McConnell & Ma 2013), even at high-z (e.g., Izumi et al. 2021). The volume density of reddened type I quasars is highly unconstrained. On the one hand, UV-luminous quasars are known to be very rare, $\sim 1000$ times rarer than these sources, as inferred by integrating the Matsuoka et al. (2018) $z = 6$

<footnote n="51">⁵¹ We adopt the flux limits reported in Civano et al. (2016) and Nandra et al. (2015) for COSMOS-Legacy and AEGIS-X, respectively. This assumes a Galactic $N_H = 2.6$ and $1.3 \times 10^{20} \, \text{cm}^{-2}$ for the two surveys, respectively, and a power-law index $\Gamma = 1.4$ for both. Given the varying flux depth across the fields, due to the varying PSF size and shape of Chandra, we adopt the 90% completeness flux limit from Nandra et al. (2015). For COSMOS, the overlapping tiling survey strategy of COSMOS-Legacy produces a relatively uniform limiting flux depth; we adopt the 50% completeness limit.</footnote>
quasar luminosity function down to $M_{UV} \sim -18$. At the same time, recent JWST/NIRspec results have revealed an abundant population of broad-line AGN in $z \sim 7$ UV-faint galaxies (Harikane et al. 2023b). However, under the red QSO interpretation, the objects presented in this work are distinct from the populations of UV-luminous quasars (which are by definition unobscured by dust) or broad-line AGN (not all of which would be expected to outshine the stellar continuum). The volume density of $z \gtrsim 7$ reddened quasars is likely somewhere in between the rare, UV-bright QSOs and the abundant broad-line AGN.

In summary, while these objects may host AGN, their measured stellar masses are robust to contributions from strong emission lines and hot dust torus emission. The major caveat is that we cannot rule out the possibility of continuum emission from dust-reddened type I quasars. However, we conclude that the dust-obscured galaxy interpretation is more likely based on the expected number densities of these classes of objects. In particular, we note that the typical star formation depletion time for $z \gtrsim 4$ DSFGs ($\sim 100-300$ Myr; Swinbank et al. 2014; Aravena et al. 2016; Williams et al. 2019; Manning et al. 2022) is significantly longer than the typical quasar lifetime of $\sim 1-10$ Myr (Marconi et al. 2004; Volonteri et al. 2015; Eilers et al. 2020). Follow-up spectroscopy will nevertheless be needed to search for AGN signatures (i.e., broadened Balmer lines or weak high-ionization lines such as N V or C III]).

## 6. Summary

In this paper, we present a search for extremely red, dust-obscured, $z > 7$ galaxies in three publicly available Cycle 1 surveys. By focusing on sources detected in JWST/NIRCam +MIRI imaging, we construct a unique selection for massive red galaxies at $z > 7$.

1. We identify two candidates, COS-z8M1 and CEERS-z7M1, which have extremely red colors ($m_{277}-m_{444} \sim 2.5$) and robust photometric redshifts of $8.5_{-0.4}^{+0.3}$ and $7.6_{-0.1}^{+0.1}$. The photometry for both sources is likely impacted by strong emission lines, particularly [O III] +Hβ in F444W and Hα+[N II] in F560W. Both candidates are significantly more dust-obscured ($A_V \sim 2-3$) than other known $z \sim 8$ galaxies.

2. We find that neither source is resolved in NIRCam/F444W, constraining the rest-frame size to $R_{\text{eff}} \lesssim 200$ pc.

3. We infer stellar masses of $\sim 10^{10} M_{\odot}$, significantly higher than the known dust continuum-detected galaxies at $z > 8$ and similar to some of the most massive $z > 8$ galaxy candidates yet identified by JWST. The inferred stellar mass density is consistent within the uncertainty with expectations from the UV-luminous population.

4. We identify a marginal, $2.9\sigma$ detection at 2 mm near the position of COS-z8M1 as part of the Ex-MORA survey. We show that this flux, if real, suggests an IR luminosity of $\sim 10^{12} L_{\odot}$, consistent with the constraints on attenuation suggested by the JWST data. There are no submillimeter constraints for CEERS-z7M1.

This work highlights the importance of long-wavelength MIR imaging for characterization of massive, dust-obscured galaxies at $z > 7$. Given the remarkable sensitivity of MIRI, almost 1 mag deeper than preflight expectations (Casey et al. 2022), it becomes possible to constrain the $z > 7$ galaxy population at rest-frame $\sim 1$ µm. Future efforts to explore the dust-obscured population at this epoch will benefit greatly from deep MIRI imaging in

multiple filters. While the completion of the full COSMOS-Web and PRIMER surveys will likely result in the detection of dozens more of these objects, spectroscopy and submillimeter follow-up will be necessary to determine the nature of the dust-obscured population in the epoch of reionization.

## Acknowledgments

Support for this work was provided by NASA through grants JWST-GO-01727 and HST-AR-15802 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. H.B.A. acknowledges the support of the UT Austin Astronomy Department and thanks the UT Austin College of Natural Sciences for support through the Harrington Graduate Fellowship. C.M.C. thanks the National Science Foundation for support through grants AST-1814034 and AST-2009577, as well as the University of Texas at Austin College of Natural Sciences for support; C.M.C. also 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.

Finally, H.B.A., C.M.C., and others at UT Austin acknowledge that they work at an institution that sits on indigenous land. The Tonkawa lived in central Texas, and the Comanche and Apache moved through this area. We pay our respects to all the American Indian and Indigenous Peoples and communities who have been or have become a part of these lands and territories in Texas.

*Facilities:* HST (ACS and WFC3), JWST (NIRCam and MIRI), Spitzer (IRAC and MIPS), Herschel (PACS and SPIRE), SCUBA-2, VLA, ALMA.

All of the JWST data used in this work can be found in MAST at doi:10.17909/4eps-dc89.

*Software:* EAzY (Brammer et al. 2008), BAGPIPES (Carnall et al. 2018), PROSPECTOR (Leja et al. 2017; Johnson et al. 2021), CIGALE (Boquien et al. 2019), IMFIT (Erwin 2015), SExtractor (Bertin & Arnouts 1996), SourceXtractorPlusPlus (Bertin et al. 2020; Kümmel et al. 2020), astropy (Astropy Collaboration 2013), matplotlib (Hunter 2007), numpy (Harris et al. 2020), photutils (Bradley et al. 2022), scipy (Virtanen et al. 2020), STScI JWST Calibration Pipeline (jwst-pipeline.readthedocs.io; Rigby et al. 2022).

## Appendix

### Additional Candidate SEDs and Photometry

Table A1 provides the photometric fluxes in the HST/ACS, JWST/NIRCam, and JWST/MIRI bands for all objects in this paper. Figure A1 shows the SEDs for COS-939 (left) and COS-3627 (right), the two candidates that fall into our color-color selection but have less robust SEDs due to redder $m_{444}-m_{770}$ colors (and generally lower fluxes). We show the fiducial EAzY fit in black and a forced $z < 6$ model in gray.

The first candidate (COS-939) has a particularly red $m_{444}-m_{770}$ color ($\sim 1.5$) and drops out blueward of F277W, yielding a redshift PDF that peaks at $z \sim 9$ but with a significant probability at $z \sim 5$ ($\sim 7\%$). We note, however, that COS-939 was initially identified as a potential close companion to COS-z8M1, located only $\sim 6''$ away ($\sim 30$ kpc at $z = 8$). This potential association lends some credence to the possibility of a $z > 7$ solution, but spectroscopic observations will be needed for confirmation. The second candidate (COS-3627) suffers from a “snowball” artifact in the F150W image; this artificially
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Figure A1.** Optical-MIR SEDs of COS-939 and COS-3627, the two candidates that fall into our color selection criteria but are fainter and thus have poorly constrained redshifts. As in Figures 2 and 3, we show cutouts in the ACS, NIRCam, and MIRI bands and the measured photometry in red. We show the best-fit EAzY SED in black and a forced $z < 6$ model in gray.</figcaption></figure>

<table><caption>**Table A1**
Photometry from HST+JWST</caption><thead><tr><th rowspan="2">ID</th><th colspan="2">HST/ACS</th><th colspan="7">JWST/NIRCam</th><th colspan="2">JWST/MIRI</th></tr><tr><th>F606W</th><th>F814W</th><th>F115W</th><th>F150W</th><th>F200W</th><th>F277W</th><th>F356W</th><th>F410M</th><th>F444W</th><th>F560W</th><th>F770W</th></tr></thead><tbody><tr><td>COS-z8M1</td><td>...</td><td>(8 ± 9)</td><td>(22 ± 25)</td><td>(6 ± 19)</td><td>...</td><td>41 ± 10</td><td>150 ± 33<sup>a</sup></td><td>...</td><td>307 ± 11</td><td>...</td><td>375 ± 73</td></tr><tr><td>CEERS-z7M1</td><td>(-1 ± 5)</td><td>(4 ± 5)</td><td>11 ± 3</td><td>(7 ± 4)</td><td>(8 ± 3)</td><td>17 ± 2</td><td>69 ± 2</td><td>158 ± 6</td><td>188 ± 5</td><td>451 ± 30</td><td>338 ± 30</td></tr><tr><td>COS-939</td><td>...</td><td>(-6 ± 9)</td><td>(-28 ± 16)</td><td>(-2 ± 13)</td><td>...</td><td>22 ± 7</td><td>(-91 ± 33)<sup>a</sup></td><td>...</td><td>118 ± 8</td><td>...</td><td>301 ± 73</td></tr><tr><td>COS-3627</td><td>...</td><td>(-3 ± 9)</td><td>(16 ± 18)</td><td>(0 ± 16)</td><td>...</td><td>(14 ± 9)</td><td>(15 ± 33)<sup>a</sup></td><td>...</td><td>68 ± 11</td><td>...</td><td>162 ± 73</td></tr></tbody></table>

Note. All fluxes are in nJy.

a Reported flux is Spitzer/IRAC [3.6].

elevated flux yields a blue $m_{150}-m_{277}$ color and drives the ultrahigh-$z$ solution as shown in Figure 1. When we mask this artifact and recompute the photometry (as shown in Figure A1), we find an incredibly broad redshift PDF with a nearly equal likelihood anywhere from $z \sim 4$ to 12; the source is too faint to achieve a robust constraint.

## ORCID iDs

Hollis B. Akins https://orcid.org/0000-0003-3596-8794

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

Natalie Allen https://orcid.org/0000-0001-9610-7950

Micaela B. Bagley https://orcid.org/0000-0002-9921-9218

Mark Dickinson https://orcid.org/0000-0001-5414-5131

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

Maximilien Franco https://orcid.org/0000-0002-3560-8599

Santosh Harish https://orcid.org/0000-0003-0129-2079

Pablo Arrabal Haro https://orcid.org/0000-0002-7959-8783

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

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

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

Daizhong Liu https://orcid.org/0000-0001-9773-7479

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

Henry Joy McCracken https://orcid.org/0000-0002-9489-7765

Louise Paquereau https://orcid.org/0000-0003-2397-0360

Casey Papovich https://orcid.org/0000-0001-7503-8482

Nor Pirzkal https://orcid.org/0000-0003-3382-5941

Jason Rhodes https://orcid.org/0000-0002-4485-8549

Brant E. Robertson https://orcid.org/0000-0002-4271-0364

Marko Shuntov https://orcid.org/0000-0002-7087-0701

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

Guang Yang https://orcid.org/0000-0001-8835-7722

Guillermo Barro https://orcid.org/0000-0001-6813-875X

Laura Bisigello https://orcid.org/0000-0003-0492-4924

Véronique Buat https://orcid.org/0000-0003-3441-903X

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

Olivia Cooper https://orcid.org/0000-0003-3881-1397

Luca Costantin https://orcid.org/0000-0001-6820-0015

Alexander de la Vega https://orcid.org/0000-0002-6219-5558

Nicole E. Drakos https://orcid.org/0000-0003-4761-2197

Andreas Faisst https://orcid.org/0000-0002-9382-9832

Adriano Fontana https://orcid.org/0000-0003-3820-2823

Steven Gillman https://orcid.org/0000-0001-9885-4589

Carlos Gómez-Guijarro https://orcid.org/0000-0002-4085-9165

Ghassem Gozaliasl https://orcid.org/0000-0002-0236-919X

Nimish P. Hathi https://orcid.org/0000-0001-6145-5090

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

Michaela Hirschmann https://orcid.org/0000-0002-3301-3321
Benne W. Holwerda https://orcid.org/0000-0002-4884-6756

Shuowen Jin https://orcid.org/0000-0002-8412-7951

Dale D. Kocevski https://orcid.org/0000-0002-8360-3880

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

Erini Lambrides https://orcid.org/0000-0003-3216-7190

Ray A. Lucas https://orcid.org/0000-0003-1581-7825

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

Benjamin Magnelli https://orcid.org/0000-0002-6777-6490

Jed McKinney https://orcid.org/0000-0002-6149-8178

Bahram Mobasher https://orcid.org/0000-0001-5846-4404

Pablo G. Pérez-González https://orcid.org/0000-0003-4528-5639

R. Michael Rich https://orcid.org/0000-0003-0427-8387

Lise-Marie Seillé https://orcid.org/0000-0001-7755-4755

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

C. Megan Urry https://orcid.org/0000-0002-0745-9792

Francesco Valentino https://orcid.org/0000-0001-6477-4011

Katherine E. Whitaker https://orcid.org/0000-0001-7160-3632

L. Y. Aaron Yung https://orcid.org/0000-0003-3466-035X

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

## References

Aihara, H., AlSayyad, Y., Ando, M., et al. 2019, PASJ, 71, 114

Algera, H. S. B., Inami, H., Oesch, P. A., et al. 2023, MNRAS, 518, 6142

Anderson, J., & King, I. R. 2000, PASP, 112, 1360

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

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

Arrabal Haro, P., Dickinson, M., Finkelstein, S. L., et al. 2023, arXiv:2303.15431

Asano, R. S., Takeuchi, T. T., Hirashita, H., & Inoue, A. K. 2013, EP&S, 65, 213

Astropy Collaboration 2013, A&A, 558, A33

Baggen, J. F. W., van Dokkum, P., Labbe, I., et al. 2023, ApJL, 955, L12

Bagley, M. B., Finkelstein, S. L., Koekemoer, A. M., et al. 2023, ApJL, 946, L12

Barro, G., Perez-Gonzalez, P. G., Kocevski, D. D., et al. 2023, arXiv:2305.14418

Barrufet, L., Oesch, P. A., Bouwens, R., et al. 2023, MNRAS, 522, 3926

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

Bertin, E., Mellier, Y., Radovich, M., et al. 2002, in ASP Conf. Ser. 281, Astronomical Data Analysis Software and Systems XI, ed. D. A. Bohlender, D. Durand, & T. H. Handley (San Francisco, CA: ASP), 228

Bertin, E., Schefer, M., Apostolakos, N., et al. 2020, in ASP Conf. Ser. 527, Astronomical Data Analysis Software and Systems XXIX, ed. E. Bertin, M. Schefer, & N. Apostolakos (San Francisco, CA: ASP), 461

Bezanson, R., Labbe, I., Whitaker, K. E., et al. 2022, arXiv:2212.04026

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

Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103

Bouwens, R., Stefanon, M., Brammer, G., et al. 2023, MNRAS, 523, 1036

Boylan-Kolchin, M. 2023, NatAs, 7, 731

Bradley, L., Sipőcz, B., Robitaille, T., et al. 2022, astropy/photutils: 1.6.0, Zenodo, doi:10.5281/zenodo.7419741

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

Brummel-Smith, C., Skinner, D., Sethuram, S. S., et al. 2023, MNRAS, 525, 4405

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

Burgarella, D., Buat, V., & Iglesias-Páramo, J. 2005, MNRAS, 360, 1413

Byler, N., Dalcanton, J. J., Conroy, C., et al. 2019, AJ, 158, 2

Calura, F., Gilli, R., Vignali, C., et al. 2014, MNRAS, 438, 2765

Calura, F., Pozzi, F., Cresci, G., et al. 2017, MNRAS, 465, 54

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Carnall, A. C., McLure, R. J., Dunlop, J. S., & Davé, R. 2018, MNRAS, 480, 4379

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

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

Casey, C. M., Kartaltepe, J. S., Drakos, N. E., et al. 2023, ApJ, 954, 31

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

Casey, C. M., Scoville, N. Z., Sanders, D. B., et al. 2014b, ApJ, 796, 95

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

Chapman, S. C., Smail, I., Windhorst, R., Muxlow, T., & Ivison, R. J. 2004, ApJ, 611, 732

Charlot, S., & Fall, S. M. 2000, ApJ, 539, 718

Chevallard, J., & Charlot, S. 2016, MNRAS, 462, 1415

Chevallard, J., Charlot, S., Wandelt, B., & Wild, V. 2013, MNRAS, 432, 2061

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

Condon, J. J. 1997, PASP, 109, 166

Condon, J. J., Huang, Z. P., Yin, Q. F., & Thuan, T. X. 1991, ApJ, 378, 65

Conroy, C., White, M., & Gunn, J. E. 2010, ApJ, 708, 58

Cooper, O. R., Casey, C. M., Zavala, J. A., et al. 2022, ApJ, 930, 32

Curtis-Lake, E., Carniani, S., Cameron, A., et al. 2023, NatAs, 7, 622

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

Dekel, A., Sarkar, K. S., Birnboim, Y., Mandelker, N., & Li, Z. 2023, MNRAS, 523, 3201

Donnan, C. T., McLeod, D. J., Dunlop, J. S., et al. 2023, MNRAS, 518, 6011

Drakos, N. E., Villasenor, B., Robertson, B. E., et al. 2022, ApJ, 926, 194

Drew, P. M., & Casey, C. M. 2022, ApJ, 930, 142

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

Dwek, E., & Cherchneff, I. 2011, ApJ, 727, 63

Eilers, A.-C., Hennawi, J. F., Decarli, R., et al. 2020, ApJ, 900, 37

Endsley, R., Stark, D. P., Whitler, L., et al. 2023a, MNRAS, 524, 2312

Endsley, R., Stark, D. P., Lyu, J., et al. 2023b, MNRAS, 520, 4609

Erwin, P. 2015, ApJ, 799, 226

Euclid Collaboration, Moneti, A., McCracken, H. J., et al. 2022, A&A, 658, A126

Faisst, A. L., Capak, P., Hsieh, B. C., et al. 2016, ApJ, 821, 122

Ferland, G. J., Porter, R. L., van Hoof, P. A. M., et al. 2013, RMxAA, 49, 137

Ferrara, A., Pallottini, A., & Dayal, P. 2023, MNRAS, 522, 3986

Ferrara, A., Sommovigo, L., Dayal, P., et al. 2022b, MNRAS, 512, 58

Finkelstein, S., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJL, 946, L13

Fu, H., Isbell, J., Casey, C. M., et al. 2017, ApJ, 844, 123

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

Fujimoto, S., Arrabal Haro, P., Dickinson, M., et al. 2023b, ApJL, 949, L25

Fujimoto, S., Brammer, G. B., Watson, D., et al. 2022a, Natur, 604, 261

Fujimoto, S., Kohno, K., Ouchi, M., et al. 2023a, arXiv:2303.01658

Fujimoto, S., Finkelstein, S. L., Burgarella, D., et al. 2022b, arXiv:2211.03896

Fujimoto, S., Oguri, M., Brammer, G., et al. 2021, ApJ, 911, 99

Furtak, L. J., Zitrin, A., Plat, A., et al. 2023, ApJ, 952, 142

Gaia Collaboration 2018, A&A, 616, A1

Glikman, E., Helfand, D. J., & White, R. L. 2006, ApJ, 640, 579

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

Gutkin, J., Charlot, S., & Bruzual, G. 2016, MNRAS, 462, 1757

Harikane, Y., Inoue, A. K., Mawatari, K., et al. 2022, ApJ, 929, 1

Harikane, Y., Ouchi, M., Oguri, M., et al. 2023a, ApJS, 265, 5

Harikane, Y., Zhang, Y., Nakajima, K., et al. 2023b, arXiv:2303.11946

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Hodge, J. A., Carilli, C. L., Walter, F., et al. 2012, ApJ, 760, 11

Holwerda, B. W., Bouwens, R., Oesch, P., et al. 2015, ApJ, 808, 6

Hopkins, P. F., Murray, N., Quataert, E., & Thompson, T. A. 2010, MNRAS Lett., 401, L19

Hunter, J. D. 2007, CSE, 9, 90

Inami, H., Algera, H. S. B., Schouws, S., et al. 2022, MNRAS, 515, 3126

Inayoshi, K., Visbal, E., & Haiman, Z. 2020, ARA&A, 58, 27

Ivison, R. J., Chapman, S. C., Faber, S. M., et al. 2007b, ApJL, 660, L77

Ivison, R. J., Greve, T. R., Dunlop, J. S., et al. 2007a, MNRAS, 380, 199

Izumi, T., Matsuoka, Y., Fujimoto, S., et al. 2021, ApJ, 914, 36

Jin, S., Daddi, E., Magdis, G. E., et al. 2022, A&A, 665, A3

Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22

Jones, A. P., Fanciullo, L., Köhler, M., et al. 2013, A&A, 558, A62

Kocevski, D. D., Onoue, M., Inayoshi, K., et al. 2023, ApJL, 954, L4

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

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

Kroupa, P. 2002, Sci, 295, 82

Kümmel, M., Bertin, E., Schefer, M., et al. 2020, in ASP Conf. Ser. 527, Astronomical Data Analysis Software and Systems XXIX, ed. R. Pizzo et al. (San Francisco, CA: ASP), 29

Labbé, I., Greene, J. E., Bezanson, R., et al. 2023b, arXiv:2306.07320

Labbé, I., van Dokkum, P., Nelson, E., et al. 2023a, Natur, 266, 266

Laird, E. S., Nandra, K., Georgakakis, A., et al. 2009, ApJS, 180, 102

Laporte, N., Ellis, R. S., Boone, F., et al. 2017, ApJL, 837, L21
Larson, R. L., Finkelstein, S. L., Kocevski, D. D., et al. 2023, ApJL, 953, L29

Larson, R. L., Hutchison, T. A., Bagley, M., et al. 2022, arXiv:2211.10035

Lau, R. M., Hankins, M. J., Han, Y., et al. 2022, NatAs, 6, 1308

Lau, R. M., Hankins, M. J., Kasliwal, M. M., et al. 2021, ApJ, 909, 113

Leja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, ApJ, 837, 170

Leśniewska, A., & Michałowski, M. J. 2019, A&A, 624, L13

Leung, G. C. K., Bagley, M. B., Finkelstein, S. L., et al. 2023, ApJL, 954, L46

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

Liu, Y., Wang, R., Momjian, E., et al. 2021, ApJ, 908, 124

Long, A. S., Casey, C. M., Lagos, C. d. P., et al. 2022, ApJ, 953, 11

Lovell, C. C., Harrison, I., Harikane, Y., Tacchella, S., & Wilkins, S. M. 2023, MNRAS, 518, 2511

Lower, S., Narayanan, D., Hu, C.-Y., & Privon, G. C. 2023, arXiv:2306.07338

Ma, J., Gonzalez, A. H., Vieira, J. D., et al. 2016, ApJ, 832, 114

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

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

Manning, S. M., Casey, C. M., Zavala, J. A., et al. 2022, ApJ, 925, 23

Marchesi, S., Civano, F., Elvis, M., et al. 2016, ApJ, 817, 34

Marconi, A., Risaliti, G., Gilli, R., et al. 2004, MNRAS, 351, 169

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

Martínez-González, S., Silich, S., & Tenorio-Tagle, G. 2021, MNRAS, 507, 1175

Matsuoka, Y., Strauss, M. A., Kashikawa, N., et al. 2018, ApJ, 869, 150

Matthee, J., Naidu, R. P., Brammer, G., et al. 2023, arXiv:2306.05448

McConnell, N. J., & Ma, C.-P. 2013, ApJ, 764, 184

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

McKinney, J., Finnerty, L., Casey, C., et al. 2023, ApJL, 946, L39

McKinney, J., Manning, S. M., Cooper, O. R., et al. 2023, arXiv:2304.07316

Menci, N., Castellano, M., Santini, P., et al. 2022, ApJL, 938, L5

Michałowski, M. J. 2015, A&A, 577, A80

Morishita, T., Roberts-Borsani, G., Treu, T., et al. 2023, ApJL, 947, L24

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

Murray, S. G., Power, C., & Robotham, A. S. G. 2013, A&C, 3, 23

Naidu, R. P., Oesch, P. A., van Dokkum, P., et al. 2022a, ApJL, 940, L14

Naidu, R. P., Oesch, P. A., Setton, D. J., et al. 2022b, arXiv:2208.02794

Nandra, K., Laird, E. S., Aird, J. A., et al. 2015, ApJS, 220, 10

Narayanan, D., Davé, R., Johnson, B. D., et al. 2018, MNRAS, 474, 1718

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

Noll, S., Burgarella, D., Giovannoli, E., et al. 2009, A&A, 507, 1793

Nozawa, T., Yoon, S.-C., Maeda, K., et al. 2014, ApJL, 787, L17

Oke, J. B. 1974, ApJS, 27, 21

Ono, Y., Harikane, Y., Ouchi, M., et al. 2023, ApJ, 951, 72

Ono, Y., Ouchi, M., Curtis-Lake, E., et al. 2013, ApJ, 777, 155

Papovich, C., Cole, J., Yang, G., et al. 2023, ApJL, 949, L18

Pérez-González, P. G., Barro, G., Annunziatella, M., et al. 2023, ApJL, 946, L16

Perrin, M. D., Sivaramakrishnan, A., Lajoie, C.-P., et al. 2014, Proc. SPIE, 9143, 91433X

Perrin, M. D., Soummer, R., Elliott, E. M., Lallo, M. D., & Sivaramakrishnan, A. 2012, Proc. SPIE, 8442, 84423D

Planck Collaboration 2020, A&A, 641, A6

Popping, G., Puglisi, A., & Norman, C. A. 2017, MNRAS, 472, 2315

Richardson, C. T., Allen, J. T., Baldwin, J. A., Hewett, P. C., & Ferland, G. J. 2014, MNRAS, 437, 2376

Rigby, J., Perrin, M., McElwain, M., et al. 2023, PASP, 135, 048001

Roberts-Borsani, G., Treu, T., Chen, W., et al. 2023, Natur, 618, 480

Robertson, B. E., Tacchella, S., Johnson, B. D., et al. 2023, NatAs, 7, 611

Rodighiero, G., Bisigello, L., Iani, E., et al. 2023, MNRAS Lett., 518, L19

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

Sarangi, A., Dwek, E., & Kazanas, D. 2019, ApJ, 885, 126

Saxena, A., Robertson, B. E., Bunker, A. J., et al. 2023, arXiv:2302.12805

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

Schouws, S., Stefanon, M., Bouwens, R., et al. 2022, ApJ, 928, 31

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

Shen, X., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2020, MNRAS, 495, 3252

Shen, Y., Wu, J., Jiang, L., et al. 2019, ApJ, 873, 35

Shim, H., Chary, R.-R., Dickinson, M., et al. 2011, ApJ, 738, 69

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

Stark, D. P., Schenker, M. A., Ellis, R., et al. 2013, ApJ, 763, 129

Stefanon, M., Bouwens, R. J., Labbé, I., et al. 2021, ApJ, 922, 29

Steidel, C. C., Giavalisco, M., Pettini, M., Dickinson, M., & Adelberger, K. L. 1996, ApJL, 462, L17

Steinhardt, C. L., Kokorev, V., Rusakov, V., Garcia, E., & Sneppen, A. 2023, ApJL, 951, L40

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

Szalay, A. S., Connolly, A. J., & Szokoly, G. P. 1999, AJ, 117, 68

Tacchella, S., Eisenstein, D. J., Hainline, K., et al. 2023, ApJ, 952, 74

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

Tamura, Y., Mawatari, K., Hashimoto, T., et al. 2019, ApJ, 874, 27

Tang, M., Stark, D. P., Chen, Z., et al. 2023, MNRAS, in press

Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709

Trump, J. R., Haro, P. A., Simons, R. C., et al. 2023, ApJ, 945, 35

Valiante, R., Schneider, R., Bianchi, S., & Andersen, A. C. 2009, MNRAS, 397, 1661

Vanden Berk, D. E., Richards, G. T., Bauer, A., et al. 2001, AJ, 122, 549

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Volonteri, M., Silk, J., & Dubus, G. 2015, ApJ, 804, 148

Wang, B., Leja, J., Bezanson, R., et al. 2023, ApJL, 944, L58

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

Weaver, J. R., Cutler, S. E., Pan, R., et al. 2023, arXiv:2301.02671

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

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

Whitler, L., Endsley, R., Stark, D. P., et al. 2023, MNRAS, 519, 157

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

Willner, S. P., Coil, A. L., Goss, W. M., et al. 2006, AJ, 132, 2159

Witsok, J., Shivaei, I., Smit, R., et al. 2023, arXiv:2302.05468

Witt, A. N., & Gordon, K. D. 1996, ApJ, 463, 681

Yang, G., Papovich, C., Bagley, M., et al. 2023, arXiv:2307.14509

Yuan, G.-W., Lei, L., Wang, Y.-Z., et al. 2023, arXiv:2303.09391

Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., Wilkins, S. M., & Gardner, J. P. 2023, arXiv:2304.04348

Zavala, J. A., Buat, V., Casey, C. M., et al. 2023, ApJL, 943, L9

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

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

Zhang, Z.-Y., Romano, D., Ivison, R. J., Papadopoulos, P. P., & Matteucci, F. 2018, Natur, 558, 260