© The Authors 2025

1. Introduction

The presence of supermassive black holes (SMBHs; M_BH > 10⁶ M_⊙) at the center of massive galaxies is well established, yet the origins of these entities remain a topic of intense investigation (Latif & Ferrara 2016; Inayoshi et al. 2020). One plausible pathway involves the formation of intermediate-mass black holes (IMBHs; 100 M_⊙ < M_BH < 10⁶ M_⊙), which serve as seeds that grow through accretion and mergers to become SMBHs (Alexander & Hickox 2012). Understanding the properties of IMBHs, particularly those formed at early epochs, will allow us to understand the evolution of SMBHs (Volonteri 2010). Unfortunately, our technological resources are still incapable of detecting the expected seeds at high redshifts (z > 6 − 12), but a census of the local low-mass black holes can provide alternative constraints. Given the relation observed between the mass of a central massive black hole and the mass of its host galaxy (Reines & Volonteri 2015) in the local universe searching for IMBHs in low stellar mass galaxies (LSMGs) in the near Universe (z < 0.15) is a plausible option to study the origins of SMBHs. An advantage of these LSMGs is that we can study the processes of black holes without the complex merge histories found in more massive galaxies (Volonteri 2010; Greene 2012), this allows us to observe the black hole evolution driven primarily by internal process and compare these findings with those in more massive merger-rich galaxies to determine the effects of different growth mechanisms. Hence, nearby LSMGs are currently the best laboratories for elucidating the initial conditions and growth mechanisms of black holes in the early Universe.

The challenge in detecting IMBHs in LSMGs comes from the faintness and weakness of their observational signatures (Reines & Volonteri 2015; Baldassare et al. 2018). Despite this, the search for active galactic nuclei (AGNs) in LSMGs has produced various results using different approaches. For instance, the first searches were in optical spectroscopic surveys using narrow emission-line (NEL) diagnostic diagrams, known as BPT diagrams (Baldwin et al. 1981), to distinguish between stellar or AGN ionization sources, together with potential detections of broad Balmer emission lines used to estimate black hole masses (e.g., Reines et al. 2013; Moran et al. 2014). However, emission line diagnostics can fail to identify AGNs in these LSMGs, as shown by photoionization models (Cann et al. 2019). For example, Birchall et al. (2020) find that among 61 dwarf galaxies (M_* ≤ 3 × 10⁹ M_⊙) that exhibit AGN X-ray activity, 85% are not classified as AGNs by BPT diagrams. Another approach is the identification of AGN activity using integral field unit spectroscopy observations in dwarf galaxies (Mezcua 2017; Mezcua & Domínguez Sánchez 2024), which is significantly more successful than searches with single-fiber spectra because it allows for the separation of the galaxy emission from the light emitted by the gas ionized by the AGN. Nevertheless, this method is limited by the expensive data required.

Another explored method to search for AGNs in LSMGs is the detection of variability from these sources; stochastic changes in brightness over time can indicate the presence of an AGN due to the dynamic processes occurring near the black hole at the center of the galaxy. For example, Martínez-Palomera et al. (2020) found 502 AGN candidates via nuclear optical variability among 12 300 galaxies with z < 0.35 and Sloan Digital Sky Survey (SDSS) legacy spectra, of which 22 were confirmed as AGN using BPT diagnostic diagrams. In another study, Burke et al. (2022) looked for the expected short variability timescales, constraining the characteristic variability timescale as log(τ/day)≤1.5), in six-year light curves acquired with the Dark Energy Camera. From a parent sample of 63 721 galaxies, they selected 706 AGN candidates including 26 LSMGs (M_* < 10^9.5 M_⊙). However, spectroscopic confirmation was achieved for only one of the five candidates with available spectra. Other authors (e.g., Baldassare et al. 2020; Kimura et al. 2020; Ward et al. 2022) have selected candidates by characterizing the variability amplitude in light curves from different photometric surveys, such as the Palomar Transient Factory (PTF), the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP), and the Zwicky Transient Facility (ZTF). However, confirmation of these candidates through X-ray observations has yielded a low success rate. In some cases, this is attributed to the depth limits of the X-ray observations, while in others, it may be due to the AGN candidate selection method. These cases are discussed in Sect. 6.2.1.

The increasing number of AGN candidates and their confirmation in large numbers is possible by the advantages offered by large surveys. Therefore, it is anticipated that upcoming observations from new large observatories and surveys will significantly advance our understanding of AGNs in faint low stellar mass galaxies. For instance, the Vera C. Rubin Large Synoptic Survey Telescope (LSST; Ivezić et al. 2019) is expected to probe weaker variability despite host contamination. The Extremely Large Telescope (ELT; Gilmozzi & Spyromilio 2007) will offer improved spatial resolution and sensitivity, and the 4-m Multi-Object Spectroscopic Telescope (4MOST; de Jong et al. 2022) will provide enhanced spectral resolution. Together, these instruments will serve as powerful tools in the search and characterization of AGNs in these galaxies.

Here we present the successful identification of AGNs in LSMGs by spectroscopic confirmation of candidates selected through variability together with color and morphology. These candidates were selected by the application of a random forest classifier to ZTF light curves, later limited to objects in the NASA-Sloan Atlas v1.0.1 of nearby galaxies¹ (NSA), with stellar mass M_* < 2 × 10¹⁰ M_⊙ and redshift z < 0.15. For spectroscopic confirmation, we used archival spectra from SDSS-DR17. Furthermore, we present the success in finding X-ray counterparts using the recent eRASSv1.1 catalog (Merloni et al. 2024).

This work is organized as follows. Section 2 describes the selection of AGN candidates; Sect. 3 presents the methods used for the spectroscopic confirmation and characterization of candidates; Sect. 4 provides the results of the characterization. In Sect. 5 we present a comparison with the X-ray catalog, and in Sect. 6 the consistency between indicators of AGN activity, and the comparison with previous similar studies. Our conclusions are summarized in Sect. 7. When needed, we assume a ΛCDM cosmology with parameters Ω_M = 0.3, Ω_λ = 0.7, and H₀ = 100 h km s⁻¹ Mpc⁻¹ with h = 70.

2. Selection of AGN candidates in low-mass galaxies

We selected AGN candidates by means of their optical variability. The variability features were measured from the multi-epoch photometry on a large portion of the sky provided by ZTF. The selections were made using hierarchical random forest (RF) algorithms.

The ZTF is a northern sky survey (Dec > −30) that has been in operation since 2018, with a typical ABmag depth limits of 20.8 and 20.6 in g- and r-band, respectively, and with a cadence of ∼3 days (Bellm et al. 2019). ZTF provides different data products, including an alert stream, data release (DR) light curves, DR images, and DR catalogs, among others. Using these, we made four different selections. The first selection was made using the ALeRCE (Automatic Learning for the Rapid Classification of Events) broker (Förster et al. 2021), specifically utilizing the classifications provided by the ALeRCE light curve classifier described in Sánchez-Sáez et al. (2021). This classifier uses variability features from the ZTF alert stream (i.e., only flux variations detected at greater than 5σ significance in a reference-subtracted image); using both g- and r-band light curves whenever possible, and single-band light curves otherwise; colors from ZTF and WISE, and a morphology score from Tachibana & Miller (2018), where values close to 1 indicate a point source and values close to 0 indicate an extended source. We included objects classified as AGN (host-dominated), QSO (core-dominated), and blazar (jet dominated) in our candidate selection. This sample covers the entire ZTF sky, namely declination larger than − 30 deg. We will refer to this set as the ‘Alerts’ set. For this set, we use data taken between March 2018 and November 2022.

The second and third sets used a similar classifier but trained on the full point spread function photometry (PSF) light curves provided by ZTF in their data releases (Masci et al. 2019). The ZTF DR11 used here contains data from March 2018 to March 2022. These light curves are built from PSF-photometry on the science images of all epochs. The benefit of these data sets is that they include photometry for all the sources present in the reference catalogs, with detection in the ZTF science images. This means that it includes objects of lower variability amplitudes and more data points per light curve than the Alerts. The drawback is that the PSF photometry on these nonreference-subtracted images is not ideal for tracking the flux of a variable point source (i.e., the AGN) superimposed on the extended image of their host galaxies. This combination results in light curves with uncertain errors and modulations caused by seeing and weather conditions. This classification was made for the southern portion of the ZTF sky, namely −30 < Dec < +7, to produce candidates for the Chilean AGN and Galaxy Evolution Survey (ChANGES; Bauer et al. 2023) project of 4MOST (de Jong et al. 2022), avoiding the Galactic plane. The classification was carried out independently for the g- and r-band light curves using data from the ZTF DR11. The classifier is described in Sánchez-Sáez et al. (2023), and our selection includes objects classified as AGNs: lowz-AGN (z ≤ 0.5), midz-AGN (0.5 < z ≤ 3), highz-AGN (z > 3), and blazar (jet dominated). We will refer to these sets as ‘DR-g’ and ‘DR-r’, respectively. Additionally, Sánchez-Sáez et al. (2023) discussed several caveats in the selection of candidates using ZTF DRs (see Sect. 7.1 in Sánchez-Sáez et al. 2023). In particular, for the r-band, there are regions of the sky with exceptionally large densities of epochs, which affect the computation of features and produce an over-density of highz-AGN candidates in the Galactic plane. Moreover, Sánchez-Sáez et al. (2023) also demonstrated that much purer candidate lists are obtained from the ZTF DRs in both g- and r-bands when the candidates are filtered by probability. Therefore, in order to have a clean sample of AGN candidates, we filtered the AGN sample from Sánchez-Sáez et al. (2023) as described below.

For ZTF DR11 g-band:

• Probability of variability pred_init_class_prob > = 0.9 and

• abs(gal_b) > = 20

For ZTF DR11 r-band:

• pred_init_class_prob > = 0.9,

• abs(gal_b) > = 20,

• Significance of the Gaia DR3 proper motion pmsig < = 3σ,

• Gaia proper motion PM < 3,

• IAR_phi > = 0.8

• GP_DRW_sigma > = 0.0001,

• GP_DRW_tau > = 5,

• number of epochs nepochs < = 400.

All the features are explained in Sánchez-Sáez et al. (2023) and Sánchez-Sáez et al. (2021) (most of the definitions are there). These filtering was made to improve the purity of the samples for the selection of targets of the 4MOST ChANGES program and will be described in detail in Bauer et al. (in prep.). In particular, a number of Galactic sources were incorrectly identified as high-redshift AGNs in the original selection, especially in the r-band selection. Applying the previously mentioned restrictions on proper motion, as well as on the structure and timescales of the variations (using IAR_phi and the DRW parameters), eliminated the majority of these misidentifications.

The fourth set is made using forced photometry light curves. However, given the large number of objects (on the order of tens of millions), the ZTF forced photometry service (Masci et al. 2023), which provides light curves on a request basis, was not adequate for our purposes. Therefore, we performed new aperture photometry on the reference-subtracted science images (i.e., difference images), provided by ZTF for all good-quality epochs (i.e., with ZTF metadata infobits = 0, maglimit > 20 mag, seeing < 4 arcsec) in the g-band only, using at most one observation per night, and using an aperture of four arcseconds. The g-band was selected because of the tendency of AGNs to show larger variability amplitude at bluer wavelengths (e.g., MacLeod et al. 2010), and the lower level of contamination from the host galaxy in bluer bands, which makes this band more sensitive to variability. The data used covered the period between March 2018 and September 2022. The photometry was forced on the location of all sources detected in the reference images. This experiment was conducted for the region of the sky with −30 < Dec < +15.5, to cover the ChANGES region. The Galactic plane area was not included due to the large number of contaminating sources in those fields. The light curves obtained from the aperture photometry were built with the purpose of selecting AGN out of the tens of millions of detectable objects in this region of the sky. The variability features were measured from the light curves and we used a similar classifier as the one used for the DR light curves. For this classifier we include some additional features, namely PANSTARRS i − z color, proper motion provided by Gaia, Mexican-Hat filtered variance at timescales of 45 and 450 days, the error on the excess variance and a flux asymmetry estimator ($lc-\mathrm{asymmetry}=\frac{N_p-N_n}{N}$, where N_p and N_n are the number of epochs with flux higher and lower than the mean flux of the light curve, and N is the number of total epochs). The processing of these light curves and the classifier with the additional features is described in Arévalo et al. (in prep.). We refer to this set as ‘Forced Photometry’.

