Tripling the Census of Dwarf AGN Candidates Using DESI Early Data

Introduction

Supermassive Black Holes (SMBHs; $M_{\rm BH} \gtrsim10^6\,M_\odot$) are found at the centers of all bulge-dominated galaxies with stellar masses $M_\star\gtrsim10^{10}\,M_\odot$ ([82, 95, 81]), yet their origin is far from understood. Theoretical studies suggest that they grow via accretion and mergers from initial “seed” black holes (BHs) formed early in the Universe ([150]). While these early BHs at high redshifts are difficult to detect, with current capabilities ($M_\star\lesssim3 \times10^9\,M_\odot$) galaxies in the nearby universe can host lower mass central BHs that can help constrain the BH seed formation models ([150, 52]).

The search for central BHs in dwarf galaxies, especially for the still elusive Intermediate-mass Black Holes (IMBHs; $M_{\rm BH} \le10^6\,M_\odot$), has been actively pursued over the past decade (see Mezcua [104]; Greene et al. [56]; Reines [128], for reviews). The accretion of matter onto these BHs powers active galactic nuclei (AGN), releasing a large amount of energy across the electromagnetic spectrum. In the optical regime, emission lines arising from ionized gas can be used to search for signatures of AGN photoionization. For instance, the BPT (“Baldwin, Phillips & Terlevich”; [11, 149]) emission-line diagnostic diagrams have been widely employed to identify hundreds of dwarf AGN candidates ([130, 114, 107, 113, 123, 135, 141, 109]). Beyond optical diagnostics, other methods such as infrared color-color diagrams ([76, 93, 87]), radio observations ([111, 129, 34]), X-ray emission ([90, 120, 105, 106, 15, 86, 23, 134]), and variability techniques ([9, 21, 152]) have been utilized to increase the census of dwarf AGN candidates. Despite this progress, questions remain regarding the prevalence of AGN in dwarf galaxies.

Optical spectroscopy also offers an avenue to infer key physical properties of AGN, including their luminosity and BH masses. The gravitational influence of the BH accelerates the gas in its vicinity, resulting in broad emission lines in the spectrum. The broad H$\alpha$ emission from single-epoch spectroscopy is frequently used to estimate BH masses in galaxies through virial techniques [54, 130, 114, 10, 131, 29, 145]. Baldassare et al. [10] identified a BH with a mass of $\approx50{,}000\ M_{\odot}$ in a nearby dwarf galaxy, which is the lowest mass BH detected in a galaxy so far. However, further efforts to identify more of such lower-mass BHs are essential for discerning the BH seed formation models.

Scaling relations between mass of the central SMBH and various galaxy properties such as the mass/luminosity of the stellar bulge [47, 97, 88, 61, 100], the velocity dispersion of stars in the bulge [47, 50, 61, 100], and the galaxy stellar mass [131, 145], have been extensively studied. The tightness of the observed correlations suggests that SMBHs and their host galaxies co-evolve [74, 81]. However, whether this co-evolution extends to lower-mass galaxies is unknown.

The low-mass end of the scaling relations can provide insight into the BH seed formation mechanisms and the efficiency of BH growth in dwarf galaxies [56]. By measuring the stellar velocity dispersion ($\sigma_{\star}$) of eight broad-line candidates from Reines et al. [130], Baldassare et al. [9] extended the $M_{\mathrm{BH}}-\sigma_{\star}$ relation down to $\log(M_{\star}/M_{\odot}) \approx8.96$. They found that the dwarf galaxies are in good agreement with the extrapolation of the existing relation [81]. In contrast, using a sample of 127 low-mass Seyfert 1 galaxies, Martín-Navarro & Mezcua [98] found a flattening at the low-mass end of the $M_{\mathrm{BH}}-\sigma_{\star}$ relation. Measuring $\sigma_{\star}$ in low-mass systems is difficult to extend to a larger sample and to farther distances due to the faintness of targets and the limitations of current telescopes and instrumentation. Therefore, even though the $M_{\mathrm{BH}}-\sigma_{\star}$ relation usually shows the tightest correlation, the $M_{\mathrm{BH}}-M_{\star}$ relation is preferred for both extending the scaling relations to lower galaxy masses and/or for comparing to higher redshifts.

Using a sample of 262 broad-line AGN and 79 galaxies with dynamical BH masses, including fourteen dwarf galaxies, Reines & Volonteri [131] measured the $M_{\mathrm{BH}}-M_{\star}$ relation down to $\log(M_{\star}/M_{\odot}) \approx9.0$. They found two distinct scaling relations for early-type and late-type galaxies, with the former having a higher normalization than the latter and the dwarf galaxies falling on the tail of the late-type galaxy relation. More recently, Mezcua et al. [110] and Mezcua et al. [109] found nineteen low-mass galaxies at $z \approx0.9-3$ with over-massive BHs that reside $\gtrsim2$ dex above the $M_{\mathrm{BH}}-M_{\star}$ relation. Several such over-massive BHs were also discovered at higher redshifts, $z \approx4-11$ using observations from the James Webb Space Telescope (JWST; [65, 96, 147]). However, these high-redshift studies are subject to observational bias, as they predominantly detect luminous AGN. Overall, these studies are constrained by limited sample sizes at the low-mass end, and the scaling relation remains poorly constrained at both low stellar and BH masses, even at low redshifts.

In this paper, we use early spectroscopic data from the Dark Energy Spectroscopic Instrument (DESI; [35] (DESI Collaboration et al. 2022)) survey to identify AGN in galaxies spanning a wide range of stellar masses ($6 \leq\log(M_{\star}/M_{\odot}) < 12$), specifically to increase the census of AGN in dwarf galaxies, and to extend the $M_{\mathrm{BH}}-M_{\star}$ scaling relation down to lower galaxy masses. The early DESI data is the largest sample of optical spectra to date, encompassing nearly 1.5 million unique galaxies at $z \leq0.45$ for our analysis.

This paper is organized as follows. § 2 describes the spectroscopic and photometric data used in the paper. We present the estimation of stellar masses and emission-line measurements of the targets in § 3. § 4 and § 5 details the AGN identification and the overall $M_{\mathrm{BH}}-M_{\star}$ scaling relation, respectively. We discuss our results in § 6 and summarize our conclusions in § 7. Throughout the paper, we assume the Chabrier initial mass function (IMF; [27]), and the Planck Collaboration et al. [122] cosmology with $H_{0}=67.4\ \mathrm{km\ s^{-1}\ Mpc^{-1}}$ and $\Omega=0.315$. All wavelengths are presented in vacuum wavelengths, and all magnitudes are given in the AB system [119].

Data

Spectroscopy

DESI is a 5000-fiber multi-object spectrograph on the 4-meter Mayall telescope at Kitt Peak National Observatory (KPNO), covering a spectral range of 3600 Å to 9800 Å. The resolving power of the instrument ranges from 2000 at 3600 Å to 5500 at 9800 Å ([140, 112] (DESI Collaboration et al. 2022)). This instrument will be used to conduct a cosmological survey over an area of $\approx14{,}000\,\mathrm{deg}^2$ over five years and will obtain optical spectra of nearly 35 million galaxies and quasars (Levi et al. 2013; DESI Collaboration et al. 2016a,b; Schlafly et al. 2023). DESI completed its survey validation (SV) observations in May 2021, and is currently carrying out its main survey. The SV data, covering $\approx1{,}400\,\mathrm{deg}^2$ of the sky, is released as the DESI Early Data Release (EDR; [38 b]). We use the entire DESI EDR and 20% of Year 1 (DA0.2) data for analysis in this paper.

DESI targets are divided into five primary classes: 1) Milky Way Survey (MWS; [31]), 2) Bright Galaxy Survey (BGS; [133, 64, 73]), 3) Luminous Red Galaxies (LRG; [158, 159]), 4) Emission Line Galaxies (ELG; [125, 126]), and 5) Quasars (QSO; [156, 28]). These are further complemented by secondary targets via spare fibers (SCND; see Appendix of DESI Collaboration et al. [38]; [33, 46]). This target selection is validated based on a thorough visual inspection of the SV spectra ([4, 83]). Depending on the selection criteria, some sources can be targeted in multiple targeting classes. The DESI targeting selection is described in detail by Myers et al. [118]. The biases related to this target selection on our measurements are briefly discussed in § 3.1 and § 4.3.

The spectra obtained from DESI are reduced using the DESI spectroscopic pipeline ([63]), and the redshifts are obtained by the Redrock redshift-fitting pipeline ([19, 5], Bailey et al., in preparation). The DESI EDR and DA0.2 spectra are reduced and released as ‘fuji’ and ‘guadalupe’ internal data releases, respectively. Combining the redshift catalogs from both these releases results in 6,059,937 spectra of 5,859,563 unique targets. For objects with multiple spectra, we choose the “best spectrum” based on the ZCAT_PRIMARY flag. This prioritizes the spectrum that has no issues with fiber and redshift measurements. If multiple spectra pass these criteria, then we choose the highest signal-to-noise ratio spectrum as the “best spectrum” for the object (see Section 3.3.3 in DESI Collaboration et al. [38] for more information). After selecting the unique spectra of all objects, we apply the following data quality cuts:

COADD_FIBERSTATUS = 0¹, ZWARN = 0 or 4². This selects 4,323,517 sources that have been observed without any fiber issues and have good redshift measurements. Of these, we select sources that are assigned spectroscopic types of GALAXY or QSO by Redrock, resulting in a total of 3,483,455 sources. Among these, 1,628,394 sources have redshifts in the range of $0.001 \le z \le0.45$. The lower redshift limit is to exclude contamination from Galactic sources. Given that we focus on identifying AGN candidates using standard emission-line ratio diagnostics (§ 4.1), the upper redshift limit is determined by the availability of [SII]$\lambda$6717,6731 emission lines in the DESI spectra (see § 3.2 for more information).

For some QSO targets, the redshifts determined by Redrock were found to be incorrect, resulting in the presence of high-redshift contaminants in our sample ([4]). To accurately determine their redshifts, the DESI team used a machine learning algorithm called QuasarNet, which is a deep convolutional neural network classifier specifically designed to identify quasars and their redshifts. Of the 1,628,394 Redrock redshifts, we find that 4,286 are classified as high-redshift ($z \gtrsim0.5$ – 4) quasars by QuasarNet. We remove these from our selection, resulting in a spectroscopic sample of 1,624,108 sources in our redshift range of interest.

Photometry

The primary DESI targets are selected based on the DESI Legacy imaging surveys ([40]). The ninth data release for this survey (DR9) comprises a set of optical and infrared imaging data covering the entire 14,000 deg$^2$ of the DESI survey (Schlegel et al., in preparation). The photometric catalog consists of optical imaging data from three different surveys: The Beijing-Arizona Sky Survey (BASS; [160]), the DECam Legacy Survey (DECaLS), and the Mayall $z$-band Legacy Survey (MzLS). Furthermore, it also includes the WISE/NEOWISE (unWISE; [84, 102, 103]) coadded images from the Wide-field Infrared Space Explorer (WISE; [154]).

Imaging from the Legacy Surveys is processed using inference modeling with the Tractor³ code ([85]). This produces a catalog that contains photometry in $g$, $r$, and $z$ bands, based on an optimal model of the source morphology. It also includes mid-IR photometry measured from the forced photometry of the unWISE coadds using the optically derived model. This detects fainter sources than the traditional approach of WISE photometry [85]. Furthermore, the morphology (TYPE) and Sérsic index (SERSIC) of each source are also computed based on this model.

The optical photometry is corrected for Galactic extinction using dust emission maps from Schlegel et al. [138]. The extinction values for the WISE filters are derived from Fitzpatrick [48], using the recommendations by Schlafly & Finkbeiner [136]. These are reported in the LS DR9 catalog in linear units of Milky Way transmission.

In addition to the primary DESI targets, the spectral data from the survey also includes secondary targets as well as targets of opportunity, which might not always have the Legacy Survey photometry. Given that our analysis requires the stellar mass calculation using this photometry, we select the 1,623,425 objects that have DR9 photometry. Additionally, we apply two quality cuts that remove false detections, extremely faint sources, and a significant fraction of fragmented sources:

We do not apply any cuts on the WISE photometry to include faint dwarf galaxies, which might not always have good available WISE detections. Given the strong cuts on optical photometry, we will instead focus on quality cuts based on stellar mass estimates (§ 3.1). The above cuts result in 1,492,821 unique sources with good spectra and valid photometry as the starting sample for our analysis.

Methods

Stellar Masses

