AGN Variability

# Are the Variations in Quasar Optical Flux Driven by Thermal Fluctuations?

Brandon C. Kelly11affiliation: bckelly@cfa.harvard.edu 22affiliation: Hubble Fellow 33affiliation: Harvard-Smithsonian Center for Astrophysics, 60 Garden St, Cambridge, MA 02138 44affiliation: Department of Astronomy, University of Arizona, Tucson, AZ 85721 , Jill Bechtold44affiliation: Department of Astronomy, University of Arizona, Tucson, AZ 85721 , Aneta Siemiginowska33affiliation: Harvard-Smithsonian Center for Astrophysics, 60 Garden St, Cambridge, MA 02138
###### Abstract

We analyze a sample of optical light curves for 100 quasars, 70 of which have black hole mass estimates. Our sample is the largest and broadest used yet for modeling quasar variability. The sources in our sample have , Å, and . We model the light curves as a continuous time stochastic process, providing a natural means of estimating the characteristic time scale and amplitude of quasar variations. We employ a Bayesian approach to estimate the characteristic time scale and amplitude of flux variations; our approach is not affected by biases introduced from discrete sampling effects. We find that the characteristic time scales stongly correlate with black hole mass and luminosity, and are consistent with disk orbital or thermal time scales. In addition, the amplitude of short time scale variations is significantly anti-correlated with black hole mass and luminosity. We interpret the optical flux fluctuations as resulting from thermal fluctuations that are driven by an underlying stochastic process, such as a turbulent magnetic field. In addition, the intranight variations in optical flux implied by our empirical model are mag, consistent with current microvariability observations of radio-quiet quasars. Our stochastic model is therefore able to unify both long and short time scale optical variations in radio-quiet quasars as resulting from the same underlying process, while radio-loud quasars have an additional variability component that operates on time scales day.

accretion, accretion disks — galaxies: active — methods: data analysis — quasars: general

## 1. Introduction

It is widely accepted that the extraordinary activity associated with quasars555Throughout this work we will use the terms quasar and AGN to refer generically to broad line AGNs. No luminosity difference between the two is assumed. involves accretion onto a supermassive black hole, with the UV/optical emission arising from a geometrically thin, optically thick cold accretion disk. Aperiodic variability across all wavebands is ubiquitous in AGN, with the most rapid variations occuring in the X-rays (for a review, see Ulrich et al., 1997). The source of quasar variability is unclear, and several models have been proposed for describing the optical variability of quasars, including accretion disk instabilities (e.g., Kawaguchi et al., 1998), supernovae (e.g., Aretxaga et al., 1997), microlensing (Hawkins, 2000), and more general Poisson process models (e.g., Cid Fernandes et al., 2000). However, recent results from reverberation mapping have shown that the broad emission lines respond to variations in the continuum emission after some time lag (e.g., Peterson et al., 2004), implying that the continuum variations are dominated by processes intrinsic to the accretion disk. If the optical/UV variations are intrinsic to the accretion disk, then thermal fluctuations appear to be a natural choice for driving the optical/UV variations, as the optical/UV emission is thought to be thermal emission from the accretion disk. The fact that quasars become bluer as they brighten is consistent with a thermal origin (e.g., Giveon et al., 1999; Trèvese et al., 2001; Geha et al., 2003). Moreover, the thermal timescale is sensitive to the disk viscosity. Variability is therefore a potentially important and powerful probe of the quasar central engine and accretion disk physics.

A considerable amount of our interpretation and understanding of AGN accretion disks, and consequently their optical/UV emission, is based on the so-called -prescription (Shakura & Syunyaev, 1973). Within the standard model, the viscosity, thought to be the source of the thermal disk emission and outward transfer of angular momentum, is parameterized as being proportional to the total pressure in the disk. Previous work on estimating from quasar variability has found a value of (e.g., Siemiginowska & Czerny, 1989; Collier & Peterson, 2001; Starling et al., 2004). However, when the disk is dominated by radiation pressure, as is thought to be the case in the inner regions, an -disk is both thermally and viscously unstable (Shakura & Sunyaev, 1976; Lightman & Eardley, 1974). In particular, for a radiation pressure dominated disk, the thermal instability is expected to grow exponentially on a time scale similar to the thermal time scale. For AGN, the thermal time scale is on the order of months to years, and in general there is no evidence for instabilities in the optical light curves of AGN which span (e.g., Collier & Peterson, 2001). However, the optical light curves for a few sources may show evidence of instability on a thermal time scale (e.g., Czerny et al., 2003; Lub & de Ruiter, 1992).

If the stress is not proportional to the total pressure, but rather some other combination of the gas and radiation pressure, then the disk becomes more stable to thermal perturbations (e.g., Stella & Rosner, 1984; Szuszkiewicz, 1990; Merloni & Fabian, 2002; Merloni, 2003). Alternatively, if part of the accretion energy is dissipated in a hot corona, then there is no longer any runaway heating in the disk, as additional cooling is provided by the corona (Svensson & Zdziarski, 1994). For example, disk models with alternative prescriptions for , as well as the addition of a corona and jet, have been able to reproduce the variability of the microquasar GRS 1915+105 (Nayakshin et al., 2000; Janiuk et al., 2000, 2002). Previous work on 3-dimensional magneto-hydrodynamic (MHD) simulations of radiation pressure-dominated AGN accretion disks have not observed a thermal instability (e.g., Turner, 2004; Hirose et al., 2008). Moreover, the most promising physical mechanism behind the viscous torque is the magneto-rotational instability (MRI, Balbus & Hawley, 1991, 1998), and recent numerical and analytical work has suggested that the prescription may be a poor representation for MRI-driven viscosity (Pessah et al., 2008). Variability studies can therefore lend observational insight into the inadequacies of the traditional -prescription, and possibly lead to a more accurate characterisation of the relationship between stress and pressure in the accretion disk.

A successful model for quasar X-ray variability describes the X-ray variations on long time scales as being the result of perturbations in the accretion rate that occur outside of the X-ray emitting region (e.g., Lyubarskii, 1997; Mayer & Pringle, 2006; Janiuk & Czerny, 2007). These accretion rate perturbations then travel inward, modulating the X-ray emitting region. It has been suggested that the origin of such perturbations is the result of a magnetic field randomly varying in time, and may be related to the appearance of an outflow (King et al., 2004). If this model for the X-ray variability is correct, we would expect to also see variations in the optical luminosity, whose origin lies in the disk at radii farther from the central source. Therefore, understanding the origin of quasar optical variations will not only lead to a better understanding of accretion disk physics, but may also lead to a better understanding of the origin of quasar X-ray variability, possibly unifying the source of variability in the two bands.

Quasars have also been observed to vary on time scales as short as hours (so-called ‘microvariability’ or ‘intranight variability’, Gopal-Krishna et al., 2003; Stalin et al., 2004, 2005; Gupta & Joshi, 2005; Carini et al., 2007). Microvariability is known to be stronger in radio-loud quasars, especially blazars (e.g., Gupta & Joshi, 2005), while microvariability in radio-quiet quasars is usually not detected above the photometric uncertainty (e.g., Gupta & Joshi, 2005; Carini et al., 2007). The enhanced level of microvariability in radio loud objects suggests that it is due to processes in a jet, while the physical cause of microvariability in radio quiet objects has remained a puzzle. Reprocessing of X-rays, disk instabilities, or a weak blazar component have been suggested (Czerny et al., 2008).

There have been numerous previous investigations of quasar optical variability. However, because of the difficulty in obtaining high quality, well sampled light curves that cover a long time span, most previous work has involved ensemble studies of quasars, or analysis of simple correlations involving variability amplitude. The most well known result from previous work is a tendency for AGN to become less variable as their luminosity increases (e.g., Hook et al., 1994; Garcia et al., 1999; Giveon et al., 1999; Geha et al., 2003; Vanden Berk et al., 2004; de Vries et al., 2005, and references therein). There have also been claims of a variability–redshift correlation (e.g., Cristiani et al., 1990; Cid Fernandes et al., 1996; Cristiani et al., 1996; Trèvese & Vagnetti, 2002), although the sign of this correlation varies between studies (e.g., see the list in Giveon et al., 1999). If real, the variability–redshift correlation is most likely caused by the fact that quasars are more variable at shorter wavelengths (e.g., Cutri et al., 1985; di Clemente et al., 1996; Helfand et al., 2001; Vanden Berk et al., 2004), corresponding to regions in the disk closer to the central black hole. Recently, a correlation between optical variability and black hole mass has been claimed (Wold et al., 2007). While many of these correlations are formally statistically significant, they often exhibit considerable scatter when one measures variability for individual objects.

