The Galaxies Missed by Hubble and ALMA: The Contribution of Extremely Red Galaxies to the Cosmic Census at 3 < z < 8

Using deep JWST imaging from JADES, JEMS, and SMILES, we characterize optically faint and extremely red galaxies at z > 3 that were previously missing from galaxy census estimates. The data indicate the existence of abundant, dusty, and poststarburst-like galaxies down to 108 M ⊙, below the sensitivity limit of Spitzer and the Atacama Large Millimeter/submillimeter Array (ALMA). Modeling the NIRCam and Hubble Space Telescope (HST) photometry of these red sources can result in extremely high values for both stellar mass and star formation rate (SFR); however, including seven MIRI filters out to 21 μm results in decreased masses (median 0.6 dex for log10(M∗/M⊙) > 10) and SFRs (median 10× for SFR > 100 M ⊙ yr−1). At z > 6, our sample includes a high fraction of “little red dots” (LRDs; NIRCam-selected dust-reddened active galactic nucleus (AGN) candidates). We significantly measure older stellar populations in the LRDs out to rest-frame 3 μm (the stellar bump) and rule out a dominant contribution from hot dust emission, a signature of AGN contamination to stellar population measurements. This allows us to measure their contribution to the cosmic census at z > 3, below the typical detection limits of ALMA (L IR < 1012 L ⊙). We find that these sources, which are overwhelmingly missed by HST and ALMA, could effectively double the obscured fraction of the star formation rate density at 4 < z < 6 compared to some estimates, showing that prior to JWST, the obscured contribution from fainter sources could be underestimated. Finally, we identify five sources with evidence for Balmer breaks and high stellar masses at 5.5 < z < 7.7. While spectroscopy is required to determine their nature, we discuss possible measurement systematics to explore with future data.


INTRODUCTION
Our picture of galaxy formation and growth during the first 2 billion years of the Universe is, currently, an incomplete one.The Hubble Space Telescope (HST) is only sensitive to the rest-frame ultra-violet (UV) light for galaxies at z > 3, which at Cosmic Noon (1 < z < 3), is known to miss the overwhelming majority of star formation activity and energy output from galaxies due to dust obscuration (e.g.Madau & Dickinson 2014).While the addition of sub-millimeter and radio measurements have greatly improved our understanding of the star formation census at Cosmic Noon (e.g.Hodge & da Cunha 2020), the difficulty in making these measurements increases dramatically at z > 3.This is in part due to a decline in the abundance of the bright, sub-millimeter galaxies that are detectable by wide-area sub-mm surveys (e.g. Simpson et al. 2014;Brisbin et al. 2017;Danielson et al. 2017;Casey et al. 2021).It is also in part because of the large uncertainty in the abundance of more common, but less-extreme galaxies, that are below the detection limits of typical multiwavelength surveys (and difficulties in determining their counterparts and redshifts, e.g.Smail et al. 2021).As the abundance of massive galaxies decline, less massive and fainter sources (e.g.Pope et al. 2017Pope et al. , 2023) ) could plausibly harbor an increasing fraction of the obscured star-formation at early times.Progress has been made by identifying fainter, dusty galaxies with ever increasing ALMA area (in particular using longer wavelength > 2 − 3 mm selections to filter out lower-redshift sources Béthermin et al. 2015;Casey et al. 2021;Cooper et al. 2022).Despite these advances, recent ALMA-based estimates at z > 4 still rely on relatively few detections, strong lensing, or extrapolations of the infrared luminosity function (e.g.Zavala et al. 2021;Algera et al. 2023; Barrufet et al. 2023a;Traina et al. 2023;Fujimoto et al. 2023a).Uniform and large-area surveys are still needed for unbiased samples, and remain challenging with ALMA's small field of view.Thus, the pre-JWST census of galaxies in the first 2 Gyr of cosmic history remained highly biased to relatively unobscured, starforming sources.A complete understanding of the evolution of the early Universe requires accurate accounting of the abundance, the growth rate, and the energy output from a complete and unbiased census of the galaxy population.
Prior to the launch of JWST (Gardner et al. 2023), observations at the limit of Spitzer Space Telescope and ALMA capabilities indicated the existence of surprisingly abundant, massive galaxies that are very red, enough to evade HST selection.Some of these massive galaxies are serendipitiously identified at submillimeter wavelengths to faint limits, thanks to ALMA's extremely sensitive receivers (otherwise lacking counterparts at shorter wavelengths; i.e. optical/near-infraredfaint or "dark").Sub-millimeter or radio sources lacking counterparts have been known for decades and hypothesized to be high redshift (e.g.Hu & Ridgway 1994;Hughes et al. 1998;Dey et al. 1999;Dunlop et al. 2004;Frayer et al. 2004).However, in recent years, ALMA's sensitivity, small field of view, and ever deeper surveys at optical/near-infrared suggest that massive, dustobscured and red galaxies were even more common at high-redshift than previously thought (Williams et al. 2019;Umehata et al. 2020;Fudamoto et al. 2021;Sun et al. 2021).The number densities and redshift distributions of these red sources are still highly uncertain, but imply that previous HST surveys might have missed up to 90% of massive galaxies (Log 10 M/M ⊙ > 10.5) at z > 3 (Wang et al. 2016(Wang et al. , 2019)).While the idea that the distant Universe may contain abundant massive galaxies waiting to be discovered is tantalizing, ALMA left many open questions regarding the properties of the stellar populations, morphologies, and potential for these missing galaxies to harbor black holes.These questions remain unanswered without deep near-to mid-infrared data.
The first year of JWST observations has now confirmed the existence of very red galaxies missed by HST.First-look papers indicated that the properties of these sources are consistent with being dusty or quiescent galaxies at z ∼ 3 − 7, with some fraction of bluer sources turned red by strong emission lines boosting long-wavelength filters (Barrufet et al. 2023b;Rodighiero et al. 2023;Pérez-González et al. 2023;Endsley et al. 2023;Smail et al. 2023;Fujimoto et al. 2023b;Barger & Cowie 2023).The new discovery of these red galaxies demonstrated conclusive evidence that HST's limited wavelength coverage had left us with a critical blind spot.
Interestingly, NIRCam imaging is revealing that distant red and optically-faint galaxies are remarkably diverse as a population.In particular, a number of galaxies with peculiar or surprisingly shaped red spectral energy distributions (SEDs) have been discovered.Some SEDs resemble strong Balmer breaks with strong red rest-optical continuum leading to high inferred stellar masses despite their early observation times (e.g.Labbé et al. 2023a).Another class exhibit similarities in SED and morphology to highly reddened QSOs (visibly appearing as "little red dots"; LRDs; Labbé et al. 2023b;Furtak et al. 2023a).A number of such LRDs have now been spectroscopically confirmed to host some form of AGN activity (Kocevski et al. 2023;Harikane et al. 2023;Greene et al. 2023;Matthee et al. 2023;Maiolino et al. 2023;Furtak et al. 2023b;Kokorev et al. 2023).Both subclasses raise important questions about how these galaxies and black holes could grow at such fast rates at early cosmic time, and motivate the identification of larger samples with expanded deep panchromatic data.
Among the broader red and massive z > 3 galaxy population, some first results have indicated that in fact their abundance has been previously significantly underestimated, with a factor of 20-25% higher stellar mass density at 3 < z < 6 and > 100% at 6 < z < 8 (Gottumukkala et al. 2023).This points to an early start to rapid stellar mass growth that is obscured by large dust columns, assuming our extrapolations from z ∼ 0 stellar population modeling are correct (e.g.Steinhardt et al. 2023;Wang et al. 2023;Woodrum et al. 2023).This could alter the picture of early galaxy growth, requiring increased efficiency of star formation than previously thought (Labbé et al. 2023a;Boylan-Kolchin 2023;Xiao et al. 2023a).However, further exploration remains to be done, since stellar mass and photometric redshift estimates at high mass and high redshift can also be sensitive to the data set used, (e.g.inclusion of near-IR medium-band filters; Desprez et al. 2023), thus motivating new exploration also using expanded wavelength coverage.
With this paper, we provide a detailed look at the SEDs and physical properties of galaxies that were missed by previous surveys using HST, and in many cases even Spitzer and ALMA, primarily due to their extremely red colors and relatively faint fluxes.The goal of this work is to determine what types of galaxies we have been missing, and, what fraction of the cosmic census has been unaccounted for as a result.This work will focus on a uniformly selected sample of extremely red sources (JWST/NIRCam 1.5µm-4.4µmcolors > 2.2) that were previously unaccounted for by HST (fainter than 27 ABmag at 1.5µm), a criteria that identifies massive and obscured sources primarily at z > 3 by bridging the Balmer break at high-redshifts (e.g.Caputi et al. 2012;Wang et al. 2016;Barrufet et al. 2023b).
In Sections 2 and 3 we overview the panchromatic data we use in the characterization of these sources.In Sections 4 and 5, we measure the physical properties using SED fitting, revealing both sources with high attenuation, and those with older stellar populations, and dustreddened AGN candidates.We also highlight how their extreme colors make it difficult to measure their properties using HST+NIRCam alone, and how the deepest MIRI plus ALMA (only available in this field) can alter the interpretation of these sources.Finally, in Section 6 we characterize their contribution to the galaxy census, identifying what fraction of star formation and stellar mass density was previously missing.We assume a ΛCDM cosmology with H 0 =70 km s −1 Mpc −1 , Ω M = 0.3, Ω Λ = 0.7, and a Kroupa 2001 initial mass function (IMF)

DATASETS
In this work, we focus on a select sample of optically faint galaxies identified using JWST/NIRCam imaging inside the Great Observatories Origins Deep Survey South (GOODS-S; Giavalisco et al. 2004), overlapping with the Hubble Ultra Deep Field (Beckwith et al. 2006).Given the exceptional ancillary data that is available in this prime region of the sky, we focus on providing a detailed look at this population to give broader insight for larger samples with more limited data.In this section we outline the various datasets we use and our photometric methods to measure fluxes in each.

