A Population of Bona Fide Intermediate-mass Black Holes Identified as Low-luminosity Active Galactic Nuclei

An AGN creates specific signatures in an optical spectrum of a galaxy (Elvis 2000): X-ray photons from the corona of an accretion disk around the BH ionize gas out to a few kiloparsecs away, which then produces easily detectable emission lines (Figure 1). The width and the flux of an allowed recombination line (e.g., hydrogen Hα) emitted from the broad-line region (BLR) in the immediate vicinity of a central BH provide a virial estimate of its mass (Greene & Ho 2005; Reines et al. 2013). We have designed an automated workflow that uses data mining in optical and X-ray astronomical data archives publicly available in the international Virtual Observatory to search for AGN signatures of IMBHs.

Zoom In Zoom Out Reset image size

The workflow automatically analyzed (Figure 1) about 1,000,000 optical spectra of galaxies and quasars from the legacy sample of the SDSS DR7 (Abazajian et al. 2009) without any preselection on host galaxy luminosity or redshift. We used a nonparametric representation of a narrow emission line profile (Chilingarian et al. 2017), which produced lower fitting residuals in Balmer lines and allowed us to detect fainter broad-line components, thus boosting the sensitivity of our analysis technique. Then, we took the resulting sample of galaxies with broad-line detections, computed emission-line flux ratios [O iii]/Hβ and [N ii]/Hα from narrow-line components, and used the Baldwin–Phillips–Terlevich (BPT) (Baldwin et al. 1981) diagnostics to reject objects where the ionization was induced by star formation, because such objects often have broad Balmer lines originating from transient stellar events (Baldassare et al. 2016) rather than from AGNs. After that, we eliminated statistically insignificant measurements by filtering the sample (see details below) on relative strengths and widths of narrow- and broad-line components, signal-to-noise ratios, and relative radial velocity offsets and assembled a list of 305 candidates in the IMBH mass range (M_BH < 2 × 10⁵ M_⊙).

This list (hereafter referred to as the parent sample) included all known nuclear IMBH candidates (Dong et al. 2007; Reines et al. 2013; Baldassare et al. 2015) except those that fell on the star-forming sequence in the BPT diagram; our BH mass estimates agreed within uncertainties with those obtained from dedicated deep spectroscopic observations (Reines et al. 2013; Baldassare et al. 2015). Then, we searched in data archives of Chandra, XMM-Newton, and Swift orbital X-ray observatories and detected 10 X-ray counterparts for candidates with virial masses as low as 4.3 × 10⁴ M_⊙. Five of them were mentioned in the literature (Reines et al. 2013; Baldassare et al. 2015), and one object had a spatially extended X-ray emission probably related to star formation, which we excluded from further analysis. All four remaining X-ray point sources were observed serendipitously.

We developed a dedicated technique to analyze optical spectra, which allows us to estimate a central BH mass by measuring a broad component in allowed emission lines (Greene & Ho 2005). As an input data set for our study we use galaxy spectra from SDSS DR7, which we reanalyzed and presented in the value-added catalog RCSED (Chilingarian et al. 2017). We extended RCSED by adding AGN spectra classified as quasars in SDSS. This input sample contains 938,487 galaxy spectra of 878,138 unique objects. We analyzed the follow-up optical spectra obtained with the 6.5 m Magellan telescope in the same fashion. We fitted and subtracted stellar continuum in SDSS spectra using the NBursts full spectrum fitting technique (Chilingarian et al. 2007a, 2007b) and then measured emission lines.

The core of our analysis method is a simultaneous fitting of all strong emission lines (Hβ, [O iii], [N ii], Hα, [S ii]) by a linear combination of a narrow-line component having a nonparametric shape and a broad Gaussian component in the Balmer lines. The broad-line component parameters, velocity dispersion (σ_BLR), and central radial velocity of the BLR component are fitted in a nonlinear minimization loop. All other parameters are fitted linearly within it at every evaluation of the function using an iterative procedure that includes the following two steps: (i) we determine fluxes of all emission-line components solving a linear problem with the non-negative constraint, and then (ii) we recover the shape of the narrow-line component in a nonparametric way by solving a linear convolution problem (see Chilingarian et al. 2017, for details). The solution is sensitive to noise in the data; therefore, we used regularization in the following form: ${\chi }^{2}+\lambda | | \mathrm{BL}| {| }^{2}$, which requires a smoothness of a solution L, where B denotes the second-order differential operator. By using a nonparametric narrow-line region (NLR) component shape, we can successfully model complex gas kinematics and avoid the degeneracy, which affects the traditionally used multiple Gaussian profile decomposition (Greene & Ho 2005; Reines et al. 2013), because Gaussian functions are not orthogonal and therefore do not form a basis. We then repeated our analysis by excluding a BLR component from the model in order to compare the χ² values between the two approaches and conclude whether adding a BLR component improved the fitting quality in a statistically significant way.

