Gaussian Process Modeling Fermi-LAT γ-Ray Blazar Variability: A Sample of Blazars with γ-Ray Quasi-periodicities

Blazars are an extreme class of active galactic nuclei (AGN) with their relativistic jets closely aligning with our line of sight, whose central engines are the super massive black holes (SMBHs) located in the cores of the host galaxies. Blazar emission is dominated by the Doppler-boosted and nonthermal emission of the powerful jet. Blazars are classified into BL Lac objects and flat spectrum radio quasars (FSRQs) according to their broad emission lines (e.g., Stickel et al. 1991; Urry & Padovani 1995). FSRQs have strong emission lines, whereas BL Lac objects have weak/no emission lines. The blazars included in the fourth catalog of AGN detected by the Fermi Large Area Telescope (LAT) consist of 38% BL Lac objects, 24% FSRQs, and 38% blazar candidates of unknown types (Ajello et al. 2020). On average, FSRQs have softer spectra and stronger variabilities in GeV γ-ray energies, compared with BL Lac objects (Ajello et al. 2020), although LAT spectra of FSRQs and BL Lac objects display different properties in flares (e.g., Williamson et al. 2014; Hayashida et al. 2015).

Flux variabilities of blazars have been detected at entire electromagnetic wavelengths from radio band to γ-rays. The variability timescales cover many orders of magnitude, from decades down to minutes. Variability analysis is a powerful tool to probe the nature of blazars (e.g., Rieger 2019).

Thanks to the successful operation of the Fermi-LAT, blazar variability study has progressed significantly in the GeV γ-ray regime. Rapid γ-ray flare on a timescale of a few minutes has been detected by LAT for the FSRQ 3C 279 (Ackermann et al. 2016). The rapid flare is considered to indicate extreme jet conditions (Petropoulou & Dermer 2016; Petropoulou et al. 2017; Asano & Hayashida 2020). Besides, with the LAT collecting data for more than 10 yr, the long-term characteristic of the flux variability becomes attractive.

In the last few years, possible quasi-periodic oscillations (QPOs) have been found in the LAT data of more than 20 blazars (e.g., Ackermann et al. 2015; Sandrinelli et al. 2016; Prokhorov & Moraghan 2017; Zhang et al. 2017b, 2017a, 2017c, 2020; Sandrinelli et al. 2018; Zhou et al. 2018; Peñil et al. 2020). The timescales of these QPOs span from one month to three years, and BL Lac objects and FSRQs seem have similar periodicity behaviors (e.g., Zhang et al. 2020). The year-type QPOs are important for the studies on accretion disks, relativistic jets, and SMBH physics (e.g., Ackermann et al. 2015).

The results of γ-ray QPOs (e.g., Ackermann et al. 2015; Zhang et al. 2020) are obtained by performing frequently used techniques such as the Lomb–Scargle periodogram (LSP; Lomb 1976; Scargle 1982) and weighted wavelet Z-transform (WWZ; Foster 1996). The above two methods calculate Fourier-like power spectral density (PSD) of γ-ray light curve, and the QPO signal appears as a peak in the PSD. However, the current LAT data only cover a few cycles of the QPOs. Therefore, the significance for these QPOs is questionable.

Some recent investigations have challenged the results of the γ-ray periodicity searchings, and it is cautioned that those γ-ray QPOs may be fake signals (e.g., Covino et al. 2019; Ait Benkhali et al. 2020). Indeed, Vaughan et al. (2016) have pointed out that the search for periodic signals among a large number of stochastic time series should be very careful. Based on the analysis with a damped random walk (DRW) model, they show that physical periodic and stochastic periodic signals are very difficult to distinguish when the time series only cover a few cycles (≲5 cycles). Covino et al. (2019) reanalyze LAT data of 10 blazars, which were reported with possible γ-ray QPOs, and they do not find any QPO signals in their sample. Ait Benkhali et al. (2020) analyze 9.5 yr LAT data of six QPO blazars, and also do not find strong evidence for any QPOs. Therefore, the traditional techniques for periodicity searching in blazar variability might be misleading, and the reported QPOs should be diagnosed with different methods.

An alternative approach is to fit the light curves using Gaussian processes in the time domain (e.g., Kelly et al. 2009; Kozłowski et al. 2010; MacLeod et al. 2010). In this case, the PSD can be calculated from the inferred autocovariance function (Kelly et al. 2014). In this field, the representative method is the continuous-time autoregressive moving average (CARMA) models developed by Kelly et al. (2014). Actually, the DRW model is the case of CARMA(1,0). CARMA models are flexible to capture the features of flux variability and to produce more accurate PSD. The CARMA method has been already applied to LAT data of blazars (e.g., Goyal et al. 2018; Ryan et al. 2019). Comparably, Foreman-Mackey et al. (2017) developed another fast and flexible Gaussian process model, celerite, for estimating the variability features of a light curve and its PSD. The capability of celerite for characterizing the variability of AGN is unknown, and it requires more investigations to find out its potentiality in the AGN variability studies.

