Evaluating Quasi-Periodic Variations in the $\gamma$-ray Lightcurves of Fermi-LAT Blazars

TL;DR

This study analyzes gamma-ray light curves of six blazars over nearly a decade to assess the presence of quasi-periodic variations, finding that apparent periodicities are not statistically significant when accounting for random variability.

Contribution

The paper provides an updated, rigorous timing analysis of multiple blazars using advanced statistical methods to evaluate the significance of quasi-periodic variations in gamma-ray light curves.

Findings

Possible year-like QPVs observed with ~3 sigma significance.

Most blazars exhibit variability consistent with log-normal distribution.

Observed periods cluster around 1.6 years, but significance is reduced globally.

Abstract

The detection of periodicities in light curves of active galacticnuclei (AGN) could have profound consequences for our understanding of the nature and radiation physics of these objects. At high energies (HE; E>100 MeV) 5 blazars (PG 1553+113,PKS 2155-304, 0426-380, 0537-441, 0301-243) have been reported to show year-like quasi-periodic variations (QPVs) with significance >3 sig. As these findings are based on few cycles only, care needs to be taken to properly account for random variations which can produce intervals of seemingly periodic behaviour. We present results of an updated timing analysis for 6 blazars (adding PKS 0447-439), utilizing suitable methods to evaluate their long term variability properties and to search for QPVs in their light curves. We generate gamma-ray light curves covering almost 10 years, study their timing properties and search for QPVs using the Lomb-Scargle Periodogram and the Wavelet Z-transform. Extended Monte Carlo simulations are used to evaluate the statistical significance. Comparing their probability density functions (PDFs), all sources (except PG 1553+113) exhibit a clear deviation from a Gaussian distribution, but are consistent with being log-normal, suggesting that the underlying variability is of a non-linear, multiplicative nature. Apart from PKS 0301-243 the power spectral density for all investigated blazars is close to flicker noise (PL slope -1). Possible QPVs with a local significance ~ 3 sig. are found in all light curves (apart from PKS 0426-380 and 0537-441), with observed periods between (1.7-2.8) yr. The evidence is strongly reduced, however, if evaluated in terms of a global significance. Our results advise caution as to the significance of reported year-like HE QPVs in blazars. Somewhat surprisingly, the putative, redshift-corrected periods are all clustering around 1.6 yr. We speculate on possible implications for QPV generation.

Figures

Tables

Table 1. Table 1: Parameters obtained for the fit of the Fermi-LAT energy spectrum of PKS 0447-439 using Power Law (PL), Log Parabola (LP), Broken Power Law (BPL) and PLSuperExpCutOff (PLExp) spectral model. The last column gives the significance, obtained using Log-Likelihood Test by comparing the log-Likelihood values for each model against those for the power-law model. Only statistical errors are shown. A BPL is preferred here. N pred. represents the number of predicted photons.

Model	Energy	Spectral Parameters	Npred.	Integral Flux	TS	Log-	Log-Likelihood
	[GeV]			${10}^{-8}\left[c{m}^{-2}{s}^{-1}\right]$		Likelihood	Ratio-Test
Power Law (PL)	$0.1$ .. $500$	$\Gamma =1.867\pm 0.009$	$16724$	$\left(7.83\pm 0.13\right)$	$29293$	-$402489.3$	$-$
Log Parabola	$0.1$ .. $500$	${\Gamma }_{\alpha }=1.742\pm 0.025$	$15726$	$\left(7.39\pm 0.15\right)$	$29234$	-$402479.8$	$19.00$
(LP)		${\Gamma }_{\beta }=0.020\pm 0.004$					$\left(4.4\sigma \right)$
Broken PowerLaw	$0.1$ .. $500$	${\Gamma }_1=1.839\pm 0.009$	$16234$	$\left(7.62\pm 0.13\right)$	$29305$	-$402462.4$	$53.80$
(BPL)		${\Gamma }_2=2.825\pm 0.171$					$\left(7.3\sigma \right)$
		$E_b=42.84\pm 1.78$
PLSuperExpCutOff	$0.1$ .. $500$	${\Gamma }_1=1.818\pm 0.012$	$16119$	$\left(7.51\pm 0.14\right)$	$29282$	-$402464.2$	$50.20$
		${\Gamma }_2=1.031\pm 0.195$					$\left(7.1\sigma \right)$
		$E_c=137.9\pm 24.3$

Table 2. Table 2: The spectral parameters obtained for the best fit of Fermi-LAT data of each sources from the energy range between 100 100 100 MeV and 500 500 500 GeV using four different spectral models PowerLaw, Broken PowerLaw, LogParabola and PLSuperExpCutoff. The forelast column gives the significance, obtained by comparing the likelihood values for each models against those for the PowerLaw model using Log likelihood ratio test. The last column summarizes the fractional variability amplitude F var calculated form monthly light curves (Vaughan et al. 2003 ; Edelson et al. 2002 ) . Error bars correspond to 1 ​ σ 1 𝜎 1\sigma error; only statistical errors are shown.

Source	Model	Energy	Spectral Parameters	Integral Flux	TS	Ratio-	Fvar
		[GeV]		${10}^{-8}\left[c{m}^{-2}{s}^{-1}\right]$		Test	[%]
PKS 0447-439	Broken PowerLaw	$0.1$ .. $500$	${\Gamma }_1=1.839\pm 0.009$	$\left(7.62\pm 0.13\right)$	$29305.1$	$7.3\sigma$	$\left(47.91\pm 1.16\right)$
			${\Gamma }_2=2.825\pm 0.171$
			$E_b=42.84\pm 1.78$
PG 1553+113	Log Parabola	$0.1$ .. $500$	${\Gamma }_{\alpha }=1.587\pm 0.013$	$\left(5.59\pm 0.15\right)$	$34701.6$	$11.01\sigma$	$\left(22.87\pm 1.14\right)$
			${\Gamma }_{\beta }=0.045\pm 0.005$
			$E_b=1.49\pm 0.95$
PKS 2155-304	PLSuperExpCutOff	$0.1$ .. $500$	${\Gamma }_1=1.749\pm 0.154$	$\left(12.82\pm 0.22\right)$	$59800.5$	$9.63\sigma$	$\left(36.62\pm 1.07\right)$
			${\Gamma }_2=0.620\pm 0.154$
			$E_c=110.57\pm 37.9$
PKS 0426-380	Log Parabola	$0.1$ .. $500$	${\Gamma }_{\alpha }=1.742\pm 0.025$	$\left(24.05\pm 0.31\right)$	$92892.5$	$21.19\sigma$	$\left(62.39\pm 0.01\right)$
			${\Gamma }_{\beta }=0.020\pm 0.004$
			$E_b=0.76\pm 0.03$
PKS 0301-243	Log Parabola	$0.1$ .. $500$	${\Gamma }_{\alpha }=1.622\pm 0.043$	$\left(3.23\pm 0.11\right)$	$9639.4$	$5.94\sigma$	$\left(82.28\pm 1.90\right)$
			${\Gamma }_{\beta }=0.039\pm 0.006$
			$E_b=0.11\pm 0.02$
