Modeling the Time Variability of SDSS Stripe 82 Quasars as a Damped Random Walk
Abstract
We model the time variability of 9,000 spectroscopically confirmed quasars in SDSS Stripe 82 as a damped random walk. Using 2.7 million photometric measurements collected over 10 years, we confirm the results of Kelly et al. (2009) and Kozłowski et al. (2010) that this model can explain quasar light curves at an impressive fidelity level (0.010.02 mag). The damped random walk model provides a simple, fast [O() for data points], and powerful statistical description of quasar light curves by a characteristic time scale () and an asymptotic rms variability on long time scales (). We searched for correlations between these two variability parameters and physical parameters such as luminosity and black hole mass, and restframe wavelength. Our analysis shows to increase with decreasing luminosity and restframe wavelength as observed previously, and without a correlation with redshift. We find a correlation between and black hole mass with a power law index of 0.180.03, independent of the anticorrelation with luminosity. We find that increases with increasing wavelength with a power law index of 0.17, remains nearly constant with redshift and luminosity, and increases with increasing black hole mass with power law index of 0.210.07. The amplitude of variability is anticorrelated with the Eddington ratio, which suggests a scenario where optical fluctuations are tied to variations in the accretion rate. However, we find an additional dependence on luminosity and/or black hole mass that cannot be explained by the trend with Eddington ratio. The radioloudest quasars have systematically larger variability amplitudes by about 30%, when corrected for the other observed trends, while the distribution of their characteristic time scale is indistinguishable from that of the full sample. We do not detect any statistically robust differences in the characteristic time scale and variability amplitude between the full sample and the small subsample of quasars detected by ROSAT. Our results provide a simple quantitative framework for generating mock quasar light curves, such as currently used in LSST image simulations.
1 Introduction
The optical variability of quasars has been recognized since they were first identified (Matthews & Sandage 1963). Indeed, most quasars are variable (90% at the 0.03 mag rms level; Sesar et al. 2007), and the variations in brightness are aperiodic and on the order of 20% on time scales of months to years (e.g., Hook et al. 1994; Vanden Berk et al. 2004). Furthermore, the smooth power spectra suggest a chaotic, or stochastic, origin for the variability. A range of models have been advanced to describe quasar variability, including supernova bursts, microlensing, and accretion disk instabilities (Aretxaga et al. 1997; Hawkins 1993; Kawaguchi et al. 1998; Trèvese & Vagnetti 2002). These models are discussed and compared in Hawkins (2007). Reverberation mapping studies (e.g., Peterson et al. 2005) show that the broad emission lines respond to continuum fluctuations, therefore providing strong evidence that the variability is intrinsic to the quasars. A number of studies have utilized standard accretion disk models to demonstrate that the opticalUV variability of quasars could be driven by a variable accretion rate (e.g., Pereyra et al. 2006; Li & Cao 2008; Liu et al. 2008). Blackburne & Kochanek (2010) find evidence in the light curves of microlensed quasars that the optical variability is caused by a change in the effective area of the accretion disk.
Recently, Kelly et al. (2009, hereafter KBS09) proposed a model where the optical variability is described by a damped random walk (a selfcorrecting term added to a random walk model that acts to push any deviations back towards the mean value). They proposed that the variability time scale might be identified with the thermal time scale of accretion disks, as also proposed by Collier & Peterson (2001). A thermal origin of the variability would explain why quasars become bluer as they brighten (e.g., Giveon et al. 1999; Trèvese et al. 2001; Geha et al. 2003).
Although the physical causes have yet to be proven, it has been established by KBS09 and Kozłowski et al. (2010a, hereafter Koz10) that a damped random walk can statistically explain the observed light curves of quasars. Using 100 wellsampled singleband light curves compiled from the literature, KBS09 show that this stochastic process is capable of modeling complex quasar light curves at an impressive fidelity level (0.010.02 mag). Koz10 applied the model to the OGLE light curves (Udalski et al. 1997; Udalski et al. 2008) of midinfraredselected quasars behind the Magellanic Clouds from Kozłowski & Kochanek (2009). Their analysis shows that this stochastic model is robust enough to efficiently select quasars from other variable sources (see Schmidt et al. 2010 for a different method of selecting quasars based on variability). The model has only three free parameters: the mean value of the light curve, the driving amplitude of the stochastic process, and the damping time scale. The predictions are only statistical, and the random nature reflects our uncertainty about the details of the physical processes.
Instead of applying a model to observed light curves for individual quasars, numerous studies have looked at the ensemble variability of quasars, particularly in samples where individual light curves are not available. Significant progress in the description of quasar variability has been made by employing the Sloan Digital Sky Survey (SDSS) data (Vanden Berk et al. 2004, hereafter VB04; Ivezić et al. 2004, hereafter I04; de Vries et al. 2005; Wilhite et al. 2005, 2006, 2008; Sesar et al. 2006). For example, the size and quality of the sample analyzed by VB04 (twoepoch photometry for 25,000 spectroscopically confirmed quasars) allowed them to constrain how quasar variability in the rest frame optical/UV regime depends upon restframe time lag (up to 2 years), luminosity, rest wavelength, redshift, the presence of radio and Xray emission, and the presence of broad absorption line systems. Using repeated SDSS photometric observations, Wilhite et al. (2008) confirmed the result of Wold et al. (2007) that variability is correlated with black hole mass, and show that this is independent of the anticorrelation between variability and luminosity established by many studies. This led them to suggest that the amplitude of variability may be driven by the quasar’s Eddington ratio, implying differences in accretion rate.
These studies typically quantify the observed optical variability of quasars using a structure function (SF) analysis (see also Hughes et al. 1992; Collier & Peterson 2001; Bauer et al. 2009; Kozłowski et al. 2010b), where the SF is the rootmeansquare (rms) magnitude difference as a function of the time lag () between measurements. This autocorrelationlike function is less sensitive to aliasing and other timesampling problems than a power spectral distribution. By studying the magnitude difference distribution for appropriately chosen subsamples with fixed values of absolute iband magnitude (), restframe time lag (, in days), and wavelength (, in Å), the mean dependence of the SF on these quantities was inferred by I04 to be
(1) 
with , , and . A qualitatively similar result was obtained by VB04. Kozłowski et al. (2010b), in the first large study of the midIR structure functions of quasars, also found lower variability for higher luminosities and longer wavelengths, but the temporal slope of the ensemble structure functions were significantly steeper than in the optical. In addition, there is evidence for a turnover in the SF on long time lags (I04; Rengstorf et al. 2006; Wold et al. 2007). Studies by de Vries et al. (2005) and Sesar et al. (2006) using SDSS combined with earlier Palomar Observatory Sky Survey measurements for 40,000 SDSS quasars constrained quasar continuum variability on time scales of 10 to 50 years in the observer’s frame. They report that the characteristic time scale, which in this context is the time lag above which the SF flattens to a constant value, is of order 1 year in the quasar rest frame. Using a shotnoise light curve model, de Vries et al. (2005) found evidence for multiple variability time scales in longterm ensemble variability measurements, while Collier & Peterson (2001) found a wide range of different time scales in their analysis of individual light curves, and even evidence for multiple time scales in a single active galactic nucleus (AGN).
These analyses of ensemble variability are based on a fundamental assumption that photometric observations at two epochs for a large number of quasars will reveal the same statistical properties as wellsampled light curves for individual objects. This assumption has been tested by MacLeod et al. (2008) using light curves for spectroscopically confirmed quasars observed roughly 50 times over 8 years in SDSS Stripe 82 (S82). They found that while the mean SF for individual sources is consistent with Eq. 1, the contribution of the mean trends to the observed dispersion in variability properties is minor compared to an intrinsic stochasticity of unknown origin. Further investigation of this stochastic behavior is one of the main goals of this study.
In order to better understand the relationship between the two types of data analyses (individual versus ensemble quasar variability), and to begin linking to physical models, we apply the damped random walk model to the light curves of 9,000 spectroscopically confirmed SDSS S82 quasars. This large sample greatly benefits from the robust, accurate, fiveband SDSS photometry. We estimate the variability parameters following Koz10, who demonstrated that their approach is more statistically powerful than the forecasting methods used by KBS09. We also note that both the Koz10 and KBS09 approaches are much faster than that used by Schmidt et al. (2010), requiring only O() rather than O() operations to determine the model parameters for a light curve with data points.
In Section 2, we describe the model, define our variability parameters, and demonstrate their relationship to those utilized in previous studies. In Section 3, we introduce the S82 data set and outline our initial light curve selection. In Section 4, we present the bestfit variability parameters for our final sample of light curves and estimate their scatter due to the limited time sampling of SDSS. We also estimate the sensitivity of our results to variations in the slope of the model power spectrum on long time scales. In Section 5, we describe the relationship between the longterm variability parameters and physical parameters such as wavelength, absolute magnitude, black hole mass, and Eddington ratio. Using these results, we also provide a prescription for simulating mock quasar light curves. In Section 6, we explore the variability properties of subsamples detected at radio and Xray wavelengths. Finally, we summarize our results in Section 7.
2 Methodology
We model the time variability of quasars as a stochastic process described by the exponential covariance matrix
(2) 
between times and . As detailed by KBS09 and Koz10,
this corresponds to a damped random walk (more specifically, an
OrnsteinUhlenbeck process) with a damping time scale , also
called the characteristic time scale, and a longterm standard
deviation of variability
2.1 Structure Function for the Damped Random Walk Model
For our analysis, we express the longterm variability in terms of the
structure function (SF) in order to relate to our previous studies and
to those of twoepoch samples.
The first order SF is (e.g., Hughes et
al. 1992), where the autocorrelation function for a damped random walk is
(Eq. A6 in KBS09). This ACF results in the structure function
(3) 
Asymptotic values of the SF at large and small are
(4) 
The form is equivalent to a power spectral distribution , where (see Appendix of Bauer et al. 2009; KBS09). The SF at small time lags is therefore equivalent to a power spectral distribution . We adopt and as our two main variability model parameters.
2.2 Model Light Curves
Equipped with a statistical description of quasar variability, we generate wellsampled light curves in order to 1) demonstrate the relationship between our variability parameters and the traditional SF analyses of many previous works, and 2) to estimate the systematic effects that the sampling rate and light curve length have on the fitted parameters. The latter is especially important because the S82 data are fairly sparse. As shown below and in Section 4.2, these indeed have a large impact.
A light curve is generated using only three input parameters: , , and the mean value of the light curve, . The magnitude at a given timestep from a previous value is drawn from a normal distribution with a mean and variance given by
(5) 
(Eqs. A4 and A5 in KBS09). The asymptotic variance of the time series is then . The top panel of Figure 1 shows a segment of a wellsampled light curve generated using days, mag, and a time sampling of . The structure function, , is computed by collecting the differences in magnitude for all points in the light curve separated by a given time lag, . The distribution of magnitude differences () is Gaussian by construction, and the rms of the distribution is the value for that time lag.
When fitting values for and for a given light curve, the length of time that it spans plays an important role. For example, in Figure 1, the SF computed for the full light curve length of (82 years) is much smoother than that computed for three equal sections, each spanning 27 years. Therefore, the determination of variability parameters for S82 quasars will be affected by their light curve lengths, which are typically . When the light curve is too short, it is easy to overestimate because there is no information on the time scale of the break (where flattens to ). The model will reproduce the observed variance in the light curve by overestimating , and therefore is the more robustly estimated model parameter when cannot be welldetermined.
2.3 Comparison with Published Work
Before interpreting the form of , we summarize the major differences between our S82 analysis and previous studies based on ensemble structure functions. The ensemble structure function is computed using only a few observations of many quasars, combining all magnitude differences to find . Using an ensemble approach is beneficial because it enables one to constrain the average variability properties when it is difficult to constrain such quantities for individual quasars. Indeed, even with wellsampled light curves, spurious breaks in the individual SFs are common (Emmanoulopoulos, McHardy, & Uttley 2010). However, in previous works (e.g., I04; de Vries et al. 2005), the characteristic time scale is defined as the time lag at which the ensemble flattens to a constant value, and thus it represents some complex average over the intrinsic distribution. In contrast, by applying a stochastic model to the individual S82 light curves, we are relatively insensitive to time sampling issues (for details, see KBS09), and we obtain a model fit for every quasar in each filter described by the parameters and .
A power law fit to has previously been a common way to describe , and even to reject certain classes of models (e.g., Kawaguchi et al. 1998; VB04; I04). However, Figure 1 shows that the bestfit power law index is extremely sensitive to the fitted range of . For example, for , is well fit by a power law with an index of 0.5, while for , we obtain a strongly biased power law index of . Furthermore, each quasar has its own values of and , and the ensemble structure function is a convolution of the individual structure functions with their distribution in parameters,
(6) 
where (Eq. 3) is the structure function at time for a quasar with variability parameters and . Hence, the ensemble structure function is only indirectly related to the structure function for any particular quasar, and results based on fitting a power law to observed ensemble structure functions should be interpreted with caution.
3 The SDSS Stripe 82 Quasar Data Set
The Sloan Digital Sky Survey (SDSS, York et al. 2000) provides homogeneous and deep () photometry in five passbands (, Fukugita et al. 1996; Gunn et al. 1998; Smith et al. 2002) accurate to 0.02 mag, of almost 12,000 deg in the Northern galactic cap, and a smaller, but deeper, survey of 290 deg in the Southern galactic hemisphere. For this 290 deg area known as Stripe 82 (S82), there are on average more than 60 available epochs of observations. These data were obtained in yearly “seasons” about 23 months long over the last decade, and the cadence effectively samples time scales from days to years. The light curve lengths are effectively shorter than the actual period of the survey because the bettersampled supernova observations begin about 5 years into the survey. Because some observations were obtained in nonphotometric conditions, improved calibration techniques have been applied to SDSS S82 data by Ivezić et al. (2007) and Sesar et al. (2007), and we use their results. For these data, photometric zeropoint errors are 0.01–0.02 mag.
We have compiled a sample of 9,275 spectroscopically confirmed quasars in S82 with recalibrated light curves (see also Bhatti et al. 2010). Most (8,974) of these are in the SDSS Data Release 5 (DR5) Quasar Catalog (Schneider et al. 2007), and the remaining are newly confirmed DR7 (Abazajian et al. 2009) quasars. Summed over all bands and epochs, the data set includes 2.7 million photometric measurements. For 41% of the sample, the random photometric errors are smaller than 0.03 mag. Only 1% have errors mag in , , and , and 2.4% have errors exceeding 0.25 mag in and filters. In MacLeod et al. (2010, in prep.), these light curves, as well as a much larger sample of quasars with two SDSS epochs selected from 12,000 deg of the sky, will be made publicly available. We adopt the Kcorrected iband absolute magnitudes from Schneider et al. (2007), and virial black hole masses and bolometric luminosities where available from Shen et al. (2008). The Shen et al. masses were estimated from emission line widths (H for , MgII for , and CIV for ). However, we note that at low spectroscopic signaltonoise, black hole masses tend to be overestimated (Denney et al. 2009).
3.1 Initial Light Curve Selection
The damped random walk model was fit to all available light curves for S82 quasars. Summed over 5 bands, there are bestfit values of the characteristic (damping) time scale and longterm structure function . For further analysis, we select light curves that satisfy the following criteria:

