How Can the Optical Variation Properties of Active Galactic Nuclei Be Unbiasedly Measured?

The variability of active galactic nuclei (AGNs) is ubiquitous, but has not yet been understood. Measuring the optical variation properties of AGNs, such as the variation timescale and amplitude, and then correlating them with their fundamental physical parameters has long served as a critical way of exploring the origin of AGN variability and the associated physics of the accretion process in AGNs. Obtaining accurate variation properties of AGNs is thus essential. It has been found that the damped random walk process can describe the AGN optical variation well, but there is a controversy over how long a minimal monitoring baseline is required to obtain unbiased variation properties. In this work, we settle the controversy by exhaustively scrutinizing the complex combination of assumed priors, adopted best-fit values, ensemble averaging methods, and fitting methods. The new proposal is then an optimized solution where unbiased variation properties of an AGN sample possessing the same variation timescale can be obtained with a minimal baseline of about 10 times their variation timescale. Finally, the new optimized solution is used to demonstrate the positive role of the time-domain surveys to be conducted by the Wide Field Survey Telescope in improving constraints on AGN variation properties.


INTRODUCTION
It is plausible that every massive galaxy contains a supermassive black hole (BH), into which gas swirls and behaves as active galactic nuclei (AGNs), including both low-luminosity AGNs and luminous quasars.Ever since the discovery of quasar, the variable nature of AGN emission has been reported.So far, its physical origin is generally explored through correlations between the variation properties, e.g., timescale (τ ) and amplitude (σ), of AGN light curves (LCs) and other physical parameters of AGNs, e.g., BH mass and Eddington ratio (Kelly et al. 2009;Burke et al. 2021).However, the physical origin of AGN variability is hitherto unclear.
Regardless of the physical origin, the optical AGN LCs with timescales of months to several years are found to be well described by a first-order continuous-time auto-regressive (CAR(1)) process or more commonly Damped Random Walk (DRW) process (Kelly et al. 2009;Zu et al. 2013).Note departures from the DRW tie@mail.ustc.edu.cn;zcai@ustc.edu.cnprocess have been reported at timescales shorter than several days (Mushotzky et al. 2011;Zu et al. 2013) and wavelengths of the rest-frame extreme ultraviolet (Zhu et al. 2016).Then, more sophisticated descriptions such as a p-order continuous-time auto-regressive q-order moving average (CARMA(p, q)) process (Kelly et al. 2014 and references therein) or simply a damped harmonic oscillator (i.e., CARMA(2,1)) model (Yu et al. 2022) have been explored.Nevertheless, parameters of the CARMA(p, q) process, more complicated than CARMA(1,0), i.e., CAR(1) or DRW, are too difficult to be physically understood, and correlations between those parameters and physical properties of AGNs are cumbersomely interpreted.Instead, the DRW is a relatively simpler model to be interpreted than the higherorder CARMA descriptions, and the two major parameters (τ and σ) of the DRW process could be, though controversial, attributed to specific physical meanings: τ , also coined by the damping timescale, may indicate the typical time required for the AGN LC to become roughly uncorrelated, while σ 2 the long-term (i.e., time interval ∆t ≫ τ ) variance of the AGN LC (Kelly et al. 2009).Therefore, the DRW process is still widely assumed to fit or simulate the AGN LCs (MacLeod et al. 2010;Koz lowski 2016Koz lowski , 2017Koz lowski , 2021;;Hu & Tak 2020;Suberlak et al. 2021;Kovačević et al. 2021;Stone et al. 2022), to select AGN candidates via variability (Koz lowski et al. 2010;Lei et al. 2022) as well as to simulate the thermal fluctuations of accretion disk (Dexter & Agol 2011;Cai et al. 2016Cai et al. , 2018Cai et al. , 2020)).
Thus, accurately measuring the variation properties of AGN LCs is essential, but many observational conditions, such as limited baseline and sparse or irregular sampling, can significantly affect measuring the variation properties of AGNs, especially τ .Therefore, DRW simulations have been extensively performed to assess to what extent we can confidently retrieve the input intrinsic variation properties of AGNs.However, confusing conclusions have been made.For example, to obtain without systematic bias an output measured timescale, τ out , of an observed AGN LC, Koz lowski (2017) suggest the observed baseline must be at least 10 times longer than the intrinsic timescale, τ in , but Suberlak et al. (2021) claim a factor of ∼ 3−5 is already sufficient.Oppositely, Koz lowski (2021) demonstrate a baseline of more than 30τ in is indispensable.
To end with these conclusions, except for the same assumption on the DRW process, there are more distinct statistical assumptions, including priors on the DRW parameters and estimators as the "best-fit" values for the DRW parameters.Therefore, to exhaustively explore the effects of all these selections and the origin of the aforementioned confusing conclusions, we in this work revisit the baseline required to accurately retrieve τ in by comparing different priors and estimators in Section 2. Effects of photometric uncertainty, sampling cadence, and season gap on retrieving the DRW parameters are considered in Section 3, before highlighting the role of the 2.5-meter Wide Field Survey Telescope (WFST; Wang et al. 2023a) in confining the variation properties of AGNs.Finally, a brief summary is presented in Section 4. The data and codes used in this work are available on Harvard Dataverse: doi:10.7910/DVN/I4XWGU.
2. DRW SIMULATION AND FITTING 2.1.Simulating AGN LCs as a DRW process Suppose X(t) is an AGN LC in magnitude as a function of epoch t, which can be described by the DRW process following the Itô differential equation (Kelly et al. 2009;Brockwell & Davis 2016) as where τ is the damping timescale, σ s the characteristic amplitude of variations per day1/2 (the short-term variances ≈ σ 2 s ∆t for ∆t ≪ τ and connected to the longterm one by σ 2 = τ σ 2 s /2), B(t) the standard Brownian motion, and cτ the mean of X(t).Practically, the value of X(t) given X(s) for s < t is where ∆t = t − s and N (0, 1) is the normal distribution with zero mean and one variance.We then use Equation ( 2) to simulate AGN LCs.In particular, the starting point of the simulated LC is X(0) ∼ N (cτ, σ 2 ).Note that the covariance function of the DRW process is Cov(∆t) = Cov(X(t), X(s)) = σ 2 e −∆t/τ (3) and the corresponding structure function is where the asymptotic rms variability on long timescales SF ∞ = √ 2σ (MacLeod et al. 2010;Zu et al. 2013).

