How Long Will the Quasar UV/Optical Flickering Be Damped?

The UV/optical light curves of Active Galactic Nuclei (AGNs) are commonly described by the Damped Random Walk (DRW) model. However, the physical interpretation of the damping timescale, a key parameter in the DRW model, remains unclear. Particularly, recent observations indicate a weak dependence of the damping timescale upon both wavelength and accretion rate, clearly being inconsistent with the accretion-disk theory. In this study, we investigate the damping timescale in the framework of the Corona Heated Accretion disk Reprocessing (CHAR) model, a physical model that describes AGN variability. We find that while the CHAR model can reproduce the observed power spectral densities of the 20 yr light curves for 190 sources from Stone et al., the observed damping timescale, as well as its weak dependence on wavelength, can also be well recovered through fitting the mock light curves with DRW. We further demonstrate that such weak dependence is artificial due to the effect of inadequate durations of light curves, which leads to best-fitting damping timescales lower than the intrinsic ones. After eliminating this effect, the CHAR model indeed yields a strong dependence of the intrinsic damping timescale on the bolometric luminosity and rest-frame wavelength. Our results highlight the demand for sufficiently long light curves in AGN variability studies and important applications of the CHAR model in such studies.


INTRODUCTION
At the centers of Active Galactic Nuclei (AGNs), growing supermassive black holes (SMBHs) are surrounded by gaseous structures such as accretion disks, broad-line regions, and dust tori.It is widely believed that most of the intrinsic properties of AGNs stem from the accretion effect of SMBHs, such as high luminosities and intense electromagnetic radiations in the entire electromagnetic spectrum.The emission of AGNs exhibits a complex and significant stochastic variability over a wide range of wavelengths, e.g., in the ultraviolet (UV) and optical bands (e.g., Ulrich et al. 1997).This inherent variability provides a unique way to understand AGNs' structures and intrinsic physical mech-anisms.It has been shown that there are relationships between the AGN UV/optical fractional variability amplitude and the AGN luminosity (e.g., MacLeod et al. 2010;Zuo et al. 2012;Morganson et al. 2014;Li et al. 2018;Sun et al. 2018;Suberlak et al. 2021), Eddington ratio (e.g., MacLeod et al. 2010;Morganson et al. 2014;Simm et al. 2016;Sun et al. 2018;De Cicco et al. 2022), rest-frame wavelength (e.g., MacLeod et al. 2010;Morganson et al. 2014;Sun et al. 2015;Simm et al. 2016;Sánchez-Sáez et al. 2018;Suberlak et al. 2021), and other relevant parameters (e.g., Ai et al. 2010;MacLeod et al. 2010;Sun et al. 2018;Kang et al. 2018).By measuring the time delay between the variability of different bands using the reverberation mapping technique (Blandford & McKee 1982), one can measure the size of the AGN accretion disk and the broad-line region, helping us estimate the virial mass of the central SMBH and test the accretion-disk theory (e.g., Fausnaugh et al. 2016;Du & Wang 2019;Cackett et al. 2021).

