Extremely Red Galaxies at z = 5 – 9 with MIRI and NIRSpec: Dusty Galaxies or Obscured Active Galactic Nuclei?

We study a new population of extremely red objects ( EROs ) recently discovered by the James Webb Space Telescope ( JWST ) based on their NIRCam colors F277W − F444W > 1.5 mag. We ﬁ nd 37 EROs in the Cosmic Evolution Early Release Science Survey ( CEERS ) ﬁ eld with F444W < 28 mag and photometric redshifts between 5 < z < 7, with median = -+ z 6.9 1.61.0 . Surprisingly, despite their red long-wavelength colors, these EROs have blue short-wavelength colors ( F150W − F200W ∼ 0 mag ) indicative of bimodal spectral energy distributions ( SEDs ) with a red, steep slope in the rest-frame optical, and a blue, ﬂ at slope in the rest-frame UV. Moreover, all these EROs are unresolved, point-like sources in all NIRCam bands. We analyze the SEDs of eight of them with MIRI and NIRSpec observations using stellar population models and active galactic nucleus ( AGN ) templates. We ﬁ nd that dusty galaxies or obscured AGNs provide similarly good SED ﬁ ts but different stellar properties: massive and dusty,   M M log ∼ 10 and A V  3 mag, or low mass and obscured,   M M log ∼ 7.5 and A V ∼ 0 mag, hosting an obscured quasi-stellar object ( QSO ) . SED modeling does not favor either scenario, but their unresolved sizes are more suggestive of AGNs. If any EROs are con ﬁ rmed to have   M M

z 6.9 1.6 1.0 .Surprisingly, despite their red long-wavelength colors, these EROs have blue short-wavelength colors (F150W − F200W ∼ 0 mag) indicative of bimodal spectral energy distributions (SEDs) with a red, steep slope in the rest-frame optical, and a blue, flat slope in the rest-frame UV.Moreover, all these EROs are unresolved, point-like sources in all NIRCam bands.We analyze the SEDs of eight of them with MIRI and NIRSpec observations using stellar population models and active galactic nucleus (AGN) templates.We find that dusty galaxies or obscured AGNs provide similarly good SED fits but different stellar properties: massive and dusty,   M M log ∼ 10 and A V  3 mag, or low mass and obscured,   M M log ∼ 7.5 and A V ∼ 0 mag, hosting an obscured quasi-stellar object (QSO).SED modeling does not favor either scenario, but their unresolved sizes are more suggestive of AGNs.If any EROs are confirmed to have   M M log  10.5, it would increase the pre-JWST number density at z > 7 by up to a factor ∼60. Similarly, if they are QSOs with luminosities in the L bol > 10 45-46 erg s −1

Introduction
The extraordinary capabilities of the James Webb Space Telescope (JWST) provide the opportunity to completely transform our understanding of the high-redshift Universe.The enhanced photometric sensitivity and spatial resolution at mid-infrared wavelengths relative to the Hubble Space Telescope (HST) or Spitzer have enabled, in the first few months of operations, a number of studies that have pushed the limits of the youngest and most distant galaxies detected in the epoch of reionization (e.g., Castellano et al. 2022;Naidu et al. 2022;Finkelstein et al. 2023a;Pérez-González et al. 2023a, 2023b;Adams et al. 2023;Whitler et al. 2023) as well as expanded our identification of more massive galaxies up to z ∼ 6 and beyond (e.g., Tacchella et al. 2022;Endsley et al. 2023;Labbé et al. 2023;Nelson et al. 2023).In the process, these papers have started to reveal the nature of the most massive galaxies that were previously undetected by HST (HST-dark) and detected only by Spitzer/IRAC, longer radio, and submillimeter wavelengths (Barrufet et al. 2023;Pérez-González et al. 2023a;Gómez-Guijarro et al. 2023;Rodighiero et al. 2023;Zavala et al. 2023), or not at all.
However, as we work our way toward a more complete census of the high-redshift Universe, there is a concern that some of these early estimates of the number density of galaxies or their (large) stellar masses could be in tension with model predictions (e.g., Boylan-Kolchin 2023; Ferrara et al. 2023;Mason et al. 2023).A potential caveat for these photometric studies is that as we probe galaxies in the first 1 Gyr of the lifetime of the Universe we might find a large number of young, low-mass galaxies with extreme emission lines and potentially large equivalent widths (EWs) of more than EW = 100-1000 Å, as suggested by early studies of faint z = 5-7 galaxies with Spitzer/IRAC (e.g., Egami et al. 2005;Eyles et al. 2007;Stark et al. 2009;González et al. 2014;Labbé et al. 2013).Such large EWs can make the Hβ, [O III], and Hα line fluxes boost the broad-and medium-band photometry in the JWST/NIRCam filters up to F444W, making them appear very red.The impact on the colors can affect both the photometric redshifts (e.g., Arrabal Haro et al. 2023) and the stellar population properties of these young, blue galaxies, introducing a bias toward older ages, more dust obscuration, and significantly larger masses ).Recent JWST-based papers have reported that emission lines with large EWs  1000 Å contaminating the NIRCam photometry are indeed a common occurrence (Endsley et al. 2021(Endsley et al. , 2023;;Matthee et al. 2023;Rinaldi et al. 2023), which may hamper the identification of true massive galaxies at z > 5. Another potential concern with massive galaxy selections based on extremely red colors is the contamination by obscured active galactic nuclei (AGNs).As shown also in IRAC-based studies, the red, power-law-like emission of an obscured AGN can also lead to very red optical to IR colors, which have been widely used to identify these galaxies in cosmological surveys (e.g., Alonso-Herrero et al. 2004;Stern et al. 2005;Lacy et al. 2007;Donley et al. 2008Donley et al. , 2012)).While the incidence of emission line or AGN contamination in color-selected samples at low to mid redshifts is only minor, the impact on JWST-based surveys is still unclear.
A way forward to overcome the degeneracy in the origin of colors in red galaxies (high-EW emission lines versus stellar or AGN continuum) is to obtain photometry in multiple bands and extend the coverage to longer wavelengths.Clear detections at wavelengths that are not severely affected by strong emission lines would be a clear confirmation of continuum emission.Likewise, long-wavelength (LW) detections probing the restframe near-infrared (NIR) of the galaxies can help distinguish between power-law AGN emission and the stellar 1.6 μm bump (Sawicki 2002;Donley et al. 2007).Observations with JWST/MIRI at λ > 5 μm help break both of these degeneracies.Similarly, JWST/NIRSpec can provide precise redshifts for these galaxies and help calibrate the impact of the emission lines in photometric observations.
In this paper, we use the first and second epochs of data from the Cosmic Evolution Early Release Science Survey (CEERS; Finkelstein et al. 2017) to identify candidates for massive dusty galaxies at z > 5 with very red colors in the LW NIRCam filters.Then, we focus on a subset of those galaxies with MIRI and NIRSpec observations to place better constraints on their redshifts and their emission at longer wavelengths, and we perform a detailed analysis of different spectral energy distribution (SED) modeling scenarios to determine the likelihood that they are blue high-EW galaxies, dusty massive galaxies, or obscured AGNs and the implications for the stellar masses and number densities of the sample in each case.
The paper is structured as follows.In Section 2, we describe the data reduction of the multiband NIRCam and MIRI imaging and the NIRSpec spectroscopy.We also describe the photometric measurements, catalog creation, and preliminary estimates of the photometric redshifts and stellar properties for the whole CEERS region.In Sections 3 and 4, we perform the extremely red object (ERO) color selection and we describe the colors, SEDs, photometric redshifts, and stellar masses of the sample selected that way.In Section 5, we perform a detailed SED modeling of a subset of eight EROs observed with MIRI and NIRSpec using a variety of SED models aimed at testing the dusty galaxy versus obscured-AGN scenarios and their implications on the stellar population properties.In Section 6, we discuss the likelihood of the different modeling scenarios based on the general properties of the EROs as well as their best-fit SEDs.Lastly, We summarize our results and discuss future prospects in Section 7.

Data
This paper is based on observations from CEERS, an early release science program (Finkelstein et al. 2017) which covers approximately 100 arcmin 2 of the Extended Groth Strip (EGS) with imaging and spectroscopy using coordinated, overlapping parallel observations by multiple JWST instruments.These images are available on the CEERS website 33 and on MAST as a High Level Science Product (doi:10.17909/z7p0-8481,Finkelstein et al. 2023b).Here we use the data acquired in 2022 June and December which comprise 10 NIRCam pointings in seven filters: three at short wavelengths (SW; F115W, F150W, and F200W), and four at LW (F277W, F356W, F410M, and F444W); and eight MIRI pointings in seven filters (F560W, F770W, F1000W, F1280W, F1500W, F1800W, and F2100W).Due to the nature of the CEERS parallel observations, some of the MIRI pointings are only observed either in the short (F560W and F770W) or long (F1000W to F2100W) wavelength filters and only six of them overlap with the NIRCam imaging.The names of these pointings in the APT observing file are 3, 6, 7, and 9, observed in F560W and F770W, and 5 and 8, observed at LW only.In addition to NIRCam imaging, pointings 3, 6, and 7 overlap with the NIRCam WFSS grism observations and two of the NIRSpec pointings named 9 and 10 in the APT file.
The NIRCam and MIRI data were calibrated using version 1.7.2 of the JWST Calibration Pipeline, reference files in pmap version 0214 (which includes a detector-to-detector-matched, improved absolute photometric calibration), with some additional modifications described in more detail in Finkelstein et al. (2023a) and Bagley et al. (2023) for NIRCam and Papovich et al. (2023) andG. Yang et al. (2024, in preparation) for MIRI.The reduced images are registered to the same world coordinate system reference frame (based on Gaia DR1.2; Gaia Collaboration et al. 2016) and coadded into single mosaics with pixel scales of 0 03 and 0 09 pixel −1 for NIRCam and MIRI, respectively.
The CEERS NIRSpec observations (P.Arrabal Haro 2024, in preparation) were processed using version 1.8.5 of the JWST Science Calibration Pipeline, with the Calibration Reference Data System mapping 1027 following similar procedures as in Fujimoto et al. (2023) and Kocevski et al. (2023).Briefly, we correct for 1/f detector noise, subtract the dark current and bias, and generate count-rate maps starting from the uncalibrated images.We apply a few additional custom steps to improve the treatment of cosmic-ray "snowballs."The resulting maps are processed with stage two of the pipeline to generate reduced 2D spectra with a rectified trace and flat slope.Custom extraction apertures are determined visually by inspecting the images for high signal-to-noise ratio continuum or emission lines.Lastly, we extract the 1D spectra boxcar apertures centered on the visually identified trace.