A few previous studies have employed spectral techniques, such as power spectra and structure functions, in the analysis of quasar optical light curves. From these studies it has been inferred that quasar optical light curves generally have variations of on timescales of months, and that the power spectra of optical light curves is well described as (Giveon et al., 1999; Collier & Peterson, 2001; Czerny et al., 2003); power spectra of this form are consistent with random walk, or more generally, autoregressive processes. In addition Collier & Peterson (2001) analyzed a sample of optical light curves from 8 low- Seyfert 1 galaxies. They found that the characteristic time scales of optical variations for the AGN in their sample are days and correlate with black hole mass, consistent with disk orbital or thermal time scales. Czerny et al. (1999) found evidence for a flattening of the optical power spectrum of NGC 5548 on time scales longer than days. Both Czerny et al. (1999) and Czerny et al. (2003) also suggested that the long time scale optical variability was due to thermal fluctuations in the accretion disks of NGC 5548 and NGC 4151, respectively.

Motivated by the potential in variability studies for increasing our understanding of the structure of quasar accretion disks, we have compiled a sample of well-sampled optical light curves from the literature. We directly model the quasar optical light curves as a stochastic process, in contrast to previous work based on more traditional Fourier (i.e., spectral) techniques. Our method allows us to describe quasar light curves with three free parameters: a characteristic time scale, amplitude of short time scale variability, and the mean value of the light curve. In addition, our method enables us to estimate the characteristic time scale of quasar variations without the windowing effects that can bias spectral approaches. Our sample consists of 100 AGN, 70 of which have black hole mass estimates. The sources in our sample have , Å, and , making this by far the largest sample yet used for this kind of study.

In this work we adopt a cosmology based on the WMAP results (, Spergel et al., 2003).

## 2. Data

In this work we analyze a sample of 100 quasar optical light curves, compiled from the literature. Our sample consists of 55 AGN from the MACHO survey (Geha et al., 2003), 37 Palomar Green (PG) quasars from the sample of Giveon et al. (1999), and 8 Seyfert galaxies from the AGN Watch database. We were able to obtain black hole mass estimates for 71 of the AGN, where has been estimated for 20 of them from reverberation mapping (Peterson et al., 2004), and is estimated for the remaining 51 AGN from the broad emission lines using standard scaling relationships (e.g., Vestergaard & Peterson, 2006). Our sample is summarized in Table 1.

### 2.1. Macho Quasars from Geha et al.(2003)

We collected band light curves from AGN selected by the MACHO survey based on their variability (Alcock et al., 1997, 1999). The motivation for the MACHO survey was to study Galactic microlensing events behind the Magellanic Clouds. However, the survey was also able to select quasars via their variability, confirmed via spectroscopic follow-up, producing 59 quasars with well-sampled light curves over a broad range in redshift () (Geha et al., 2003). The quasar light curves span years, and have a typical sampling interval of days, although much longer gaps exist for some light curves. In general, the MACHO light curves are the highest quality in our sample, often being frequently and regularly sampled. However, we note that because the MACHO sources were selected based on their variability, this sample is biased toward more variable objects. See Geha et al. (2003) for more details of the MACHO quasar catalogue.

Spectra for the MACHO quasars are presented in Geha et al. (2003), and spectra for 27 of these quasars were kindly provided to us by Marla Geha. We calculated black hole mass estimates from the source luminosity and the of the H, Mg II, or C IV broad emission line using standard scaling relationships (e.g., Vestergaard & Peterson, 2006, Vestergaard 2008). However, the spectra are not flux-calibrated, so the luminosities at and Å were estimated from the photometric data, assuming a power-law continuum with (Richards et al., 2001). Typical uncertainties on the broad line mass estimates are dex (e.g., Vestergaard & Peterson, 2006). When combined with the measurement error in the measurements, and the uncertainty due to the intrinsic scatter in quasar spectral slopes (, Richards et al., 2001), our typical adopted uncertainty in was dex for the MACHO quasars.

### 2.2. PG Quasars from Giveon et al.(1999)

Giveon et al. (1999) report and photometric light curves for 42 PG quasars spanning years with a typical sampling interval of days. These quasars are all bright ( mag) and nearby (). Peterson et al. (2004) report estimates of and Å calculated from reverberation mapping for 12 of these quasars. Estimates of and Å were taken from Vestergaard & Peterson (2006) for the remaining Giveon et al. (1999) quasars. See Giveon et al. (1999) for further details.

### 2.3. Seyfert Galaxies from AGN Watch Database

The light curves for the remaining 8 AGN in our sample are from the AGN watch project. The optical light curves for the Seyfert Galaxies AKN 564, Fairall 9, MRK 279, MRK 509, NGC 3783, NGC 4051, NGC 4151, NGC 5548, and NGC 7469 were taken from the AGN watch website. In general, we used the light curves at 5100Å, with the exception of AKN 564, for which we used the -band light curve. We excluded 3C390.3 from our analysis because it is a classified as an optically-violent variable, and thus may experience a different variability mechanism than the other AGN in our sample. Black hole masses and 5100Å luminosities were taken from Peterson et al. (2004). Further details may be found in the references listed in Table 1. Although the AGN Watch sample is small compared to the other two, these sources are particularly important in our analysis because they help to anchor the low- and low- end of the correlations analyzed in § 4.

## 3. The Statistical Model: Continuous Autoregressive Process

The previous analysis of Giveon et al. (1999) and Collier & Peterson (2001) suggests that the optical variability can be well-described by a power-law power spectrum with slope . The power spectra of the MACHO quasars have not been analyzed for individual objects, although Hawkins (2007) investigated the power spectra of the ensemble of objects. In Figure 1 we show the geometric mean -band power spectra for the MACHO quasars, along with the 90% confidence region on the geometric mean. The power spectra for the MACHO sources are well described by a power law of slope , consistent with previous work. The lack of any peaks in the power spectra, as well as the aperiodic and noisy appearance of quasar light curves, suggests that quasar light curves are stochastic or chaotic in nature.

In this work we model quasar light curves as a continuous time first order autoregressive process (CAR(1)). Power spectra of the form are consistent with a first order autoregressive (AR(1)) process. We model this stochastic process in continuous time, both because the actual physical processes in the accretion disk are continuous, and because doing so allows a natural way of handling the irregular sampling of our light curves. Moreover, this process is well-studied, only has three parameters, and provides a natural and consistent way of estimating a characteristic time scale and variance of quasar light curves.

The CAR(1) process is described by the following stochastic differential equation222Strictly speaking, the stochastic differential equation is complicated by the fact that white noise does not exist as a derivative in the usual sense. However, we ignore the mathematical technicalities for ease of interpretation of Equation (1) (e.g., Brockwell & Davis, 2002):

 dX(t)=−1τX(t)dt+σ√dtϵ(t)+b dt,   τ,σ,t>0. (1)

Here, is called the ‘relaxation time’ of the process , and is a white noise process with zero mean and variance equal to one. Within the context of this work, is the quasar flux. We assume that the white noise process is also Gaussian. The mean value of is and the variance is . Further details on the CAR(1) process are described in the Appendix.

The relaxation time, , can be interpreted as the time required for the time series to become roughly uncorrelated, and can be interpreted as describing the variability of the time series on time scales short compared to . Within the context of this work, is the quasar light curve. It is tempting to associate with a characteristic time scale, such as the time required for diffusion to smooth out local accretion rate perturbations, and to represent the variability resulting from local random deviations in the accretion disk structure, such as caused by turbulence and other random magneto-hydrodynamic (MHD) effects.

The power spectrum of a CAR(1) process is

 PX(f)=2σ2τ21+(2πτf)2. (2)

From Equation (2) we infer that there are two important regimes for : for and for . Therefore, the CAR(1) process has a power spectrum that falls off as at time scales short compared to the relaxation time, and flattens to white noise at time scales long compared to the relaxation time. Because ‘characteristic’ time scales of quasar light curves are often defined by a break in the power spectrum, this is an additional justification of associating with a characteristic time scale. In addition, because the power spectra of quasar optical light curves are well described by , it suggests that a CAR(1) process should provide a good description of the light curves, with being on the order of the length of the light curves or longer.

