Estimating AGN Black Hole Masses via Continuum Reverberation Mapping in the Era of LSST

Spectroscopic reverberation mapping (RM) is a direct approach widely used to estimate the mass of black holes (BHs) in active galactic nuclei (AGNs). However, it is very time consuming and difficult to apply to a large AGN sample. The empirical relation between the broad-line region size and luminosity (H$\beta$ $R_{\rm BLR}\unicode{x2013}L$) provides a practical alternative yet is subject to large scatter and systematic bias. Based on the relation between the continuum emitting region size and luminosity ($R_{\rm CER}\unicode{x2013}L$) reported by Netzer (2022), we present a new BH mass estimator via continuum RM (CRM) by comparing $R_{\rm CER}$ and $R_{\rm BLR}$, assuming that the continuum lags are dominated by the diffuse continuum emission. Using a sample of 21 AGNs, we find a tight $R_{\rm BLR}\unicode{x2013}R_{\rm CER}$ relation (scatter$\sim$0.28 dex) and that $R_{\rm BLR}$ is larger than $R_{\rm CER}$ at 5100 \r{A} by an average factor of 8.1. This tight relation enables the BH mass estimation based on the CRM combined with the velocity information. Applying the relation to rest objects in our CRM sample, we demonstrate that the predicted $R_{\rm BLR,CRM}$ follows the existing H$\beta$ $R_{\rm BLR}\unicode{x2013}L$ relation well and the estimated CRM BH masses are consistent with the RM/single-epoch BH masses using H$\beta$. This method will provide significant applications for BH mass estimation thanks to the short continuum lags and the easily accessible high-cadence, large-area photometric data, especially in the era of Legacy Survey of Space and Time.