PKS 0537-441	PLSuperExpCutOff	$0.1$ .. $500$	${\Gamma }_1=1.662\pm 0.026$	$\left(19.38\pm 0.20\right)$	$63570.9$	$16.78\sigma$	$\left(85.13\pm 0.08\right)$
			${\Gamma }_2=0.267\pm 0.009$
			$E_c=0.24\pm 0.06$

Table 3. Table 3: The table shows the best Fit parameters for the log-normal and normal (Gaussian) distribution, respectively; the AD Test statistics and χ r ​ e ​ d . 2 subscript superscript 𝜒 2 𝑟 𝑒 𝑑 \chi^{2}_{red.} are the Anderson Darling test and the reduced χ 2 superscript 𝜒 2 \chi^{2} , respectively. [*] in units of photon/cm s − 1 2 superscript superscript 𝑠 1 2 {}^{2}s^{-1} .

Source	Log-normal				Gaussian
Name	Center*	Width	${\chi }_{red.}^2$	AD Test	Center*	Width	${\chi }_{red.}^2$	AD Test
	[ $\times {10}^{-7}$ ]	[ $\times {10}^{-8}$ ]		statistics	[ $\times {10}^{-8}$ ]	[ $\times {10}^{-8}$ ]		statistics
PG 1553+113	$\left(1.41\pm 0.02\right)$	$\left(2.13\pm 0.13\right)$	$14.59$	$2.57$	$\left(1.40\pm 0.02\right)$	$\left(3.21\pm 0.15\right)$	$12.29$	$1.83$
PKS 0447-439	$\left(0.97\pm 0.01\right)$	$\left(2.48\pm 0.15\right)$	$5.70$	$0.18$	$\left(0.93\pm 0.02\right)$	$\left(3.66\pm 0.16\right)$	$9.94$	$6.61$
PKS 2155-304	$\left(1.65\pm 0.02\right)$	$\left(3.65\pm 0.23\right)$	$8.00$	$0.35$	$\left(1.58\pm 0.02\right)$	$\left(5.14\pm 0.22\right)$	$14.38$	$8.76$
PKS 0426-380	$\left(2.19\pm 0.09\right)$	$\left(8.99\pm 0.58\right)$	$13.29$	$0.79$	$\left(1.91\pm 0.08\right)$	$\left(12.42\pm 0.89\right)$	$14.02$	$7.84$
PKS 0301-243	$\left(0.65\pm 0.01\right)$	$\left(1.57\pm 0.09\right)$	$6.56$	$0.77$	$\left(0.62\pm 0.01\right)$	$\left(2.37\pm 0.08\right)$	$17.52$	$22.52$
PKS 0537-441	$\left(1.28\pm 0.04\right)$	$\left(4.62\pm 0.29\right)$	$16.30$	$11.13$	$\left(1.02\pm 0.09\right)$	$\left(7.82\pm 0.94\right)$	$35.99$	$26.41$

Table 4. Table 4: D’Agostino’s K 2 Test tests the null hypothesis that a sample comes from a normal distribution. It is based on D’Agostino and Pearson’s test (D’Agostino 1970 ; D’Agostino & Pearson 1973 ) and combines skew and kurtosis to produce an omnibus test of normality. Similarly, the Shapiro-Wilk test quantifies how likely it is that the data was drawn from a Gaussian distribution (Shapiro & Wilk 1965 )

Source	K2	Pvalue	Shapi.-	Pvalue
	Stati.	(${\chi }^2$ Prob.)	Wilk	(Prob.)
PG 1553+113	$56.8$	$4.6\times {10}^{-13}$	$0.972$	$4.7\times {10}^{-8}$
PKS 0447-439	$101.3$	$1.0\times {10}^{-22}$	$0.935$	$9.9\times {10}^{-14}$
PKS 2155-304	$222.1$	$5.8\times {10}^{-49}$	$0.893$	$6.4\times {10}^{-18}$
PKS 0426-380	$90.7$	$2.0\times {10}^{-20}$	$0.927$	$1.2\times {10}^{-14}$
PKS 0301-243	$721.4$	$2.2\times {10}^{-157}$	$0.607$	$8.1\times {10}^{-32}$
PKS 0537-441	$185.4$	$5.6\times {10}^{-41}$	$0.815$	$2.9\times {10}^{-23}$

Table 5. Table 5: Periodicity results derived from the GLSP method for the monthly and weekly γ 𝛾 \gamma -ray light curves, respectively, and including the intrinsic (redshift-corrected) period. [*] in units of days.

Source	Energy	Redshift	Observed Period		Intrinsic Period
	[ GeV ]	$z$	Monthly bins*	Weekly bins*	Monthly bins*	Weekly bins*
PG 1553+113	$0.1$ .. $500$	$0.360$	($821.9\pm 74.5$)	($811.9\pm 58.4$)	($604.3\pm 54.8$)	($596.9\pm 42.9$)
PKS 0447-439	$0.1$ .. $500$	$0.343$	($929.0\pm 88.4$)	($935.4\pm 81.0$)	($691.7\pm 65.8$)	($696.5\pm 60.3$)
PKS 2155-304	$0.1$ .. $500$	$0.116$	($619.7\pm 40.8$)	($626.7\pm 35.6$)	($555.3\pm 35.8$)	($561.5\pm 31.9$)
PKS 0426-380	$0.1$ .. $500$	$1.111$	($1218.0\pm 158.7$)	($1225.8\pm 146.9$)	($576.9\pm 75.2$)	($580.7\pm 69.6$)
PKS 0301-243	$0.1$ .. $500$	$0.260$	($755.1\pm 58.4$)	($760.5\pm 90.5$)	($599.3\pm 46.3$)	($603.6\pm 71.8$)

Table 6. Table 6: Summary results of the strongest periods for the Weighted Wavelet Z-transform method WWZ based on the monthly γ 𝛾 \gamma -ray light curves, including observed and intrinsic period.

Source Name	$z$	Observed Period	Intrinsic Period
		[ days ]	[ days ]
PG 1553+113	$0.360$	$807.8\pm 102.4$	$593.9\pm 75.3$
PKS 0447-439	$0.343$	$937.3\pm 152.7$	$697.9\pm 113.7$
PKS 2155-304	$0.116$	$624.7\pm 69.5$	$559.8\pm 62.3$
PKS 0426-380	$1.111$	$1257.9\pm 280.4$	$595.9\pm 132.8$
PKS 0301-243	$0.260$	$752.4\pm 89.2$	$597.1\pm 70.8$