To illustrate the CAR(1) process, we simulate four CAR(1) light curves. The light curves were simulated by first simulating a random variable from a normal distribution with mean and variance ; note that this is the mean and variance of the CAR(1) process. Then, from this random initial value, we simulated the rest of the light curve using Equations (A4) and (A5) in the Appendix. These simulated light curves span a length of 7 yrs and are sampled every 5 days. The simulated light curves span a period in time similar to the quasar light curves analyzed in this work, but are better sampled than most of the quasar light curves. Three characteristic time scales of interest for quasars are the light crossing time, the gas orbital time scale, and the accretion disk thermal time scale. These time scales are

 tlc = 1.1×(MBH108M⊙)(R100RS) days (3) torb = 104×(MBH108M⊙)(R100RS)3/2 days (4) tth = 4.6×(α0.01)−1(MBH108M⊙)(R100RS)3/2 yrs, (5)

where is the mass of the black hole, is the emission distance from the central black hole, is the Schwarzschild radius, and is the standard disk viscosity parameter. For the simulated quasar light curves, we use and , and set equal to each of these three time scales. In addition, we use and . Note that the assumed value of only affects the thermal time scale, and a higher value of results in a shorter time scale.

The simulated light curves are shown in Figure 2, and their corresponding power spectra are shown in Figure 3. The increased amount of variation on long time scales with increasing is apparent. In addition, because is shorter than the time sampling, the first simulated light curve is only sampling frequencies on the flat part of the power spectrum, giving it the appearance of white noise. In contrast, the two simulated light curves with the longest time scales are sampled on the part of the power spectrum, giving them more of a ‘red noise’ appearance. In addition, ‘red noise’ leak affects the estimated power spectrum of the light curve with yrs, evidenced by the constant offset between the true power spectrum and the estimated one. Red noise leak occurs when power from time scales longer than the span of the time series ‘leaks’ into the shorter time scales, biasing the power spectrum when estimated as the modulus of the discrete Fourier transform (e.g., van der Klis, 1997).

### 3.1. Estimating the Parameters of a CAR(1) Process

The parameters for a CAR(1) process are commonly estimated by maximum-likelihood directly from the observed time series. This is an advantage over non-parameteric approaches, such as the discrete power spectrum or the structure function. The observed power spectrum and structure function can both suffer from windowing effects caused by the finite duration and sampling of the light curve, whereby power from high frequencies can leak to low frequencies (aliasing), and power at low frequencies can leak to high frequencies (e.g., red noise leak). For ground based optical observations, an additional complication is the periodicity enforced in the sampling caused by the Earth’s rotation around the sun, as objects are only observable during certain times of the year. For example, this periodic sampling can be seen in the light curve for the MACHO source shown in Figure 4. All of these effects can bias the power spectrum or structure function when estimated directly from the light curve in a non-parameteric fashion. Similarly, these effects can bias parameter estimates when fitting a model to the observed power spectrum or structure function.

In contrast, estimating a ‘characteristic’ time scale and variance directly from the observed time series, instead of from the observed power spectrum or structure function, has the advantage of being free of windowing effects, giving unbiased estimates of and . Of course, this requires one to assume a parameteric model for the time series (or power spectrum), but as we will show below the CAR(1) process provides a good description of most of the AGN light curves analyzed in this work. Furthermore, higher-order terms can be added to Equation (1) to allow additional flexibility (e.g., Brockwell & Davis, 2002), but this is beyond the scope of the current work.

When the light curve is measured with error, the likelihood function can be calculated using a ‘state-space’ representation of the time series (e.g., Brockwell & Davis, 2002). Denoting the measured fluxes as , observed at times with measurement error variances , the likelihood function, , is a product of Gaussian functions:

 p(x1,…,xn|b,σ,τ) = n∏i=1[2π(Ωi+σ2i)]−1/2 (6) ×exp{−12(^xi−x∗i)2Ωi+σ2i} x∗i = xi−bτ (7) ^x0 = 0 (8) Ω0 = τσ22 (9) ^xi = ai^xi−1 (10) + aiΩi−1Ωi−1+σ2i−1(x∗i−1−^xi−1) Ωi = Ω0(1−a2i) (11) + a2iΩi−1(1−Ωi−1Ωi−1+σ2i−1) ai = e−(ti−ti−1)/τ. (12)

The maximum-likelihood estimate is then found by maximizing Equation (6) with respect to and .

In this work we employ a Bayesian approach in order to directly compute the probability distribution of and , given our observed light curves. The probability distribution of the parameters, given the observed data (i.e., the posterior distribution), is calculated as the product of the likelihood function with a prior probability distribution. In this work we assume a uniform prior on and . For deriving a prior on , we note that when the data are regularly sampled, the CAR(1) process reduces to the AR(1) process described by Equation (A1) in the Appendix with , where is the time sampling interval. In this work we assume a uniform prior on . For an AR(1) process, gives the correlation between and . Therefore, we consider it reasonable to assume that any value of is a priori likely, i.e., we do not assume anything a priori about the correlations between subsequent data points, and therefore we assume a uniform prior on from 0 to 1. This prior is non-informative in the sense that all of the information on comes from the data.

The accuracy of the fit can be assessed by comparing the residuals of the light curve with the values expected under the assumption of a CAR(1) process. Equation (6) implies that if the CAR(1) process provides a good model of the observed data, then the residuals should be uncorrelated and follow a normal distribution:

 χ≡x∗i−^xi√Ωi+σ2i∼N(0,1). (13)

Here, the notation means that is distributed according to a normal distribution with mean equal to zero and variance equal to one. The goodness of fit can then be assessed by inspecting a plot of the residuals with time to ensure that they are uncorrelated, and by comparing a histogram of the residuals with the expected standard normal distribution.

### 3.2. Fitting the Quasar Light Curves

In this work we model the logarithm of the flux as following a CAR(1) process, or equivalently the apparent magnitudes. We do this because the assumption of a Gaussian white noise process in Equation (1) produces both positive and negative values of , while flux is a strictly positive quantity. The logarithm function maps a strictly positive quantity to the interval , and therefore Equation (1) is likely to be a better description of the light curve for the source apparent magnitude, as opposed to the source flux. Uttley et al. (2005) have also argued for describing accreting black holes as a Gaussian stochastic process in the logarithm of the flux.

An additional correction is needed to correct for cosmological time dilation. Because the quasars are fit in the observed frame, and time scales decrease as , spurious correlations may arise if one does not correct to the quasar rest frame. This is particularly problematic when dealing with flux limited samples, which create an artificial correlation between and both luminosity and . Noting that , Equation (1) can be expressed in the forms

 dX(t) = −1τobsX(t)dtobs+σobs√dtobsϵ(t)+bobs dtobs = −1+zτobsX(t)dtrest+σobs√(1+z)dtrestϵ(t) + (1+z)bobs dtrest = −1τrestX(t)dtrest+σrest√dtrestϵ(t) + brest dtrest. (16)

From Equations (3.2)–(16) it is apparent that the observed and rest frame parameters are related as

 τrest = (1+z)−1τobs (17) σrest = (1+z)1/2σobs (18) brest = (1+z)bobs. (19)

Noting that the mean of the CAR(1) process is , and that the variance is , Equations (17)–(19) imply that the mean and variance of a CAR(1) process are unaffected by cosmological time dilation. However, the variance observed from a light curve with a finite duration and sampling is still affected by time dilation, as the observed variance over a time interval is the integral of Equation (2) over that time interval. In what follows, the quantities and will always refer to the quasar rest frame quantities, unless specified otherwise.