Zhou et al.
Although much progress has been achieved in the studies of AGN variability, the primary physical mechanism causing variability is still an unresolved issue.Many theories have been proposed, such as the global variation of the accretion rate (Lyubarskii 1997;Li & Cao 2008;Liu et al. 2016), the local temperature fluctuation (Kelly et al. 2009;Dexter & Agol 2011;Cai et al. 2016), and the effect of large-scale fluctuations on local temperature fluctuations (Cai et al. 2018;Neustadt & Kochanek 2022;Secunda et al. 2023).However, these theories failed to fully account for observational results, e.g., AGN power spectral densities (PSDs).More sophisticated theoretical models and better observational data are required to investigate AGN variability further.
The Damped Random Walk (DRW) model is an effective statistical model for fitting AGN light curves (e.g., Kelly et al. 2009;Koz lowski et al. 2010;MacLeod et al. 2010;Zu et al. 2013;Suberlak et al. 2021).Its power spectral density is described by a power-law function f −2 at high frequencies, which transits to white noise at low frequencies.The transition frequency is f 0 = 1/(2πτ DRW ), where τ DRW is the damping timescale that characterizes the variability.Then, τ DRW should be relevant to characteristic timescales of the accretion disk, such as the thermal timescale.According to the static standard accretion disk theory (hereafter SSD; Shakura & Sunyaev 1973), for a given wavelength λ, if one assumes the corresponding emission-region size R λ merits the criteria k B T (R λ ) = hc/λ, the emissionregion size R λ /R S ∝ M −1/3 BH ṁ1/3 λ 4/3 and the thermal timescale τ th ∼ α −1 Ω −1 k ∝ α −1 λ 2 Ṁ 0.5 , where k B , T , h, c, R S ≡ 2GM BH /c 2 , α, and Ω k are the Boltzmann constant, the effective temperature, the Planck constant, the speed of light, the Schwarzschild radius, the viscosity parameter and the Keplerian angular velocity, respectively; ṁ = Ṁ / ṀEdd = 0.1c 2 Ṁ /L Edd is a dimensionless accretion rate, which is the ratio of the accretion rate Ṁ to the Eddington accretion rate ṀEdd , and L Edd is the Eddington luminosity.If τ DRW is relevant to the thermal timescale, τ DRW should also scale as λ 2 or Ṁ 0.5 .
Several studies have established empirical relationships between the best-fitting τ DRW and λ and the physical properties of AGNs.MacLeod et al. (2010) obtained τ DRW ∝ λ 0.17 , τ DRW ∝ Ṁ −0.075 , and τ DRW ∝ M 0.21 BH based on 10-year variability data for 9000 quasars in SSDS Stripe 82; Suberlak et al. (2021) obtained τ DRW ∝ Ṁ −0.088 and τ DRW ∝ M 0.14 BH based on 15-year variability data for 9000 quasars in SSDS Stripe 82; Burke et al. (2021) obtained τ DRW ∝ M 0.38  BH and τ DRW ∝ Ṁ 0.33 based on the variability data of 67 AGNs, and they even proposed to use the τ DRW − M BH relation to estimate SMBH masses; very recently, Stone et al. (2022) (hereafter S22; also see Stone et al. 2023 for erratum) obtained τ DRW ∝ λ 0.20 based on 20-year variability data for 190 quasars in SSDS Stripe 82 (Z.Stone, private communication).These observations demonstrate a weak dependence between the best-fitting τ DRW and λ and Ṁ (but see Sun et al. 2018;Arevalo et al. 2023).One possible explanation is that the relationship between the damping timescale and the thermal timescale is not a simple linear one.In summary, these observations seem to disfavor the static SSD model strongly.
It is worth pointing out that when using the DRW model to fit light curves, the durations of light curves (referred to as the baselines) must be much longer than the intrinsic damping timescale; otherwise, the damping timescale will be significantly underestimated.Koz lowski (2017) pointed out that if the intrinsic damping timescale is larger than 10% of the baseline, the lowfrequency white noise cannot be identified accurately, leading to an underestimation of the damping timescale.Suberlak et al. (2021) followed the simulation methodologies of Koz lowski (2017) and revisited the criteria for baselines to obtain unbiased damping timescales; they claim that the best-fitting damping timescale is an unbiased estimator of the intrinsic damping timescale if the best-fitting one is less than 20% of the baseline.Recently, Koz lowski (2021) stressed that the baseline should be at least 30 times the intrinsic damping timescale.Hu et al. (2023) further suggested that the deviation from the best-fitting damping timescale to the intrinsic value also depends upon the statistical assumptions, including the priors for the DRW parameters and the estimators of the best-fitting damping timescale.In observational studies, it is often argued that if the bestfitting damping timescale is less than 10% (or 20%) of the baseline, the best-fitting damping timescale is unbiased (e.g., Suberlak et al. 2021;Burke et al. 2021).We stress that this criteria is incorrect.Even if the intrinsic damping timescale is larger than 10% of the baseline, the best-fitting damping timescale can still be much smaller than 10% of the baseline in some random realizations of AGN variability (especially for observations with poor cadences).To verify whether the best-fitting damping timescale is unbiased, one must know the intrinsic damping timescale!S22 studied a sample of 190 quasars with a 20-year baseline in the observed frame using the DRW model.There are only 27 sources that merit the criteria of having a best-fitting damping timescale of less than 20% of the baseline.In addition, the best-fitting damping timescale also has large uncertainties.In our opinion, for the 27 sources, it is still unclear whether their intrinsic damping timescales are less than 20% of the baseline.
The same argument also holds for other observational studies (e.g., Suberlak et al. 2021;Burke et al. 2021).In summary, current observational studies of damping timescales may still be seriously limited by the baseline.
Theoretically, the AGN variability study should not be based on the static SSD model.Instead, timedependent physical models should be proposed to explain AGN variability.The Corona Heated Accretion disc Reprocessing (CHAR; Sun et al. 2020a) model established a new connection between disk fluctuations and observed variability; this model can reproduce the SDSS quasar variability (Sun et al. 2020b;Sun 2023) and the larger-than-expected UV/optical time lags (Li et al. 2021).In this model, the black hole accretion disk and the corona are coupled by magnetic fields.When the magnetic field in the corona perturbates, not only does the X-ray luminosity of the corona changes, but also the heating rate of the accretion disk fluctuates.Thus, the temperature of the accretion disk and the UV/optical luminosity vary on the thermal timescales.The CHAR model takes into account the time-dependent evolution of the SSD model to describe the temperature fluctuations, albeit it can be extended to include other disk models.The CHAR model may be effective for understanding the observational results of AGN variability.Is the CHAR model able to account for the observed damping timescales?With future wide-field time-domain surveys, e.g., the Wide Field Survey Telescope (WFST; Wang et al. 2023a) and the Legacy Survey of Space and Time (LSST; Ivezić et al. 2019), a large amount of AGN data covering multi-band time domains will be available.Extracting information beneficial to AGN studies from a large amount of variability data is crucial.In the future, longer AGN variability data will provide more accurate damping timescales.Hence, it is of great importance to understand the physical nature of the damping timescale.
The main objectives of this paper are to test the CHAR model with S22 observations, to explain the physical nature of the damping timescale in the DRW model, and to propose a new relation between the intrinsic damping timescale and the AGN properties.
The manuscript is organized as follows.In Section 2, the CHAR model is compared with the sample of S22; in Section 3, we offer new relationships between the intrinsic damping timescales and AGN properties in the CHAR model simulations; in Section 4, we study the PSD shapes; and Section 5 summarizes the main conclusions.

CONSISTENCY BETWEEN THE CHAR MODEL AND REAL OBSERVATIONS
In this section, we use the CHAR model to reproduce the DRW fitting results of S22.In Section 2.1, we configure the parameters of the CHAR model.In Sections 2.2 and 2.3, we use the CHAR model to reproduce the damping timescale dependence on wavelength and the PSDs of the S22 sample.

