A Census of Photometrically Selected Little Red Dots at 4 < z < 9 in JWST Blank Fields

Observations with the James Webb Space Telescope (JWST) have uncovered numerous faint active galactic nuclei (AGN) at z ∼ 5 and beyond. These objects are key to our understanding of the formation of supermassive black holes (SMBHs), their coevolution with host galaxies, as well as the role of AGN in cosmic reionization. Using photometric colors and size measurements, we perform a search for compact red objects in an array of blank deep JWST/NIRCam fields totaling ∼640 arcmin2. Our careful selection yields 260 reddened AGN candidates at 4 < z phot < 9, dominated by a point-source-like central component (〈r eff〉 < 130 pc) and displaying a dichotomy in their rest-frame colors (blue UV and red optical slopes). Quasar model fitting reveals our objects to be moderately dust-extincted (A V ∼ 1.6), which is reflected in their inferred bolometric luminosities of L bol = 1044–47 erg s−1 and fainter UV magnitudes M UV ≃ −17 to −22. Thanks to the large areas explored, we extend the existing dusty AGN luminosity functions to both fainter and brighter magnitudes, estimating their number densities to be ×100 higher than for UV-selected quasars of similar magnitudes. At the same time, they constitute only a small fraction of all UV-selected galaxies at similar redshifts, but this percentage rises to ∼10% for M UV ∼ − 22 at z ∼ 7. Finally, assuming a conservative case of accretion at the Eddington rate, we place a lower limit on the SMBH mass function at z ∼ 5, finding it to be consistent with both theory and previous JWST observations.


INTRODUCTION
The remarkable sensitivity and angular resolution of the James Webb Space Telescope (JWST ) at infrared wavelengths is enabling us to explore the distant Universe like never before.This allows for an exceptionally detailed examination of the characteristics of known high-z sources (e.g.Bunker et al. 2023;Maiolino et al. 2023a) and, at the same time, reveals the presence of more and farther galaxies (e.g.Adams et al. 2023;Atek et al. 2023;Austin et al. 2023;Bradley et al. 2023;Casey et al. 2023;Finkelstein et al. 2023;Naidu et al. 2022;Robertson et al. 2023), some of them spectroscopically confirmed beyond z > 13 (Curtis-Lake et al. 2023;Wang et al. 2023).
What truly tests our models and preconceived vision of galaxy evolution is not how early we can see these objects, but the questions they raise regarding the balance between their mass, UV luminosity and age.The excess of high-z galaxies at the bright end (M UV ≤ −20) of the UV luminosity function is in tension with current theoretical frameworks (Behroozi & Silk 2015;Dayal et al. 2017;Yung et al. 2019Yung et al. , 2020;;Behroozi et al. 2019Behroozi et al. , 2020;;Davé et al. 2019;Wilkins et al. 2022;Kannan et al. 2023;Mason et al. 2023;Mauerhofer & Dayal 2023), which suggests exotic initial mass functions, little to no dust attenuation, or a higher than anticipated density of galaxies undergoing active galactic nuclei (AGN) phenomena (e.g.Finkelstein & Bagley 2022;Pacucci et al. 2022;Boylan-Kolchin 2023;Ferrara et al. 2023;Fujimoto et al. 2023a;Lovell et al. 2023;Steinhardt et al. 2023;Sun et al. 2023).
Although early hints also existed in prior works (Morishita et al. 2020;Fujimoto et al. 2022;Endsley et al. 2023), one of the most intriguing discoveries from early JWST imaging is that of compact red sources with a "v-shaped" spectral energy distribution (SED), namely a blue UV continuum and a steep red slope in the restframe optical (Labbé et al. 2023a,b;Furtak et al. 2023a).While the first photometric selections of these objects included spatially resolved targets that could be early massive compact galaxies (Barro et al. 2023), spectra revealed clear evidence for broad Hα and/or Hβ emission indicative of actively accreting supermassive black holes (SMBH; Furtak et al. 2023b;Fujimoto et al. 2023a;Greene et al. 2023;Killi et al. 2023;Kocevski et al. 2023;Kokorev et al. 2023a;Matthee et al. 2023;Übler et al. 2023).
Dubbed "little red dots" (LRDs), these sources have SEDs characterized by a unique "v-shaped" continuum combined with their point source morphology (Labbé et al. 2023a,b;Furtak et al. 2023a).However, what truly makes the LRDs stand out is their high number densities.It appears that LRDs may account for a few percent of the galaxy population at z > 5, and are far more numerous than the lowest luminosity known UV-selected quasars.Likewise, they appear to account for ∼ 20% of broad-line selected active galactic nuclei (AGN) at z ∼ 5−6 (Greene et al. 2023;Harikane et al. 2023;Labbé et al. 2023b;Maiolino et al. 2023b), which is higher than the fraction of dusty red quasars at z < 2 (Banerji et al. 2015;Glikman et al. 2015).These red dots are generally observed at z ∼ 5 (Labbé et al. 2023b), but can potentially be found even at z > 9 (Leung et al. 2023).However, these initial LRD studies were performed with limited spectroscopic samples and/or small areas of the sky, covering only ∼ 20 -40 arcmin 2 .The numbers of compact red objects could therefore be further affected by cosmic variance, which makes it quite difficult to assess their real importance and diversity.
Extending the selection of this compact red population of low-luminosity broad-line AGN candidates to larger areas would thus be necessary to study their complete demographics, limiting the effects of cosmic variance.In addition, this would provide us with a sufficient level of detail toward a better understanding of the total number densities of obscured AGN at high-z as well as the potential role that these sources play in cosmic reionization (e.g.see Grazian et al. 2018;Mitra et al. 2018;Dayal et al. 2020;Trebitsch et al. 2023;Dayal et al. 2024).
In this work we present a carefully selected sample of 260 reddened AGN candidates in the ∼ 640 arcmin 2 area covering some of the deepest blank extragalactic JWST fields.Examining such a large area will ensure that we are reducing the effects of cosmic variance to a minimum, while our focus on blank fields lessens the selection biases and avoids volume uncertainties arising from lensing magnification.