HST Imaging
Figure 1.F150W-F444W color vs F444W magnitude illustrating our selection from the JADES catalog with S/N > 20 in F444W (black points).We identify very red sources by their F150W-F444W> 2.2 color (all colored points) color coded by their EAZY photometric redshift.We further identify the optically fainter sample whose colors and magnitudes are consistent with z > 3 massive galaxies that were absent from earlier HST selected samples (i.e. have F150W < 27 ABmag: the black diagonal line).In this work, we focus on the optically faint sources behind the black line marked as colored squares.Existing samples of optically faint objects (maroon open circles) are collected from the literature using pre-JWST datasets (Wang et al. 2016(Wang et al. , 2019;;Franco et al. 2018;Alcalde Pampliega et al. 2019;Xiao et al. 2023b).We note the detection limit of the Spitzer GREATS program (dotted line) indicating newly discovered sources by JWST.
In our analysis we include deep optical and near-IR HST imaging from the Advanced Camera for Surveys (ACS) and Wide Field Camera 3 (WFC3) as compiled by use of the Hubble Legacy Field (HLF; Illingworth et al. 2016;Whitaker et al. 2019, and references therein) imaging in GOODS-S from HST.The HLF represents the deepest composite imaging including nine filters between 0.4-1.6µmwavelength.We include data from ACS (F435W, F606W, F775W, F814W, F850LP) and WFC3/IR (F105W, F125W, F140W, F160W).

Photometric catalogs
Also following the methods outlined in the JADES first data release (Rieke et al. 2023b) we measure photometry jointly for the HST and NIRCam imaging based on a detection image constructed as an inverse varianceweighted stack of the NIRCam longwave filters we have across the entire GOODS-S footprint: F277W, F335M, F356W, F410M, and F444W.Photometry is measured via a series of steps that perform source detection with advanced deblending algorithms (see Rieke et al. 2023b, Robertson, in prep) that utilize photutils, astropy, scipy, cupy and sextractor packages.We use the photutils package to identify source detections, and measure forced circular aperture photometry (convolved to the resolution of the F444W filter) for all detected objects.For compact or unresolved sources we use a 0.5" diameter aperture, which we validate is the appropriate size by visual inspection.However, by eye we identified a fraction of sources are more extended than this aperture, and in those cases we instead use an 0.7" diameter aperture.To estimate photometric uncertainties, we combine in quadrature the Poisson noise and the noise estimated from the root mean squared (RMS) of 100,000 random apertures (following methodology outlined in Labbé et al. 2005;Quadri et al. 2007;Whitaker et al. 2011).We generate aperture corrections to account for flux lost from the fixed apertures using encircled energy curves constructed using WebbPSF (Perrin et al. 2014) as described in Ji et al. (2023).

Forced photometry
Once we identify our sample (as will be discussed in Section 3), we also include forced photometry on the following data measured based on the NIRCam positions of our galaxies.

JWST MIRI Imaging
We include wide-field JWST/MIRI (Wright et al. 2023) imaging at 5-25.5µm from the Systematic Midinfrared Instrument Legacy Extragalactic Survey program (SMILES; PID 1207; PI Rieke), which covers 34 square arcminutes of the GOODS-S/HUDF region with a nearly complete overlap with JADES and JEMS imaging.For each of the 15 pointings, ∼ 2.2 hr of sci-ence time was spread between 8 MIRI filters1 , reaching 5σ point source sensitivities of 0.20, 0.19, 0.38, 0.59, 0.68, 1.7, 2.8 and 16 µJy for F560W, F770W, F1000W, F1280W, F1500W, F1800W, F2100W, and F2550W respectively, in apertures containing 65% of the encircled energy for each PSF (see below).Survey design, data reduction, and source extraction are described in Lyu et al. (2023) 2 .Briefly, data reduction was performed using the JWST Calibration Pipeline v1.10.0 (Bushouse et al. 2023) with custom external background subtraction (see Álvarez-Márquez et al. 2023) and astrometry correction matched to the JADES astrometric solution.Flux calibration is updated based on the new reference files released in CRDS jwst 1130.pmap, which include the time-dependent count rate loss (Gordon, in prep.).
The blind MIRI photometric catalog is made using the same pipeline as the NIRCam catalog (Rieke et al. 2023b) with the differences that F560W and F770W are stacked for the detection image and the source detection and deblending are optimized for MIRI, given its noise properties and source density.Aperture corrections are applied based on model PSFs from WebbPSF (Perrin et al. 2014) for F1000W-F2550W and with empirical PSFs for F560W and F770W to account for the "cruciform" detector artifact (Gáspár et al. 2021, see Lyu et al. 2023 for more details).The JADES NIRCam sources that also have a blind SMILES MIRI counterpart are then matched within 0.15 ′′ .Given MIRI's remarkable sensitivity and well-behaved noise, we also include forced photometric measurements using seven MIRI filters out to 21µm at the NIRCam positions for sources that are not formally detected in the blind MIRI catalog, in order to enable better constraints on the observed mid-IR SEDs (restframe near-infrared).To perform forced photometry, we use photutils to measure flux in circular apertures of size 0.7" arcsecond diameter (comparable to the FWHM at F2100W, which balances the higher resolution MIRI shortwave bands against minimizing the aperture corrections in the longwave bands).Our forced photometry is verified by repeating this process on detected sources and comparing to the catalog fluxes, which are found to be in good agreement at the few percent level.To assess the uncertainties, we measured the median flux in boxes at random sourceless locations across the map to get the RMS in MJy/sr and then to estimate the RMS in the given aperture, we scale by the square root of the aperture area.

ALMA imaging
Extensive ALMA ∼1 mm imaging exists within the JADES region of the GOODS-S field, including some of the deepest contiguous ALMA maps on the sky to date.We take advantage of this ancillary data in order to constrain the far-infrared emission.This critical constraint allows a more robust estimate of the dustobscured star formation, given the relatively featureless power-law SEDs in the rest optical and faintness in the rest-UV.In cases where galaxies are significantly detected at 1-mm in either the deep ASAGAO imaging (Hatsukade et al. 2018) or the wider-area GOODS-ALMA program (Franco et al. 2018) we use the published total fluxes in our SED modeling.
For cases where galaxies are not detected in any of the available ALMA imaging, the deep data still provides strong constraints on the limiting far-infrared luminosity associated with galaxies.To incorporate these deep ALMA limits as constraints in our SED-modeling, we perform forced photometry at the locations of our NIRCam detections using the deepest available map at the sources' position.The data we use for our analysis include the deep combined map built by the ASAGAO team (Hatsukade et al. 2018) combined with data from the deep HUDF ALMA footprint (Dunlop et al. 2017) and overlapping coverage from the shallower, but much wider GOODS-ALMA program (Franco et al. 2018;Gómez-Guijarro et al. 2022).We prioritize the combined map from Hatsukade et al. (2018) if objects appear in multiple datasets, and provide GOODS-ALMA limits (map σ rms ∼ 180µJy/beam) for the objects outside the ASAGAO coverage.
The resolution of the combined ALMA data is relatively high (beam size ∼ 0.5" × 0.5" FWHM) similar in size to our NIRCam photometric aperture (Hatsukade et al. 2018).Without knowledge of the actual size of our sources at far-infrared wavelengths, to ensure we include any extended flux, we perform forced aperture photometry at the location of the NIRCam centroid, using a 1.0" diameter.To perform forced photometry on ALMA, we broadly follow the methodology of Sun et al. (2021); Shivaei et al. (2022); Betti et al. (2019).In short, we multiply the map units (flux/beam) by the number of pixels per beam, and use photutils (Bradley et al. 2022) to measure the flux within a circular aperture of 1.0" diameter from the primary beam corrected map.We assess the uncertainty as the standard deviation of the fluxes within 100 randomly placed apertures within 30" distance from the NIRCam source.Where the source is within 30" of the edge of the mosaic, we decrease this distance to 15", and do not assess forced photometry for sources closer than 15" to the edge of the map.

JWST spectroscopy
For this work, we will cross-match the samples we select in Section 3 with available spectroscopy from JWST.We make use of the public spectroscopic confirmations released by the JADES-NIRSpec program (Bunker et al. 2023;Eisenstein et al. 2023b).We also use data from the First Reionization Epoch Spectroscopic COmplete Survey (FRESCO; PI: Oesch; PID: 1895; Oesch et al. 2023) that obtained 4µm NIRCam wide-field slitless spectroscopy.We will discuss these further in Section 4.1.

COLOR SELECTION OF TARGETS
In this paper we target extremely red sources at z > 3 that would have been missed by typical pre-JWST surveys given wavelength coverage and detection limits.To build an inclusive sample of red sources, we adopt a selection based on the HST/WFC3 F160W -Spitzer 4.5µm colors (H-[4.5]> 2.3; e.g.Caputi et al. 2012;Wang et al. 2016Wang et al. , 2019)).Due to the redshifting of the Balmer break beyond 1.5µm at z > 3, this selection targets massive and red (including both dust obscured and quiescent) high-redshift galaxy candidates.At typical detection limits of HST, such galaxies were referred to as "H-band dropouts" (i.e. because they are red, and often not detected in HST; see Wang et al. 2019).For this work, we adopt a comparable selection based on JWST/NIRCam colors F150W-F444W > 2.2, slightly less restrictive in order to include existing objects that were identified pre-JWST based on data with larger photometric scatter.Since the first year of GOODS-S imaging from JADES yielded smaller area coverage from F150W than in F444W, we additionally check for additional H-dropouts using the F160W imaging for which we have the other NIRCam filters.We identify two additional sources which happen to fall off of our F150W footprint, but that otherwise have nearly complete NIR-Cam and MIRI imaging.
We present our red sources in F150W-F444W vs F444W magnitude space in Figure 1, color coded by their initial photometric redshifts measured with EAZY (see Section 4.1), along with all significant JADES NIR-Cam detections (black).This figure highlights the relative rarity of galaxies with this extreme red color compared to the full JADES sample, especially at the fainter F444W magnitudes (which can be seen to roughly correspond to higher-redshift sources based on the EAZY color coding).After visual inspection to remove imaging artifacts and hot pixels, our selection identifies a parent sample of 240 red galaxies accross redshifts, that are both F150W-F444W > 2.2 and F444W S/N > 20 (all colored points in Figure 1).
Following Barrufet et al. (2023b), we retain the F150W-fainter sample for the primary analysis of this paper, which roughly corresponds to a redshift selection (see their Figure 1) to identify the primary targets of study: those with F150W magnitude fainter than 27 (the dashed diagonal line), which are sources that are below the typical detection limit of HST/F160W in GOODS-S.We will refer to this primary sample following convention in the literature as "HST-dark", of which we identify 66 candidate objects (colored squares in Figure 1).Additionally, we note that ∼20% of those sources are also below the detection limit of the deepest Spitzer 3.6 and 4.5 µm imaging to date from the GREATS program (Labbé et al. 2015) which reached a 3σ limiting depth of 27.2 (AB magnitude).H-dropouts rightward of this limit (dotted vertical line) represent galaxies that are newly discovered by JWST.
Our HST-dark sources are listed in Tables 1 and 2. NIRCam cutouts of our sample, along with the aperture used to measure their photometry, are shown in Figures 11.The corresponding HST/F160W imaging is also shown.
The sample we target are very red in the restframe optical, exhibiting significant overlap in color space with very late-type stars.Thus, we check whether any of our objects are both unresolved, and have F115W, F150W, F277W and F444W colors that look like red stellar sources (e.g.brown dwarfs with effective temperatures < 1500K) and are poorly fit with stellar population synthesis models, yet are well fit by Sonora (2018) substellar atmosphere models (Marley et al. 2021;Hainline et al. 2023a).Only one source (ID 190413) is a highly probable brown dwarf candidate based on its colors.This source also has a 0.62¨positional offset measured between GREATS 3.6µm Spitzer imaging and our JWST data (2σ confidence, given its low S/N in the Spitzer imaging of 5σ, suggesting the positional uncertainty is ∼1.2" / 5 = 0.2-0.3";see Hainline et al. 2023a).This source is additionally poorly fit by stellar population models with a high redshift and improbably high mass solution.Thus, based on the evidence for proper motion, we catalog its properties in Table 1 but remove it from our analysis and census measurements in Sections 5 and 6.The color similarity indicates the uncertainty involved in distinguishing red galaxies at these redshifts from late-type stars (see e.g.Hainline et al. 2023a;Burgasser et al. 2023).
Our sample has some overlap with a number of red optically faint galaxy populations studied to date with JWST.For context, we note what fraction of objects in JADES identified by their criteria also meet our criteria of F150W-F444W>2.2 and F150W<27 ABmag.We find that only 12% of massive double-break candidates as in Labbé et al. (2023a) overlap with our HST-dark sample, likely because of their requirement of a relatively blue SED between F150W-F277W<0.7.We will discuss this overlap further in Section 6.2.Our sources are also redder than 50% of sources selected based on F150W-F356W>1.5 in Pérez-González et al. (2023), owing primarily to the bluer color chosen to include sources at z > 2. Thus in comparison to these other samples, our objects tend to be those with consistently red SEDs from the rest-frame UV to near-IR, and, our sources are typically redder across the NIRCam wavelength baseline.