Setting the parameters for the CHAR model
The CHAR model uses the static SSD as the initial conditions.Hence, this model requires only three parameters to simulate the light curves: the dimensionless viscosity parameter α, the black hole mass M BH , and the dimensionless accretion ratio ṁ(= L bol /((1 + k)L Edd ))1 .For all simulations in this paper, we adopt α = 0.4 (King et al. 2007;Sun 2023).For the 190 quasars in the S22 sample, we use their M BH and the bolometric luminosity L bol as the input parameters of the CHAR model (the base parameterization).The black-hole mass is estimated via the single-epoch virial mass estimators, which have substantial uncertainties (∼ 0.5 dex; for a review, see, e.g., Shen 2013).We, therefore, also alter the black hole masses by 0.5 dex but leave other parameters (e.g., L bol ) unchanged, and the resulting damping timescales and their relation to wavelengths are almost unaltered.This is because the damping timescales of the CHAR model are independent of M BH for fixing L bol (see Section 3.2).There is an uncertainty of 0.2 − 0.3 dex in the observationally determined L bol (e.g., Netzer 2019).To assess the luminosity uncertainties, we consider an extreme case, in which L bol in the S22 sample are systematically reduced by 0.2 dex, leaving the other parameters unchanged (the "faint" parameterization).In summary, there is no free parameter in the simulations.
S22 uses the g, r, and i light curves.To facilitate calculations and comparisons, we use the CHAR model to calculate the observed frame (according to each source's redshift) 4500−5500 Å, 5500−6500 Å, and 7000−8000 Å emission to represent the g, r, and i bands, respectively.

Dependence of τ DRW on wavelength
The baseline of the light curves simulated by the CHAR model is 20 years in the observed frame (i.e., identical to S22).We adopt two sampling patterns for the light-curve simulations: uniform sampling with a cadence of 10 days (hereafter the uniform sampling) and irregular sampling with realistic cadences (hereafter the real sampling).For each quasar with the base and "faint" parameterizations, we use the CHAR model to generate light curves, encompassing both uniform and real sampling.We fit the light curves with the DRW model and derive the two DRW parameters (i.e., the damping timescale and the variability amplitude) using the taufit code of Burke et al. (2021) 2 which is based on the celerite Gaussian-process package (Foreman-Mackey et al. 2017).The celerite package fits light curves to a given kernel function using the Gaussian Process regression.The DRW kernel function in the taufit code is where t ij = |t i − t j | is the time interval between two measurements in the light curve, 2σ 2 DRW is the long-term variance of variability, τ DRW is the damping timescale, and σ 2 i δ ij denotes an excess white noise term from the measurement errors, with σ i being the excess white noise amplitude and δ ij being the Kronecker δ function.The taufit code utilizes the Markov Chain Monte Carlo (MCMC) code of emcee (Foreman-Mackey et al. 2013) with uniform priors to obtain the posterior distributions for the DRW parameters.We take the posterior medians as the best-fitting parameters and the 16 th to 84 th percentiles of the posterior as the 1σ uncertainties for the measured parameters.The simulation process is repeated 1, 000 times.Then, we take the medians and 16 th to 84 th percentiles of the one thousand best-fitting damping timescales obtained from the simulations, and the results are shown in Table 1.All simulation results are consistent with S22 within 1σ uncertainties.The sampling slightly affects the best-fitting τ DRW , and the effect increases with wavelengths.
We then investigate the dependence of the best-fitting τ DRW on wavelength obtained from the CHAR model and compare them with the results in S22.Previous studies have shown that the best-fitting τ DRW is significantly underestimated if the baseline is not longer than ten (or five) times the intrinsic damping timescale (Koz lowski 2017; Suberlak et al. 2021;Hu et al. 2023).S22 selects a subsample with the bestfitting damping timescales less than 20% of the baseline containing 27 sources.As mentioned in Section 1, their source selection criteria cannot eliminate the ef-fects of baseline inadequacy.The real observations in S22 yield τ DRW,S22 ∝ λ 0.20±0.20 for the subsample and τ DRW,S22 ∝ λ 0.30±0.13for the full sample.
, where σ CHAR is the 1σ uncertainty of the best-fitting τ DRW,CHAR .The best-fitting results for m and n are taken as the posterior medians, and their 1σ uncertainties are taken as 16 th to 84 th percentiles of the posterior distribution.
For different sampling methods and parameterizations, we obtain the probabilities of having the CHARmodel slope m agree with S22 within 1σ uncertainty through one thousand repeated CHAR model simulations, and the results are shown in Table 2. Thus, we cannot statistically reject the hypothesis that the dependence of the best-fitting τ DRW on wavelength obtained by the CHAR model is statistically consistent with S22.
Real observations show that the best-fitting τ DRW exhibits a weak dependence on wavelength.One possible reason is that the baseline is not long enough, leading to an underestimation of the intrinsic τ DRW (see Section 3.1 for details).We stress that, as we mentioned in Section 1, this bias still exists even if one only selects sources with the best-fitting damping timescale less than 10% (or 20%) of the baseline.Indeed, if we increase the simulation baseline, the best-fitting damping timescale also increases, and the timescale-wavelength relation has a larger slope.Another possible reason is that the observed τ DRW is the average of thermal timescales at different radii of the accretion disk, hence the relationship between τ DRW and wavelength may not be very simple (see Section 3.3 for details).

Power spectral density
Power spectral density is one of the most essential tools for studying AGN variability.We use the CHAR model to generate light curves for the full sample in S22 3.20 +0.17   4) are the CHAR model simulation results with uniform sampling for the base and "faint" parameterizations, respectively; Columns ( 5) and ( 6) are the CHAR model simulation results with real sampling.In summary, our studies above suggest that the CHAR model can reproduce the ensemble PSDs of the sample in S22.We then can use the CHAR model to predict the intrinsic damping timescale for other AGNs and its dependence upon AGN properties after properly eliminating the biases due to the limited baseline.
In the observational point of view, the DRW (a.k.a., the CARMA(1, 0) model) modeling is still an efficient way to understand AGN UV/optical variability.First, the real light curves are often very sparse.Several sophisticated statistical models are proposed to fit AGN UV/optical light curves, e.g., the damped harmonic oscillator (DHO, a.k.a., CARMA(2, 1); Moreno et al. 2019) model and other high-order CARMA models (Kelly et al. 2014).Compared with the DRW model, these sophisticated statistical models have more free parameters.For most AGN UV/optical light curves, these free parameters cannot be simultaneously well constrained.Second, the damping timescale of the DRW model is related to the timescales of the DHO model (see the lower-right panel of Figure 14 in Yu et al. 2022).Kasliwal et al. (2017) compared the PSDs of the DRW and the DHO models and found no difference between the two PSDs on timescales longer than ∼ 10 days.Hence, while the ensemble PSD analysis suggests that AGN UV/optical variability is not consistent with the DRW model (Figures 2, 3 and 4), one can still use this model to measure the long breaking timescales.We stress that, as pointed out by Vio et al. (1992), there is no one-to-one relationship between statistical models and intrinsic dynamics.Ideally, one should directly fit real AGN UV/optical light curves with physical models that have been tested (e.g., the CHAR model), which is beyond the scope of this manuscript.