Table 7. Table 7: The table shows the slopes β 𝛽 \beta obtained from fitting the γ 𝛾 \gamma -ray power spectra with simple power-law model, and the correspondents success fraction calculated based on the PSRESP method (Uttley et al. 2002 ) (see subsection 4.1 4.1 4.1 ) (the success fractions indicate the goodness of fit obtained from the PSRESP method). The β 𝛽 \beta values have been determined using both monthly and weekly light curves, as well as the Timmer & Koenig (1995) and Emmanoulopoulos et al. (2013) simulation methods. The errors are results from fitting Gaussians to each of the slope profile.

Source	Simulation Method		Monthly		Weekly
	PSRESP	Simulation	PSD-Slope	Success Fraction	PSD-Slope	Success Fraction
	Model	Algorithm	${\beta }_m$	[ % ]	${\beta }_w$	[ % ]
PG 1553+113	Power Law	Timmer & Koenig	$0.88\pm 0.13$	$96.0\pm 3.2$	$0.66\pm 0.11$	$82.4\pm 3.5$
	Power Law	Emmanoulopoulos	$0.90\pm 0.14$	$93.0\pm 2.9$	$0.67\pm 0.11$	$84.8\pm 3.6$
PKS 0447-439	Power Law	Timmer & Koenig	$1.06\pm 0.15$	$98.4\pm 3.0$	$0.68\pm 0.10$	$88.3\pm 3.8$
	Power Law	Emmanoulopoulos	$1.07\pm 0.13$	$83.5\pm 5.2$	$0.69\pm 0.10$	$86.6\pm 4.1$
PKS 2155-304	Power Law	Timmer & Koenig	$0.95\pm 0.13$	$97.4\pm 3.1$	$0.88\pm 0.09$	$86.9\pm 3.3$
	Power Law	Emmanoulopoulos	$0.99\pm 0.14$	$94.4\pm 2.7$	$0.89\pm 0.09$	$90.5\pm 3.6$
PKS 0426-380	Power Law	Timmer & Koenig	$0.98\pm 0.22$	$95.1\pm 5.0$	$1.08\pm 0.08$	$96.5\pm 3.6$
	Power Law	Emmanoulopoulos	$1.02\pm 0.23$	$88.6\pm 4.1$	$1.09\pm 0.09$	$94.2\pm 3.7$
PKS 0301-243	Power Law	Timmer & Koenig	$0.32\pm 0.20$	$96.9\pm 3.6$	$0.53\pm 0.14$	$94.8\pm 5.3$
	Power Law	Emmanoulopoulos	$0.36\pm 0.22$	$95.1\pm 2.9$	$0.58\pm 0.15$	$93.9\pm 3.2$

Table 8. Table 8: The table shows the significance values of the observed periods obtained for the two method local and global significance evaluated using the Timmer & Koenig (1995) and the Emmanoulopoulos et al. (2013) and based on on the monthly and weekly γ 𝛾 \gamma -ray light curves. the slopes β 𝛽 \beta are obtained from fitting the γ 𝛾 \gamma -ray power spectra with simple power-law model, and the corresponding success fraction are calculated based on the PSRESP method (Uttley et al. 2002 ) .

Source	Simulation	Monthly			Weekly
	Algorithm	PSD-Slope	local	global	PSD-Slope	local	global
		$⟨{\beta }_m⟩$	Significance	Significance	$⟨{\beta }_w⟩$	Significance	Significance
PG 1553+113	Timmer & Koenig	$0.88$	$3.89$ $\sigma$	$2.14$ $\sigma$	$0.66$	$3.70$ $\sigma$	$2.07$ $\sigma$
	Emmanoulopoulos	$0.90$	$3.82$ $\sigma$	$1.73$ $\sigma$	$0.67$	$3.83$ $\sigma$	$2.16$ $\sigma$
PKS 0447-439	Timmer & Koenig	$1.06$	$2.87$ $\sigma$	$1.00$ $\sigma$	$0.68$	$3.26$ $\sigma$	$1.49$ $\sigma$
	Emmanoulopoulos	$1.07$	$2.94$ $\sigma$	$1.16$ $\sigma$	$0.69$	$3.16$ $\sigma$	$1.46$ $\sigma$
PKS 2155-304	Timmer & Koenig	$0.95$	$3.62$ $\sigma$	$1.61$ $\sigma$	$0.89$	$3.32$ $\sigma$	$0.84$ $\sigma$
	Emmanoulopoulos	$0.99$	$3.72$ $\sigma$	$1.64$ $\sigma$	$0.88$	$3.16$ $\sigma$	$0.77$ $\sigma$
PKS 0426-380	Timmer & Koenig	$0.98$	$2.73$ $\sigma$	$1.04$ $\sigma$	$1.08$	$2.17$ $\sigma$	$0.52$ $\sigma$
	Emmanoulopoulos	$1.02$	$2.72$ $\sigma$	$0.96$ $\sigma$	$1.09$	$2.35$ $\sigma$	$0.18$ $\sigma$
PKS 0301-243	Timmer & Koenig	$0.32$	$3.54$ $\sigma$	$1.77$ $\sigma$	$0.53$	$3.39$ $\sigma$	$1.02$ $\sigma$
	Emmanoulopoulos	$0.36$	$3.72$ $\sigma$	$2.17$ $\sigma$	$0.58$	$3.43$ $\sigma$	$1.43$ $\sigma$

Table 9. Table 9: The table summarizes the significance values of the observed periods for the local and global method based on the monthly and weekly γ 𝛾 \gamma -ray light curves of PG 1553+113 and PKS 0447-489, respectively. The significance levels are evaluated via MC-Simulation using Timmer & Koenig 1995 and Emmanoulopoulos et al 2013 algorithms with 50000 trials. The mean slopes ⟨ β m ⟩ delimited-⟨⟩ subscript 𝛽 𝑚 \langle\beta_{m}\rangle are obtained from fitting the γ 𝛾 \gamma -ray power spectra with simple power-law model, and the corresponding success fraction are calculated based on the PSRESP method (Uttley et al. 2002 ) . β m ​ i ​ n subscript 𝛽 𝑚 𝑖 𝑛 \beta_{min} and β m ​ i ​ n subscript 𝛽 𝑚 𝑖 𝑛 \beta_{min} represents the minimum and maximum of the 1 ​ σ 1 𝜎 1\sigma error bars resulting from fitting Gaussians to each of the slope profile.

Source	Simulation	Monthly			Weekly
	Algorithm	PSD-Slope	local	global	PSD-Slope	local	global
		${\beta }_m$	Significance	Significance	${\beta }_w$	Significance	Significance