OBSERVATIONS AND DATA
In this work we use JWST data from the following programs/fields -CEERS (# 1345, PI: S. Finkelstein, Bagley et al. 2022) in the EGS, PRIMER (# 1837, PI: J. Dunlop) in COSMOS and UDS.For the GOODS-S we combine the available data from multiple broad and medium band programs -FRESCO (# 1895, PI: P. Oesch, Oesch et al. 2023), JADES (# 1180, 1210, 1286, 1287PIs: D. Eisenstein, N. Luetzgendorf, Eisenstein et al. 2023a,b) and JEMS (# 1963, PI: C. Williams, Williams et al. 2023b).We provide a general overview of these 4 fields in Table 1.More detailed information, including specific filters, depths and survey designs can be found in overview papers for each data release.

JW ST Imaging Data Reduction
We homogeneously processed all the publicly available JWST imaging obtained with the NIRCam and MIRI in a variety of public JWST fields, presented in Table 1.The images have all been reduced with the grizli pipeline (Brammer 2023), using the jwst 1084.pmap, and follow the same methodology of (multiple) previous studies (e.g., Jin et al. 2023;Kokorev et al. 2023b;Valentino et al. 2023).Compared to the standard pipeline, we incorporate additional corrections to account for cosmic rays and stray light (see e.g., Bradley et al. 2023), 1/f noise, detector level artifacts ("wisps" and "snowballs") and bias in individual exposures (see e.g., Rigby et al. 2023).For the PRIMER data, we introduce an additional procedure that alleviates the detrimental effects of the diagonal striping seen in some exposures as was done in Valentino et al. (2023).Finally, our mosaics include the updated sky flats for all NIR-Cam filters.These reductions are publicly available as a part of the DAWN JWST Archive (DJA 1 ).
These data-sets are further complemented by including all available optical and near-infrared data from the Complete Hubble Archive for Galaxy Evolution (CHArGE, Kokorev et al. 2022).Individual JWST and HST exposures were aligned to the same astrometric 1 dawn-cph.github.io/dja/reference frame by using the Gaia DR3 (Gaia Collaboration et al. 2021), then co-added and drizzled (Fruchter & Hook 2002) to a 0. ′′ 04 pixel scale for all the JWST and HST filters.
Some of the fields we examine in this work have also been observed with MIRI, in one or more filters, sampling mostly the rest-frame near-infrared (NIR) at z ≳ 4.These data are, however, not uniform in the wavelength coverage, depth and area.In fact, only about a third of the objects in areas we examine have public MIRI data and even fewer are actually detected.While the inclusion of the MIRI photometry can assist in further identifying the presence (or absence) of dusty, power law-like AGN component in galaxies (e.g.see Yang et al. 2023;Williams et al. 2023a), doing so appro-priately within a context of a population study requires a degree of uniformity which the current MIRI data do not possess.Therefore, we have opted to exclude MIRI photometry from our current analysis to maintain consistency across various fields.

Source Extraction
The initial JWST catalog was constructed by utilizing a detection image combined from all noise weighted "wide" (W) NIRCam Long Wavelength (LW) filters available, which includes F277W, F356W and F444W.A similar detection method was already successfully employed in several works (see e.g.Jin et al. 2022;Kokorev et al. 2023b;Valentino et al. 2023;Weaver et al. 2023a).To extract the sources and produce a segmentation map, we used sep (Barbary 2016), a Python version of SExtractor (Bertin & Arnouts 1996).Photometry was extracted in circular apertures of increasing size.Correction from the aperture to the "total" values was performed by using the flux auto column output by sep, which is equivalent to the MAG AUTO from SExtractor, ensuring that for each source only flux belonging to its segment is taken into account.This method was shown to apply to both point-like and extended objects (Weaver et al. 2022(Weaver et al. , 2023b)), so we believe it to be adequate for our sources.
Additionally, we introduce a correction to account for the missing flux outside the Kron aperture (Kron 1980), by utilizing a method similar to the one used in Whitaker et al. (2011) and Weaver et al. (2023a).In short, this procedure involves computing the fraction of the missing light outside the circularized Kron radius by analyzing curves of growth of the point spread functions (PSF), which were obtained empirically, by stacking stars in these various fields.This correction is then applied to the flux auto values for each source.However, since our work focuses on compact (sub NIRCam PSF size) AGN candidates this additional correction does not strongly influence the derived flux densities.For the same reason, we use the total fluxes, computed from D=0. ′′ 36 apertures, unless specified otherwise.

IDENTIFYING COMPACT RED OBJECTS
The data from CEERS, PRIMER and various programs covering GOODS-S are well-suited for a photometric search for compact obscured AGN candidates.The available HST and JWST photometry covers a complete wavelength range from 0.4 -5 µm, in at least 7 broad and medium bands, reaching a median 5σ depth of 28.3 AB mag in F444W filter (Table 1).In our search, we explore blank fields covering a large area of ∼ 640 arcmin 2 , which are also completely independent.In return this will significantly limit the impact of cosmic variance and enable us to avoid dealing with the cosmic volume uncertainties introduced by the lensing magnification.
Colors alone would end up selecting both LRDs and extended red galaxies (see e.g.Labbé et al. 2023b;Williams et al. 2023a,b), so we introduce a further "compactness" cut to only select sources with high central flux concentration.To do that we use the ratio between the total flux in F444W between 0. ′′ 4 and 0. ′′ 2 apertures.Since roughly 17% of the LRD candidates followed up with NIRSpec turned out to be brown dwarfs (Burgasser et al. 2023), we would also like to minimize the incidence of these objects in our sample.To do that we adopt the brown dwarf removal criterion from Greene et al. (2023), based on the LRD spectra from NIRSPec/MSA.Finally, we also require our sources to be significantly (> 14σ) detected in F444W, and be brighter than 27.7 AB mags, to be consistent with the UNCOVER selection.The im-posed color cuts are then: which are effectively selecting our low (z < 6) and high (z > 6) redshift samples, respectively.The compactness is given by: compact = f f444w (0. ′′ 4)/f f444w (0. ′′ 2) < 1.7.
To limit the number of brown dwarfs in the sample we also adopt: bd removal = F115W − F200W > −0.5.
The final selection then becomes (red Applying the color criteria also means that every object has to be detected (> 3σ) in at least one band per color to make the selection meaningful.In case of a non-detection we use the 2σ upper limits, but only if the "brighter" band in the color is detected.Out of ∼ 408 000 objects covering 4 fields of interest, we end up selecting 334.Most importantly, we note that no information about photometric redshifts and underlying galaxy/AGN SEDs is used at this stage to avoid being biased by models.We discuss our photometric redshift estimate and its agreement with spec-z for sub-samples in the next sub section.

Size Measurements
While the compactness cut alone already successfully manages to select PSF-dominated point sources, we would like to provide a further fine-tuning to provide a fully quantitative rather than qualitative assessment.Figure 3. Observed MUV compared to the MUV expected from the dust-corrected L bol .We derive the expected MUV values by following the relation from Shen et al. (2020).We also show the data for red dots from Greene et al. (2023) (blue squares) and Matthee et al. (2023) (orange circles), as well as broad-line AGN at z > 4 from Maiolino et al. (2023b) (gray triangles).The gray dotted line shows the 1:1 relation.The observed MUV we derive is fainter than the expected value, with AUV varying from ∼ 0.6 − 4.2.While extreme, the UV attenuation is more than 5 magnitudes smaller when compared to values expected given our best-fit AV.
Miller 2023) in the F444W band.The primary goal of this is to ensure that the source is dominated by the PSF component in the reddest, least dust obscured, band as was done in Labbé et al. (2023b).We focus on the F444W band for this analysis, as the galactic origin of the rest-UV can not be ruled out with current photometric (or even spectroscopic) observations.Moreover, if an object is dominated by a single star-forming region, it could appear compact in rest-UV bands, but still be extended in the redder filters, making the F444W band the most physically constraining for our type of study.
Taking the PSF into account is imperative when measuring sizes of unresolved objects.We generate our F444W PSFs empirically for each field by following the methodology described in Skelton et al. (2014), Whitaker et al. (2019) and Weaver et al. (2023a).In brief, we identify non-saturated stars in every field by considering objects on the stellar locus, that are brighter than 24 AB mags and extract these candidates in stamps.These stamps are then centered and normalized to unity.The final PSFs are derived by averaging the weighted stamps, and are then normalized to match the enclosed energies of the expected JWST calibration levels within 4" diameter apertures2 .For more detail see the Appendix in Weaver et al. (2023a).
The light is modeled with a single Sérsic (Sérsic 1963) profile with the center, brightness, effective radius, Sérsic index, and axis ratio as free parameters.The prior for the index is uniform between 0.65 − 6 and the effective radius uniform between 0.25−5 pixels (0. ′′ 01 -0.′′ 2 ).For each source we create a 3" square cutout (75 pixels by 75 pixels) and mask any additional sources within the stamp.Parameter values and uncertainties are calculated using the Laplace approximation, assuming that the posterior is Gaussian.We exclude fits where the resulting χ 2 per pixel is greater than 2 or the best fit flux differs from the catalog value by more than 2 AB magnitudes.This excludes 15 sources from our sample which by visual inspection we find are untrustworthy due to contamination of bright nearby objects.
A source can be considered to be point-like if its effective radius in F444W band is lower than the empirical PSF FWHM (∼ 0. ′′ 15).It appears that our compact criterion is extremely effective at identifying PSF-like objects as none of the 319/334 sources with reliable fits exceed a diameter of 0. ′′ 08, when considering the 95 % size upper limits, corroborating the effectiveness of the compactness criterion described in Section 3.After carefully considering both colors and morphology when selecting our sample of AGN candidates, we are now able to proceed directly to the SED fitting.

Photometric Redshifts
To calculate photometric redshifts (z phot ) for our objects, we use the Python version of EAZY (Brammer et al. 2008).We choose the blue sfhz 13 model subset3 that contains redshift-dependent SFHs, and dust attenuation values.More specifically, the linear combinations of log-normal SFHs included in the template set are not allowed to exceed redshifts that start earlier than the age of the Universe (for more detail see Blanton & Roweis 2007).These models are further complemented by a blue galaxy template, derived from a JWST spectrum of a z = 8.50 galaxy with extreme line equivalent widths (ID4590; Carnall et al. 2022).
While it might seem counter-intuitive to use galaxy templates for what we believe to be AGN candidates, similar efforts presented in Labbé et al. (2023b) report a good agreement between deriving z phot with stellar templates alone, as opposed to stellar+AGN models, finding a very good agreement between the two.This is not surprising, as when it comes to photometric redshift fitting, the key deciding factors are the positions of the Lyman (∼ 912 Å) and Balmer (∼ 4000 Å) breaks.
For the LRDs, a general absence of significant stellar contribution in the rest-frame optical (e.g.Greene et al. 2023) would result in a lack of a noticeable Balmer break, however the trough of the "v-shape" in the restframe SEDs of observed LRDs is also located at roughly 4000 Å (e.g.Furtak et al. 2023b;Kokorev et al. 2023a).Indeed the existence of such a feature in LRDs has resulted in their misidentification as dusty star-forming galaxies, leading to stellar mass estimates which are in tension with ΛCDM, if all the light is attributed to starformation alone (e.g.see discussion in Boylan-Kolchin 2023; Kocevski et al. 2023;Labbé et al. 2023a;Steinhardt et al. 2023).
Spectroscopic follow-up of red compact objects hosting AGN BL emission has in fact shown a remarkable agreement between the z phot derived with EAZY (or similar routines) and z spec .For example, in GOODS-S, Matthee et al. (2023) report an average σ z = |∆z|/(1 + z spec ) = 0.01, and UNCOVER LRDs presented in Greene et al. (2023) have shown σ z ∼ 0.04.Similar consistency was also found between the initial photometric source selection and final spectra in JADES and CEERS fields (Maiolino et al. 2023b;Kocevski et al. 2023;Andika et al. 2024).As such we consider that utilizing EAZY to derive redshifts is adequate for our sample.
We fit all the available photometry, and upper limits from HST/F435W (λ obs ∼ 0.4 µm) to JWST/F444W (λ obs ∼ 4.4 µm) filters for our sample of 319 LRDs, limiting the redshift grid between 0.01 < z < 20.From the best fit EAZY SEDs we only derive photometric redshifts, delegating the estimation of physical parameters to a different template set discussed in the next section.The uncertainties on the photometric redshift are computed from the 16th and 84th percentiles of the redshift probability distributions -p(z).The availability of HST photometry allows us to securely constrain the presence of the Lyman break either through diminishing flux where a given filter overlaps with the break, or via upper limits for 40 % sources in our sample.Notably, however, the presence of the Lyman break is not required to securely constrain redshift for high-z LRDs, as the break in the optical part of the SED already places these objects in a unique color-color space, as initially shown in Labbé et al. (2023b) and then further confirmed in Greene et al. (2023).In addition, access to at least one NIRCam medium band further enhances redshift quality by allowing us to identify emission lines in broad band photometry.As a result none of our objects have double peaked redshift solutions.Despite that, since our sample is still only identified photometrically, appropriately taking into account z phot uncertainties is crucial when deriving the physical parameters and luminosity functions in the upcoming sections.

Quasar Template Fitting
While the origin of the rest-frame UV light in LRDs remains elusive, growing samples of JWST spectra consistently show either a complete absence or a lack of a significant contribution from the host galaxy to the total flux in the rest-frame optical (λ obs ≳ 2 µm) (Furtak et al. 2023b; Greene et al. 2023;Kokorev et al. 2023a).This is generally evidenced by comparing the expected L 5100 from broad Balmer series lines (generally Hβ and/or Hα) to the observed values.For example Greene et al. (2023) find that Hα-derived and observed L 5100 agree within a factor of two for the objects which have Hα PRISM coverage.Supporting this, Furtak et al. (2023b) and Kokorev et al. (2023a) also find that black holes masses (M BH ) derived via broad Hβ line luminosity and continuum are identical, given the scatter of the relations derived from AGN reverberation mapping (see e.g.Kaspi et al. 2000;Greene & Ho 2005), hinting at negligible stellar components.Furthermore none of the currently known spectroscopically confirmed LRDs in Abell 2744 are detected in ALMA at 1.2 mm down to < 70 µJy (2σ), which strongly limits the contribution of obscured star formation (see e.g.Labbé et al. 2013Labbé et al. , 2023b)), unless the dust is either very cold, very hot or diffuse.Indeed, when both JWST data and ALMA upper limits (Fujimoto et al. 2023b,c) are considered in a joint AGN+galaxy template fitting for the objects described in Furtak et al. and Kokorev et al. the contribution of the galaxy model to the total rest-frame optical light is negligible.Finally, robust measurements of effective radii for all UNCOVER LRDs, while also taking into account the empirically derived PSFs (see Weaver et al. 2023a), find no strong evidence for extended emission associated with the host galaxy in the F444W band.
Unfortunately, a lack of deep and uniform ALMA coverage for our objects prevents us from carrying out joint AGN+galaxy template fitting to ascertain the amount of AGN contribution to the optical SED.While it is possible to do it with only JWST photometry, such a fit would be too degenerate given the available number of bands and the number of models required.However, objects in our work were specifically selected with the color and compactness criteria largely mirroring those used to identify broad-line AGN in UNCOVER.It is reasonable therefore to assume that given similarly red (f λ ∝ λ 0−2 at 3100 -5200 Å rest) slopes, the rest optical continuum in our sources is also dominated by AGN light.
In terms of luminosity, the dust-obscured component is dominating the light from LRDs, and must be substantially attenuated (A V ∼ 1 − 2) in order to fit the observed red slope.Given that, the rest-UV light should not be visible at all (A UV > 10).From our photometry, however, we see that while the blue component is weak (only a few percent of red component), it is not reddened.This emission can be interpreted as either scattered light from the AGN itself, or the host galaxy (see discussion in Labbé et al. 2023b;Greene et al. 2023).However, even when spectra are available (Greene et al. 2023), given the similarities between the UV slopes of quasars and young star-forming galaxies, these two models are equally good representations of the observed light.Our available data also do not allow us to make a clear distinction between these two possibilities, therefore to avoid over-interpreting the origins of the rest-UV emission, we will assume the scattered light (unreddened) only template in our modeling.We caution the reader that as a result of the unknown origin of the blue light, the rest-UV properties derived in this paper do not necessarily represent physical conditions of the potential AGN our LRDs might host.Due to the aforementioned similarity between UV slopes in quasars and SFGs, the M UV values derived from both galaxy and quasar fits are thus nearly identical.
Following galaxy-only fits presented in Section 3.3 and keeping the above considerations in mind, we now would like to explore an AGN-only scenario where we model the observed light with a two component AGN model.The first one is the empirical model based on a composite of 2200 SDSS quasar spectra (Vanden Berk et al. 2001), and the second is derived from 27 near-infrared quasar spectra by Glikman et al. (2006).We then combine and renormalize both templates, allowing us to cover the full range from rest-UV to the near-infrared.
The same approach was already successfully employed in Labbé et al. (2023b) for a photometrically selected sample of red dots, and then later for PRISM spectra of 14 such objects in Greene et al. (2023) and Kokorev et al. (2023a).We fit the unreddened AGN component together with the Small Magellanic Cloud (SMC) law (Gordon et al. 2003) attenuated (A V = 0.1 -4) version of the same composite template.With the photometric redshift being fixed, we are fitting for a total of three free parameters.
We find the AGN-only fits to be a marginally better representation of the observed photometry, when compared to galaxy-only EAZY fits, with ⟨χ 2 ν ⟩ = 3.0 +3.7 −1.8 for the former and ⟨χ 2 ν ⟩ = 4.2 +6.6 −1.9 for the latter, with a difference of approximately ⟨∆χ 2 ν ⟩ ∼ 1.Similar findings were also presented in Labbé et al. (2023b), even without ALMA photometry, and Barro et al. (2023), where no significant χ 2 difference exists between dusty star-formation and reddened AGN models.