Simulation setting
Globally following Koz lowski (2017) and Suberlak et al. (2021), we fix an observed baseline of T = 8 yr and revisit how well the input intrinsic τ in relative to T can be retrieved.For this purpose, we consider 61 different τ in evenly distributed from 0.001T to T in steps of 0.05 dex.For each τ in , we simulate 10 4 different LCs with fixed SF ∞ = 0.2 mag (or σ = 0.14 mag), and we also consider two kinds of cadences: SDSS-like and OGLE-like with N = 60 and N = 445 epochs, respectively.To mimic the real observation, simulated epochs are randomly distributed at night and, otherwise specified, only one epoch is allowed within 3 hours before and after midnight.For SDSS-like and OGLE-like observations, we consider fixed typical mean magnitudes, ⟨m⟩, for AGNs with ⟨r SDSS ⟩ = 17 mag and ⟨I OGLE ⟩ = 18 mag, but magnitude-dependent photometric uncertainties as Suberlak et al. (2021): Hereafter, we refer to an SDSS-observed AGN with conditions of ⟨r SDSS ⟩ = 17 mag, T = 8 yr, N = 60, and σ e = σ SDSS ≃ 0.013 mag, while an OGLE-observed one with ⟨I OGLE ⟩ = 18 mag, T = 8 yr, N = 445, and σ e = σ OGLE ≃ 0.025 mag.Note here σ OGLE ≃ 2σ SDSS .
Specifically, the simulated LC is y(t) = s(t) + n(t), where s(t) is the DRW process around ⟨m⟩ and the observational uncertainty, n(t), is randomly drawn from a normal distribution, N (0, σ 2 e (t)).

Fitting DRW parameters
Several libraries are available for fitting the DRW parameters, such as javelin 1 (Zu et al. 2011) (Yu & Richards 2022;Yu et al. 2022).Among them, the javelin utilizes the deterministic amoeba method to maximize the so-called PRH likelihood (see Appendix A), while the others are Monte Carlo Markov Chain (MCMC) samplers over likelihood functions.The carma pack is an MCMC sampler for inferring parameters of the CARMA process, the celerite is capable of fast modeling different kinds of one-dimensional Gaussian processes, including the DRW process, and the EzTao is an easier CARMA modeling based on celerite.
In this work, we assume the DRW kernel (Equation 3) and fiducially compute the likelihood functions for the DRW parameters with celerite, complemented with 104 MCMC samplings by emcee5 .We use the default "stretch move" algorithm in emcee for the actual MCMC sampling and confirm that using different MCMC sampling algorithms, such as the "differential evolution" or "walk" ones contained in emcee, does not alter our conclusions.

Priors on and estimators for the DRW parameters
Assessing the quality of the fit and determining the "best-fit" value for the DRW parameters are challenging tasks.Given an observed or simulated AGN LC, we can directly estimate a potential "best-fit" value for its DRW parameters from the corresponding likelihood functions using the Maximum Likelihood Estimation (MLE) without priors on the DRW parameters.Instead, if priors are assumed, there are three potential Bayesian estimators (Wei 2016), i.e., Maximum A Posterior (MAP), Posterior Expectancy (PE), and Posterior Median (PM), for the "best-fit" values assessed from the posterior probability distribution which is the product of the likelihood function with a prior probability distribution of the DRW parameter.Note the MAP, PE, and PM measure the peak value, the mean value (or the one-order moment), and the 50% percentile of the posterior distribution, respectively.
Here, the MLE and MAP are obtained by maximizing the likelihood function and the posterior probability distribution, respectively, by virtue of the L-BFGS-B algorithm (Byrd et al. 1995;Zhu et al. 1997), while both PE and PM are derived from the posterior probability distributions constructed using 10 4 MCMC samplings.Compared to fully explore the parameter space on fixed fine grid for the posterior probability distribution, the introduction of MCMC sampling is computational cheaper.
Previous studies have employed distinct priors and estimators for the DRW parameters and reached quite confusing conclusions.For example, Koz lowski (2017) assumes priors of P (τ, σ) = 1/σ √ τ (or equivalently P (τ, σ s ) = 1/σ s τ ; hereafter K17 prior; see also MacLeod et al. 2010;Koz lowski et al. 2010) and uses the MAP estimator, suggesting that a baseline longer than at least 10 times intrinsic τ in is required for accurately retrieving τ in .With the same prior as Koz lowski (2017), even a baseline as long as 30 τ in has been suggested to be indispensable for retrieving unbiased DRW parameters (Koz lowski 2021).However, Suberlak et al. (2021) specify priors of P (τ, σ) = 1/στ (hereafter S21 prior) and adopt the PE estimator, claiming that a shorter baseline with only 3 − 5 times τ in is sufficient for retrieving unbiased DRW parameters.
Note that both Koz lowski (2017) and Suberlak et al. (2021) take the ensemble median value of hundreds of recovered timescales from simulated LCs as the final "bestfit" timescale for an input τ in .
The left panel of Figure 1 intuitively illustrates the prominent differences among aforesaid four estimators for the "best-fit" values of τ given a typical simulated LC of an SDSS-observed AGN with τ in = T /10.Therefore, for the purpose of this work, we exhaustively explore all combinations of the two priors and four estimators, as well as both the ensemble mean and median, to determine which estimator is more applicable to the assumed DRW process.

