PHOTOMETRIC REVERBERATION MAPPING OF THE BROAD EMISSION LINE REGION IN QUASARS

Quasar light curves are reminiscent of red noise, following a power spectrum P(ω) ∼ ω^α (ω is the angular frequency) with α ≃ −2 (Giveon et al. 1999, and references therein). Here we model the pure continuum light curve in some band X as a sum of discrete Fourier components, each with a random phase ϕ_i:

$$\begin{equation} f_X^c(t)=\sum _{i=1}^\infty A(\omega _i){\rm sin}(\omega _it+\phi _i). \end{equation} \tag{ 1 }$$

The amplitude of each Fourier mode, A(ω_i)∝ω^α/2_i and t is time.⁴ We define the normalized variability measure, $\chi =\ssty\sqrt{\sigma _f^2-\delta ^2}/\left< f \right>$ (see Kaspi et al. 2000, who express this in percent), where <f > is the mean flux of the light curve, σ_f is its standard deviation, and δ is the mean measurement uncertainty of the light curve. For the PG sample of quasars (Schmidt & Green 1983), the continuum normalized variability measure χ_c ∼ 0.1–0.2 (Kaspi et al. 2000) and, typically, δ ≳ 0.01 (see Table 1). A light-curve realization, with δ ≪ σ_f and χ_c ∼ 0.2, is shown in Figure 1.

Zoom In Zoom Out Reset image size

Quasar emission occurs over a broad range of photon energies from radio to hard X-ray wavelengths and originates from spatially distinct regions in the quasar. In the optical bands, there is some evidence for red continua to lag behind blue continua (Krolik et al. 1991; Collier et al. 1998; Sergeev et al. 2005; Bachev 2009). Without loss of generality, we associate the Y band with the redder parts of the spectrum so that its instantaneous continuum flux may be expressed as a function of the flux in the X band:

$$\begin{equation} f_Y^c(t)= f_X^c * \psi ^c \equiv \int _{-\infty }^\infty d\tau f_X^c(\tau) \psi ^c(t-\tau), \end{equation} \tag{ 2 }$$

where ψ^c(δt) is the continuum transfer function.⁵ Presently, ψ^c is poorly constrained by observations and only the inter-band continuum time delay, t^c_lag, is loosely determined (Bachev 2009, and references therein). Nevertheless, so long as R_BLR/c≫ t^c_lag, the particular form of ψ^c is immaterial to our discussion (see Section 3.2.3). For simplicity, we shall therefore take ψ^c to be a Gaussian with a mean of t^c_lag and a standard deviation t^c_lag/2.

Emission by photoionized (BLR) gas is naturally sensitive to the continuum level of the ionizing source. In particular, for typical BLR densities and distances from the SMBH, the response time of the photoionized gas to continuum variations is much shorter than the light travel time across the region, and so may be considered instantaneous. Therefore, the flux in some line i, which is emitted by the BLR, can be written, in the linear approximation, as

$$\begin{equation} f_i^l= f_X^c * \psi _i^l, \end{equation} \tag{ 3 }$$

where the transfer function ψ^l_i depends, essentially, on the geometry of the BLR region as seen from the direction of the observer (Horne et al. 2004, and references therein). Our present knowledge of ψ^l(δt) for the quasar population as a whole is, however, rather limited (see Goad et al. 1993 for a few relevant forms for ψ^l). In what follows we shall consider two forms for the transfer function: (1) a Gaussian kernel with a mean t^l_lag and a standard deviation t^l_lag/4, and (2) a transfer function with a flat core in the range 0 ⩽ δt ⩽ 1.36t^l_lag, with a following linear decline to zero at δt = 2.72t^l_lag (such a transfer function also results in a time lag being ≃ t^l_lag). The latter transfer function is motivated by the observations of the Seyfert 1 galaxy NGC 4151 (Maoz et al. 1991) and corresponds to a case where a geometrically thick, spherical shell of broad emission line gas isotropically encompasses the continuum emitting region. The different transfer functions defined above will henceforth be termed (1) the Gaussian and (2) the thick-shell BLR models. The transfer functions are depicted in Figure 2. The main qualitative difference between the two transfer functions, as far as light curve behavior is concerned, is that while the Gaussian one is relatively localized in time, the thick-shell BLR model is much broader and results in longer range correlations between the light curves of the two bands. Unless otherwise stated, our results are presented for the Gaussian model. Cases for which we find a qualitative difference between the two transfer functions are explicitly treated (e.g., Section 3.2.5).

Zoom In Zoom Out Reset image size

