Characterizing the $\gamma$-Ray Variability of Active Galactic Nuclei with Stochastic Process Method

The $\gamma$-ray astronomy in time domain has been by now progressed further as the variabilities of Active Galactic Nuclei (AGNs) on different timescales have been reported a lot. We study the $\gamma$-ray variabilities of 23 jetted AGNs through applying a stochastic process method to the ~12.7 yr long-term light curve (LC) obtained by Fermi-Large Area Telescope (Fermi-LAT). In this method, the stochastically driven damped simple harmonic oscillator (SHO) and the damped random walk (DRW) models are used to model the long-term LCs. Our results show that the long-term variabilities of 23 AGNs can be characterized well by both SHO and DRW models. However, the SHO model is restricted in the over-damped mode and the parameters are poorly constrained. The SHO power spectral densities (PSDs) are same as the typical DRW PSD. In the plot of the rest-frame timescale that corresponds to the broken frequency in the PSD versus black hole mass, the intrinsic $\gamma$-ray characteristic timescales of 23 AGNs occupy almost the same space with the optical variability timescales obtained from the accretion disk emission. This suggests a connection between the jet and the accretion disk. Same as the optical variability of AGN accretion disk, the $\gamma$-ray timescale is also consistent with the thermal timescale caused by the thermal instability in the standard accretion disk of AGN.


INTRODUCTION High energy (HE;
100 MeV) γ-ray observations suggest that the emissions from Active Galactic Nuclei (AGNs) dominate the extragalactic γ-ray sky (Abdollahi et al. 2020).The strongly Doppler-boosted blazars, an extreme class of AGNs of which emissions are mainly from the nonthermal relativistic jets, are dominant in these powerful emitters.Blazars are classified into BL Lac objects (BL Lacs) and flat spectrum radio quasars (FSRQs), according to the strength of their optical emission lines (Ajello et al. 2020).FSRQs have strong, broad emission lines, while BL Lacs have weak, narrow, or no such lines.
AGN variability has already been detected at entire electromagnetic wavelengths with timescales covering from decades down to minutes.Radio-loud AGNs are highly variable γ-ray emitters.This not only applies to the blazars which have strong and incessant flux variability, but also to misaligned jet sources such as radio galaxies (MAGIC Collaboration et al. 2018;Ait Benkhali et al. 2019).The underlying physical processes can be investigated by characterizing the variabilities (e.g., H. E. S. S. Collaboration et al. 2017;Yan et al. 2018;Rieger 2019;Bhatta & Dhital 2020).
The observations of Fermi-Large Area Telescope (Fermi-LAT) have advanced the studies in HE time-domain.An attractive phenomenon is the γ-ray quasi-periodic oscillations (QPOs) detected in LAT data of blazars (e.g., Ackermann et al. 2015;Sandrinelli et al. 2016;Zhou et al. 2018;Peñil et al. 2020;Zhang et al. 2020).However, the reliability of these QPOs is always questionable (e.g., Covino et al. 2019;Ait Benkhali et al. 2020;Yang et al. 2021).Another interesting phenomena is the fast γ-ray flares on timescale of a few minutes detected in the LAT data of FSRQs (Ackermann et al. 2016;Meyer et al. 2019;Shukla & Mannheim 2020).Besides, statistical characteristics of the γ-ray variability in AGNs have been extensively investigated.The commonly used methods are the analyses of power spectral density (PSD) and flux distribution (e.g., Abdo et al. 2010;Shah et al. 2018;Meyer et al. 2019).Abdo et al. (2010) presented γ-ray variabilities of the 106 AGNs systematically, using the first 11 months of the Fermi survey.They reported that more than 50 % of the sources are found to be variable with a power-law (PL) PSD, and a random walk underlying mechanism was reflected in some blazars.Nakagawa & Mori (2013) reported a characteristic timescale of ∼7.9 days in the PSD of 3C 454.3 by analyzing the first four years Fermi-LAT data, and they used an internal shock model to interpret this timescale.Meyer et al. (2019) presented a detailed analysis of LAT LCs of bright γ-ray FSRQs, and put strong constraints on blazar jet physics accordingly.
A stochastic process model has been wildly used to describe optical variability of AGN accretion disk (e.g., Kelly et al. 2009;MacLeod et al. 2010;Zu et al. 2013;Rakshit & Stalin 2017;Li & Wang 2018;Zhang et al. 2018;Burke et al. 2021).Generally, the damped random walk (DRW) model can provide successful fit to the long-term variability of AGN accretion disk.It is proved that such a stochastic process model provides a powerful tool to extract information from AGN variability (e.g., Kasliwal et al. 2017;Burke et al. 2021).In the past few years, the stochastic process model has been applied to γ-ray variabilities of AGNs in several papers (Sobolewska et al. 2014;Goyal et al. 2018;Ryan et al. 2019;Tarnopolski et al. 2020;Covino et al. 2020;Yang et al. 2021;Zhang et al. 2021).Based on the stochastic model developed by Kelly et al. (2009), Sobolewska et al. (2014) modeled the γ-ray LCs of 13 blazars observed during the first four years of the Fermi sky survey with the OrnsteinUhlenbeck (OU) process (also called DRW) and a linear superposition of the OU processes (sup-OU).They showed that 10 of 13 blazars prefer the sup-OU process over the OU process.The continuous-time autoregressive moving average (CARMA) method (Kelly et al. 2014), a generalized stochastic model which can be applied to astronomical time series, is flexible to capture the features of flux variability and to produce more accurate PSD.Applying this method to the 9.5 yr LAT data of the same 13 blazars in Sobolewska et al. (2014), Ryan et al. (2019) reported that the DRW model is good enough to describe the γ-ray variability of the 13 blazars.In addition to CARMA, celerite is a newly developed method for modeling LC with the stochastic process model (Foreman-Mackey et al. 2017).It was applied to the LAT data of AGNs to examine the significance of the γ-ray QPOs (Yang et al. 2021;Zhang et al. 2021).
So far, the Fermi-LAT with collecting data for more than 12 yr has provided an excellent opportunity to study the long-term γ-ray variability in AGNs.In this paper, we apply the celerite model to 12.7 yr Fermi-LAT LCs of 23 bright LAT AGNs including 10 BL Lac objects, 12 FSRQ objects and one radio galaxy.We aim to investigate the γ-ray variabilities of AGNs on long-term timescales.The format of this paper is as follows.In Section 2, we briefly introduce the Fermi-LAT data processing method.In Section 3, the stochastic process models are briefly described.In Section 4, we show the modeling results for the LCs of 23 AGNs with celerite method.In Section 5, we focus on the variability characteristic timescales in the jets, and compare them with optical results obtained from AGN accretion disk emissions.In Section 6, we discuss our results.Finally, a summary is presented in the Section 7.