In this paper, we carefully study the γ-ray light curves of 27 Fermi blazars, which have possible QPO features, using the different Gaussian process methods: CARMA and celerite. We aim to primarily figure out two points: (i) whether or not the γ-ray light curves can be fitted well by the Gaussian processes and (ii) whether or not the reported QPO signals still appear in the derived PSDs. Additionally, we also test the efficiency of the celerite model for the AGN variability studies. The format of this paper is as follows. In Section 2, we briefly introduce the basic concepts of CARMA and celerite. In Section 3, we briefly describe the analysis procedure of Fermi data covering ∼11 yr. In Section 4 we give the modeling results of the light curves with CARMA and celerite. Finally, we discuss the results in Section 5.

If a light curve y(t) is considered as a zero-mean CARMA(a, b) process, it can be obtained by solve the stochastic differential equation (Kelly et al. 2014):

$$\begin{eqnarray}&&\begin{array}{l}{\alpha }_{a}\displaystyle \frac{{d}^{a}y(t)}{{{dt}}^{a}}+{\alpha }_{a-1}\displaystyle \frac{{d}^{a-1}y(t)}{{{dt}}^{a-1}}+\cdots +{\alpha }_{0}y(t)\\ \ \ ={\beta }_{b}\displaystyle \frac{{d}^{b}\epsilon (t)}{{{dt}}^{b}}+{\beta }_{b-1}\displaystyle \frac{{d}^{b-1}\epsilon (t)}{{{dt}}^{b-1}}+\cdots +{\beta }_{0}\epsilon (t),\end{array}\end{eqnarray} \tag{ 1 }$$

where (t) is a white-noise process with zero mean and variance σ², and α_a and β₀ are defined to equal 1. The PSD of the stationary CARMA(a, b) process can be written as

$$\begin{eqnarray}&&P(f)={\sigma }^{2}\displaystyle \frac{{\left|{\displaystyle \sum }_{j=0}^{b}{\beta }_{j}{\left(2\pi {if}\right)}^{j}\right|}^{2}}{{\left|{\displaystyle \sum }_{k=0}^{a}{\alpha }_{k}{\left(2\pi {if}\right)}^{k}\right|}^{2}}.\end{eqnarray} \tag{ 2 }$$

In practice, the parameters σ², α₀, ⋯ ,α_a−1, β₁, ⋯ , and β_b can be obtained by fitting y(t) to an observed light curve with determined a and b, and then it is convenient to calculate the PSD using Equation (2). Any characteristic timescale including QPO could manifest a structure on the PSD. For more details about CARMA, the reader is referred to Kelly et al. (2014). The software we used in this paper is carma_pack³ released by Kelly et al. (2014).

Celerite (Foreman-Mackey et al. 2017) is also a Gaussian process method for modeling light curves, and shares many features with CARMA. According to its original kernel function (i.e., covariance function), it can be used to calculate any CARMA model theoretically. The official software celerite⁴ is designed to be restrictive and it requires a specific kernel function chosen or provided by the user. In other words, users are free to analyze their time series with stochastic or physical models. Nevertheless, we adopt the stochastically driven damped simple harmonic oscillator (SHO) as the model for our condition, following the recommendation of the developers. An SHO system has the differential equation:

$$\begin{eqnarray}&&\left[\displaystyle \frac{{d}^{2}}{{{dt}}^{2}}+\displaystyle \frac{{\omega }_{0}}{Q}\displaystyle \frac{d}{{dt}}+{\omega }_{0}^{2}\right]y(t)=\epsilon (t),\end{eqnarray} \tag{ 3 }$$

where ω₀ is the frequency of the undamped oscillator, Q is the quality factor of the oscillator, and (t) is assumed to be white noise. Then the PSD can be expressed as

$$\begin{eqnarray}&&P(\omega )=\sqrt{\displaystyle \frac{2}{\pi }}\displaystyle \frac{{P}_{0}{\omega }_{0}^{4}}{{\left({\omega }^{2}-{\omega }_{0}^{2}\right)}^{2}+{\omega }_{0}^{2}{\omega }^{2}/{Q}^{2}},\end{eqnarray} \tag{ 4 }$$

where P₀ is proportional to the power at ω = ω₀. As CARMA does, a specific model for a time series can be built up with a mixture of determined number of SHO terms. The number of the oscillators should be determined using model selection criteria.

We collect 27 blazars with possible γ-ray QPOs, including 17 BL Lac objects and 10 FSRQs. The reported period and relevant information of each blazar are demonstrated in Table 1.

We use the Fermi-LAT photon data of these sources covering ∼11 yr (from 2008 August 4 to 2019 September 5). We keep only SOURCE class events (evclass=128) and the three event types (evtype=3). The region of interest (ROI) for each source is a circle centered at the corresponding position with 15° radius. To avoid contamination from the Earth's limb, the maximum zenith angle is set to be 90°. The good time intervals are selected with the recommended filter expression (DATA_QUAL>0)&&(LAT_CONFIG==1). Additionally, to avoid the impact from the Sun emission, we exclude the time intervals when the distances from the Sun directions to the sources are less than 15°, with the expression angsep(src_ra, src_dec, RA_SUN, DEC_SUN)>15. We use the script make4FGLxml and the source catalog gll_psc_v19 to generate our model file for the likelihood analysis. The Galactic and extragalactic diffuse and isotropic emission models are gll_iem_v07 and iso_P8R3_SOURCE_V2_v1, respectively. The instrument response function P8R3_SOURCE_V2 and the unbinned likelihood analysis method are used. The whole procedure follows the official unbinned likelihood tutorial,⁵ and is performed in the Fermitools-conda (version 1.0.10) environment.