Extreme Equivalent Width of Emission Lines
Before focusing on the final sample of reddened AGN candidates we would like to conduct one final test which concerns the potential presence of strong emission lines, particularly Hα in the spectra of LRDs.Empirical quasar templates presented in Vanden Berk et al. ( 2001), which we used to fit our objects, generally contain bright AGN with a rest-frame Hα EW∼ 190 Å.Conversely, the recent literature results which analyze LRD spectra (Killi et al. 2023;Matthee et al. 2023) have found that the EW of Hα can reach and even exceed 500 Å.Such strong emission lines can contribute to the flux observed in the medium and even broad-band JWST filters in a non-negligible way, making the observed colors redder.In return, if such strong emission lines are not present in the templates, the value of the A V , and subsequently other physical properties dependent on it (e.g.L bol ) can be overestimated.
To test the significance of this effect we do the following.Starting with the original combined Vanden Berk et al. ( 2001) and Glikman et al. (2006) template set, we isolate the regions that cover the Hβ+[OIII] and Hα lines, and use a spline function to fit the continuum, while masking out the regions containing line complexes.While doing this we successfully verify that the measured rest-frame EW of these lines is exactly as the one reported in Vanden Berk et al. (2001).Finally, we uniformly boost the continuum subtracted spectrum to a point where the EW of the Hα line measures at ∼ 500 Å, and add back the continuum.We then re-fit all of our sources, following the same considerations as described in Section 3.4.
Using models with boosted emission line strengths, we indeed find the best-fit A V values to be systematically lower, albeit only by ∼ 0.1 mag, on average, compared to the original templates.This offset is well within our quoted uncertainty on the A V from the SED fitting.We thus conclude that even if some of our AGN candidates indeed contained very high EW Hα emission lines, the physical properties derived with the original Vanden Berk et al. (2001) template set should still remain valid.Despite being small, this offset is systematic, so we still incorporate it into out uncertainties when computing number densities in the subsequent sections.