A simulated light curve of an emission line that is driven by continuum variations (assuming t^l_lag = 200 days) is shown in Figure 1 (blue curves). Time-delay and smearing effects, with respect to the continuum light curve, are clearly visible. Here too, we define a normalized variability measure for the line, χ_l, and consider χ_l ≃ 0.75χ_c (Kaspi et al. 2000, see their Table 5). The sum of continuum and lines' emission to the total flux in band Y is therefore given by

$$\begin{equation} f_Y(t)= f_Y^c(t)+f_Y^l(t) = f_Y^c(t)+\sum _i f_i^l(t), \end{equation} \tag{ 4 }$$

where f^l_Y is the total flux in the emission lines. While a similar expression applies, in principle, also for the total flux in band X, unless otherwise stated, we shall assume that f_X = f^c_X, i.e., that emission line contribution to band X is either negligible or that the line variability associated with it operates on timescales that are much too short to be resolved by the cadence of the experiment, or too long to be detectable in the time series. For simplicity, we shall limit our discussion to a single emission line contributing to the Y band (this restriction is relaxed in Section 3.2.4).

The contribution of the emission line(s) to the total flux in the Y band requires knowledge of the equivalent width of the broad emission line, W (typically ≲ 100 Å in the optical band for strong lines; Vanden Berk et al. 2001), the transmission curve of the broadband filter, as well as the redshift of the quasar. For example, for filters with sharp features in their response function, even small redshift changes could lead to substantial differences in the relative contribution of the emission line to the total flux in the band, η. For a flat filter response function over a wavelength range W_b (typically, W_b ∼ 10³ Å for broadband filters), η = W/W_b. (A more accurate treatment of the line contribution to the broadband flux for the case of Johnson–Cousins filters is given in Section 4.) Example light curves for the X and Y bands, assuming η = 0.133, are shown in Figure 1. In the following calculations, we shall work with normalized light curves that have a zero mean and a unit variance, i.e., f → (f − <f>)/σ_f.

In the upcoming era of large surveys, statistical information may be gathered for many objects. For example, the LSST is expected to provide photometric light curves for >10⁶ quasars having a broad range of luminosities and redshifts. With such large sets of data at one's disposal, sub-samples of quasars may be considered having a prescribed set of properties (e.g., luminosity and redshift), and the line to continuum time-lag measurement may be quantified for the ensembles. In this case, it is possible to measure the characteristic line-to-continuum time delay in two ways: by measuring the peak of ξ_CA(δt) or ξ_AA(δt) on a case-by-case basis and then constructing time-lag distributions to measure their moments, or by defining the ensemble averages of the statistical estimators. Specifically, we define <ξ_CA(δt)>, <ξ_AA(δt)> as the arithmetic averages⁹ of the corresponding statistical estimators at every δt. (Alternative definitions for the mean do not seem to be particularly advantageous and are not considered here.) The average statistical estimators for ensembles of varying sizes are shown in Figure 4. As expected, averaging leads to a more robust measurement of the peak, as the various ξ_CA(δt) "constructively interfere" at ∼t^l_lag. The ensemble size required to reach a time-lag measurement of a given accuracy naturally depends on the noise level.

We find that measuring the time lag from the peak of <ξ_CA> (or <ξ_AA>) is more accurate than measuring the time lags for individual objects in the ensemble, and then quantifying the moments of the distribution (Figure 7).¹⁰ This is particularly important when noisy data are concerned, in which case reliable time-lag measurements for individual objects are difficult if not impossible to obtain. To see this, we note the rising uncertainty on the deduced time lag, σ(t^l_{lag, measured}), as the number of visits per object is reduced (and despite the fact that a correspondingly larger sample of objects is considered; see the right panel of Figure 7 for a case in which δ/<f> = 0.1). In contrast, when working with <ξ_CA>-statistics, the uncertainty monotonically declines so long as the sampling of the transfer function is adequate.

The above analysis highlights the importance of using large samples of objects (even with less than ideal sampling) in addition to studying individual objects having relatively well-sampled light curve. In particular, by using large enough samples of quasars, the effects of noise, and even sparse sampling, may be overcome.

Current studies suggest that there is a scatter of 30%–40% in the BLR size versus quasar luminosity relation (Bentz et al. 2009; Kaspi et al. 2000). To assess the effect of such a scatter on <ξ_CA> (similar results apply also for <ξ_AA>), we averaged the results for an ensemble of simulated objects whose time lags follow a Gaussian distribution with a mean t_lag and standard deviation σ(t_lag). We modify this Gaussian such that the input lag is always positive definite: for σ(t_lag)/t_lag ≪1 the distribution is Gaussian while for σ(t_lag)/t_lag ∼ 1 there is an excess of systems with zero lag. <ξ_CA(δt)> is shown in Figure 8: evidently, the average statistical estimator peaks at the mean of the distribution so long as σ(t_lag)/t_lag <0.8. For larger ratios, the distribution considerably deviates from Gaussian and its mean is>t_lag. In this case, <ξ_CA(δt)> still peaks around the mean but its shape is highly asymmetric. This is partly due to the time-lag distribution, but also so due the fact that the transfer functions with the largest t^l_lag are not adequately sampled by the light curve (t_tot/t^l_lag = 4 was assumed). To conclude, photometric reverberation should yield accurate mean BLR sizes for a sample of quasars, given the observed scatter in the BLR properties.