Having obtained the flux and the width of the Hα broad component in all 938,487 spectra of the input sample, we use the conservative empirical calibration to estimate a virial BH mass (Reines et al. 2013):

$$\begin{eqnarray}{M}_{\mathrm{BH}} & = & 3.72\times {10}^{6}\,{({\mathrm{FWHM}}_{{\rm{H}}\alpha }/{10}^{3}\mathrm{km}{{\rm{s}}}^{-1})}^{2.06}\\ & & \times \ {({L}_{{\rm{H}}\alpha }/{10}^{42}\mathrm{erg}{{\rm{s}}}^{-1})}^{0.47}\,{M}_{\odot }.\end{eqnarray} \tag{ 1 }$$

By comparing broad Hα-based virial estimates with other BH mass measurement techniques (e.g., reverberation, stellar dynamics), it was demonstrated (Xiao et al. 2011; Dong et al. 2012b) that they agree within 0.3 dex, i.e., a factor of 2. One, however, has to keep in mind that this uncertainty also includes statistical and systematic errors of BH mass measurements used in the calibration, which might significantly contribute to the 0.3 dex error budget. We use this as a rough estimate of the systematic uncertainty of our method, which also defines the mass range of our search: M_BH < 2 × 10⁵ M_⊙.

We then apply multiple selection criteria in order to filter reliable IMBH candidates from the input sample and eliminate spurious broad-line detections:

• 1.  
  ${M}_{\mathrm{BH}}\lt 2\times {10}^{5}\,{M}_{\odot }$ in order to select objects in the IMBH mass range given the assumed 0.3 dex systematic uncertainty.
• 2.  
  No night-sky airglow lines falling in the regions around Hα+[N ii], [O iii] λ5007, and Hβ λ4861, which we use for the spectral line profile fitting and decomposition.
• 3.  
  An empirical constraint that the width of the broad-line component is at least $\sqrt{5}$ times larger than that of the narrow-line component. Because of the quadratic nature of the second moment of a distribution and also because ${M}_{\mathrm{BH}}\propto {\sigma }_{\mathrm{BLR}}^{2}$, this constraint leaves only those galaxies with the relative BLR-to-NLR contributions to the M_BH measurement exceeding 5:1.
• 4.  
  Signal-to-noise ratio exceeding 3 in every emission line used in the BPT (Baldwin et al. 1981) classification (Hβ, [O iii], [N ii], Hα), which ensures its reliability.
• 5.  
  The BPT classification is "AGN" or "composite" (Kewley et al. 2006), which discards pure star-forming galaxies because broad-line components in them are often transient (Baldassare et al. 2016).
• 6.  
  The Hα/Hβ Balmer decrement for both narrow- and broad-line components <4. According to Dong et al. (2008), only ∼1% of type I AGNs exhibit the BLR Hα/Hβ > 4; however, this constraint allows us to eliminate spurious BLR detections originating, e.g., from a complex Hα profile formed by internal kinematics of a star-forming galaxy without Hβ counterparts because of strong internal extinction.
• 7.  
  Statistical error on M_BH better than 33%.
• 8.  
  $| {v}_{\mathrm{BLR}}-{v}_{\mathrm{NLR}}| \lt 3{\sigma }_{\mathrm{NLR}}$ to reject strongly asymmetric BLR profiles. While this constraint may potentially exclude some legit IMBHs that have not yet settled to the center of the potential well of a host galaxy after a recent merger, this allows us to get rid of underestimated M_BH measurements originating from the approximation of an asymmetric profile by a Gaussian.

These criteria, joined together with Boolean and, form our main selection filter. It leaves a sample of 305 IMBH candidates out of nearly 1 million input spectra. We call it the parent sample.