To produce the light curve, we first select the data between 100 MeV and 300 GeV in the time span given above, and perform the binned likelihood analysis to get the best-fit parameters. Then we set these best-fit parameters as initial parameters in the input model file, and use the unbinned likelihood method to generate the light curves. All the parameters of sources at the positions greater than 10° from the center of ROI are fixed. The normalizations of the sources in the region of ≤10° are set to be free. For each source, two light curves are built: one 30 day binning light curve and one 14 day binning light curve.

First, we would like to stress that the use of the Gaussian methods is limited for faint sources, since one also needs to include the nonrobust detections in the analysis to have a sufficient number of data points required for these models. Consequently, the results given for such sources should be treated with caution and tested by means of the alternative models for which a relatively low number of detections (only the robust ones) is required.

Therefore, we exclude the data points with Test Statistic (TS) value of  < 25 or the predicted photon number of N_pred < 10 in the light curves in order to get reliable results (see, e.g., Raiteri et al. 2013; Kapanadze et al. 2018, 2020). Unfortunately, there are three sources, GB6 J0043+3426, MG2 J130304+2434, and 87 GB 164812.2+524023, whose light curves cannot be selected by the criterion given above, because they are too faint to meet the criterion. We then apply a loose exclusion criterion of TS  < 4 to the light curves of the three faint sources.

Note that, as shown in Kelly et al. (2014) and Foreman-Mackey et al. (2017), CARMA and celerite models have the capability of handling with the irregular sampling, and they both consider the errors of the data points in the modeling.

We perform CARMA model selections to the light curves with an extended parameter space of a = 1, ⋯ ,7, b = 0, ⋯ , a − 1. In modeling a light curve, a Markov Chain Monte Carlo (MCMC) technique is used to sample the parameters, and a Bayesian criterion deviance information criterion (DIC) is calculated (Kelly et al. 2014). The model with the minimum DIC is considered to be the best one.

For celerite, a more complicated model is needed for fitting the real data. Therefore, DRW is added to the SHO model, as the DRW is used to characterize the aperiodic component in variability (e.g., Liu et al. 2018; Covino et al. 2020). Namely, the model can be written as DRW+SHO × n. In this work, n is set to be  ≤ 4. We perform the corrected Akaike information criterion (AIC) to do the model selection. The AIC is defined as

$$\begin{eqnarray}&&\mathrm{AICc}=2k-2\mathrm{log}{ \mathcal L }+\displaystyle \frac{2k(k+1)}{n-k-1},\end{eqnarray} \tag{ 5 }$$

where k is the number of parameters, n is the number of data points of light curves, and ${ \mathcal L }$ is the maximum likelihood. We improve the fitting procedure in celerite by adopting the MCMC sampler emcee⁶ (Foreman-Mackey et al. 2013). We run 100 optimizations with random starting values of parameters to avoid the effect caused by the instability of algorithm L-BFGS-B. We calculate the maximum likelihood for each optimization, and the maximum value among the output is used to calculate the AIC. We perform this procedure on each model, and the best model is the one with the minimum AIC.

In our analysis, we run the MCMC sampler for 50,000 iterations in the CARMA model. The first 25,000 iterations are taken as burn-in sampling, which is not involved in the posterior analysis. In the celerite model, we run the emcee sampler using 32 parallel walkers for 10,000 steps as burn-in and 20,000 steps for MCMC sampling.

4.1. Modeling Light Curves

The γ-ray light curves in the 0.1–300 GeV energy range of the 27 γ-ray blazars are, respectively, modeled by CARMA and celerite. The LAT light curves and modeled light curves are shown in Figure 1 (for the 24 bright sources) and Figure 2 (for the three faint sources). The modeling results are given in Table 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (24 panels)

Tarball of figure set images

Zoom In Zoom Out Reset image size

Table 2.  Modeling Results for 27 Blazars

Source	Method	30 Days			14 Days
		Model	p-value	Slope	Model	p-value	Slope
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
TXS 0059+581	CARMA	(7,6)	0.459	$-{1.02}_{-0.33}^{+0.31}$	(6,5)	0.080	$-{1.03}_{-0.33}^{+0.30}$
	celerite	S × 3	0.044	$-{1.36}_{-0.27}^{+0.25}$	S × 3	0.137	$-{1.27}_{-0.30}^{+0.21}$
PKS 0208-512	CARMA	(3,0)	0.020	$-{1.45}_{-0.22}^{+0.25}$	(7,0)	0.524	$-{1.32}_{-0.19}^{+0.17}$
	celerite	S × 2	0.003	$-{1.33}_{-0.20}^{+0.18}$	S × 4	0.209	$-{1.12}_{-0.22}^{+0.19}$