Random draws of and from the posterior probability distribution are obtained using a Metropolis-Hastings algorithm. We calculated an estimate of each parameter as the median of the posterior, and the posterior medians, standard deviations, and confidence intervals are computed using these random draws. In general the posterior median values were not significantly different from a maximum likelihood fit. The goodness of fit for the light curves was determined by examining a histogram of the residuals and a plot of the residuals against time, as described in § 3.1. Occasionally outlying values of the flux are present for the light curves with more data points, possibly due to unidentified systematic error. These outlying data points were removed and the light curves were refit. In Figure 4 we show the light curve and best fit CAR(1) model for the most densely sampled object in our sample, NGC 5548, for a representative light curve from the MACHO sample, and for a representative light curve from the PG sample of Giveon et al. (1999). In general, the CAR(1) model provided a good fit to the quasar light curves analyzed in this work, and we only flagged 9 out of 109 light curves as having a bad fit. The best fit CAR(1) parameters, along with their uncertainties, for the 100 AGN in our sample are listed in Table 2. The 9 objects for which the fit was deemed unacceptable are not used in the regression analysis.

## 4. Results

The distribution of and light curve standard deviations for our sample are shown in Figure 5. The light curve standard deviation is calculated as the square root of the light curve variance, . The best fit quasar relaxation times have a median value of 540 days and a dispersion of 0.64 dex, and show typical long time scale optical variations of per cent. However, the uncertainties on are large and make a considerable contribution to the observed scatter. Correcting for the contribution from the uncertainties implies an intrinsic dispersion in relaxation times of dex.

In order to look for correlations of quasar variability properties on luminosity, redshift, black hole mass, and Eddington ratio, we used the linear regression method of Kelly (2007). The method of Kelly (2007) takes a Bayesian approach to linear regression, and accurately accounts for intrinsic scatter in the regression relationships, as well as measurement errors in both the dependent and independent variables. Measurement errors can be large for both the estimates of and , and thus can have a significant effect on the observed correlations (e.g., Kelly & Bechtold, 2007). Therefore, it is necessary to correct for the measurement errors when attempting to recover any underlying trends.

### 4.1. Dependence of Quasar Variability on Luminosity and Redshift

In order to investigate whether quasar variability properties depend on luminosity, redshift, or both, we performed a regression analysis. Throughout this section, the luminosity will always be taken to be at 5100Å. There has been considerable debate over whether quasar variability is correlated with luminosity or redshift, and the artificial correlation between the two has made it difficult for previous work to uncover the true intrinsic correlation. However, because we perform a linear regresion of variability properties on both and simultaneously, we are able to break the degeneracy between and . This is because the multiple linear regression describes how variability depends on at a given , and likewise for at a given .

In Figure 6 we show the relaxation time as a function of luminosity and redshift, and in Figure 7 we show as a function of luminosity and redshift. The results of the regressions for are

 logτ = (−10.29±3.76)+(0.29±0.08)logλLλ  [days] (20) logτ = (2.32±0.10)+(1.12±0.41)log(1+z)  [days] (21) logτ = (−8.13±0.12)+(0.24±0.12)logλLλ (22) + (0.34±0.58)log(1+z)  [days],

and the results of the regression for are

 logσ2 = (4.73±2.34)−(0.19±0.05)logλLλ  [R mag2/day] (24) logσ2 = (−3.84±0.06) (25) − (0.32±0.25)log(1+z)  [R mag2/day] logσ2 = (8.00±3.29)−(0.27±0.07)logλLλ (26) + (0.47±0.33)log(1+z)  [R mag2/day].

There is a statistically significant correlation between and both and , where the light curve relaxation time increases with increasing and . In addition, there is a statistically significant anti-correlation between and , implying that the short time scale variance decreases with increasing . However, the multiple regression results show that the luminosity trends are the dominant ones, and that there is no significant trend between either or and , at a given . This therefore implies that the observed trends with redshift are caused by the artificial correlation between and resulting from selection effects.

We also looked for trends of the light curve variance, , with and and found statistically significant evidence for a correlation between the variance of the light curve and . Both the Spearman and Kendall rank correlation statistic was significant at , although there is considerable scatter in the correlation. This correlation is most likely a reflection of the well-known fact that quasar emission at shorter wavelengths is more variable (e.g., Vanden Berk et al., 2004)

### 4.2. Dependence of Quasar Characteristic Time Scale and Variability on Mbh

We also looked for trends in quasar variability properties with and the Eddington ratio, . In this work we assume a constant bolometric correction of to the luminosity at 5100Å(Kaspi et al., 2000). However, we stress that a constant bolometric correction can introduce significant error in the Eddington ratio (e.g., Vasudevan & Fabian, 2007; Kelly et al., 2008), and that our use of is only suggestive. Strictly speaking, what is being used in the following regressions is the ratio Å, and in general this will not equal the true Eddington ratio. In Figure 6 we also show as a function of and , and in Figure 7 we also show as a function of and .

The results of the regressions for are

 logτ = (−2.29±1.17)+(0.56±0.14)logMBH  [days] (28) logτ = (2.50±0.24)−(0.06±0.27)logL/Lbol  [days]. (29)

and the results of the regression for are

 logσ2 = (0.33±0.73) (30) − (0.52±0.08)logMBH  [R mag2/day] logσ2 = (−4.33±0.19) (31) − (0.25±0.22)logL/LEdd  [R mag2/day].

Based on these regression results, there is a statistically significant correlation between and , where the relaxation time increases with increasing . In addition, there is significant evidence that decreases with increasing . However, there is no evidence for a dependence of or on the Eddington ratio.

In this work we have employed the linear regression technique described in Kelly (2007). However, this technique assumes symmetric errors in , but in reality the error distribution is asymmetric. In order to assess whether this has any effects on our results, we modifed the technique of Kelly (2007) to handle asymmetric errors in , and used our modifed code to recompute the dependencies involving . However, we did not notice any difference in the slopes by assuming symmetric or asymmetric error bars. This is most likely because of the assumption of a linear relationship between and the other quantities. The asymmetric error bars for the high-mass quasars imply that the variability time scales for these sources are consistent with time scales much longer than the span of the lightcurve, which would lead to a steepening in the slope. However, this steepening of the slope would be inconsistent with the (better determined) time scales for the low-mass AGN, since a steeper line would predict shorter time scales than are observed. So, the assumption of a linear relationship, in combination with the well determined time scales at the low-mass end, counteracts the effects of the asymmetric error bars at the high-mass end. However, if we were to assume a nonlinear form for the dependence of on , such that the slope increases for the high mass objects, then the asymmetric error bars would likely have a noticeable effect, since the low-mass objects would no longer convey a lot of ’information’ about the regression relationship at the high mass end.

Due to the correlation between and , it is unclear whether the observed dependency of and on these quantities is real for both and , or whether one correlation is simply a reflection of the other. Similar to breaking the degeneracy, we can investigate which correlation is the fundamental one, or if both are, by performing a multiple regression of and on and . The results are

 logσ2 = (−3.83±0.17)−(0.09±0.19)log(λLλ1045 erg s−1) (32) − (0.25±0.24)log(MBH108M⊙)   [R mag2/day] τ = (80.4+66.9−35.8)(λLλ1045 erg s−1)−0.42±0.28 (33) × (MBH108M⊙)1.03±0.38   [days]

Here, we have expressed as a function of and instead of for more direct comparison with the characteristic time scales described by Equations (3)–(5). The joint probability distributions of the slopes are shown in Figure 8 for both regressions. It is unclear whether depends on solely , solely , or both and , although the data favor a dependence on over one on . However, there is significant evidence that the relaxation time scale depends on at least , and possibly on as well. In fact, the dependence of on has steepened, and the relationship described by Equation (33) is similar to that for the orbital or thermal time scale, assuming a viscosity parameter of and location of the emitting region to be .

## 5. Discussion

### 5.1. Comparison with Previous Work

Most previous work on quasar optical variability has been based on analysis of structure functions or power spectra, either of individual quasars or an ensemble of quasars. From these studies, an anti-correlation between variability and luminosity has often emerged (e.g., Hook et al., 1994; Garcia et al., 1999; Vanden Berk et al., 2004; Wilhite et al., 2008), while results on a variability–redshift correlation have remained mixed. Our result that long-term quasar variability is uncorrelated with luminosity may appear to be in conflict with previous work. However, the evidence for a variability–luminosity correlation is considerably weaker in the studies that have computed variability measures for individual objects. Indeed, studies that have compared variability with luminosity for individual objects have noted the significant scatter in the relationship, producing a very weak correlation, often leading to a detection of ‘moderate’ statistical significance at best. Ensemble studies, on the other hand, cannot investigate the scatter in the relationship and therefore cannot assess the strength of the correlation. Instead, ensemble studies simply look for an average trend in variability properties, and, given enough quasars, are able to detect even a weak trend of variability with luminosity. Furthermore, quasars exhibit a range in characteristic time scale and variability amplitude at a given luminosity or black hole mass, and it is unclear how this affects an ensemble structure function or power spectrum.