We fit the spectral energy distributions (SED) using Code Investigating GALaxy Emission (CIGALE v.22.1; [18]) to estimate the stellar masses of our sample while accounting for AGN contribution. We use Bruzual & Charlot [20] single stellar population models assuming a Chabrier IMF (Chabrier [27]) and solar metallicity. Furthermore, we consider a delayed star formation history (SFH) model with an optional exponential burst to generate the grid of models. We include the standard nebular emission model from Inoue [69] and the dust attenuation model using the Calzetti et al. [25] attenuation curve. The reprocessed dust emission is modeled adopting the dust emission models of Dale et al. [32], whereas the AGN emission is modeled using the models from Fritz et al. [49]. CIGALE

simultaneously fits the AGN and the galaxy component using all available photometry ($g$, $r$, $z$, W1, W2, W3, and W4) at the spectroscopic redshift of each galaxy. It returns the estimates for the main galaxy properties, including stellar mass, star formation rate, and the relative contribution of the dusty AGN torus to the total infrared luminosity. In cases where the fit is consistent with no or negligible AGN contribution, CIGALE outputs the AGN fraction to be zero. The quality of fit is expressed by reduced $\chi^2$, and the estimates and errors of physical properties are the likelihood-weighted mean and standard deviations of the probability distribution function, respectively ([18]). The physical properties of galaxies and AGN for both fuji and guadalupe releases are generated as value-added catalogs (VACs), and the former is now publicly available with DESI EDR⁴ ([142]).

CIGALE performs better when mid-IR photometry is available, but it can still converge to a reliable solution using confident optical photometry (in $grz$ bands) along with any available WISE photometry ([142]). We therefore do not place any constraints on the WISE photometry of sources (§ 2.2), and instead focus on the following quality cuts to select sources with confident stellar masses:

CIGALE computes two estimates of the different parameters based on the best-fit model (best) and the likelihood-weighted mean measured from the probability density function marginalized over all the parameters (bayes). The comparison between the two measurements is parametrized by FLAG_MASSPDF, which is defined as $\log M_{\mathrm{best}}/\log M_{\mathrm{bayes}}$. To select accurate stellar mass measurements, we further apply the following cut as recommended by the CIGALE VAC (see [115, 142]).

This results in 1,358,177 sources with stellar masses ranging from $\log(M_\star/M_\odot)=6.0-12.3$, with a median of $\log(M_\star/M_\odot)=10.3$. For consistency with similar studies on dwarf AGN candidates, we define galaxies with $\log(M_\star/M_\odot) \le9.5$ as dwarf galaxies throughout this paper. This threshold is approximately equivalent to the stellar mass of LMC ([148]), the largest satellite galaxy of the Milky Way. The left panel of Figure 1 shows the 2D distribution of all 1,358,177 sources in the $\log(M_\star/M_\odot)-z$ space. The source density peaks at $\log(M_\star/M_\odot) \approx10.5$ and $z \approx0.2$. We find that the stellar masses at $z \approx0.45$ extend down to $\log(M_\star/M_\odot) \approx7.5$. There are two overlapping distributions varying with redshift, with the higher-density sample extending from low-redshift dwarf galaxies to high-redshift massive galaxies. The lower-density sample hosts the majority of dwarf galaxies in our sample. This is primarily due to the DESI targeting selection (see § 2.1) as shown on the right panel of the Figure. The top-right panel shows the redshift distribution of sources, divided by their DESI target type. All the target types are uniformly distributed across our redshift range of interest, except the LRG targets that have a redshift distribution skewed toward the high end of the redshift range. The MWS sources are intended to be Milky Way stars by the targeting criteria, but Redrock recovers them as GALAXY or QSO. We also see that BGS targets dominate across all redshifts. This is expected because BGS is designed to focus on low redshifts, while LRG, ELG, and QSO targets focus on $z > 0.5$. The bottom-right panel shows the stellar mass distribution of the sources, distributed between different target types. The dotted vertical line denotes the separation of dwarf galaxies from high-mass galaxies that we consider in this paper. Of the 1,385,177 total sources, we have 222,680 dwarf galaxies (with $\log(M_\star/M_\odot) \le9.5$) and 1,135,497 high-mass galaxies (with $\log(M_\star/M_\odot) > 9.5$). The BGS targets constitute the majority of both dwarf galaxies and high-mass galaxies in our sample. The rest of the dwarf galaxies come from SCND, ELG, and QSO, with <1% sources from LRG and MWS. On the other hand, the remaining high-mass galaxies are primarily LRG and SCND targets and <1% are from ELG, QSO, and MWS. We discuss the impact of these targeting selections on AGN identification in § 4.3.

Emission-Line Measurements

To accurately determine the emission-line ratios for AGN selection and to estimate the BH masses using virial techniques, we need careful measurements of the flux and width ($\sigma$) values of various emission lines. For this purpose, we developed a Python-based emission-line fitting code (EmFit), which uses a non-linear least squares fitting algorithm, for fitting H$\beta$ $\lambda$4861, [O III]$\lambda$5007, [N II]$\lambda$6548,6584, H$\alpha$ $\lambda$6563, and [S II]$\lambda$6716,6731 emission lines. The code includes testing for extra components for [O III], [S II], and H$\alpha$. We run EmFit on all 1,624,108 sources with good spectra, irrespective of photometry and stellar mass measurements. We describe the fitting procedure in this subsection.

For each object, we correct the DESI spectrum for galactic reddening using the Fitzpatrick (1999) dust extinction law. We then obtain the stellar continuum from FastSpecFit⁵ ([116], J. Moustakas et al., in preparation). We use FastSpecFit v3.2 for fuji and v3.1 for guadalupe; the VAC for DESI-EDR is publicly available⁶. FastSpecFit models the stellar continuum using 168 composite stellar population templates of varying age, stellar metallicity, and dust attenuation. The spectra are then corrected for any residual, unmodeled flux by constructing a “smooth” continuum using a sliding median with iterative outlier-clipping. We subtract this total continuum from the observed spectrum and convert the resulting emission-line spectrum to rest-frame.

We have two different fitting modes for the rest-frame spectra: the “default” mode and the “Extreme Broad Line (EBL)” mode. The EBL sources are defined here as sources where the broad component of H$\alpha$ extends up to the [SII] region. For selecting such sources, we compare the blue side of [SII] emission-lines (6600 − 6670 Å) to the red side (6700 − 6900 Å). If the difference between median flux densities from both sides exceeds 0.5 × $10^{-17}$ ergs s$^{-1}$ cm$^{-2}$, and if the flux density on the blue side is at least five times greater than the red side, we flag the source as an EBL candidate. Of the 1,624,108 sources, 3,457 sources (0.21%) are fit via the EBL mode, while the rest are fit via the default mode.

For both modes, we model each emission line using one or more Gaussian functions, which are constrained to have positive flux. We also include a zero-order polynomial to each fit, as a proxy for any residual continuum around the emission lines. Adding an extra Gaussian component to a fit typically improves the fit and decreases its reduced $\chi^2$. To evaluate the statistical significance of this improvement, we compute the $\chi^2$ difference, $\Delta\chi^2$, with and without the extra component. By Wilk’s theorem, the $\Delta\chi^2$ asymptotically approaches a $\chi^2$ distribution under the null hypothesis of a simple fit with a single Gaussian (Wilks [153]). Thus, we convert the $\Delta\chi^2$ into a tail probability ($p$-value), and we reject the null hypothesis when $p$-value is less than 3 × $10^{-7}$. This provides at least a 5$\sigma$ confidence for the existence of an extra component.

We first describe our default fitting method, followed by the EBL fitting method for the flagged sources.

Default Fitting Mode

For the sources that are not flagged as EBL sources, we divide the rest-frame emission-line spectrum into four windows: 1) 4700 Å $\le\lambda_{\mathrm{rest}} \le$ 4930 Å for the H$\beta$ region; 2) 4900 Å $\le\lambda_{\mathrm{rest}} \le$ 5100 Å for the [O III] region; 3) 6300 Å $\le\lambda_{\mathrm{rest}} \le$ 6700 Å for the [N II]+ H$\alpha$ region; and 4) 6630 Å $\le\lambda_{\mathrm{rest}} \le$ 6900 Å for the [S II] region.

We first fit the [S II] doublet with a pair of Gaussian profiles tied to have the same redshift and same intrinsic linewidth in velocity space. In addition to the one-component fit, we also test for a two-component fit to the [S II] doublet, such that the amplitude ratio of the two components is equal for the lines (following Reines et al. [130]). When the spectrum is noisy, there is a risk that the model includes a low-amplitude broad component to fit the continuum. To avoid this unphysical scenario, the width of the second component is allowed to vary up to 600 km s$^{-1}$. However, if in such a scenario, a single component has a width $> 450$ km s$^{-1}$, we revert back to a two-component model. We apply quality cuts based on this scenario to remove unphysical fits (see § 3.2.4). The two-component model is accepted if the $p$-value from the $\Delta\chi^2$ provides a $> 5\sigma$ confidence. In such a case, we label the narrower component as the primary component and the broader component as the secondary component. The second component can either be an “outflow” component or a second narrow peak in double-peaked emission line sources (see Appendix § 8). Only 11,000 ($\approx0.8\%$) sources require a second component. The [S II] emission line profile has been shown to generally match other narrow emission lines (Greene & Ho [53]; Reines et al. [130]). We, therefore, use the final [S II] model as a template for the [N II] and H$\alpha$ emission lines.

We then fit the [N II]+ H$\alpha$ complex, where the redshifts of the [N II] and H$\alpha$ lines are tied to the [S II] lines. The amplitude ratio of [N II]$\lambda$6583/[N II]$\lambda$6548 is fixed to the theoretical value of 2.96 (Acker et al. [1]). When [S II] is best-fit with a single component model, the width of [N II] lines is fixed to the width of [S II] lines in velocity space, and the width of the narrow H$\alpha$ component is initially allowed to increase up to 30% of the [S II] width (similar to Reines et al. [130]). If H$\alpha$ does not satisfy this constraint, the width of the narrow H$\alpha$ is also held fixed to that of the [S II] component. For the two-component [S II] model, all the three emission lines - [N II] doublet and H$\alpha$ - are also considered to have two components such that the width of the respective components is fixed to the [S II] components in velocity space. Additionally, the separation between the two components is also fixed as measured from the [S II] template. Having the same second component as [S II] lines is important to achieve a physically probable broad component after removing the contribution from the second component (see Figure 2 and Appendix § 9). For both one and two-component fits, we further test for a broad component (H$\alpha$;b) with a minimum FWHM (H$\alpha$;b) $> 300$ km s$^{-1}$, based on the $p$-value from the $\Delta\chi^2$ with and without a broad component. The broad component is accepted when it delivers $> 5\sigma$ confidence. However, this component is rejected if this component is narrower than the second component. We use the resulting H$\alpha$ fit as a template for the H$\beta$ fit and the number of components of H$\beta$ is fixed to be the same as H$\alpha$. The redshifts and widths of all the available H$\beta$ components (primary, second, and broad) are fixed to the widths of the respective H$\alpha$ components in velocity space. Only the amplitude of the components is allowed to vary.

Given that the [O III] profile can exhibit a blue-shifted outflow component and does not match the profile of other narrow emission lines, we fit the [O III] doublet separately. We fix the relative redshifts of the peaks, fix the amplitude ratio of [O III]$\lambda$5007/[O III]$\lambda$4959 to 2.98 (Dimitrijević et al. [41]), and add the condition that the widths of the two lines are equal in velocity space. Both one and two-component models are tested and the two-component fit is accepted if the $p$-value from the $\Delta\chi^2$ satisfies the $5\sigma$ confidence cut. Similar to the [S II] model, the narrower component is taken as the primary component, and the broader component can either be an outflow component or a second component (see Appendix § 8).

Figure 2 shows an example of a dwarf galaxy ($\log(M_\star/M_\odot) \approx9.4$) spectrum fit via the default mode. The top panel shows the best-fit of the stellar continuum (in blue) from FastSpecFit overplotted on the DESI spectrum (in gray). The resulting emission line spectrum (in black) and the best-fit models (in red) are shown in the bottom panel. We see that all the emission lines contain a second component (a possible outflow in this case). The primary and outflow Gaussian components are shown in orange, and the reduced $\chi^2$ values are shown in each of the panels. This galaxy is a dwarf AGN candidate (see § 4) and shows the importance of adding a second component in the narrow emission lines.

Extreme Broad Line Fitting Mode

For the sources that are flagged as EBL sources, we divide the rest-frame emission-line spectrum into two windows: 1) 4700 Å $\le\lambda_{\mathrm{rest}} \le$ 5100 Å for the H$\beta$ + [O III] region; and 2) 6300 Å $\le\lambda_{\mathrm{rest}} \le$ 6900 Å for the [N II]+ H$\alpha$ + [S II] region.