PG 1553+113	Timmer & Koenig	${\beta }_{min}=0.75$	$5.33$ $\sigma$	$2.46$ $\sigma$	${\beta }_{min}=0.55$	$5.32$ $\sigma$	$3.06$ $\sigma$
		$⟨{\beta }_m⟩=0.88$	$3.89$ $\sigma$	$2.14$ $\sigma$	$⟨{\beta }_w⟩=0.66$	$3.70$ $\sigma$	$2.07$ $\sigma$
		${\beta }_{max}=1.01$	$3.72$ $\sigma$	$1.86$ $\sigma$	${\beta }_{max}=0.77$	$2.90$ $\sigma$	$0.93$ $\sigma$
	Emmanoulopoulos	${\beta }_{min}=0.76$	$5.33$ $\sigma$	$2.40$ $\sigma$	${\beta }_{min}=0.56$	$5.32$ $\sigma$	$3.10$ $\sigma$
		$⟨{\beta }_m⟩=0.90$	$3.89$ $\sigma$	$2.06$ $\sigma$	$⟨{\beta }_w⟩=0.67$	$3.83$ $\sigma$	$2.16$ $\sigma$
		${\beta }_{max}=1.04$	$3.48$ $\sigma$	$1.75$ $\sigma$	${\beta }_{max}=0.78$	$2.79$ $\sigma$	$0.79$ $\sigma$
PKS 0447-439	Timmer & Koenig	${\beta }_{min}=0.91$	$3.12$ $\sigma$	$1.06$ $\sigma$	${\beta }_{min}=0.58$	$3.48$ $\sigma$	$2.01$ $\sigma$
		$⟨{\beta }_m⟩=1.06$	$2.87$ $\sigma$	$1.00$ $\sigma$	$⟨{\beta }_w⟩=0.68$	$3.26$ $\sigma$	$1.49$ $\sigma$
		${\beta }_{max}=1.21$	$2.67$ $\sigma$	$0.46$ $\sigma$	${\beta }_{max}=0.78$	$2.78$ $\sigma$	$1.02$ $\sigma$
	Emmanoulopoulos	${\beta }_{min}=0.94$	$3.11$ $\sigma$	$1.04$ $\sigma$	${\beta }_{min}=0.59$	$3.54$ $\sigma$	$1.99$ $\sigma$
		$⟨{\beta }_m⟩=1.07$	$2.94$ $\sigma$	$1.16$ $\sigma$	$⟨{\beta }_w⟩=0.69$	$3.16$ $\sigma$	$1.46$ $\sigma$
		${\beta }_{max}=1.20$	$2.71$ $\sigma$	$0.54$ $\sigma$	${\beta }_{max}=0.79$	$2.80$ $\sigma$	$1.00$ $\sigma$

Equations

$$\displaystyle L(\phi\textbf{ }|\textbf{ }\mu,\sigma)=\frac{1}{\sqrt{2\pi}\sigma\phi}exp\textbf{ }\bigg{(}-\frac{\big{(}\log(\phi)-\mu\big{)}^{2}}{2\sigma^{2}}\bigg{)}$$

$$\displaystyle G(\phi\textbf{ }|\textbf{ }\mu,\sigma)=\frac{1}{\sqrt{2\pi}\sigma}exp\textbf{ }\bigg{(}-\frac{\big{(}\phi-\mu\big{)}^{2}}{2\sigma^{2}}\bigg{)}\,,$$

$$P(\omega )$$

$$\tan (\tau )$$

