Issue 
A&A
Volume 634, February 2020



Article Number  A120  
Number of page(s)  11  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/201935117  
Published online  20 February 2020 
Evaluating quasiperiodic variations in the γray light curves of FermiLAT blazars
^{1}
MaxPlanckInstitut für Kernphysik, Saupfercheckweg 1, 69117 Heidelberg, Germany
email: fbenkhal@mpihd.mpg.de
^{2}
ZAH, Institut für Theoretische Astrophysik, Universität Heidelberg, Philosophenweg 12, 69120 Heidelberg, Germany
^{3}
Data Assimilation Research Centre, University of Reading, Whiteknights Rd, Reading RG6 6BB, UK
Received:
23
January
2019
Accepted:
13
March
2019
Context. The detection of periodicities in the light curves of active galactic nuclei (AGNs) could have profound consequences for our understanding of the nature and radiation physics of these objects. At high energies (HE; E > 100 MeV), five blazars (PG 1553+113, PKS 2155−304, PKS 0426−380, PKS 0537−441 and PKS 0301−243) have been reported to show yearlike quasiperiodic variations (QPVs) with significance > 3σ. As these findings are based on only a few cycles, care needs to be taken to properly account for random variations that can produce intervals of seemingly periodic behavior.
Aims. We present results of an updated timing analysis for six blazars (adding PKS 0447−439 to the above), using suitable methods to evaluate their longterm variability properties and to search for QPVs in their light curves.
Methods. We generate γray light curves covering almost ten years, study their timing properties, and search for QPVs using the LombScargle Periodogram and the Wavelet Ztransform. Extended Monte Carlo simulations are used to evaluate the statistical significance.
Results. (1) Comparing their probability density functions, all sources (except PG 1553+113) exhibit a clear deviation from a Gaussian distribution, but are consistent with being lognormal, suggesting that the underlying variability is of a nonlinear, multiplicative nature. (2) Apart from PKS 0301−243, the power spectral density for all investigated blazars is close to flicker noise (powerlaw slope −1). (3) Possible QPVs with a local significance ≳3σ are found in all light curves (apart from PKS 0426−380 and PKS 0537−441), with observed periods in the range (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 yearlike HE QPVs in blazars. Somewhat surprisingly, the putative redhiftcorrected period all cluster around ∼1.6 yr. We speculate on possible implications for QPV generation.
Key words: γ rays: galaxies / BL Lacertae objects: general / galaxies: jets / radiation mechanisms: nonthermal / galaxies: active
© F. Ait Benkhali et al. 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Open Access funding provided by Max Planck Society.
1. Introduction
Blazars belong to the most luminous and variable extragalactic sources in the Universe. They represent a special subclass of active galactic nuclei (AGNs) characterized by a relativistic plasma jet oriented very close to the line of sight at angles θ ≤ 1/Γ, where Γ is the bulk Lorentz factor (e.g., Urry & Padovani 1995; Ackermann et al. 2015a). Blazar sources are known to be variable across all wavelengths from radio frequencies to teraelectronvolt γray energies, and on a wide range of timescales from subminute to several years (e.g., Böttcher & Chiang 2002; Aharonian et al. 2007). Their multiwavelength spectral energy distributions (SED) often exhibit two bumps. The first, lowenergy bump (peaking at infrared to Xray frequencies) is usually interpreted as synchrotron emission from highly relativistic electrons within the jet. The second bump on the other hand, peaking in the HE range, has frequently been related to inverse Compton emission (IC) upscattering 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 bands (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 ∼12 and 29.6 years being reported, respectively (Valtonen et al. 2006; Li et al. 2009). In the highenergy gammaray 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 to exhibit longterm quasiperiodic variability (these blazars are therefore referred to as QPVs) in their γray fluxes with significance higher than 3σ (Ackermann et al. 2015b; Sandrinelli et al. 2014, 2016; Zhang et al. 2017a,b; 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 behavior at different wavelengths and on various timescales to the motion of two supermassive black holes in a binary system, general helical jet structures, shocks, or instabilities of the disk or jetplasma flow (e.g., Rieger 2004; Wiita 2011; Sobacchi et al. 2017; Caproni et al. 2017; Tavani et al. 2018; Yan et al. 2018).
In the present paper we use FermiLAT γray data between 2008 August 4 and 2018 March 15 to reevaluate the five blazars mentioned above and to investigate an additional source, PKS 0449−439. Results of a periodicity search based on the generalized LombScargle periodogram (GLSP) and the Weighted Wavelet Ztransform (WWZ) are presented. The signifiance of the inferred periods is then estimated on the basis of Bootstrap resampling and lightcurve simulations using the Timmer & Koenig (1995) and Emmanoulopoulos algorithms (Timmer & Koenig 1995; Emmanoulopoulos et al. 2013)
The paper is organized as follows. Section 2 provides an exemplary illustration of the performed FermiLAT data analysis procedures to generate γray light curves based on a binned likelihood method and also presents the aperture photometry technique. In Sect. 3 we describe the GLSP and WWZ techniques used here to search for periodicity. Significance estimation and variability characterization is presented in Sect. 4. Finally, Sect. 5 provides a summary and discussion of the results.
2. Observations and analysis
2.1. FermiLAT likelihood analysis
FermiLAT on board the Fermi satellite is an electronpositron pairconversion γray detector sensitive to photon in the energy range between ∼20 MeV and 500 GeV. The LAT has a wide field of view of ∼2.4 sr and observes the entire sky every two orbits. It has a point spread function (PSF) < 0.8° and the largest effective area of approximately 8000 cm^{2} above 1 GeV. FermiLAT data were downloaded from the FermiLAT Data Server^{1} and events were selected between 2008 August 4 and 2018 March 15 (MET from 239557418.0 to 542821500.0), covering about ∼9.5 years with energy range from 100 MeV to 500 GeV in a circular region of interest (ROI) of 15° 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° 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 (IRFs). 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 binnedlikelihood analysis with a bin size of 0.1° was performed. The spectral indices were allowed to vary for sources located within a radius of 5° around the position of the investigated source. The Test Statistic value (TS) defined as TS = −2ln(L_{0}/L_{1}) was used to determine the source detection significance, with the threshold set to TS = 25 (∼5σ). The significance of a source detection is given by (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 ≃ 29305 corresponding to 171σ. Both energy and temporal bins with TS < 9 (or ∼3σ) are set as upper limits throughout the paper.
The appearance of new sources within the ROI could in principle strongly influence the γray spectrum and light curve of the investigated source. In order to investigate this we also searched for possible new γray sources within a FOV of 15°. Our analysis of 9.5 years of FermiLAT 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 (Power law) and maximize the likelihood as a function of its flux.
We produced the γray SED and light curves of each source through the binned maximum likelihoodfitting technique, 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; to do so, 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 ten energy bins. We performed a maximum likelihood analyses, exploring different spectral forms for the whole energy range, namely single powerlaw (PL), broken powerlaw (BPL), logparabola (LP), and powerlaw super exponential cutoff (PLSEC). In the case of PKS 0447−439, a likelihood ratio test comparison yields TS = −2log(L_{PL}/L_{BPL})≃53.8, thus preferring BPL over PL at a significance level of 7.3σ. The best fit indices are Γ_{1} = 1.84 ± 0.01 and Γ_{2} = 2.82 ± 0.17 with a break at an energy E_{b} = (42.8 ± 1.8) GeV. The obtained results of the maximum likelihood analysis are summarized for the different models in the Table 1.
Fig. 1. SED of PKS 0447−439 in the energy band of 100 MeV to 500 GeV as extracted from the complete data set (2008−2018) along with a Broken Power Law fit (solid red line). The regions show the 1σ confidence intervals resulting from the fit with BPL. The resulting break energy is E_{b} = 42.84 ± 1.78. 
Parameters obtained for the fit of the FermiLAT energy spectrum of PKS 0447−439 using power law (PL), log parabola (LP), broken power law (BPL) and PLSuperExpCutOff (PLSEC) spectral models.
2.2. FermiLAT γray light curves
We generate γray light curves for all our six sources using the individual bestfit 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 γ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 γevents are selected in a circular ROI of 15° 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 γray flux of ⟨ϕ⟩=(7.62 ± 0.13)×10^{−8} cm^{2} s^{−1}. The source exhibits flux variability during the whole observational period, with flux levels of peaktopeak oscillations changing by a factor four. The calculated fractional variability index is F_{var} = (47.91 ± 1.16)% (Vaughan et al. 2003; Edelson et al. 2002).
Fig. 2. γray light curve of PKS 0447−439 from August 2008 to December 2017 in the energy range between 100 MeV and 500 GeV in intervals of 30 days. The dotted line gives a fit with a constant flux. There is clear evidence for variability. Vertical error bars indicate the onesigma error bars. The gray shaded regions show the 1σ, 3σ, and 5σ confidence intervals resulting from the fit with a constant function. 
Spectral parameters obtained for the best fit of FermiLAT data for each source from the energy range between 100 MeV and 500 GeV using four different spectral models: PL, BPL, LP and PLSEC.
The aperture photometry method on the other hand provides a modelindependent measure of the γ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° to exclude most background γevents and to focus on the events that are most likely associated with the source itself. We produce γ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 γray light curve for PKS 0447−439 obtained with this method.
Fig. 3. γray light curve of PKS 0447−439 from 2008 to 2018 in the energy range (0.1 − 500) GeV obtained using the aperture photometry technique for a radius of 1° around the source. A binning of one week is employed. Error bars correspond to 1σ error bars. The gray shaded regions show the 1σ, 3σ, and 5σ confidence intervals resulting from the fit with a constant function. 
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 FermiLAT analysis. 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 ρ ∼ 0.936.
Fig. 4. Aperture photometry flux in comparison with the binned likelihood flux calculated for PKS 0447−439 (logarithmic scale). The red solid line represents a fit of the data with a linear function and the gray shaded regions show the resulting 3σ and 5σ confidence intervals. Error bars indicate the 1σ error interval. No background subtraction is used in the aperture photometry method. 
2.3. Gaussian versus lognormal 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 longterm γray emission reveals any preference for a Gaussian (normal) or a lognormal flux distribution. It is expected that such a feature offers important clues for understanding the central engine of a blazer (e.g., Uttley et al. 2005; Shah et al. 2018; Rieger 2019). A lognormal distribution of γray fluxes, for example, would suggest that the mechanism driving the variability is multiplicative in nature rather than additive. Such multiplicative behavior could possibly be related to accretion disk fluctuations (Lyubarskii 1997; Arévalo & Uttley 2006) or the particle acceleration process itself (Sinha et al. 2018). As a possible caveat we note that the 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 fluxhistograms. For each source the weekly γray fluxes were distributed in a histogram of fluxes. We fit all histograms in logscale, with Gaussian G(ϕ) and lognormal L(ϕ) distribution functions given by
and
respectively, where σ and μ 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 Fig. 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 lognormal distribution. The AD allows one to determine whether one can accept or reject the null hypothesis that a sample is drawn from a population that follows a particular distribution; in our case a normal distribution. The results obtained are shown in Table 3. The obtained (AD test) values provide strong evidence in favor of rejecting the null hypothesis of an underlying normal distribution as the AD test value is greater than the relevant critical value (e.g., Stephens 1974).
Fig. 5. Normalized histograms of γray photon fluxes fitted with a lognormal (solid red line) and Gaussian (blue solid line), respectively. Apart from PG 1553+113, all sources show a clear preference for a lognormal distribution. 
Best fit parameters for the lognormal and normal (Gaussian) distribution, respectively; the AD test statistics and are the Anderson Darling test and the reduced χ^{2}, respectively, [*] in units of photon cm^{−2} s^{−1}.
Except for PG 1553+113, all sources exhibit a clear deviation from a Gaussian distribution. The flux distributions are instead compatible with being lognormal, which suggests that the underlying variability is of a nonlinear, multiplicative origin. Similar results are obtained from the D’Agostino’s K^{2} and the Shapiro–Wilk test; see Table 4.
K^{2} Test of D’Agostino tests the null hypothesis that a sample comes from a normal distribution.
3. Searching for periodicity
In order to search for periodicity in blazars we apply two widely used methods to light curves: the LombScargle periodogram (LSP) and the WWZ.
3.1. LombScargle periodogram
The LSP is a commonly used algorithm for detecting and characterizing periodicity in unevenly sampled light curves (Lomb 1976; Scargle 1982). The standard normalized LSP is equivalent to fitting sine waves of the form y(t) = A cos(ωt)+B sin(ωt), and is defined for an uneven, simple time series (t_{i}, y_{i} ) as
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 GLSP (Zechmeister & Kürster 2009). As an example, the calculated GLSP for PKS 0447−439 using monthly γray light curves (generated based on the binned likelihood) is shown in Fig. 8. A strong peak at around a period of 945 ± 40 days is apparent in this figure. We have estimated the uncertainty of the calculated period based on the halfwidths at halfmaxima (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 Sect. 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 shown in Figs. 6 and 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.
Fig. 6. Smoothed count map (logarithmic scale) of the PKS 0447−439 region in energy range from 100 MeV up to 500 GeV as seen by FermiLAT based on data between August 2008 and December 2017. The color bar has units of counts per pixel and the pixel dimensions are 0.1 × 0.1°. The locations of the background region and surrounding sources (labeled P1 = new source beyond the 3FGL catalog, P2 = 3FGL J0438.8−4519, P3 = 3FGL J0437.2−4713 and P4 = 3FGL J0455.7−4617) are shown as circles. The magenta circle shows the location of PKS 0447−439 and the yellow circles represent the background region. 
Fig. 7. Generalized LSPs obtained based on light curves generated using aperture photometry for PKS 0447−439 and surrounding sources (P1,P2,P3,P4) shown as circles in Fig. 6. 
Periodicity results derived from the GLSP method for the monthly and weekly γray light curves, respectively, and including the intrinsic (redshiftcorrected) period. [*] in units of days.
The observed period P_{obs.} for a blazar at redshift z is related to its intrinsic period by P_{int.} = P_{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. Yeartype HE periodicity is found for all five sources, with the results based on aperture photometry or binned likelihood method being in good agreement. For the blazar PKS 0537−441 no significant period was found from a search of the entire light curve. Furthermore, the source is estimated to have a rednoise power spectrum (Covino et al. 2019), which suggests a greater possibility of detecting a fake peak. This source is therefore not included in the following. We confirm however that during the initial approximately threeyear high state, a peak can be seen on a periodicity timescale of ∼280 d (see e.g., Sandrinelli et al. 2016). However, a longer observation window showing a larger number of cycles is needed, given the timescale of periodicity (Vaughan et al. 2016), for this result to be confirmed.
3.2. Weighted wavelet Ztransform
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 quasiperiodic oscillations may develop that vary significantly in frequency and amplitude over a specified period of time. In such cases, the 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 in 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 quasiperiodic oscillations (QPOs; Torrence & Compo 1998; Bravo et al. 2014).
The WWZ is based on a weighed projection onto three trial functions, y(t) = ∑_{i}y_{i}ϕ_{i}(t),
where the “bestfit” 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
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 finally be defined as
with V_{x} and V_{y} as weighted variations of the data and the model function, respectively, and N_{eff} representing the effective number of data points (for more details, see Foster 1996).
The WWZ powers for PKS 0447−439 are shown 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 quasiperiodic modulation in the data. The WWZ indicates a characteristic period at (937.3 ± 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 quasiperiodicity is probably driven by a physical process that is stable over the duration of the observation. The results of the WWZ method are shown in Table 6 and the results agree well with those from the GLSP.
Fig. 8. GLPS obtained using the monthly light curves of PKS 0447−439 is shown in cyan (circles); the red solid line, green line, and orange line represent the 3σ, 2.5σ, and 2σ confidence levels, respectively, calculated based on simulations of 50 000 light curves using the method of Timmer & Koenig (1995). 
Fig. 9. Twodimensional WWZ of the γray light curve of PKS 0447−439 based on FermiLAT data from August 2008 to December 2017. The color scale represents the Zstatistics of the WWZ power of a certain period at a given time. The panel shows the signature of a quasiperiodic oscillation at a (observed) period of ∼2.6 years without any significant changes over time. 
Summary results of the strongest periods for the WWZ method based on the monthly γray light curves, including observed and intrinsic period.
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 situations lead to the generation of false (artificial) periods in the GLSP that could be mistaken as a real periodic signal of the source. For this reason, it is important to take such effects into account when searching for periodicities and to analyze their impact on the GLSP carefully.
The variability of AGN light curves often exhibits a colorednoiselike behavior with PSD characterized by a simple power law of the form PSD(ν) ∼ ν^{−β} where ν is the (temporal) frequency and β the powerlaw index.
4.1. Power spectrum response method
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 (PSRESP; Uttley et al. 2002) method, 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 determine which model maximizes the probability that the 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 follows.
We consider a simple powerlaw model for the underlying power spectrum of the form
where A is the amplitude of the model at the reference frequency ν_{0}, β corresponds to the powerspectral 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 resample these with the same sampling interval of the observed light curves. We then calculate the periodogram of the observed light curve (P(ν)_{obs.}) and that of each of the simulated light curves (P(ν)_{sim, i} where i = 1, …N). The statistic is calculated from the underlying model average and observed PSD of each light curve, with
where ν_{min} and ν_{max} are the minimum and maximum frequencies measured by P(ν)_{obs.}, respectively. We calculate the for each simulated PSD P(ν)_{sim, i} over all frequencies between ν_{min} and ν_{max} with respect to the data. We then compare with the distribution and count the number m of for which is smaller than . Finally, we calculate the success fraction m/M (goodness of fit), which gives the probabilities of a model being accepted, for a range of β from 0.0 to 2.0 with a step size of 0.05. The obtained results according to the PSRESP method are summarized in Table 7 and shown in Fig. 10.
Fig. 10. Success fraction vs. slope (β) for all three PSDs (FermiLAT). The success fractions indicate the goodness of fit obtained from the PSRESP method (see text). (a) Weekly, (b) monthly light curves. 
Slopes β obtained from fitting the γray power spectra with a simple powerlaw model, and the corresponding success fraction calculated based on the PSRESP method (Uttley et al. 2002) (see Sect. 4.1) (the success fractions indicate the goodness of fit obtained from the PSRESP method).
4.2. MCsimulations of colorednoise light curves
The Timmer and Koenig (TK) algorithm (Timmer & Koenig 1995) is a commonly used method to produce artificial light curves. This technique allows for the generation of nondeterministic (stochastic Gaussiandistributed) 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 burstlike behavior).
Given the preference for lognormality, we also simulate 5 × 10^{4} light curves using the method proposed by Emmanoulopoulos et al. (2013, E13) to obtain bestfitting results from the PSRESP method. The latter E13 method is able to account for a general PDF, that is, it can match both the PSD and the PDF of an observed light curve, thus relaxing restrictions of the TK method. This does in fact better comply with our previous findings of nonGaussianity in Sect. 2.3.
The results are shown in Fig. 8 using the (main value of the) best fit slope β. In many cases the detected periods are close to or above 3σ. Given available data there is rather little difference between the results based on the Timmer & Koenig (1995) and Emmanoulopoulos et al. (2013) methods. The outcome is however obviously dependent on the slope β. Within the PSRESP inferred range quite different results can be obtained as shown in Table 9, suggesting that a narrowingdown of the PSD slope will be most relevant for assessing the real significance of the detected periods.
Significance values of the observed periods for the local and the global method using the Timmer & Koenig (1995) and Emmanoulopoulos et al. (2013) algorithm, based on the monthly and weekly γray light curves.
Significance values of the observed periods using either the local or the global method with respect to the monthly and weekly γray light curves of PG 1553+113 and PKS 0447−489, respectively.
4.3. Statistical confidence of the GLSPdetected period
4.3.1. Local significance
In the presence of an a 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 socalled “local significance” as referenced in several papers (e.g., Bell et al. 2011). Therefore, using the local method, the frequency channels corresponding to the GLSPdetected period were sought and the power spectra in which the peak power at this period exceeded the observed value were recorded. Finally, the fraction of occurrences with power greater than the detected period is the probability of a false positive signal 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 gammarays, we do not have strong a priori 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 socalled “lookelsewhere” effect (or also “multiple testing problem” in statistics) can be quantified in terms of a trial factor (e.g., Lyons 2008; Gross & Vitells 2010). This is defined as the ratio of 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 lookelsewhere effect has been explored and factored in for detection of resonant peaks in particle physics and indeed other fields including astronomy. We therefore reevaluated the probability that the power of any observed peak is equal to or greater than a selected value somewhere in the periodogram and calculated the global significance (see e.g., 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σ 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 periodicities of years in duration in the γray light curves of FermiLATdetected blazars. Given the still limited duration of the light curves (≲10 yr), care needs to be taken however to properly assess the significance of these periods against a colorednoise 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 longterm QPO indications (local significance ≳3σ) in their FermiLAT 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 bestfit PSD slope, the (local) significance can encompass a range from ∼2.6σ to > 5σ (method 1); see Table 8. The QPO significance is however strongly reduced in terms of global significance, suggesting that claims of longterm periods should be treated with caution.
While all sources except PG 1553+113 show clear indications for lognormality 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, that is, 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 over the possible physical mechanisms capable of causing these longer periodicities. It is interesting, for example, that the intrinsic periods appear quite similar, clustering around P ≈ 1.6 yr for sources at different redshifts (with the caveat of some possible redshift uncertainty for PKS 0447−439).
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 radioloud AGNs (hosted by elliptical galaxies) at some stage. However, the periods inferred here appear too short to plausibly relate them to orbital SBBH motion. Gravitation radiation would lead to coalescence on a characteristic timescale of only t_{grav} ≃ 7 × 10^{3} (1 + q)^{1/3}/[q (M_{BH}/5 × 10^{8} M_{⊙})^{5/3}] (P/2 yr)^{8/3} yr, meaning that for typical systems with q := m/M ≥ 0.05 the presumed source state would be highly unlikely. In fact, one alternatively expects orbital periods of the order of ∼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 suggests that caution should be taken when directly relating periodicities of the order of P ∼ 1 yr to such systems.
Alternatively, QPOs might be related to quasiperiodic changes in accretion flow conditions that are effectively transmitted to the jet, modulating its nonthermal emission properties. Timedependent modulations of the transition radius r_{t} between an outer coolingdominated (standard) disk 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_{ad}(r_{t})∼r_{g} (0.5α c)^{−1}(r_{t}/r_{g})^{3/2} to be (at most) comparable to P, with α = 0.25 being the viscosity coefficient and r_{g} = GM_{BH}/c^{2} the gravitational radius, this would place the transition radius at a characteristic scale of r_{t} ≲ 200 r_{g}(P/1.6 yr)^{2/3} (5 × 10^{8} M_{⊙}/M_{BH})^{2/3} for a reference mass of M_{BH} = 5 × 10^{8} M_{⊙}. 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 massscale for the systems investigated here.
On the other hand, yeartype QPVs 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 diskrelated jet is initially rotating, leading to a helical trajectory for a component injected on scales r_{0} ∼ 10 r_{L} beyond the light cylinder r_{L} ∼ 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πr_{L}(r_{0}/r_{L})^{2}/c ≲ 2 (r_{0}/20r_{L})^{2}(M_{BH}/5 × 10^{8} M_{⊙}) yr. While collimation might occur earlier (e.g., Fendt 1997), thereby reducing the jet radius, slightly changing footpoint radii could possibly compensate for this. The lighthouse model was originally designed to account for observed QPOs with P_{obs} of less than or equal to a few weeks by taking into account travel time effects with respect to an outwardly moving, (single) flaring component. It seems likely however that the fundamental period P might be visible even if the flux were quickly suppressed. 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 typical black hole mass range.
The fact that the intrinsic periods seem to be clustering around P ∼ 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 cycles of data (i.e., ∼2 − 3 yr) is expected to significantly improve the situation and to help to clarify their putative presence and physical implication.
Acknowledgments
Useful comments by an anonymous referee are gratefully acknowledged. We would also like to thank J. Graham, A. Sinha, and S. Wagner (LSW Heidelberg) for helpful comments and suggestions. FMR acknowledges support by a DFG Heisenberg Fellowship RI 1187/61.
References
 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJS, 187, 460 [NASA ADS] [CrossRef] [Google Scholar]
 Acero, F., Ackermann, M., Ajello, M., et al. 2015, ApJS, 218, 23 [Google Scholar]
 Ackermann, M., Ajello, M., Atwood, W. B., et al. 2015a, ApJ, 810, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Ackermann, M., Ajello, M., Albert, A., et al. 2015b, ApJ, 813, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Aharonian, F., Akhperjanian, A. G., BazerBachi, A. R., et al. 2007, ApJ, 664, L71 [Google Scholar]
 Anderson, T. W., & Darling, D. A. 1954, J. Am. Stat. Assoc., 49, 765 [Google Scholar]
 Arévalo, P., & Uttley, P. 2006, MNRAS, 367, 801 [NASA ADS] [CrossRef] [Google Scholar]
 Atwood, W., Albert, A., Baldini, L., et al. 2013, 2012 Fermi Symposium Proceedings – eConf C121028 [Google Scholar]
 Bell, M. E., Tzioumis, T., Uttley, P., et al. 2011, MNRAS, 411, 402 [NASA ADS] [CrossRef] [Google Scholar]
 Bhatta, G., Zola, S., Stawarz, Ł., et al. 2016, ApJ, 832, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Böttcher, M., & Chiang, J. 2002, ApJ, 581, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Bravo, J. P., Roque, S., Estrela, R., Leão, I. C., & De Medeiros, J. R. 2014, A&A, 568, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Camenzind, M., & Krockenberger, M. 1992, A&A, 255, 59 [NASA ADS] [Google Scholar]
 Cao, X. 2003, ApJ, 599, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Caproni, A., Abraham, Z., Motter, J. C., & Monteiro, H. 2017, ApJ, 851, L39 [NASA ADS] [CrossRef] [Google Scholar]
 Chatterjee, R., Jorstad, S. G., Marscher, A. P., et al. 2008, ApJ, 689, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Covino, S., Sandrinelli, A., & Treves, A. 2019, MNRAS, 482, 1270 [NASA ADS] [CrossRef] [Google Scholar]
 D’Agostino, R. B. 1970, Biometrika, 57, 679 [Google Scholar]
 D’Agostino, R., & Pearson, E. S. 1973, Biometrika, 60, 613 [Google Scholar]
 Edelson, R., Turner, T. J., Pounds, K., et al. 2002, ApJ, 568, 610 [NASA ADS] [CrossRef] [Google Scholar]
 Emmanoulopoulos, D., McHardy, I. M., & Papadakis, I. E. 2013, MNRAS, 433, 907 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, J. H., Lin, R. G., Xie, G. Z., et al. 2002, A&A, 381, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fendt, C. 1997, A&A, 319, 1025 [NASA ADS] [Google Scholar]
 Foster, G. 1996, AJ, 112, 1709 [NASA ADS] [CrossRef] [Google Scholar]
 Gracia, J., Peitz, J., Keller, C., & Camenzind, M. 2003, MNRAS, 344, 468 [NASA ADS] [CrossRef] [Google Scholar]
 Gross, E., & Vitells, O. 2010, Eur. Phys. J. C, 70, 525 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Han, J., & van der Baan, M. 2013, Geophysics, 78, O9 [CrossRef] [Google Scholar]
 Holgado, A. M., Sesana, A., Sandrinelli, A., et al. 2018, MNRAS, 481, L74 [NASA ADS] [CrossRef] [Google Scholar]
 Kadler, M., Hughes, P. A., Ros, E., Aller, M. F., & Aller, H. D. 2006, A&A, 456, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 King, O. G., Hovatta, T., MaxMoerbeck, W., et al. 2013, MNRAS, 436, L114 [NASA ADS] [CrossRef] [Google Scholar]
 Li, H. Z., Xie, G. Z., Chen, L. E., et al. 2009, PASP, 121, 1172 [NASA ADS] [CrossRef] [Google Scholar]
 Lomb, N. R. 1976, Ap&SS, 39, 447 [NASA ADS] [CrossRef] [Google Scholar]
 Lyons, L. 2008, Ann. Appl. Stat., 2, 887 [CrossRef] [MathSciNet] [Google Scholar]
 Lyubarskii, Y. E. 1997, MNRAS, 292, 679 [NASA ADS] [CrossRef] [Google Scholar]
 Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396 [NASA ADS] [CrossRef] [Google Scholar]
 Mohan, P., & Mangalam, A. 2015, ApJ, 805, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Muriel, H., Donzelli, C., Rovero, A. C., & Pichel, A. 2015, A&A, 574, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Prokhorov, D. A., & Moraghan, A. 2017, MNRAS, 471, 3036 [NASA ADS] [CrossRef] [Google Scholar]
 Rieger, F. M. 2004, ApJ, 615, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Rieger, F. M. 2007, Ap&SS, 309, 271 [NASA ADS] [CrossRef] [Google Scholar]
 Rieger, F. 2019, Galaxies, 7, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Sandrinelli, A., Covino, S., & Treves, A. 2014, ApJ, 793, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Sandrinelli, A., Covino, S., & Treves, A. 2016, ApJ, 820, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Scargle, J. D. 1982, ApJ, 263, 835 [NASA ADS] [CrossRef] [Google Scholar]
 Shah, Z., Mankuzhiyil, N., Sinha, A., et al. 2018, Res. Astron. Astrophys., 18, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Shapiro, S. S., & Wilk, M. B. 1965, Biometrika, 52, 591 [Google Scholar]
 Sinha, A., Khatoon, R., Misra, R., et al. 2018, MNRAS, 480, L116 [NASA ADS] [CrossRef] [Google Scholar]
 Sobacchi, E., Sormani, M. C., & Stamerra, A. 2017, MNRAS, 465, 161 [NASA ADS] [CrossRef] [Google Scholar]
 Stephens, M. A. 1974, J. Am. Stat. Assoc., 69, 730 [Google Scholar]
 Tavani, M., Cavaliere, A., MunarAdrover, P., & Argan, A. 2018, ApJ, 854, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Timmer, J., & Koenig, M. 1995, A&A, 300, 707 [NASA ADS] [Google Scholar]
 Torrence, C., & Compo, G. P. 1998, Bull. Am. Meteorol. Soc., 79, 61 [Google Scholar]
 Urry, C. M., & Padovani, P. 1995, PASP, 107, 803 [NASA ADS] [CrossRef] [Google Scholar]
 Uttley, P., McHardy, I. M., & Papadakis, I. E. 2002, MNRAS, 332, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Uttley, P., McHardy, I. M., & Vaughan, S. 2005, MNRAS, 359, 345 [NASA ADS] [CrossRef] [Google Scholar]
 Valtonen, M. J., Lehto, H. J., Sillanpää, A., et al. 2006, ApJ, 646, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Vaughan, S. 2010, MNRAS, 402, 307 [NASA ADS] [CrossRef] [Google Scholar]
 Vaughan, S., Edelson, R., Warwick, R. S., & Uttley, P. 2003, MNRAS, 345, 1271 [NASA ADS] [CrossRef] [Google Scholar]
 Vaughan, S., Uttley, P., Markowitz, A. G., et al. 2016, MNRAS, 461, 3145 [NASA ADS] [CrossRef] [Google Scholar]
 Wiita, P. J. 2011, J. Astrophys. Astron., 32, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Xie, Z. H., Hao, J. M., Du, L. M., Zhang, X., & Jia, Z. L. 2008, PASP, 120, 477 [NASA ADS] [CrossRef] [Google Scholar]
 Yan, D., Zhou, J., Zhang, P., Zhu, Q., & Wang, J. 2018, ApJ, 867, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zhang, P.F., Yan, D.H., Zhou, J.N., et al. 2017a, ApJ, 845, 82 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, P., Yan, D., Liao, N., et al. 2017b, ApJ, 842, 10 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Parameters obtained for the fit of the FermiLAT energy spectrum of PKS 0447−439 using power law (PL), log parabola (LP), broken power law (BPL) and PLSuperExpCutOff (PLSEC) spectral models.
Spectral parameters obtained for the best fit of FermiLAT data for each source from the energy range between 100 MeV and 500 GeV using four different spectral models: PL, BPL, LP and PLSEC.
Best fit parameters for the lognormal and normal (Gaussian) distribution, respectively; the AD test statistics and are the Anderson Darling test and the reduced χ^{2}, respectively, [*] in units of photon cm^{−2} s^{−1}.
K^{2} Test of D’Agostino tests the null hypothesis that a sample comes from a normal distribution.
Periodicity results derived from the GLSP method for the monthly and weekly γray light curves, respectively, and including the intrinsic (redshiftcorrected) period. [*] in units of days.
Summary results of the strongest periods for the WWZ method based on the monthly γray light curves, including observed and intrinsic period.
Slopes β obtained from fitting the γray power spectra with a simple powerlaw model, and the corresponding success fraction calculated based on the PSRESP method (Uttley et al. 2002) (see Sect. 4.1) (the success fractions indicate the goodness of fit obtained from the PSRESP method).
Significance values of the observed periods for the local and the global method using the Timmer & Koenig (1995) and Emmanoulopoulos et al. (2013) algorithm, based on the monthly and weekly γray light curves.
Significance values of the observed periods using either the local or the global method with respect to the monthly and weekly γray light curves of PG 1553+113 and PKS 0447−489, respectively.
All Figures
Fig. 1. SED of PKS 0447−439 in the energy band of 100 MeV to 500 GeV as extracted from the complete data set (2008−2018) along with a Broken Power Law fit (solid red line). The regions show the 1σ confidence intervals resulting from the fit with BPL. The resulting break energy is E_{b} = 42.84 ± 1.78. 

In the text 
Fig. 2. γray light curve of PKS 0447−439 from August 2008 to December 2017 in the energy range between 100 MeV and 500 GeV in intervals of 30 days. The dotted line gives a fit with a constant flux. There is clear evidence for variability. Vertical error bars indicate the onesigma error bars. The gray shaded regions show the 1σ, 3σ, and 5σ confidence intervals resulting from the fit with a constant function. 

In the text 
Fig. 3. γray light curve of PKS 0447−439 from 2008 to 2018 in the energy range (0.1 − 500) GeV obtained using the aperture photometry technique for a radius of 1° around the source. A binning of one week is employed. Error bars correspond to 1σ error bars. The gray shaded regions show the 1σ, 3σ, and 5σ confidence intervals resulting from the fit with a constant function. 

In the text 
Fig. 4. Aperture photometry flux in comparison with the binned likelihood flux calculated for PKS 0447−439 (logarithmic scale). The red solid line represents a fit of the data with a linear function and the gray shaded regions show the resulting 3σ and 5σ confidence intervals. Error bars indicate the 1σ error interval. No background subtraction is used in the aperture photometry method. 

In the text 
Fig. 5. Normalized histograms of γray photon fluxes fitted with a lognormal (solid red line) and Gaussian (blue solid line), respectively. Apart from PG 1553+113, all sources show a clear preference for a lognormal distribution. 

In the text 
Fig. 6. Smoothed count map (logarithmic scale) of the PKS 0447−439 region in energy range from 100 MeV up to 500 GeV as seen by FermiLAT based on data between August 2008 and December 2017. The color bar has units of counts per pixel and the pixel dimensions are 0.1 × 0.1°. The locations of the background region and surrounding sources (labeled P1 = new source beyond the 3FGL catalog, P2 = 3FGL J0438.8−4519, P3 = 3FGL J0437.2−4713 and P4 = 3FGL J0455.7−4617) are shown as circles. The magenta circle shows the location of PKS 0447−439 and the yellow circles represent the background region. 

In the text 
Fig. 7. Generalized LSPs obtained based on light curves generated using aperture photometry for PKS 0447−439 and surrounding sources (P1,P2,P3,P4) shown as circles in Fig. 6. 

In the text 
Fig. 8. GLPS obtained using the monthly light curves of PKS 0447−439 is shown in cyan (circles); the red solid line, green line, and orange line represent the 3σ, 2.5σ, and 2σ confidence levels, respectively, calculated based on simulations of 50 000 light curves using the method of Timmer & Koenig (1995). 

In the text 
Fig. 9. Twodimensional WWZ of the γray light curve of PKS 0447−439 based on FermiLAT data from August 2008 to December 2017. The color scale represents the Zstatistics of the WWZ power of a certain period at a given time. The panel shows the signature of a quasiperiodic oscillation at a (observed) period of ∼2.6 years without any significant changes over time. 

In the text 
Fig. 10. Success fraction vs. slope (β) for all three PSDs (FermiLAT). The success fractions indicate the goodness of fit obtained from the PSRESP method (see text). (a) Weekly, (b) monthly light curves. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.