FERMI-LAT DATA ANALYSIS
We collect 23 bright γ-ray AGNs with significant variability in the LAT data.The information of these sources are listed in the Table 1.
All data analyzed here come from the LAT 8 yr Source Catalog (4FGL; Abdollahi et al. 2020), spanning the time range of MJD 54682-59332 which gives in total ∼ 12.7 yr of data in the energy range of 0.1-500 GeV.We consider only SOURCE class events (evclass=128) and event type three (evtype=3) from the region of interest (ROI) at 15 • for each source.The maximum zenith angle is set to be 90 • to avoid contaminating from Earth's limb.DATAQUAL >0 and LATCONFIG == 1 options are chosen to ensure the good data quality and the proper time intervals.The instrument response function P8R3 SOURCE V3 is applied in the analysis.We use the Galactic (gll iem v07.fits) and the extragalactic (iso P8R3 SOURCE V3 v1.txt) diffuse emission models, which are the latest Pass 8 background models.We use the binned maximum likelihood analysis1 (Abdo et al. 2009), which is the preferred method for most types of LAT analysis.It is a three-dimensional maximum likelihood algorithm, i.e., events are binned into channels of energy and position in sky (Abdo et al. 2009), and a maximum likelihood optimization technique is Note-(1)(4) source name, (2)(3) RA and Dec (J2000), (5) redshift, (6) source type, (7) Eddington ratio from Xiong et al. (2015), (8) black hole masses and their uncertainties (in solar mass) collected from the references in the last column.The uncertainty in the relation between stellar velocity dispersion and black hole mass (MBH − σ) is 0.21 dex (Tremaine et al. 2002), and we use 0.2 dex.The uncertainty on the zero point of the line width-luminosity-mass relation is approximately 0.5 dex (Gebhardt et al. 2000;Ferrarese et al. 2001) performed to determine the best-fit parameters and the Test Statistic TS = 2∆log(likelihood) between models with and without the source (Mattox et al. 1996).In the fitting, we use the spectral model of the LogParabola (LP) form

