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 JWSTNIRCam fields totaling 640 arcminSUP2SUP. Our careful selection yields 260 reddened AGN candidates at 4 lt z SUBphotSUB lt 9 dominated by a pointsourcelike central component r SUBeffSUB lt 130 pc and displaying a dichotomy in their restframe colors blue UV and red optical slopes. Quasar model fitting reveals our objects to be moderately dustextincted A SUBVSUB  1.6 which is reflected in their inferred bolometric luminosities of L SUBbolSUB  10SUP4447SUP erg sSUP1SUP and fainter UV magnitudes M SUBUVSUB  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 UVselected quasars of similar magnitudes. At the same time they constitute only a small fraction of all UVselected galaxies at similar redshifts but this percentage rises to 10 for M SUBUVSUB   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.
['Galaxies', 'High-redshift galaxies', 'Active galaxies', 'Active galactic nuclei', '573', '734', '17', '16', 'Astrophysics - Astrophysics of Galaxies']
A Census of Photometrically Selected Little Red Dots at 4 lt z lt 9 in JWST Blank Fields
{'No Header': '7', 'A Census of Photometrically Selected Little Red Dots at 4 < z < 9 in JWST Blank Fields': "Vasily Kokorev, 1 Karina I. Caputi, 1, 2 Jenny E. Greene, 3 Pratika Dayal, 1 Maxime Trebitsch, 1 Sam E. Cutler, 4 Seiji Fujimoto, 5, 2 Ivo Labb'e, 6 Tim B. Miller, 7 Edoardo Iani, 1 Rafael Navarro-Carrera, 1 and Pierluigi Rinaldi 1  \n1  \nKapteyn Astronomical Institute, University of Groningen, 9700 AV Groningen, The Netherlands 2 Cosmic Dawn Center (DAWN), Niels Bohr Institute, University of Copenhagen, Jagtvej 128, København N, DK-2200, Denmark 3 Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544 4 Department of Astronomy, University of Massachusetts, Amherst, MA 01003, USA 5 Department of Astronomy, The University of Texas at Austin, Austin, TX 78712, USA 6 Centre for Astrophysics and Supercomputing, Swinburne University of Technology, Melbourne, VIC 3122, Australia Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), Northwestern University, 1800 Sherman Ave, Evanston IL 60201, USA  \n(Received n/a; Revised n/a; Accepted n/a)  \nSubmitted to ApJ", 'ABSTRACT': '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 co-evolution 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 arcmin 2 . 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 = 10 44 -47 erg/s, 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.  \nKeywords: Active galactic nuclei (16), High-redshift galaxies (734), Early universe (435)', '1. 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 highz sources (e.g. Bunker et al. 2023; Maiolino et al.  \nCorresponding author: Vasily Kokorev  \nkokorev@astro.rug.nl  \n2023a) 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).  \nWhat 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 highz galaxies at the bright end ( M UV ≤ -20) of the  \nUV luminosity function is in tension with current theoretical frameworks (Behroozi & Silk 2015; Dayal et al. 2017; Yung et al. 2019, 2020; Behroozi et al. 2019, 2020; Dav'e 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).  \nAlthough 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'e 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; Ubler et al. 2023).  \nDubbed 'little red dots' (LRDs), these sources have SEDs characterized by a unique 'v-shaped' continuum combined with their point source morphology (Labb'e 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'e 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'e 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.  \nExtending the selection of this compact red population of low-luminosity broad-line AGN candidates to  \nlarger 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 highz 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).  \nIn 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.  \nThroughout this work we assume a flat ΛCDM cosmology (e.g. Planck Collaboration et al. 2020) with Ω m , 0 = 0 . 3, Ω Λ , 0 = 0 . 7 and H 0 = 70 km s -1 Mpc -1 , and a Chabrier (2003) initial mass function (IMF) between 0 . 1 -100 M ⊙ . All magnitudes are expressed in the AB system (Oke 1974).", '2. 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 GOODSS 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, 1287 PIs: 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.', '2.1. JWST 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  \nFigure 1. 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 from Labb'e 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.  \n<!-- image -->  \ne.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 NIRCam filters. These reductions are publicly available as a part of the DAWN JWST Archive (DJA 1 ).  \nThese 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  \nreference 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.  \nSome 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-  \niately 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.", '2.2. 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, 2023b), so we believe it to be adequate for our sources.  \nAdditionally, 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.", '3. 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  \nFigure 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), highz 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 M BH = 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 A V ∼ 1 . 6 +1 . 1 -1 . 0 . On the side of each panel we also show histograms highlighting the redshift and L bol distributions.  \n<!-- image -->  \nvariance and enable us to avoid dealing with the cosmic volume uncertainties introduced by the lensing magnification.', '3.1. Color and Morphology Selection': "Recently, Labb'e et al. (2023b) published a large sample of photometrically identified compact red sources from the Cycle 1 JWST UNCOVER program (PIs: I. Labb'e, R. Bezanson; Bezanson et al. 2022). Subsequent follow-up of 17 such objects with NIRSPec/MSA  \nTable 1. Properties of the observed fields with JWST /NIRCam observations.  \nNote -NIRCam depths: expressed as 5 σ within the 0 . '' 36 apertures used for the photometric extraction in the area covered by F115W/F150W/F200W/F277W/ F356W/F444W.  \nPRISM has resulted in a remarkable success rate of 83 %, with 14/17 photometrically selected targets confirmed as broad-line (BL) AGN at 4 < z < 8 . 5 (Greene et al. 2023), and 3/17 as brown dwarfs (Burgasser et al. 2023). In brief, the color cuts introduced in Labb'e et al. (2023b) are designed to catch the break between the red continuum slope in rest-frame optical, and the blue rest-UV emission ( λ rest ∼ 4000 ˚ A). This color selection requires that the red continuum slope is rising in more than one adjacent filter pair, to avoid selecting galaxies with strong emission lines. Indeed, currently available spectra of LRDs (e.g. Fujimoto et al. 2023a; Furtak et al. 2023b; Kocevski et al. 2023; Kokorev et al. 2023a; Greene et al. 2023; Matthee et al. 2023) display a remarkable dichotomy in their observed spectral shapes. In particular the SEDs at 1 -2 µ m (1000 - 2000 ˚ A rest) are blue ( f λ ∝ λ -2 ) and red ( f λ ∝ λ 0 -2 ) at 3 -5 µ m (3100 - 5200 ˚ A rest). As such we keep the Labb'e et al. color criteria largely unchanged, only introducing some further adjustments based on the UNCOVER spectra of LRDs, namely to limit the contamination of our sample by brown dwarfs as was suggested in Greene et al. (2023).  \nColors alone would end up selecting both LRDs and extended red galaxies (see e.g. Labb'e 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-  \nr cuts are then:  \nred1 = F115W -F150W < 0 . 8 & F200W -F277W > 0 . 7 & F200W -F356W > 1 . 0  \nor red2 = F150W -F200W < 0 . 8 & F277W -F356W > 0 . 6 & F277W -F444W > 0 . 7 ,  \nwhich are effectively selecting our low ( z < 6) and high ( z > 6) redshift samples, respectively. The compactness is given by:  \ncompact = f f444w (0 . '' 4) /f f444w (0 . '' 2) < 1 . 7 .  \nTo limit the number of brown dwarfs in the sample we also adopt:  \nbd removal = F115W -F200W > -0 . 5 .  \nThe final selection then becomes ( red1 | red2 ) & compact & bd removal . 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 specz for sub-samples in the next sub section.", '3.2. 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. To do that we fit our sources with pysersic (Pasha &  \nFigure 3. Observed M UV compared to the M UV expected from the dust-corrected L bol . We derive the expected M UV 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 M UV we derive is fainter than the expected value, with A UV 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 A V .  \n<!-- image -->  \nMiller 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'e 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.  \nTaking 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  \nlevels within 4' diameter apertures 2 . For more detail see the Appendix in Weaver et al. (2023a).  \nThe light is modeled with a single S'ersic (S'ersic 1963) profile with the center, brightness, effective radius, S'ersic 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.  \nA 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.", '3.3. 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 subset 3 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).  \nWhile it might seem counter-intuitive to use galaxy templates for what we believe to be AGN candidates, similar efforts presented in Labb'e et al. (2023b) report a good agreement between deriving z phot with stellar templates alone, as opposed to stellar+AGN models, finding  \na 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 ˚ A) and Balmer ( ∼ 4000 ˚ A) breaks.  \nFor 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 ˚ A (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'e et al. 2023a; Steinhardt et al. 2023).  \nSpectroscopic 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.  \nWe 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 highz 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'e 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  \nbroad 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.", '3.4. 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'e et al. 2013, 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.  \nUnfortunately, 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  \n( f λ ∝ λ 0 -2 at 3100 - 5200 ˚ A rest) slopes, the rest optical continuum in our sources is also dominated by AGN light.  \nIn 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'e 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.  \nFollowing 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.  \nThe same approach was already successfully employed in Labb'e 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.  \nWe 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  \nfor 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'e 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.", '3.5. 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 ˚ A. 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 ˚ A. 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.  \nTo 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 ˚ A, and add back the continuum. We then re-fit all of our sources, following the same considerations as described in Section 3.4.  \nUsing 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.', '3.6. 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 highz 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 ˚ A. 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 ˚ A is essential to accurately derive M UV , and the 5100 ˚ A 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 ˚ A 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.", '3.7. 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'e 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.  \nUsing 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 ˚ A continuum, measured directly from  \nbest-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'e 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.  \nIn 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.  \nThe final table which contains photometry, sizes and the physical parameters we derive for our sample is available in full online 4 . We show an excerpt of the full table in the Appendix.", '4.1. 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'e 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.  \nFigure 4. The UV luminosity function (UVLF) for the LRDs in our sample in the z ∈ [4 . 5 , 6 . 5] (left) and z ∈ [6 . 5 , 8 . 5] (right) bins (blue), 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 M UV 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. (2021) (solid blue line), extrapolated quasar UVLF relations from Niida et al. (2020) at z ∼ 5, as well as an upper bound provided by Kulkarni et al. (2019) (green lines). At z ∼ 7, in green, we show the UVLF derived from bright quasars from Matsuoka et al. (2023). We highlight the spectroscopically identified LRDs from Greene et al. (2023) (blue squares), Matthee et al. (2023) (orange circles) and Kocevski et al. (2023) (blue pentagon). Furthermore, we show densities of BL AGN quasars from Maiolino et al. (2023b) (gray triangles) and Harikane et al. (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.  \n<!-- image -->  \nFocusing 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:  \nΦ( x ) = 1 ∆ x ∑ i V max , i ( A,z min , z max ) -1 ,  \nwhere ∆ 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.  \nWe 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  \nto 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.  \nAccounting 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'e 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.", '4.2. 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 ˚ A 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.  \nAt 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.  \nBefore 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).  \nTaking 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'e et al. 2023b; Matthee et al. 2023), at least at M UV ∼ -19 and brighter. Confirming the initial findings for the UNCOVER red-dots presented in Greene et al. (2023) and Labb'e et al. (2023b), we also find that our sample accounts for ∼ 10 - 30 % of total BL AGN populations at highz (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.  \nWhen 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 UNCOVER 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 ˚ A 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.  \nFollowing 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  \nreport 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.  \nProvided that our color and morphology selection is comparably successful at identifying reddened AGN as was previously shown (Labb'e 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 preJWST 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.", '4.3. Bolometric Luminosity Function': "Our SED fitting results show that the fraction of the UV light contributing to the total luminosity is small as  \na 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 .  \nUnderstanding 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.  \nOur 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'e 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-  \nFigure 5. Bolometric luminosity functions in the z ∈ [4 . 5 , 6 . 5] (left) and z ∈ [6 . 5 , 8 . 5] (right) bins, derived from L 5100 , assuming rest-frame optical continuum is AGN dominated. The number densities have been V max 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 from Greene et al. (2023). In addition in the lowest redshift bin we show the NIRCam grism result of Matthee et al. (2023) (open circle). Dashed lines show the preJWST L bol relation derived in Shen 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.  \n<!-- image -->  \ne, 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.  \nWe additionally compare our bolometric luminosity function to the latest version of the semi-analytic Delphi models (Dayal et al. 2019, 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.  \nWe 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  \nobjects 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).  \nFinally, 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 photoz (Marchesini et al. 2009) uncertainties.", '4.4. The z ∼ 5 SMBH Mass Function': 'With the data we have obtained, we will now derive and describe the measurement for the supermassive BH mass function which our compact objects potentially host. Computing the mass of the central black hole generally requires knowledge of the width of the broad  \nFigure 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, 2019, 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.  \n<!-- image -->  \nlines (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 ˚ A (see e.g. Kaspi et al. 2000; Greene & Ho 2005).  \nWhile 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.  \nTable 2. Bolometric and UV ( λ rest =1450 ˚ A) luminosity functions, as well as a SMBH mass function for our sample of LRDs.  \nWecalculate 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  \nrange. 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.  \nWe 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.  \nFurthermore, 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..  \nNaturally, 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.  \nBefore 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, 2019, 2020, 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  \nTable 3. Best fit Schechter parameters for the rest-frame UVLF at λ rest = 1450 ˚ A, across blank JWST fields.  \nerg/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.  \nExamining 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 BHpopulation 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.', '5.1. Abundance of bright compact sources': "Previously limited to UV -selected samples at z ≲ 6 (Kashikawa et al. 2015; Ba˜nados 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.  \n2023; Matthee et al. 2023; Ubler 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'e et al. 2023b; Matthee et al. 2023).  \nWhile 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'e 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.  \nMotivated 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.  \nUsing 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 highz 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  \npinned 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).  \nAssuming 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.", '5.2. 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'e 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.  \nFirstly, 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  \nassessment 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 ˚ A), the CIV, SiIV or HeII ( λ rest ∼ 1840 ˚ A) 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.  \nSecondly, 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 highz . 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.  \nThirdly, it is crucial to note that a substantial proportion of massive SMBHs (with M BH > 10 8 M ⊙ ) at high redshifts (highz ) 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.  \nFinally, as was already shown recently by Williams et al. (2023a) and P'erez-Gonz'alez 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.  \nEarly 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 highz in an effort to better understand their properties and  \nabundance. 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.", 'ACKNOWLEDGMENTS': "We thank the anonymous referee for a number of constructive suggestions, which helped to greatly improve the quality of this manuscript. We are grateful to Dale Kocevski and Kohei Inayoshi for their patience with helping us spot and correct minor inconsistencies in the manuscript. The authors would like to thank Sarah Bosman for insightful discussions about UV-bright quasars at high redshift. We also thank Mauro Giavalisco, Hollis Akins and Meghana Killi for useful discussion regarding the nature of dust obscured AGN. VK and KIC acknowledge funding from the Dutch Research Council (NWO) through the award of the Vici Grant VI.C.212.036. J.E.G. acknowledges support from NSF/AAG grant# 1007094, and also support from NSF/AAG grant # 1007052. PD & MT acknowledge support from the NWO grant 016.VIDI.189.162 ('ODIN'). PD also acknowledges support from the European Commission's and University of Groningen COFUND Rosalind Franklin program. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. All the JWST and HST data used in this paper can be found on MAST: 10.17909/de9v-7893. Some of the data products presented herein were retrieved from the Dawn JWST Archive (DJA). DJA is an initiative of the Cosmic Dawn Center, which is funded by the Danish National Research Foundation under grant No. 140. TBM was supported by a CIERA Postdoctoral Fellowship. This work used computing resources provided by Northwestern University and the Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA). This research was supported in part through the computational resources and staff contributions provided for the Quest high performance computing facility at Northwestern University which is jointly supported by the Office of the Provost, the Office for Research, and Northwestern University Information Technology.  \n```\nSoftware: EAZY(Brammer et al. 2008), FSPS (Conroy et al. 2009), pysersic (Pasha & Miller 2023), grizli (Brammer 2023), msaexp (Brammer 2022). Facilities: JWST , HST\n```", 'APPENDIX': "Table 4. An example of the table containing all properties of our sources. A full version of this table is available in the electronic format.  \nNote -Sizes are measured in F444W band on the 0 . '' 04 images. The FWHM of the F444W PSF is 3.45 pixels.", 'REFERENCES': "Adams, N. J., Conselice, C. J., Ferreira, L., et al. 2023, MNRAS, 518, 4755, doi: 10.1093/mnras/stac3347 Aird, J., Coil, A. L., Georgakakis, A., et al. 2015, MNRAS, 451, 1892, doi: 10.1093/mnras/stv1062 Akins, H. B., Casey, C. M., Allen, N., et al. 2023, ApJ, 956, 61, doi: 10.3847/1538-4357/acef21 Andika, I. T., Jahnke, K., Onoue, M., et al. 2024, arXiv e-prints, arXiv:2401.11826. https://arxiv.org/abs/2401.11826 Atek, H., Chemerynska, I., Wang, B., et al. 2023, MNRAS, 524, 5486, doi: 10.1093/mnras/stad1998 Austin, D., Adams, N., Conselice, C. J., et al. 2023, ApJL, 952, L7, doi: 10.3847/2041-8213/ace18d Ba˜nados, E., Carilli, C., Walter, F., et al. 2018, ApJL, 861,  \nL14, doi: 10.3847/2041-8213/aac511  \nBaggen, J. F. W., van Dokkum, P., Labb'e, I., et al. 2023, ApJL, 955, L12, doi: 10.3847/2041-8213/acf5ef  \nBagley, M. B., Finkelstein, S. L., Koekemoer, A. M., et al.  \n2022, arXiv e-prints, arXiv:2211.02495.  \nhttps://arxiv.org/abs/2211.02495  \nBanerji, M., Alaghband-Zadeh, S., Hewett, P. C., & McMahon, R. G. 2015, MNRAS, 447, 3368,  \ndoi: 10.1093/mnras/stu2649  \nBarbary, K. 2016, Journal of Open Source Software, 1, 58, doi: 10.21105/joss.00058  \nBarro, G., Perez-Gonzalez, P. G., Kocevski, D. D., et al. 2023, arXiv e-prints, arXiv:2305.14418,  \ndoi: 10.48550/arXiv.2305.14418  \nBehroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143, doi: 10.1093/mnras/stz1182  \nFurtak, L. J., Zitrin, A., Weaver, J. R., et al. 2023c, MNRAS, 523, 4568, doi: 10.1093/mnras/stad1627  \nGaia Collaboration, Brown, A. G. A., Vallenari, A., et al.  \n2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657  \nGehrels, N. 1986, ApJ, 303, 336, doi: 10.1086/164079  \n- Giallongo, E., Grazian, A., Fiore, F., et al. 2019, ApJ, 884, 19, doi: 10.3847/1538-4357/ab39e1\n- Glikman, E., Helfand, D. J., & White, R. L. 2006, ApJ, 640, 579, doi: 10.1086/500098\n- Glikman, E., Simmons, B., Mailly, M., et al. 2015, ApJ, 806, 218, doi: 10.1088/0004-637X/806/2/218\n- Gordon, K. D., Clayton, G. C., Misselt, K. A., Landolt, A. U., & Wolff, M. J. 2003, ApJ, 594, 279,  \ndoi: 10.1086/376774  \nGoulding, A. D., Greene, J. E., Setton, D. J., et al. 2023,  \narXiv e-prints, arXiv:2308.02750,  \ndoi: 10.48550/arXiv.2308.02750  \n- Grazian, A., Giallongo, E., Boutsia, K., et al. 2018, A&A, 613, A44, doi: 10.1051/0004-6361/201732385\n- Greene, J. E., & Ho, L. C. 2005, ApJ, 630, 122,  \ndoi: 10.1086/431897  \n- Greene, J. E., Labbe, I., Goulding, A. D., et al. 2023, arXiv e-prints, arXiv:2309.05714,  \ndoi: 10.48550/arXiv.2309.05714  \nHarikane, Y., Zhang, Y., Nakajima, K., et al. 2023, arXiv e-prints, arXiv:2303.11946,  \ndoi: 10.48550/arXiv.2303.11946  \nHe, W., Akiyama, M., Enoki, M., et al. 2023, arXiv e-prints, arXiv:2311.08922, doi: 10.48550/arXiv.2311.08922  \nInayoshi, K., Visbal, E., & Haiman, Z. 2020, ARA&A, 58, 27, doi: 10.1146/annurev-astro-120419-014455  \n- Jin, S., Sillassen, N. B., Magdis, G. E., et al. 2022, arXiv e-prints, arXiv:2212.09372.\n- https://arxiv.org/abs/2212.09372\n- -. 2023, A&A, 670, L11,  \ndoi: 10.1051/0004-6361/202245724  \nKannan, R., Springel, V., Hernquist, L., et al. 2023,  \nMNRAS, 524, 2594, doi: 10.1093/mnras/stac3743  \nKartaltepe, J. S., Rose, C., Vanderhoof, B. N., et al. 2023, ApJL, 946, L15, doi: 10.3847/2041-8213/acad01  \nKashikawa, N., Ishizaki, Y., Willott, C. J., et al. 2015, ApJ, 798, 28, doi: 10.1088/0004-637X/798/1/28  \nKaspi, S., Smith, P. S., Netzer, H., et al. 2000, ApJ, 533, 631, doi: 10.1086/308704  \nKato, N., Matsuoka, Y., Onoue, M., et al. 2020, PASJ, 72, 84, doi: 10.1093/pasj/psaa074  \nKilli, M., Watson, D., Brammer, G., et al. 2023, arXiv e-prints, arXiv:2312.03065,  \ndoi: 10.48550/arXiv.2312.03065  \nWilliams, C. C., Alberts, S., Ji, Z., et al. 2023a, arXiv e-prints, arXiv:2311.07483, doi: 10.48550/arXiv.2311.07483 Williams, H., Kelly, P. L., Chen, W., et al. 2023b, Science, 380, 416, doi: 10.1126/science.adf5307 Yang, G., Caputi, K. I., Papovich, C., et al. 2023, ApJL,  \n950, L5, doi: 10.3847/2041-8213/acd639  \nYang, J., Wang, F., Fan, X., et al. 2021, ApJ, 923, 262, doi: 10.3847/1538-4357/ac2b32 Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., Popping, G., & Dav'e, R. 2019, MNRAS, 483, 2983, doi: 10.1093/mnras/sty3241 Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., et al. 2020, MNRAS, 496, 4574, doi: 10.1093/mnras/staa1800"}