ZTF provides unique object IDs for unique combinations of RA, Dec, filter, field, CCD, and CCD quadrant. Thus, a single object can have multiple IDs in a given band. Since the flux calibration is done on a quadrant basis, there can be small offsets between the mean fluxes of observations with different IDs, which could be mistaken for rapid fluctuations. Thus, for the DR light curves as well as for the Forced Photometry ones, we kept only the light curve (i.e., unique object ID) with the largest number of epochs associated with a single object. This can result in a reduced number of epochs for objects in the overlap regions between pointings. However, it does not adversely reduce the total baseline of the light curve or the sampling rate for the intermediate and long timescale fluctuations.

The sets described above compose our variability-selected AGN samples. To search for massive BHs, we cross-matched these samples to the population of galaxies with low stellar masses retrieved from the NASA-Sloan Atlas v1.0.1² (NSA) catalog. The NSA v1.0.1 contains a catalog of all low-redshift galaxies (z < 0.15) inside the footprint of SDSS DR8, identified in images of this survey. The catalog provides, among many other derived quantities, measurements of the total, K-corrected stellar mass, which we used to select our parent sample. For our work, we selected galaxies by SERSIC-MASS, which is listed in units of M_⊙ h⁻². We adopt the value of h = 0.7 for the conversion into masses in M_⊙. We set a limit to search for AGN candidates in LSMGs at SERSIC − MASS < 10¹⁰ M_⊙ h⁻² (i.e., M_* < 2 × 10¹⁰ M_⊙). The NSA catalog contains 220, 830 objects that meet the mass condition. All these objects lie within the sky region covered by the Alerts set (Dec > −30). The smaller Forced Photometry sky area (−30 < Dec < +15.5) contains 103 054 (47%) of the NSA galaxies, while the DR-g and DR-r sky region (−30 < Dec < +7) contains 73 964 (33%) of the NSA objects. We cross-matched the positions of the four variable AGN samples with the LSMGs using a radius of 1.5 arcseconds. The total number of AGN candidates matched to low-mass galaxies is 383 in the Alerts set, 215 in the Forced Photometry set, 71 in the DR-g set, and 69 in the DR-r set. Since one candidate can be part of one or more sets, counting unique objects the match produced a total of 506 AGN variability-selected candidates.

Figure 1 shows the sky locations of the matched sample. The Alerts set is marked in orange, the Forced Photometry set in blue, the DR-g in green, and DR-r in red. Evidently, most of the sky of the NSA sample is only covered by our Alerts set. The Forced Photometry set covers the area below Dec = +15.5 and the DR-g and DR-r sets below Dec = +7. We note that in the sky area covered by both sets, most of the Alerts candidates also appear in the Forced Photometry set (96 out of 137), while a majority of Forced Photometry candidates are not selected in the Alerts set (120 out of 216 are not in the Alerts set). More strikingly, almost all DR-g candidates are also selected in the Forced Photometry set, but about half of the Forced Photometry candidates are not selected in the DR-g set, and similarly for the DR-r set. These differences happen, in the case of the DR sets, because the DR light curves of extended objects are too noisy to allow the detection of low-amplitude variations, and for the case of the alert light curves, because the 5-sigma threshold prevents the identification of low-amplitude variations with respect to the reference image. The number of sources in each set and the coincidences in the regions of overlap between Forced Photometry and the other sets are summarized in Table 1. Additionally, Fig. 2 shows a Venn diagram for the four sets to visualize the number of elements within each intersection between sets.

Fig. 1. Sky distribution of AGN candidates in low stellar mass galaxies. The different selection sets are marked by different colors, whereas AGN candidates selected in multiple sets have their symbols overlaid.

Fig. 2. Venn diagram showing the number of elements within each intersection between the different sets. The number of objects included only in the Alerts set is higher than in the other sets due to the larger sky region coverage of the Alerts set.

Additionally, we present some properties of the selected objects using the data from the NSA v1.0.1 catalog. Figure 3 shows the distribution of the redshift for each set, which appear similar, with median redshift values of 0.076 for the Alerts, 0.099 for the Forced Photometry, 0.089 for DR-g, and 0.092 for DR-r. As we note in the next section, we find some objects with different redshift to the ones listed in the NSA catalog. However, this number is small and does not affect the redshift distribution and the magnitude and color distributions shown below.

Fig. 3. NSA v1.0.1 catalog redshift distributions for all the objects in each of the four sets: orange for Alerts, blue for Forced Photometry, green for ZTFDR11 g-band, and red for ZTFDR11 r-band.

Figure 4 shows the normalized distribution of the apparent rest-frame g^* and r^* magnitudes using the NSA v1.0.1 Sersic absolute magnitude (these magnitudes are computed using the flux from a Sersic profile modelling in each SDSS bands; we use the ^* symbol to note the rest-frame bands) and redshift. The median g*-band AB magnitudes for the Alerts, Forced Photometry, DR-g and DR-r sets are 19.73, 20.15, 19.72, and 19.86, respectively, while for the r^*-band they are 18.67, 19.09, 18.71, and 18.82 for the same sets. The Forced Photometry set contains objects with higher median g^*-band and r^*-band magnitudes than the other sets. The color g^* − r^* distributions are also shown in Fig. 4, with median ABmag colors of 1.02, 0.88, 0.82, 0.96 for the Alerts, Forced Photometry, DR-g, and DR-r sets.

Fig. 4. Apparent rest-frame magnitude distribution of AGN candidates in low stellar mass galaxies for the g* (top) and r* (middle) bands. The bottom panel shows the color g* − r* distributions. Each set of variability-selected AGN candidates is plotted in different colors: orange for Alerts, blue for Forced Photometry, green for ZTFDR11 g-band, and red for ZTFDR11 r-band.

3. Selection and analysis of spectroscopic data

3.1. Selection of archival SDSS spectra

Spectroscopic characteristics in AGN candidates, like broad emission lines (BELs), are strong confirmation of nuclear activity. As a consequence, AGN candidates selected by different approaches are normally confirmed from their optical spectra (e.g., Sánchez-Sáez et al. 2019; Hviding et al. 2024). Here, we searched for optical spectra of our variability selected sets using the tools available in the SDSS webpage³ using the ZTF RA and Dec coordinates from the Data Release 17 catalog (DR17; Abdurro’uf et al. 2022) with a matching radius of 1.5 arcsec. From the total sample, we found 454 objects with SDSS spectra, of which 450 were of good quality (i.e., SN-median-all > 2). In Sect. 3.2 we show the main differences between the samples with and without spectra.

The NSA parent sample used spectroscopic redshifts from SDSS, NASA Extragalactic Database, the Six-degree Field Galaxy Redshift Survey, the Two-degree Field Galaxy Redshift Survey, and ZCAT (the CfA Redshift Survey), to limit the study to the local universe. However, we found that 16 of our matched AGN candidates have a higher spectroscopic redshift according to the SDSS pipeline and subsequent visual verification, in cases where the redshifts listed in the NSA catalog are from a source different of SDSS. Fifteen of these 16 objects are classified as QSO by SDSS spectral classification. The SDSS classification uses a least-squares minimization performed by the comparison of each spectrum to a full range of templates spanning galaxies, quasars, and stars, all in a range of redshifts (see Bolton et al. 2012 for full description). These objects are listed in Table A.1 and have spectroscopic redshifts ranging from 0.225 to 2.427, whereas the NSA catalog redshifts were all below 0.15.

Additionally, visual inspection showed that a further seven candidates present BELs consistent with an AGN, but lie at higher redshifts than the value reported by the SDSS pipeline. All these higher redshift AGN were classified as such in the SDSS quasar catalog with visual inspections produced by Lyke et al. (2020) and are listed in Table A.2, showing spectroscopic redshifts in the range z = 0.47942–3.650. In all these cases, although the AGN classification was correct, the stellar mass estimate is most likely wrong. Visual inspection of the remaining spectra revealed that six appear like galaxies without emission lines, three have blue, featureless continua consistent with either blazars or stars, and two were not readily identifiable as AGN and could correspond to broad absorption line (BAL) quasars. The spectra of these last five objects are shown in Appendix B. Finally, one spectrum shows no data in the H_α region, flat, featureless continuum, and bad pixels in the H_β region, preventing us from obtaining a useful fit. All the aforementioned 35 candidates with spectra were excluded from further consideration.

Whit this reduction, we are left with 415 AGN candidates at low redshift (z < 0.15) with good spectra, we will refer to these objects as the visually cleaned spectra (VCS) sample. For this sample, spectral fitting was performed following the methodology described below in Sec. 3.3.

3.2. Assessment of biases in the spectroscopic sample

Here, we investigate how the missing candidates differ from those with spectra. We note that of the 506 AGN candidates selected by variability features, only 454 have archival spectra of any quality from SDSS. The top panel in Fig. 5 shows the distribution in the sky of all candidates with and without available SDSS spectra. Some candidates fall outside the main areas of SDSS spectroscopic surveys, in particular targets below Dec = −10. This explains why the DR-g and DR-r sets have lower fractions of candidates with spectra than the other two sets, since they are restricted to regions of the sky with lower declination and therefore also lower SDSS coverage (see Fig. 1). The bottom panel in Fig. 5 shows the distribution of redshifts taken from the NSA catalog, for the candidates with and without spectra, which is the quantity where the samples differ more strongly. Evidently, a larger fraction of the lowest redshift (and also lowest-mass) galaxies do not have an archival spectrum. A consequence of this bias in the availability of spectra is that some of the lowest-mass black holes will not be included in the sample with mass measurements, which limits the ability of our method to confirm low-mass black holes. We note, however, that the missing spectra only represent about 10% of the AGN candidate sample and hence this limitation does not strongly bias the mass distribution of confirmed AGN candidates and other properties reported in the next sections. In any case, the low-mass AGN candidates without spectra should be confirmed or rejected with future spectroscopic campaigns to refine the selection process.

Fig. 5. Comparison of AGN candidates with SDSS spectra (blue) and without spectra (orange). The top panel shows the distribution of these samples in the sky and the bottom panel shows the distribution of their redshifts from the NSA catalog.

3.3. Spectral fitting

For the spectral fitting, we used the pPXF software (Cappellari & Emsellem 2004; Cappellari 2017) with a model that includes the stellar populations from the E-MILES library (Vazdekis et al. 2016), a set of Gaussian profiles for the NELs with a wavelength shorter than 6300 Å and another set for NELs with wavelength longer than 6300 Å, a set of Gaussians with four Gauss-Hermite moments for the BELs, templates for the Balmer higher order emission and Balmer continuum, templates for the FeII pseudo-continuum, and a set of power laws for the accretion disc continuum emission. The power laws are defined as $f_{\lambda }={(\frac{\lambda }{N})}^{\alpha }$, with λ being the wavelength, N the normalization factor, and α the slope of the power law, which can take values in the range −3 ≤ α ≤ 0. We will refer to this combination of components as Model 1. Each component has its own radial velocity and velocity dispersion. For the set of NELs with wavelengths shorter and longer than 6300 Å, the kinematic moments were not tied to allow for deviations in wavelength calibrations. However, the differences between them, as determined from the best-fit model, are small and within the range of the spectral resolution, except for a few specific objects as discussed below. We highlight this here but did not conduct further analysis.

We note that 46 objects show more complex spectra profiles for the broad H_α and H_β, and/or offset wings on the [O III]5007 emission lines. For this reason, a second model was fitted (Model 2) for 39 of the 46 objects. This model adds to Model 1 an extra Gaussian profile for both broad H_α and H_β emission lines, as well as an additional Gaussian profile for the [OIII]5007 doublet, to model possible winds. For the remaining seven objects, Model 2 was used with an adjustment: the H_β narrow emission line Gaussian had its radial velocity and velocity dispersion allowed to vary independently, leading to an improved fit. This modified approach is referred to as Model 3.

