# Two massive, compact, and dust-obscured candidate
$\boldsymbol{z\simeq 8}$ galaxies discovered by JWST

Hollis B. Akins Department of Astronomy, The University of Texas at Austin,
2515 Speedway Blvd Stop C1400, Austin, TX 78712, USA  Caitlin M. Casey
Department of Astronomy, The University of Texas at Austin, 2515 Speedway Blvd
Stop C1400, Austin, TX 78712, USA  Cosmic Dawn Center (DAWN), Denmark  Natalie
Allen Cosmic Dawn Center (DAWN), Denmark  Niels Bohr Institute, University of
Copenhagen, Jagtvej 128, DK-2200, Copenhagen N, Denmark  Micaela B. Bagley
Department of Astronomy, The University of Texas at Austin, 2515 Speedway Blvd
Stop C1400, Austin, TX 78712, USA  Mark Dickinson NSF’s National Optical-
Infrared Astronomy Research Laboratory, 950 N. Cherry Ave., Tucson, AZ 85719,
USA  Steven L. Finkelstein Department of Astronomy, The University of Texas at
Austin, 2515 Speedway Blvd Stop C1400, Austin, TX 78712, USA  Maximilien
Franco Department of Astronomy, The University of Texas at Austin, 2515
Speedway Blvd Stop C1400, Austin, TX 78712, USA  Santosh Harish Laboratory for
Multiwavelength Astrophysics, School of Physics and Astronomy, Rochester
Institute of Technology, 84 Lomb Memorial Drive, Rochester, NY 14623, USA
Pablo Arrabal Haro NSF’s National Optical-Infrared Astronomy Research
Laboratory, 950 N. Cherry Ave., Tucson, AZ 85719, USA  Olivier Ilbert Aix
Marseille Université, CNRS, CNES, LAM, Marseille, France  Jeyhan S. Kartaltepe
Laboratory for Multiwavelength Astrophysics, School of Physics and Astronomy,
Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, NY
14623, USA  Anton M. Koekemoer Space Telescope Science Institute, 3700 San
Martin Drive, Baltimore, MD 21218, USA  Daizhong Liu Max-Planck-Institut für
Extraterrestrische Physik (MPE), Giessenbachstr. 1, D-85748 Garching, Germany
Arianna S. Long Hubble Fellow  Department of Astronomy, The University of
Texas at Austin, 2515 Speedway Blvd Stop C1400, Austin, TX 78712, USA  Henry
Joy McCracken Institut d’Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne
Université, 98 bis boulevard Arago, F-75014 Paris, France  Louise Paquereau
Institut d’Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Université, 98
bis boulevard Arago, F-75014 Paris, France  Casey Papovich 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  Nor
Pirzkal ESA/AURA Space Telescope Science Institute  Jason Rhodes Jet
Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove
Drive, Pasadena, CA 91001, USA  Brant E. Robertson Department of Astronomy and
Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa
Cruz, CA 95064, USA  Marko Shuntov Cosmic Dawn Center (DAWN), Denmark  Niels
Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200, Copenhagen N,
Denmark  Sune Toft Cosmic Dawn Center (DAWN), Denmark  Niels Bohr Institute,
University of Copenhagen, Jagtvej 128, DK-2200, Copenhagen N, Denmark  Guang
Yang 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  Guillermo Barro
Department of Physics, University of the Pacific, Stockton, CA 90340, USA
Laura Bisigello 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  Véronique Buat Aix Marseille Université, CNRS, CNES, LAM, Marseille,
France  Jaclyn B. Champagne Steward Observatory, University of Arizona, 933 N
Cherry Ave, Tucson, AZ 85721, USA  Olivia Cooper NSF Graduate Research Fellow
Department of Astronomy, The University of Texas at Austin, 2515 Speedway Blvd
Stop C1400, Austin, TX 78712, USA  Luca Costantin Centro de Astrobiología
(CAB/CSIC-INTA), Ctra. de Ajalvir km 4, Torrejón de Ardoz, E-28850 Madrid,
Spain  Alexander de la Vega Department of Physics and Astronomy, University of
California Riverside, 900 University Ave, Riverside, CA 92521, USA  Nicole E.
Drakos Department of Astronomy and Astrophysics, University of California,
Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA  Andreas Faisst
Caltech/IPAC, 1200 E. California Blvd., Pasadena, CA 91125, USA  Adriano
Fontana INAF, Osservatorio Astronomico di Roma, via di Frascati 33, 00078
Monte Porzio Catone, Italy  Seiji Fujimoto Hubble Fellow  Department of
Astronomy, The University of Texas at Austin, 2515 Speedway Blvd Stop C1400,
Austin, TX 78712, USA  Steven Gillman Cosmic Dawn Center (DAWN), Denmark  DTU-
Space, Technical University of Denmark, Elektrovej 327, DK-2800 Kgs. Lyngby,
Denmark  Carlos Gómez-Guijarro Université Paris-Saclay, Université Paris Cité,
CEA, CNRS, AIM, 91191, Gif-sur-Yvette, France  Ghassem Gozaliasl Department of
Physics, University of Helsinki, P.O. Box 64, FI-00014 Helsinki, Finland
Nimish P. Hathi Space Telescope Science Institute, 3700 San Martin Drive,
Baltimore, MD 21218, USA  Christopher C. Hayward Center for Computational
Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA
Michaela Hirschmann 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, 34131 Trieste, Italy  Benne W.
Holwerda Physics & Astronomy Department, University of Louisville, 40292 KY,
Louisville, USA  Shuowen Jin Cosmic Dawn Center (DAWN), Denmark  DTU-Space,
Technical University of Denmark, Elektrovej 327, DK-2800 Kgs. Lyngby, Denmark
Dale D. Kocevski Department of Physics and Astronomy, Colby College,
Waterville, ME 04901, USA  Vasily Kokorev Kapteyn Astronomical Institute,
University of Groningen, P.O. Box 800, 9700 AV Groningen, The Netherlands
Erini Lambrides NASA Goddard Space Flight Center, 8800 Greenbelt Rd,
Greenbelt, MD 20771, USA  NPP Fellow  Ray A. Lucas Space Telescope Science
Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  Georgios E. Magdis
Cosmic Dawn Center (DAWN), Denmark  Niels Bohr Institute, University of
Copenhagen, Jagtvej 128, DK-2200, Copenhagen N, Denmark  DTU-Space, Technical
University of Denmark, Elektrovej 327, DK-2800 Kgs. Lyngby, Denmark  Benjamin
Magnelli Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM,
91191, Gif-sur-Yvette, France  Jed McKinney Department of Astronomy, The
University of Texas at Austin, 2515 Speedway Blvd Stop C1400, Austin, TX
78712, USA  Bahram Mobasher Department of Physics and Astronomy, University of
California Riverside, 900 University Ave, Riverside, CA 92521, USA  Pablo G.
Pérez-González Centro de Astrobiología (CAB/CSIC-INTA), Ctra. de Ajalvir km 4,
Torrejón de Ardoz, E-28850 Madrid, Spain  R. Michael Rich Department of
Physics and Astronomy, UCLA, PAB 430 Portola Plaza, Box 951547, Los Angeles,
CA 90095-1547  Lise-Marie Seillé Aix Marseille Université, CNRS, CNES, LAM,
Marseille, France  Margherita Talia 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  C. Megan Urry Physics Department and
Yale Center for Astronomy & Astrophysics, Yale University, PO Box 208120, CT
06520, USA  Francesco Valentino European Southern Observatory, Karl-
Schwarzschild-Strasse 2, D-85748, Garching bei München, Germany  Cosmic Dawn
Center (DAWN), Denmark  Katherine E. Whitaker Department of Astronomy,
University of Massachusetts, Amherst, MA 01003, USA  Cosmic Dawn Center
(DAWN), Denmark  L. Y. Aaron Yung NASA Goddard Space Flight Center, 8800
Greenbelt Rd, Greenbelt, MD 20771, USA  NASA Postdoctoral Fellow  Jorge Zavala
National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo
181-8588, Japan

###### Abstract

We present a search for extremely red, dust-obscured, $z>7$ galaxies with
JWST/NIRCam+MIRI imaging over the first 20 arcmin2 of publicly-available Cycle
1 data from the COSMOS-Web, CEERS, and PRIMER surveys. Based on their red
color in F277W$-$F444W ($\sim 2.5$ mag) and detection in MIRI/F770W ($\sim$25
mag), we identify two galaxies—COS-z8M1 and CEERS-z7M1—which have best-fit
photometric redshifts of $z=8.5^{+0.3}_{-0.4}$ and $z=7.6^{+0.1}_{-0.1}$,
respectively. We perform SED 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_{\rm eff}\lesssim 200$
pc), highly obscured ($A_{V}\sim 1.5$–$2.5$), and, at our best-fit redshift
estimates, likely have strong $[{\rm O}\,\textsc{iii}]$+H$\beta$ 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 AGN. We identify a marginal (2.9$\sigma$) ALMA
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 modestly dust-
obscured galaxies at cosmic dawn.

††journal: (Affiliations can be found after the references)††facilities: HST
(ACS and WFC3), JWST (NIRCam and MIRI), Spitzer (IRAC and MIPS), Herschel
(PACS and SPIRE), SCUBA-2, VLA, ALMA. All the JWST data used in this work can
be found in MAST: . (catalog 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),
SourceExtractor (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).