Recently Wold et al. (2007) have reported a correlation between optical variability and on time scales days, but no correlation is seen on shorter time scales. In contrast, we observed an anti-correlation between and short time scale variability, but no correlation between and long time scale variability. The Wold et al. (2007) result is unexpected, since and are correlated (e.g., Peterson et al., 2004), and previous studies have found that variability is anti-correlated with , even on long time scales. The Wold et al. (2007) result is based on an ensemble structure function, and the uncertainties on the structure function for the high bins are large. Furthermore, the time sampling of the Wold et al. (2007) sample for the high bins is worse than for the low bins, and windowing effects due to the finite length of the time series may be at work here. When considering variability measurement of individual sources, Wold et al. (2007) still find a positive correlation between variability and . However, while this correlation is statistically significant, it is very weak and exhibits considerable scatter. We performed a Kendall and Spearman rank correlation test between long time scale variability and the estimated black hole mass, and also find a marginally significant correlation, but the significance disappeared when we accounted for the uncertainty in the mass estimates. In addition, because the long time scale variability for a CAR(1) process is , Equations (28) and (30) imply that the long time scale variability depends only weakly on , if at all.

Collier & Peterson (2001) calculated structure functions of optical light curves for 12 low- AGN with reverberation mapping data. Consistent with our work, they find a correlation between and characteristic time scale, where a characteristic time scale was defined as the location of a break in the structure function. The time scales found by Collier & Peterson (2001) were consistent with dynamical or disk thermal time scales, although they tended to be somewhat shorter than those derived in this work. This difference may be explained by somewhat different definitions of ‘characteristic’ time scale between their work and ours, and systematic errors in the estimated structure functions caused by time sampling effects.

Two of the sources in our sample, NGC 5548 and NGC 4151, were studied by Czerny et al. (1999) and Czerny et al. (2003), respectively, using spectral techniques. These authors fit a functional form to the optical power spectra similar to that implied by the CAR(1) model, with the exception that they allow the high-frequency slope to vary. In contrast, the CAR(1) model we employ fixes this slope to be , i.e., . For NGC 5548, Czerny et al. (1999) find evidence for a flattening of the optical power spectrum on time scales days. This is longer than our observed characteristic time scale of at significance. This difference is only marginally significant, and may be explained by the slightly different parameteric forms assumed for the power spectra, the longer light curve used by Czerny et al. (1999), the different photometric bands, and errors introduced by the window function and smoothing of the power spectra in the Czerny et al. (1999) analysis.

For NGC 4151, Czerny et al. (2003) did not find any evidence for a flattening to white noise of the -band optical power spectra over the time scales probed by their 90 yr light curve. However, they find some evidence that the power spectra flattens below time scales of days, in agreement with the characteristic time scale we find of days. Furthermore, we note that our analysis constrains the characteristic time scale for NGC 4151 to be yr at confidence. If the power spectra does flatten on time scales yr, it would be difficult to detect in the spectral analysis employed by Czerny et al. (2003), as such a turn-over would occur in the few lowest frequency bins. This is especially true when one considers that the lowest frequency bins of the power spectra are also among the bins most biased by windowing effects and smoothing.

### 5.2. Connection with Accretion Physics

In this work we have found that both the characteristic time scale of quasar light curves, and the magnitude of the short time scale variations, depend on black hole mass. This, in combination with the evidence from reverberation mapping, strongly argues that the source of quasar optical variability is intrinsic to the accretion disk. In addition, we have found that the characteristic time scales of quasar light curves are similar to what would be expected for disk dynamical or thermal time scales, assuming a viscosity parameter of . Recent MHD simulations of radiation-dominated AGN accretion disks have found that the thermal time scale is shorter than that implied by the standard -prescription (Turner, 2004), making the association of with thermal time scales more consistent. Our results imply that on time scales shorter than or , the accretion disk has difficulty generating variations in optical flux in response to random variations of some input process, such as, for example, a time varying magnetic field. Instead, these short time scale variations get ‘smoothed out’, creating a power spectrum for frequencies higher than or .

While the quasar characteristic time scales are consistent with both orbital and thermal time scales, we find it more appropriate to associate these time scales with , as we might expect some sort of periodic activity in the light curves if the flux variations were driven by orbital motion. In addition, quasar optical emission is thought to be thermal emission from an optically thick accretion disk (e.g., Krolik, 1999; Frank, King, & Raine, 2002), and quasars tend to be bluer as they brighten (e.g., Giveon et al., 1999; Trèvese et al., 2001), suggesting that the flux variations are due to thermal variations. An association with orbital time scales cannot be ruled out, and indeed, the characteristic time scale of the turbulence driving the magnetic energy and kinetic energy density in the disk is (e.g., Hirose et al., 2008). However, if the radiation energy density in the disk is due to dissipation of magnetic and kinetic energy density, then it is expected to respond to fluctuations in the magnetic and kinetic energy densities with a characteristic time scale of order the thermal time.

In order to interpret the CAR(1) process in terms of accretion disk physics, we rewrite Equation (1) as

 dlogL(t)=−1τ(logL(t)−μ)dt+σ[B(t+dt)−B(t)]. (34)

Here, denotes the luminosity of the quasar at time , is the mean value of the quasar light curve, and denotes Brownian motion. In Equation (34) we have used the fact that the derivative of Brownian motion is white noise, i.e., . Brownian motion is a non-stationary random walk process that has a power spectrum , and is described by Equation (1) in the limit . In addition, Equation (34) implicitly assumes that the variance in the random variable is . The notation should not be confused with the common physics notation of denoting the amplitude of a magnetic field; although in this work we suggest that may be associated with a turbulent magnetic field, the term should be understood as referring to a Brownian motion.

Writing the Equation for a CAR(1) process as Equation (34) reveals a number of interesting properties of this process. First, we note that the first term on the right side is what keeps the time series stationary. Considering only this term (i.e., ), is negative when the value of is brighter than the mean, and is positive when is fainter than the mean. Therefore, the first term on the right side stabilizes the process by always driving toward its mean value, while the second term generates random perturbations to that cause to deviate from its expected path. For highly accreting objects like quasars, the accretion disks are expected to be radiation pressure dominated in the regions emitting the optical and UV flux. Under the standard -prescription for the viscosity, where the viscous torque is assumed to be proportional to the total pressure, a radiation pressure dominated disk is unstable to perturbations in the heating rate (e.g., Shakura & Sunyaev, 1976; Krolik, 1999). The fact that the quasar light curves in our sample are described well by a CAR(1) process with relaxation times similar to disk thermal time scales rules out instabilities in the disk that grow as , consistent with results obtained from 3-dimensional MHD simulations (e.g., Turner, 2004; Hirose et al., 2008).

From Equation (34) it is apparent that the stochastic input into the differential equation, which drives the random variations in , is itself a stochastic process. In particular, a random deviation in the input process, , over a time interval causes a random perturbation to the change in expected over the interval , with controlling how sensitive is to . For example, the input process could be due to variations in accretion rate, or perturbations caused by a time-varying magnetic field. Indeed, recently some authors have modeled quasar X-ray variability as being driven by variations in a magnetic field, with the magnetic field density being modelled as an AR(1) process (e.g., King et al., 2004; Mayer & Pringle, 2006; Janiuk & Czerny, 2007).

Hirose et al. (2008) have argued for a similar interpretation of quasar optical variability, based on 3-dimensional MHD simulations. They argue that fluctuations in the magnetic field are dissipated, transferring energy from the magnetic field to heat in the plasma, creating thermal fluctuations in the accretion disk. The fluctuations in the heat content of the disk then create fluctuations in ; however, because the disk radiation cannot react to changes in heat content on time scales less than the thermal time scale, the shorter time scale fluctuations in flux are smoothed out and damped. Short time scale variations in radiation are correlated because the radiation energy density of the disk has not had time to completely react to the change in heat content induced by the change in magnetic energy density. However, on time scales , the disk has had time to adjust to the heat content variations, thereby ‘forgetting’ about the previous perturbations in heat content. The result is a red noise power spectrum on time scales , and a white noise power spectrum on time scales .