INTRODUCTION
Black hole (BH) mass is a key parameter in active galactic nucleus (AGN) studies, and is vital for understanding the coevolution of BHs and host galaxies (e.g., Kormendy & Ho 2013), the early formation history of supermassive BHs (e.g., Inayoshi et al. 2020), and accretion physics (e.g., Pringle 1981). To accurately estimate the BH mass, one of the primary methods is observing the motion of stars or gas around the BHs to derive the dynamical mass, which is the most accurate approach for nearby BHs whose influence of sphere can be spatially resolved (e.g., Ghez et al. 2008;Hicks & Malkan 2008;GRAVITY Collaboration et al. 2019). However, this approach is not accessible for more distant galaxies and is difficult to apply to a large sample.
Reverberation mapping (RM, Blandford & McKee 1982) of AGN provides a unique chance to estimate the BH mass in distant galaxies. It measures the time delay (τ ) between the variability of the continuum and the broad-line emission to infer the broad-line region (BLR) size (R BLR = cτ ). Combined with the broad-line width (∆V ) which is a proxy of the virial velocity of gas clouds in the BLR, we can estimate the BH mass following: where G is the gravitational constant and f is the virial factor accounting for the unknown BLR orientation, kinematics, structure, and etc.. The traditional spectroscopic RM requires intensive observational resources and is usually very time-consuming, especially for luminous quasars at high redshift (e.g. Kaspi et al. 2000;Grier et al. 2017;Woo et al. 2019b;Malik et al. 2023). A scaling relation between the BLR size and the continuum luminosity (R BLR -L relation) is established using Hβ based on ∼ 40 local AGNs (e.g., Kaspi et al. 2000;Bentz et al. 2013). Its best-fit slope is 0.533 +0.035 −0.033 (Bentz et al. 2013), fully consistent with the expectation from simple photoionization models. According to the definition of the ionization parameter U H = Q(H) 4πR 2 nHc (Davidson 1972), where Q(H) is the production rate of arXiv:2302.05261v2 [astro-ph.GA] 28 May 2023 hydrogen-ionizing photons that is proportional to the L, and n(H) is the hydrogen density of the gas, the size R that is responsible for the line emission will be proportional to √ U H n H L. The observed Hβ R BLR -L relation is a natural consequence if U H and n H of Hβ-emitting regions are similar among different AGNs.
The R BLR -L relation enables the BH mass estimation using single-epoch (SE) luminosity and broad-line width, which can be easily applied to a large sample of AGNs (e.g., Shen et al. 2011;Liu et al. 2019;Rakshit et al. 2020;Wu & Shen 2022). However, the current R BLR -L relation becomes much more complex and less tight than the canonical one (Bentz et al. 2013). AGNs with super-Eddington accretion rates have significantly smaller BLR size (e.g., Du et al. 2015Du et al. , 2016Du et al. , 2018 than the expectation from canonical R BLR -L relation, which results in overestimation of their SE BH mass using traditional estimators (Du & Wang 2019).
Apart from the BLR, the more extended dusty torus is also found to obey a similar R-L relation (Suganuma et al. 2006;Koshida et al. 2014;Lyu et al. 2019). Recently, Gravity Collaboration et al. (2023) confirmed that the size of hot dust continuum R dust measured using optical/near-infrared interferometry is tightly correlated with the R BLR , with an intrinsic scatter of 0.25 dex. This suggests that the BH mass can be estimated via the hot dust continuum size combined with a broadline width, although the slope of R dust − L relation is found slightly flatter than 0.5. This method is currently limited to the brightest objects with K band brighter than 11 mag. A large sample study will be feasible using the next-generation instrument with upgraded sensitivity and sky coverage (GRAVITY+ Collaboration et al. 2022).
The optical continuum RM (CRM) programs have also made significant progress in resolving the continuum emitting region (CER) sizes of a number of local AGNs with observations through X-ray to UV/optical. (e.g., Fausnaugh et al. 2016;Edelson et al. 2017;Cackett et al. 2018Cackett et al. , 2020Hernández Santisteban et al. 2020;Vincentelli et al. 2021;Kara et al. 2021). Several other works successfully measured the optical continuum lags using modern photometric time-domain surveys (Jiang et al. 2017;Mudd et al. 2018;Homayouni et al. 2019;Yu et al. 2020;Jha et al. 2022;Guo H. et al. 2022a;Guo W. et al. 2022). A key result from these CRM campaigns is that the observed CER sizes (R CER ) are about 3 times larger than the prediction of the standard thin disk model (SSD, Shakura & Sunyaev 1973). One of the popular explanations suggests that the diffuse continuum (DC) emission from the BLR dominates the observed continuum lags (e.g., Korista & Goad 2001Lawther et al. 2018;Netzer 2022). This idea is strongly supported by the U/u band lag excess in the lag spectrum as a function of wavelength (e.g., NGC 4593, Cackett et al. 2018), although not every object present this feature (e.g., Mrk 817, Kara et al. 2021).
Based on several local AGNs, Netzer (2022) proposed the R CER -L relation and suggested that the R CER -L relation is likely a down-scale version of the Hβ R BLR -L relation. This result is quickly confirmed by Guo H. et al. (2022a, hereafter G22) with a much larger sample of AGNs using the light curves from Zwicky Transient Facility (ZTF) and is also extended to the low mass regime (Montano et al. 2022).
In this Letter, we present the tight scaling relation between R CER and R BLR and propose a new method to estimate the BH mass via CRM. We describe the sample in §2 and investigate the feasibility of CRM BH mass estimation in §3. We discuss the advantages and limitations of this new approach, as well as future prospects in the era of Legacy Survey of Space and Time (LSST) in §4. Finally, we draw our conclusions in §5. Throughout this paper, we use the ΛCDM cosmology, with H 0 = 72.0, and Ω m = 0.3.

DATA AND SAMPLE
The CRM sample is mainly from G22 who performed a uniform lag analysis for a large sample of AGNs using gri -band light curves from ZTF (Bellm et al. 2019). ZTF is a time-domain survey that utilizes a 1.2 m telescope at the Palomar Observatory to scan the entire visible sky (DEC ≥ −30 • ) with a cadence of ∼ 3 days since 2018. G22 cross-matched the Million Quasar Catalog (Flesch 2021) with ZTF Data Release 7 that contains the PSFbased photometric light curves from March 2018 to June 2021. A total of 455 spectroscopically confirmed type-I AGNs were selected at redshift z < 0.8 to have wellsampled light curves (N epoch > 20 for each of the three bands) and show good cross-correlation between g and r band light curves (see section 2.1.1 in G22 for details).
The lags of these 455 AGNs were measured using both the interpolated cross-correlation function (ICCF, Peterson et al. 1998b;Sun et al. 2018) and JAVELIN (Zu et al. 2011). The R CER were estimated by fitting a power-law function R CER = R CER,0 λ β to the inter-band (g-r, g-i) lag-wavelength relation. The power-law slope β was fixed to β = 4/3 as suggested by previous AGN CRM studies (Fausnaugh et al. 2016;Jiang et al. 2017;Yu et al. 2020;Homayouni et al. 2019Homayouni et al. , 2022. The R CER at rest-frame 2500Å are available in the Table 4 of G22. We use the high-quality sub-samples selected by G22 based on lag uncertainties, consistency between two lag   Kaspi et al. (2000) . For objects with multiple BLR RM measurements, we use the averaged R BLR,Hβ and logL5100 from the Table 1 by Du & Wang (2019). measuring approaches, lag reliability, etc. A total of 94 objects with reasonably good lag quality were selected as the parent sample, from which 38 objects with the most confident lag measurements were labeled as the core sample. Note that in this work, we use the parent sample to refer to the 56 objects excluding the core sample objects (94 − 38 = 56) and the core sample still represents the best 38 objects. We adopt R CER based on ICCF measurements. Using JAVELIN results will not change our conclusion since the lag difference between ICCF and JAVELIN is required to be small for these two samples.
We supplement G22 CRM sample with 12 local AGNs (Fausnaugh et al. 2016(Fausnaugh et al. , 2018Edelson et al. 2017;Cack-ett et al. 2018Cack-ett et al. , 2020Hernández Santisteban et al. 2020;Vincentelli et al. 2021;Kara et al. 2021;Montano et al. 2022), which is denoted as the local sample. Most of these AGNs are extensively monitored at multi-wavelengths ranging from X-ray to near-infrared, providing the best constraints on the R CER to date. Their R CER at rest-frame 2500Å are summarized in the Table 3 of G22. Finally, our CRM sample consists of 94+12 =106 objects.
In total there are 21 objects successfully cross-matched, among which there are 4, 5, and 12 objects from the core, parent sample, and local sample, respectively. For each object, we collect the Hβ BLR size, optical luminosity λL λ (5100Å). Note that if there are multiple measurements of individual object, we use the average BLR size and luminosity provided by Du & Wang (2019). The luminosity is transferred into the same cosmology if necessary. Table 1 lists the properties for all the crossmatched objects.

RESULTS
The upper panels in Figure 1 compare the BLR and CER R-L relation. The BLR sizes are characterized by the most well-studied emission line Hβ, except for NGC 4395 where only Hα is available. In this section, we will investigate the connection including and excluding NGC 4395, respectively. The R CER is characterized at rest-frame 5100Å which is transferred from 2500Å assuming the lag-wavelength power-law index β = 4/3 (i.e., R CER,5100 = R CER,2500 (5100/2500) On the other hand, the R CER -L relation has a slope of 0.48±0.04, established using 49 objects (see G22, 38 objects from the core sample, and 12 objects from the local sample). The two slopes are fully consistent with each other but the CER size is substantially smaller. Another interesting fact is that both the R BLR and R CER of NGC 4395 show substantial offset to the best-fit relation of more luminous AGNs.
The similar dependency of R CER and R BLR on luminosity can be interpreted as a significant contribution from the BLR DC in the optical continuum band (e.g., Li et al. 2021;Netzer 2022, G22). If this scenario is the real case, this similarity provides an opportunity to use the CRM derived R CER as a surrogate of the spectroscopic RM measured R BLR to estimate the BH mass.
To study this feasibility, we plot the direct comparison between R BLR and R CER in the lower panels in Figure  1 and fit the relation using the following equation: log (R BLR /light-day) = α log (R CER /R CER,0 ) + K (2) where α and K are the slope and intercept, respectively. R CER,0 is the reference point, which is set to logR CER,0 = 0.5, close to the median of our crossmatched sample. We perform linear regression using emcee (Foreman-Mackey et al. 2013) with the likelihood expressed by: where y i is the ith observed logR BLR , m i is the model prediction based on ith logR CER through equation 2, and s 2 i = (αx err,i ) 2 + y 2 err,i + σ 2 int involves the uncertainties of R CER , R BLR and the intrinsic scatter σ int . We perform the fitting with and without NGC 4395, respectively. The best-fit slopes and intercepts as well as their uncertainties are summarized in Table 2.
As shown by the lower panels of Figure 1, the correlation between R BLR and R CER has a slope of α = 0.984 +0.101 −0.101 if NGC 4395 is included. Otherwise, it shows a slightly flatter slope of α = 0.759 +0.150 −0.157 , which is still consistent with linear relation in normal space (α = 1.0) within 2σ uncertainty. We adopt the fitting result with NGC 4395 as the fiducial relation because he inclusion of NGC 4395 provides much larger dynamical range thus better constraints of the slope.
The best-fit intercept of the fiducial case is K = 1.410 +0.501 −0.348 at R CER,0 = 0.5, which reveals that the R BLR is 8.1 times (0.910 dex) larger than the R CER . In the case of excluding NGC 4395, the result is 8.6 times (0.934 dex). These results are consistent with the predicted value from the radiation pressure confined (RPC) cloud model (Netzer 2022) where the R BLR /R CER ∼ 8.7 (0.939 dex) assuming that continuum lags are fully contributed by the DC component and the BLR covers 20% percent of the sky of the central disk. However, the broad-line RM usually assumes that the continuum is emitted from a very compact region, whose size is considered negligible when compared to R BLR . However, if the R CER is actually ∼1/8 of the size of Hβ emitting region, it implies that the previous RM may have underestimated the BH mass by an average of 12.5% using 5100Å continuum. This underestimation can be even larger if the continuum is close to the Balmer limit, e.g., using B or g band as the reference continuum for objects at redshift z ∼ 0.2 to z ∼ 0.6.
The intrinsic scatter of the fiducial case is σ int = 0.284 dex, comparable to the current scatter of R BLR -L relation (∼0.22 dex Dalla Bontà et al. 2020;Malik et al. 2023). In the case of excluding NGC 4395, the intrinsic scatter is even lower (σ int = 0.265 dex). This scatter can be partly resulted from different contributions of BLR DC relative to the disk component among different objects if we assume DC is the main reason leading to this relation. The AGN variability between the BLR RM and   Notes. These two fitting results represent the relation shown in the lower left and right panel of Figure 1, respectively. We adopt the fitting result with NGC 4395 as the fiducial relation.
CRM campaigns can also introduce a large scatter. In our sample, all R CER were measured within recent ∼ 7 years, while the R BLR of 7 out of 21 objects include measurements more than 20 years ago (e.g. Peterson et al. 1998a;Kaspi et al. 2000). However, a direct comparison between the luminosity used in CRM studies and that reported by Hβ RM campaign shows no significant offset, indicating no dramatic variability for objects in our sample. Other factors can also contribute to the scatter. For example, G22 measures the R CER using three optical bands but the local sample uses more bands, especially in the UV wavelengths. A simple combination can lead to extra scatter.
Based on the R BLR -R CER correlation, we can estimate the R BLR from CRM, i.e., R BLR,CRM . Using the fiducial relation, we predict the R BLR for the rest objects in the core and parent sample without BLR size measurements. As shown in Figure 2, the objects in the core sample (filled and empty diamonds) perfectly follow the trends of the existing Hβ R BLR -L relation, exhibiting a similar scatter as direct Hβ RM measurements (0.22 dex, Dalla Bontà et al. 2020;Malik et al. 2023). On the other hand, some objects in the parent sample (filled and empty squares) show smaller BLR size than the prediction from the canonical R BLR -L relation (Dalla Bontà et al. 2020), which is inherited from G22. It can be explained by the fact that the parent sample in G22 includes some small lag measurements that are consistent with zero within 1σ (less reliable). In addition, the generally lower fractional variability (e.g., weaker variability with relatively larger photometric noise) of the parent sample can lead to an underestimation of the lag as found in the two-night observation of NGC 4395 (Montano et al. 2022).
Next, we estimate the CRM BH mass combining the R BLR,CRM and the FWHM from SE spectra in the literature (e.g., Shangguan et al. 2018;Du & Wang 2019;Liu et al. 2019). We compare the CRM BH mass with RM/SE mass using Hβ in the right panel of Figure 2 and find that they are in good agreement with each other. The median and 1σ scatter of the ratio log (M BH,CRM /M BH,SE ) is −0.03 and 0.23 dex for core+local sample, respectively, confirming the feasibility of CRM BH mass estimation. Compared with the ∼0.28 dex intrinsic scatter of R BLR -R CER relation, the ∼0.4 dex intrinsic scatter of virial factor f (Woo et al. 2010(Woo et al. , 2013Ho & Kim 2014) is still the primary source of the CRM BH mass uncertainty, which includes the intrinsic scatter of M BH -σ * relation (Ferrarese & Merritt 2000;Gebhardt et al. 2000) as well as the diversity of BLR orientation and dynamic structures. The f factor varies between different bulge types (Ho & Kim 2014) and is suggested to correlate with BLR orientation (Mejía-Restrepo et al. 2018;Yu et al. 2019). The use of SE line width can introduce additional uncertainties due to asynchronous measurements of line width and the R CER as well as the difference between SE/mean spectrum and rms spectrum (Peterson et al. 2004;Collin et al. 2006;Wang et al. 2019;Yu et al. 2019). Without determining f factor for individual object as done in dynamical modelling using spectroscopic RM data (Pancoast et al. 2014;Williams et al. 2018;Villafaña et al. 2022), the CRM BH mass uncertainty is not likely to be lower than 0.4 dex.
Next, we explore the BH mass of NGC 4395. If adopting the fiducial relation, the estimated R BLR of NGC 4395 from CRM is 166 ± 6 minutes with a 0.284 dex intrinsic scatter, and the CRM BH mass is (3.5 ± 0.4)×10 4 M ⊙ with a ∼0.5 dex intrinsic scatter. This BH mass is more consistent with the RM BH mass (Woo et al. 2019a;Cho et al. 2021) relative to the dynamical BH mass (den Brok et al. 2015).
Finally, we investigate the affection of accretion properties on the relation, calculating the Eddington ratios (λ Edd ) as well as the dimensionless accretion rates (Ṁ , Wang et al. 2014;Du et al. 2015) of our sample. None objects in our sample exhibit an Eddington ratio λ Edd > 1 (see Table 1) while if usingṀ ∼ 3 as the criterion to distinguish sub/super-Eddington accretion, six objects belong to super-Eddington accretion category representing roughly 30% of our sample. We find no apparent difference between these two categories in the R BLR -R CER relation within our sample. However, due to the limited sample size and quality of lag measurement, the dependence on accretion properties need to be studied more thoroughly with future observations of a larger sample. • This approach provides a promising opportunity to measure BH mass for high redshift AGNs since CRM requires a much shorter monitoring baseline. For example, the g-i lags of the 9 AGNs (logL 5100 ∼44.5 erg s −1 , z∼0.1) from core and parent sample range from 3 to 13 days in the observedframe, roughly 1/5 to 1/14 of Hβ lags. An AGN at z = 2 with similar luminosity will have an expected observed-frame continuum lag of less than 40 days. Thus, a 120-day baseline will be sufficient to recover its continuum lag, much shorter than the requirement of broad emission-line RM (e.g., Grier et al. 2019). It is also possible to measure the CRM BH mass of an AGN at z = 6 using near-infrared photometric monitoring if the variability is not fully smoothed by time dilation. As estimated, the inter-band lag between J and K is ∼ 80 days in the observed frame for a quasar with logL bol ∼47.0 erg s −1 , which can be recovered with a ≳ 180 days monitoring using ≳ 6m telescope.
• It also allows small telescopes (aperture size 1 ∼ 2 m) to measure the mass of small BHs. For example, NGC 4395 contains an intermediate-mass BH, and exhibits low continuum luminosity and very weak broad emission lines. It needs 8m-level telescopes for spectroscopic RM campaign (Cho et al. 2021), while its CRM BH mass is still accessible using a 2m telescope as in the case of the successful measurement of minute-level continuum lags by Montano et al. (2022). Such analysis can be extended to a large variability selected IMBH sample (e.g., Baldassare et al. 2020;Ward et al. 2022;Burke et al. 2022;Shin et al. 2022;Treiber et al. 2022) with well-sampled light curves. In our work, we find that the ratio R BLR /R CER of NGC 4395 seems consistent with higher luminosity AGNs, which further paves the way to extend the CRM method to the IMBH regime.
• The BH mass obtained from the R BLR -R CER relation can be potentially less biased than the SE BH mass based on the R BLR -L relation which is affected by the accretion rate and the host light. On the other hand, the R CER is a direct size measurement and both the diffuse continuum and broad emission-line originate from the BLR sharing similar kinematics. The tight R BLR -R CER relation presented in this work demonstrates the feasibility of direct BH mass estimation through CRM despite using a SE line width. Moreover, the intrinsic scatter of the relation can be further reduced with simultaneous R BLR and R CER measurements. These expectations can be tested with a larger sample.
On the other hand, CRM BH mass also has several caveats. As noticed in §3, the predicted R BLR of some objects in the parent sample shows large offsets. It may suggest that the current quality of the light curves, i.e., cadence and number of bands, are not enough to accurately determine the R CER,BLR . In addition, the DC contribution to the continuum lag depends on the continuum luminosity (G22). This correlation probably caused by the Baldwin effect of the DC  will also systematically affect the observed R CER -L relation. The study of such secondary dependence on different AGN properties is beyond the scope of this work. Apart from these, it is known that some objects do not show clear u/U band excess (Kara et al. 2021;McHardy et al. 2023), indicating a small diffuse continuum contribution to the continuum lags. At this point, its influence on our R CER -L relation is still unclear as separating the DC component from the disk is not easy (Guo H. et al. 2022b).
The next-generation time-domain survey, i.e., LSST, will help solve some of these problems and provide more accurate BH mass estimation. Its main survey, the Wide-Fast-Deep survey (WFD), will have six filters (ugrizy) covering from optical to near-infrared wavelength range. Any footprint in its ∼18000 square degree sky coverage will be observed ∼ 1000 times in the baseline of 10 years. The single visit can reach a depth of 24.5 mag in r band, and the expected number of AGNs monitored by WFD will be 6.2 × 10 6 (De Cicco et al. 2021). In addition to WFD, LSST will allocate ∼ 10 % of its observing time to five Deep-Drilling-Field (DDF) for a high-cadence (∼ 14000 visits) deep monitoring survey down to a coadded depth of ugrizy deeper than 26.5 ∼ 28.5 mag (Brandt et al. 2018). The accuracy of the lags can be 5% or 15% with a 2 or 5-day cadence, respectively (Pozo Nuñez et al. 2023). As studied in Kovačević et al. (2022), if the observational strategy for these fields has a cadence of 1 ∼ 5 days and duration > 9 years, the expected number of sources, whose continuum lags can be successfully measured, is >1000 for each DDF (10 square degree) in any filter. Therefore, the DDFs will become ideal laboratories to measure CRM BH mass as well as investigate the possibility of disentangling the R BLR (Czerny et al. 2023) and R CER through the photometric light curves.
Finally, it needs to mention that the CRM approach needs a SE spectrum to provide the velocity information. Using non-simultaneous observation of spectrum and CRM will provide additional uncertainties of the BH mass. Therefore, large-area spectroscopic surveys that will accompany the LSST, e.g., the 4-meter Multi-Object Spectroscopic Telescope (4MOST) (de Jong et al. 2019) and the Manuakea Spectroscopic Explorer (MSE) (The MSE Science Team et al. 2019), will benefit our method and provide more accurate BH mass estimations.

CONCLUSIONS
In this Letter we present a new BH mass estimator via optical continuum RM. Based on the similarity between the BLR and the CER size-luminosity relations (Netzer 2022, G22), we find a tight scaling relation between the size of BLR and CER using a sample of 21 AGNs. We find that the BLR size is about 8.1 larger than that of CERs at 5100Å (Figure 1), consistent with the model prediction from Netzer (2022). The intrinsic scatter of this relation is about 0.28 dex, comparable to the Hβ R BLR − L relation (0.22 dex). We apply this relation to the objects with reliable CER size measurements but without BLR size measurements to estimate the BLR size and BH mass (Figure 2). We find that the predicted R BLR of these objects follow the existing Hβ R BLR -L relation well, and the estimated BH mass is also consistent with the RM/SE BH mass using Hβ. Our proposed continuum RM BH mass approach will play an important role in the era of LSST. With the advent of more accurate and high-cadence light curves, the new continuum R-L relation will allow us to estimate the BH mass for a large sample of AGNs, especially for the IMBHs and high-redshift AGNs.