Zoom In Zoom Out Reset image size

Biases are present in the time-lag measurement from noisy data. These are negligible for noise levels δ/<f> <η but may surface otherwise. Specifically, for a (Gaussian) noise level of ∼50%, and a relative line contribution to the band of 10%, both photometric estimators overestimate t^l_lag by as much as 20% with ξ_AA being more biased than ξ_CA (Figure 9). This bias could therefore be present in situations where large and noisy ensembles are considered, or when weak emission lines are targeted. Quantifying the bias requires good understanding of the noise properties.

Zoom In Zoom Out Reset image size

A related bias exists when the noise level in the two bands differs substantially. More specifically, when |δ_X/<f_X> −δ_Y/<f_Y> | ∼ η (δ_X/δ_Y is the noise level in the X/Y bands), our simulations show that t_{lag, measured} may differ from t_lag by as much as 25%. Whether an overestimate or an underestimate of t_lag is obtained depends on the particular noise properties, as mentioned above. Observationally, one should therefore aim to work with data having comparable noise levels for the different bands to obtain a bias-free result. As noted above, for relatively strong lines or low noise levels, such biases are negligible.

Varying seeing conditions leading to aperture losses, as well as improperly accounting for stellar variability in the field of the quasar, could lead to correlated errors between the bands. This in turn could lead to artificial zero-lag peaks in the correlation function of the light curves (Welsh 1999), and several solutions have been devised to overcome this problem (Alexander 1997; Zu et al. 2011, and references therein). Nevertheless, the formalism developed here subtracts off, by construction, the contribution of similar processes (e.g., correlated noise) to both light curves. In fact, as far as the algorithm is concerned, the continuum contribution to the Y band is merely an effective correlated noise to be corrected for, which it does. Therefore, photometric reverberation mapping, as proposed here, is considerably less susceptible to the effects of correlated noise than the spectroscopic method. That being said, if the noise between the bands correlates over sufficiently long (e.g.,≫day) timescales, a secure detection of the line-to-continuum time delay may be more challenging (for an effectively similar case see Section 3.2.3).

Sampling is rarely uniform. In fact, seasonal gaps could be of the order of the time lag in some objects (Kaspi et al. 2000). The full treatment of sampling related issues is beyond the scope of this paper, and we restrict our discussion to light curve gaps (real data are treated in Section 4). To check for potential systematics, we resampled our simulated light curves with equal "on" and "off" observing periods with durations t_gap. We identify a bias when t_gap ≃ t^l_lag. Nevertheless, both the spectroscopic and photometric methods are equally affected by it so that their results remain consistent. As our aim is to focus on cases in which the photometric approach leads to new or different biases, we do not further quantify this bias.¹¹

Thus far, t^c_lag ≪ t_lag^l has been assumed. While the spectroscopic method is unaffected by inter-band continuum delays (the continuum in the immediate vicinity of the emission line is considered), our definition of the continuum is that which corresponds to the filter having the least contribution of emission lines to its flux, i.e., the X band.

For a large set of simulation runs (obtained by varying both η and t^c_lag/t_lag^l; not shown), we find that biases become non-negligible for η⁻¹|t^c_lag|/t_lag^l ≳ 1 (t^c_lag> 0 implies that continuum emission in the Y band lags behind the X band). Here, we consider particular cases in which |t^c_lag|/t_lag^l ⩽ 0.3 (η ≃ 0.13 is assumed) and find that |t^c_lag|/t_lag^l ≳ 0.1 results in substantial,>20% bias (Figure 9) for the following reason: although modest values of |t^c_lag|/t_lag^l are considered, the contribution of continuum processes to both light curves by far exceeds that of the emission line. Specifically, when the continuum in the Y band lags behind the X band, the deduced time lag reflects more on t^c_lag than on t^l_lag. When the X band lags behind the Y band, and t^c_lag/t_lag^l is significant, the peak in both statistical estimators is strongly suppressed and time-lag measurements become less reliable. Realistically, inter-band continuum time delays seem to be of order 10⁻²t^l_lag (Bachev 2009), and it is therefore likely that biases of the type discussed here would be negligible for most strong broad emission lines.