$$\text{ $y(t)=\sum_{i}y_{i}\phi_{i}(t)$ :}\qquad\left\{\begin{aligned} \phi_{1}(t)=1(t)=1\\ \phi_{2}(t)=cos\textbf{ }\big{(}\omega(t-\tau)\big{)}\\ \phi_{3}(t)=sin\textbf{ }\big{(}\omega(t-\tau)\big{)}\end{aligned}\right.$$

$${\omega }_i=e^{-c\omega (t_i-\tau {)}^2},$$

$$WWZ=\frac{(N_{eff}-3)V_y}{2(V_x-V_y)},$$

$$P(\nu)=A\bigg{(}\frac{\nu}{\nu_{0}}\bigg{)}^{-\beta}+C_{noise}\,,$$

$$\displaystyle{\chi^{2}_{dist,i}}={\sum\limits_{\nu=\nu_{min}}^{\nu_{max}}}\frac{\big{[}P(\nu)_{sim,i}-\langle P(\nu)_{sim}\rangle\big{]}^{2}}{\langle\Delta P(\nu)_{sim}\rangle^{2}}\,,$$

$$\displaystyle\ {\chi^{2}_{obs}}={\sum\limits_{\nu=\nu_{min}}^{\nu_{max}}}\frac{\big{[}P(\nu)_{obs}-\langle P(\nu)_{sim}\rangle\big{]}^{2}}{\langle\Delta P(\nu)_{sim}\rangle^{2}}\,,$$

Full text

*Aims. *We present results of an updated timing analysis for six blazars (adding PKS 0447-439 to the above), utilizing suitable methods to evaluate their longterm variability properties and to search for QPVs in their light curves.

*Methods. *We generate $\gamma$-ray light curves covering almost ten years, study their timing properties and search for QPVs using the Lomb-Scargle Periodogram and the Wavelet Z-transform. Extended Monte Carlo simulations are used to evaluate the statistical significance.

*Results. * (1) Comparing their probability density functions (PDFs), all sources (except PG 1553+113) exhibit a clear deviation from a Gaussian distribution, but are consistent with being log-normal, suggesting that the underlying variability is of a non-linear, multiplicative nature. (2) Apart from PKS 0301-243 the power spectral density (PSD) for all investigated blazars is close to flicker noise (power-law slope $-1$). (3) Possible QPVs with a local significance \lower 2.15277pt\hbox{;\buildrel>\over{\sim};}3\sigma are found in all light curves (apart from PKS 0426-380 and PKS 0537-441), with observed periods between $(1.7-2.8)$ yr. The evidence is strongly reduced, however, if evaluated in terms of a global significance.

*Conclusions. * Our results advise caution as to the significance of reported year-like HE QPVs in blazars. Somewhat surprisingly, the putative, redhift-corrected periods are all clustering around $\sim 1.6$ yr. We speculate on possible implications for QPV generation.

1 Introduction

Blazars belong to the most luminous and variable extragalactic sources in the Universe. They represent a special sub-class of AGN characterized by a relativistic plasma jet oriented very close to the line of sight at angles $\theta\leq 1/\Gamma$, where $\Gamma$ is the bulk Lorentz factor ( (e.g., Urry & Padovani 1995; Ackermann et al. 2015). Blazar sources are known to be variable across all wavelengths from radio frequencies to TeV $\gamma$-ray energies, and on a wide range of timescales from sub-minute to several years (e.g. Böttcher & Chiang 2002; Aharonian et al. 2007). Their multi-wavelength spectral energy distribution (SED) often exhibits two bumps. The first, low-energy bump (peaking at infrared to X-ray frequencies) is usually interpreted as synchrotron emission from highly relativistic electrons within in the jet. The second bump, on the other hand, peaking in the HE range has been frequently related to inverse Compton emission (IC) up-scattering of various soft photon fields (mostly synchrotron or external thermal radiation), although hadronic interactions may contribute as well.

Periodic variability in the light curves of blazars has been investigated extensively in the radio and optical band (e.g., Fan et al. 2002; Kadler et al. 2006; King et al. 2013; Bhatta et al. 2016). In this context the two prominent sources OJ 287 and 3C 279 are worth mentioning, with longterm optical periods of $\sim 12$ and $29.6$ years being reported, respectively (Valtonen et al. 2006; Li et al. 2009). In the high-energy gamma-ray regime only five blazars (PG 1553+113, PKS 2155-304, PKS 0301-243, PKS 0426-380 and PKS 0537-441) have been reported so far exhibiting longterm quasi-periodic variability in their $\gamma$-ray fluxes with significance higher than $3\sigma$ (Ackermann et al. 2015; Sandrinelli et al. 2014; Zhang et al. 2017, 2017; Sandrinelli et al. 2016; Prokhorov & Moraghan 2017). Apart from PKS 0537-441, these QPVs have been seen over the whole length of the light curve. In the majority of these cases, attempts have been made to relate the detected periodic behaviour at different wavelengths and on various timescales to the motion of two SMBHs in a binary black holes system, general helical jet structures, shocks or instabilities of the disk or jet-plasma flow (e.g., Rieger 2004; Wiita 2011; Sobacchi et al. 2017).

In the present paper we use Fermi-LAT $\gamma$-ray data between $08.2008$ and $12.2017$ to re-evaluate the five blazars mentioned above and to investigate an additional source, PKS 0449-439. Results of a periodicity search based on the generalized Lomb-Scargle periodogram (GLSP) and the Weighted Wavelet Z-transform (WWZ) are presented. The signifiance of the inferred periods is then estimated on the basis of Bootstrap resampling and light curves simulations using the Timmer$\&$ Koenig as well as the Emmanoulopoulos algorithm (Timmer & Koenig 1995; Emmanoulopoulos et al. 2013)

The paper is organized as follows. Section 2 provides an exemplary illustration of the performed Fermi-LAT data analysis procedures to generate $\gamma$-ray light curves based on a binned likelihood method and aperture photometry technique, respectively. In Sec. 3 we describe the GLSP and WWZ technique used here to search for periodicity. Signifiance estimation and variability characterization is presented in Sec. 4. Finally, Sec. 5 provides a summary and discussion of the results.

2 Observations and Analysis

2.1 Fermi-LAT likelihood analysis

Fermi-LAT on board the Fermi satellite is a electron-positron pair-conversion $\gamma$-ray detector sensitive to photon in energy range between $\sim 20$ MeV and $500$ GeV. The LAT has a wide field of view of $\sim 2.4$ sr and observes the entire sky every 2 orbits. It has a point spread function (PSF) $<0.8^{\circ}$ and the largest effective area of approximately $8000$ cm2 above 1 GeV. Fermi-LAT data were downloaded from the Fermi-LAT Data Server111https://fermi.gsfc.nasa.gov/ssc/data/access/ and events were selected between $2008$ August $4$ and $2018$ April $1$ (MET from $239557418.0$ to $541555205.0$), covering about $\sim 9.5$ years with energy range from $100$ MeV to $500$ GeV in a circular region of interest (ROI) of $15$ degree centered on the position of each source (Atwood et al. 2013).

A binned maximum likelihood analysis (Mattox et al. 1996) is performed using the Pass 8 (P8R2) algorithms and employing the standard Fermi Science Tools v10r0p5 software package. We follow the standard procedure provided by the Fermi Science Support Center (FSSC) to reduce the data. In addition, photons coming from zenith angles larger than $90$ degree were all rejected to reduce the background from gamma rays produced in the atmosphere of the Earth (albedo) with P8R2_SOURCE_V6 instrument response functions (IRF). We performed standard quality cuts in accordance with the Pass8 data analysis criteria. The background emission was modeled using the Galactic and isotropic diffuse emission gll_iem_v06_iso and P8R2_SOURCE_V6_v06 files two. All sources from the third LAT source catalog (3FGL; (Acero et al. 2015) within the ROI are included in the model to ensure a satisfactory background modeling.

A binned-likelihood analysis with a bin size of $0.1^{\circ}$ was performed. The spectral indices were allowed to vary for sources located within a radius of $5^{\circ}$ around the position of the investigated source. The Test Statistic value (TS), defined as TS$=-2\ln(L_{0}/L_{1})$ was used to determine the source detection significance, with threshold set to TS$=25$ ($\sim 5\sigma$). The significance of a source detection is given by $\sim\sqrt[]{TS}\sigma$ (Abdo et al. 2010). In the case of PKS 0447-439 for example, the likelihood analysis reveals a point source with a high statistical significance TS $\simeq 29305$ corresponding to $171\sigma$. Both energy and temporal bins with TS$<9$ (or $\sim 3\sigma$) are set as upper limits throughout the paper.

The appearance of new sources within the ROI could in principle strongly influence the $\gamma$-ray spectrum and light curve of the investigated source. In order to investigate this we also searched for possible new $\gamma$-ray sources within a FOV of $15^{\circ}$. Our analysis of $9.5$ year of Fermi-LAT data indicates additional transient sources beyond the 3FGL catalog; these were appropriately taken into account in the full analysis, i.e. for each additional source we sequentially add a new point source with a standard spectral definition (PowerLaw) and maximize the likelihood as a function of its flux.

We produced the $\gamma$-ray SED and light curves of each sources through the binned maximum likelihood fitting technique, respectively with gtlike to determine the flux and TS value for each time bin. The effect of energy dispersion below $300$ MeV is also accounted for in the analysis. For that we enabled the energy dispersion correction.

As an example, the SED of PKS 0447-439 is shown in Fig. 1 for the full time and energy range using 10 energy bins. We performed maximum likelihood analyses, exploring different spectral forms for the whole energy range, namely Single Power-Law (PL), Broken Power-Law (BPL), Log-Parabola (LP) and Power-Law Super Exponential Cutoff (PLExpCutOff). In the case of PKS 0447-439 a likelihood ratio test comparison yields TS$=-2\log(L_{PL}/L_{BPL})\simeq 53.8$, thus preferring BPL over PL at a significance level of $7.3\ \sigma$. The best fit indices are $\Gamma_{1}=1.84\pm 0.01$ and $\Gamma_{2}=2.82\pm 0.17$ with a break at an energy $E_{b}=(42.8\pm 1.8)$ GeV. The obtained results of the maximum likelihood analysis are summarized for the different models in the Table 1.

2.2 Fermi-LAT $\gamma$-ray light curves

We generate $\gamma$-ray light curves for all our six sources using the individual best fit results (i.e., a BPL model in the case of PKS 0447-439), see Table 2 for details. For comparison two methods are employed: the maximum likelihood optimization and the aperture photometry. In the analysis with the first method the $\gamma$-ray light curves are produced by using separate equal time-bins of one month (30 days). We apply a maximum likelihood fitting technique by running the ScienceTools gtlike to extract the flux and TS value for each time bin. In this computationally more intensive procedure, the $\gamma$-events are selected in a circular ROI of $15^{\circ}$ radius centered on the position of the source. The resulting light curve for PKS 0447-439 is shown in Fig. 2 and has an average $\gamma$-ray flux of $\langle{\phi}\rangle=(7.62\pm 0.13)\times 10^{-8}$cm2s*-1*. The source exhibits flux variability during the whole observational period, with flux levels of peak-to-peak oscillations changing by a factor four. The calculated fractional variability index is F${}_{var}=(47.91\pm 1.16)\%$ (Vaughan et al. 2003; Edelson et al. 2002).

The aperture photometry method on the other hand, provides a model-independent measure of the $\gamma$-ray flux and is less computationally demanding. It also enables the use of short time bins whereas the maximum likelihood technique requires that time bins contain sufficient photons for the analysis. We select a very small aperture radius of $1^{\circ}$ to exclude most background $\gamma$-events and to focus on the events that are most likely associated with the source itself. We produce $\gamma$-ray light curves in the range between $100$ MeV and $500$ GeV using a weekly binning. The exposure of each time bin is determined with the ScienceTools gtexposure. Figure  3 shows the resultant $\gamma$-ray light curve for PKS 0447-439 obtained with this method.

In order to verify that the aperture photometry (AP) approach gives results which are consistent with the results obtained using the binned likelihood (BL) analysis, we performed diagnostic tests on each of the investigated blazars. The BL method is the preferred procedure for most types of Fermi-LAT analyses. Taking for granted that BL flux values are indeed more accurately measured than the AP flux, we examine the correlation between the AP and BL calculated flux to assess the performance of the AP method. This correlation is plotted for PKS 0447-439 in Fig. 4 using monthly binned light curves and the calculated value of the Pearson correlation coefficient is $\rho\sim 0.936$.

2.3 Gaussian vs Log-normal Flux Distributions

The almost continuous detection of the investigated sources over time bins of seven days provides a good opportunity to investigate whether the probability density function (PDF) for the long-term $\gamma$-ray emission reveals any preference for a gaussian (normal) or a log-normal flux distribution. It is expected that such a feature offers an important clue for understanding the central engine of a blazer (e.g., Uttley et al. 2005; Shah et al. 2018). A log-normal distribution of $\gamma$-ray fluxes, for example, would suggest that the mechanism driving the variability is multiplicative in nature rather than additive. Such multiplicative behaviour might possibly be related to accretion disk fluctuations (Lyubarskii 1997; Arévalo & Uttley 2006) or the particle acceleration process itself (Shah et al. 2018). As a possible caveat we note that PDF of stochastic processes with power spectral density (PSD) steeper than index $-1$ can show deviations from normality owing to divergence of power at low frequencies (e.g., Vaughan et al. 2003).

We explore the PDF and quantify this in terms of fluxes using flux-histograms. For each source the weekly $\gamma$-ray fluxes were distributed in a histogram of fluxes. We fit all histograms in log-scale, with Gaussian G($\phi$) and log-normal L($\phi$) distribution functions given by:

[TABLE]

and

[TABLE]

respectively, where $\sigma$ and $\mu$ are the standard deviation and the mean of the distribution, respectively. For illustration the obtained flux histograms for PKS 0447-439, PG 1553+113 and PKS 2155-304 are shown in Figure  5. We compute the Anderson Darling test (AD) statistics (Anderson & Darling 1954) for each of the light curves as coming from a Gaussian or a log-normal distribution. The AD allows to test the null hypothesis that a sample is drawn from a population that follows a particular distribution, e.g. in our case a normal distribution. The results obtained are shown in Table 3. The obtained (AD Test) values provide strong evidence to reject the null hypothesis of an underlying normal distribution as the AD Test statistics is greater than the relevant critical value (e.g., Stephens 1974).

Except for PG 1553+113 all sources exhibit a clear deviation from a Gaussian distribution. The flux distributions are instead compatible with being log-normal, which suggests that the underlying variability is of a non-linear, multiplicative origin. Similar results are obtained from the D’Agostino’s K2 and the Shapiro-Wilk test, see Table 3.

3 Searching for periodicity

In order to search for periodicity in blazars we apply two widely used methods, the Lomb-Scargle Periodogram and the Weighted Wavelet Z-transform to light curves.

3.1 Lomb-Scargle Periodogram

The Lomb-Scargle periodogram (LSP) is a commonly used algorithm for detecting and characterizing periodicity in unevenly sampled light curves (Lomb 1976; Scargle 1982). The standard normalized Lomb-Scargle Periodogram (LSP) is equivalent to fitting sine waves of the form $y(t)=A\textbf{ }\cos(\omega t)+B\textbf{ }\sin(\omega t)$, and is defined for a uneven simpled time series ($t_{i},y_{i}$ ) as

[TABLE]

The standard LSP, however, has some limitations. On the one hand, it does not take measurement errors into account; on the other hand, in the calculation the mean of the flux is subtracted, which assumes that the mean of the data and the mean of the fitted sine function are the same. To overcome those limitations we instead use the generalized LSP (GLSP) (Zechmeister & Kürster 2009). As an example the calculated GLSP for PKS 0447-439 using monthly $\gamma$-ray light curves (generated based on the binned likelihood) is shown in Fig. 8. A strong peak at around a period of $945\pm 40$ days is apparent in this figure. We have estimated the uncertainty of the calculated period based on the half-widths at half-maxima (HWHM) of Gaussian fits to the profile at the position of the highest peak. An evaluation of confidence levels determined from Monte Carlo simulations of colored noise is explained in detail in Sec. 4.

For comparison the GLSP method has also been applied to the light curves based on aperture photometry. The possible interference of nearby sources (including background contribution) is evaluated for each source as exemplary shown in Fig. 6 and Fig. 7 for the case of PKS 0447-439. The results agree well with those obtained based on the binned likelihood method, as evident from Table 5.

The observed period Pobs. for a blazar at redshift $z$ is related to its intrinsic period by $P_{\rm int.}=P_{\rm obs.}/(1+z)$. Table  5 shows the GLSP results for all sources, except PKS 0537-441, based on both, monthly (generated based on the binned likelihood) and weekly (aperture photometry) binned light curves. Year-type HE periodicity is found for all five sources, with the results based on aperture photometry and binned likelihood method being in good agreement. For the blazar PKS 0537-441 no significant period has been found when searched for the whole light curve. Furthermore the source is estimated to have a red-noise power spectrum (Covino et al. 2019), which suggests a greater possibility of detecting a fake peak. This source has thus not been included in the following. We confirm, however, that during the initial $\sim 3$-year high state a peak can be seen at periodicity timescale of $\sim 280$ d (see e.g., Sandrinelli et al. 2016). However, for this to be a robust result, a longer observation window showing a larger number of cycles is needed given the timescale of periodicity (Vaughan et al. 2016).

3.2 Weighted Wavelet Z-transform

The GLSP provides an excellent tool for the periodicity analysis of light curves with unevenly-spaced sampling. Nevertheless it does not account for the possibility that in some astrophysical systems quasi-periodic oscillations may develop that vary significantly in frequency and amplitude over a specified period of time. In such cases, the Weighted Wavelet Z-transform (WWZ) method turns out to be a more convenient technique for detecting and quantifying such variations (Foster 1996; Han & van der Baan 2013), and has been applied for the timing analysis of AGN light curves at different wavelengths (e.g. Bhatta et al. 2016; Mohan & Mangalam 2015). The method is based on a similar concept like LSP, where sinusoidal functions are used to fit the data; however, the waves can now be localized in both time and frequency domain to account for the possible transient nature of the QPOs (Torrence & Compo 1998; Bravo et al. 2014).

The WWZ is based on a weighed projection onto three trial functions, $y(t)=\sum_{i}y_{i}\phi_{i}(t)$

[TABLE]

where the ’best-fit’ coefficients $y_{i}$ are the coefficients for which the model function $y(t)$ best fits the data. For each projection the statistic weight is given by

[TABLE]

where $c$ is a constant that determines how rapidly the Morlet wavelet decays, and is usually chosen to be close to $0.0125$. The WWZ power can then be finally defined as

[TABLE]

with Vx and Vy as weighted variations of the data and the model function, respectively and Neff representing the effective number of the data points (for more details, see (Foster 1996).

The WWZ powers for PKS 0447-439 are showed in Fig. 9 as a function of both observing time (horizontal axis) and period (vertical axis). The peaks in the power characterize the strength and duration of a possible quasi-periodic modulation in the data. The WWZ indicates a characteristic period at $(937.3\pm 152.7)$ days (around $2.6$ years) with peak power of about $18.5$ and no significant changes over time. The constancy of the period indicates that the quasi-periodicity is probably driven by a physical process that is stable over duration of observation. The results of the WWZ method are shown in Table 6. The results agree well with the GLSP results.

4 Significance and Uncertainty Estimation

The effects related to irregular sampling of a light curve and the noisy nature of the periodogram can in some situation lead to the generation of false (artificial) periods in the GLSP that could be mistaken as real periodic signal of the source. For this reason, it is important to take such effects into accounts when searching for periodicities, and to analyze their impact on GLSP carefully.

The variability of AGN light curves often exhibits a colored-noise-like behavior with power spectral density (PSD) characterized by a simple power-law of the form PSD($\nu$) $\sim\nu^{-\beta}$ where $\nu$ is the (temporal) frequency and $\beta$ the power-law index.

4.1 Power Spectrum Response Method (PSRESP)

In order to determine the appropriate form of the underlying colored noise needed as an input for simulating the light curves, we first applied the Power Spectrum Response method (PSRESP; Uttley et al. (2002)), which is a widely used technique for the characterization of AGN power spectra (e.g., Chatterjee et al. 2008; Bhatta et al. 2016). The method attempts to fit the binned periodogram with different realisations of a given PSD model in order to estimate the model which maximizes the probability that observed PSD can be reproduced. Our implementation of the PSRESP method is described in full detail in Chatterjee et al. (2008). Selected details are as follow:

We consider a simple power-law model for the underlying power spectrum of the form

[TABLE]

where $A$ is the amplitude of the model at the reference frequency $\nu_{0}$, $\beta$ correspond to the power-spectral slope and $C_{noise}$ is a constant which is fixed at the Poisson noise level for the light curve. We simulate $N=1000$ light curves starting from the underlying model and using the Timmer & Koenig algorithm (see below), and re-sample these with the same sampling interval of the observed light curves. We then calculate the periodogram of the observed light curve (P($\nu$)obs.) and that of each of the simulated light curves (P($\nu$)sim,i where $i=1,N$). The $\chi^{2}_{dist,i}$ statistic is calculated from the underlying model average and observed PSD of each light curve, with

[TABLE]

where $\nu_{min}$ and $\nu_{max}$ are the minimum and maximum frequencies measured by P($\nu$)obs., respectively. We calculate the $\chi^{2}_{dist}$ for each simulated PSD P($\nu$)sim,i over all frequencies between $\nu_{min}$ and $\nu_{max}$ with respect to the data. We then compare $\chi^{2}_{dist}$ with the $\chi^{2}_{obs}$ distribution and count the number $m$ of $\chi^{2}_{dist,i}$ for which $\chi^{2}_{obs}$ is smaller than $\chi^{2}_{dist,i}$. Finally, we calculate the success fraction m/M (goodness of fit), which gives the probabilities of a model being accepted, for a range of $\beta$ from $0.0$ to $2.0$ with step size of $0.05$. The obtained results according to the PSRESP method are summarized in Table  7 and shown in Fig.  10

4.2 MC-simulations of colored-noise light curves

The Timmer and Koenig (TK) algorithm (Timmer & Koenig 1995) is a commonly used method to produce artificial light curves. This technique allows to generate non-deterministic (stochastic Gaussian-distributed) time series from a given underlying PSD model by randomizing both the phase and the amplitude of the Fourier components. Limitations could arise, however, for light curves that exhibit strong deviations from Gaussian distributions (e.g., a burst-like behaviour).

Given the preference for log-normality, we also simulate $5\times 10^{4}$ light curves using the method proposed by (Emmanoulopoulos et al. 2013, E13) to obtain best-fitting results from the PSRESP method. The latter E13 method is able to account for a general PDF, i.e. to match both the PSD and the probability density function (PDF) of an observed light curve, thus relaxing restrictions of the TK method. This does in fact better comply with our previous findings of non-Gaussianity in Sec. 2.3.

The results are shown in Fig. 8 using the (main value of the) best fit slope $\beta$. In many cases the detected periods are close to or above $3\sigma$. Given available data there is rather little difference between the results based on the Timmer & Koenig (1995) and the Emmanoulopoulos et al. (2013) methods. The outcome is, however, obviously dependent on the slope $\beta$. Within the PSRESP inferred range quite different results can be obtained as shown in Table 9, suggesting that a narrowing-down of the PSD slope will be most relevant for assessing the real significance of the detected periods.

4.3 Statistical Confidence of GLSP-detected Period

4.3.1 Local significance

In presence of an priori expectation of periodicity either from theory or from entirely independent observations, one can compute the statistical significance of the detected period at that position. This is the so-called ”local significance” as referenced in several papers (e.g., Bell et al. 2011). Thus, in the local method, the frequency channels corresponding to the GLSP-detected period were searched and the power spectra were recorded in which the peak power at this period exceeds the observed value. Finally, the fraction of occurrences with greater power than detected period is the probability of false positive, resulting from random red noise in the observed light curve.

4.3.2 Global significance

In the majority of the cases of reported periodicities, particularly in gamma-rays, we do not have strong apriori indications of the expected periodicity. In such a circumstance, it is statistically more rigorous to evaluate the ”global” rather than the local significance (e.g. Bell et al. 2011). This constitutes computing the fraction of occurrences of larger peaks at any period within a reasonable range of timescales (dependent on the cadence properties of the light curve), relative to the detected one. This ensures that we factor in false-positives from this larger range of timescales rather than a specific period. This, so called ”look elsewhere effect” (or also Multiple Testing Problem in Statistics), can be quantified in terms of a trial factor (cf. Lyons 2008; Gross & Vitells 2010). This is defined as the ratio between the probability of detecting a peak (or period/excess) at some fixed frequency, to the probability of detecting it anywhere in the (tested) range. The ”look elsewhere effect” has been explored and factored in for detection of resonant peaks in particle physics and indeed other fields including astronomy. We therefore re-evaluated the probability that the power of any observed peak is equal to or greater than a selected value somewhere in the periodogram and calculated in this wise the global significance (cf. Vaughan 2010). The results are shown in Table. 8 and reveal that in the absence of other physical reasons for restricting the period range the QPV evidence is strongly reduced, with none of the sources reaching $3\sigma$ significance. These findings are in line with similar indications in (Covino et al. 2019).

5 Discussion and Conclusions

There are an increasing number of reports suggesting the presence of year-type periodicities in the $\gamma$-ray light curves of Fermi-LAT detected blazars. Given the still limited duration of the light curves (\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10 yr), care needs to be taken, however, to properly assess the significance of these periods against a colored-noise background. Vaughan et al. (2016) have shown, for example, that clear phantom periodicities can be found even in pure noise data, with typical periods corresponding to $(1.5-2.5)$ cycles over the available data.

Our systematic investigation performed in this study reveals that four out of the six investigated blazars show long-term QPO indications (local significance \lower 2.15277pt\hbox{;\buildrel>\over{\sim};}3\sigma) in their Fermi-LAT light curves with an intrinsic period around $1.6$ yr. As we have shown there is some uncertainty as to the significance of these periodicities. Depending on the inferred best-fit PSD slope, the (local) significance can encompass a range from $\sim 2.6\sigma$ to $>5\sigma$ (method 1), cf. Table. 8. The QPO significance is however strongly reduced in terms of a global significance, suggesting that longterm period claims should be treated with caution.

While all sources, except PG 1553+113, show clear indications for log-normality in their distribution of fluxes, incorporation of the appropriate simulation method (Timmer & Koenig 1995 or Emmanoulopoulos et al. 2013) does not have a strong impact on the significance evaluation. Improving the PSD characterization, i.e. by narrowing down the range for the PSD slope will instead be more relevant to better assess the significance.

Though the inferred periods are only tentative, one could speculate about physical mechanisms capable of accounting for year-type periodicity. It seems in fact surprising that the intrinsic periods appear quite similar, clustering around $P\sim 1.6$ yr for sources at different redshifts.

Perhaps one of the most natural scenarios for the origin of longterm periodic variability is related to the orbital motion in a supermassive binary black hole system (SBBHS). In principle a SBBHS phase is likely to occur in radio-loud AGN (being hosted by elliptically galaxies) at some stage. However, the periods inferred here appear rather short to plausibly relate them to orbital SBBH motion. Gravitation radiation would lead to coalescence on a characteristic timescale $t_{grav}\simeq 7\times 10^{3}\,(1+q)^{1/3}/[q\,(M_{BH}/5\times 10^{8}M_{\odot})^{5/3}]\,(P/2\leavevmode\nobreak\ \mathrm{yr})^{8/3}$ yr only, so that for typical systems with $q:=m/M\geq 0.05$ the presumed source state would be highly unlikely. In fact, one rather expects orbital periods of the order $\sim 10$ yr for SBBHSs in blazars (Rieger 2007). This would then also be compatible with pulsar timing constraints on the inferred gravitational background (Holgado et al. 2018). As orbital motion usually allows for the shortest periodic driving (e.g., Rieger 2004), a direct SBBH origin of the observed periods appears less likely. This does not argue against the presence of SBBHSs in blazars in general, but simply cautions to directly relate periodicities of the order of $P\sim 1$ yr to such systems.

Alternatively, QPOs might be related to quasi-periodic changes in accretion flow conditions that are effectively transmitted to the jet, modulating its non-thermal emission properties. Time-dependent modulations of the transition radius $r_{\rm t}$ between an outer cooling-dominated (standard) disc and an inner radiatively inefficient flow (ADAF), for example, could lead to periodic mass flux variations (e.g., Gracia et al. 2003). If one requires the advective timescale $t_{\rm ad}(r_{\rm t})\sim r_{g}\,(0.5\alpha\,c)^{-1}(r_{\rm t}/r_{g})^{3/2}$, with $\alpha=0.25$ the viscosity coefficient and $r_{g}=GM_{\rm BH}/c^{2}$ the gravitational radius, to be (at most) comparable to $P$, this would place the transition radius at a characteristic scale of $r_{\rm t}\lesssim 200\leavevmode\nobreak\ r_{g}(P/1.6\leavevmode\nobreak\ \mathrm{yr})^{2/3}\,\ (5\times 10^{8}M_{\odot}/M_{\rm BH})^{2/3}$ for a reference mass of $M_{\rm BH}=5\times 10^{8}M_{\odot}$. This seems compatible with estimates for the transition radius in BL Lacs (e.g., Cao 2003; Xie et al. 2008). This would then suggest a similar black hole mass-scale for the systems investigated here.

On the other hand, year-type QPO could perhaps also trace plasma motion in the jet close to its outer jet radius $r_{0}$. In the lighthouse model (Camenzind & Krockenberger 1992), for example, the disk-related jet is initially rotating, leading to a helical trajectory for a component injected on scales $r_{0}\sim 10\,r_{L}$ beyond the light cylinder $r_{L}\sim 10\,r_{g}$ of the innermost part of the disk magnetosphere. Angular momentum conservation would imply a characteristic intrinsic period for such a component of $P=2\pi r_{L}(r_{0}/r_{L})^{2}/c\lesssim 2\;(r_{0}/20r_{L})^{2}(M_{\rm BH}/5\times 10^{8}M_{\odot})$ yr. While collimation might occur earlier (e.g., Fendt 1997), thereby reducing the jet radius, sligthly changing footpoint radii could possibly compensate for this. The lighthouse model was originally designed to account for observed QPOs with $P_{obs}\leq$ few weeks by the taking travel time effects with respect to an outwardly moving, (single) flaring component into account. It seems likely however, that the fundamental period $P$ might be visible even if the flux would be suppressed quickly. It will be interesting to probe this with a larger source sample, as intrinsic periods in this case are expected to be less than $2$ yr for a typcial black hole mass range.

The fact that the intrinsic periods seems to be clustering around $P\sim 1.6$ yr remains particularly interesting and suggestive of a common physical mechanism. The inferred periods are, however, only tentative. Adding one or two more cycle of data (i.e., $\sim 2-3$ yr) is expected to significantly improve the situation and to help to clarify their putative presence and physical implication.

Acknowledgements.

We would like to thank Stefan Wagner from Landessternwarte Heidelberg (LSW) for comments and suggestions. FMR acknowledges support by a DFG Heisenberg Fellowship RI 1187/6-1.