Final Sample of Little Red Dots
Following the initial object selection and SED fitting we are now in a position to define our final sample of "little red dots".The primary goal of this work is to explore the photometrically selected dusty AGN candidates in the high-z Universe, compare these results to robust samples of spectroscopically identified BL AGN, and potentially extend these examinations to fainter UV magnitudes and bolometric luminosities.The accurate determination of these parameters is contingent upon good coverage of the spectral break between the blue and red components at ∼ 4000 Å.This is crucial to confirm that the selected objects indeed exhibit the characteristic features of LRDs.Furthermore, a thorough sampling of the rest-frame UV around ∼ 1450 Å is essential to accurately derive M UV , and the 5100 Å rest-frame optical continuum is needed for determining the bolometric luminosity -L bol .With the exception of CEERS, all of our fields benefit from full NIRCam filter coverage, spanning from F090W to F444W, which will cover the rest-frame UV at z ≳ 4. On the other hand CEERS has extremely deep (∼ 29.6 mag at 5σ) HST /ACS F814W coverage instead, which will also allow us to adequately compute M UV at 1450 Å in the same redshift range.We thus limit our exploration only to objects which have z > 4.
To do that we take into account the p(z) and ensure that the 16th percentile, rather than just the median of the p(z) lies above our redshift threshold (e.g.see Valentino et al. 2023).This final selection leaves us with a total of 260 red dots.

Physical Parameters
The physical sizes of objects in our final sample are extremely compact, with a median effective radius of r eff < 130 pc (95 % upper limit).This is much smaller when compared to the typical rest-optical sizes of starforming galaxies measured at z > 5 (e.g.see Kartaltepe et al. 2023;Ormerod et al. 2024), but is similar to the extremely compact red objects presented in Labbé et al. (2023a,b); Baggen et al. (2023) and LRDs spectroscopically confirmed as BL AGN (Furtak et al. 2023b;Kokorev et al. 2023a).Curiously, dusty galaxies at z > 7 explored in Akins et al. (2023) also show a lack of extended bright component (r eff < 200 pc), similar to LRDs.Although not as faint or centrally concentrated as our objects or other LRDs at these redshifts, some of these similarities might imply that these dusty objects can act as potential AGN hosts.
Using the standard relations, with the scatter, presented in Greene & Ho (2005) and taking into account our best-fit A V (∼ 0.6 -3.7 mags), we derive the L bol from the 5100 Å continuum, measured directly from best-fit SEDs.While this is not ideal, and assumes that the red continuum is AGN dominated, the SED model-dependent values represent our best guess for the intrinsic AGN luminosities.The inferred bolometric luminosities for the compact red objects from our sample thus range from L bol ≃ 10 43.5 − 10 46.5 erg/s.This range is slightly brighter than that derived in Labbé et al. (2023b) as we are not including any lensed fields, and thus likely fail to detect intrinsically fainter LRDs.We show the dust-corrected and observed L bol values in Figure 2.
In Figure 3 we explore how the observed M UV values of our LRDs compare to the expectations derived from the dust corrected bolometric luminosity (Shen et al. 2020).Given our median ⟨A V ⟩ ∼ 1.6, we expect the UV extinction to be large with A UV ∼ 9, however what we find is ⟨A UV ⟩ ∼ 2.5 (similar to e.g.Greene et al. 2023;Maiolino et al. 2023b;Matthee et al. 2023), more than six magnitudes difference.Adding to this, the shape of the rest-UV spectrum, while faint, does not hint at any dust extinction.This suggests that a second component, different from a reddened AGN spectrum is present in LRDs, however with our current data its origin can not be determined.
The final table which contains photometry, sizes and the physical parameters we derive for our sample is available in full online4 .We show an excerpt of the full table in the Appendix.