Within this interpretation, the parameter describes the amplitude of variations in flux caused by variations in the magnetic energy density. Therefore, the anti-correlation between and may imply two things. First, it may imply that the magnitude of the stochastic input process, , possibly associated with a turbulent magnetic field, decreases with increasing . Or, it may imply that the sensitivity of changes in radiation energy density to changes in the turbulent magnetic field depends on , where fluctuations in the magnetic field create smaller fluctuations in as increases. Alternatively, it could imply that both effects are at work.

Modelling quasar variability as a stochastic process provides an opportunity to unify both short and long time scale variability as the result of a single process. The source of variations in radio-quiet quasar optical luminosity over time scales of hours (so-called ‘microvariability’ or ‘intranight variability’) has remained a puzzle, although reprocessing of X-rays or a weak blazar component have been suggested (Czerny et al., 2008). In general, microvariability in radio-quiet quasars is not detected above the photometric uncertainty (e.g., Gupta & Joshi, 2005; Carini et al., 2007); however, for those sources for which it is detected the standard deviation in the variability over the course of a night is mag (Gopal-Krishna et al., 2003; Stalin et al., 2004, 2005; Gupta & Joshi, 2005). As noted in the Appendix, for a CAR(1) process the standard deviation of short time scale variations is . As can be seen from figure 5, our best-fitting CAR(1) processes for quasar light curves predict variations of mag over 8 hours, consistent with what has been observed for radio-quiet quasars. Furthermore, it implies that the amplitude of microvariability should decrease with increasing black hole mass.

To further investigate whether our best-fitting CAR(1) processes are consistent with observed microvariability of radio-quiet quasars, we calculate the expected distribution of microvariability amplitude for three different values of black hole mass: and . We simulate light curves using the values of and calculated from Equations (28) and (30) for each value of . Each simulated light curve had 40 data points, regularly sampled over the course of 5 hrs; these values were chosen to be consistent with recent observations of microvariability (e.g., Stalin et al., 2005; Gupta & Joshi, 2005). For simplicity we neglect observational error. Following Romero et al. (1999) we calculate the amplitude of microvariability from the difference in the maximum and minimum observed flux values, expressed as a per cent:

 ψ=100¯¯¯¯¯¯¯¯¯¯L(t)(Lmax(t)−Lmin(t)). (35)

Here, is the average value of the light curve. The results are shown in Figure 9. As can be seen, the amplitudes of microvariability expected from our best fit CAR(1) processes are consistent with what is observed, at least for radio-quiet quasars (e.g., Gopal-Krishna et al., 2003; Stalin et al., 2004, 2005; Gupta & Joshi, 2005).

Because the magnitude of microvariability predicted by our best-fit CAR(1) process is consistent with observations of microvariability in radio-quiet quasars, it is not necessary to invoke an additional physical mechanism to explain the short time scale variations in these objects. Instead, the CAR(1) process is able to explain both short and long time scale variations as being driven by the same process, possibly thermal fluctuations being driven by a turbulent magnetic field. Czerny et al. (2008) investigated three different models for microvariability in radio-quiet quasars and concluded that the most promising source of the microvariability was a weak blazar-like jet. Furthermore, Czerny et al. (2008) concluded that luminosity fluctuations driven by the magneto-rotational instability are not strong enough to create microvariability. Our empirical results are inconsistent with the conclusions of Czerny et al. (2008) in that we do not find any evidence for an additional process that drives the microvariability. This is not to say that our results imply that a weak jet does not exist in radio-quiet quasars, but rather that if a weak jet does exist then it’s variability amplitude is small compared to that caused by the process that also drives the longer time scale variations, such as a turbulent magnetic field.

It is unclear how sensitive the conclusions of Czerny et al. (2008) are to the specifics of their model, and the investigation of Czerny et al. (2008) does not necessarily rule a turbulent magnetic field as being the source of microvariability in radio-quiet quasars. Of particular interest is the manner in which the time dependence of the magnetic field is modeled. Czerny et al. (2008), as well as King et al. (2004), Mayer & Pringle (2006), and Janiuk & Czerny (2007), have modeled the time series of the magnetic field as an AR(1) process:

 ui+1=αARui+ϵi, (36)

where the amplitude of the local magnetic field at the time step is proportional to , and is a random variable uniformly distributed between . Because the characteristic time scale of the magnetic field development is , (Czerny et al., 2008) used a time step equal to the orbital time. However, Equation (4), in combination with Figure 10 of (Czerny et al., 2008), implies that the orbital time scale should be days, depending on the mass of the black hole. This results in a time sampling of the magnetic field that is much longer than the time scales over which microvariability are measured, therefore artificially suppressing the amplitude of microvariability in the simulated light curves. Moreover, the amplitude of microvariability can also be increased by increasing the variance of the random variable .

Following King et al. (2004), Czerny et al. (2008) assumed a value of . As described in the appendix, a value of corresponds to a characteristic time scale for the equivalent continuous process of , i.e., a value of . A more accurate simulated light curve can be obtained by using a sampling interval in Equation (36) that is small compared to that probed by microvariability, such as minutes. Then, the value of can be related to the characteristic time scale of the magnetic field development as . In this case, the only free parameters are and the variance of . Further improvement can be obtained by employing 3-dimensional MHD simulations to test whether a turbulent magnetic field can create microvariability with fluctuations .

## 6. Summary

In this work we have modeled quasar light curves as a type of stochastic process called a first-order continuous autoregressive process. This statistical model has three free parameters: the characteristic time scale for the process to ‘forget’ about itself, , the magnitude of the small time scale variations, , and the mean of the time series, . We used this model to fit 100 quasar light curves at , including 70 quasars with black hole mass estimates. Our conclusions are summarized as follows:

• Quasar optical light curves are often well described by a continuous autoregressive process (CAR(1)). For the quasars in our sample, the amplitude of variations on time scales long compare to their characteristic time scales (i.e., decades) are typically , consistent with previous work. The characteristic time scales of the quasar light curves vary between days and yrs and are consistent with accretion disk orbital or thermal time scales, assuming a viscosity parameter of and location of the emitting region to be . In addition, the short time scale variations are .

• The characteristic time scales of quasar optical light curves are correlated with and luminosity, while the magnitude of the short time scale variations are anti-correlated with and luminosity. We did not find any evidence for an additional redshift correlation. A multiple regression analysis suggested that the primary correlation is with . At a given luminosity, the characteristic time scales depend on as

 τ = (80.4+66.9−35.8)(λLλ1045 erg s−1)−0.42±0.28 × (MBH108M⊙)1.03±0.38   [days]

where the errors are quoted at confidence ().

• Within the CAR(1) process model of quasar light curves, the random perturbations to are caused by a stochastic process. This stochastic input process may be a turbulent magnetic field that creates fluctuations in the accretion disk radiation energy density. Variations in the input process over an interval create variations in which are smoothed out on time scales shorter than an orbital or thermal time scale.

• The fact that radio-quiet quasar optical light curves can be well fit by a CAR(1) model suggests that it is not necessary to invoke an additional physical mechanism to describe short time scale variations in these objects. Instead, within the CAR(1) model, both short and long time scale variations are driven by an underlying stochastic process that causes random perturbations to . This is supported by the fact that the intranight variations predicted by our best fit CAR(1) processes are in the -band, consistent with observations of intranight variability. Furthermore, our best-fit CAR(1) parameters predict an anti-correlation between radio-loudness and microvariability, inconsistent with what has been observed. This suggests that radio-loud quasars have an additional optical variability mechanism that operates on time scales shorter than those probed by our study, days.

Before concluding, we stress that the CAR(1) model is a statistical model and not a physical model. Quasar light curves are stochastic in nature (e.g., Czerny & Lehto, 1997), and the dependence of luminosity on the time-varying properties of the disk is complex. In this sense, the randomness in the stochastic model is not due to that fact that the physical processes themselves are not deterministic, but rather is a reflection of our lack of knowledge of the complex physical processes that generate variations in flux. While a physical model is needed in order to interpret the stochastic model in terms of accretion disk physics, and thus lead to a proper understanding of quasar light curves, the stochastic model is sufficient for modeling the data, given our current knowledge. Furthermore, much of the mathematical formalism of accretion physics is in the language of differential equations, suggesting that stochastic differential equations are a natural choice for modeling quasar light curves.