Source Extraction and Photometry
The source extraction and multiband photometric measurements were performed following the same methods as for the first epoch data described in detail in Finkelstein et al. (2023a).Briefly, photometry was computed on point-spread function (PSF)-matched images using SExtractor (Bertin & Arnouts 1996) v2.25.0 in two-image mode, with an inverse-variance weighted combination of the PSF-matched F277W and F356W images as the detection image.Photometry was measured in all seven of the NIRCam bands observed by CEERS, as well as the F606W, F814W, F105W, F125W, F140W, and F160W HST bands using data obtained by the CANDELS and 3D-HST surveys (Brammer et al. 2011;Grogin et al. 2011;Koekemoer et al. 2011).

Circular Aperture Photometry
We recompute the photometry of the subsample of objects studied in Section 5 using smaller circular apertures to improve the precision in the photometric errors and to avoid potential photometric contamination by nearby sources or background subtraction problems.Given that the nature of our galaxies is very homogeneous, and all sources analyzed in this paper are barely resolved or unresolved (see Section 4), photometric apertures with a 0 4 diameter were the most adequate to obtain the most precise and reliable SEDs.Photometry was measured in original and PSF-matched images, and after applying aperture corrections for point-like sources for the former, we arrived at consistent colors within at least half the value of the photometric corrections.

Photometric Redshifts and Stellar Population Properties
We estimate photometric redshifts for the whole parent catalog by fitting the multiband SEDs using the code EAZYpy (Brammer et al. 2008).The code fits nonnegative linear combinations of templates to the observed data to derive probability distribution functions (PDFs).Here we use the default template set "tweak fsps QSF 12 v3" which consists of a set of 12 templates derived from the stellar population synthesis code FSPS (Conroy et al. 2010).As a result, in addition to the photometric redshift the code also provides an estimate of the stellar mass as well as the dust attenuation.In addition, we also estimate stellar population properties by fitting the optical and NIR SEDs using FAST (Kriek et al. 2009), assuming Bruzual & Charlot (2003) stellar population synthesis models, following a Chabrier (2003) initial mass function (IMF), a delayed exponential star formation history (SFH), and the Calzetti et al. (2000) dust law with attenuation 0 < A V < 4 mag.

ERO Color Criterion
We identify extremely red galaxies at high redshift using a single color cut of F277W − F444W > 1.5 mag.This method is similar to the traditional ERO (R − K; e.g., McCarthy et al. 2004) or IERO (K − [4.5]; e.g., Wang et al. 2012;Caputi 2013;Stefanon et al. 2013) selections, which use red optical to NIR colors to find massive, dusty, or quiescent galaxies with strong Balmer or 4000 Å breaks at z  3.With the arrival of JWST, this technique has been extended to fainter magnitudes and higher redshifts by using filters at longer wavelengths, for example, F150W − F444W in Barrufet et al. (2023), or F150W − F356W in Pérez- González et al. (2023a).Recently, Labbé et al. (2023) used a threshold of F277W − F444W > 1 mag to identify candidates to massive galaxies at z > 7.Here we use a slightly redder color and we drop the additional color constraints to lower the selection redshift to z  5.A redder color threshold can also reduce the contamination by galaxies with high-EW emission lines (>1000 Å) boosting the NIR fluxes of blue galaxies with a relatively shallow stellar continuum.For example, Endsley et al. (2023) find red colors, F277W − F444W  1 and F277W − F356W  1, in a sample of low-mass, Lyman-break galaxy candidates at z = 6.5-8 which were largely driven by high-EW [O III]/Hβ lines boosting the flux in F444W.Such strong lines have also been spectroscopically confirmed by recent NIRCam/WFSS surveys at slightly lower redshifts of z > 5.3 (Matthee et al. 2023).
The left panel of Figure 1 illustrates the sample selection in a color-magnitude diagram compared to the overall distribution of galaxies in the CEERS catalog, color coded by different properties, and a subsample of F150W EROS (F150W − F444W > 2 and F444W < 25 mag; red circles).The 13 galaxies from Labbé et al. (2023) are shown with black squares.All of them except the four with colors F277W − F444W < 1.5 mag are included in our sample.The color code in the CEERS sample highlights the trend of increasing NIR colors with extinction (and similarly with stellar mass and redshift in the other panels).As discussed above, galaxies redder than the color threshold (dashed line) are candidates for massive galaxies with red, dusty, or quiescent SEDs and possibly some galaxies with high-EW emission lines.Interestingly, there are some differences between the sample of F150W EROs and F277W EROs.First, F277W EROs are fainter, with a median magnitude of F444W = -+ 25.9 1.1 0.8 mag, whereas F150W EROs span a broader range in magnitude starting at F444W  20 mag, which is consistent with the notion that by selecting in a redder band, F277W EROs lean more toward the higher-redshift tail of the massive galaxy selection.Second, F150W EROs are typically selected within a brighter limiting magnitude to restrict the number of galaxies in the lower-mass end of the selection criteria ∼ 10 (e.g., F444W  25-26; Alcalde Pampliega et al. 2019;Gómez-Guijarro et al. 2023).However, using a fainter limiting magnitude increases the overlap between the two ERO samples, as shown for example in Pérez- González et al. (2023a).Nevertheless, we find that, even within a similar magnitude range, the F150W selection misses some F277W EROs because they have bluer colors in F150W − F444W =2.2 - + 0.5 0.7 mag.The reason for this key difference is highlighted in the central panel of Figure 1, which shows that all the F277W EROs are surprisingly blue in F150W − F200W ∼ 0 mag, which probes the rest-frame UV at z > 5. Consequently, these EROs populate a very different region of the color-color diagram far from the loci of the F150W EROs, and all other massive galaxies which typically have red colors, F150W − F200W = 0.5-1.5 mag.This means that, unlike the majority of other massive galaxies, which are red across their whole SEDs, the F277W EROs are blue in the rest-frame UV and red in the rest-frame optical.Such peculiar colors indicate that these EROs have bimodal blue-red SEDs (L shaped or V shaped in f λ ), as noted by Labbé et al. (2023).The right panel shows that the goal of the second color threshold (F150W − F270W < 0.7 mag) in the selection method of Labbé et al. (2023) is to remove galaxies at z < 7 from the sample.The F150W − F270W color acts as a pseudoredshift because the F277W filter shifts from the steep optical side of the SED to the flat UV with increasing redshift.Consequently, the color quickly declines toward F150W − F270W ∼ 0 for galaxies at z  7.For the same reason, the primary selection in F277W − F444W might start missing galaxies of this type at z  9 when the F444W filter starts to shift out the steep restframe optical range.Lastly, we note that the selection in F277W − F444W > 1.5 mag is surprisingly clean as it only  Labbé et al. (2023).The left and central panels show the general trends toward redder colors with increasing mass and dust attenuation (arrows), which suggest that F277W EROs are massive and dusty galaxies.However, the central panel reveals that F277W EROs have surprisingly blue colors at SW, F150W − F200W ∼ 0 mag, very different from those of F150W EROs and massive dusty galaxies in general.The red square shows a massive, dusty, submillimeter galaxy at z = 5.1 from Zavala et al. (2023) which is also red in all bands.This implies that F277W EROs have bimodal SEDs with blue SW colors and red LW colors.The right panel shows the correlation between photometric redshift and F150W − F277W color for the F277W EROs.As the F277W filter shifts from the steep, rest-frame optical range to the flat rest-frame UV range with increasing redshift, the color declines to F150W − F277W ∼ 0 mag.identifies these peculiar EROs with bimodal, blue-red SEDs with no contamination from typical EROs (i.e., red across their whole SEDs).
We identify 37 EROs with the color criterion described above.We visually inspect all the candidates and we remove some unreliable detections (e.g., hot pixels or fake objects extracted near the diffraction spikes of bright stars).Their average magnitudes in F444W, F536W, F277W, and F150W are -+ 25.9 1.1 0.8 , -+ 26.8 1.2 0.9 , -+ 27.6 1.4 0.9 , and -+ 28.2 1.3 1.0 mag, respectively, which are consistent with the color selection criterion.Their very faint magnitudes in F150W imply that these objects are all HST/WFC3 dropouts at the depth of the CANDELS data in the CEERS region.

Brown Dwarf Contamination
The bimodal SEDs of the EROs are very different from the SEDs of typical massive, dusty galaxies at any redshift.However, they do exhibit some similarities with the 1-5 μm SEDs of cool, brown dwarfs in the Milky Way (e.g., Wilkins et al. 2014).Indeed, the SW NIRCam colors of brown dwarfs are quite blue and their fluxes drop abruptly in F090W, which can be misinterpreted as a Lyman break in a high-z (z  7) galaxy.At the same time, brown dwarfs also exhibit an upturn in their SEDs starting around 3.0 μm and peaking at ∼4.5 μm, which leads to red LW colors.Brown dwarf candidates have already been identified photometrically due to their peculiar blue-red SEDs using color-color thresholds similar to the ERO selection criteria (e.g., Hainline et al. 2023;Holwerda et al. 2023;Nonino et al. 2023;Wang et al. 2023), and recent NIRSpec observations have confirmed the stellar nature of handful of them (Langeroodi et al. 2023;Burgasser et al. 2024).
To investigate the likelihood of brown dwarf contamination in our sample of EROs we study the overlap between the two populations in color-color space.Figure 2, shows the F115W -F200W versus F277W -F444W colors for the bulk of the CEERS sample (gray scale) and the EROs (colored circles) identified in the previous section (central panel of Figure 1).Simultaneously, the black lines depict the colorcolor tracks for brown dwarfs within a small range of temperatures (T = 500-1500 K) and metallicities (log(Z/Z e ) = −1 and 0), computed using LOWZ stellar atmosphere templates (Meisner et al. 2021).The figure highlights again the dual blue-red nature of the EROs relative to the bulk of the galaxies, although, over this longer color baseline, some of the EROs have slightly redder colors in rest-frame UV (average F115W -F200W ∼ 0.25 versus F150W − F200W ∼ 0).Interestingly, all the brown dwarfs with red F277W -F444W > 1.5 colors are much bluer than any of the EROs (i.e., they have blue slopes rather than relatively flat).Based on this distribution we expect the contamination to be essentially nonexistent for EROs with F115W − F200W < −0.5.We identify only two potential brown dwarf contaminants in our sample of 37 EROs based on their much bluer F115W -F200W colors.We flag these objects in Table 1.
In addition, we use the brown dwarf stellar templates to study their typical colors at MIRI wavelengths.Since their SEDs have a maximum at around 4.5 μm, the LW colors quickly turn blue relative to F444W.On average, we find blue MIRI F444W -F560W = −0.8 and F444W -F777W = −1.5 colors, which contrast with the red NIRCam colors F356W -F444W = 0.5.As described in the following section, the EROs appear to have red colors in the MIRI bands, continuing the steep SED trend of the NIRCAM bands.Consequently, the MIRI colors provide additional leverage to distinguish galaxies from brown dwarfs.