We also note that the method implemented by pPXF is sensitive to the initial values of velocity and velocity dispersion. In some cases, the best fit was only achieved when one or both of these parameters were fixed. For 53 objects, when all parameters were left free, visual inspection revealed unsatisfactory fits, particularly in the velocity shifts of the narrow line templates. As a result, for these 53 cases, the most problematic parameter was manually adjusted and fixed for subsequent fitting. The fixed parameters, only for the narrow line templates, were either the line-of-sight velocity (nV) or the velocity dispersion (nVd), or both. The broad line kinematic parameters and the flux normalization for all templates were always free to vary. The specific model used, along with any fixed parameters (i.e., nV and nVd), are detailed for each object in the complete online version of Table 4.

In order to account for the errors on the fitted parameters, we followed the advice in Cappellari & Emsellem (2004) and performed Monte Carlo simulations. We created simulated spectra by combining the best-fit model for each object with different realizations of the observational noise. Below we describe the process we followed. First, we calculated the standard deviation of the residuals from the best-fitting model. In the next step, we used the standard deviation of the residuals to scale a Gaussian deviation to produce random values for each wavelength bin, thereby generating the simulated noise. This provides an empirical measure of the error in the flux. Finally, we combined the best-fit model with the simulated noise and then used the same models, and the same fixed parameters if needed, as explained above to fit this simulated spectra, with the difference that in this case, we included only the stellar templates selected in the best-fitting model. This approach assumes that the stellar template is well-determined in the best-fit model of the observed spectrum. The highest uncertainty is assumed to come from the normalization of the stellar templates in the final fit, and this normalization is always free during the fitting process. This procedure is repeated 100 times for each spectrum. Finally, we computed the 16th and 84th percentiles of the fitted parameters to report the lower and upper estimations of each parameter. Note that use the residuals to generate the simulated noise instead of errors provided by SDSS, which allows us to test the robustness of the fit by evaluating not only the impact of systematic errors but also how specific deviations between the best-fit model and the observed spectrum influence the final results. As the average deviations are slightly larger than the SDSS errors, the use of residuals is a more conservative approach for the estimation of errors on the fitted parameters.

3.4. Estimating black hole masses and Eddington ratios

Estimates of the black hole masses were made using the relation obtained by Mejía-Restrepo et al. (2016),

$$\begin{array}{c}\hfill M_{\mathrm{BH}}=K{(L_{\lambda })}^{\alpha }{\mathrm{FWHM}}^2,\end{array}$$(1)

where the values of K and α depend on combinations of the monochromatic luminosity of the AGN continuum L_λ and full width at half maximum (FWHM) of the BEL used (see Table 7 of Mejía-Restrepo et al. 2016). For L_λ = L₅₁₀₀ and the FWHM of H_β, the values are K = 10^6.864 and α = 0.568. For L_λ = L₅₁₀₀ and the FWHM of H_α, the values are K = 10^6.958 and α = 0.569. The monochromatic luminosity is defined as L_λ = 4πr²F_λ, where F_λ is the monochromatic flux at λ, and r = z × c/H₀, with z being the redshift, c the speed of light, and H₀ = 70 km s⁻¹ Mpc⁻¹. For the Eddington ratio, we used the relation

$$\begin{array}{c}\hfill L_{\mathrm{REdd}}=C_{\mathrm{Bol}}\frac{L_{5100}}{1.5\times {10}^{38}(M_{\mathrm{BH}}/M_{\odot })},\end{array}$$(2)

with C_Bol = 9.26 (see Shen et al. 2008) and M_BH being the mass calculated with the previous equation. To estimate the lower limits in Mass and Eddington ratio, we calculated these values for each Monte Carlo simulated spectra. Then, we obtain the 16th and 84th percentiles, for the lower and upper limits of the mass and Eddington ratio. To report the errors for black hole mass and Eddington ratio we add in quadrature the uncertainty associated with Eqs. (1) and (2). For the mass the uncertainty is 0.19 dex (Mejía-Restrepo et al. 2016) and for the C_Bol value 0.1 dex (Shen et al. 2008).

4. Results

4.1. Confirmation of type I AGN

Type I AGNs are characterized by the presence of BELs in their spectra. We used the criteria of the equivalent width of the broad H_α emission line to be EW_Hα > 5 Å and the signal-to-noise ratio (S/N) of broad H_α flux to be larger than three to consider a BEL detection on the fitting of the spectra. The S/N is defined using the simulation results as the ratio between the 50th percentile (p50) and the difference between the 50th percentile and the 16th percentile (p16), expressed as S/N = $\frac{p50}{p50-p16}$. Additionally, for objects at the low end of the FWHM (FWHM = 940km/s) and high end (FWHM > 9400km/s), we visually inspect the reliability of the measured BEL. As a result, of the 415 AGN candidates in the VCS sample, the fits returned 355 objects that met the EW_Hα and broad H_α flux criteria, of which 323 also had EW_Hβ > 5 Å. Figure 6 shows the distribution of the S/N of H_α flux for all the candidates and for those with EW_Hα > 5 Å. We additionally include two objects, where the spectra show no data in the H_α range, but show EW_Hβ > 5 Å and S/N of broad H_β flux larger than 3. Therefore, the fits returned 357 objects that met the criteria for the presence of BELs, either in H_α or in H_β.

Fig. 6. Distribution of the S/N of the flux of the broad component of Hα. All the variability-selected AGN candidates in low-mass galaxies are plotted in orange; of these, only the ones with EWHα > 5 Å are plotted in blue.

Table 2 shows the number of candidates in each set that had available good spectra and the number of these spectra classified as AGN, either by the presence of BELs described here for the LSMG in column LSMG-EW, or by the higher-redshift QSO classification criteria described in the previous section. We also include a “suspected AGN” category after a further visual inspection of the spectra that were not selected for fitting. These objects have mainly either a featureless blue continuum, which might correspond to a blazar classification, or an uncommon spectrum that might correspond to broad absorption line quasars. The number of suspected AGNs objects are included in Table 2. The spectra of these objects are described in Appendix B and shown in Fig. B.1. Given their uncertainty, we do not include them in the further analysis, and note that their small numbers do not affect the percentage of confirmed candidates.

For 54 candidates in the VCS sample, a BEL was not detected in the spectra above EW_Hα = 5 Å, all of these arise from the Alerts set. In addition, four objects were classified as unreliable detections after visual inspection of the best-fit model, as the resulting EW_Hα > 5 was due to a large FWHM and relatively low peak flux, indicating that the feature was not a true BEL, two of them in the Alerts set and two in the Forced Photometry set. Furthermore, five objects in the Alerts set and two in the DR-r set, which have spectra but are not in the VCS sample (see Sec. 3.1), are also not confirmed or suspected AGN. In total, 65 objects, representing 12.8% of the 506 variability-selected candidates, are not confirmed as AGNs according to their optical spectra.

We visually inspected the light curves and image stamps for the 65 candidates described above, using the ALeRCE ZTF explorer⁴ with the options “Apparent Magnitude” and “Toggle DR” activated, which allows the user to visualize the alert light curve together with the ZTF DR 6 light curve. We found that 15 of the 65 corresponded to a bad subtraction in the ZTF difference images⁵, 30 had a nuclear transient in the ZTF template⁶ two candidates corresponded to bright and seven to weak nuclear transients⁷ including potential tidal disruption events (TDEs)⁸. The last 11 candidates correspond to sources with an unclear nature, but whose alert light curves seem to show real stochastic variations⁹. The footnotes provide one example of each mentioned case. We summarize all these cases in Table 3. The misclassifications due to bad subtractions or transients in the templates are hard to correct from the point of view of a light curve classifier. For the case of nuclear transients, a model that includes the class nuclear transient would solve the issue. Currently, the ALeRCE broker team is working on a new model that includes the class TDE, which will prevent issues like this in the future (F. Förster, priv. comm.).

On the other hand, considering the candidates in the VCS sample as reliable LSMGs, we identified 258 out of 314 objects in the Alerts set with significant BELs. In the Forced Photometry set, 168 out of 170 objects in the VCS exhibit the same characteristic, while all 40 objects in the DR-g and 40 in DR-r subsets display significant BELs. Additionally, high-redshift objects identified through their spectra are also consistent with AGN BELs, with 13 in the Alerts set, 14 in the Forced Photometry set, 11 in the DR-g set, and 1 in the DR-r set.

For all the fitted objects (objects in the VCS sample) the distribution of H_α equivalent widths is shown as a histogram in the top panel in Fig. 7, where we have plotted each variability set separately, even though many objects belong to more than one set. The equivalent widths cluster around 100 − 200 Å, with the Alerts set peaking at a lower EW than the DR-g and DR-r sets and the Forced Photometry set peaking at even lower values. This shift demonstrates the ability of the Forced Photometry selection (and to a lesser extent the Alerts selection) to detect weaker AGN, embedded in stellar continua than the data release light curves. Only the Alerts set shows a tail toward very low equivalent widths, below our threshold of 5 Å. As described above, these AGN candidates are probably produced by spurious variations such as bad image subtractions and transient events appearing in the reference images, which affect the alerts-based classification much more strongly than those based on full light curves.

Fig. 7. Distribution of the equivalent width of the broad component of Hα (top) and of black hole masses derived using Eq. (1) (bottom) for objects with EWHα > 5 Å. In both histograms, each set of variability-selected AGN candidates is plotted in different colors: orange for Alerts, blue for Forced Photometry, green for ZTFDR11 g-band, and red for ZTFDR11 r-band. In the top panel, the vertical line indicates EWHα = 5 Å.

4.2. Black hole masses and Eddington ratio

We estimated black hole masses for all objects with sufficiently significant broad Balmer lines, selecting only galaxies where the EW_Hα > 5 Å and the S/N of broad H_α flux is larger than 3. The fitted widths, equivalent widths, fluxes, and mass estimates of this sample are presented in Table 4.

The bottom panel in Fig. 7 shows the distribution of the black hole masses derived using Eq. (1) for all cases where EW_Hα > 5 Å, except for the object J102530.29+140207.3. This object shows an absorbed (red) spectrum profile so the fitted model only used the stellar populations to fit the continuum and it is therefore not possible to estimate L₅₁₀₀ for the AGN continuum. Similarly to the top panel in this figure, the distribution is shown independently for each variability set. The mean black hole mass for the Forced Photometry set is 1.7 × 10⁷ M_⊙, for the Alerts set is 2.2 × 10⁷ M_⊙, for DR-g is 2.7 × 10⁷ M_⊙, and for DR-r is 2.8 × 10⁷ M_⊙. Almost all black hole masses below 3 × 10⁶ M_⊙ are only found in the Alerts and Forced Photometry sets. To estimate the significance of the differences in the black hole mass distributions of the four sets we performed an Anderson-Darling test¹⁰. Based on the p-values obtained, the DR-r, DR-g, and Alerts sets are all consistent with each other (p-value ≥ 0.17–0.25). The Forced Photometry set, however, has a different distribution compared to any of the other sets (p-value < 0.001). Table 5 summarizes the distribution of black hole masses obtained using Eq. (1) for each set of AGN candidates. The lowest black hole mass obtained corresponds to the well-known low-mass AGN NGC 4395 included in the Alerts set, and our result of $M_{\mathrm{BH}}=3.9_{-1.2}^{+0.9}\times {10}^4{M}_{\odot }$ is in agreement with the reverberation mapping mass M_BH = (4.9 ± 2.6)×10⁴ M_⊙ estimated by Edri et al. (2012). We note that the masses derived using the FWHM of H_β in Eq. (1) are systematically offset and 30% lower compared to those derived using the FWHM of H_α. Since they are measured using the same continuum flux and the fitted width of H_α and H_β are the same, the difference is only set by the parameters used in Eq. (1).