STOCHASTIC PROCESS METHOD
We use the stochastic process method implemented in celerite package (Foreman-Mackey et al. 2017).In celerite package, a specific and stationary kernel function (i.e., covariance function) is required, which can be defined by users.

DRW Model
The DRW process is described by a first order stochastic differential equation (see details in Kelly et al. 2009;Moreno et al. 2019).It represents a competition between a process trying to maintain an equilibrium state and a perturbation  (3) model, (4) AICc, (5)(6)(7) the results and the reduced χ 2 of normal distribution fitting to the standard residuals.making the system out of stability.It is sometimes also written as a Langevin equation of the form where τ DRW is the damping timescale and σ DRW is the amplitude of random perturbations.Following the setting in celerite package, the covariance function for the DRW is written as where t nm = |t n − t m | is the time lag between measurements m and n, a = 2σ2 DRW , and c = 1/τ DRW .The PSD is written as (Foreman-Mackey et al. 2017) (3) The DRW PSD is a broken power-law form, and the index changes from 0 at low frequencies to -2 at high frequencies.

SHO Model
The dynamics of a stochastically driven damped simple harmonic oscillator (SHO) provides a physically motivated model, as it can describe the variability driven by noisy physical processes, which grows most strongly at the characteristic timescale but is also damped owing to dissipation in the system (Foreman-Mackey et al. 2017).The differential equation for this system is with the frequency of the undamped oscillator ω 0 , the quality factor of the oscillator Q, and a stochastic driving force (t).When the (t) is white noise, the PSD of this process is where S 0 is proportional to the power at ω = ω 0 .The SHO PSD is complex.For the over-damped SHO (Q < 0.5), it is also a broken power-law form at low frequencies, very similar to the DRW PSD, while at high frequencies, the index can be as small as ∼ −4 (Kasliwal et al. 2017;Moreno et al. 2019).For the under-damped SHO (Q > 0.5), a Lorentzian appears in the PSD, i.e., a QPO signal (Foreman-Mackey et al. 2017;Moreno et al. 2019).

Model Selection
Akaike information criterion (AIC) estimates the relevant information that is lost when models are used to represent the underlying processes that generate the data.It is an estimator of the relative quality of models for a given set of data: the less information a model loses, the higher the quality of that model.We use the corrected AIC (AIC c ) to perform model selection, which is given by where k is the number of model parameters, L is the maximum likelihood, and n is the number of data points.A preferred model is one that minimizes AIC c .It is accepted that ∆(AIC c ) 10 is a difference substantial enough to prefer the model with smaller AIC c (Burnham & Anderson 2004;Sobolewska et al. 2014).

RESULTS
In the fittings, Markov chain Monte Carlo (MCMC) implemented in the package emcee 2 (Foreman-Mackey et al. 2013) is used to sample the posterior probability density in our analysis.The priors for the parameters are assumed to be flat.We run the MCMC sampler for 50,000 steps of which the first 20,000 steps are taken as burn-in.We calculate the maximum likelihood for optimization which is executed 100 times to resolve possible instability of the algorithm L-BFGS-B, and then calculate AIC c by using the maximum value among the 100 values.
For the data from MJD 54682 to 59332, the 15-day binning LCs of the 23 AGNs are structured by performing the binned likelihood method for each time bin.The time bins having TS value of ≥ 25 are selected here, in order to get reliable and high signal-to-noise ratio results (e.g., Kapanadze et al. 2020).In Table 2, we report the mean cadence for each LC.
Each LC is fitted with the SHO and DRW models, respectively.The goodness of the fit is assessed by analyzing the probability densities of the standardized residuals and the auto-correlation function (ACF) of the standardized residuals.The distribution of the standardized residuals is fitted by a normal distribution.The parameters and the reduced χ 2 (χ 2 red ) are given in Table 2.It is shown that the distribution of the standardized residuals is in good agreement with the normal distribution with the mean value close to zero and the standard deviation less than one.In Figure 1, we show the fitting results for 3C 454.3 and 3C 279 for an example.The ACFs of the residuals are randomly distributed around zero, and are almost inside the 95% confidence limits of the white noise.It indicates that the model has captured the characteristics of the time series.It is notable that the standardized residuals corresponding to the two or three highest flux are large.The DRW and SHO models likely fail to describe the brightest flares.When we fit the LC excluded the two or three highest flux points, the modeling results are unchanged.
We give the posterior probability density distribution of parameters resulting from the SHO and DRW modelings in Figure 2. The values of the model parameters are given in Table 3.One can see that the model parameters in DRW model are constrained well, although the two parameters are degenerate.While in the SHO model, the strong degeneracy between ω 0 and Q leads to large uncertainties on the two parameters.In some cases, the large uncertainties on Q cause the unilateral distribution of ω 0 , and the upper limit for ω 0 is meaningless.For all sources, the SHO model is constrained in the over-damped mode (Q < 0.5).
In Figure 3, we show the PSDs for 3C 454.3 and 3C 279 constructed from the modeling results with SHO and DRW models.The PSDs for SHO and DRW are almost the same.The index changes from 0 at frequencies below f b to -2 at frequencies above f b .
The difference between the AIC c for SHO and DRW are small (Table 2), indicating that the fittings of the two models are comparable.However, the poor constrained ω 0 and Q suggest that the DRW model is preferred over the SHO model.
There is no significant difference among the forms of the PSDs from different types of AGNs (Figure 4).The PSDs for the 23 sources are typical DRW PSD, and the broken frequencies are between 0.001 day −1 and 0.01 day −1 .