MIRI Detection and NIRSpec Spectroscopy of the EROs
We search for counterparts of the 37 EROs in the CEERS MIRI and NIRSpec observations.Unfortunately, the MIRI coverage of the CEERS/NIRCam mosaic is quite limited (less than ∼8% of the area) and none of the pointings have simultaneous observations in the SW and LW bands.Overall, only four of the MIRI pointings in F560W and F770W and two of the pointings observed in F1000W and onward overlap with the NIRCam coverage.Surprisingly, we find clear detections for all four of the 37 EROs that lie within the MIRI-observed area.Three of them are detected in F560W and F770W with an average magnitude of 25.3 0.1 0.2 mag, respectively, and one of them is weakly but clearly detected in F1000W at 24.6 mag.While it is difficult to extrapolate from such a small sample, the high recovery fraction of observed objects, as well as the very red, power-law-like slope of the SED in the LW NIRcam bands, suggests that follow-up observations of similarly selected EROs in other fields with denser MIRI coverage is likely to yield a significant number of detections.Given the median magnitude of these objects, F444W = 25.9 - + 1.1 0.8 mag, we would expect detections in F560W and F777W in the 25-26 mag range, which is clearly within the 5σ limit for surveys similar to CEERS, like PRIMER or COSMOS-Web.Note also that above z > 7, the Hα emission line shifts into the F560W filter (see discussion in Section 4.3),    which might further enhance the flux and facilitate the detection.
In addition to the MIRI detections, four other EROs have been observed as part of the CEERS NIRSpec survey.All of them have clear emission lines that provide a robust estimate of their redshifts.Two of them have already been presented and discussed in Fujimoto et al. (2023) and Kocevski et al. (2023), nircam3-2232 and nircam3-3210, respectively.The two galaxies at z > 7 exhibit only Hβ and [O III] detections, while the other two at z  6 show Hα as well.The galaxy discussed in Kocevski et al. (2023), at z = 5.62, is the only one that has a continuum detection and exhibits a broad-line Hα emission, which confirms that it is an AGN.All galaxies have relatively low [O III]/Hβ ratios, however, as noted Kocevski et al. (2023), the narrow emission-line ratios are very similar to those of star-forming galaxies (SFGs) observed at similar redshifts, which means that the line-ratio AGN diagnostic might not be particularly effective at z  5.

Photometric Redshifts and Stellar Masses
The left panel of Figure 3 shows the overall distribution of the F277W EROs in photometric redshift and stellar mass compared to the bulk of the CEERS sample (green density map) and the sample of F150W EROs from Figure 1 (red).Overall, the F277W EROs are relatively massive and dusty with median values of mag, similar to those reported in Labbé et al. (2023) for the z > 7 population.The redshift distribution ranges between 5 < z  9 with a median of = -+ z 6.9 1.7 1.0 .This indicates that nearly half the sample is at redshifts 5 < z < 7, as suggested by Pérez- González et al. (2023a).We remove a single object at z  5 for homogeneity, but the color selection is, overall, very effective at identifying galaxies at z > 5.As expected from the color selection, the EROs tend to be among the most massive galaxies at their redshift (i.e., relative to the green map).Compared to the other F150W EROs at lower redshift, the F277W EROs tend to follow the expected decline in the number of very massive galaxies, .5, with redshift.However, we find a handful of galaxies with large masses, .5, even at z > 7, which, if confirmed, would be hard to reconcile with the observed stellar mass functions (SMFs) as well as models of galaxy evolution (e.g., see the discussion in Boylan-Kolchin 2023).Not surprisingly, these galaxies are also among the brightest in F444W by nearly 1 or 2 mag relative to the median of the sample.The reliability of the stellar mass estimates is indeed one of the fundamental questions about these EROs with unusual SEDs.The values discussed in this section are computed with FAST based on typical modeling assumptions (see Section 2.3), which work well for most galaxies at low to mid redshifts.However, this method might have limitations for these EROs (e.g., because of strong emission lines or extreme obscurations).In Section 5 we analyze in detail the impact of using different codes and modeling assumptions on the inferred stellar masses.

Sizes and Morphologies
The right panel of Figure 3 shows the distribution of the F277W EROs in a stellar mass versus size diagram compared to F150W EROs at z > 3, and the overall distribution of galaxies in the CANDELS F160W catalog in the overlapping area with CEERS (green density map).The CANDELS measurements are derived from Stefanon et al. (2015) and van der Wel et al. (2014).Sizes are represented by the effective radius, r e , of the Sérsic (Schmidt 1968) profile fit performed with GALFIT v3.0.5 (Peng et al. 2002) in the F356W band.The code was run on the background-subtracted images with sizes 2.5 times the Kron radius.An array, which includes background sky, Poisson, and read noise, was used as the input noise map.Empirical PSFs were constructed using stars in all CEERS pointings.All galaxies in the image cutout within M M log ∼10, and dusty, A V ∼ 3 mag and they span the redshift range 5 < z < 9. Overall, F277W EROs are among the most massive galaxies at their redshift, but less massive than the F150W EROs at lower redshift, following the expected decline in the number of very massive galaxies with redshift.However, a few of them are much more massive (   M M log  10.5), suggesting that there might be limitations in the fitting of their bimodal SEDs or perhaps that their continua are not stellar, but AGN dominated (see the discussion in Section 5).Right: stellar mass vs. F356W effective radius for the same galaxies.For the Labbé et al. (2023) galaxies we use our effective radii estimates.The blue and red lines show the mass-size relations for star-forming and quiescent galaxies from van der Wel et al. (2014).The dashed lines indicate the approximate resolution limit from the half-width at half-maximum (HWHM) of the PSF in F356W (FHWM = 0 15) at z = 5 and z = 9.Remarkably, all the F277W EROs appear to be unresolved point-like sources in contrast with the typical spread of F150W EROs and other massive galaxies.We find similar results in the other NIRCam bands, suggesting that the EROs are unresolved at all wavelengths.The panels on the right show the best fits to a PSF in F444W for the four galaxies with MIRI detections, which show negligible residuals.
3 mag of the primary source were fit simultaneously.All other sources were masked out during the fitting.The fitting parameters were allowed to vary within the following reasonable bounds: Sérsic index (0.2 n 8.0), effective radius (0.3 r e 400 pixels), axis ratio (0.01 q 1), magnitude (±3 mag from the initial value), and position (±3 pixels from the initial value).
Overall, we find that while the F150W EROs tend to overlap with the bulk of the galaxy sample, scattered in between the expected mass-size relations for star-forming and quiescent galaxies (blue and red lines from van der Wel et al. 2014 at z = 3), all the F277W EROs are extremely small, systematically under the resolution limits regardless of their stellar masses.The best-fit GALFIT r e ∼ 0 009 (0.3 pixels) returns in most cases the absolute lower limit set for the fitting, suggesting that the galaxies are not resolved.The dashed lines indicate the approximate minimum sizes measurable as the HWHM of the PSF (0 07) at z = 5 and z = 9, roughly r e ∼ 0.3-0.4kpc.We further explore the size measurements of the EROs in F200W, F277W, and F444W, obtaining similar results which suggest that they are unresolved in all the observed wavelengths.Note that the EROs are typically very faint (∼27-28 mag) in all the SW NIRcam bands and, in most cases, they have only a handful of bright pixels for the fitting.Lastly, we also fit the profiles of the eight EROs with MIRI and NIRSpec detections using point-like PSFs and we find excellent agreement with negligible residuals (right panel of Figure 3), indicating that they are indeed unresolved.