We estimated the luminosity at 5100 Å from the measured flux from the best-fit power-law continuum component in erg s⁻¹, and a bolometric correction factor of C_Bol = 9.26 as mentioned in Sec. 3.4. Combining these luminosities with the previously estimated mass, we calculated the Eddington ratio R_Edd, and show them in the same table (Table 4). The distribution of the Eddington ratio is shown in Fig. 8. We found that, unlike the mass distribution, the Eddington ratio for all sets is similarly distributed. The mean Eddington ratio in each set are 7.1 × 10⁻² for the Forced Photometry, 5.7 × 10⁻² for the Alerts, 6.4 × 10⁻² for DR-g and 6.5 × 10⁻² for DR-r. Performing an Anderson-Darling test to estimate the difference between the Eddington ratio distributions, we find a p-value of 0.23 for a simultaneous comparison of the four sets, confirming their consistency.

Fig. 8. Distribution of the Eddington ratios derived for each set of AGN candidates. Each set is plotted in different colors: orange for Alerts, blue for Forced Photometry, green for ZTFDR11 g-band, and red for ZTFDR11 r-band.

4.3. Galaxy classification by narrow emission-line ratios

The spectral fitting described above also produced fluxes for the narrow permitted and forbidden emission lines. The relative fluxes of these lines can differentiate between different excitation mechanisms, such as star formation (SF), Seyfert activity, and LINERs (i.e., Baldwin, Phillips & Terlevich (BPT) diagnostics, Baldwin et al. 1981). We state whether each object falls in the SF, Seyfert, or LINER portions of each BPT diagram, following the relations in Kewley et al. (2001), Kauffmann et al. (2003a), Kewley et al. (2006) and Schawinski et al. (2007).

Figure 9 shows the BPT diagrams for all the objects where the emission lines were measured. The color bar indicates the black hole mass on a logarithmic scale. Objects with no BH mass estimate, mainly because they lacked BELs, are plotted with black triangles. In the top left of each panel we show the average error bars of both line ratios. These are estimated from the percentiles of the distributions of line ratios for each object resulting from the Monte Carlo simulations described in Sec. 3.3. In these diagrams objects with lower black hole mass are preferably found in the star-forming region but the distribution of masses largely overlaps. For example, the mean, median, and standard deviation of the mass in logarithm scale and solar masses (log(M/M_⊙)) for objects in the star-forming class are 6.83, 6.90, and 0.49 respectively; for the composite class are 7.01, 7.05, and 0.50; and for the Seyfert class are 7.08, 7.09, and 0.48. From the 415 objects in the VCS sample, 20 have no BPT classification on the [OIII]/H_β – [NII]/H_α diagram: for three of them, the spectrum has no data in the H_α region (note that these objects are different from the one excluded in Sect. 3.1 because of no data in the H_α region, as we were able to fit the H_β region for the three mentioned here, and in two cases, we found a BEL); in three cases the best-fit model does not fit a narrow H_β line; and in 14 cases the best-fit model could not separate the overlapped [NII] and H_α narrow emission lines. All these 14 objects are classified as QSO by the SDSS pipeline and all have a significant broad H_α detection.

Fig. 9. BPT diagrams from the measured flux of narrow emission lines. The color bar indicates the black hole mass on a logarithmic scale. The black triangles represent the objects with no mass estimation. In the top left corner, the black lines represent the mean errors. There are four objects classified as Seyfert (LINER) in all three BPT diagrams with no mass estimation (green diamonds). For these objects the WHAN classification is presented in Table 6. Furthermore, the object represented as a green square has an X-ray counterpart in the eRASSv1.1 catalog.

In Fig. 9, the objects represented as green diamonds, J124438.47+061804.6, J120141.43+382821.5, J131305.85+232733.6, J215055.73−010654.1, and a square (differentiated because this object has an X-ray counterpart; see Secs. 5 and 6.1), J121736.78+293628.8, fall on the AGN region. All of them show an EW of broad H_α and H_β lower than 5 Å, which we classify as having no BELs. For these objects, we used the equivalent width of the narrow H_α emission line versus the [NII]/H_α (WHAN) criteria introduced by Cid Fernandes et al. (2011). The WHAN diagnostic helps to distinguish between LINER spectra excited by an AGN and similar spectra produced instead by old stellar populations. Table 6 shows the classification according to the WHAN diagram. As shown in Fig. 10 the five objects fall on the AGN classification.

Fig. 10. WHAN diagram for the five objects resulting in the AGN region on BPT diagrams, but showing no BELs. The blue markers show the position of these objects in the diagram. Three of them fall in the S-AGN region and two in the W-AGN region.

There are 121 objects classified as star-forming by the BPT diagnostic diagram [OIII]/H_β versus [N II]/H_α. From these, 75 have a significant EW of the broad H_α emission line, with 55 objects in the Alerts set, 37 in the Forced Photometry, 13 in the DR-g, and nine in the DR-r set. We note again that one object can be in one or more sets. Furthermore, 62 of the 70 objects that were classified as composite also show significant EW of broad H_α. Of these 62, 40 are in the Alerts set, 30 in the Forced Photometry set, four in the DR-g and three in the DR-r set.

We used the publicly available measurements provided by SDSS to compare our BPT classification with those obtained by other groups. We utilized the catalog from the MPA-JHU group, whose method for spectral fitting and measurements is described on the SDSS webpage¹¹ and based on Brinchmann et al. (2004), Kauffmann et al. (2003b) and Tremonti et al. (2004). We found a total of 410 objects after cross-matching the VCS sample with the MPA-JHU catalog, using the coordinates provided in the NSA catalog, with a radius of 1.5 arcsec. In addition, we used the Portsmouth Stellar Kinematics and Emission Line Fluxes Value Added Catalog from SDSS for the same comparison. The description and methodology are also in the SDSS webpage¹² and in Thomas et al. (2013). After a cross-match between the Portsmouth catalog and our VCS sample, we found 116 objects, from which 17 are also in the MPA-JHU catalog. In Table 7 we show the comparison of the different BPT classes between this work and the MPA-JHU and Portsmouth catalogs. We note that our results for BPT classification are more consistent with the results obtained by the Portsmouth group. This could be a consequence of the differences in the models used to fit the spectra. In the case of the MPA-JHU method, the goal had been to separate star-forming galaxies from AGNs, so the model does not distinguish broad from narrow emission lines. As a consequence, the fluxes of the narrow H_α and H_β lines are overestimated in the spectra that contain a broad component. This explains why many objects classified as Seyfert and composite in our work appear as star-forming in the MPA-JHU catalog, while almost all the star-forming classifications in our work are also classified as star-forming in the MPA-JHU catalog (see Table 7). We explore these flux differences due to the differences in model fitting in Appendix C.

4.4. Contribution of the AGN component to the continuum

The continuum emission in the spectra of low-redshift (z < 0.15) AGNs at optical wavelengths is mainly composed of the stellar light of the host galaxy and the continuum from the accretion disk; the relative contribution of these two components is dependent on the properties of the host galaxy and the AGN. For objects where the EW_Hα > 5 Å, we measured the relative contribution of the AGN component at 5100 Å as the ratio of the monochromatic flux of the AGN to the total continuum flux (f_5100, AGN/f_{5100, cont.}). To compare the emissivity profile of the AGN component obtained from pPXF, we fit this component as one power law using the model $f_{\lambda }={(\frac{\lambda }{N})}^{\gamma }$, and in Fig. 11 we compare the slope γ and the fractional contribution of the AGN. For bluer (meaning more negative values of γ) AGN continuum, the relative contribution of the AGN to the total continuum shows a tendency to increase.

Fig. 11. Here we present the AGN-continuum slope and the AGN contribution to the continuum. Center: comparison of the AGN-continuum slope from the best fit with the ratio of the AGN component and total continuum at 5100 Å. The different colors represent the different sets: orange for Alerts, blue for Forced Photometry, green for DR-g, and red for DR-r. The black square and error bars represent the typical error in the AGN-continuum slope constrained by the 16th and 84th percentile values from the simulations. In the top panel, the distribution of the AGN relative contribution (f5100, AGN/f5100, cont.) is plotted for each set. Additionally, in the right panel, the power-law slope (γ) distribution of each set is shown.

There are 108 objects with −0.1 < γ ≤ 0.0. In most cases, this flat AGN component is expected due to the limit in the α parameter in the model introduced in Sect. 3.3 ($f_{\lambda }={(\frac{\lambda }{N})}^{\alpha }$, with −3 ≤ α ≤ 0) when fitting the spectra that show flat profiles, possibly dominated by the host galaxy. However, we identified 10 objects with γ = 0 that exhibit spectra dominated by the AGN component, with 0.51 < f_5100, AGN/f_{5100, cont.} < 0.76; the spectra of these objects show flat profiles and do not display very strong stellar absorption features. For these cases, the normalization of the stellar component and the AGN component could be degenerate (see spectral decompositions in Appendix D). However, note that the black hole mass measured using the AGN continuum (Eq. (1)) and using the H_α luminosity are in agreement as shown in Sec. 4.5. In addition, in the top panel of Fig. 11 we show the histograms of the AGN relative contribution of the different sample sets. The distributions show that Alerts and Forced Photometry sets find more objects with lower (< 0.4) relative AGN contribution than the DR. In the right panel of Fig. 11 we show the distribution of the power-law slope, where the Alerts and Forced Photometry sets tend to find more objects with redder AGN continua or host dominated.

When comparing the distribution of the AGN relative contribution and power-law slope γ for different BPT classes, namely star-forming, composite, and Seyfert/LINER (Fig. 12), we do not see strong distinctions; this means that objects classified as star-forming by the BPT diagrams show characteristics of their continuum emission similar to those classified as AGN. The difference between star-forming and AGN classes comes from the strength of the narrow emission lines. As shown by Mezcua (2017), and in more recent work by Mezcua & Domínguez Sánchez (2024), using IFU data, dwarf galaxies can have narrow-line characteristics of AGN in small spatial regions, while the rest of the galaxy is dominated by star formation. In our case, it is possible that the broad-line AGN with star-forming classification are similarly affected by contamination from the rest of the galaxy within the fiber aperture. This contamination appears to influence the line emission and continuum differently, potentially depending on the level of star formation activity, which alters the narrow-line fluxes relative to the stellar continuum. Consequently, we do not find a correlation between BPT type and AGN continuum dominance.

Fig. 12. Here we present the distribution of the AGN contribution to the continuum and the AGN-continuum slope for the different BPT classes. Top: distribution of the AGN relative contribution (f5100, AGN/f5100, cont.). For different BPT classes, according to the [OIII]/Hβ vs. [N II]/Hα diagram. The distributions for different classes, star-forming (purple), composite (pink), and Seyfert/LINER (brown) are similar to the Seyfert/LINER objects concentrated between 0.3 < f5100, AGN/f5100, cont. < 0.75. Bottom: Distribution of the fitted power slope for the same BPT classes.

4.5. Black hole versus stellar mass

In the top panel in Fig. 13 we show a histogram of the ratio between our estimate of the black hole mass and the stellar mass of its host galaxy given in the NSA catalog. For all samples, the ratio is distributed mainly between values of M_*/M_BH = 300 and 3000, with the DR sets distributed almost uniformly in this range and the Alerts and Forced Photometry sets peaking around M_*/M_BH = 1000. The median value of the ratios M_*/M_BH for the Alerts, Forced Photometry, DR-g, and DR-r sets are 906, 1252, 592, and 787, respectively. An Anderson-Darling test results in a p-value lower than 0.05 for the comparison between the Forced Photometry set with the other sets. While for the comparison between any pair of the Alerts, DR-g, and DR-r sets the p-value is always higher than 0.05. This shows that the mass ratio distribution of the Forced Photometry set is significantly different from the other sets.