Determining the "best-fit" value for the DRW parameters
For each input τ in , we simulate 10 4 LCs.Given the i-th simulated LC, we may find several potential estimators for the "best-fit" value of the DRW parameters, such as τ i type with type = MLE, MAP, PM, or PE for the DRW τ parameter.As illustrated in the left panel of Figure 1, the four estimators are quite different and, for this specific LC, τ PM happens to be consistent with τ in .
Owing to the randomness of AGN variability, values of the same estimator would be different from one LC to another.Therefore, we treat the ensemble mean (median) value of {τ i=1,...,10 4 type } as the output τ mean out (τ median out ) to be compared with the input τ in .For example, the right panel of Figure 1 illustrates the probability distribution of τ PM retrieved from 10 4 simulated LCs assuming the K17 prior.In this case, τ mean out is found to be in stable agreement with τ in , while τ median out is always smaller.Note that there are individual τ PM much smaller (larger) than τ in .It is found that they are typically retrieved from LCs with smaller (larger) variation amplitude, confirming the bias due to the limited baseline reported by Koz lowski (2021).
Furthermore, for a series of τ in , Figures 2 and 3 illustrate to what accuracy can we measure the DRW τ parameter for SDSS-observed and OGLE-observed AGNs, respectively, by considering all combinations of two priors, four estimators, and two ensemble average methods.Globally, estimators inferred from the S21 prior are more or less smaller than those from the K17 prior.Roughly, we find τ PE > τ PM ≳ τ MLE > τ MAP , while differences among them complicatedly depends on T /τ in .For the same prior and estimator, τ mean out is somewhat larger than τ median out , especially around T /τ in ∼ 10.Comparing τ out to τ in , there are prominent departures at either too long average cadence or too short baseline relative to τ in .The effect of the former can be found in Figure 2 for the SDSS-observed AGNs, whose average cadence ∆T is only ∼ T /60, resulting in ∆T /τ in > 1 when T /τ in > 60.Without priors, τ MLE /τ in increases monotonically with increasing T /τ in since the information of the short-term variation is lost more significantly for larger ∆T given a fixed number of epochs.Instead, once a prior on τ is assumed, a competition between the likelihood without prior (preferring larger timescales) and the prior (preferring smaller timescales for those discussed here) results in the special dependence (like a valley6 ) of τ MAP , τ PE , and τ PM on T /τ in when LCs are too sparsely sampled at T /τ in > 10 2 (Figure 2).Although by increasing the number of epochs or shortening cadence the effect of the special dependence becomes marginal (Figure 3), unavoidable is the significant departure at too short baselines with T /τ in < 10 (Fig- ures 2 and 3).Therefore, in terms of Figure 3, we look for the best combinations that can recover τ in as much as possible.
As a rule of thumb, we find that there are two best combinations: K17 prior plus ensemble mean of τ PM (hereafter K17PMm solution) and S21 prior plus ensemble median of τ PE (hereafter S21PE50 solution).Both of them can quite well retrieve τ in down to T /τ in ≃ 10, within uncertainties of 2% and 4% for the K17PMm and S21PE50 solutions, respectively (Figure 3).For other combinations, T /τ in as large as ≃ 100 is a prerequisite for achieving comparable accuracy as the two best combinations.
Employing the MAP estimator to determine the DRW parameters, Koz lowski (2021) point out that underestimated variances in shorter AGN LCs with T /τ in ≲ 30 lead to underestimated timescales as compared to τ in .Contrarily, utilizing the K17PMm and S21PE50 solutions, timescales can be retrieved quiet well down to T /τ in ≃ 10.Since the variance and timescale are correlated, we confirm in Figure 4 that the MLE-estimated variances of our simulated LCs are nicely consistent with the prediction of DRW LCs (Koz lowski et al. 2010), and that variances retrieved using the K17PMm and S21PE50 solutions are not significantly biased down to T /τ in ≃ 10, especially for the K17PMm solution.Therefore, timescales down to T /τ in ≃ 10 can be well retrieved with the K17PMm solution.
Besides less unbiased variances retrieved using the K17PMm solution down to T /τ in ≃ 10, we also find that, for both the timescale and variance, the mean square error (MSE) of the K17PMm solution is at least half of that of the S21PE50 solution, thus the K17PMm solution is recommended as the final best way of determin- ing the "best-fit" value for the DRW parameters down to T /τ in ≃ 10.

Ours versus Koz lowski (2017)
To determine the DRW parameters, Koz lowski (2017) assume the K17 prior plus ensemble median of τ MAP (hereafter, K17MAP50 solution) and find that a minimum of T /τ in ≃ 10 is required for meaningfully measuring the timescale of AGN LC.However, according to the dot-dashed curve presented in the top-right panel of Figure 3, our K17MAP50 solution does not perform well at T /τ in ∼ 10, contrary to the statement of Koz lowski (2017).Note, rather than using celerite we adopted here, Koz lowski (2017) utilize the so-called PRH (Press et al. 1992) method summarized in Koz lowski et al. (2010).This fitting method is also implemented in javelin (S.Koz lowski 2023, personal communication; see also Appendix A).
To understand the difference between us and them, we repeat the same simulation but determine the DRW parameters using the PRH method (S.Koz lowski 2023, personal communication for his PRH code).As illustrated in the left panel of Figure 5, the K17MAP50 solution involving the PRH method indeed performs much better down to T /τ in ≃ 10, consistent with Koz lowski (2017).The slight underestimation of τ out at T /τ in ∼ 10 in Figure 2 of Koz lowski ( 2017) is not so prominent owing to the fact that their τ out /T is presented in logarithmic scale (see also Suberlak et al. 2021).Instead, presenting τ out /τ in in linear scale as we choose here can sensitively highlight the difference between τ out and τ in .In the left panel of Figure 5, from T /τ in ≃ 1 through 10 to 100, τ out derived from the PRH method is larger than that from celerite by a factor of ≃ 2.7 through ≃ 1.3 to ≃ 1.1, as a result of the prominent difference in the likelihood functions and larger MAP implied by the PRH method (See Appendix A).Therefore, we conclude that, in the sense of well retrieving the DRW parameters, the K17MAP50 solution is likely more compatible with the PRH fitting method but not celerite.

