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

We present a search for extremely red, dust-obscured, z > 7 galaxies with JWST / NIRCam + MIRI imaging over the ﬁ rst 20 arcmin 2 of publicly available Cycle 1 data from the COSMOS-Web, CEERS, and PRIMER surveys. Based on their red color in F277W − F444W ( ∼ 2.5 mag ) and detection in MIRI / F770W ( ∼ 25 mag ) , we identify two galaxies, COS-z8M1 and CEERS-z7M1, that have best-ﬁ t photometric redshifts of


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 ∼ 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 restframe 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 ∼ 2), where they domi-nate the star-formation rate density (SFRD) of the Universe (Casey, Narayanan, & Cooray 2014a;Madau & Dickinson 2014).While DSFGs have been detected out to z ∼ 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 ∼ 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 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 IR 10 12-13 L ) DSFGs, as confirming lowerluminosity 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). 1 It remains unclear whether the perceived rarity of DSFGs in the first ∼ 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 ∼ 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 lowski 2015;Leśniewska & Micha lowski 2019).
In this paper, we present a search for extremely red, z 7 galaxies across the overlapping NIR-Cam+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 ∼ 8.Both are among the reddest galaxies identified in the entirety of their respective surveys, and detected in 7.7 µm MIRI imaging which provides robust constraints on their stellar mass.
The paper is organized as follows.In Section 2 we describe an overview of the imaging data used and the construction of photometric catalogs.In Section 3 we describe our sample selection and vetting of individual candidates.In Section 4 we present and characterize the two robust candidates, including their photometric redshifts ( §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 Λ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).

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

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 deg 2 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 arcmin 2 (∼ 4% of the overall area), with 8.7 arcmin 2 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σ 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) 2 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. 03/pixel.Astrometric calibration is conducted via the JWST tweakreg procedure, with a reference catalog based on a HST /F814W 0. 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. 06/pixel and aligned with ancillary HST /F814W imaging of the region.

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 arcmin 2 of overlap between MIRI/F770W and NIRCam imaging, at an approximate 5σ depth in F770W of ∼ 27 mag.We utilize the NIRCam and MIRI reductions produced by the CEERS team.The NIRCam reduction is described in detail in Bagley et al. (2022), and the MIRI reduction will be described in Yang et al. (in prep).

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 arcmin 2 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.

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 Y JHK s imaging from the UltraVISTA survey (McCracken et al. 2012).We additionally utilize the HST /ACS F814W imaging covering the entire COS-MOS 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 UltraV-ISTA data which is updated to the newest public data release 5.

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 groundbased data in COSMOS as well as high-resolution nearinfrared JWST imaging.We perform source detection on a χ 2 detection image constructed from all four NIR-Cam bands; priors for the source centroid positions are determined from this detection image.For each detected source in the χ 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(Perrin et al. , 2014)).The Sérsic model parameters (centroid position, n, R 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 1.The SE++ catalogs yield 1789 (585) objects detected with both NIRCam and MIRI for COSMOS-Web (PRIMER-COSMOS).

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 inversevariance-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 × 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 ∼ 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.

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.We show in grey two model tracks generated from bagpipes, for young (10 Myr) stellar populations with AV = 2.5 and 3.5 from z ∼ 1-9.We additionally show in red a QSO model with AV = 4, which, at z > 3, is redder in m444 − m770 than the galaxy models.The colorcolor selection shown here is optimized for selecting z 6-7 obscured galaxies with young ages and bright emission lines but rejecting objects with steeply rising MIR SEDs (e.g.obscured AGN).

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 χ 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 ∆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 observedframe colors and redshift and outline the sample selection of massive, dusty, and red z > 7 galaxies in the context of the full catalog.
We identify high-redshift, dusty galaxy candidates based on their colors in m 444 − m 770 and m 277 − m 444 .In particular, we design our selection to specifically target galaxies that are red in m 277 − m 444 but not as red in m 444 − m 770 ; this helps mitigate contamination from lower-redshift, extremely obscured sources or 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 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 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 aper-

P(z)
Figure 2. The optical through mid-infrared SED of COS-z8M1.The top panels show 1.5 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σ upper limits) from all available JWST bands (shown in dark red) as well as HST and groundbased 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.
ture 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.
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-τ 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 ∼ 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-β 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.Based on the sample selection and candidate vetting described in §3, we identify two robust z ∼ 8 dusty galaxies, with very red colors (m 277 − m 444 ∼ 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 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 χ 2 is similar to the other codes, and forcing a z < 7 solution with EAzY yields a higher χ 2 by > 10.
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 ∼ 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.Similarly, for CEERS-z7M1, the redshift solution is significantly constrained by the elevated flux in NIR-Cam/F410M+F444W and MIRI/F560W, likely due to bright [O iii] and Hα 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 ∼ 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 λ < 2 µm data gives z ∼ 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 /M > 9.5), a steeply rising SFH (with SFR 10 Myr /SFR 100 Myr ∼ 10), and a high dust attenuation (A V ∼ 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.

Rest-frame Optical Sizes
We use the 2D image fitting code imfit3 (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 χ 2 statistic, the resulting Sérsic parameters are poorly constrained.In order to provide a constraint on R eff , we fix the Sérsic index n = 1, the axis ratio b/a = 1 and the position angle θ = 0 (i.e. an exponential disk profile).We fit a series of models with R eff ∼ 0. 01-0.07 (50-350 pc at z = 8).We find that the resulting residuals are significant at the 3σ level for R eff 200 pc, but consistent with the background rms for R eff 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.

Stellar masses
We compare the stellar masses derived from SED fitting for these two sources to estimates for similar objects in the literature.Figure 5 shows the 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. ) 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 ∼ 20-30%.Both COS-z8M1 and CEERS-z7M1 appear to be among the most massive dust-obscured galaxies at this epoch.
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 dustobscured high-z population, with stellar masses 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 ∼ 7-9, across 20 arcmin 2 of combined NIRCam+MIRI imaging, suggests an approximate volume density of ∼ 2 × 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 halo associated with a volume density n(> M halo ) ≈ 2 × 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, .The orange stripe in Figure 5 shows the expected stellar masses for ∼ 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.

Far-Infrared Constraints and Dust Mass
Given that both COS-z8M1 and CEERS-z7M1 are very red, and our SED fits imply significant dust obscuration, they should be luminous in the FIR.Both the COSMOS and EGS fields have coverage with Spitzer /MIPS, Herschel /PACS & 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 highredshift sub-mm galaxies (an extension of the original MORA survey presented in Casey et al. 2021;Zavala et al. 2021). 4While Ex-MORA is still relatively shallow (5σ ∼ 1 mJy), it is deeper than the available SCUBA-2 data and has a smaller beam size (θ beam ∼ 1.5 ).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σ and flux density S 2 mm = 0.19 ± 0.07 mJy offset by ∼ 0.5 from the NIRCam source.Figure 6 shows contours of this marginal detection overlaid on the NIRCam/F444W image.While FIR positional offsets of ∼ 0.5-1.0"are not uncommon in high-redshift DSFGs (e.g.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 ∼ 0.5 (Condon 1997;Ivison et al. 2007a) due to the low S/N and ∼ 1.5 beam size, making it feasibly associated with the NIRCam source; as such, we do not consider the offset significant.
Though it requries follow-up verification, we examine the implications of this marginal detection for the FIR luminosity L IR .The only constraining power in the FIR SED comes from the marginal 2 mm Ex-MORA detection and the 850 µm upper limit from SCUBA-2 (3σ = 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 IR , T dust , and β.In the absence of significant FIR constraints beyond our one 2 mm data point, we fix α MIR = 2.3 and λ 0 = 200 µm (following the recommendations in Drew & Casey 2022); we allow T dust to vary from 26 K (≈ T CMB at z = 8.5) to 90 K and β to vary from 1.5 to 2.4.Due to the very negative Kcorrection 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 dust or β, but rather estimate L IR marginalized over the uncertainty in the other SED parameters.
This fit yields an IR luminosity of log L IR /L ∼ 11.9 +0.4  −0.3 .Based on the Murphy et al. (2011) calibration, this corresponds to an obscured SFR of ∼ 110 +160 −60 M yr −1 .This implies a large fraction of obscured star-formation (∼ 99%), as the upper limits in the rest-frame UV (uncorrected for dust) imply SFR UV 1 M yr −1 .Moreover, we apply the methodology outlined in Scoville et al. (2016, see also Casey et al. 2019 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 β = 1.75), and the massweighted 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 dust /M = 8.2 +0.5 −0.4 , for β ∼ 1.5-2.4. 5 This dust mass is consistent with the inferred attenuation A V ∼ 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.

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 phot > 8 candidates from early CEERS imaging which exhibited "double breaks," both the Lyman and Balmer break.They interpret the red 2-4 µm color as tracing the 4000 Å break and derive very large masses of M 10 10 M , potentially in excess of limits from Λ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 ∼ 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 35 pc, providing strong evidence for a lowluminosity 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 λ 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.

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 farinfrared spectroscopic follow-up has been limited to the brightest DSFGs (with L IR ∼ 10 13 L 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 ∼ 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.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 /M = 10), obscured (AV = 4) galaxy combined with a low-mass (log M /M = 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.
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 ∼ 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 ∼ 8-10.Indeed, many z > 9 galaxies being confirmed by JWST are ultra-compact with effective radii 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 ∼ 7.5-8.5, which may not be the best tracer of the morphology of the stellar continuum.
The composite blue+red SED of CEERS-z7M1 could be due to a patchy distribution of dust in an overall very dust-obscured galaxy.In fact, many observed submillimeter galaxies (SMGs) at lower-redshift (z 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 /M ∼ 7, A V = 0) emits alongside a dustobscured massive galaxy (log M /M ∼ 10, A V ∼ 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 dustobscured galaxies, which are expected to have complex star-dust geometry (Ma et al. 2019).

Implications for dust production at z > 8
The interpretation of these candidates as high-z dustobscured galaxies would imply an early buildup of dust.In particular, the marginal 2 mm flux from the Ex-MORA survey, if real and associated with COS-z8M1, suggests a high IR luminosity (log L IR /L = 11.9 +0.3 −0.4 ) and high dust mass (log M dust /M = 8.2 +0.5 −0.4 ).These estimates are highly uncertain and require direct submm 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 dust /M ≈ 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(Calura et al. , 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 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 lowski (2015), we estimate N AGB and N SN by integrating the IMF from 3-8 M and 8-40 M , 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 per SN6 and 0.04 M per AGB star (Micha lowski 2015).Based on the derived stellar masses of COS-z8M1 and CEERS-z7M1 (Table 2), we derive maximum dust masses of ∼ 10 8 M from SNe production and ∼ 10 7 M 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 (∼ 0.3 Z ), ISM grain growth can dominate over stellar dust production and can form ∼ 10 7 M of dust in ∼ 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(Lau et al. , 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 highmass 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 ∼ 8 The apparent ubiquity of massive (M > 10 10 M ) galaxies identified at z > 8 with JWST has produced tensions with ΛCDM (Boylan-Kolchin 2022; Menci et al. 2022).In the dust-obscured galaxy interpretation, the observed photometry implies stellar masses of ∼ 10 10 M 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, ∼ 25%.This is consistent with results of Inayoshi et al. (2020), who derive constraints on the starformation 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.
To explore the implications of these large stellar masses, we derive the cumulative stellar mass density ρ (> M ) vs. M 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 ∼ 7-9 and scaled to the total survey area of 20.6 arcmin 2 (Table 1).Error bars include the cosmic variance uncertainty, which we compute from halo number counts in the DREaM simulation (Drakos et al. 2022) as well as Poisson counting uncertainty.We note that the Poisson uncertainty dominates over cos- mic variance (∼ 70-90% vs. ∼ 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 ρ in the same manner using their full z ∼ 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 /M ∼ 10.9 candidate at z ∼ 7.5; we show in light blue the ρ estimates with this candidate removed, just to highlight its impact on ρ .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 UVluminous, IRAC-detected galaxies.To produce a curve comparable to our derived ρ , we integrate the ∼ 7, 8, and 9 SMFs and take a weighted mean.We weight each curve by the volume associated with redshift bins from z ∼ 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 ∼ 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 /M ∼ 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.

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 ∼ 7-10 star-forming galaxies (Ono et al. 2013;Holwerda et al. 2015), we might expect R eff ∼ 800 pc at M ∼ 10 10 M , significantly larger than the upper limit inferred from the F444W imaging (∼ 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 (∼ 200-300 pc), and (as discussed in §5.1) compact star-formation appears common at z ∼ 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 ∼ 2 µm, which is unconstrained by the present depth and filter coverage.Indeed, our prospector-β model includes emission from the dusty torus, but yields highly unconstrained values of f AGN and τ 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α 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 ∼ 5 but require unphysically high emission line EWs to match the observed photometry ([O iii]+Hβ EW ∼ 22000 Å) and even still achieve poor fits (χ 2 ν ∼ 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 (∼ 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 µ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 ∼ 10 10 M 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 BH ∼ 10 7 M (based on the L bol -M BH relation) with a comparable host galaxy stellar mass of M ∼ 10 7 M , 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, ∼ 1000 times rarer than these sources as inferred by integrating the Matsuoka et al. (2018) z = 6 quasar luminosity function down to M UV ∼ −18.At the same time, recent JWST /NIRspec results have revealed an abundant population of broad-line AGN in z ∼ 7 UVfaint 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 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 4 DSFGs (∼ 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 ∼ 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]).

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.
2. We find that neither source is resolved in NIR-Cam/F444W, constraining the rest-frame size to R eff 200 pc.
3. We infer stellar masses of ∼ 10 10 M , significantly higher than known dust-continuumdetected galaxies at z > 8 and similar to some of the most massive z > 8 galaxy candidates yet identified by JWST.The inferred stellar mass density is consistent within the uncertainty with expectations from the UV-luminous population.
4. We identify a marginal, 2.9σ 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 ∼ 10 12 L , 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 ∼ 1 µm.Future efforts to explore the dust-obscured population at this epoch will benefit greatly from deep MIRI imaging in multiple filters.While the completion of the full COSMOS-Web and PRIMER surveys will likely result in the detection of dozens more of these objects, spectroscopy and sub-mm follow-up will be necessary to determine the nature of the dust-obscured population in the epoch of reionization.

Figure 1 .
Figure 1.The m277 − m444 vs. m444 − m770 color-color diagram indicating selection of z 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 7.We show in grey two model tracks generated from bagpipes, for young (10 Myr) stellar populations with AV = 2.5 and 3.5 from z ∼ 1-9.We additionally show in red a QSO model with AV = 4, which, at z > 3, is redder in m444 − m770 than the galaxy models.The colorcolor selection shown here is optimized for selecting z 6-7 obscured galaxies with young ages and bright emission lines but rejecting objects with steeply rising MIR SEDs (e.g.obscured AGN).

Figure 3 .
Figure 3. Same as Figure2, but for CEERS-z7M1.The NIRCam SW bands reveal a blue UV slope (β ∼ −2.5) at 1-2 µm despite the very red continuum at > 2 µ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.

Figure 4 .
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 eff 200 pc, as discussed in the text.

M 3 Figure 5 .
Figure5.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 starforming 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 ∼ 20-30%.Both COS-z8M1 and CEERS-z7M1 appear to be among the most massive dust-obscured galaxies at this epoch.

Figure 6 .
Figure 6.Contours of the marginal 2 mm emission from Ex-MORA overlaid on the JWST /F444W image.The location of the 2.9σ peak is marked with an X, offset ∼ 0.5 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 LIR/L ∼ 11.9 +0.4 −0.3 , comparable to local ULIRGs.

Figure 7 .
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 /M = 10), obscured (AV = 4) galaxy combined with a low-mass (log M /M = 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.

Figure 8 .
Figure 8. Cumulative stellar mass density ρ vs. M .The red points indicate the cumulative stellar mass density inferred from the detection of COS-z8M1 and CEERS-z7M1.We additionally plot the result for z ∼ 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 ∼ 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 ∼ 7-9 SMF, though the uncertainties are large given the identification of just two sources.

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

Table 2 .
Kroupa (2002)erties derived from SED fitting for the two dust-obscured, z ∼ 8 sources.Note-Stellar masses and SFRs from beagle and cigale are corrected by a factor of 1.12 to be consistent with our assumption of aKroupa (2002)IMF.
(Fujimoto et al. 2023a;Saxena et al. 2023;Trump et al. 2023ger-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 Figure2.At 7 < z < 9, [O iii] 5007 Å emission falls into F444W while Hα falls blueward of F770W, yielding an extreme color differential.A high [O iii]+Hβ EW of ∼ 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).