We first fit the [N II]+ H$\alpha$ + [S II] emission-lines with a single Gaussian component for each of the narrow lines and an extra broad Gaussian component. In the case of EBL sources, multiple broad components might be needed for properly fitting the complex emission line profile ([54, 92]). However, we approximate the complex spectrum with a single broad Gaussian component for this work. The amplitude ratio of [NII]$\lambda$6583/[NII]$\lambda$6548 is again fixed to the theoretical value of 2.96 ([1]). The redshifts of all the narrow components ([NII] doublet, H$\alpha$, and [SII] doublet) are tied together and their widths are constrained to be equal in velocity space. In some sources, a noisy spectrum is selected as an EBL source. Their fits using this mode lead to an unrealistic broad H$\alpha$ component. Using visual inspection, we find that such sources are poorly fit with broad H$\alpha$ widths $< 1000$ km s$^{-1}$. We therefore use this criterion ($\sigma$(H$\alpha$; b) $< 1000$ km s$^{-1}$) to pick such sources and rerun their fitting using the default mode. The median FWHM of broad H$\alpha$ components for the sources that are fit with the EBL mode is $\approx2000$ km s$^{-1}$.

For the H$\beta$ + [OIII] region, we use the H$\alpha$ model as a template for the H$\beta$ model and the [OIII] lines are fit independently. We fix the redshifts and widths of the H$\beta$ components to scale with the H$\alpha$ components. The redshifts, amplitudes, and widths of the [OIII] components are tied as mentioned in the default mode. Both one- and two-component models are tested for the [OIII] lines, and the two-component fit is allowed if the $p$-value from the $\Delta\chi^2$ is less than the $5\sigma$-confidence cut.

Figure 3 shows an example of the best-fit of a high-mass EBL source ($\log(M_\star/M_\odot) \approx11.6$). The panel shows the total continuum from FastSpecFit in blue, overplotted on the spectrum. The bottom panels show the best-fits in the H$\beta$ + [OIII] and [NII] + H$\alpha$ + [SII] regions. The narrow and outflow components are shown in orange, while the broad Balmer components are shown in blue. The [OIII] lines in this galaxy clearly show evidence of outflows.

Measurements and Uncertainties

We compile the amplitude, mean, and standard deviation of all the detected Gaussian components from the best-fit model. The fluxes and widths ($\sigma$) are computed using these parameters and the resulting width is corrected for instrumental resolution by subtracting the resolution element in quadrature. To estimate the fit uncertainties, we consider a Gaussian error distribution with standard deviation derived from the inverse variance at each pixel of the spectrum. We vary the spectrum within this distribution and repeat the fitting procedure a hundred times. The uncertainties are estimated as the standard deviation of all measurements from these iterations. We assume that the pixels are uncorrelated and the errors are therefore underestimated. In addition to the flux and $\sigma$ measurements, we also measure the noise around the emission lines as the root-mean-square (rms) of the continuum around the lines. We also make note of PROB_BROAD as the fraction of iterations in which the broad H$\alpha$ component is measured with non-zero flux. This is especially important for selecting broad-line candidates when the spectra are noisy.

The emission line ratio measurements require total fluxes from the narrow component of the emission lines (§ 4). As mentioned earlier, some of the second components in [S II], [O III], [N II], and H$\alpha$ appear as narrow components for double-peaked emission line sources. We add the fluxes from the two components for such sources. The method of picking these double-peaked sources is described in Appendix § 8.

Applying Quality Cuts

Even after the detailed careful fitting process, the fits to some noisy spectra lead to poor representations of the observations and unphysical flux and sigma values. One such instance is when the stellar continuum subtraction is erroneous and creates a non-physical jump between the continuum on either side of the emission lines. This is especially prominent in fitting [S II] lines, leading to broad second components. We remove these sources by applying the following cut:

where $\sigma_{\mathrm{[SII]},2}$ is the width of the second [S II] component. Another issue is related to the forced fitting of a negative continuum and an extremely high-velocity width narrow component to account for the noisy spectrum. These affect all the fitting windows and we remove them by applying an upper limit on the width of the primary components:

These quality cuts result in 1,350,112 sources across the stellar mass and redshift range, which we use for our analysis.

Selecting Line-Emitting Galaxies

We explore two complementary ways of detecting AGN candidates in this paper: emission-line ratios and broad H$\alpha$ detection (§ 4). This depends on the accurate detection of four emission lines: H$\beta$, [O III]$\lambda$5007, H$\alpha$, and [N II]$\lambda$6583. We apply the following SNR and AoN (amplitude-over-noise) cuts on these four lines for selecting the sample of line-emitting galaxies:

It is often found that H$\beta$ lacks significant detection in galaxies and setting a higher threshold often removes good candidates that would have been classified as AGN [30]. We therefore use a lower threshold for H$\beta$ compared to the other emission lines. For double-peaked sources, we consider the fluxes as the sum of the fluxes of the two components, and their errors are added in quadrature for computing the SNR values. We still consider the AoN of the primary H$\beta$ component for these sources as a way of ensuring that the components are detected over the noise of the spectrum. We have 459,208 sources that satisfy the SNR criteria for [O III], H$\alpha$, and [N II], of which 48,451 sources are rejected based on the H$\beta$ cuts. This results in 410,757 line-emitting galaxies, of which 114,496 are dwarf galaxies (51.6%), and the rest of 296,261 (26.3%) are high-mass galaxies. We observe broad H$\alpha$b detection in 13,526 of these sources. To select confident H$\alpha$b candidates, we apply the following cuts:

Finally, we apply a cut on PROB_BROAD for selecting statistically likely broad component detections:

These cuts result in 6,185 broad-line (BL) candidates, and the remaining 404,572 line-emitting galaxies are considered narrow-line (NL) sources.

Incidence of AGN in Line-Emitting Galaxies

In this section, we search for AGN signatures in our sample of line-emitting high-mass and dwarf galaxies (§ 4.1). We further compare our identified dwarf AGN candidates with results from prior single-fiber spectroscopic surveys (§ 4.2). Lastly, we explore the variation of the observed AGN fraction with stellar mass and redshift (§ 4.3).

AGN Candidates from the [Nii]-BPT Diagram

Two primary processes contribute to the photoionization of the interstellar medium (ISM) in galaxies: star formation and AGN activity. AGN produce a harder spectrum (i.e., a larger fraction of high-energy photons) than hot stars, causing identifiable signatures in the resulting emission line spectrum. Two-dimensional emission-line diagnostic diagrams, commonly known as the BPT diagrams [11, 149], can be used to distinguish the primary

source of ionization in galaxies and identify AGN. In this paper, we focus on the most widely used BPT diagram with $[\mathrm{N\,II}]/\mathrm{H}\alpha$ vs. $[\mathrm{O\,III}]/\mathrm{H}\beta$ (hereafter $[\mathrm{N\,II}]$-BPT Diagram). We explain our choice to use only this specific diagnostic in § 6.1.

Using pure stellar photoionization models, Kewley et al. [79] suggested a “maximum starburst line” (dashed line in all panels of Figure 4), above which a non-stellar ionization source, likely an AGN, is necessary to account for the emission-line ratios. The solid line in all the panels of Figure 4 separates the sources with AGN contribution from the empirical branch attributed to star-forming galaxies (Kauffmann et al. [75]). Galaxies located between these two lines contain contributions from both AGN and star formation and are typically referred to as “composite” galaxies. For this study, we classify all sources in the AGN-dominated and composite regions as AGN candidates.

The $[\mathrm{N\,II}]$-BPT diagnostic is metallicity-sensitive, with lower metallicity systems occupying regions with higher $[\mathrm{O\,III}]/\mathrm{H}\beta$ and lower $[\mathrm{N\,II}]/\mathrm{H}\alpha$ ratios ([144, 26]). Consequently, low metallicity AGN can occupy the same regions as low metallicity starbursts (Groves et al. [60]).

Initial BPT Selection

We categorize the line-emitting galaxies in our sample into two groups: those exhibiting only narrow components (NL candidates) and those with a secure broad H$\alpha$ detection (BL candidates). While our primary goal is to identify AGN candidates, we also aim to separate the BL candidates where the presence of a broad H$\alpha$ component is indicative of the broad-line region (BLR) associated with an AGN. This will facilitate the estimation of their BH masses (§ 5.1) and the study of the $M_{\mathrm{BH}} - M_{\star}$ scaling relation (Section 5.3). Thus, we focus on identifying AGN and composite sources in both NL and BL candidates (§ 3.3) via the BPT diagnostic diagram.

Figure 4 shows the $[\mathrm{N\,II}]$-BPT diagram of high-mass galaxies (Top Panel) and dwarf galaxies (Bottom Panel) in our sample. The left panels show the BPT diagram for the NL candidates, while the right panels show the BPT diagram for the BL candidates constructed using the line ratios of the narrow components of the emission lines.

The shape of the bivariate distributions on the BPT diagram is different for dwarf and high-mass galaxies. The high-mass galaxies are characterized by a thin star-forming branch and a well-extended AGN branch for both NL and BL candidates. On the other hand, low-mass galaxies have a dense star-forming branch that extends to lower $[\mathrm{N\,II}]/\mathrm{H}\alpha$ ratios ($< -1.4$) than their counterparts. The dwarf galaxy diagrams also have sparsely populated AGN branches for both NL and BL candidates. This is not surprising as dwarf galaxies are typically star-forming galaxies and their star-forming signatures can therefore dominate the emission-line signatures. In addition, the low metallicity of dwarf galaxies leads to low $[\mathrm{N\,II}]$ and high $[\mathrm{O\,III}]$ values ([94, 67]), which explains the concentration of dwarf galaxies at the top-left corner of the star-forming branch ([144, 26]).

Of the 296,261 high-mass line-emitting galaxies (both NL and BL candidates) in our sample, we find that 24,227 (8.2%) sources are AGN-dominated, 51,701 (17.4%) are composites, and 220,333 (74.4%) are star-forming galaxies. In contrast, we find 405/114,496 (0.3%), 2,039/114,496 (1.8%), and 112,052 (97.9%) of the dwarf line-emitting galaxies (both NL and BL candidates) reside on the AGN-dominated, composite, and star-forming branch, respectively. Table 1 describes the number of NL and BL candidates in each of the BPT regions from Figure 4.

	Dwarf ($\log(M_\star/M_\odot) \le9.5$) Galaxies		High-Mass ($\log(M_\star/M_\odot) > 9.5$) Galaxies
All Sources	221,778		1,128,334
	NL Candidates	BL Candidates	NL Candidates	BL Candidates
Before Visual Inspection$^{a}$
Line Emitting Galaxies	113,655	841	290,917	5,344
AGN Dominated Sources	350 (0.3%)	55 (6.5%)	21,322 (7.3%)	2,905 (54.4%)
Composites	1,933 (1.7%)	106 (12.6%)	50,330 (17.3%)	1,371 (25.6%)
Star-Forming	111,372 (98.0%)	680 (80.9%)	219,265 (75.4%)	1,068 (20.0%)
After Visual Inspection$^{b}$
AGN Dominated Sources	365 (0.3%)	40 (4.8%)	21,405 (7.4%)	2,822 (53.6%)
Composites	1,951 (1.7%)	88 (10.7%)	50,470 (17.3%)	1,231 (23.6%)

Visual Inspection of BL Candidates

Of the total 6,185 BL candidates in our sample (both high-mass and dwarf galaxies), 1,748 (28.2%) lie on the star-forming branch, 1,477 (23.9%) are composites, and 2,960 (47.9%) are AGN-dominated galaxies (see Table 1). The broad H$\alpha$ emission in the BL candidates can originate from several processes within the galaxy. For instance, infalling gas in the broad-line region (BLR) of the AGN can produce this emission. Additionally, dynamics in the accretion disk surrounding the BH, outflows driven by AGN activity, and/or star formation can contribute extra components to the H$\alpha$ emission line. Stellar processes, such as Type II supernovae and Luminous Blue Variables (LBVs), may also generate broad H$\alpha$ emission ([143, 62, 121]). To distinguish between the star-formation-driven and AGN-driven broad H$\alpha$ sources, we focus exclusively on the BL candidates identified as AGN and composites in the $[\mathrm{N\,II}]$-BPT diagram.

