An Alternative Approach To Measuring Reverberation Lags in Active Galactic Nuclei
aa r X i v : . [ a s t r o - ph . C O ] M a y Draft version May 29, 2018
Preprint typeset using L A TEX style emulateapj v. 11/10/09
AN ALTERNATIVE APPROACH TO MEASURING REVERBERATION LAGS IN ACTIVE GALACTICNUCLEI
Ying Zu, C.S. Kochanek, and Bradley M. Peterson
Draft version May 29, 2018
ABSTRACTMotivated by recent progress in the statistical modeling of quasar variability, we develop a new ap-proach to measuring emission-line reverberation lags to estimate the size of broad-line regions (BLRs)in active galactic nuclei. Assuming that all emission-line light curves are scaled, smoothed, anddisplaced versions of the continuum, this alternative approach fits the light curves directly using adamped random walk model and aligns them to recover the time lag and its statistical confidencelimits. We introduce the mathematical formalism of this approach and demonstrate its ability to copewith some of the problems for traditional methods, such as irregular sampling, correlated errors, andseasonal gaps. We redetermine the lags for 87 emission lines in 31 quasars and reassess the BLRsize–luminosity relationship using 60 H β lags. We confirm the general results from the traditionalcross-correlation methods, with a few exceptions. Our method, however, also supports a broad rangeof extensions. In particular, it can simultaneously fit multiple lines and continuum light curves whichimproves the lag estimate for the lines and provides estimates of the error correlations between them.Determining these correlations is of particular importance for interpreting emission-line velocity–delaymaps. We can also include parameters for luminosity-dependent lags or line responses. We use thisto detect the scaling of the BLR size with continuum luminosity in NGC 5548. Subject headings: galaxies: active — galaxies: nuclei — galaxies: Seyfert — quasars: general INTRODUCTION
While it is widely accepted that the enormous lu-minosities of active galactic nuclei (AGNs) are at-tributable to accretion of matter onto supermassive blackholes (BH), detailed studies are extremely challengingon account of the small angular scales of the regions in-volved in the accretion process. Direct probes of thesub-microarcsecond structure of AGNs has been there-fore limited to VLBI studies of the radio-emitting re-gions, gravitational microlensing studies of the accretiondisk (see review by Wambsganss 2006) and reverberationmapping of the broad-line regions (Blandford & McKee1982; Peterson 1993). The technique of reverberationmapping (a.k.a. echo mapping) exploits the light traveltime between the central engine and the broad-line re-gion (BLR) to deduce the structure of the BLR (seePeterson 2001 for a tutorial). The continuum radiationfrom the accretion disk photoionizes gas clouds near theAGN to produce broad emission lines, thus encoding thegeometry and kinematics of the clouds (Osterbrock 1989;Peterson 1997; Krolik 1999). The physical ansatz for re-verberation mapping is straightforward:1. The continuum emission of the quasar shows(stochastic) variability that drives emission-linevariations after a light travel-time delay.2. The unobservable ionizing UV continuum thatdrives the emission lines is simply related to theobservable satellite UV or optical continuum (i.e., Department of Astronomy, The Ohio State Uni-versity, 140 West 18th Avenue, Columbus, OH 43210;[email protected]. Center for Cosmology and AstroParticle Physics, The OhioState University, 191 West Woodruff Avenue, Columbus, OH43210 the pattern and phase of variations are closely cor-related).3. The light-travel time is the most important timescale; specifically, the local emission-line responsetime to continuum changes is assumed to be instan-taneous and the dynamical time scale of the BLRis much larger than the light-travel time across it.The relationship between the observables, continuumlight curve s c ( t ) and the emission-line light curve s l ( t, V )where V is the line-of-sight velocity, is taken to be s l ( t, V ) = Z dτ Ψ( τ, V ) s c ( t − τ ) , (1)where Ψ( τ, V ) is known as the “transfer function” or“velocity–delay map.” In reality, the relationship be-tween the continuum and emission-line variations can benon-linear, but the amplitude of variation on reverbera-tion time scales is sufficiently small that the linear ap-proximation seems to be justified. Inspection of Equation(1) shows that Ψ( τ, V ) is the observed response of thebroad emission-line region to a delta-function continuumoutburst, mapped into the observable quantities time de-lay τ and line-of-sight velocity V . The data requirementsfor successful recovery of the transfer function are quitedemanding (Horne et al. 2004) and consequently mostefforts to date have concentrated on measuring only the total emission-line response to continuum variations. Thetransfer equation (eq. 1) then becomes s l ( t ) = Z dτ Ψ( τ ) s c ( t − τ ) , (2)where Ψ( t ) = R Ψ( τ, V ) dV is variously known as the“one-dimensional transfer function” (so that Ψ( t, V ) is Zu et al.the “two-dimensional transfer function” or the “delaymap”). For the remainder of this paper, we will refer toΨ( t ) simply as the “transfer function.” In most investi-gations to date, it is the mean response time or “lag” h τ i that one tries to measure, generally by cross-correlationof the continuum and emission-line light curves, as wediscuss further below. The importance of measuring theemission-line lag is two-fold: first, h τ i yields a charac-teristic physical scale for emission of a particular line, R = c h τ i , and this can be combined with some measureof the emission-line Doppler width ∆ V to obtain an es-timate of the central black hole mass. Assuming thatgravity is the dominant force on the line-emitting gas,the virial equation for the central black hole mass is M BH = f ∆ V RG , (3)where G is the gravitational constant and f is a dimen-sionless factor of order unity that depends on the geome-try, velocity field, and inclination of the BLR. We note inpassing that there is currently an active debate about therelative importance of radiation pressure on the BLR gasand how this affects reverberation-based mass measure-ments (Marconi et al. 2008; Netzer 2009; Marconi et al.2009; Netzer & Marziani 2010). While the possible roleof radiation pressure in measurement of black hole massesis an important issue, it has no direct bearing on thepresent discussion, which is about measuring time delays.Similarly, there is still active discussion about the meanvalue of the scaling factor h f i (e.g., Onken et al. 2004;Labita et al. 2006; Woo et al. 2010; Graham et al. 2011)that is beyond the scope of this contribution. Second,reverberation studies have established a tight empiricalrelationship between the BLR size and the AGN con-tinuum luminosity (Kaspi et al. 1996; Kaspi et al. 2000,2005; Bentz et al. 2006, 2009) that allows us to use theluminosity as a surrogate for the BLR size in eq. (3) andthus estimate the masses of black holes in AGNs from in-dividual spectra (Wandel et al. 1999; Vestergaard 2002;McLure & Jarvis 2002, 2004; Vestergaard & Peterson2006; Kollmeier et al. 2006; Shen et al. 2008). Thisallows us to explore BH properties and evolutionwith redshift (e.g., Kollmeier et al. 2006; Peng et al.2006; Hopkins & Hernquist 2006; Shankar et al. 2009;Steinhardt & Elvis 2010; Kelly et al. 2010), thus pro-viding valuable insights into the mystery of black holegrowth and its connection to galaxy evolution at highredshift, where the quasar population is evolving dra-matically. The potential for obtaining simple esti-mates for the masses of black holes in quasars pro-vide a means of exploring the correlations between theBH mass and global properties of their host galax-ies such as the bulge luminosity ( M BH – L bulge rela-tionship; Kormendy & Richstone 1995; Magorrian et al.1998; Bentz et al. 2009) and bulge stellar velocitydispersion ( M BH – σ ⋆ relationship; Ferrarese & Merritt2000 Gebhardt et al. 2000a,b; Tremaine et al. 2002;Ferrarese et al. 2001; Onken et al. 2004; Nelson et al.2004; G¨ultekin et al. 2009) both locally and potentiallyover cosmic time.As a practical problem in aperiodic time-series dataanalysis, reverberation mapping requires high-fidelityspectroscopic monitoring of the continuum and emission- line variations for a duration long compared to theemission-line lag (Horne et al. 2004), which is observedto range from hours to a year or more, depending onthe luminosity of the AGN and the cosmic time dila-tion at its redshift. Emission-line lags have been mea-sured for more than three dozen AGNs by cross corre-lation of the continuum and emission-line light curves.The particular challenge of dealing with reverberationtime series is that they are generally irregularly sam-pled for various reasons, including unfavorable weatherand, for higher-luminosity objects with larger lags, an-nual conjunctions with Sun that cause seasonal gaps inthe observed time series. In practice, two methodolo-gies have been widely employed to deal with unevenlysampled data. The first method is to interpolate be-tween real data points to obtain a regular samplinggrid for computation of the cross-correlation function(CCF) as a function of time delay τ (Gaskell & Sparke1986; Gaskell & Peterson 1987; White & Peterson 1994;Peterson et al. 1998; Welsh 1999; Peterson et al. 2004).The second method, the discrete correlation function(DCF) method (Edelson & Krolik 1988), bins the dataover discrete time intervals on which the data are rea-sonably well-sampled and a correlation coefficient iscomputed for the time delays between each pair ofcontinuum/emission-line time bins. A variant on this isthe Z -transformed DCF (Alexander 1997) which variesthe width of the time bins to better distribute the datapoints among the time bins. White & Peterson (1994)show that when common assumptions and normaliza-tions are used, the interpolation CCF method and theDCF method give similar results. However, as the time-sampling becomes sparser, the interpolation method sig-nificantly outperforms the DCF method as long as inter-polation of the light curves (usually linear in practice)remains a reasonable assumption .Presumably even more accurate lags could be mea-sured given more realistic modeling of the continuumbehavior between real measurements of the continuumand emission-line fluxes. This now seems to be a realpossibility given the recent work of Kelly et al. (2009),who find that quasar variability can be well describedby a damped random walk. By applying the variabilitymodel to the light curves of known quasars and com-paring them to other variable sources, Koz lowski et al.(2010) show that quasars occupy a very distinctive re-gion in the model parameter space of time scale and vari-ability amplitude. MacLeod et al. (2010) then apply themodel to ∼ Consider, as an example, the UV and optical monitioringcampaign on NGC 5548 undertaken with the
International Ul-traviolet Explorer (Clavel et al. 1991) and ground-based telescopes(Peterson et al. 1991) in 1989. The UV data were sampled at ap-proximately regular 4-day intervals. Analysis of these data usingthe interpolation CCF (Peterson & Wandel 1999) revealed for thefirst time a “virial relationship” between the emission-line lags andline widths (i.e., h τ i ∝ ∆ V − ), thus providing an empirical jus-tification for using Equation (3) to estimate the black hole mass.Analysis of these same data with the DCF method (Krolik et al.1991) obscured this result at least in part because of “discretizationnoise” introduced by the DCF: because of the regular 4-day sam-pling, the smallest usable time bin for the DCF was also four days,resulting in lag measurements that were integer multiples of 4 days,significantly reducing the time resolution of the lag measurementsand smearing out the h τ i –∆ V anticorrelation. easuring Reverberation Lags 3scribe quasar variability well and they explore the corre-lation of the variability parameters with other propertiesof quasars such as wavelength, luminosity, BH mass, andEddington ratios in detail. Importantly for reverbera-tion studies, the formalism is able to statistically predictthe value of light curve at an unmeasured time basedon the overall statistical properties of the light curve. Itprovides a well-defined statistical model for interpolatinglight curves and can do appropriate statistical averagesover the uncertainties in the model predictions.Given a complete statistical framework for describingthe continuum variability, and the overall ansatz thatemission-line variability is a scaled and smoothed versionof the continuum, we can build an alternative approachto measuring reverberation lags, aspects of which werepreviously noted by Rybicki & Kleyna (1994). Amongthe advantages of this approach are:1. It not only interpolates between data points, butalso self-consistently estimates and includes the un-certainties in the interpolation.2. It can separate light curve means, trends, and sys-tematic errors in flux calibration from variabilitysignals and meansurement noise in a self-consistentway.3. Correlated errors can be treated naturally.4. Lags of multiple emission lines and their covari-ances can be derived simultaneously.5. It provides statistical confidence limits on the lagestimates as well as other parameters.We describe the methodology of our approach in detailin §
2. In §
3, we present the statistical process model forthe continuum light curves. We briefly describe our dataset and apply this method to the estimate of H β lagsin §
4. We further show how the method can address theproblem of correlated errors in § §
6. We also fit the R BLR – L relationship using H β lags determined by our method in §
7. In §
8, we add aluminosity dependence to the lag and solve for the lag–luminosity relationship of NGC 5548. We summarizeour main findings and discuss future applications andexpansions of our approach in § METHODOLOGY
Press et al. (1992) and Rybicki & Press (1992) devel-oped a method to statistically analyze irregularly sam-pled light curves, and Rybicki & Kleyna (1994) appliedthe variant we now consider to four seasons of optical re-verberation data on NGC 5548. Here we reintroduce thisapproach, which we have named “Stochastic Process Es-timation for AGN Reverberation (SPEAR ),” with sev-eral modest changes in algorithm and a broad range ofnew applications.Except for the transfer function Ψ( τ ), our notationis chosen for comparison with Rybicki & Kleyna (1994). We start with a model process driving the continuum s c ( t ) that has a covariance between times t i and t j of h s c ( t i ) s c ( t j ) i = σ exp( −| t i − t j | /τ d ) . (4)We adopt here an exponential covariance matrix forconcreteness, since we know from Kelly et al. (2009),Koz lowski et al. (2010) and MacLeod et al. (2010) thatquasar light curves are well modeled by this process.Physically, the model corresponds to a random walk de-scribed by an amplitude σ = ˆ σ τ d / τ d , where ˆ σ and τ d are used as our model parameters. Rybicki & Press(1992) estimated the covariance matrix based on thestructure function of the continuum light curve, whilehere we adopt a specific parametrized model that will beoptimized as part of the analysis.Slightly rewriting Equation (2) for convenience and tofacilitate comparison with Rybicki & Kleyna (1994), thelight curve of a line is s l ( t ) ≡ Z dt ′ Ψ( t − t ′ ) s c ( t ) . (5)Since the lines and continuum are related by the transferfunction, we can also determine the covariance betweenthe line and continuum h s l ( t i ) s c ( t j ) i = Z dt ′ Ψ( t i − t ′ ) h s c ( t ′ ) s c ( t j ) i , (6)between the line and itself h s l ( t i ) s l ( t j ) i = Z dt ′ dt ′′ Ψ( t i − t ′ )Ψ( t j − t ′′ ) h s c ( t ′ ) s c ( t ′′ ) i , (7)and between two different lines h s l ( t i ) s ′ l ( t j ) i = Z dt ′ dt ′′ Ψ( t i − t ′ )Ψ ′ ( t j − t ′′ ) h s c ( t ′ ) s c ( t ′′ ) i . (8)If the light curve of the line is divided into velocitybins δV , then there is a transfer function for each binΨ( t − t ′ , V ) and we can compute all the expected covari-ances between the light curves. For convenience, let s be a vector comprised of all the light curves, both lineand continuum, and S = h ss i be the covariance matrixbetween all the elements of s . By definition, in Gaussianstatistics the probability of the light curve is simply P ( s ) ∝ | S | − / exp (cid:18) − s T S − s (cid:19) . (9)We do not measure the actual light curve, but some re-alization of it, y = s + n + L q , in which there is mea-surement error n , whose probability distribution is P ( n ) ∝ | N | − / exp (cid:18) − n T N − n (cid:19) . (10)where N = h nn i is the covariance matrix of the noise.Note that nothing requires N to be diagonal, so there isno formal difficulty to including covariances in the noisebetween the line and continuum.In defining y , we have also allowed for the simultaneousfitting of a general trend defined by a response matrix L and a set of linear coefficients q . In particular, we use thisto fit and remove separate means from the light curves. Zu et al. Fig. 1.—
Distribution of χ per degree of freedom for the con-tinuum fits. The solid histogram is the χ /dof distribution of ourstochastic model, while the dashed one shows the distribution ex-pected for models with correctly estimated Gaussian uncertainties.The dotted histogram is the χ /dof distribution of the joint modelof the continuum and H β light curves from § In this application to a model with two light curves, L isa 2 × K matrix with entries of (1 ,
0) for the continuumdata points and (0 ,
1) for the line data points, where K is the total number of data points. The linear parame-ters are a very general tool. For example, separate lineartrends would be removed with a 4 × K matrix with en-tries of (1 , t i , ,
0) for continuum epoch t i and (0 , , , t j )for line epoch t j . Two sources of data with potentiallydifferent but constant levels of contamination from thehost galaxy can be reconciled by using different meansfor each line and continuum data source, correspondingto a 4 × K matrix with entries of (1 , , ,
0) for the firstcontinuum source, (0 , , ,
0) for the second continuumsource, (0 , , ,
0) for the first line source and (0 , , , y given the linear coefficients q , the intrinsic light curves s , and any other parameters of the model p is P (cid:0) y (cid:12)(cid:12) q , s , p (cid:1) ∝ | SN | − / Z d n n d n s δ ( y − ( s + n + L q )) exp (cid:18) − s T S − s + n T N − n (cid:19) . (11)After evaluating the Dirac delta function, we “completethe squares” in the exponential with respect to both theunknown intrinsic source variability s and the linear co-efficients q . This exercise determines our best estimatefor the intrinsic variability ˆs = SC − ( y − L ˆq ) (12)and the linear coefficients ˆq = ( L T C − L ) − L T C − y ≡ C q L T C − y (13)where C = S + N is the overall covariance matrix of the data and C q = ( L T C − L ) − . With these definitions wecan factor the argument of the exponential into P (cid:0) y (cid:12)(cid:12) q , s , p (cid:1) ∝ | SN | − / exp − ∆s T ( S − + N − ) ∆s − ∆q T C − q ∆q − y T C − ⊥ y ! (14)where C − ⊥ = C − − C − LC q L T C − (15)is the component of C that is orthogonal to the fittedlinear functions, the variances in the linear parametersare h ∆q i = ( L T C − L ) − ≡ C q , (16) ∆s = s − ˆs and ∆q = q − ˆq . We can marginalize theprobability over the light curve s and the linear param-eters q under the assumption of uniform priors for thesevariables to find that P (cid:0) y (cid:12)(cid:12) p (cid:1) ∝ L≡ | S + N | − / (cid:12)(cid:12) L T C − L (cid:12)(cid:12) − / exp (cid:18) − y T C − ⊥ y (cid:19) (17)where L represents the likelihood function we are to max-imize, and the remaining parameters p are those de-scribing the process (Equation 4) and the transfer func-tions. The term in the exponent, y T C − ⊥ y , is the gener-alized χ that we present throughout the paper. Whilethis treatment of linear parameters was included byRybicki & Press (1992), Rybicki & Kleyna (1994) choseto subtract fixed means rather than marginalizing overthem as part of the analysis as we do here. The variancein the estimate for the mean light curve is h ∆s i = S − S T C ⊥ S. (18)We can estimate the light curve s ( t ) at any unmeasuredtime using the same formalism. The simplest means ofdoing so is simply to pad the data vector y d with ad-ditional fake points y f that have infinite measurementuncertainties in the sense that N − → ˆs f = S fd ( S dd + N dd ) − y d (19)with variance relative to the true light curve of h ∆s f i = S ff − S fd ( S dd + N dd ) − S df . (20)where S dd , S ff , S fd and S df are the data-data, fake-fake, fake-data and data-fake covariance matrices of theprocess and N dd is the noise matrix of the data. Theinclusion of the fake points has no effect on the expectedresults for the measured data points.Just to re-emphasize the point, this formalism wasfirst outlined by Rybicki & Kleyna (1994) based onPress et al. (1992) and Rybicki & Press (1992). We haverefined it slightly to use a specific process model, to op-timize the parameters of that model and to include themeans of the light curves as parameters that are auto-matically marginalized. Unfortunately, we will not beable to use the fast implementation of this method foreasuring Reverberation Lags 5 Fig. 2.—
Continuum models. The solid line shows the expected mean source light curve ˆs (Equation 12) and the dashed line shows theexpected spread (Equation 18) of light curves about the mean consistent with the data. An individual light curve realization consistentwith the data (see Equation 21) will show more structure than this mean light curve and have excursions outside the dashed line consistentwith the estimated variance. exponential covariance matrices from Rybicki & Press(1995), because the inclusion of the transfer functionsmeans that S is not a simple exponential covariance ma-trix and hence does not have a simple, tridiagonal inversefor the fast method.We can, however, use the fast methods for generatingsimulated light curves. In particular, we are interested inlight curves constrained to resemble the continuum lightcurve. As discussed by Rybicki & Press (1992), sucha light curve is simply the estimated mean light curvegiven by Equation (12) with an added random compo-nent that has the covariance matrix Q = ( S − + N − ) − .Rybicki & Press (1992) suggest determining the eigen-modes of Q which are then the independent “normal”modes that can be added to the mean light curve to pro-duce a random realization constrained by the continuumlight curve. This is computationally expensive. Instead,we note that if we Cholesky decompose Q = M T M , where M is an upper triangular matrix, and define therandom component of the light curve by u = M r where r is a vector of zero-mean, unit-dispersion Gaussian ran-dom deviates, that h uu T i = M h rr T i M T = M M T = Q T = Q (21)since the covariance matrix h rr T i of the Gaussian devi-ates is simply the identity matrix and Q is symmetric.Since Q − is a tridiagonal matrix given the exponentialcovariance matrix and a diagonal noise matrix, we cangenerate very high dimension u that can be convolvedwith the transfer function to produce a simulated linelight curve in O ( K ) operations rather than the O ( K )needed following the eigenmode approach.The original application of the method by Press et al.(1992) was to cross correlate the light curves of two im-ages of a lensed quasar in order to estimate the timedelay between them. While this was not discussed in Zu et al.terms of transfer functions, it does correspond to a trans-fer function of the form Ψ( t i − t j ) = δ ( t i − ( t j + ∆ t )),making the second light curve a lagged version of thefirst. Press et al. (1992) also treated the parameters cor-responding the process as fixed parameters, derived byfitting a power law to the structure functions of the lightcurve. It is likely that some combination of neglectinguncertainties in the process model or covariances in theerrors of the light curves led Press et al. (1992) to ob-tain an incorrect estimate of the time delay despite theelegance of the approach.In Rybicki & Kleyna’s (1994) expansion of the methodto reverberation mapping, they used rising and fallingsawtooth and isosceles triangle transfer functions, findinglittle difference between the results or ability to discrimi-nate between them. Thus, for this initial reconnaissance,we will simply use a top-hat (rectangular function) forthe transfer function,Ψ( t − t ′ ) = A ( t − t ) − for t ≤ t − t ′ ≤ t (22)which has a mean lag of h τ i = ( t + t ) / τ = t − t . The necessary integrals for Equa-tions (6), (7), and (8) are all analytic (see the Appendix)and the model includes the limits of a delta function as∆ τ → t →
0. The scal-ing coefficient A determines the line response for a givenchange in the continuum (i.e., the responsivity of BLRclouds), but for present purposes we will largely view itas a nuisance variable.We use the amoeba minimization method (Press et al.1992) to optimize the solution and then either a MonteCarlo Markov Chain (MCMC, Metropolis et al. 1953;Hastings 1970) or optimization over a grid to estimateparameter uncertainties. We carry out the analysis intwo phases. We first analyze the continuum light curveon its own, using logarithmic priors for τ d and ˆ σ to deter-mine the range of the variability process parameters con-sistent with the continuum light curve. The logarithmicprior on τ d essentially penalizes values that deviate fromthe median sampling intervals to avoid both unphysi-cally large τ d and a second class of solutions of τ d → σ to broaden the uncertainties until obtain-ing an acceptable fit. Then we do the joint analysis ofthe continuum and the lines using Gaussian priors for τ d and ˆ σ determined from the analysis of the continuum inisolation. In detail, we take the results of the MCMCanalysis of the continuum and used uncorrelated priorson ln τ d and ln ˆ σ (which is conservative), where the priorfor each variable was centered at the median value withthe Gaussian width chosen to match the upper and lower1 σ confidence regions. We then used uniform priors for A , t and t .The reason for using the continuum to define a strongerGaussian prior on the process variable before carryingout the joint analysis is to eliminate the aforementionedsecond class of solutions of τ d → Fig. 3.—
Rest-frame damping timescale τ d of the continuum lightcurves as a function of optical luminosity. The uncertainties in τ d are the ± σ range. multaneously at the wrong lag, the optimal solution willbe to let τ d → τ d andˆ σ that were statistically consistent with the continuumvariability. THE STATISTICAL PROCESS MODEL OF THECONTINUUM
This approach depends on using a statistical model forthe variability process of the continuum in order to opti-mally model the underlying light curve of the continuum.Here we use the exponential covariance matrix suggestedby Kelly et al. (2009), although it was also introducedby Rybicki & Press (1995) to enable a fast version of theSPEAR approach. Physically, the exponential covariancematrix in Equation 4 corresponds to a damped randomwalk with an amplitude scale ˆ σ and a damping time scale τ d . On long time scales the variance of the light curve isˆ σ ( τ d / / and on short time scales it is ˆ σ √ t .Kelly et al. (2009) use this to model the light curvesof 100 quasars, including some of the objects we willconsider here, using a light curve forecasting approachto estimate the process parameters. Koz lowski et al.(2010) show how the Kelly et al. (2009) approach canbe derived from the SPEAR approach and demonstratedthat forecasting is less statistically optimal for parame-ter estimation than using the complete light curve mod-eling method of SPEAR, and then applied the pro-cess model and the SPEAR method to the OGLE-III (Udalski et al. 2008) light curves of ∼ τ d –ˆ σ parameter space. This is further confirmedby MacLeod et al. (2010), who used this approach tomodel 9,000 SDSS quasars to examine the correlationsof σ and τ d with other quasar properties.Unlike the previous papers, we fit flux rather than mag-nitude light curves because the line flux is more closelyeasuring Reverberation Lags 7 Fig. 4.—
Comparison of rest-frame H β time lags from the CCF and the SPEAR methods. Green triangles, red squares and blue circleswere used for the PG, NGC and other objects, respectively, with ± σ error bars indicated on both estimates. The labeled points linked bydashed vertical lines are objects having multiple lag solutions and the filled symbol is the higher likelihood solution. The two intersectingstripes indicate the region where the solutions from both methods may be false due to the seasonal gap (140–200 days, with time dilation). related to the continuum flux than to the continuummagnitude. Thus, we start by examining how well thedamped random-walk process models the 60 continuumflux light curves for the 31 systems we consider in § χ per degreeof freedom for the best-fit models of all the continuumlight curves we consider. Since half of the continuumlight curves in our sample have less than 50 data points,the expected χ /dof distribution is broader than thatof the OGLE light curves ( ∼
500 points) considered byKoz lowski et al. (2010). Nevertheless, the χ /dof distri-bution indicates that the statistical process model pro-vides a reasonable fit to the light curves. The fact thatthe distribution is narrower than expected for correctlyestimated Gaussian uncertainties suggests that the re-ported photometric errors are somewhat larger than thetrue uncertainties, or that there has been some pruningof outliers from the light curves. Figure 2 shows three examples of modeled continuumlight curves interpolated and extrapolated from Equa-tion (12) and their uncertainties from Equation (18), aswell as the observed light curve. The estimated lightcurve at time t is in essence a weighted average over datapoints within the damping time | t − t ′ | < ∼ τ d that bal-ances the variance expected on those time scales due tothe process against the uncertainties in the data point todetermine how closely the model light curve approachesa particular data point. Far from any data points, themodel returns to the light curve mean on the time scale τ d . Remember, however, that Equation (12) is an esti-mate for the average of all possible light curves that couldbe drawn from the process that would be consistent withthe data — a particular realization of such a light curvewould show additional structure (see Rybicki & Press1992). The “error snake” surrounding the model lightcurve is the variance in these possible light curves. Near Zu et al. Fig. 5.—
Sensitivity of the lag estimates to the noise correlation coefficient r between the H β and the continuum light curves of PG 0844.The left top panel shows the dependence of the lag on the correlation coefficient r . The left bottom panel shows the corresponding changein the likelihood function with r at the best-fit lag. In these panels, the blue triangle, green circle, and red square mark the results for r = − , , +1, respectively, and the dotted line indicates the 3 σ limit of the likelihood function. The right panel compares the lag likelihooddistribution for these 3 cases: r = − data points, its width approaches that of the measure-ment errors and then grows as the distance ∆ t from anydata point increases. The variance from the process ini-tially increases as ˆ σ | ∆ t | / , but then saturates at theoverall process variance once | ∆ t | ≫ τ d . Thus, in theextrapolated regions we see the model light curve be-comes a constant and the error snake expands and thenbecomes constant.The three objects shown in Figure 2 represent threetypical levels of light-curve sampling quality for the ob-jects we consider. Generally, the light curves of thePalomar–Green (PG) quasars obtained by Kaspi et al.(2000) were sampled every 1–4 months over a baseline aslong as 7.5 yr, as opposed to most of the low-luminositySeyfert 1 AGNs that were more densely sampled overshorter baselines. The rest of the sample mainly con-sists of nearby bright Seyfert galaxies (Peterson et al.1998) whose light curves are sparsely sampled over ashort baseline. Peterson et al. (2004) discuss the datain detail. In addition, we also include new light curvesfrom a recent high sampling rate, multi-month rever-beration mapping campaign on six local Seyfert galax-ies (Denney et al. 2010).We expect the damping timescales τ d to show cor-relations with the physical characteristics of the accre-tion disk such as the mass of the central black hole,and the AGN luminosity (Peterson 2008). Kelly et al.(2009) demonstrated this scaling relationship between τ d and L AGN by performing a linear regression of τ d on L AGN , while Collier & Peterson (2001) also found apositive correlation between the characteristic timescaleand black hole masses. Their characteristic timescale,which is defined by the timescale where the structurefunctions flattened, is roughly equivalent to τ d . Fig-ure 3 shows that the more luminous central engines havelonger exponential damping timescales, as we would ex-pect from Kelly et al. (2009), up to any minor differ- ences from fitting fluxes rather than magnitudes. Notethat MacLeod et al. (2010) argue that the dependenceof τ d on black hole mass M BH is the real driver of thecorrelation between τ d and luminosity. We can use thesecorrelations to estimate τ d for sources lacking sufficientlygood light curves. ESTIMATING EMISSION-LINE LAGS
As our first application of the SPEAR method werecompute the lags of 101 emission-line light curvesfor 31 objects in the literature (the compilation ofPeterson et al. 2004 with the addition of data fromBentz et al. 2006, Grier et al. 2008, and Denney et al.2006, 2010). We carried out the analysis in three stages.First, as discussed in §
3, we modeled the continuum aloneto determine the range of process parameters ( τ d , ˆ σ ) con-sistent with the continuum light curve. We use this dis-tribution of estimated τ d and ˆ σ as a prior for the jointmodels of the continuum and line light curves in order toavoid the secondary solutions with τ d → §
2. Second, for each joint model, we find the best-fittop hat transfer function (Equation 22) which maximizesthe model likelihood calculated by Equation (17), alongwith an updated set of process parameters. Finally, weran an MCMC analysis on each joint model to calculatethe statistical confidence limits on each best-fit param-eter found by global optimization on a grid, especiallythe time lag. We then compare these estimates to thosederived from previous CCF analyses. We refer to thesemodels as the “single-line” fits since they are solving for asingle top-hat transfer function. The dotted histogram inFigure 1 shows the χ /dof distribution of the single-linemodel. It has a similar shape to the χ /dof distribu-tion of the stochastic model for only the continuum lightcurve, and confirms that the statistical model provides agood fit to the quasar variability, as well as the overall ansatz that the H β variability is a scaled and smoothedeasuring Reverberation Lags 9 -20-15-10-50 l n ( / m a x ) H Lag (Days) LL -20-15-10-500 10 20 30 40 50 60 H Lag (Days) H L a g ( D a y s ) H Lag (Days)
H Lag (Days) H L a g ( D a y s ) H Lag (Days) (a) single-linetwo-line (b) single-linetwo-line(c) single-line (d) two-line (e) single-linetwo-line
Fig. 6.—
Comparison between independent (single-line) and joint (two-line) fits to the H α and H β light curves of NGC 3516. Thered solid lines are the estimate from the single-line fits, while the blue dashed lines are those from the two-line fits. The top left (right)panel compares the likelihood distributions of the two fits for the H α (H β ) line. The interval between the two dotted lines correspondsto a 3 σ range in the likelihood, while the two blocks above indicate the ± σ range of the CCF peak analysis (upper) and CC centroiddistribution (lower), where the central lines mark the τ peak and τ cent values, respectively. The two bottom left panels shows the color-codedcovariance map between the two lags for the single-line and two-line fits, respectively. The contours in the bottom right panel comparethe confidence levels calculated from MCMC for the H α /H β lags near the peak (black boxes inside left two panels). Working outward, thethree contour curves are for 1 σ , 2 σ and 3 σ levels, respectively. Note that those are all observed-frame lags and the H β light curve here isthe older of the two we have for NGC 3516. version of the continuum. The χ /dof distribution ofthe single-line model is somewhat worse than for fittingthe continuum alone, but still reasonably consistent withstatistical expectations.For the sake of uniformity of emission-line species inthe comparison between the SPEAR and the CCF meth-ods, and to avoid confusion in the figures for sources withmultiple line observations, we will focus on the the 66H β light curves in our subsequent analyses, and tabulateall the other emission-line lags we successfully computedwith SPEAR in Table 1.Figure 4 shows the comparison between CCF centroidtime lags τ CCF and our lags τ SPEAR for all the H β lines.The range of uncertainties for τ CCF contains 68 .
3% ofMonte Carlo realizations in the cross-correlation centroiddistribution (CCCD), while our estimated error bound-aries are defined by the 68 .
3% (ln L / L max = 0 .
5) con-fidence levels that encloses the best-fit lags (i.e., ± σ errors if the probability distribution is Gaussian in bothcases). Based on the structure of the lag probability dis-tribution, we can divide the “single-line” fits into fivequality groups:(I) In most of the cases (43 of 66), the likelihood dis-tribution for the lags has a single peak and there isan unambiguous H β lag.(II) In 9 cases, the likelihood distribution has multi- ple peaks with significant ( > σ ) likelihood differ-ences. This occurred for one season of Akn 120(JD49980–50175) , Mrk 110 (JD48953–49149), andMrk 590 (JD49183–49338); two seasons of Mrk 79(JD48193–48393 and JD49996–50220); NGC 4051,PG 0844 , PG 1411, and PG 1617. Compared toour estimate, the CCF analysis picks a lower like-lihood peak or aliases several peaks into one broadpeak. Generally, the two peaks are so close that thedifferences between the results from the two meth-ods are insignificant compared to the uncertainties.(III) In 7 cases, the likelihood distribution has multiplepeaks of comparable significance ( ≤ σ : one seriesof NGC 3516 (JD47894–48047), Fairall 9, PG 0026,PG 0052, PG 1211, PG 1226, and PG 1307). Theyare shown in Figure 4 as the objects with a dashedline connecting the possible solutions. The tradi-tional CCF method seems to find one broad peakfor these sources, rather than multiple peaks, lead-ing to large reported uncertainties for the estimateof τ CCF . These degeneracies are largely caused For brevity, we retain only the five least significant digits ofthe Julian Date. For brevity, we truncate the PG coordinate identifiers to rightascension only since this introduces no ambiguity in the presentsmall sample. -20-15-10-50 l n ( / m a x )
50 100 150 200
H Lag (Days) LL -20-15-10-5050 100 150 200 250 H Lag (Days) H L a g ( D a y s )
50 100 150 200
H Lag (Days)
50 100 150 200 250
H Lag (Days) H L a g ( D a y s )
90 100 110 120
H Lag (Days) (a) single-linetwo-line (b) single-linetwo-line(c) single-line (d) two-line (e) single-linetwo-line
Fig. 7.—
Comparison between independent (single-line) and joint (two-line) fits to the H α and H β light curves of PG 0026. The formatis the same as in Figure 6. by poor light curve sampling that allows the lightcurve of the emission line to be mapped into thesampling gaps of the continuum. This problem isworst for the PG objects, which have many “sea-sonal gaps” over the long observing baselines ( ∼ ∼
40 days in the observed frame. Al-though we find a sub-peak at 40 days, the modelfavors another peak of much higher significance at259 days. For PG 1613 we obtain a lag of ∼ ∼
50 day CCF estimate. In both cases, the longerlag is favored because it minimizes the data over-lap — 259 and 575 days put most of the line datain the seasonal gaps, and many points also lie be-fore the start of the continuum light curve. This is essentially an aliasing problem in our method.We also note that the continuum flux varied by upto 50% over the 7 year span of the light curves.We know empirically that the scaling coefficient A in the transfer function is inversely correlated withionizing continuum flux (see the right panel of Fig-ure 10 and the discussion in § A asa constant parameter in each individual fit. Thismay create problems for light curves with the sig-nificant long term trends observed for these objects.Allowing A to vary and adopting a prior that pe-nalizes large lags that minimize light curve overlapwould likely solve these problems. MODEL TEST FOR CORRELATED ERRORS
Correlated errors have long been viewed as a prob-lem in traditional CCF analysis. Observations made ata common epoch are inevitably correlated by the pro-cesses required for calibration, light curve extraction,broad/narrow line modeling and removal of host or Fe ii contamination. Because no assumption about the prop-erties of the noise matrix N was made in §
2, it is easyto include the effects of correlated errors within our ap-proach. While we did not make an extensive survey ofour ability to model noise correlations between the con-tinuum and lines, we did carry out some experiments forobjects noted as potentially having strong covariances byPeterson et al. (2004).The simplest test is to introduce a covariance factor − ≤ r ≤ N for line and continuum points measured at thesame epoch of N cl ( t, t ) = N lc ( t, t ) = rσ l ( t ) σ c ( t ) in orderto examine the sensitivity of the lag estimates to corre-easuring Reverberation Lags 11lated noise between the line and the continuum measuredat each epoch. This should be present in the data at somelevel because of the challenge of consistently subtractingthe contribution of the host galaxy to the line and thecontinuum in the presence of variable observational con-ditions.Figure 5 illustrates the effects of adding off-diagonalcorrelated noise terms on the H β lag estimate of PG 0844.The shift in the estimated lag (left top panel) induced by r varying from − r > r = +1, suggesting that theerrors in the two light curves are positively correlated.The lag likelihood distribution (right panel) changes ifwe assume different levels of correlations r between thelight curves. While the overall lag likelihood is greatlydepressed in the r = − r = 0 and r = +1 cases. How-ever, the peaks near the best lag estimate ( ∼
12 days)are slightly more significant in the r = +1 case than inthe r = 0 case. We explored this issue for several othersystems, and generally the impact on the estimated lag isnegligible, although different levels of (anti-)correlationsare detected. JOINT ANALYSIS OF MULTIPLE LINES In §
4, we found that poor light curve sampling was asignificant problem in many systems, particularly in ob-jects with observed-frame lags on time scales similar tothe seasonal gap spacing. However, if multiple lines havebeen measured, then we have significant, additional datato better sample the light curves under our overall ansatz that all light curves are scaled, smoothed, and displacedversions of the continuum. Simultaneous fits also de-termine the covariance between the lags of the differentlines. In this section, we explore simultaneously fittingthe continuum and two emission-line light curves (here-after “two-line” fits, as opposed to the “single-line” fits in §
4, as we are now fitting two top-hat transfer functions).Figure 6 summarizes the significant improvement inestimating the H β time lag of NGC 3516 of theWanders et al. (1993) data (JD47894–48047) after in-cluding the H α light curve (two-line) compared to us-ing the H β line alone (single-line). NGC 3516 is a casewhere the single-line H β fits shows a secondary peak at ∼
42 days whose likelihood relative to the main peakat ∼ L nd / L max ) = − . α fit does not show such a secondarypeak (panel a). When we fit both simultaneously, the H α light curve together with its well-determined lag adds ex-tra information to the continuum light curve, and thusbetter constrains the H β lag. The second H β peak issuppressed and there is a single unambiguous H β lag forthe two-line fit (dotted curve in panel b). The improve-ment is most clearly seen in the structure of the H α /H β lag likelihood plane (panel c and d). If we zoom in on theremaining peak and run a MCMC chain using a flat prioron lags in the zoomed region, we can see that the two-line fits not only suppress the secondary peaks but alsoshrink the uncertainties in the primary peak to producebetter results for both lines (panel e).The joint analysis of multiple lines is especially useful for the PG objects, whose light curves show observationalgaps of period ∼
180 days in the observed-frame. In thesingle line fits, the model would always show (sub)peaksfor lags ∼
180 days because of the seasonal aliases (theseasonal stripes in Figure 4). It is not possible, however,to do this for 2 lines simultaneously, so the two-line fitslargely eliminate seasonal aliasing. Figure 7 illustratesthis for PG 0026. In particular, the broad H β likelihooddistribution shrinks significantly and the maximum like-lihood lag drops from ∼
160 days to ∼
106 days ( ∼
140 to ∼
93 in the rest frame) and is in better agreement withthe H α results. Although the traditional CCF methodmakes similar estimates (green and blue bands in twotop panels), it yields significantly larger uncertainties byaliasing several peaks into one broad CCF centroid dis-tribution.We performed similar joint analyses for the 21 sourcesfor which we have multiple emission line light curves andrecompile the results for the H β lags, as shown in Fig-ure 8. Fortunately, all the sources whose H β lags werefound to be ambiguous in the single-line fits (i.e., the7 H β lags from groups III in §
4) are improved by thetwo-line fits, although the degree of improvement varies.We also dropped lag estimates that were either flaggedas unreliable or believed to be wrong (i.e., the 6 H β lagsfrom groups IV and V in §
4) and keep only those objectsdeemed to give robust estimates of lag by our method(i.e., the 60 H β lags from groups I, II and III). To illus-trate the quality of the final result for each source, wedivide all 60 remaining sources into 4 new groups basedon the results of both the single-line fits in § § β lag. Seven of the objects havelight curves of lines other than H β to carry outtwo-line fits, but they provided little gain when thesingle-line fits already provided good lag estimates.(B) The 10 group II sources from § β lag estimate but potentially larger uncertaintiesdue to the presence of low significance ( > σ ) sub-peaks in the lag likelihood distribution. Most ofthose sources do not have the multiple line lightcurves needed to carry out two-line fits.(C) The four group III sources (NGC 3516, PG 0026,PG 0052, and PG 1307) from § § β light curves following the discussion in § ambiguity Fig. 8.—
Comparison of the H β time lags from CCF analysis and the SPEAR method, similar to Figure 4, updated where we have usedthe two-line fits and dropping the 6 unreliable sources. Four types of symbols are used to indicate our estimate for increasing levels ofambiguity in the lag estimate. Objects with inconsistent lag estimates between the two methods are labeled. red triangles correspond to sources of group A, B, C, andD, respectively. There is general agreement between thetwo methods, but also several discrepancies, as 7 of ourH β lag estimates are inconsistent with the CCF resultsgiven their error estimates. We marked these sources inFigure 8 and now discuss each case individually, NGC 7469. — We estimate an H β lag of 11.7 +0 . − . days,as opposed to τ CCF =4.7 +0 . − . . However, if we use aDirac delta function for the transfer function instead ofa tophat, the estimated time lag changes to 4.3 days,in agreement with the CCF result. Thus, the discrep-ancy originates from the improvement of fit with a tophatsmoothing kernel. The continuum of NGC 7469 was in-tensively monitored to search for time lags between theUV and optical continuum (Collier et al. 1998), so itscontinuum light curve is densely sampled while the H β light curve is much less so. The model has to smooththe continuum light curve heavily (i.e., a broad tophat width) to obtain a good fit, which at the same time shiftsthe time lag estimate to a longer value than it would bewith a zero width (i.e., a delta function). This is sugges-tive of the continuum errors being underestimated, or amore realistic transfer function is required. Mrk 79 (years 2 and 4) . — In both cases, we estimatelarger time lags than the CCF results, although there aresub-peaks which correspond to the CCF lags. For year2 (JD47838–48044), while the CCF centroid gives a lagof 16.4 +6 . − . days, the CCF peak estimate is 19 +11 − days,more consistent with our estimate of 30.9 +1 . − . days. Foryear 4 (JD49996–50220), Peterson et al. (2004) flaggedit as “unreliable” for the poor light curve sampling. Ourmethod shows a dense array of sub-peaks in the lag like-lihood distribution, but the most significant peak is at43.6 +1 . − . days.easuring Reverberation Lags 13 Fig. 9.—
The R BLR – L relation for H β . The luminosity is λL λ (5100 ˚A) and the BLR radius is equivalent to the lag in units of lightdays. The open symbols and gray solid circles indicate the measurement from SPEAR method and from CCF method for the same set ofsources, respectively. The gray solid curve is the fit to the CCF R BLR – L relation, while the rest of the curves are the fits to the SPEAR R BLR – L relation, using four subsets of the sources (see Table 2 for details of each fit). The slope of the fit to the SPEAR R BLR – L relation α is steeper than the CCF relation, but the two are consistent within the uncertainties. σ rms is the rms scatter of each fit. PG 0844. — As discussed in §
5, the CCF estimate ofthe H β lag (34.4 +14 . − . ) for PG 0844 is likely susceptibleto correlated errors, while our method estimates a lag of12.2 +5 . − . days regardless of the value of correlation coef-ficient r . PG 0052. — We estimate an H β lag of 149.3 +4 . − . days,as opposed to τ CCF =103 +28 . − . . The single-line fit showsmultiple peaks and usually one would be inclined to mis-trust a peak at the seasonal alias (a rest-frame lag of150 days corresponds to 170 days in the observed-frame).However, the joint H α /H β fit clearly reinforced this so-lution. PG 1307. — We estimate an H β lag of 188.8 +5 . − . days,as opposed to τ CCF = 121.9 +41 . − . . The joint H α /H β fitsuppressed the false peak which corresponds to the τ CCF lag, favoring a longer lag that is more consistent withlags of the other Balmer lines.
PG 1426. — We estimate an H β lag of 161.6 +6 . − . days,as opposed to τ CCF = 103.2 +32 . − . . Similar to PG 0052and PG 1307, the joint H α /H β fit reinforced a solutionwhich is otherwise susceptible to the seasonal gap effect.We carried out a similar analysis for each data set, in-cluding all emission lines besides H β , as summarized inTable 1. Note that in the table we only include 87 lightcurves for which we have successfully computed lags. Theobject is identified in column (1). The emission lineand its light curve Heliocentric Julian Date range arelisted in columns (2) and (3), respectively. Column (4)gives the rest-frame time lag estimate from the SPEARmethod, while column (5) indicates the associated “am-biguity” (i.e., the group membership) defined above.4 Zu et al. Fig. 10.—
Lag (left) and scaling coefficient (right) of the H β transfer function as a function of continuum luminosity from 14 years ofNGC 5548 data. The best-fit slopes are also reported for each panel and shown by the black solid lines. The black dashed line in the leftpanel is the best fit with a fixed slope of 0.50. The number inside each solid circle indicates the year of observation for each light curvestarting from Dec. 1988. Note that in the left panel the point for year 12 is hidden under that of year 13. THE R BLR – L RELATION FOR H β With the revised set of H β lags, and the starlight-corrected optical luminosity of each AGN fromBentz et al. (2009), we have calculated the fit to the R BLR – L relationship for our sample R BLR = C · (cid:18) L . ergs s − (cid:19) α (23)and compared it to that based on CCF lags in Figure 9.We obtained a slope α = 0 . ± .
010 for all the SPEARlags regardless of the level of “ambiguity” at which weprobe the slope. This slope is slightly steeper than pre-vious estimates, and only marginally consistent with thena¨ıve theoretical prediction of α = 0 .
5. Compared to theCCF-based R BLR – L fit of the same sample of AGNs (bluefilled circles), our R BLR – L fit has a steeper slope butcomparable rms scatter σ rms , which grows smaller as weuse more reliable lags. Table 2 gives the results from thedifferent fits using the 4 combinations of groups indicatedby column (1). Column (2) gives the number of datapoints used in each fit. We fit each combinatorial dataset using lag estimates from both the SPEAR (columns3–6) and CCF methods (columns 7–10). Two param-eters in Equation (23) are listed in columns (3) and(4) for SPEAR method, and in column (7) and (8) forCCF method, respectively. Column (5) and (6) give the χ /dof and rms scatter for our fit, while column (9) and(10) give these statistics correspondingly for the CCFmethod.Our R BLR – L fits have a larger χ /dof than the CCFones. This does not necessarily mean they are poorer fits,because our lag estimates generally have tighter error-bars than the CCF estimates. It could indicate that ourapproach underestimates uncertainties, that the CCFmethod overestimates uncertainties, or that we are nottaking into account intrinsic scatter in the R BLR – L re-lationship. Since most of the group C and D sources are high-redshift luminous PG objects, the rms scatterfor our method decreases from 0.229 dex to 0.196 dexafter dropping them from the fit. Three outliers fromthe CCF R BLR – L relation (NGC 7469, years 1 and 4of Mrk 79) are also the sources where our lag estimatesare inconsistent with CCF results. When we use our lagestimates, these three CCF outliers lie on the R BLR – L re-lation, which reduces the rms scatter near L opt ∼ . ergs s − . Note that there is significant scatter in the R BLR – L relation even for multiple estimates for a sin-gle source, as shown in the left panel of Figure 10 forNGC 5548. R BLR – L RELATION OF NGC 5548 REVISITED
So far, we have carried out our calculations assumingthat the parameters are constant during a season. Thisis likely true for the underlying variability process. If wemodel either the full continuum light curve or the individ-ual seasons, we find estimates for the process parameters τ d and ˆ σ that are statistically consistent. We do observelags that vary from season to season, and these are ar-guably correlated with luminosity. If so, they should alsobe varying within seasons, and we have not accountedfor this. Similarly, we assume the scaling between thecontinuum and line fluxes does not vary over a season,although we do observe it to vary between seasons.The nearby Seyfert 1 galaxy NGC 5548, with its manycontinuous years of monitoring data, serves as an idealexample of an AGN changing its variability levels fromseason to season. Figure 10 illustrates the continuumflux dependence of both the H β lag h τ i and the scal-ing coefficient A for 14 seasons of NGC 5548 data. Weclearly see trends that the lag increases with luminosityand the amplitude of the response diminishes. If we fitthe lag, we find a steep slope, h τ i ∝ L . ± . that isinconsistent with the expected h τ i ∝ L . . However, thepoor fit ( χ /dof = 4 .
17) suggests that either the uncer-tainties are underestimated or intrinsic scatter dominateseasuring Reverberation Lags 15
TABLE 1Rest-Frame Lag Estimates
Julian Dates τ SPEAR
Object Line ( − β . +2 . − . A3C 390.3 Ly α . +34 . − . D3C 390.3 C iv λ . +2 . − . CAkn 120 H β . +6 . − . AAkn 120 H β . +3 . − . BFairall 9 H β . +42 . − . DFairall 9 He ii λ . +0 . − . CFairall 9 Ly α . +0 . − . CMrk 79 H β . +2 . − . BMrk 79 H β . +1 . − . AMrk 79 H β . +7 . − . BMrk 79 H β . +1 . − . AMrk 110 H β . +2 . − . BMrk 110 H β . +6 . − . AMrk 110 H β . +2 . − . AMrk 279 H β . +1 . − . AMrk 290 H β . +0 . − . AMrk 335 H β . +3 . − . AMrk 335 H β . +3 . − . AMrk 509 H β . +0 . − . AMrk 509 He ii λ . +0 . − . DMrk 590 H β . +1 . − . AMrk 590 H β . +2 . − . AMrk 590 H β . +3 . − . BMrk 590 H β . +2 . − . AMrk 817 H β . +2 . − . AMrk 817 H β . +1 . − . AMrk 817 H β . +4 . − . AMrk 817 H β . +1 . − . ANGC 3227 H β . +6 . − . ANGC 3227 H β . +0 . − . ANGC 3516 H α . +0 . − . ANGC 3516 H β . +0 . − . CNGC 3516 H β . +1 . − . ANGC 3783 H β . +0 . − . ANGC 4051 H β . +0 . − . BNGC 4151 H β . +0 . − . ANGC 4593 H β . +0 . − . ANGC 7469 H β . +0 . − . ANGC 7469 Si iv λ . +0 . − . ANGC 7469 C iv λ . +0 . − . ANGC 7469 He ii λ . +0 . − . APG 0026+129 H α . +1 . − . BPG 0026+129 H β . +7 . − . CPG 0052+251 H α . +2 . − . APG 0052+251 H β . +4 . − . CPG 0052+251 H γ . +1 . − . CPG 0804+761 H α . +8 . − . CPG 0804+761 H β . +2 . − . APG 0804+761 H γ . +46 . − . DPG 0844+349 H α . +0 . − . APG 0844+349 H β . +5 . − . BPG 0844+349 H γ . +2 . − . CPG 0953+414 H β . +3 . − . APG 0953+414 H γ . +3 . − . DPG 1211+143 H α . +0 . − . CPG 1211+143 H β . +0 . − . DPG 1211+143 H γ . +15 . − . APG 1226+023 H α . +40 . − . BPG 1226+023 H β . +284 . − . DPG 1226+023 H γ . +9 . − . APG 1229+204 H α . +2 . − . APG 1229+204 H β . +2 . − . A TABLE 1—
Continued
Julian Dates τ SPEAR
Object Line ( − α . +4 . − . APG 1307+085 H β . +5 . − . CPG 1307+085 H γ . +7 . − . DPG 1411+442 H α . +10 . − . APG 1411+442 H β . +13 . − . BPG 1426+015 H β . +6 . − . APG 1617+175 H α . +9 . − . BPG 1617+175 H β . +31 . − . BPG 2130+099 H β . +4 . − . APG 2130+099 He ii λ . +3 . − . ANGC 5548 H β . +0 . − . ANGC 5548 H β . +0 . − . ANGC 5548 H β . +2 . − . ANGC 5548 H β . +1 . − . ANGC 5548 H β . +1 . − . ANGC 5548 H β . +1 . − . ANGC 5548 H β . +1 . − . ANGC 5548 H β . +0 . − . ANGC 5548 H β . +0 . − . ANGC 5548 H β . +1 . − . ANGC 5548 H β . +3 . − . ANGC 5548 H β . +1 . − . ANGC 5548 H β . +0 . − . BNGC 5548 H β . +1 . − . A Note . — Lag estimates and confidence limits for Groups A and Bare calculated by the single-line fits, while those for Groups C andD are from the two-line fits. the goodness-of-fit. If we rescale the uncertainties so thatthe best-fit model has χ /dof ≡
1, the flatter L . slopeis not ruled out, with a ∆ χ of only 0 . h τ i ∝ L α instead of shifting the entire lightcurve by the same t lag , and then optimize the fits overthe additional parameter α . Unfortunately, we cannotfit the full NGC 5548 light curve because the resultingmatrix dimensions are impractically high ( K =3085 datapoints). We instead estimate the normalized likelihooddistribution for α in each season and then combined thelikelihoods, as shown in Figure 11 (we did not includeyear 13, which was part of the less reliable group B).This “breathing” effect is clearly detected, and the log-arithmic slope estimate of α = 0 . +0 . − . is consistentwith the na¨ıve expectation α = 1 / R BLR – L re-lation in Figure 9. Using almost the same set of lightcurves from NGC 5548 (we add year 14 and excludeyear 13), Cackett & Horne (2006) find a much shallowerslope (0 . .
46) with a luminosity-dependent delay map,in better agreement with the prediction of photoioniza-tion models ( ∼ .
23; Korista & Goad 2004). However,their small correction for the host galaxy starlight mayartificially flatten their estimate of the slope (Bentz et al.2007). Note that for this experiment we did not makethe line-to-continuum scaling coefficient A a function ofcontinuum luminosity in the fit. Such a full scale calcu-lation should be carried out using the complete data set.6 Zu et al. TABLE 2BLR Size-Luminosity Relation
Groups
N C
SPEAR α SPEAR χ σ SPEARrms C CCF α CCF χ σ CCFrms
Included (lt-days) (dex) (lt-days) (dex)(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)A,B,C,D 60 1 . ± .
01 0 . ± .
010 8.13 0.229 1 . ± .
01 0 . ± .
017 3.50 0.203A,B,C 57 1 . ± .
01 0 . ± .
011 8.52 0.207 1 . ± .
01 0 . ± .
019 3.54 0.205A,B 53 1 . ± .
01 0 . ± .
013 8.39 0.206 1 . ± .
01 0 . ± .
020 3.81 0.213A 43 1 . ± .
01 0 . ± .
013 9.44 0.196 1 . ± .
01 0 . ± .
020 4.22 0.211
Fig. 11.—
Likelihood distribution of α for 13 years of NGC 5548light curves. The normalized log-likelihood is calculated by addingthe likelihood distribution functions for the 13 individual yearstogether. DISCUSSION
We have demonstrated that direct fitting of contin-uum and line light curves is a viable approach to mea-suring reverberation lags, confirming the initial study ofRybicki & Kleyna (1994). It provides a full statisticalframework for determining time lags and estimating theiruncertainties, including the full contributions from cor-related noise, de-trending and interpolation. In essence,the lags are determined using a weighted average of allstatistically acceptable models for interpolating the un-derlying true light curve. While we used the assump-tion that the underlying variable process had an expo-nential correlation function corresponding to a dampedrandom walk, any other statistical process could be sub-stituted. We note, however, that Kelly et al. (2009),Koz lowski et al. (2010) and MacLeod et al. (2010) havefound the exponential correlation function to be an excel-lent model of quasar light curves, just as we have foundhere, although we modeled the light curves in flux ratherthan magnitude.Because we are explicitly modeling the light curves, wemust include an explicit model of the transfer function.Here we used a top hat for simplicity. It includes the sim-ple limits of a delta function and a uniform thin shell, andis likely a reasonable model for any single-peaked transferfunction given the available data (see Rybicki & Kleyna1994). As with the model for the variability process, us-ing an alternative transfer function simply requires com-puting the appropriate terms of the covariance matrix. Aside from the case of NGC 7469 where it seemed toaffect the lag estimation, we did not discuss the tophatwidth. In general, there is a relatively strong degeneracybetween ∆ t , the width, and A , the scaling between thecontinuum and line light curves. When ∆ t is large andthe continuum is heavily smoothed, the model will try toincrease the variability amplitude by artificially boosting A to re-align the continuum and line light curves. How-ever, the degeneracy does not seem to lead to problemsin estimating the mean lag unless the line light curve isvery poorly sampled. The traditional CCF method doesnot implicitly assume a shape for the transfer functionbut calculates the lag as either the barycenter ( τ cent )or the peak ( τ peak ) of underlying transfer function (con-volved with data). The difference between the two some-times can be large and hard to reconcile unless the trans-fer function can be modeled explicitly. For future high-fidelity datasets, our approach should also have no diffi-culty constraining the shape of transfer functions.The most important future path for this method is tosimultaneously fit multiple line components, whether dif-ferent lines (e.g., H β , H α , etc.), velocity sub-componentsof individual lines or multiple continuum bands. As longas the overall ansatz that all light curves are scaled andsmoothed versions of the continuum holds, combiningmany light curves with differing lags means that the lagestimate for any given light curve is now derived froma better sampled estimate of the continuum variabil-ity. A second advantage, particularly for attempts tostudy the velocity structure of a particular broad line,is that such joint analyses will correctly infer the co-variances between the individual lags. Current velocity-dependent lags have uncertainties comparable to theirdifferences (Bentz et al. 2008, 2010; Denney et al. 2009,2010), but it may be true that these differences actu-ally have a strong covariances, so that the differencesare far more significant than estimates from analyzingthe light curves in isolation. The method can also allowfor luminosity-dependent lags or line-continuum scalingfactors. Also note that while we only use the linear pa-rameters of the model to remove the light curve means, itis a very flexible tool for de-trending or cross-calibratinglight curves whose model uncertainties will be fully in-cluded in lag estimate.The most important observational implication of thisapproach is the value of measuring multiple lines, es-pecially those with high ionization potentials. In ourapproach, multiple lines with differing lags allow oneto overcome many of the sampling problems inherentto cross-correlation methodologies. At its simplest, onelight curve can be aliased into a (seasonal) sampling gap,but two cannot be unless the transfer functions are sim-ilar (i.e., the lines have similar lags). Given the radialeasuring Reverberation Lags 17ionization stratification of the BLR (Peterson 1993), thelag difference between two lines is proportional to thedifference in their ionization levels. In this paper, how-ever, the lines we used for two-line fits are mostly pairsof two Balmer lines, which have similarly low ionizationlevels. Thus, the observational goal should be to obtaindata for multiple lines with a broad range of ionizationpotentials. Indeed, with a wide variety of emission-linelines, it is in principle possible to combine the reverber-ation results with photoionization equilibrium modelingto highly constrain the geometry and physics of the BLR(Horne, Korista, & Goad 2003).The only significant algorithmic challenge comes fromthe O ( K ) scaling of the computational cost with thenumber of data points K . Unfortunately, the reverber-ation mapping problem is very different from simply us-ing the damped random walk to model the continuumlight curves, where we can take advantage of the partic- ular structure of the covariance matrix to calculate thenecessary matrix inversions in O ( K ) operations. Sincethe expensive matrix inversion is required for each like-lihood calculation, it becomes difficult to analyze largedata sets, particularly if the number of parameters alsoincreases greatly as in a full simultaneous model of lagsas a function of line velocity. These problems can be ad-dressed using hyper–threaded or parallel versions of theunderlying algorithm. ACKNOWLEDGEMENTS
We thank Kelly D. Denney and Catherine J. Grier forkindly providing some of the light curves. Thanks also toMisty C. Bentz for her starlight corrected AGN luminosi-ties. CSK is supported by NSF Grant AST-0708082 andAST-1009756. BMP is supported by NSF Grant AST-1008882 and YZ is supported by an OSU DistinguishedUniversity Fellowship.
APPENDIX
COVARIANCE MATRIX OF THE CORRELATION FUNCTIONS
The expressions for the covariance matrices used in this paper and the accompanying code assume that the transferfunction is a simple top hat, Ψ( t − t ′ ) = A ( t − t ) − for t ≤ t − t ′ ≤ t . (A1)For this transfer function, we can analytically calculate the correlation functions in Equation (6), (7) and (8), respec-tively. The Covariance Matrix Between the Continuum and One Line
The covariance between continuum s c ( t ) at t j and line s l ( t ) at t i with transfer function defined as in Equation (A1)is h s c ( t j ) s l ( t i ) i = τ d σ A e − t L /τ d − e − t H /τ d if t L > e t H /τ d − e t L /τ d if t H < − e t L /τ d − e − t H /τ d if t L ≤ ≤ t H , (A2)where t L ≡ t i − t j − t and t H ≡ t i − t j − t . The Covariance Matrix Between Two Lines
Consider the case when the first line s l ( t ) has transfer function Ψ( t − t ′ ) as defined in Equation (A1) and the otherline s ′ l ( t ) has transfer function Ψ ′ ( t − t ′ )Ψ ′ ( t − t ′ ) = B ( t − t ) − for t ≤ t − t ′ ≤ t , (A3)where t − t ≤ t − t . The covariance between line s l ( t ) at time t i and line s ′ l ( t ) at time t j (Equation 8) is h s l ( t i ) s ′ l ( t j ) i = τ σ A B e −| t L | /τ d + e −| t H | /τ d − e −| t M | /τ d − e −| t M | /τ d + ( t H /τ d if t M ≤ < t H t − t ) /τ d if t M ≤ < t H − t L /τ d if t L ≤ < t M e −| t L | /τ d + e −| t H | /τ d − e −| t M | /τ d − e −| t M | /τ d if t L > t H < , (A4)where t L ≡ ( t i − t j ) − ( t − t ) ,t M ≡ ( t i − t j ) − ( t − t ) ,t M ≡ ( t i − t j ) − ( t − t ) , and t H ≡ ( t i − t j ) − ( t − t ) . (A5)By definition, the covariance for the autocorrelation of line s l ( t ) between time t i and t j (Equation 7) can be obtainedby equating Ψ ′ ( t − t ′ ) with Ψ( t − t ′ ) so that B ≡ A , t ≡ t and t ≡ t . REFERENCESAlexander, T. 1997, in Astronomical Time Series, ed. Maoz, D.,Sternberg, A., & Leibowitz, E. M. (Dordrecht: Kluwer), p. 163 Bentz, M. C., Peterson, B. M., Pogge, R. W., Vestergaard, M., &Onken, C. A. 2006, ApJ, 644, 1338 Zu et al.