Overall SEDs and Possible Modeling Scenarios
The right panel of Figure 4 shows the stacked SED of all the EROs normalized to the median of the relatively flat rest-frame UV continuum traced by F115W, F150W, and F200W, divided into two groups at redshifts below and above z = 7 with purple and red markers, respectively.Both groups exhibit the distinctive, bimodal SEDs discussed in Section 3, which consist of extremely red colors at λ > 2 μm, with a relatively constant power-law slope ∼3.5 ± 0.5 μJy μm −1 , and a flat SED at shorter wavelengths.The red, power-law-like emission is typically associated with large amounts of dust attenuation.However, as discussed in Section 3.1, it is also possible that the flux in some of the LW filters is partially boosted by strong emission lines, making the colors redder than the underlying stellar continuum.The right panel of Figure 4 highlights the location of some of the strongest lines that can boost the emission in different filters as a function of redshift.At 5 < z < 7, the Hα and [O III] lines can contaminate the F444W and F356W filters while F277W probes the continuum redward of the 4000 Å break.At z > 7, the same lines shift into F444W and MIRI/F560W while F356W probes the red continuum.The average, stacked fluxes in F277W and F356W for the lowand high-redshift groups are both clearly above the flat continuum in the rest-frame UV, suggesting that there is at least some continuum emission redward of 4000 Å.Furthermore, it would be difficult to reproduce a constant power-law slope spanning both the NIRCam and MIRI bands with relatively normal, low-EW (∼100 Å) emission lines since typically at least one, but probably several bands, should not be affected by the most prominent emission lines.
Nonetheless, the very pronounced change in the slope from the blue to the red spectral regions is also difficult to model in terms of a single stellar continuum.Indeed, the best-fit templates from EAZYpy at z = 5.5 and z = 7.5 shown in the left panel of Figure 4 are often composites of two templates with very different stellar ages, masses, and dust attenuations: on the one hand, a young, low-mass, low-attenuation galaxy (i.e., a typical Lyman-break galaxy) and, on the other, a more massive and dusty galaxy.As a consequence, the inferred stellar mass and extinction of the composite is usually quite large, because it is dominated by the larger mass-to-light ratio of the older galaxy.
Recently, other works (e.g., Endsley et al. 2023;Furtak et al. 2023) discussed the possibility that the SEDs of some of these EROs could be explained partially, or completely, by very strong, AGN-driven emission lines.The presence of high-EW (>1000 Å) emission lines can boost the flux in all the filters since these are not restricted to just the brightest emission lines due to star formation.Similarly, the peculiar SEDs can also be explained in terms of continuum emission from an AGN which outshines the galaxy host in different spectral regions.This possibility was recently explored in Kocevski et al. (2023) for one of the EROs at z = 5.62 with NIRSpec observations, which is also included in our sample (nircam3-3210).This galaxy was also discussed in Labbé et al. (2023) but the estimated photo-z was much higher, z ∼ 8.This highlights again the potential pitfalls in the SED modeling of these galaxies.Kocevski et al. (2023) proposed some AGN-dominated scenarios where the SED could be explained by: (1) a heavily obscured quasi-stellar object (QSO) dominating the LW fluxes and a small percentage of scattered light from the broad-line component causing the blue, SW emission (e.g., as in the Polletta et al. 2006 torus template); (2) a heavily obscured QSO dominating the LW fluxes plus a blue, low-mass galaxy host, which dominates the SW fluxes; or (3) a blue, type-1 QSO dominating the SW fluxes in a dusty starburst galaxy, which in turn dominates the LW emission.The latter is also similar to the red QSO scenario in Fujimoto et al. (2022).
Crucially, many of these different scenarios can be confirmed or ruled out with additional observations such as the NIRSpec spectroscopy in Kocevski et al. (2023) or with additional photometry at longer wavelengths from JWST/ MIRI.For example, Papovich et al. (2023) and Rinaldi et al. (2023), have recently shown that many of the blue, low-mass Lyman-break galaxies at z > 7 with emission-line-driven excesses in F444W have clear detections in MIRI at F560W and F777W that can trace the continuum in a spectral region without prominent emission lines.For these EROs, MIRI detections in the rest-frame optical continuum can distinguish between scenarios where the red optical colors are primarily driven by high-EW emission line versus any kind of continuum-dominated emission by a red, dusty galaxy or a QSO.In the Section 5, we study the likelihood and implications of the different scenarios outlined above from a detailed analysis of the SED modeling of the four galaxies with additional photometric constraints from MIRI and the four galaxies with spectroscopic redshifts from NIRSpec.In Section 6 we use those results to inform the discussion on what would be the most likely scenario for the whole population of EROs.

Modeling Codes
In this section, we perform more detailed SED modeling of the eight EROs with MIRI and NIRSpec observations using the SEDs derived from the circular aperture photometry described in Section 2.2 and a variety of codes aimed at exploring the likelihood of the different dusty galaxy versus obscured-AGN scenarios outlined in the previous section.A detailed description of the modeling assumptions adopted for each code is provided in Appendix A. Briefly, we use EAZYpy (Brammer et al. 2008), Synthesizer (Pérez- González et al. 2008a), Prospector (Johnson et al. 2021), and a custom code to perform a hybrid fit of the stellar population models from Prospector with the AGN templates of Polletta et al. (2006).The EAZYpy fits are based on the same default template set used in Section 2.3.The Synthesizer run uses parametric SFHs, following a delayed-τ function characterized with the Bruzual & Charlot (2003) stellar population models, a Calzetti et al. (2000) attenuation law, and nebular emission following Ferland et al. (1998).With Prospector, we use three different options: (1) a fiducial model with a parametric delayed-τ SFH and Calzetti et al. (2000) attenuation law; (2) a nonparametric SFH based on the continuity priors of Prospector-α (e.g., Leja et al. 2019or Tacchella et al. 2022) but with a maximum age of 100 Myr and using a Calzetti et al. (2000) attenuation law; and (3) a similar nonparametric SFH with a more complex dust attenuation model based on Charlot & Fall (2000) and Kriek & Conroy (2013).All three options are based on FSPS models (Conroy et al. 2009) and include nebular emission from young stars.They also have a number of other modeling assumptions in common (gas and stellar metallicity, ionization parameter, etc.) described in the Appendix.The first two options are aimed at exploring the impact of using parametric/nonparametric SFHs and different stellar population models with respect to Synthesizer, while the third focuses on the impact of the dust attenuation law.The last SED model is a hybrid of a galaxy and a dustobscured QSO.Here we assume that the emission in the LW NIRCam and MIRI bands is largely dominated by an obscured QSO modeled after the QSO2 from Polletta et al. (2006), while the flux in the SW fluxes comes from the galaxy host.We also show the fits to an intrinsically blue QSO template, QSO1 from Polletta et al. (2006), with a large A V = 3-4 based again on a Calzetti attenuation law.While this template fits worse than the QSO2 one, it is useful to illustrate the differences and it provides a way to estimate the bolometric luminosity of the QSO from the unobscured emission.We fit the QSO model in three steps.First, we do a coarse fit of the QSO2 template to the LW fluxes, then we fit all the photometry subtracting the bestfit fluxes from the QSO template with Prospector delayed-τ models, and lastly, we perform a simultaneous fit with the QSO template to galaxy SEDs drawn from the posterior of the Prospector fit.The results from this method are similar to those obtained with the modified version of FAST (Aird et al. 2018) used in Kocevski et al. (2023).The advantage of the Prospector fit is that it includes emission lines that can help shore up the limitations of the obscured QSO template, which has a fixed set of emission lines.While this is not a fully selfconsistent AGN method, it helps to account for the contribution of emission lines to the photometry.

Photometric Redshifts of the Four EROs with MIRI Detections
The peculiar SEDs of the EROs and the high chances that some of the fluxes are at least partially boosted by emission lines make the photometric redshift estimates one of the key parameters and potentially one of the most problematic.For that reason, we run EAZYpy twice, first using the default modeling assumptions and a second time using the recently updated models which include a blue galaxy template with strong, high-EW emission lines similar to those observed in recent NIRSpec spectra of z > 7 galaxies.We also include in the analysis the redshift probability distributions (PDFz), shown in Figure 5, the values computed in Finkelstein et al. (2023a) using the original version of EAZY with an updated template set optimized for high redshift presented in Larson et al. (2023).The latter fits do not include the MIRI fluxes and thus allow us to gauge the impact of the additional photometry in the redshift likelihood.Lastly, we also include the PDFz estimate from the fiducial Prospector fit described in the previous section.
nircam5-5815.The primary EAZYpy and Prospector solutions agree on a value of z ∼ 5 for which strong Hα emission would boost the fluxes in F410M and F444W.There is a secondary solution at z ∼ 9 for which the red F277W − F444W color is caused by a strong Balmer break.However, at that redshift, the galaxy should be an F150W dropout, and the galaxy is clearly detected at >5σ.Therefore, we adopt the lower redshift solution as the primary.
nircam5-9553.The photometric redshift distributions from EAZYpy and prospector are quite consistent, peaking around z ∼ 5.8.At this redshift, the [O III]/Hβ and Hα lines can contribute to the flux in F356W and F444W but not in F410M (or at least not significantly).There is a secondary peak at z = 8.7 which also produces a relatively good fit.However, like in the previous galaxy, this would require the F150W flux to be a dropout, and the galaxy is faint but clearly detected in that band.Therefore we consider the low-redshift solution as the primary.
nircam6-7042.This is the only galaxy observed in the LW MIRI bands.It has a faint but clear detection in F1000W but is not detected in F1500W.Similarly to the galaxies above, the PDFz exhibits a primary peak at z = 6.4 and a secondary peak at z ∼ 8.5, which is closer to the value presented in Labbé et al. (2023), z = 8.11.The two different solutions try to fit an excess in F444W relative to F410M with a strong emission line, either Hα or [O III] at low and high z respectively.We notice however that the F277W flux for this source is above the relatively flat continuum delineated by the SW bands, suggesting that it might be sampling the continuum redward of the 4000 Å break and therefore favoring the low-z solution.The F277W photometry in Labbé et al. (2023) appears to be fainter and closer to the bluer bands, which might favor the high-z solution.At the redshift of the two possible solutions, the F1000W detection (and the upper limit in F1500W) still probes rest-frame wavelengths shorter than the 1.6 μm bump and thus cannot help discriminate between them.
nircam5-6746.This galaxy presents a PDFz centered around z = 7-8 with no secondary peaks at significantly different redshifts.The brighter MIRI flux in F560W relative to F770W also favors a redshift of z = 7.5, suggesting that strong Hα emission is boosting the flux in F560W and similarly [O III] in F410M and F444W.This galaxy is also discussed in Akins et al. (2023) with a similar photometric redshift and consistent stellar population fits.

