How Far Is Quasar UV/Optical Variability from a Damped Random Walk at Low Frequency?

SDSS Stripe 82, lying along the celestial equator in the southern Galactic hemisphere, covers an area of ∼290 square degrees, and has been scanned ∼60 times in the ugriz bands by the SDSS imaging survey. MacLeod et al. (2012) presented recalibrated ∼10 year long five-band light curves for 9275 spectroscopically confirmed quasars in the field.⁶ Measurements of black hole mass, absolute magnitude (K-corrected) and bolometric luminosity for most of them are available, taken from Shen et al. (2008). With these long light curves of such a large sample of quasars, we can extensively study the scatter of quasar variability parameters.

We model the SDSS light curves using the DRW process. The DRW as a stochastic process is well described by the exponential covariance matrix of the signal, with its structure function expressed as (e.g., Hughes et al. 1992; MacLeod et al. 2010)

$$\begin{eqnarray}&&\mathrm{SF}({\rm{\Delta }}t)={\mathrm{SF}}_{\infty }{(1-{e}^{-| {\rm{\Delta }}t| /\tau })}^{1/2},\end{eqnarray} \tag{ 1 }$$

with ${\mathrm{SF}}_{\infty }=\sqrt{2}\sigma$, τ the characteristic timescale, and ${\rm{\Delta }}t$ the time lag. Clearly,

$$\begin{eqnarray}&&\mathrm{SF}({\rm{\Delta }}t\gg \tau )={\mathrm{SF}}_{\infty }=\sqrt{2}\sigma ,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\mathrm{SF}({\rm{\Delta }}t\ll \tau )=\sigma \sqrt{\displaystyle \frac{2| {\rm{\Delta }}t| }{\tau }}={\mathrm{SF}}_{\infty }\sqrt{\displaystyle \frac{| {\rm{\Delta }}t| }{\tau }}=\hat{\sigma }\sqrt{{\rm{\Delta }}t}.\end{eqnarray} \tag{ 3 }$$

In reality, measuring τ and ${\mathrm{SF}}_{\infty }$ is often difficult due to the limited light curve length. Instead, $\hat{\sigma }$ (=${\mathrm{SF}}_{\infty }/\sqrt{\tau }$) can be measured with much better accuracy, even when τ and ${\mathrm{SF}}_{\infty }$ are poorly constrained, and can be used as a good proxy to evaluate the properties of quasar variabilities. According to the continuous-time first-order autoregressive process (CAR(1), Kelly et al. 2009), the PSD of a DRW process is given by

$$\begin{eqnarray}&&P(f)=\displaystyle \frac{4{\sigma }^{2}\tau }{1+{(2\pi \tau f)}^{2}},\end{eqnarray} \tag{ 4 }$$

which infers that $P(f)\propto {f}^{-2}$ when $f\gt {(2\pi \tau )}^{-1}$ and P(f) yields a constant when $f\lt {(2\pi \tau )}^{-1}$. In our paper, the nomenclature "DRW" only refers to the CAR(1) process, and "non-DRW" models refer to several special non-CAR(1) models, which are included in a more general model (continuous-time autoregressive moving average, CARMA; Kelly et al. 2014).

In this work, to make a direct comparison with the results of MacLeod et al. (2012), we elect to use the same r-band light curves, which have the best photometries among the five bands for the quasars. We also adopt similar sample selections: (1) only quasars with ≥40 epochs (to ensure reasonable DRW fitting) and measurements of black hole mass and absolute magnitude are retained; (2) ${\rm{\Delta }}{L}_{\mathrm{noise}}=\mathrm{ln}({L}_{\mathrm{best}}/{L}_{\mathrm{noise}})\gt 2$ and ${\rm{\Delta }}{L}_{\infty }=\mathrm{ln}({L}_{\mathrm{best}}/{L}_{\infty })\gt 0.05$, where ${\rm{\Delta }}{L}_{\mathrm{noise}}$, ${\rm{\Delta }}{L}_{\infty }$, ${\rm{\Delta }}{L}_{\mathrm{best}}$ are the likelihoods for τ = 0, $\tau =\infty$, and the best-fit τ (Kozłowski et al. 2010; Zu et al. 2013). The former equation is equivalent to setting a noise limit to request the obtained characteristic timescale larger than the average cadence. The latter eliminates the fitting results that cannot distinguish the best τ and $\tau =\infty$ based on the insufficient light curve length.