Realistic quasar spectra could result in more than one emission line contributing to the Y band, or even different emission lines contributing to both bands. While such complications are irrelevant for the spectroscopic method (note the small statistical scatter around unity in Figure 9), the photometric method is only sensitive to the combined contribution of all emission processes in some waveband. Here we study the biases introduced by including the contribution of a weaker secondary emission line to either of the bands whose time lag, t^l_{lag, s} ≠ t_lag^l.

We find potentially considerable biases when the relative contribution of the weaker line to the flux, η_s, satisfies η_s/η ≲ 1. However, the magnitude of the bias is also affected by (t^l_{lag, s} − t_lag^l)/t^l_lag. We do not attempt to fully cover the parameter space in this work, and focus on particular cases. Figure 9 shows an example with t^l_{lag, s}/t_lag^l = 0.5 where we find ≳ 20% biases for η_s/η ≳ 0.5: when the weaker line contributes to the Y band, t^l_{lag, measured} <t_lag^l while the opposite occurs when the weaker line contributes to the X band. Comparing to a case in which t^l_{lag, s} = t_lag^l/10 and η_s/η = 0.8, we find similar trends but an overall smaller(<5%) bias. Assuming instead t^l_{lag, s} = 2t_lag^l and η_s/η = 0.8, we find a bias of order 20% (10%) when the weaker line contributes to the Y (X) band. While the magnitude of the bias is clearly not solely determined by η/η_s, we find that the bias is, generally, small for η_s/η ≪ 1. Clearly, if t^l_{lag, s} ≪ t_sam or t^l_{lag, s}≫ t_tot then the time series is insensitive to the presence of the secondary line, and the bias would be negligible regardless of η_s/η.

Further quantifying the bias as a function of the line properties (e.g., transfer functions and time delays), filter response functions, quasar redshifts, and observing patterns is beyond the scope of this paper and should be carried out on a case by case basis.

It is well known that the ratio of the total duration of the light curve to the time lag, t_tot/t^l_lag, needs to be greater than some value to be able to adequately sample the transfer function and reliably measure the time lag using spectroscopic means. Our calculations demonstrate that the same holds also for the photometric method. In particular, for the Gaussian (thick shell) transfer function, reliable measurements require that t_tot/t^l_lag ≳ 3(9); see Figure 10.

Zoom In Zoom Out Reset image size

Differences between the photometric and the spectroscopic approach arise for t_tot/t^l_lag≫ 1: while the spectroscopic method gradually converges to t^l_lag (assuming sufficient data and adequate sampling; see Figure 10), the photometrically deduced time lag, t^l_{lag, measured}, develops a bias with respect to t^l_lag with increasing t_tot/t^l_lag. Generally, the bias is comparable (but not identical) for <ξ_CA> and <ξ_AA>, and overestimates the true lag. We quantify the bias for the two line transfer functions, ψ^l, defined above (Section 2.1.2), and find it to be different in each case: a broader ψ^l that extends over larger δt intervals leads to more biased results. For example, while the thick shell BLR model shows a factor ∼2 bias for t_tot/t^l_{lag, measured} = 50, the corresponding bias for the Gaussian BLR is only ∼40% (see Figure 10). The bias is seen to grow logarithmically with t_tot/t^l_{lag, measured}, and so would be most prominent when analyzing long time series of low-luminosity quasars. Consider, for example, ∼10 years of LSST data for a low-z quasar with a thick shell BLR whose size is ∼100 light days. In this case, t_tot/t^l_lag ∼ 30, which would result in t^l_{lag, measured}/t_lag^l ∼ 2. Not correcting for this bias would lead to an overestimation of the BH mass, by a similar factor, with considerable implications for BH formation and evolution.

Consider the bias in <ξ_CA>: as its value is sensitive to the line transfer function, its origin may be traced to the f_Y*f_X term in Equation (5), which is akin to the bias discussed in Welsh (1999). In a nutshell, with the quasar's power spectrum being red, most of the power is associated with the lowest frequency modes. Nevertheless, those are the same modes that are not adequately sampled by any finite time series. As such, the data appear to originate in a non-stationary process where, for example, the mean flux changes substantially between the extreme ends of the light curve. This and the fact that emission lines' contribution to the flux is always additive result in a positive bias. More extended line transfer functions, resulting in longer-term correlations between the light curves, are therefore more prone to such biases. With <ξ_AA(δt)> ≃ <ξ_CA(δt)>, similar biases are also encountered there that trace back to the f_Y*f_Y term.