Fig. 13. Stellar mass and black hole mass comparison. Top: distribution of the ratio of our estimates of the black hole mass to the stellar mass of the host galaxies listed in the NSA catalog, separated by selection set. Bottom: black hole mass estimated using Eq. (1) vs. the stellar mass of the host galaxy. The results from Reines & Volonteri (2015) are shown in gray for comparison. The purple line represents the MBH = 0.001M* relation (Kormendy & Ho 2013), while the black dashed line indicates the relation obtained in Reines & Volonteri (2015). The mean errors are shown with error bars in the bottom right corner.

Reines & Volonteri (2015) studied the relation between black hole mass and total stellar mass for a sample of 341 galaxies with stellar masses in the range 10⁸ − 10¹² M_⊙, although only about 40 of them with masses below 10¹⁰ M_⊙. They find a linear relation between black hole mass and total stellar mass with a ratio M_*/M_BH ∼ 4000 (i.e., with black hole masses four times less massive than those found here). In the bottom panel in Fig. 13 we show the relation between black hole mass and total stellar mass in our sample, together with the best-fitting line obtained by Reines & Volonteri (2015) of M_BH = 0.00025M_*. Although some of our sources reach this level, the majority of our sample lies above. We also plot for reference the relation M_BH = 0.001M_* (Kormendy & Ho 2013), found for generally more massive elliptical galaxies, which is closer to our data.

One key difference to examine is how black hole mass is estimated. In the Eq. (1) in this work we use a single-epoch mass estimation from Mejía-Restrepo et al. (2016), that uses the FWHM of the BELs and the AGN luminosity at 5100 Å. While in Reines & Volonteri (2015) (Eq. (1)) the mass is estimated, also using a single-epoch, but considering only the H_α properties (FWHM, Luminosity). Figure 14 shows the comparison of the masses estimated by the two methods for our sample, which give very consistent results, with only a slight discrepancy towards higher masses at the level of ∼0.3 dex. As a consequence, this cannot account for the difference in the M_*/M_BH.

Fig. 14. Comparison of the black hole mass computed by Eq. (1) in this work (i.e., using the FWHM of BELs and AGN luminosity at 5100 Å) vs. the mass computed using Eq. (1) of Reines & Volonteri (2015) (i.e., using the Hα FWHM and its luminosity). The purple line shows the 1:1 relation in logarithmic scale. The dashed purple lines indicate the 0.5 dex uncertainty from Eq. (1) of Reines & Volonteri (2015). The mean errors are shown in the top left with error bars. Since both mass estimates use the FWHM of Hα, the errors are partly correlated

Another key difference to investigate is the adopted samples. Although the parent sample in Reines & Volonteri (2015) and this work is the same (i.e., the NSA catalog of local galaxies), the criteria for the selection of AGN differs. In our work, we select 506 candidates based on the optical variability features, of which 357 have a significant broad Balmer line detection (i.e., a flux S/N > 3 and EW ≥ 5 Å in broad H_α or H_β). In Reines & Volonteri (2015), 341 AGN candidates were selected based on the simultaneous identification of broad Balmer lines and a Seyfert or AGN classification in all three BPT diagnostic diagrams using the optical spectra; by contrast, a broad line was considered detected in their work if it led to a 50% improvement in the χ² of the fitted model. As shown in Secs. 4.3 and later in Sect. 6.1, 169 of our variability-selected AGN are not classified as AGN in all three BPT diagrams in Fig. 9, despite having significant broad Balmer lines and in many cases an X-ray detection as well. The 222 galaxies classified as star-forming or composite (including mixed BPT classifications) in our sample appear in all stellar mass bins, but are more frequent for lower masses and comprise almost all the AGN candidates in galaxies with stellar masses below 10⁹ M_⊙. Therefore, the removal of candidates with non-AGN BPT classification in Reines & Volonteri (2015) at least partially explains why our sample has many more galaxies with stellar masses below 10¹⁰ M_⊙ (i.e., 161 AGN in galaxies with stellar mass below 10¹⁰ M_⊙ in our work versus 41 in Reines & Volonteri 2015, with 15 objects in common in this mass range).

The difference in the confirmation of AGN criteria, however, does not explain why our objects have, on average, larger black hole masses for a given stellar mass. The difference is probably produced by the limitations of the variability selection: although it can find AGN in lower mass galaxies, it preferentially detects the galaxies with the largest black holes. This supposition is supported by the median of the mass ratios detected by the different sets, where the more sensitive sets (i.e., Alerts and Forced Photometry) detect AGN up to higher ratios of stellar mass to black hole masses than the DR sets (see the top panel in Fig. 13). This limitation might be reduced with even more sensitive, higher S/N light curves that will be produced in the future by the Vera C. Rubin LSST (Ivezić et al. 2019).

Finally, we note that the stellar masses calculated via SED fitting, such as those in the NSA catalog, can be overestimated in active galaxies, as demonstrated in Buchner et al. (2024). These authors show that the bias factor is a function of AGN bolometric luminosity and stellar mass. For the majority of the bolometric luminosities (i.e., log L_bol = 41 − 44) and stellar masses (i.e., log M_* = 9 − 10) estimated for our sample, the mean and standard deviation values of this bias are in the range 0.1 ± 0.2 to 0.2 ± 0.4, which should not affect considerably our results.

5. Frequency of X-ray counterparts of the variability-selected AGN in low-mass galaxies

AGN activity often results in the emission of X-rays and conversely, sufficiently high X-ray luminosity in galaxies is often considered proof of black hole activity (Pounds et al. 1994). We explored the incidence of X-ray detections in our sample of variability-selected AGN candidates in the interest of establishing whether X-ray emission is also ubiquitous in low-mass galaxies and/or low mass SMBHs, in particular those selected by optical variability. Arcodia et al. (2024) carried out a similar analysis on several literature samples, which we discuss in Sec. 6.2.1.

We used the recent publication of the first data release (DR1) of the X-ray source catalog of the SRG/eROSITA all-sky survey from the German Consortium (Merloni et al. 2024), to search for soft X-ray counterparts of our variability-selected AGN sets. We will refer to the soft X-ray counterparts simply as X-ray counterparts in the following text.

The eRASSv1.1 catalog covers half the sky, specifically the Galactic longitudes in the range 180–360 deg, to a depth of about 5 × 10⁻¹⁴ erg s⁻¹ cm⁻² in the 0.2–2.3 keV band. Of our 506 AGN candidates, 230 fall in this region of the sky. Performing a simple cross-match by sky coordinates between the two samples with a matching radius of 5 (10) ″ results in 123 (150) AGN candidates with an X-ray counterpart (i.e., 53(65)% of the candidates have a match). We note that the source density of the eRASSv1.1 catalog in this area of the sky is about 30 points per deg² or equivalently 0.00073 sources in a circle of 10″ radius. Therefore, even using the larger cross-matching radius the probability of chance alignments is small. Increasing the matching radius to 15″ only results in an addition of five matches, from which we conclude that the counterparts are most likely related to the AGN candidates and not to chance alignments.

The discussion above includes all variability-selected AGN in the low stellar mass sample of the NSA catalog. As discussed in Sec. 3.1, some of these galaxies are indeed QSOs but have higher redshifts than expected for the NSA sample, so their stellar masses are incorrect. Therefore, their spectra were not modeled. If we consider only the 415 VCS sample of galaxies with good S/N spectra, that are indeed LSMGs at z < 0.15 and show at least narrow emission lines, we find 195 in the portion of the sky included in the eRASSv1.1 catalog. Of these, 130 have an X-ray counterpart within 10″, or 67% of the sample, a similar fraction as in the total variability-selected AGN sample in the eROSITA-DE sky.

Table 8 shows the number of AGN candidates in the eROSITA-DE sky in each set, and the number of those candidates with an X-ray counterpart in the eRASSv1 catalog. We note that the ratio of matched candidates in each individual set is larger than the ratio of matches for the combined sets in this area of the sky (150/230). This counter-intuitive result is caused by the larger probability of X-ray matches for sources that appear in multiple sets, while sources that appear in only one set have a lower ratio of matches. The total number of unmatched candidates in the combined set is 80, while the sum of the unmatched candidates in each individual set is 100, showing that most of the unmatched candidates only appear in one set while only a few appear in two or more sets.

In Fig. 15 we show the distribution of equivalent widths of broad H_α obtained for all the AGN candidates in bona fide LSMGs that are in the eRASSv1.1 sky (orange), the ones with X-ray counterparts (blue) and, of these, the ones where the BPT diagnostic [OIII]/H_β – [N II]/H_α returns a star-forming type (red). We added 0.6 to objects with EW_Hα = 0Å (12 in total), so they can be shown in the plot with a logarithmic scale. There are 56 AGN candidates with fitted spectra in this region of the sky which have a star-forming classification according to the [O III]/H_β – [NII]/H_α BPT diagnostic. From these, 27 have an X-ray counterpart and all have equivalent widths of H_α > 5 Å.

Fig. 15. Distribution of Hα equivalent widths for the AGN candidates with fitted spectra in the eROSITA-DE sky (gray), of these the ones with X-ray counterparts (blue), and of these the ones classified as star-forming by the BPT diagnostics (red). Objects with EWHα = 0 Å, 12 in total, are counted using an added EW value of 0.6.

We note that, only one of the 25 objects with an equivalent width EW_Hα < 5 Å in the eRASS sky has an X-ray counterpart. The object J121736.78+293628.8, which has an X-ray counterpart, shows a spectrum profile dominated by the stellar light and is classified as Seyfert in the [OIII]/H_β – [N II]/H_α BPT diagram and it is also classified as W-AGN in the WHAN diagram. From these numbers, we conclude that the majority of AGN candidates with EW_Hα > 5 Å have X-ray counterparts, even if they are classified as star-forming by the BPT diagnostics.

The relationship between X-ray emission and H_α luminosity has been studied by various authors (e.g., Panessa et al. 2006; Ho et al. 2001; Shi et al. 2010). Panessa et al. (2006) demonstrates the relationship between X-ray luminosity in the 2–10 keV band and H_α luminosity for Type 1 Seyferts, Type 2 Seyferts, mixed Seyferts, Compton-thick candidates, and low-redshift quasars. This relationship is consistent with the findings of Ho et al. (2001), where the best-fit equation is given by log(L_X) = (1.11 ± 0.054)log L_Hα − (3.50 ± 2.27). For comparison, in Fig. 16, we plot the X-ray luminosity in the 2–8 keV band in the observed frame of reference from the eRASSv1.1 catalog and the H_α luminosity measured from our best-fit models, including both the broad and narrow components, for objects with EW_Hα > 5 Å. We find that the data are well distributed around the relationship from Ho et al. (2001) (purple line in the plot). To contrast, in Fig. 16, we also plot the relation between X-ray luminosity in the 2–10 keV band and H_α luminosity for low-mass star-forming galaxies, given by log(L_X) = log L_Hα − 1.40 as presented in Rosa González et al. (2009) (red line in the plot).

Fig. 16. X-ray luminosity (2–8 keV band) and Hα (narrow+broad) luminosity comparison. The purple line plots the relation found by Ho et al. (2001) for AGN host galaxies, while the red line plots the relation for star-forming galaxies presented in Rosa González et al. (2009). The mean errors are shown in the bottom right with error bars.

As expected, objects without X-ray detection in eRASSv1.1 have a lower value of H_α (narrow + broad) flux. For these objects, the median is Flux_Hα = 988 erg cm⁻² s⁻¹, while for those with an X-ray counterpart, it is Flux_Hα = 2395 erg cm⁻² s⁻¹. We intend to conduct further comparisons and analyses of X-ray properties in future work, as this is beyond the scope of this paper.

6. Discussion

6.1. Consistency between different indicators of activity

We consider the NEL diagnostics, the presence of broad permitted lines, and the detection of an X-ray counterpart as indicators of black hole activity. In this section we quantify the candidates that fulfill one or more of the criteria, limiting the discussion to the visually cleaned 415 objects with good-quality optical spectra that are indeed low-redshift low stellar mass galaxies, and that show at least narrow emission lines.