##  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 (Naidu et al., 2022a; Finkelstein et
al., 2022; Donnan et al., 2023), and several have been spectroscopically
confirmed with JWST/NIRSpec (e.g. Fujimoto et al., 2023a; Larson et al., 2023;
Saxena et al., 2023; Curtis-Lake et al., 2022; Arrabal Haro 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 which 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; Harikane et al.,
2023a; Ferrara et al., 2022a; Menci et al., 2022; 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 are dusty star-forming
galaxies (DSFGs), which are ubiquitous at cosmic noon ($z\sim 2$), where they
dominate the star-formation rate density (SFRD) of the Universe (Casey,
Narayanan, & Cooray, 2014a; Madau & Dickinson, 2014). While DSFGs have been
detected out to $z\sim 6$–$7$, with obscured star-formation rates in excess of
$1000$ M⊙ yr-1 (Zavala et al., 2018; Marrone 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. Zavala et al., 2018; Marrone 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_{\rm
IR}\gtrsim 10^{12\mbox{--}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, b). At
present, most $z>7$ dust-continuum detections come from follow-up of UV-
selected LBGs (e.g. Watson et al., 2015; Bouwens et al., 2022; Fujimoto et
al., 2022a) or serendipitous detections (e.g. Fudamoto et al., 2021) with only
one robust dust-continuum-detected object at $z>8$: MACS0416-Y1 at $z=8.31$
(Tamura et al., 2019).111A2744-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). 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
constraints at this epoch are limited, the near and mid-IR has 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 (e.g. Labbé et al., 2023; Pérez-González et al., 2022; Rodighiero et
al., 2023). In particular, the mid-infrared coverage of JWST/MIRI allows
constraints on the full rest-frame optical 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 (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), CEERS (Finkelstein et al., 2022), and PRIMER-
COSMOS (P.I. J. Dunlop, GO#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\,\mu$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 (§4.1), sizes (§4.2), stellar masses (§4.3), and
far-infrared emission (§4.4). Finally, in Section 5 we discuss the
implications of these discoveries for galaxy formation within a $\Lambda{\rm
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~{}\mu$m JWST/NIRCam+MIRI imaging from three publicly
available Cycle 1 programs to identify candidate massive, dusty galaxies at
$z\gtrsim 8$. Table 2.2 provides the approximate $5\sigma$ depth for point
sources and the effective area (i.e. the NIRCam+MIRI overlapping area) for
each survey.

###  2.1 JWST Observations & Data Reduction

####  2.1.1 COSMOS-Web

COSMOS-Web is a large Cycle 1 treasury program imaging a contiguous $0.54$
deg2 in the COSMOS field with NIRCam and $0.2$ deg2 with MIRI in parallel
(Casey et al., 2022). As of this writing, 6 of the 152 visits have been
completed during observations executed in January 2023, constituting a
contiguous 77 arcmin2 ($\sim 4\%$ of the overall area), with 8.7 arcmin2 of
overlap between the NIRCam and MIRI coverage. The COSMOS-Web imaging includes
four NIRCam filters—F115W, F150W, F277W, and F444W—and one MIRI filter, F770W,
at an approximate $5\sigma$ depth of 26 AB mag.

The full details on the NIRCam and MIRI reduction process will be presented in
upcoming papers (M. Franco et al. in prep and S. Harish et al. in prep,
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 (CRDS)222jwst-
crds.stsci.edu pmap 0989 which corresponds to NIRCam instrument mapping imap
0232. The final mosaics are created in Stage 3 of the pipeline with a pixel
size of $0\farcs 03$/pixel. Astrometric calibration is conducted via the JWST
tweakreg procedure, with a reference catalog based on a HST/F814W $0\farcs
03$/pixel mosaic in the COSMOS field with astrometry tied to Gaia-EDR3 (Gaia
Collaboration, 2018). The median offset in RA and Dec 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$\farcs$06/pixel and aligned with ancillary
HST/F814W imaging of the region.

####  2.1.2 CEERS

The Cosmic Evolution Early Release Science survey (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 9 MIRI pointings in the
CANDELS Extended Groth Strip (EGS) field, alongside spectroscopy with NIRSpec
and NIRCam WFSS. Each CEERS/NIRCam pointing includes 7 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 arcmin2 of
overlap between MIRI/F770W and NIRCam imaging, at an approximate $5\sigma$
depth in F770W of $\sim 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.
(in prep).

####  2.1.3 PRIMER-COSMOS

The Public Release Imaging for Extragalactic Research (PRIMER) survey (P.I. J.
Dunlop, GO#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 has been observed,
constituting just $4.1$ arcmin2 of overlap between NIRCam+MIRI. PRIMER imaging
includes 8 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 and 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 January 2024.

###  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 Subaru/Hyper Suprime-Cam (HSC; Aihara et
al., 2019), and $YJHK_{s}$ imaging from the UltraVISTA survey (McCracken et
al., 2012). We additionally utilize the HST/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 is used in Weaver et al. (2022) with the exception of the UltraVISTA data
which is 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).

Table 1: Combined NIRCam+MIRI area and approximate 5$\sigma$ depths for the
three surveys included in this work.

Survey | Area | F277W | F444W | F770W
---|---|---|---|---
(arcmin2) | ($5\sigma$) | ($5\sigma$) | ($5\sigma$)
COSMOS-Web† | 8.7 | 28.3 | 28.2 | 26.0
CEERS | 7.8 | 29.2 | 28.6 | 27.1
PRIMER-COSMOS† | 4.1 | 28.9 | 28.7 | 26.0

††footnotetext: Table reflects data available as of March 2023. Total
NIRCam+MIRI overlap upon survey completion will be $650$ arcmin2 for COSMOS-
Web and $140$ arcmin2 for PRIMER-COSMOS.

References. — COSMOS-Web (Casey et al., 2022, GO#1727), CEERS (Finkelstein et
al., 2022, DD-ERS#1345), PRIMER (P.I. J. Dunlop; GO#1837)

Note. — Depths are computed in $0.3^{\prime\prime}$ ($0.6^{\prime\prime}$)
diameter apertures for NIRCam (MIRI).

###  2.3 Multiwavelength Catalogs

####  2.3.1 COSMOS-Web & PRIMER-COSMOS

For both COSMOS-Web and PRIMER-COSMOS we conduct source detection, perform
model-based photometry, and construct multi-band catalogs using
SourceXtractor++ (Bertin et al., 2020; Kümmel et al., 2020, hereafter SE++),
an updated version of the popular SExtractor package (Bertin & Arnouts, 1996).
The use of SE++ 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. We perform source detection on a
$\chi^{2}$ detection image constructed from all four NIRCam bands; priors for
the source centroid positions are determined from this detection image. For
each detected source in the $\chi^{2}$ image, SE++ then fits a Sérsic model
convolved with the filter-specific PSF in each of the measurement band. Here
we use model PSFs from WebbPSF (Perrin et al., 2012, 2014). The Sérsic model
parameters (centroid position, $n$, $R_{\rm eff}$, $b/a$) are fit jointly
between all bands, weighted by their respective S/N, while the total flux is
fit independently in each band. Details on the SE++ catalogs for COSMOS-Web
will be presented in M. Shuntov et al. (in prep.).

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 in circular apertures with
a diameter of 0.3” (for ACS/NIRCam), 0.6” (for MIRI), and 2” (for ground-based
data). We adopt the measured depths from Weaver et al. (2022) and Casey et al.
(2022); we report depths in relevant JWST filters in Table 2.2. 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 multi-band 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 source detection on the MIRI F770W
images using astropy/photutils. Before performing detection we convolve the
image with a $5\times 5$ pixel Gaussian smoothing kernel with FWHM of 2 pixels
in order to better identify faint objects. We set a detection threshold at 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 signal-to-noise, 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 $\sim 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

We construct empirical point spread functions (ePSFs) in each 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 which
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 slightly broader than
the corresponding WebbPSF models due to the nature of the mosaicking process.

##  3 Sample selection & vetting

We determine initial photometric redshift estimates for the entire NIRCam+MIRI
samples across all three imaging surveys using EAzY (Brammer et al., 2008).
EAzY computes linear combinations of pre-defined 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 multi-wavelength 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 6 templates from Larson et al. (2022). We allow the redshift to vary
from 0 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.

Figure 1: The $m_{277}-m_{444}$ vs. $m_{444}-m_{770}$ color-color diagram
indicating selection of $z\gtrsim 7$ dusty galaxies. Points are colored by
their best-fit photometric redshift from first-pass EAzY runs. The dashed
lines indicate our proposed color selection criteria, which captures 4
objects, all with photometric redshifts $\gtrsim 7$. We show in grey two model
tracks generated from bagpipes, for young (10 Myr) stellar populations with
$A_{V}=2.5$ and $3.5$ from $z\sim 1$–$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). Figure
2: The optical through mid-infrared SED of COS-z8M1. The top panels show
$1.5^{\prime\prime}$ square 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$\sigma$ upper limits) from all available
JWST bands (shown in dark red) as well as HST and ground-based bands (in light
red). We additionally show maximum a posteriori (MAP) model SEDs from cigale,
beagle, bagpipes, and prospector (for clarity we do not show the EAzY SED).
The inset panel 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.

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 dusty AGN whose SED continues to
rise in the mid-infrared 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 $A_{V}=3.5$ (dashed line). We also plot a QSO
model by reddening the 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 a rough color-color selection criteria to encompass the reddest
galaxies in our sample, which all have EAzY photometric redshifts $>7$.
Specifically, we use

$\displaystyle m_{277}-$ $\displaystyle m_{444}>1.8\qquad{\rm and}$ (1)
$\displaystyle m_{277}-$ $\displaystyle m_{444}>1.5(m_{444}-m_{770})-0.6$ (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
mid-infrared 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 this 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 which may be present in the catalog
photometry used in Figure 1. In particular, we use the astropy/photutils
package to perform aperture photometry with aperture diameters ranging from
$0.3$–$0.5$”. We correct the aperture flux based on the PSF flux falling
outside the circular aperture.

With this photometry, we re-run 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 Appendix A,
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.

Figure 3: Same as Figure 2, but for CEERS-z7M1. The NIRCam SW bands reveal a
blue UV slope ($\beta\sim-2.5$) at $1$–$2\,\mu$m despite the very red
continuum at $>2\,\mu$m. This could be from SF in an unobscured line of sight,
or an unobscured AGN component. The additional filters constrain the redshift
precisely to $z=7.6$.

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 (Gutkin et al., 2016; Chevallard & Charlot, 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; Whitler et al., 2023; Furtak et al., 2022). We
additionally include, in all three codes, a late starburst (with age $\sim
10$–$100$ Myr) to allow for bright emission lines from H ii regions. We allow
$A_{V}$ to vary from 0 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 non-parametric star-formation history and informed,
joint priors encoding empirical constraints on the redshifts, stellar masses,
and SFHs of observed galaxies (see Wang et al., 2023, Table 1). Finally, we
include nebular emission in all our fits. beagle uses nebular emission
templates from Gutkin et al. (2016) which 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 Boquien et al. (2019). While each of these fits use
slightly different physical assumptions, the 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 §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^{\prime\prime}$ square 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 higher $\chi^{2}$ by $>10$.

Table 2: Physical properties derived from SED fitting for the two dust-
obscured, $z\sim 8$ sources.

ID | R.A., Decl. | Photo-$z$ | Photometric | $\log M_{\star}$ | SFR${}_{10\,\mathrm{Myr}}$ | SFR${}_{100\,\mathrm{Myr}}$ | $A_{V}$
---|---|---|---|---|---|---|---
[J2000] | Code | Redshift | [${\rm M}_{\odot}$] | [${\rm M}_{\odot}~{}{\rm yr}^{-1}$] | [${\rm M}_{\odot}~{}{\rm yr}^{-1}$] | [mag]
COS-z8M1 |  09h59m40.18s, +02d17m32.19s |  prospector-$\beta$ | $8.5^{+0.3}_{-0.4}$ | $9.8^{+0.1}_{-0.1}$ | $520^{+110}_{-100}$ | $59^{+22}_{-14}$ | $1.6^{+0.3}_{-0.2}$
bagpipes | $8.4^{+0.4}_{-0.5}$ | $9.9^{+0.2}_{-0.2}$ | $320^{+140}_{-90}$ | $73^{+51}_{-30}$ | $2.4^{+0.2}_{-0.2}$
beagle | $8.2^{+0.4}_{-0.5}$ | $9.5^{+0.2}_{-0.2}$ | $270^{+120}_{-90}$ | $27^{+13}_{-8}$ | $1.9^{+0.3}_{-0.3}$
cigale | $8.7^{+0.2}_{-0.4}$ | $10.3^{+0.4}_{-0.3}$ | $480^{+290}_{-460}$ | $57^{+35}_{-24}$ | $1.1^{+0.3}_{-0.1}$
eazy | $8.4^{+0.3}_{-0.3}$ | … | … | … | …
CEERS-z7M1 |  14h19m43.08s, +52d53m16.48s |  prospector-$\beta$ | $7.60^{+0.07}_{-0.07}$ | $10.1^{+0.1}_{-0.1}$ | $1010^{+270}_{-320}$ | $120^{+36}_{-33}$ | $2.6^{+0.1}_{-0.2}$
bagpipes | $7.69^{+0.05}_{-0.06}$ | $9.7^{+0.2}_{-0.2}$ | $360^{+100}_{-150}$ | $43^{+30}_{-20}$ | $3.2^{+0.2}_{-0.1}$
beagle | $7.69^{+0.02}_{-0.02}$ | $9.8^{+0.1}_{-0.1}$ | $600^{+120}_{-100}$ | $61^{+12}_{-10}$ | $3.1^{+0.1}_{-0.1}$
cigale | $8.26^{+0.12}_{-0.64}$ | $10.2^{+0.3}_{-0.1}$ | $810^{+190}_{-270}$ | $86^{+21}_{-16}$ | $1.5^{+0.4}_{-0.1}$
eazy | $8.22^{+0.31}_{-0.70}$ | … | … | … | …

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.

COS-z8M1 (Figure 2), is only detected in F277W, F444W, F770W, and in
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., 2023), which would suggest a maximally old
stellar population which 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 broad-band photometry in high-redshift galaxies (e.g. Shim
et al., 2011; Stark et al., 2013; Faisst et al., 2016; McKinney et al., 2022;
Fujimoto et al., 2022b; Naidu et al., 2022b; Zavala et al., 2023; Arrabal Haro
et al., 2023), and longer-wavelength MIRI data is critical to constrain
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$ 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 LW 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 is 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 high stellar mass for this epoch ($\log M_{\star}/{\rm
M}_{\odot}>9.5$), a steeply rising SFH (with
SFR${}_{10\,\mathrm{Myr}}$/SFR${}_{100\,\mathrm{Myr}}\sim 10$), and a high
dust attenuation ($A_{V}\sim 1.5$–$2.5$ mags) for both candidates. 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
imfit333https://www.mpe.mpg.de/~erwin/code/imfit/ (Erwin, 2015) to
characterize the rest-frame optical sizes of COS-z8M1 and CEERS-z7M1. In
particular, we fit PSF-convolved models to the NIRCam/F444W images. While
F444W has the largest PSF of all 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_{\rm 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_{\rm
eff}\sim 0\farcs 01$–$0\farcs 07$ ($50$–$350$ pc at $z=8$). We find that the
resulting residuals are significant at the 3$\sigma$ level for $R_{\rm
eff}\gtrsim 200$ pc, but consistent with the background rms for $R_{\rm
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 (Robertson et al., 2022; Roberts-Borsani et al., 2022; Ono
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.

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_{\rm eff}\lesssim 200$ pc, as discussed in the text.

###  4.3 Stellar masses

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 sources of this
rarity, 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.

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 shows
stellar mass vs. 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$ arcmin2 of
combined NIRCam+MIRI imaging, suggests an approximate volume density of $\sim
2\times 10^{-5}$ 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) parametrization, modified with the redshift-dependent parameters
from Rodríguez-Puebla et al. (2016) to be consistent with a Planck cosmology
(see also discussion in Yung et al., 2023). We compute the halo mass $M_{\rm
halo}$ associated with a volume density $n(>M_{\rm halo})\approx 2\times
10^{-5}$ Mpc-3. This halo mass is then converted to a stellar mass assuming
the 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 don’t include uncertainties (e.g. from cosmic
variance) into 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 Far-Infrared Constraints and Dust Mass

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 an X,
offset $\sim 0.5^{\prime\prime}$ from the NIRCam source but consistent within
the positional uncertainty given the large beam size. The Ex-MORA flux, if
real, suggests a high IR luminosity of $\log L_{\rm IR}/L_{\odot}\sim
11.9^{+0.4}_{-0.3}$, comparable to local ULIRGs.

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 & 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 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. in prep), a
blind 2 mm ALMA survey of the COSMOS field designed to identify high-redshift
sub-mm galaxies (an extension of the original MORA survey presented in Casey
et al., 2021; Zavala et al., 2021).444CEERS-z7M1 is not covered by any
archival NOEMA pointing in the EGS field. 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_{\rm beam}\sim 1.5^{\prime\prime}$). At 2
mm, such a dataset is optimal for 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\,{\rm mm}}=0.19\pm 0.07$ mJy offset by
$\sim 0.5^{\prime\prime}$ from the NIRCam source. Figure 6 shows contours of
this marginal detection overlaid on the NIRCam/F444W image. While FIR
positional offsets of $\sim 0.5$–$1.0$” are not uncommon in high-redshift
DSFGs (e.g. Biggs & Ivison, 2008; Chapman et al., 2004; Hodge et al., 2012;
Inami et al., 2022), these are typically observed between the rest-frame UV
and the FIR (rather than rest-frame optical and FIR) and in more extended
sources. However, we note that the 2 mm centroid position is uncertain to
$\sim 0.5^{\prime\prime}$ (Condon, 1997; Ivison et al., 2007a) due to the low
S/N and $\sim 1.5^{\prime\prime}$ beam size, making it feasibly associated
with the NIRCam source; as such, we do not consider the offset significant.

Though it requries follow-up verification, we examine the implications of this
marginal detection for the FIR luminosity $L_{\rm 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 a MIR power law and a FIR modified blackbody (as
in Casey, 2012; Drew & Casey, 2022). We use a custom Markov Chain Monte Carlo
(MCMC) routine (based upon MCIRSED; Drew & Casey, 2022) with flat priors on
$L_{\rm IR}$, $T_{\rm dust}$, and $\beta$. In the absence of significant FIR
constraints beyond our one 2 mm data point, we fix $\alpha_{\rm MIR}=2.3$ and
$\lambda_{0}$ = 200 $\mu$m (following the recommendations in Drew & Casey,
2022); we allow $T_{\rm dust}$ to vary from $26$ K ($\approx T_{\rm 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 sub-mm (Casey, Narayanan, & Cooray, 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 CMB following da Cunha et al. (2013). We
note that this MCMC fitting is not intended to constrain $T_{\rm dust}$ or
$\beta$, but rather estimate $L_{\rm IR}$ marginalized over the uncertainty in
the other SED parameters.

This fit yields an IR luminosity of $\log L_{\rm IR}/L_{\odot}\sim
11.9^{+0.4}_{-0.3}$. Based on the Murphy et al. (2011) calibration, this
corresponds to an obscured SFR of $\sim 110^{+160}_{-60}\,{\rm
M}_{\odot}\,{\rm yr}^{-1}$. This implies a large fraction of obscured star-
formation ($\sim 99$%), as the upper limits in the rest-frame UV (uncorrected
for dust) imply SFR${}_{\rm UV}\lesssim 1~{}M_{\odot}~{}{\rm yr}^{-1}$.
Moreover, we apply the methodology outlined in Scoville et al. (2016, see also
section 3.3) to estimate the dust mass in COS-z8M1. This method depends on the
monochromatic flux at some wavelength (here 2 mm), the emissivity index (we
assume $\beta=1.75$), and the mass-weighted dust temperature, which is not the
same as the luminosity-weighted temperature that can be derived in SED
fitting. We adopt a temperature of 25 K, consistent with Scoville et al.
(2016). This yields a dust mass of $\log M_{\rm
dust}/M_{\odot}=8.2^{+0.5}_{-0.4}$, for $\beta\sim 1.5$–$2.4$.555We note that
this method assumes the dust is optically thin, which is likely not the case
at 2 mm (rest-frame $\sim 200\,\mu$m); however, in this case the dust mass
would be underestimated. This dust mass is consistent with the inferred
attenuation $A_{V}\sim 1.5$–$2.5$, and 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

JWST imaging 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. (2023) presented a sample of red, $z_{\rm phot}>8$
candidates from early CEERS imaging which exhibited “double breaks,” both the
Lyman and Balmer break. They interpret the red $2$–$4\,\mu$m color as tracing
the 4000 Å break and derive 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 $z>10$ candidates may be overestimated
due to differences in the IMF (e.g. Steinhardt et al., 2022) or contamination
from strong emission lines, e.g. from AGN (as discussed in Labbé et al., 2023;
Endsley et al., 2022a). Along those lines, Kocevski et al. (2023) presented
NIRSpec spectroscopy confirming one of the candidates in the Labbé et al.
(2023) work as a $z=5.6$ broad-line AGN, which reduces their estimated volume
density of massive $z>8$ galaxies. Furtak et al. (2022) identify an extremely
red, compact object at $z\sim 7.7$, similar to those in Labbé et al. (2023),
in deep JWST/NIRCam imaging as part of the UNCOVER survey (Bezanson et al.,
2022). Aided by lensing magnification, they constrain the effective radius to
$\lesssim 35$ pc, providing strong evidence for a low-luminosity quasar with
strong emission lines driving the red colors, though they do not rule out the
possibility of highly dust-obscured star-formation.

These results raise the question of whether these red, compact, high-$z$
objects are massive dust-obscured galaxies, reddened quasars, or both. Early
massive black holes have indeed been invoked as a possible explanation for the
high stellar masses inferred for the $z>9$ galaxy candidates being discovered
by JWST (Yuan et al., 2023; Brummel-Smith et al., 2023). Moreover, the objects
presented in Kocevski et al. (2023) and Furtak et al. (2022) show a red
continuum at $\lambda_{\rm rest}>3000$ Å but a blue UV slope at shorter
wavelengths, which has been interpreted as a composite galaxy+AGN signature.
This SED shape is similar to CEERS-z7M1, and consistent with COS-z8M1 given
the shallower depth of the NIRCam data. Even with spectroscopic data, Kocevski
et al. (2023) do not come to a conclusion as to the origin of the two
components. 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 vs.
quasar) for the candidates presented in this work.

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: a galaxy-only composite model with a massive ($\log
M_{\star}/M_{\odot}=10$), obscured ($A_{V}=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: a 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: a
galaxy+quasar model with a low-mass, unobscured galaxy combined with a bright
but heavily dust-reddened Type I quasar.

###  5.1 The high-z extreme of the DSFG population?

The $z>7$ regime is notoriously difficult for DSFG identification (Casey et
al., 2018a, b) and therefore far-infrared spectroscopic follow-up has been
limited to the brightest DSFGs (with $L_{\rm IR}\sim 10^{13}\,L_{\odot}$
Marrone et al., 2018; Endsley et al., 2022b; Fujimoto et al., 2022a). At the
same time, far-IR follow-up of UV-luminous galaxies, which are far more
numerous and more 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). 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.
Zavala et al., 2018; Marrone 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}\,$Mpc-3, is comparable to the (albeit poorly constrained) volume
density of $z\sim 7$ luminous infrared galaxies (LIRGS; $L_{\rm
IR}>10^{11}\,L_{\odot}$), for which integrating the IR luminosity function
yields $n\sim 0.5$–$10\times 10^{-5}$ 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 morphology 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 sub-millimeter galaxies (Condon et al., 1991; Ma et
al., 2016; Jin et al., 2022). Theoretical work has suggested that massive
galaxies at ultra-high-redshift may be able to form efficiently via feedback-
free starbursts (Dekel et al., 2023), which predicts compact, massive objects
by $z\sim 8$–$10$. Indeed, many $z>9$ galaxies being confirmed by JWST are
ultra-compact with effective radii $\lesssim 200$ pc (Tacchella et al., 2023;
Robertson et al., 2022; Roberts-Borsani et al., 2022; Ono et al., 2022). We
also note 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 sub-
millimeter galaxies (SMGs) 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$). 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 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_{\rm IR}/L_{\odot}=11.9^{+0.3}_{-0.4}$) and high dust
mass ($\log M_{\rm dust}/M_{\odot}=8.2^{+0.5}_{-0.4}$). These estimates are
highly uncertain and require direct sub-mm 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_{\rm
dust}/M_{\star}\approx 0.03$ is significantly higher than 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 alone 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 which 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_{\rm AGB}$ and $N_{\rm SN}$ by integrating the IMF from
$3$–$8\,M_{\odot}$ and $8$–$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 SN666This
assumes no dust destruction, likely an unphysical assumption. A more realistic
yield of $\sim 0.1$–$0.15\,M_{\odot}$ per SN (Leśniewska & Michałowski, 2019)
would give a lower dust mass by a factor of $\sim 10$. 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 SNe 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 can 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), in the wake around Wolf-Rayet stars (Lau et
al., 2021, 2022), in winds around an AGN accretion disk (Sarangi et al.,
2019), or in 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
SNe rate may also drive accelerated grain growth in the ISM by the
contribution additional seed metals at $z>10$. Taken altogether, these results
highlight the importance of deep sub-mm follow-up of these objects and future
samples of similarly-selected objects to constrain the dust masses and the
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_{\star}>10^{10}\,M_{\odot}$) galaxies
identified at $z>8$ with JWST has produced tensions with $\Lambda$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 derive constraints on the star-formation efficiency based on the
$z>10$ candidates identified in early JWST imaging (Harikane et al., 2022;
Naidu et al., 2022a; Finkelstein et al., 2022; Donnan et al., 2023; Harikane
et al., 2023a). Inayoshi et al. note 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.