To verify the source of the bias, we repeated the calculations using a truncated form for the power spectrum where all frequencies below a minimum frequency, ω_min, have no power associated with them (we have experimented with 2π/t_lag≫ ω_min≫ 2π/t_tot). The resulting <ξ_CA(δt)> is indeed unbiased and peaks at t^l_lag (not shown). For intermediate cases where, for example, the power spectrum flattens at low frequencies (Czerny et al. 1999), the results will be less biased but may not be entirely bias free.

Having identified the source of the bias, it is now possible to correct for it using one of several methods, which we term de-trending, sub-sampling, and bootstrapping. De-trending has been shown to be useful to correct for biases in the spectroscopic method (Welsh 1999). Basically, one filters out long-term trends in the light curve by fitting (done here in the least-squares sense) and subtracting off a polynomial of some degree, d, from the data. As a polynomial of degree d has d zeros, subtracting a higher degree polynomial from the data results in a suppressed variance at frequencies ω ≲ d(2π/t_tot). Conversely, the highest degree polynomial that can be subtracted from the data, while not suppressing the reverberation signal, is d ≲ ⌊t_tot/t_lag⌋. The effect of de-trending is shown in Figure 10 for a range of d-values and assuming t_tot/t^l_lag = 50. While subtracting a d = 2 polynomial already decreases the bias substantially, a bias-free result is recovered by de-trending with a d = 20(= 0.4t_tot/t^l_lag) polynomial. Further increasing d suppresses the sought after signal to insignificant levels by removing modes with frequencies ∼2π/t^l_lag. The dependence of the remaining bias on the de-trending polynomial degree, d, is shown in the inset of Figure 10.

An alternative method to remove the bias is to analyze continuous sub-samples of the time series, thereby reducing the effective t_tot and decreasing the bias. This method can be implemented in cases where one is not data-limited.

The third method for removing the bias involves bootstrapping: given an observed continuum light curve (i.e., f_X), and an assumed ψ^l, the light curve of the line-rich band is simulated. Upon sampling it using the cadence of the real data, an artificial f_Y light curve is obtained. By calculating the statistical estimators for the semi-artificial data, it is possible to compare t^l_{lag, measured} to the input lag and correct for the bias. Conversely, should the bias be known by other means, inferences concerning the nature of the line transfer function may be obtained.

In analyzing real data for individual objects in the PG sample (Table 1), it is important to verify that the detected signal is real and is not some artifact of non-even sampling or the particular light curve observed. There are several methods in the literature for cross-correlating realistic data sets: the interpolated cross-correlation function (Gaskell & Sparke 1986, ICCF; used in Section 3) and the discrete correlation function (Edelson & Krolik 1988, DCF). While the ICCF may be somewhat more sensitive to low-level signals, it is also more prone to sampling artifacts (Peterson 1993). To verify the significance of the ICCF signal, two algorithms are often used to estimate the uncertainty in the deduced correlation function: the flux randomization (FR) and the random subset selection (RSS). In a nutshell, the first algorithm adds noise according to the measurement errors and repeats the analysis many times to estimate the uncertainty associated with the results, while the latter method resamples the light curves (while losing ∼37% of the data) to reach the same goal. Quite often both tests are combined to form the FR–RSS method, which is our method of choice in the present analysis. When considering the DCF approach, the z-transformed DCF (ZDCF) algorithm by Alexander (1997) often produces the best results and is the approach adopted here. While all the aforementioned algorithms have been thoroughly tested with respect to spectroscopic reverberation mapping, their adequacy for the photometric reverberation mapping method is yet to be examined.

Our implementation of the ICCF (FR and RSS algorithms alike) is standard (Peterson et al. 1998). Nevertheless, because photometric reverberation mapping requires the subtraction of two cross-correlation functions, and to avoid interpolations of any kind,¹⁴ our implementation of the ZDCF algorithm requires that both the f_Y*f_X and f_X*f_X terms in Equation (5) are evaluated over an identical time grid, and that the bin center is identified with the time tag for that bin, i.e., δt. Similar considerations apply also when evaluating ξ_AA. Furthermore, we require that the number of independent pairs contributing to each bin in the correlation functions is ⩾11 (Alexander 1997). For simplicity, we use equally spaced bins to evaluate ξ_CA and ξ_AA, and repeat the analysis many times by randomly selecting independent pairs (Alexander 1997). This provides us with a distribution of ξ_CA and ξ_AA values at each time bin, from which the mean and its uncertainty (associated with one standard deviation) are evaluated.

Spectroscopic reverberation mapping often employs two distinct approaches to estimate the time lag: by measuring the position of the peak of the correlation function or by measuring its centroid. While the results are comparable, they are not identical and many works have discussed the pros and cons associated with each approach (Kaspi et al. 2000, and references therein). Similar algorithms are implemented here.