Ours versus Suberlak et al. (2021)
As discussed in the last subsection, the S21PE50 solution also works quite well on retrieving τ in down to T /τ in ≃ 10.The same solution is adopted by Suberlak et al. (2021), but they suggest that T /τ in ∼ 3−5, shorter than what we find here, is sufficient for retrieving the timescale.Note that Suberlak et al. ( 2021) state, without quantitative assessment, that as long as T /τ in ≳ 3 they can recover the timescale without substantial bias.However, this is likely not true as we should demonstrate in the following.
Except adopting the same fitting method (i.e., celerite) and S21PE50 solution, there are still three major differences between their simulation7 and ours: 1.The starting points, s(0), of simulated LCs are all fixed to the mean magnitude, ⟨m⟩, by Suberlak et al. ( 2021), while we free s(0) by drawing it from a normal distribution, N (⟨m⟩, σ 2 ), mimicking a long enough "burn-in".
2. When estimating τ PE , Suberlak et al. (2021) only consider a fixed sparse grid (3600 points, that is, 60 points for τ from 1 day to 5000 days and 60 points for σ from 0.02 to 0.7 mag, both of them are evenly spaced in logarithm) for the two DRW parameters, while we adopt the MCMC approach (Foreman-Mackey et al. 2013) to approximate the posterior distributions for the two DRW parameters.
3. For each τ in , τ out is the median of τ PE retrieved from only 100 simulated LCs, while we simulate 10 4 LCs, pursuing a stable estimation for the ensemble median.
As illustrated in the right panel of Figure 5, even both Suberlak et al. (2021) and us have adopted the same fitting method and S21PE50 solution, differences on simulating LCs and exploring the parameter space indeed induce some prominent diversities.Using a fixed sparse grid for the DRW parameters would overestimate τ out as compared to using the MCMC approach, while fixing the start point would give rise to smaller retrieved τ out .Thus, considering both a fixed grid for the DRW parameters and fixing the start point as done by Suberlak et al. (2021) surprisingly results in τ out comparable to what we find, so we do not confirm their finding that T /τ in ∼ 3 − 5 is sufficient for retrieving the timescale.Instead, we show that T /τ in ≃ 10 or even larger is essential under the same approach as Suberlak et al. (2021).
In the right panel of Figure 5, we further demonstrate that achieving unbiased fitting with T /τ in ∼ 3 − 5 is likely attributed to the randomness of the ensemble median of τ out , which are based on only 100 simulated LCs as done by Suberlak et al. (2021).

DISCUSSIONS
As introduced in Section 2 and summarized in Table 1, we have demonstrated that the K17PMm solution works quite well on retrieving the DRW parameters of an AGN sample down to T /τ in ≃ 10.Here, adopting the K17PMm solution, we further consider effects of the photometric uncertainty (σ e ), the sampling cadence (with ⟨∆t⟩ for the average cadence), and the season gap on the retrieved DRW parameters.Then we address how long the baseline is necessary in order to determine the DRW parameters for a single AGN down to a desired accuracy.Finally, time domain surveys to be conducted by the WFST (Wang et al. 2023a) are briefly introduced and the role of WFST surveys in improving constraints on the DRW parameters is highlighted.
In the following comparisons, we consider an AGN sample observed with a baseline of T ≳ 10τ in and an average cadence ⟨∆t⟩.In other words, an AGN would be observed at N epochs, where N ≃ T /⟨∆t⟩.These observed epochs are randomly distributed, rather than evenly distributed, within the baseline because we confirm that for random epochs the DRW parameters can still be well retrieved even with a quite large ⟨∆t⟩ since (5) The sampling method for the model parameters: deterministic, fixed sparse grid, or MCMC sampling using emcee; (6) The estimator for the "best-fit" model parameter of a single LC: MAP, PE, or PM; (7) The method averaging the "best-fit" model parameters retrieved from an ensemble LCs: either median or mean; (8) The abbreviation for the solution mentioned in this work; (9) The minimal T /τin required for achieving a 2% accuracy for the retrieved τ and are estimated using the orange solid line in the left panel of Figure 5, the black solid line in the right panel of Figure 5, and the green solid line in the top-left panel of Figure 3 for the Koz lowski (2017), Suberlak et al. (2021), and our solutions, respectively; (10) Same as the ninth column but for achieving a 10% accuracy.
short-term variations can be kept by random observations.To reduce the randomness, we simulate 10 4 LCs for each set of conditions, i.e., {τ in , σ, T , ⟨∆t⟩, σ e }.