Estimating Effective Volumes
One of the key motivations for our work is to conduct an unbiased search for LRDs in some of the deepest blank fields observed with JWST.Our goal is to extend the existing luminosity functions that have been deduced from spectroscopic samples, for a larger sample covering a wider area.By applying the color and size criteria that have been adopted in recent LRD studies, such as those discussed in Greene et al. (2023) and Labbé et al. (2023b), we aim to exploit the large area in our blank fields to get better statistics, particularly at the bright and faint ends.This approach should allow us to examine how much these objects contribute to the observed L bol and M UV number densities.However, as we are working with a photometrically selected sample our analysis will be focused on the aggregate characteristics of the LRDs, rather than on detailed examinations of individual objects.2023) (open squares).Green stars show the UV number densities of the X-ray detected quasars at z ∼ 5 from Giallongo et al. (2019).Finally, light blue octagons represent the UVLF derived for galaxies hosting Type 2 AGN from Scholtz et al. (2023).We offset some of the literature points by ±0.05 dex horizontally for visualization purposes.Note that our measured UV luminosities do not decompose the AGN emission from the potential galaxy light.
Focusing only on the blank fields allows us to estimate the effective volumes for our objects in a rather simple way.In order to measure the observed number densities of our sample, we follow the standard V max method (Schmidt 1968).The 1/V max estimator has the advantage of simplicity and does not require prior assumptions on the functional form for the luminosity distribution, ideal for LRDs since their intrinsic luminosity/mass distributions are unknown.To compute the number density for some property -x, we can then say: where ∆x is the width of the bin and V max,i is the maximum volume over which a source can be detected.In return, V max,i depends on the effective survey area -A, lower redshift bin boundary -z min and maximum observable redshift -z max .The latter is computed empirically from the detection limits of the survey, given the selection criteria, and cannot exceed the maximum redshift of the bin.We obtain the total survey areas by adding up all the non-masked pixels in our detection images, as presented in Table 1.Given how bright we require our objects to be (F444W< 27.7 mags at SN> 14) it might seem that z max would always exceed the maximum redshift of the bin, however this does not take into account the fact that our objects have to be detected in at least four bands (at > 3σ) to make color selection robust.We choose to remain conservative with our volume corrections, by only requiring one band per color combination to be detected.The z max values for each object are then estimated by considering our color selection laid out in Section 3. Uncertainties on our number densities are then derived in the following way.We consider the standard errors arising from Poisson statistics and compute them as prescribed in Gehrels (1986).Given that we only consider a photometric sample in our work, the uncertainty on the photometric redshift has to be taken into account appropriately in order to derive realistic errors on the physical parameters and number densities.
To do that we follow the approach described in Marchesini et al. (2009).Briefly, for each object we use Monte Carlo simulations to determine whether the objects fall into the redshift bin by considering their p(z).The final uncertainties are then a quadrature sum of the Poisson and p(z) errors.
Accounting for magnitude incompleteness effects as it is normally done for galaxy luminosity functions is not possible in our case, since it relies on making assumptions regarding the intrinsic source distributions.However, as Labbé et al. (2023b) already note, the requirement for objects to be bright in the detection band should lessen, but not eliminate altogether, the effect of magnitude incompleteness.Given that our sources are compact, we also expect that all of them will be detected above the brightness limit, diminishing the need to consider the incompleteness as a function of surface brightness.Despite the complex selection function, it is still possible to define a limit beyond which the derived number densities, for observed quantities, are expected to become incomplete.We will discuss this in the next sections.

UV Luminosity Function
In Figure 4 we present the UV luminosity functions in two redshift bins, at z ∼ 5 and z ∼ 7, derived from the continuum luminosity at rest frame 1450 Å as normally done for blue quasars.We list the number counts alongside the uncertainties in Table 2.The widths of our redshift bins were chosen to best align with the current literature results, for ease of comparison, as well as to ensure that photometric redshift uncertainties have a minimal impact on the luminosity functions.
At z ∼ 5, we find that the number densities of our red color selected AGN are ∼ 2 dex higher compared to the UV-selected quasars at similar magnitudes, depending on the extrapolation (Niida et al. 2020).As an upper limit on number density of quasars at z ∼ 5, we also compare to the results presented in Kulkarni et al. (2019), which combine both UV-bright quasars (M UV < −24) and UV-faint X-ray detected AGN (Giallongo et al. 2019) in their UVLF.
Before comparing to current observational results (Kocevski et al. 2023;Greene et al. 2023;Matthee et al. 2023), we note that it is difficult to accurately define the selection function for spectroscopically observed samples, and therefore derive the V max corrections.As such the number densities computed in these works should be treated as lower-limits.In our case the sample is selected via photometry and we derive our V max correction based on the selection criteria alone.This is done to avoid significantly over-estimating the number counts, thus misrepresenting the true abundance of red dots.It is also unlikely that this difference is a result of brown dwarfs contaminating our sample here, since we introduce an additional color cut from Greene et al. (2023) based on the spectra from Burgasser et al. (2023).
Taking the uncertainties into account, we find that our UV number counts are consistent with JWST -selected red BL AGN samples (Greene et al. 2023;Labbé et al. 2023b;Matthee et al. 2023), at least at M UV ∼ −19 and brighter.Confirming the initial findings for the UN-COVER red-dots presented in Greene et al. (2023) and Labbé et al. (2023b), we also find that our sample accounts for ∼ 10 -30 % of total BL AGN populations at high-z (Harikane et al. 2023;Maiolino et al. 2023b) and is largely comparable to the X-ray selected quasars from Giallongo et al. (2019), although in the case of the latter we infer higher number densities at fainter UV magnitudes.However it is worth noting that differences between the resolution of Chandra X-ray data and optical light from HST can lead to uncertainties when associating X-ray emission to the galaxies being present in the same patch of the sky.Curiously enough, the recovered scarcity of compact red sources compared to galaxies is in stark contrast to the density of Type 2 AGN hosts inferred from the recent JADES spectra (Scholtz et al. 2023) which report as much as a 20% contribution to the galaxy luminosity functions at z ∼ 5.
When moving to the z ∼ 7 bin, the results for the UVLF at both bright and faint luminosities are inconclusive, due to the limited number of objects and the uncertainty on the photometric redshifts.However, we are again consistent with the number densities of UN-COVER BL AGN from Greene et al.. Comparing to the luminosity functions of UV selected quasars from Matsuoka et al. (2023) at z ∼ 7, and extrapolating to fainter magnitudes, we find a 2 − 3 dex offset between the number densities at M UV > −22, roughly a factor of ten larger than in the lower redshift bin.Alongside our UVLF we also highlight the median M UV 5σ completeness limits.This is derived by considering the depths of filters covering the rest frame ∼ 1450 Å at a given redshift, and whether a source of a given M UV would be detected at a S/N> 5.As such we should be complete down to M UV ∼ −18.5 at z = 5, and M UV ∼ −19.0 at z = 7.
Following Bouwens et al. (2015), we fit our observed UV number densities with a Schechter (Schechter 1976) function, allowing all parameters to be free.We only fit data brighter than M UV = −18.5 at z ≃ 5 and M UV = −19.0 at z ≃ 7 as our number densities indicate that we are likely becoming incomplete at such faint magnitudes.The best-fit is shown in Figure 4 and the parameters are listed in Table 3.In both redshift bins we find that our red compact objects constitute roughly 3 − 5% of the total star-forming galaxy populations (Bouwens et al. 2021), consistent with the spectroscopic samples of red-dots (Greene et al. 2023).We also report shallower faint-end slopes compared to SF galaxies, however it is likely that the observed flattening of the UVLF for LRDs is induced by the incompleteness of our sample at fainter UV magnitudes, rather than any lack of compact red sources at fainter UV magnitudes.Deeper surveys would be required to robustly constrain the faint-end slope of LRDs.It appears that the LRD luminosities start to become comparable or even outshine galaxies at brighter (M UV ∼ −23) magnitudes, which is particularly prominent in the z ∼ 7 bin.This might be an expected consequence of the assembly of increasingly massive black holes with cosmic time (e.g.see Piana et al. 2022), or selection effects (see Volonteri et al. 2017), however we note that our number counts for the brightest objects are uncertain due to a limited amount of detections available.
Provided that our color and morphology selection is comparably successful at identifying reddened AGN as was previously shown (Labbé et al. 2023b;Greene et al. 2023), it appears that the compact red sources identified in blank JWST fields are ∼ 1 − 2 dex more numerous compared to the pre-JWST studies of known UVselected faint quasars (M UV > −21).While this trend has been consistently re-emerging in the new JWST results (e.g.Furtak et al. 2023b;Kokorev et al. 2023a;Maiolino et al. 2023b;Pacucci et al. 2023), it is worth noting that earlier works have already hinted that the number density of UV faint, dusty active black holes could have been much higher than previously thought (Laporte et al. 2017;Morishita et al. 2020;Fujimoto et al. 2022).For example, both Fujimoto et al. (2022) and Morishita et al. (2020) find that the less-luminous red quasar population could be anywhere from 10 to 100 times more common at z ∼ 7 − 8, compared to quasars luminosity functions at z ∼ 6, constructed from ground-based datasets (e.g.Matsuoka et al. 2018;Kato et al. 2020;Niida et al. 2020).The results of this work, together with the recent efforts to study compact red sources, therefore imply that these faint quasar populations, missed by previous surveys, are now being uncovered by the deep and rich multi-wavelength photometry and spectra from JWST.It is also important to highlight that if we extrapolate our UVLF to brighter magnitudes, the number density of LRDs becomes comparable to and then drops below the density of UV-selected quasars.Currently, however, it is not possible to speculate whether this is a real physical effect, or simply a consequence of insufficient volumes sampled.