The field of stochastic processes is a rich field with well-developed methodology, predominantly because of its importance in financial and economic modeling. We have utilized the CAR(1) model because of its simplicity, and because it allows us to perform statistical inference without having our results biased by the irregular sampling, measurement errors, and finite span of the time series. However, the CAR(1) model is the simplest of stationary continuous autoregressive processes, and additional flexibility may be achieved through the addition of higher order derivatives to Equation (1). This provides a rich and flexible method of modeling the power spectra of quasar light curves without suffering from the windowing effects that can bias traditional Fourier and structure function techniques. For example, quasi-periodic oscillations can be modeled through the addition of second order derivatives to Equation (1), as has been done in the analysis of the frequency of sun spot numbers (Phadke & Wu, 1974). In addition, Equation (1) can be generalized to a vector form, allowing the simultaneous modeling of quasar light curves across multiple observing wavelengths, and thus introducing additional constraints on physical models of quasar variability. Both the use of higher order terms and multiwavelength modelling of quasar light curves will be the subject of future research.

A computer routine for fitting the CAR(1) model to a time series is available on request from B. Kelly.

We would like to thank Tim Axelrod for supplying the MACHO light curves, and Marla Geha for providing optical spectra for a number of the MACHO sources and providing helpful comments on an early draft of this paper. We would also like to thank the referee, Bożena Czerny, for a careful reading of the manuscript and comments that led to its improvement. BK acknowledges support from NASA through Hubble Fellowship grant #HF-01220.01 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. AS acknowledges support from NASA contract NAS 8-39073.

## Appendix A Description of Autoregressive Processes

The first order autoregressive process (AR(1)) process is a well-studied stochastic process (e.g., see Scargle, 1981, and references therein), generated according to

 xi=αARxi−1+ϵi, (A1)

where is a normally-distributed random variable with zero mean and variance , and the data are observed at regular time intervals. The parameters of the AR(1) model are and , and is usually constrained as in order to ensure stationarity; a time series is said to be stationary when its mean and covariance do not vary with time222Technically, this is referred to as ’weakly stationary’, while stationarity in the strictest sense implies that the probability distribution of the time series does not change with time. However, because we are only concerned with Gaussian processes in this work, the two definitions are equivalent. The case = 1 corresponds to a random walk. For quasar light curves, the observed data points correspond to the observed fluxes at times , for .

The discrete AR(1) process is only defined for regularly sampled time series. However, astronomical time series are rarely regularly sampled, and often large gaps in time can exist. Furthermore, the AR(1) process is a discrete process, but the underlying physical process that gives rise to the observed flux is continuous. Because of these two considerations, we instead model the quasar light curves as a first order continuous autoregressive (CAR(1)) process. As noted in § 3, the CAR(1) process is described by the following stochastic differential equation (e.g., Brockwell & Davis, 2002):

 dX(t)=−1τX(t)dt+σ√dtϵ(t)+b dt,   τ,σ,t>0, (A2)

where, is called the ‘relaxation time’ of the process , and is a white noise process with zero mean and variance equal to one. Throughout this work we have assumed that the white noise process is also Gaussian. In the physics literature, Equation (A2) is often referred to as an Ornstein-Uhlenbeck (O-U) process, and plays a central role in the mathemematics of Brownian motion; see Gillespie (1995) for a review of the O-U process.

The solution to Equation (A2) is

 X(t)=e−t/τX(0)+bτ(1−e−t/τ)+σ∫t0e−(t−s)/τdB(s). (A3)

Here, is a random variable describing the initial value of the time series, is a temporally uncorrelated normally distributed random variable with zero mean and variance . Strictly speaking, is an interval of Brownian motion, and the integral on the right side represents the stochastic component of the time series. From Equation (A3) it can be seen that if , i.e., if there is no stochastic component, and if represents a random perturbation, relaxes to it’s mean value with an -folding time scale ; hence the identification of as the relaxation time. If , then the path that takes will vary randomly about the expected exponential relaxation.

The expected value of given for is

 E(X(t)|X(s))=e−Δt/τX(s)+bτ(1−e−Δt/τ) (A4)

and the variance in given is

 Var(X(t)|X(s))=τσ22[1−e−2Δt/τ] (A5)

where . If , then Equation (A5) implies that the variance on time scales much shorter than the relaxation time is . Therefore, can be interpreted as representing the variance in the light curve on short time scales. In addition, one can show that when the time sampling is regular with , then the CAR(1) process reduces to an AR(1) process with and .

In astronomical time series analysis it is common to interpret a light curve in terms of its autocorrelation function and power spectrum. The autocovariance function at time is defined to be the expected value of the product of and , and the autocorrelation function is calculated by dividing the autocovariance function by the variance of the time series. The autocorrelation function of the CAR(1) process is

 ACF(t′)=e−t′/τ. (A6)

Equation (A6) states that the correlations in CAR(1) light curve fall off exponentially with lag , with an -folding time equal to the relaxation time, . Following Gillespie (1995), the power spectrum of a process is computed from the autocovariance function as

 PX(f)=4∫∞0rX(t′)cos(2πft′) dt′,   f≥0. (A7)

For a CAR(1) process, , from which it follows that the power spectrum of a CAR(1) process is

 PX(f)=2σ2τ21+(2πτf)2. (A8)

The power spectrum of the CAR(1) process falls off as on time scales short compared to the relaxation time, and flattens to white noise at time scales long compared to the relaxation time.

## References