Figure 8: Cumulative stellar mass density $\rho_{\star}$ vs. $M_{\star}$. 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. (2023) 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. (2023) sample. The light blue points show
the Labbé et al. (2023) 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$ SMF, though the uncertainties are large given the
identification of just two sources.

To explore the implications of these large stellar masses, we derive the
cumulative stellar mass density $\rho_{\star}(>M_{\star})$ vs. $M_{\star}$
based on the two candidates presented in this work. We estimate the effective
volume as the differential comoving volume integrated over a redshift bin from
$z\sim 7$–$9$ and scaled to the total survey area of $20.6$ arcmin2 (Table
2.2). 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 Poisson counting uncertainty. We note that the Poisson uncertainty
dominates over cosmic variance ($\sim 70$–$90\%$ vs. $\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. (2023); we compute $\rho_{\star}$ 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. (2023) is driven by the one $\log M_{\star}/M_{\odot}\sim 10.9$
candidate at $z\sim 7.5$; we show in light blue the $\rho_{\star}$ estimates
with this candidate removed, just to highlight its impact on $\rho_{\star}$.
The line and shaded region shows the stellar mass density corresponding to the
Schechter fits 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_{\star}$, we integrate the $\sim 7$, $8$, and $9$ SMFs
and take a weighted mean. We weight each curve by the volume associated with
redshift bins from $z\sim 7$–$7.5$, $7.5$–$8.5$, and $8.5$–$9$, respectively.

We find that the stellar mass density inferred from this work is consistent
with the $z\sim 7$–$9$ stellar mass function inferred from the UV-luminous,
IRAC-detected population. 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. (2023) sample is also fully consistent with the
Stefanon et al. (2021) SMF if the $\log M_{\star}/M_{\odot}\sim 10.9$
candidate were in fact 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_{\rm eff}\sim 800$ pc at
$M_{\star}\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 §5.1) compact
star-formation appears common at $z\sim 8$ (e.g. Ono et al., 2022).