VARIABILITY CHARACTERISTIC TIMESCALES IN RELATIVISTIC JETS
Characteristic timescale is an important parameter in source variability.Our results show that DRW model can describe the γ-ray variability of 23 targets successfully, and the higher-order model SHO is unnecessary.From the fittings with the DRW model, we can obtain the characteristic timescales for the 23 sources.The characteristic (damping) timescales with errors in 1σ are given in Table 3.One can see that the timescales are between 20 days and 250 days.Note that the measurement of the damping timescale can be biased by the insufficient length of the LC (Koz lowski 2017; Suberlak et al. 2021).Moreover, Burke et al. (2021) found that the measurement of damping timescale is reliable when it is larger than the typical cadence of LC.We use the following criteria to select the reliable measurements of the damping timescale: (1) The length of the LC should be ten times of the timescale at least; (2) The derived timescale should be larger than the mean cadence of the LC.It is found that all the γ-ray characteristic timescales for the 23 sources are reliable.The LCs are fitted in the observed frame, and the values of the damping timescales in Table 3 are in the observed frame.In order to get the timescale in the rest frame (τ rest damping ), the timescale should be corrected by cosmological time dilation and and Doppler beaming effect, where δ D is the Doppler factor.The Doppler factor for the γ-ray emission region in AGN jet is difficult to measure.It is estimated by different methods based on, for example, modeling the broadband spectral energy distributions (Chen 2018;Pei et al. 2020), opacity of the γ-rays, and the brightness temperature of radio flare (Liodakis et al. 2017).In Table 4, we list the average values of δ D for different types of AGNs estimated by recent three papers (Liodakis et al. 2017;Chen 2018;Pei et al. 2020).One can see that the uncertainties on the Doppler factor are very large.We use the Figure 5. Variability damping timescale (in the rest frame) as a function of black hole mass.The gray data, lines and area represent optical results taken from Burke et al. (2021).The data in color are our results from the γ-ray LCs of AGNs.Data points in red, blue and black represent FSRQ, BLL and RDG, respectively.average results of the three papers to correct the timescale.The average δ D for blazars is 10.It is found that τ rest damping is between 100 days and 1500 days, and the average τ rest damping is ≈ 510 days.The characteristic timescale in the optical variability of AGN accretion disk has been extensively studied (e.g., Collier & Peterson 2001;Kelly et al. 2009;MacLeod et al. 2010;Simm et al. 2016;Suberlak et al. 2021).Recently, Burke et al. (2021) reported a correlation between the optical characteristic timescale and the black hole mass.We also show our results in the plot of τ rest damping − M BH (Figure 5), together with the results in Burke et al. (2021).It can be seen that the γ-ray variability timescales of AGNs occupy the same space with the optical variability timescales.Namely, in the same range of the black hole mass, the γ-ray variability timescales are consistent with the optical variability timescales within the errors.There is no correlation between the γ-ray variability timescale of AGN and the black hole mass.This is probably due to the small dynamic range in the black hole mass in the sample of ∼ 10 8 − 10 9 M (with the exception of NGC 1275) making any correlation difficult to be identified.We have searched the γ-ray AGNs with smaller black hole mass.However, they are not bright enough to perform the variability analysis.In Appendix A, we use the γ-ray flares from Crab Nebula to extend the mass of the central engine to much smaller range.
The γ-ray τ rest damping values slightly and systematically lie above the optical relation of Burke et al. (2021).In addition, we fit our γ-ray results and the optical results together, resulting in the best-fit relation, with an intrinsic scatter of 0.23 ± 0.03 dex and Pearson correlation coefficient r = 0.80.It is similar to the optical result.
We get the Eddington ratio (the ratio between accretion disk luminosity and the Eddington luminosity) for 15 out of 23 AGNs from Xiong et al. (2015), which are listed in Table 1.Except for Mrk 421, BL Lac and OJ 287, the rest of 12 sources have Eddington ratio between 0.01 and 1.From Figure 6, one can see that the 12 sources have similar Eddington ratio with the normal quasars in Burke et al. (2021), and the characteristic timescale is independent on the Eddington ratio.The stochastic process models are increasingly used to model γ-ray variability of blazars, providing an effective tool to study the statistical properties of the variability.We model the γ-ray variability of 23 bright LAT AGNs with the DRW and SHO models by performing the celerite package.It is found that the DRW model with two parameters can fit the γ-ray variability of 23 AGNs successfully, and the fits with the SHO model with three parameters for the 23 sources are not improved.The SHO is constrained in the over-damped mode (Q < 0.5), and the strong degeneracy between ω 0 and Q leads to poor constraints on the two parameters.The PSDs of the 23 sources are typical DRW form with f b between 0.001 day −1 and 0.01 day −1 .
The brightest flares ( 10 −5 ph cm −2 s −1 ) in 3C 279 and 3C 454.3 are poorly fitted by both DRW and SHO models.Indeed, the extreme γ-ray flares in 3C 279 seem special (e.g., Shukla & Mannheim 2020).The γ-ray photon index in the extreme flare on Dec 16 2012 is 1.7 (Hayashida et al. 2015), significantly smaller than the typical γ-ray photon index of 2.4 for FSRQ (Abdollahi et al. 2020).Minute-scale GeV γ-ray variability from 3C 279 was observed in an extreme flare on June 15 2016 (Ackermann et al. 2016).Nalewajko (2013) studied the individual γ-ray flares of 3C 454.3 with the flux above 0.71×10 −5 ph cm −2 s −1 , and found that the γ-ray flares of 3C 454.3 have more complex light curves than other blazars.The extreme flares may have a different physical mechanism than the underlying long-term stochastic variability.
These brightest flares are expected to have impact on the slope of PSD at high frequencies (Ryan et al. 2019).We have examined that the brightest flares in 3C 279 and 3C 454.3 cannot affect the modeling results for the long-term variabilities.A further and careful study on the brightest flares is worthy of performing by using an adaptive binning algorithm.
The theoretical PSD expected by the one-zone leptonic emission model has been investigated (Finke & Becker 2014, 2015;Chen et al. 2016).Thiersen et al. (2022) simulated multi-wavelength variability of blazars from a purely numerical approach by using a time-dependent one-zone leptonic emission model.They showed that a power-law PSD for the emission variability is produced by introducing stochastic variations for model parameters in the emission region, and the PSD is similar to the underlying power law of the model parameter variation.No spectral break is found in their produced PSDs.The results of Thiersen et al. (2022) indicate that in the frame of one-zone emission model, the physical processes associated with electron cooling, light crossing, and electron escape would not produce a break in the PSD.The broken frequencies we obtained are between 10 −8 Hz and 10 −7 Hz.The corresponding intrinsic timescale is several hundred days at least, which cannot be the timescale corresponding to electron cooling or acceleration process.
The γ-ray timescales of AGNs we obtained are very close to the optical timescales obtained from modeling AGN accretion disk emissions in Burke et al. (2021).Burke et al. (2021) speculated that the optical timescales could be associated with the thermal timescales 3 expected in the AGN standard accretion disk theory, and the optical variability may be driven by the thermal instability of the accretion disk.The similarity between the γ-ray and optical characteristic timescales could imply a connection between jet and accretion disk.The thermal instability may also causes the γ-ray variability in the jet.However, the detailed mechanism that connects the accretion disk and the jet is unclear.The γ-ray timescales are slightly larger than the optical timescales of normal quasars.This may be due to that the distance from the γ-ray emission region to accretion disk extends the intrinsic timescale from accretion disk.Ruan et al. (2012) modeled the nonthermal optical variabilities of 51 γ-ray blazars, and found that blazar optical τ rest damping peaks at ∼1000 days (assuming a typical Doppler factor of 10), which is systematically larger than the γ-ray τ rest damping in this work.The discrepancy between blazar γ-ray and optical τ rest damping may imply that the γ-ray and optical emissions are produced in different regions.The γ-ray emission region is closer to the accretion disk than the optical emission region.Ruan et al. (2012) found that blazar nonthermal optical characteristic timescales are ∼4 times smaller than normal quasars.They considered that the discrepancy between the optical characteristic timescales for blazars and normal quasars could be caused by the Doppler effect, if the jet variability and accretion disk variability have the same origin.Combining with our γ-ray results, we suppose that the discrepancy between the characteristic timescales for blazars and normal quasars is not only caused by the Doppler effect, but also related to the location of the jet emission region (the distance from the accretion disk).The jet long-term variability may be the convolution of the accretion disk variability with a transfer function which is related to Doppler factor and the distance from the jet emission region to the accretion disk at least.