Photometric uncertainty
Relative to the long-term variance (σ) of AGN LC, we compare three different photometric uncertainties with σ e = 0, 0.1σ, and 0.5σ.For an AGN sample observed with T = 10τ in and fine enough ⟨∆t⟩ (i.e., ⟨∆t⟩ < 0.2τ in ; see the left panel of Figure 6), unbiased timescales can be retrieved for all σ e even though a slightly larger dispersion is found for very large σ e .Interestingly, decreasing σ e from 0.1σ to 0 does not further decrease the dispersion of the retrieved timescales, suggesting that for small photometric uncertainty celerite can take good care of the scatter imposed by the photometric uncertainty.
In addition, the measured errors on the LCs could not be completely accurate.Thus, an extra white noise term has been introduced to account for an additional source of photometric error (e.g., Burke et al. 2021;Wang et al. 2023b).However, we find this extra noise term has little effect on our conclusions (see Appendix B).

Sampling cadence
The minor effect of the average cadence on the retrieved timescale is clearly illustrated in Figure 6.For T = 10τ in , the timescales can be well retrieved as long as ⟨∆t⟩ < 0.2τ in , with which the dispersion of retrieved timescales is independent of ⟨∆t⟩ (left panel of Figure 6).Instead, the dispersion can only be reduced by increasing the baseline, e.g., from T = 10τ in to 100τ in (right panel of Figure 6).
Moreover, when the baseline is 100τ in , unbiased timescales can be obtained even with ⟨∆t⟩ ≃ 1τ in .This is because when the baseline increases from 10τ in to 100τ in and the unbiased timescale is obtained at ⟨∆t⟩ increasing from ≃ 0.2τ in to ≃ 1τ in , the number of observed epochs indeed increases from 50 to 100 such that the short-term variations can still be somewhat covered by random observations.Note that ⟨∆t⟩ ≃ 1τ in does not mean that there are no samplings with cadences smaller than τ in .Instead, there are ∼ 65% samplings with cadences smaller than τ in .For comparison, there are ∼ 99% samplings with cadences smaller than τ in for ⟨∆t⟩ ≃ 0.2τ in .

Season gap
Real observations are affected by the season gap, which means that some AGNs can not be continuously observed throughout the year.To investigate the effect of season gap on the retrieved timescale, we simulate an AGN sample observed by 10 years, but within each year the observable duration is limited to 3 months and 6 months.A quite sparse cadence of ⟨∆t⟩ = 12 days is selected such that AGNs are observed at ≃ 8, 15, and 30 epochs for 3, 6, and 12 months duration, respectively.For comparison, the well-known SDSS Stripe 82 is mapped 8 times on average in a 2-to-3 months duration per year and the average cadence is about 6 to 7 days (MacLeod et al. 2010), while a 3-day cadences is reached by the subsequent time domain surveys, such as Pan-STARRS1 (PS1; only in the median deep fields; Chambers et al. 2016) and ZTF (Masci et al. 2019).The upcoming deep drilling fields of LSST will have a 2-day cadence (Brandt et al. 2018), and a 1-day cadence is proposed for the deep high-cadence survey of WFST (Wang et al. 2023a).out /τin ≃ 1) can be retrieved with the average cadence ⟨∆t⟩ shorter than 0.2τin, while the dispersion of the retrieved timescales (i.e., 16%-84% percentile range based on 10 4 LCs) does not significantly depend on the photometric uncertainty.The 16%-84% percentile range for the case with σe ≃ 0.1σ is indicated by the green region, while those for the other two cases are indicated by the corresponding thinner lines with the same color and line style.Right panel: similar to the left panel but for an AGN sample observed with σe ≃ 0.1σ and T = 10τin and 100τin.With a longer baseline, unbiased timescales can be retrieved with larger average cadence ⟨∆t⟩ and the corresponding 16%-84% dispersion is narrower.As illustrated in Figure 7, unbiased timescales can be well retrieved for τ in ≲ 300 days or again ≲ 0.1T , regardless of the presence of season gap.Indeed, larger dispersion is induced by the season gap.Owing to the given correlation between variance and timescale of AGN LC (Koz lowski 2021), the retrieved timescales are slightly smaller than τ in as a result of the smaller variations implied by LCs with larger season gap.

Can the timescale of individual AGN be accurately determined?
Unbiased timescales of an AGN sample can already be well retrieved given T ≃ 10τ in and ⟨∆t⟩ ≃ 0.2τ in , but with large dispersion (left panel of Figure 6).In other words, for individual AGN the retrieved timescale could be different from τ in , e.g., by at most ∼ 50% for a 68% probability.
In terms of above analyses, the only effective way to reduce the dispersion of retrieved timescales is by increasing the baseline.In Figure 8, fixing σ e ≃ 0.1σ, dis- persion of retrieved timescales is illustrated as a function of the baseline T from 10τ in to 1000τ in for two cadences of ⟨∆t⟩ = 0.2τ in and 1τ in .For ⟨∆t⟩ = 0.2τ in , the dispersion is just gradually decreasing with increasing baseline.The 16%-84% dispersion of τ out /τ in decreases from ∼ 50% through ∼ 15% to ∼ 5% for baseline increasing from 10τ in through 100τ in to 1000τ in , respectively.
Even for ⟨∆t⟩ = 1τ in , unbiased timescale can be obtained with T ≳ 100τ in and the 16%-84% dispersion of τ out /τ in can decrease from ∼ 20% to ∼ 8% for baseline increasing from 100τ in to 1000τ in (Figure 8).This has a practical implication that, for intermediate-mass BHs whose τ in are likely several days (Burke et al. 2021), continual monitoring up to several decades with ⟨∆t⟩ ∼ τ in would be sufficient for obtaining τ in with ∼ 20% accuracy.
3.5.The role of WFST in constraining the optical variation properties of AGNs