Since we push to lower widths of the broad H$\alpha$ component than any prior studies ([130, 131, 145, 135]), we aim to ensure the confident identification of such low-width broad components. For this purpose, we visually inspect 556 AGN candidates that have $\mathrm{FWHM}_{\mathrm{H}\alpha;\mathrm{b}} < 1000\ \mathrm{km\,s}^{-1}$ (Appendix § 9).

After a thorough visual inspection, we identify 146 sources exhibiting a clear broad H$\alpha$ component. We consider these, along with the rest of the BL AGN-dominated and composite candidates with FWHM$_{\mathrm{H}\alpha;\mathrm{b}} > 1000$ km s$^{-1}$ as confident BL-AGN candidates. Additionally, we classify 154 sources with a statistical broad H$\alpha$ component, as tentative BL-AGN candidates, although their broad components are less clear visually.

Of the remaining 256 visually inspected BL candidates, we identify 252 sources that display a visual second/outflow component in the [N II] line but are not detected by EmFit. We suspect this explains the additional H$\alpha$ component rather than a true broad component. Additionally, four sources exhibit particularly complex emission line profiles around the [N II] + H$\alpha$ emission lines, with multiple overlapping components, preventing us from constraining the presence of a broad component. Consequently, we do not classify these 256 galaxies as BL-AGN and relabel them as NL-AGN. Even though we do not visually inspect any BL candidates that are present on the star-forming branch, it does not affect the percentages reported in § 4.1.1 or any subsequent analysis in this paper. We present details of our visual inspection in Appendix § 9.

Final BPT Selection

Of the 296,261 high-mass line-emitting galaxies, we have a final sample of 71,875 NL-AGN and 4,053 BL-AGN candidates (AGN Fraction $\approx25.6\%$). These numbers are influenced by the selection criterion of line-emitting galaxies (§ 3.3). Given that robust detection of emission lines is required to appear on the BPT diagram, some high-mass quiescent galaxies hosting an AGN may be missed by this diagnostic. The spectra of such galaxies are typically dominated by strong absorption lines from stars, weakening the emission-line contribution from the AGN and resulting in a drop of high-mass AGN candidates detected via the [N II]-BPT diagram.

Of the 114,496 dwarf line-emitting galaxies, we have a final sample of 2,316 NL-AGN and 128 BL-AGN candidates (AGN Fraction $\approx2.1\%$). The majority of dwarf galaxies reside on the star-forming branch. This does not necessarily mean that no AGN is present in these sources. The high star formation in low-mass galaxies can dilute the contribution from AGN ([146]). There might be a few dwarf AGN candidates on the star-forming branch that may be detected via other diagnostics ([113, 135]). In fact, Birchall et al. [15] found that 85% of their X-ray-confirmed dwarf AGN candidates fail the BPT criterion for AGN. However, the dwarf AGN candidates that are detected from the [N II]-BPT diagram have been observed to be robust AGN candidates ([130, 9, 135]). We note that we are attempting to detect confident AGN candidates in our sample of galaxies, but do not claim to be complete. Our estimates therefore provide lower limits to the incidence of AGNs in galaxies. Table 1 provides more detailed information regarding our emission-line classification.

AGN in Line Emitting Dwarf Galaxies

Our result of 2,444 dwarf AGN candidates more than triples the existing census of optically selected dwarf AGN candidates to date, and is based on using just the SV data and 20% of DESI Year 1 data. By extrapolating to the entire DESI footprint, we expect to conservatively find $\gtrsim$10,000 dwarf galaxies with detectable AGN signatures within the redshift range of 0.001 $\le z \le$ 0.45, after the completion of the five-year survey.