Given the lack of spectroscopy or rest-frame MIR coverage, the present data
does not allow a robust constraint on any AGN contribution to the photometry.
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, our prospector-$\beta$ model
includes emission from the dusty torus, but yields highly unconstrained values
of $f_{\rm AGN}$ and $\tau_{\rm AGN}$. At this epoch, constraining the rest-
frame MIR SED will be difficult and require ultra-deep MIRI imaging in the
reddest wavelengths, which are also the least sensitive. Alternatively, deep
X-ray or radio imaging could indicate the presence of an obscured AGN.
However, neither source is detected in existing VLA radio images at $1.4$,
$3$, and $5$ GHz (Schinnerer et al., 2007; Smolčić et al., 2017; Ivison et
al., 2007b; Willner et al., 2006) or X-ray imaging from Chandra (Laird et al.,
2009; Nandra et al., 2015; Civano et al., 2016; Marchesi et al., 2016).

A second scenario involves strong emission lines from an AGN contaminating the
broadband photometry (e.g. Furtak et al., 2022; 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 for 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 22000$ Å) and even still achieve poor fits ($\chi_{\nu}^{2}\sim 1.8$).

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$ 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,
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 from 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 dust-enshrouded 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_{\star}\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 (shown in the
rightmost panel in Figure 7, and explored further in Barro et al. in prep.)
would significantly impact the derived physical properties of the galaxy. The
continuum in this case would be dominated by thermal emission from a dust-
obscured accretion disk, with little to no contribution from a host galaxy,
which would be completely sub-dominant. Unfortunately, constraining the
relative contribution from stellar vs. quasar continuum will require rest-
frame mid-infrared diagnostics or deep radio/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_{\rm BH}\sim 10^{7}\,M_{\odot}$ (based
on the $L_{\rm bol}$-$M_{\rm BH}$ relation) with a comparable host galaxy
stellar mass of $M_{\star}\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$ quasar luminosity function down to $M_{\rm
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 dominate over 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 contribution 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 SF depletion time for $z\gtrsim 4$ DSFGs ($\sim 100$–$300$
Myr; Swinbank et al., 2014; Aravena et al., 2016; Manning et al., 2022;
Williams et al., 2019) 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 presented 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. 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.3}_{-0.4}$ and $7.6^{+0.1}_{-0.1}$. The photometry for both sources
is likely impacted by strong emission lines, particularly [O iii]+H$\beta$ in
F444W and H$\alpha$+[N ii] in F560W. Both candidates are significantly more
dust obscured ($A_{V}\sim 2$–$3$) than other known $z\sim 8$ galaxies.

  2. 2.

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

  3. 3.