Survey Strategy of WFST
Brief introductions to the WFST surveys and the relevant AGN sciences are presented in the following, while readers are referred to the WFST white paper for a panchromatic view (Wang et al. 2023a).The 2.5-meter WFST with a field of view of 6.5 square degrees is designated to quickly survey the northern sky in four optical bands (u, g, r, and i).There will be two planned key programs across 6 years: a deep high-cadence u-band survey (DHS) and a wide field survey (WFS).
The DHS program tends to cover ∼ 720 deg 2 surrounding the equator and reaches a depth of ≃ 23 and ≃ 24 mag (AB) in u and g bands in a 90-second exposure, respectively.Two separate DHS fields of ∼ 360 deg 2 would be continuously monitored for 6 months per year in two to three bands (probably more u and less i).At least 1-day cadence is allocated to each band (probably except u around the full moon and i around the new moon).These quasi-simultaneous observations are dedicated to systematically unveiling the multi-band continuum lags of AGNs (Z.B. Su et al. 2023, in preparation).
The WFS program would cover ∼ 8000 deg 2 in the northern sky and reach a depth of ≃ 22.2 and ≃ 23.3 mag (AB) in u and g bands in a 30-second exposure, respectively.Four separate WFS fields of ∼ 2000 deg 2 would be continuously monitored for 3 months per year in the four bands.On average, a 6-day cadence is allocated to each band, or ≃ 15 visits per 3 months per band.
Both DHS and WFS are valuable for constraining the variability properties of AGNs, especially the u band.According to the WFST schedule and the first light on September 17, 2023, we hereafter assume that the formal WFST scientific observation would start in spring 2024 and consider three WFST-extended baselines of 1, 6, and 10 years, coined as W1, W6, and W10, respectively.