The observed dwarf AGN fraction in this work ($\approx$2.1%) is higher than any previous study using single-fiber integrated optical spectra (AGN fraction $<$1%; [130, 135] or using infrared and X-ray searches (AGN Fraction $<$1%; [106, 93, 15, 87, 16, 23]. In this subsection, we compare our dwarf AGN candidates with those detected using the BPT diagnostics from prior single-fiber spectroscopic surveys i.e., the Sloan Digital Sky Survey (SDSS; York et al. [157]) and Galaxy and Mass Assembly (GAMA; Driver et al. [43]) survey.

Figure 5 shows the distribution of dwarf AGN candidates in our sample (blue circles), alongside candidates identified by Reines et al. [130] (with SDSS; pink triangles) and Salehirad et al. [135] (with GAMA; black squares) in the $\log(M_{\star}/M_{\odot}) - z$ space. The top and right panels display the distribution of redshift and stellar masses of these samples, respectively. All stellar masses have been estimated using the Chabrier IMF.

Reines et al. [130] conducted the first systematic search for AGN in dwarf galaxies using the SDSS DR8 spectroscopic catalog (Aihara et al. [2]), covering an area of $\approx9200\ \mathrm{deg}^{-2}$. Using the [N ii]-BPT emission-line diagnostic, they found 136 candidates with AGN signatures down to $\log(M_{\star}/M_{\odot}) = 8.1$ and out to $z = 0.055$, with a median $\log(M_{\star}/M_{\odot}) = 9.3$ and a median redshift of $z = 0.028$. More recently, Salehirad et al. [135] applied various spectroscopic diagnostics on galaxies from the GAMA DR4 (Driver et al. [44]), spread over an area of $\approx250\ \mathrm{deg}^{-2}$. They found 35 dwarf AGN candidates down to $\log(M_{\star}/M_{\odot}) = 7.99$ and out to $z = 0.2$. These candidates have a slightly higher median redshift ($z = 0.077$) compared to Reines et al. [130], but have a similar median $\log(M_{\star}/M_{\odot}) = 9.3$.

DESI has greatly expanded the discovery space of optical dwarf AGN candidates. Our candidates extend down to $\log(M_{\star}/M_{\odot}) \approx6.12$, almost 1.5 dex lower than both SDSS and GAMA surveys, and also extend out to $z \approx0.45$. The early DESI data covers a sky area of $\approx4,400\ \mathrm{deg}^{-2}$, nearly half of the SDSS DR8 coverage ($\approx9,200\ \mathrm{deg}^{-2}$). The larger number of dwarf AGN candidates from our study is partly due to a higher target density of DESI compared to the previous spectroscopic sky surveys. We explore the selection effects compared to Reines et al. [130] in more detail in § 6.2.

AGN Fraction in Galaxies

The fraction of dwarf galaxies hosting a BH, i.e., their BH occupation fraction, is a key diagnostic in understanding the BH seed formation mechanisms ([150, 104]). Some of these BHs can be detected as an AGN via different multi-wavelength methods. Therefore, the AGN fraction acts as a lower limit for the BH occupation fraction. We estimate the AGN fraction in DESI sources based on the fraction of galaxies identified as BPT-AGN candidates. As is common in previous optical studies of dwarf galaxies ([130, 135, 107, 109]), we do not apply any completeness corrections. However, we will consider the effect of accounting for the galaxies that fail the emission line detection criteria from § 3.3.

To investigate trends with galaxy properties, we start by computing the AGN fraction in line-emitting galaxies in bins of stellar mass and redshift (Figure 6). The size of the squares shown in the figure is proportional to the number of galaxies within those bins, while their color reflects the AGN fraction as shown on the color bar.

We find that the BPT-AGN fraction, on average, increases with increasing stellar mass, across the entire redshift range, which is consistent with previous optical emission-line studies (e.g., Juneau et al. [71]). The BPT-AGN fraction is $\lesssim10\%$ at low stellar masses ($\log(M_{\star}/M_{\odot}) \lesssim9$), and reaches $\approx100\%$ at high stellar masses ($\log(M_{\star}/M_{\odot}) \gtrsim11$). At $\log(M_{\star}/M_{\odot}) \approx10$, the BPT-AGN fraction appears to decrease with increasing redshift. The change with redshift is unclear at lower and higher masses than this region. We discuss the variation of the BPT-AGN fraction with stellar mass and redshift, along with various selection effects below.

Variation with Stellar Mass

In Figure 7, we show the observed BPT-AGN fraction as a function of stellar mass. We arrange all the sources in order of increasing stellar masses and divide them into bins of nearly equal number of sources ($\approx54,000$ sources in each bin). We compute two different BPT-AGN fraction estimates: the fraction of BPT-AGN candidates in line-emitting galaxies (shown as red squares in Figure 7); and the fraction of BPT-AGN candidates in all galaxies, including the galaxies without line emission (shown as black circles in Figure 7). Even though some non-line-emitting galaxies may host a BH, we assume that they do not host an AGN. Therefore, we are reporting the lower limits of the total BPT-AGN fraction. Within each bin, we compute the median and standard deviation of stellar masses as the representative stellar mass and error. The right panel of Figure 7 shows the variation of the BPT-AGN fraction in the dwarf galaxy regime. In this case, we divide all the dwarf galaxies into bins of $\approx27,700$ sources and compute both versions of the BPT-AGN fraction. The horizontal dashed-dotted red line denotes the observed average dwarf AGN fraction in line-emitting galaxies (§ 4.2). This fraction is nearly 4 times higher than the pink dashed line that denotes the average dwarf AGN fraction from Reines et al. [130].

Focusing on line-emitting galaxies, we observe the BPT-AGN fraction sharply rises from $\approx1.2\%$ at $\log(M_{\star}/M_{\odot}) \approx8.0$ to $\approx93.5\%$ at $\log(M_{\star}/M_{\odot}) \approx11.4$. While it is well established that nearly all massive galaxies host BHs (Kormendy & Ho [81]), only a subset of these BHs are actively accreting at any given time. Furthermore, strong selection effects are present in this regime, as many sources fail to be detected in all four BPT lines. As a result, there is an increasing bias towards the detection of emission lines linked to the presence of AGN with increasing stellar mass.

At the low-mass end, whether BHs are present at the centers of all dwarf galaxies has not yet been established. The BPT-AGN fraction in dwarf galaxies (right panel of Figure 7) increases from $\approx1.2\%$ at $\log(M_{\star}/M_{\odot}) \approx7.6$ to $\approx3.6\%$ at $\log(M_{\star}/M_{\odot}) \approx9.4$. It drops to $<2.1\%$ (red horizontal line) by $\log(M_{\star}/M_{\odot}) \approx9$, which is the average dwarf AGN fraction. This suggests that the majority of dwarf AGN candidates from the BPT selection criterion are dwarf galaxies with $\log(M_{\star}/M_{\odot}) > 9$, which is what we see in Figure 6. An important selection bias in this mass regime is that for lower-mass systems, the detected AGN sources are likely accreting at higher Eddington ratios, and hence only the more extreme sources would be selected in dwarf galaxies compared to massive galaxies.

When we include non-line-emitting galaxies that do not have significant emission-line detection (black points in Figure 7), we observe that the BPT-AGN fraction increases from $\approx0.3\%$ at $\log(M_{\star}/M_{\odot}) \approx7.8$, peaks at $\approx9.7\%$ at $\log(M_{\star}/M_{\odot}) \approx10.5$, turnovers and decreases to $\approx2\%$ at $\log(M_{\star}/M_{\odot}) \approx11.4$. The discrepancy between the two estimates of the BPT-AGN fraction is due to the fraction of line-emitting galaxies selected at different stellar masses (see top panel of Figure 8). In the stellar mass range of $9 \leq\log(M_{\star}/M_{\odot}) \leq10$, the difference is small, and $>60\%$ of the sources have strong emission-line detection. However, there is an increase in the divergence between the line-emitting and total BPT-AGN fraction (red and black curves) at higher stellar masses. This divergence occurs because massive galaxies tend to be metal-rich and have weak [OIII] emission (Kauffmann et al. [75]), often failing the selection criterion. In contrast, the differences in the low-mass regime are primarily driven by the [NII] selection criterion, which is expected due to the low metallicity of galaxies in this regime.

The increase in the AGN fraction with increasing stellar mass is similar to some previous studies focusing on X-ray-detected AGN ([106, 15]). However, other X-ray and multiwavelength studies have found weak to no variation with stellar mass ([72, 16] (Birchall et al. 2023)), suggesting that different AGN selection criteria have different selection biases.

Variation with Redshift

We found that the BPT-AGN fraction strongly depends on stellar mass. To examine the variation of BPT-AGN fraction with redshift, we group all the dwarf galaxies and sub-divide high-mass galaxies into stellar mass bins of 0.5 dex. We then organize these sources in increasing order of redshift and calculate the fraction of line-emitting galaxies, as well as the line-emitting and total BPT-AGN fractions, within redshift bins containing nearly equal number of sources. The top panel of Figure 8 illustrates the fraction of line-emitting galaxies across different stellar mass selections as a function of redshift, while the middle and bottom panels display the variations of line-emitting and total BPT-AGN fractions with redshift for these groups, respectively. The horizontal red line in the middle panel is the observed AGN fraction in dwarf galaxies.

A noticeable dip in the estimated fractions is evident in all the panels around $z \approx0.18$. This dip arises from two issues concerning the [O III] emission line. In this redshift range, the [O III] line coincides with the overlap region of the B and R cameras of the DESI spectrograph (shown by grey region in Figure 8), and also one of the skylines (the Na I doublet). As a result, the detection of this emission line becomes very noisy, often resulting in the SNR ([O III]) failing to meet the selection criterion in this redshift range.

In the dwarf galaxy regime (black line), the fraction of line-emitting galaxies remains relatively constant at $\approx50\%$ up to $z \approx0.2$, before declining to $\approx20\%$ at $z \approx0.37$. The BPT-AGN fraction for line-emitting dwarf galaxies fluctuates around the overall dwarf AGN fraction ($\approx2.1\%$; dash-dotted line) and remains nearly constant throughout the entire redshift range. When considering all dwarf galaxies and assuming that the non-line-emitting galaxies do not host an AGN, the BPT-AGN fraction drops to $\lesssim1\%$ but remains nearly flat with redshift. Thus, the dwarf AGN fraction is almost a factor of two higher when considering only line-emitting galaxies versus all galaxies but with no significant redshift dependence.

When we focus on massive galaxies (colored lines), we find a clear dependence of the trends on stellar mass. As stellar mass increases, the fraction of line-emitting galaxies decreases across all redshifts. Conversely, the BPT-AGN fraction for line-emitting galaxies (middle panel) exhibits an opposite trend and increases with stellar mass across all redshifts. When considering all galaxies (bottom panel), there is no longer a monotonic trend of AGN fraction with stellar mass. However, the trend with redshift behaves similarly between the mass bins: we find that the total BPT-AGN fraction initially decreases, but is nearly constant at $z \gtrsim0.2$.

Overall, the observed trends are influenced by various selection effects including the emission-line detection at different stellar masses. Moreover, the DESI targeting selection also introduces significant biases (§ 2.1). The BGS sources are uniformly distributed across all stellar masses and redshifts. However, the rest of the high-mass galaxies are predominantly LRG targets, which primarily target high-redshift sources. The SCND targets also include a majority of high-mass and dwarf galaxies in our sample, but their distribution across redshifts is irregular due to different selection criteria (see Figure 1). These variations in targeting complicate the line-emitting galaxy fraction and BPT-AGN fraction as a function of redshift. Detailed analyses involving selection functions and incompleteness corrections are necessary for a thorough physical interpretation, which we will address in a future publication.

Black Hole Mass - Stellar Mass Scaling Relation

BH-galaxy scaling relations hold clues about BH seed formation mechanisms and also provide a way to study the co-evolution of BHs and their host galaxies. In this section, we study the $M_{\mathrm{BH}} - M_{\star}$ scaling relation of the BL-AGN candidates and compare it with previous studies.

Black Hole Masses

The broad H$\alpha$ emission line can be used to compute BH masses using single-epoch virial techniques. Under the assumption that the BLR gas around the BH is virialized, the kinematics of the gas can be used as a dynamical tracer of BH mass ($M_{\mathrm{BH}}$):

where $R$ is the radius of the BLR around the BH, $V$ is the average velocity of the gas, and $G$ is the gravitational constant. The average velocity of the gas can be inferred using the width of the broad H$\alpha$ emission line, while the radius of the BLR is approximated using the radius-luminosity relation.

Using the same approach as previous studies (Reines et al. [130]; Reines & Volonteri [131]; Salehirad et al. [135]), we use the Greene & Ho [54] formula for estimating the BH masses, with the modified radius-luminosity relationship of Bentz et al. [13] as derived by Reines et al. [130]:

where $L_{\mathrm{H}\alpha;\mathrm{b}}$ and $\mathrm{FWHM}_{\mathrm{H}\alpha;\mathrm{b}}$ are the luminosity and FWHM of the broad H$\alpha$ component and $\epsilon$ is the scale factor, which spans a range of $\sim0.75$–$1.4$ (Greene & Ho 2007; Grier et al. [59, 58]). Here, we assume $\epsilon= 1$ to compare with previous results (Reines & Volonteri [131]).

We compute BH masses based on the EmFit measurements of the broad H$\alpha$ flux and line width (§ 3.2.3) and propagate their uncertainties. The median error for our sample of BL-AGN candidates is $\approx0.02$ dex. However, virial BH mass measurements from single-epoch spectroscopy can have several systematic uncertainties on the order of 0.5 dex (McGill et al. [101]; Shen [139]). We therefore add 0.5 dex to our estimated uncertainties in quadrature to get the final uncertainties in $M_{\mathrm{BH}}$.

Figure 9 presents the distribution of BH masses for the 4,027 confident BL-AGN candidates (in black) and 154 tentative BL-AGN candidates (in hatched magenta) from our sample. The BH masses of the confident BL-AGN sources range from $\log(M_{\mathrm{BH}}/M_{\odot}) = 4.4$ to 9.2, with a median of 7.2. On the other hand, the tentative BL-AGN candidates are concentrated at the lower end of the BH mass range, with values between $\log(M_{\mathrm{BH}}/M_{\odot}) = 4.3$ and 6.3, with a median of 5.3.

Among the BH candidates, 151 confident and 147 tentative BL-AGN candidates have BH masses of $\log(M_{\rm BH}/M_\odot) \le6$, identifying them as potential candidates for IMBHs⁹. This constitutes the largest sample of IMBHs identified to date, almost doubling their existing census ([29, 51]).

Dependency on Galaxy Morphology

Combining BH masses (§ 5.1) and host galaxy stellar masses (§ 3.1), we can now build the $M_{\rm BH} - M_\star$ scaling relation. First, we assess whether this relation depends on host galaxy morphology by dividing our sample of BL-AGN based on the TYPE and SERSIC columns from the LS DR9 catalog (§ 2.2). Point sources are defined based on MORPHTYPE = PSF, while the late-type sources are selected if the MORPHTYPE is either a round-exponential (REX) or an exponential (EXP) or a Sèrsic (SER) with a Sèrsic index $\le2$. The rest of the candidates that do not have PSF or late-type morphologies are considered early-type galaxies. We show the bivariate distribution of these sources in the $M_{\rm BH} - M_\star$ space in Figure 10, along with the empirical fits by Reines & Volonteri [131] and Suh et al. [145].

Reines & Volonteri [131] studied the $M_{\rm BH} - M_\star$ scaling relation by considering 262 BL-AGNs from SDSS DR8 spectroscopic catalog, 15 reverberation-mapped AGNs from Bentz & Katz [12], and 79 inactive galaxies with dynamical BH masses from Kormendy & Ho [81] at $z \le0.055$. They identified two distinct fits with different slopes for active and inactive galaxies, which correlate with galaxy morphology. Specifically, the late-type galaxies are located around the fit for active galaxies, while the early-type galaxies align with the fit for inactive galaxies. Building on this, Suh et al. [145] combined the local AGN sample from Reines & Volonteri [131] with an additional 100 BL-AGNs from the COSMOS field, extending their analysis to $z \approx2.5$. They find a comparatively steeper fit, with their high-redshift sample predominantly consisting of high stellar mass ($\log(M_\star/M_\odot) \gtrsim11$) and high BH mass ($\log(M_{\rm BH}/M_\odot) \gtrsim7$) sources.

From Figure 10, we find that the late-type and early-type galaxies occupy similar regions in the $M_{\rm BH} - M_\star$ space, while the point sources lie $\approx1$ dex above them. Both the distributions for $\log(M_\star/M_\odot) \ge10$ are consistent with active galaxy fit (late-type fit) from Reines & Volonteri [131], with the sources with $\log(M_{\rm BH}/M_\odot) \le6$ falling below it.

Galaxies with a point source morphology lie toward higher values of BH masses and do not appear to follow the late-type (active) fit from Reines & Volonteri [131] nor the Suh et al. [145] fit. However, the early-type galaxy fit from Reines & Volonteri [131] overlaps with the distribution of these sources. Visual inspection of their spectra reveals quasar-like characteristics, featuring extremely broad emission lines and weak to negligible stellar continuum. This makes their stellar mass estimates highly uncertain. Keeping this limitation in mind, we find that 43 dwarf galaxies with PSF morphologies host BHs with $\log(M_{\rm BH}/M_\odot) \gtrsim5$. If their stellar masses were accurate or possibly over-estimated, these sources could represent low-redshift counterparts to the over-massive BHs identified by JWST ([65, 96, 147]). We discuss this further in § 6.3.

To analyze and find empirical fits to the $M_{\rm BH} - M_\star$ scaling relation, we exclude the 518 point sources due to the uncertainty with their stellar mass estimates. Given that the rest of the BL-AGN candidates are consistent with the Reines & Volonteri [131] fit based on active galaxies, we will focus solely on comparing our empirical fit to this late-type (active) fit in the following subsection.

Empirical Fit to the $M_{\mathrm{BH}}$- $M_\star$ Scaling Relation

We use the 3,633 BL-AGN candidates (3,517 confident and 146 tentative), excluding the point sources, to study the $M_{\mathrm{BH}} - M_\star$ scaling relation. Figure 11 shows the bivariate distribution of the confident BL-AGN (left panel) and all BL-AGN sources (including tentative candidates; right panel) in the $M_{\mathrm{BH}} - M_\star$ space.

Focusing on confident BL-AGN candidates, we observe a significant spread in BH masses for a given stellar mass (nearly 3 dex at $\log(M_\star/M_\odot) \approx10.5$). Compared to previous local measurements of the relation Reines & Volonteri [131]; Greene et al. [56], the scaling relation now extends to lower galaxy and BH masses, down to $\log(M_\star/M_\odot) \approx8.5$ and $\log(M_{\mathrm{BH}}/M_\odot) \approx4.4$. When we include the tentative candidates, we also see a similar spread in BH masses at lower mass galaxies ($\approx3$ dex at $\log(M_\star/M_\odot) \approx9.2$). We find a total of 86 dwarf BL-AGN (64 confident and 22 tentative) candidates, extending up to $z \le0.45$. These together represent 5 times more candidates than the 14 low-mass candidates presented in Reines & Volonteri [131].

To study the correlation between $M_{\mathrm{BH}}$ and $M_\star$, and to compare with previous studies, we fit a linear model for the two samples separately, using the Bayesian approach by Kelly [78].¹⁰ We parametrize the relation similar to Reines & Volonteri [131], as follows:

Considering only the confident BL-AGN candidates, we find:

When we consider all the BL-AGN candidates, we find:

We overplot these fits as solid and dashed red lines, respectively, in Figure 11. We see that including the tentative candidates leads to a slightly steeper slope but with a similar normalization compared to the fit with only confident candidates. We find that our fits exhibit a similar slope but a higher normalization compared to the Reines & Volonteri [131] fit (black line). In contrast, the fit from Suh et al. [145] (orange line) has a steeper slope relative to both our fit and the Reines & Volonteri [131] fit. The primary divergence arises at the high-mass end, which is influenced by selection biases. Specifically, Suh et al. [145] focuses on massive BHs selected based on their strong X-ray emission and extremely broad H$\alpha$ lines, while Reines & Volonteri [131] do not include such sources in their sample. The closest counterparts of such sources in our sample are the point sources that have high BH masses (left panel of Figure 10) and have been excluded from our best-fit.¹¹ Excluding the point sources removes candidates with potentially over-massive BHs that are being ubiquitously discovered at high redshifts by Harikane et al. [65]; Maiolino et al. [96]; Greene et al. [57]. However, we find potential under-massive BHs near $9 \leq\log(M_{\star}/M_{\odot}) \leq11$ and $4.4 \leq\log(M_{\mathrm{BH}}/M_{\odot}) \leq6$ (right panel of Figure 11). This scenario is possible if the BH growth is inefficient or delayed in these galaxies. The presence of star formation in low-mass galaxies can stunt the growth of BHs due to supernova feedback heating the gas and inhibiting accretion. This feedback can further hinder the gas from getting to the nuclear regions as a result of the low gravitational potential of these sources (Dubois et al. [45]; Byrne et al. [24]). A possibility for such under-massive BHs in massive galaxies may be associated with recent mergers, where an SMBH is merging with a lower-mass BH (Kelley [77]). In such instances, the lower-mass BH could be undergoing accretion (D’Orazio & Charisi [42]), causing us to observe the broad H$\alpha$ component from its BLR kinematics. The main uncertainty in these sources comes from the intrinsic difficulty in differentiating whether the line widths are due to AGN or outflows (Appendix § 9). However, excluding or including them does not significantly affect our fit.

While the extension of the $M_{\mathrm{BH}}-M_{\star}$ scaling relation to lower galaxy masses may provide clues regarding the origin of BH seeds (Mezcua [104]; Volonteri [150]), it can also point to the variations in the modes of accretion in galaxies (Ricarte & Natarajan [132]). We explore this connection of the scaling relation with theoretical models in § 6.3.

Discussion

In this section, we examine how our main results on the AGN fraction in dwarf galaxies, BH masses, and the $M_{\mathrm{BH}}-M_{\star}$ scaling relation may be impacted by our choices and assumptions. We describe our reasoning for using only the [NII]-BPT diagram in § 6.1. We investigate the higher dwarf AGN fraction observed from DESI compared to SDSS in § 6.2. Finally, we compare our $M_{\mathrm{BH}}-M_{\star}$ scaling relation and its relation to the BH seed formation models in § 6.3.

[SII]-BPT Diagram

The [NII]-BPT diagram is a valuable diagnostic for selecting AGN-dominated sources, but it is metallicity sensitive (Storchi-Bergmann et al. [144]; Carvalho et al. [26]). This may be problematic for dwarf galaxies, which generally have low metallicities. As a result, low-metallicity star-forming and AGN dwarfs can occupy the same region as low-mass starbursts.

Apart from the [NII]-BPT diagram, two other emission-line diagnostics are often used in identifying AGN signatures in galaxies: [SII]/H$\alpha$ vs. [OIII]/H$\beta$ ([SII]-BPT Diagram) and [OI]/H$\alpha$ vs. [OIII]/H$\beta$ ([OI]-BPT Diagram). These two diagnostics reduce the effect of metallicity and separate the star-forming galaxies from AGN-dominated Seyferts and low-ionization nuclear-emitting regions (LINERs). Polimera et al. [123] revised the AGN selection criterion in dwarf galaxies using all three BPT diagrams. They included a new class of dwarf galaxies, the star-forming AGN (SF-AGN) that lie on the star-forming branch in the [NII]-BPT diagram, but are Seyferts/LINERs from the [SII] or [OI]-BPT diagrams. Using this new classification scheme for dwarf galaxies, they find an AGN fraction, ranging from $\approx3\%$ to $\approx16\%$ depending on the catalog used. More recently, Mezcua & Domínguez Sánchez [109] also found a higher AGN fraction $\approx20\%$ when combining all three diagnostic diagrams. However, the validity of the current demarcation lines is unclear and there have been several efforts to improve them (Ji & Yan (2020); Law et al. [89]). In this subsection, we examine the demarcation lines used for the [SII]-BPT diagrams.

We start with the 1,350,112 sources in our sample with robust stellar masses and emission-line measurements in the redshift range of interest (§ 3). We apply the following SNR and AoN cuts on the four emission lines required for the [SII]-BPT diagram:

These cuts result in 449,295 sources with $\log(M_{\star}/M_{\odot}) \geq6.0$. Figure 12 shows the [SII]-BPT diagram of all these sources. The solid black lines in the figure denote the demarcation line for separating star-forming galaxies from the AGN-dominated Seyferts and LINERs (Kewley et al. [79]). The dashed black line separates Seyferts (left of the line) and LINERs (right of the line) (Kewley et al. [80]). These lines are traditionally used in AGN-related studies (Polimera et al. [123]; Salehirad et al. [135]; Mezcua & Domínguez Sánchez [109]). A new demarcation line put forth by Law et al. [89] is shown by the dotted black line. They used line-of-sight velocity distributions of MaNGA sources to constrain the dynamical properties of the ionized gas and provided new demarcation lines for all three BPT diagrams.

From Figure 12, we see that the density of sources in the star-forming branch (dense black region) extends beyond the solid black line. In particular, the density contours on the top-left corner that overlap the star-forming branch extend above the demarcation line by Kewley et al. [79] and are instead more consistent with the Law et al. [89] demarcation line. This finding suggests that the demarcation line needs to be revised for the $[\mathrm{SII}]$-BPT diagnostic.

By considering the Kewley et al. [79] line, we find 55,208/449,295 AGN candidates (AGN Fraction $\approx$ 12.3%). However, if we consider the Law et al. [89] line, we identify 27,953/449,295 AGN candidates (AGN Fraction $\approx$ 6.2%), nearly half of the traditional diagnostic. When we focus on dwarf galaxies, the AGN fraction

changes dramatically from 17.4% with the Kewley et al. [79] line to just 4.6% with the Law et al. [89] line. This stark difference highlights that the AGN fraction is sensitive to the chosen diagnostic line. This could explain the unusually high dwarf AGN fraction in dwarf galaxies reported by Polimera et al. [123] compared to other single-fiber spectroscopic studies. Some galaxies near these demarcation lines may exhibit both star formation and AGN activity. Therefore, careful analysis is essential for accurately estimating the AGN fraction using these diagnostic methods.

To identify AGN signatures in galaxies and facilitate comparisons with similar studies, especially in the dwarf galaxy regime, we utilize the $[\mathrm{NII}]$-BPT diagram in this paper, while deferring the analysis of the other two diagnostics for future publications. Even with just this single diagnostic, we identify a robust sample of dwarf AGN candidates (§ 4.1). Incorporating the other two diagnostics will likely reveal an even higher dwarf AGN fraction than what we observe in this study.

Dwarf AGN Fraction: DESI vs. SDSS

Using the $[\mathrm{NII}]$-BPT diagnostic, Reines et al. [130] reported a dwarf AGN fraction of $\approx$0.5% based on SDSS DR8 spectroscopic data, which utilized the original SDSS spectrograph. Applying the same diagnostic on the early DESI data, our study reveals a dwarf AGN fraction of $\approx$2.1% (§ 4.2), nearly four times higher than this previous estimate. We investigate possible explanations for this increased fraction, such as differences in target selection or instrumental characteristics (e.g., aperture size, spectral resolution). We begin by comparing the targets from both surveys, followed by building a matched sample between SDSS and DESI to isolate the impact of the instrument design.

The DESI targets probe fainter sources than the SDSS targets. Figure 13 compares dwarf AGN candidates from our study (shown as blue circles) with those identified by Reines et al. [130] (shown as pink triangles) in the $r-z$ space. Our candidates reach approximately 2.5–5 magnitudes fainter and extend about 10 times further in redshift compared to the Reines et al. [130] sample.

To facilitate a direct comparison, we create matched samples of galaxies with similar redshift and $r$ magnitudes to compare DESI vs. SDSS spectra. For this purpose, we first recreate the starting sample of Reines et al. [130] using the NASA-Sloan Atlas (NSA). By applying the same selection criteria outlined in their paper, we identify 24,966 dwarf galaxies¹². We cross-match these galaxies with the LS DR9 catalog to establish a common photometric selection, focusing on sources brighter than the SDSS magnitude limit of $r < 17.7$ mag. We obtain 21,468 dwarf galaxies including 123 candidate AGN, resulting in an AGN fraction of $\approx0.6\%$ based on the SDSS spectra, in agreement with Reines et al. [130].

We check the available sample size spanning the same parameter space of $r < 17.7$ mag and $z \le0.055$. The number of DESI targets (4,999 sources) within those limits is smaller than the SDSS targets (21,468 sources). Consequently, we identify up to three potential SDSS sources for each DESI target by locating the nearest neighbors in the $r-z$ space. Some DESI sources have only one closest neighbor, resulting in 4,798 unique DESI sources with a total of 10,478 SDSS neighbors.

To validate our matched sample, we perform the Kolmogorov–Smirnov (K-S test) on $r$ and $z$ distributions of these sub-samples, yielding a $p$-value of 0.99 and 1.0, respectively. These high values indicate that the matched samples have similar distributions in both $r$ and $z$. We identify 189/4,798 dwarf AGN candidates (AGN Fraction $\approx3.9\%$) from DESI and 64/10,478 dwarf AGN candidates (AGN Fraction $\approx0.6\%$) from SDSS. This difference in AGN fraction is even more pronounced, increasing from a factor of 4 to nearly 8 for the matched DESI sample compared to the SDSS-matched

sample. This suggests that the DESI target selection does not account for the observed higher AGN fraction.

We now focus on the differences in instrument properties. The DESI spectrograph offers higher spectral resolution (60–150 $\mathrm{km\,s}^{-1}$) compared to the SDSS spectrograph (100–190 $\mathrm{km\,s}^{-1}$). We note that the SDSS resolution is already sufficient to isolate the emission lines used in the BPT diagram and to resolve the profile of H$\alpha$ in case of broad components. Nevertheless, it is plausible that a higher spectral resolution can lead to a more efficient detection of faint emission lines. Additionally, DESI fibers have a smaller diameter (1.5″) compared to SDSS fibers (3″). Since AGNs are typically located at the centers of galaxies¹³, the smaller fiber size of DESI captures light primarily from this central AGN, thereby reducing contamination from the surrounding galaxy and minimizing dilution from star formation (Trump et al. [146]). We examine these scenarios by comparing the spectra obtained for the same galaxies from the two surveys.

We perform a positional cross-match between the 21,468 dwarf galaxies from SDSS and our starting sample of dwarf galaxies and find 1,725 common candidates. Among these, we find 9 dwarf AGN candidates from Reines et al. [130], while our analysis reveals a total of 24 dwarf AGN candidates. Of the 9 AGN candidates from SDSS, five are also classified as AGN by our selection, two are classified as star-forming, and two do not meet our SNR criteria. In addition to the five common dwarf AGN candidates, we identify 19 additional candidates not selected as AGN by Reines et al. [130]. We visually inspect the DESI spectra and fits of these 19 candidates and confirm that they are valid AGN candidates. Given that the SDSS spectra still detect the emission lines (even if the sources lie on the star-forming branch of the BPT), we conclude that spectral resolution or depth is unlikely to be the primary factor. Instead, we suspect that the smaller fiber size of DESI plays a significant role in detecting a higher number of dwarf AGN candidates compared to SDSS, thereby contributing to the observed increase in the dwarf AGN fraction.

In fact, Moustakas et al. [117] demonstrated that the fraction of AGN in a galaxy sample is strongly influenced by the integrated light collected within the spectroscopic aperture (see Table 5 and Figure 5). More recently, Albán & Wylezalek [3] also found that the number of detected AGN candidates decreases as the aperture size increases. If our supposition is correct, we further expect that a smaller aperture would facilitate the detection of less luminous AGNs by effectively isolating their emission from those of the surrounding star-forming regions.

We estimate the bolometric luminosities ($L_{bol}$) of all dwarf AGN candidates from the matched samples using the [O iii] luminosity, following the formula ($L_{bol} = L_{\mathrm{[OIII]}} \times1000$) from Moran et al. [114]. For the 64 candidates from SDSS, the $\log(L_{bol}/\mathrm{ergs\ s}^{-1})$ values range from 40.0 − 44.2, with a median of 42.2. In contrast, the $\log(L_{bol}/\mathrm{ergs\ s}^{-1})$ values for the 189 DESI dwarf AGN candidates span from 39.3 − 44.0, with a median of 41.6, ≈0.6 dex lower than the SDSS candidates.

Figure 14 shows the distribution of the $\log(L_{bol}/\mathrm{ergs\ s}^{-1})$ values for DESI (in blue) and SDSS (in pink). The lower number of DESI sources at higher luminosities can be attributed to the 1:3 selection ratio of the matched sample. However, as the figure indicates, DESI recovers a greater overall number of AGN candidates compared to SDSS. We perform the K-S test on these luminosities and obtain a $p$-value of ≈ $1.4 \times10^{-5}$, suggesting that the luminosity distributions for the two samples are significantly different. This indicates that DESI effectively probes lower luminosity AGN candidates compared to SDSS within the same $r-z$ range, further supporting our hypothesis that the aperture size is indeed the dominant effect driving the difference in the dwarf AGN fraction.

BH Seed-Formation Mechanisms

The presence of AGN in dwarf galaxies and the Intermediate-Mass Black Holes (IMBHs) that power them hold clues to the origin of SMBHs in the early Universe. In this subsection, we explore our results on the $M_{\mathrm{BH}}-M_{\star}$ scaling relation in the context of available literature regarding BH seed formation models.

Early BH seeds can be broadly divided into two types depending on their initial mass ([104]; [56]): 1) Light seeds ($M_{\mathrm{BH}} \approx10^{2}-10^{3}\ M_{\odot}$) that form as end products of Population III stars in the early Universe (“Pop-III scenario”); 2) Heavy seeds ($M_{\mathrm{BH}} \approx10^{4}-10^{6}\ M_{\odot}$) that form when pristine gas in the early Universe directly collapses into a BH instead of undergoing fragmentation to form stars (“Direct Collapse scenario”). As we cannot directly observe their formation with our current telescopes, we focus on their imprints in the local universe. Most of these early BHs grow via accretion and mergers into the SMBHs we see today. But, a fraction of these might be left over in isolated galaxies and we expect them to be observed as IMBHs in local dwarf galaxies.