We fit SDSS r-band light curves with the Javelin code⁷ (Zu et al. 2011) to measure the damping timescale τ and asymptotic magnitude ${\mathrm{SF}}_{\infty }$ ⁸ (see the upper panel of Figure 1 for a demonstration). According to the simulations of Kozłowski (2017a), only if the characteristic timescale (in the observed frame, the same hereafter unless otherwise stated) is less than ∼10% of the light curve length, the input DRW parameters ${\tau }_{\mathrm{in}}$ can be reasonably recovered by DRW fitting, although biased, depending on different priors and fitting method. To address this issue, we only include quasars at lower redshift ($z\lt 1.2$) and with low black hole mass (Log ${M}_{\mathrm{BH}}\lt 9$, see Figure 2) for which the damping timescale is considerably small compared with the light curve length (e.g., Kozłowski 2016b). The final sample includes 1678 quasars. The median value of their observed ${\tau }_{\mathrm{obs}}$ is 301 days, and the typical length of the Stripe 82 light curves is ∼2800 days. The ratio of the typical damping timescale to the light curve duration is 11%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 3 we plot the distributions of the measured ${\tau }_{\mathrm{obs}}$, ${\mathrm{SF}}_{\infty ,\mathrm{ratio}}$, ${\hat{\sigma }}_{\mathrm{ratio}}$, and K_ratio for our quasar sample, where $K=\tau \sqrt{{\mathrm{SF}}_{\infty }}$ is a variable orthogonal to $\hat{\sigma }$ in log space. All of them show clear scatters, but slightly smaller than those of MacLeod et al. (2010) as we only include quasars with small black hole mass and redshift, for which τ, ${\mathrm{SF}}_{\infty }$ can be more reliably constrained.

Zoom In Zoom Out Reset image size

We revisited the empirical relations given by MacLeod et al. (2010) between the measured DRW parameters (τ, ${\mathrm{SF}}_{\infty }$ and $\hat{\sigma }$), and physical parameters including the i-band absolute magnitude (${M}_{i}$), rest frame wavelength (${\lambda }_{\mathrm{RF}}$), black hole mass (${M}_{\mathrm{BH}}$), and redshift (z). Following their method, we first determined the rest wavelength dependence of the variability parameters (τ, ${\mathrm{SF}}_{\infty }$, and $\hat{\sigma }$) using various SDSS bands, since each individual quasar has other physical parameters fixed. Note we only utilize the $u,g,r$ bands here, and exclude the $i,z$ bands since the weaker intrinsic variabilities and longer intrinsic damping timescales of the latter make the measurements of τ and ${\mathrm{SF}}_{\infty }$ less reliable than in the bluer bands. We then use power laws to fit other physical parameters simultaneously with fixed wavelength coefficients for r band:

$$\begin{eqnarray}{\rm{log}}\,f & = & A+B{\rm{log}}\left(\displaystyle \frac{{\lambda }_{{\rm{RF}}}}{4000\mathring{{\rm{A}}}}\right)+C({M}_{i}+23)\\ & & +D{\rm{log}}\left(\displaystyle \frac{{M}_{{\rm{BH}}}}{{10}^{9}{M}_{\odot }}\right)+E{\rm{log}}(1+z)\end{eqnarray} \tag{ 5 }$$

Using our redefined quasar sample, we obtain A = −0.50, B = −0.48, C = 0.13, and D = 0.18 for $f={\mathrm{SF}}_{\infty }$ (mag), A = 2.1, B = 0.66, C = −0.19, and D = 0.14 for f = τ (days, in the rest frame), and A = −1.72, B = −0.86, C = 0.22, and D = 0.12 for $f=\hat{\sigma }$ ($\mathrm{mag}/{\mathrm{day}}^{1/2}$, in the rest frame). All E are fixed to 0, assuming the variability has no evolution with redshift. The typical errors of these coefficients are 0.05. The correlations imply that not only the long/short-term variability but also the damping timescale strongly depend on the rest frame wavelength, luminosity, and black hole mass. Note in Equation (5) that the dependence on ${M}_{i}$ and ${M}_{\mathrm{BH}}$ can be coupled as both quantities correlate with redshift. Caution is thus needed when comparing the coefficients with those reported in other studies (e.g., Kozłowski 2016b, 2017b). Obtaining decoupled coefficients requires samples spanning larger ranges in luminosity and ${M}_{\mathrm{BH}}$, and is beyond the scope of this work. For this study the coefficients we obtained using our sample are sufficient to derive the residual scatters in the DRW parameters (excluding the dependence on the physical parameters listed in Equation (5), see next paragraph).

We then obtain the expected ${\tau }_{\exp }$, ${\mathrm{SF}}_{\infty ,\exp }$, ${\hat{\sigma }}_{\exp }$, and K_exp for each quasar based on the empirical relations. In Figure 3, we plot the distribution of ${\tau }_{\mathrm{obs}}/{\tau }_{\exp }$, ${\mathrm{SF}}_{\infty ,\mathrm{obs}}/{\mathrm{SF}}_{\infty ,\exp }$, ${\hat{\sigma }}_{\mathrm{obs}}/{\hat{\sigma }}_{\exp }$, and K_obs/K_exp respectively, where we can see that the residual scatters (the blue dotted line refers to the scatter after excluding the dependence on the physical parameters listed in Equation (5)), in all four quantities are still prominent. Note that the uncertainties in the measurement of black hole mass contribute little to these residual scatters. For instance, based on Equation (5), a 0.3 dex uncertainty in BH mass can produce a scatter of 0.036 in log($\hat{\sigma }$), which is only 7% of the residual scatter of $\hat{\sigma }$ (${0.036}^{2}/{0.14}^{2}$).

The residual scatters can be at least partially attributed to the uncertainties in the measurement of the DRW parameters, due to limited light curve length, sparse time sampling, and photometric errors. Following MacLeod et al. (2010, 2012), we perform Monte Carlo simulations to quantify such uncertainties. Using the observed ${\tau }_{\mathrm{obs}}$ and ${{\rm{SF}}}_{{\rm{\infty }},{\rm{obs}}}$, we generate a 100 year long artificial DRW light curve spaced every 1 day for each quasar after eliminating the burn-in part (Kelly et al. 2009). A 10 year segment is randomly cut to cross-match with the observed light curve, with the Stripe 82 time sampling and photometric errors imposed at the same time (see the lower panel of Figure 1). We generate one artificial light curve for each quasar in our sample, then fit the artificial light curves again to obtain the output DRW parameter. The ratios of the output to input DRW parameters are over-plotted in Figure 3, from which we can see that the uncertainties in the measurement of the DRW parameters can only account for 41%, 68%, 33%, and 44% of the residual scatter (${\sigma }_{\mathrm{red}}^{2}/{\sigma }_{\mathrm{blue}}^{2}$), for τ, ${\mathrm{SF}}_{\infty }$, $\hat{\sigma }$, and K, respectively. Below we focus on the $\hat{\sigma }$ for which DRW fitting to the simulated light curves recovers the input values with the smallest scatter (red line in the lower left panel of Figure 3), and the largest fraction of the residual scatter (67%) remains unaccounted for.

To better demonstrate how the input the DRW parameters were recovered through fitting the artificial light curves, in Figure 4 we plot the contour distributions of output τ and ${\mathrm{SF}}_{\infty }$ for four different pairs of input values. The input values of τ are 100, 300, 1000, and 3000 days, and the input ${\mathrm{SF}}_{\infty }$ was adjusted to keep $\hat{\sigma }$ unchanged. For each pair of input τ and ${\mathrm{SF}}_{\infty }$, we generate ∼2000 artificial light curves and measure output values for each of them. The contour plots show that the output τ and ${\mathrm{SF}}_{\infty }$ (for each pair of input values) are clearly coupled, following the direction of constant $\hat{\sigma }$. This indicates that the input $\hat{\sigma }$ can be recovered with remarkably high accuracy and small scatter, comparing with τ and ${\mathrm{SF}}_{\infty }$ (see also the red lines in Figure 3). This holds even for those light curves with intrinsic large τ (∼1000 days), for which the ${\tau }_{\mathrm{out}}$ and ${\mathrm{SF}}_{\infty ,\mathrm{out}}$ are poorly constrained with huge scatters. This is also one of the key reasons that we opted to focus on $\hat{\sigma }$ (but not ${\tau }_{\mathrm{out}}$ or ${\mathrm{SF}}_{\infty }$) in this work. Offsets of the median output values from input ones are seen for τ, ${\mathrm{SF}}_{\infty }$ and $\hat{\sigma }$. Such offsets are dependent on the sampling and priors adopted in the DRW fitting. Correcting such offsets will not change the results we present in this work.

Zoom In Zoom Out Reset image size