Beneficial complements to archive surveys
The DHS footprint is planned to entirely cover the well-monitored SDSS Stripe 82 (S82) region (≃ 290 deg 2 ; MacLeod et al. 2010), while the other DHS and WFS footprints would be enclosed by that of ZTF.The legacy surveys of WFST will provide time domain optical data with a cadence denser than SDSS S82 and a waveband shorter than ZTF.
In the following, we perform simulations by mocking WFST observations on AGNs to demonstrate how well can the extension of WFST baseline improve measuring variation properties of AGNs.We consider a total baseline of 24, 30, and 34 years, starting from the SDSS observations at around 2000 to 2024, 2030, and 2034 and including the 1, 6, and 10 years of WFST observations, respectively.On the basis of 9258 SDSS quasars in S82 region (MacLeod et al. 2012), we search for their PS1 counterparts within one arcsec, resulting in 9254 quasars with 16 PS1 epochs complement to 62 SDSS epochs on average.Then we simulate 9254 LCs for any set of timescale (τ in = 10 -1000 days) and amplitude (σ = 0.14 mag), and implement the observed r-band epochs and photometric errors of 9254 quasars to the simulated 9254 LCs.Two relative photometric uncertainties (σ W = 0.1σ and 0.5σ) and two typical even cadences (∆t W = 1 day and 6 days for DHS and WFS, respectively) are considered for comparison.Finally, variation properties of AGNs are retrieved using the afore-proposed optimized solution (i.e., K17PMm; Section 2).Note at the moment the similarity in r-band transmission curves among SDSS, PS1, and WFST leads to neglect the small corrections from SDSS or PS1 to WFST photometry.
As illustrated in Figure 9, unbiased mean timescales of quasars in the S82 region can be retrieved as long as ≃ 400 − 600 days for the fiducial 6-year WFST observation on the S82 region with a total baseline of ∼ 30 years, regardless of the relative photometric uncertainties and cadences.Compared to the 6-year WFST observation, the 1-year WFST observation is already valuable, as a result of nearly doubled baseline increasing from ∼ 13 years (SDSS+PS1) to ∼ 24 years (SDSS+PS1+W1), while the help of extending 6-year to 10-year WFST observations is small because the final baseline is further increased by only ∼ 13%.Finding low-luminosity AGNs with intermediate-mass BHs (IMBHs) in dwarf galaxies is a challenge because the weak AGN signal is overwhelmed by the star-forming activity when conventional imaging and spectroscopic methods are adopted.Instead, variability is probably a promising avenue for selecting IMBH candidates (Baldassare et al. 2018;Greene et al. 2020) and the measured timescale could be correlated with the BH mass Figure 9. Adopting the real r-band cadences and photometry uncertainties of 9254 quasars in the S82 region with both SDSS and PS1 observations, we predict how well the retrieved timescales can be achieved by extending WFST observations, compared to solely considering SDSS+PS1 observations (dot-dashed lines).For 9254 AGNs with σ = 0.14 mag and the fiducial WFST baseline of 6 years (W6; or a SDSS+PS1+W6 baseline of 30 years), we consider two WFST photometry uncertainties (i.e., σW ≃ 0.1σ and 0.5σ; top-left panel) and two even cadences (i.e., ∆tW = 1 day and 6 days; top-right panel).Although dispersions (16%-84% percentile; green region or enclosed by two thin lines with the same style) of the retrieved timescales are larger for larger WFST photometry uncertainty, unbiased mean timescales can still be retrieved as long as ≃ 400 − 600 days for extending 6-year WFST observation (+W6).Compared to +W6, extending the 1-year WFST observation (+W1) is already valuable, while the help of extending the 10-year WFST observation (+W10) is small, in retrieving unbiased mean timescales (bottom panel).
of IMBH, though still subject to the lack of a sample of IMBHs with well measured BH masses in order to calibrate the intrinsic scatter of the M BH − τ relation at low luminosities (Burke et al. 2021;Wang et al. 2023b).
Nevertheless, to assess to what extent can the WFST DHS be used to retrieve unbiased timescales for IMBHs with ∼ 1 − 10 days, probably corresponding to M BH ∼ 10 4 − 10 5 M ⊙ (Burke et al. 2021), we mock 6-monthper-year WFST DHS observations for IMBHs with τ in ∼ 0.1 − 10 days and fixed σ = 0.14 mag.Although some IMBHs can have variation amplitudes as large as ∼ 0.14 mag, the typical variation amplitudes of IMBHs are likely ∼ 0.03 mag or even smaller (Baldassare et al. 2018).Therefore, we consider two WFST photome-try uncertainties, σ W , relative to the fixed σ, that is, σ W = 0.1σ and 0.5σ.
The left panel of Figure 10 shows, with the 6-month WFST DHS observations in the first year survey, unbiased timescales of IMBHs from ∼ 1 day to 10 days can be well retrieved, even for the case of σ W ≃ 0.5σ.If IMBHs within the WFST DHS fields are continually monitored for 6 years (the right panel of Figure 10), dispersions of retrieved timescales will be smaller than ∼ 10% and the minimal unbiased timescale decreases to ≃ 0.5 day.Therefore, we expect the WFST DHS would be helpful in retrieving accurate and unbiased timescales for individual IMBHs likely glowing as low-luminosity AGNs in dwarf galaxies.4. SUMMARY In this work, assuming the optical AGN variability is DRW-like, we examine how can the variation properties of AGNs be unbiasedly measured using the mainstream celerite.We assess the effects of different priors and distinct "best-fit" values and find two best-combined estimators, that is, the K17 prior (1/σ √ τ ) plus the ensemble mean of PM values (K17PMm solution) and the S21 prior (1/τ σ) plus the ensemble median of PE values (S21PE50 solution).Since the MSE of the former is smaller, we propose the K17PMm solution as an optimized method to infer variation properties of AGNs and demonstrate that a minimal baseline 10 times longer than the input variation timescale is essential.Although the proposed solution may not perform well for the non-DRW-like AGN variability, our procedure in unveiling the corresponding optimized solution can be applied.
After scrutinizing our findings to those of Koz lowski (2017), Suberlak et al. (2021), and Koz lowski (2021), we find different combinations of priors, "best-fit" values, and fitting methods indeed result in distinct conclusions on the minimal baselines required to retrieve unbiased variation timescales.Koz lowski (2017) use the PRH fitting method and also report a minimal baseline of T /τ in ∼ 10.We find the PRH fitting method should only be combined with another specific estimator, that is, the K17 prior plus the ensemble median of MAP values (K17MAP50 solution).Instead, our solution performs somewhat better than theirs and the utility of the fast celerite fitting method has a promising application in the time domain astronomy.Suberlak et al. (2021) suggest a shorter requested minimal baseline of T /τ in ∼ 3 − 5, but we can not confirm their finding by freeing the start points of simulated 10 4 LCs.Therefore, we conclude the shorter minimal baseline suggested by them is attributed to the bias and randomness as a result of fixing the start point of simulated 100 LCs.
Koz lowski (2021) claims a longer requested minimal baseline of T /τ in ∼ 30 is necessary and attributes the underestimation of timescales to the smaller variances of short LCs.Instead, our solution performs better on estimating the true long-term variability amplitude and so can shorten the requested minimal baseline to T /τ in ∼ 10.
Furthermore, utilizing the new optimized solution, we examine the impacts of several observational factors, including baseline T , number of epochs N (cadence ∆t), photometric uncertainty, and seasonal gap.Our findings are listed as follows: 1. Larger photometric uncertainty only slightly increases the dispersion of retrieved timescales.
2. The number of epochs should be sufficient, but excessive epochs are not very helpful.For T = 10τ in , timescales can be well retrieved as long as the average cadence ⟨∆t⟩ = 0.2τ in .For a longer baseline T = 100τ in , a sparse average cadence ⟨∆t⟩ = τ in is sufficient.Further decreasing cadences does not significantly reduce the dispersion of retrieved timescales.
3. Seasonal gaps do not affect the result obviously, except for adding little dispersion and slightly depressing the ensemble mean of retrieved timescales as a result of smaller variance typically implied by LCs with larger season gaps.
4. For sufficient small cadences, the retrieved timescale of individual AGN can reach an accuracy of ∼ 50% (1σ confidence level) for T = 10τ in , while extending baseline to 100τ in would increase the accuracy to ∼ 15% (1σ confidence level).
Finally, we apply our solution to evaluate constraints on the variation properties of AGNs to be provided by the WFST DHS.Complementing archive surveys in the SDSS S82 region, unbiased mean timescales of quasars covered by the WFST DHS can be retrieved as long as ≃ 400 − 600 days for the fiducial 6-year WFST DHS observation.Moreover, for the 1-year WFST observation, unbiased timescales as short as ≃ 1 day for IMBHs can be retrieved, while ≃ 0.5 day for the 6-year WFST DHS observation.Therefore, we expect the WFST DHS would be helpful in improving constraints on AGN variability.σ * (Figure 12).The best-fit σ * is also very small with a ratio of σ * to the average σ OGLE only several times 10 −3 .Therefore, we conclude that adding an extra noise term is unnecessary for us and has little effect on our conclusion.

Figure 1 .
Figure 1.Left panel: given a typical LC of an SDSS-observed AGN with τin = T /10 (i.e., τin = 292 days for given T = 8 years), the likelihood function for its DRW τ parameter is compared to the corresponding marginalized posterior probability distribution assuming the K17 prior, i.e., P (τ, σ) = 1/σ √ τ .Four distinct estimators (τMLE for the MLE estimator from the likelihood function, while τMAP, τPM, and τPE for the MAP, PM, and PE estimators from the posterior probability distribution, respectively) are compared to the input intrinsic τin.Right panel: shown is the distribution of τPM from 10 4 LCs simulated in the same way as presented in the left panel, while the mean and median values of the distribution are compared to τin.