Bolometric Luminosity Function
Our SED fitting results show that the fraction of the UV light contributing to the total luminosity is small as a result of significant dust reddening (A V = 0.6-2.7) in these objects.Even with spectra in hand (e.g.Greene et al. 2023), it is not easy to establish the origins of the rest UV light, which could be AGN light, either scattered or transmitted through patchy dust clouds, or unobscured light from star-formation in the host galaxy.As such, while we put our LRDs in the context of their observed UV luminosities, this does not explicitly describe the physics of potential AGN these compact objects host.Due to that, and also to carry out a comparison with existing spectroscopic bolometric luminosity functions of dusty BL AGN, we also present bolometric luminosity functions in Figure 5 and Table 2.While dust attenuation, estimated from SED fitting, can be an important source of uncertainty, we note that even if all our A V values were grossly overestimated, this would only change the dust-corrected L bol by a factor of ×5, given the average A V ∼ 1.6.This would thus only impact the number densities by ∼ √ 5 on average, which is insignificant when compared to the Poisson and z phot errors.Finally, to account for the potential presence of emission lines with high EW, we incorporate the additional ∼ 0.1 systematic shift in A V and apply it to the uncertainties on the L bol .
Understanding where the bolometric luminosity functions start to become incomplete is less straightforward compared to the observed quantities like M UV , as the former also relies on dust correction derived via SED modeling, and the assumptions made regarding the AGN contribution to the rest-frame optical emission.For this reason we do not define a completeness cut, like we do for M UV , however for each bin of bolometric luminosity with a width of 1 dex, we also compute the V max correction as was described in Section 4.1.
Our number densities again confirm that the redcompact AGN candidates are roughly 100 times more abundant compared to the UV-selected AGN at similar intrinsic luminosities (Shen et al. 2020) at z ∼ 5.The L bol number densities which we recover are comparable to the previous results for these objects derived in Greene et al. (2023), Labbé et al. (2023b) and Matthee et al. (2023) for L bol −10 45−46 erg/s.Curiously however, we find a factor of ten more LRDs compared to Greene et al. (2023) at L bol ∼ 10 44 erg/s.Nominally, the median NIRSpec depth at 4 µm of the UNCOVER followup of Abell 2744 is shallower compared to the fields we examine, as such it is perhaps unsurprising that we can recover a large fraction of intrinsically faint objects.However since the L bol is not an observed quantity and depends on SED modeling to calculate the dust correction it is difficult to ascertain whether the higher number densities we recover are indeed caused by the depth dif- ference, or simply the bias caused by the spectroscopiconly sample selection and lensed volumes in UNCOVER.Moreover, the mask design of the UNCOVER NIRSpec observations in the Abell 2744 field was also driven by optimizing the MSA coverage to include other targets of interest and was not just limited to LRDs.This, in return, induces selection effects which would not be possible to trace back and correct for.We additionally compare our bolometric luminosity function to the latest version of the semi-analytic Delphi models (Dayal et al. 2019(Dayal et al. , 2024)).In brief, these models follow the seeding and growth of BHs from z ∼ 40 down to z ∼ 4.5.Included are also all the key processes of merger-and accretion driven assembly of dark matter halos and their baryonic component (including black holes).The model also follows star formation and black hole growth and their respective feedbacks in determining the assembly of these early systems.Finally, Delphi models also include key dust processes to yield dust-to-stellar mass ratios, which with a baseline constructed against the latest ALMA observations (Dayal et al. 2022;Mauerhofer & Dayal 2023).All of this was specifically done to ensure Delphi can reproduce both the intrinsically faint and reddened sources in the recent literature, i.e. the LRDs.
We find that while our observations are comparable to Delphi results at L bol < 10 47 erg/s at z ∼ 7, these models fail to reproduce the high number density of bright objects we report.At z ∼ 5 on the other hand, our densities consistently fall 1 dex below Delphi predictions.This in return could suggest that the fraction of dusty AGN is diminishing toward later times, as they potentially transition to unobscured quasars (Fu et al. 2017;Fujimoto et al. 2022).
Finally, we also see a higher prevalence of intrinsically brighter objects at ∼ L bol − 10 47 , which is likely consequence of larger volumes sampled in our analysis.As already mentioned in Greene et al. (2023) however, it is worth noting that the uncertainties on the L bol − L 5100 relation, dust correction and assuming that these objects are dominated by AGN light at rest-frame optical could cause objects to scatter upwards into the high luminosity bins.We only recover a single object above L bol = 10 47 erg/s at z ∼ 5 and below L bol = 10 44 erg/s at z ∼ 7, respectively.As this is insufficient to properly compute luminosity functions, these are shown as upper (lower) limits in Figure 5 derived by combining the Poisson (Gehrels 1986) and photo-z (Marchesini et al. 2009) uncertainties.The solid and dashed blue lines show the result from Delphi (Dayal et al. 2014(Dayal et al. , 2019(Dayal et al. , 2020) ) simulations for all and bright (L bol ≳ 10 44 erg/s) black holes, respectively.Measured number densities of our LRDs agree well with the spectroscopic sample from (Matthee et al. 2023) and simulations given a λ edd ∼ 1.
lines (e.g.Hα, Hβ in rest-optical or Mg II in rest-UV) coupled with the luminosity of their broad components or the luminosity derived from the AGN continuum at λ rest = 5100 Å (see e.g.Kaspi et al. 2000;Greene & Ho 2005).
While secure determination of the black hole mass in our compact objects is not possible due the photometric nature of the sample, we can still place a lower limit on black hole masses (M BH ), by making a set of conservative assumptions.To do that, we adopt a scenario where all our AGN candidates accrete at Eddington rate (the physical limit at which outward radiation pressure balances inward gravitational force), such that L bol ∼ L edd , where L edd is directly proportional to M BH .While in the literature the Eddington rate (λ edd ) for confirmed AGN in LRDs was found to vary between 10 − 40 % (Furtak et al. 2023b;Greene et al. 2023;Kokorev et al. 2023a), we would like to remain conservative and compute a lower limit on the M BH .It is also worth noting that, given this range of λ edd in the literature, this is still small compared to other sources of uncertainty in our work.We calculate the M BH directly from the dust-corrected L bol and compute the SMBH mass function as described in the previous sections.We present our SMBH mass function in Figure 6 and Table 2, binned into 0.5 dex intervals to allow a direct comparison with the existing observational and theoretical results in this redshift range.We limit this investigation to the z ∼ 5 range only.As before, we note that the effect of the A V uncertainty on our number densities is expected to be at most ∼ √ 5, largely overshadowed by the Poisson and redshift errors.
We are now in a position to compare our mass function to the existing samples of both bright and faint quasars at z ∼ 5. We start with the latest ground based examination of the quasar mass function at z ∼ 4 from He et al. (2023).The authors focus on a sample of ∼ 1500 faint broad-line AGN, from a combined Hyper Suprime Cam (HSC) and SDSS dataset, allowing them to extend their examination to a low mass range we are most interested in (M BH ≃ 10 7−8 M ⊙ ).We find that, while our result is consistent with the ground based mass function in the high mass regime M BH > 10 8 M ⊙ , our number densities diverge below that mass and continue to rise up to ∼ 10 −4 cMpc −3 at M BH ≃ 10 6 M ⊙ .Barring the color selection, it is possible that this effect is purely observational, as the SDSS/HSC detection limits in the rest-UV are much shallower compared to the JWST fields we explore.
Furthermore, we contrast our result to the BH mass function at z ∼ 5 based on a sample of LRDs from the slitless JWST derived in Matthee et al. (2023).Given the uncertainties, we find our results to be consistent within 1σ, although we do not find a sharp drop-off in number densities at M BH < 10 7 M ⊙ , likely driven by low mass incompleteness of grism data as mentioned in Matthee et al..
Naturally, the fact that both our result and Matthee et al. (2023) find more low mass black holes compared to He et al. (2023) is unsurprising, given the depth and wavelength coverage of JWST data.Quite curiously however, it appears that the mass function derived from dusty compact LRDs seems to nicely continue the rising trend of ground-based data, and extend the SMBH mass functions toward M BH ∼ 10 6 M ⊙ .In this redshift range, the maximum volume sampled by our multi-field investigation is roughly equal to ∼ 3.0×10 6 cMpc 3 .Therefore, taking into account the results of He et al. (2023), we should expect only one object with M BH ∼ 10 8.5 M ⊙ in our images, which indeed is the case.Detection of AGN hosting black hole with masses larger than that, would, however, require survey sizes ten to twenty times larger.
Before drawing conclusions, we would like to conduct a final check, and compare our result to the hydrodynamical simulation EAGLE (Rosas-Guevara et al. 2016) and semi-analytic Delphi (Dayal et al. 2014(Dayal et al. , 2019(Dayal et al. , 2020(Dayal et al. , 2024) ) simulations describing masses of SMBHs in the same redshift range.We limit our examination of the Delphi models to the bright (L bol ≳ 10 44 erg/s) regime to match the same luminosity range covered by our objects.In the intermediate, to low mass end (M BH < 10 7.5 M ⊙ ) our results are in broad (∼ 2σ) agreement with both EAGLE and Delphi, however in both cases we start to see a significant difference in number densities as we move to higher masses.Perhaps a worthwhile question to ask in this case, is whether more of these high-mass BHs would be found in larger areas, we will discuss this in the later section.
Examining both the UV and bolometric luminosity functions we note that LRDs only represent ∼ 25 % of the total type I (broad-line) AGN population as inferred by Harikane et al. (2023); Maiolino et al. (2023b), even less so compared to the most recent examination of type II AGN hosts from JADES (Scholtz et al. 2023), where LRDs are 30-40 times less numerous.Taking this into account, we can conclude that, at least at z ∼ 5, LRDs appear to represent at most 1 % of the total accreting BH population over the L bol ∼ 10 44−47 erg/s range.The fact that LRDs are truly a distinct population of dusty broad-line AGN can therefore explain the observed ∼ 2 dex disparity between our results and simulations.Finally we would like to reiterate that our investigation of the BH mass function relies on assuming the most conservative case of accretion at exactly the Eddington rate, as we do not want to erroneously overestimate the number of high-mass black holes.Keeping that in mind, in Figure 6 we also show how our mass functions would change if we were to assume an Eddington ratio of 10 % instead.In this case we find that while our number densities compared to UV samples are still high, we now more closely match the abundance of high mass SMBHs predicted by Delphi.However, until broad emission line observations for all our sources are available, the value of λ edd will remain uncertain.

Abundance of bright compact sources
Previously limited to UV−selected samples at z ≲ 6 (Kashikawa et al. 2015;Bañados et al. 2018;Matsuoka et al. 2018;Inayoshi et al. 2020;Wang et al. 2021;Fan et al. 2022) we are now able to use JWST to reveal the presence of AGN during (e.g.Kocevski et al. 2023;Matthee et al. 2023;Übler et al. 2023) and even beyond the epoch of reionization, (e.g.Bogdan et al. 2023;Furtak et al. 2023c;Goulding et al. 2023;Kokorev et al. 2023a;Larson et al. 2023;Lambrides et al. 2023;Maiolino et al. 2023a) only hundreds of millions of years after the Big Bang.Standing out among these early studies of active black holes, is the population of reddened type I AGN, the so called "little red dots" (Greene et al. 2023;Labbé et al. 2023b;Matthee et al. 2023).
While the study of this unique population has been mostly limited to small spectroscopic samples, most recent efforts focused on the expansive Abell 2744 JWST data-set (Labbé et al. 2023b) have shown great promise at using a combination of NIRCam colors and morphology to identify reddened AGN.This initial photometric selection was shown to be remarkably successful with ∼ 80 % of targets indeed confirmed as z > 5 dusty broad-line AGN (Fujimoto et al. 2023a;Furtak et al. 2023c;Greene et al. 2023;Kokorev et al. 2023a).It is clear that these objects play an important role in the story of black hole growth at early times, however so far a systematic review of these enigmatic AGN across multiple fields has not been undertaken.
Motivated by the success of this photometric selection, we present a sample of 260 reddened BL AGN candidates in the 4 < z < 9 redshift range, covering 4 separate blank JWST fields with a total area of ∼ 640 arcmin 2 .We uniformly reduce the NIRCam JWST data from a variety of public programs, complementing our photometric coverage with archival HST observations.We perform a color and morphology selection to identify the most promising compact objects which display a dichotomy in their observed SED shapes, namely a blue rest-UV continuum, and a red power lawlike rest-optical.
Using model fitting, we derive photometric redshifts as well as a range of physical parameters including A V , M UV and dust corrected L bol .We split our objects into two redshifts bins at z ∼ 5 and z ∼ 7 and explore their contribution to the UV and bolometric luminosity functions of star-forming galaxies, as well as UV-selected quasars.Consistent with the previous works (Greene et al. 2023;Matthee et al. 2023;Maiolino et al. 2023b) exploring high-z BL AGN, we find that number densities of these objects at z > 5 are surprisingly high, in excess of ×100 compared to faint UV selected quasars (e.g.Niida et al. 2020;He et al. 2023), while also accounting for ∼ 20% of the total BL AGN population (Harikane et al. 2023;Maiolino et al. 2023b), and ∼ 1 − 2 % of UVselected star forming galaxies (e.g.Bouwens et al. 2021).Moreover, while some of these objects were potentially pinned down as potential sources of reionization in their local environment (Fujimoto et al. 2023a), it appears that their UV luminosities are still insufficient to contribute to reionization in a significant way (Dayal et al. 2024).
Assuming accretion at the Eddington rate, we also place a lower limit on the M BH of our objects, finding that some of these can already be very massive (M BH > 10 7 M ⊙ ) only a few hundred of years after the Big Bang.Using these masses we were also able to construct our prediction for the SMBH mass function, and for the first time, extend it to the low-mass (< 10 7 M ⊙ ) regime.We find that our mass function results are completely consistent with the number densities derived for faint dusty AGN from Matthee et al. (2023) at intermediate masses, and are comparable to those from UV-selected samples at high mass (He et al. 2023).We note however that while their number densities are similar, the sample presented in He et al. consists of unobscured quasars, and not LRDs, which are thought to be dust obscured AGN.We find that both hydrodynamical and semi-analytic predictions for the number of black holes at this redshift match our observations below ∼ 10 7.5 M ⊙ , however start to disagree by almost 2 dex at higher masses.These massive and bright black holes are likely to be heavily obscured in the rest-frame UV, and thus are not selected as LRDs due to a lack of a clear blue component.

Final Remarks
Using observed NIR colors to pick out active black holes in extragalactic fields is by no means a novel endeavor and has already been successfully done with the IRAC instrument, onboard the Spitzer Space Telescope (Lacy et al. 2004;Stern et al. 2005;Donley et al. 2012).This method, however, is still in its infancy when it comes to JWST (Labbé et al. 2023b;Andika et al. 2024).As already pointed out by Matthee et al. (2023), very few JWST programs that detect dusty AGN, were actually designed with AGN in mind, implying that there is still more that we can to address a growing number of questions about this population.
Firstly, the physical mechanisms that govern BH formation and growth in these systems are still poorly understood, however there exists already a growing body of works which try to decipher this enigmatic population (Greene et al. 2023;Silk et al. 2024).One such puzzle, is the origin of the blue light found in LRDs.The similarity between blue slopes of low-metallicity star-forming galaxies and quasars does not allow us to make a clear assessment of whether the rest-UV light originates from the AGN itself or the compact host galaxy surrounding it from the continuum alone.One way to solve this is to target the Mg II doublet (λ rest ∼ 2800 Å), the CIV, SiIV or HeII (λ rest ∼ 1840 Å) lines (e.g.see Maiolino et al. 2023a) and confirm whether these are broadened or not.This however would require longer integration times with NIRSpec as medium or even high resolution gratings would be required to achieve the necessary spectral fidelity.Moreover, while in principle detection of broad UV emission lines could hint at the AGN being responsible for some UV light, this would not necessarily mean that UV continuum also originates from the same source.
Secondly, there are no models as of yet, which can adequately describe the light we see emerging from these objects.So far, we mainly have had to rely on combinations of dust-free and dust attenuated empirical models of local quasars, which might be adequately describing AGN at high-z.Moreover, the lingering uncertainty on the A V correction can introduce some biases in our estimates of L bol and the M BH .One solution to alleviate this is to stack spectra of known LRDs, to define sets of reliable models describing these populations and aiding with further photometric selection.
Thirdly, it is crucial to note that a substantial proportion of massive SMBHs (with M BH > 10 8 M ⊙ ) at high redshifts (high-z) can be heavily obscured, as implied by Delphi.Similar conclusions were already reached from the X-ray luminosity functions at z > 5, both from simulations (Ni et al. 2020), and observations (Aird et al. 2015;Vito et al. 2018).In addition Trebitsch et al. (2019) have shown that accreting SMBHs in Lyman break galaxies are rarely UV-bright.With this in mind, selecting these massive AGN as LRDs would thus not be possible, as a combination of very deep rest-UV imaging and large areas are required.Despite that, these objects should still appear bright in near-infrared, which opens up the possibility of effectively identifying them through large or parallel surveys using MIRI.
Finally, as was already shown recently by Williams et al. (2023a) and Pérez-González et al. (2024), MIRI can also assist with clarifying true numbers of AGN among LRDs as some of these could be dusty progenitors of compact ellipticals.
Early results from JWST have already provided us with quite unexpected and remarkable results regarding number densities of early AGN, leading to a shift in our understanding of their formation and growth in the early Universe.Our results highlight the potential of using NIRCam alone to select reddened AGN at high-z in an effort to better understand their properties and abundance.While some limitations to this technique exist, as we already discuss in our work, this provides a crucial set of next steps in order to bridge the gap between UV bright quasars and faint SMBHs.However, it is already evident that the importance of faint, reddened AGN at early times can not be overlooked.

Figure 1 .
Figure1.Selection and analysis of LRD candidates.Top: Sample selection criteria.The left and central panels show modified "red 1" (z ≲ 6) and "red 2" (z ≳ 6) color-color cuts fromLabbé et al. (2023b).The right panel shows the compactness cut of our sample.Selected objects are highlighted as maroon circles, while grayscale hexbins show the full catalog.The compact red sources are clear outliers in color-color-compactness space.Colorbar is shared between all plots.Bottom: An example of best-fit SEDs to the photometry of LRD candidates with the dust-free (blue) and dusty (red) AGN templates(Vanden Berk et al. 2001;Glikman et al. 2006) at representative redshifts of z ∼ 6 and z ∼ 8.The combined model is shown in black.Detections (> 3σ) are shown as red circles, while upper limits (primarily from HST ) are shown as downward arrows.On the right of each SED we show 1. ′′ 5 color composite cutouts in the short (F115W/F150W/F200W) and long (F277W/F356W/F444W) NIRCam filters.

Figure 2 .
Figure 2. Top: Distribution of the L bol assuming an AGN dominated rest-frame optical continuum, with best-fit dust correction applied.Our compact sources span a wide range of luminosities across the redshift range of interest.We contrast our results to the BL AGN from Greene et al. (2023, blue squares), Matthee et al. (2023, orange circles), z > 6.5 quasars from Yang et al. (2021, gray crosses), high-z AGN from Maiolino et al. (2023b, magenta circles) and finally objects hosting Type 2 AGN from Scholtz et al. (2023, light blue circles).Assuming λ edd = 1, we show what L bol would correspond to MBH = 10 6−8 M⊙ (dashed lines).Bottom: Same as before, but we show the distribution of L bol without correcting for the best-fit dust-attenuation.Our final sample has a median AV ∼ 1.6 +1.1 −1.0 .On the side of each panel we also show histograms highlighting the redshift and L bol distributions.

Figure 5 .
Figure5.Bolometric luminosity functions in the z∈ [4.5, 6.5] (left) and z∈ [6.5, 8.5] (right) bins, derived from L5100, assuming rest-frame optical continuum is AGN dominated.The number densities have been Vmax and completeness corrected.Uncertainties are derived from Poisson noise(Gehrels 1986).Arrows show upper limits on the derived number densities.Blue squares show upper limits derived for spectroscopically confirmed "little red dots" in the UNCOVER data fromGreene et al. (2023).In addition in the lowest redshift bin we show the NIRCam grism result ofMatthee et al. (2023) (open circle).Dashed lines show the pre-JWST L bol relation derived inShen et al. (2020).Finally, the blue lines show the luminosity function from the Delphi semi-analytic models(Dayal et al. 2019) that grow SMBHs from seeds.

4. 4 .Figure 6 .
Figure 6.The SMBH mass function, assuming λ edd = 1, of our sample in the 4.5 < z < 6.5 range.Red arrows show how our mass function would change, if we assumed a lower Eddington ratio of 10 % .We overlay the SMBH mass function from Matthee et al. (2023) at z ∼ 5 in orange and HSC+SDSS derived BH mass function from (He et al. 2023) in magenta.The maroon line shows the results from the EAGLE simulation at z ∼ 5 (Rosas-Guevara et al. 2016).The solid and dashed blue lines show the result from Delphi(Dayal et al. 2014(Dayal et al. , 2019(Dayal et al. , 2020) ) simulations for all and bright (L bol ≳ 10 44 erg/s) black holes, respectively.Measured number densities of our LRDs agree well with the spectroscopic sample from(Matthee et al. 2023) and simulations given a λ edd ∼ 1.

Table 1 .
Properties of the observed fields with JWST /NIRCam observations.
To do that we fit our sources with pysersic(Pasha & , derived from observed rest-UV light.Upper limits are shown is downward pointing arrows.The dashed red line and the shaded area correspond to our best-fit Schechter function and its 68 % confidence interval, respectively.Vertical maroon lines highlight the MUV completeness limit calculated from the average depth of F814W/F090W bands.We compare our derived number densities to the luminosity functions of Lyman break galaxies from Bouwens et al.

Table 2 .
Bolometric and UV (λrest =1450 Å ) luminosity functions, as well as a SMBH mass function for our sample of LRDs.