Theoretical studies propose that the low-mass end of the $M_{\mathrm{BH}}-\sigma_{\star}$ relation can help differentiate between these BH seed formation models ([151]; [150]). While these studies specifically examine the correlation between the BH mass and stellar velocity dispersion, the expected correlations are likely similar to the relationship between BH mass and stellar mass. If the initial BHs are heavy seeds, the scaling relation is predicted to flatten at the low-mass end. Conversely, if the initial seed BHs are light seeds, the scaling relation from higher-mass galaxies is expected to extend down to lower masses. This is based on the assumption that the observed black holes in local dwarf galaxies represent remnant seed black holes from the early universe, which experienced minimal evolution.

From Figure 11, we observe that our scaling relation extends down to $\log(M_{\star}/M_{\odot}) \approx8.5$, with no evident flattening in the relation. The spread among low-mass galaxies appears similar to that of high-mass galaxies. Given that our candidates reach down to $\log(M_{\mathrm{BH}}/M_{\odot}) \approx4.4$, with most residing below the established relation, these results may suggest a Population-III seed-formation scenario. Despite their uncertain stellar masses, incorporating the point sources from Figure 10 could provide additional insights into this relation. If the stellar masses of these sources are accurate, their distribution is consistent with our confident candidates. Also, if their stellar masses are underestimated, they would shift to the right on the plot, likely aligning with our empirical relation or the one proposed by Suh et al. [145]. Conversely, if their stellar masses are over-estimated, they might represent over-massive BHs in dwarf galaxies. In this case, it could lead to a flattening of the $M_{\rm BH}-M_{\star}$ relation, suggesting a direct collapse seed-formation scenario.