SUMMARY
We have applied a stochastic process method to the ∼12.7 yr Fermi-LAT LCs of 23 jetted AGNs in order to investigate the γ-ray variability properties.The SHO and DRW models are both used to model the long-term LCs.Our main results are as follows.
(i) The long-term variability of 23 sources in our sample can be described well by both SHO and DRW models.However, the modelings with the SHO are not improved, and the parameters ω 0 and Q are poorly constrained.This suggests that the DRW model is preferred over the SHO model for the γ-ray long-term variability of AGNs.The PSDs for the 23 sources are the typical DRW PSD form.
(ii) The intrinsic characteristic timescale of AGNs extracted from modeling the γ-ray variability is between 100 days to 1500 days.Such a long timescale cannot be produced in a one-zone leptonic emission model within the typical parameter space.In the plot of τ rest damping − M BH , the γ-ray timescales obtained from jet emissions occupy almost the same space with the optical timescales obtained from the accretion disk emissions.Both the γ-ray and optical In conclusion, our results suggest that the origin of the γ-ray variability could be related to the thermal instability in the accretion disk, however the detailed process that drives the variability is unclear.

Figure 1 .
Figure1.Fitting results of 3C 454.3 and 3C 279 for example.For each source, the left column presents results obtained from DRW model and we give the LAT LC (black points) and the modeled LC (orange line) in the top panel.In the middle panel, we show the standardized residuals (black points), the probability density of standardized residuals (orange histogram) as well as the best-fit normal distribution (green solid line).The reduced χ 2 is labeled in the figure.The ACF of residuals with the 95% confidence limits of the white noise (the gray region) are shown in the bottom panel.The results obtained by the SHO model are given in the right column.The symbols and lines are the same as the left column but with different colors.