Estimating the uncertainty on the measured time lag is done in the following way: we use the set of Monte Carlo simulations described above and calculate ξ_CA, ξ_AA for each realization. The time lag is identified either with the peak or with the centroid of the statistical estimators,¹⁵ and a distribution of time lags is obtained. The uncertainty on the time lag is associated with one standard deviation of that distribution.

In addition to the above algorithms, aimed at testing the robustness of the signal, we provide an additional test to identify systematic effects of the type discussed in Section 3.2.5 and to assist in screening against problematic data sets that may lead to spurious signals. To this end we use the bootstrapping method discussed in Section 3.2.5, verifying that the mean contribution of the simulated emission line to the flux of the line-rich band is that given in Figure 11, and that the variance of the line light curve is consistent with the Kaspi et al. (2000, see their Table 5) results.

The B and R light curves for PG 0804+761 (henceforth termed PG 0804), at z = 0.1, are shown in Figure 12 where large non-monotonic flux variations are evident. This object has the best-sampled light curves in the Kaspi et al. (2000) sample, with a mean uncertainty on the flux measurement in both bands being of order 1% (with maximum uncertainties on individual measurements being ∼2%). The normalized variability measure is among the largest in the sample ∼17%/14% for the B/R bands (Table 1). As such, it is a promising candidate for photometric reverberation mapping.

Zoom In Zoom Out Reset image size

Figure 11 indicates that, at its redshift, the R band is relatively line-rich (the net contribution of the lines to the bands is ∼6%), with Hα being the dominant emission line. We calculated ξ_CA and ξ_AA using the ICCF and ZDCF algorithms. Unsurprisingly, the ZDCF method is somewhat more noisy than the ICCF methods (Kaspi et al. 2000, see their Figure 4). That said, we find the results of all methods of analyses to be consistent given the uncertainties. This means that interpolation is not a major source of uncertainty for this object, which is consistent with the relative regular sampling of its light curves.

Our analysis shows the expected features: a hump in both ξ_CA and ξ_AA, which peaks at similar (δt ∼ 200 days) timescales. For comparison, we used a published light curve for the broad Hα emission line (Kaspi et al. 2000),¹⁷ and calculated ξ_lc using the ZDCF algorithm. Clearly, both statistical estimators defined here trace the line-to-continuum cross-correlation function with a varying degree of significance: ξ_CA shows a ≳ 2σ peak (σ being the uncertainty around the maximum), while the peak in ξ_AA is only ≳ 1σ significant.

Using the Monte Carlo simulations described in Section 4.2.1, peak or centroid time distributions were obtained and are shown in Figure 13. Clearly, time lag estimates qualitatively agree between the different methods. Quantitatively, however, there are some differences: the centroid method results in somewhat larger values for the time lag. This results from the non-symmetric nature of ξ_CA, ξ_AA showing a more extended wing toward longer time delays (Figure 12). Results based on the ZDCF approach result in errors, which are typically smaller than those obtained using the ICCF method. This results from the fact that the ZDCF peak is largely determined by a few individual bins (note the two bins around 200 days that lie above most other ZDCF points, as well as above the ICCF curve; Figure 12). In fact, the resulting uncertainty on the ZDCF, being of order 20% in this case, is smaller than predicted based on our Monte Carlo simulations presented in Figure 7. We conclude that estimating the uncertainty on the time lag using the ZDCF results may be less reliable.

Zoom In Zoom Out Reset image size

Using the ICCF time-lag distributions, we determine the observed time lag for PG 0804 to be 250 ± 100 days, which is statistically consistent with the spectroscopically measured value of 193⁺²⁰_{− 17} days (Kaspi et al. 2000). As far as systematics is concerned, we do not expect a bias in the measured lag to exceed ∼40% given that t_tot/t^l_{lag, measured} ≳ 10 (Figure 10).

Using the numerical bootstrapping scheme defined in Section 4.2.2, we find that an input time lag of 100–200 days is qualitatively consistent with the data (Figure 14). In particular, much shorter or much longer time lags are inconsistent with our findings in terms of signal amplitude and shape. Interestingly, a weak yet marginally significant hump at around 200–300 days is implied by the data even when much shorter artificial time lags are used as input. Therefore, one should be careful when interpreting marginally significant signals. For the specific case considered here we can conclude that the data cannot be used to detect short time lags, which is of no surprise given that the cadence of the light curves is at 20–30 day intervals.

Zoom In Zoom Out Reset image size

To conclude, the time lag estimated here by the photometric reverberation mapping approach is ≃ 200 ± 100 days and is therefore consistent with the value deduced by Kaspi et al. (2000), albeit with an expectedly larger uncertainty.