High-redshift observations using JWST have identified $\gtrsim20$ over-massive BHs pointing to a direct collapse scenario [65, 96, 147, 57]. However, these observations are affected by selection biases favoring the brightest AGN and are likely not detecting the faint sources hosting light seed BHs. Moreover, using a new suite of cosmological simulations, Bhowmick et al. [14] suggest that it may be possible to grow these high-redshift over-massive BHs from light seeds with efficient mergers.

Ricarte & Natarajan [132] argued that the low-mass end of the $M_{\rm BH}-\sigma_{\star}$ scaling relation is driven primarily by the accretion modes of BHs and not the seed-formation models. Furthermore, they show that a hybrid seed-formation model with both light and heavy seeds is almost indistinguishable from a light seed model when studying the scaling relations. Regan & Volonteri [127] also emphasized that both seed formation scenarios can occur simultaneously. They suggest that rather than a bimodal distribution of light and heavy seeds, there exists a continuum of BH seeds with the light seeds more abundant and the heavier seeds becoming rarer.

The increased census of dwarf AGN candidates from this study over a spread of BH masses provides an avenue to study the accretion modes in these sources. This number will increase with future DESI releases and will begin the era of the statistical study of dwarf galaxies and their central BHs. Multi-wavelength observations, including future releases of eROSITA X-ray data [124, 134], will be useful in constraining the low-mass end of the galaxy-BH scaling relations and their connection to the BH seed formation models.

Conclusions

Using early DESI spectroscopic data from DESI EDR and 20% of Year 1 (DA0.2), we identify AGN in a sample of 1,385,177 galaxies split into dwarf ($\log(M_{\star}/M_{\odot}) \le9.5$) and high-mass ($\log(M_{\star}/M_{\odot}) > 9.5$) galaxies. We also constrain the $M_{\rm BH}-M_{\star}$ scaling relation to lower galaxy and BH masses compared to previous studies. We use the photometric data from LS DR9 to estimate stellar masses with the SED fitting code CIGALE (§ 3.1, Siudek et al. [142]). We develop a Python-based emission-line fitting code, EmFit, to measure the fluxes and widths of narrow and broad components of various emission lines. By selecting line-emitting galaxies at $z \le0.45$ across the stellar mass range, we conclude the following:

• Using the optical emission-line [NII]-BPT diagnostic, we find 75,928/296,261 ($\approx$25.6%) high-mass AGN candidates and 2,444/114,496 ($\approx$2.1%) dwarf AGN candidates. With these sources, we have more than tripled the existing census of optical dwarf AGN candidates (§ 4).

• DESI has significantly expanded the discovery space for dwarf AGN candidates. We have extended the search to $\approx1.5$ dex lower galaxy masses and up to 10 times higher redshifts compared to previous spectroscopic surveys (Figure 5). The full DESI survey, set to be completed by 2026, will extend the sample size to $\approx10{,}000$ dwarf AGN candidates.

• Our estimation of the dwarf BPT-AGN fraction ($\approx$2.1%) is nearly four times higher than that from a comparable systematic search using SDSS (Reines et al. [130]). This increase can be primarily attributed to the smaller fiber size of DESI compared to SDSS, which aids with the identification of lower luminosity AGN within the same magnitude and redshift range (§ 6.2).

• We find that the BPT-AGN fraction in line emitting galaxies increases with increasing stellar mass, from $\approx$1.2% at $\log(M_{\star}/M_{\odot}) \approx8.0$ to $\approx$93.5% at $\log(M_{\star}/M_{\odot}) \approx11.4$. On average, the BPT-AGN fraction slightly decreases with increasing redshift. However, these trends are affected by selection effects including stellar mass variations, emission line detection limits, and DESI targeting algorithms (§ 4.3).

• Of the 410,757 line-emitting galaxies in our sample, 6,185 ($\approx$1.5%) sources show a broad H$\alpha$ component. Using the [NII]-BPT diagnostics, we have a sample of 4,181 BL-AGN candidates that are likely powered by an AGN (§ 4.1).

• We estimate the BH masses of the BL-AGN candidates using the flux and width measurements of the broad H$\alpha$ component. The BH masses extend down to $\log(M_{\rm BH}/M_{\odot}) \approx4.4$ for confident BL-AGN candidates and to $\log(M_{\rm BH}/M_{\odot}) \approx4.3$ for tentative candidates. Among these, we find 151 confident (and 147 tentative) BL-AGNs have

$M_{\mathrm{BH}} \le10^{6}\,M_{\odot}$, making this the largest sample of IMBH candidates to date (§ 5.1).

• We extend the $M_{\mathrm{BH}} - M_{\star}$ scaling relation down to $\log(M_{\star}/M_{\odot}) \approx8.5$ and $\log(M_{\mathrm{BH}}/M_{\odot}) \approx4.4$ (Figure 11). The empirical fit from these sources has a similar slope to Reines & Volonteri [131], but a slightly higher normalization (Section 5.3).

• The majority of the tentative candidates lie below the empirical fit, suggesting that these sources host under-massive BHs assuming that their broad lines from BLR kinematics. This could be due to the low efficiency of BH growth in galaxies (Section 5.3).

The anticipated increase in the sample of dwarf AGN candidates over the next five years with DESI will accelerate studies of AGN in dwarf galaxies. This expansion will enhance the identification of AGN candidates at the low-mass end of the $M_{\mathrm{BH}} - M_{\star}$ scaling relation, thereby providing insights into BH seed formation. The statistical sample of dwarf AGN candidates will be invaluable for addressing several key questions related to galaxy evolution on the smallest scales, including accretion modes in low-mass galaxies and the co-evolution of galaxies and their central BHs.

DATA AVAILABILITY

All the data from the figures are available in machine-readable form at https://doi.org/10.5281/zenodo.14009453.

ACKNOWLEDGEMENTS

The authors thank the anonymous reviewer for their valuable comments and suggestions that helped improve this paper. R.P. thanks Gurtina Besla, Xiaohui Fan, Rob Kennicutt, and Anil Seth for their valuable comments on this project.

R.P. is currently supported by the University of Utah, and was also previously supported by the University of Arizona and in part by NSF NOIRLab. The research of S.J. and A.D. is supported by the U.S. NSF NOIRLab, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. M.S. acknowledges support by the Polish National Agency for Academic Exchange (Bekker grant BPN/BEK/2021/1/00298/DEC/1), the State Research Agency of the Spanish Ministry of Science and Innovation under the grants ‘Galaxy Evolution with Artificial Intelligence’ (PGC2018-100852-A-I00) and ‘BASALT’ (PID2021-126838NB-I00). This work was partially supported by the European Union’s Horizon 2020 Research and Innovation program under the Maria Sklodowska-Curie grant agreement (No. 754510). M.M. acknowledges support from the Spanish Ministry of Science and Innovation through the project PID2021-124243NBC22. This work was partially supported by the program Unidad de Excelencia María de Maeztu CEX2020-001058-M.

This material is based upon work supported by the U.S. Department of Energy (DOE), Office of Science, Office of High-Energy Physics, under Contract No. DE-AC02-05CH11231, and by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract. Additional support for DESI was provided by the U.S. National Science Foundation (NSF), Division of Astronomical Sciences under Contract No. AST-0950945 to the NSF’s National Optical-Infrared Astronomy Research Laboratory; the Science and Technology Facilities Council of the United Kingdom; the Gordon and Betty Moore Foundation; the Heising-Simons Foundation; the French Alternative Energies and Atomic Energy Commission (CEA); the National Council of Humanities, Science, and Technology of Mexico (CONAHCYT); the Ministry of Science, Innovation, and Universities of Spain (MICIU/AEI/10.13039/501100011033), and by the DESI Member Institutions: https://www.desi.lbl.gov/collaborating-institutions. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the U. S. National Science Foundation, the U. S. Department of Energy, or any of the listed funding agencies.

The authors are honored to be permitted to conduct scientific research on Iolkam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation.

The DESI Legacy Imaging Surveys consist of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS), the Beijing-Arizona Sky Survey (BASS), and the Mayall $z$-band Legacy Survey (MzLS). DECaLS, BASS and MzLS together include data obtained, respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, NSF’s NOIRLab; the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab. NOIRLab is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. Pipeline processing and analyses of the data were supported by NOIRLab and the Lawrence Berkeley National Laboratory (LBNL). Legacy Surveys also uses data products from the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), a project of the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. Legacy Surveys was supported by: the Director, Office of Science, Office of High Energy Physics of the U.S. Department of Energy; the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility; the U.S. National Science Foundation, Division of Astronomical Sciences; the National Astronomical Observatories of China, the Chinese Academy of Sciences and the Chinese National Natural Science Foundation. LBNL is managed by the Regents of the University of California under contract to the U.S. Department of Energy. The complete acknowledgments can be found at https://www.legacysurvey.org/acknowledgment/.

This research uses services and data provided by the Astro Data Lab, which is part of the Community Science and Data Center (CSDC) program at NSF NOIRLab.

Facilities: DESI, Mayall

Software: Astropy [6, 7] (Astropy Collaboration et al. 2022), CIGALE [18], FastSpecFit [116], Matplotlib [68], NumPy [66], QuasarNet [22], Redrock [19, 5].

Double-Peaked Emission Line Sources

When fitting narrow emission lines ($[\mathrm{S\,II}]$ and $[\mathrm{O\,III}]$), we evaluate the presence of additional components as described in § 3.2. These extra components may correspond to an outflow component (see Figure 2) or an additional narrow component. Whenever an extra component is detected in $[\mathrm{S\,II}]$, we apply the same profile to $[\mathrm{N\,II}]$, $H\alpha$, and $H\beta$ fits.

Figure 15 shows an example of a double-peaked emission line candidate that requires two components to fit all the narrow lines: $H\beta$, $[\mathrm{O\,III}]$, $[\mathrm{N\,II}]$, $H\alpha$, and $[\mathrm{S\,II}]$. Such dual-peaks in emission lines can be a result of complex gas kinematics (e.g., double nuclei, accretion disk kinematics, ejecta, bar) in the galaxy (e.g., [99]). However, not all galaxies exhibit a resolved second component across all narrow lines simultaneously. In fact, of the 26,819 galaxies in our sample that show a detectable second component in at least one of the four emission lines used in the $[\mathrm{N\,II}]$-BPT diagram, only 2,803 galaxies have a non-zero flux in the second component across all four emission lines. To ensure accurate AGN classification, it is therefore important to identify these sources and carefully account for their emission line fluxes.