Best-fit Properties and SEDs
Figures 6 and 7 show the multiband images, NIRSpec spectra, and SEDs for the eight MIRI-and NIRSpec-detected galaxies jointly with the best-fit models obtained with the different codes outlined in the previous section.From left to right, the panels show the stellar population fits with Prospector (τ-model and nonparametric) and Synthesizer, the composite stellar populations with EAZYpy (middle), and the hybrid galaxy + AGN models (right).
MIRI fluxes and the high-EW emission-line scenario.The four galaxies with MIRI detections exhibit F560W and F777W fluxes that continue the red power-law trend outlined by the NIRCam LW bands.For three of them, the MIRI bands probe a spectral region redward of Hα, which does not have any prominent emission lines.The exception is nc5-6746 at z = 7.5, which seems to have an excess in F560W due to a strong Hα line, but not in F777W, which also continues the same trend of increasingly larger fluxes as the other three galaxies.Therefore, the MIRI detections strongly suggest the presence of red continuum emission in these galaxies, which disfavors the scenario where the red optical fluxes originate in a blue galaxy with very high-EW emission lines masquerading as a red continuum.Nonetheless, we note that the best-fit SEDs for these EROs show strong emission lines and even emissionline-driven excess in one or two of the LW NIRCam bands.However, these lines have relatively normal EWs for a massive SFG (∼100 Å) due to the presence of a red stellar continuum.
Prospector-τ, -np, and Synthesizer.Overall, these models based on different SFHs but using the same Calzetti et al. (2000) attenuation law provide a relatively good fit to the Figure 5. Photometric redshift distributions (PDFz) for the four MIRI-detected EROs computed using EAZY, EAZYpy, and Prospector.The PDFzs derived with the default and blue versions of the EAZYpy templates agree well with one another and with the Prospector estimates for all the galaxies.For the three galaxies at z < 7, the PDFzs based on the templates with very high-EW lines (blue) suggest a secondary peak at higher redshift that is not supported by the detections in F150W.The key difference between the lowand high-z peaks is typically an emission-line-driven excess in F444W which could be attributed to Hα or [O III], respectively (see also Figure 4).majority of the LW NIRCam bands and the MIRI fluxes.However, they all fail to reproduce the rest-frame UV fluxes probed by F115W, F150W, and, in some cases, F200W, regardless of the SFH.Both Prospector fits yield systematically lower fluxes in the rest-frame UV, while Synthesizer sometimes finds a trade-off between improving the fit to the UV bands at the expense of a worse fit to the optical bands.The reason behind this systematic issue for all the models is that the large dust attenuations required to reproduce the extremely red optical colors lead to even larger attenuations in the UV which completely suppress the predicted emission regardless of the stellar population parameters or SFHs; i.e., even nonparametric SFHs having substantial star formation rates (SFRs) in the last 5-10 Myr still yield very red colors in the rest-frame UV.This problem is unavoidable for the typical attenuation laws such as Calzetti (A 2500 /A V ∼ 2), and it would  (2006).The left panels illustrate that fits based on a single stellar population component provide a good fit to the overall LW NIRCam and MIRI photometry (black and green squares) but they systematically fail to reproduce the rest-frame UV probed by the SW NIRCam bands.The middle panels show that a composite model consisting of two (or more) stellar populations provides an excellent fit to all the bands by combining a red, massive, and dusty galaxy that fits the LW bands and a blue, low-mass galaxy that fits the SW bands but has little impact on the stellar mass.The right panels show that the hybrid galaxy + QSO model (QSO1 and QSO2, orange and red, respectively) provides an equally good (or better) fit to the SED than the other models.Here, a dust-obscured QSO dominates the LW photometry but does not contribute to the stellar mass of a blue unobscured host, and consequently leads to total stellar masses ∼two orders of magnitude smaller than the other scenarios.The two stellar templates (gray) illustrate the 16%-84% confidence range in stellar mass for the galaxy component.The SEDs exhibit similar UV emission but increasingly larger optical emission with mass.
be worse for steeper attenuation laws such as an SMC type (A 2500 /A V ∼ 2.6) or a Milky Way type with a UV bump at 2175 Å.However, a shallower, grayer attenuation law, resulting perhaps from a more patchy distribution of the dust in the galaxy, could alleviate this problem.
Prospector-np-cf.Indeed, the best-fit SED models derived with Prospector using nonparametric SFHs and a more complex, two-component dust attenuation model based on Charlot & Fall (2000) and Kriek & Conroy (2013) provide a better match to the UV fluxes with varying degrees of improvement.In this model, the diffuse attenuation is multiplied by a power law with index n that increases/lowers the slope of the attenuation law relative to Calzetti (i.e., for n = 0 it becomes Calzetti).The models that fit the UV fluxes best (e.g., nc5-5815, nc6-7042, or nc1-2441) all have similar attenuation laws which lean heavily toward the shallowest (grayest) possible attenuation law allowed by the priors (n = 0.4 and A 2500 /A V ∼ 1.4); i.e., the posterior is not evenly sampled but rather skewed to the maximum value.The models without a significant improvement of the UV fit still return a better χ 2 than the Calzetti-based fits.For these galaxies differential attenuation between the stellar continuum and the emission lines introduced by the two-component Charlot & Fall (2000) prescription appears to allow stronger emission lines that improve the fit to the bands with emission-line excesses.EAZYpy.These models provide a good match to both the rest-frame UV and optical SEDs.The difference with respect to the Prospector and Synthesizer fits is that EAZYpy uses composite models that are linear combinations of templates with different ages, SFRs, and, crucially, dust attenuations.Consequently, the composite SED is not necessarily bounded by the same dust attenuation across the whole spectral range.The best-fit models for all the EROs are always a combination of at least two templates with very different properties: a young, blue galaxy with low dust attenuation that fits the relatively flat rest-frame UV emission and an older galaxy with large dust attenuation that fits the red optical emission.
Hybrid galaxy + red QSO.The rightmost panels of Figures 6  and 7 show the fits to the hybrid model of a blue galaxy and a dust-obscured QSO (QSO2 template in red).This model shows an excellent fit to the overall SED including the rest-frame UV and the MIRI fluxes.In this scenario, the continuum emission from the obscured QSO dominates the SED redward of F277W while the galaxy component dominates the rest-frame UV emission.Consequently, the best-fit galaxy model is a blue, low-extinction galaxy similar to the blue component in the EAZYpy composite.The gray and magenta lines in the fits illustrate the 1σ range in the stellar masses which are, in all cases, very small,   M M log = 7-8.The main difference in the best-fit SEDs of QSO-dominated versus galaxy-dominated scenarios is that in the latter, the stellar continuum typically exhibits a peak around ∼1.6 μm, whereas the QSO emission increases continuously toward the rest-frame mid-infrared.Unfortunately, at z > 5 the MIRI detections in F560W and F777W still probe rest-frame wavelengths shorter than 1.6 μm, and even for the one galaxy detected in F1000W, the rest-frame flux is still too close to 1.6 μm.Detections at longer wavelengths are clearly necessary to distinguish conclusively between a declining stellar continuum and rising QSO emission.The panels show also the fits using the blue QSO1 template with very large attenuations (A V  3, orange).These are generally a worse fit to the MIRI data because they have more steeply rising SEDs, but they help provide an order-ofmagnitude estimate of the QSO bolometric luminosity.
Hybrid blue QSO + dusty galaxy and pure QSO + torus.Figure 8 shows the best-fit SED of nc3-3210 (the broad-lined AGN with NIRSpec studied in Kocevski et al. 2023) to the other two possible scenarios involving a QSO: QSO-torus emission and a hybrid model consisting of a blue QSO and a red, dusty galaxy.In the torus model, the SED is completely QSO dominated at all wavelengths (e.g., scattered UV light, attenuated optical emission, and mid-to far-infrared reemission by dust).Here we use the torus template from Polletta et al. (2006) to fit the observed SED and we find that while the intrinsic shape of the torus SED template is, to some extent, similar to the bimodal SED of the EROs, a single template is not flexible enough to obtain a better fit than any of the other scenarios discussed above.This is likely a limitation of our approach based on a single template, and it is possible that a more comprehensive AGN modeling code can fully reproduce the observed SED with higher accuracy.The scenario involving a galaxy plus a blue QSO is, to some extent, similar to the EAZYpy model.In both of them, the LW NIRCam bands are largely dominated by the emission of a red, dusty galaxy while the SW bands are dominated by a blue, low-extinction galaxy or QSO.Consequently, the inferred stellar masses and dust attenuations for the bulk of the galaxy are also very similar, since none of the blue components contribute significantly to the mass.These two scenarios are not discussed in detail for the other objects because the blue QSO model leads to similar results for the stellar properties as the other galaxy-dominated scenarios, and the torus model does not provide constraints on the stellar mass of the host or the luminosity of the QSO.

Stellar Masses and Attenuations
Figure 9 shows the ranges of stellar masses and dust attenuations for the eight EROs obtained with the different SED modeling codes.We also include the stellar masses and attenuations computed with FAST and use it as a benchmark Figure 8.Additional SED modeling scenarios involving a QSO.Left: a hybrid of a dusty-galaxy-dominated SED with a blue, low-extinction QSO contributing only to the rest-frame UV emission.This scenario is similar to the EAZYpy fits replacing the blue galaxy with a blue QSO with a minimal impact on the stellar mass of the composite.Right: a pure QSO-dominated model based on the torus template from Polletta et al. (2006) where the emission from the QSO outshines the galaxy host at all wavelengths.The intrinsic shape of the torus SED is very similar to the bimodal SED of the EROs.However, we find that using a single template limits the flexibility of the fits and it leads to generally worse agreement (χ 2 ) with the data.
model for the comparisons to study systematic effects.FAST has been widely tested in typical galaxies at low to mid redshift with accurate results, but it is critical to understand if there are potential issues modeling these high-z galaxies with peculiar SEDs.
The stellar masses computed with FAST and EAZYpy tend to be the largest, and they are very similar, with a median difference and scatter of Δ Although the SED fits with FAST do not reproduce the UV fluxes like the composite SEDs with EAZYpy, the effect on the stellar mass is very minor.This is because the red, dusty component in the EAZYpy fit, which is similar to the overall FAST fit, dominates the stellar mass over the young, blue component, which has a much lower mass-to-light ratio.
Interestingly, the median difference with respect to the stellar masses computed with the fiducial Prospector fits (τ-model with Calzetti attenuation) is relatively small, with a larger scatter Δ   M M log (FAST -Prospector-τ) = −0.16± 0.49 dex.This means that despite the more flexible modeling of key parameters like emission-line strength or metallicity, the stellar mass is mostly driven by the need to fit the red optical slope with high dust attenuation.In fact, the cases where the Prospector fits obtain the largest stellar masses are typically those where the extinction is a maximum A V ∼ 4. Note also that the extinction values from FAST and Prospector are typically the largest, ranging between A V = 3-4.The EAZYpy fits have lower extinctions A V = 2-3 in part because of the combination with a blue template (A V = 0), but sometimes because it includes a red, quiescent template that also has a low attenuation but a large mass-to-light ratio, which, in turn, leads to larger stellar masses (e.g., as in nc3-2232).
The Prospector fits with nonparametric SFHs capped at a maximum formation age of 100 Myr and a Calzetti attenuation law leads to systematically lower stellar masses than the fiducial Prospector-τ, with Δ   M M log (τnp) = −0.23 ± 0.21 dex.This is because the fiducial model has a maximally old start of the SFH at 90% of the age of the Universe at the redshift of the galaxy and, consequently, tends to form more stars over a longer period of time.Consequently, the masses are even smaller relative to FAST Δ The Prospector fits with nonparametric SFHs and a more flexible attenuation law based on Charlot & Fall (2000), which provides the best SED fits, exhibit an interesting behavior.For the four galaxies with MIRI detections, the stellar masses are significantly lower with Δ   M M log (FAST -Prospectornp-cf) = −0.68 ± 0.28 dex, but for the four galaxies with NIRSpec data the difference is nearly zero, Δ   M M log (FAST -Prospector-np) = −0.01 ± 0.16 dex.The reason for this difference is clearly visible in the SED fits shown in Figures 6 and 7. Without MIRI data to constrain the continuum emission beyond F444W, Prospector favors solutions with a stronger continuum (i.e., more massive) and lower EWs for the lines.For example, in nc3-3210 or nc1-9410, the best-fit models with Prospector-τ versus Prospector-np-cf would exhibit differences in the predicted F560W and F770W fluxes of the order of 1-1.5 mag.
The fits with Synthesizer provide the smallest stellar mass estimates, nearly 1 dex smaller than FAST, Δ   M M log (FAST -Synthesizer) = −0.98 ± 0.33 dex.As discussed in the previous section, these SED fits are, overall, less accurate than the other codes, but tend to fit the UV region a bit better at the expense of a worse fit to the optical.As a result, they have lower attenuations of A V ∼ 2 mag and, consequently, lower stellar masses.
Lastly, in the hybrid galaxy plus obscured QSO scenario, the latter completely dominates the bulk of the emission in the LW bands.However, it does not contribute to the stellar mass, which depends exclusively on the faint blue galaxy host.Consequently, inferred stellar masses are ∼two orders of magnitude,   M M log = 7-8, smaller than in any of the scenarios where the bright LW continuum originates in a dusty massive galaxy.
In summary, the commonly used methods based on τ-models and Calzetti attenuation, or variants of EAZY with the default templates (including the reddest dusty/old templates) are likely to obtain the largest stellar masses.Nonparametric or similar SFHs that limit the age of the galaxy to relatively young values (100 Myr) lead to lead to lower stellar masses by 0.4 dex.The addition of a more flexible dust modeling to allow grayer attenuation curves can lead to stellar masses up to 0.7 dex smaller.However, without MIRI data, the stellar masses can also be as high as for the fiducial τ-models.6. Discussion