First, we remove light curves with fewer than 10 observations. The topleft panel of Figure 2 shows the distribution of the number of observations () per light curve before this cut, which reduces our sample to . At a given , the distribution of the ratio of light curve length to is similar to that at all other values of . Therefore, any systematic effects of on derived parameter distributions should be small.

We then require that the stochastic model must provide a better fit than uncorrelated noise. Following Koz10, we select light curves with a likelihood improvement over simply broadening the measurement errors of , where is the likelihood of the stochastic model and is that for the noise solution where . The topright panel of Figure 2 shows the distribution of before this cut, which removes 14% of and light curves, whose photometric errors are larger. About 7% of our light curves (82% of which are or band) are removed in this step
^{20} , reducing our total to . 
Finally, we remove cases where is merely a lower limit due to the length of the light curve. We define , where is the likelihood that , indicating that the light curve length is too short to accurately measure . The bottomleft panel shows a peak at , and we exclude these objects by requiring that . Most (95%) of the rejected light curves have lengths ; 76% have days, and 64% have mag. The latter is due to the fact that as the value becomes long and uncertain, the model will necessarily overestimate in order to keep the overall light curve rms fixed (see Section 2.2).
The rejected light curves tend to be higher redshift quasars because stronger time dilation leads to shorter restframe light curve lengths, making it increasingly difficult to constrain long restframe . This criterion removes 22% of our sample, leaving a total of values of and . Because we are limited by the duration of the S82 survey, this is a significant loss, and therefore our final distribution is biased low, but the bias should not be significant considering our results in Section 4.2.
The resulting light curves are wellfit by the stochastic model, as can be seen from the distribution of shown in the bottomright panel in Figure 2, where is the number of degrees of freedom. The expected Gaussian distribution with rms is also shown in the panel, where we have averaged over the distribution of the light curves. The observed distribution is centered at . This difference is some combination of errors in the estimated errors, outliers in the light curves, and any poorly modeled physics. Koz10 noted a similar difference in their analysis of OGLE light curves. Only 5% of the light curves have , confirming that most quasars are variable (at the level), and that a damped random walk is a good description of quasar variability.
4 Variability Parameters for Stripe 82 Quasars
4.1 Observed Distributions
We assume the variability to be intrinsic to the quasars and convert the time scales to the rest frame (dividing by ) before further analysis. The longterm structure function should be independent of redshift other than through evolution in physical parameters, and this view is confirmed by KBS09. Figure 3 shows the distributions in and restframe found for the S82 quasars. If we consider only the brighter () quasars, based on the bestfit mean magnitude, the distributions are generally similar but biased towards lower asymptotic amplitudes (peaked at 0.12 mag). The distributions are consistent with what was found in KBS09 and Koz10; however, in these studies, there are not as many objects with runaway () time scales because of the improved time sampling of the OGLE survey and many of the light curves in the KBS09 sample. Also, the observedframe distribution lacks many of the short time scales observed in Koz10; this is likely due to either the better time sampling of the OGLE light curves or unrecognized stellar contamination in their sample.
The bottom panel of Figure 3 shows that the bestfit variability parameters are highly correlated with each other, indicating that quasars with larger asymptotic amplitudes of variability also have longer characteristic time scales. We fit a power law slope of dex/dex for all 33,112 data points, and this trend persists within each band as well. Note that a correlation between and is expected even if is independent of the driving amplitude of shortterm variations, . Since the power law slope of is steeper than that expected if and are independent (0.2, see below), the time scale must also be correlated with the amplitude of shortterm variability. We confirm that this correlation is intrinsic, rather than an artifact of short light curve lengths, in Section 4.2.
The approximate values for and can often be guessed from light curves by simple visual inspection. Figure 4 shows representative light curves from several regions of space, indicated by stars in Figure 3. The weighted average of all model light curves consistent with the data is also shown, and the “error snake” is the range of those light curves about this mean. The top panel shows a light curve with a relatively short – this is due to the large amount of variability within each season (i.e. each grouping of points), as compared to that in the second panel, which shows a curve with a much larger . The third panel shows a light curve with a relatively low , while the bottom panel shows one with a larger , as can be seen by the larger difference in median band magnitude between about 10 and 11 years. Note the outlier in the third panel: outliers such as these are generally responsible for higher values of . The (datamodel) brightness difference provides a convenient way to identify outliers, and they may be an indicator of variability behavior not captured by the damped random walk model. However, they may also be caused by occasional nonGaussian photometric errors and thus their analysis requires a detailed and careful study (e.g., by utilizing control samples of appropriately chosen nearby nonvariable stars). Since the model provides satisfactory fits for the overwhelming majority of objects (see Section 3.1), we do not further investigate such outliers in this work. We have also searched for periodic signals in observed light curves and did not find any convincing cases. This analysis is summarized in the Appendix.
4.2 The Effect of S82 Time Sampling and Estimate of Fitting Errors
There are three contributions to the scatter in the bestfit variability parameters and . First, the fitting errors, including those due to insufficient time sampling and light curve length, will introduce some scatter. Second, trends with physical parameters (see Section 5) result in a certain distribution width. There may also be scatter due to other sources of variability that are not captured by the model, such as flares, or other activity related to radio emission, for example. Here, we carry out two tests in order to understand how the bestfit parameters are affected by the limited data sampling of S82.
Since the correlation between and seen in Figure 3 might be expected if the light curves are not sampled over sufficiently long periods of time (see Section 2.2), we first test whether the correlation is real or simply an artifact. We generated light curves as described in Section 2.2 using the time sampling and photometric uncertainties of the S82 data. For each object, the input parameters and are randomly drawn from a uniform distribution in and limited by the dotted rectangle in the left panel of Figure 5. These artificial light curves are then fit to obtain and , and the resulting distribution is shown by the contours in Figure 5. The open circles show data points that do not satisfy (3% of all input points) – these are concentrated at small and . The closed circles show those that do not satisfy (21% of all input points) – these have been smeared to large values of and . After omitting points with and , the distribution of the output estimates is similar to the input distribution and shows no strong correlation between and . This suggests that the correlation in Figure 3 is largely real, and not an artifact of sampling and fitting procedures. These results also justify the selection cuts outlined in Section 3.1. The right panel of Figure 5 shows the expected correlation that the larger the overestimate of , the larger the overestimate of , with a slope between them following that expected by Eq. 4 (for ). We repeated the test using a uniform input distribution in (the driving amplitude of shortterm variations) rather than in , where . In this case, we find the output and are correlated with a power law slope of . Since this slope is smaller than that for the observed distribution in Figure 3, the and must be intrinsically correlated for the S82 sample as well.
For our second test, we used the bestfit , , and for the S82 sample to generate new light curves with the same time sampling and photometric uncertainties as the S82 data. By comparing the output and input parameter distributions, we can estimate how much the intrinsic stochasticity and the time sampling issues affect the results. Figure 6 compares an observed and “regenerated” light curve, where the differences are due to the stochastic nature of the process. The fit parameters for the “regenerated” light curve can be very different because of how the particular realization is affected by the time sampling, as illustrated in Figure 6. Figure 7 shows the ratio of output to input distributions for both and . The input distributions normalized by their median values are also shown to illustrate their dynamic range. These two distributions are compared to each other in order to estimate the effect of fitting errors (the ratio of output to input should be a delta function centered at 1 for perfect time sampling). The bottomright panel shows that the correlation between and becomes slightly weaker than that seen in Figure 3. Based on these results, we conclude that the uncertainties due to sparse sampling and limited lengths of the S82 light curves can account for 71% of the spread in and 57% of that in . As shown in Section 2.2, very long light curve lengths are needed to estimate accurate time scales and asymptotic amplitudes. Nevertheless, the observed distributions indicate that the underlying intrinsic distributions of and have finite widths that are similar to the observed widths.
4.3 Relationship between the Individual SFs and the Ensemble SF
The distribution of is important to consider when relating the ensemble SF, such as those determined using twoepoch datasets, to the SFs for individual light curves. The analysis based on twoepoch data measures the mean value of the SF, but provides no information about the SF variance among individual objects. To measure the latter, individual light curves must be available.
MacLeod et al. (2008) analyzed individual light curves for S82 quasars in order to test the common assumption that photometric observations at two epochs for a large number of quasars will reveal the same statistical properties seen in light curves for individual objects. They found that the dependence of the mean SF computed using SFs for individual light curves on luminosity, restframe wavelength and time lag is indeed qualitatively and quantitatively similar to that derived from twoepoch observations of a much larger sample. However, they also found that the scatter in the lightcurve based SFs for fixed values of , and is very large, and in fact, similar to the scatter for the whole sample (see Figure 8). This large scatter was attributed to an intrinsic stochasticity of unknown origin. The new modelbased analysis discussed here allows us to explain this puzzling result as a consequence of the finite width of the distribution.
An important piece of information missing from the MacLeod et al. (2008) analysis is the existence of a characteristic time scale . In MacLeod et al. (2008), the individual structure functions () were computed following the standard approach with a fixed observer’s frame time lag of 1 year (and with restframe time lags spanning 100 to 300 days). However, Eq. 4 indicates that the SF at small time lags should be proportional to , indicating that the observed structure functions will vary between quasars even if they have similar . Indeed, the distribution of has an almost identical shape and width as that of (see Figure 8). We therefore conclude that variations in are responsible for most of the scatter in for quasars with similar luminosity, rest wavelength, and time lag. In addition, it is likely that the distribution is responsible for the exponential tails of the magnitude difference distribution for quasars reported by I04 and Sesar et al. (2006; see MacLeod et al. 2010, in prep., for further discussion). Therefore, the published structure function results based on twoepoch datasets can only be interpreted in the context of Eq. 6.
4.4 Dependence on the Underlying PSD
The analysis throughout this paper assumes that the variability is described by a damped random walk, which has a power spectral distribution (PSD) described by at frequencies , flattening to a constant at lower frequencies. In this Section, we investigate the sensitivity of our resulting parameter distributions to the possibility that the low frequency part of the PSD is not flat. For example, the Xray variability of Seyfert galaxies is welldescribed by a broken power law with a slope of at high frequencies, breaking to a shallower slope () at low frequencies (Arévalo et al. 2008a). Since the optical and Xray variations are correlated (Arévalo et al. 2008b), it is plausible that the optical variability might have a similar underlying PSD. On the other hand, as reviewed by McHardy (2010), the optical and Xray fluctuations are only correlated on time scales which are shorter than the typical characteristic optical time scale (i.e., on time scales corresponding to the part of the optical power spectra). If the damped random walk is a good description for the optical variability, then we might expect that the optical and Xray fluctuations are no longer correlated on time scales longer than , as the optical fluctuations resemble white noise on these time scales.
We consider three cases for a broken, or bending, PSD described by at high frequencies and at low frequencies. The first case is where , which is a damped random walk. The second case is where , and the third is where , which is nearly a constant powerlaw slope. In each case, 7000 light curves are simulated using the algorithm from Timmer & Koenig (1995) with the chosen PSD. For each realization, the break frequency is set to and the total is fixed, using the observed band and values for the S82 quasars that satisfy the selection criteria in Section 3.1. Therefore, the observed distributions in Figure 3 (for the band) are the “input” values. Figure 9 shows 3 example light curves simulated using (red), (blue), or (green). The lines in the last panel show each PSD from which the light curves are generated, with the black dotted line indicating the break frequency. As a check, the PSD was computed for each simulated light curve, and as seen by the colored dots, the shapes match the input.
Each light curve was simulated over 100 years and then truncated to the 30–40 year segment in order to account for the additional variability and bias in the inferred , which may result from a red noise leak or an incorrect specification of the PSD model (see Uttley et al. 2002). The simulated light curves were then modeled as a damped random walk to obtain the (“output”) parameter distributions shown in Figures 10 and 11. For Figure 10, the simulated observations are spaced every 5 days over 10 years with typical errors of 0.01 mag, and for Figure 11, the S82 window function is imposed (i.e., all light curves have the S82 time sampling and photometric accuracy). The filled histograms show the input distributions and are the same for both Figures.
It is clear from comparing the input and output parameter distributions that and map simply onto the values of the powerspectral break time scale and amplitude. Indeed, the correlation between and is preserved in the output distributions, as seen in the bottomright panels of Figures 10 and 11. However, as seen in the top and bottomleft panels, as the PSD slope at low steepens to , the number of “runaway” time scales (where saturates at days) increase due to that fact that can no longer be constrained. Whereas for the S82 data, 20% of light curves were rejected as runaway cases (see step 3 of Section 3.1), for an PSD, the runaway fraction is 30%. This significant increase rules out the PSD as the correct model, as we do not see this large runaway fraction in the data. The fractions are similar for the and cases, suggesting that the correct model has . The increasing fraction of runaway with decreasing allows one to distinguish between each PSD in the wellsampled case, as seen in the topright panel of Figure 10. In this panel, the difference between the for the bestfit damped random walk and that for a solution, , is shown for each input PSD. The total distribution shows two peaks; that on the left corresponds to cases where the model is able to constrain . For wellsampled light curves (Figure 10), the case with (dashed line) yields an overall lower and therefore a more likely damped random walk fit, as expected. However, for the S82 sampling, it is difficult to distinguish between the red and blue lines ( and ), and the (input) observed distribution is similar to both. Therefore, within the limited S82 sampling, we are unable to distinguish reliably between a damped random walk and a PSD on long time scales.
5 Dependence of Variability Parameters on Luminosity, Wavelength, Redshift, and Black Hole Mass
Next, we discuss correlations between each of the variability parameters, and , and the four available physical parameters: restframe wavelength (), redshift (), absolute iband magnitude (), and black hole mass (). It is important to fit a multiple regression to each of the physical parameters, functions of the form , because of the correlations between physical parameters. For example, when searching for a correlation between and , one must take into account that a more luminous quasar hosts a more massive black hole (e.g., Kollmeier et al. 2006), or else a trend with luminosity may be mistaken for a trend with black hole mass. A similar example is the wellknown luminosityredshift (–) degeneracy seen in fluxlimited samples of quasars. Magnitude limits result in the illusion that only the most luminous quasars are seen at high , and therefore the observed increases with independent of reality. Having a large sample size helps to alleviate these degeneracies. For example, one can look for trends using 2dimensional grids in any two physical parameters of interest. The numbers of data points per twodimensional bin of and redshift, and similarly for and , are shown in the bottom panels of Figure 12. The two left panels of Figure 12 show the selection effect that quasars at higher redshift must have higher luminosities to be included in the survey. From the twodimensional distribution in the topleft panel, we can see that this – degeneracy is essentially independent of . Furthermore, shorter restframe wavelengths are probed at higher redshifts for two reasons: first, quasars emitting at shorter rest wavelengths must be at higher redshifts in order to be observed within the filters, and second, quasars at high redshifts are closer to their Eddington limit as a result of the – degeneracy, and possibly cosmological downsizing (e.g., Kollmeier et al. 2006). Therefore, any dependence of variability on wavelength must be accounted for when considering a dependence with redshift (or luminosity).
When fitting power laws to a large number of data points throughout this paper, we fit to the median values in each bin of the independent variable, where all bins have the same number of data points . The errors in the medians are computed as 0.93(IQR), where IQR is the 25% – 75% interquartile range, and these constrain our formal errors in the power law slopes.
5.1 Trends with Restframe Wavelength
We start by examining the wavelength dependence of the variability parameters. Since there are multiple bands for each quasar, this dependence can be determined for individual sources for which , , and are fixed. The restframe wavelength is found by dividing the observed bandpasses (, , , , and Å for , respectively) by . We fit a power law to the estimates of and for every quasar observed in at least two filters (8,000 quasars). The median values are and for and , respectively. We searched for significant correlations between and other physical parameters and did not find any. We use these median values to correct for the wavelength dependence of and before searching for their correlations with other physical parameters in Sections 5.2, 5.3, and 6.
This method of fixing before investigating other correlations naturally eliminates any degeneracies between rest wavelength and the other physical parameters. This is especially important in the case of the variability amplitude, . Figure 13 shows how the two variability parameters vary with in each of the filters. The median power law from the above analysis, shown as a straight line stretching across all wavelengths, accurately traces the overall trend. However, when the data are fit to ensembles of quasars in each band separately, the slopes for , shown by the short lines for each filter, are very different (0.4). This difference is a consequence of the correlation between and : for a given band, shorter corresponds to higher , but at higher , quasars have higher and thus smaller (see below), creating a bias in the inferred slope of the wavelength dependence.
5.2 Trends with Luminosity, Redshift, and Black Hole Mass
In the upper panels of Figure 14, the median values of and are shown as a function of absolute magnitude () and redshift. The structure function parameters are normalized to a fixed rest wavelength of 4000Å using the fitted power law dependencies of with and 0.17 for and , respectively, from Section 5.1, before finding the median in each pixel. For , the anticorrelation with luminosity clearly dominates any trend with redshift. The distribution also shows negligible correlation with redshift. The bottom panels show the dependence on . Using a grid of vs. allows one to search for trends in the variability parameters while accounting for the selection effect that more luminous galaxies host more massive black holes. A positive correlation between and is apparent, independent of the correlation with . This is in agreement with the result from Wilhite et al. (2008), who used ensemble structure functions. The parameter shows a clear correlation with , which dominates any trend with .
Motivated by these qualitative results, we fit a power law of the form
(7) 
to all and data points in each band separately, keeping fixed to and , respectively, in order to avoid the bias discussed in Section 5.1. While there is a lot of scatter in the variability parameters for fixed , , and (see Section 4.2), Eq. 7 describes the mean trends well. The bestfit coefficients, averaged over the five bands, are reported in the first and sixth rows of Table 1. The bestfit coefficients when simultaneously fitting all bands are consistent with these averages. The reported error bars are computed from the variation of the bestfit parameters over the five bands.
The dependence of on redshift and is only marginally detected, and can be attributed to the – degeneracy: as redshift increases, the bestfit decreases, while the coefficient for indicates that as luminosity increases, increases. When is fixed to zero, the dependence on has low significance, while the correlation with remains the same (seventh row of Table 1). For illustration, when is fixed to 1 [so , third and eighth rows], a spurious dependence on emerges. Therefore, the controlling variable for the characteristic time scale is clearly , suggesting that more massive black holes vary on longer time scales. Similarly, with fixed to 0 for (second row), the () and () coefficients remain largely unchanged, confirming that there is no significant correlation between amplitude and redshift. Therefore, we force the final adopted model to have no redshift dependence (). Table 1 also provides fits for when is fixed to zero so that and can be estimated in the absence of black hole mass information (fourth and ninth rows).
Of the analyzed parameters, the probable errors for black hole mass estimates are the largest: of order 0.2 – 0.4 dex (Marconi et al. 2008; Vestergaard & Peterson 2006). When ignored, these large statistical uncertainties result in underestimated values for (Kelly 2007). With an assumed random uncertainty of 0.2 dex in the black hole mass measurements, we find the bestfit coefficients reported in the fifth and tenth rows of Table 1 (continuing to assume that ). It can be seen that including mass uncertainties of 0.2 dex increases from 0.1 to 0.2 for both and (a bias of about 23 formal statistical fitting errors). The coefficients in the fifth and tenth rows represent our bestfit model for the variability parameters, and we assume these values in the remaining Sections.
The bestfit model shows a smaller slope for the correlation between and of than the value of found by KBS09. The two results are still marginally statistically consistent given the large uncertainties in KBS09. Moreover, their sample contains relatively lower mass and luminosity quasars than those analyzed in this study, and this difference might result in additional biases. Koz10 also noted that the KBS09 results may be affected by contamination from host galaxy emission, a problem which will be far smaller for the generally more luminous SDSS quasars.
()  ()  ()  ()  