Figure 2 .
Figure 2. Posterior probability densities of model parameters for SHO (left) and DRW (right).The top panels are for 3C 454.3, and the bottom panels are for 3C 279.The vertical dotted lines mark the median value and 68% confidence intervals of the parameters distribution.

Figure 6 .
Figure6.Plot of the rest-frame timescale versus the Eddington ratio.The gray data points with error bars are optical results taken fromBurke et al. (2021).The data points in red denote the γ-ray results of blazars in our samples.

3
The thermal timescale reads t th = 4.6 × α yrs, where R is the emission distance on the accretion disk from the central black hole, R S = 2GM BH /c 2 is the Schwarzschild radius, and α is the standard disk viscosity parameter.timescalesare consistent with the thermal timescale expected by the AGN standard accretion disk.It may indicate a connection between the jet and the accretion disk.

Figure 7 .
Figure 7. DRW modeling results of γ-ray flare (MJD 55654.65-55678.65)from the Crab Nebula.The symbols and lines are the same as those in Figure 1.

Figure 8 .
Figure 8. Left panel: posterior probability densities of DRW parameters for the Crab Nebula.The symbols and lines are the same as those in Figure 2. Right panel: DRW PSD of the γ-ray LC of the Crab Nebula.The corresponding color region denotes 1σ confidence interval.

Figure 9 .
Figure 9. Variability damping timescale (in the rest frame) as a function of the mass of the central engine.The gray data, lines, area as well as the crosses represent optical results taken from Burke et al. (2021).The data in color are our results from the γ-ray LCs of AGNs and the Crab Nebula.The orange line and shaded band are the best-fit relation and 1σ uncertainty for 23 AGNs and the Crab Nebula.

Table 4 .
Doppler factor for AGNs