Likelihood of the Dusty-galaxy Scenario
In the previous sections, we discussed three possible scenarios where a dusty, SFG can fit the overall SEDs of the EROs dominating the emission in the rest-frame optical: with a flat, gray attenuation law or with a secondary component which is either a blue, low-extinction galaxy or a blue QSO, that fits the rest-frame UV.
Looking at these possibilities in the light of the point-like, unresolved sizes for all these galaxies, the scenario with two distinct stellar components seems quite unlikely.Such a model would make more sense for an extended galaxy with clearly differentiated regions (e.g., clumps, or a bulge).On the other hand, a compact size might help explain the very gray attenuation law in terms of the geometry and distribution of dust in a high-density environment.For example, rather than a dust-shell scenario we might have a mixed star-dust distribution (probably clumpy) which produces gray attenuation laws including huge extinctions (A V  20 mag or more), but also significant scattering resulting in much lower and grayer total attenuations and, consequently, bluer UV colors (Witt & Gordon 2000).
The scenario involving a blue, low-luminosity QSO is also plausible as it can help explain why the colors of these EROs are very different from those of F150W EROs and other dusty galaxies at higher redshift recently identified with JWST (e.g., Pérez- González et al. 2023a;Zavala et al. 2023), which are red in all the NIRCam bands.As discussed in Kocevski et al. (2023; see also Fujimoto et al. 2022) this scenario could be a transitional phase in the evolution of a dust-obscured starburst that is clearing up the dust and leading the way to an unobscured QSO.Note that while bluer UV colors have been reported in dusty SFGs at z  3 with large IR luminosities (e.g., Casey et al. 2014), these EROs are very blue, with relatively flat UV continua in f ν , which imply very steep UV slopes of β  −2 for the high attenuations implied by the SED modeling, i.e., A V > 3 mag.
Taken together, the different colors and morphologies of these EROs relative to the other massive dusty galaxies might be an indication that they are a distinct population, perhaps undergoing a strong nuclear starburst phase as seen for example in some of the radio/submillimeter-detected galaxies at z > 3 (Barro et al. 2017;Tadaki et al. 2017).To some degree, this scenario might be similar to that of the compact SFGs at z  2-3 which are also small (but resolved, r e ∼ 1 kpc), massive, and dusty (e.g., Barro et al. 2013;Nelson et al. 2014;Williams et al. 2014;van Dokkum et al. 2015), and exhibit a large fraction of X-ray AGN detections (Kocevski et al. 2017).Indeed, galaxy formation models suggest that the progenitors of those compact SFGs could be even smaller at higher redshift due to the larger gas reservoirs leading to wet-compaction events that result in the formation of a very dense core (e.g., Wellons et al. 2015;Zolotov et al. 2015;Tacchella et al. 2016).
Nevertheless, it seems odd that all these EROs at z = 5-9 are unresolved.If they were to evolve into compact SFGs at z  3 we would expect some of them to be transitioning from purely unresolved to the characteristic mass-size relation that compact SFGs follow at z  2-3 (Barro et al. 2017).Furthermore, we note that if the intrinsic sizes of these galaxies are under the 200-300 pc half-light radius limit (or even 150 pc; e.g., Baggen et al. 2023) the implied stellar mass densities for the most massive EROs, over > 10, would exceed even the stellar mass densities observed even in the most massive galaxies at z = 0 (Σ M  10 11 M e kpc −2 ; Bezanson et al. 2009;Hopkins et al. 2010).

Likelihood of the Obscured-AGN Scenarios
An alternative scenario to the dusty SFGs where we expect unresolved, point-like sources and peculiar, nonstellar SEDs is in AGNs where a bright QSO can outshine the emission of its host in different spectral ranges from the UV to the midinfrared.For example, hybrid galaxy + AGN SEDs where the latter dominates the near-to mid-infrared emission are a relatively common occurrence in galaxy surveys at mid to high redshifts (Stern et al. 2005;Lacy et al. 2007;Donley et al. 2012Donley et al. , 2018)).In the previous sections, we discussed two possible scenarios where an obscured AGN can fit the overall SED of the EROs dominating the red optical emission: (1) combined with a blue, low-mass galaxy host or (2) in a pure AGN model where the emission from the QSO dominates at all wavelengths.
The first scenario would imply that all these EROs are lowmass galaxies whose optical to IR fluxes are completely outshined by the emission of an obscured QSO.The limiting factor in this scenario is the bolometric luminosity and implied black hole mass of the QSOs, which should be at least one or two orders of magnitude lower than the stellar masses of the hosts (e.g., Kormendy & Ho 2013).The stellar masses of the blue, low-extinction hosts inferred in the previous section range between . Therefore, we would expect black hole masses of the order of   M M log = 6-7 and, based on the typical luminosity-black hole mass relation (Greene & Ho 2007), QSO bolometric luminosities of L bol ∼ 10 44-45 erg s −1 or smaller, since this is the value at the higher end of the accretion rate, L bol /L Edd = 1.
Unfortunately, the estimate of the bolometric luminosity of an obscured QSO requires X-ray, UV, or bolometric luminosities, none of which can be easily computed for these galaxies.For intrinsically blue QSOs the total luminosities can be estimated from monochromatic luminosities using bolometric corrections (e.g., Richards et al. 2006).However, for obscured AGNs, the total luminosities are usually inferred from rest-frame IR luminosities of the total IR luminosity (e.g., Donley et al. 2012;Runnoe et al. 2012), which for these galaxies would require MIRI fluxes at the longest wavelengths.Thus, the only alternative to estimate a luminosity is to fit the SED with a blue QSO template heavily obscured with a Calzetti attenuation law and then transform the dust-corrected UV luminosity into L bol (e.g., L bol = 5.15 L 3000 ; Richards et al. 2006).The values obtained for the EROs with this method range between L bol ∼ 10 46-47 erg s −1 , which are 1-2 dex larger than the expectation from typical low-redshift black hole mass to stellar mass ratios (i.e., they would be very luminous QSOs).We caution however, that this estimate is a large oversimplification since, as shown in Section 5, the red QSO SED (QSO2) differs from the attenuated blue QSO SED (QSO1 + Calzetti).The true attenuation law of an obscured QSO depends on multiple factors such as the geometry and distribution of dust in the torus or the line-of-sight inclination.In compact galaxies at high z, it might even depend on the galaxy-wide conditions (gas/dust fractions; e.g., Gilli et al. 2014).Consequently, the bolometric luminosities of obscured QSOs are probably lower than the values estimated with a blue QSO template, which should be considered upper limits.
In the second scenario, the obscured QSO completely outshines the galaxy host emission across the whole spectral range; i.e., both the SW and LW NIRCam fluxes arise from the QSO.This scenario would be the most plausible based on the unresolved sizes of these galaxies in all the NIRCam bands.Unfortunately, a more detailed characterization of the bolometric luminosity in this type of scenario requires more complex modeling of the extinction and scattering of the QSO emission that is beyond the scope of this paper.Interestingly, in this scenario, the constraints on the bolometric luminosity of the QSO might be less strict since the galaxy host does not have to be detected in the UV.Therefore, a slightly more massive and dusty galaxy, , can perhaps hide under the bright red continuum of the QSO without having a significant impact on the observed SED.