Note. – In each row, the coefficient was determined and fixed before fitting a multiple regression in all other parameters (see Section 5.1). The cosmology used for determining is , , and (Schneider et al. 2007), whereas that used in the estimates is , , and (Shen et al. 2008). This difference should only have a 1% effect on the bestfit coefficients.
5.3 The Eddington Ratio as the Driver of Variability?
Since the Eddington ratio, , is dependent on luminosity and black hole mass, the trends for in Table 1 might be explained if is simply driven by . To test this, we estimated as the ratio of the bolometric luminosities from Shen et al. (2008) to the Eddington luminosity, erg/s. Figure 15 shows as a function of and , and demonstrates that lines of constant are similar in slope to those of constant (see Figure 14). The median versus the median for each bin is shown in the right panel, where we find a power law slope of . This significant anticorrelation was also found by Wilhite et al. (2008), Bauer et al. (2009), and Ai et al. (2010). KBS09 did not report such an anticorrelation; however, they did not compare with but rather with , the driving amplitude of short term variability. If is the sole driver of the quasar variability amplitude , we would expect that the coefficients for and are related as . However, we find , 4.7 away from the presumed slope of 2.5. This suggests an additional source of variability, such as a dependence on luminosity or on in addition to the dependence on . Moreover, if we ignore this additional source, the result that depends on supports a scenario where the amplitude of the optical variability is determined by the accretion rate (see discussion in Wilhite et al. 2008).
If the Eddington ratio is a proxy for the quasar’s age (e.g., Martini & Schneider 2003), then a lower Eddington ratio, and thus a larger amplitude of variability, could indicate a dwindling fuel supply and a more variable rate at which it is supplied to the black hole. However, it is unlikely that the observed variability is due to fluctuations in the external fuel supply to the disk, because the fluctuation time scales of 10 – 10,000 days are very short compared to the viscous time scales of – days that should control large scale changes in the accretion rate. Moreover, these fluctuations in the fueling rate will also be smoothed out and damped as they travel inward, and will likely have effectively disappeared by the time they reach the optical emitting region (e.g., see discussion in Churazov et al. 2001). Instead, the origin of the fluctuations is probably more local.
Another possibility is that the dependence on may simply be a reflection of the dependence on wavelength, which in turn depends on the disk radius. If a higher means a hotter disk, then the optical flux originates at a larger radius. Assuming that longer wavelengths are emitted further out in the disk, the lower at longer wavelengths would lead to the anticorrelation between and . In thin disk theory (e.g., see Frank, King, & Raine 2002), the characteristic radius for emission at wavelength scales as , and the thermal time scale is related to the orbital time scale as . Therefore, under the assumption that is related to the thermal time scale, and that variability at wavelength is dominated by the scale , scales with , , and as:
(8) 
Since , this means scales as
(9) 
which does not match the observed scaling of . If we assume that is related to the viscous time scale,
(10) 
still in conflict with the measured values.
The variability time scale simply does not show the strong dependence on , , and expected from these simple scalings. Therefore, using this naive scaling, we are not able to relate the observed to either a thermal or viscous time scale of the radius associated with the wavelength of the variability. However, in reality there is a range of radii, which corresponds to a range in time scales, contributing to the observed flux in each band, and this will cause some degree of smoothing. Also, the radial regions might overlap for each band, causing a single radius to contribute flux in multiple bandpasses, and this might attenuate some of the dependencies on , , and . Finally, our observed power law indices may be biased due to the uncertainty in the bolometric correction.
5.4 A “Recipe” for Generating Mock Light Curves
Given the success of the damped random walk model in explaining a large body of observations, it represents a simple quantitative framework for generating mock quasar light curves. The optical light curves can be simulated for quasars at different redshifts and for a wide range of luminosity and black hole mass, which provides the basis for quantitative modeling and optimization of quasar variability surveys. The model presented has already been implemented as a part of the simulation effort in support of the Large Synoptic Survey Telescope (Ivezić et al. 2008).
At a given redshift, a simulated value for absolute magnitude can be readily drawn from an adopted luminosity function (e.g., Croom et al. 2009 and references therein). An estimate for the black hole mass is also required in order to apply our model, and it can be generated using the adopted absolute magnitude as follows. According to the Shen et al. (2008) results, the quasar luminosity and black hole mass are strongly correlated (see the bottomright panel of Figure 12). Using Shen et al. (2008) values, we quantify this correlation in Figure 16. At a given , the distribution has a finite width due to the distribution of Eddington ratios and measurement errors. Both effects can be welldescribed using a Gaussian distribution:
(11) 
where and (the black hole mass is expressed in solar mass units).
With adopted values for and , the variability amplitude and the characteristic time scale can be estimated using Eq. 7. We note that choosing a time scale from an adopted value for by utilizing the correlation seen in Figure 3 (and drawing from a Gaussian distribution similar to Eq. 5.4) is not as accurate as using Eq. 7 and the adopted values of and . Finally, given a mean magnitude and wavelength, a quasar light curve can be easily generated using Eq. 5.
6 Variability Properties of Radio and Xray–detected Quasars
Previous studies, based mostly on twoepoch datasets, found that radioloud quasars are marginally more variable than radioquiet quasars for restframe time lags in the range 50–400 days (e.g., VB04 and references therein). It was also found that Xray–detected quasars are significantly more variable than Xray–undetected quasars at restframe time lags up to 250 days.
The large size and the availability of light curves for our sample allows us to revisit these results with more statistical power. In addition, the damped random walk model and our bestfit results from Table 1 provide a convenient method to account for various selection effects. For example, it is possible that subsamples selected by various means, such as the requirement of radio or Xray detections, have different distributions of luminosity and black hole mass. In this case, they would display different variability behavior not because of intrinsic differences in the variability mechanism, but rather because of the trends captured by Eq. 7. A comparison of the ratio of observed and predicted variability parameters for any subsample and the full sample will automatically take into account these sampling effects, and this is the main statistical method we use in this section. We first discuss a radio subsample, and then analyze a sample of quasars with Xray detections.
6.1 Radiodetected Subsample
We use the unified radio catalog of Kimball & Ivezić (2008) to access the FIRST and NVSS data (20 cm continuum data) for our objects. We refer the reader to this paper and references therein for details about these radio surveys and object association. Figure 17 shows the radiodetected fraction of quasars in S82 as a function of magnitude. The overall fraction (5%) is considerably lower than that in the Northern galactic cap footprint of the SDSS for the magnitude range (White et al. 1997; Ivezić et al. 2002). This is due to a difference in targeting algorithms between the two surveys, and results in nearly 300 quasars with both optical light curves and radio data. When available, we use NVSS values for radio flux, but otherwise we adopt the FIRST integrated flux. The radioloudness parameter is calculated exactly as the in Ivezić et al. (2002), but using magnitudes instead. We assume a spectral index of for the optical and 0 for the radio (because this sample is dominated by the core radio emission; see Kimball & Ivezić 2008) when computing the Kcorrection.
Table 2 compares the mean properties for various subsamples detected in the radio to the radioundetected subsample (essentially the full sample due to small fraction of radiodetected objects). For each subsample, the mean and are listed, as well as the mean ratio of the observed values to those predicted using Eq. 7 and the measured values of , , and . The mean values are computed iteratively with outliers excluded. Errors are reported as the (clipped) rms divided by , where N is the number of data points in the distribution. Numbers of objects are lower for columns involving model quantities because black hole mass estimates are not available for all objects.
Table 2 shows that only the variability amplitude of the radioloudest quasar subsample () is significantly different ( deviation) from the behavior of the full sample. The radioloudest quasars have systematically larger variability amplitudes, when corrected for trends described by Eq. 7, by about 30%, compared to the full sample dominated by radioquiet objects, in agreement with VB04.
6.2 Xray Detected Subsample
For the analysis of variability properties of quasars detected at Xray wavelengths, we use data from the ROSAT AllSky Survey (RASS; Voges et al. 1999). The Xray subsample consists of 82 quasars with RASS fullband count rates greater than ct/s, taken from the Schneider et al. (2007) catalog. As can be seen from Table 2, the variability properties of this subsample are statistically indistinguishable from the full sample. VB04 detected a significant increase in structure function at rest time lags below 250 days for their Xray subsample, but on long time scales there was no significant difference. It is plausible that their result for short time scales was influenced by small sample size and effects discussed in Section 4.3.
no radio  6989  4467  7067  4508  
no radio  1519  1279  1527  1287  
radio  277  154  283  154  
radio  105  73  108  74  
135  68  137  68  
211  133  219  140  
30  19  30  19  
resolved 
111  63  111  63  
unresolved 
167  91  172  91  
xray  81  59  82  58  
no xray  6950  4559  7020  4598 
Note. – and are the observed time scales and asymptotic amplitudes of optical variability (Section 4), while and refer to those predicted from Eq. 7, using the coefficients in the fifth and the tenth rows of Table 1. Errors are reported as the rms divided by , where N is the number of data points, listed after each column. Numbers of objects are lower for columns involving model estimates because black hole mass estimates are not available for all objects. Only the variability amplitude for the radioloudest quasar subsample () is significantly different ( deviation) from the behavior of the full sample.
7 Summary and Conclusions
We have used the damped random walk model of KBS09 and Koz10 to model the optical/UV variability of 9000 SDSS Stripe 82 quasars with the light curves. The dataset includes 2.7 million photometric measurements collected over 10 years. We confirm that this is a good model of quasar variability, and quantify the dependence of two variability parameters, the longterm rms variability and the damping time scale , on physical parameters such as wavelength, luminosity, black hole mass and Eddington ratio. Our main results are the following:

A stochastic process with an exponential covariance function characterized by an amplitude and time scale provides a good fit to observed quasar light curves, as shown by KBS09 and Koz10, using smaller samples with less wavelength coverage, but better time sampling.

The longterm rms variability, , has a mode at 0.2 mag and characteristic time scales, , are roughly 200 days in the rest frame, as found previously by KBS09 and Koz10. These time scales are consistent with thermal time scales, but simple accretion disk models fail to reproduce the observed scaling of with physical parameters.

Quasars with similar physical parameters can have different characteristic time scales for variability. It is now clear that the distribution of accounts for most of the scatter in the structure function on short time scales for quasars with similar luminosity, rest wavelength, and time lag, which explains the puzzling results from MacLeod et al. (2008). Results from fitting a power law to observed ensemble structure functions should be interpreted with caution.

The variability time scale is correlated with the longterm rms variability with a slope of dex/dex. Quasars that have large longterm variability amplitudes generally vary on longer characteristic time scales. The amplitude of short term variations is also correlated with . This conclusion is unaffected by any time sampling issues in the S82 dataset.

The damped random walk corresponds to a PSD proportional to at frequencies , flattening to a constant at lower frequencies. At large , the data are in great agreement with . In terms of the structure function, this means that . Whereas previous analyses of the SF obtained a power law slope of , here we demonstrated that this would be a consequence of fitting the data around the “knee” (turnover) of the SF. Our constraints for small are much weaker. As discussed in Section 4.4, due to a lack of sufficient longtime scale information, we are unable to distinguish between a or a PSD at frequencies using the data and computational technique described here.

The restframe variability parameters show a negligible trend with redshift, suggesting that they are intrinsic to the quasars, and these properties do not evolve over cosmic time for fixed physical parameters of the quasar (, , and ).

For fixed luminosity and black hole mass, increases with increasing restframe wavelength with a power law index of , and decreases with a power law index of . The latter result is similar to previous findings (e.g., MacLeod et al. 2008, VB04). Koz10 also observed that the variability increases to shorter , but they kept fixed in their fits. If wavelength is a proxy for radius in the accretion disk, this implies that the characteristic time scales are longer and the variability amplitudes are smaller in the outer regions than in the inner regions.

The longterm variability is strongly anticorrelated with luminosity as found in previous studies such as VB04, Wilhite et al. (2008), and references therein. By studying the median in the plane of absolute magnitude and black hole mass, we can separate the anticorrelation of amplitude with luminosity from the positive correlation with black hole mass. As suggested in Wilhite et al. (2008), these trends may be largely explained if the amplitude of variability is tied to changes in the accretion rate in the disk, and is simply related to the Eddington ratio. However, despite the strong anticorrelation between and (which accounts for most of the dependence on and ), the exact dependence with and is not consistent with as the sole driver of quasar variability.

The damping time scale appears to be nearly independent of luminosity and correlated with with a power law index of . The mild discrepancy with the KBS09 result (1.00.4) may be due to the different range of sampled luminosity and black hole mass, as well as contamination by host galaxy emission in many of the very low luminosity systems they consider (see Koz10).

While the mean variability parameters can be related to physical parameters, for fixed values of , , and , there is still a large scatter around the mean values, similar to the variance of the observed distributions. Some of that scatter can be attributed to measurement and fitting errors (60% for and 70% for ), but there is enough evidence for residual stochastic nature of quasar variability. Therefore, it cannot be assumed that quasars with similar , , and will necessarily have similar variability properties.

The radioloudest quasars have systematically larger variability amplitude by about 30%, while the distribution of their characteristic time scale is indistinguishable from that of the full sample. There are no statistically robust differences in the characteristic time scale and variability amplitude between the full sample and a small subsample of quasars detected by ROSAT.
With this paper, and results from KBS09 and Koz10, the ability of the damped random walk model to quantitatively describe quasar variability is wellestablished. As emphasized by Koz10, this means that variability studies can become fully quantitative because the entire process of identifying and assigning parameters to quasars can be simulated to allow estimates of completeness and parameter biases. In particular, an important next step is to determine the variability equivalents of luminosity functions, i.e., the intrinsic distributions of the variability parameters. While our results represent a good first step in this direction, we caution that a nonnegligible fraction of the S82 quasars have indeterminately long time scales. If, following KBS09, we identify the characteristic time scale with the thermal time scale, then the next question is whether there are additional time scales (such as the dynamical or viscous time scale), or sources of variability. For example, Blackburne & Kochanek (2010) recently found evidence for changes in disk size with changes in luminosity using gravitational microlensing. Such searches for additional sources of variability will likely require better sampled light curves and over a longer time scale.
The prospect of advancing these studies of quasar variability now faces a bottleneck. The S82 quasars have the advantage of sample size, wavelength coverage, and spectroscopy, but the light curves have poor sampling and modest overall lengths, leading to significant problems for accurately estimating when it is long. The quasars behind the Magellanic Clouds (Kozłowski & Kochanek 2009) are a smaller sample without spectroscopic confirmation, but have superb, long term light curves that continue to be extended because of the continuing microlensing projects. Improving on our present results in the short term depends on either reviving the monitoring of S82 or spectroscopically confirming the Magellanic Cloud quasars.
Resuming the monitoring of Stripe 82 is challenging with the decommissioning of the SDSS imaging system. The best shortterm prospects are the PanSTARRS project (Kaiser et al. 2002), the PalomarQUEST project (Schweitzer et al. 2006), or using the DECam being built for the Dark Energy Survey (Honscheid et al. 2008). Since the challenge is to constrain long time scales, the presence of a multipleyear gap is mainly a complication for ensuring that problems in matching photometric bands are not interpreted as a form of longterm variability. Obtaining spectra of the Magellanic Cloud quasars is in some ways easier, because the quasar magnitudes and densities are wellsuited to the AAOmega fiber spectrograph on the AAT and, to a lesser degree, the IMACS spectrograph on Magellan.
The more general, Northern monitoring projects such as PanSTARRS and PalomarQUEST will slowly build light curves for essentially all the SDSS quasars, but at present their cadences are not ideal (see Schmidt et al. 2010) and it will take nearly a decade to build the long duration light curves needed for the analysis. In the long term, observations will be significantly improved with the advent of nextgeneration sky surveys. Most notably, the Large Synoptic Survey Telescope (LSST, Ivezić et al. 2008) will obtain accurate, wellsampled light curves for millions of AGN. The observed distribution of restframe characteristic time scales for S82 quasars spans the range from about 10 days to 1000 days (c.f. Figure 3). To probe the time scales as short as , and assuming a characteristic redshift of 2, the light curves should be sampled every 3 days in the observer’s frame, which is in good agreement with the baseline cadence of LSST. With a 10 yearlong survey, the length of the light curves will be in the range (1200). A combination of the SDSS, PanSTARRS, DES, and LSST data for 10,000 Stripe 82 quasars would span well over two decades, with multiband photometry obtained for hundreds of epochs, and would represent the best sample to date for studying the optical continuum variability of quasars. In particular, such a dataset would enable a robust measurement of the lowfrequency behavior of their PSD (c.f. point 5 above). For illustration, the LSST photometric errors in the band will be mag for , and there are roughly 23 million AGN with in the 20,000 sq. deg. covered by the main LSST survey (see Table 10.2 in the LSST Science Book; Abell et al. 2009). Each of these objects will be observed about 1000 times, yielding a database of over 2 billion photometric measurements. This data set, roughly a thousand times larger than that analyzed here, will enable a significant improvement in our understanding of quasar variability.
Appendix: Search for Periodic Light Curves
Although a stochastic process has proven to be an accurate statistical description of quasar light curves, any discovery of periodic behavior (even of a single source) would have interesting physical implications. However, a periodogram can only be used as an indicator of significant periodicity in a signal as compared to pure white noise (i.e., having no signal at all). Since quasars are genuinely variable, as described by a damped random walk (i.e., a red noise process; see KBS09), one can only evaluate the significance of the periodicity when also allowing for the signal covariance that is also present (see Markwardt et al. 2009; Gotz et al. 2009; Cenko et al. 2009, and references therein). Nevertheless, as a quick test for outstanding cases of periodicity among the S82 light curves, we analyzed 8,863 light curves for evidence of periodicity using the LombScargle periodogram (Lomb 1976; Scargle 1982). The threshold for considering the strongest peak in the periodogram as candidate evidence for periodicity was set following Horne & Baliunas (1986), with an adopted false alarm probability of 0.05. However, this is the threshold for ruling out white noise in favor of periodicity, and therefore provides no information on how a periodic description compares to a colored noise process such as a damped random walk. To determine the latter, a different threshold is needed (see Koz10). When adopting the threshold for ruling out white noise in favor of periodic variability, we identify 88 light curves as good candidates.
A close inspection of the period distribution for the full sample and the selected 88 candidates shows important differences: while the full sample displays a fairly flat distribution ranging from 100 days to values exceeding 10,000 days, the period distribution for 88 candidates is bimodal. The first peak with 22 objects corresponds to aliasing at roughly oneyear sampling cadence, while the second peak is centered on periods of about 67 years, similar to the total length of observations. We have visually inspected the light curves and phased light curves for all 88 candidates. It turns out that candidates with proposed periods of the order one year have light curves consistent with aliasing, while those with longer periods typically have only one observed “oscillation” that might not be used as robust evidence for periodicity (six examples are shown in Figure 18). Therefore, our search for periodic light curves in the S82 quasar sample has not yielded any convincing cases. Again, it is important to note that in cases such as in Figure 18, the periodogram indicates that these sources with long are likely to be periodic, not because they show true periodicity, but because it computes the likelihood relative to the wrong null hypothesis (white noise rather than colored noise). Moreover, the fact that 5% of the light curves exceed the threshold for periodicity being a better description than pure white noise is a further indication that the periodogram is not a very powerful statistic for poorly sampled light curves, since it is clear from the previous Sections that quasar light curves are not white noise, but rather are well described by a damped random walk. In fact, the single “oscillation” observed for the best candidate long period objects is entirely expected from colored noise processes such as a damped random walk (see right panels of Figure 18). Therefore, this analysis is further evidence that a damped random walk is a good description of quasar light curves.
Footnotes
 affiliation: cmacleod@astro.washington.edu
 affiliation: Department of Astronomy, University of Washington, Box 351580, Seattle, WA 98195
 affiliation: Department of Astronomy, University of Washington, Box 351580, Seattle, WA 98195
 affiliation: Department of Astronomy, The Ohio State University, 140 West 18th Avenue, Columbus, OH 43210
 affiliation: The Center for Cosmology and Astroparticle Physics, The Ohio State University, 191 West Woodruff Avenue, Columbus, OH 43210
 affiliation: Department of Astronomy, The Ohio State University, 140 West 18th Avenue, Columbus, OH 43210
 affiliation: Hubble Fellow
 affiliation: HarvardSmithsonian Center for Astrophysics, 60 Garden St, Cambridge, MA 02138
 affiliation: Department of Astronomy, University of Washington, Box 351580, Seattle, WA 98195
 affiliation: Department of Astronomy, University of Washington, Box 351580, Seattle, WA 98195
 affiliation: Department of Astronomy, University of Washington, Box 351580, Seattle, WA 98195
 affiliation: Department of Astronomy, University of Washington, Box 351580, Seattle, WA 98195
 affiliation: James Cook Univ., Centre for Astronomy, Townsville, QLD 4811, Australia
 affiliation: Department of Astronomy, University of Washington, Box 351580, Seattle, WA 98195
 affiliation: Department of Astronomy, University of Washington, Box 351580, Seattle, WA 98195
 affiliation: Department of Astronomy, University of Washington, Box 351580, Seattle, WA 98195
 affiliation: University of California, One Shields Ave, Davis, CA 95616
 The used here is related to the used in KBS09 () and the parameter used in Koz10 as .
 The functional form of SF fit to the longterm SDSSPOSS data in Sesar et al. (2006), see their Eq. 5, is similar but not identical to the functional form given by Eq. 3.
 A detailed analysis of the use of this variability model to select candidate quasars will be presented elsewhere (Brooks et al., in prep.).
 Based on Kcorrected magnitudes. Without the Kcorrection, this coefficent changes to .
 Measurement errors in of 0.2 dex have been included in the fitting. These coefficients are recommended when applying the model.
 Based on Kcorrected magnitudes. Without the Kcorrection, this coefficent changes to .
 Measurement errors in of 0.2 dex have been included in the fitting. These coefficients are recommended when applying the model.
 The morphological radio classes are defined using the integrated and peak FIRST fluxes as in Kimball & Ivezić (2008).
 The morphological radio classes are defined using the integrated and peak FIRST fluxes as in Kimball & Ivezić (2008).
