On the Performance of Quasar Reverberation Mapping in the Era of TimeDomain Photometric Surveys
Abstract
We quantitatively assess, by means of comprehensive numerical simulations, the ability of broadband photometric surveys to recover the broad emission line region (BLR) size in quasars under various observing conditions and for a wide range of object properties. Focusing on the general characteristics of the Large Synoptic Survey Telescope (LSST), we find that the slope of the sizeluminosity relation for the BLR in quasars can be determined with unprecedented accuracy, of order a few percent, over a broad luminosity range and out to . In particular, major emission lines for which the BLR size can be reliably measured with LSST include H, Mg ii , C iii] , C iv , and Ly, amounting to a total of timedelay measurements for all transitions. Combined with an estimate for the emission line velocity dispersion, upcoming photometric surveys will facilitate the estimation of black hole masses in AGN over a broad range of luminosities and redshifts, allow for refined calibrations of BLR sizeluminosityredshift relations in different transitions, as well as lead to more reliable crosscalibration with other black hole mass estimation techniques.
Subject headings:
galaxies: active — methods: data analysis — quasars: emission lines — quasars: supermassive black holes — techniques: photometric1. Introduction
Mass estimation of supermassive black holes (SMBHs) in active galactic nuclei (AGN) relies on locallyestablished relations between the SMBH mass and various source observables, such as the luminosity and the velocity dispersion of the broad emission line region, BLR (Kaspi et al., 2000; Peterson et al., 2004; Vestergaard & Peterson, 2006; Bentz et al., 2009; Denney et al., 2009). While the adequacy of such relations has been demonstrated for local samples of objects where several independent means for estimating BH masses exist (e.g., via the reverberation mapping technique and the stellar velocity dispersion in the inner regions of the host; Onken et al. 2004), it is not clear that this approach is warranted also at higher , where one probes earlier cosmic times, and is sensitive only to the most luminous sources at those epochs (Netzer, 2003). Further, as quasar^{1}^{1}1The terms AGN and quasars are used here interchangeably. observables in a given spectral band are redshiftdependent, it is not clear that intercalibration of various phenomenological relations (e.g., between luminosity, the particular emission line probed, and SMBH mass), even if justified, is free of biases (Denney et al., 2009). In this respect, the challenge in SMBH mass estimates at high is akin to the cosmological distance scale problem, with more data and better control of systematics required to place them on firmer ground.
With upcoming (photometric) surveys that will monitor a fair fraction of the sky to unprecedented depth and photometric accuracy, and with good cadence, the timedomain field of quasar study, and specifically reverberation mapping of the BLR, is expected to undergo a major revision (Chelouche & Daniel, 2012). As has been demonstrated in several previous works (Haas et al., 2011; Chelouche & Daniel, 2012; Chelouche et al., 2012; Chelouche, 2013; Chelouche & Zucker, 2013; Edri et al., 2012; Pozo Nuñez et al., 2012, 2013; Zu et al., 2013), while individual emission lines cannot be resolved using photometric means, their signal can be recovered at the light curve level (see also Rafter et al., 2013), and reverberation mapping is, in principle, possible.
In this work we build upon the formalism of Chelouche & Zucker (2013), and provide more realistic benchmarks for reverberation mapping using broadband photometric surveys, having in mind the characteristics of the Large Synoptic Survey Telescope (LSST)^{2}^{2}2http://www.lsst.org/lsst/scibook experiment. Specifically, we are interested in the ability of photometric campaigns in determining BLRassociated time delays in individual sources as well as statistically, in subsamples of sources. This paper is organized as follows: section 2 describes the model used here for constructing mock AGN light curves in different bands, as well as the analysis technique employed to recover the input linetocontinuum timedelays. Section 3 provides benchmarks for the measurement of the linetocontinuum timedelay under various assumptions concerning the sampling, redshift determination accuracy, filter choice, object properties, and the priors used in the analysis. The discussion, with particular emphasis on LSSTenabled science, follows in section 4, and the summary in section 5.
2. Light Curve Simulations and their Analysis
Below we specify our approach for simulating quasar light curves. For demonstrative purposes, we use the characteristics of the LSST and standard cosmology^{3}^{3}3.
2.1. Continuum light curves
We model the continuum light curve, , using the method of Timmer & Koenig (1995) so that its Fourier transform is of a powerlaw form, , where we choose to be the frequency powerlaw slope (note that our definition for differs from that used in, e.g., Giveon et al., 1999, who report instead slopes for the power spectrum). Our choice of results in a power spectrum that is somewhat redder than the one which characterizes continuous firstorder autoregressive [CAR(1)] processes, for which (Kelly et al., 2009). This is in better agreement with the slopes found by highcadence observations of AGN using Kepler data (Mushotzky et al., 2011), and may be more characteristic of the variability of high sources given the cadence and lifetime of the LSST experiment (Ivezic et al., 2008). We explore departures from our chosen value of in section 3.
2.1.1 The light curve variance
The normalized variability measure (i.e., the standard deviation of the light curve normalized by the mean flux level, and ignoring measurement uncertainties) in some filter , which is characterized by a central wavelength, , is , and depends, in principle, on the quasar luminosity, , its redshift, , and the wavelength range probed. Other quasar properties, such as radio loudness, black hole mass, and its Eddington ratio may also be important but are not taken into account in the present treatment of the problem. In particular, we do not consider blazars or sources where the jet contributes significantly to the variability.
Our model for is of the form
(1) 
where is the rest frame normalized variability, where is the restframe optical luminosity of the AGN, and are powerlaw indices whose values are discussed below, and are assumed to be independent of redshift and AGN characteristics.
Motivated by the fact that the variance of the light curve is (here we assumed rednoise spectra and ), we note that, unless the physical bound on is observable (e.g., when a spectral break lies in the frequency range probed by the time series), the reduced variability measure of the light curve depends on the duration of the experiment, , so that . For an experiment of a fixed duration over the entire sky (e.g., LSST) we get that because of this time dilation effect alone. Nevertheless, a further effect determines the final value of : observationally, it is well established that longer (rest) wavelength data have smaller variations such that their reduced variability measure hence (Garcia et al., 1999; Giveon et al., 1999; Meusinger et al., 2011). Therefore, observing quasars over a broad range of redshifts, using a particular band, would lead to an independent redshift dependence of the form . Combining the two sources for redshift dependence, we obtain that, for , . In this work we do not consider the possibility of quasar evolution, i.e., an intrinsic dependence of quasar variability on the cosmological epoch, an assumption which appears to be consistent with the data (Meusinger et al., 2011). If present, this could imply enhanced (suppressed) high variability of quasars than assumed here, thereby facilitating (restricting) reverberation mapping.
The luminosity dependence of the normalized variability measure is determined by searching for a reasonable agreement between the model predictions (equation 1), and the combined variability dataset for Seyfert galaxies (Sergeev et al., 2005) and PalomarGreen (PG) quasars (Giveon et al., 1999; Kaspi et al., 2000) at low . The model is then scaled up to higher using equation 1, and its predictions are checked against data for luminous high objects (Kaspi et al., 2007), which cover the relevant timescales. The normalization constant, and are then simultaneously determined yet we note that the solution may not be unique, and we make no attempt to cover the plausible solution phase space in this work. We find that and provide an adequate description for the median reduced variability measure of all data sets. Nevertheless, we caution that (a) the normalized variability measure for some Seyfert galaxies is somewhat underestimated due to host contamination in the Sergeev et al. (2005) dataset (Cackett et al., 2007), and that (b) some of the high objects in the Kaspi et al. (2007) sample are extremely radioloud objects, blazars, whose variability may not be characteristic of the bulk of the high radioquiet population.
Lastly, individual objects show considerable scatter around the typical (median) normalized variability measure at a given luminosity range. We find that a model where equation 1 is multiplied by a factor , where is a Gaussian random variable with a standard deviation of and a zero mean, qualitatively agrees with the observations (see Fig. 1). Overall, our model provides a fair statistical description of quasar variability, which is in line with the qualitative statistical framework considered in this work. However, it is not designed to account for the full richness of quasar variability (e.g., with a possible dependence also on the BH mass, or the Eddington ratio), nor for the particular properties of extreme AGN types (e.g., the narrow line objects), and the possibility of truly nonstationary light curve behavior.
2.1.2 Time delays
It is known that the continua light curves of quasars in different bands are highly correlated (Ulrich et al., 1997). Nevertheless, there are a few examples of lowluminosity Seyfert 1 galaxies showing discernible lags between the continuum light curves in different bands (Collier et al., 1998, 2001; Chelouche & Zucker, 2013; Chelouche, 2013). These are consistent with a physical picture in which the inner accretion disk irradiates its outer parts resulting in longer wavelength reprocessed emission lagging the short wavelength one (see however Korista & Goad 2001 for a different interpretation of the data). We note, however, that such time delays have not been established for highluminosity sources, nor in high quasars.
In cases where continuum transfer effects across the accretion flow are relevant, the light curve in some band , may be obtained from the light curve in band by means of a convolution with an appropriate transfer function, . For simplicity, we shall assume a rectangular function such that , and is defined in the range , where is the time delay between the bands (for more theoretically motivated transfer functions see Cackett et al. 2007, Chelouche 2013, and references therein). We use theoretically expected time delays from standard thin accretiondisk theory, which seem to be supported by the current data, such that
(2) 
where we implicitly assumed that , and normalized according to Chelouche (2013). The redshift dependence results from the combined effect of cosmological time dilation, and the smaller sizes of the accretion disk probed by shorter wavelength emission that enters the optical bands at higher . The continuum light curve in band is then set once is determined by
(3) 
where ’’ denotes convolution. If Shakura & Sunyaev (1973) disks provide a fair description of the accretion physics in AGN, then this relation is expected to hold except for objects with the lowest/highest accretion rates. In those physical regimes, advection becomes relevant (Abramowicz et al., 1988; Narayan & Yi, 1994).
2.2. Line and nonIonizing Continuum Emission
ID  days  