PHYSICAL INTERPRETATION FOR THE DAMPING TIMESCALE
We now discuss the intrinsic damping timescale in the CHAR model.The simulation steps are introduced in Section 3.1; the relationship between the intrinsic damping timescale and M BH , ṁ and rest-frame wavelength λ rest is given in Section 3.2; the relationships between the intrinsic damping timescale and a series of physical properties are presented in Section 3.3.

Simulation steps
The simulation parameter settings and key points are as follows: 1.The simulation parameters M BH and ṁ used for the CHAR model are shown in Table 3, including 27 cases with the bolometric luminosity larger than 10 44 erg s −1 .
2. The simulation bands are shown in Table 4 which cover rest-frame UV-to-optical wavelengths.
3. The light curves of the integrated thermal emission from the whole disk with rest-frame 20-year and 40-year baselines are simulated using the CHAR model with a cadence of 10 days.The light curves are then fitted using the DRW model, and the best-fitting damping timescale is denoted as τ DRW,disk .The simulation is repeated unless at least one hundred sets of light curves are obtained for each case in Table 3.
4. For the rest-frame 20-year baseline, the light curves at different radii of the accretion disk are also fitted using the DRW model, and their best-fitting damping timescales are denoted as τ DRW,radius .
If the baseline is not longer than ten (or five) times the intrinsic damping timescale, τ DRW,disk will be strongly underestimated (e.g., Koz lowski 2017; Suberlak et al. 2021;Hu et al. 2023).Under such circumstances, the fitted τ DRW,disk should increase with the baseline.Figure 5 shows the relationship between the median τ DRW,disk over one hundred simulations and L bol at different bands, for both the 20-year and 40-year baselines.Note that, here we only consider simulations with τ DRW,disk /baseline < 20% for all bands.At the high bolometric luminosity end, there are clear differences between τ DRW,disk obtained for the 20-year baseline and the 40-year baseline, with the results for the 20-year baseline being significantly smaller than the 40-year baseline.The difference is more significant at longer wavelength bands.Hence, even if one only keeps the simulations with the best-fitting damping timescale less than 20% of the baseline, the intrinsic damping timescale cannot be obtained for the high luminosity cases.Likewise, in real observations, the best-fitting damping timescale is also biased even if it is less than 20% of the observed baseline.
As we mentioned in Section 1, the requirement for recovering the intrinsic damping timescale is that the intrinsic damping timescale (rather than the best-fitting one) should be less than 10% (or 20%) of the baseline.
The results of Hu et al. (2023) show that if the intrinsic τ DRW is 10% of the baseline, although the statistical expectation of the output best-fitting τ DRW is the same as the intrinsic τ DRW , its 1σ dispersion is as large as 50%.If one only requires that the best-fitting τ DRW is less than 10% (or 20%) of the baseline, the statistical expectation of τ DRW can still be significantly biased.This is because a DRW model with a large intrinsic τ DRW (larger than 10% or 20% of the baseline) can accidentally have a small best-fitting τ DRW because of statistical fluctuations.
If τ DRW,disk does not increase with baseline, one can regard τ DRW,disk as the intrinsic one.Practically speaking, the difference (∆log 10 τ DRW,disk ) between the 20year and 40-year baseline results is less than 0.1 at the longest wavelength range we probed ([7000, 8000] Å).Thus, for the 20-year baseline, τ DRW,disk is reliable (i.e., the best-fitting τ DRW,disk equal to the intrinsic value) only when L bol < 10 45.1 erg s −1 (cases marked with * in Table 3).For such sources, the best-fitting τ DRW,disk for the 20-year baseline is the same as that for the 40year baseline.For the L bol < 10 45.1 erg s −1 cases, their τ DRW,disk is less than 10% of baseline, which is consistent with Koz lowski (2017) and Hu et al. (2023).We only take the cases with L bol < 10 45.1 erg s −1 in Table 3 for subsequent analysis, but our conclusions are generalizable as long as the baseline is long enough.Previous studies have suggested that the best-fitting damping timescale may be related to the SMBH mass or the AGN luminosity and established empirical relationships between them based on observations (e.g., MacLeod et al. 2010;Burke et al. 2021;Wang et al. 2023b).However, real observations cannot properly eliminate the effects of inadequate baselines.We consider only the cases with L bol < 10 45.1 erg s −1 in Table 3 and use the medians of one hundred simulations of τ DRW,disk obtained from step 4 in Section 3.1 to establish the relationship between τ DRW,disk and M BH , ṁ, and rest-frame wavelength λ rest .λ rest is taken as the medians of the different bands in Table 4.The fitted equation is