This low-z, low-luminosity object has a potentially prominent (∼8%) line contribution to the R band from the Balmer emission lines. With ≲ 70 points in its light curve, and what appears to be significant variations over 200 days timescales, our analysis reveals a marginal (≳ 1σ using the FR–RSS method and ≳ 2 using ZDCF statistics) peak at ∼500 days (Figure 15). The spectroscopically measured time lag for this object is of order 50 days and agrees with the value expected from the BLR size versus luminosity relation (Table 1). Nevertheless, the data do not show a corresponding peak at those timescales. In fact, our simulations indicate that only a time lag ≳ 100 days could have been detected with some certainty. Interestingly, simulations with an input lag of 500 days (Section 4.2.2) are able to qualitatively account for the observed signal (dashed magenta curve in Figure 15).

Zoom In Zoom Out Reset image size

It is not clear whether a 500 day lag, if real, is associated with a long-term reverberation signal of the BLR, since there is no corresponding signal in the Kaspi et al. (2000) analysis of the spectroscopic data. We note, however, that it is also possible that the measurement uncertainties are somewhat underestimated for this object, hence the significance of our deduced signal overestimated. That this might be the case is hinted at by the normalized variability measure in both the B and R bands being similar (contrary to naive theoretical expectations) and the fact that contamination by the host galaxy in this low-z quasar may be important (see below).

This quasar, being at low redshift (z = 0.064), has one of the most prominent (≲ 8%) Balmer line contributions to the R band among the PG objects in our sample. The analysis detects a significant trough at ≳ 200 days, in both statistical estimators, using both the ICCF and ZDCF methods. This is unexpected since a peak rather than a trough is predicted by the simulations for an expected lag of ∼34 days (Figure 15). Simulating cases with time lags of order 200 days also produces a peak, in contrast to the observed signal. The interpretation that the B rather than the R band is the line-rich band is not supported by the optical spectrum showing rather typical quasar emission lines (Kaspi et al. 2000).

We cannot positively identify the origin of the unexpected signal but note that, at its low redshift, and being a rather faint low-luminosity quasar, the host galaxy is likely to substantially contribute to the flux in the bands. This has a dual result: (1) the actual contribution of the emission lines to the bands is smaller than that predicted here and (2) the data may be substantially affected by varying seeing conditions between the nights, thereby contributing to the variance in the light curves and masking the true signal. Interestingly, the reduced variability measure in both bands is similar for this object, which may indicate the presence of an additional source of noise not accounted for by the reported measurement uncertainties.

Analyzing the data for this object, a marginally significant signal (at the 1σ–2σ level) is detected in both statistical estimators, at ≲ 300 days (Figure 15). Formally, ξ_CA statistics implies that the time lag is 300 ± 130 days. At its redshift, the net contribution of emission lines is of order 5% and is therefore larger than the typical flux measurement uncertainty in this object (≲ 2%). Simulations with an artificial lag of ∼300 days are able to qualitatively account for the observed signal. Much shorter time lags (e.g., ∼40 days, as spectroscopically deduced by Kaspi et al. 2000) do not seem to be supported by our findings. Interestingly, for this object, spectroscopic time-lag measurements using the peak distribution function give very large errors for the Hα line, and given its luminosity, PG 1613+658 lies considerably below the BLR size versus luminosity relation (Kaspi et al. 2000, see also our Table 1). Using our deduced lag, a better agreement is obtained with the Kaspi et al. (2000) size–luminosity relation.

This object shows a statistically significant (2σ–3σ), wide peak extending over the range 150–450 days. Assuming that Hβ is the main line contributor to the R band's flux, our simulations indicate that the signal is consistent with a time lag of ∼430 days (formally, ξ_CA statistics yields a 300 ± 160 day' lag). This value is statistically consistent with the time-delay measurement of Kaspi et al. (2000), but our deduced lag is in somewhat better agreement with the expected BLR size given the source luminosity (Table 1). Our simulations show that time lags considerably shorter (≪200 days) or longer (≫800 days) are less favored by the data and cannot fully account for the observed signal (Figure 15).

This object, having ∼40 visits in each band, shows a marginal (≳ 1σ) hump at 200–300 days. Formally, the deduced time lag is 400 ± 200 days and is consistent with the values of Kaspi et al. (2000). At z = 0.371, the net contribution of emission lines to the line-rich R band is only about twice as large as the mean measurement uncertainty. This and the small number of data points account for the low significance of the observed signal. Further, simulations show that a signal should indeed be marginally detectable at times corresponding roughly to the spectroscopically determined lag.