Regarding classification based on NEL ratios, we note that four candidates without BELs are classified as Seyfert/AGN or LINER types in all three BPT diagrams. One additional candidate without BELs has the same classification in two BPT diagrams but lacks classification in the [OIII]/H_β – [OI]/H_α diagram due to the nondetection of the [OI] line. Of these, three are classified as S-AGN and two as W-AGN in the WHAN diagram, and hence consistent with black hole activity, of either Seyfert II or LINER type. Three of the five are in the eROSITA-DE sky and one has an X-ray counterpart in eRASSv1.1. The BPT and WHAN classification of the object with X-ray counterpart is highlighted by using a different (square) marker in Figs. 9 and 10. The other two of these galaxies were re-observed after they showed alerts in ZTF and therefore became candidates of changing state AGN (CSAGN): SDSS J215055.73−010654.1/ZTF18abtizze by López-Navas et al. (2022) and J120141.43+382821.5/ZTF18aaqjyon by López-Navas et al. (2023), but neither had developed BELs. Therefore, they represent either true type II AGN or transient flaring events in AGN that have switched off. Apart from these five galaxies, all objects classified as AGN/Seyfert by all three BPT diagnostics have significant BELs.

Of the 415 galaxies with fitted spectra, 185 have an AGN/Seyfert, composite, or LINER classification in all three BPT diagrams. Of these, 181 also have significant BELs. Of the 185 galaxies with consistent AGN classification, 86 broad-line and two narrow-line AGN fall in the eROSITA-DE sky. Of the broad-line objects, 68 have an X-ray counterpart (i.e., 68/86 = 79%) while only one of the two narrow-line objects has an X-ray counterpart. We note here that the spectrum of the source with an X-ray counterpart was taken in 2006 and the X-ray observation published in the eRASSv1.1 started on December 13, 2019 (Predehl et al. 2021). Thus, considering the existence of CSAGNs (e.g., López-Navas et al. 2023), new spectroscopic observations can help to unveil the presence or absence of BELs together with the X-ray emission.

At the other extreme, there are 104 galaxies classified as SF in all three BPT diagrams. Of these, 46 do not have significant broad H_α. Of the 58 galaxies with consistent SF classification that do show BELs, 29 are in the eROSITA-DE sky and 23 have X-ray counterparts (i.e., 23/29 = 79%), while none of the ones without BELs have X-ray counterparts. For completeness, we note that there are 118 galaxies with mixed (meaning different classification in one or two BPT diagrams, or no classification in the [O III]/H_β – [NII]/H_α diagram) BPT classifications, and of these 111 have BELs.

Therefore, the likelihood of finding broad lines is much higher for consistent AGN/Seyfert/composite classifications (181/185 = 98%) than for consistent SF classifications (58/104 = 56%), while mixed types also have a higher likelihood (111/118 = 94%). The probability of finding X-ray counterparts is much higher for galaxies that show BELs: of 169 galaxies with significant BELs in the eROSITA-DE sky 129 have X-ray counterparts (129/169 = 76%), while for the 26 galaxies in this region of the sky that do not have BELs only one has an X-ray counterpart (1/26 = 4%).

Finally, the probability of finding X-ray counterparts in SF versus AGN/Seyfert/composite BPT types only depends on whether they have BELs or not. For galaxies with BELs, the fraction of X-ray matches is 79% for both (galaxies with consistent SF or with consistent AGN/Seyfert/composite classifications), while almost all galaxies without broad lines have no X-ray counterparts, evidently regardless of where they appear in the BPTs.

In summary, Fig. 17 presents a Venn diagram illustrating the distribution of 195 VCS objects that are located within the eROSITA-DE sky. These candidates are grouped based on the following criteria: objects with BEL detections, classified as either ‘Seyfert’ or ‘LINER’ in at least two BPT diagrams (this criterion was chosen due to the number of candidates showing mixed BPT classifications but displaying BELs); objects with X-ray counterparts; and the various intersections between the categories.

Fig. 17. Venn diagram of the 195 candidates that are in the VCS sample and that lie in the eROSITA-DE sky. The groups are objects with BEL detections (purple); objects classified as AGN (Seyfert or LINER) in at least two BPT diagrams (green); objects with X-ray counterparts (red); and the intersections. “AGN variability only” refers to the candidates in the eROSITA-DE sky falling outside the three mentioned groups.

6.2. Comparison with previous variability-based selections

6.2.1. Comparison by optical properties

Several works have used optical variability to search for AGN in LSMGs, finding candidates with different methods and validation criteria. Here, we compare AGN candidates selected by variability, confirmed either through the detection of BELs or by classification as AGNs in BPT diagrams, with large optical samples. Specifically, we relate our findings to the results from Baldassare et al. (2018, 2020), and Ward et al. (2022).

In Baldassare et al. (2018), AGN candidates were selected from the SDSS spectroscopy within Stripe 82, covering approximately 300 deg². To identify AGN candidates, the full sample of around 28 000 objects was cross-matched with the NSA v1.0.1 catalog, and difference-image photometry light curves were used to assess variability, yielding 135 AGN candidates with redshifts z < 0.15 and galaxy stellar masses in the range of 2 × 10⁸ < M_* < 5 × 10¹¹ M_⊙. When considering only objects with stellar mass M_* < 10¹⁰ M_⊙, the number of candidates reduces to 35, of which 10 (29%) are classified as AGN, 7 (14%) as composite and 16 (46%) as star-forming. Of these 35 candidates, 16 (46%) exhibit H_α BELs, with estimated black hole masses ranging from 10^6.1 M_⊙ to 10^7.9 M_⊙. In our study, when limiting the host galaxy mass to M_* < 10¹⁰ M_⊙, the VCS sample includes 161 candidates. Of these, 60 (37%) are classified as AGN based on BPT diagrams, 24 (15%) as composite, and 69 (43%) as star-forming. Among the 161 candidates, 127 show BELs, and the estimated black hole masses range from 10^5.3 M_⊙ to 10^8.0 M_⊙. In this comparison, the success rate of confirmed AGNs is higher in our study when considering the detection of BELs in LSMGs (M_* < 10¹⁰ M_⊙): 46% in Baldassare et al. (2018) versus 79% in the VCS sample. Further, restricting our sample to the SDSS Stripe 82 sky area we obtain 13 candidates with M_* < 10¹⁰ M_⊙. Of these, 12 (92%) exhibit BELs, 7 (54%) are classified as AGN by BPT diagrams, and 6 (46%) as star-forming.

Our sample is more easily compared to the low-redshift sample of Baldassare et al. (2020), who investigated variability-selected AGNs in low-mass galaxies by analyzing low-amplitude variability in the r-band light curves from the Palomar Transient Factory (PTF). They identified 417 candidates using the NSA v0 catalog, selecting objects with z < 0.055 and stellar mass M_* < 10¹² M_⊙. When constraining the candidates to M_* < 10¹⁰ M_⊙, they have 237 candidates, all with SDSS spectral data. Of these 25% are AGN/composite, and 12 (5%) exhibit BELs. If we restrict the redshift of our VCS sample to z < 0.055, and the stellar mass to M_* < 10¹⁰ M_⊙, we have 55 candidates, of these 25 (45%) are AGN/composite, and 30 (54%) exhibit BELs.

The third set for comparison is the one studied by Ward et al. (2022). These authors identified 44 AGN candidates selected by optical variability from ZTF light curves in LSMGs from the NSA v1.0.1 catalog, with a stellar mass limit set at M_Petrossian < 3 × 10⁹ M_⊙/h². Because of computational limitations, only about a third of the available galaxies were studied (25 714 out of 81 462). These 44 candidates have z < 0.15 and M_* < 10^9.75 M_⊙, with 4 (9%) classified as AGN by BPT diagrams, 4 (9%) as composite, and 35 (81%) as star-forming. Six (14%) of the 44 candidates exhibited BELs and the estimated black hole masses range from 10^6.3 M_⊙ to 10^7.5 M_⊙. Applying the same stellar mass restriction, M_* < 10^9.75 M_⊙, to our VCS sample, we obtain 57 candidates, of which 17 (30%) are classified as AGN, 5 (9%) as composite, and 33 (58%) as star-forming in BPT diagrams, while 37 (65%) exhibit BELs, with black hole masses estimated to range from 10^5.3 M_⊙ to 10^7.7 M_⊙. Comparatively, our success rate in confirming AGNs in LSMGs is significantly higher than that reported in Ward et al. (2022), particularly regarding the detection of BELs (i.e., 6/44 in Ward et al. 2022, 37/57 in the present work).

Since the Ward et al. (2022) selection was also based on ZTF light curves in the g band, we directly compare our samples further. Cross-matching the 44 AGN candidates selected from optical variability by Ward et al. (2022) to our AGN candidate list, we find only five matches. These all belong to the Alerts set and the only one of the five that falls inside our forced photometry footprint is also in the Forced Photometry set. All five have significant BELs in our analysis and are classified as QSO/broad-line by the SDSS pipeline, with four of the five carrying the additional label of starburst. Of the 44 AGN candidates identified by Ward et al. (2022), 18 have Dec < 15.5 deg and we had produced forced photometry light curves for 11 of these. Our forced photometry classifier returned a classification of AGN for one of them (the one included also in the Alerts set as discussed above), a transient class for another one, and a nonvariable galaxy-class for the remaining nine. We checked the SDSS spectra of these 11 galaxies, finding significant BELs in only two, the ones classified as AGN and transient above, while all the others only show narrow emission lines and BPT classifications of either star-forming type (8/9) or composite (1/9). In summary, the AGN candidates identified by both our classifier and Ward et al. (2022) are spectrally confirmed, while the ones identified by Ward et al. (2022) and where our classification based on forced photometry light curves produces a nonvariable galaxy-class are not confirmed as AGN by their spectra. The one case where Ward et al. (2022) identified an AGN candidate that we classified as transient and not as AGN with our forced photometry classification is confirmed as an AGN by its spectrum. Therefore, it is possible that our classification has missed some AGN by assigning them a transient class rather than an AGN class.

6.2.2. Comparison by X-ray detection

We now consider the detection of X-rays in the AGN candidates as an independent validation criterion, following the work presented by Arcodia et al. (2024). This study reported very low X-ray counterpart match rates among optical variability-selected AGN candidates in LSMGs, notably lower than the rates for our sample (Secs. 5 and 6.1). We compare our selection to the largest optical sets studied by Arcodia et al. (2024), namely Baldassare et al. (2020), Kimura et al. (2020), Burke et al. (2022) and the optically-selected sample by Ward et al. (2022). As we show below, the lack of X-ray counterparts in some cases is caused by shallow X-ray catalog sensitivity limits, and differences in selection in other cases.

Burke et al. (2022) searched for variable galaxies in Dark Energy Camera images, with a ∼7-day cadence and 6-year baseline in three different fields. These data are much deeper and also much more restricted in sky coverage than our sample. They selected low-mass black holes on the characteristic timescales of Damped Random Walk models fitted to the light curves, with shorter timescales expected for low-mass black holes. Of 46 AGN candidates in the eROSITA-DE sky with short variability timescales (i.e., log(τ/day) <  = 1.5) Arcodia et al. (2024) only found one X-ray counterpart at the depth of eRASSv1.1. These AGN candidates are considerably dimmer than our sample, with a median g-band magnitude of 21.3 and a standard deviation of 1.5. Therefore, it is likely that the lack of X-ray counterparts is a product of the limited depth of the X-ray catalog. The number of matches in Arcodia et al. (2024) rises to four when the deeper, stacked X-ray catalog is used.

The sample by Kimura et al. (2020) is even more extreme, with fewer but deeper epochs restricted only to the COSMOS field. In their work, the identification of their variability-selected AGN with X-ray counterparts from deep observations of the Chandra COSMOS Legacy survey is 90%, although no cross-matches are reported for the eROSITA data, pointing again to a limitation in the depth of eRASSv1.1.