We infer stellar masses of $\sim 10^{10}\,M_{\odot}$, significantly higher
than 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. 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 $\sim 10^{12}\,L_{\odot}$, consistent with the
constraints on attenuation suggested by the JWST data. There are no sub-mm
constraints for CEERS-z7M1.

This work highlights the importance of long wavelength MIRI imaging for
characterization of massive, dust-obscured galaxies at $z>7$. Given the
remarkable sensitivity of MIRI, almost $1$ mag deeper than pre-flight
expectations (Casey et al., 2022), it becomes possible to constrain the $z>7$
galaxy population at rest-frame $\sim 1\,\mu$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 sub-mm follow-up will be necessary to determine the
nature of the dust-obscured population in the epoch of reionization.

## Acknowledgements

Support for this work was provided by NASA through grant 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.

## References

  * Aihara et al. (2019) Aihara, H., AlSayyad, Y., Ando, M., et al. 2019, PASJ, 71, 114, doi: 10.1093/pasj/psz103
  * Algera et al. (2023) Algera, H. S. B., Inami, H., Oesch, P. A., et al. 2023, MNRAS, 518, 6142, doi: 10.1093/mnras/stac3195
  * Anderson & King (2000) Anderson, J., & King, I. R. 2000, PASP, 112, 1360, doi: 10.1086/316632
  * Aravena et al. (2016) Aravena, M., Decarli, R., Walter, F., et al. 2016, ApJ, 833, 68, doi: 10.3847/1538-4357/833/1/68
  * Aretxaga et al. (2011) Aretxaga, I., Wilson, G. W., Aguilar, E., et al. 2011, MNRAS, 415, 3831, doi: 10.1111/j.1365-2966.2011.18989.x
  * Arrabal Haro et al. (2023) Arrabal Haro, P., Dickinson, M., Finkelstein, S. L., et al. 2023, arXiv, 2303.15431, doi: 10.48550/arXiv.2303.15431
  * Asano et al. (2013) Asano, R. S., Takeuchi, T. T., Hirashita, H., & Inoue, A. K. 2013, Earth, Planets and Space, 65, 213, doi: 10.5047/eps.2012.04.014
  * Astropy Collaboration (2013) Astropy Collaboration. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Bagley et al. (2022) Bagley, M. B., Finkelstein, S. L., Koekemoer, A. M., et al. 2022, arXiv, 2211.02495, doi: 10.48550/arXiv.2211.02495
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, Astronomy and Astrophysics Supplement, v.117, p.393-404, 117, 393, doi: 10.1051/aas:1996164
  * Bertin et al. (2020) Bertin, E., Schefer, M., Apostolakos, N., et al. 2020, Proceedings of Astronomical Data Analysis Software and Systems XXIX, 527, 461
  * Bezanson et al. (2022) Bezanson, R., Labbe, I., Whitaker, K. E., et al. 2022, arXiv, 2212.04026, doi: 10.48550/arXiv.2212.04026
  * Biggs & Ivison (2008) Biggs, A. D., & Ivison, R. J. 2008, MNRAS, 385, 893, doi: 10.1111/j.1365-2966.2008.12869.x
  * Boquien et al. (2019) Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103, doi: 10.1051/0004-6361/201834156
  * Bouwens et al. (2022) Bouwens, R., Stefanon, M., Brammer, G., et al. 2022, arXiv, 2211.02607, doi: 10.48550/arXiv.2211.02607
  * Boylan-Kolchin (2022) Boylan-Kolchin, M. 2022, arXiv, 2208.01611, doi: 10.48550/arXiv.2208.01611
  * Bradley et al. (2022) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2022, Zenodo, doi: 10.5281/zenodo.7419741
  * Brammer et al. (2008) Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503, doi: 10.1086/591786
  * Brummel-Smith et al. (2023) Brummel-Smith, C., Skinner, D., Sethuram, S. S., et al. 2023, arXiv, 2302.04882, doi: 10.48550/arXiv.2302.04882
  * Bruzual & Charlot (2003) Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000, doi: 10.1046/j.1365-8711.2003.06897.x
  * Burgarella et al. (2005) Burgarella, D., Buat, V., & Iglesias-Páramo, J. 2005, MNRAS, 360, 1413, doi: 10.1111/j.1365-2966.2005.09131.x
  * Byler et al. (2019) Byler, N., Dalcanton, J. J., Conroy, C., et al. 2019, ApJ, 158, 2, doi: 10.3847/1538-3881/ab1b70
  * Calura et al. (2014) Calura, F., Gilli, R., Vignali, C., et al. 2014, MNRAS, 438, 2765, doi: 10.1093/mnras/stt2329
  * Calura et al. (2017) Calura, F., Pozzi, F., Cresci, G., et al. 2017, MNRAS, 465, 54, doi: 10.1093/mnras/stw2749
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682, doi: 10.1086/308692
  * Carnall et al. (2018) Carnall, A. C., McLure, R. J., Dunlop, J. S., & Davé, R. 2018, MNRAS, 480, 4379, doi: 10.1093/mnras/sty2169
  * Casey (2012) Casey, C. M. 2012, MNRAS, 425, 3094, doi: 10.1111/j.1365-2966.2012.21455.x
  * Casey et al. (2018a) Casey, C. M., Hodge, J., Zavala, J. A., et al. 2018a, ApJ, 862, 78, doi: 10.3847/1538-4357/aacd11
  * Casey et al. (2014a) Casey, C. M., Narayanan, D., & Cooray, A. 2014a, Phys. Rep., 541, 45, doi: 10.1016/j.physrep.2014.02.009
  * Casey et al. (2014b) Casey, C. M., Scoville, N. Z., Sanders, D. B., et al. 2014b, ApJ, 796, 95, doi: 10.1088/0004-637X/796/2/95
  * Casey et al. (2018b) Casey, C. M., Zavala, J. A., Spilker, J., et al. 2018b, ApJ, 862, 77, doi: 10.3847/1538-4357/aac82d
  * Casey et al. (2019) Casey, C. M., Zavala, J. A., Aravena, M., et al. 2019, ApJ, 887, 55, doi: 10.3847/1538-4357/ab52ff
  * Casey et al. (2021) Casey, C. M., Zavala, J. A., Manning, S. M., et al. 2021, ApJ, 923, 215, doi: 10.3847/1538-4357/ac2eb4
  * Casey et al. (2022) Casey, C. M., Kartaltepe, J. S., Drakos, N. E., et al. 2022, arXiv, 2211.07865, doi: 10.48550/arXiv.2211.07865
  * Chapman et al. (2005) Chapman, S. C., Blain, A. W., Smail, I., & Ivison, R. J. 2005, ApJ, 622, 772, doi: 10.1086/428082
  * Chapman et al. (2004) Chapman, S. C., Smail, I., Windhorst, R., Muxlow, T., & Ivison, R. J. 2004, ApJ, 611, 732, doi: 10.1086/422383
  * Charlot & Fall (2000) Charlot, S., & Fall, S. M. 2000, ApJ, 539, 718, doi: 10.1086/309250
  * Chevallard & Charlot (2016) Chevallard, J., & Charlot, S. 2016, MNRAS, 462, 1415, doi: 10.1093/mnras/stw1756
  * Chevallard et al. (2013) Chevallard, J., Charlot, S., Wandelt, B., & Wild, V. 2013, MNRAS, 432, 2061, doi: 10.1093/mnras/stt523
  * Civano et al. (2016) Civano, F., Marchesi, S., Comastri, A., et al. 2016, ApJ, 819, 62, doi: 10.3847/0004-637X/819/1/62
  * Condon (1997) Condon, J. J. 1997, PASP, 109, 166, doi: 10.1086/133871
  * Condon et al. (1991) Condon, J. J., Huang, Z. P., Yin, Q. F., & Thuan, T. X. 1991, ApJ, 378, 65, doi: 10.1086/170407
  * Conroy et al. (2010) Conroy, C., White, M., & Gunn, J. E. 2010, ApJ, 708, 58, doi: 10.1088/0004-637X/708/1/58
  * Cooper et al. (2022) Cooper, O. R., Casey, C. M., Zavala, J. A., et al. 2022, ApJ, 930, 32, doi: 10.3847/1538-4357/ac616d
  * Curtis-Lake et al. (2022) Curtis-Lake, E., Carniani, S., Cameron, A., et al. 2022, arXiv, 2212.04568, doi: 10.48550/arXiv.2212.04568
  * da Cunha et al. (2013) da Cunha, E., Groves, B., Walter, F., et al. 2013, ApJ, 766, 13, doi: 10.1088/0004-637X/766/1/13
  * Dekel et al. (2023) Dekel, A., Sarkar, K. S., Birnboim, Y., Mandelker, N., & Li, Z. 2023, arXiv, 2303.04827, doi: 10.48550/arXiv.2303.04827
  * Donnan et al. (2023) Donnan, C. T., McLeod, D. J., Dunlop, J. S., et al. 2023, MNRAS, 518, 6011, doi: 10.1093/mnras/stac3472
  * Drakos et al. (2022) Drakos, N. E., Villasenor, B., Robertson, B. E., et al. 2022, ApJ, 926, 194, doi: 10.3847/1538-4357/ac46fb
  * Drew & Casey (2022) Drew, P. M., & Casey, C. M. 2022, ApJ, 930, 142, doi: 10.3847/1538-4357/ac6270
  * Dunne et al. (2011) Dunne, L., Gomez, H. L., da Cunha, E., et al. 2011, MNRAS, 417, 1510, doi: 10.1111/j.1365-2966.2011.19363.x
  * Dwek & Cherchneff (2011) Dwek, E., & Cherchneff, I. 2011, ApJ, 727, 63, doi: 10.1088/0004-637X/727/2/63
  * Eilers et al. (2020) Eilers, A.-C., Hennawi, J. F., Decarli, R., et al. 2020, ApJ, 900, 37, doi: 10.3847/1538-4357/aba52e
  * Endsley et al. (2022a) Endsley, R., Stark, D. P., Whitler, L., et al. 2022a, arXiv, 2208.14999, doi: 10.48550/arXiv.2208.14999
  * Endsley et al. (2022b) Endsley, R., Stark, D. P., Lyu, J., et al. 2022b, arXiv, 2206.00018, doi: 10.48550/arXiv.2206.00018
  * Erwin (2015) Erwin, P. 2015, ApJ, 799, 226, doi: 10.1088/0004-637X/799/2/226
  * Euclid Collaboration et al. (2022) Euclid Collaboration, Moneti, A., McCracken, H. J., et al. 2022, A&A, 658, A126, doi: 10.1051/0004-6361/202142361
  * Faisst et al. (2016) Faisst, A. L., Capak, P., Hsieh, B. C., et al. 2016, ApJ, 821, 122, doi: 10.3847/0004-637X/821/2/122
  * Ferland et al. (2013) Ferland, G. J., Porter, R. L., van Hoof, P. a. M., et al. 2013, Revista Mexicana de Astronomía y Astrofísica Vol. 49, pp. 137-163 (2013), 49, 137, doi: 10.48550/arXiv.1302.4485
  * Ferrara et al. (2022a) Ferrara, A., Pallottini, A., & Dayal, P. 2022a, arXiv, 2208.00720, doi: 10.48550/arXiv.2208.00720
  * Ferrara et al. (2022b) Ferrara, A., Sommovigo, L., Dayal, P., et al. 2022b, MNRAS, 512, 58, doi: 10.1093/mnras/stac460
  * Finkelstein et al. (2022) Finkelstein, S., Bagley, M. B., Ferguson, H. C., et al. 2022, arXiv, 2211.05792, doi: 10.48550/arXiv.2211.05792
  * Fu et al. (2017) Fu, H., Isbell, J., Casey, C. M., et al. 2017, ApJ, 844, 123, doi: 10.3847/1538-4357/aa7c63
  * Fudamoto et al. (2021) Fudamoto, Y., Oesch, P. A., Schouws, S., et al. 2021, Nature, 597, 489, doi: 10.1038/s41586-021-03846-z
  * Fujimoto et al. (2021) Fujimoto, S., Oguri, M., Brammer, G., et al. 2021, arXiv e-prints, 2101, arXiv:2101.01937
  * Fujimoto et al. (2022a) Fujimoto, S., Brammer, G. B., Watson, D., et al. 2022a, Nature, 604, 261, doi: 10.1038/s41586-022-04454-1
  * Fujimoto et al. (2022b) Fujimoto, S., Finkelstein, S. L., Burgarella, D., et al. 2022b, arXiv, doi: 10.48550/arXiv.2211.03896
  * Fujimoto et al. (2023a) Fujimoto, S., Arrabal Haro, P., Dickinson, M., et al. 2023a, arXiv, 2301.09482, doi: 10.48550/arXiv.2301.09482
  * Fujimoto et al. (2023b) Fujimoto, S., Kohno, K., Ouchi, M., et al. 2023b, arXiv, 2303.01658, doi: 10.48550/arXiv.2303.01658
  * Furtak et al. (2022) Furtak, L. J., Zitrin, A., Plat, A., et al. 2022, arXiv, 2212.10531, doi: 10.48550/arXiv.2212.10531
  * Gaia Collaboration (2018) Gaia Collaboration. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Glikman et al. (2006) Glikman, E., Helfand, D. J., & White, R. L. 2006, ApJ, 640, 579, doi: 10.1086/500098
  * Grogin et al. (2011) Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35, doi: 10.1088/0067-0049/197/2/35
  * Gutkin et al. (2016) Gutkin, J., Charlot, S., & Bruzual, G. 2016, MNRAS, 462, 1757, doi: 10.1093/mnras/stw1716
  * Harikane et al. (2022) Harikane, Y., Inoue, A. K., Mawatari, K., et al. 2022, ApJ, 929, 1, doi: 10.3847/1538-4357/ac53a9
  * Harikane et al. (2023a) Harikane, Y., Ouchi, M., Oguri, M., et al. 2023a, ApJS, 265, 5, doi: 10.3847/1538-4365/acaaa9
  * Harikane et al. (2023b) Harikane, Y., Zhang, Y., Nakajima, K., et al. 2023b, arXiv, 2303.11946, doi: 10.48550/arXiv.2303.11946
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hodge et al. (2012) Hodge, J. A., Carilli, C. L., Walter, F., et al. 2012, ApJ, 760, 11, doi: 10.1088/0004-637X/760/1/11
  * Holwerda et al. (2015) Holwerda, B. W., Bouwens, R., Oesch, P., et al. 2015, ApJ, 808, 6, doi: 10.1088/0004-637X/808/1/6
  * Hunter (2007) Hunter, J. D. 2007, Comput Sci Eng, 9, 90, doi: 10.1109/MCSE.2007.55
  * Inami et al. (2022) Inami, H., Algera, H. S. B., Schouws, S., et al. 2022, MNRAS, 515, 3126, doi: 10.1093/mnras/stac1779
  * Inayoshi et al. (2020) Inayoshi, K., Visbal, E., & Haiman, Z. 2020, ARA&A, 58, 27, doi: 10.1146/annurev-astro-120419-014455
  * Ivison et al. (2007a) Ivison, R. J., Greve, T. R., Dunlop, J. S., et al. 2007a, MNRAS, 380, 199, doi: 10.1111/j.1365-2966.2007.12044.x
  * Ivison et al. (2007b) Ivison, R. J., Chapman, S. C., Faber, S. M., et al. 2007b, ApJ, 660, L77, doi: 10.1086/517917
  * Izumi et al. (2021) Izumi, T., Matsuoka, Y., Fujimoto, S., et al. 2021, ApJ, 914, 36, doi: 10.3847/1538-4357/abf6dc
  * Jin et al. (2022) Jin, S., Daddi, E., Magdis, G. E., et al. 2022, A&A, 665, A3, doi: 10.1051/0004-6361/202243341
  * Johnson et al. (2021) Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22, doi: 10.3847/1538-4365/abef67
  * Jones et al. (2013) Jones, A. P., Fanciullo, L., Köhler, M., et al. 2013, A&A, 558, A62, doi: 10.1051/0004-6361/201321686
  * Kocevski et al. (2023) Kocevski, D. D., Onoue, M., Inayoshi, K., et al. 2023, arXiv, doi: 10.48550/arXiv.2302.00012
  * Koekemoer et al. (2007) Koekemoer, A. M., Aussel, H., Calzetti, D., et al. 2007, ApJS, 172, 196, doi: 10.1086/520086
  * Koekemoer et al. (2011) Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36, doi: 10.1088/0067-0049/197/2/36
  * Kroupa (2002) Kroupa, P. 2002, Science, 295, 82, doi: 10.1126/science.1067524
  * Kümmel et al. (2020) Kümmel, M., Bertin, E., Schefer, M., et al. 2020, 527, 29
  * Labbé et al. (2023) Labbé, I., van Dokkum, P., Nelson, E., et al. 2023, Nature, 1, doi: 10.1038/s41586-023-05786-2
  * Laird et al. (2009) Laird, E. S., Nandra, K., Georgakakis, A., et al. 2009, ApJS, 180, 102, doi: 10.1088/0067-0049/180/1/102
  * Laporte et al. (2017) Laporte, N., Ellis, R. S., Boone, F., et al. 2017, ApJ, 837, L21, doi: 10.3847/2041-8213/aa62aa
  * Larson et al. (2022) Larson, R. L., Hutchison, T. A., Bagley, M., et al. 2022, arXiv, 2211.10035, doi: 10.48550/arXiv.2211.10035
  * Larson et al. (2023) Larson, R. L., Finkelstein, S. L., Kocevski, D. D., et al. 2023, arXiv, 2303.08918, arXiv:2303.08918, doi: 10.48550/arXiv.2303.08918
  * Lau et al. (2021) Lau, R. M., Hankins, M. J., Kasliwal, M. M., et al. 2021, ApJ, 909, 113, doi: 10.3847/1538-4357/abd8cd
  * Lau et al. (2022) Lau, R. M., Hankins, M. J., Han, Y., et al. 2022, Nat Astron, 6, 1308, doi: 10.1038/s41550-022-01812-x
  * Leja et al. (2017) Leja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, ApJ, 837, 170, doi: 10.3847/1538-4357/aa5ffe
  * Leśniewska & Michałowski (2019) Leśniewska, A., & Michałowski, M. J. 2019, arXiv, 1904.11185, doi: 10.48550/ARXIV.1904.11185
  * Long et al. (2022) Long, A. S., Casey, C. M., Lagos, C. d. P., et al. 2022, arXiv, 2211.02072
  * Lovell et al. (2023) Lovell, C. C., Harrison, I., Harikane, Y., Tacchella, S., & Wilkins, S. M. 2023, MNRAS, 518, 2511, doi: 10.1093/mnras/stac3224
  * Ma et al. (2016) Ma, J., Gonzalez, Anthony. H., Vieira, J. D., et al. 2016, ApJ, 832, 114, doi: 10.3847/0004-637X/832/2/114
  * Ma et al. (2019) Ma, X., Hayward, C. C., Casey, C. M., et al. 2019, MNRAS, 487, 1844, doi: 10.1093/mnras/stz1324
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415, doi: 10.1146/annurev-astro-081811-125615
  * Manning et al. (2022) Manning, S. M., Casey, C. M., Zavala, J. A., et al. 2022, ApJ, 925, 23, doi: 10.3847/1538-4357/ac366a
  * Marchesi et al. (2016) Marchesi, S., Civano, F., Elvis, M., et al. 2016, ApJ, 817, 34, doi: 10.3847/0004-637X/817/1/34
  * Marconi et al. (2004) Marconi, A., Risaliti, G., Gilli, R., et al. 2004, MNRAS, 351, 169, doi: 10.1111/j.1365-2966.2004.07765.x
  * Marrone et al. (2018) Marrone, D. P., Spilker, J. S., Hayward, C. C., et al. 2018, Nature, 553, 51, doi: 10.1038/nature24629
  * Martínez-González et al. (2021) Martínez-González, S., Silich, S., & Tenorio-Tagle, G. 2021, Monthly Notices of the Royal Astronomical Society, 507, 1175, doi: 10.1093/mnras/stab2190
  * Matsuoka et al. (2018) Matsuoka, Y., Strauss, M. A., Kashikawa, N., et al. 2018, ApJ, 869, 150, doi: 10.3847/1538-4357/aaee7a
  * McConnell & Ma (2013) McConnell, N. J., & Ma, C.-P. 2013, ApJ, 764, 184, doi: 10.1088/0004-637X/764/2/184
  * McCracken et al. (2012) McCracken, H. J., Milvang-Jensen, B., Dunlop, J., et al. 2012, A&A, 544, A156, doi: 10.1051/0004-6361/201219507
  * McKinney et al. (2022) McKinney, J., Finnerty, L., Casey, C., et al. 2022, arXiv, doi: 10.48550/arXiv.2301.00017
  * McKinney et al. (2023) McKinney, J., Manning, S. M., Cooper, O. R., et al. 2023, arXiv, 2304.07316, doi: 10.48550/arXiv.2304.07316
  * Menci et al. (2022) Menci, N., Castellano, M., Santini, P., et al. 2022, ApJ, 938, L5, doi: 10.3847/2041-8213/ac96e9
  * Michałowski (2015) Michałowski, M. J. 2015, A&A, 577, A80, doi: 10.1051/0004-6361/201525644
  * Morishita et al. (2022) Morishita, T., Roberts-Borsani, G., Treu, T., et al. 2022, arXiv, 2211.09097, doi: 10.48550/arXiv.2211.09097
  * Murphy et al. (2011) Murphy, E. J., Condon, J. J., Schinnerer, E., et al. 2011, ApJ, 737, 67, doi: 10.1088/0004-637X/737/2/67
  * Murray et al. (2013) Murray, S. G., Power, C., & Robotham, A. S. G. 2013, Astronomy and Computing, Volume 3, article id. 23., 3, 23, doi: 10.1016/j.ascom.2013.11.001
  * Naidu et al. (2022a) Naidu, R. P., Oesch, P. A., van Dokkum, P., et al. 2022a, ApJ, 940, L14, doi: 10.3847/2041-8213/ac9b22
  * Naidu et al. (2022b) Naidu, R. P., Oesch, P. A., Setton, D. J., et al. 2022b, arXiv, 2208.02794, doi: 10.48550/arXiv.2208.02794
  * Nandra et al. (2015) Nandra, K., Laird, E. S., Aird, J. A., et al. 2015, ApJS, 220, 10, doi: 10.1088/0067-0049/220/1/10
  * Narayanan et al. (2018) Narayanan, D., Davé, R., Johnson, B. D., et al. 2018, MNRAS, 474, 1718, doi: 10.1093/mnras/stx2860
  * Narayanan et al. (2015) Narayanan, D., Turk, M., Feldmann, R., et al. 2015, Nature, 525, 496, doi: 10.1038/nature15383
  * Noll et al. (2009) Noll, S., Burgarella, D., Giovannoli, E., et al. 2009, A&A, 507, 1793, doi: 10.1051/0004-6361/200912497
  * Nozawa et al. (2014) Nozawa, T., Yoon, S.-C., Maeda, K., et al. 2014, ApJ, 787, L17, doi: 10.1088/2041-8205/787/2/L17
  * Oke (1974) Oke, J. B. 1974, The Astrophysical Journal Supplement Series, 27, 21, doi: 10.1086/190287
  * Ono et al. (2013) Ono, Y., Ouchi, M., Curtis-Lake, E., et al. 2013, ApJ, 777, 155, doi: 10.1088/0004-637X/777/2/155
  * Ono et al. (2022) Ono, Y., Harikane, Y., Ouchi, M., et al. 2022, arXiv, 2208.13582, doi: 10.48550/arXiv.2208.13582
  * Papovich et al. (2022) Papovich, C., Cole, J., Yang, G., et al. 2022, arXiv, doi: 10.48550/arXiv.2301.00027
  * Pérez-González et al. (2022) Pérez-González, P. G., Barro, G., Annunziatella, M., et al. 2022, arXiv, 2211.00045, doi: 10.48550/arXiv.2211.00045
  * Perrin et al. (2014) Perrin, M. D., Sivaramakrishnan, A., Lajoie, C.-P., et al. 2014, in Proceedings of the SPIE, Vol. 9143, 91433X, doi: 10.1117/12.2056689
  * Perrin et al. (2012) Perrin, M. D., Soummer, R., Elliott, E. M., Lallo, M. D., & Sivaramakrishnan, A. 2012, in Proceedings of the SPIE, Vol. 8442, 84423D, doi: 10.1117/12.925230
  * Planck Collaboration (2020) Planck Collaboration. 2020, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Popping et al. (2017) Popping, G., Puglisi, A., & Norman, C. A. 2017, MNRAS, 472, 2315, doi: 10.1093/mnras/stx2202
  * Richardson et al. (2014) Richardson, C. T., Allen, J. T., Baldwin, J. A., Hewett, P. C., & Ferland, G. J. 2014, Monthly Notices of the Royal Astronomical Society, 437, 2376, doi: 10.1093/mnras/stt2056
  * Rigby et al. (2022) Rigby, J., Perrin, M., McElwain, M., et al. 2022, arXiv, 2207.05632, doi: 10.48550/arXiv.2207.05632
  * Roberts-Borsani et al. (2022) Roberts-Borsani, G., Treu, T., Chen, W., et al. 2022, arXiv, 2210.15639, doi: 10.48550/arXiv.2210.15639
  * Robertson et al. (2022) Robertson, B. E., Tacchella, S., Johnson, B. D., et al. 2022, arXiv, 2212.04480, arXiv:2212.04480, doi: 10.48550/arXiv.2212.04480
  * Rodighiero et al. (2023) Rodighiero, G., Bisigello, L., Iani, E., et al. 2023, Monthly Notices of the Royal Astronomical Society, 518, L19, doi: 10.1093/mnrasl/slac115
  * Rodríguez-Puebla et al. (2016) Rodríguez-Puebla, A., Behroozi, P., Primack, J., et al. 2016, MNRAS, 462, 893, doi: 10.1093/mnras/stw1705
  * Sarangi et al. (2019) Sarangi, A., Dwek, E., & Kazanas, D. 2019, ApJ, 885, 126, doi: 10.3847/1538-4357/ab46a9
  * Saxena et al. (2023) Saxena, A., Robertson, B. E., Bunker, A. J., et al. 2023, arXiv, 2302.12805, doi: 10.48550/arXiv.2302.12805
  * Schinnerer et al. (2007) Schinnerer, E., Smolčić, V., Carilli, C. L., et al. 2007, ApJS, 172, 46, doi: 10.1086/516587
  * Schouws et al. (2022) Schouws, S., Stefanon, M., Bouwens, R., et al. 2022, ApJ, 928, 31, doi: 10.3847/1538-4357/ac4605
  * Scoville et al. (2016) Scoville, N., Sheth, K., Aussel, H., et al. 2016, ApJ, 820, 83, doi: 10.3847/0004-637X/820/2/83
  * Shim et al. (2011) Shim, H., Chary, R.-R., Dickinson, M., et al. 2011, ApJ, 738, 69, doi: 10.1088/0004-637X/738/1/69
  * Smolčić et al. (2012) Smolčić, V., Aravena, M., Navarrete, F., et al. 2012, A&A, 548, A4, doi: 10.1051/0004-6361/201219368
  * Smolčić et al. (2017) Smolčić, V., Novak, M., Bondi, M., et al. 2017, A&A, 602, A1, doi: 10.1051/0004-6361/201628704
  * Stark et al. (2013) Stark, D. P., Schenker, M. A., Ellis, R., et al. 2013, ApJ, 763, 129, doi: 10.1088/0004-637X/763/2/129
  * Stefanon et al. (2021) Stefanon, M., Bouwens, R. J., Labbé, I., et al. 2021, ApJ, 922, 29, doi: 10.3847/1538-4357/ac1bb6
  * Steidel et al. (1996) Steidel, C. C., Giavalisco, M., Pettini, M., Dickinson, M., & Adelberger, K. L. 1996, ApJ, 462, L17, doi: 10.1086/310029
  * Steinhardt et al. (2022) Steinhardt, C. L., Kokorev, V., Rusakov, V., Garcia, E., & Sneppen, A. 2022, arXiv, 2208.07879, doi: 10.48550/arXiv.2208.07879
  * Swinbank et al. (2014) Swinbank, A. M., Simpson, J. M., Smail, I., et al. 2014, MNRAS, 438, 1267, doi: 10.1093/mnras/stt2273
  * Tacchella et al. (2023) Tacchella, S., Eisenstein, D. J., Hainline, K., et al. 2023, arXiv, 2302.07234, doi: 10.48550/arXiv.2302.07234
  * Talia et al. (2021) Talia, M., Cimatti, A., Giulietti, M., et al. 2021, ApJ, 909, 23, doi: 10.3847/1538-4357/abd6e3
  * Tamura et al. (2019) Tamura, Y., Mawatari, K., Hashimoto, T., et al. 2019, ApJ, 874, 27, doi: 10.3847/1538-4357/ab0374
  * Tang et al. (2023) Tang, M., Stark, D. P., Chen, Z., et al. 2023, arXiv, 2301.07072, doi: 10.48550/arXiv.2301.07072
  * Tinker et al. (2008) Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709, doi: 10.1086/591439
  * Trump et al. (2023) Trump, J. R., Haro, P. A., Simons, R. C., et al. 2023, ApJ, 945, 35, doi: 10.3847/1538-4357/acba8a
  * Valiante et al. (2009) Valiante, R., Schneider, R., Bianchi, S., & Andersen, A. C. 2009, MNRAS, 397, 1661, doi: 10.1111/j.1365-2966.2009.15076.x
  * Vanden Berk et al. (2001) Vanden Berk, D. E., Richards, G. T., Bauer, A., et al. 2001, AJ, 122, 549, doi: 10.1086/321167
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Volonteri et al. (2015) Volonteri, M., Silk, J., & Dubus, G. 2015, ApJ, 804, 148, doi: 10.1088/0004-637X/804/2/148
  * Wang et al. (2023) Wang, B., Leja, J., Bezanson, R., et al. 2023, arXiv, 2302.08486, doi: 10.48550/arXiv.2302.08486
  * Watson et al. (2015) Watson, D., Christensen, L., Knudsen, K. K., et al. 2015, Nature, 519, 327, doi: 10.1038/nature14164
  * Weaver et al. (2022) Weaver, J. R., Kauffmann, O. B., Ilbert, O., et al. 2022, ApJS, 258, 11, doi: 10.3847/1538-4365/ac3078
  * Whitaker et al. (2017) Whitaker, K. E., Pope, A., Cybulski, R., et al. 2017, ApJ, 850, 208, doi: 10.3847/1538-4357/aa94ce
  * Whitler et al. (2023) Whitler, L., Endsley, R., Stark, D. P., et al. 2023, Monthly Notices of the Royal Astronomical Society, 519, 157, doi: 10.1093/mnras/stac3535
  * Williams et al. (2019) Williams, C. C., Labbe, I., Spilker, J., et al. 2019, ApJ, 884, 154, doi: 10.3847/1538-4357/ab44aa
  * Willner et al. (2006) Willner, S. P., Coil, A. L., Goss, W. M., et al. 2006, AJ, 132, 2159, doi: 10.1086/508202
  * Witstok et al. (2023) Witstok, J., Shivaei, I., Smit, R., et al. 2023, arXiv, 2302.05468, doi: 10.48550/arXiv.2302.05468
  * Yuan et al. (2023) Yuan, G.-W., Lei, L., Wang, Y.-Z., et al. 2023, arXiv, 2303.09391, doi: 10.48550/arXiv.2303.09391
  * Yung et al. (2023) Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., Wilkins, S. M., & Gardner, J. P. 2023, arXiv, 2304.04348, doi: 10.48550/arXiv.2304.04348
  * Zavala et al. (2018) Zavala, J. A., Montaña, A., Hughes, D. H., et al. 2018, Nature Astronomy, 2, 56, doi: 10.1038/s41550-017-0297-8
  * Zavala et al. (2021) Zavala, J. A., Casey, C. M., Manning, S. M., et al. 2021, ApJ, 909, 165, doi: 10.3847/1538-4357/abdb27
  * Zavala et al. (2023) Zavala, J. A., Buat, V., Casey, C. M., et al. 2023, ApJ, 943, L9, doi: 10.3847/2041-8213/acacfe
  * Zhang et al. (2018) Zhang, Z.-Y., Romano, D., Ivison, R. J., Papadopoulos, P. P., & Matteucci, F. 2018, Nature, 558, 260, doi: 10.1038/s41586-018-0196-x