The dependencies of τ
The MCMC code emcee (Foreman-Mackey et al. 2013) is adopted to sample the posterior distributions of the fitting parameters with the model likelihood and uniform priors.The logarithmic likelihood function is ln L = −0.5   .19, i.e., τ DRW,disk is strongly related to the bolometric luminosity L bol and the rest-frame wavelength λ rest , with little or no correlation with M BH , which ensures the feasibility of using the damping timescale to probe the cosmological time dilation (Lewis & Brewer 2023).
For the cases in Table 3 but with L bol > 10 45.1 erg s −1 (gray dots in Figure 6), their best-fitting damping timescales are strongly underestimated for the rest-frame 20-year baseline.Hence, they fall below the relation of Equation 2.
Our relationship (Equation 2) also holds for luminosity ranges lower than all cases in Table 3.To demonstrate this, we consider a low-luminosity case with M BH = 10 7.0 M ⊙ and ṁ = 0.01.We set a cadence of one day rather than ten days because the expected damping timescales (Equation 2) can be as short as ten days.Same as in Section 3.1, we only consider simulations with τ DRW,disk /baseline < 20% for all bands.The simulations are repeated 100 times.The purple-open dots in Figure 6 are the median values of the hundred simulations.These damping timescales are in good agreement with Equation 2.
In Equation 2, τ DRW,disk ∝ λ 1.19 , whereas the relationships between the best-fitting damping timescale and wavelength obtained in both Figure 1 and real observations of S22 are much weaker than this relationship.This is because, for most of the targets in S22, the 20-year baseline in the observed frame does not  2).In both panels, the purple-filled and grey dots represent the CHAR model prediction for cases in Table 3 with L bol < 10 45.1 erg s −1 (whose τ DRW,disk is unbiased for the 20-year baseline) and L bol > 10 45.1 erg s −1 (whose τ DRW,disk is strongly biased for the 20-year baseline), respectively.5. Darker colors indicate smaller τ DRW,cal -to-baseline ratios in the rest-frame.The lower panels show the relationship between ∆log10τDRW = log10(τ DRW,obs /τ DRW,cal ) and τ DRW,cal -to-baseline ratios in the rest-frame.∆log10τDRW decreases with τ DRW,cal -to-baseline ratios, which strongly suggests that the observationally-determined damping timescales are significantly underestimated.
yield unbiased damping timescales.For the subsample in S22, the selection criterion is that the best-fitting damping timescale rather than the intrinsic one is less than 20% of the baseline; hence, the best-fitting damping timescales in the subsample are also biased to smaller values.The bias should increase with wavelength if the intrinsic damping timescale positively correlates with the wavelength.As a result, the dependence of bestfitting damping timescale and wavelength obtained by S22 is weaker than Equation 2. Lewis & Brewer (2023) uses the damping timescales of the S22 sample to probe the cosmological time dilation.Their conclusion might also be affected by the same bias.The Equation 2 can predict a given AGN's intrinsic damping timescale and justify whether or not the best-fitting damping timescales in observational studies are biased.Figure 6 compares the best-fitting damping timescales from real observations (hereafter τ DRW,obs ) with Equation 2(hereafter τ DRW,cal ).The left panel presents the subsample of S22 with the bestfitting damping timescales less than 20% of the base-line.For almost all sources, the τ DRW,cal -to-baseline ratio is significantly larger than 10%, and the bestfitting damping timescales should be strongly biased.Hence, τ DRW,obs is always smaller than τ DRW,cal , just like the gray dots in Figure 6.Moreover, targets with smaller ∆log 10 τ DRW = log 10 (τ DRW,obs /τ DRW,cal ) values tend to have larger τ DRW,cal -to-baseline ratios.This anti-correlation again strongly supports the idea that the best-fitting damping timescales in S22 are probably underestimated.The right panel shows the 67 AGNs in Burke et al. (2021).Again, we find that, for sources with large τ DRW,cal -to-baseline ratios (i.e., > 10%), their τ DRW,obs are lower than the predictions from Equation 2. Interestingly, five sources4 in Burke et al. (2021) have small τ DRW,cal -to-baseline ratios (i.e., < 10%), and their best-fitting damping timescales agree well with Equation 2. Burke et al. (2021) combined the timescale mea- surements from AGNs and white dwarfs.The white dwarfs have short damping timescales of ≲ 0.01 day and can be unbiasedly measured.Then, they found that τ DRW ∝ M 0.5 BH .Given that the sources in their sample roughly have similar ṁ, their result also suggests that τ DRW ∝ L 0.5 bol , which is close to our Equation 2. To further test Equation 2, we measure the damping timescale for three local AGNs listed in Table 5.These sources have relatively long baselines and small luminosities or black-hole masses.According to Equation 2, these sources should have small intrinsic damping timescales and can be unbiasedly measured.Their light curves are obtained from literature as indicated in Table 5.We again use taufit to measure their bestfitting damping timescales and 1σ uncertainties (diamonds in the right panel of Figure 6).The results are again roughly consistent with Equation 2. Interestingly, two sources with relatively large deviations, i.e., NGC 4151 and NGC 3516, are both changing-look AGNs.
For NGC 3516, we separately measured their damping timescales in the high-flux state (from 1996 to 2007) and low-flux state (from 2018 to 2021).In the high-flux state, the measured damping timescale is consistent with our relation within 2σ; in the low-flux state, the measured damping timescale is almost identical to τ DRW,cal .Hence, our results suggest that the thermal structure of changing-look AGN's accretion disks changes as the line appears or disappears.
Having very long light curves is vital to obtain unbiased damping timescales.According to Equation 2, to obtain unbiased damping timescales, luminous targets with long-wavelength emission have expected baselines that are decades-long.In contrast, faint counterparts with short-wavelength emission have expected baselines of only a few days to a few years.Therefore, in the case of finite observational baselines, choosing targets with low bolometric luminosity or short rest-frame wavelength emission can improve the reliability of the bestfitting damping timescales.The LSST (Brandt et al. 2018) and WFST (Wang et al. 2023a) will provide a vast amount of AGN optical variability data in the southern and northern celestial hemispheres, respectively.These programs will extend the observational baselines and expand the variability data, helping one to obtain unbiased damping timescales.
Equation 2 also provides a new method for calculating the absolute accretion rate Ṁ , which ∝ M BH ṁ.While ensuring that the best-fitting damping timescale is un- Note-Column (1) is object designation; column (2) is the redshift; column (3) is the black-hole mass; column (4) is the bolometric luminosity; column (5) is the rest-frame wavelength; column (6) is the baseline in rest-frame; column (7) is the rest-frame best-fitting damping timescale τ DRW,obs from light curves; column (8) is τ DRW,obs -to-baseline ratio; column (9) is the damping timescale τ DRW,cal calculated by Equation 2; column (10) is τ DRW,cal -to-baseline ratio; column ( 11 (3)

Variability radius R var
If the damping timescale is related to the characteristic timescales of the accretion disk, such as the thermal timescale, then according to the static SSD model, the relationship between τ DRW and wavelength should be τ DRW ∝ λ 2 .The observations in S22 concluded τ DRW ∝ λ 0.20 , which is a biased result due to the baseline limitation.In Section 3.2, we obtain τ DRW ∝ λ 1.19 after eliminating the effect of baseline inadequacy by simulations, and the dependence of τ DRW on wavelength is still weaker than expected from the static SSD model.This is because the best-fitting damping timescale is an average of the radius-dependent characteristic timescales at different radii of the accretion disk.Figure 7 shows the 5100 Å flux variations of a given wavelength at different radii (hereafter, the single-radius light curves) of the same accretion disk.On short timescales, the observed variability is dominated by contributions from inner regions and vice versa.Hence, the best-fitting damping timescale is different from the thermal timescale at which k B T (R λ ) = hc/λ.
We are now interested in finding a characteristic radius (hereafter R var ) at which its single-radius light curve has a local damping timescale (i.e., τ DRW,radius ) equaling the damping timescale (i.e., τ DRW,disk ) of the integrated disk light curve.For each case in Table 3 with L bol < 10 45.1 erg s −1 , in step 4 in Section 3.1, we generate one hundred sets of single-radius light curves.We fit every single-radius light curve with the DRW model for each band and obtain the corresponding local damping timescale, τ DRW,radius .Figure 8  M BH = 10 8.0 M ⊙ and ṁ = 0.05.We define the variability radius R var to be the characteristic radius at which |τ DRW,radius − τ DRW,disk | < 0.05τ DRW,disk .Then, we find R var for each case and each band.

The dependencies of τ DRW,disk upon Rvar/RS
To investigate the relationship between τ DRW,disk and R var , we fit τ DRW,disk and R var at different bands for a specific M BH and ṁ using the following equation, log 10 τ DRW,disk = αlog 10 (R var /R S ) + β. (4) For each case in Table 3 with L bol < 10 45.1 erg s −1 , we repeat the fitting of Equation 4 one hundred times and adopt the medians as the parameters α and β corresponding to each case.Figure 9 represents a fitting   result for M BH = 10 8.0 M ⊙ and ṁ = 0.05.The slope differs from the scaling law of ∼ R 3/2 .We can obtain α and β for each case in Table 3 with L bol < 10 45.1 erg s −1 and find that they depend upon M BH and ṁ.Hence, we aim to find the relationships between the parameters α and β and M BH and ṁ (Fig-     For a given wavelength λ, the emission-region size at which k B T (R λ ) = hc/λ for the CHAR model is (Sun et al. 2020a) where k B , σ, G, η, and h denote the Boltzmann constant, the Stefan-Boltzmann constant, the gravitational constant, the radiation efficiency (fixed at 0.1), and the Planck constant, respectively.We try to find the relationship between R var and R λ .We fit R var /R S and R λ /R S using the following relation: Figure 12 shows the relationship between R var /R S and R λ /R S for M BH = 10 8.0 M ⊙ and ṁ = 0.05.The relationship between R var /R S and R λ /R S is non-linear since A is less than unity.Figure 13 shows the values of parameters A and B for the different cases.It is evident that A and B depend weakly upon M BH or ṁ.We cannot find a simple function to describe the relationship between A (or B) and M BH or ṁ.
3.3.4.The dependencies of τ DRW,disk upon R λ /RS We also aim to establish a relationship between the directly measurable quantity τ DRW,disk and R λ /R S .The fitted equation is Figure 14 is the relationship between τ DRW,disk and R λ /R S for M BH = 10 8.0 M ⊙ and ṁ = 0.05.Again, K 1 and K 2 depend upon M BH or ṁ (Figure 15), and the best-fitting relationships are K 1 = 0.47 (12) We stress that K 1 is a function of M BH and ṁ, rather than being fixed to 3/2 as expected from the thermal timescale.

POWER SPECTRAL DENSITIES AT DIFFERENT RADII
We also use the FFT method to calculate the PSDs for the light curves at each radius.Figure 16 shows the PSDs at different radii for a typycal case of M BH = 10 7.5 M ⊙ and ṁ = 0.2.It is obvious that the PSD is steeper than that of a DRW on short timescales.This is because, in the CHAR model, the temperature fluctuations have a PSD of ∝ f −3 on short timescales (Sun et al. 2020a).As a result, the PSD of each single-radius light curve is inconsistent with the DRW model.The PSD of the integrated disk light curve is a superposition of various non-DRW PSDs at various disk radii and can resemble the DRW PSD on timescales from months to years (see Figure 4 of Sun et al. 2020a).We compare the PSD corresponding to R var with that of the whole disk.At high frequencies, the PSDs corresponding to R var are smaller than that of the whole disk.The two PSDs are consistent at low frequencies.Hence, the damping timescales from the DRW fitting and PSD analysis should be generally consistent.The PSD at the characteristic radius R var can roughly represent the PSD of the whole disk.

SUMMARY
We have used the CHAR model to reproduce the observations of S22.In addition, we have obtained a new scaling relation for the intrinsic damping timescale and its connection to the AGN properties through CHAR model simulations.The main conclusions are as follows: 1.The CHAR model can reproduce the DRW fitting results for the S22 sample, including the bestfitting damping timescales (see Table 1), the dependence of the best-fitting damping timescale on wavelength (see Figure 1, Table 2, Section 2.2) and the ensemble PSDs of the sample (see Figures 2, 3, and 4; Section 2.3).
2. The observational baselines for most luminous AGNs are not long enough to recover the intrinsic damping timescale.The damping timescale may be biased even if the best-fitting damping timescale is less than 20% (or 10%) of the observed baseline (Figure 5, Section 3.1).
3. We have obtained the relationship between the intrinsic damping timescale and M BH , ṁ, λ rest after eliminating the effect of baseline inadequacy by simulations.The main factors affecting the intrinsic damping timescale are the bolometric luminosity and rest-frame wavelength (Equation 2, Figure 6, Section 3.2).Furthermore, we argue that our Zhou et al.   results can be used to estimate the absolute accretion rate (see Equation 3) based on the intrinsic damping timescale.
4. The observed light curves are a summation of variable emission from various disk radii, and the damping timescale should be some averages of radius-dependent characteristic timescales (see Figure 7).This leads to a weaker dependence between the damping timescale and wavelength than the static SSD model.We have obtained the relationship between the directly measurable quantity τ 5.The PSDs corresponding to R var are consistent with that of the whole disk at low frequencies; the damping timescales from DRW fitting and PSD analysis should be generally consistent (see Figure 16, Section 4).
The upcoming LSST and WFST are expected to provide tremendous AGN variability data that will enlarge the light-curve baselines significantly, further validating our results.

3
Figure 1 shows a typical realization of the dependence of bestfitting rest-frame τ DRW,CHAR on wavelength obtained from the CHAR model: the left panel includes only the subsample in S22, while the right panel shows the result for their full sample.We fit the τ DRW,CHAR − λ rest relation with log 10 τ DRW,CHAR = mlog 10 λ rest +n.The posterior distributions of m and n are obtained by the MCMC code emcee (Foreman-Mackey et al. 2013) with uniform priors, and the logarithmic likelihood function is ln L = −0.5 (log 10 τ DRW,CHAR −(m log 10 (1) represents the wavelength ranges in the observed frame; Column (2) is the results obtained by S22; Columns (3) and (

Figure 1 .
Figure1.A typical realization of the dependence of the best-fitting τDRW on wavelength obtained from the CHAR model.MBH, L bol , and redshifts are taken from the S22 sample.In each panel, the red error bar in the upper right corner shows the median 1σ uncertainty for all best-fitting τDRW.The Spearman's rank correlation coefficient (ρ), p-value, and the best-fitting slope (m) are shown in the lower right corner.The left panel is for the same subsample in S22 (i.e., with 27 quasars) and the correlation is statistically insignificant.The right panel is for the full sample of 190 quasars in S22; in this sample, the correlation is statistically significant.The solid lines and the shaded areas are the best-fitting lines and 1σ uncertainties.
Figures 3 and 4 show the ensemble PSDs of the subsamples grouped by M BH and the subsamples grouped by L bol and redshift in λ obs = 4500 − 5500 Å band, respectively.The grouping strategies are the same as Figures 16 and 18 of S22.The model ensemble PSDs of the subsamples exhibit similar behaviors to those of S22: the high-frequency breaks of quasars with smaller SMBHs or lower luminosities tend to occur on shorter timescales and vice versa.

Figure 2 .Figure 3 .
Figure 2. Comparisons of the CHAR model and real observations of rest-frame ensemble PSDs for the full sample in different bands.The purple curves are the real observations from S22.The blue, red, and yellow curves and shaded areas are the simulation results of the CHAR model and their 1σ confidence intervals.The dashed grey lines represent the f −2 power-law function.The lower panels are the differences between the CHAR model and the real observations on a logarithmic scale.The ensemble PSDs of the CHAR model are consistent with those of the S22 observations within the 1σ confidence intervals.

Figure 4 .
Figure 4.The ensemble PSDs of the CHAR model in the λ obs = 4500 − 5500 Å band for different L bol and redshift groups.The grouping strategy is the same as Figure 16 in S22.The dashed grey lines represent the f −2 power-law function.The dependence of the ensemble PSDs upon L bol and redshift is similar to that of S22.

Figure 5 .
Figure 5.The relationship between the best-fitting damping timescale of the integrated disk light curve τ DRW,disk and L bol .Different panels indicate different bands.The yellow inverted triangles and green squares are the medians of the one hundred simulations of the CHAR model for the 20-year and 40-year baselines, respectively.The error bars are 1σ uncertainties.The red dotted and red dashed lines represent 20% of the 20-year and 40-year baselines, respectively.The blue dotted lines represent 10% of the 20-year baseline.When L bol > 10 45.1 erg s −1 , τ DRW,disk is biased since its values for the 20-year baseline are significantly smaller than the 40-year baseline.andf model,i , and σ i are the measured damping timescale, the model timescale, the 1σ uncertainty of f i , respectively.The best-fitting values for the parameters are taken as the posterior medians, and their 1σ uncertainties are taken as 16 th to 84 th percentiles of the posterior distribution, as are all subsequent fits.Figure6shows the relationship between τ DRW,disk and M BH , ṁ, and λ rest (purple-filled dots).The bestfitting parameters and their 1σ uncertainties are a = 0.65 +0.01 −0.01 , b = 0.65 +0.01 −0.01 , c = 1.19 +0.01 −0.01 , d = −6.04+0.05 −0.05 .The result demonstrates that τ DRW,disk ∝ L 0.65 bol λ 1.19 , i.e., τ DRW,disk is strongly related to the bolometric luminosity L bol and the rest-frame wavelength λ rest , with little or no correlation with M BH , which ensures the feasibility of using the damping timescale to probe the cosmological time dilation(Lewis & Brewer 2023).For the cases in Table3but with L bol > 10 45.1 erg s −1 (gray dots in Figure6), their best-fitting damping timescales are strongly underestimated for the rest-

Figure 6 .
Figure6.The relationship between τ DRW,disk and MBH, ṁ, and λrest (see Equation2).In both panels, the purple-filled and grey dots represent the CHAR model prediction for cases in Table3with L bol < 10 45.1 erg s −1 (whose τ DRW,disk is unbiased for the 20-year baseline) and L bol > 10 45.1 erg s −1 (whose τ DRW,disk is strongly biased for the 20-year baseline), respectively.The purple-open dots are the CHAR model calculations for a low-luminosity case with MBH = 10 7.0 M⊙ and ṁ = 0.01.The gray dashed lines indicate the one-to-one relation.For comparison purposes, we also include real observations.The 27 sources in the left panel are taken from S22, where τ DRW,obs is less than 20% of the baseline.Dots of different colors indicate different bands.In the right panel, the orange dots indicate the 67 sources inBurke et al. (2021), and the color diamonds are model results for the cases in Table5.Darker colors indicate smaller τ DRW,cal -to-baseline ratios in the rest-frame.The lower panels show the relationship between ∆log10τDRW = log10(τ DRW,obs /τ DRW,cal ) and τ DRW,cal -to-baseline ratios in the rest-frame.∆log10τDRW decreases with τ DRW,cal -to-baseline ratios, which strongly suggests that the observationally-determined damping timescales are significantly underestimated.
Figure6.The relationship between τ DRW,disk and MBH, ṁ, and λrest (see Equation2).In both panels, the purple-filled and grey dots represent the CHAR model prediction for cases in Table3with L bol < 10 45.1 erg s −1 (whose τ DRW,disk is unbiased for the 20-year baseline) and L bol > 10 45.1 erg s −1 (whose τ DRW,disk is strongly biased for the 20-year baseline), respectively.The purple-open dots are the CHAR model calculations for a low-luminosity case with MBH = 10 7.0 M⊙ and ṁ = 0.01.The gray dashed lines indicate the one-to-one relation.For comparison purposes, we also include real observations.The 27 sources in the left panel are taken from S22, where τ DRW,obs is less than 20% of the baseline.Dots of different colors indicate different bands.In the right panel, the orange dots indicate the 67 sources inBurke et al. (2021), and the color diamonds are model results for the cases in Table5.Darker colors indicate smaller τ DRW,cal -to-baseline ratios in the rest-frame.The lower panels show the relationship between ∆log10τDRW = log10(τ DRW,obs /τ DRW,cal ) and τ DRW,cal -to-baseline ratios in the rest-frame.∆log10τDRW decreases with τ DRW,cal -to-baseline ratios, which strongly suggests that the observationally-determined damping timescales are significantly underestimated.

Figure 7 .
Figure7.The 5100 Å flux variations at different radii of the same accretion disk.The x-axis is time, and the y-axis is the ratio of the luminosity at a given radius to the luminosity of the whole disk at the wavelength of 5100 Å. Darker colors indicate larger radii.Smaller-radii light curves contain more short-term variations than larger-radii ones.

Figure 8 .
Figure 8. Relation between τ DRW,radius and R/RS for MBH = 10 8.0 M⊙ and ṁ = 0.05.The six panels represent six bands.Red dashed lines indicate 20% of the baseline.Horizontal gray lines and shaded gray areas are τ DRW,disk and its narrow range, which is 10% of τ DRW,disk .Vertical lines and shaded areas indicate a narrow radius range in which τ DRW,radius is the same as τ DRW,disk , and the average of this range is noted as Rvar.

Figure 9 .
Figure 9.The relationship between τ DRW,disk and Rvar for MBH = 10 8.0 M⊙ and ṁ = 0.05.The solid line is the bestfitting line and the shaded area is 1σ confidence intervals.

Figure 10 .
Figure10.Fitting results (data points) for parameters α and β in Equation4.Different panels represent different ṁ.We fit all data points with two linear relations for α and β, respectively.The light-and dark-shaded regions represent the 1σ and 2σ confidence intervals of the best-fitting linear relations, respectively.

Figure 12 .
Figure 12.The relationship between Rvar/RS and R λ /RS for MBH = 10 8.0 M⊙ and ṁ = 0.05.The dots represent our calculations.The solid line and shaded areas are the bestfitting result and 1σ confidence intervals, respectively.

Figure 13 .
Figure 13.Parameters A and B for different MBH and ṁ in Equation 9. Different colors represent different ṁ.

Figure 14 .
Figure 14.The relationship between τ DRW,disk and R λ /RS for MBH = 10 8.0 M⊙ and ṁ = 0.05.The solid line and shaded areas are the best-fitting relation and 1σ confidence intervals, respectively.

Table 1 .
The best-fitting logarithmic τDRW of S22 observations and the CAHR model simulations.

Table 2 .
(Mushotzky et al. 2011;Kasliwal et al. 2015;Smith et al. 2018) the ensemble PSDs of both the S22 observations and the CHAR model simulations are steeper than the f −2 power-law function, suggesting that the DRW model overestimates the variability on short timescales, which is consistent with previous Kepler studies(Mushotzky et al. 2011;Kasliwal et al. 2015;Smith et al. 2018).As in S22, we also study the model ensemble PSDs after grouping the full sample by different methods.
rest-frame ensemble PSDs for the full sample in different bands.The ensemble PSDs of the CHAR model are consistent with S22 observations within the 1σ confidence interval.
DRW,disk upon M BH , ṁ, and λ rest Zhou et al.

Table 3 .
Model parameters for the CHAR model simulation

Table 5 .
Additional low-luminosity targets with long light curves.
Burke et al. 2021)In previous studies (e.g.,Burke et al. 2021), it is often argued that the damping timescale should be related to the thermal timescale at R λ , which scales as R The relationship between Rvar/RS and R λ /RS Figure 11.The relationship between Rvar and MBH, ṁ and λrest.The gray dashed line indicates the one-to-one relation.