As mentioned, our sample is more readily comparable to the low-redshift sample of Baldassare et al. (2020). As in our work, the parent sample was selected from the NSA catalog of local galaxies, and they found 237 AGN candidates in galaxies with stellar mass below 10¹⁰ M_⊙. Of the 52 candidates in that sample that fall inside the eROSITA-DE sky, Arcodia et al. (2024) only found counterparts for three sources (i.e., 6%) at the depth of eRASSv1.1. To make a fairer comparison with the sample in Baldassare et al. (2020), we note that their parent sample was the v0 of the NSA catalog, reaching up to z = 0.055, whereas we used version 1.0.1, including galaxies up to z = 0.15. Applying these same cuts, we have 55 AGN candidates. Of these, 42 fall in the eROSITA-DE sky and 16 have X-ray counterparts within 10″ (i.e., 38%), a lower fraction than our total sample but still much higher than the fraction of X-ray matches found for the variability-selected PTF-NSA AGN sample studied by Arcodia et al. (2024). For completeness, limiting our sample by redshift to z < 0.055 but not by mass (i.e., allowing candidates up to M_* < 2 × 10¹⁰ M_⊙), results in 56 galaxies in the eROSITA-DE sky, of which 26 have X-ray counterparts (i.e., 46%). This small increase in stellar mass changes the fraction of counterparts noticeably. In any case, the random forest algorithms used here, possibly together with the better sampled and higher cadence light curves from ZTF as compared to PTF, are significantly more successful at selecting AGN that do have X-ray counterparts.

Finally, Ward et al. (2022) found 44 AGN candidates selected by optical variability from ZTF light curves in LSMGs of the NSA catalog. Of these 44 candidates, seven were included in the analysis by Arcodia et al. (2024), with no X-ray counterparts found. As mentioned, cross-matching the 44 AGN candidates selected based on optical variability by Ward et al. (2022) with our AGN candidate list yielded only five matches. Unfortunately, none of these matches fall within the eROSITA-DE footprint. However, based on their optical spectral classification, we regard them as correct AGN identifications.

7. Conclusions

In this work we presented the performance of a variability-, color-, and morphology-based random-forest classifier applied to ZTF data to search for IMBH candidates in low stellar mass galaxies (LSMGs) limited by M_* < 2 × 10¹⁰ M_⊙ and redshift z < 0.15 from the NSA v1.0.1 catalog. Here, four different data sets were used to select candidates, while archival spectra from SDSS were used to confirm the AGN activity in the sample. The analysis above shows that the selection based on the variability properties of type I AGN is effective in the majority of candidates selected from well-sampled light curves. The different selections used different light curves from the ZTF data: the alert stream, the complete light curves available through data releases, in g- and r-band separately, and our custom-made forced photometry on the reference-subtracted images provided by ZTF. The Alerts and Forced Photometry light curves can be more accurate since they are extracted from difference images, while the data release and Forced Photometry light curves are more complete than the Alerts set because they include flux measurements for all epochs regardless of the variability amplitude. Our main conclusions are summarized as follows:

1. The cross-match between the variability selection and the NSA v1.0.1 catalog produced 506 AGN candidates in LSMGs. From these, 450 have good-quality SDSS archival spectra. After visual inspection, 415 objects were characterized by fitting the spectrum. Subsequently, applying the criteria of an EW_Hα > 5 Å and the S/N of H_α flux > 3, a total of 357 candidates exhibited a significant BEL. The 357 objects include the two selected by the same criteria, but for the H_β emission line, since the H_α region was missing in the spectra. Additionally, five objects that do not show BELs were classified as AGNs by the BPT and WHAN diagrams, bringing the total number of confirmed AGNs in the VCS sample to 362 out of 415. This means that 87% of the LSMG candidates were confirmed as AGNs. Among the 35 objects with spectra that were not fitted, 22 are classified as AGN by the SDSS pipeline or are included in the quasar catalog of Lyke et al. (2020), but are at a higher redshift than expected from the NSA catalog.

2. From the spectral analysis, 58 candidates in the VCS sample were not classified as AGN: 56 are from the Alerts set and two are from the Forced Photometry set. Additionally, eight objects that were not included in the candidates with good spectra, were also not classified as AGN. After visual inspection of the misclassified light curves, it was found that a few are the product of bad image subtraction, and most are produced by transients in the data. These problems in the candidate selection can be addressed by this analysis with the inclusion of the nuclear transient class in the light curve classifier. It is noteworthy that sets using complete light curves are minimally affected by transients or bad image subtraction. When comparing the performance of Forced Photometry and DR, both are equally pure, but Forced Photometry is capable of finding twice as many AGNs in LSMGs in the same sky region, selecting more objects with lower black hole mass.

3. In the Forced Photometry set, 97% of candidates with good-quality spectra show significant BELs. The spectral analysis of these candidates results in equivalent widths between 20 and 800 Åand FWHM > 1500 km/s. This means that all the broad lines are strongly evident and distinguishable from the narrow emission lines in the spectral profile. This strongly implies that the identification of lower-mass black holes (M_⊙ < 10⁶) is not limited by spectroscopic identification, but by the variability selection from the available photometry data. New facilities with more sensitive light curves, such as the Vera C. Rubin LSST, will allow the discovery of more low-mass black holes from variability selection.

4. The black hole mass estimates were performed for 355 objects using Eq. (1), with mean values of 1.6 × 10⁷ M_⊙ for the Forced Photometry set, 2.1 × 10⁷ M_⊙ for the Alerts set, 2.7 × 10⁷ M_⊙ for DR-g, and 2.8 × 10⁷ M_⊙ for DR-r. Almost all the black hole masses below 3 × 10⁶ M_⊙ are in the Alerts and Forced Photometry sets, showing that the photometry performed on reference-subtracted images is more sensitive to lower black hole masses. This is clear when we note that 10% of the masses in the Forced Photometry set are below 1.8 × 1⁶ M_⊙, while for the DR this limit is at 6.4 − 6.4 × 10⁶ M_⊙.

5. The Eddington ratios for the different sets are mainly distributed within the range of 10⁻² and 10⁻¹, showing that this is the range where the selection method is most sensitive.

6. The classification of the 415 visually cleaned candidates by NEL ratios (i.e., BPT diagrams) results in 185 classified as AGN/Seyfert, composite, or LINER in all three diagrams, while 104 are classified as star-forming in all three BPT diagrams. The remaining show mixed classifications as AGN/Seyfert, composite, or LINER among one or more diagrams. Among the objects consistently classified as star-forming, 59 objects show a significant BEL, proving their type I AGN nature. As a consequence, methods searching for AGN candidates using the BPT AGN class as an obligatory criterion will miss these objects. The comparison with other methods for measuring the NEL ratios suggests that a more careful fitting provides better constraints for the BPT classification and can also improve the methods of AGN candidate selection.

7. Cross-matching the eRASSv1.1 catalog of X-ray sources with the 415 variability-selected AGN candidates with good spectra and z < 0.15, we find that 67% of the latter have X-ray counterparts when restricting to the eROSITA-DE sky. Candidates that appear in more than one set are more likely to have an X-ray counterpart. Furthermore, we found that 75% of candidates in the eROSITA-DE sky that show significant BELs have an X-ray counterpart. This result considers all BPT classifications, suggesting that the probability of finding an X-ray counterpart only depends on the presence of BELs and not on the BPT class. This also indicates that the SF classification of objects with X-ray counterparts is likely due to the dominant contamination of the host galaxy emission in the optical spectrum. In addition, as demonstrated by other studies (e.g., Mezcua 2017; Mezcua & Domínguez Sánchez 2020, 2024), observations with higher spatial resolution have been shown to effectively resolve the narrow emission line region of the AGN.

8. Variability-selected AGN candidates with higher AGN relative contributions tend to have bluer AGN continua. However, some objects (ten in total) with flat spectral profiles and evident BELs also show AGN relative contributions between 0.51 and 0.76. When comparing different sets, the Alerts and Forced Photometry find more objects with AGN relative contributions lower than 0.4. Additionally, we find no differences in the AGN relative contribution among different BPT classes. This suggests that emission lines from the host galaxy can contaminate the AGN emission lines even when the AGN has a high relative contribution to the total continuum. Therefore, as mentioned in the previous point, observations with higher spatial resolution are needed to reveal the AGN emission lines component.

9. The ratio of the host galaxy stellar mass to the black hole mass for our sample shows a distribution closer to the relation found for massive elliptical galaxies (M_BH = 0.001M_*) in strong contrast to the results of Reines & Volonteri (2015), for the same stellar mass range. We attribute this to the limitations of the variability selection in identifying black holes with the lowest masses for a given stellar mass host since the more sensitive selection sets (i.e., Alerts and Forced Photometry) find AGNs with lower M_BH/M_* ratios than the other sets. We expect that this limitation might be mitigated with more sensitive and higher S/N data by new observatories like the Vera C. Rubin LSST.

10. A comparison of results from different studies utilizing optical variability to find AGN in LSMGs reveals that the method used here produced a much higher fraction of AGN candidates with eROSITA-DE X-ray counterparts. We argue that, at least for works with similar input data, the main difference is the selection method, in our case using RF algorithms that include several classes of variable and nonvariable objects.

The results presented in this work, aimed at searching for IMBH in LSMGs showed that only a few were consistent with a mass lower than 10⁶ M_⊙. However, the number of confirmed candidates selected by variability demonstrates the high performance of the applied selection method and offers insights on how to improve it. Notably, the number of candidates studied was made possible by the advantages of photometric and spectroscopic surveys. Therefore, future surveys and new generations of instruments will provide better and significantly more data for similar studies. These advancements, together with the refinement of selection techniques, will increase the number of IMBH detections, helping to reveal the history of AGNs in low stellar mass galaxies.

Data availability

Full Table 4 is available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/694/A127

Acknowledgments

Part of this work was supported by the European Southern Observatory (ESO) project SSDF 28/23D. The authors acknowledge the National Agency for Research and Development (ANID) grants: Programa de Becas/Doctorado Nacional 21212344 (SB); Millennium Science Initiative Programs NCN2023_002 (SB, PA, PL), ICN12_009 (FEB) and AIM23-0001 (FEB, PA); and Fondecyt Regular 1241422 (PA, FEB, PL), 1240875 (PL, PA), and 1241005 (FEB, PA), and CATA-BASAL FB210003 (FEB). PA acknowledges support from Centro de Astrofísica de Valparaíso – CAV, CIDI N. 21 (Universidad de Valparaíso, Chile). Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is www.sdss4.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) /University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University. Based on observations obtained with the Samuel Oschin 48-inch Telescope at the Palomar Observatory as part of the Zwicky Transient Facility project. ZTF is supported by the National Science Foundation under Grant No. AST-1440341 and a collaboration including Caltech, IPAC, the Weizmann Institute for Science, the Oskar Klein Center at Stockholm University, the University of Maryland, the University of Washington, Deutsches Elektronen-Synchrotron and Humboldt University, Los Alamos National Laboratories, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, and Lawrence Berkeley National Laboratories. Operations are conducted by COO, IPAC, and UW. Based on observations obtained with the Samuel Oschin Telescope 48-inch and the 60-inch Telescope at the Palomar Observatory as part of the Zwicky Transient Facility project. ZTF is supported by the National Science Foundation under Grants No. AST-1440341 and AST-2034437 and a collaboration including current partners Caltech, IPAC, the Weizmann Institute for Science, the Oskar Klein Center at Stockholm University, the University of Maryland, Deutsches Elektronen-Synchrotron and Humboldt University, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, Trinity College Dublin, Lawrence Livermore National Laboratories, IN2P3, University of Warwick, Ruhr University Bochum, Northwestern University and former partners the University of Washington, Los Alamos National Laboratories, and Lawrence Berkeley National Laboratories. Operations are conducted by COO, IPAC, and UW. Part of this work is based on data products from eROSITA aboard SRG, a joint Russian-German science mission supported by the Russian Space Agency (Roskosmos), in the interests of the Russian Academy of Sciences represented by its Space Research Institute (IKI), and the Deutsches Zentrum für Luft- und Raumfahrt (DLR). The SRG spacecraft was built by Lavochkin Association (NPOL) and its subcontractors, and is operated by NPOL with support from the Max Planck Institute for Extraterrestrial Physics (MPE). The development and construction of the eROSITA X-ray instrument was led by MPE, with contributions from the Dr. Karl Remeis Observatory Bamberg & ECAP (FAU Erlangen-Nuernberg), the University of Hamburg Observatory, the Leibniz Institute for Astrophysics Potsdam (AIP), and the Institute for Astronomy and Astrophysics of the University of Tübingen, with the support of DLR and the Max Planck Society. The Argelander Institute for Astronomy of the University of Bonn and the Ludwig Maximilians Universität Munich also participated in the science preparation for eROSITA.