Comparison with previously identified optically faint galaxies in GOODS-S
Our sample includes the majority of known optically faint or HST-dark galaxies that were identified in previous searches in this field (see open circles in Figure 1 Wang et al. 2016Wang et al. , 2019;;Alcalde Pampliega et al. 2019;Franco et al. 2018;Xiao et al. 2023b).Five out of seven of the H-dropouts identified by Wang et al. (2016Wang et al. ( , 2019) ) that fully overlap with our imaging are identified by our parent sample H-dropout color selection, (all of which also overlap with Alcalde Pampliega et al. 2019).Five optically dark sources identified by pre-JWST data were not red enough to meet our selection (which were identified using different methods Franco et al. 2018;Wang et al. 2019;Xiao et al. 2023b) and a further four are brighter than our F150W magnitude limit.
We note that our sample includes the optically dark object AGS11 identified in Franco et al. (2018, our ID 279678) which had been identified based on a blind ALMA detection with no optical counterpart.JWST now demonstrates that it has one of the most extreme red colors among our sample (F150W-F444W>4), remains still undetected in all shortwave filters (despite our ultra deep imaging), and is also one of the largest, reaching nearly one arcsecond in extent.This source was previously hypothesized to reside in a confirmed overdensity using ALMA spectroscopy to be at z ∼ 3.4 (Zhou et al. 2020), although the source itself was not detected in the ALMA spectroscopy.Based on our photometry, we find a higher redshift of z ∼ 4.7 for this source (which likely explains their finding of no spectroscopic confirmation).Based on this comparison, our selection identifies a relatively unexplored sample in this field, below the detection limits of other multiwavelength selections (including ALMA).

Redshift estimations
Our sample exhibits very red and sometimes featureless SEDs.To assist in our more detailed SED modeling to infer the physical properties, we first measure preliminary photometric redshifts based on the HST and JWST/NIRCam photometry using EAZY (Brammer et al. 2008) as presented in Hainline et al. (2023b).We measure the redshift as the probability weighted average peak of the photometric redshift distribution without any priors.This preliminary photometric redshift measurement from EAZY is used to set a prior on redshift for our more detailed photometric modeling that we discuss in Section 4.2.
A small number of our sample also have spectroscopy.These include 9 sources with one emission line with S/N > 3 in FRESCO data (Oesch et al. 2023), using spectral extractions presented in Sun et al. (in preparation) combining JADES and FRESCO.Two sources that received slits from the JADES NIRSpec program also have one emission line with S/N>3 (Bunker et al. 2023;Eisenstein et al. 2023b).Below we describe our procedure for visually inspecting the spectroscopic redshifts in conjunction with the photometric PDF(z) measured using EAZY, and our process for deciding how to set the redshift priors that we use in our SED modeling in the case of uncertain spectroscopic redshift solutions.
For the case of the FRESCO sources, all nine sources have only one line, resulting in ambiguity in the redshift solution.For four out of the nine, the one line is weakly detected with 3 < S/N < 5.The five brighter, significantly detected FRESCO sources also only show a single emission line.Thus for all, the solution is heavily dependent on the most probable photometric redshift measured with EAZY.We find that this often leads to a degeneracy between Hα at z ∼ 5-6 and [OIII]λ5007 at z ∼ 7-8 (assumed to be the detected line, given its 3× brighter than the other [OIII] doublet).This is because our sources are all very red, and the Lyman break is often faint and poorly constrained, thus photometric boosting by strong rest-optical line emission tends to have a strong influence on the photometric redshift solution.Thus, EAZY sometimes yields comparable probability for both Hα or Hβ+[OIII] solutions.In these cases, we visually inspect the EAZY χ 2 surface in conjunction with the SED shape, the detected line's wavelength, and we also consider whether the physical parameters derived from the SED modeling are reasonable (see next Section).
After this iterative process, we find that for four sources with confidently detected lines (219000,154428,184838,204851) the photometric evidence clearly agrees with the FRESCO redshift solutions.For a fifth confidently detected single line source, 217926, we find that in fact the photometry (primarily near the Lyman break) supports an altered redshift solution at z = 5.04 rather than z = 7.6, assuming the securely detected line is Hα, not [OIII]λ5007.For three of the less-securely identified objects, (IDs 90354, 120484, 104849) the marginal emission lines are more difficult to interpret so close to the limiting signal to noise, although the candidate lines do have solutions that are consistent with the EAZY photometric redshifts.We decide to explore the SEDmodeling solutions that are retrieved for both the case where redshift is a free parameter, and also while fix-ing to the tentative redshift.We get consistent results within the uncertainties either using the tentative spectroscopic redshifts or leaving the redshift free.We therefore consider these spectroscopic redshifts as robust, but in any case this choice does not impact our results.
For the last marginal case (ID 126594) we find a lowconfidence marginal emission line which, assuming it is the stronger of the [OIII] doublet, puts this source at z = 7.9.In this case, the redshift lines up with one of 3 plausible (narrow) photometric redshift solutions based on the presence of emission line boosting to the photometry, although this redshift solution is not the one most favored by the photometry (which prefers a higher photometric redshift of z = 9.9).However, at the redshift preferred by EAZY (which is also the redshift preferred by our SED-modeling when redshift is left as a free parameter), the galaxy is quite bright and in excess of expected stellar masses given our small survey area (based on our SED-modeling procedure outlined in the next Section).Thus, we take the approach of comparing stellar masses measured for both cases.We find that at the z = 7.9 redshift solution, which is in agreement with the tentative line detection in FRESCO data, the SED modeling yields a more realistic stellar mass given our small area (see Section 6.2).Thus, we opt for this more conservative redshift constraint and fix to the marginal spectroscopic redshift.
In addition, two of our sources were observed as part of the JADES NIRSpec campaign (Bunker et al. 2023;Eisenstein et al. 2023b).ID 198459 was spectroscopically confirmed at z = 3.588, which is consistent with the photometric redshift we measured using EAZY (z = 3.65).For the second source ID 132229 however, the solution is less obvious.While the redshift is tentative (based on detection of the [OIII]λ5007 line at z = 7.247), and consistent with the EAZY redshift (z ∼ 7.5), we found that leaving the redshift as a free parameter yielded an inconsistently high photometric redshift of z ∼ 8.1.If modeled at such a high redshift, we find that the inferred stellar mass is Log 10 M * /M ⊙ ∼10.7, unphysically high for a z ∼ 8 source in a small area (which we will discuss further in Section 6.2), lending some credibility to the lower redshift solution.Therefore, for this source we chose to accept the tentative spectroscopic redshift at z = 7.247 in the modeling and results.

SED Modeling
To measure the more detailed physical properties, we use the prospector Bayesian code to model the SEDs (Johnson et al. 2021) using the Flexible Stellar Population Synthesis (FSPS) models (Conroy et al. 2009), MIST stellar isochrone libraries (Choi et al. 2016;Dotter 2016) and the stellar spectral libraries MILES (Falcón-Barroso et al. 2011).We use the MCMC sampling code dynesty (Speagle 2020), adopting the nested sampling procedure (Skilling 2004).Our fiducial prospector setup broadly follows that outlined in Ji et al. (2023) with a few alterations.We adopt the Madau (1995) IGM absorption model.We briefly summarize the other model assumptions and priors used here.
We use a non-parametric star formation history (SFH) composed of nine time bins with a constant star formation rate in each bin.We fix the first two bins to be at 0 − 30 and 30 − 100 Myr.Throughout this work, we will refer to the SFR as modeled by our SED-fitting as that inferred in the most recent 30 Myr (the latest time bin).The last time bin is assumed to be 0.85t H − t H where t H is the Hubble Time at the time of observation; and the remaining 6 bins are evenly spaced in logarithmic space between 100 Myr − 0.85t H .It has been shown that the recovered physical properties are largely insensitive to the number of bins used, when it is greater than 5 (Leja et al. 2019).We further adopt the continuity prior (to weight for physically plausible SFH forms, thus mitigating overfitting the data), which has been demonstrated to work well across various galaxy types (Leja et al. 2019).
We adopt the Byler et al. (2017) nebular continuum and line emission model.We set both the stellar metallicity and gas phase metallicities as free parameters and assume flat priors in logarithmic space (with log(Z * /Z ⊙ ) ∈ (−2, 0.19) and log(Z gas /Z ⊙ ) ∈ (−2, 0.5)).Ionization parameter U is also left as free parameter using a flat prior with log U ∈ (−4, 1).
We adopt the (Draine & Li 2007) dust emission model with priors as defined in (Williams et al. 2019) to allow for more flexibly hotter dust temperatures which may be prevalent at higher redshift (da Cunha et al. 2013).These include flat priors on the starlight intensity on dust grains U min ∈ (1, 25), and the faction of stars at U min , γ ∈ (0.01, 0.99).These parameters are related to T dust ∼ 18× < U > 1/6 K as in Draine et al. (2014).We also adopt flat priors on the polycyclic aromatic hydrocarbon (PAH) mass fraction, qpah ∈ (0.5, 4).
We assume a two-component dust attenuation model where the dust attenuation of nebular emission and young stellar populations, and of old stellar populations, are treated differently (Charlot & Fall 2000).For stellar populations older than 10 Myr, we assume the dust attenuation using the parametrization from Noll et al. 2009(i.e. a modified Calzetti et al. (2000) dust attenua-  2.84 +0.17 1.64 +0.17 −0.17 1.73 +0.17 tion law).Stellar populations younger than 10 Myr are assumed to have the same dust attenuation law as for the nebular emission (for further details on the various dust parameter priors, dependencies, and prior ranges, see Tacchella et al. 2022a;Ji et al. 2023).The dust model priors are set such that the V-band attenuation (Av) can vary between A V ∈ (0, 10) with a flat prior.We also include AGN dust torus templates from Nenkova et al. (2008a) and Nenkova et al. (2008b), with flat priors in logarithmic space for both the ratio of bolo-metric luminosity from the galaxy divided by that from the AGN (f AGN ∈ (10 −5 , 3)), and the optical depth of clumps in AGN dust torus at 5500 Å (τ AGN ∈ (5, 150)).
We use the photometric redshift measured using EAZY in the last section as a photometric redshift prior for prospector modeling (the mean of a Gaussian prior width ± 0.5).For the cases where the EAZY fit resulted in a redshift z > 8, we instead use a flat prior on the redshift to allow the possibility of lower-redshift solutions.For the sources with spectroscopic redshift con-straints, we fix to the spectroscopic redshift that was identified in the last Section.We further limit the S/N of any photometric point, which is capped at 20 (minimum 5% uncertainty, reflecting uncertainties in relative photometric calibration between filters).

Impact of MIRI+ALMA data
To date, photometric studies of similarly red galaxies at z > 3 have been restricted to NIRCam and HST data, with some limited wavelength coverage from 1-2 MIRI bands (Barrufet et al. 2023b;Pérez-González et al. 2023;Labbé et al. 2023a,b;Barro et al. 2023;Rodighiero et al. 2023;Akins et al. 2023;McKinney et al. 2023;Endsley et al. 2023).To explore the impact on recovered properties when using this more limited wavelength coverage, we run our prospector modeling for the 29 sources inside the MIRI footprint using only the HST+NIRCam data, and, again including also the MIRI+ALMA data.
In Figure 2 we present a comparison of the inferred best-fit parameters: photometric redshift, stellar mass, and SFR using HST+NIRCam only, with the result obtained when including the MIRI+ALMA data.We find that including photometry from both the seven MIRI filters and ALMA 1-mm data overall results in consistent redshifts as HST+NIRCam (with a few outliers; left panel of Figure 2).However, including MIRI+ALMA significantly alters other key parameters recovered using SED modeling.We find that the addition of MIRI+ALMA data serves to lower both the stellar masses and SFRs of galaxies compared with using HST+NIRCam alone.In particular, we find a systematic reduction in stellar mass of (median decrease of 0.6 dex, and as large as 1 dex) for galaxies with HST+NIRCam-measured Log 10 M * /M ⊙ >10.We similarly find that for SFRs in excess of ∼100 M ⊙ /yr as measured by HST+NIRCam alone we calculate a median factor of 10 decrease in SFR inferred when including MIRI+ALMA.These findings indicate that studies based on the more limited data sets are likely to overestimate both the SFRD and the stellar mass density, at high redshift significantly, in particular for galaxies where HST+NIRCam infer high masses and SFRs.This seems to be the case even despite the excellent medium band coverage of JADES.
We note that this finding may not be representative of the general galaxy population, since our sample has much more extreme red colors than typical galaxies at these redshifts.A systematic assessment of the impact of multiple MIRI bands on mass and SFR across redshifts has not yet been undertaken (although see simulations based on mock galaxies in Kauffmann et al. 2020;Kemp et al. 2019;Bisigello et al. 2019).However, analyses with more limited filters (F560W and/or F770W) suggest mixed results.A similar impact to ours was noted in Papovich et al. (2023) using a more representative sample of galaxies at similar redshifts, (although comparing to modeling with a much more limited set of only HST plus IRAC data).Those authors find that, on average, the stellar mass decreases by ⟨∆logM * ⟩ = 0.25 dex at 4 < z < 6 and 0.37 dex at 6 < z < 9.Those authors also find a systematic reduction in SFRs by ⟨∆logSF R⟩= 0.14 dex at 4 < z < 6 and 0.27 dex at 6 < z < 9. MIRI likely has a large impact in this case owing to a lack of NIRCam photometry, which can break degeneracies between continuum and emission lines, which translates to a less robust constraint on rest-optical continuum at high-redshifts.While their result is qualitatively in line with what we find for the most extreme of our red sources, accumulating evidence indicates that including some MIRI data (in particular, F770W, covering rest-frame ∼1-2µm) may not be essential to achieving more accurate estimates for typical galaxies.Helton et al (in preparation) finds that for 7 < z < 9 galaxies (using our same NIRCam filter set), that the stellar population models with and without MIRI data are similar.Additionally, Alberts et al (in preparation) find that stellar masses of Log 10 M * /M ⊙ >9 galaxies at 3 < z < 6 also demonstrate no systematic improvement when F770W is included.
Regardless, since in this study we have both multiple NIRCam medium bands plus multiple MIRI filters (mitigating uncertainties in constraining the stellar mass from both angles), our measurements are likely to be relatively reliable despite the exceptional character of our very red sources.However, this exercise demonstrates that the MIRI and ALMA data are of particular importance for sources with extreme red colors.Thus, caution should be exercised when interpreting the SED modeling of sources with red (and often featureless) SEDs without mid and far infrared wavelength coverage.Further, the uncertainties that arise from lacking the fuller wavelength coverage are not reflected in errorbars of mass and SFR measured with HST+NIRCam alone.

Properties of Optically faint (HST-dark) galaxies
In this section, we characterize the properties of our sample of very red sources, in particular highlighting the diversity of our sample.Our prospector modeling indicates that our sources range from 3 < z < 8, revising all of the higher-redshift solutions that were measured using EAZY (owing to non-detections in the observed optical and near-infrared; Figure 1).We also find that our sample includes galaxies at a range of stellar masses from Log 10 M * /M ⊙ ∼ 8.2 − 10.8, with moderate median SFR (∼ 50 M ⊙ /yr), high attenuation (A V ∼ 2), moderately evolved stellar populations (mass-weighted age ∼ 250 Myr).These are similar to initial JWST explorations (Barrufet et al. 2023b;Rodighiero et al. 2023;Pérez-González et al. 2023).
As shown in Figure 11, our sample of H-dropouts have incredibly diverse morphologies, ranging from large extended disks (e.g.Gibson et al submitted; Nelson et al. 2023), to potentially merging clumps or groups, while a large fraction are remarkably compact, close to the resolution limit of F444W.This diversity is in line with expectations from simulations that 1) whether dusty galaxies are "dark" is a strong function of viewing angle for a range of galaxy types, and 2) the most massive galaxies will pass through this phase early in the Universe's history due to prodigious dust production (Cochrane et al. 2023).In the subsections below, we review the sub-categories of objects that we identify in our sample.

Dusty star forming galaxies in our sample
A minority of our sources exhibit detectable ALMA emission, indicative of dust obscured star formation at the level typical of dust obscured galaxies (DSFGs; SFR>100 M ⊙ /yr; Casey et al. 2014).Four sources between 3.6 < z < 5 show significant ALMA detections (> 4σ) in the range 0.5-1.0mJy/beam, well below those of prototypical submillimeter galaxies (their ALMA properties have been studied elsewhere; Franco et al. 2018;Hatsukade et al. 2018;Xiao et al. 2023b).The NIRCam and MIRI cutouts (in Figure 11 of the Appendix) show that 3/4 ALMA sources are the most extended disk-like objects in our sample, while the 4th (201793) is compact but clearly resolved in the shortwave filters.These sources are qualitatively similar to ALMA-only objects that have been identified by 2-3mm ALMA imaging with z > 4, total infrared luminosity log 10 L IR /L ⊙ ∼ 12 − 12.5, Log 10 M * /M ⊙ ∼ 10.5 − 11, SFR∼ 200−300 and ∼ 25 ABmag at 4µm (e.g.Williams et al. 2019;Manning et al. 2022, among others).Noting that our MIRI area is quite small, we find a similar abundance of ∼ 0.1 square arcmin −1 .
However, the overwhelming majority of our sources are not detected in any of the 1-mm ALMA data (25 out of 29) and thus have only upper limits to their total infrared luminosity.These upper limits are log 10 L IR /L ⊙ ≲12, based on integrating the 8-1100µm (restframe) best-fit SED from prospector (see Figure 7).While this means we can't robustly constrain how much lower is the exact level of obscured star formation individually, we stack the 1.1mm ALMA images for all sources that are not individually detected (calculating the inverse variance weighted average) from the GOODS-ALMA program (σ rms ∼ 180µJy/beam; Franco et al. 2018) for which we have coverage of our entire MIRI footprint.We do not find a significant detection in the stacked data at 3 < z < 4, or 4 < z < 5, and only a 2-sigma detection from the combined redshift range.However, a majority of the sources sit inside the deeper ASAGAO footprint.We perform the stack again for sources inside ASAGAO, finding again a marginal detection for 3 < z < 5 sources in our sample with 55±24 µJy/beam.The ASAGAO results from inside the MIRI footprint are shown in Figure 3.
It is possible that our limited sample size is too small for a significant detection.Therefore, we repeat the stacking while including our broader NIRCam-only sample within the GOODS-ALMA footprint, which roughly doubles the sample size.We also show the GOODS-ALMA stack from the broader NIRCam-only sample in Figure 3.We find that this test does yield a more significant (3σ) stacked detection in ALMA in the lowest redshift bin, 3 < z < 4. Collectively, we interpret these stacking experiments to indicate that our sample is likely dominated by faint sources with some small level of dust obscured star formation that are (individually) below the detection limit of ALMA at the low redshift end, z < 5.
The NIRCam imaging is sufficiently deep to have detected lower-luminosity analogs of DSFGs at even higher-redshifts z ∼ 6−7 (e.g. the serendipitious sources identified in Fudamoto et al. 2021).These ALMA-only sources were thought to have ∼25-26 ABmag at 4µm, with obscured SFRs in the range ∼40-70 M ⊙ /yr (based on f 1mm ∼ 110 − 190µJy) and Log 10 M * /M ⊙ ≲ 10.3.We identify nine sources with inferred properties that are consistent with these, noting that the stellar population modeling for three of those indicate substantially lower sSFRs to Fudamoto et al. (2021) To determine whether on average our sources at 6 < z < 8 may be comparable, we repeat the ALMA stacking for sources in this redshift range, and we find no stacked detection in GOODS-ALMA to a limit of f 1.1mm = 38 ± 41 µJy/beam (Figure 3; in ASAGAO, we find f 1.1mm = −9 ± 28µJy/beam, 1σ upper limits that are a factor of 4-7 below the 1-mm continuum detections in that work).This result suggests that our 6 < z < 8 sample is not actually dominated by similar lower luminosity DSFGs.
In Figure 4 we explore further whether galaxies are red due to age, dust, or both in the restframe U-V and V-J colors of our sources.All sources exhibit substantial dust obscuration (Av>0.7)except one.Although we do not see a strong dependence in Av with redshift, we note  (Franco et al. 2018), which covers our entire sample.Bottom panel: ALMA 1.1 mm stacks of a subset of our galaxies that are inside the footprint of the deeper ASAGAO imaging (individual detections removed).Only one source is inside the footprint at 3 < z < 4, so we instead show a combined 3 < z < 5 redshift bin.Our sample below z < 5 exhibit possible cold dust emission (2.2 − 3σ for both maps respectively), consistent with low luminosity DSFGs.Sources at z > 5 (including all LRDs) do not show detectable dust emission.Far right panels shows the ALMA stack of all 9 LRDs, which reaches a non-detection upper limit of 32µJy/beam.the most extreme sources with Av>2.5 are all at z < 6.5.We also note that by stacking the ALMA data in redshift bins, we find that at 3 < z < 4 and 4 < z < 5, the stacks show weak but more significant detections.This is consistent with a majority of the sources we identified at z < 5 being faint dust obscured sources.
While the stacked non-detection in ALMA at z > 6 could mean that the level of obscured star formation traced by cold dust is relatively low, (i.e.below the ALMA detection limit, log 10 L IR /L ⊙ <12), this is not definitive, and we unfortunately have limited data to further constrain this on an individual-galaxy basis.A lack of ALMA detection could also imply that our sources either 1) contain obscured star formation, but the dust is hotter than is typically assumed at lower redshift (see e.g.De Rossi et al. 2018), 2) contain hot dust heated by AGN activity, or, 3) the galaxies have primarily evolved or older stellar populations.In the following sub-sections we now explore these possible scenarios.

Possible evidence of AGN among red galaxies
To look for AGN evidence, we have matched our sample to the pre-JWST AGN catalog built by Lyu et al. (2022) that has integrated the Chandra X-ray, HST optical to near-IR, Spitzer mid-IR and JVLA radio data for a comprehensive search of AGN in the GOODS-S field.In total, we found only two matches among the more extended NIRCam-only sample (Table 2): JADES 171973 and 284527; both of them are identified as AGN by their high X-ray luminosity and the X-ray to radio ratio (see details in Lyu et al. 2022).As pointed out in Lyu et al. 2022, the AGN selection is complicated by the survey depth, wavelength range and object variations and many AGNs are still likely missed.
With the improved sensitivity and wavelength coverage of JWST data, significant progress has been made to identify AGN.Based on semi-empirical SED analysis of MIRI-detected sources with JADES NIRCam and SMILES MIRI photometry, Lyu et al. (2023) has drastically improved the AGN census in the central regions of GOODS-S.For our sample, three new AGN candidates have been revealed from that study: JADES 57356, 106502 and 204851.Notably, 204851 has been confirmed to be a broad-line AGN at z=5.48 in FRESCO data (Matthee et al. 2023).
Meanwhile, several groups have demonstrated the existence of fainter broad-line AGNs by selecting sources that feature as Little Red Dots (LRDs; Matthee et al. 2023;Labbé et al. 2023b;Greene et al. 2023, among other references) -objects with strong red continuum and compact morphologies in the NIRCam bands.Although the nature of these objects is still debated (e.g., Barro et al. 2023), the success rate of AGN search by this selection has been high (e.g., 9/12 in the sample of Greene et al. (2023) show evidence of broad Hα).While confirmed broad-line AGN may be prevalent among LRDs, it remains unclear whether the AGN dominates the host galaxy, and what its contribution is to the restframe UV and restframe optical continua.We now apply such selections to our sample and discuss their nature via the color/SED analysis with the addition of MIRI data points at longer wavelengths.Matthee et al. (2023) describe the LRD selection criteria as relatively flat or blue at observed 1 -2 micron, with a (very) red continuum from 2 -4 micron.We cross-check our sample with the following selection criteria for LRD: −0.5 <F115W-F200W< 1 and F277W-F444W>1.6 (Greene et al. 2023).Those criteria, based on spectroscopic confirmation of broad line AGN, should contain an estimated 80% AGN fraction (Greene et al. 2023).We find that 13 out of 66 objects in the NIR-Cam footprint are candidate LRDs by this criterion (9 of which lie in our MIRI coverage, excluding the probable brown dwarf candidate 190413).Two of the 13 include sources with evidence of an AGN, including 204851 (Lyu et al. 2023;Matthee et al. 2023), as well as 154428, which show possible evidence of broadened Hα and a weak narrow component in FRESCO data (Sun et al. in prep).While we do not fold in the explicit compact-ness cut of Greene et al. (2023) to identify LRDs, we note that all candidates identified by the LRD colors are visibly unresolved, or consistent with point sources, in F444W (including 204851, which in the single-band cutouts of Figure 11 appears blended with 2 neighbors).In general, our F150W-F444W selection does not pick up LRDs systematically, or, similarly, the extremely red object (ERO) selection that was used in Barro et al. 2023 (F277W-F444W>1.5).This is because while the LRDs and EROs are red in the long-wave NIRCam filters, some fraction are bluer in the shorter wavelength ones due to the rising rest-frame UV SEDs.The net result is that the bluer, more "V" SED-shaped sources are preferentially excluded by our H-dropout selection unless combined with a very red rest-optical continuum.

Nature of the LRDs in Our Sample
Our LRD subset has distinct SEDs that appear different from most of the H-dropouts.Some LRDs are also poorly fit by e.g.single component dusty or quiescent stellar populations.In the literature, this has prompted a number of explorations into differing origins for their blue rest-UV SEDs (Barro et al. 2023;Labbé et al. 2023b;Endsley et al. 2023;Furtak et al. 2023a;Greene et al. 2023).On the one hand, the rest-frame UV could be unobscured star formation.Alternatively, it could be scattered light from either hot stars or from an AGN accretion disk (depending on geometry).Similar ambiguity exists over the origin of their very red rest-optical SEDs, which could be driven by either obscured AGN continuum, or dust obscured stellar emission.Recently, deep ALMA non-detections for LRDs have provided evidence against red continuum produced by dust-obscured star formation, because the expected amount of re-processed dust emission in the far-IR under typical assumptions (e.g.dust temperatures in the range T dust ∼ 20 − 60 K) would be dramatically in excess of the deep ALMA limits (Labbé et al. 2023b).We find similar results for our LRD sample (see right panels in Figure 3).This could point to an AGN-dominant SED with hot dust emission.However, recent studies of compact dusty star forming galaxies at high redshift indicate that the dust temperature can be significantly higher due to the higher density of star formation (De Rossi et al. 2018;Sommovigo et al. 2020).So, this argument in favor of AGN dominance is not definitive.
Our MIRI photometry enables us to explore whether the LRD/AGN candidates from our sample are plausibly sources whose red rest-optical emission is dominated by a heavily reddened AGN (e.g.A V ∼ 1 − 4, based on the shape of the rest-optical continuum Greene et al. 2023).To explore this hypothesis, we plot our sources with MIRI coverage on two NIRCam-MIRI color-color diagrams in Figure 6 along with the redshift evolution in colors for a template dominated by an obscured AGN (i.e.excludes contribution from stars Lyu & Rieke 2018).From the empirical AGN SED library built in Lyu & Rieke (2018), we choose a model template for a typical (normal) AGN obscured by an extended dust distribution featuring large grains with an optical depth τ V = 3.This template matches the typical SEDs of lower-redshift "extremely red quasars" (ERQ; Ross et al. 2015;Hamann et al. 2017) that have similarly red F277W-F444W color as our LRD/ERO sample.
Looking at the observed colors, we find that all of our LRD candidates exhibit a turnover in their SED redward of F444W, between restframe 0.5-3µm.This turnover results in blue or flat F444W-F2100W color that is inconsistent with obscured AGN templates, which show a steeply rising shape at longer wavelengths (e.g.Lyu & Rieke 2022).These colors instead favor a scenario where the rest-optical emission between 0.5 − 3µm restframe is dominated by continuum from the stellar populations.This is potentially the spectral signature of the stellar bump at restframe 1.6µm (caused by a minimum in the H-opacity in the atmospheres of cool and low-mass stars that dominate the near-infrared spectrum of galaxies at age >10 Myr; e.g.Sawicki 2002).Reddened AGN, in contrast, typically exhibit a steeply rising red continuum redward of rest-frame 1.6µm, tracing hot dust emission from a torus (Alexander & Hickox 2012;Lyu & Rieke 2022).
While the MIRI data is deep enough to rule out rising red continuum from a dominant AGN for individual sources, the majority of the LRD subset are not detected in the longer-wavelength filters.To obtain a stronger constraint on the MIRI colors (on average) we median stack the MIRI imaging for the LRD subset.For this experiment we also include the 25µm band (F2550W) which otherwise is shallow relative to the other filters, but on average, can provide a meaningful constraint at restframe >3µm.The MIRI stacks are shown in the top panel Figure 5.
We measure average aperture photometry using the stacked MIRI images following the same procedure as Section 4. On average, the sample remains undetected at wavelengths >18µm, even in the deeper stacked image.For an accurate comparison to the average NIR-Cam SED of LRDs, we also calculate the median and inter-quartile ranges of the LRD photometric points (in lieu of stacking the NIRCam images, since all sources are already strongly detected in the NIRCam F444W).The full median-stacked SED for the sample of LRDs is shown in the bottom panel of Figure 5.For comparison, we also plot two example SEDs: one for a reddened quasar and its host galaxy, placed at the median redshift of our LRD sample (Mrk 231; Polletta et al. 2007), and, a prospector best fit model that is representative in shape of a number of our LRDs among the H-dropout sample (ID 121710).It is immediately clear that the turnover in the median MIRI SED strongly disfavors a reddened AGN shape, and instead agrees better with the moderately dusty and moderately old stellar model describing 121710 (mass weighted age ∼500 Myr, Av∼1).Notably, we find that the limits in the color from the median stack indicate that the SEDs of our sources must be quite flat between F444W and F2100W, with color F444W-F2100W∼ 0 (and, even consistent with blue color, i.e. a turnover, given the lack of significant stacked detection in F2100W).
Emission line boosting in F444W could also generate red F277W-F444W with blue F444W-F2100W colors, either from [OIII]+Hβ at z > 7 or Hα at z∼5.To explore this possibility, we also plot the rest-optical colors vs F444W-F770W (right panel of Figure 6) which should be very blue if F444W is boosted by emission lines (note that we choose F770W because F560W can also be contaminated by Hα emission when [OIII]+Hβ is in F444W).We find that the majority of the LRDs have weakly red or flat F444W-F770W colors, suggesting that the F277W-F444W color is not red due to significant line boosting in F444W.We additionally find that the median SED color of F444W-F770W=0.7, which is not consistent with the idea that the F444W-F2100W is blue or flat in spite of a red continuum because of strong emission line boosting F444W.These colors also suggest that the red rest-optical colors are dominated by stellar emission and not obscured AGN emission (which would predict a continual red rise into the mid-IR from AGN continuum rather than turnover from the stellar bump).However, the colors cannot rule out an obscured mid-infrared-dominant AGN whose continuum begins to dominate the SED at restframe wavelengths >3µm.
We note that the theoretical AGN templates provided with prospector may not correspond fully to reality; an alternative approach uses empirically based templates (see review by Lyu & Rieke 2022).To explore this approach, we run prospector on our full sample while using the modified and more realistic mid-infrared AGN template set (Lyu & Rieke 2018).We follow the prospector setup as outlined in Lyu et al. 2023.We find that, in line with our exploration of the NIRCam-MIRI colors, the contribution to the restframe 0.5-3µm continuum of LRDs is not obviously dominated by an obscured AGN.In reality, individual sources may exhibit a broad variety of behavior near restframe 3µm (where the presence of a mid-IR AGN is expected to be most obvious among typical star forming galaxies) and analysis is complicated by the effect of strong emission lines and low S/N of the longest wavelength MIRI data.
Further, we do a similar exploration using CIGALE (Boquien et al. 2019), which for samples in the literature has found that the steep slopes of the restoptical continuum for LRDs has a preferred origin from AGN continuum, based on HST+NIRCam photometry alone.For this experiment we use the SKIRTOR AGN model, a clumpy two-phase torus model from Stalevski et al. 2012Stalevski et al. , 2016.We find that, in nearly all cases (for the LRDs) CIGALE prefers fits where AGN continuum dominates the SED between restframe 0.5-0.8µm,but these AGN dominant models significantly overpredict the MIRI observed SED at F1800W-F2100W by a median/typical factor of 25-40 (well in excess of the photometric uncertainties).We also find that CIGALE overpredicts the MIRI photometry for the parent H-dropout sample, although to a lesser degree for non-LRDs (factor 2-3).This is likely driven by the redder rest-optical (F277W-F444W) color of the our subset of LRDs compared to the parent sample of H-dropouts (see left panel of Figure 6).
Thus, we would find similar conclusions based on HST+NIRCam photometry alone as the AGN-dominant solution identified by other studies in the literature.This highlights that while LRDs may be redder in the rest-frame optical than the non-LRD H-dropout sample, the MIRI data demonstrates that these SEDs are in fact mostly quite similar in the restframe near-infrared, suggesting a stellar origin for both subsets of H-dropouts (see right panel of Figure 6).
Given the consistency of the rest-optical and nearinfrared SED with dust obscured stellar emission, it remains puzzling why this does not translate to brighter far-infrared emission from re-processed energy by dust.
To explore this further, we stack the 1.1mm ALMA imaging covering our 9 LRDs from the GOODS-ALMA program (Franco et al. 2018) to obtain an average farinfrared flux (see Figure 3).We find that the 1.1mm flux is not detected (3±37µJy/beam).Similarly low stacked ALMA limits disfavoring dust-obscured star formation were found in Labbé et al. (2023b).It remains plausible that compact star formation at high redshift and low metallicity could heat dust well above typical expectations, and we discuss this scenario further in Section 6.1.

Evidence for older stellar populations in our sample
Such a low far-infrared flux measured in the ALMA stacks is consistent with a primary result from the SEDmodeling, which is that a large fraction of LRDs in particular are preferred to have older stellar populations over high levels of dust-obscured star formation.This is presumably, at least in part, a result of prospector trying to account for these deep upper limits from ALMA (given the fixed low dust temperatures assumed by the modeling).This possibility of quiescent SED solutions for LRDs was also explored in Labbé et al. (2023b), however that work determined that the extremely red rest-optical continuum slope disfavored a purely quiescent stellar population.However, we find that an older and evolved stellar population (mass weighted age > 200Myr), in combination with significant dust attenuation (Av > 1), adequately fits the SEDs of a number of H-dropouts, as well as a large fraction of our LRD subset.This is demonstrated in Figure 4, which shows that a number of our LRDs reside near the UVJ quiescent box (Williams et al. 2009;Muzzin et al. 2013).The typical level of dust attenuation that we infer is relatively high.However, owing to their faintness, we note that the uncertainties in the restframe colors near both the restframe U and J part of the SED are quite large.Thus, we caution against overinterpretation of the location of our sources in the UVJ diagram.
Regardless of the large errors in rest-frame U and J, we note that a number of the LRD SEDs with red U-V and V-J colors also exhibit clear evidence at high significance for strong Balmer breaks (e.g.IDs 121710, 132229, 219000 and perhaps also 154428) as does one non-LRD, 200576.The Balmer break evidence is clear in part due to the NIRCam medium band images, which at F182M, F210M and F335M are finely sampling the spectral region near 4000 Å across our entire redshift range.These strong Balmer breaks (all at early cosmic time z > 5.5) are remarkable, in particular because we simultaneously can rule out that the breaks are degenerate with strong line emission or AGN contamination with the 4µm medium bands and MIRI data.We will discuss their implication later in Section 6.2.
We note that for 154428 in particular, in addition to having evidence of a Balmer break, the FRESCO detection indicates possible evidence of broadened Hα with a weak narrow component, which could indicate a weak optical AGN (e.g.similar to the broad Hα in the quiescent galaxy confirmed at z = 4.6 in Carnall et al. 2023).Unfortunately, we are unable to robustly rule out the presence of an broad-line region in the other sources with FRESCO coverage, since any broad components may be dust obscured and below the FRESCO detection limit.
To summarize, the MIRI colors of AGN-dominated objects should be red, however, since we don't find evidence for that, we suspect that the light is instead dominated by stars.However, the conclusion that the LRDs in our sample are dominated by stellar emission is not definitive.In Section 6, we will show results including and excluding these sources from the sample, and discuss these results in the context of the assumption that these galaxies are dominated by stellar emission.

Excess UV emission
As a final note, we discuss the apparent presence of a restframe UV "excess" in a handful of our sources: flux which is not easily modeled without a secondary SED component (either by unobscured SF, or, scattered UV light from an AGN).We find clear evidence that composite SEDs are needed to explain the flat UV slopes in a number of objects (90354, 120484, 203749; and potentially also 81400, 132229, 183348).Our selection based on very red F150W-F444W colors may have rejected a number of LRDs with more obvious needs for composite SEDs.The slopes of the UV SEDs are at low S/N, but nonetheless consistent with typical slopes of either AGN or young stars (see discussion in Greene et al. 2023).Thus we cannot rule out models that have been proposed for similar sources with our data (Labbé et al. 2023b;Matthee et al. 2023;Barro et al. 2023).Spectroscopy has now confirmed that the continuum slope alone cannot easily differentiate between a SF or AGN origin using similar samples (Greene et al. 2023) nor their intrinsic luminosity (which is degenerate with fraction of light scattered).Thus we can't hope to do better with photometry, and take the presence of excess UV emission as an indication that scattered light from either an AGN or SF may contribute.We note that, if the origin of the UV flux is indeed from a frosting of unobscured star formation, the amount of SF is very small (i.e.these objects are barely detected in the very deep JADES imaging, and we measured the typical implied SFR(UV) based on the flux is <1 M ⊙ /yr).Thus any Figure 7. log10LIR/L⊙vs redshift of our sources with MIRI coverage that go into our estimate of the cosmic SFRD.For comparison, the average RMS limit of the ASAGAO imaging used is 60 µJy and GOODS-ALMA is 180 µJy.The majority of sources are inferred to be upper limits to the log10LIR/L⊙based on their non-detection upper limits from the ALMA imaging.
assumption about the origin of the UV flux will not impact our results in Section 6.1.Pursuing an explanation for the UV emission is outside the scope of this paper.
6. DISCUSSION: WHAT FRACTION OF THE GALAXY CENSUS WAS MISSED DUE TO "DARK" GALAXIES?
In this section, we explore how these previouslymissed galaxies may contribute to the star formation and stellar mass budget of the early Universe from 3 < z < 8, a regime that was previously incomplete among infrared measurements, and which could not be probed uniformly with earlier data.Given the results of Section 5.1, for this section we only consider sources for which we have MIRI and ALMA constraints, since we find major uncertainties that can change the mass and SFR estimates by up to 1 dex if we have only HST+NIRCam photometry.

The cosmic star formation rate density
Here we estimate the cosmic star formation rate density (SFRD) contribution from these sources (a previous estimate has been made for a similar sample by Barrufet et al. 2023b; we now include sub-mm and MIRI data that improves SFR constraints compared to HST+NIRCam data; see Section 5.1).Given that our sources are relatively bright in the detection band compared to the very deep JADES imaging (F444W S/N > 20) and the relative depths of our F444W and F150W imaging (Eisenstein et al. 2023a) we find we are sensitive to colors F150W-F444W<3.2 even for the fainter  (Zavala et al. 2021), along with the UV-based compilation in that work (unobscured, uncorrected for dust; blue shaded region).We also plot the SFRD contribution from our sources with log10LIR/L⊙<12 (maroon points).Although our MIRI data suggests our LRD subset are stellar-dominated, we show the SFRD when removing them since their nature is uncertain (open red circles).We also include an estimate based on the average ALMA flux from stacking the 1.1 mm imaging scaling with a hotter dust template (peach squares).sources that are not detected at F150W.Thus, we expect that selection of H-dropouts is complete for sources brighter than our limiting F444W AB magnitude of 29.4.
We do a simple estimate of the total SFRD by summing the total SFR among our sample (i.e. the SFR during the most recent 30 Myr, as derived from our SED modeling) divided by the cosmic volume within several redshift bins: z ∼3.5, 4.5, 5.5, and 7 (with ∆z ∼ 1, except for the highest redshift bin which has width ∆z ∼ 2).Without completeness our measurements should be considered lower limits, although we pick a high S/N where we are likely complete in both magnitude and color.
These results are shown as red circles in Figure 8, and are compared to the obscured SFRD and unobscured (and un-dust-corrected) SFRD from the MORA survey (orange and blue curves; Zavala et al. 2021;Casey et al. 2021).First, we find that at 3 < z < 4, optically faint galaxies make up a relatively small fraction of the obscured contribution to the SFRD.This is similar to cosmic noon where the SFRD is still dominated by brighter sources such as sub-millimeter galaxies (e.g.Dudzevičiūtė et al. 2021), and is consistent with earlier findings (e.g.Wang et al. 2019;Sun et al. 2021).However, at 4 < z < 6, we find that the missing population identified by our H-dropout selection likely contributes non-negligibly to the obscured fraction of the SFRD.We find that in these two redshift bins, our full sample is comparable to the total obscured fraction of cosmic SFRD characterized using existing ALMA and far-infrared observations (orange line).
Since ALMA and far-infrared detections become sparse at z > 3, the orange region representing the obscured SFRD is measured by combining individual bright detections along with an extrapolation of the infrared luminosity function to faint infrared luminosities (Zavala et al. 2021).Since the shape of the infrared luminosity function at such early times is relatively un-certain (in particular at the faint end, log 10 L IR /L ⊙ <11, where a number of studies report significant differences, e.g.Koprowski et al. 2017;Gruppioni et al. 2020;Fujimoto et al. 2023a;Barrufet et al. 2023a;Zavala et al. 2021;Traina et al. 2023), our data provides an opportunity to compare the contribution of populations below ALMA detection limits to that typically extrapolated from the luminosity function.
To make this comparison, we put our sample in context of existing far-infrared measurements at z > 3 and estimate their total infrared luminosity by integrating the maximum likelihood prospector model between restframe 8-1100µm.We then convert this to SFR using the indicator based on total infrared luminosity of Kennicutt & Evans (2012).We find that almost all of our sources are below the detection limits of earlier submm surveys at z > 3, with the majority of our sample having upper limits to L IR of 9 <log 10 L IR /L ⊙ < 11.8 (see Figure 7).We have only 4 detected sources with confirmed ALMA flux consistent with log 10 L IR /L ⊙ >12.Since many of our objects are non-detections (and thus upper limits), we are likely reaching the "extrapolation" regime (9 <log 10 L IR /L ⊙ < 11 Lsun) of the dustobscured SFRD at z > 3 (Zavala et al. 2021).
We focus more specifically on sources with log 10 L IR /L ⊙ <12, to make a direct comparison with the contribution from similar-luminosity populations in Zavala et al. (2021).That work determined that from 4 < z < 7, the fraction of the obscured SFRD contributed by log 10 L IR /L ⊙ <12 sources was relatively small (20%) and this fraction was flat with redshift.To compare, we also plot the SFRD contribution from our subset of sources with log 10 L IR /L ⊙ <12 (maroon points in Figure 8).We find that our measurements indicate that log 10 L IR /L ⊙ <12 sources make up a higher fraction of the obscured SFRD (26%) at 4 < z < 5, compared to 20% estimated by Zavala et al (a relatively minor factor of 1.3 increase).However, at z > 5, the fractional contribution from sources with log 10 L IR /L ⊙ <12 is much larger, and is comparable with the total obscured SFRD previously estimated.Compared to the earlier estimate of flat 20% fraction from log 10 L IR /L ⊙ <12, this suggests that the contribution of log 10 L IR /L ⊙ < 12 sources could be underestimated a factor of ∼ 5 at z > 5.In comparison to the total obscured contribution from Zavala et al. (2021, orange curve), a factor of 5 increase in the contribution from log 10 L IR /L ⊙ < 12 sources could double the obscured fraction of the SFRD at these redshifts.
While we note that major uncertainties exist in our measurements (see next Section), our data indicates that the overall census of dust obscured SF from earlier studies could be underestimated.Our findings are consistent with and similar to the estimates based on pre-JWST datasets that find higher SFRD at early times (Fujimoto et al. 2023a;Algera et al. 2023;Traina et al. 2023), in particular when including estimates based on existence of optically faint or or various "dark" sources (Williams et al. 2019;Gruppioni et al. 2020;Talia et al. 2021).
Based on our data, it may also be the case that the dust obscured star formation is underestimated at z > 6.We find that at z ∼ 7 the SFRD derived using the SEDmodeling is log 10 ρ SF R = −2.44 +0.12 −0.18 M ⊙ yr −1 Mpc −3 .This is consistent with Algera et al. 2023, but in excess of estimates in both Zavala et al. 2021 andBarrufet et al. 2023a.However, the fraction of our sources which are classified as candidate LRDs at these redshifts by color selection is high.While our analysis in Section 5.2.2 indicates that the NIRCam and MIRI photometry used to measure their stellar populations is not dominated by the light from the AGN, as a conservative estimate we also calculate the SFRD assuming the star formation contribution of these candidate AGN cannot be robustly determined.Thus we plot a second estimate with these LRD sources removed (open red circles).We find that in doing so, at z > 6 the obscured SFRD is substantially lower than implied by previous studies using NIRCam and HST data alone (Barrufet et al. 2023b), and that the obscured SFRD contributed by our own sample of dark galaxies would be substantially lower.We thus caution against overinterpretation of this measurement, and note that a complete assessment of the SFRD at z > 6 from dust obscured star formation using JWST selected samples will likely remain uncertain until we understand the true nature of LRDs.
6.1.1.Impact of modeling assumptions on the cosmic SFRD Despite the panchromatic data and limits used in our SED-fitting, our prospector-based SFRs may still be underestimating the intrinsic amount of SFR due to the model assumptions used to infer obscured star formation (namely, that the dust emission model assumed by prospector inherently prefer cold dust temperatures at the default prior settings).However, a wealth of evidence now points to hotter dust temperatures among compact, low metallicity systems that make up a higher fraction of the population at high-redshift (T dust ∼ 40-60 K; e.g.De Rossi et al. 2018;Sommovigo et al. 2020Sommovigo et al. , 2022;;Bakx et al. 2020;Schreiber et al. 2018;Faisst et al. 2017;Behrens et al. 2018).The effect of assuming colder dust temperature is to underestimate the SFR, since at fixed 1.1mm flux, the obscured SFR can be higher if dust is hotter.While we have adjusted our priors to allow higher dust temperatures (our prospector fits yielded a typical T dust ∼ 37K), this is still below some empirical constraints from dusty galaxies at high redshifts.
Thus, to test the impact of hotter dust on our estimates, we return to our ALMA image stacking that we used to estimate the average 1.1mm flux in the same redshift bins (see Section 5.2.1).Now including the detected sources in the stacks in order to assess the true average flux per bin, we find that in the above redshift bins, the stacked fluxes are 83±44, 213±39, 2±48, and 25±38 µJy/beam.However, to interpret the stacked 1.1 mm flux, we now assume a hotter dust SED template.To estimate the corresponding average L IR , we scale a far-infrared template for Haro-11, an analog for a high-redshift dusty SFG that factors in elevated dust temperatures that more realistically describe compact SFGs at z > 4 (T dust ∼ 47K, emissivity index β ∼ 1.9; Lyu et al. 2016).To obtain the L IR and SFR we integrate the template scaled to the ALMA stacked flux and integrate between restframe 8-1100µm and convert to SFR as earlier.
We find that our stacked ALMA fluxes are consistent with log 10 L IR /L ⊙ = 11.6, 12, < 10, < 11 for redshifts z ∼3.5, 4.5, 5.5, 7, which correspond to average SFRs that are 69±37, 152±28, 1±31, 15±23 M ⊙ /yr per redshift bin.For comparison, the average SFR per bin based on the prospector SED modeling is similar within 2σ (although systematically underestimated at z < 5).Using the number of galaxies per redshift bin, we translate these average estimates derived from the ALMA stack to an estimate of the total SFR per redshift bin contributed by this sample, scaled to the same signal-to-noise as reflected in the image stacked flux (see peach squares in Figure 8).For the highest two redshift bins where stacked flux is not detected, we instead plot 1σ upper limits.Generally, we find that at z < 5 the hotter dust template does increase our measured SFRD, while at z > 5 the stacks point to lower SFRD than found with the SED-modeling.More data sampling the far infrared SED to constrain the dust temperature would be needed to resolve the discrepancies between the different measurements.
One potential downstream impact of prospector not allowing for hotter dust in the modeling is that it likely forces a more quiescent (redder) stellar population solution (which, due to higher M/L ratios of older stars, could result in higher mass solutions at fixed luminosity).We note that redder stars is an easy solution to justify for the sources where we identified clear Balmer breaks indicative of older stellar populations.However, in cases without clear breaks, we may need dust emission modeling with hotter dust more typical of dusty galaxies observed at similar redshifts.However, we note the Figure 9.The stellar masses vs redshift of our sample of galaxies with MIRI coverage (34 sq arcmin area).For context we include the limiting stellar mass expected for 100% baryon conversion efficiency (solid line) and 20% efficiency (dashed line).For six galaxies between z ∼ 5.5 − 8 we measure stellar masses that likely exceed expectations for a small area, even after the stellar masses have decreased by up to an order of magnitude by including MIRI+ALMA data (see Figure 2).This comparison to the expectation curves from halo abundance and star formation efficiency show these measurements are likely still overestimated.
dust temperatures required may be even hotter (e.g.De Rossi et al. 2018), given the extremely compact nature of the objects in our z > 6 sample.
In fact this parameter space (young compact star forming regions driving hot dust temperature) is seen in the nearby Universe, and could potentially be analogous to these sources (e.g.Hainline et al. 2016).Thus, there is likely to be a significant bias at high redshift (where our sizes are also the most compact) that would mean that these estimates of SFRs are biased low for z > 5.This is because galaxies are more difficult to detect based on their 1.1 mm emission than they would be if local (colder) far infrared SEDs are used (see also Shivaei et al. 2022).Thus, our estimated contributions to the cosmic SFRD may also additionally be underestimated.Additional higher-frequency dust continuum imaging that would constrain the dust temperature would be needed to investigate this possibility further.

Stellar mass census: the abundance of massive galaxies
Despite the findings in Section 4 that MIRI+ALMA overall result in lower stellar masses, we still identify a remarkably large number of galaxies ( 16) above z > 3.5 with Log 10 M * /M ⊙ >10.The MIRI data also covers a relatively small area of 34 sq arcmin, making the identifi-cation of so many massive galaxies unlikely.In fact, 80% of our objects above z > 7 have Log 10 M * /M ⊙ >10.It has already been pointed out by Narayanan et al. (2023) that stellar masses are essentially unconstrained above z > 7 (and can dramatically over or underestimate stellar mass) owing to outshining of older stars by young low metallicity stars complicating the reconstruction of the SFH.However, we note that our high-redshift sources are not particularly young (as prospector preferred to model them with older stellar populations, either to allow low dust content to fit the ALMA data, or to accomodate clear Balmer breaks in a number of sources).However, even at redshifts below this problematic epoch noted by Narayanan et al. (2023), the number of high mass objects are also surprisingly, and problematically, high.
To demonstrate this, we show the stellar mass vs redshift for our MIRI+ALMA sample in Figure 9.For context, we also plot the expected stellar mass limit (i.e.mass where we only expect one halo, given our MIRI survey area), based on the halo mass function evolution with redshift.To estimate this we use the halo mass function calculator HMF published by Murray et al. (2013) and assume the halo mass function of Behroozi et al. (2013).We use the limiting halo mass to convert to a limiting expected stellar mass by assuming a fiducial baryon conversion efficiency into stars (0.2, dashed line), and for the limit where 100% of baryons converted in to stars (solid line).While our sample size is small, these curves help to indicate the number of sources whose mass is probably physically unlikely, under the typical assumptions of ΛCDM.We find one source at z = 5.5 and five sources at z ≳ 6.5 that are either improbably high mass, or, imply extremely efficient baryon conversion (>20%) compared to typical assumptions.
In Sections 2.3.3 and 4 we noted that a number of the photometric and spectroscopic solutions allowed revision to lower redshift upon detailed inspection (thus also allowing further reduction in the stellar masses even after inclusion of MIRI+ALMA data).However, in particular for high mass objects 90354, 121710, 200576 and 219000, we were unable to justify the possibility of lower redshift solutions.Further, while we identified a tentative emission line for 132229 (where the spectroscopic solution is ∆z ∼ 1 lower in redshift than the photometric redshift inferred with prospector), this galaxy still remains at a problematically high mass.We note that the last photometric candidate with unprobably high mass, 203749, has an uncertain redshift solution (we identified two comparable solutions at z=2.41 and the z ∼ 7 one we use for the analysis; we note that its LRD colors and unresolved morphology would make it an outlier among known z < 3 galaxies).Both redshift solutions are poor fits to the restframe UV photometry due to the presence of UV excess (Section 5.2.5).While we do not understand the nature of this object, we conclude it is likely to have alternative explanations besides a high-redshift massive object and do not consider it in the following discussion.
For the other 5 objects with less redshift ambiguity however (all are LRDs except 200576, and are relatively red) we find that under typical SED-modeling assumptions, these galaxies remain well above stellar mass expectations for their redshifts.None of them are well fit with the Sonora-Cholla brown dwarf atmospheric models, as their 2 -3µm colors are redder than what is observed in ultracool dwarfs.Three of these sources also meet the double-break criteria used to identify extremely massive candidates in Labbé et al. (2023a), and exhibit clear evidence of Balmer breaks (as discussed in Section 5.2.4).In fact, they all have very similar restframe SED shapes, with a strong balmer break, clear turnover at 1.6µm due to the stellar bump, and with deep ALMA limits, are best fit by moderately aged (mass weighted age ∼500 Gyr) and dusty (Av∼1) SEDs (see Figure 10).Due to our excellent wavelength sampling including at least four NIRCam medium bands and seven MIRI bands, it is extremely unlikely we would not be able to account for emission line boosting.All except 90543 exhibit significant MIRI detections even out to 10µm, making it less possible that we still overestimated the restframe near-infrared continuum due to having only upper limits from MIRI.
While there exists the possibility that these stellar masses are accurate, and these objects trace a population of galaxies that demonstrate highly efficient mass growth (e.g.Labbé et al. 2023a;Boylan-Kolchin 2023;Xiao et al. 2023a), we offer a few other possible explanations that would require additional data to disentangle.We note that five out of the six high mass sources are LRDs, and there remains large ambiguity about the nature of such objects.However, despite the evidence presented in Section 5.2.2 that the rest-optical SED is not dominated by a rising AGN continuum, it remains possible that the redder restframe optical SEDs of LRDs are (still) poorly described by stars alone.There is also the obvious possibility that various stellar population assumptions are incorrect at higher redshift (i.e.locallycalibrated models do not apply to high redshift phenomena where the stellar properties and environmental conditions may be dramatically different).Below we outline a few hypotheses based on the potential impact of AGN, and modeling assumptions, that could be tested with future high resolution spectroscopy.
Figure 10.The normalized restframe SEDs of the five most confident candidate massive galaxies in Figure 9 which include four LRDs, plus a Balmer break non-LRD, at z ∼ 5.5 − 8 (gray points; 203749 is excluded).All sources exhibit similar double-break SEDs exhibiting well constrained Lyman and Balmer breaks, evidence for a turnover at 1.6µm, and are consistent with high mass weighted age (∼ 500 Myr) and moderate dust (Av∼1).For comparison are a dustreddened AGN (yellow), a dusty star forming galaxy (teal; Polletta et al. 2007).While all SED types look very similar in the rest-frame optical (probed by HST+NIRCam) the SEDs diverge at λrest > 1µm, where MIRI shows a clear flattening due to the stellar bump, and ALMA rules out cold dust emission.

Sub-dominant AGN contamination
One possibility is that some sub-dominant (but perhaps still impactful) fraction of AGN flux in the restoptical and near-infrared is still driving up the stellar mass estimates of these four LRDs (ignoring for now ID 203749 with ambiguous redshift).We showed in Section 5.2.2 that the CIGALE modeling of the rest-optical SED of LRDs predicted much larger discrepancies with the MIRI data, primarily due to their redder NIRCam longwave photometry.The AGN model that overpredicted the MIRI flux included two components: the black body emission from the accretion disk in the restoptical (e.g.Big Blue Bump), plus hot+warm dust emission in the near and mid infrared.While we have shown that the strongly rising continuum from the hot+warm dust emission component does not agree with our MIRI data, this does not directly constrain the contribution from an accretion disk, unless the relative contributions are physically linked.However, some examples where the relative contribution of an accretion disk could be larger than we assumed (e.g.hot or warm dust deficient AGN; Lyu et al. 2017).In that case, the black body emission from an accretion disk could still contribute a small fraction of light to the rest-optical without being so visible in the restframe near-to-mid-infrared.This could be plausible (and could be an explanation given the high confirmation fraction of broad line AGN among LRDs with these colors; Greene et al. 2023) but we do not have the data or evidence to determine whether this is the case.AGN contribution from an accretion disk is not definitively ruled out by our relatively flat MIRI SEDs.This could in part explain the large inferred stellar masses.We are not aware of analogous sources at lower-redshifts, except for dust-deficient quasars.We note that while strong rest-optical emission lines driven by AGN can impact the interpretation of high stellar mass (Endsley et al. 2023;Kocevski et al. 2023), our inclusion of at least four NIRCam medium bands plus the longer-wavelength MIRI data to anchor the SED beyond the wavelength range with the most contamination is a strong mitigator of this uncertainty.
Future spectroscopy of this sample (in particular to measure restframe optical Hα, Hβ, [OIII] equivalent widths and line profile shapes) could potentially reveal the origin of the rest-optical continuum and validate the photometric measurements.In addition to confirming the likely presence of an AGN accretion disk via the broad lines, Greene et al. (2023) use the relatively low Hα EWs to argue that the continuum is not likely dominated by dust-obscured young stars.However, low EWs could also be explained by older stars and low level star formation over the past 10 Myr (which is consistent with our SED modeling results).Given that these sources typically have high attenuation (Av∼ 1 − 3), deeper spectroscopy than was obtained with FRESCO is likely required to adequately measure both broad and narrow line components (if they exist).High spectral resolution will also improve the differentiation between emission line broadening due to outflowing gas vs AGN.

Potentially errant modeling assumptions
We note that we use mostly conservative assumptions in our modeling, including a flexible-slope attenuation curve that allows extra attenuation in the UV to avoid "hiding" stellar mass with a fixed flat slope model (noting that the assumed attenuation curve slope can impact the recovered mass by almost an order of magnitude; Lo Faro et al. 2017;Williams et al. 2019).Additionally, while non-parametric SFHs have been shown to raise stellar mass by 0.3 dex on average based on representative galaxies at lower redshifts, (e.g.Leja et al. 2019) these objects would still remain systematically too high (unless the typical systematic offset is not representative for such extreme and red galaxies).
A number of works now report that the choice of prior on the SFH can significantly impact the ages and masses inferred by the non-parametric SFH (e.g.Leja et al. 2019;Lower et al. 2020;Tacchella et al. 2022b;Ji & Giavalisco 2022).Therefore, we test whether our choice of continuity prior for the SFH has weighted against bursty SFH solutions, which could result in lower mass solutions by enabling us to explain the SEDs with a larger fraction of stars at younger stellar ages.We re-run our prospector modeling for these six high mass sources instead using the Dirichlet prior (Leja et al. 2017), which allows for sudden and extreme changes in SFR in adjacent time bins (and is a weaker prior on the inferred shape of the SFH compared to the continuity).We find that the differences in stellar mass measured with the burstier Dirichlet prior does not cause a systematic shift in mass of our high mass sample, and further, are all consistent within the uncertainties of the stellar masses based on the continuity prior.We do find that the best SFH shapes do change with the Dirichlet solution and appear more stochastic, (indicating that robust stellar ages will require spectroscopy).This typical difference to the SFR in the most recent 30 Myr of only ∼1 M ⊙ , but we note that our most active source among the massive sample (219000) sees a 30% decrease to its recent SFR with the Dirichlet.However, our stellar masses, and thus primary conclusions, do not rely heavily on our prior choice.
Importantly, while the stellar masses from the Dirichlet prior are still consistent within their uncertainties with those inferred using the continuity prior (typical difference is < 0.1dex) we note that the uncertainties in mass derived from the posteriors of each set of modeling does not marginalize over these assumptions.Unfortunately, this is also the case for a number of assumptions that have gone into our modeling, including IMF shape (see e.g.Wang et al. 2023;Woodrum et al. 2023).Thus, the true uncertainty in stellar mass is larger than implied in Figure 9.We conclude that more advanced priors or an alternative IMF, among other assumptions, may be able to help lower the stellar masses and account for this tension.Otherwise, these targets are likely excellent testbeds for studying modeling systematics, or the potential for more exotic explanations in the future with near-infrared spectroscopy.

CONCLUSIONS
We study a sample of optically faint sources at z > 3 that are below the detection limits of the deepest HST and ALMA surveys to date, and have previously been missed from the galaxy census at 3 < z < 8.We find that these sources are relatively abundant within the JADES survey (66) and study a subset of those (29) for which deep multi-wavelength MIRI is also available.Our findings include: • The population of red optically faint galaxies are diverse in morphology (including both extremely extended sources as well as compact unresolved sources) and in SED shape, including sources resembling dust obscured star forming galaxies, post starburst galaxies, some objects exhibiting evidence of strong Balmer breaks despite high redshifts.
• We find that stellar population modeling for sources using HST+NIRCam data alone can result large masses and SFRs.When MIRI+ALMA are included, we find a median decrease in 0.6 dex in stellar mass and median decrease of 10×, for sources sources where HST+NIRCAM alone infer Log 10 M * /M ⊙ >10 and SFR>100 M ⊙ /yr, respectively.Thus, caution should be exercised when interpreting the SEDs of very red sources from HST+NIRCam data alone.
• Our sample includes ∼30% candidate AGN (selected as LRDs) and the fraction of LRDs is 100% among red galaxies above z > 6.5.Novel measurements with MIRI out to 25µm for this population confidently rule out that their very red rest-optical continuum primarily originates as obscured AGN continuum.Instead, evidence for a turnover in the SED between restframe 1-3µm suggests we are seeing the stellar bump, and the red rest-optical continuum is stellar in origin.We cannot rule out the presence of an AGN that becomes dominant in the restframe mid-infrared SED, which requires longer wavelength data.
• Noting that AGN are not likely dominating the restframe optical emission, we use our stellar population modeling to assess the contribution to the galaxy census of this previously-hidden population of galaxies.We estimate lower limits to the cosmic SFRD and find that galaxies at log 10 L IR /L ⊙ <12 and 4 < z < 6 may contribute 5× more to the obscured fraction of the SFR than previously estimated based on extrapolation of the infrared luminosity function, which could effectively double the obscured SFRD in this redshift range.
• We also assess the stellar masses we measure in the context of the limited area we probe in our survey, finding that five sources between z ∼ 5.5 − 8 have very high mass for our small survey area, despite the revision to lower stellar mass provided by the MIRI data.These sources have strong Balmer breaks and SED turnovers consistent with the stellar bump in the restframe near infrared, well described by moderately old and dusty SEDs (age∼500 Myr, Av∼1).We discuss plausible reasons for overestimated stellar masses, based on existing assumptions in stellar population and AGN modeling, motivating future work to characterize the physical properties of very red high-redshift galaxies.

Figure 2 .
Figure 2.For sources with MIRI and ALMA coverage, we compare inferred photometric redshift, stellar mass and SFR from our prospector modeling of HST+NIRCam+MIRI+ALMA vs just the HST+NIRCam photometry.Left panel shows the inferred photometric redshift (sources with spectroscopic redshifts are excluded from this panel).Middle panel shows the stellar mass, and right panel shows the SFR.At the very high mass and SFR end, we find that including the MIRI+ALMA data lowers the extreme masses and SFRs that are otherwise inferred (median decrease of 0.6 dex for sources above Log10M * /M⊙>10, and median 10× less for SFR > 100 M⊙/yr).This figure shows that stellar masses and SFRs can be wrong for extreme and red sources using NIRCam data without longer wavelength constraints.

0
Note-Properties of our H-dropout sample, measured including MIRI+ALMA data: (a) JADES DR1 ID (b) Photometric redshift measured by prospector, or from EAZY in the case of a spectroscopic redshift (c) Log10 of stellar mass in units M⊙ (d) SFR measured by prospector using the most recent 30 Myr time bin of the SFH in units M⊙/yr (e) V-band Attenuation (f) Mass weighted age in units Gyr (g) the ratio of bolometric luminosity from the galaxy divided by that from the AGN (h) flag (1) indicating the source meets our LRD color selection (j) confirmed in JADES NIRSpec data (k) confirmed in FRESCO data.

Figure 3 .
Figure3.Top panel: ALMA 1.1 mm stacks of our entire sample galaxies (including NIRCam-only sources, with individual detections removed) using the GOODS-ALMA data from(Franco et al. 2018), which covers our entire sample.Bottom panel: ALMA 1.1 mm stacks of a subset of our galaxies that are inside the footprint of the deeper ASAGAO imaging (individual detections removed).Only one source is inside the footprint at 3 < z < 4, so we instead show a combined 3 < z < 5 redshift bin.Our sample below z < 5 exhibit possible cold dust emission (2.2 − 3σ for both maps respectively), consistent with low luminosity DSFGs.Sources at z > 5 (including all LRDs) do not show detectable dust emission.Far right panels shows the ALMA stack of all 9 LRDs, which reaches a non-detection upper limit of 32µJy/beam.

Figure 4 .
Figure 4. UVJ diagram of our MIRI sample.Points are color coded by their inferred Av and symbol shape indicates old vs young mass weighted ages.LRDs are flagged with black dots.We omit the 1-σ uncertainties from the modeling posteriors for clarity, noting the uncertainties are large (sources have low S/N near both the rest frame U and J) enough to scatter away from the UVJ box.SEDs are consistent with our sample including primarily star forming sources, some with extreme dust obscuration (red in V-J), and some candidate older/post-starburst-like SEDs (Balmer breaks are also visible in the SEDs).

Figure 5 .
Figure 5. Top panel: Stacked MIRI images from 5.6-25µm of the LRD subset of H-drops (3.6" on a side).Bottom panel: stacked observed SED of the LRD subset, where NIRCam points are the median and interquartile range of the measured photometry, and MIRI points are measured from the median stacked imaging in top panel.The stacks are consistent with the flattening of the restframe near-infrared SED that was observed in the forced photometry of individual sources.They are not consistent with rising mid-infrared flux from a dust obscured AGN (e.g.Mrk 231 Polletta et al. 2007) and are a better match to stellar models of moderate dust and age (red).

Figure 6 .
Figure 6.NIRCam-MIRI color color diagrams for our sample (squares) with the LRD flagged (black dots).Objects are color-coded by redshift.Red star indicates the average colors obtained from a median-stacked LRD SED in Figure 5. Left: F277W-F444W vs F444W-F2100W colors of our sample.For comparison we plot a heavily obscured AGN template from Lyu & Rieke (2018, excluding stellar host) shown to match the red SED of so-called Extremely Red Quasars (ERQ) which at z ∼ 3 − 8 are extremely red in the restframe near-infrared.In comparison, LRDs have bluer or flat colors in F444W-F2100Wthan SEDs dominated by obscured AGN, more similar to SEDs dominated by stellar emission (tracing the stellar bump).Right: F277W-F444W vs F444W-F770W colors of our sample.Most LRDs are red or flat in F444W-F770W color, inconsistent with an interpretation that emission line boosting in F444W produces the red in F277W-F444W color (which would be blue in F444W-F770W).These colors are more consistent with flat or moderately-rising red continuum driven by starlight (the dashed-line box, specifically diagonal line, defines the color region at F444W-F770W>1.5 that excludes heavily reddened AGN continuum at z > 3;Akins et al. 2023).The NIRCam-MIRI colors in both plots favor SEDs that are stellar in origin.

Figure 8 .
Figure8.The cosmic SFRD of our full H-dropout sample in the MIRI footprint (red points).For comparison we show far-infrared measurements (obscured SFRD; orange shaded region) based on(Zavala et al. 2021), along with the UV-based compilation in that work (unobscured, uncorrected for dust; blue shaded region).We also plot the SFRD contribution from our sources with log10LIR/L⊙<12 (maroon points).Although our MIRI data suggests our LRD subset are stellar-dominated, we show the SFRD when removing them since their nature is uncertain (open red circles).We also include an estimate based on the average ALMA flux from stacking the 1.1 mm imaging scaling with a hotter dust template (peach squares).

Table 1 .
Properties of H-dropouts that are identified inside the MIRI SMILES footprint.

Table 2 .
Properties of H-dropouts without MIRI data.Spec-z Stellar Mass c SFR