References
 Abazajian, K. N., et al. 2009, ApJS, 182, 543
 Abell, P. A. et al. 2009, arXiv:0912.0201; http://www.lsst.org/lsst/scibook
 Ai, Y. L., Yuan, W., Zhou, H. Y., Wang, T. G., Dong, X. ., Wang, J. G., & Lu, H. L. 2010, arXiv:1005.0901
 Aretxaga, I., Cid Fernandes, R., & Terlevich, R. J. 1997, MNRAS, 286, 271
 Arévalo, P., Uttley, P., Kaspi, S., Breedt, E., Lira, P., & McHardy, I. M. 2008b, MNRAS, 389, 1479
 Arévalo, P., McHardy, I. M., & Summons, D. P. 2008a, MNRAS, 388, 211
 Bauer, A., Baltay, C., Coppi, P., Ellman, N., Jerke, J., Rabinowitz, D., & Scalzo, R. 2009, ApJ, 696, 1241
 Bhatti, W. A., Richmond, M. W., Ford, H. C., & Petro, L. D. 2010, ApJS, 186, 233
 Blackburne, J. A., & Kochanek, C. S. 2010, arXiv:1002.3126
 Cenko, S. B., et al. 2009, arXiv:0911.3150
 Churazov, E., Gilfanov, M., & Revnivtsev, M. 2001, MNRAS, 321, 759
 Collier, S., & Peterson, B. M. 2001, ApJ, 555, 775
 Croom, S.M., Richards, G.T., Shanks, T., et al. 2009, MNRAS, 399, 1755
 Denney, K. D., Peterson, B. M., Dietrich, M., Vestergaard, M., & Bentz, M. C. 2009, ApJ, 692, 246
 de Vries, W. H., Becker, R. H., White, R. L., & Loomis, C. 2005, AJ, 129, 615
 Emmanoulopoulos, D., McHardy, I. M., & Uttley, P. 2010, arXiv:1001.2045
 Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics (Cambridge, UK: Cambridge University Press)
 Fukugita, M., Ichikawa, T., Gunn, J.E., Doi, M., Shimasaku, K., & Schneider, D.P. 1996, AJ, 111, 1748
 Geha, M., et al. 2003, AJ, 125, 1
 Gotz, D., Mereghetti, S., von Kienlin, A., & Beck, M. 2009, GRB Coordinates Network, 9649, 1
 Gunn, J.E., et al. 1998, AJ, 116, 3040
 Giveon, U., Maoz, D., Kaspi, S., Netzer, H., & Smith, P. S. 1999, MNRAS, 306, 637
 Hawkins, M. R. S. 1993, Nature, 366, 242
 Hawkins, M. R. S. 2007, A&A, 462, 581
 Honscheid, K., DePoy, D. L., & for the DES Collaboration 2008, arXiv:0810.3600
 Hook, I. M., McMahon, R. G., Boyle, B. J., & Irwin, M. J. 1994, MNRAS, 268, 305
 Horne, J.H., & Baliunas, S.L. 1986, ApJ, 302, 757
 Hughes, P. A., Aller, H. D., & Aller, M. F. 1992, ApJ, 396, 469
 Ivezić, Ž., et al. 2002, AJ, 124, 2364
 Ivezić, Ž., et al. 2004, The Interplay Among Black Holes, Stars and ISM in Galactic Nuclei, 222, 525 [I04]
 Ivezić, Ž., et al. 2007, AJ, 134, 973
 Ivezić, Ž., Tyson, J. A., Allsman, R., Andrew, J., Angel, R., & for the LSST Collaboration 2008, arXiv:0805.2366
 Kaiser, N., et al. 2002, Proc. SPIE, 4836, 154
 Kawaguchi, T., Mineshige, S., Umemura, M., & Turner, E. L. 1998, ApJ, 504, 671
 Kelly, B. C. 2007, ApJ, 665, 1489
 Kelly, B. C., Bechtold, J., & Siemiginowska, A. 2009, ApJ, 698, 895 [KBS09]
 Kimball, A. E., & Ivezić, Ž. 2008, AJ, 136, 684
 Kollmeier, J. A., et al. 2006, ApJ, 648, 128
 Kozłowski, S., & Kochanek, C. S. 2009, ApJ, 701, 508
 Kozłowski, S., et al. 2010a, ApJ, 708, 927 [Koz10]
 Kozłowski, S., et al. 2010b, arXiv:1002.3365
 Li, S.L., & Cao, X. 2008, MNRAS, 387, L41
 Liu, H. T., Bai, J. M., Zhao, X. H., & Ma, L. 2008, ApJ, 677, 884
 Lomb, N.R. 1976, Ap&SS, 39, 447
 MacLeod, C., Ivezić, Ž., de Vries, W., Sesar, B., & Becker, A. 2008, American Institute of Physics Conference Series, 1082, 282
 Marconi, A., Axon, D. J., Maiolino, R., Nagao, T., Pastorini, G., Pietrini, P., Robinson, A., & Torricelli, G. 2008, ApJ, 678, 693
 Markwardt, C. B., Gavriil, F. P., Palmer, D. M., Baumgartner, W. H., & Barthelmy, S. D. 2009, GRB Coordinates Network, 9645, 1
 Martini, P., & Schneider, D. P. 2003, ApJ, 597, L109
 Matthews, T. A., & Sandage, A. R. 1963, ApJ, 138, 30
 McHardy, I. 2010, Lecture Notes in Physics, Berlin Springer Verlag, 794, 203
 Pereyra, N. A., Vanden Berk, D. E., Turnshek, D. A., Hillier, D. J., Wilhite, B. C., Kron, R. G., Schneider, D. P., & Brinkmann, J. 2006, ApJ, 642, 87
 Peterson, B. M., et al. 2005, ApJ, 632, 799
 Press, W. H., Rybicki, G. B., & Hewitt, J. N. 1992, ApJ, 385, 404
 Rengstorf, A. W., Brunner, R. J., & Wilhite, B. C. 2006, AJ, 131, 1923
 Rybicki, G. B., & Press, W. H. 1992, ApJ, 398, 169
 Rybicki, G. B., & Press, W. H. 1994, Computer, 5004
 Scargle, J.D. 1982, ApJ, 263, 835
 Schmidt, K. B., Marshall, P. J., Rix, H.W., Jester, S., Hennawi, J. F., & Dobler, G. 2010, arXiv:1002.2642
 Schneider, D. P., et al. 2007, AJ, 134, 102
 Schweitzer, M., et al. 2006, ApJ, 649, 79
 Sesar, B., et al. 2006, AJ, 131, 2801
 Sesar, B., et al. 2007, AJ, 134, 2236
 Shen, Y., Greene, J. E., Strauss, M. A., Richards, G. T., & Schneider, D. P. 2008, ApJ, 680, 169
 Smith, J. A., et al. 2002, AJ, 123, 2121
 Timmer, J., & Koenig, M. 1995, A&A, 300, 707
 Trèvese, D., Kron, R. G., & Bunone, A. 2001, ApJ, 551, 103
 Trèvese, D., & Vagnetti, F. 2002, ApJ, 564, 624
 Udalski et al. 1997, Acta Astron 47 319
 Udalski et al. 2008, Acta Astron 58 69
 Uttley, P., McHardy, I. M., & Papadakis, I. E. 2002, MNRAS, 332, 231
 Vanden Berk, D. E., et al. 2004, ApJ, 601, 692 [VB04]
 Vestergaard, M., & Peterson, B. M. 2006, ApJ, 641, 689
 Voges, W., et al. 1999, A&A, 349, 389
 White, R. L., Becker, R. H., Helfand, D. J., & Gregg, M. D. 1997, ApJ, 475, 479
 Wilhite, B. C., Vanden Berk, D. E., Kron, R. G., Schneider, D. P., Pereyra, N., Brunner, R. J., Richards, G. T., & Brinkmann, J. V. 2005, ApJ, 633, 638
 Wilhite, B. C., Vanden Berk, D. E., Brunner, R. J., & Brinkmann, J. V. 2006, ApJ, 641, 78
 Wilhite, B. C., Brunner, R. J., Grier, C. J., Schneider, D. P., & Vanden Berk, D. E. 2008, MNRAS, 383, 1232
 Wold, M., Brotherton, M. S., & Shang, Z. 2007, MNRAS, 375, 989
 York, D.G. et al. 2000, AJ, 120, 1579