• Alcock et al. (1997) Alcock, C., et al. 1997, ApJ, 486, 697
• Alcock et al. (1999) Alcock, C., et al. 1999, PASP, 111, 1539
• Aretxaga et al. (1997) Aretxaga, I., Cid Fernandes, R., & Terlevich, R. J. 1997, MNRAS, 286, 271
• Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
• Balbus & Hawley (1998) Balbus, S. A., & Hawley, J. F. 1998, Reviews of Modern Physics, 70, 1
• Brockwell & Davis (2002) Brockwell, P. .J., & Davis, R. A. 2002 Introduction to Time Series and Forecasting (2nd Ed.; New York, NY: Springer)
• Carini et al. (2007) Carini, M. T., Noble, J. C., Taylor, R., & Culler, R. 2007, AJ, 133, 303
• Carone et al. (1996) Carone, T. E., et al. 1996, ApJ, 471, 737
• Cid Fernandes et al. (1996) Cid Fernandes, R. J., Aretxaga, I., & Terlevich, R. 1996, MNRAS, 282, 1191
• Cid Fernandes et al. (2000) Cid Fernandes, R., Sodré, L., Jr., & Vieira da Silva, L., Jr. 2000, ApJ, 544, 123
• Collier et al. (1998) Collier, S. J., et al. 1998, ApJ, 500, 162
• Collier & Peterson (2001) Collier, S., & Peterson, B. M. 2001, ApJ, 555, 775
• Cristiani et al. (1990) Cristiani, S., Vio, R., & Andreani, P. 1990, AJ, 100, 56
• Cristiani et al. (1996) Cristiani, S., Trentini, S., La Franca, F., Aretxaga, I., Andreani, P., Vio, R., & Gemmo, A. 1996, A&A, 306, 395
• Cutri et al. (1985) Cutri, R. M., Wisniewski, W. Z., Rieke, G. H., & Lebofsky, M. J. 1985, ApJ, 296, 423
• Czerny et al. (2003) Czerny, B., Doroshenko, V. T., Nikołajuk, M., Schwarzenberg-Czerny, A., Loska, Z., & Madejski, G. 2003, MNRAS, 342, 1222
• Czerny & Lehto (1997) Czerny, B., & Lehto, H. J. 1997, MNRAS, 285, 365
• Czerny et al. (1999) Czerny, B., Schwarzenberg-Czerny, A., & Loska, Z. 1999, MNRAS, 303, 148
• Czerny et al. (2008) Czerny, B., Siemiginowska, A., Janiuk, A., & Gupta, A. C. 2008, MNRAS, 386, 1557
• de Vries et al. (2005) de Vries, W. H., Becker, R. H., White, R. L., & Loomis, C. 2005, AJ, 129, 615
• di Clemente et al. (1996) di Clemente, A., Giallongo, E., Natali, G., Trevese, D., & Vagnetti, F. 1996, ApJ, 463, 466
• Frank, King, & Raine (2002) Frank, J., King, A., & Raine, D. 2002, Accretion Power in Astrophysics (3rd Edition;Cambridge, UK:Cambridge Univ. Press)
• Garcia et al. (1999) Garcia, A., Sodré, L., Jablonski, F. J., & Terlevich, R. J. 1999, MNRAS, 309, 803
• Geha et al. (2003) Geha, M., et al. 2003, AJ, 125, 1
• Gillespie (1995) Gillespie, D. T. 1996, Am. J. Phys., 64, 225
• Giveon et al. (1999) Giveon, U., Maoz, D., Kaspi, S., Netzer, H., & Smith, P. S. 1999, MNRAS, 306, 637
• Gopal-Krishna et al. (2003) Gopal-Krishna, Stalin, C. S., Sagar, R., & Wiita, P. J. 2003, ApJ, 586, L25
• Gupta & Joshi (2005) Gupta, A. C., & Joshi, U. C. 2005, A&A, 440, 855
• Hawkins (2000) Hawkins, M. R. S. 2000, A&AS, 143, 465
• Hawkins (2007) Hawkins, M. R. S. 2007, A&A, 462, 581
• Helfand et al. (2001) Helfand, D. J., Stone, R. P. S., Willman, B., White, R. L., Becker, R. H., Price, T., Gregg, M. D., & McMahon, R. G. 2001, AJ, 121, 1872
• Hirose et al. (2008) Hirose, S., Krolik, J. H., & Blaes, O. 2008, in press at ApJ(arXiv:0809.1708)
• Hook et al. (1994) Hook, I. M., McMahon, R. G., Boyle, B. J., & Irwin, M. J. 1994, MNRAS, 268, 305
• Janiuk et al. (2000) Janiuk, A., Czerny, B., & Siemiginowska, A. 2000, ApJ, 542, L33
• Janiuk et al. (2002) Janiuk, A., Czerny, B., & Siemiginowska, A. 2002, ApJ, 576, 908
• Janiuk & Czerny (2007) Janiuk, A., & Czerny, B. 2007, A&A, 466, 793
• Kaspi et al. (1996) Kaspi, S., et al. 1996, ApJ, 470, 336
• Kaspi et al. (2000) Kaspi, S., Smith, P. S., Netzer, H., Maoz, D., Jannuzi, B. T., & Giveon, U. 2000, ApJ, 533, 631
• Kawaguchi et al. (1998) Kawaguchi, T., Mineshige, S., Umemura, M., & Turner, E. L. 1998, ApJ, 504, 671
• Kellermann et al. (1989) Kellermann, K. I., Sramek, R., Schmidt, M., Shaffer, D. B., & Green, R. 1989, AJ, 98, 1195
• Kelly (2007) Kelly, B. C. 2007, ApJ, 665, 1489
• Kelly & Bechtold (2007) Kelly, B. C., & Bechtold, J. 2007, ApJS, 168, 1
• Kelly et al. (2008) Kelly, B. C., Bechtold, J., Trump, J. R., Vestergaard, M., & Siemiginowska, A. 2008, ApJS, 176, 355
• King et al. (2004) King, A. R., Pringle, J. E., West, R. G., & Livio, M. 2004, MNRAS, 348, 111
• Krolik (1999) Krolik, J. H. 1999, Active Galactic Nuclei: From the Central Engine to the Galactic Environment (Princeton, NJ:Princeton Univ. Press)
• Lightman & Eardley (1974) Lightman, A. P., & Eardley, D. M. 1974, ApJ, 187, L1
• Lub & de Ruiter (1992) Lub, J., & de Ruiter, H. R. 1992, A&A, 256, 33
• Lyubarskii (1997) Lyubarskii, Y. E. 1997, MNRAS, 292, 679
• Manmoto et al. (1996) Manmoto, T., Takeuchi, M., Mineshige, S., Matsumoto, R., & Negoro, H. 1996, ApJ, 464, L135
• Mayer & Pringle (2006) Mayer, M., & Pringle, J. E. 2006, MNRAS, 368, 379
• Merloni (2003) Merloni, A. 2003, MNRAS, 341, 1051
• Merloni & Fabian (2002) Merloni, A., & Fabian, A. C. 2002, MNRAS, 332, 165
• Nayakshin et al. (2000) Nayakshin, S., Rappaport, S., & Melia, F. 2000, ApJ, 535, 798
• Pessah et al. (2008) Pessah, M. E., Chan, C.-K., & Psaltis, D. 2008, MNRAS, 383, 683
• Peterson et al. (2000) Peterson, B. M., et al. 2000, ApJ, 542, 161
• Peterson et al. (2002) Peterson, B. M., et al. 2002, ApJ, 581, 197
• Peterson et al. (2004) Peterson, B. M., et al. 2004, ApJ, 613, 682
• Phadke & Wu (1974) Phadke, M. S., & Wu, S. M. 1974, J. Amer. Statist. Assoc., 69, 325
• Richards et al. (2001) Richards, G. T., et al. 2001, AJ, 121, 2308
• Romero et al. (1999) Romero, G. E., Cellone, S. A., & Combi, J. A. 1999, A&AS, 135, 477
• Santos-Lleó et al. (2001) Santos-Lleó, M., et al. 2001, A&A, 369, 57
• Scargle (1981) Scargle, J. D. 1981, ApJS, 45, 1
• Shakura & Syunyaev (1973) Shakura, N. I., & Syunyaev, R. A. 1973, A&A, 24, 337
• Shakura & Sunyaev (1976) Shakura, N. I., & Sunyaev, R. A. 1976, MNRAS, 175, 613
• Shemmer et al. (2001) Shemmer, O., et al. 2001, ApJ, 561, 162
• Siemiginowska & Czerny (1989) Siemiginowska, A., & Czerny, B. 1989, MNRAS, 239, 289
• Spergel et al. (2003) Spergel, D. N., et al. 2003, ApJS, 148, 175
• Stalin et al. (2004) Stalin, C. S., Gopal-Krishna, Sagar, R., & Wiita, P. J. 2004, MNRAS, 350, 175
• Stalin et al. (2005) Stalin, C. S., Gupta, A. C., Gopal-Krishna, Wiita, P. J., & Sagar, R. 2005, MNRAS, 356, 607
• Starling et al. (2004) Starling, R. L. C., Siemiginowska, A., Uttley, P., & Soria, R. 2004, MNRAS, 347, 67
• Stella & Rosner (1984) Stella, L., & Rosner, R. 1984, ApJ, 277, 312
• Stirpe et al. (1994) Stirpe, G. M., et al. 1994, ApJ, 425, 609
• Svensson & Zdziarski (1994) Svensson, R., & Zdziarski, A. A. 1994, ApJ, 436, 599
• Szuszkiewicz (1990) Szuszkiewicz, E. 1990, MNRAS, 244, 377
• Trèvese et al. (2001) Trèvese, D., Kron, R. G., & Bunone, A. 2001, ApJ, 551, 103
• Trèvese & Vagnetti (2002) Trèvese, D., & Vagnetti, F. 2002, ApJ, 564, 624
• Turner (2004) Turner, N. J. 2004, ApJ, 605, L45
• Ulrich et al. (1997) Ulrich, M.-H., Maraschi, L., & Urry, C. M. 1997, ARA&A, 35, 445
• Uttley et al. (2005) Uttley, P., McHardy, I. M., & Vaughan, S. 2005, MNRAS, 359, 345
• van der Klis (1997) van der Klis, M. 1997, Statistical Challenges in Modern Astronomy II, 321
• Vanden Berk et al. (2004) Vanden Berk, D. E., et al. 2004, ApJ, 601, 692
• Vasudevan & Fabian (2007) Vasudevan, R. V., & Fabian, A. C. 2007, MNRAS, 381, 1235
• Vestergaard & Peterson (2006) Vestergaard, M., & Peterson, B. M. 2006, ApJ, 641, 689
• Wilhite et al. (2008) Wilhite, B. C., Brunner, R. J., Grier, C. J., Schneider, D. P., & vanden Berk, D. E. 2008, MNRAS, 383, 1232
• Wold et al. (2007) Wold, M., Brotherton, M. S., & Shang, Z. 2007, MNRAS, 375, 989
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters