A degeneracy in DRW modelling of AGN light curves

Abstract

Individual light curves of active galactic nuclei (AGNs) are nowadays successfully modelled with the damped random walk (DRW) stochastic process, characterized by the power exponential covariance matrix of the signal, with the power β = 1. By Monte Carlo simulation means, we generate mock AGN light curves described by non-DRW stochastic processes (0.5 ≤ β ≤ 1.5 and β ≠ 1) and show they can be successfully and well modelled as a single DRW process, obtaining comparable goodness of fits. A good DRW fit, in fact, may not mean that DRW is the true underlying process leading to variability and it cannot be used as a proof for it. When comparing the input (non-DRW) and measured (DRW) process parameters, the recovered time-scale (amplitude) increases (decreases) with the increasing input β. In practice, this means that the recovered DRW parameters may lead to biased (or even non-existing) correlations of the variability and physical parameters of AGNs if the true AGN variability is caused by non-DRW stochastic processes. The proper way of identifying the processes leading to variability are model-independent structure functions and/or power spectral densities and then using such information on the covariance matrix of the signal in light-curve modelling.

accretion, accretion discs, methods: data analysis, galaxies: active, quasars: general

1 INTRODUCTION

The damped random walk (DRW) model is an increasingly successful method of quantifying the variability of active galactic nuclei (AGNs; Kelly, Bechtold & Siemiginowska 2009; Kozłowski et al. 2010; MacLeod et al. 2010; Zu, Kochanek & Peterson 2011; Zu et al. 2013). Kelly et al. (2009) introduce DRW as an underlying stochastic process leading to AGN variability, also known as the first order continuous-time autoregressive process [CAR(1)] or Ornstein–Uhlenbeck process (Uhlenbeck & Ornstein 1930). The model has two parameters, the time-scale τ after which the light curve becomes uncorrelated and the modified amplitude |$\hat{\sigma }$| (Kozłowski et al. 2010), or asymptotic amplitude SF_∞ (MacLeod et al. 2010). These two parameters show correlations with the physical parameters of AGNs, such as the black hole mass, luminosity, Eddington ratio, and rest-frame wavelength. For example, Kelly et al. (2009) report that the time-scale τ is correlated with the black hole mass and luminosity, while the amplitude is anticorrelated with these parameters. MacLeod et al. (2010) study ∼9000 AGN from Stripe 82 of the Sloan Digital Sky Survey (SDSS) and find that the time-scale τ is correlated with the rest-frame wavelength and the black hole mass, and does not depend on redshift or luminosity. The asymptotic variability SF_∞ is anticorrelated with the luminosity, rest-frame wavelength, and the Eddington ratio. Kozłowski (2016) reanalysed the same set of SDSS AGN light curves with the ‘sub-ensemble’ structure function (SF) analysis that is model-independent and essentially confirms these correlations, albeit with a minute differences in these relations. He noticed, however, that the SF power-law slope γ steepens from β ≡ 2γ ≈ 1 for the fainter AGNs to about β ≈ 1.2 for the brightest AGNs and is independent of the black hole mass. Such a change means a departure from DRW that is paralleled by the DRW time-scale increase obtained from light-curve modelling (but a bulk of this is the true correlation with the black hole mass).

Can a non-DRW stochastic process be successfully and well modelled as DRW, and return correct variability parameters? Or will it rather return biased parameters, for example, longer time-scales for steeper SFs as in the SF analysis from Kozłowski (2016)? If the latter is the case, then it may have profound implications for the reported correlations of variability with the physical parameters of AGNs. In this paper, we are interested in the modelling of simulated AGN light curves as the DRW process, that are caused by other than DRW underlying processes, in order to find answers to the above questions.

In Section 2, we present the methodology of simulations and modelling of the quasar light curves, while in Section 3, we discuss our findings. The paper is concluded in Section 4.

2 METHODOLOGY

We simulate AGN light curves as a single stochastic process with the power exponential covariance matrix of the signal

\begin{equation} {\rm cov}(\Delta t) = \sigma _{\rm s}^2 e^{-\left(\frac{|\Delta t|}{\tau }\right)^\beta }, \end{equation}

(1)

where τ > 0 is the decorrelation time-scale, |$\sigma _{\rm s}^2$| is the signal variance, Δt = t_i − t_j is the rest-frame time difference (or time lag) between ith and jth data points, and 0 < β < 2, where β = 1 corresponds to DRW. To simulate a light curve with N points, first, the (N × N) covariance matrix of the signal

\begin{equation} C_{ij} = {\left(\begin{array}{c@{\,\,\,}c@{\,\,}c@{\,\,\,}c}\sigma _{\rm s}^2 & \sigma _{\rm s}^2 e^{-\left(\frac{|t_1-t_2|}{\tau }\right)^\beta } & \cdots & \sigma _{\rm s}^2 e^{-\left(\frac{|t_1-t_N|}{\tau }\right)^\beta } \\ \sigma _{\rm s}^2 e^{-\left(\frac{|t_2-t_1|}{\tau }\right)^\beta } & \sigma _{\rm s}^2& \cdots & \sigma _{\rm s}^2 e^{-\left(\frac{|t_2-t_N|}{\tau }\right)^\beta } \\ \vdots & \vdots & \ddots & \vdots \\ \sigma _{\rm s}^2 e^{-\left(\frac{|t_N-t_1|}{\tau }\right)^\beta } & \sigma _{\rm s}^2 e^{-\left(\frac{|t_N-t_2|}{\tau }\right)^\beta } & \cdots & \sigma _{\rm s}^2 \end{array}\right)} \end{equation}

(2)

must be Cholesky decomposed as |$\boldsymbol {C} = \boldsymbol {L}^{\rm T}\boldsymbol {L}$| (e.g. Press et al. 1992), where |$\boldsymbol {L}$| is the upper triangular matrix. Next, the light curve is obtained from |$\boldsymbol {y}=\boldsymbol {Lr}$|⁠, where |$\boldsymbol {r}$| is a vector of Gaussian deviations with the variance of unity (e.g. Zu et al. 2011). Finally, we add the photometric noise n_i(y_i) dependent on the magnitude y_i, y_i = y_i + G(n_i(y_i)), that is drawn from a Gaussian (G) distribution of the true SDSS photometric noise (Ivezić et al. 2007).

We simulate sets of a 1000 light curves in a range 0.5 ≤ β ≤ 1.5, spaced every 0.05. With the exception of the case with β = 1, they are regarded as the non-DRW stochastic processes. The light curves have the mean magnitude 〈r〉 = 17 mag and the noise properties of the SDSS Stripe82 quasars (Ivezić et al. 2007). The input time-scale is τ = 500 d, the asymptotic variability amplitude is SF_∞ = 0.18 mag, the light-curve length is 8000 d, and the cadence is 20 d (hence 400 points). Exemplary simulated light curves for β = 0.5, 1.0, and 1.5, are shown in Fig. 1.

Figure 1.

Typical simulated AGN light curves with 〈r〉 = 17 mag, τ = 500 d, SF_∞ = 0.18 mag, 400 data points and the length of 8000 d, with β = 1.5 (top), β = 1.0 (DRW; middle), and β = 0.5 (bottom). The top and bottom light curves are shifted by ±0.6 mag for clarity.

Subsequently, the light curves are modelled with DRW [see appendix in Kozłowski et al. 2010 for fast (only O(N) operations for a light curve with N points) modelling with DRW, also Kelly et al. 2009; MacLeod et al. 2010; Zu et al. 2011, 2013], i.e. with fixed β = 1.

3 DISCUSSION

In Fig. 2, we present the results of modelling AGN light curves as the DRW process, for which the underlying processes were set to be, with the exception of β = 1, non-DRW. For each 0.5 ≤ β ≤ 1.5, spaced every 0.05, we model a 1000 light curves and calculate the median-measured parameters along with the dispersions, measured as 0.74 of the interquartile range of these distributions (see MacLeod et al. 2012; Kozłowski 2016). Because the recovered parameters are also a function of the ratio of the time-scale τ to the experiment length (Kozłowski in preparation), we normalize the returned parameters to be unity for β = 1 (τ is divided by 0.86 and |$\hat{\sigma }$| is divided by 0.96).

Figure 2.