The intrinsic τ for some of the quasars in our sample may still have been significantly underestimated due to the limited length of the light curves. We perform further simulations to check whether such effect may alter our results. We enlarge the input τ value for each source in our sample by a factor of 2–10, and adjust ${\mathrm{SF}}_{\infty }$ accordingly to keep input $\hat{\sigma }$ unchanged. However, we see no significant difference (<5%) in the scatter of ${\hat{\sigma }}_{\mathrm{out}}/{\hat{\sigma }}_{\mathrm{in}}$, demonstrating that our results are not affected by possible underestimations of τ in our sample. It can also be seen in Figure 4 that, although larger input τ yield much stronger scatter in ${\tau }_{\mathrm{out}}$, the scatter of ${\hat{\sigma }}_{\mathrm{out}}$ appears insensitive to input τ. We also double-check this issue by dividing our sample into two equal-size subsamples according to the BH mass. We find that adopting only the lower-mass sample ($\mathrm{Log}\ {M}_{\mathrm{BH}}\lesssim 8.5$), for which the intrinsic τ should be even smaller, does not alter the results presented above.

Non-DRW light curves are simulated from a broken power-law PSD with the public Python code pyLCSIM,⁹ based on the algorithm of Timmer & Koenig (1995). Two key variables are needed in this approach: one is the break frequency ${(2\pi \tau )}^{-1}$, where τ is the damping timescale; the other is the rms variability amplitude ($\mathrm{rms}=\sigma ={\mathrm{SF}}_{\infty }/\sqrt{2}$ for DRW). As a starting point, we use the observed ${(2\pi {\tau }_{\mathrm{obs}})}^{-1}$ and ${\mathrm{SF}}_{\infty ,\mathrm{obs}}/\sqrt{2}$ (where ${\tau }_{\mathrm{obs}}$ and ${\mathrm{SF}}_{\infty ,\mathrm{obs}}$ are the best-fit DRW parameters) of each quasar as input values to generate 100 year long artificial light curves spaced every 1 day. The slopes of the broken power-law model can be arbitrarily decided. We investigate four non-DRW models comparing with the standard DRW model: setting low-frequency slopes ${\alpha }_{l}=-1$ (model A), −1.9 (model B) with fixed high-frequency slopes at ${\alpha }_{h}=-2$; and ${\alpha }_{h}=-1.8$ (model C), −3 (model D) with fixed ${\alpha }_{l}=0$ (see Figure 5).

Zoom In Zoom Out Reset image size

Figure 6 briefly illustrates the deviation of non-DRW models from the DRW one. The first three panels show three 100 year light curves produced by three models with the same input values (rms = 0.1 mag, τ = 10 days ), with ${\alpha }_{h}$ fixed at −2 and slopes ${\alpha }_{l}=0$ (DRW), −1 (model A), and −1.9 (model B), respectively. In the last panel, the corresponding power density spectra are plotted. Note the pyLCSIM code generates light curves with rms (measured through the whole light curve) simply fixed at the input value. Short segments of such light curves have smaller rms for red noise PSD, but with dispersion in rms, reflecting red noise leakage. We find that selecting a 10 year segment out of a 100 year long simulated light curve is sufficient to reproduce the effect of red noise leakage (for the PSD we chose), and simulating even longer light curves is not needed. For the DRW process, the variation is white noise at timescales longer than τ, thus there is no red noise leakage. On the other hand, utilizing the DRW model to fit the red noise PSD will misestimate the τ and ${\mathrm{SF}}_{\infty }$, as shown in the last panel of Figure 6 (for the case of ${\alpha }_{l}=-1$).

Zoom In Zoom Out Reset image size

In Figure 7 we demonstrate the discrepancies of the output τ and ${\mathrm{SF}}_{\infty }$ from the input ones used to simulate artificial light curves based on various models. For each set of input parameters (shown in Figure 7 with fixed τ = 300 days and ${\mathrm{SF}}_{\infty }$ changing from 0.01 to 0.64 mag, or fixed ${\mathrm{SF}}_{\infty }=0.18$ mag and τ changing from 50 to 3200 days), we generate 10,000 artificial resampled light curves and plot the averaged output parameters versus the input ones. From Figure 7 we see that for the DRW model, the input parameters are well recovered (with only slight offsets), except that ${\tau }_{\mathrm{out}}$ is significantly underestimated at ${\tau }_{\mathrm{in}}$ > 1000 days (due to the limited length of the light curves), consistent with the results of Kozłowski (2017a). For the non-DRW models, only models C and D yield ${\mathrm{SF}}_{\infty ,\mathrm{out}}$ consistent with input values. Significant deviations from input τ are seen for all non-DRW models. For model B, particularly, the output τ barely correlates with the input value. This is because the PSD break in model B is so weak that modeling it with the DRW process somehow makes no sense (imaging modeling pure red noise PSD slope of −2 with a DRW process).

