COSMOS-Web: The over-abundance and physical nature of “little red dots”—Implications for early galaxy and SMBH assembly

Akins et al.

Hollis B. Akins

0000-0003-3596-8794]Hollis B. Akins
NSF Graduate Research Fellow
Department of Astronomy, The University of Texas at Austin, 2515 Speedway Blvd Stop C1400, Austin, TX 78712, USA

0000-0002-0930-6466]Caitlin M. Casey
Department of Astronomy, The University of Texas at Austin, 2515 Speedway Blvd Stop C1400, Austin, TX 78712, USA
Cosmic Dawn Center (DAWN), Denmark

0000-0003-3216-7190]Erini Lambrides
NASA Goddard Space Flight Center, 8800 Greenbelt Rd, Greenbelt, MD 20771, USA
NPP Fellow

0000-0001-9610-7950]Natalie Allen
Cosmic Dawn Center (DAWN), Denmark
Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200, Copenhagen N, Denmark

0000-0001-6102-9526]Irham T. Andika
Technical University of Munich, TUM School of Natural Sciences, Department of Physics, James-Franck-Str. 1, 85748 Garching, Germany
Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany

0000-0002-0245-6365]Malte Brinch
Cosmic Dawn Center (DAWN), Denmark
DTU-Space, Technical University of Denmark, Elektrovej 327, DK-2800 Kgs. Lyngby, Denmark

0000-0002-6184-9097]Jaclyn B. Champagne
Steward Observatory, University of Arizona, 933 N Cherry Ave, Tucson, AZ 85721, USA

0000-0003-3881-1397]Olivia Cooper
NSF Graduate Research Fellow
Department of Astronomy, The University of Texas at Austin, 2515 Speedway Blvd Stop C1400, Austin, TX 78712, USA

0000-0002-0786-7307]Xuheng Ding
Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan

0000-0003-4761-2197]Nicole E. Drakos
Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA

0000-0002-9382-9832]Andreas Faisst
Caltech/IPAC, 1200 E. California Blvd., Pasadena, CA 91125, USA

0000-0001-8519-1130]Steven L. Finkelstein
Department of Astronomy, The University of Texas at Austin, 2515 Speedway Blvd Stop C1400, Austin, TX 78712, USA

0000-0002-3560-8599]Maximilien Franco
Department of Astronomy, The University of Texas at Austin, 2515 Speedway Blvd Stop C1400, Austin, TX 78712, USA

0000-0002-7530-8857]Seiji Fujimoto
Hubble Fellow
Department of Astronomy, The University of Texas at Austin, 2515 Speedway Blvd Stop C1400, Austin, TX 78712, USA

0000-0002-8008-9871]Fabrizio Gentile
University of Bologna, Department of Physics and Astronomy (DIFA), via Gobetti 93/2, I-40129, Bologna, Italy
INAF – Osservatorio di Astrofisica e Scienza dello Spazio, via Gobetti 93/3, I-40129, Bologna, Italy

0000-0001-9885-4589]Steven Gillman
Cosmic Dawn Center (DAWN), Denmark
DTU-Space, Technical University of Denmark, Elektrovej 327, DK-2800 Kgs. Lyngby, Denmark

0000-0002-0236-919X]Ghassem Gozaliasl
Department of Physics, University of Helsinki, P.O. Box 64, FI-00014 Helsinki, Finland

0000-0003-0129-2079]Santosh Harish
Laboratory for Multiwavelength Astrophysics, School of Physics and Astronomy, Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, NY 14623, USA

0000-0003-4073-3236]Christopher C. Hayward
Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA
0000-0002-3301-3321]Michaela Hirschmann
Institute for Physics, Laboratory for Galaxy Evolution and Spectral modelling, Ecole Polytechnique Federale de Lausanne, Observatoire de Sauverny, Chemin Pegasi 51, 1290 Versoix, Switzerland
INAF, Osservatorio Astronomico di Trieste, Via Tiepolo 11, 34131 Trieste, Italy

0000-0002-7303-4397]Olivier Ilbert
Aix Marseille Université, CNRS, CNES, LAM, Marseille, France

0000-0001-9187-3605]Jeyhan S. Kartaltepe
Laboratory for Multiwavelength Astrophysics, School of Physics and Astronomy, Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, NY 14623, USA

0000-0002-8360-3880]Dale D. Kocevski
Department of Physics and Astronomy, Colby College, Waterville, ME 04901, USA

0000-0002-6610-2048]Anton M. Koekemoer
Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

0000-0002-5588-9156]Vasily Kokorev
Department of Astronomy, The University of Texas at Austin, 2515 Speedway Blvd Stop C1400, Austin, TX 78712, USA

0000-0001-9773-7479]Daizhong Liu
Max-Planck-Institut für Extraterrestrische Physik (MPE), Giessenbachstr. 1, D-85748 Garching, Germany

0000-0002-7530-8857]Arianna S. Long
Hubble Fellow
Department of Astronomy, The University of Texas at Austin, 2515 Speedway Blvd Stop C1400, Austin, TX 78712, USA

0000-0002-9489-7765]Henry Joy McCracken
Institut d’Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Université, 98 bis boulevard Arago, F-75014 Paris, France

0000-0002-6149-8178]Jed McKinney
Hubble Fellow
Department of Astronomy, The University of Texas at Austin, 2515 Speedway Blvd Stop C1400, Austin, TX 78712, USA

0000-0003-2984-6803]Masafusa Onoue
Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan
Center for Data-Driven Discovery, Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan
Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, People's Republic of China

0000-0003-2397-0360]Louise Paquereau
Institut d’Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Université, 98 bis boulevard Arago, F-75014 Paris, France

0000-0002-7093-7355]Alvio Renzini
INAF - Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, I-35122 Padova, Italy

0000-0002-4485-8549]Jason Rhodes
Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91001, USA

0000-0002-4271-0364]Brant E. Robertson
Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA

0000-0002-7087-0701]Marko Shuntov
Cosmic Dawn Center (DAWN), Denmark
Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200, Copenhagen N, Denmark

0000-0002-0000-6977]John D. Silverman
Department of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo, 113-0033, Japan
Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan
Center for Data-Driven Discovery, Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan
Center for Astrophysical Sciences, Department of Physics & Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA

0009-0003-4742-7060]Takumi S. Tanaka
Department of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo, 113-0033, Japan
Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan
Center for Data-Driven Discovery, Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan

0000-0003-3631-7176]Sune Toft
Cosmic Dawn Center (DAWN), Denmark
Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200, Copenhagen N, Denmark

0000-0002-3683-7297]Benny Trakhtenbrot
School of Physics and Astronomy, Tel Aviv University, Tel Aviv 69978, Israel

0000-0001-6477-4011]Francesco Valentino
European Southern Observatory, Karl-Schwarzschild-Strasse 2, D-85748, Garching bei München, Germany

0000-0002-7051-1100]Jorge Zavala
National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

(Affiliations can be found after the references)

has revealed a population of compact and extremely red galaxies at $z\gtrsim4$, which likely host active galactic nuclei (AGN). We present a sample of 434 “little red dots” (LRDs), selected from the 0.54 deg$^2$ COSMOS-Web survey. We fit galaxy and AGN SED models to derive redshifts and physical properties; the sample spans $z\sim5$–$9$ after removing brown dwarf contaminants. We consider two extreme physical scenarios: either LRDs are all AGN, and their continuum emission is dominated by the accretion disk, or they are all compact star-forming galaxies, and their continuum is dominated by stars. If LRDs are AGN-dominated, our sample exhibits bolometric luminosities $\sim10^{45-47}$ erg s$^{-1}$, spanning the gap between  AGN in the literature and bright, rare quasars. We derive a bolometric luminosity function (LF) $\sim100$ times the (UV-selected) quasar LF, implying a non-evolving black hole accretion density of $\sim10^{-4}$ M$_\odot$ yr$^{-1}$ Mpc$^{-3}$ from $z\sim2$–$9$. By contrast, if LRDs are dominated by star formation, we derive stellar masses $\sim10^{8.5-10}\,M_\odot$. MIRI/F770W is key to deriving accurate stellar masses; without it, we derive a mass function inconsistent with $\Lambda$CDM. The median stellar mass profile is broadly consistent with the maximal stellar mass surface densities seen in the nearby universe, though the most massive $\sim50$% of objects exceed this limit, requiring substantial AGN contribution to the continuum. Nevertheless, stacking all available X-ray, mid-IR, far-IR/sub-mm, and radio data yields non-detections. Whether dominated by dusty AGN, compact star-formation, or both, the high masses/luminosities and remarkable abundance of LRDs implies a dominant mode of early galaxy/SMBH growth.

§ INTRODUCTION

The launch of  has revealed new classes of objects unseen before; in particular, an abundant population of compact and extremely red objects, so-called “little red dots” <cit.>.
These objects are characterized by point-like morphology, red colors from observed-frame $\sim 2$–$5$ , and often exhibit flat or even blue colors from $1$–$2$ .
The population appears to be ubiquitous from $z\gtrsim 3$ <cit.> to $z\sim 9$ <cit.>, spanning several orders of magnitude in luminosity.
At these redshifts, the distinctive colors imply a blue UV slope with a steep/red optical continuum.
Initial studies interpreted these objects as high-mass ($\gtrsim 10^{10}\,\Msun$), $z\gtrsim 6$ galaxy candidates <cit.>, with either high equivalent width (EW) emission lines or strong Balmer breaks and likely with significant dust attenuation driving the red optical color.
Reconciling these high masses with the point-like morphology requires rapid starbursts in compact cores <cit.>, and implies a high stellar mass density in the early Universe, which for some time appeared to be in tension with expectations for the halo mass function in $\Lambda$CDM <cit.>.
While more recent mass estimates have brought down the exceedingly large stellar masses, relieving the possible tension with $\Lambda$CDM <cit.>, the sheer abundance of these objects seems to be at odds with a massive galaxy interpretation.

An alternative interpretation of these “little red dots” is that they are dominated by emission from an active galactic nucleus (AGN).
Indeed, spectroscopy has been obtained for $\gtrsim 50$ LRDs, in $\sim 80\%$ of cases confirming the high redshifts ($z>5$) and identifying broad Balmer lines <cit.>.
The point-like morphology and broad Balmer emission lines are signatures of unobscured (Type I) AGN; i.e., we have a direct view to the broad-line region (BLR) and the accretion disk, but with foreground dust attenuation either from a dusty interstellar medium (ISM) or polar dust, nuclear to the AGN <cit.>.
This configuration is similar to the “red quasars” at $1<z<3$ <cit.>, for which the dust attenuation is thought to originate from the galaxy-scale ISM <cit.>, and which have also been found out to high-$z$ <cit.>.
Although at face value, the LRDs appear to populate the faint end of the AGN UV luminosity function, when accounting for dust obscuration, their intrinsic luminosities can be quite large <cit.>.
Consistent with these high luminoisities, many authors have derived large black hole masses relative to the host galaxy mass, implying exotic black hole seed/growth mechanisms <cit.>.
However, the origin of the red continuum emission—whether it is dominated by AGN or host galaxy—remains unclear, but has major implications for our interpretation of these objects as massive galaxies, luminous quasars, or some combination thereof.

There exists some evidence that the continuum emission from the LRDs is actually dominated by stars.
For one, multi-band /MIRI observations from $5$–$25$  have shown a remarkably flat SED in the rest-frame mid-IR for a number of LRDs, placing strong upper limits on the contribution from hot dust in the AGN torus, in tension with the implied luminosities of the quasars <cit.>.
Moreover, several broad-line objects exhibit clear Balmer breaks in their spectra, implying a dominant contribution from evolved stars in the rest-frame optical, though the precise contribution is degenerate, with stellar mass estimates ranging over 2 orders of magnitude <cit.>.
Even still, not all LRDs with spectroscopic coverage exhibit clear broad lines <cit.>, perhaps indicating that the population may not be homogeneous.
It remains unclear whether the LRDs can be explained entirely by AGN, with unique geometry/dust properties, or also represent a significant population of compact/massive galaxies in the early Universe.

In this work, we present a large sample of 434 $z\gtrsim 5$ LRDs selected from the COSMOS-Web  survey.
With the large sample size, nearly doubling the number of known LRDs <cit.>, and large on-sky area,  minimizing the impact of cosmic variance, we provide strong constraints on the volume density of this population.
In Section <ref> we describe the COSMOS-Web survey, our data reduction methodology, and ancillary multi-wavelength data.
In Section <ref> we describe our sample selection and SED fitting methodology; we fit each object to both galaxy and QSO models to characterize the two alternative scenarios.
In Section <ref> we consider the interpretation of the LRDs as AGN; in particular, we derive their contribution to the bolometric luminosity function and black hole accretion density.
In Section <ref>, we consider the interpretation of the LRDs as instead compact star-forming galaxies; in particular, we present the derived stellar mass function and examine the feasibility of such dense stellar systems.
In Section <ref> we stack the available multi-wavelngth data over the entire sample, providing constraints on the panchromatic SED of the typical LRD.
In Section <ref> we discuss this sample in context with other  results, with particular focus on the physical processes that may be responsible for the unique characteristics of the LRD population.
Throughout this paper, we adopt a <cit.> initial mass function and a cosmology consistent with the <cit.> results ($H_0=67.66$ km s$^{-1}$ Mpc$^{-1}$, $\Omega_{m,0} = 0.31$).
All magnitudes are quoted in the AB system <cit.>.

Color-magnitude and color-color diagrams illustrating the selection of LRDs in this paper. Left: The $m_{277}-m_{444}$ vs. $m_{444}$ color-magnitude diagram. We select all sources with $S/N_{444} > 12$ (roughly  $m_{444}$ brighter than $27.2$) and $m_{277}-m_{444}>1.5$. Dark red points indicate the subsample which is covered by MIRI imaging, while light red points indicate those with only NIRCam coverage. COS-756434, which is spectroscopically-confirmed at $z=6.99$, is marked with a star. Blue points indicate objects for which a brown dwarf solution is preferred via SED fitting.Center: The $m_{277}-m_{444}$ color vs. compactness $C_{444}$, defined as the ratio of the F444W fluxes in $0\farcs2$ and $0\farcs5$ diameter apertures. Objects with $C_{444}<0.5$ are generally extended, while those with $C_{444}>0.7$ are generally imaging artifacts, i.e. hot pixels. Right: The $m_{277}-m_{444}$ vs. $m_{150}-m_{277}$ color-color diagram. While we do not select for blue SW colors (either $m_{115}-m_{150}$ or $m_{150}-m_{277}$), as the imaging is too shallow to constrain SW colors at high SNR, we demonstrate here that the NIRCam colors can serve to filter out brown dwarfs, which tend to have bluer SW colors than galaxies at temperatures $T_{\rm BD} < 1200$ K.

§ DATA

§.§ COSMOS-Web

COSMOS-Web is a large  Cycle 1 treasury program imaging a contiguous 0.54 deg$^2$ with NIRCam and 0.19 deg$^2$ with MIRI <cit.> in the COSMOS field <cit.>.

The COSMOS-Web imaging includes four NIRCam filters (F115W, F150W, F277W, and F444W) at approximate $5\sigma$ depths of $27.4$–$28.2$ AB mag, and one MIRI filter (F770W) at a $5\sigma$ depth of 26 AB mag.
The full details on the NIRCam and MIRI reduction process will be presented in upcoming papers (M. Franco et al. in prep and S. Harish et al. in prep, respectively) but are briefly described here.
The raw NIRCam imaging was reduced by the JWST Calibration Pipeline version 1.12.1, with the addition of several custom modifications (as has also been done for other JWST studies, e.g. ) including the subtraction of $1/f$ noise and sky background.
We use the Calibration Reference Data System (CRDS)[<jwst-crds.stsci.edu>] pmap 1170 which corresponds to NIRCam instrument mapping imap 0273.
The final mosaics are created in Stage 3 of the pipeline with a pixel size of 003/pixel.
Astrometric calibration is conducted via the / alignment tool (JHAT; ), with a reference catalog based on an /F814W 003/pixel mosaic in the COSMOS field <cit.> with astrometry tied to Gaia-EDR3 <cit.>.
The median offset in RA and Dec between our reference catalog and the NIRCam mosaic is less than 5 mas.
The MIRI/F770W observations were reduced using version 1.8.4 of the
JWST Calibration pipeline, along with additional steps for background subtraction that was necessary to mitigate instrumental effects.
The resulting mosaic was resampled onto a common output grid with a pixel scale of 003/pixel and aligned with ancillary /F814W imaging of the region.

§.§ PRIMER-COSMOS

The Public Release Imaging for Extragalactic Research (PRIMER) survey (P.I. J. Dunlop, GO#1837) is a large Cycle 1 Treasury Program to image two HST CANDELS Legacy Fields (COSMOS and UDS) with NIRCam+MIRI <cit.>.
The PRIMER-COSMOS field comprises $\sim 130$ sq. arcmin and is contained entirely within the COSMOS-Web footprint (see  for details on the field layout), providing additional depth and filter coverage over $\sim 8\%$ of the field.
PRIMER imaging includes 8 NIRCam bands (adding F090W, F200W, F356W, and F410M to the COSMOS-Web filter set), plus two MIRI bands (adding F1800W).
We have reduced the PRIMER data following the same procedure as with the COSMOS-Web data.

§.§ Ancillary Multi-wavelength Data

While the COSMOS field is covered with substantial ground-based (Subaru and UltraVISTA) and Spitzer/IRAC imaging, we do not use these additional photometric data in our SED fitting.
The lower-resolution of these maps, relative to the HST+JWST imaging, mandates careful deblending, which can suffer from significant uncertainties.
However, the deep ground-based imaging of the COSMOS field can be useful to constrain photometric redshifts and rule out degenerate low-redshift solutions.

Even beyond the optical/near-infrared, the COSMOS field is rich with ancillary multiwavelength data in the X-ray, MIR, FIR/submm, and radio.
We make use of Chandra X-ray imaging from the COSMOS Legacy Survey <cit.>, Spitzer/MIPS 24+70 µm imaging from S-COSMOS <cit.>, Herschel/PACS+SPIRE imaging from PEP and HerMES <cit.>, SCUBA-2 850 µm imaging from S2COSMOS <cit.>, ALMA 1.2mm imaging from CHAMPS (over 0.18 deg$^2$; Faisst et al. in prep), ALMA 2 mm imaging from MORA <cit.>, and radio imaging from MeerKAT at 1.4 GHz <cit.> and VLA at 3 GHz <cit.>.

§.§ Catalog construction

Details on the photometric catalogs for COSMOS-Web will be presented in a future paper (M. Shuntov, L. Paquereau et al. in prep.); here, we provide a brief summary of the photometric catalog used in this work.
We conduct source detection and aperture photometry using SEP <cit.>, the python implementation of SourceExtractor <cit.>.
We perform source detection on a $\chi^2$ detection image constructed from all four NIRCam bands <cit.>.
First, F115W, F150W, and F277W are PSF-homogenized to F444W and converted to noise-equalized images by multiplying by the square root of the weight map.
We then fit the negative tail of the pixel distribution to a Gaussian function to measure the per-pixel rms of the noise-equalized map.
The maps are normalized by this rms and truncated at negative values to avoid false positives arising from multiple colocated negative noise peaks.
The detection image is then computed by adding these PSF-homogenized, truncated signal-to-noise maps in quadrature.
The resulting $\sqrt{\chi^2}$ image can be directly related to the probability of an individual pixel value being drawn from the sky noise distribution <cit.>.

We adopt a detection threshold of $\sqrt{\chi^2} = 3.698$ (with no background subtraction), equivalent to a 2.33$\sigma$ detection, or a 1% probability of drawing a given pixel value from the noise distribution.
We then perform aperture photometry on the available /ACS, /NIRCam, and /MIRI imaging.
Similar to <cit.>, we utilize small elliptical apertures  (with a Kron factor $k=1.1$) and correct to the total flux based on the default Kron aperture flux in F444W ($k=2.5$).
ACS and NIRCam images are PSF-matched to F444W (as with the detection stage), and photometry is corrected for the F444W PSF flux falling outside of the large Kron aperture.
For the MIRI bands, we adopt the same aperture sizes and apply a correction for the fractional flux loss due to the larger PSF.
In particular, we compare the aperture flux in the native resolution F444W image to the flux after PSF-matching F444W to the MIRI F770W PSF, and apply this correction to the measured MIRI flux.

Photometric uncertainties are adjusted using a random aperture method to account for correlated noise in the mosaics.
For each band, and for a range of aperture sizes, we place 200,000 random circular apertures across the field, avoiding areas with detected sources, and measure the standard deviation of the measured fluxes by fitting the negative tail of the Gaussian distribution.
Similar to <cit.> and <cit.>, we fit a power law of the form $\sigma_N = \alpha N^{\beta/2}$, where $N$ is the number of pixels in the aperture.
All photometric uncertainties are then added in quadrature with this random aperture noise measurement $\sigma_N$ for the corresponding aperture size.

§ SAMPLE SELECTION & CHARACTERIZATION

From the COSMOS-Web data, we construct our sample of LRDs, fit them each to a series of SED models, and characterize the typical properties and volume density of the population.
Before describing our sample selection and SED fitting methodology, we present a brief note of clarification on the motivation behind the strategy we adopt.
While, ideally, we would be able to fit each object with a composite galaxy/AGN SED model, decomposing the two components, the limited information available (largely due to the lack of spectroscopic redshifts) makes this kind of analysis prone to overfitting.
Instead, we would like to limit the number of free parameters;
rather than attempting to decompose the SEDs into galaxy/AGN components, we explore two extreme scenarios:

* QSO models: the LRD population is primarily composed of dust-reddened AGN, with the optical continuum dominated by the direct thermal emission from the accretion disk ($\ll 1$ pc scales, temperatures $\gtrsim 10^5$ K) <cit.>. Note that in this scenario, the host galaxy may dominate in the rest-frame UV, but is negligible in the rest-frame optical.

* Galaxy models: the LRD population is primarily composed of compact/dusty starbursts, with the optical continuum dominated by the emission from young stellar populations <cit.>. Note that this scenario does not preclude the presence of AGN in the population, but just requires that the red continuum emission in the optical is dominated by the stellar component.

These two scenarios may be thought of as “edge cases”; while neither likely represents the true nature of the entire population, exploring the implications of either scenario yields valuable physical insight.
In both cases, we assume that the characteristic blue UV slope of the LRDs is simply leakage or scattering through a clumpy obscuring medium, similar to what is seen in some lower-redshift DSFGs <cit.> and red quasars <cit.>, and we note that we do not select directly for this blue component.

§.§ Sample selection

Figure <ref> outlines our color-compactness sample selection procedure.
From our master photometric catalog described in <ref>, we first downselect our catalog to objects which appear compact in F444W, defined via the compactness metric

\begin{equation}
C_{444} = f_{444}(d=0\farcs2)/f_{444}(d=0\farcs5)
\end{equation}

We use F444W, despite the larger PSF relative to shorter-wavelength bands, since the LRDs are brightest in this band (and it is sometimes the only band with sufficient S/N to yield a reliable measure of compactness).
Through a series of source injection simulations, we found that point sources reliably have $0.5<C_{444}<0.7$ provided the signal-to-noise is $\gtrsim 20$.
We select objects with $C_{444} > 0.5$, meaning that $\gtrsim 50\%$ of the F444W broadband emission within a $0\farcs5$ diameter aperture is contained within a $0\farcs2$ diameter aperture.
We exclude objects with $C_{444} > 0.7$, which tend to be imaging artifacts (hot pixels), which don't follow the curve-of-growth associated with the PSF.
In Section <ref> we fit these objects to 2D profiles to verify their compact/point-like nature.

We then select objects with red F277W–F444W colors,

\begin{equation}
m_{277}-m_{444} > 1.5.
\end{equation}

This is the same threshold that was adopted by <cit.>, who found that almost all objects satisfying this color threshold appear point-like and blue from $1$–$2$ µm.
Similarly, while <cit.> selected for broad-line AGN based on the H$\alpha$ emission, the majority of their objects show similarly red colors from $2$ to $\sim 3$–$5$ µm.
We note that this color threshold is more stringent than several recent works which select LRDs via F277W$-$F444W $> 0.7$ or F200W$-$F356W $>1.0$ <cit.>.
We focus on the reddest subset of the LRD population, which may bias our sample towards more dust-obscured and/or higher-redshift objects.
Figure <ref> shows our sample compared to other literature samples in color-redshift space; here, we adopt the redshifts from QSO SED models as described in section <ref>.
We also show in Figure <ref> the F277W$-$F444W color vs. redshift derived by redshifting the high SNR NIRSpec/PRISM spectra of two lower-redshift LRDs (J0647_1045 at $z=4.532$ from <cit.> and RUBIES-BLAGN-1 at $z=3.104$ from <cit.>).
Indeed, by selecting for the reddest objects, we bias our sample to $z\gtrsim 5$, where F277W$-$F444W covers the rest-frame optical range.
Moreover, the somewhat bimodal redshift distribution is a result of either  (from $z\sim 5$–$7$) or  $\lambda 5007$ ($z\sim 7$–$9$) falling into F444W and contributing to the extremely red color.

Comparison of our sample to literature LRD samples in F277W$-$F444W color vs. redshift space. Our sample is shown in red, as in Figure <ref>, and literature photometric/spectroscopic samples are shown in grey open/filled points, respectively <cit.>. The blue and green points/lines highlight two lower-redshift LRDs with high-SNR PRISM spectra (J0647_1045 at $z=4.5$ from  and RUBIES-BLAGN-1 at $z=3.1$ from ) and their F277W$-$F444W color when redshifted.
The dark orange line shows a redshifted SDSS extreme emission line galaxy (EELG) template <cit.>.
We see that selecting for the reddest objects, ($m_{277}-m_{444}>1.5$, as indicated by the grey line) biases our sample to $z\gtrsim 5$ due to the shifting rest-frame wavelengths probed by F277W and F444W.
The bimodal redshift distribution of our sample is a result of either  (from $z\sim 5$–$7$) or  $\lambda 5007$ ($z\sim 7$–$9$) falling into F444W and contributing to the extremely red color, though the selection for $m_{277}-m_{444}>1.5$ mitigates contamination from intrinsically blue, strong emission line objects.

We note that our selection is perhaps better described as a selection for extremely red objects <cit.>, as LRDs are often selected based on the composite blue+red spectral shape <cit.>.
However, as the physical origin of the blue UV component is unclear, requiring it in selection may bias the sample towards objects with unobscured host galaxies or unique dust geometries.
Because of this, and due to the shallow depth of the F115W+F150W imaging in COSMOS-Web ($5\sigma \sim 27.1$–$27.6$ AB mag), we opt to use only the red F277W–F444W criterion.
We also note that other LRD studies <cit.> often require red colors in multiple filter pairs to rule out contamination from extreme emission line galaxies (EELGs) which may exhibit a very red F277W–F444W color despite a blue continuum.
While our selection has the potential to yield such contaminants, we note that by focusing on the reddest objects ($\sim 0.5$ mag redder than other selections), we largely mitigate this effect; in fact, $\sim 90\%$ of our sample that is covered by the PRIMER survey ($N=37$) satisfy the multi-band color selection described in <cit.>.
If we instead adopted a less stringent single color criteria, F277W$-$F444W$>0.7$, this fraction drops to $65\%$.

The simple selection criteria described above yields a sample of 533 sources across the 0.54 sq. degree COSMOS-Web field.
Of these, 163 sources are covered with MIRI F770W imaging (32 of which also have F1800W imaging from PRIMER), and 39 are covered with PRIMER NIRCam imaging.

§.§ Removal of Likely Brown Dwarfs

In order to remove contaminanting brown dwarf stars in the Milky Way <cit.>, we fit each candidate with a grid of low-temperature stellar atmosphere models.
We construct a grid of models that spans a large parameter space, including cloudy, chemical equilibrium models <cit.>, cloud-free models in both chemical equilibrium (Sonora-Bobcat, $T\sim 200$–$1300$ K; ) and disequilibrium (ElfOwl, $T\sim 300$–$1000$ K; ), and low-metallicity models (LOWZ, $[{\rm M}/{\rm H}]=-1$; ).
We explore a range of temperatures $T\sim 200$–$2400$ K, surface gravities $g=100$ and $3160$, and metallicities $[{\rm M}/{\rm H}] = -1$, $-0.5$, and $0$.
We convolve the model SEDs with the HST+JWST filter curves and perform a simple grid-fitting routine, scaling the fluxes of each model to minimize the $\chi^2$.
We adopt the model with the lowest $\chi^2$ as the best-fitting stellar model, and compare with the best-fit galaxy/quasar SEDs described later in Section <ref>.
104 objects in our sample have brown dwarf model $\chi^2$ less than the minimum galaxy/quasar model $\chi^2$ and are removed.
The fraction of sources estimated to be brown dwarfs is $\sim 9\%$ for the MIRI subsample, and $\sim 24\%$ for the NIRCam-only sample.
This difference is due to the presence of strong absorption bands in cold brown dwarfs at $\sim 8$ µm, allowing MIRI/F770W imaging to break degeneracies between galaxy/QSO and brown dwarf models, particularly that our NIRCam data have only four bands.

After the brown dwarf removal procedure, our final LRD sample contains 434 sources, 148 of which are covered with MIRI F770W imaging and 29 of which are covered with PRIMER NIRCam imaging.

§.§ Morphological Fitting

We perform 2D image fitting to the /NIRCam imaging to validate the simplified “compactness” selection and determine whether any sources are marginally resolved.
We do this using galsim, a galaxy image forward-modeling code designed to produce realistic simulated images based on analytical profiles <cit.>.
We opt for this approach, rather than the commonly-used image fitting tools (such as galfit; ), in order to conduct Bayesian inference; specifically, we utilize the nested sampling algorithim MultiNest <cit.> to fully explore the parameter space and obtain robust uncertainties on the derived sizes for each object.
We fit the F444W imaging for each object to a PSF-convolved 2D Sérsic profile with the effective radius allowed to vary from $0\farcs001$ to $0\farcs3$, the Sérsic index $n$ from $0.5$ to $4$, and the axis ratio $q$ from $0.7$ to $1$.
Figure <ref> shows the measured F444W effective radius vs. the F444W magnitude for our sample.
The half-width at half-maximum for the F444W PSF is shown in the black dotted line, as an approximate resolution limit for the observations.
We can constrain the sizes even below the PSF size, for high SNR sources, as even slightly extended objects would show some residuals after PSF subtraction.
To quantify this effect, we performed our Sérsic profile fitting to PSF stars identified via their position in the half-light radius vs. magnitude plane.
To avoid including marginally resolved sources in this selection, at $m_{444}>24$, we adopt the brown dwarf sample (identified as contaminants in Section <ref>) to estimate the limiting size.
We then fit a curve to the 90th percentile of the sizes measured in bins of F444W magnitude, deriving effective resolution limit shown in the dashed line in Figure <ref>.

Nearly all of our candidates are unresolved, which is to be expected, as we are selecting for compact sources based on the F444W aperture photometry.
The unresolved nature of the LRDs places an upper limit on their physical sizes, which ranges from $R_{\rm eff} \lesssim 0\farcs01$ at $m_{444} \approx 22$ to $\lesssim 0\farcs04$ at $m_{444} \approx 27$.
At $z=6$, this corresponds to effective radii $\lesssim 100$–$300$ pc.
A few sources are marginally resolved, particularly the brightest subset ($m_{\rm 444} \lesssim 24$), perhaps indicating some host galaxy contribution.
We keep these marginally resolved sources in the sample, as they satisfy the original compactness criteria and have effective radii $\lesssim 200$ pc.

Effective radius vs. F444W magnitude for our sample. The dotted line indicates the half-width at half-maximum (HWHM) for the F444W PSF. The dashed line indicates the effective resolution limit based on the measured sizes of PSF stars. While the majority of our sample is unresolved, a few sources appear marginally resolved, perhaps indicating host galaxy contribution.

§.§ Spectral Energy Distribution Fitting

We fit both galaxy and QSO models to all objects in our sample.
However, we do not intend to use these results to classify objects as either galaxy or AGN-dominated based on goodness-of-fit—there are, in general more free parameters in the model than we have data points.
Instead, as we are testing the validity of either extreme scenario (the “edge cases”), intend to examine the results from the model fitting under the assumption of either scenario.

692866_sed.pdf 845096_sed.pdf
Cutouts and SEDs for six particularly bright ($m_{444} \lesssim 24$) LRDs. We show 1.5” square cutouts in /ACS F814W, /NIRCam F115W, F150W, F277W, and F444W, and /MIRI F770W (where available).
Measured photometry is shown in red, or 2$\sigma$ upper limits for non-detections; the best-fitting QSO model is shown in black,  in blue, and brown dwarf in light grey.
The inset panels show the redshift probability distribution with the maximum likelihood value marked with a dashed line.
The complete figure set (434 images) is available at <https://github.com/hollisakins/akins24_cw>.

Same as Figure <ref>, but for four bright LRDs with MIRI F1800W coverage from PRIMER. Two sources (above) are only in COSMOS-Web NIRCam coverage, but fall within the PRIMER MIRI coverage; two sources (below) have both PRIMER NIRCam and MIRI coverage. In all cases, the MIRI F1800W imaging yields marginal or non-detections, consistent with a flat SED in $f_\nu$ from rest-frame $\sim 1$–$3$ µm.

§.§.§ QSO models

We first fit each LRD in the sample to Type I AGN models.
In particular, we adopt broad-line AGN templates, consistent with the high fraction of broad-line AGN confirmed in LRDs observed with  spectroscopy <cit.>.
While existing SED fitting codes such as CIGALE <cit.> or prospector <cit.> do include AGN models, these are primarily designed to decompose host galaxy/AGN components, particularly in the rest-frame mid-IR, and thus are not ideal for this study.
Moreover, these codes do not typically include strong/broad emission lines from the AGN, which may contribute significantly to the observed photometry.
We therefore employ a custom SED fitting code, assuming that the emission is dominated by AGN, and employing flexibile dust attenuation and emission line models.

We model the intrinsic accretion disk emission as a single power law with a slope $\beta$ ($f_\lambda \propto \lambda^\beta$).
The normalization is set by the intrinsic UV magnitude $M_{\rm UV,int}$, prior to any dust attenuation. We allow $M_{\rm UV,int}$ to vary from $-18$ to $-27$ and $\beta$ to vary from $-1.5$ to $-2.8$, both with flat priors.

We employ the <cit.> QSO emission line template with , corresponding to the average emission line template at $z=2$, though we employ a few adjustments.
In particular, we allow the overall emission line scaling to vary from $1$ to $2$ times the nominal value, which is calibrated to a reference QSO continuum.
This is motivated by the higher broad-line EWs observed for -selected AGN compared to local samples <cit.>.
We also disable Lyman-$\alpha$ emission due to the opacity of the IGM at the redshifts relevant to our study (though we note that several LRDs exhibit Ly$\alpha$ emission, even at $z>7$, e.g. ).
Even at $z<7$, this choice shouldn't impact the results significantly, given the relatively low EWs typical of Ly$\alpha$ emission <cit.> compared to the filter bandwidths ($\sim 3000$ Å).

We additionally use a flexible dust attenuation prescription to account for the blue UV slope characteristic of this population.
We adopt the QSO attenuation law from <cit.>, which is similar to the SMC extinction curve (except at $\lambda \lesssim 1700$ Å, where it is significantly shallower).
We add an additional parameter, a “scattering” fraction, or a fraction of the intrinsic spectrum that escapes unattenuated. We allow this scattering fraction to vary from $10^{-5}$ to $10^{-1}$ with a log-uniform prior.
The effect of this is similar to employing an even greyer attenuation law, but adds a bit more

flexibility and aids in the interpretation of the blue slope as coming from scattering around the AGN torus <cit.> or leakage through holes in the obscuring medium, perhaps created by outflows <cit.>.
We additionally include a parameter, $\eta_{\rm NLR}$, representing the scaling of the dust attenuation in the narrow line region (NLR).
While the BLR and accretion disk are subject to the nominal $A_V$, the more spatially-extended NLR is attenuated by $\eta_{\rm NLR} A_V$.
We adopt a fixed $\eta_{\rm NLR} = 0.3$, consistent with the LRD Balmer decrement measurements from <cit.>.

We fit these QSO models using the MultiNest nested sampling algorithm <cit.>.
As a test, we apply our SED fitting code to the reported /NIRCam photometry for the two confirmed $z>7$ broad-line AGN reported in <cit.> and <cit.>.
We recover the spectroscopic redshifts within $1\sigma$, and our estimate of the bolometric luminosity (from the continuum) is consistent with the estimate derived from the H$\beta$ line for the source in <cit.>.
We note that we overestimate the bolometric luminosity for the source in <cit.> by a factor of $\sim 5$, though this is consistent with the continuum-based $L_{\rm bol}$ estimate from <cit.>.

§.§.§ Galaxy models

To test the other “edge case,” that the population is dominated by compact star-forming galaxies, we fit each object to a galaxy SED model.
Here we assume that all the light originates from the galaxy, including the blue/UV and red/optical components.
We use bagpipes, a Bayesian SED fitting code <cit.> designed to fit galaxy photometry and spectra.
We use the BPASS library of stellar SED models <cit.> and a non-parametric SFH model.
The SFH is parametrized by the $\Delta \log({\rm SFR})$ in adjacent time bins; we adopt the continuity prior <cit.>, i.e. the prior on $\Delta\log({\rm SFR})$ is a $t$-distribution with $\sigma=0.3$ and $\nu=2$ degrees of freedom.
We adopt five age bins for the SFH, with the first three at fixed ages from $0$–$10$ Myr, $10$–$50$ Myr, and $50$–$200$ Myr and two more logarithmically-spaced from 200 Myr to $z=20$.

We adopt log-uniform priors on the stellar mass (from $10^{6}$ to $10^{13}\,M_\odot$) and metallicity (from $10^{-3}$ to $1.5\,Z_\odot$).
Nebular emission is implemented via cloudy, specifically the updated grids presented in <cit.>, with $\log U_{\rm ion}$ allowed to vary from $-4$ to $-1$.
We adopt an SMC dust law with $A_V$ allowed to vary from 0 to 6; to be consistent with our QSO modeling, we modify the dust attenuation prescription to include a similar “scattering fraction” parameter, which we allow to vary with a log-uniform prior from $10^{-5}$ to $10^{-1}$.
We note, however, that in the galaxy case, this is better interpreted as representing “leakage” through holes in the ISM dust screen, as is found in some lower-redshift DSFGs <cit.>.

§.§ Individual Objects

We highlight a few individual objects which stand out from our sample.
Figure <ref> shows the cutouts and SEDs for six particularly bright LRDs ($m_{444} \sim 22$–$24$).
In each panel, the best-fit QSO model is shown in black, and the best-fit bagpipes model is shown in blue.
For sources which fall in the PRIMER coverage, we plot the PRIMER photometry as well.
For the brightest LRDs in our sample, we derive bolometric luminosities $10^{46-47.5}$ erg s$^{-1}$ from the QSO modeling, or stellar masses $10^{10-12}$ from the bagpipes modeling.
This is due to their steep, red optical continuum slopes, anchored by bright F444W/F770W magnitudes, implying significant dust obscuration.

/NIRSpec PRISM spectrum of COS-756434. The source is clearly detected in continuum and exhibits strong +H$\beta$ and H$\alpha$ emission, confirming the redshift at $z=6.9993\pm 0.0001$, consistent with the photometric redshift $z_{\rm phot} \sim 7.05^{+0.09}_{-0.05}$.
Broad components are visible in both the H$\beta$ and H$\alpha$ lines, with a FWHM of $2750\pm 40$ km s$^{-1}$, a clear signature of AGN activity.
The top panels highlight the broad+narrow line decomposition for H$\beta$ and H$\alpha$, as well as the deblending of $\lambda 4363$ and H$\gamma$.
The source also exhibits several high-ionization emission lines, including [Ne iv]$\lambda 2423$ and [Ne iii]$\lambda 3869$, consistent with the AGN interpretation.

Several bright LRDs in our sample also fall into the PRIMER MIRI coverage, which includes F1800W.
We highlight four of these objects in Figure <ref>.
The only object in our sample detected in F1800W is COS-105829, which is also the brightest LRD in the F1800W coverage ($m_{770} \sim 22.7$) but is not in the PRIMER NIRCam coverage.
No other source is detected, despite several (pictured in Figure <ref>) being similarly bright in F444W and F770W.
The MIRI limits imply flat SED (in $f_\nu$) between observed-frame $7.7$–$18$ µm, or rest-frame $\sim 1$–$3$ µm at the redshifts of these objects.
This may indicate marginal contribution from a hot dust torus, which we return to in Section <ref>.

We provide the SED-derived properties of our sample in Table <ref>.

§.§ Spectroscopic Confirmation of COS-756434

We identify one source in our sample, COS-756434, with public /NIRSpec observations.
COS-756434 was observed with $\sim 6$ hr PRISM spectroscopy as part of a DDT program (#6585, P.I. D. Coulter) targeting high-redshift supernovae.
We reduce the NIRSpec data using the standard  pipeline (version 1.14.0); the automatically extracted 1D spectrum is shown in Figure <ref>.
The spectrum exhibits strong  $\lambda\lambda4959,5007$$+$H$\beta$ and H$\alpha$ emission, confirming the redshift at $z_{\rm spec}=6.9993$, consistent with the photometric redshift of $z_{\rm phot} \sim 7.05^{+0.09}_{-0.05}$.
We fit the spectrum with a simple model including narrow+broad H$\beta$ and H$\alpha$ emission, narrow emission for other detected lines, and a piecewise power-law continuum.
We detect the H$\gamma$ and  $\lambda4363$ lines, $[{\rm Ne}\,\textsc{iii}]\,\lambda\lambda 3869,3967$ and even the high-ionization line $[{\rm Ne}\,\textsc{iv}]\,\lambda2423$, consistent with the AGN interpretation.
We measure a broad-line FWHM $= 2920^{+40}_{-50}$ km/s, corresponding to a black hole mass $\log M_{\rm BH} \sim 7.8 \pm 0.5$ using the calibration from <cit.>.
Though we do not discuss the spectrum of this object in detail, we report the properties measured from the spectrum in Table <ref> and adopt its spectroscopic redshift in the remainder of the paper.

Spectroscopic measurements for COS-756434
Property     Units     Value
$z_{\rm spec}$     …    6.9993±0.0001

R.A.     hms     10:00:25.6561

Decl.     dms     +02:21:36.165

$f_{\lambda,5100}$     $10^{-21}$ erg s$^{-1}$ cm$^{-2}$ Å$^{-1}$     3.20_-0.03^+0.04

$A_{V{\rm ,BLR}}$     AB mag     2.9^+0.2_-0.2

$\log M_{\rm BH}$     $M_{\odot}$     7.8 ±0.5

$\log L_{{\rm bol,H}\alpha}$     erg s$^{-1}$     45.2 ±0.4

$\lambda_{\rm Edd}$     …    0.19^+0.30_-0.11

FWHM$_{\rm narrow}$     km s$^{-1}$     494^+5_-11

FWHM$_{\rm broad}$     km s$^{-1}$     2920^+40_-50

$F_{[{\rm O}\,\textsc{iii}]\lambda4363}$     $10^{-20}$ erg s$^{-1}$ cm$^{-2}$     7.07_-0.67^+0.63

$F_{{\rm H}\beta{\rm ,narrow}}$     $10^{-19}$ erg s$^{-1}$ cm$^{-2}$     1.18_-0.13^+0.13

$F_{{\rm H}\beta{\rm ,broad}}$     $10^{-19}$ erg s$^{-1}$ cm$^{-2}$     0.23_-0.02^+0.02

$F_{[{\rm O}\,\textsc{iii}]\lambda4959}$     $10^{-19}$ erg s$^{-1}$ cm$^{-2}$     2.77_-0.07^+0.07

$F_{[{\rm O}\,\textsc{iii}]\lambda5007}$     $10^{-19}$ erg s$^{-1}$ cm$^{-2}$     7.98_-0.08^+0.08

$F_{{\rm H}\alpha{\rm ,narrow}}$     $10^{-18}$ erg s$^{-1}$ cm$^{-2}$     1.11_-0.03^+0.02

$F_{{\rm H}\alpha{\rm ,broad}}$     $10^{-18}$ erg s$^{-1}$ cm$^{-2}$     2.44_-0.03^+0.03

§ TWO INTERPRETATIONS OF LRDS

We now consider the two alternative interpretations of LRDs as either reddened quasars/dusty AGN or compact star-forming galaxies.
In the following sections we explore the implications of either scenario, in particular the resulting mass/luminosity functions.

§.§ Dusty AGN/reddened quasars

First, we discuss the interpretation of the LRDs as a population of reddened AGN, as has been posited <cit.>.
We characterize the LRDs in terms of their bolometric luminosity, i.e. the total AGN luminosity after correcting for dust attenuation; we consider the bolometric luminosity rather than the UV luminosity as these sources are quite obscured.
We compute the bolometric luminosity from the intrinsic model SED (i.e., before any dust attenuation), using the monochromatic luminosity at $3000$ Å and a bolometric correction of 5.15 <cit.>.

Figure <ref> shows the COSMOS-Web LRDs in the AGN bolometric luminosity vs. redshift plane.
Literature AGN with  spectroscopy are shown in purple (reddened AGN) and blue (broad-line AGN).
We also show a compilation of bright quasars from the literature <cit.>.
The COSMOS-Web LRDs span the range between the literature  samples and the bright quasars, with a handful reaching similar bolometric luminosities $\sim 10^{47}$ erg s$^{-1}$.
These luminosities are generally much higher than other LRD samples, due to the wide/shallow nature of COSMOS-Web, but we note that similarly luminous LRDs have been observed at $z\sim 3$–$5$ (e.g. UNCOVER-45924 at $z=4.66$; , and RUBIES-BLAGN-1 at $z=3.1$; ).
However, the high luminosities and abundance of these objects is surprising given the relative scarcity of luminous quasars.

Bolometric luminosity vs. redshift.
AGN from the literature with  spectroscopy are shown in purple <cit.> and blue <cit.>. Our sample from COSMOS-Web spans the range of luminosities between the literature  AGN and the UV-bright quasars discovered in wide-field ground-based surveys <cit.>.

To quantify the excess of luminous AGN in our sample, we compute the bolometric luminosity function.
We note that here, and in our estimates of the stellar mass function in  Section <ref>, we compute the volume of the survey as the differential comoving volume scaled to the survey area of 0.54 deg$^2$ and integrated over the redshift bin.
In order to capture the large uncertainties on the redshifts and luminosities of our objects, we marginalize over the joint posterior distribution of $L_{\rm bol}$ and redshift for each source; that is, individual objects can fractionally span multiple redshift and luminosity bins.
We additionally marginalize over the bin size, computing the luminosity function in bins centered at $\log L_{\rm bol} \sim 43$, $43.2$, $\dots$, $48$ with bin widths randomly drawn from $0.3$ to $1.0$ dex.
This avoids biases imparted by the choice of bin size.

We additionally employ a simple completeness correction based on our photometric catalog and source selection procedure.
In short, we inject 10,000 mock sources, constructed to be compact and extremely red, into the NIRCam data and repeat our cataloging procedure and selection.
We then estimate the completeness as a function of $m_{444}$, for a color distribution representative of our full sample.
We translate this to be a function of $L_{\rm bol}$ using the correlation between $m_{444}$ and the median $L_{\rm bol}$ derived from our SED fitting for each source.
Uncertainty on the completeness is computed from the same correlation using the 16th and 84th percentile $L_{\rm bol}$; this is added in quadrature with the Poisson uncertainty from the luminosity function.

AGN bolometric luminosity function at $z\sim 5$–$7$ (top) and $z\sim 7$–$9$ (bottom). Our results are shown in red, for the MIRI subsample (solid) and the full sample (dashed). The bolometric luminosity function limits from <cit.> and <cit.>, based on spectroscopically-confirmed LRDs, are shown in purple/blue. Photometric estimates from <cit.>, derived from other  imaging surveys, are shown in dark blue. The $z\sim 6$ quasar luminosity function from <cit.> is shown in grey. Our results imply a factor of $\sim 100$ overabundance of luminous AGN relative to pre- estimates from UV-bright quasars.

Binned AGN bolometric luminosity function
2*Redshift     $\log L_{\rm bol}$     2c$\log \phi~[\#\,{\rm dex}^{-1}\,{\rm Mpc}^{-3}]$

$[{\rm erg\,s}^{-1}]$     Full sample     MIRI subsample
5–7     45     -3.99^+1.31_-0.81     -3.89^+1.33_-0.84

5–7     45.5     -4.39^+0.65_-0.34     -4.24^+0.64_-0.35

5–7     46     -4.54^+0.31_-0.18     -4.50^+0.32_-0.18

5–7     46.5     -4.80^+0.22_-0.14     -4.84^+0.23_-0.16

5–7     47     -5.32^+0.19_-0.16     -5.39^+0.25_-0.28

7–9     45.5     -4.19^+1.30_-0.59     -4.17^+1.33_-0.65

7–9     46     -4.42^+0.58_-0.28     -4.42^+0.58_-0.28

7–9     46.5     -4.67^+0.26_-0.17     -4.73^+0.27_-0.20

7–9     47     -5.32^+0.29_-0.35     -5.32^+0.29_-0.35

The black hole accretion rate density as a function of redshift. Our results are shown in red, for our full sample (light red, smaller points) and our MIRI subsample (dark red, larger points).
Results from X-ray AGN surveys <cit.> are shown in black, including extrapolation at $z>5$.
Results from /MIRI are shown in green <cit.>.
If the red dots are universally AGN, our results imply a non-evolving black hole accretion density of $\sim 0.5$–$1\times 10^{-4}$ M$_\odot$ yr$^{-1}$ Mpc$^{-3}$ from $z\sim 1$ to $z\sim 9$.

Figure <ref> shows the AGN bolometric luminosity function inferred from our photometric candidate sample at $z\sim 5$–$7$ (top) and $z\sim 7$–$9$ (bottom).
Our sample is shown in red, with the solid line corresponding to the MIRI subsample.
We also show lower limits on the bolometric LF from LRDs selected in spectroscopic surveys (UNCOVER, , and FRESCO, ) and estimates from other photometric surveys <cit.>.
Our results are consistent with spectroscopic samples, and suggest an overall  excess of red AGN relative to the quasar bolometric luminosity function measured from rest-UV selected quasars (shown in grey, ).
While the faint-end slope of the quasar bolometric luminosity function remains uncertain, our results suggest that the knee of the LF could be dominated by these dust-reddened AGN, even up to $L_{\rm bol} \sim 10^{47}$ erg s$^{-1}$.
In particular, we find a factor of $\sim 100$ over-abundance of the LRDs relative to UV-bright quasars, at fixed $L_{\rm bol}$.
This could also be explained if we were overestimating the bolometric luminosities of the LRDs, though we note that this would also require a factor of $\sim 100$ in luminosity to be consistent with the QLF.
Our estimated bolometric luminosity function is tabulated in 0.5 dex bins of luminosity in Table <ref>

How does the abundance of AGN at $z\gtrsim 5$ compare to the population at lower redshift?
Figure <ref> shows the evolution of the black hole accretion rate density (BHAD) with redshift.
Results from X-ray AGN surveys <cit.> are shown in black, extrapolated above $z\sim 5$.
We also show the SFRD <cit.> as black-hole galaxy co-evolution scenarios would suggest that the BHAD should follow a similar evolution to the SFRD.
Results from /MIRI observations of $z\sim 1$–$5$ obscured AGN in CEERS are shown in blue <cit.>.
We compute the BHAD by integrating the bolometric luminosity function and assuming a radiative efficiency factor $\epsilon = 0.1$.
As Figure <ref> shows, if the LRD are universally AGN, our results imply a non-evolving black hole accretion density of $\sim  10^{-4}$ M$_\odot$ yr$^{-1}$ Mpc$^{-3}$ from $z\sim 1$ to $z\sim 9$.
This is an order of magnitude above the scaled SFRD, and several orders of magnitude above the X-ray extrapolations.

§.§ Compact SFGs

Next we consider the interpretation of the “little red dots” as compact dust-obscured star-forming galaxies, as has been suggested in the literature <cit.>.
Under this interpretation, the continuum emission would originate primarily from starlight, with the red color either coming from an old stellar population with a strong Balmer break or high EW emission lines combined with strong dust attenuation.
The lack of mid-IR detections for these sources may further support their star-forming nature <cit.>, and the relative abundance is more consistent with galaxies than quasars.

Figure <ref> shows the derived stellar mass vs. redshift for our LRD sample.
We indicate the region of parameter space that would be disfavored or disallowed in a $\Lambda$CDM framework in dark grey, corresponding to the stellar mass of the most massive halo in a given area assuming the cosmic baryon fraction of $f_b = 0.158$ and stellar baryon fraction $\epsilon = 1$.
Here we compute the halo mass function using the  python package <cit.> assuming an SMT fitting function <cit.>.
We plot the relevant threshold for the full sky, the COSMOS-Web area ($0.5$ deg$^2$) and the MIRI area ($0.17$ deg$^2$) in successively lighter shades of grey.
We note that a more reasonable upper limit on the observed stellar mass would use $\epsilon \approx 0.2$, the maximum inferred from the peak of the stellar mass to halo mass relation across a range of redshifts <cit.>.

A few of our LRDs exceed the maximum thresholds expected/allowed in $\Lambda$CDM, particularly at $z\gtrsim 7$.
This may be seen as reason to favor the AGN interpretation: if they are dominated by stellar light, they are simply too massive.
However, we note that the subsample with MIRI coverage does not violate these thresholds; the stellar masses derived for the MIRI subsample are generally lower than those for the NIRCam-only sample.
To test this further, we perform a series of  runs on the MIRI sample, both including and omitting the F770W data point.
We find systematic offsets of $\sim 1$–$2$ dex, towards larger masses, when not including the MIRI data.
Similar results have been found in other surveys with MIRI imaging (e.g. CEERS, ; SMILES, ).
In our case, this effect is largely due to the impact of strong emission lines boosting the photometry in the F444W bands.
With the additional data point on the red end, the modeling tends to favor younger stellar populations with stronger emission lines and lower masses—i.e., a compact starburst.
We note, however, that the total contribution of the LRDs to the cosmic star-formation rate density is still $\lesssim 5\%$, based on the sample with MIRI data.
We therefore conclude that the LRDs do not individually violate $\Lambda$CDM limits, but that MIRI imaging is critical to derive accurate masses, especially in the case of limited broadband coverage.

Stellar mass vs. redshift from galaxy SED fitting with .
The dashed lines indicate the maximum stellar mass we would expect to find in a given volume based on the halo mass function (assuming a global star-formation efficiency $\epsilon=1$, or $\epsilon=0.2$ for a more realistic assumption).
We plot the stellar masses inferred for “little red dots” in JADES from <cit.>, as well as the three ultra-massive objects identified in <cit.>.

Stellar mass functions at $z\sim 5$–$7$ (top) and $z\sim 7$–$9$ (bottom). As in Figure <ref>, our sample is shown in red, for the MIRI subsample (solid) and the full sample (dashed). The blue and green lines show the pre- estimates of the SMF from <cit.> and <cit.>. The grey shaded region shows the parameter space prohibited by $\Lambda$CDM, scaling the the halo mass function by the cosmic baryon fraction $f_b$ and a maximal “efficiency” $\epsilon=1$. The grey dashed line shows a more realistic efficiency $\epsilon=0.2$. While the SMF for the full sample violates the limits imposed by $\Lambda$CDM, particularly at $z>7$, the impact of MIRI on the stellar masses brings the SMF into agreement with literature estimates.

However, even without individual objects violating $\Lambda$CDM, the ubiquity of these sources suggests a high number density of massive galaxies.
Figure <ref> shows the stellar mass function (SMF) derived from our LRD sample.
We show the pre- estimates of the SMF from CANDELS <cit.> and COSMOS <cit.> in blue and green, respectively.
As discussed previously, our full sample shows a strong excess relative to pre- estimates, and with several candidates in the region prohibited by $\Lambda$CDM.
However, the subsample with MIRI coverage does not show as strong of an excess, and in fact shows a bright-end slope more consistent with the steep slope derived in previous studies.
This suggests that these objects could be dominated by star-formation.
Importantly, though, this implies that the massive galaxy population at $z\gtrsim 5$ is dominated by these compact starbursts, which has important implications for our understanding of early stellar mass growth.

A stronger constraint can be placed on the stellar content of these objects based on their point-like morphology.
There exists a maximum stellar mass surface density observed in dense stellar systems <cit.>, which remains consistent across $z\sim 0$–$3$ and over $\sim 9$ dex in stellar mass.
In other words, local and $z\lesssim 3$ stellar systems show similar maximum stellar densities $\sim 10^{11}$ M$_\odot$ kpc$^{-2}$ in their cores.
Figure <ref> shows the stellar mass profiles derived for our sample, in comparison to the average profiles for various local and high-redshift populations from <cit.>.
Here we include only the MIRI subsample, for which we can derive reliable stellar masses.
We find a median central stellar density of $\log \Sigma_\star \sim 11.6$, higher than—but broadly consistent with—that observed in the lower-redshift universe.
Approximately half of the sample, however, has stellar mass profiles significantly in excess of this limit.
On the one hand, these central densities may be reduced over time due to $N$-body relaxation, the timescale for which can be $\lesssim 1$ Gyr for particularly compact stellar systems <cit.>.
On the other hand, this may imply that the stellar masses are simply overestimated, requiring significant AGN contribution to the continuum to reconcile the derived central densities.

Median stellar mass surface density profile for the COSMOS-Web LRDs, compared to the compilation of dense stellar systems in <cit.>. Individual profiles are shown in the thin, light red lines. The different colors/symbols correspond to the different classes of objects compiled by <cit.>, ranging from local globular clusters/nuclear disks to massive ellipticals at $z\gtrsim 2$; all exhibit similar central densities, roughly consistent with the median LRD in our sample.

§ MULTIWAVELENGTH STACKING

Median-stacked panchromatic SED for the COSMOS-Web LRDs. Top: Stacked images in all / bands available, as well as Chandra X-ray, Spitzer/MIPS 24 µm, JCMT/SCUBA-2 850 µm, ALMA 1.2 and 2.0 mm imaging, and MeerKAT/VLA 1.3 and 3GHz imaging. Each stack is annotated with the number of objects contributing to the stack. For /, we show in light red the apertures used to extract photometry directly from the stacks. For the remaining bands, we show a 5” square to guide the eye given the varying resolution of the images.
Bottom: full panchromatic SED for the median COSMOS-Web LRD. Photometry measured from stacks (or $5\sigma$ upper limits from the non-detections) are shown in black.
We plot the maximal allowed SED model in red.
The faint blue lines in the X-ray and UV/optical show the intrinsic (i.e., unattenuated) SEDs for the corona and the accretion disk.
We plot two alternative X-ray SEDs consistent with the data, a Compton-thin model with an anomalously low $\alpha_{\rm OX}$ (solid) and a Compton-thick model with a normal $\alpha_{\rm OX}$ given the optical luminosity (dashed).
In the far-IR, we show several single-temperature blackbodies at different normalizations, ranging from to cold ISM dust the hottest possible dust, at the sublimation temperature.
The maximal SED model includes a total dust luminosity of $2\times 10^{12}\,L_{\odot}$ with a Gaussian temperature distribution $\sim 200\pm 75$ K.
Finally, in the radio, we adopt a slope $\alpha=-1.5$ consistent with radio-loud AGN <cit.>.
For reference, we also plot the best-fit model SEDs for a $z=5.85$ DSFG (MAMBO-9; Akins et al. in prep) and a $z=6.85$ obscured AGN <cit.>.

The COSMOS field is rich in multiwavelength coverage from the X-ray, MIR, FIR/submm, and radio.
As described in Section <ref>, we inspect cutouts around each source in Chandra X-ray imaging, Spitzer/MIPS 24 µm, SCUBA-2 and ALMA mm/sub-mm imaging, and VLA/MeerKAT radio data.
No source in our sample is individually detected in any of these data; this is somewhat surprising, given the large sample size and inclusion of several remarkably luminous objects in the sample.
In particular, under the AGN interpretation, objects with $L_{\rm bol} \gtrsim 10^{46}$ erg/s at $z\gtrsim 4$ are regularly detected in the X-ray <cit.>, mid-infrared (from the dusty torus), and radio <cit.>.

To push the multiwavelength constraints deeper, we stack the data for the entire sample.
Even without individual detections, if these objects exhibit ubiquitous, faint flux, it would show up in stacking.
We perform a simple median stack of the cutouts in each map.
We note that for Chandra we compute the sum, rather than the median, as the median is always 0 given the low X-ray counts.
We additionally perform stacking on Spitzer/MIPS 70 µm and Herschel imaging from PACS <cit.> and SPIRE <cit.>; while SPIRE imaging is severely confusion-limited, stacking on maps with a mode equal to 0 removes flux boosting from adjacent sources <cit.>.
Due to the widely varying PSF size of the various datasets, we stack relatively large cutouts (ranging from $5''$ for / to $10'$ for Herschel/SCUBA data).

The stacked cutouts are shown in the top panels of Figure <ref>.
No flux is detected in our stacks in any band aside from /.
We measure fluxes directly from the stacked cutouts for +, yielding the median (observed-frame) SED for the sample of LRDs; note that we find marginal flux in F1800W ($2.9\sigma$), but formally take this as an upper limit.

To quantify the flux limits implied by the non-detections in the various multi-wavelength datasets, we plot in the lower panel of Figure <ref> the median panchromatic SED for the COSMOS-Web LRDs.
We derive the upper limits for each band in the following manner.
For the ALMA data, which comes from the CHAMPS (Faisst et al. in prep) and MORA surveys <cit.>, and SCUBA-2, we add in quadrature the RMS at the position of each source, and divide by $N$.
For all other maps, which have relatively uniform depth, we adopt the reported depths, scaled down by $\sqrt{N}$.
For Chandra, we use the limiting depths reported in <cit.> for the $0.5$–$2$ and $2$–$7$ keV bands.
For Spitzer/MIPS and Herschel/PACS and SPIRE we use the depths reported in <cit.>, <cit.>, and <cit.>, respectively.

Do the stacked limits provide useful constraints on the physical properties of the LRDs?
To answer this, we show in the dark red line in Figure <ref> the maximal SED allowed by the limits.
The UV/optical/NIR SED is constructed from our QSO model (see <ref>) to match the photometry, with $A_V=3$ and $f_{\rm scat} = 1\%$ at $z=6$.
The intrinsic (unattenuated) SED is shown in light blue, with a UV magnitude of $M_{\rm UV} = -21.8$ and a 3000 Å continuum luminosity of $L_{3000} = 5.5\times10^{44}$ erg/s.
Note that these properties are very similar to the average SED-derived properties for the sample as a whole.

For the X-ray SED, we adopt a photon index of $\Gamma=1.8$ and an upper energy cut at 300 keV <cit.>.
The normalization is set by the (unattenuated) optical SED, using the conversion parameter $\alpha_{\rm OX} \equiv -0.3838 \log (L_{\nu,2500}/L_{\rm 2\,keV})$.
We then attenuate the SED assuming a column density $N_{\rm H}$, using the photoelectric absorption cross sections provided in <cit.>.
We show in the darker, solid line the X-ray SED assuming Compton-thin absorption $\log N_{\rm H}/{\rm cm}^{-2} = 23$.
This column density is consistent with the moderate dust attenuation ($A_V\approx 3$) based on empirical relations between dust attenuation and gas absorption in red quasars <cit.>.
With this column density, the X-ray limits require $\alpha_{\rm OX} < -2$; by contrast, we expect $\alpha_{\rm OX}\approx -1.45$ ($\pm 0.1$) from the $\alpha_{\rm OX}$–$L_{2500}$ relation <cit.>.
An alternative model adopts $\alpha_{\rm OX} = -1.45$ but a much higher column density $\log N_{\rm H}/{\rm cm}^{-2} > 24.5$.
In this Compton-thick model, the X-ray SED falls just below the limits for the shallower hard-band, and is entirely attenuated in the soft band.
Several other authors have reported X-ray measurements for LRDs that imply lower masses/luminosities than the broad lines in the optical <cit.>, which may indicate a fundamental difference in the nuclear structure in these objects.

Turning to the far-IR, the stacked limits provide a constraint on the dust SED of the LRDs.
We plot four single-temperature modified blackbody models, from cold ISM dust ($T=30$ K) to the hottest dust in the AGN torus, at the sublimation temperature ($T=1500$ K).
The ALMA/SCUBA constraints limit the cold dust component to $\lesssim 10^{11}\,$L$_\odot$, and the MIRI/MIPS constraints limit the hottest dust components to $\lesssim 10^{10}$ L$_\odot$.
The shallow Herschel limits, however, provide little constraints on warm dust $\sim 100$–$200$ K.
We include in the dark red model SED a dust component with a Gaussian distribution of temperatures, with a mean of 200 K and a standard deviation of 75 K.
We do not impose any energy balance with the rest-frame optical, instead adopting a luminosity of $2\times 10^{12}$ L$_\odot$, which falls just below the stacked $5\sigma$ limits from Herschel.
With the available limits, in the maximal SED model, the vast majority of the bolometric luminosity could be emitted in the far-IR and still result in non-detections.
We note that the bolometric luminosity inferred from $L_{3000}$, $L_{\rm bol} \approx 7.4\times 10^{11}$ L$_\odot$ <cit.> is still less than the maximal luminosity allowed in the FIR, though only by a factor of $\sim 3$—the current FIR limits are consistent with a standard bolometric correction for LRDs.

Finally, in the radio, we adopt a power-law slope of $\alpha = -1.5$, consistent with radio-loud AGN at high-redshift <cit.>.
The normalization is set by the luminosity density at 1.4 GHz, $L_{\nu,1.4\,{\rm GHz}} = 5.2\times 10^{23}$ W Hz$^{-1}$.
This is quite a bit lower than we would expect for radio-loud AGN; in terms of the ratio of the ratio of the 1.4 GHz luminosity to the optical $B$-band luminosity, $R^*_{1.4} = L_{\nu,1.4\,{\rm GHz}}/L_{\nu,4400} = 30$ defines the boundary between radio-loud and radio-quiet AGN <cit.>.
Our radio limits constrain $R^*_{1.4}<10$, at the $5\sigma$ level, suggesting that the average LRD is inconsistent with a radio-loud AGN.

Observed-frame stacking results
Band     Wavelength     $N$     Units     $F_\nu$
Chandra 2–10 keV$^\dagger$     0.28 nm     434     erg s$^{-1}$ cm$^{-2}$     $<7.2\times 10^{-17}$

Chandra 0.5–2 keV$^\dagger$     1.24 nm     434     erg s$^{-1}$ cm$^{-2}$     $<1.1\times 10^{-17}$

/ACS F606W     0.60 µm     29     nJy     $<3.85$

/ACS F814W     0.81 µm     434     nJy     4.11 ±0.52

/NIRCam F090W     0.90 µm     29     nJy     7.37 ±1.22

/NIRCam F115W     1.16 µm     434     nJy     15.76 ±0.62

/NIRCam F150W      1.50 µm     434     nJy     16.46 ±0.51

/NIRCam F200W     2.00 µm     29     nJy     17.66 ±0.77

/NIRCam F277W     2.77 µm     434     nJy     30.51 ±0.23

/NIRCam F356W      3.58 µm     29     nJy     83.30 ±0.56

/NIRCam F410M      4.08 µm     29     nJy     113.56 ±1.03

/NIRCam F444W     4.42 µm     434     nJy     173.00 ±0.33

/MIRI F770W      7.66 µm     148     nJy     277.99 ±2.94

/MIRI F1800W     18.00 µm     24     nJy     (435.2 ±150.9)

Spitzer/MIPS     24 µm     434     µJy     $<3.43 $

Spitzer/MIPS     70 µm     434     mJy     $<0.36 $

Herschel/PACS     100 µm     434     mJy     $<0.40$

Herschel/PACS     160 µm     434     mJy     $<0.82$

Herschel/SPIRE     250 µm     434     mJy     $<0.65$

Herschel/SPIRE     350 µm     434     mJy     $<0.86$

Herschel/SPIRE     500 µm     434     mJy     $<1.24$

SCUBA-2     850 µm     434     µJy     $<239$

ALMA     1.2 mm     30     µJy     $<157$

ALMA     2.0 mm     128     µJy     $<35.4 $

VLA 3 GHz     10 cm     434     uJy     $<0.56 $

MeerKAT 1.28 GHz     23 cm     434     uJy     $<1.80$
$^\dagger$For Chandra, we quote the absolute flux limit (rather than flux density) to be consistent with X-ray surveys. Uncertainties are quoted as $1\sigma$; for non-detections, we quote the $5\sigma$ upper limits.

§ DISCUSSION

Our results do not necessarily favor one interpretation over the other for the origin of the LRDs, and present challenges to our understanding of early galaxy/SMBH growth in either case.
Under the AGN interpretation, our estimate of the bolometric luminosity function is consistent with the lower limits provided by recent spectroscopic observations, with a $\sim 2$ dex excess relative to blue quasars across $\log L_{\rm bol}\sim 45$–$47$.
Under the galaxy interpretation, our stellar mass function estimate is consistent with $\Lambda$CDM limits, but implies that the LRDs would dominate the high-mass end of the SMF, and have central stellar mass densities significantly greater than the maximum observed at lower-redshift.
In reality, neither interpretation is likely 100% true; both galaxy and AGN components will contribute to the SED.
Moreover, the mechanisms that trigger compact nuclear starbursts (i.e. those in local LIRGs or URLIRGs) are thought to be the same mechanisms that trigger quasars; specifically, gas-rich major mergers or disk instabilities <cit.>.
Distinguishing between these two possibilities will inevitably require more data, both to examine sub-populations and decompose individual SEDs.

§.§ What dominates the continuum?

The major remaining uncertainty surrounding the LRDs is whether the red optical continuum emission is dominated by stellar light or the AGN accretion disk.
This directly impacts the measured stellar masses and AGN bolometric luminosities, and, to a lesser degree, the black hole masses.
On the one hand, the broad line EWs measured for LRDs with  spectroscopy have been found to be largely consistent with local calibrations, i.e. the continuum is consistent with being dominated by the accretion disk <cit.>.
However, it is possible that the broad lines are unusually strong relative to the accretion disk continuum in these objects.
This could be due to different AGN physics (such as a higher BLR covering factor, e.g. ), differential attenuation between a spatially-extended BLR and the accretion disk <cit.>, or kinematically broadened lines in a multiphase ISM, even without an AGN (e.g. N. Roy et al. in prep).
While the likelihood of these scenarios is unclear, and will need to be examined in more depth in future work, each would imply that the rest-optical continuum emission is not necessarily dominated by the accretion disk, resulting in lower bolometric luminosity estimates.
Lower AGN luminosities would be more consistent with the mid-IR limits from MIRI imaging surveys <cit.> and X-ray limits from Chandra <cit.>.
This may also bring the luminosity function into closer agreement with the quasar luminosity function, though, as we note in Section <ref>, it would require a factor of $\sim 100$ reduction in order to match the bolometric luminosity function from UV-selected quasars.
The bolometric LF could also come down if the bolometric correction were found to be lower than the nominal value, e.g. due to the lack of hot dust in this population.
We note, though, that the the current FIR data are not sufficient to constrain this; our stack gives an upper limit on the 3000 Å bolometric correction of $\lesssim 15$, entirely consistent with the nominal value of 5.15 <cit.>

While reducing the AGN contribution to the continuum may alleviate the tension with AGN models and the quasar LF, if the rest-optical continuum is instead dominated by stellar light, it implies a significant contribution of these objects to the stellar mass function.
As we show in Section <ref>, the LRDs could dominate the stellar mass function at the high-mass end, particularly at $z>7$.
This is true even with the MIRI data generally preferring solutions in which a recent, young starburst dominates the mass, which effectively represent a lower-limits on the mass.
Recent /NIRSpec results have identified LRDs with clear Balmer breaks, indicating the presence of old stellar populations alongside broad emission lines <cit.>; if ubiquitous in LRDs, this would raise their contribution to the stellar mass function.
However, even with high SNR spectroscopy from $1$–$5$ µm, the relative contribution from a Balmer break as opposed to a reddened power-law continuum is degenerate, with derived stellar masses ranging from $10^9$ to $10^{11}~M_\odot$ depending on the breakdown of the two components <cit.>.
MIRI imaging may help to break some of these degeneracies, as we find in our sample, and will be crucial to carefully decompose the continuum contribution from galaxies/AGN in the future.

§.§ What is the origin of the dust obscuration?

Whether the LRDs represent star-forming galaxies or AGN, their dust-reddening provides a view into the physics of AGN/ISM dust at high-$z$.
Despite SED-derived attenuation $A_V\gtrsim 2$, the LRDs entirely lack submillimeter or millimeter detections with SCUBA, ALMA, or otherwise, even in our deep stacks of hundreds of objects.
Their sub-mm faintness has been interpreted as evidence in favor of an AGN interpretation, as the ALMA limits are inconsistent with energy-balance models in which the attenuated starlight is reprocessed in the sub-mm <cit.>.
However, in a given bandpass, sub-mm observations are extremely sensitive to the temperature ($S_\nu \propto T^{-3.5}$ for fixed $L_{\rm IR}$), such that the difference between $T=25$ K and $T=75$ K represents $\sim 100$ times fainter emission <cit.>.
The typical assumption of a cold temperature ($T\sim 25$–$35$ K) may not apply in these systems; in fact, hotter dust is in fact expected at high redshift, as the increased CMB temperature, $T_{\rm CMB} = 2.73\,{\rm K}\times(1+z)$, sets an absolute floor on the dust temperature.
Moreover, a simple application of the Stefan-Boltzmann law ($L_{\rm IR} = 4\pi R_e^2\sigma T^4$) would predict higher dust temperatures for a more compact IR-emitting region <cit.>.
As we show in Figure <ref>, even a substantial luminosity of warm dust ($T\sim 150$ K) would be completely missed at the current limits of our samples, but would be expected if the obscuration originates on small scales.
This warm dust could be heated by the AGN, either directly or by successive absorption/re-emission <cit.>, or from intense star-formation, similar to what is seen in local ULIRGS <cit.>.

At the same time, the lack of mid-IR detections in LRDs appear at odds with the AGN interpretation.
The one object in our sample individually-detected in MIRI/F1800W (COS-105829) exhibits a flat SED from $\sim 7$–$18$ µm, as does the median stack.
Similar results have been reported for little red dots covered by the SMILES and JADES programs, with deep MIRI imaging in multiple bands <cit.>.
As discussed in Section <ref>, this may indicate minimal contribution from hot dust in the torus, which in turn may imply that the AGN contribution is overestimated.
However, an alternative explanation is that significant torus dust is not common in this sample, despite being ubiquitous in lower-redshift AGN.
Dust-deficient quasars are known to exist both in the local universe <cit.> and at $z\sim 6$ <cit.>, though these objects are typically blue in the rest-frame optical, lacking significant dust attenuation.
Given that the LRDs exhibit significant attenuation, the dust deficit is rather a hot dust deficit—a lack of dust at temperatures $\gtrsim 500$ K.

Perhaps a natural explanation is that the remarkable abundance and characteristic dust-reddening of the LRDs is directly related to their missing mid-IR emission.
In the canonical picture of obscured AGN, the hot dust emitting in the mid-IR is also the material obscuring the AGN, with $A_V \gtrsim 10$ mag across the lines of sight blocked by the torus <cit.>.
In the LRDs, this obscuring material may instead be located in clouds further from the central source, with colder dust.
These may originate from nuclear starburst disks <cit.>, or from “polar dust” in the NLR, ejected via winds <cit.>.
The hottest dust, in the clouds at the sublimation radius, may itself be optically thick, or be heated inefficiently due to obscuration by the (dust-free) BLR clouds <cit.>.
Regardless, if the colder obscuring clouds have a high covering factor, i.e. a roughly isotropic distribution, it could explain the ubiquity of reddened AGN at a given luminosity.

Such a geometry may be expected if these objects are in a nascent evolutionary stage.
At $z<2$, red quasars are thought to be in the blowout phase; as the “effective Eddington ratio” is lower for dusty gas, as opposed to pure ionized Hydrogen, moderate column densities of dusty gas around the accretion disk can be easily blown out via radiative feedback <cit.>.
In this case, the AGN effectively goes super-Eddington, blowing out its obscuring material, clearing lines of sight directly to the accretion disk.
In addition to radiatively-driven outflows, MHD winds driven off the accretion disk may present the right conditions for dust formation, and may even be linked to the formation of the dust torus <cit.>.
In fact, even the simplest models for torus formation begin with a larger, isotropic distribution of clouds, which form into a torus via anisotropic radiatively-driven winds <cit.>.
The mid-IR deficit in the LRDs may therefore indicate that dynamic feedback/winds are at play, regulating the dust content at the sublimation radius, or that the torus has not fully formed at this epoch.

§.§ Where there are AGN, there are galaxies

Finally, we want to emphasize that despite the focus of this paper on the two extreme physical interpretations of the LRDs, neither is likely completely true.
The physical mechanisms that trigger starbursts and AGN are similar, and we would expect that where there is early AGN growth, there is also significant star-formation <cit.>.
Already, we have clear evidence for AGN in the LRDs in the form of ubiquitous broad lines <cit.>, and also evidence for evolved stellar populations with clear Balmer breaks <cit.>.
Where these objects fit into the picture of black-hole galaxy co-evolution is yet to be determined.

It is widely believed that sub-millimeter galaxies at $z>2$ (i.e., 850 µm-selected dusty star-forming galaxies) represent the progenitors of massive $z\sim 1$–$2$ elliptical galaxies <cit.>.
Stellar age gradients in quiescent galaxies imply an inside-out formation scenario, in which the central cores formed in early, rapid starbursts, consistent with the compact sizes of SMGs <cit.>.

Extending this framework to higher redshift, the LRDs at $z\sim 5$–$9$ could represent the progenitors of at least some of the  first compact quiescent galaxies at $z\sim 3$–$4$.
Recent  observations find evidence for an over-abundance of massive quiescent galaxies at $z\sim 3$–$5$ <cit.>, including several which are quite compact <cit.> or maximally old, with formation redshifts $\gtrsim 10$ <cit.>.
The overall volume density of the LRDs in our sample is $n\sim 2\times 10^{-5}$ Mpc$^{-3}$, similar to the overall volume density of massive ($\log M_\star/M_\odot > 10.6$) quiescent galaxies presented in <cit.>.
Given also their compact sizes and similar formation redshifts, the LRDs may therefore be a natural progenitor population for these early quiescent galaxies.
Simulations predict that AGN feedback drives quenching <cit.>, and <cit.> find significant broad H$\alpha$ emission in the spectrum of an otherwise quiescent galaxy at $z=4.5$.
However, a number of complicating factors muddy any conclusions regarding the connections between LRDs and quiescent galaxies.
Samples of quiescent galaxies and LRDs are not mass-matched (nor are the stellar masses of LRDs well-constrained), and cover different redshift ranges.
Moreover, the duty cycle of the LRD phase would have to be of order unity if they are the direct progenitors of massive quiescent galaxies.
Understanding any potential evolutionary connection between LRDs and lower-redshift compact quiescent galaxies will require systematic follow-up and careful analysis to constrain the stellar contribution to the continuum.

§ CONCLUSIONS

We have presented a large sample of compact, extremely red objects/little red dots in COSMOS-Web.

* We present a sample of 434 LRDs selected based on their F444W compactness and F277W–F444W color $>1.5$, after removing likely brown dwarfs. By selection, we focus on the reddest subset of the LRD population, biasing our sample towards higher redshifts ($z\gtrsim 5$) but mitigating contamination from EELGs.

* We fit each candidate to both galaxy and quasar SED models to explore the two “edge cases,” or the two alternatively physical interpretations of the LRDs: either they are dominated by galaxy light from a compact, dusty starburst, or AGN light from a reddened quasar.

* We report the spectroscopic confirmation of one source, COS-756434, at $z=6.9993\pm 0.0001$, from public /NIRSpec PRISM data. The source exhibits clear broad lines in H$\alpha$ and H$\beta$, indicating strong AGN contribution.

* Under the interpretation of the LRDs as red quasars, we measure bolometric luminosities $\sim 10^{45-47}$ erg s$^{-1}$, spanning the gap between other -selected AGN and low-luminosity UV-selected quasars. We infer a bolometric luminosity function $\sim 100$ times higher than that for UV-selected quasars, consistent with lower limits from spectroscopic observations, implying a black hole accretion rate density at $z\sim 5$–$9$ rivaling that at cosmic noon.

* Under the alternative interpretation of the LRDs as massive/compact galaxies, perhaps hosting AGN but not dominated by their light, we measure stellar masses $\sim 10^{9-11}~M_\odot$. While many of the LRDs have masses too large to be found in this volume in a $\Lambda$CDM cosmology, when focusing on the subset of the sample with MIRI/F770W coverage, no source violates the limits imposed by $\Lambda$CDM. However, they dominate the high-mass end of the stellar mass function, and based on their compact sizes, must approach or exceed the maximal stellar mass surface densities observed in the local universe.

* We search all available X-ray, mid-IR, far-IR/sub-mm, and radio data in the COSMOS field to examine the multiwavelength properties of the LRDs. No source is individually detected, and stacks of the full sample yield non-detections in all bands except /. We provide the median-stacked limits and a model for the maximal SED consistent with the stacked constraints. This maximal SED model requires an anomalously low X-ray normalization ($\alpha_{\rm OX} \lesssim -2$) or Compton-thick absorption ($\log N_{\rm H} \gtrsim 24.5$) and is radio-quiet ($R^*_{1.4} \lesssim 10$). The model also allows for a substantial warm dust component with $T_{\rm dust} \sim 200\pm 75$, K despite non-detections in the mid-IR and millimeter wavelengths.

The interpretation of the “little red dots” as a population of reddened AGN has important implications for our understanding of black hole growth in the early universe.
Future large spectroscopic surveys will be needed to conduct a comprehensive census of LRDs, to elucidate their true nature at a population level, and deeper observations in the X-ray, mid/far-IR, and radio will be needed to determine their true properties.

HST (ACS), JWST (NIRCam, NIRSpec, and MIRI)

<cit.>,  <cit.>, SourceExtractor <cit.>,  <cit.>,  <cit.>,  <cit.>,  <cit.>,  <cit.>, STScI JWST Calibration Pipeline <cit.>.

§ ACKNOWLEDGEMENTS

The authors thank Gene Leung, Ryan Endsley, Jenny Greene, and Ignas Juodžbalis for useful discussions on the nature of the “little red dots.”
Support for this work was provided by NASA through grant JWST-GO-01727 and HST-AR-15802 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.
H.B.A. acknowledges the support of the UT Austin Astronomy Department and the UT Austin College of Natural Sciences through Harrington Graduate Fellowship, as well as the National Science Foundation for support through the NSF Graduate Research Fellowship Program.
C.M.C. thanks the National Science Foundation for support through grants AST-1814034 and AST-2009577 as well as the University of Texas at Austin College of Natural Sciences for support; C.M.C. also acknowledges support from the Research Corporation for Science Advancement from a 2019 Cottrell Scholar Award sponsored by IF/THEN, an initiative of Lyda Hill Philanthropies.
FG acknowledges the support from grant PRIN MIUR 2017-20173ML3WW_001. 'Opening the ALMA window on the cosmic evolution of gas, stars, and supermassive black holes.'
H.B.A., C.M.C., and others at UT-Austin acknowledge that they work at an institution that sits on indigenous land.
The Tonkawa lived in central Texas, and the Comanche and Apache moved through this area.
We pay our respects to all the American Indian and Indigenous Peoples and communities who have been or have become a part of these lands and territories in Texas.