Implications for Number Densities and Mass and
Luminosity Functions 6.3.1.If the EROs Are Massive, Dusty Galaxies As discussed in Labbé et al. (2023), if all of these EROs are dusty galaxies with relatively large masses ∼ 11 for some of the most extreme objects, their number densities can lead to some tension with the observed SMFs and would imply higher than expected star formation efficiencies (Boylan-Kolchin 2023).We review the number density estimates using the full sample of 37 EROs selected over the larger area of the full CEERS survey and spanning a broader redshift range from z = 5-9.Figure 10 shows the redshift evolution in the number density of galaxies with stellar masses larger than   M M log = 10 (9.5 and 10.5 in dashed lines) derived from pre-JWST SMFs in the literature (Muzzin et al. 2013;Grazian et al. 2015;Stefanon et al. 2015;Song et al. 2016;Stefanon et al. 2021).The orange lines show a similar prediction from mock catalogs based on the Santa Cruz semianalytic models (Somerville et al. 2015;Yung et al. 2019b;Somerville et al. 2021;Yung et al. 2022), which shows the median and 84th and 16th percentiles from 100 CEERSsized fields subsampled from a 2 deg 2 light cone (Yung et al. 2023) to illustrate the effect of cosmic variance.These results have been shown to agree well with observed luminosity functions and other observations in this redshift range (Yung et al. 2019a(Yung et al. , 2019b)).The purple, red, and green markers show the number densities of EROs with   M M log > 10 at redshifts z = 5-7 and z = 7-9, estimated with FAST, Prospector-np, and Synthesizer, respectively.The error bars indicate Poissonian errors.As discussed in the previous sections, these values generally bracket the largest to smallest stellar mass estimates and therefore provide a way to estimate the impact of the SED modeling choices on the number densities.
The densities of EROs with   M M log >10 at z = 5-7 are all slightly under the values from the literature, which still allow additional, non-ERO, massive galaxies to exist at this redshift (e.g., Zavala et al. 2023) without tension with the literature.At z ∼ 7, the expected number of galaxies with > 10 in the area of CEERS is roughly one (with large errors) while the densities of EROs, inferred from the different stellar mass estimates, range between 2 and 10.Nonetheless, these differences are still within the range of the uncertainties in the SED modeling and cosmic variance.Furthermore, it is possible that pre-JWST SMFs missed > 10 if they are massive dusty galaxies.The purple, green, and red colors show the predictions based on three SED modeling scenarios which typically encompass the minimum/maximum stellar mass estimates (see Section 5.4).The orange lines show a similar prediction derived from the median and percentiles of 100 CEERS-sized draws of a 2 deg 2 light cone based on the Santa Cruz semianalytic models.The scatter around the median illustrates the impact of cosmic variance in an area of the size of the CEERS footprint.Similarly, the cyan lines show the predictions from mock light cones with larger baryon conversion efficiencies (ò = 0.5 and 1).At z = 5-7 the density of EROs is lower than the values from the literature, however, at z > 7 the density can be up to a factor of ∼10 larger for some of the estimates with the largest stellar masses.While this difference can still be reconciled with the large uncertainties from the SMFs and the SED modeling variations, the discrepancy at larger masses,   M M log > 10.5, is much larger.We expect one galaxy in an area 10-20 times larger than CEERS and we find three.Right: galaxy number densities with stellar masses > 8 from the literature (gray squares) as a function of redshift.The red stars show the density of EROs if they are a hybrid of an obscured AGN and an unobscured, blue galaxy host (see Section 5.3).In this scenario, the hosts have significantly lower stellar masses than in the dusty-galaxy scenario by up to 2 dex.As a result, the EROs only make up for a small fraction (10%) of the much more abundant low-mass galaxy population.some of the most massive galaxies.At the largest masses,   M M log > 10.5, however, the difference increases up to a factor of ∼60.We should not detect any such galaxies in the area of CEERS (or even an area 10 times larger).Labbé et al. (2023) reported one of those galaxies in their sample.Here we identify five galaxies (including Labbé et al. 2023ʼs) with masses above .5 by at least two out of the three estimates, two of them at z = 5-7, and another three at z = 7-9 (nc1-10084, nc5-3637, and nc8-13596, see Table 2).These galaxies are also among the brightest in F444W ∼ 23 mag (up to 2 mag brighter than the median of all EROs), which indicates that they are different in some way and perhaps they are the ones that are AGNs.Nonetheless, an individual analysis of these sources carefully characterizing their photo-z values and masses is required to clarify the strong discrepancy with respect to the expected densities.In summary, while the number densities of EROs exhibit some tension with the predictions of pre-JWST SMFs if they are all dusty galaxies, the numbers match relatively well if all the masses are closer to , as predicted by some of the SED modeling scenarios.Nonetheless, even if just a few of them are confirmed to be very .5, the discrepancy with the SMFs would be very large.
Comparing to simulations, the cyan lines in Figure 10 show the predictions based on mock light cones presented in Yung et al. (2023), for which dark matter halos have been extracted from N-body simulations in a standard Lambda cold dark matter (ΛCDM) cosmology.Each line indicates the density of objects that would result if each halo is able to convert a different fraction of its baryon content into stars (i.e., m * = òf b M halo ), where ò = 0.5 and 1, respectively.A value of ò = 1 is expected to yield an extreme upper limit since the fractions in the local Universe are typically less than ò ∼ 0.2.Models based on similarly low efficiencies, such as the Santa Cruz semianalytic model (orange lines), yield good agreement with the density of EROs at z = 5-7.However, the implied masses of the EROs in the z = 7-9 bin, if they are primarily powered by stars, would imply significantly higher than expected values of this baryon conversion efficiency ò, although not in fundamental tension with ΛCDM.

If the EROs Are Obscured AGNs
As discussed in Section 5.4, if the bright optical emission of the EROs is dominated by an obscured AGN, but the faint UV emission arises from an unobscured galaxy host, the implied stellar masses of the EROs can be up to two orders of magnitude smaller than in the dusty-galaxy scenario.We apply the SED modeling described in Section 5.3 for the hybrid galaxy + red QSO case to the full sample of 37 EROs and we obtain similarly low stellar masses with a median and 1σ percentiles of 8.10 0.67 0.61 .In this scenario, the resulting number densities of EROs are one to two orders of magnitude lower than the typical densities for low-mass galaxies inferred from the SMFs in the literature (right panel in Figure 10).This implies that, as opposed to the dusty-galaxy scenario, the fraction of obscured AGNs in blue, low-mass hosts would be very small compared to the parent sample of galaxies with similar stellar masses.
In this scenario, the contribution of the EROs to the overall AGN population at high redshift is difficult to quantify because of their very obscured nature.Ideally, one would like to use an intrinsic property of the AGNs, such as bolometric luminosity or black hole mass; however, these measurements are not always available for AGN samples selected with different methods (e.g., X-rays or broad emission lines).Instead, comparisons are often done in the context of UV luminosity functions (at λ = 1450 Å), which allow for a direct comparison to the luminosity functions of galaxies bright, blue QSOs, which, until the advent of JWST, was the largest population of high-z AGNs.
Unfortunately, this comparison is potentially misleading for obscured AGNs, which are UV faint but intrinsically bright and massive.Figure 11 shows the UV luminosity function for the 37 EROs in two redshift bins, z = 5-7 and z = 7-9, compared to other samples of bright QSOs (Matsuoka et al. 2018), X-ray selected AGNs (Giallongo et al. 2019), broad-emission-line AGNs (Harikane et al. 2023;Kocevski et al. 2023), and galaxies (Bouwens et al. 2021) at the same redshifts.The filled red circles show the density of EROs based on their observed (i.e., obscured) UV luminosities.At face value, their density is much higher than the UV-bright QSO population but is comparable to those of the X-ray or broad-lined AGNs, which exhibit similar luminosities.We note, however, that so far these EROs have not been X-ray detected, most likely due to their high obscuration (e.g., Kocevski et al. 2023).Compared to the density of galaxies, the EROs still make up only a small fraction of the parent galaxy population, although, as noted in Harikane et al. (2023), this number is larger than in the local Universe (1-2%; Stern & Laor 2012).
If instead we compare the densities derived from the intrinsic, dust-corrected UV luminosities of the AGNs (open red circles), the ERO population overlaps in luminosity with the tail of the UV-bright QSO distribution, and they are nearly three orders of magnitude more abundant.This would imply a surprisingly large increase in the number of relatively luminous and massive AGNs.To emphasize the key distinction between the observed/intrinsic UV luminosity it is worth noting that if the EROs were instead a combination of massive, dusty galaxies and faint, blue QSOs dominating the UV emission (as in the SED fits shown in the left panel of Figure 8) we would obtain exactly the same UV luminosity function (filled red circles) but, in that case, the AGNs would be intrinsically faint with low bolometric luminosities of L bol ∼ 10 44 erg s −1 .
The median and 1σ bolometric luminosity of the 37 EROs is = -+ L 10 bol 46 0.5 0.3 erg s −1 , which implies black hole masses of the order of logM BH  8, assuming Eddington ratios of L bol /L edd ∼ 1, or larger for sub-Eddington regimes.The comparison of these black hole masses with the stellar masses of their hosts, derived with the hybrid blue galaxy + red QSO model, would imply large, unprecedented mass fractions, M BH /M galaxy = 0.81 - + 0.48 0.78 , much larger than the typical ratio of ∼0.01 seen locally (Kormendy & Ho 2013;Reines & Volonteri 2015).
In summary, if the EROs are obscured AGNs, the significant increase in the number of bright, massive BHs would mean that we are beginning to unveil a key era of very quick, heavily obscured, black hole growth with short-duty cycles occurring in the first 1 Gyr of the Universe.

Prospects for Revealing the Nature of These EROs
These sources are complex to interpret.Even though the SW photometry with MIRI and the spectroscopy with NIRSpec help place better constraints on the presence of a red continuum or the redshift of these sources, they are not enough to break the degeneracies in the possible modeling scenarios.Additional MIRI photometry at longer wavelengths can distinguish between the rising continuum of an obscured QSO and the decline in the stellar SED past 1.6 μm.Deeper NIRSpec spectroscopy of their rest-frame UV or rest-frame optical can reveal high-excitation emission lines (e.g., C II and Mg II, or He II and [Ne V]), indicative of AGNs or reach the stellar continuum showing absorption lines that would confirm the presence of an underlying stellar population.