MG1 J021114+1051	CARMA	(7,6)	0.023	$-{0.85}_{-0.37}^{+0.35}$	(3,2)	0.277	$-{0.76}_{-0.41}^{+0.38}$
	celerite	D + S	0.090	$-{0.84}_{-0.21}^{+0.17}$	S × 2	0.085	$-{0.63}_{-0.21}^{+0.16}$
PKS 0250-225	CARMA	(1,0)	0.011	$-{0.72}_{-0.17}^{+0.13}$	(1,0)	0.072	$-{0.83}_{-0.16}^{+0.13}$
	celerite	D + S	0.024	−0.72 ± 0.16	S	0.002	$-{0.71}_{-0.14}^{+0.12}$
PKS 0301-243	CARMA	(4,2)	0.308	$-{0.65}_{-0.30}^{+0.48}$	(3,1)	0.379	$-{0.76}_{-0.45}^{+0.37}$
	celerite	D + S	0.000	$-{0.49}_{-0.32}^{+0.15}$	D + S × 2	0.000	$-{0.58}_{-0.16}^{+0.12}$
PKS 0426-380	CARMA	(7,5)	0.013	$-{1.16}_{-0.22}^{+0.21}$	(6,4)	0.181	$-{1.20}_{-0.19}^{+0.18}$
	celerite	D + S	0.008	$-{1.29}_{-0.16}^{+0.14}$	S × 2	0.001	$-{1.12}_{-0.17}^{+0.14}$
PKS 0447-439	CARMA	(2,1)	0.461	$-{1.04}_{-0.29}^{+0.23}$	(3,1)	0.808	$-{1.02}_{-0.26}^{+0.22}$
	celerite	S	0.260	$-{1.00}_{-0.25}^{+0.21}$	D + S	0.431	$-{0.95}_{-0.22}^{+0.17}$
PKS 0454-234	CARMA	(3,1)	0.090	$-{1.27}_{-0.29}^{+0.24}$	(4,2)	0.077	$-{1.21}_{-0.20}^{+0.15}$
	celerite	D + S	0.139	$-{1.35}_{-0.17}^{+0.15}$	D + S × 2	0.006	$-{1.04}_{-0.13}^{+0.10}$
S3 0458-02	CARMA	(1,0)	0.001	$-{0.69}_{-0.18}^{+0.13}$	(1,0)	0.001	$-{0.55}_{-0.10}^{+0.09}$
	celerite	D + S	0.001	$-{0.71}_{-0.17}^{+0.16}$	S	0.000	$-{0.50}_{-0.11}^{+0.09}$
TXS 0518+211	CARMA	(3,1)	0.455	$-{1.31}_{-0.28}^{+0.25}$	(2,1)	0.599	$-{1.22}_{-0.32}^{+0.26}$
	celerite	S	0.671	$-{1.23}_{-0.28}^{+0.22}$	S	0.422	$-{1.16}_{-0.27}^{+0.21}$
PKS 0537-441	CARMA	(5,1)	0.100	$-{1.67}_{-0.24}^{+0.20}$	(5,2)	0.224	$-{1.72}_{-0.27}^{+0.23}$
	celerite	D + S	0.173	$-{1.67}_{-0.12}^{+0.14}$	D + S × 2	0.097	$-{1.70}_{-0.15}^{+0.14}$
S5 0716+714	CARMA	(3,0)	0.133	$-{0.90}_{-0.22}^{+0.21}$	(5,0)	0.002	$-{0.85}_{-0.23}^{+0.21}$
	celerite	S	0.183	$-{0.94}_{-0.13}^{+0.12}$	S	0.000	−0.90 ± 0.21
PKS 0805-077	CARMA	(5,4)	0.119	$-{0.91}_{-0.25}^{+0.23}$	(4,1)	0.054	$-{1.19}_{-0.18}^{+0.16}$
	celerite	D + S	0.034	$-{0.99}_{-0.22}^{+0.18}$	D + S × 2	0.005	$-{0.55}_{-0.14}^{+0.12}$
OJ 014	CARMA	(2, 1)	0.972	$-{1.23}_{-0.41}^{+0.45}$	(2, 1)	0.539	$-{0.47}_{-0.55}^{+0.47}$
	celerite	S	0.794	$-{1.35}_{-0.50}^{+0.48}$	S	0.491	$-{0.64}_{-0.58}^{+0.44}$
S4 0814+42	CARMA	(2, 1)	0.686	$-{1.48}_{-0.26}^{+0.33}$	(2, 0)	0.702	$-{1.50}_{-0.24}^{+0.39}$
	celerite	S	0.782	$-{1.48}_{-0.42}^{+0.38}$	D + S × 2	0.754	$-{1.64}_{-0.32}^{+0.43}$
4C +01.28	CARMA	(7, 6)	0.819	$-{0.81}_{-0.28}^{+0.22}$	(7, 6)	0.515	$-{0.82}_{-0.28}^{+0.24}$
	celerite	D + S	0.230	$-{0.80}_{-0.21}^{+0.17}$	D + S	0.049	$-{0.68}_{-0.20}^{+0.15}$