References

Appendix A: Higher redshift AGN

AGN candidates at higher redshift that are included in the NSA low stellar mass sample because their redshift was incorrect are listed in Table A.1, together with the spectroscopic redshifts and classification provided by the SDSS pipeline.

Another group of six AGN candidates also have spectra consistent with higher redshift QSO but have been miss-classified by the SDSS pipeline and assigned lower redshifts. These are J032620.06-064221.7, J120948.84+224432.5, J023259.60+004801.7, J163844.20+242601.9, J144540.13+113521.0, all classified as GALAXY by the SDSS pipeline and J101336.37+561536.3, classified as low-redshift QSO. The corresponding spectra are plotted below. These targets appear in the SDSS DR16Q quasar catalog with visual classifications (Lyke et al. 2020), the redshifts and broad type classifications from this table are summarized in Table A.2. The objects in these two tables are considered correct AGN identifications for purposes of establishing the success rate of the classifications but are not included in the study of black hole masses, since their host-galaxy stellar masses are under-predicted.

Fig. A.1. AGN candidates with low-redshift GALAXY or QSO classification by the SDSS pipeline, that we identify as higher redshift AGN.

Fig. A.1. Continued.

Appendix B: Suspected AGN

Here we show the spectra of suspected AGN objects. The list of objects and the corresponding set are presented in Table B.1. The visual inspection of light curves and spectra reveals variability and profiles that can potentially be attributed to an AGN. The spectra profiles are presented in Fig. B.1.

Fig. B.1. Objects classified as suspected AGN. The spectral profiles show BELs and/or blue continuum.

Fig. B.1. Continued.

Appendix C: Comparison with the BPT classification results of the MPA-JHU group

In Table 7 we compare the BPT classification obtained in this study with the Value added catalogs of SDSS. We found that the major difference is for 117 objects that we classify as Seyfert and that are classified as SF in the results of the MPA-JHU group. To explore this difference we compare the flux measured for the different emission-lines used to calculate the ratios, meaning, H_α, H_β, [OIII]λ5007, [NII]λ6584. In Fig. C.1 it is evident that flux measurements are different for emission-lines where a BEL is present. In the comparison of [OIII]λ5007, the agreement is clear with the exception of one object (J000926.40+001932.1). For this and the objects that show more difference in the fluxes, we display the quality of our fits in Fig. C.2. This comparison and the quality of fits leads us to interpret the difference in classification as a product of the difference in the method used to measure the flux of the H_α and H_β. In particular, the method used by the MPA-JHU group does not distinguish between NELs and BELs, resulting in an increase in the flux of the narrow H_α and H_β emission lines.

Fig. C.1. Comparison of fluxes, between this work and the MPA-JHU group, for different emission-lines. The fluxes obtained by this study are in the x-axes and the ones obtained by the MPA-JHU group are in the y-axes. We show the comparison with only the narrow component of our model (flux narrow), and for the cases of Hβ, Hα and Hα complex we also compare with the total narrow+broad components (flux narrow+broad). Flux units in all cases are 10−17ergcm−2s−1. The gray lines indicate the one-to-one relation.

Fig. C.1. Continued.

Fig. C.2. Best-fit models for objects that show more difference in the comparison of Hβ flux between our fits and the MPA-JHU group results.

Fig. C.2. Continued.

Fig. C.2. Continued.

Appendix D: Flat spectra with relative high AGN contribution

In Sect. 4.4 we compare the relative contribution of the AGN continuum. From the analysis, we found objects that have a flat profile and where the best-fit model results in a power-law slope equal to zero (α = 0). For these objects, we show the best-fit model in Fig. D.1. We find that the stellar absorption lines are well-fitted, this indicates that the decomposition of the continuum is highly accurate.

Fig. D.1. Best-fit models for objects that show flat spectra and high AGN relative contribution. The dashed gray line correspond to zero flux.

Fig. D.1. Continued.

Fig. D.1. Continued.

All Tables

All Figures

	Fig. 1. Sky distribution of AGN candidates in low stellar mass galaxies. The different selection sets are marked by different colors, whereas AGN candidates selected in multiple sets have their symbols overlaid.
In the text

	Fig. 2. Venn diagram showing the number of elements within each intersection between the different sets. The number of objects included only in the Alerts set is higher than in the other sets due to the larger sky region coverage of the Alerts set.
In the text

	Fig. 3. NSA v1.0.1 catalog redshift distributions for all the objects in each of the four sets: orange for Alerts, blue for Forced Photometry, green for ZTFDR11 g-band, and red for ZTFDR11 r-band.
In the text

	Fig. 4. Apparent rest-frame magnitude distribution of AGN candidates in low stellar mass galaxies for the g* (top) and r* (middle) bands. The bottom panel shows the color g* − r* distributions. Each set of variability-selected AGN candidates is plotted in different colors: orange for Alerts, blue for Forced Photometry, green for ZTFDR11 g-band, and red for ZTFDR11 r-band.
In the text

	Fig. 5. Comparison of AGN candidates with SDSS spectra (blue) and without spectra (orange). The top panel shows the distribution of these samples in the sky and the bottom panel shows the distribution of their redshifts from the NSA catalog.
In the text

	Fig. 6. Distribution of the S/N of the flux of the broad component of Hα. All the variability-selected AGN candidates in low-mass galaxies are plotted in orange; of these, only the ones with EWHα > 5 Å are plotted in blue.
In the text

Fig. 7. Distribution of the equivalent width of the broad component of Hα (top) and of black hole masses derived using Eq. (1) (bottom) for objects with EWHα > 5 Å. In both histograms, each set of variability-selected AGN candidates is plotted in different colors: orange for Alerts, blue for Forced Photometry, green for ZTFDR11 g-band, and red for ZTFDR11 r-band. In the top panel, the vertical line indicates EWHα = 5 Å.

In the text

	Fig. 8. Distribution of the Eddington ratios derived for each set of AGN candidates. Each set is plotted in different colors: orange for Alerts, blue for Forced Photometry, green for ZTFDR11 g-band, and red for ZTFDR11 r-band.
In the text

Fig. 9. BPT diagrams from the measured flux of narrow emission lines. The color bar indicates the black hole mass on a logarithmic scale. The black triangles represent the objects with no mass estimation. In the top left corner, the black lines represent the mean errors. There are four objects classified as Seyfert (LINER) in all three BPT diagrams with no mass estimation (green diamonds). For these objects the WHAN classification is presented in Table 6. Furthermore, the object represented as a green square has an X-ray counterpart in the eRASSv1.1 catalog.

In the text

	Fig. 10. WHAN diagram for the five objects resulting in the AGN region on BPT diagrams, but showing no BELs. The blue markers show the position of these objects in the diagram. Three of them fall in the S-AGN region and two in the W-AGN region.
In the text

Fig. 11. Here we present the AGN-continuum slope and the AGN contribution to the continuum. Center: comparison of the AGN-continuum slope from the best fit with the ratio of the AGN component and total continuum at 5100 Å. The different colors represent the different sets: orange for Alerts, blue for Forced Photometry, green for DR-g, and red for DR-r. The black square and error bars represent the typical error in the AGN-continuum slope constrained by the 16th and 84th percentile values from the simulations. In the top panel, the distribution of the AGN relative contribution (f5100, AGN/f5100, cont.) is plotted for each set. Additionally, in the right panel, the power-law slope (γ) distribution of each set is shown.

In the text

Fig. 12. Here we present the distribution of the AGN contribution to the continuum and the AGN-continuum slope for the different BPT classes. Top: distribution of the AGN relative contribution (f5100, AGN/f5100, cont.). For different BPT classes, according to the [OIII]/Hβ vs. [N II]/Hα diagram. The distributions for different classes, star-forming (purple), composite (pink), and Seyfert/LINER (brown) are similar to the Seyfert/LINER objects concentrated between 0.3 < f5100, AGN/f5100, cont. < 0.75. Bottom: Distribution of the fitted power slope for the same BPT classes.

In the text

Fig. 13. Stellar mass and black hole mass comparison. Top: distribution of the ratio of our estimates of the black hole mass to the stellar mass of the host galaxies listed in the NSA catalog, separated by selection set. Bottom: black hole mass estimated using Eq. (1) vs. the stellar mass of the host galaxy. The results from Reines & Volonteri (2015) are shown in gray for comparison. The purple line represents the MBH = 0.001M* relation (Kormendy & Ho 2013), while the black dashed line indicates the relation obtained in Reines & Volonteri (2015). The mean errors are shown with error bars in the bottom right corner.

In the text

Fig. 14. Comparison of the black hole mass computed by Eq. (1) in this work (i.e., using the FWHM of BELs and AGN luminosity at 5100 Å) vs. the mass computed using Eq. (1) of Reines & Volonteri (2015) (i.e., using the Hα FWHM and its luminosity). The purple line shows the 1:1 relation in logarithmic scale. The dashed purple lines indicate the 0.5 dex uncertainty from Eq. (1) of Reines & Volonteri (2015). The mean errors are shown in the top left with error bars. Since both mass estimates use the FWHM of Hα, the errors are partly correlated

In the text

	Fig. 15. Distribution of Hα equivalent widths for the AGN candidates with fitted spectra in the eROSITA-DE sky (gray), of these the ones with X-ray counterparts (blue), and of these the ones classified as star-forming by the BPT diagnostics (red). Objects with EWHα = 0 Å, 12 in total, are counted using an added EW value of 0.6.
In the text

	Fig. 16. X-ray luminosity (2–8 keV band) and Hα (narrow+broad) luminosity comparison. The purple line plots the relation found by Ho et al. (2001) for AGN host galaxies, while the red line plots the relation for star-forming galaxies presented in Rosa González et al. (2009). The mean errors are shown in the bottom right with error bars.
In the text

Fig. 17. Venn diagram of the 195 candidates that are in the VCS sample and that lie in the eROSITA-DE sky. The groups are objects with BEL detections (purple); objects classified as AGN (Seyfert or LINER) in at least two BPT diagrams (green); objects with X-ray counterparts (red); and the intersections. “AGN variability only” refers to the candidates in the eROSITA-DE sky falling outside the three mentioned groups.

In the text

	Fig. A.1. AGN candidates with low-redshift GALAXY or QSO classification by the SDSS pipeline, that we identify as higher redshift AGN.
In the text

	Fig. A.1. Continued.
In the text

	Fig. B.1. Objects classified as suspected AGN. The spectral profiles show BELs and/or blue continuum.
In the text

	Fig. B.1. Continued.
In the text

Fig. C.1. Comparison of fluxes, between this work and the MPA-JHU group, for different emission-lines. The fluxes obtained by this study are in the x-axes and the ones obtained by the MPA-JHU group are in the y-axes. We show the comparison with only the narrow component of our model (flux narrow), and for the cases of Hβ, Hα and Hα complex we also compare with the total narrow+broad components (flux narrow+broad). Flux units in all cases are 10−17ergcm−2s−1. The gray lines indicate the one-to-one relation.

In the text

	Fig. C.1. Continued.
In the text

	Fig. C.2. Best-fit models for objects that show more difference in the comparison of Hβ flux between our fits and the MPA-JHU group results.
In the text

	Fig. C.2. Continued.
In the text

	Fig. C.2. Continued.
In the text

	Fig. D.1. Best-fit models for objects that show flat spectra and high AGN relative contribution. The dashed gray line correspond to zero flux.
In the text

	Fig. D.1. Continued.
In the text

	Fig. D.1. Continued.
In the text