Summary
We identify 37 EROs in the CEERS field with NIRCam colors F277W − F444W > 1.5 mag, down to a limiting magnitude of F444W < 28 mag.These are candidate massive dusty galaxies at z > 5.
1.A key defining feature of these EROs is that all of them have blue colors in the SW NIRCam bands (F150W − F277W ∼ 0).The color difference in the SW and LW bands indicates that these galaxies have bimodal SEDs consisting of a red, power-law slope (α ν > 3) in the rest-frame optical, and a blue, flat slope in the rest-frame UV.These colors and SEDs are very different from those of other EROs or massive dusty galaxies at lower or similar redshifts.2. Another key feature is that all of them are remarkably compact and featureless.The light profile fits with GALFIT indicate that they are unresolved, point-like sources in all the NIRCam bands.This differs again from the typical spread in stellar mass and size of other EROs or massive galaxies at similar redshifts.3. Their photometric redshifts, stellar masses, and dust extinctions derived with standard SED fitting codes EAZYpy and FAST range from 5 < z < 9 with median values á ñ = -+ z 6.9 1.6 1.0 , 10.2 0.4 0.5 , and = .3, respectively.However, if the red colors are not due to stellar continuum emission in a dusty galaxy, these values might be overestimated.Alternative scenarios include: emission lines with extreme EWs > 1000 Å from a galaxy or an AGN boosting the LW fluxes, a hybrid of a galaxy and a dusty QSO with the latter dominating the LW continuum, or an AGN dominating the whole SED. 4. Four of these EROs within the limited MIRI views overlap with the CEERS/NIRCam mosaic are clearly detected, showing that the extremely red colors extend to longer wavelengths.Another four EROs were observed with NIRSpec and they exhibit [O III] and Hα emission lines which confirm spectroscopic redshifts in the z = 5-9 range.The MIRI detections at rest wavelengths redward of the most prominent emission lines indicate the presence of a continuum and disfavors a scenario where these EROs are intrinsically blue galaxies with high-EW emission lines masquerading as a red continuum.5. We investigate the likelihood and implications of the different modeling scenarios using the eight MIRI-and NIRSpec-detected EROs to test a variety of codes with flexible options to characterize the stellar continuum, emission lines, dust attenuation, SFH, etc.For scenarios where the LW bands are dominated by a dusty galaxy, we find: (1) SED models based on either parametric or nonparametric SFHs and a Calzetti attenuation law fail to reproduce the blue, rest-frame UV emission regardless of the modeling assumptions (age, metallicity, etc.) and often lead to the largest stellar masses,   M M log > 10; (2) models with a flatter, gray attenuation law provide a better fit the UV region and lower stellar masses; and (3) composite SEDs with a dusty galaxy and either a blue galaxy or a blue QSO dominating the SW bands also provide a good overall fit to SED and similar masses with scenario (1).For scenarios where the LW bands are dominated by an obscured AGN, we find that models based on an obscured QSO plus a blue galaxy dominating the SW bands, or pure AGN models, where a combination of obscured and scattered emission by the torus dominates the whole SED, provide a good fit to the overall SED and lead to stellar masses for the galaxy host that are two orders of magnitude lower than in the dustygalaxy-dominated scenarios, with   M M log = 7-8.6.The unresolved, point-like sizes of all the EROs are more suggestive of an AGN-dominated scenario and disfavor a galaxy-dominated scenario where the blue and red SEDs are caused by different stellar populations in distinct regions of the galaxy.7. The NIRCam colors are not enough to break the SED model degeneracies and meaningfully distinguish between galaxy-or AGN-dominated scenarios.Additional MIRI photometry redward of F1000W, probing the rest-frame SED between 1 and 3 μm, will be able to answer this question definitely.8.The number densities do not favor any of the galaxy versus AGN scenarios either since both have potentially problematic implications if the extreme properties of some of these EROs are confirmed.The dusty-galaxy scenario would imply an increase in the number density of very massive galaxies,   M M log > 10.5, at z > 7 of up to a factor of ∼60, relative to the pre-JWST estimates even if just a handful of them are confirmed to be that massive.Similarly, in the obscured-AGN scenario, the large, dust-corrected UV luminosities would imply an unexpectedly large number, ∼10 −5 Mpc −3 , of obscured but luminous, L bol = 10 45-46 erg s −1 , QSOs at z > 7, more than three orders of magnitude larger than the observed density of unobscured QSOs with similar M UV ∼ −24 mag.

Figure 1 .
Figure 1.Color-magnitude and color-color diagrams showing the selection threshold for F277W EROs (circles; F277W − F444W > 1.5 mag), relative to the bulk of the CEERS galaxy catalog, color coded by stellar mass and A V , and a subset of F150W EROs (F150W − F444W > 2 mag).The blue and purple markers indicate the EROs observed with MIRI and NIRSpec, respectively.The black squares show the EROs from Labbé et al.(2023).The left and central panels show the general trends toward redder colors with increasing mass and dust attenuation (arrows), which suggest that F277W EROs are massive and dusty galaxies.However, the central panel reveals that F277W EROs have surprisingly blue colors at SW, F150W − F200W ∼ 0 mag, very different from those of F150W EROs and massive dusty galaxies in general.The red square shows a massive, dusty, submillimeter galaxy at z = 5.1 fromZavala et al. (2023) which is also red in all bands.This implies that F277W EROs have bimodal SEDs with blue SW colors and red LW colors.The right panel shows the correlation between photometric redshift and F150W − F277W color for the F277W EROs.As the F277W filter shifts from the steep, rest-frame optical range to the flat rest-frame UV range with increasing redshift, the color declines to F150W − F277W ∼ 0 mag.

Figure 2 .
Figure 2. NIRCam color-color, F115W -F200W vs. F277W -F444W diagram showing the bulk of the CEERS galaxy population (gray scale) and the EROs (circles) selected in Section 3.1 based on their characteristic blue-red colors in the SW and LW bands.The colors are the same as in Figure 1.The solid and dashed black lines depict the color tracks as a function of temperature (T = 500-1500 K) and metallicity (log(Z/Z e ) = −1 and 0) derived from LOWZ brown dwarf stellar templates (Meisner et al. 2021).While brown dwarfs also appear to have blue-red SEDs in the NIRCam bands, their SW colors are typically bluer, F115W -F200W < −0.5, than the ERO population at similar LW colors.Based on this distinction we identify two potential brown dwarf candidates in our sample of 37 EROs.

Figure 3 .
Figure 3. Left: photometric redshift vs. FAST stellar mass diagram for the F277W EROs (circles, color coded as in Figure 1), the CEERS galaxy sample (green density map), and the F150W EROs (red).For comparison, we also show the galaxies from Labbé et al. (2023) using their redshifts and stellar masses (squares).The F277W EROs are relatively massive,  

Figure 4 .
Figure 4. Left: stacked SED of the 37 EROs (gray squares) divided into two groups below and above z = 7, shown in purple and red, respectively.The MIRI photometry is shown in green.All galaxies exhibit a characteristic bimodal SED.Representative best-fit SEDs with EAZYPy at z = 5.5 and z = 7.5 (purple and red solid lines, respectively) show that this peculiar SED shape is typically reproduced by a composite SED with a blue, flat continuum in the rest-frame UV and red, steep continuum in the optical.Indeed, the best-fit power law to the fluxes redward of F277W is quite large (α ν ∼ 3-4), indicative of a heavily reddened continuum.The stacked SEDs also highlight the difference in F277W as the bands shift from the steep to the flat slope with increasing redshift.Right (top): 2 5 × 2 5 cutouts of EROs in the two redshift bins showing their similar compact and featureless visual appearances.Right (bottom): list of some of the strongest emission lines that can potentially cause emission-line-driven excesses in the NIRCAM and MIRI photometry at different redshifts.The locations of the strongest Hα and [O III] lines are also indicated in the left panel.

Figure 6 .
Figure6.Multiband 2 5 × 2 5 cutouts of the MIRI-detected EROs and best-fit SED models computed with EAZYpy, Prospector, Synthesizer, and a hybrid of a galaxy plus a red QSO template (either QSO1 or QSO2) fromPolletta et al. (2006).The left panels illustrate that fits based on a single stellar population component provide a good fit to the overall LW NIRCam and MIRI photometry (black and green squares) but they systematically fail to reproduce the rest-frame UV probed by the SW NIRCam bands.The middle panels show that a composite model consisting of two (or more) stellar populations provides an excellent fit to all the bands by combining a red, massive, and dusty galaxy that fits the LW bands and a blue, low-mass galaxy that fits the SW bands but has little impact on the stellar mass.The right panels show that the hybrid galaxy + QSO model (QSO1 and QSO2, orange and red, respectively) provides an equally good (or better) fit to the SED than the other models.Here, a dust-obscured QSO dominates the LW photometry but does not contribute to the stellar mass of a blue unobscured host, and consequently leads to total stellar masses ∼two orders of magnitude smaller than the other scenarios.The two stellar templates (gray) illustrate the 16%-84% confidence range in stellar mass for the galaxy component.The SEDs exhibit similar UV emission but increasingly larger optical emission with mass.

Figure 7 .
Figure 7. Multiband 2 5 × 2 5 cutouts and 2D/1D NIRSpec spectra of the NIRSpec-detected EROs.The best-fit SED models computed with EAZYpy, Prospector, Synthesizer, and a hybrid of galaxy plus the QSO2 template from Polletta et al. (2006) are the same as in Figure 6 but fixed to the spectroscopic redshift.

Figure 9 .
Figure 9. Ranges of stellar masses and A V values obtained with different SED modeling assumptions for the eight EROs with MIRI and NIRSpec observations.Overall, the values derived with the commonly used EAZYpy and FAST methods provide similar estimates as the fiducial Prospector-τ model, and they are typically the largest (red, gray, and light blue markers, respectively).The Prospector-np nonparametric model (dark blue) leads to smaller stellar masses by 0.4 dex, on average.A nonparametric model with a gray attenuation law, Prospector-np-cf (black), can lead to even smaller masses by 0.7 dex when MIRI fluxes are available, but it obtains similar values to the fiducial models where they are not.The stellar masses from Synthesizer (green) are the smallest by ∼1 dex relative to the fiducial values, but the accuracy of the fit is worse.The values obtained with the hybrid galaxy plus obscured QSO model (not shown) are much smaller,   M M log = 7-8 because the QSO dominates the SED without contributing to the stellar mass of the blue, low-mass host.

Figure 10 .
Figure 10.Left: galaxy number densities with stellar masses above

Figure 11 .
Figure 11.UV luminosity functions at redshifts z = 5-7 and z = 7-9.The gray markers show the density of galaxies and AGNs from the literature identified with different criteria (QSO, X-ray, and broad-line detection).The red circles indicate the density of the full sample of EROs if they are obscured QSOs.The filled and open markers show the different densities computed either from the observed or the dust-corrected UV luminosities.While the observed density and luminosities of the EROs are roughly similar to those of the X-ray population at z ∼ 6, if they are obscured AGNs, their intrinsic luminosities are much larger, L bol ∼ 10 46 erg s −1 , comparable to the faint end of the bright QSO population but nearly three orders of magnitude more abundant.Such luminosities would also imply that the EROs have black hole masses of the order of  M log 8 BH , nearly as large as the estimated stellar masses of their blue galaxy hosts, which would lead to unexpectedly large mass ratios of M BH /M galaxy = 0.8.

Figure 12 .
Figure 12.Color composite (F277W + F356W + F444W) 2 5 × 2 5 cutouts of the 29 other EROs in the color-selected sample.Similar to the eight primary sources in Figures 6 and 7, these objects are also very red and remarkably homogeneous and compact.

Table 1
MIRI and NIRSpec EROs at 5