S4 1144+40	CARMA	(6, 5)	0.179	$-{1.34}_{-0.21}^{+0.19}$	(5, 2)	0.564	−1.17 ± 0.22
	celerite	D + S	0.437	$-{1.15}_{-0.26}^{+0.18}$	D + S	0.030	$-{0.90}_{-0.23}^{+0.16}$
PG 1246+586	CARMA	(4, 1)	0.569	$-{0.68}_{-0.70}^{+0.69}$	(7, 6)	0.325	$-{0.45}_{-0.68}^{+0.56}$
	celerite	S	0.952	$-{0.79}_{-0.43}^{+0.29}$	S × 2	0.826	$-{0.88}_{-0.57}^{+0.45}$
PG 1553+113	CARMA	(4, 1)	0.776	$-{1.05}_{-0.47}^{+0.45}$	(3, 0)	0.401	$-{1.12}_{-0.39}^{+0.32}$
	celerite	S × 4	0.436	$-{0.70}_{-0.57}^{+0.41}$	S × 4	0.565	$-{0.72}_{-0.45}^{+0.49}$
TXS 1902+556	CARMA	(2, 1)	0.814	$-{0.78}_{-0.56}^{+0.44}$	(6, 4)	0.086	$-{0.41}_{-0.63}^{+0.35}$
	celerite	S × 2	0.952	$-{0.81}_{-0.66}^{+0.43}$	S × 2	0.300	$-{0.84}_{-0.75}^{+0.47}$
PKS 2052-477	CARMA	(5, 0)	0.021	$-{1.15}_{-0.26}^{-0.32}$	(6, 4)	0.075	$-{1.02}_{-0.24}^{+0.21}$
	celerite	D + S × 2	0.090	$-{1.09}_{-0.38}^{+0.25}$	S × 4	0.249	$-{0.91}_{-0.23}^{+0.14}$
PKS 2155-304	CARMA	(1, 0)	0.050	$-{1.00}_{-0.24}^{+0.18}$	(5, 0)	0.119	$-{0.95}_{-0.25}^{+0.21}$
	celerite	S	0.026	$-{0.81}_{-0.17}^{+0.13}$	S × 2	0.265	$-{0.79}_{-0.15}^{+0.13}$
BL Lacertae	CARMA	(2, 1)	0.332	$-{1.30}_{-0.20}^{+0.25}$	(2, 1)	0.039	$-{1.23}_{-0.14}^{+0.13}$
	celerite	D + S	0.341	$-{1.21}_{-0.23}^{+0.16}$	D + S	0.005	$-{1.20}_{-0.21}^{+0.18}$
PKS 2255-282	CARMA	(2, 1)	0.165	−1.43 ± 0.22	(7, 4)	0.558	$-{1.21}_{-0.29}^{+0.28}$
	celerite	S × 2	0.302	$-{1.33}_{-0.28}^{+0.23}$	D + S	0.165	$-{1.01}_{-0.23}^{+0.19}$
Three faint sources
GB6 J0043+3426	CARMA	(7, 5)	0.446	$-{0.98}_{-0.56}^{+0.48}$	(6, 3)	0.028	$-{0.77}_{-0.52}^{+0.50}$
	celerite	S × 2	0.805	$-{1.11}_{-0.58}^{+0.46}$	S × 2	0.404	$-{1.10}_{-0.65}^{+0.64}$
MG2 J130304+2434	CARMA	(1, 0)	0.466	$-{0.77}_{-0.27}^{-0.20}$	(3, 2)	0.059	$-{0.49}_{-0.28}^{+0.23}$
	celerite	D + S	0.078	$-{0.61}_{-0.22}^{+0.17}$	S × 3	0.142	$-{0.50}_{-0.21}^{+0.15}$
87 GB 164812.2+524023	CARMA	(2, 1)	0.352	$-{0.68}_{-0.62}^{+0.67}$	(7, 3)	0.048	$-{0.30}_{-0.86}^{+0.40}$
	celerite	S	0.921	$-{0.89}_{-0.72}^{+0.69}$	D + S	0.373	$-{1.51}_{-0.45}^{+0.76}$

Note. The results of the three faint sources are listed at the end of the table. (1) Source name, (2) CARMA/celerite method, (3, 6) the best CARMA/celerite model we selected for the modeling, (4, 7) KS test p-value on the normality of the standardized residuals for the fit, (5, 8) slope of the resulting PSD, and the uncertainties represent 1σ confidence intervals.

Download table as:  ASCIITypeset image

To assess the goodness of the fitting with Gaussian models, we perform a standardized residual analysis. In the Kalman filter approach that is used in CARMA, the standardized residuals, χ, can be calculated as

$$\begin{eqnarray}&&{\chi }_{i}=\displaystyle \frac{{y}_{i}-{\hat{E}}_{i}}{\sqrt{{Var}({\hat{E}}_{i})}},\end{eqnarray} \tag{ 6 }$$

where y_i is the observation, ${\hat{E}}_{i}$ is the expectation from the Kalman filtering, ${Var}({\hat{E}}_{i})$ is the variance of ${\hat{E}}_{i}$.

To link the Gaussian process of celerite to moving averages, the observed data can be expressed as