Figure 2 .
Figure 2.For each given τin relative to the baseline T (x-axis), the ensemble mean/median timescales τ mean/median out relative to τin (y-axis) for four different estimators (i.e., τPE, τPM, τMAP, and τMLE; see the left panel of Figure 1) estimated from 10 4 LCs of SDSS-observed AGNs are shown in the left/right panels.When estimating τPE, τPM and τMAP, two priors are considered, that is, K17/S21 prior in the top/bottom panels.

Figure 3 .
Figure 3. Same as Figure 2, but for OGLE-observed AGNs.For the two combinations that can well retrieve τin down to T /τin ≃ 10, the 16%-84% percentile range of 10 4 recovered timescales is shown as the green (blue) region for the K17 prior plus mean of τPM (the S21 prior plus median of τPE) in the top-left (bottom-right) panel.

Figure 4 .
Figure 4.For OGLE-observed AGNs, σout retrieved using the K17PMm and S21PE50 solutions are compared to the mean of σMLE, retrieved using MLE without prior, which is found to nicely agree with the analytical formula worked out by Koz lowski et al. (2010) for DRW LCs, where x = T /τin.

Figure 5 .
Figure5.Left panel: considering the same simulations but different fitting methods for OGLE-observed AGNs, our ensemble median of τMAP retrieved using celerite as a function of T /τin is compared with that retrieved using the PRH method by Koz lowski (2017).Right panel: considering the same celerite but different ways of simulating LCs (i.e., fixed or free starting point to the mean of LC) and of estimating the DRW parameters (i.e., fixed grid or MCMC searching for the best parameters) for OGLE-observed AGNs, our ensemble median (based on 10 4 LCs) of τPE obtained using "free starting point + MCMC" is compared with that obtained using "fixed starting point + grid" as performed bySuberlak et al. (2021) and that obtained using "free starting point + grid".For the latter two results, we have adopted the code on GitHub provided bySuberlak et al. (2021).SinceSuberlak et al. (2021) take the median of τPE of only 100 LCs, we also illustrate 100 τ median out by taking the median of τPE of 100 LCs, randomly selected from the corresponding 10 4 LCs, simulated and fitted in the same way asSuberlak et al. (2021).

Figure 6 .
Figure6.Left panel: adopting the K17PMm solution to retrieve the DRW timescales for an AGN sample observed with T = 10τin and different photometric uncertainties (i.e., σe ≃ 0, 0.1σ, and 0.5σ, where σ is the long-term variance of AGN LC), unbiased timescales (i.e., τ mean out /τin ≃ 1) can be retrieved with the average cadence ⟨∆t⟩ shorter than 0.2τin, while the dispersion of the retrieved timescales (i.e., 16%-84% percentile range based on 10 4 LCs) does not significantly depend on the photometric uncertainty.The 16%-84% percentile range for the case with σe ≃ 0.1σ is indicated by the green region, while those for the other two cases are indicated by the corresponding thinner lines with the same color and line style.Right panel: similar to the left panel but for an AGN sample observed with σe ≃ 0.1σ and T = 10τin and 100τin.With a longer baseline, unbiased timescales can be retrieved with larger average cadence ⟨∆t⟩ and the corresponding 16%-84% dispersion is narrower.

Figure 7 .
Figure 7. Similar to Figure 6, but for an AGN sample observed with σe ≃ 0.1σ and a fixed baseline of T = 10 year, unbiased timescales can be well retrieved for τin ≲ 0.1T ≃ 300 days, regardless of the quite sparse cadence ⟨∆t⟩ = 12 days and of the limited observable duration of one year (i.e., 3 months and 6 months compared to the whole year).

Figure 8 .
Figure 8. Similar to Figure 6, but for an AGN sample observed with σe ≃ 0.1σ and ⟨∆t⟩ = 0.2τin and 1τin, the retrieved timescale and the 16%-84% dispersion as a function of the baseline T relative to τin.
3.5.3.Retrieving unbiased timescales as short as about one day for intermediate-mass BHs in the WFST DHS

Figure 10 .
Figure10.Left panel: for the WFST DHS with even cadence of ∆tW = 1 day, unbiased timescales as short as ≃ 1 day can be retrieved, regardless of the photometry uncertainties, for the 1-year WFST DHS observation with a baseline of 6 months.Dispersions of retrieved timescales (16%-84% percentile) are only somewhat larger for larger relative photometry uncertainty.Right panel: similar to the left panel but for 6-month-per-year WFST DHS observations across 6 years.Accumulating WFST DHS observations from 1 year to 6 years, not only do the dispersions of the retrieved timescales get smaller but also the minimal unbiased timescale retrieved decreases to ≃ 0.5 day.

Figure 11 .
Figure 11.Given a typical simulated LC of an OGLE-observed AGN with τin = T /10 = 292 day for T = 8 years and σin = 0.14 mag (open star for the input parameters), the likelihood surfaces (i.e., red, green, and blue contours for ∆ ln L = −0.5, -2.0, and -4.5 or, equivalently, 1, 2, and 3σ confidence levels) implied by the celerite (solid contours) and PRH (dotted contours) likelihood functions are compared together with their maximum likelihood points, i.e., filled and open circles, respectively.

Figure 12 .
Figure 12.Left panel: same as the top-left panel of Figure 3, but only the variation timescale retrieved by the K17PMm solution (green solid line) is shown and compared to that retrieved including σ * (grey dot-dash line).Middle panel: same as the left panel but for the retrieved variation amplitude σ.Right panel: adopting the same K17PMm solution, illustrated is the retrieved σ * relative to the average OGLE photometric uncertainty (i.e., ⟨σOGLE⟩ ≃ 0.025 mag) as a function of T /τin.

Table 1 .
Settings of the three fitting methods and the corresponding performances