(1) Ly  941  4.7  42 
(2) C iii  977  6.6  50 
(3) Ly  1033  9.8  42 
(4) Fe iii  1117  3.7  50 
(5) Ly  1216  93  42 
(6) N v  1240  1  18 
(7) Si iv  1397  8  50 
(8) C iv  1549  24  40 
(9) C iii  1909  21  120 
(10) Fe iii  2077  2.5  50 
(11) Fe ii  2324  3.6  50 
(12) Mg ii  2799  32  70 
(13) Fe ii  2964  5  50 
(14) He i  3189  1  100 
(15) Fe ii  3498  1.4  50 
(16) H  4103  5  36 
(17) Fe ii  4160  1  50 
(18) H  4342  13  60 
(19) Fe ii  4564  20  200 
(20) He ii  4687  1  40 
(21) H  4863  46  100 
(22) Fe ii  5305  22  200 
(23) He i  5877  5  100 
(24) H  6565  195  130 
(25) He i  7067  3  100 
(26) O i  8457  10  200 
(27) Fe ii  9202  4  200 
(28) Pa  9545  7  150 
(29) Pa  10049  21  150 
(30) He i  10830  36  100 
(31) Pa  10941  7  150 
(32) Balmer cont.  34004000  100  30 
(33) Paschen cont.  6500 8300  300  150 
(34) Hot dust  1000020000  10000  500 
BLR emission features (i.e., lines, line blends, and continua) included in this work, with their parameters ( and ) given in the restframe; see table 2 in Vanden Berk et al. (2001) for wavelength definitions. Quoted lags correspond to an quasar (see text).
Approximate wavelength are used. For the recombination continua, allowance is made for the combined contribution of higherorder emission lines redward of the edge ionization energy.
The flux in this component is poorly constrained and is inspired by the theoretical calculations of Kwan & Krolik (1981); Korista & Goad (2001).
Once the wavelengthdependent (ionizing) continuum emission from the accretion disk has been defined for all bands, the contribution of emission lines and line blends, as well as variable nonionizing continuum emission from the outer regions (e.g., Balmer and Paschen continuum emission from the BLR, Korista & Goad 2001, and heated dust emission from its far outskirts, Netzer & Laor 1993) must be included. To this end, we compiled a list of major emission components, and estimated their relative contribution and time delays with respect to the adjacent continuum level from a mean quasar composite spectrum (Vanden Berk et al., 2001), and from various reverberation mapping campaigns (Kaspi et al., 2000; Peterson et al., 2004; Metzroth et al., 2006; Suganuma et al., 2006; Kaspi et al., 2007; Bentz et al., 2010; Barth et al., 2013; Chelouche et al., 2014)^{4}^{4}4We resort to educated guesses and theoretical works  e.g., Netzer & Laor (1993)  to estimate the lags in lines and blends for which reverberation results are currently unavailable.. We note that we consider luminosityindependent rest equivalent width values in this work, which is a clear simplification given the scatter observed (Shen et al., 2011) and the presence of the Baldwin effect, i.e., the anti correlation between line equivalent width and AGN luminosity (Baldwin, 1977; Bian et al., 2012). As our emission line properties are deduced from SDSS data for relatively bright objects (Vanden Berk et al., 2001), we expect the median equivalent width for LSST objects to be comparable to, or somewhat larger than, that assumed here^{5}^{5}5It turns out that the luminosity of the bulk of the quasars for which reverberation mapping would be feasible with LSST is , hence consistent with the SDSS composite (Vanden Berk et al., 2001).. Throughout this work, we assume that the time delays for a given emission line (or blend) scale as (Bentz et al., 2009; Chelouche et al., 2014, but see Kaspi et al. 2005), and that the variance in the line light curves is comparable to that in the continuum (Kaspi et al., 2000; Woo, 2008). Table 1 lists the various emission components used in this work, to which we refer to henceforth as emission lines.
The total AGN light curve in band is then given by:
(4) 
where the summation is over all entries in Table 1. is the line transfer function for the ’th emission feature, and is assumed to be of a rectangular form (see above), which starts at zero lag, and extends out to ( is the time delay). The telescope throughput curve for band , , is shown in Fig. 2 for the proposed set of LSST filters (including telescope efficiency and atmospheric effects). The rest equivalent width values for the relevant transitions, , are given in Table 1. We define the mean efficiency for an emission component centered on as,
(5) 
where is the profile of the emission component. For simplicity, we assume rectangular profiles with a width Å for the emission lines (this corresponds to velocities of , typical of broad line widths). For large scale continuum emission (e.g., the Balmer bump), we assume a rectangular profile which spans the wavelength range given in Table 1 (rows 32–34).
Examples for the light curves in different LSST bands for a AGN are shown in Figure 2. For clarity, we do not consider in this example the decreasing variance with wavelength (Eq. 1) and assume measurement uncertainties are negligible (those effects are included in the full calculations presented in sections 3, 4). Clearly, light curves in all bands are, to zeroth order, similar due to the dominant contribution of correlated continuum processes. Nevertheless, finite differences are noticeable: at progressively redder bands, the light curves are less spiky due to transfer effects across the accretion disk, which suppress small scale fluctuations (e.g., compare the light curve in the band and the driving continuum signal). In addition, the finite contribution of emission lines and blends to the broadband flux leads to longer trend deviations, and further suppresses small scale fluctuations. For example, the (normalized) band light curve lies above the driving continuum level at days, but below the driving continuum light curve at days, which is due to the lagging contribution of the H line to the flux. A more prominent effect is seen in the band, where a considerably delayed largeamplitude contribution from dust emission is apparent. Such differences between light curves in different bands can be detected in good quality data, allowing for the measurement of time delays associated with major emission lines (Chelouche & Daniel, 2012; Chelouche & Zucker, 2013).
2.3. Sampling and Measurement Noise
With LSST being the largest planned photometric survey in the coming decades, we consider the following observing patterns, which are inspired by the current LSST operations simulator (Opsim^{6}^{6}6https://www.lsstcorp.org/opsim/home): a) deep drilling fields (DDF) cadence, and b) universal sampling (UNS). For the DDF, daily contemporaneous acquisition of data in all bands is considered, but with two variants: with 120 days seasonal gaps (characterizing a fair fraction of the observable sky), or in a continuous viewing mode (e.g., for fields aimed toward the pole). For the UNS, we take the representative number of visits in each band over the sky, as given by the current Opsim, so that bands have visits over the entire tenyear duration of the survey. We note, however, that these are mere averages, and that some directions over the sky will likely benefit from better sampling while others will be characterized by sparser light curves. Four months seasonal gaps are assumed for all UNS sources but otherwise regular sampling is considered per band. We do not simulate data gaps due to downtime of the telescope, nor do we consider more realistic sampling patterns for the LSST, as those are, currently, publicly unavailable.
To estimate the noise level of a given observation, in each filter, we follow the results of the LSST exposure time calculator (ETC^{7}^{7}7http://dls.physics.ucdavis.edu/etc/), which provides approximate (average) values for the background noise, dark current noise, and readout noise for a 15 sec on (point) source exposure. The effect of varying background resulting from object declination and time of observation, as well as weather conditions, is not considered here. We assume Poisson noise associated with the target signal, and normalize according to the ETC for each band. The total measurement noise of an observation is obtained by summing all sources of noise in quadrature.
Representative lightcurves for a Seyfert 1 galaxy, and a quasar in the and bands are shown in figure 3 (we do not include the effect of intervening systems, such as the Ly forest, on the flux level of high sources). In the examples shown, the light curve of the low Seyfert galaxy is considerably noisier than the quasar’s, with the latter being characterized by a smaller variance. The effect of seasonal gaps is also evident.
2.3.1 Host Contribution
Finite aperture observations result in some of the host light contributing to the signal, and may lead to additional noise under varying seeing conditions^{8}^{8}8Note that a constant host contribution to the light curves is immaterial as reverberation mapping algorithms are insensitive to it (Welsh, 1999; Chelouche & Zucker, 2013).. For high sources, being the bulk of the AGN population in upcoming surveys (see Section 4), the contribution of the host decreases sharply because the surface brightness scales as and galaxies are restUV faint relative to the nuclear emission. To mitigate host contamination problems in faint low objects, largeaperture photometry or imagesubtraction techniques (Alard & Lupton, 1998) provide viable solutions. For the LSST, imagesubtraction techniques will be part of the standard reduction pipeline, and are expected to reduce varying seeing effects down to photon noise levels. Simulating secondorder host galaxy subtraction effects, which could affect a small fraction of the quasar population, is beyond the scope of the present work and will not alter the main conclusions of this paper.
2.4. Analysis
We apply the photometric reverberation mapping technique of Chelouche & Zucker (2013) to our mock photometric light curves. Briefly, one considers a pair of filters, one of which is relatively line rich, while the other one is line poor, and calculates the multivariate correlation function (MCF) to deduce a lag. In cases where the noncontinuum emission contribution to the linerich band is dominated by a single emission line, while the other (linepoor) band is devoid of emission components other than the primary ionizing continuum, the lag is the linetocontinuum time delay to a good approximation (Chelouche & Zucker, 2013). Denoting the linepoor (fluxed) light curve as , and the linerich (fluxed) light curve as , then, to firstorder approximation, the model for is
(6) 
where is the relative emission line contribution to the signal in the band, and is the linetocontinuum lag. One may deduce both and by searching for the peak of the MCF defined by
(7) 
where is the standard deviation of the time series, and is the number of points in the sum (see Chelouche & Zucker, 2013, and references therein for further details). Unlike Chelouche & Zucker (2013), we do not carry out Monte Carlo simulations on a casebycase basis to test for the significance of individual solutions since these would result in prohibitivelylong execution times on current hardware. As we shall show below (section 3.1), it is possible to effectively screen against unphysical solutions^{9}^{9}9The degree to which improved lag statistics for large samples of AGN may be obtained by considering the significance of individual measurements is worth exploring in future studies. This awaits better errorestimation algorithms for reverberation mapping, and further improvements in computer hardware..
While the MCF method has been shown to work well in regions of the spectrum with relatively sparse emission line density (e.g., for lowredshift objects in the redder bands), it is yet to be established whether it can also be applied at higher, where the optical wavelength bands are populated by many, possibly overlapping emission lines and blends, each of which is characterized by its own time delay and relative contribution to the flux. In such instances the aforementioned formalism is not strictly valid, yet could still serve as a fair approximation to the ideal experiment. The degree to which this statement applies is quantitatively explored in this work.
It also needs to be demonstrated that the MCF formalism of Chelouche & Zucker (2013) can work under more realistic observing conditions, which may be characterized by sparse, noncontemporaneous sampling in different bands. Further, redshift determination by photometric means alone has finite accuracy (Richards et al., 2001; Brescia et al., 2013, and note that the expected photometric redshift accuracy for LSST quasars is at the % level hence of the order of the typical emission line width^{10}^{10}10http://www.lsst.org/files/docs/sciencebook/SB_10.pdf), and so the choice of linerich and linepoor bands may be erroneous for certain redshift intervals, and for a fraction of all objects. The degree to which such issues affect linetocontinuum timedelay measurements is explored below.
2.4.1 Implementation
Given light curves in six filters per object, we first identify linepoor and linerich bands so that the above model (Eq. 6) is justified and the analysis meaningful.
Typical emission line contribution to the bands is of order 10% (Chelouche & Daniel, 2012), but can reach 30% for sources where the H and the Paschen bump contribution are dominant (Fig. 4). Line emission may even dominate the flux in the (relatively narrow) and bands at redshift intervals where the strongest hydrogen emission lines, Ly and H, are relevant.
While no band may be strictly free from BLR emission, we identify with the band characterized by the smallest contribution of BLR emission features to its flux (e.g., the band for sources using our emission model; see Fig. 4). It should be emphasized that the suitability of any given band to serve as depends not only on the relative emission line contribution to the flux, but also on the host galaxy contribution to it and the effect of seeing in nearby faint AGN, as well as on the number of visits, and how far in wavelength the band is removed from the linerich band in question (see more §3.4). Once is defined, timedelay measurements proceed via the calculation of the MCF with respect to each of the linericher bands, resulting in up to five timedelay measurements per object.
In this work we take the following approach to benchmark photometric reverberation mapping: the deduced time delay is compared to the (input) time delay of the emission component which contributes most to the linerich band. This is certainly a reasonable interpretation when one dominant emission line exists, and the degree to which it is generally useful is critically examined below (section 3.3).
3. Results
3.1. Single Emission Line
We first consider the simplest case, where a single emission line, H, is present in the spectrum, and attempt to recover its delay from photometric light curves that are characterized by DDF sampling in continuous viewing mode. For simplicity, all quasars are placed at (i.e., resulting in high signaltonoise light curves), and span a broad range of luminosity^{11}^{11}11We do not consider luminosityfunction statistics in this section (see §4).; in this case, is the band light curve, and we arbitrarily assign the band light curve to . Unless otherwise stated, of order objects were simulated, and the MCF was calculated for each, resulting in the timedelay statistics provided below.
Figure 5 shows the recovered radiusluminosity relation^{12}^{12}12Realistically, the BLR sizes for a given emission line in quasars of a given set of properties, such as luminosity, follow a physical distribution of some form rather than reflect on a single value. The resulting distribution functions for the recovered lags could therefore be somewhat broader than calculated here. Currently higher moments of the BLR sizeluminosity relations are poorly understood hence not included in our model. for H. On average, the slope of the lagluminosity relation is well traced over the intermediate luminosity range. In particular, the deduced slope of the sizeluminosity relation for objects with is , which is in excellent agreement with the input value^{13}^{13}13The uncertainty was determined by considering sub samples of the data, and obtaining fit statistics by means of bootstrapping.. At the lowluminosity end, the effects of sampling discretization are evident as timelags approach the sampling period, and the transfer function is under sampled. At the highluminosity end, clear deviations from the input sizeluminosity relation are evident as the survey’s lifetime (10 years for the first phase of the LSST) is comparable to the lag, and the full extent of the transfer function cannot be observed. In particular, the recovered lag distribution for objects is considerably skewed to unphysically short lags. As we discuss below, reliable timedelay measurements are possible provided additional information exists.
Considering the ratio between the output (i.e., recovered) and input quantities (O/I) for quasars of all luminosities, we see that most points cluster, as expected, around (Fig. 6). In particular, lag measurements that considerably deviate from the input lag are very likely to result in an offset solution for : e.g., lag solutions for the most luminous sources in figure 5 are skewed to shorter lags and to low values, which do not match the input values. Therefore, provided independent information on is available  for example, through spectroscopic observations, or from statistical knowledge of the quasar population as a whole  it is possible to screen against nonphysical solutions by rejecting highly discrepant values. Generally, allowing for smaller uncertainty margins on results in a smaller scatter in the lag O/I statistics, and also with fewer objects contributing to the sample (not shown).
As ”valid solutions” in this work, we consider values that deviate from the input value by no more than . This criterion seems to provide satisfactory results for a large range of object redshifts and luminosities (see section 4), and is consistent with the results of Shen et al. (2011) who show that % of all H rest equivalent width measurements for SDSS quasars lie within a similar interval. Valid solutions under our interval definition are relatively well clustered around the input values as can be seen in Figure 5. In particular, even for the most luminous sources, whose line transfer function is poorly sampled by the data, it is possible to identify reliable solutions in a fraction of all objects (note how the green points follow the input relation up to high luminosities in figure 5).
There is a statistically significant, and consistent, overestimation of the median lag with being underestimated by a similar factor (see the bulk of points in Fig. 6). This small but significant effect can only be detected with good enough statistics and reliable independent lag measurements, hence could not have been determined by Chelouche & Zucker (2013) in their analysis of Sergeev et al. (2005) data. The underlying reason for these offsets has to do with the MCF formalism being effectively blind to the structure of the line transfer function on short timescales. In particular, any emission line contribution to the light curve at short times (e.g., due to clouds which lie closer to our sightline and whose reprocessed radiation arrives nearly simultaneously with the continuum light curve), may be regarded as pure continuum emission. This is not surprising since line and continuum emission cannot be independently traced by photometric data. Considering narrower emissionline transfer functions, whose centroids are identical yet rise at finite times  i.e., the reverberating signal is sufficiently removed from zero lags  results in a reduced bias (Fig. 7; see also Zu et al. 2013 who use a narrow Gaussian transfer function). Similarly, the small offset in also vanishes, indicating its common origin; by not being able to trace part of the transfer function, the relative contribution of the line to the photometric light curve is effectively diminished. This in turn means that some information about the transfer function of emission lines is lurking also in broadband photometric data.
Lastly, we wish to study the effect of the powerdensity spectral shape, , on the ability of the MCF algorithm to reliably recover time delays. In particular, a considerable range of powerlaw indices have been found to characterize quasar lightcurves, with slopes, (Giveon et al., 1999; Mushotzky et al., 2011). To this end we consider a narrow range of intermediate luminosity objects (so that sampling issues are irrelevant), and plot the statistics of the recovered signal in Figure 8. As shown, the range covered by most AGN lies around a sweet spot in the parameter space, where the scatter in individual time delays around the median value is minimal, and the fraction of objects for which reliable timelag solutions are obtained peaks (as before, we use the constraint that the recovered deviates by less than % from the input value). For objects whose light curves show considerable variability only at the longest timescales (i.e., ), the recovered lags show the greatest dispersion around the median lag: while the latter quantity is robust, individual lag measurements are less reliable. Similarly, the fraction of unphysical solutions is somewhat increased for small values. At the other extreme, where AGN light curves gradually resemble white noise () , the fraction of objects with physical solutions drops significantly, and the dispersion in the lag increases so that lag measurements become unreliable. Moreover, even statistical averages become meaningless as the algorithm prefers to fit white noise with white noise, and the median time delay goes to zero. The slight increase in the fraction of ”reliable” physical solutions that satisfy our criterion for (note the last bin in Fig. 8) is purely artificial and has to do with the fact that the recovered in this limit is drawn from a uniform distribution in the range . We note that the aforementioned (transferfunction dependent) offset in the median lag depends relatively little on for the parameter range typical of quasars (Giveon et al., 1999; Mushotzky et al., 2011).
3.2. Sampling
Realistic experiments rarely benefit from strictly uniform sampling since seasonal gaps, telescope down times, bad weather, and scheduling constraints are all limiting factors. The effect of different sampling patterns on the sizeluminosity relation for H (again, neglecting other lines), are shown in Figure 9. Generally, with more timescales characterizing the observing pattern, the noisier the measurement becomes. In particular, introducing 120 day seasonal gaps results in noisier lag measurements around those timescales, which act as effective solution attractors. Interestingly, there is a tendency for the algorithm to somewhat overestimate the lag around the gap timescale (Fig. 9). Again, erroneous lag measurements can be discarded using additional information, if available, on : screening against discrepant solutions, some % of all measurements are discarded, and the remaining solutions match well the input sizeluminosity relation of the BLR. Specifically, the cluster of overestimated timedelay solutions around seasonal gap timescales is naturally excluded by this filtering scheme.
Experimenting instead with UNSlike sampling for the LSST, and restricting our analysis to light curves with visits in each band, we find that the ability to recover short delays (e.g., in faint enough AGN or in highionization emission lines) is compromised. In particular, the sampling period in this case is days, which is of the order of the minimal lag that can be reliably recovered. As before, using prior knowledge on reduces the scatter, especially at the lowluminosity end of the quasar population, which is characterized by the shortest lags, and allows for a better characterization of the sizeluminosity relation.
Lastly, we consider UNS, with one of the bands having visits, and the other band a mere visits^{14}^{14}14We do not find significant differences between cases in which the sampling is interchanged between the linerich and linepoor bands. (e.g., the and filters given the current version of the Opsim). In this case, the scatter is considerably increased (note the highly broadened lag distributions in Figure 9), as is the minimal time delay that can be reliably recovered; e.g., measuring the H lag for lowluminosity AGN is impossible. Discarding nonphysical solutions using our current screening algorithm is effective over a relatively narrow luminosity range () but even then individual measurements contain little information, with only large ensembles being able to characterize the typical time delay up to a % bias. Our calculations show that further screening the sample for objects showing the greatest variance in their light curves may be useful to further narrow down the lag distributions yet it is not clear that such a selection would provide an unbiased view of the quasar population as a whole, hence it is not considered here. Additional tests for the MCF solution may prove beneficial (see Chelouche & Zucker, 2013, for one possible version of such a test) in further weeding out erroneous solutions. These, however, are beyond the scope of the present work due to their prohibitive computational demands (see also section 2.4).
Simulations show that more reliable solutions under sparse sampling may be obtained when the emission line contribution to the band is larger. This may occur in cases where effectively narrower photometric bands are considered or when stronger emission lines are probed (see Fig. 9 and section 4 discussing the results for the Ly line). In particular, for the MCF reduces to standard crosscorrelation techniques with their advantages and shortcomings (Netzer & Peterson, 1997; Welsh, 1999, and references therein).
We note that, in the limit of truly random sampling, infinite characteristic timescales of the time series exist, and the solution phase space is uniformly covered (not shown). In such cases it is practically impossible to deduce the true lag, regardless of whether additional information on exists.
3.3. Multiple Emission Lines
Here we include all emission lines and blends (as well as the Paschen and Balmer bumps, and dust emission) listed in Table 1, and examine the ability of the MCF algorithm to recover the linetocontinuum time delay. More specifically, we first identify the linepoor band and associate its light curve with . We then determine the (effective) and from the MCF for each filter combination. The solutions are then compared to the input lags and the relative contributions to the flux of the emission features that contribute the most to each of the linerich bands. That is, whether a particular emission component contributes solely to the flux in the linerich band, or whether its contribution is only slightly larger than that of the other emission components contributing to the same band, does not alter our interpretation of the signal.
Results for individual objects having a luminosity of , over the redshift range , assuming DDF sampling with seasonal gaps, are shown in Figure 10 for the different bands. For the chosen luminosity, the algorithm is able to determine a lag up to . Specifically, as the S/N decreases, the recovered timedelay distribution extends to shorter recovered lags whence becoming less reliable. Nonphysical solutions may be discarded, as before, using prior information on , leading to solutions which hover around the input value, with a typical scatter of 50% (see also upper panels of figure 11).
Interestingly, some filter combinations give rise to more biased results. For example, time delay measurements in the band for objects (here, the band is the linepoor band), show a factor overestimation of the lag (see Figs. 10,11). The primary reason for that has to do with the fact that, at such redshifts, Fe ii and H have comparable contributions to the flux in the band, with the former transition being characterized by longer timedelays (Barth et al., 2013; Chelouche et al., 2014, see our Table 1). We also note that, at , the contribution of the emission lines to the band is only slightly larger than that to the band (see Fig. 4) and caution is advised when interpreting the signal.
It is also interesting to note the increased scatter in the recovered lags at redshifts where the identity of the linepoor filter is replaced, or around redshift intervals where the linerich band gradually becomes line poor, hence the BLR signal weakens (e.g., the case of the band at ; Fig. 11).
3.3.1 Redshift uncertainties
Redshift inaccuracies may lead to an erroneous identification of the linepoor and linerich bands, and to poor estimation of the relative contribution of emission lines to the flux, hence to an improper interpretation of the signal. We note that typical redshift uncertainties in photometric surveys are of the order of a few percent, and incorporated such uncertainties in our simulations pipeline (a Gaussian distribution with a standard deviation, was assumed, and a hard lower bound of enforced).
Figure 11 (middle panels) portrays the effect of redshift uncertainties, where the standard deviation of the recovered lag distribution rapidly increases at bandspecific redshift ranges. As expected, such intervals correspond to redshifts where the identity of the most prominent emission line at any given band is changed over narrow redshift intervals, or in cases where the identity of the linepoor band is switched. Most importantly, while individual measurements are less reliable at such redshift intervals, the median of the recovered lag distribution does not show significant deviation from the ideal case where the redshift is precisely known. Clearly, provided systematic effects are well controlled, large number statistics wins.
Lastly, we mention a subtle effect having to do with the emission line contribution to the band across filter edges. Specifically, for relatively boxy filter throughput curves it is possible for only one of the extended line wings to contribute to the flux in the band. As has been shown in several previous works, the time delay associated with emission line wings could differ from that which characterizes the bulk of the emission line (Grier et al., 2013, and references therein). Such effects could lead to small biases in the median lag close to certain redshift intervals. However, under most circumstances pertaining to broadband data, the emission line contribution at the filter edge will be overwhelmed by that of another transition close to the peak sensitivity curve of the band. The full treatment of such cases is beyond the scope of the present work.
3.4. Continuum Time Delays
The presence of time delays between the continuum emission in different wavebands (equation 2) could considerably bias linetocontinuum timedelay measurements, as already discussed in Chelouche & Daniel (2012). Figure 11 (bottom panels) shows this effect when considering a quasar with an optical luminosity of over a broad redshift range.
We find that, for much of the parameter space covered by our simulations, the effect of continuum time delays on the deduced (linetocontinuum) lag is modest. Nevertheless, there are particular filter combinations and redshift intervals where the deduced (median) lag can considerably deviate from the input value, which cannot be corrected by using priors on . For example, the emissionline timedelays associated with the band appear to be highly biased (by an order of magnitude) compared to the input delay for sources. This is due to the fact that the band is used as the linepoor band at this redshift interval, and due to the continua time delay between the band and the band being days in our model, i.e., of the linetocontinuum time delay in this case (see also Chelouche & Daniel 2012). A similar bias, although somewhat smaller in amplitude, is noticeable when considering the and bands at .
There are two straightforward ways to mitigate the aforementioned complications: (1) to use a linepoor band that is as adjacent as possible to the linerich band probed, so that interband continuum time delays are minimized (equation 2), or (2) to use a more comprehensive formalism, which is able to handle situations where two lagging emission components are present in the problem. The full treatment of such situations is naturally more complicated and may also lead to degenerate results (Chelouche & Zucker, 2013).
Lastly, we mention the possibility of nonnegligible diffuse continuum emission from the BLR (Korista & Goad, 2001, see their Figs. 2, 4). While the quantitative treatment of this scenario is beyond the scope of the present work, and is subject to model uncertainties, at the basic level, the problem is akin to that of the continua timedelays. Specifically, one should seek linepoor and linerich filter combinations for which diffuse continuum emission from the BLR has similar characteristics (e.g., lag and relative contribution to the flux).
4. Discussion
We have shown that photometric surveys with quasiregular sampling in several bands, and with characteristics similar to the LSST, can be used to measure the linetocontinuum timedelay in major emission lines, line blends, and nonionizing continua from the BLR. Lag measurements show enhanced scatter when the time delay is a fair fraction of the survey lifetime, or when the lag is comparable to the sampling period. Nevertheless, even in such cases, it is possible to reduce the scatter and reach statistically robust results if additional constraints on the value of (e.g., from independent spectroscopic measurements, or from knowledge of the quasar population as a whole) are incorporated in the MCF analysis. Probing linetocontinuum timedelays on subsampling timescales is unreliable, as photometry cannot disentangle line and continuum light curves on short timescales; an exception to this rule is in cases where the lagging signal dominates the flux in the band (Chelouche & Zucker, 2013; Chelouche, 2013), and provided the light curves are well sampled.
Time series whose sampling pattern exhibits many different characteristic timescales lead to considerable scatter in individual lag measurements, especially on timescales comparable to the sampling periods and their harmonics. Here too, it is possible to reject nonphysical solutions by employing priors on so that a statistical measurement of the median lag, in properly defined quasar samples, provides a good estimator for the true lag. The degree to which such filtering can be effectively used depends on the strength of the emission feature probed. Generally, as with spectroscopic reverberation mapping campaigns, regular sampling should be sought as one attempts to reduce the number of sampling timescales in the problem.
Good redshift determination leads to more robust lag measurement on a casebycase basis since preferable filter combinations may be reliably selected, and meaningful priors on may be set. This also exemplifies the advantage of having followup (singleepoch) spectroscopy which can secure the object identification and precisely determine its redshift. Moreover, spectroscopic followup can determine the velocity dispersion of the relevant emission lines, thereby providing critical information required for SMBH mass estimates.
We find that, even under ideal observing conditions, the recovered median lag may overestimate the true lag, i.e. the centroid of the line transfer function, by %. This is a direct consequence of the inability of photometric data to disentangle line and continuum light curves on short timescales, whose sum makes up the total signal. The exact value of the bias depends primarily on the line transfer function, which is rather loosely constrained by observations: for transfer functions with a large amplitude at zero time delays (e.g., when line emissivity close to our sightline is considerable, as in edgeon configurations) the bias is more significant. Consequently, the bias is smaller for faceon systems, which may be more relevant to typeI quasars (Maiolino et al., 2001). This implies that some statistical information concerning the line transfer function may be obtained using photometric means.
Our simulations indicate that photometric reverberation mapping is especially advantageous when large samples are concerned, as the median time delay is less susceptible to samplinginduced noise, redshift uncertainties, and the underlying properties of the powerdensity spectrum of the quasar. Therefore, good control of systematic effects is essential, especially as far as redshift determination is concerned. Provided large enough samples exist, statistical averaging would lead to BLR sizeluminosity relations with unprecedented accuracy for various emission lines, and over a broad luminosity and redshift ranges (see also Chelouche & Daniel 2012; Zu et al. 2013, and below).
It is worthwhile to consider more specific predictions for the LSST. We first consider the case where the signal per visit is obtained by combining the data of four consecutive (to within 30 minutes) exposures of 15 sec each, which should characterize the majority of LSST data^{15}^{15}15See http://www.lsst.org/files/docs/sciencebook/SB_2.pdf. Thus, we effectively neglect quasar variations on hour timescales, which is reasonable given the red powerdensity spectra of quasars, the suppressed variability of luminous high sources (Fig. 1), and the considerably longer timescales associated with the BLR.
Figure 12 divides the luminosityredshift space into small [] segments, each typically having simulated sources, and shows the statistical properties of the obtained lag solutions for Mg ii . The timedelay statistics includes objects for which the recovered is within % of the input value (section 3). The relevant parameter space for measuring the size of the Mg ii emission region is welldefined and forms an envelope that is bounded at lowluminosities by minimal S/N requirements, and by the finite duration of the survey at the highluminosity end. Within this envelope the recovered lag is within dex of the input value, with a median value of over the relevant phase space. In particular, several trends are observed: in regimes of effective lowS/N (either low source flux and/or small contribution of the emission line to the flux in the band), there is a tendency for a biased lag measurement, by up to 60%; c.f., the second column of Figure 11. There is also a tendency to somewhat underestimate the lag, by typically 30%, in cases where it is comparable to the lifetime of the experiment. Discarding those regions near the envelope’s rims, typical median lag determination is at the % level for the Mg ii line out to . Taking into account the quasar luminosity function and the LSST selection criteria, some lags may be determined over the LSST lifetime (see also Table 2) suggesting that Mg ii may have a crucial role in constructing more reliable singleepoch BH mass estimations by crosscalibrating scaling relations used at different redshift ranges.
Our calculations imply that the expected scatter in individual Mg ii lag measurements is, typically, %, with a tendency for an increased scatter under low S/N conditions (see the middle panels of Fig. 12). As noted above, these statistics were obtained after screening against erroneous solutions that, typically, results in only 10% of the objects being used (right panel of Fig. 12). Nevertheless, there are regions in the parameter space ( and ) that are characterized by good S/N, and wellsampled light curves with respect to the BLR extent, for which % of the sources lead to robust solutions. We note that it may be possible to further narrow the scatter in individual measurements by testing the robustness of the solutions (Chelouche & Zucker, 2013). This, however, is beyond the scope of the present work due to the prohibitively long computation time involved.
Reducing the exposure time per visit to 15 sec instead of 60 sec, shrinks the parameter space over which reliable Mg ii lag measurements are obtained while maintaining similar lag statistics in regions of the parameter space where such measurements are possible. Specifically, with the S/N reduced by a factor 2, Mg ii lag measurements are feasible up to redshift of , instead of , and only for the brighter sources, which reduces the number of reliable lag measurements by an order of magnitude for this transition (Table 2).
Figure 13 shows an analysis similar to the above for several additional major emission lines in quasar spectra, and examines some of their statistical properties using DDF and UNS samplings (60 sec exposure times are assumed throughout). Clearly, there exists a nonnegligible region in the redshiftluminosity plane in which time delays may be determined with good accuracy. For example, in lowluminosity, low quasars, the lagluminosity relation for H may be probed over orders of magnitude in luminosity with a total expected number of H lag measurements of (Table 2). For the C iv emission line, the parameter range probed extends to lower luminosities than probed by Kaspi et al. (2007), and is dominated by radioquiet objects rather than radio loud ones, as in their work. The small blue bump (Balmer continuum) has a similar role to the Mg ii line in the sense that its lag may be quantified over a broad range of quasar luminosity and redshift. Obviously, this spectral feature is of limited use for directly estimating BH masses as there is no kinematic information associated with it. Interestingly, we find that Ly lag measurements at the peak of quasar activity are facilitated by the effectively narrow and spiky throughput curve of the band, and are feasible out to .
UNS  DDF  DDF  

ID  (60 sec)  (60 sec)  (15 sec) 
Ly  64000  95000  7000 
Si iv  150  300  30 
C iv  8000  16000  1200 
C iii  21000  45000  6500 
Mg ii  46000  48000  4300 
H  2000  6000  300 
Fe ii  9400  16000  900 
H  47000  110000  19000 
Fe ii  6500  27000  4800 
He i  0  700  40 
H  28000  73000  25000 
Balmer cont.  310000  310000  50000 
Paschen cont.  14000  15000  4000 
Total for all BLR features:  560000  760000  120000 

Number of timedelays measured for individual emission features in quasars within the LSST footprint for DDF (using 15 sec and 60 sec exposure times) and UNS sampling (60 sec exposure time). Features for which lag measurements may be obtained are included in the table, and measurements for which the measured lag agrees to better than 0.5 dex with the input value. Typical uncertainties on the quoted figures are at the 10% level given the model assumptions (but may be of order unity for those transitions with detection statistics).
Full coverage of the LSST footprint by DDFs exceeds the LSST resources (see text).
Approximate wavelengths are quoted for blends (Table 1).
The fact that the lags of different emission lines may be determined in overlapping luminosityredshift ranges means that a more reliable BLRsize ladder (in analogy with the cosmological distance ladder) may be obtained using LSST, allowing the cross calibration of different prescriptions for SMBH mass estimation over cosmic time. In this respect, it is interesting to note that for an exponentiallysmall (but finite) fraction of all quasars, simultaneous lag determination for several transitions may be possible. The number statistics in this case is less secure and may depend on the exponential tails of currently poorlydetermined quasar property distributions (e.g., the fraction of objects with extreme variability amplitude). Current predictions for the statistics of multi timelag measurements within the LSST footprint for different levels of accuracy, and for UNS and DDF sampling (for the latter either 15 sec or 60 sec exposure times are considered) are shown in Figure 14. Clearly, an exponentially small (but finite) number of sources will have up to 5 timedelays simultaneously determined. The number of such sources is very sensitive to sampling: e.g., the number of quasars with five lag measurements drops by more than an order of magnitude when switching from DDF to UNS sampling. S/N has a smaller effect on the relative number of objects with multi timedelay measurements.
Lastly, we consider the total number of reliable timelag detections as a function of the number of DDFs covered assuming each field consumes 1% of the LSST resources^{16}^{16}16Assuming our definition of a DDF, each field, observed daily, will require four 15 sec exposures, each with an overhead of 2 sec of reading time, in six bands. Slewing between adjacent fields will take additional 5 sec. Together, this amounts to sec per DDF, which is of order 1% of a full observing night.. The total sky coverage in DDFs is then sqdeg, and the remaining UNS coverage is then sqdeg. The total number of time lag measurements as a function of is shown in the bottomright panel of Figure 14 demonstrating that sheer number statistics prefers UNS over DDF sampling with up to time delay measurements possible for full UNS coverage of the sky. Note, however, that the sparser UNS will result in the loss of short timedelay information concerning, e.g., accretion disks or highionization optical emission lines, such as He i 5877 (see Table 2), especially in low, lowluminosity sources.
5. Summary
We show that large timedomain photometric surveys, such as the LSST, are expected to transform the field of AGN reverberation mapping by increasing the number of objects with measured timedelays by several orders of magnitude and by including sources at the epoch of peak quasar activity, at . Specifically, the LSST is expected to yield a total of timedelay measurements including all major emission lines, blends, and pseudocontinuum features out to . For a given exposure time per visit and finite survey resources, covering a larger area of the sky at the expense of sampling would result in additional timelag measurements while inhibiting the study of short timescale phenomena. Spectroscopic followup of photometric quasars will lead to improved BLR lagdetermination via the incorporation of priors in the MCF analysis. It will also facilitate BH mass estimations thereby improving our understanding of black hole demography at high with implications for galaxy evolution, black hole seeds, and gravitational wave detection.
References
 Abramowicz et al. (1988) Abramowicz, M. A., Czerny, B., Lasota, J. P., & Szuszkiewicz, E. 1988, ApJ, 332, 646
 Alard & Lupton (1998) Alard, C., & Lupton, R. H. 1998, ApJ, 503, 325
 Baldwin (1977) Baldwin, J. A. 1977, ApJ, 214, 679
 Barth et al. (2013) Barth, A. J., Pancoast, A., Bennert, V. N., et al. 2013, ApJ, 769, 128
 Bentz et al. (2009) Bentz, M. C., Peterson, B. M., Netzer, H., Pogge, R. W., & Vestergaard, M. 2009, ApJ, 697, 160
 Bentz et al. (2010) Bentz, M. C., Walsh, J. L., Barth, A. J., et al. 2010, ApJ, 716, 993
 Bian et al. (2012) Bian, W.H., Fang, L.L., Huang, K.L., & Wang, J.M. 2012, MNRAS, 427, 2881
 Brescia et al. (2013) Brescia, M., Cavuoti, S., D’Abrusco, R., Longo, G., & Mercurio, A. 2013, ApJ, 772, 140
 Vanden Berk et al. (2001) Vanden Berk, D. E., et al. 2001, AJ, 122, 549
 Cackett et al. (2007) Cackett, E. M., Horne, K., & Winkler, H. 2007, MNRAS, 380, 669
 Chelouche (2013) Chelouche, D. 2013, ApJ, 772, 9
 Chelouche & Daniel (2012) Chelouche, D., & Daniel E. 2012, ApJ, 747, 62
 Chelouche et al. (2012) Chelouche, D., Daniel, E., & Kaspi, S. 2012, ApJ, 750, L43
 Chelouche et al. (2014) Chelouche, D., Rafter, S. E., Cotlier, G. I., Kaspi, S., & Barth, A. J. 2014, ApJ, 783, L34
 Chelouche & Zucker (2013) Chelouche, D., & Zucker S., 2013, ApJ, 769, 124
 Collier et al. (1998) Collier, S. J., Horne, K., Kaspi, S., et al. 1998, ApJ, 500, 162
 Collier et al. (2001) Collier, S., Crenshaw, D. M., Peterson, B. M., et al. 2001, ApJ, 561, 146
 Denney et al. (2009) Denney, K. D., Peterson, B. M., Dietrich, M., Vestergaard, M., & Bentz, M. C. 2009, ApJ, 692, 246
 Edri et al. (2012) Edri, H., Rafter, S. E., Chelouche, D., Kaspi, S., & Behar, E. 2012, ApJ, 756, 73
 Garcia et al. (1999) Garcia, A., Sodré, L., Jablonski, F. J., & Terlevich, R. J. 1999, MNRAS, 309, 803
 Giveon et al. (1999) Giveon, U., Maoz, D., Kaspi, S., Netzer, H., & Smith, P. S. 1999, MNRAS, 306, 637
 Glikman et al. (2006) Glikman, E., Helfand, D. J., & White, R. L. 2006, ApJ, 640, 579
 Grier et al. (2013) Grier, C. J., Peterson, B. M., Horne, K., et al. 2013, ApJ, 764, 47
 Haas et al. (2011) Haas, M., Chini, R., Ramolla, M., et al. 2011, A&A, 535, A73
 Ivezic et al. (2008) Ivezic, Z., Tyson, J. A., Acosta, E., et al. 2008, arXiv:0805.2366
 Kaspi et al. (2000) Kaspi, S., Smith, P. S., Netzer, H., et al. 2000, ApJ, 533, 631
 Kaspi et al. (2005) Kaspi, S., Maoz, D., Netzer, H., et al. 2005, ApJ, 629, 61
 Kaspi et al. (2007) Kaspi, S., Brandt, W. N., Maoz, D., et al. 2007, ApJ, 659, 997
 Kelly et al. (2009) Kelly, B. C., Bechtold, J., & Siemiginowska, A. 2009, ApJ, 698, 895
 Korista & Goad (2001) Korista, K. T., & Goad, M. R. 2001, ApJ, 553, 695
 Kwan & Krolik (1981) Kwan, J., & Krolik, J. H. 1981, ApJ, 250, 478
 Maiolino et al. (2001) Maiolino, R., Salvati, M., Marconi, A., & Antonucci, R. R. J. 2001, A&A, 375, 25
 Metzroth et al. (2006) Metzroth, K. G., Onken, C. A., & Peterson, B. M. 2006, ApJ, 647, 901
 Meusinger et al. (2011) Meusinger, H., Hinze, A., & de Hoon, A. 2011, A&A, 525, A37
 Mushotzky et al. (2011) Mushotzky, R. F., Edelson, R., Baumgartner, W., & Gandhi, P. 2011, ApJ, 743, L12
 Narayan & Yi (1994) Narayan, R., & Yi, I. 1994, ApJ, 428, L13
 Netzer (2003) Netzer, H. 2003, ApJ, 583, L5
 Netzer & Laor (1993) Netzer, H., & Laor, A. 1993, ApJ, 404, L51
 Netzer & Peterson (1997) Netzer, H., & Peterson, B. M. 1997, Astronomical Time Series, 218, 85
 Onken et al. (2004) Onken, C. A., Ferrarese, L., Merritt, D., et al. 2004, ApJ, 615, 645
 Peterson et al. (2004) Peterson, B. M., Ferrarese, L., Gilbert, K. M., et al. 2004, ApJ, 613, 682
 Pozo Nuñez et al. (2012) Pozo Nuñez, F., Ramolla, M., Westhues, C., et al. 2012, A&A, 545, A84
 Pozo Nuñez et al. (2013) Pozo Nuñez, F., Westhues, C., Ramolla, M., et al. 2013, A&A, 552, A1
 Rafter et al. (2013) Rafter, S. E., Kaspi, S., Chelouche, D., et al. 2013, ApJ, 773, 24
 Richards et al. (2001) Richards, G. T., Weinstein, M. A., Schneider, D. P., et al. 2001, AJ, 122, 1151
 Sergeev et al. (2005) Sergeev, S. G., Doroshenko, V. T., Golubinskiy, Y. V., Merkulova, N. I., & Sergeeva, E. A. 2005, ApJ, 622, 129
 Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
 Shen et al. (2011) Shen, Y., Richards, G. T., Strauss, M. A., et al. 2011, ApJS, 194, 45
 Suganuma et al. (2006) Suganuma, M., Yoshii, Y., Kobayashi, Y., et al. 2006, ApJ, 639, 46
 Timmer & Koenig (1995) Timmer, J., & Koenig, M. 1995, A&A, 300, 707
 Ulrich et al. (1997) Ulrich, M.H., Maraschi, L., & Urry, C. M. 1997, ARA&A, 35, 445
 Vestergaard & Peterson (2006) Vestergaard, M., & Peterson, B. M. 2006, ApJ, 641, 689
 Welsh (1999) Welsh, W. F. 1999, PASP, 111, 1347
 Woo (2008) Woo, J.H. 2008, AJ, 135, 1849
 Zu et al. (2013) Zu, Y., Kochanek, C. S., Kozłowski, S., & Peterson, B. M. 2013, arXiv:1310.6774