For the double-peaked emission lines, we find that either the amplitudes or the widths of the two components are similar. This is in contrast with the outflows for which the width and amplitude of the two components are distinct. We define the amplitude ratios and $\sigma$ ratios of the two components for each of the emission lines as follows:

where ‘em’ denotes the different emission lines: [SII], [NII], H$\alpha$, and [OIII]. The second component (em; out) is defined such that its width is always greater than the width of the primary component, leading to $\mathrm{SR}_{\mathrm{em}}$ being always greater than 1.

We identify 10,934 sources that show an extra component in [SII] emission line fit, where the same profile is then used to fit [NII], H$\alpha$, and H$\beta$ emission lines. In Figure 16, we show the distribution of these sources in the $\log(\mathrm{AR}) - \log(\mathrm{SR})$ space for [SII], [NII], and H$\alpha$ lines. We see that there are two different distributions of sources in these three panels. The cases where the amplitude of the extra component is greater than that of the primary component ($\log(\mathrm{AR}) > 0$) are sources without an outflow component. We visually inspect the sources where either the amplitude or $\sigma$ values of the two components are nearly equal. We focus on [NII] and H$\alpha$ lines as they are critical for the [NII]-BPT diagnostic. Based on the visual inspection, we find the following optimum cut for separating the double-peaked emission line sources from the outflow candidates:

OR

This criterion is shown as solid red lines in the left three panels of Figure 16. All the sources above these solid red lines (7,789/10,934) are the cases where the extra component is not considered as an outflow component.

We repeat the visual inspection for the [OIII] emission line. The right-most panel of Figure 16 shows the distribution of 18,725 sources that have an extra [OIII] component, in the $\log(\mathrm{AR})-\log(\mathrm{SR})$ space. From visual inspection, we find the following optimum cut for the [OIII] components:

The solid red lines in the [OIII] panel of Figure 16 mark this criterion and the sources that are above or left of these lines (4,247/18,725) are the cases where the extra component is not considered as an outflow component.

We classify these sources selected via the above criteria as double-peaked emission-line candidates, and they have several complexities associated with their emission-line flux measurements. The dual peaks are only marginally resolved and are not consistently present across all emission lines. In several sources, there is also some ambiguity in identifying the primary component of the emission line. For a robust analysis of their ionization signatures, we sum the fluxes of the individual components to derive the total emission-line flux. The uncertainties in these fluxes are estimated by adding the uncertainties of the individual components in quadrature.

To illustrate the importance of this approach, we analyze the [NII]-BPT diagram of the 6,230 galaxies with double peaks in [NII], H$\alpha$, and H$\beta$ emission lines, but no extra component in [OIII]. We compare the results when using only the primary component (left panel of Figure 17) versus when summing both the components (right panel of Figure 17). We find that the star-forming and AGN branches are unclear when using the primary components alone, but are visible when using the sum of the two components. The left panel of Figure 17 results in 3,255 sources ($\approx52\%$) that lie on the AGN-dominated or composite region of the BPT diagram. On the other hand, the right panel results in only 1,288 sources ($\approx21\%$) lying on these regions of the BPT diagram. This suggests that relying solely on the primary components leads to an over-detection of AGN candidates. We, therefore, adopt this procedure for all the emission line-related studies mentioned in this paper (§ 3 and § 4).

!Two-panel [NII]-BPT diagram showing primary components only and the sum of two components bbox=0.13,0.09,0.90,0.49

Visual Inspection of Broad-Line Candidates

Given that we are extending the search for broad components of H$\alpha$ to narrower widths (FWHM (H$\alpha$;b) $>$ 300 km s$^{-1}$), we want to ensure that we are including robust BL-AGN candidates in our analysis. We therefore visually inspect the spectra and fits of 556 BL sources with (FWHM (H$\alpha$;b) $<$ 1000 km s$^{-1}$) and which are BPT-AGN or BPT-composite candidates, and assign them one of the four VI Flags. These are not to be confused with the standard VI Flags used for DESI Survey Validation (Alexander et al. [4]).

• VI Flag = 0: This flag denotes good fits with clear broad components. We consider these as confident BL-AGN candidates. Of the 556 sources, we have 146 sources that satisfy this visual confidence. Figure 18 shows example cases of four such candidates.

• VI Flag = 1: This flag denotes cases where a broad Hα component is detected statistically but is not obvious in visual inspection. Of the 556 sources, we have 154 such cases and we call them tentative BL-AGN candidates and include them in our analysis (§ 5). Figure 19 shows four example candidates with this flag.

!Example emission-line spectra and fits of the [N II] + Hα region of four candidates with VI Flag = 1 bbox=0.10,0.08,0.91,0.21

• VI Flag = 2: This flag is for candidates where the fitting code fails to identify a second or outflow component. This issue arises when the [S II] region is noisy and the fitting does not detect the second or outflow component, but is visible in the [N II]λ6583 line. Consequently, the broad Hα component in these cases may correspond to the missed second or outflow component, rather than being attributed to the kinematics of the BLR. As a result, the estimated BH masses for these candidates are not accurate, and we exclude 252/556 sources with this VI flag. Four examples of such cases are illustrated in Figure 20.

!Example emission-line spectra and fits of the [N II] + Hα region of four candidates with VI Flag = 2 bbox=0.10,0.45,0.91,0.59

• VI Flag = 3: This flag denotes cases with uncertain fits due to complex kinematics in the H$\alpha$ and [NII] region, with only four candidates falling into this category (see Figure 21). All of these sources exhibit strong radio emission, as observed in VLA Sky Survey images available through the LS Sky Viewer¹⁴. While this suggests the presence of an AGN, the broad H$\alpha$ components lead to estimates of under-massive black holes. The kinematics in these cases may be complex, making it challenging for EmFit to accurately separate the emission components. We therefore exclude these candidates from our analysis.

!Example emission-line spectra and fits of the [N II] + Hα region of four candidates with VI Flag = 3 bbox=0.10,0.70,0.91,0.83

In Figure 22, we show the location of these visually inspected BL-AGN candidates on the $M_{\mathrm{BH}}-M_{\star}$ space based on their VI flag. The confident (VI Flag = 0) and tentative candidates (VI Flag = 1) are shown as cyan circles and blue squares, respectively. The sources where the outflow components are possibly selected as broad components (VI Flag = 2) are shown as red crosses. Even though the assumption of BLR kinematics is not valid in such a case, we estimate their BH masses using this extra component. The sources with possible outflows are located in the same region as tentative candidates and are removed from our bestfit estimates of the empirical relation in Section 5.3.

All Authors and Affiliations

Ragadeepika Pucha,^{1, 2} S. Juneau,^{3} Arjun Dey,^{3} M. Siudek,^{4, 5} M. Mezcua,^{6, 4} J. Moustakas,^{7} S. BenZvi,^{8} K. Hainline,^{2} R. Hviding,^{9, 2} Yao-Yuan Mao,^{1} D. M. Alexander,^{10, 11} R. Alfarsy,^{12} C. Circosta,^{13, 14} Wei-Jian Guo,^{15} V. Manwadkar,^{16, 17} P. Martini,^{18, 19, 20} B. A. Weaver,^{3} J. Aguilar,^{21} S. Ahlen,^{22} D. Bianchi,^{23} D. Brooks,^{14} R. Canning,^{12} T. Claybaugh,^{21} K. Dawson,^{24} A. de la Macorra,^{25} Biprateep Dey,^{26, 27} P. Doel,^{14} A. Font-Ribera,^{14, 28} J. E. Forero-Romero,^{29, 30} E. Gaztañaga,^{6, 12, 4} S. Gontcho A Gontcho,^{21} G. Gutierrez,^{31} K. Honscheid,^{18, 32, 20} R. Kehoe,^{33} S. E. Koposov,^{34, 35} A. Lambert,^{21} M. Landriau,^{21} L. Le Guillou,^{36} A. Meisner,^{3} R. Miquel,^{37, 28} F. Prada,^{38} G. Rossi,^{39} E. Sanchez,^{40} D. Schlegel,^{21} M. Schubnell,^{41, 42} H. Seo,^{43} D. Sprayberry,^{3} G. Tarlé,^{42} and H. Zou^{15}

$^{1}$ Department of Physics and Astronomy, University of Utah, 115 South 1400 East, Salt Lake City, UT 84112, USA

$^{2}$ Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA

$^{3}$ NSF NOIRLab, 950 N. Cherry Ave., Tucson, AZ 85719, USA

$^{4}$ Institute of Space Sciences, ICE-CSIC, Campus UAB, Carrer de Can Magrans s/n, 08913 Bellaterra, Barcelona, Spain

$^{5}$ Instituto de Astrofísica de Canarias, C/ Vía Láctea, s/n, E-38205 La Laguna, Tenerife, Spain

$^{6}$ Institut d’Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Spain

$^{7}$ Department of Physics and Astronomy, Siena College, 515 Loudon Road, Loudonville, NY 12211, USA

$^{8}$ Department of Physics & Astronomy, University of Rochester, 206 Bausch and Lomb Hall, Rochester, NY 14627-0171, USA

$^{9}$ Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany

$^{10}$ Centre for Extragalactic Astronomy, Department of Physics, Durham University, South Road, Durham, DH1 3LE, UK

$^{11}$ Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK

$^{12}$ Institute of Cosmology and Gravitation, University of Portsmouth, Dennis Sciama Building, Portsmouth, PO1 3FX, UK

$^{13}$ European Space Agency (ESA), European Space Astronomy Centre (ESAC), Camino bajo del Castillo s/n, Villanueva de la Canada, E28692 Madrid, Spain

$^{14}$ Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, UK

$^{15}$ National Astronomical Observatories, Chinese Academy of Sciences, A20 Datun Rd., Chaoyang District, Beijing, 100012, P.R. China

$^{16}$ Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA 94305, USA

$^{17}$ Kavli Institute for Particle Astrophysics and Cosmology, SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA

$^{18}$ Center for Cosmology and AstroParticle Physics, The Ohio State University, 191 West Woodruff Avenue, Columbus, OH 43210, USA

$^{19}$ Department of Astronomy, The Ohio State University, 4055 McPherson Laboratory, 140 W 18th Avenue, Columbus, OH 43210, USA

$^{20}$ The Ohio State University, Columbus, 43210 OH, USA

$^{21}$ Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA

$^{22}$ Physics Dept., Boston University, 590 Commonwealth Avenue, Boston, MA 02215, USA

$^{23}$ Dipartimento di Fisica “Aldo Pontremoli”, Università degli Studi di Milano, Via Celoria 16, I-20133 Milano, Italy

$^{24}$ Department of Physics and Astronomy, The University of Utah, 115 South 1400 East, Salt Lake City, UT 84112, USA

$^{25}$ Instituto de Física, Universidad Nacional Autónoma de México, Cd. de México C.P. 04510, México

$^{26}$ Department of Astronomy & Astrophysics, University of Toronto, Toronto, ON M5S 3H4, Canada

$^{27}$ Department of Physics & Astronomy and Pittsburgh Particle Physics, Astrophysics, and Cosmology Center (PITT PACC), University of Pittsburgh, 3941 O’Hara Street, Pittsburgh, PA 15260, USA

$^{28}$ Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra Barcelona, Spain

$^{29}$ Departamento de Física, Universidad de los Andes, Cra. 1 No. 18A-10, Edificio Ip, CP 111711, Bogotá, Colombia

$^{30}$ Observatorio Astronómico, Universidad de los Andes, Cra. 1 No. 18A-10, Edificio H, CP 111711 Bogotá, Colombia

$^{31}$ Fermi National Accelerator Laboratory, PO Box 500, Batavia, IL 60510, USA

$^{32}$ Department of Physics, The Ohio State University, 191 West Woodruff Avenue, Columbus, OH 43210, USA

$^{33}$ Department of Physics, Southern Methodist University, 3215 Daniel Avenue, Dallas, TX 75275, USA

$^{34}$ Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK

$^{35}$ Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

$^{36}$ Sorbonne Université, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies (LPNHE), FR-75005 Paris, France

$^{37}$ Institució Catalana de Recerca i Estudis Avançats, Passeig de Lluís Companys, 23, 08010 Barcelona, Spain

$^{38}$ Instituto de Astrofísica de Andalucía (CSIC), Glorieta de la Astronomía, s/n, E-18008 Granada, Spain

$^{39}$ Department of Physics and Astronomy, Sejong University, Seoul, 143-747, Korea

$^{40}$ CIEMAT, Avenida Complutense 40, E-28040 Madrid, Spain

$^{41}$ Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA

$^{42}$ University of Michigan, Ann Arbor, MI 48109, USA

$^{43}$ Department of Physics & Astronomy, Ohio University, Athens, OH 45701, USA