We carried out follow-up imaging and spectroscopic observations of several IMBH candidates and their host galaxies in the optical and near-infrared domains using the 6.5 m Magellan Baade telescope, Las Campanas Observatory, Chile.

Our primary goal was to obtain quasi-simultaneous optical spectroscopy of the IMBH galaxies selected for follow-up X-ray observations using Chandra and XMM-Newton within a period of 2–6 weeks between observations. Our secondary goal was to obtain the second spectroscopic epoch for several prominent X-ray-confirmed IMBHs and get independent BH mass estimates. Finally, we aimed to take advantage of superior seeing conditions at the Magellan telescope to obtain near-infrared images of several IMBH host galaxies with the spatial resolution 05–07  crucial for the analysis of structural properties, which is 2–3 times better than the resolution of SDSS images. The complete log of our follow-up observations for confirmed IMBHs is provided in Table 1.

For our spectroscopic observations, we used the Magellan Echellette spectrograph (Marshall et al. 2008) with the 10'' long by 07 wide slit that provides cross-dispersion spectroscopy with spectral resolving power λ/Δλ = 6500 or σ_inst = 20 km s⁻¹ in 14 spectral orders covering the wavelength range 0.3 μm < λ < 1.0 μm. Each object was integrated for 40–80 minutes in individual 20-minute-long exposures along either the major or minor axis of its host galaxy. Objects that were small enough to fit in the slit (<5'') were observed along the minor axis, and the sky model was constructed from the "empty" part of the slit. Larger galaxies that did not fit in a 10 slit were observed along the major axis; then, we used offset sky observations of 5 minutes renormalized in flux to match science observations.

We reduced the data using a dedicated MagE data reduction pipeline, which we developed. The pipeline builds a wavelength solution with uncertainties as small as 2 km s⁻¹, merges echelle orders, and creates a flux-calibrated, sky-subtracted merged 2D spectrum, which is then fed to the NBursts spectral fitting procedure to subtract the stellar continuum and then to the emission-line analysis procedure.

The pipeline uses standard stars observed shortly before or after a science source to perform flux calibration and telluric correction. However, in order to perform an extra check and eliminate possible systematic flux calibration errors, we compare our reduced spectra to SDSS and use stellar continuum to perform independent flux calibration. We first extract a spectrum in a 3 × 07  box and use a published light profile of each galaxy (Simard et al. 2011) in order to calculate the expected flux difference between the circular 3'' SDSS fiber aperture and the box extraction. Then, we calibrate a reduced 2D spectrum using that flux ratio, and we perform an optimal extraction of a nuclear point source using the value of the image quality reported by the guider in order to estimate an extraction profile, because the BLR in an AGN is supposed to stay unresolved. At the end, we apply a flux correction computed for a point source observed through a 07 wide slit to the extracted spectrum. This approach yields a flux-calibrated spectrum of the galaxy nucleus that can be directly compared to SDSS.

For our imaging observations we used the FourStar camera (Persson et al. 2013), which covers a field of view of 11 × 11'' with a mosaic containing four Hawaii2-RG detectors. We observed each IMBH host galaxy from a subsample selected for imaging in the K_s band with a total on-source integration of 8 minutes. We used the Poisson random dithering pattern with a box size of 52'' in order to provide enough background sampling for flat-fielding and background subtraction. We reduced FourStar images using the fsred data reduction pipeline, which performs preprocessing of raw near-IR images obtained in the fowler2 mode, that is, two read-outs in the beginning and two at the end of each exposure; flat-fielding; background subtraction; and flux calibration using the Two Micron All Sky Survey (Skrutskie et al. 2006) sources inside the field of view. The final result of the pipeline is a sky-subtracted, flux-calibrated image and its flux uncertainties.