Zoom In Zoom Out Reset image size

To correct the aforementioned biases (mismatches between input and output τ and ${\mathrm{SF}}_{\infty }$), we perform simulations for a grid of input parameters (τ ranging from 1 to 10⁵ days and ${\mathrm{SF}}_{\infty }$ from 0.01 to 0.8 $\mathrm{mag}$). The output parameters (averaged over a large number of simulations with the same input parameters) build a ${\tau }_{\mathrm{out}}\mbox{--}{\mathrm{SF}}_{\infty ,\mathrm{out}}$ surface. Finding the position of each quasar on the ${\tau }_{\mathrm{out}}\mbox{--}{\mathrm{SF}}_{\infty ,\mathrm{out}}$ plane, we derive the corresponding ${\tau }_{\mathrm{in}}$ and ${\mathrm{SF}}_{\infty ,\mathrm{in}}$ to be input for simulations.

In this section we examine whether non-DRW models could account for the residual scatter in $\hat{\sigma }$. Comparing with the DRW, models A and B (deviation from the DRW at low frequency: ${\alpha }_{l}=-1$, −1.9 instead of 0, with ${\alpha }_{h}$ fixed at −2) produce obviously larger scatters in ${\hat{\sigma }}_{\mathrm{out}}/{\hat{\sigma }}_{\mathrm{in}}$ (Figure 8). This is because of red noise leakage which refers to the variability power transferred from low to high frequencies (Uttley et al. 2002; Zhu & Xue 2016). This effect is more significant for PSDs with steeper slopes at low frequencies. An offset from unity of the ${\hat{\sigma }}_{\mathrm{out}}$/${\hat{\sigma }}_{\mathrm{in}}$ distribution is seen for model B. As explained in Section 3.1, this is because the PSD break in model B is so weak that modeling it with the DRW yields an output τ insensitive to the input value, thus our corrections in Section 3.1 are not sufficient. Nevertheless, this does not affect the major conclusions of this work.

Zoom In Zoom Out Reset image size

We repeat our simulations for different low-frequency slopes. In Figure 9, we plot the reproduced scatter in ${\hat{\sigma }}_{\mathrm{out}}$/${\hat{\sigma }}_{\mathrm{in}}$ versus the low-frequency PSD slope, where we see a clear trend that a steeper lower-frequency PSD yields a stronger scatter in $\hat{\sigma }$. Simulating and modeling quasar light curves with the DRW provides a much lower scatter in $\hat{\sigma }$ ($0.08\ \mathrm{mag}/{\mathrm{yr}}^{1/2}$) as compared to empirically measured values (0.14 $\mathrm{mag}/{\mathrm{yr}}^{1/2}$). A non-DRW model with low-frequency slope ${\alpha }_{l}\simeq -1.3$ appears just sufficient to explain the observed residual scatter in $\hat{\sigma }$.

Zoom In Zoom Out Reset image size

When modeling the SDSS stripe 82 light curves, we obtained a runaway fraction (RF) of ${\mathrm{RF}}_{\mathrm{obs}}=19$%, i.e., for 19% of the light curves τ cannot be well constrained (output τ hits the boundary at τ = 10⁵ days). The RF was also measured for our simulations: ${11}_{-2}^{+2}$% for the DRW model, ${16}_{-2}^{+2}$% for model A, and ${40}_{-3}^{+3}$% for model B, ${6}_{-1}^{+1}$% for model C and ${28}_{-2}^{+2}$% for model D. The errors in RF are measured as scatters from different sets of simulations. The observed RF appears more consistent with model A, also suggesting the low-frequency PSD slope is likely ${\alpha }_{l}\sim -1$, deviating from the DRW.

For model C and D (${\alpha }_{h}=-1.8$, −3 instead of −2, with ${\alpha }_{l}$ fixed at 0), as there is no red noise leakage; both models yield ${\hat{\sigma }}_{\mathrm{out}}/{\hat{\sigma }}_{\mathrm{in}}$ scatter similar to that of the DRW model (see Figure 10).

Zoom In Zoom Out Reset image size