$$\begin{eqnarray}&&{\boldsymbol{y}}* \,=\,\mathrm{chol}{\left({\boldsymbol{K}}\right)}^{\top }{\boldsymbol{w}},\end{eqnarray} \tag{ 7 }$$

where y* is zero-mean vector of observation data, $\mathrm{chol}\left({\boldsymbol{K}}\right)$ is the upper triangular Cholesky factorization of the covariance matrix K, and w is the white noise which can be considered as the standardized residuals described above. Thus w can be calculated as

$$\begin{eqnarray}&&{\boldsymbol{w}}=\mathrm{chol}{\left({\boldsymbol{K}}\right)}^{-\top }{\boldsymbol{y}}* .\end{eqnarray} \tag{ 8 }$$

The expected ${\hat{E}}_{i}$ and K are appropriately obtained using the maximum-likelihood estimates.

If the model is correct, the χ or w, which are referred to as standardized residuals here, should be approximately a normal distribution with mean zero and standard deviation of one. The standardized residuals should follow a Gaussian white-noise sequence, which can be assessed through the auto-correlation function (ACF) of the standardized residuals.

With the modeling results, we calculate the standardized residuals and analyze their probability densities. We first examine the deviation of the distribution of the standardized residuals from a normal distribution using the Kolmogorov–Smirnov (KS) test. The p-values of the tests are given in Table 2. We then calculate the ACF of the residual and squared residual sequences to assess their behaviors (see Figure 1).

Looking at the distribution of the standardized residuals, one can find that the residuals of most of the fits are consistent with the expected normal distribution with p > 0.05. The 30 day light curves have smaller errors and their rapid flares are easier to smooth, compared with the 14 day light curves. In general, both CARMA and celerite can well describe the light curves of the two time bins. The ACFs of the residuals and squared residuals are almost inside the 95% confidence limits, which indicates that the models have captured the correlation structures and no significant nonlinear behaviors appear in the time series.

However, we also notice that the distribution of the standardized residuals of several sources (e.g., PKS 0250−225, PKS 0426−380, PKS 0454−234, S3 0458−02, S5 0716+714, and PKS 0805−077) deviate from the expected normal distribution. The distributions of the standardized residuals of these sources are concentrated around zero, suggestive of overfitting. The overfitting is likely due to the errors of the data.

Generally speaking, a simple DRW (or Ornstein–Uhlenbeck process) model fails to model the light curves, and high-order CARMA and celerite models are needed to characterize the γ-ray variabilities (Table 2). This is consistent with the results of previous works (Sobolewska et al. 2014; Goyal et al. 2018; Ryan et al. 2019) obtained from different samples.

It is also worth noting that both CARMA and celerite models cannot characterize large amplitude flares (see the results of PKS 0301-243), despite this failure being statistically meaningless.

No significant difference is found between the modeling results of CARMA and celerite.

4.2. Building PSDs

The PSD can be constructed by calculating Equations (2) and (4) with the modeling results. The resulting PSDs for the 24 bright sources are shown in Figure 3, and the PSDs for the three faint sources are shown in Figure 4. We also give the noise level caused by measurement errors, which is estimated by (assuming Gaussian errors)

$$\begin{eqnarray}&&{P}_{\mathrm{noise}}=2\times {\rm{\Delta }}t\times \overline{{\sigma }_{\mathrm{err}}^{2}},\end{eqnarray} \tag{ 9 }$$

where Δt is the median of sampling time and σ_err is the measurement error. It should be noted that such a noise has been deducted from each PSD.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We estimate the slope of each PSD sample by fitting it with a linear function in log–log space. The structure at low frequencies ( ∼ 10⁻³ day⁻¹) is not taken into account in the fitting. The average slopes and the corresponding confidence intervals are then calculated. The results are given in Table 2.

From Table 2, one can find that the PSDs produced by CARMA and celerite are consistent for all sources. Considering the uncertainties, there is no divergence between the PSDs of the light curves in different time bins. For the PSDs above 100 MeV, the slope is between −0.6 and −1.7, as shown in Figure 5. FSRQs and BL Lac objects have similar distributions of the PSD slopes (Figure 5). It should be stressed that the errors on the slopes are large for all sources.

Zoom In Zoom Out Reset image size

No QPO feature is found in all celerite and CARMA PSDs, except for PKS 0537−441 and PG 1553+113. We further analyze the posterior distributions of the QPO features in PKS 0537−441 and PG 1553+113 (Table 3 and Figures 6–7). According to the Nyquist criterion, the lower limit of the posterior distribution of the period is 60 days when the time bin is 30 days, and is 28 days when the time bin is 14 days. The calculation of the uncertainties on the period considers the whole posterior probability density distribution.

Zoom In Zoom Out Reset image size

Table 3.  Period, Quality Factor Q, and the Corresponding Significance Estimated by the Posterior Analyses of CARMA and Celerite Modeling for PKS 0537-441 and PG 1553+113

Source	Method	30 Days				14 Days
		Period	Q	L. S.	G. S.	Period	Q	L. S.	G. S.
		(days)		(%)	(%)	(days)		(%)	(%)
