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

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 $\dot{m}(={L}_{\mathrm{bol}}/((1+k){L}_{\mathrm{Edd}}))$. ⁴ 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.

The baseline of the light curves simulated by the CHAR model is 20 yr 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) ⁵ 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

$$\begin{eqnarray}&&k({t}_{{ij}})=2{\sigma }_{\mathrm{DRW}}^{2}{e}^{{t}_{{ij}}/{\tau }_{\mathrm{DRW}}}+{\sigma }_{i}^{2}{\delta }_{{ij}},\end{eqnarray} \tag{ 1 }$$

where t_ij = ∣t_i − t_j ∣ is the time interval between two measurements in the light curve, $2{\sigma }_{\mathrm{DRW}}^{2}$ is the long-term variance of variability, τ_DRW is the damping timescale, and ${\sigma }_{i}^{2}{\delta }_{{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 16th–84th percentiles of the posterior as the 1σ uncertainties for the measured parameters. The simulation process is repeated 1000 times. Then, we take the medians and 16th–84th percentiles of the 1000 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.

Table 1. The Best-fitting Logarithmic τ_DRW of S22 Observations and the CAHR Model Simulations

Bands	Observations	Uniform sampling		Real sampling
		base	"faint"	base	"faint"
(1)	(2)	(3)	(4)	(5)	(6)
4500–5500 Å (g band)	${2.88}_{-0.28}^{+0.56}$	${3.20}_{-0.29}^{+0.17}$	${3.13}_{-0.30}^{+0.18}$	${3.19}_{-0.15}^{+0.11}$	${3.15}_{-0.18}^{+0.12}$
5500–6500 Å (r band)	${3.05}_{-0.34}^{+0.61}$	${3.33}_{-0.26}^{+0.16}$	${3.25}_{-0.30}^{+0.17}$	${3.23}_{-0.13}^{+0.10}$	${3.20}_{-0.14}^{+0.11}$
7000–8000 Å (i band)	${3.11}_{-0.35}^{+0.58}$	${3.47}_{-0.25}^{+0.16}$	${3.40}_{-0.25}^{+0.16}$	${3.30}_{-0.10}^{+0.08}$	${3.27}_{-0.12}^{+0.09}$

Note. Column (1) represents the wavelength ranges in the observed frame; column (2) is the results obtained by S22; columns (3) and (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.

Download table as:  ASCIITypeset image

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 10 (or five) times the intrinsic damping timescale (Kozłowski 2017; Suberlak et al. 2021; Hu et al. 2024). S22 selects a subsample with the best-fitting damping timescales <20% of the baseline containing 27 sources. As mentioned in Section 1, their source-selection criteria cannot eliminate the effects of baseline inadequacy. The real observations in S22 yield τ_DRW,S22 ∝ λ^0.20±0.20 for the subsample and τ_DRW,S22 ∝ λ^0.30±0.13 for the full sample. ⁶ Figure 1 shows a typical realization of the dependence of best-fitting 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 ${\mathrm{log}}_{10}{\tau }_{\mathrm{DRW},\mathrm{CHAR}}=m{\mathrm{log}}_{10}{\lambda }_{\mathrm{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 $\mathrm{ln}{ \mathcal L }=-0.5\sum \{({\mathrm{log}}_{10}{\tau }_{\mathrm{DRW,CHAR}}$ $-{(m{\mathrm{log}}_{10}{\lambda }_{\mathrm{rest}}+n))}^{2}/{\sigma }_{\mathrm{CHAR}}^{2}+\mathrm{ln}{\sigma }_{\mathrm{CHAR}}^{2}\}$, 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 16th–84th percentiles of the posterior distribution.

Zoom In Zoom Out Reset image size

For different sampling methods and parameterizations, we obtain the probabilities of having the CHAR model slope m agree with S22 within 1σ uncertainty through 1000 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 <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 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 with a 20 yr baseline in the observed frame and a cadence of 1 day. Because the light curves generated by the CHAR model are uniformly sampled, we use the Fast Fourier Transform (FFT) to obtain the model ensemble PSDs for the full sample. We repeat the above process 200 times and then take the average of these simulations to obtain the PSDs. Whereas the light curves of the observations are not uniformly sampled, S22 used the continuous autoregressive with moving-average (CARMA; Kelly et al. 2014) model, which is a high-order model, to obtain the ensemble PSDs. Figure 2 compares the CHAR model with real observations of 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. At high frequencies, the ensemble PSDs of both the S22 observations and the CHAR model simulations are steeper than the f⁻² 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).

Zoom In Zoom Out Reset image size

As in S22, we also study the model ensemble PSDs after grouping the full sample by different methods. 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In summary, our studies above suggest that the CHAR model can reproduce the ensemble PSDs of the sample in S22. We can then 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.

The simulation parameter settings and key points are as follows:

• 1.  
  The simulation parameters M_BH and $\dot{m}$ used for the CHAR model are shown in Table 3, including 27 cases with the bolometric luminosity >10⁴⁴ erg s⁻¹.
• 2.  
  The simulation bands that cover rest-frame UV-to-optical wavelengths are shown in Table 4.
• 3.  
  The light curves of the integrated thermal emission from the whole disk with rest-frame 20 and 40 yr 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 100 sets of light curves are obtained for each case in Table 3.
• 4.  
  For the rest-frame 20 yr 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.

Table 3. Model Parameters for the CHAR Model Simulation

No.	${\mathrm{log}}_{10}({M}_{\mathrm{BH}}/{M}_{\odot })$	$\dot{m}$	${\mathrm{log}}_{10}({L}_{\mathrm{bol}}/(\mathrm{erg}\,{{\rm{s}}}^{-1}))$
1*	8.0	0.01	44.24
2*	8.5	0.01	44.74
3	9.0	0.01	45.24
4*	7.5	0.05	44.44
5*	8.0	0.05	44.94
6	8.5	0.05	45.44
7	9.0	0.05	45.94
8*	7.0	0.1	44.24
9*	7.5	0.1	44.74
10	8.0	0.1	45.24
11	8.5	0.1	45.74
12	9.0	0.1	46.24
13*	7.0	0.15	44.41
14*	7.5	0.15	44.91
15	8.0	0.15	45.41
16	8.5	0.15	45.91
17	9.0	0.15	46.41
18*	7.0	0.2	44.54
19*	7.5	0.2	45.04
20	8.0	0.2	45.54
21	8.5	0.2	46.04
22	9.0	0.2	46.54
23*	7.0	0.5	44.94
24	7.5	0.5	45.44
25	8.0	0.5	45.94
26	8.5	0.5	46.44
27	9.0	0.5	46.94

Note. Column (1) shows the case number, those with * are cases with L_bol < 10^45.1 erg s⁻¹; column (2) represents the logarithmic black hole mass; column (3) indicates the dimensionless accretion ratio $\dot{m};$ and column (4) shows the logarithmic bolometric luminosity.

Download table as:  ASCIITypeset image

If the baseline is not longer than 10 (or five) times the intrinsic damping timescale, τ_DRW,disk will be strongly underestimated (e.g., Kozłowski 2017; Suberlak et al. 2021; Hu et al. 2024). Under such circumstances, the fitted τ_DRW,disk should increase with the baseline. Figure 5 shows the relationship between the median τ_DRW,disk over 100 simulations and L_bol at different bands, for both the 20 yr and 40 yr 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 yr baseline and the 40 yr baseline, with the results for the 20 yr baseline being significantly smaller than the 40 yr baseline. The difference is more significant at longer wavelength bands. Hence, even if one only keeps the simulations with the best-fitting damping timescale <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 <20% of the observed baseline.

Zoom In Zoom Out Reset image size

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 <10% (or 20%) of the baseline. The results of Hu et al. (2024) 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 <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 (>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 (${\rm{\Delta }}{\mathrm{log}}_{10}{\tau }_{\mathrm{DRW},\mathrm{disk}}$) between the 20 and 40 yr baseline results is <0.1 at the longest wavelength range we probed ([7000, 8000] Å). Thus, for the 20 yr 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⁻¹ (cases marked with * in Table 3). For such sources, the best-fitting τ_DRW,disk for the 20 yr baseline is the same as that for the 40 yr baseline. For the L_bol < 10^45.1 erg s⁻¹ cases, their τ_DRW,disk is <10% of baseline, which is consistent with Kozłowski (2017) and Hu et al. (2024). We only take the cases with L_bol < 10^45.1 erg s⁻¹ 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⁻¹ in Table 3 and use the medians of 100 simulations of τ_DRW,disk obtained from step 4 in Section 3.1 to establish the relationship between τ_DRW,disk and M_BH, $\dot{m}$, and rest-frame wavelength λ_rest. The parameter λ_rest is taken as the medians of the different bands in Table 4. The fitted equation is

$$\begin{eqnarray}\begin{array}{rcl}{\mathrm{log}}_{10}{\tau }_{\mathrm{DRW},\mathrm{disk}} & = & a{\mathrm{log}}_{10}({M}_{\mathrm{BH}}/{M}_{\odot })+b{\mathrm{log}}_{10}\dot{m}\\ & & +c{\mathrm{log}}_{10}({\lambda }_{\mathrm{rest}}/{\rm{\mathring{\rm A} }})+d,\end{array}\end{eqnarray} \tag{ 2 }$$

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 $\mathrm{ln}{ \mathcal L }\,=-0.5\sum \left\{{({f}_{{\rm{i}}}-{f}_{\mathrm{model},{\rm{i}}})}^{2}/{\sigma }_{{\rm{i}}}^{2}+\mathrm{ln}{\sigma }_{{\rm{i}}}^{2}\right\}$, where f_i, f_model,i, and σ_i are the measured damping timescale, the model timescale, and 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 16th–84th percentiles of the posterior distribution, as are all subsequent fits. Figure 6 shows the relationship between τ_DRW,disk and M_BH, $\dot{m}$, and λ_rest (purple-filled dots). The best-fitting 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 ${\tau }_{\mathrm{DRW},\mathrm{disk}}\propto {L}_{\mathrm{bol}}^{0.65}{\lambda }^{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).

Zoom In Zoom Out Reset image size

For the cases in Table 3 but with L_bol > 10^45.1 erg s⁻¹ (gray dots in Figure 6), their best-fitting damping timescales are strongly underestimated for the rest-frame 20 yr 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 $\dot{m}=0.01$. We set a cadence of one day rather than 10 days because the expected damping timescales (Equation (2)) can be as short as 10 days. 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 100 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 yr baseline in the observed frame does not 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 <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 best-fitting 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.

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 best-fitting damping timescales <20% of the baseline. For almost all sources, the τ_DRW,cal-to-baseline ratio is significantly >10%, and the best-fitting 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 ${\rm{\Delta }}{\mathrm{log}}_{10}{\tau }_{\mathrm{DRW}}={\mathrm{log}}_{10}({\tau }_{\mathrm{DRW},\mathrm{obs}}/{\tau }_{\mathrm{DRW},\mathrm{cal}})$ values tend to have larger τ_DRW,cal-to-baseline ratios. This anticorrelation 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 sources ⁷ 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 measurements 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 ${\tau }_{\mathrm{DRW}}\propto {M}_{\mathrm{BH}}^{0.5}$. Given that the sources in their sample roughly have similar $\dot{m}$, their result also suggests that ${\tau }_{\mathrm{DRW}}\propto {L}_{\mathrm{bol}}^{0.5}$, 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 best-fitting 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.

Table 5. Additional Low-luminosity Targets with Long Light Curves

Object	z	MBH	Lbol	λrest	Baseline	${\mathrm{log}}_{10}({\tau }_{\mathrm{DRW},\mathrm{obs}}/[\mathrm{days}])$	$\tfrac{{\tau }_{\mathrm{DRW},\mathrm{obs}}}{\mathrm{baseline}}$	${\mathrm{log}}_{10}({\tau }_{\mathrm{DRW},\mathrm{cal}}/[\mathrm{days}])$	$\tfrac{{\tau }_{\mathrm{DRW},\mathrm{cal}}}{\mathrm{baseline}}$	References
		(107 M⊙)	(1044 erg s−1)	(Å)	(day)		(%)		(%)
NGC 4151 a	0.0032	5.10	1.15 ± 0.13	4754	8118	${2.73}_{-0.10}^{+0.26}$	6.64	2.12	1.62	(1)
NGC 7469	0.0163	1.10	4.82 ± 0.74	5100	6868	${2.40}_{-0.16}^{+0.22}$	3.66	2.56	5.31	(2)
NGC 3516 b	0.0088	4.73	1.2	5100	4112	${2.62}_{-0.26}^{+0.42}$	5.25	2.17	3.58	(3, 4)
NGC 3516 c	0.0088	4.73	0.27	4728	1085	${1.75}_{-0.12}^{+0.16}$	5.15	1.71	4.69	(4, 5)

Notes. 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; and column (11) is the references: Ref. (1) Chen et al. (2023), Ref. (2) Shapovalova et al. (2017), Ref. (3) Shapovalova et al. (2019), Ref. (4) Mehdipour et al. (2022), Ref. (5) The g-band data from Zwicky Transient Facility (ZTF; Masci et al. 2019).

^a NGC 4151 is a changing-look AGN, and L_bol is taken from the high-flux state. We remove data points with the signal-to-noise ratios <20. ^b NGC 3516 is also a changing-look AGN, and here are the results in the high-flux state (from 1996 to 2007; Popović et al. 2023). ^c The results of NGC 3516 in the low-flux state (from 2018 to 2021; Popović et al. 2023).

Download table as:  ASCIITypeset image

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 best-fitting 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 $\dot{M}$, which $\propto {M}_{\mathrm{BH}}\dot{m}$. While ensuring that the best-fitting damping timescale is unbiased, $\dot{M}$ should satisfy the following equation,

$$\begin{eqnarray}&&\dot{M}=43.51{\tau }_{\mathrm{DRW},\mathrm{disk}}^{1.54}{\lambda }_{\mathrm{rest}}^{-1.83}\ [{M}_{\odot }{\mathrm{yr}}^{-1}].\end{eqnarray} \tag{ 3 }$$

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 ∝ λ². 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/λ.

Zoom In Zoom Out Reset image size

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⁻¹, in step 4 in Section 3.1, we generate 100 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 represents the relation between τ_DRW,radius and R/R_S for M_BH = 10^8.0 M_⊙ and $\dot{m}=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.

Zoom In Zoom Out Reset image size

3.3.1. The Dependencies of τ_DRW,disk upon R_var/R_S

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 $\dot{m}$ using the following equation:

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{\tau }_{\mathrm{DRW},\mathrm{disk}}=\alpha {\mathrm{log}}_{10}({R}_{\mathrm{var}}/{R}_{{\rm{S}}})+\beta .\end{eqnarray} \tag{ 4 }$$

For each case in Table 3 with L_bol < 10^45.1 erg s⁻¹, we repeat the fitting of Equation (4) 100 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 $\dot{m}=0.05$. The slope differs from the scaling law of ∼R^3/2.

Zoom In Zoom Out Reset image size

We can obtain α and β for each case in Table 3 with L_bol < 10^45.1 erg s⁻¹ and find that they depend upon M_BH and $\dot{m}$. Hence, we aim to find the relationships between the parameters α and β and M_BH and $\dot{m}$ (Figure 10) by simultaneously fitting these cases. The best-fitting results are

$$\begin{eqnarray}&&\alpha ={0.49}_{-0.01}^{+0.01}{\mathrm{log}}_{10}\displaystyle \frac{{M}_{\mathrm{BH}}}{{M}_{\odot }}+{0.52}_{-0.02}^{+0.02}{\mathrm{log}}_{10}\dot{m}-{2.26}_{-0.10}^{+0.10},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\beta =\left\{\begin{array}{l}{0.51}_{-0.03}^{+0.03}{\mathrm{log}}_{10}\displaystyle \frac{{M}_{\mathrm{BH}}}{{M}_{\odot }}-{0.11}_{-0.04}^{+0.04}{\mathrm{log}}_{10}\dot{m}-{3.73}_{-0.22}^{+0.22}\\ \dot{m}\leqslant 0.1,\\ -{0.04}_{-0.05}^{+0.05}{\mathrm{log}}_{10}\displaystyle \frac{{M}_{\mathrm{BH}}}{{M}_{\odot }}-{0.85}_{-0.08}^{+0.08}{\mathrm{log}}_{10}\dot{m}-{0.33}_{-0.36}^{+0.36}\\ \dot{m}\gt 0.1.\end{array}\right.\end{eqnarray} \tag{ 6 }$$

In summary, the relation between the damping timescale and R_var depends upon the black hole mass M_BH and dimensionless accretion ratio $\dot{m}$.

Zoom In Zoom Out Reset image size

3.3.2. The Dependencies of R_var/R_S upon M_BH, $\dot{m}$, and λ_rest

We also establish the relationship between R_var/R_S and M_BH, $\dot{m}$, and rest-frame wavelength λ_rest. The fitting equation is

$$\begin{eqnarray}\begin{array}{rcl}{\mathrm{log}}_{10}({R}_{\mathrm{var}}/{R}_{{\rm{S}}}) & = & u{\mathrm{log}}_{10}({M}_{\mathrm{BH}}/{M}_{\odot })+v{\mathrm{log}}_{10}\dot{m}\\ & & +s{\mathrm{log}}_{10}({\lambda }_{\mathrm{rest}}/{\rm{\mathring{\rm A} }})+\gamma .\end{array}\end{eqnarray} \tag{ 7 }$$

Figure 11 illustrates the fitting results of Equation (7). The best-fitting parameters are $u=-{0.63}_{-0.01}^{+0.01},v={0.01}_{-0.01}^{+0.01},s={1.06}_{-0.01}^{+0.01},\gamma ={2.99}_{-0.06}^{+0.06}.$ According to Equation (7), ${R}_{\mathrm{var}}\sim {\lambda }_{\mathrm{rest}}^{1.06}$, which is less steep than ${R}_{\lambda }\sim {\lambda }_{\mathrm{rest}}^{4/3}$ based on k_B T(R_λ ) = hc/λ. 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}_{\lambda }^{3/2}\sim {\lambda }^{2}$. In the CHAR model, the damping timescale scales as ${R}_{\mathrm{var}}^{\alpha }\sim {\lambda }^{s\alpha }$, which is less steep than the scaling relation of λ².

Zoom In Zoom Out Reset image size

3.3.3. The Relationship between R_var/R_S and R_λ /R_S

For a given wavelength λ, the emission-region size at which k_B T(R_λ ) = hc/λ for the CHAR model is (Sun et al. 2020a)

$$\begin{eqnarray}&&{R}_{\lambda }/{R}_{{\rm{S}}}={\left(\displaystyle \frac{3(1+k){k}_{{\rm{B}}}^{4}{L}_{\mathrm{Edd}}}{64\pi \sigma \eta {G}^{2}{M}_{\odot }^{2}{h}^{4}}{\lambda }^{4}{\left(\displaystyle \frac{{M}_{\mathrm{BH}}}{{M}_{\odot }}\right)}^{-2}\dot{m}\right)}^{1/3},\end{eqnarray} \tag{ 8 }$$

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:

$$\begin{eqnarray}&&{\mathrm{log}}_{10}({R}_{\mathrm{var}}/{R}_{{\rm{S}}})=A{\mathrm{log}}_{10}({R}_{\lambda }/{R}_{{\rm{S}}})+B.\end{eqnarray} \tag{ 9 }$$

Figure 12 shows the relationship between R_var/R_S and R_λ /R_S for M_BH = 10^8.0 M_⊙ and $\dot{m}=0.05$. The relationship between R_var/R_S and R_λ /R_S is nonlinear since A is less than unity.

Zoom In Zoom Out Reset image size

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 $\dot{m}$. We cannot find a simple function to describe the relationship between A (or B) and M_BH or $\dot{m}$.

Zoom In Zoom Out Reset image size

3.3.4. The Dependencies of τ_DRW,disk Upon R_λ /R_S

We also aim to establish a relationship between the directly measurable quantity τ_DRW,disk and R_λ /R_S. The fitted equation is

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{\tau }_{\mathrm{DRW},\mathrm{disk}}={K}_{1}{\mathrm{log}}_{10}({R}_{\lambda }/{R}_{{\rm{S}}})+{K}_{2}.\end{eqnarray} \tag{ 10 }$$

Figure 14 is the relationship between τ_DRW,disk and R_λ /R_S for M_BH = 10^8.0 M_⊙ and $\dot{m}=0.05$.

Zoom In Zoom Out Reset image size

Again, K₁ and K₂ depend upon M_BH or $\dot{m}$ (Figure 15), and the best-fitting relationships are

$$\begin{eqnarray}&&{K}_{1}={0.47}_{-0.02}^{+0.02}{\mathrm{log}}_{10}\displaystyle \frac{{M}_{\mathrm{BH}}}{{M}_{\odot }}+{0.48}_{-0.02}^{+0.02}{\mathrm{log}}_{10}\dot{m}-{2.22}_{-0.16}^{+0.16},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{K}_{2}=-{0.92}_{-0.06}^{+0.06}{\mathrm{log}}_{10}\displaystyle \frac{{M}_{\mathrm{BH}}}{{M}_{\odot }}-{0.36}_{-0.06}^{+0.06}{\mathrm{log}}_{10}\dot{m}+{1.57}_{-0.43}^{+0.43}.\end{eqnarray} \tag{ 12 }$$

We stress that K₁ is a function of M_BH and $\dot{m}$, rather than being fixed to 3/2 as expected from the thermal timescale.

Zoom In Zoom Out Reset image size