##  Appendix A Additional candidate SEDs and photometry

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

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^{\prime\prime}$ 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 will
be needed for confirmation. The second candidate (COS-3627) suffers from a
“snowball” artifact in the F150W image; this artificially elevated flux yields
a blue $m_{150}-m_{277}$ color and drives the ultra high-$z$ solution as shown
in Figure 1. When we mask this artifact and recompute photometry (as is shown
in Fig. A.1), we find an incredibly broad redshift PDF with nearly equal
likelihood anywhere from $z\sim 4$ to $z\sim 12$—the source is too faint to
achieve a robust constraint.

Figure A.1: The optical-MIR SEDs of COS-939 and COS-3627, the two candidates
that fall in 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 grey. Table A.1: Photometry
from HST+JWST. All fluxes are in nJy.

ID |  HST/ACS |  JWST/NIRCam |  JWST/MIRI
---|---|---|---
F606W | F814W | F115W | F150W | F200W | F277W | F356W | F410M | F444W | F560W | F770W
COS-z8M1 | … | $(8\pm 9)$ | $(22\pm 25)$ | $(6\pm 19)$ | … | $41\pm 10$ | $150\pm 33^{\dagger}$ | … | $307\pm 11$ | … | $375\pm 73$
CEERS-z7M1 | $(-1\pm 5)$ | $(4\pm 5)$ | $11\pm 3$ | $(7\pm 4)$ | $(8\pm 3)$ | $17\pm 2$ | $69\pm 2$ | $158\pm 6$ | $188\pm 5$ | $451\pm 30$ | $338\pm 30$
COS-939 | … | $(-6\pm 9)$ | $(-28\pm 16)$ | $(-2\pm 13)$ | … | $22\pm 7$ | $(-91\pm 33)^{\dagger}$ | … | $118\pm 8$ | … | $301\pm 73$
COS-3627 | … | $(-3\pm 9)$ | $(16\pm 18)$ | $(0\pm 16)$ | … | $(14\pm 9)$ | $(15\pm 33)^{\dagger}$ | … | $68\pm 11$ | … | $162\pm 73$

††footnotetext: Reported flux is Spitzer/IRAC [3.6].