PKS 0537-441	CARMA	${233.42}_{-23.01}^{+27.21}$	${2.77}_{-1.10}^{+1.62}$	99.337	96.685	${280.57}_{-40.43}^{+63.44}$	${2.03}_{-0.79}^{+0.87}$	98.453	95.907
	celerite	${273.10}_{-13.52}^{+13.32}$	${9.20}_{-5.45}^{+25.27}$	99.911	99.603	${270.79}_{-198.32}^{+282.42}$	${6.70}_{-6.70}^{+225.70}$	93.608	49.264
PG 1553+113	CARMA	${840.25}_{-164.38}^{+337.60}$	${1.46}_{-0.77}^{+1.43}$	94.538	90.772	${101.94}_{-45.39}^{+134.36}$	${1.42}_{-0.82}^{+1.36}$	⋯	⋯
	celerite	${791.98}_{-289.75}^{+19.47}$	${62.47}_{-62.47}^{+64016.75}$	96.419	48.005	${780.81}_{-675.87}^{+30.37}$	${299.14}_{-299.14}^{+102933.42}$	98.161	52.380

Note. L. S. and G. S are the abbreviations of local significance and global significance, respectively. The uncertainties represent 1σ confidence intervals.

Download table as:  ASCIITypeset image

The period of PKS 0537−441 is constrained well at $\sim {233}_{-23}^{+27}$ days by CARMA in the 30 day light curve (Table 3 and Figure 6). Similarly, the period in the 14 day light curve is at $\sim {280}_{-40}^{+63}$ days. The two results are consistent within their errors. Celerite captures a QPO signal at $\sim {273}_{-14}^{+13}$ days in the 30 day light curve. For the 14 day light curve, a sign of a period of ∼270 days appears (Table 3 and Figure 6).

For PG 1553+113, a period of $\sim {840}_{-164}^{+338}$ days is found by CARMA in the 30 day binning light curve. Note that the uncertainties of the period are large. No QPO is captured in the 14 day light curve by CARMA (Figure 7 and Table 3). In the celerite model, the probability densities of the period obtained from two light curves both show a peak at ∼800 days (Figure 7); however, in both cases the period is not constrained (Table 3).

Zoom In Zoom Out Reset image size

4.3. Estimating the Significance of the QPOs in PKS 0537−441 and PG 1553+113

Some blazars' γ-ray light curves continuously monitored by Fermi-LAT are mainly behaving as stochastic processes (e.g., Sobolewska et al. 2014; Ryan et al. 2019), which can be explained by the internal shock model or turbulence in the jet (e.g., Rachen et al. 2010; Marscher 2014). However, the year-type periodic behaviors in the γ-ray band were recently claimed for blazars. If the periodicities are confirmed, these blazars will be an interesting sample for the investigation of accretion physics (e.g., McKinney et al. 2012), jet geometry (e.g., Rieger 2004), and the gravitational wave induced by the SMBH binary (e.g., Komossa 2006; Rieger 2007; Burke-Spolaor et al. 2019). However, the γ-ray QPOs of blazars are questioned (e.g., Covino et al. 2019; Ait Benkhali et al. 2020).

Stochastic models have been developed to characterize astronomical variability (e.g., Zu et al. 2013; Kelly et al. 2014; Foreman-Mackey et al. 2017; Li & Wang 2018). This kind of model assumes that the light curve is a realization of a Gaussian process. CARMA and celerite are two popular stochastic models for modeling astronomical variability. CARMA has been applied to γ-ray variabilities of blazars (Sobolewska et al. 2014; Goyal et al. 2018; Ryan et al. 2019).⁷ Celerite is a new tool in the study of AGN variability.

Sobolewska et al. (2014) applied their stochastic models to the 4 yr Fermi-LAT light curves of 13 blazars, and found that the γ-ray activity of the blazar is consistent with stochastic processes. Ryan et al. (2019) modeled the 9.5 yr Fermi-LAT light curves of the same 13 blazars with CARMA, and claimed that CARMA(2,1) models provide adequate descriptions of the variability. Covino et al. (2020) investigated the γ-ray QPOs in PG 1553+113 and PKS 2155−304 with the Gaussian process, and confirmed the QPO in 1553+113.

So far, possible year-type γ-ray QPOs are reported in 27 blazars. We analyze the 11 yr Fermi-LAT data of the 27 blazars and build their light curves. CARMA and celerite are applied to these light curves. Generally speaking, the γ-ray light curves can be well characterized by high-order CARMA and celerite models, which suggest that the γ-ray activity of the blazars is essentially a Gaussian process. The overfitting problem in several blazars can be resolved by modeling the logarithmic flux (see the Appendix).

PSDs are produced from the modeling results. The PSDs recovered by CARMA and celerite are consistent for the 27 sources with slopes from −0.6 to −1.7. For all sources, there is no significant difference between the PSDs of the light curves in different time bins. However, it is noted that the uncertainties on the PSDs are large.