There is a hint for a localized (∼1σ) peak around 200–300 days in both statistical estimators (ξ_CA statistics yields a lag of 240 ± 100 days) using both the ICCF and ZDCF methods. At face value, this timescale is consistent with the results of Kaspi et al. (2000) who find a time lag of ∼200 days. Nevertheless, as discussed by Grier et al. (2008), the cadence of the Kaspi et al. (2000) observations may not be adequate to uncover the true time lag in this object. Specifically, Grier et al. (2008) find a time lag of order 30 days using better sampled light curves on short timescales. The fact that we do not find a significant peak on short timescales is not surprising since we are essentially using the Kaspi et al. (2000) data set.

Like PG 0804, PG 0026+129 (z = 0.142) is among the best-sampled quasars in our sample with ∼70 points per light curve (see Figure 16). While the normalized variability measure is comparable to that of PG 0804, most of the power is associated with long-term variability. In particular, when subtracting a linear trend from the light curve, the normalized variability measure drops by a factor of ∼2 in both bands. At its redshift, the R band is the line-rich (Figure 11) with a net relative line contribution of ∼3% due to Hα (or ∼7% if iron is important).

Zoom In Zoom Out Reset image size

The statistical estimators calculated here do not yield a significant time lag. Specifically, there is no discernible peak detected in either estimator using either cross-correlation method (i.e., ICCF or ZDCF). While de-trending seems to somewhat improve the signal (there is a 1σ hump at around 100 days), it falls short of revealing a significant result. We attribute our failure to securely identify a signal to a combination of factors related to the smaller variability measure on the relevant BLR-size timescales, as well as to the lower net contribution of emission lines to the filters, compared to the case of PG 0804.

PG 0052+251 has 72 points in each band, and a normalized variability measure of 26%/20% in the B/R bands, partly due to a sharp feature extending over ∼200 day scales. The typical measurement uncertainty is ≲ 3% in both bands. At z = 0.155, the net line contribution is of order 4%, i.e., comparable to the noise level. Therefore, despite the relatively large number of visits, no clear signal is detected, in either statistical estimator, using either numerical scheme.

The remaining quasars in the Kaspi et al. (2000) PG sample also do not yield a significant signal. This is however expected given the small number of points in the photometric light curve of these quasars being ∼40 (see Figure 7). The fact that such objects do have a spectroscopically determined time lag is also not surprising: the photometric approach to reverberation mapping is less sensitive than the spectroscopic one.

To overcome the limitations plaguing the analysis of individual objects, ensemble averages may be considered. As demonstrated in Section 3, combining the results for many objects allows the time lag to be determined with good accuracy and was a main motivation on our part to develop this approach to reverberation mapping. Clearly, large ensembles are required to reach a robust result when less than favorable light curves are available on an individual case basis. Nevertheless, with upcoming photometric surveys, statistics is unlikely to be a limiting factor, and defining ensemble averages over narrow luminosity and redshift bins will be possible.

As our small sample of 17 PG quasars consists of a range of redshifts and quasar luminosities, we use the following transformation to put the results for different objects on an equal footing: δt → δt/((1 + z)L^γ₄₅), where L₄₅ ≡ λL_λ(5100 Å)/10⁴⁵ erg s⁻¹. This transformation corrects for both cosmological time dilation, and the fact that the BLR size scales roughly as L^γ. Following Kaspi et al. (2000), we take γ = 0.7 but note that similar results are obtained for 0.5 ≲ γ ≲ 0.7 (Bentz et al. 2009, and references therein). Using the above transformation, all PG quasars are effectively transformed to a z = 0 quasar with L₄₅ = 1.

Figure 17 shows the ensemble-averaged results for ξ_CA, ξ_AA for the entire PG sample (Table 1). We deliberately include all objects regardless of data quality issues (as in the case of, e.g., PG 1229). A peak in both (mean) statistical estimators is seen around 200 days. This value is similar to a time lag of 170 days, which is the value predicted by the Kaspi et al. (2000) BLR size–luminosity relation for an L₄₅ = 1 quasar. Also plotted are the median results for both statistical estimators. Both the mean and the median provide qualitatively consistent results, although quantitative differences do exist. Given the small nature of the sample, we do not present estimates for the uncertainty on the mean/median values.

Zoom In Zoom Out Reset image size

To conclude, averaging the results for the PG quasar sample, we obtain a similar peak in both statistical estimators, occurring at times which are qualitatively consistent with the BLR sizes implied by the Kaspi et al. (2000) results. A larger sample of objects will allow to better quantify the mean time lag for sub-samples of quasars over narrow luminosity and redshift intervals. Alternatively, better quality data (higher S/N and better-sampled light curves) could be used to better constrain the time lag for individual objects in the PG sample of quasars.