Ratios of the median recovered to input parameters (top and middle panels) and χ²/dof (bottom panel) as a function of the input power β. Each point corresponds to a 1000 simulated light curves with τ = 500 d, SF_∞ = 0.18 mag, 400 data points and the length of 8000 d, that are modelled as DRW (β = 1). Light curves with input β < 1 (β > 1) modelled as DRW, have underestimated time-scales and overestimated amplitudes (overestimated time-scales and underestimated amplitudes) while the goodness of fit weakly improves with increasing β. The error bars are in fact dispersions of the measured values calculated as 0.74 interquartile range (IQR) of the recovered parameter distributions.

In top panel of Fig. 2, we show the ratio of the median of the measured time-scales τ to the input value (τ = 500 d) as a function of the input parameter β. We see that the returned time-scale is correlated with β, but also the goodness of fit improves with increasing β (bottom panel of Fig. 2). Because in Kozłowski (2016), sub-ensemble SFs steepen from β ≈ 1.0 for fainter AGNs to about β ≈ 1.2 for the brightest ones, it means that DRW should return longer time-scales for the latter sources, even if the true time-scales were identical. The middle panel of Fig. 2 presents the dependence of the recovered modified amplitudes as a function of β. It is obvious that the two parameters are anticorrelated.

Finding the exact (parametric) form of these biases is not the goal of this paper, because they will depend on the photometric quality and length of data being analysed. The goal here is simply to make the point about their existence, their possible implications on (mis)understanding of AGN physics, and to provide a solution to avoid them. Because in Kozłowski (2016), β ≈ 1 and weakly changes (to ∼1.2 for the brightest AGNs), the positive side is that modest deviations from the DRW model seem to be nearly unimportant for the estimated variability parameters and they weakly affect the correlations with the physical AGN parameters in MacLeod et al. (2010). The increase of the input (or true) β from 1.0 to 1.2 leads to the overestimation of the time-scale τ by a factor of 1.5 (Fig. 2). The negative side is that typical AGN light curves are not good enough to notice the deviations from DRW and so one may misinterpret parameters. Because SFs or power spectral densities are a model-independent means of estimating the shape of the covariance function of the signal (e.g. Kozłowski 2016), one should rather estimate β this way, and then use it as input parameter in direct light-curve modelling to obtain correct model parameters. Zu et al. (2013) discusses modelling of light curves with additional parameters to that from DRW.

4 CONCLUSIONS

In this paper, we have been interested if AGN variability caused by non-DRW stochastic processes can be well modelled with a single DRW process, nowadays frequently considered in AGN variability studies. By simulation means, we have tested the implications of modelling non-DRW processes on the measured DRW model parameters and found that they are biased, where the time-scale increases and the amplitude decreases with the increasing input parameter β (the power of the power exponential covariance matrix of the signal). Equally important finding here is the goodness of fit being unable to recognize what process is being modelled, hence, the word ‘degeneracy’ in the title. A good DRW fit should not and cannot be used as a proof for DRW as the true underlying process leading to variability. Instead, the covariance matrix of the signal should be obtained from model-independent measures of variability such as the SFs or power spectral densities, and serve as input for the covariance matrix used in direct light-curve modelling.

An answer to a question if DRW was a good model describing the AGN variability would be yes. Yes, because both DRW and non-DRW processes described by the power exponential covariance matrix of the signal are very well modelled by a single DRW process. And yes, because Kozłowski (2016) based on model-independent SFs shows that β ≈ 1 for 9000 SDSS AGNs, consistent with DRW. The caveat is, however, that some of the underlying processes may be non-DRW, as indicated by Mushotzky et al. (2011) and Kasliwal, Vogeley & Richards (2015) based on steeper than DRW power spectral distributions and SFs of Kepler AGNs (the conversion between the power spectral density modelled as a single power law with the slope α is β = −0.5α, so we have explored in this paper −1 < α < −3, where α ≈ −3 was reported by Mushotzky et al. (2011)). Then the DRW light-curve modelling will not be able to identify a non-DRW process and will return biased DRW variability parameters.

As a matter of fact, in Kozłowski et al. (2010), we already modelled deterministic processes (non-stochastic) such as periodic variable stars, and found that DRW modelling provides a good description, where the time-scale τ is identified with the variability period.

I am grateful to Chris Kochanek for discussions of this topic. This work has been supported by the Polish National Science Centre OPUS grant number 2014/15/B/ST9/00093 and MAESTRO grant number 2014/14/A/ST9/00121.