The QPO feature only appears in the PSDs of PKS 0537−441 and PG 1553+113. Further posterior analyses are performed to examine the two periodicities. A period of ∼250 days for PKS 0537−441 is well constrained in the CARMA model. The modeling result of celerite gives a period of ∼280 days. It is noted that celerite result is only constrained well in one light curve. Our results are close to the period of ∼280 days claimed by Sandrinelli et al. (2016). The global significance of the period in our analysis is  >95% in the CARMA model. It is  >99% given by the celerite model with the 30 day binning light curve, while a low global significance is obtained in celerite modeling the 14 day binning light curve. The corresponding local significance is higher. In particular, the local significance is  >99% given by both CARMA and celerite models with the 30 day binning light curve. However, Covino et al. (2019) and Ait Benkhali et al. (2020) did not find QPO signal for this source. For PKS 0537−441, the slope of the PSD we obtained is consistent with that given by Covino et al. (2019). We note that a peak also appears at ∼280 days in the PSD of PKS 0537−441 given by Covino et al. (2019), although the signal is not significant. The estimate for the significance of such a QPO signal should be investigated more deeply .

In the modeling results of PG 1553+113, possible evidence for the period of ∼800 days appears, with local significance of  ≳ 95%. However, the global significance is 50%-90%. In the analyses of Covino et al. (2019) and Ait Benkhali et al. (2020), they did not find significant γ-ray QPO for PG 1553+113. While Covino et al. (2020) claimed to confirmed the period of the PG 1553+113 reported by Ackermann et al. (2015) with their Gaussian modeling method. The Gaussian model used in Covino et al. (2020) is similar to the celerite model. For PG 1553+113, a peak at ∼800 days can also be found in the posterior probability density functions given by the celerite model (Figure 7), but its global significance is considerably weakened by the distribution below 100 days. It seems that a broader posterior probability density distribution of period considered here causes different results from Covino et al. (2020).

Here, it is concluded that possible evidence for the γ-ray QPOs in PKS 0537−441 and PG 1553+113 are found in our Gaussian process analyses.

The γ-ray QPO of the blazar may be caused by a periodic modulation in Doppler factor (Yan et al. 2018), which could be the result of jet wobbling motion or jet procession. Jet wobbling motion is found for PG 1553+113 in the recent study of the radio emission from the jet on parsec scales (Lico et al. 2020). However, there is no QPO found in its radio emission, suggesting that the connection between jet motion and its emission is rather complex. It is hard to determine the origin of the γ-ray QPO.

The γ-ray QPO in PKS 0537−441 appears in the early years of the Fermi monitoring. There is an outburst along with the oscillation (Figure 1). Zhou et al. (2018) reported the detection of a γ-ray QPO after a γ-ray outburst in the blazar PKS 2247−13, and the QPO vanished after about six period cycles. Zhou et al. (2018) explained the QPO with periodic changes of Doppler factor. This scenario is also applicable to the origin of the QPO in PKS 0537−441. The method proposed by Yan et al. (2018) can be used to test this scenario, which will be carried out in a future work. On the other hand, the intrinsic origins for such a γ-ray QPO cannot be excluded, for instance, a periodic change in the particle acceleration rate in a binary SMBH system (e.g., Tavani et al. 2018).

Finally, we would like to stress that celerite models have been used to successfully characterize the γ-ray light curves of blazars, suggesting that celerite models have strong potential for studying AGN variability. Some issues are also found in the celerite modeling results, for example, the cause for the peaks at 40–50 days in the posterior distribution of period and the unconstrained Q, which likely cause the much larger uncertainties on the periods. A deep investigation is needed to resolve these issues, and to test numerical stability of celerite.

We thank the anonymous reviewer for constructive suggestions and Yan-Rong Li (IHEP) for valuable discussions. We acknowledge financial support from National Key R&D Program of China under grant No. 2018YFA0404204, the National Science Foundation of China (U1738124, 11803081, U1531131, and U1738211), and the Science Foundation of Yunnan Province (NO. 2018FA004). D.H.Y. acknowledges financial support from the joint foundation of Department of Science and Technology of Yunnan Province and Yunnan University [2018FY001(-003)]. The work of D.H.Y. is also supported by the CAS Youth Innovation Promotion Association and Basic research Program of Yunnan Province (202001AW070013).

Facility: Fermi (LAT) -

Software: Fermitools-conda, carma_pack (Kelly et al. 2014), celerite (Foreman-Mackey et al. 2017), emcee (Foreman-Mackey et al. 2013), NumPy, Matplotlib (Hunter 2007).

The overfitting problem in our results (Section 4.1) may be mitigated by reducing the weight of the flares in the whole variability. Here, we reanalyze the six sources (PKS 0250-225, PKS 0426-380, PKS 0454-234, S3 0458-02, S5 0716+714, and PKS 0805-077) using CARMA and celerite after taking logarithm for flux. The model selection and fitting analysis procedures are the same as that described in Section 4. The fitting results compared to Figure 1 are shown in Figure 8. One can see that all of the standardized residuals are consistent with the normal distribution. The PSDs of the six sources are also given in Figure 9. The PSDs are slightly different from those of nonlogarithmic flux at high frequencies. This is because modeling the light curve with logarithmic flux weakens the power of the flares and then change, the PSD shape at high frequencies (Ryan et al. 2019).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size