For one galaxy, SDSS J1714+5849, we used archival HST images in the F814W photometric band downloaded from the Hubble Legacy Archive (http://hla.stsci.edu/; data set JA2S0M010). For SDSS J1227+0757, which we did not observe with FourStar, we used Pan-STARRS archival images (Chambers et al. 2016) with subarcsecond seeing quality.

In order to check the position of all our candidates on the ${M}_{\mathrm{BH}}\mbox{--}{M}_{\mathrm{bulge}}$ scaling relation, we analyzed imaging data for the candidate IMBH host galaxies. For all of them but one we have used a two-dimensional decomposition using the galfit v.3 software (Peng et al. 2010). For one galaxy, SDSS J1714+5849, which harbors a strong stellar bar, we used instead a one-dimensional decomposition of a light profile extracted using the ellipse task in noao iraf. For SDSS J1342+0530, SDSS J1107+1347, and SDSS J0228−0901 the follow-up imaging data taken on FourStar have been used. We used psfextractor (Bertin 2011) to extract the point-spread function convolution kernel for every image. Usually we found the best-fitting solution with the simple photometric model "bulge+disk." In case of compact bulges being limited by the atmospheric seeing quality, we had to model a bulge using a point source. Then, apparent magnitudes of bulges were translated into luminosities using the WMAP9 cosmology (Hinshaw et al. 2013) and later converted into stellar masses using stellar population properties provided in the RCSED (Chilingarian et al. 2017).

We performed X-ray observations of two objects with Chandra (observations 20114 and 20115) and two objects with XMM-Newton (observations 0795711301 and 0795711401) using director's discretionary time quasi-simultaneously with optical observations.

Both Chandra observations were carried out with the Advanced CCD Imaging Spectrometer (ACIS) detector in the faint data mode with 10,000 s long exposures. Target galaxies were always placed on-axis of the back-illuminated S3 chip of ACIS-S. The data were reduced and analyzed with the ciao 4.9 package following standard recipes.

In Chandra observation 20114 a single bright point-like X-ray source was detected at the position of the optical center of SDSS J1107+1347 galaxy. We performed its aperture photometry using the srcflux task and detected 518 net counts in the 0.5–7 keV band, which corresponds to the observed flux (4.5 ± 0.4) × 10⁻¹³ $\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$.

No source was detected at the position of SDSS J135750.71+223100.8 in the Chandra observation 20115. We estimate a 3σ detection limit of this observation as 8.0 × 10⁻¹⁵ $\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$.

Two XMM-Newton observations, each 10,000 s long, were performed with the EPIC detector in FullFrame mode with the Thin filter. The data were reduced and analyzed with the common XMM-Newton analysis threads, with the sas 16.1.0 software package running in a virtual machine.

We were not able to detect any source at a position of SDSS J161251.77+110621.6 in the XMM-Newton observation 0795711301 up to the limiting flux level of 5.0 ×10⁻¹⁵ $\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$.

However, we marginally detected a faint source coincident within positional errors with the nucleus of SDSS J1440+1155 in the XMM-Newton observation 0795711401. We estimate its flux in the standard XMM-Newton 0.2–10 keV band as (5 ± 2) × 10⁻¹⁵ $\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$.

For other sources in this study (including previously known objects) we used X-ray data from XMM-Newton's 3XMM-DR7 catalog (Rosen et al. 2016), accessible through the catalog website (Zolotukhin et al. 2017, http://xmm-catalog.irap.omp.eu), the Chandra Source Catalog Release 1.1 (Evans et al. 2010), and the Swift 1SXPS catalog (Evans et al. 2014). We performed a cross-match with a point nonspurious subset of sources in those catalogs using their 3σ X-ray positional uncertainties.

Using data mining in wide-field sky surveys and applying dedicated analysis to archival and follow-up optical spectra, we identified a sample of 305 IMBH candidates having masses 3 × 10⁴ M_⊙ < M_BH < 2 × 10⁵ M_⊙, which reside in galaxy centers and are accreting gas that creates characteristic signatures of a type I AGN. Out of these 305 candidates, we confirmed the AGN nature of 10 sources, including five previously described in the literature (Dong et al. 2007; Reines et al. 2013; Baldassare et al. 2015), by detecting the X-ray emission from their accretion disks, thus defining the first bona fide sample of IMBHs in galactic nuclei. The very existence of nuclear IMBHs supports the stellar-mass seed scenario of the massive BH formation.

In Table 2 we present main properties of 10 IMBHs confirmed as AGNs by the X-ray identification and their host galaxies. Every object is identified by the IAU designation, which includes its J2000 coordinates. For every source we present a central BH virial mass estimate, flux and width of the broad Hα component, redshift, X-ray flux, and an estimated stellar mass of a spheroidal component. For host galaxies of five newly detected sources presented in the top part of the table, we also provide estimates of the absolute magnitude of the bulge or spheroid obtained from the photometric decomposition of their direct images. For the confirmed sources from the literature (bottom part of the table) we also provide published BH mass estimates. In Figures 3–6 we present SDSS and MagE (when available) spectra and line profile decomposition results for these 10 objects.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Here we briefly describe properties of bona fide IMBHs detected in X-ray for the first time:

• 1.  
  J122732.18+075747.7: The least massive IMBH (M_BH = 3.6 × 10⁴ M_⊙) detected by our workflow hosted in a barred spiral galaxy with a star-forming ring; the X-ray counterpart is very faint. The BPT diagnostics places a galaxy in the composite region because the AGN emission is heavily contaminated by star formation in the inner ring, which becomes less of a problem in MagE data, where it is spatially resolved.
• 2.  
  J134244.41+053056.1: A particularly interesting source, which we matched with the Swift source 1SXPS J134244.6+053052. Yang et al. (2013) classified it as a TDE candidate based on the variability of highly ionized iron lines in its optical spectrum. The claim is that the TDE must have happened close to the SDSS spectrum epoch (2002 April 9), which is, however, in clear contradiction with the hypothesis that its X-ray emission with a luminosity of 1.3 × 10⁴¹ erg s⁻¹ observed by Swift 7 yr later on 2009 May 15 is connected to the TDE. Therefore, we attribute the detected X-ray source to the AGN activity in SDSS J1342+0530.
• 3.  
  J171409.04+584906.2: Hosted in a barred spiral with a compact bulge well resolved in archival HST images, this IMBH is another example of a source falling into the composite region in the BPT diagram. Because of high declination, we were unable to obtain the second-epoch spectroscopy.
• 4.  
  J111552.01−000436.1: Located in a nearly edge-on spiral galaxy with a compact bulge, this is another example of a weak AGN whose signature is contaminated by star formation in its host galaxy.
• 5.  
  J110731.23+134712.8: This IMBH located in a low-luminosity disk galaxy with a very compact bulge is the only one of five falling in the AGN region of the BPT diagram despite its contamination by star formation. This object has a very bright X-ray counterpart detected in our Chandra data set, which corresponds to the X-ray luminosity alone over 10% of the Eddington limit for a 70,000 M_⊙ BH, which suggests that the bolometric luminosity should be close to the Eddington limit.

We also mention one object previously described in the literature, J022849.51−090153.8 (Greene & Ho 2007). Its BH mass estimate from the follow-up spectroscopy with MagE, (3.7 ± 0.3) × 10⁵ M_⊙, puts it above the IMBH mass threshold adopted in this work; however, similarly to J110731.23+134712.8, it also exhibits very bright X-ray emission, which corresponds to the bolometric luminosity close to the Eddington limit for its mass.

Table 3 contains properties of all 305 sources selected as IMBH candidates (the parent sample) regardless of the availability of X-ray data. The columns are similar to Table 2, with the exception of X-ray identification and literature data. We used bulge luminosities in the r band reported in the photometric catalog by Simard et al. (2011) in the "Sérsic+disk" decomposition table (pure exponential disk and n = 4 de Vaucouleurs bulge), which we converted into stellar masses using mass-to-light ratios of SDSS galaxies presented in Saulder et al. (2016). Our parent sample includes five objects: J091424.76+115625.6, J093401.24+245342.5, J130141.56+100100.1, J144850.09+160803.2, andJ160531.85+174826.2, which have recently been identified as candidate low-mass BHs by Liu et al. (2018). In Figure 7 we display the BPT diagnostics diagram for our parent sample and IMBHs confirmed in X-ray. The vast majority of the 305 objects fall into the composite region of the diagram, which is not surprising because they are mostly disk galaxies with ongoing star formation.

Zoom In Zoom Out Reset image size

Table 3.  List of 305 Candidate IMBHs Identified as AGNs Based on the SDSS Archival Data

Object	z	MBH	${\sigma }_{\mathrm{BLR}}$	LBLR Hα	${M}_{\mathrm{abs}}^{\mathrm{sph}}$
		(103 M⊙)	(km s−1)	(1039 erg s−1)	(109 M⊙)
J111835.82+002511.2	0.025	138 ± 20	230 ± 13	13.24 ± 2.07	1.87 ± 0.04
J112209.97+010114.8	0.075	99 ± 19	216 ± 14	86.3 ± 2.51	1.94 ± 0.3
J141215.60−003759.0	0.025	62 ± 17	269 ± 30	1.24 ± 0.39	0.46 ± 0.05
J094733.06+001302.9	0.063	181 ± 39	322 ± 30	5.47 ± 1.15	10.65 ± 2.50
J003826.70+000536.8	0.071	100 ± 22	257 ± 25	4.06 ± 0.87	n/a

(This table is available in its entirety in FITS format.)

Download table as:  ASCIITypeset image

We studied the behavior of our nonparametric emission-line analysis algorithm with Monte Carlo simulations. For each object in our final sample of 10 X-ray-confirmed IMBHs we generated a set of 90,000 mock emission spectra, which included all strong optical emission lines (Hβ, [O iii], [N ii], Hα, [S ii]). To generate synthetic forbidden lines in these spectra, we took their model profiles derived from the narrow-line line-of-sight velocity distribution, which was recovered at the first pass of the algorithm, and added random noise with the distribution derived from the observed spectrum noise in the vicinity of a line. For allowed lines we also added broad Gaussian components with a random Balmer decrement in the range 2.8–3.2, whose widths σ and luminosities L_Hα were distributed in a grid to cover the region of interest in the parameter space. At each point of this grid we generated 100 spectra with random noise realizations that were fed to the nonparametric emission-line analysis algorithm. M_BH recovered by the algorithm was then compared to the true input value used. We considered an individual trial successful in case the recovered BH mass fell within 0.3 dex of the input value, and we computed the recovery rate at each point of the (σ, L_Hα) grid as a fraction of successful to total number of trials.

The maps of recovery rate for MagE and SDSS spectra of objects from the final sample of 10 IMBHs are presented in Figure 8. In almost all cases the derived BH masses lie in the region with a reliable recovery rate. This modeling also shows that under favorable circumstances our nonparametric emission-line analysis algorithm could recover from SDSS spectra BH masses as low as 10⁴ M_⊙.

Zoom In Zoom Out Reset image size

We estimate the contamination of the parent sample of IMBH candidates produced by our method using several approaches. By contamination here we mean the fraction of sources in the sample that are not actual IMBHs, that is, they do not exhibit all the required observed IMBH properties: persistent broad Hα emission, X-ray emission from an accretion disk, and the AGN or composite emission-line ratio in the BPT diagram. By construction, our parent sample contains sources with single-epoch spectroscopic mass estimate without X-ray confirmation for most of them. Hence, it inevitably contains sources that, e.g., showed IMBH features at some point in time and then changed their appearance. The contaminating sources likely have diverse origins. These could be supernovae, transient stellar processes (Reines et al. 2013; Baldassare et al. 2016), or spurious detections caused by the imperfection of our spectral analysis.

It is generally very hard to precisely estimate the contamination, so here for simplicity we derive an upper limit of the contamination, i.e., the most pessimistic estimate of the quality of our IMBH search method. It requires that bona fide IMBH candidates satisfy the most stringent observing constraints: they must have an X-ray detection and consistent multiepoch spectroscopic mass estimates (more precisely, we require that broad Hα emission is detected at different epochs, possibly with different instruments, and all its detections satisfy our quality criteria, the mass estimates at different epochs are consistent within 0.3 dex, and the BPT classification does not change).

First, we checked our parent sample against the 3XMM-DR5 upper limit server (http://www.ledas.ac.uk/flix/flix.html) and found 14 objects that serendipitously fell in the footprint of archival XMM-Newton observations but were not detected in them. We compared detection limits of these observations with the fluxes expected from the 14 IMBH candidates given the existing L_X–L_{[O iii]} correlation (Heckman et al. 2005). None of these observations were deep enough to exclude X-ray emission at the level of L_X–L_{[O iii]} correlation minus its 1σ scatter. We therefore cannot reject the accreting IMBH hypothesis for these objects. Given that the objects of our interest are all low-mass AGNs with small luminosities, we expect that other X-ray archives are unlikely to contain many deep-enough observations of our IMBH candidates.

Hence, out of 305 IMBH candidates in our parent sample, only 18 possess sufficiently deep X-ray observations to confirm or rule out the accreting IMBH hypothesis: 14 from X-ray archives and 4 IMBH candidates observed with dedicated X-ray observations by Chandra and XMM-Newton in this work. Out of those 18, 2 sources (SDSS J161251.77+110621.6, observed by XMM-Newton, and SDSS J135750.71+223100.8, observed by Chandra) were not detected in X-ray at a level below that expected from the L_X–L_{[O iii]} correlation minus its 1σ scatter, which secures their non-IMBH nature. One source, SDSS J144005.82+115508.7, while detected in X-ray with XMM-Newton in our observation, did not display broad Hα emission at the second spectroscopic epoch when observed with Magellan/MagE. In addition to this, we discarded five sources without performing the second-epoch spectroscopic follow-up observations, considering them spurious detections. This happens, for example, when our automated emission-line decomposition algorithm underestimates the broad Hα emission flux, and after more careful emission-line decomposition with manually adjusted constraints, a BH mass estimate exceeds 2 × 10⁵ M_⊙. Thus, out of 18 sources with deep X-ray data, we discard 8 sources in total. The 10 remaining sources are those presented in Table 2. For six sources from Table 2 we have both X-ray and second-epoch spectroscopic confirmation (objects found in this work: SDSS J1227+0757, SDSS J1342+0530, SDSS J1115−0004, SDSS J1107+1347; previously known objects: SDSS J1523+1145, SDSS J0228−0901). Four remaining sources from Table 2 (object found in this work: SDSS J1714+5849; previously known objects: SDSS J1534+0408, SDSS J1605+1748, SDSS J1123+6711) are detected in X-ray and have reliable single-epoch detections of broad Hα emission but do not possess second-epoch spectroscopic observations and therefore cannot be used for the contamination estimate.

This leaves us with a sample of 14 sources that have enough evidence to tell whether they pass all required tests (multiepoch spectroscopy and deep enough X-ray observations) or not: 6 sources successfully pass them all, and 8 sources fail in at least one test. The upper limit on the contamination of our parent sample can be estimated as 8/14 = 57%. A more realistic estimate should lower this value. In particular, it was shown (Baldassare et al. 2016) that 100% of objects classified as AGNs on the BPT diagram, which at the same time show broad Hα emission, exhibit the same properties in the second spectroscopic epoch. Out of four objects without the second spectroscopic epoch, one (SDSS J1605+1748; Dong et al. 2007) is classified as an AGN on the BPT diagram. It is then natural to anticipate that it is a true IMBH, which would lower the contaminating fraction in our sample to 53%.

Therefore, we have all evidence to expect that at least 0.43 × 305 = 131 sources in our parent sample of IMBH candidates are real IMBHs in a sense that they will satisfy the most stringent observing criteria once the necessary follow-up observations have been performed. If we assume no strong dependence of the contamination level on the BH mass, we find that at least 42 of our IMBH candidates from Table 3 with M_BH < 10⁵ M_⊙ must be actual IMBHs.

In Figure 9 we compare IMBH masses and host galaxy properties to the recent compilations of bulge/spheroid masses of host galaxies of massive BHs (Graham & Scott 2015; Graham et al. 2015) and to upper limits for central BH masses obtained using stellar dynamics in four low-luminosity galaxies presented in Nguyen et al. (2017, 2018). We also display the "bulgeless" spiral galaxy NGC 4395 hosting a 360,000 M_⊙ BH (Peterson et al. 2005) and measurements of BH masses in tidally stripped compact elliptical (cE; Kormendy et al. 1997; Nguyen et al. 2018) and ultracompact dwarf (UCD) galaxies (Seth et al. 2014; Ahn et al. 2017, 2018; Afanasiev et al. 2018). All IMBH hosts have stellar masses of their bulges between 4 × 10⁸ M_⊙ and 4 × 10⁹ M_⊙ and reside on the low-mass extension of the ${M}_{\mathrm{BH}}\mbox{--}{M}_{\mathrm{bulge}}$ scaling relation established by more massive BHs and their host galaxies (Laor 2001; Graham & Scott 2015), thus filling a sparsely populated region of the parameter space. This argues for the validity of the search approach that looks for AGN signatures created by IMBHs and also supports the connection between the BH mass growth and the growth of their host galaxy bulges via mergers down to the IMBH regime.

Zoom In Zoom Out Reset image size

Galaxy mergers were frequent when the universe was younger (redshifts z > 1; Conselice et al. 2003; Bell et al. 2006; Lotz et al. 2011). They are thought to be responsible for the growth of bulges (Aguerri et al. 2001; Boylan-Kolchin et al. 2006), hence suggesting the coevolution of central BHs and their hosts (Kormendy & Ho 2013), and explaining the observed correlations between central BH masses and bulge properties, stellar velocity dispersion (Ferrarese & Merritt 2000; Gebhardt et al. 2000; van den Bosch 2016), and stellar mass (Marconi & Hunt 2003; Häring & Rix 2004; Gültekin et al. 2009). Therefore, because their bulges are small, IMBH host galaxies must have experienced very few if any major mergers over their lifetime.

The upper limits for BH masses in low-luminosity galaxies displayed in Figure 9 follow the scaling relation established by bulges of massive BHs and extended toward lower masses by IMBHs, hence suggesting that it is only a question of improving the sensitivity of the stellar dynamics approach by a factor of 2–3 before they can be measured and confirmed. On the other hand, cEs and UCDs are clearly offset to lower bulge masses for their M_BH values that reflect their evolutionary status: the host galaxies lost 90%–98% of the stellar bulge mass (Chilingarian et al. 2009; Seth et al. 2014) by tidal stripping. About 4% of all quiescent low-luminosity early-type galaxies in the mass range $5\times {10}^{8}\,{M}_{\odot }\lt {M}_{* }\lt 5\times {10}^{9}\,{M}_{\odot }$ are cE galaxies (Chilingarian & Zolotukhin 2015), which must have lost about 90% of their stars. Most of them likely host "overweight" SMBHs, which will populate the upper part of the ${M}_{\mathrm{bulge}}\mbox{--}{M}_{\mathrm{BH}}$ scaling relation and increase its scatter. One should, therefore, be careful to identify galaxies, which underwent substantial stellar mass loss from their spheroidal components when interpreting this data in the context of the BH–bulge coevolution.

From the multiepoch optical spectroscopy and X-ray observations, we estimate that our IMBH candidate sample may include up to 57% of transient broad lines or spurious detections (see the detailed discussion above), even though, keeping in mind that virial masses are uncertain to a factor of 2, it should contain at least 42 objects with masses smaller than 10⁵ M_⊙. These objects are the relics of the early SMBH formation that survived through the cosmic time almost intact, and their host galaxies must have had very poor merger histories. Their existence suggests that at least some SMBHs did not originate from massive (>10⁵ M_⊙) seeds but from stellar-mass BHs. The efficiency of mass growth via super-Eddington accretion is questionable because of radiation-driven outflows (King & Muldrew 2016). Therefore, stellar-mass BH mergers must have played an important role in the SMBH assembly. We would also like to mention that alternative theories exist, which involve phase transitions (Rubin et al. 2001; Dolgov & Blinnikov 2014) or the modified Affleck–Dine mechanism of baryogenesis in the early universe (Dolgov & Postnov 2017) and result in a continuous mass spectrum of primordial BHs (e.g., lognormal) ranging from a few to a few thousand solar masses. The modified baryogenesis will result in specific peculiarities in the chemical composition of stars, which, in principle, can be detected observationally.

Our sample likely represents the tip of the iceberg of the IMBH population. The sphere of influence of an IMBH is too small to significantly affect the stellar dynamics of its host galaxy and cannot be detected beyond the Local Group with currently available instruments; therefore, we can find only IMBHs in the actively accreting phase, which requires the gas supply onto the galaxy center. On the other hand, most non-star-forming galaxies with small bulges reside in galaxy clusters, which does not favor AGNs, because they lost their gas completely owing to environmental effects. Therefore, while the exact fraction of actively accreting IMBHs is unknown, it is likely smaller than that of more massive BHs (≲1%).

The main limitations of our technique are the relatively shallow flux limit of the SDSS spectroscopy and the lack of wide-field X-ray surveys reaching the flux limit (5 ×10⁻¹⁵ erg cm⁻² s⁻¹) of a snapshot Chandra or XMM-Newton observation: four serendipitously detected sources reside in ≃2% of the sky observed by both SDSS and X-ray satellites. The future deep eRosita X-ray survey may confirm several dozen IMBH candidates from our current sample of 305. A targeted optical spectroscopic probe of nearby galaxies with small bulges deeper than SDSS will likely bring new IMBH identifications at even lower masses.