Lowluminosity active galactic nuclei AGNs with lowmass black holes BHs in the early universe are fundamental to understanding the BH growth and their coevolution with the host galaxies. Utilizing JWST NIRCam Wide Field Slitless Spectroscopy we perform a systematic search for broadline H emitters BHAEs at z  45 in 25 fields of the A SPectroscopic survey of biased halos In the Reionization Era ASPIRE project covering a total area of 275 arcminSUP2SUP. We identify 16 BHAEs with FWHM of the broad components spanning from 1000 to 3000 km sSUP1SUP. Assuming that the broad line widths arise as a result of Doppler broadening around BHs the implied BH masses range from 10SUP7SUP to 10SUP8SUP M SUBSUB with broad Hconverted bolometric luminosities of 10SUP44.5SUP10SUP45.5SUP erg sSUP1SUP and Eddington ratios of 0.070.47. The spatially extended structure of the F200W stacked image may trace the stellar light from the host galaxies. The H luminosity function indicates an increasing AGN fraction toward the higher H luminosities. We find possible evidence for clustering of BHAEs two sources are at the same redshift with a projected separation of 519 kpc one BHAE appears as a composite system residing in an overdense region with three close companion H emitters. Three BHAEs exhibit blueshifted absorption troughs indicative of the presence of high column density gas. We find that the broadlineselected and photometrically selected BHAE samples exhibit different distributions in the optical continuum slopes which can be attributed to their different selection methods. The ASPIRE broadline H sample provides a good database for future studies of faint AGN populations at high redshift.
['High-redshift galaxies', 'Active galactic nuclei', '734', '16', 'Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Cosmology and Nongalactic Astrophysics']
A SPectroscopic Survey of Biased Halos In the Reionization Era ASPIRE Broadline AGN at z  45 Revealed by JWSTNIRCam WFSS
{'A SPectroscopic survey of biased halos In the Reionization Era (ASPIRE): Broad-line AGN at z = 4 -5 revealed by JWST/NIRCam WFSS': "Xiaojing Lin , 1, 2 Feige Wang , 2 Xiaohui Fan , 2 Zheng Cai , 1 Jaclyn B. Champagne , 2 Fengwu Sun , 2, 3 Marta Volonteri , 4 Jinyi Yang , 2 Joseph F. Hennawi , 5, 6 Eduardo Ba˜nados , 7 Aaron Barth , 8 Anna-Christina Eilers , 9 Emanuele Paolo Farina, 10 Weizhe Liu , 2 Xiangyu Jin, 2 Hyunsung D. Jun , 11 Alessandro Lupi , 12, 13, 14 Koki Kakiichi, 15, 16 Chiara Mazzucchelli , 17 Masafusa Onoue , 18, 19, 20 Zhiwei Pan, 21, 22 Elia Pizzati, 23 Sof'ıa Rojas-Ruiz , 24 Jan-Torge Schindler , 25 Benny Trakhtenbrot, 26 Yue Shen , 27, 28 Maxime Trebitsch, 29 Ming-Yang Zhuang , 30 Ryan Endsley, 31 Romain A. Meyer, 7 Zihao Li, 1 Mingyu Li, 1 Maria Pudoka, 2 Wei Leong Tee, 2 Yunjing Wu, 1 and Haowen Zhang 2  \n1 Department of Astronomy, Tsinghua University, Beijing 100084, China  \n2  \nSteward Observatory, University of Arizona, 933 N Cherry Ave, Tucson, AZ 85721, USA  \nCenter for Astrophysics  \n|  \nHarvard & Smithsonian, 60 Garden St., Cambridge MA 02138 USA  \n3  \n- 4 Institut d'Astrophysique de Paris, Sorbonne Universit'e, CNRS, UMR 7095, 98 bis bd Arago, 75014 Paris, France\n- 5 Department of Physics, Broida Hall, University of California, Santa Barbara, CA 93106-9530, USA  \n6  \nLeiden Observatory, Leiden University, P.O. Box 9513, NL-2300 RA Leiden, The Netherlands  \n- 7 Max Planck Institut fur Astronomie, Konigstuhl 17, D-69117, Heidelberg, Germany\n- 8 Department of Physics and Astronomy, 4129 Frederick Reines Hall, University of California, Irvine, CA 92697-4575, USA 9 MIT Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA 10 Gemini Observatory, NSF's NOIRLab, 670 North A'ohoku Place, Hilo, HI 96720, USA\n- 11 Department of Physics, Northwestern College, 101 7th St SW, Orange City, IA 51041, USA  \n12  \nDipartimento di Scienza e Alta Tecnologia, Universit'a degli Studi dell'Insubria, via Valleggio 11, I-22100, Como, Italy  \n13  \nINFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy  \n- 14 Dipartimento di Fisica 'G. Occhialini', Universit'a degli Studi di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy 15 Cosmic Dawn Center (DAWN), Denmark  \n16 Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen N, Denmark  \n- 17 Instituto de Estudios Astrof'ısicos, Facultad de Ingenier'ıa y Ciencias, Universidad Diego Portales, Avenida Ejercito Libertador 441, Santiago, Chile.\n- 18 Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU, WPI), The University of Tokyo Institutes for Advanced Study, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan\n- 19 Center for Data-Driven Discovery, Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan 20 Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, P.R.China\n- 21 Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, People's Republic of China  \n22  \nDepartment of Astronomy, School of Physics, Peking University, Beijing 100871, People's Republic of China  \n23  \nLeiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands  \nDepartment of Physics and Astronomy, University of California, Los Angeles, 430 Portola Plaza, Los Angeles, CA 90095, USA  \n25  \nHamburger Sternwarte, Universitat Hamburg, Gojenbergsweg 112, D-21029 Hamburg, Germany  \n26 School of Physics and Astronomy, Tel Aviv University, Tel Aviv 69978, Israel  \n27 Department of Astronomy, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA  \nNational Center for Supercomputing Applications, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA  \n29 Kapteyn Astronomical Institute, University of Groningen, 9700 AV Groningen, The Netherlands  \n30  \nDepartment of Astronomy, University of Illinois Urbana-Champaign, Urbana, IL 61801, USA  \n31 Department of Astronomy, University of Texas, Austin, TX 78712, USA  \n(Received XXX; Revised YYY; Accepted ZZZ)  \nSubmitted to AJ  \nlinxj21@mails.tsinghua.edu.cn  \nCorresponding author: Xiaojing Lin  \n24  \n28", 'ABSTRACT': 'Low-luminosity AGNs with low-mass black holes (BHs) in the early universe are fundamental to understanding the BH growth and their co-evolution with the host galaxies. Utilizing JWST NIRCam Wide Field Slitless Spectroscopy (WFSS), we perform a systematic search for broad-line H α emitters (BHAEs) at z ≈ 4 -5 in 25 fields of the ASPIRE (A SPectroscopic survey of biased halos In the Reionization Era) project, covering a total area of 275 arcmin 2 . We identify 16 BHAEs with FWHM of the broad components spanning from ∼ 1000 km s -1 to 3000 km s -1 . Assuming the broad linewidths arise due to Doppler broadening around BHs, the implied BH masses range from 10 7 to 10 8 M ⊙ , with broad H α -converted bolometric luminosity of 10 44 . 5 -10 45 . 5 erg s -1 and Eddington ratios of 0 . 07 -0 . 47. The spatially extended structure of the F200W stacked image may trace the stellar light from the host galaxies. The H α luminosity function indicates an increasing AGN fraction towards the higher H α luminosities. We find possible evidence for clustering of BHAEs: two sources are at the same redshift with a projected separation of 519 kpc; one BHAE appears as a composite system residing in an overdense region with three close companion H α emitters. Three BHAEs exhibit blueshifted absorption troughs indicative of the presence of high-column-density gas. We find the broad-line and photometrically selected BHAE samples exhibit different distributions in the optical continuum slopes, which can be attributed to their different selection methods. The ASPIRE broad-line H α sample provides a good database for future studies of faint AGN populations at high redshift.  \nKeywords: high-redshift - AGN - Black hole assembly', '1. INTRODUCTION': "Tremendous efforts have been made in the past decades to search for and characterize supermassive black holes (SMBHs) in distant bright quasars and active galactic nuclei (AGNs; Fan et al. 2023). Luminous quasars at z > 5, and BH mass of ∼ 10 9 M ⊙ , pose challenges to BH formation theories, with rapid mass assembly within a short amount of cosmic time (e.g., Ba˜nados et al. 2018; Wang et al. 2021; Yang et al. 2021). Such a rapid growth with sustained episodes would erase the imprints of initial BH seeds (Volonteri 2010), necessitating studies of high redshift low-mass BHs and low-luminosity AGNs. However, the knowledge of early low-mass BHs is very limited before the launch of James Webb Space Telescope (JWST). With the unprecedented infrared capabilities of JWST, we are now able to their early, low-mass, low-luminosity systems, closer to the early stage of BH growth and the predicted seed BH population (Inayoshi et al. 2022). Faint AGNs at high redshift with smaller BH provide crucial constraints to the BH seeding, growth, and coevolution with galaxies (Pacucci et al. 2023; Volonteri et al. 2023; Ding et al. 2023; Li et al. 2024).  \nIn pursuit of the low-luminosity AGNs (e.g., Onoue et al. 2023; Ubler et al. 2023; Furtak et al. 2023; Kocevski et al. 2023), JWST reveals an abundant population of red compact sources at z > 4 (e.g., Barro et al. 2023; Labbe et al. 2023), nicknamed as 'little red dots' (Matthee et al. 2023, LRDs). Their number density is 10 -100 × higher than the extrapolation from quasars  \nluminosity functions (e.g., Kocevski et al. 2023). The nature of LRDs remains a mystery. JWST spectroscopy detects moderately luminous and broad H α or Hβ lines in the LRDs, suggesting SMBHs with masses ranging from ∼ 10 6 to ≲ 10 9 M ⊙ (Kokorev et al. 2023; Kocevski et al. 2023; Furtak et al. 2023; Maiolino et al. 2023; Matthee et al. 2023; Greene et al. 2024; Harikane et al. 2023). Most of them have SMBH masses 1-2 dex lower than those of UV bright quasars at similar redshifts identified before the launch of JWST (Onoue et al. 2019; Yang et al. 2023). LRDs typically exhibit V-shape spectral energy distribution (SED, Labbe et al. 2023; Furtak et al. 2023; Greene et al. 2024; Killi et al. 2023), with the optical continuum redder than that in the UV in the f λ space (typically the UV continuum slope β UV < 0, the optical continuum slope β opt > 1). Whether the UV emission of LRDs originates from host galaxies or the scattered light of AGNs is still under debate (Labbe et al. 2023; Kocevski et al. 2023; P'erez-Gonz'alez et al. 2024).  \nHowever, before a full understanding of the physical nature of LRDs, one of the main areas of uncertainty is the selection criteria. Spectroscopic selections (Kocevski et al. 2023; Matthee et al. 2023; Maiolino et al. 2023; Harikane et al. 2023) are based on the broad H α or H β emission lines, indicative of broad-line regions around BHs. On the other hand, photometric selections are usually based on LRDs' V-shape SED and their compactness and redness (Labbe et al. 2023; Greene et al. 2024; Kokorev et al. 2024; P'erez-Gonz'alez et al. 2024). Some photometrically selected LRD candidates have been con-  \ned as brown dwarfs spectroscopically (Greene et al. 2024; Hainline et al. 2023). Without showing broadened Balmer emission lines or other AGN line diagnostics (Scholtz et al. 2023), we also cannot rule out the possibility of dusty starburst galaxies. JWST /MIRI studies on photometrically selected LRDs show flattened spectral shape in the rest-frame near-infrared ( ≳ 1 µ m), inconsistent with emission from AGN dust tori (Williams et al. 2023; P'erez-Gonz'alez et al. 2024). P'erez-Gonz'alez et al. (2024) claimed that these photometrically selected LRDs are more likely to be extremely intense dusty starburst galaxies with obscured AGNs contributing in the mid-infrared. It is unclear whether the photometrically selected LRDs and those broad-line selected ones through WFSS or NIRSpec share the same physical properties. Furthermore, it also deserves further exploration to confirm if these LRDs are indeed AGNs.  \nTo better constrain the nature of low-luminosity AGNs, including the so-called LRDs, and study BH mass assembly at high redshifts, we need a large sample of broad-line selected faint AGNs. In this work, we perform a systemic search for broad H α line emitters in the ASPIRE program (Wang et al. 2023; Yang et al. 2023). ASPIRE is the JWST Cycle 1 grism program covering a large sky area, totaling ≈ 65 hours over ≈ 275 arcmin 2 , over 25 differnt fields. Using JWST /NIRCam WFSS (Greene et al. 2016) in F356W ( R ≈ 1600), we can build up a flux-limited broad H α emitter sample at z = 4 -5. ASPIRE's large survey area over multiple fields allows quantifying the number density and AGN fraction, effectively mitigating the cosmic variance. To avoid confusion, we refer to our broad-line selected sample as broad H α emitters throughout this paper (BHAEs for short), rather than LRDs.  \nThis paper is organized as follows. We describe the JWST observations, data reduction, and sample selection in § 2. In § 3 we present the properties of our selected broad H α emitter sample, including their spectral properties ( § 3.1), BH masses and bolometric luminosities ( § 3.2), UV and optical continuum slopes ( § 3.3 and morphology ( § 3.4). We analyze the H α luminosity function in § 3.5. We discuss interesting individual sources in our sample in § 4. Our results are summarized in § 5. Throughout this work, a flat ΛCDM cosmology is assumed, with H 0 = 70 km s -1 Mpc -1 , Ω Λ , 0 = 0 . 7 and Ω m, 0 = 0 . 3. Distances are provided in the comoving frame if not specified.", '2.1. JWST observations': "The ASPIRE (A SPectroscopic survey of biased halos In the Reionization Era) program (GO 2078; PI  \nF. Wang) targets 25 quasars at 6 . 5 < z < 6 . 8 with JWST /NIRCam imaging and wide field slitless spectroscopic (WFSS) data, with a total survey area of ∼ 11 × 25 arcmin 2 . For each field, ASPIRE obtains 2834.5s exposures in NIRCam/F356W Grism R, accompanied with direct imaging of 1417.3s in NIRCam/F115W and F356W, and 2800s in F200W. The data is reduced in the same way as described in Wang et al. (2023). We summarize the main procedure below.  \nThe NIRCam images are reduced with the standard pipeline v1.10.2 1 . We utilize the reference files ( jwst 1080.pmap ) from version 11.16.21 of the standard Calibration Reference Data System 2 (CRDS) for calibration. The 1/f noise is subtracted on a row-by-row basis and column-by-column basis after CALWEBB Stage 1. After Stage 2, a master median background is generated using all available exposures from ASPIRE for each band. The master backgrounds are then scaled and subtracted for individual exposures. All the images are mosaicked with a pixel scale of 0.031 '' and pixfrac=0.8 . Their WCS are registered using the Gaia Data Release 2 catalog (Gaia Collaboration et al. 2016, 2018). we construct the photometric catalogs with SourceExtractor++(Bertin et al. 2020), using F356W as the detection images. We first measure the flux in elliptical apertures with a Kron factor of k = 1 . 2 (Kron 1980). Then we performed aperture correction for each source, multiplying the ratio of F356W flux within k = 2 . 5 v.s. k = 1 . 2 apertures. SourceExtractor++ also provides the F356W half-light radii measurements. More details about the photometric catalogs are presented in Champagne et al. in prep.  \nFor NIRCam WFSS data, we subtract the 1/f noise only along columns. The median background models for WFSS are constructed using all available ASPIRE observations obtained at similar times and are scaled and subtracted from individual WFSS exposures. The WCS of individual exposures is first aligned using the CALWEBB pipeline. We then measure the astrometric offsets between the direct images and the fully calibrated F356W-band mosaic, which are added to the spectral tracing model as described below. The spectral tracing models, dispersion models, and sensitivity functions are constructed as described in Sun et al. (2023). The spectral tracing model and the sensitivity functions used in this work are updated using the spectral traces of point sources observed in the LMC field (PID #1076) and Cycle-1 calibration programs (PIDs  \n#1536, #1537, #1538) respectively. we construct and stack the 2D spectra for each source, using 2D spectra from individual exposures that are registered to a common wavelength and spatial grid following the algorithms in PypeIt (Prochaska et al. 2020). The pixel scales for the 2D spectra and 1D spectra described below are both resampled to be 10 ˚ A per pixel.  \nWe first extract 1D spectra using the optimal extraction algorithms and remove the background and continuum using median filters (see Wang et al. 2023 for details). These continuum-removed 1D spectra are used for the initial emission line identification. After identifying the H α emitter candidates, we re-extract the 1D spectra for our analysis, without continuum removal. For those with contaminated 2D spectra, we mask regions of ± 5000 km s -1 around the H α emission to avoid over-subtraction around the broad wing. We then run median filters iteratively along each row of the 2D spectra to construct continuum maps. We re-extract 1D spectra from these continuum-removed 2D spectra using the optimal extraction algorithms. For samples whose 2D grism spectra are not contaminated by nearby bright sources, we re-extract the 1D spectra based on their original grism data.", '2.2. Sample selection': 'We systematically search for broad-line H α emitters in the 25 ASPIRE fields over a total area of 275 arcmin 2 . We select compact red sources first with the following morphology and color cut:  \n- (1) F356W half-light radius < 0.2 arcsec;\n- (2) F200W -F356W > 0.75 mag.  \nThe first requirement for compactness ensures that the source morphology does not broaden the H α line profiles during the dispersion of grism spectra. The second criterion describes the F356W excess boosted by the strong H α emission line. We define the color cut based on the broad-H α sample in Matthee et al. (2023). We then perform an emission line search for all the photometrically selected red compact sources and select objects with strong emission lines at a signal-to-noise ratio (S/N) of integrated line fluxes greater than five. Careful visual inspection is performed to remove artifacts in the imaging and grism data, false detection by the search algorithms, and misidentified lines from contamination of other sources along the dispersion direction. We then fit the line profiles of the emitter candidates with single and two-component Gaussian models. Broad H α emitters are identified in two cases: (i) A single Gaussian model can well describe the line profile and it has an FWHM > 1000 km s -1 ; or (ii) the two-component  \nmodel with an FWHM > 1000 km s -1 broad component performs better, resulting in | χ 2 red -1 | at least 0.5 smaller than that of single Gaussian models. The spectroscopic selection can be summarized as:  \n- (3) Emission lines with S/N of integrated line fluxes > 5;\n- (4) Robust broad components with FWHM broad > 1000 km s -1 .  \nTo test our selection criteria, we generate a mock spectrum assuming a source with an effective radius of 0.2 arcsec and a single emission line with an intrinsic FWHM of 500 km s -1 . The mock emission line in the dispersed spectrum, broadened by both the line spread function of NIRCam WFSS and the source morphology, displays a Gaussian profile with FWHM of 715 km s -1 . This suggests that with flux radii < 0.2 arcsec and FWHMs of H α broad components > 1000 km s -1 , the H α of our selected BHAE should have intrinsically broad lines.  \nWe finally identify a total of 16 BHAEs in the 25 ASPIRE fields. Figure 1 presents our selection criteria and the 16 selected BHAEs. Their basic properties are listed in Table 1. We note that our selection does not rely on the spectral energy distribution (SED) shape of sources (e.g., Greene et al. 2024). We specifically select compact, point-source-like BHAEs, rather than selecting all BHAEs from a complete H α emitter (HAE) sample. The constraint of compactness maximizes the purity of our BHAE sample.', '3. ANALYSES AND RESULTS': 'In this section, we present our analyses and results on a sample of 16 ASPIRE BHAEs. We describe the properties of the broad H α emission lines in § 3.1 and derive the black hole masses and bolometric luminosities in § 3.2. The UV and optical continuum slopes are discussed in § 3.3. Finally, we study the morphology and H α luminosity function of BHAEs in § 3.4 and § 3.5.', '3.1. Broad H α emission': "As described in § 2.2, we fit the extracted 1D spectra with two-component Gaussian models for all 16 BHAEs. We do not impose any limits on the FWHM of the narrow-component Gaussian. Figure 2 and 3 display the 16 BHAEs and their best-fit line profiles. The measurements are presented in Table 2. Among them, J0430M1445-BHAE-1 and J0923P0402BHAE-1 have two broad components both with FWHM > 1000 km s -1 . The broad-line H α emission is often not adequately modeled as a single Gaussian (Shen  \nTable 1. The coordinate, redshift, and AB magnitude in F356W of the 16 LRDs in the ASPIRE fields. BHAE refers to Broad H α emitters throughout the paper. The redshift uncertainties are determined from the uncertainties in Gaussian fitting of H α . The uncertainties in F356W magnitude are based on measurements. However, we apply an uncertainty floor of ∆ m = 0 . 05 when calculating the continuum slope.  \nTable 2. The H α properties of LRDs in the ASPIRE fields. FWHM broad and FWHM narrow is the full-width-half-maximum measured for the broad and narrow H α components respectively. L H α, broad is the H α luminosity of the broad component and L H α, total is the H α luminosity of the total H α emission without absorption. EW H α, broad is the rest-frame equivalent widths of the broad H α emission. The narrow H α component of J0218P0007-BHAE-1 is unresolved and thus we do not show its FWHM narrow . J0430M1445-BHAE-1 and J0923P0402-BHAE-1 exhibit two Gaussian components both with FWHMs > 1000 km s -1 , so we measure the FWHM and luminosity of the entire H α emission lines as the FWHM broad and L broad .  \n<!-- image -->  \nFigure 1. The selection criteria of ASPIRE BHAE (red stars) among all photometric sources in the 25 ASPIRE fields (gray dots). We label the applied color cut and size cut as red lines.  \n<!-- image -->  \nTable 3. The black hole and galaxy properties of LRDs in the ASPIRE fields. M UV is the absolute UV magnitude at the restframe 1500 ˚ A. β UV is the UV slope derived from the F115W and F200W magnitudes. β opt is the optical slope derived from the F200W and F356W magnitudes, where we subtract the H α emission flux in the F356W to obtain the continuum flux. The black hole mass M BH and bolometric luminosity. The black hole mass and bolometric luminosity values have been calculated from the observed H α emission without dust correction, hence their intrinsic values might be higher.  \nFigure 2. The broad-H α emitters identified in 25 ASPIRE fields. For each source, the upper left panel shows a 2 '' × 2 '' RGB thumbnail composed of JWST /NIRCam F356W, F200W, and F115W images. The top panel shows the 2D grism spectrum with continuum and background removed. The bottom panel is the optimally extracted 1D spectrum (black) and the corresponding error spectrum (gray-shaded region). We show the total best-fit line profile as the red solid line, the narrow component as the orange dashed line, and the broad component as the green dashed line. For J0923P0402-BHAE-1, J1526M2050-BHAE-1, and J1526M2050-BHAE-3, we also show the H α profiles without the absorption features as the red dashed lines. J0923P0402BHAE-1 and J0430M1445-BHAE-1 present two broad H α components with FWHM > 1000 km s -1 . We label the FWHMs of their total H α emission, which are used to estimate the BH masses.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 3. Same as Figure 2. For J2232P2930-BHAE-1, we also show the [N II ] λλ 6548 , 6583 emission lines in blue and [S II ] λλ 6718 , 6732 in cyan. In the 2D spectrum of J2232P2930-BHAE-1, the emission on its top right originates from a companion HAE at the same redshift (companion-2, see § 4.2.2).  \n<!-- image -->  \nFigure 4 presents the properties of broad H α emission of ASPIRE BHAEs, together with literature highredshift BHAEs (Greene et al. 2024; Matthee et al. 2023; Harikane et al. 2023; Maiolino et al. 2023) and lowredshift type-1 quasars observed by SDSS (Wu & Shen 2022). We only present SDSS quasars at z < 0 . 6 whose H α lines are observable in SDSS spectra. As shown in the left panel of Figure 4, the FWHMs and luminosities of the broad H α emission of BHAEs are distributed similarly to those of type-1 quasars. It implies that BHAEs likely have broad-line regions similar to those of typical low-redshift type-1 quasars. On the other hand, the fraction of H α flux in the broad components, with a median of 0 . 70, as shown in the middle panel of Figure 4, is generally lower than that  \n<!-- image -->  \nl  \na  \nt  \no  \nt  \n,  \nH  \nL  \n/  \nd  \na  \no  \nr  \nb  \n,  \nH  \nL  \n<!-- image -->  \n<!-- image -->  \nL  \n/  \nL  \nFigure 4. The properties of broad H α emission lines of BHAEs. In the left panel, we show the relation between the FWHM and the luminosity of the broad H α components. We show the fraction of H α luminosity in the broad components as a function of the broad H α FWHMs in the middle panel, and as a function of the broad H α luminosity in the right panel. The ASPIRE BHAEs are denoted as the red stars. We present the measurements of literature BHAEs in Greene et al. (2024) (4 . 5 < z < 8 . 5, orange circle), Harikane et al. (2023) (4 < z < 7, orange square), Maiolino et al. (2023) (4 < z < 7, orange diamond) and Matthee et al. (2023) (4 . 2 < z < 5 . 5, orange triangle). The blue-shaded region denotes the distribution of z < 0 . 6 quasars with H α emission lines in SDSS spectra (Wu & Shen 2022).  \net al. 2011). We measure the FWHMs and luminosity of their entire H α emission for BH mass estimates. J0923P0402-BHAE-1, J1526M2050-BHAE-1, and J1526M2050-BHAE-3 show significantly blueshifted narrow absorption features of the H α emission peaks (see the third to last and second to last panel of Figure 2, and the first panel of 3). Similar features were also found in Matthee et al. (2023) and are interpreted as narrow H α absorption. Whether the features originate from H α absorption or just complexes of multi-component emission lines remains uncertain. In this work, we assume they are caused by H α absorption and fit them as a single Gaussian. We adopt the total H α luminosity values without absorption as the intrinsically emitted H α luminosity for L H α, total . We note that the broad H α fit is robust under the inclusion of the absorption, and exclude this component when measuring the broad H α components.  \nin most of low-redshift quasars (median 0.93). It indicates that the contributions from the narrow-line regions (NLRs) or host galaxies in these BHAEs to the total H α luminosity (median L H α, narrow /L H α, total = 0 . 30) are higher than those of low-redshift type-1 quasars (median L H α, narrow /L H α, total = 0 . 07). If the narrow H α components largely originate from the host galaxies, the high values of L H α, narrow /L H α, total suggest a higher star formation rate in the host galaxies of these BHAEs than those in typical low-redshift quasars. The right panel of Figure 4 shows that L H α, broad /L H α, total generally increases with L H α, broad . It is likely a natural outcome of dust obscuration within BLRs, where BLRs of BHAEs with lower L H α, broad are subject to greater attenuation.", '3.2. Black Hole mass and bolometric luminosity': 'We estimate the bolometric luminosity ( L bol ) and black hole mass ( M BH ) based on the broad H α emission line as done in recent works (e.g., Matthee et al. 2023; Harikane et al. 2023). Following Reines et al. (2013), we estimate M BH as:  \nlog 10 ( M BH /M ⊙ ) = 6 . 57 + log 10 ( ϵ ) +0 . 47 log 10 ( L H α, broad / 10 42 erg s -1 ) +2 . 06 log 10 ( FWHM broad / 10 3 km s -1 ) . (1)  \nThe geometric correction factor related to the broad line region, ϵ , is assumed to be 1.075 (Reines & Volonteri 2015). The BH estimator has an intrinsic uncertainty of 0.5 dex (Reines & Volonteri 2015). Note that, limited by our wavelength coverage, we derive the estimated M BH values here without dust attenuation correction, as done in Matthee et al. (2023); the M BH values could be even higher after correcting for the dust attenuation.  \nFigure 5. Black hole mass ( M BH ) versus bolometric luminosity ( L bol ). The red stars show the estimates in this work. The lengths of arrows indicate the changes of M BH and L bol when A V varies from 0 to 2. The orange circles, diamonds, triangles, and squares show the estimates from the literature BHAEs in Kocevski et al. (2023); Greene et al. (2024); Maiolino et al. (2023); Matthee et al. (2023); Harikane et al. (2023). Note that the L bol estimator used in Harikane et al. (2023) is different as discussed in § 3.2. We present the distribution of low-redshift ( z < 0 . 6) quasars from Wu & Shen (2022) as the blue-shaded region, with the depth demonstrating their number density. The blue-edged circles are measurements for z > 6 . 5 quasars, deriving M BH from H β emission lines and L bol from L 5100 . We draw equal Eddington ratio lines (gray dashed) of λ Edd = 0 . 01 , 0 . 1 , 1 , 10.  \n<!-- image -->  \nWe derive L bol by adopting the relation between the rest-frame 5100 ˚ A luminosity ( L 5100 ) and the AGNinduced H α luminosity in Greene & Ho (2005), with the bolometric correction in Richards et al. (2006):  \nL 5100 = 10 44 ( L H α / 5 . 25 × 10 42 erg s -1 ) 1 / 1 . 157 erg s -1 L bol = 10 . 33 × L 5100 . (2)  \nwhere L H α in Equation 2 should, include H α from both AGN narrow line and broad line regions in principle. Since it is still not clear whether the narrow H α components of BHAEs originate from AGN narrow line regions or the host galaxies, we use L H α, broad as the input. The L bol could be higher if the dust attenuation correction is applied or the contribution from AGN narrow-line regions is included. We discuss the caveat of the L bol estimator in § 4.1. The Eddington ratio λ Edd is defined as L bol /L Edd , where L Edd is the theoretical maximum luminosity achievable when radiation pressure and gravity are balanced in a spherical geometry.. We calculate  \nλ Edd following Trakhtenbrot et al. (2011):  \nλ Edd = L bol / erg s -1 1 . 5 × 10 38 M BH /M ⊙ . (3)  \nWe list the derived physical properties of BHAEs in Table 3. Figure 5 shows the estimated M BH and L bol of ASPIRE BHAEs together with literature BHAEs 3 and SDSS quasars. To display the effect of dust attenuation, we also show arrows associated with each BHAE. The lengths of the arrows indicate changes in M BH and L bol with A V varying from 0 to 2. The M BH of ASPIRE BHAEs ranges from 10 7 to 10 8 M ⊙ , comparable to those of low-redshift quasars. The H α converted L bol of BHAEs implies that the BHs are accreting at λ Edd ≈ 0 . 07 -0 . 47, with a median λ Edd of 0.17. The λ Edd of BHAEs are comparable to or slightly lower than those of z > 5 quasars (median values ∼ 0 . 3 -0 . 8, Shen et al. 2019; Farina et al. 2022; Mazzucchelli et al. 2023).', '3.3. UV and optical slopes': "The large difference in UV and optical slopes ( β UV , β opt , defined as f λ = λ β ) is one of the primary puzzles about the nature of so-called 'little red dots' (e.g., Matthee et al. 2023; Killi et al. 2023). We approximate β UV for all ASPIRE BHAEs by fitting the F115W and F200W photometry with a power law, which spans over rest-frame 2000 ˚ A - 3600 ˚ A at z ∼ 4 . 5. We assume that there are no strong emission lines in the two filters. The values of β opt are estimated using F200W and emission-line subtracted F356W photometry. We estimated the uncertainties of β UV and β opt through Monte Carlo (MC) method. This involved generating mock photometry based on the measurements and uncertainties and then recalculating the beta values 10 4 times. An uncertainty floor of 0.05 mag is applied. The derived slopes are listed in Table 3.  \nFigure 6 displays the relation between β opt and β UV , the fraction of H α flux in the broad component L H α, broad /L H α, total and M BH . Most of the BHAEs show redder β opt than β UV except J2232P2930-BHAE1. J2232P2930-BHAE-1 is the most extended source in our sample with bright and diffuse H α emission (see § 4.2.2 for more details). We do not see a correlation between β opt and L H α, broad /L H α, total in ASPIRE BHAEs. The Kendall's τ correlation analysis  \nFigure 7 shows the stacked images of BHAEs and the corresponding radial profiles. With the longest exposure time in F200W, both the F200W stack of original images and PSF-subtracted images present extended morphology. The radial profiles of the two F200W stacked images both deviate significantly from the shape of PSFs, with σ ≈ 3 . 5 at 0.1 arcsec and σ ≈ 1 at 0 . 4 arcsec. The  \n<!-- image -->  \n<!-- image -->  \nH ,total  \nFigure 6. The relation between the optical slope ( β opt ) and UV slope ( β UV ), the fraction of H α flux in the broad components ( L H α, broad /L H α, total ) and black hole mass ( M BH ). The measurements for ASPIRE BHAEs are denoted as red stars. The literature BHAEs in Greene et al. (2024) and Matthee et al. (2023) are shown as the orange circles and triangles. We label the source J2232P2930-BHAE-1 using a dashed square. More details about it are presented in § 4.2.2.  \n<!-- image -->  \non ASPIRE BHAEs yields a correlation coefficient of r = 0 . 27 ± 0 . 22 and p -value of 0 . 2 +0 . 4 -0 . 2 . Likewise, no strong correlation can be concluded from the β opt -M BH relation in ASPIRE BHAEs. The Kendall's τ correlation analysis between β opt and log M BH of ASPIRE BHAEs shows r = 0 . 19 +0 . 18 -0 . 17 and p -value of 0 . 3 +0 . 4 -0 . 3 . Our results are not consistent with the tight correlation seen in Matthee et al. (2023), who interprets the correlations as transitions from star formation-dominated phase with small M BH to dusty AGN-dominated phase with massive BHs. Our results suggest that the evolutionary track of these high-redshift AGNs might be non-linear or not identical among individual objects. However, we caution that the β UV and β opt values in both this work and Matthee et al. (2023) are calculated based on only two broad bands. The UV and optical continuum shapes may not be simple power laws, so the estimates with two broad bands can be highly uncertain. Future JWST /NIRSpec Prism and Grating spectra of larger BHAE samples are required to test the evolutionary scenario depicted in Matthee et al. (2023).", '3.4. Morphology': "Limited by the depth of ASPIRE observations, it is difficult to identify the host or companion galaxies of individual BHAEs. To detect possible extended emission, we stack the images of 14 BHAEs in the ASPIRE sample. Among the total 16 ASPIRE BHAEs, J2232P2930-BHAE-1 is not included, because it is a composite system presenting significant extended H α emission ( § 4.2.2). J2002M3013-BHAE-1 is not included, because its image is severely contaminated by bright spikes of saturated stars nearby. We generate 3 '' × 3 '' cutouts of the three bands for each broad H α emitter and normalize the images to the value of the central brightest pixels. We stack the normalized cutouts and  \nobtain the 3 σ clipped mean. We use the median value of the outskirt 10 pixels in the stack of each band as the background level and subtract it from the final stack. For cross-checking, we also perform PSF subtraction for each AGN using Galfit (Peng et al. 2002). The PSF models in each band are built by stacking gaia stars in the 23 ASPIRE fields (see more details in Champagne et al. in prep). We subtract the PSF models within the central 8-pixel × 8-pixel regions of the cutouts. The PSF-subtracted images are stacked in the same way, normalized using the central brightest pixels of the original cutouts. We calculate the radial profiles of the two stack versions. The calculation of PSF-subtracted radial profiles starts at r = 1 pixel (0 . 031 arcsec), to avoid the 2-pixel × 2-pixel over-subtracted region in the centers. To compare the radial profiles of stacked images with the PSF profiles, we normalize all the profiles to their maximum values, so that the maximum values are unity at the minimum radius. To quantify the deviation of stacked images from the PSF shape, we define the deviation as  \nσ ( r ) = S ( r ) -S PSF ( r ) δS ( r ) (4)  \nwhere σ ( r ) is the extent of the S ( r ) profile deviation from the PSF profile S PSF ( r ) with respect to the uncertainty δS ( r ). For instance, σ = 1 means that the radial profile S ( r ) deviates from the PSF shape at a level of one sigma.  \nstacked F200W extends to ∼ 2 . 6 physical kpc (pkpc) at z = 4 . 5. The F115W and F356W stacks are also tentatively extended with σ ≳ 1 within r < 0 . 3 -0 . 4 arcsec. The presence of extended structure in F200W, corresponding to ∼ restframe 3600 ˚ A at z ∼ 4 . 5, is likely caused by the stellar emission from the host galaxies.", '3.5. Luminosity function': "The large survey area of ASPIRE program on 275 arcmin 2 (25 ASPIRE fields) allows us to measure the number density of broad H α line-selected BHAEs at the bright end and reduce the effect of cosmic variance. We do not calculate the UV luminosity function (LF) of BHAEs since it is still uncertain whether the UV components of BHAEs originate from the host galaxies or the scattered light from AGNs. Instead, we focus on the H α LF using the total H α flux (narrow + broad) directly measured from the observations.  \nTo calculate the total survey volume, we first compute the effective survey area for each BHAE, where the noise levels enable us to reliably identify the broad wing of the H α emission lines. To do this, we need to ascertain the maximum root mean square (RMS) value for the identification of the broad H α lines. We generate a series of mock 2D grism spectra. We convolve the twocomponent Gaussian model of each ASPIRE AGN with the line spread function of NIRCam WFSS along the wavelength axis and the F356W PSF along the spatial axis. We then add noise to the 2D mock spectrum and extract the 1D spectrum using the boxcar algorithm. We fit the extracted 1D spectrum using a single Gaussian model and a two-component Gaussian model. If the two-component Gaussian model yields a S/N > 3 and | χ 2 red -1 | value that is 0.1 smaller than that of the single Gaussian model, the broad H α line is successfully identified. We progressively increase the noise levels added to the mock 2D spectra until the broad H α components cannot be characterized. The noise values at this stage are the maximum RMS for the broad H α identification. The survey volume is computed following the procedure in Sun et al. (2023). We define a series of redshifts spanning from 3.6 to 5.1, with H α emission lines residing at 3-4 µ m. We compute the effective sky area based on the spectral tracing and grism dispersion models for a specific redshift and construct the RMS maps using continuum-removed WFSS cal files. The maximum sky area at this redshift is determined as the area on the RMS map with RMS noise smaller than the maximum RMS for H α broad line identification. The maximum survey volume ( V max ) for each BHAE is integrated from the maximum sky area across z = 3 . 6 -5 . 1  \nover all the 25 ASPIRE fields. The median V max for our AGN sample is 911561 cMpc 3 .  \nThe H α LF of BHAEs is then computed using the direct 1/ V max method (Schmidt 1968):  \nΦ( L ) = 1 d log L ∑ i 1 C i V max ,i , (5)  \nwhere L is the H α luminosity of each LF bin, C i accounts for the completeness correction, and V max ,i is the maximum survey volume for the i -th source. To measure completeness, we randomly draw noise values from the effective RMS map as described above. We generate 2D mock grism spectra for each BHAE based on the two Gaussian models, adding the selected RMS values. We then re-extract the 1D spectrum following the procedure in Section 2. We run 1000 mock realizations for each BHAE, re-fitting the two-component Gaussian model to the extracted spectra. The completeness is defined as the fraction of the realizations in which the broad H α components are successfully characterized according to the criteria mentioned above. Additionally, we perform 1000 MC experiments for each luminosity bin to calculate the uncertainties, propagating the uncertainties of line fluxes to the LF. In each MC run, we account for the Poisson noise associated with small number statistics following the prescription in Gehrels (1986). We also take into account the potential incompleteness caused by the color cut (see Section 2.2). We generate 10 4 mock BHAE optical spectra assuming uniformly distributed L broad in log space with log( L broad / erg s -1 from 41.5 to 44, L broad /L tot from 0.1 to 1, β opt frin -2.5 to 2.5, FWHM broad from 1000 to 3000 km s -1 , FWHM narrow from 500 to 900 km s -1 , F200W magnitude from 25.5 to 27 mag, and redshift within z = 4 -5. We measure the F200W-F356W color excess for these mock BHAEs. We find that, for the LF of the total H α emission, the color cut only excludes 8.6% of BHAEs in the faintest H α luminosity bin ( L tot ≈ 10 42 . 5 erg s -1 ), and is nearly complete in more luminous bins. For the broad H α LF, the fraction is 12.6% in the faintest bin ( L broad ≈ 10 42 . 3 erg s -1 ) and also complete in more luminous bins. This minor incompleteness is insignificant when considering the substantial uncertainty and the completeness correction due to the background RMS. We do not correct the incompleteness caused by the color cut because the intrinsic distribution of the physical parameters remains unclear.  \nFigure 8 shows the H α LF for ASPIRE BHAEs and the number density of each H α luminosity bin is listed in Table 4. The number density of BHAEs are much lower than star-forming galaxies at similar redshifts. The fraction of BHAEs at each luminosity bin is computed as the  \nNormalized  \nFigure 7. Left: stack of images of 13 BHAEs in NIRCam F115W, F200W and F356W, normalized to the value of central pixels. The top panel shows the 3 '' × 3 '' stack of each band. In the middle panel, the radial profiles of the AGN stacks are shown as orange squares, and the radial profiles of PSFs are shown as red dashed lines. The bottom panel shows the deviations of the stacked AGN radial profiles from the corresponding PSFs as defined in Equation 4. Right: the stack of PSF-subtracted images of 13 BHAEs, also normalized by the value of central pixels in the original images. The radial profiles of the AGN stacks and the deviation from the PSFs are shown in the middle and bottom panels. The calculation of PSF-subtracted radial profiles starts at r = 1 pixel (0 . 031 arcsec), to avoid the 2-pixel × 2-pixel over-subtracted region in the centers. The r < 1 pixel ranges are grey-shaded in the middle and bottom panels.  \n<!-- image -->  \nnumber density of BHAEs over the number density of Lyman break galaxies as derived in Bollo et al. (2023). The fractions increase strongly with L H α , from 2% at log( L H α / erg s -1 ) ≈ 42 . 5 to 17% at log( L H α / erg s -1 ) ≈ 43 . 3. The most luminous H α emitters include a large portion of broad-line AGNs at z ≈ 4 -5. The rising trend of AGN fraction towards the higher luminosity end is in agreement with those observed at z ≈ 0 . 8 -2 . 2 (Stott et al. 2013; Sobral et al. 2016). To more precisely quantify the AGN fraction at high redshift, deeper spectroscopic observations are required to resolve the broad components in fainter H α emitters.  \nTo compare the number density of BHAEs with quasars and AGNs at similar redshifts, we estimate the bolometric LF of BHAEs by converting their L H α, broad into L bol following Equation 2. The bolometric LF is shown in Figure 9. For comparison, we also convert the broad H α LF of Matthee et al. (2023) into L bol . At the most luminous bin ( L bol ≈ 10 45 . 5 erg s -1 ), our result is consistent with the extrapolation of bolometric LF of quasars at z ∼ 4 but are ∼ 2 × higher than those of quasars at z ∼ 5, despite the large uncertainty. At the faintest end ( L bol ≈ 10 44 . 7 erg s -1 ) the number density of BHAEs is approximately 1 . 6 × that of z ∼ 4 quasars and 3 . 7 × that of z ∼ 5 quasars. The underlying BHAE population could have a higher number density with fainter L H α, broad . Our measured number density of  \nBHAEs is generally consistent with Greene et al. (2024) and Matthee et al. (2023), considering the uncertainties induced by dust attenuation, different selection criteria, and the cosmic variance.  \nHowever, it is unclear whether it is reasonable to compare the bolometric LF of quasars with that of BHAEs. The conversion between H α and bolometric luminosity in Equation 2 was calibrated based on the mid-infrared and optical properties of type-1 quasars at low redshifts. These high-redshift BHAEs may exhibit very different SEDs from those of typical quasars (e.g., P'erez-Gonz'alez et al. 2024). The estimated L bol are subject to significant systematic uncertainties. The different bolometric correction recipes also introduce uncertainties (Runnoe et al. 2012; Duras et al. 2020). Multiwavelength observations and comprehensive analysis of BHAE SEDs are essential to investigate their bolometric properties further.", '4.1. High number density of BHAEs': "Figure 9 suggests a high AGN number density at the lower end of the bolometric luminosity, albeit with significant potential systematic. We note that we do not apply any reddening correction in Figure 9. If we correct L bol using a typical value of A V = 2, L bol would increase by a factor of 2.3. Applying A V = 4 would  \nFigure 8. The H α (broad + narrow) luminosity function (LF) of the BHAEs (red stars). The H α LF of Lyman Break galaxies at z ∼ 4 . 5 is presented as the blue dasheddotted line and the prediction of JAGUAR mocks (Williams et al. 2018) for star-forming galaxies is presented as the blue dashed line. We show the fraction of BHAEs in each luminosity bin. The AGN fractions reported in Matthee et al. (2023), Harikane et al. (2023), and Maiolino et al. (2023) are labeled as cyan, blue, and purple-shaded regions. For comparison, we show the AGN fraction from the HiZELs sample at z ∼ 0 . 8 -2 . 23: the AGN fraction based on BPT diagram (Stott et al. 2013) as the green square, the AGN fraction identified using Spitzer/IRAC colors (Sobral et al. 2016) as the green diamonds and the AGN fraction by broad H α components (Sobral et al. 2016) as the green circles.  \n<!-- image -->  \nTable 4. The H α LF, BHAE fraction, and bolometric LF as shown in Figure 8 and 9. The number densities and fractions are all lower limits since we do not correct for the completeness.  \nincrease L bol by 6.9. Consequently, the bolometric LF would shift to the bright end by 0.4-0.8 dex, increasing the excess of BHAE number density over that of type-1 AGNs. As discussed in Greene et al. (2024), while the  \nFigure 9. The bolometric luminosity function of BHAEs (red stars). The bolometric luminosities of ASPIRE BHAEs are estimated based on the broad H α luminosities. We also show the broad H α -converted bolometric LF of z ≈ 4 . 5 -6 . 5 AGNs in Greene et al. (2024) and the L 5100 -converted bolometric LF in Kokorev et al. (2024) as the blue diamonds and orange squares, and convert the broad H α LF of Matthee et al. (2023), which were not corrected for dust attenuation either, into the bolometric LF. The bolometric LFs of quasars at z ∼ 4 and z ∼ 5 are shown as the gray dashed and dashdotted lines (Shen et al. 2020). We also present the predicted AGN bolometric LFs of Delphi at z = 5 (Dayal et al. 2019), CAT at z = 4 (Trinca et al. 2022), Ricarte & Natarajan (2018) at z = 4, and the range of AGN bolometric LFs at z = 4 from different hydrodynamical cosmological simulations concluded in Habouzit et al. (2022).  \n<!-- image -->  \nnumber density of AGNs is larger than what would be expected from the extrapolation of the quasar luminosity function towards faint luminosities, a large number density of faint quasars is generally consistent with predictions from theoretical models, at least at z < 6 . 5, while there may be some tension at the bright end for z > 6 . 5. Indeed, the overprediction of the faint end of the luminosity function appeared to be common, as shown for instance in the compilation of results from large-volume simulations of Habouzit et al. (2022). The updated AGN number densities is in much better agreement with the predicted large number of high-redshift faint AGNs. Most of these large-volume simulations, however, include simplistic predictions for black hole formation and growth, since they cannot resolve the subparsec scales relevant for black hole physics.  \nSemi-analytical models tend to include more specific recipes for black hole formation (e.g., Ricarte & Natara-  \njan 2018; Dayal et al. 2019; Trinca et al. 2022). Many of these models are based on bimodal black hole seed models: 'light' seeds with mass < 10 3 M ⊙ remnants of the first stars, or 'heavy' seeds with mass > 10 5 M ⊙ formed from the collapse of supermassive stars formed under very specific conditions (see Inayoshi et al. 2022; Volonteri et al. 2023, for reviews). They usually find that the number density of black hole seeds decreases with the seed mass, but the ability to grow instead increases with the seed mass. Interestingly, recent simulations and theoretical models of black hole seed formation suggest a less extremely bimodal situation (e.g., Regan et al. 2020; Schleicher et al. 2022), with numerous seeds at mass ∼ 10 3 M ⊙ , more consistent with a widespread population of relatively bright AGNs presented in this paper.  \nFinally, as discussed in § 3.5, the observed high-redshift BHAE population is not yet fully characterized. The bolometric luminosity of BHAEs might not be well described by Equation 2. A valid comparison of BHAE number density between simulations and observations relies on thoroughly understanding BHAEs' nature, bolometric properties, selection function, and completeness through future observations.", '4.2. Evidence of clustering of BHAEs': 'In this section we present early evidence of strong clustering of BHAEs.', '4.2.1. A pair of BHAEs': "J1526M2050-BHAE-2 and J1526M2050-BHAE-3, both at z = 4 . 87, are 1.36 arcmin apart from each other, corresponding to a projected separation of 519 pkpc. The velocity offset between the two is ≈ 222 km s -1 . There might be physical connections between these two BHAEs, which could reside in the same overdensity or large-scale structure. A systematic search for galaxies at a similar redshift is needed for further detailed studies.  \n4.2.2. J2232P2930-BHAE-1: an extended H α emitter with close neighbors  \nAmong the 15 ASPIRE BHAEs presented in this paper, J2232P2930-BHAE-1 at z = 4 . 135 has the brightest F356W magnitude ( m F356W = 23 . 39) and very high H α luminosity ( L H α, tot = 2 × 10 43 erg s -1 ). Its F356W morphology is more extended than that in the F200W band, indicative of diffuse extended H α emission. We present its F356W image, PSF-homogenized F200W contour, and the corresponding radial profiles in the left and middle panels of Figure 10. Its F115W and F200W images reveal four marginally resolved components as shown  \nin the right panel of Figure 10. It also shows significant spatial variation in β UV . The upper component (( x, y ) = (0 , 0 . 15 '' )) has a red UV slope ( β UV ≳ 1) and the left component is comparatively blue ( β UV ≲ -0 . 5). It is not clear if there is β UV variation in the two central F200W contour peaks, limited by the spatial resolution of F200W. The variation of β UV might be caused by the inhomogeneous distribution of dust content. As a composite system with multiple components, the broad line profile of J2232P2930-BHAE-1 could be a result of blended H α emission from individual substructures, and not necessarily corresponding to an AGN broad line region. Future observations with JWST /NIRSpec IFU are required to investigate its internal structure.  \nMoreover, J2232P2930-BHAE-1 has three close companion galaxies at the same redshift. As shown in Figure 11, the three companion galaxies, named Companion-1, Companion-2, and Companion-3, are 2.1 arcsec, 1.6 arcsec, and 3.85 arcsec away from J2232P2930-BHAE-1, corresponding to 14.5 pkpc, 10.8 pkpc and 26.4 pkpc at z = 4 . 13. The Companion-1, Companion-2, and Companion-3 have H α luminosity of (1 . 7 ± 0 . 7) × 10 42 erg s -1 , (1 . 9 ± 0 . 2) × 10 42 erg s -1 and (5 . 2 ± 0 . 5) × 10 41 erg s -1 , respectively. It implies that J2232P2930-BHAE-1 may reside in an overdense region. There is a tentative detection of diffuse H α emission between Companion-2 and J2232P2930-BHAE-1, as shown in both the RGB image (left panel of Figure 11) and the 2D grism spectrum (middle panel of Figure 11), indicating possible interaction between the two sources. Detailed studies on the environment of J2232P2930-BHAE-1 will be performed in future works with a systematic search for H α emitters at the same redshift.", '4.3. Blue-shifted H α absorption': 'J0923P0402-BHAE-1, J1526M2050-BHAE-1, and J1526M2050-BHAE-3 present blueshifted absorption features on top of their H α emission lines. The velocity offsets between the absorption and the H α emission peak for J0923P0402-BHAE-1, J1526M2050-BHAE-1, and J1526M2050-BHAE-3 are 586 km s -1 , 357 km s -1 , and 266 km s -1 , respectively.  \nIf these features are H α absorption of high-columndensity gas, the fraction is high (19% in ASPIRE BHAEs, 10% in the sample of Matthee et al. 2023, 15% in the sample of Kocevski et al. 2024) compared to lowredshift type-1 AGNs ( < 0 . 1%, e.g., Aoki et al. 2006; Zhang et al. 2018), as shown in Figure 12. The prevalence of H α absorption, if confirmed, suggests that these high-redshift BHAEs could differ significantly from lowredshift type-1 AGNs. High-redshift BHAEs may contain abundant outflowing dense gas, which could be re-  \nFigure 11. J2232P2930-BHAE-1 and its three neighboring galaxies. The left panel shows the RGB image composed of F115W, F200W, and F356W. We label the three companion galaxies using white plus markers. In the middle panel, we present the 2D WFSS spectra of the three galaxies, and the extracted 1D spectra are shown in the right panel. The 2D and 1D spectra are both displayed in the velocity frame relative to J2232P2930-BHAE-1. We can see tentative diffuse H α gas between Companion-2 and J2232P2930-BHAE-1 in the RGB image and the 2D grism, suggesting possible interaction between the two sources.  \n<!-- image -->  \n<!-- image -->  \nFigure 10. Left: JWST /NIRCam F356W image of J2232P2930-BHAE-1. The blue lines denote the 2 σ and 7 σ contours of the F356W image. The green lines are the 2 σ and 7 σ contours of the PSF-homogenized F200W image, matched to the F356W PSF. Middle: The surface brightness radial profile of J2232P290-BHAE-1. The blue squares show the F356W surface brightness profile of J2232P2930-BHAE-1 and the green circles show the profile of PSF-homogenized F200W. We show the radial profile of the F356W PSF as the black dashed line. Right: The UV slope map of J2232P2930-BHAE-1, produced using F115W and F200W images, PSF-matched to F200W. We smooth the map using a 1-pixel Gaussian kernel. The magenta and cyan lines denote the 2 σ , 5 σ , and 10 σ contours of the F115W and F200W images, respectively.  \n<!-- image -->  \n<!-- image -->  \nCompanion-1  \nCompanion-2  \nCompanion-3  \n-2000  \n0  \n2000  \nv  \n(km s  \n1  \n)  \n2000  \n0  \n2000  \nv  \n(km s  \n1  \n)  \nlated to either BH activity or strong galactic feedback in the host galaxies.', '4.4. Comparison with literature LRDs and selection effect': "We compare photometrically and broad-line selected samples in the literature in Figure 13. The broad-line selected samples (e.g., ASPIRE, Matthee et al. 2023, and Harikane et al. 2023), which imposes no constraints on the SED shapes, contain BHAEs with β opt ≲ 0, as illustrated in the left panel of Figure 13. On the other hand, the samples photometrically selected based on V-shaped SEDs (Greene et al. 2024 and Kocevski et al. 2024) do not include such objects. Meanwhile, the broad-line se-  \nted BHAEs lack UV-luminous but red objects (with M UV ≲ -19 . 5 and β opt > 0).  \nThe photometric selection requiring V-shaped SEDs is biased towards sources with reddened rest-frame optical spectra and against those with relatively blue optical continua ( β opt < 0). The flux-limited broad-line selection fails to detect objects with weak broad H α flux, as their H α broad wings could be overwhelmed by background noise. As shown in the middle panel of Figure 13, the UV-luminous and red sources in Greene et al. (2024) exhibit broad H α luminosity lower than those of most JWST /WFSS-selected BHAEs. Their broad H α components may not be reliably identified with the shallow grism data. We note that among the NIRSpec-  \nz  \n=4.135  \nS/N=2.3  \nz  \n=4.138  \nS/N=12.6  \nz  \n=4.136  \nS/N=9.2  \n1.0  \n0.5  \n1.0  \n0.5  \n1.0  \n0.5  \n)  \n1  \nÅ  \n2  \nm  \nc  \n1  \ns  \ng  \nr  \ne  \n8  \n1  \n0  \n1  \n(  \nf  \n100%  \nRedshift  \nFigure 12. The fraction of AGNs with Balmer absorption across the cosmic time. The errorbars in the redshift axis indicate the redshift range of the AGN sample, while the uncertainties in the fraction on the y -axis are calculated assuming a Poisson distribution.  \n<!-- image -->  \nselected sample from Harikane et al. (2023), the two sources with reddened rest-frame optical continua also have the faintest M UV ( M UV > -18). Since their broadline selection relies on archive data, the prior source selection for the NIRSpec MSA is crucial but can be intricate.  \nThe distribution of grism-selected BHAEs on the β opt -M UV and M UV -EW H α, broad planes (the left and right panels of Figure 13) is attributed to the fluxlimited selection effect. The absence of UV-faint and optical-blue objects ( M UV ≳ -18 and β opt < 0) can be attributed to their weak broad H α wings falling below the detection limit, as they might have lower-mass BHs and host galaxies following the co-evolution of the BH-galaxy. Although the UV-luminous and optically red BHAEs ( M UV > -19 and β opt > 0) missing from grism-selected sample all have low L H α, broad , similar objects with high L H α, broad might exist but are rare due to the rapid decline towards the bright end of the LFs. On the other hand, it remains unclear whether BHAEs with β opt > 0 and β opt < 0 belong to the same population. A plausible interpretation is that β opt < 0 BHAEs may closely resemble type-1 quasars. It is also possible that BHAEs with β opt < 0 are counterparts of those with β opt > 0, but with less obscured AGN components. In this case, if the UV emission originates from the host galaxies, the stellar light from UV-luminous star-forming galaxies may contribute significantly to the  \nrest-frame optical spectra, potentially leading to a bluer β opt .  \nWe notice that the two optically red BHAEs identified in Harikane et al. (2023) are both compact and pointsource-like, in contrast to the other objects they selected, which are either extended or have bright companions. The depth of ASPIRE images does not allow us to resolve the morphology of individual BHAEs. Instead, we examine the BHAEs in Matthee et al. (2023). Among the nine BHAEs with β opt < 0 in Matthee et al. (2023), seven have extended F200W (or F182M + F210M) morphologies, and only two are compact and point-sourcelike. Among the 11 β opt > 0 BHAEs, four show significant extended F200W morphologies and the remaining seven are compact. The 'compactness rate' is 63% for β opt > 0 BHAEs compared to 22% for β opt < 0 BHAEs. This can also be due to a selection effect. As discussed above, the grism-selected BHAE sample misses UV-luminous but optically red objects, and as shown in the left panel of Figure 13, the grism-selected β opt > 0 BHAEs are all UV faint. Their UV faintness makes it challenging to resolve the extended components in shortwavelength images. To conclusively determine if there is any potential correlation between optical redness and morphologies, a larger and unbiased sample is essential.  \nIn summary, different selection criteria introduce various biases that must be carefully considered to understand the nature of this newly discovered BHAE population. A systematic classification is essential for understanding their properties. To explore a broader range of samples in the future, a more inclusive selection criterion that considers both AGN and galaxy contributions in the rest-frame optical band is necessary.", '5. SUMMARY': 'We report the discovery of 16 broad H α emitters (BHAEs) at z = 4 -5 selected from the JWST ASPIRE program. These BHAEs are primarily selected due to their compactness and redness. Their broad H α emission lines are identified utilizing the JWST /NIRCam WFSS of F356W. Our main results are summarized below.  \n- (1) The ASPIRE BHAEs exhibit broad H α emission components with FWHM (FWHM broad ) ranging from 1000 to 3000 km s -1 and luminosity ( L H α, broad ) from 10 42 to ≳ 10 43 erg s -1 . The fraction of H α flux in the broad components ( L H α, broad /L H α, total ) spans a wide range from ∼ 0 . 4 up to 1 with a median value of 0.7. These BHAEs are distributed similarly in the FWHM broad -L H α, broad diagram compared to lowredshift ( z < 0 . 6) SDSS quasars. On the other  \nFigure 13. Left panel: the β opt versus M UV distribution of JWST broad-line and photometrically selected LRDs. The broadline selected samples include the ASPIRE BHAEs (red) and the BHAEs in Matthee et al. (2023) (orange). The photometrically selected samples include spectroscopically confirmed BHAEs in Greene et al. (2024) (cyan) and LRDs in Kocevski et al. (2024) (gray). We show the density distribution of β opt for the Kocevski et al. 2024 sample. Middle panel: the broad H α luminosity versus M UV . The filled markers denote the BHAEs with reddened optical continua ( β opt > 0) while white ones have β opt < 0. Right panel: the rest-frame equivalent width of the broad H α versus M UV , color-coded by β opt .  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nUV  \nhand, BHAEs have lower L H α, broad /L H α, total than most of the SDSS quasars, suggesting larger contributions from the host galaxies or narrow line regions to the total H α luminosity.  \n- (2) The black hole masses ( M BH ) of ASPIRE BHAEs range from ∼ 10 7 to ∼ 10 8 M ⊙ , comparable to those of BHs in low-redshift SDSS quasars. The H α converted bolometric luminosity ( L bol ) suggests that these BHAEs are accreting at eddington ratios of 0.07-0.47, with a median value of 0.17.\n- (3) In general, the ASPIRE BHAEs have UV continuum slopes ( β UV ) bluer than the optical continuum slopes ( β opt ). In contrast to Matthee et al. (2023), we do not see significant correlations between β opt and L H α, broad /L H α, total or M BH .\n- (4) The stacked images of ASPIRE BHAE show extended components in F200W. The F200W stack extends to ∼ 0 . 4 arcsec, corresponding to 2.6 physical kpc at z ∼ 4 . 5. The extended morphology in the F115W and F356W are tentative, and limited by the shallow depth of observations. The extended structure suggests the presence of host galaxies.\n- (5) The H α LF of ASPIRE BHAEs suggests that the AGN fraction increases towards higher H α luminosities, consistent with low-redshift H α emitters. The BHAE fraction spans from 2% at total H α luminosity of log( L H α / erg s -1 ) ≈ 42 . 5 to 17% at log( L H α / erg s -1 ) ≈ 43 . 3.\n- (6) We find a pair of BHAEs at the same redshift ( z ≈ 4 . 87) with a separation of 519 physical kpc.  \nAnother isolated BHAEs, J2232P2930-BHAE-1, as a composite system with bright extended H α emission and significant β UV variation, resides in an overdense region with three close companion galaxies. These, for the first time, provide tentative evidence for the strong clustering of BHAEs.  \n- (7) We find three BHAEs with blueshifted H α absorption, potentially indicating the presence of outflowing dense gas. The prevalence of H α absorption (19% in ASPIRE BHAEs) might be a unique feature of this high-redshift AGN population.\n- (8) We compare the broad-line and photometrically selected BHAE samples from the literature. The photometrically selected sample, based on Vshaped SED, does not contain optically blue ( β opt < 0) objects, while the flux-limited grismselected sample lacks UV-luminous but optically red ( M UV ≲ -19 . 5 and β opt > 0) objects. The optically red BHAEs tend to be more compact in rest-frame UV bands. However, these differences can be explained by the selection effect. A more inclusive selection method is needed to study a broader range of BHAEs.  \nThe ASPIRE BHAE sample with spectroscopically identified broad H α emission provides a good database for further studies on low-luminosity AGNs in the early Universe. To further investigate the nature of highredshift low-luminosity AGNs and their roles on the evolutionary tracks of early BH assembly, future followups in both the restframe UV and infrared wavelengths are required. Deep spectroscopic and IFU observations covering the UV continua and emission lines can provide constraints on the origin of the extended blue  \ncomponents. Observations in the infrared band with JWST /MIRI are crucial to determine the full SED shapes and to understand the roles of AGNs and starbursts in these broad-line selected samples.', 'ACKNOWLEDGMENTS': "X.L. and Z.C. are supported by the National Key R&D Program of China (grant no. 2023YFA1605600), the National Science Foundation of China (grant no. 12073014), the science research grants from the China Manned Space Project with No. CMS-CSST2021-A05, and Tsinghua University Initiative Scientific Research Program (No. 20223080023). F.W. and X. F. acknowledge support from NSF Grant AST2308258. F.S. acknowledges JWST/NIRCam contract to the University of Arizona NAS5-02015. JBC acknowledges funding from the JWST Arizona/Steward Postdoc in Early Galaxies and Reionization (JASPER) Scholar contract at the University of Arizona. AL ac-  \nknowledges support by the PRIN MUR '2022935STW'. MT acknowledges support from the NWO grant 016.VIDI.189.162 ('ODIN'). JTS is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Project number 518006966.  \nThis 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. These observations are associated with program #2078. The specific observations analyzed can be accessed via 10.17909/vt74-kd84. Support for this program was given through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127.  \nFacility: JWST  \nSoftware: JWST calibration pipeline (v1.10.2 Bushouse et al. 2023) SciPy (Virtanen et al. 2020) SourceExtractor++ (Bertin et al. 2020)", 'REFERENCES': "Farina, E. P., Schindler, J.-T., Walter, F., et al. 2022, ApJ, 941, 106, doi: 10.3847/1538-4357/ac9626 Furtak, L. J., Labb'e, I., Zitrin, A., et al. 2023, arXiv e-prints, arXiv:2308.05735,  \n```\ndoi: 10.48550/arXiv.2308.05735 Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272 Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051 Gehrels, N. 1986, ApJ, 303, 336, doi: 10.1086/164079 Greene, J. E., & Ho, L. C. 2005, ApJ, 630, 122, doi: 10.1086/431897\n```  \nGreene, J. E., Labbe, I., Goulding, A. D., et al. 2024, ApJ, 964, 39, doi: 10.3847/1538-4357/ad1e5f  \nGreene, T. P., Chu, L., Egami, E., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9904, Space Telescopes and Instrumentation 2016: Optical, Infrared, and Millimeter Wave, ed. H. A. MacEwen, G. G. Fazio, M. Lystrup, N. Batalha, N. Siegler, & E. C. Tong, 99040E, doi: 10.1117/12.2231347  \nHabouzit, M., Somerville, R. S., Li, Y., et al. 2022, MNRAS, 509, 3015, doi: 10.1093/mnras/stab3147  \n```\nWilliams, C. C., Alberts, S., Ji, Z., et al. 2023, arXiv e-prints, arXiv:2311.07483, doi: 10.48550/arXiv.2311.07483 Wu, Q., & Shen, Y. 2022, ApJS, 263, 42, doi: 10.3847/1538-4365/ac9ead\n```  \n```\nYang, J., Wang, F., Fan, X., et al. 2021, ApJ, 923, 262, doi: 10.3847/1538-4357/ac2b32 Yang, J., Fan, X., Gupta, A., et al. 2023, ApJS, 269, 27, doi: 10.3847/1538-4365/acf99b Zhang, S., Zhou, H., Shi, X., et al. 2018, AJ, 156, 4, doi: 10.3847/1538-3881/aac433\n```"}