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

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.


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 θ ≤ 1/Γ, where Γ 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 γ-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 multiwavelength 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 ∼ 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 andPKS 0537-441) have been reported so far exhibiting longterm quasi-periodic variability in their γ-ray fluxes with significance higher than 3σ Sandrinelli et al. 2014;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 γ-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 Ztransform (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 Article number, page 1 of 12 arXiv:1901.10246v1 [astro-ph.HE] 29 Jan 2019 A&A proofs: manuscript no. aanda 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 γ-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.

Fermi-LAT likelihood analysis
Fermi-LAT on board the Fermi satellite is a electron-positron pair-conversion γ-ray detector sensitive to photon in 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 2 orbits. It has a point spread function (PSF) < 0.8 • and the largest effective area of approximately 8000 cm 2 above 1 GeV. Fermi-LAT data were downloaded from the Fermi-LAT Data Server 1 and events were selected between 2008 August 4 and 2018 April 1 (MET from 239557418.0 to 541555205.0), covering about ∼ 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 • 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= −2 ln(L 0 /L 1 ) was used to determine the source detection significance, with threshold set to TS= 25 (∼ 5σ). The significance of a source detection is given by ∼ √ T S σ (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 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 1 https://fermi.gsfc.nasa.gov/ssc/data/access/ 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 γ-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 ) 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.

Fermi-LAT γ-ray light curves
We generate γ-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 γ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 peak-to-peak 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). The aperture photometry method on the other hand, provides a model-independent 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.
In order to verify that the aperture photometry (AP) approach gives results which are consistent with the results obtained us- ing 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 ρ ∼ 0.936.

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 γ-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 γ-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 γ-ray fluxes were distributed in a histogram of fluxes. We fit all histograms in logscale, with Gaussian G(φ) and log-normal 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 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 lognormal 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 K 2 and the Shapiro-Wilk test, see Table 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.

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 cos(ωt) + B sin(ωt), and is defined for a uneven simpled 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 generalized LSP (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 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 in-terference 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 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. 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 ∼ 3-year high state a peak can be seen at periodicity timescale of ∼ 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

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 The WWZ is based on a weighed projection onto three trial functions, y(t) = i y i φ i (t) 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 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 Article number, page 7 of 12  Table 6. Summary results of the strongest periods for the Weighted Wavelet Z-transform method WWZ based on the monthly γ-ray light curves, including observed and intrinsic period.
with V x and V y as weighted variations of the data and the model function, respectively and N e f f 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 ± 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.

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 colorednoise-like behavior with power spectral density (PSD) charac-terized by a simple power-law of the form PSD(ν) ∼ ν −β where ν is the (temporal) frequency and β the power-law index.

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 where A is the amplitude of the model at the reference frequency ν 0 , β 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(ν) obs. ) and that of each of the simulated light curves (P(ν) sim,i where i = 1, N). The χ 2 dist,i 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 χ 2 dist for each simulated PSD P(ν) sim,i over all frequencies between ν min and ν max with respect to the data. We then compare χ 2 dist with the χ 2 obs distribution and count the number m of χ 2 dist,i for which χ 2 obs is smaller than χ 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 β 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 Article number, page 8 of 12

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 × 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 β. 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 the 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 narrowing-down of the PSD slope will be most relevant for assessing the real significance of the detected periods.

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.

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 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 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σ significance. These findings are in line with similar indications in (Covino et al. 2019).

Discussion and Conclusions
There are an increasing number of reports suggesting the presence of year-type periodicities in the γ-ray light curves of Fermi-LAT detected blazars. Given the still limited duration of the light curves ( < ∼ 10 yr), care needs to be taken, however, to properly  Table 7. The table shows the slopes β obtained from fitting the γ-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) (the success fractions indicate the goodness of fit obtained from the PSRESP method). The β values have been determined using both monthly and weekly light curves, as well as the Timmer & Koenig (1995) and Emmanoulopoulos et al. (2013) Timmer & Koenig (1995) and the Emmanoulopoulos et al. (2013) and based on on the monthly and weekly γ-ray light curves. the slopes β are obtained from fitting the γ-ray power spectra with simple power-law model, and the corresponding success fraction are calculated based on the PSRESP method (Uttley et al. 2002).
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 > ∼ 3σ) 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 ∼ 2.6σ to > 5σ (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 1995or 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 ∼ 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 7×10 3 (1+q) 1/3 /[q (M BH /5× 10 8 M ) 5/3 ] (P/2 yr) 8/3 yr only, so that for typical systems with q := m/M ≥ 0.05 the presumed source state would be highly unlikely. In fact, one rather expects orbital periods of the order The mean slopes β m are obtained from fitting the γ-ray power spectra with simple power-law model, and the corresponding success fraction are calculated based on the PSRESP method (Uttley et al. 2002). β min and β min represents the minimum and maximum of the 1σ error bars resulting from fitting Gaussians to each of the slope profile.
∼ 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 ∼ 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 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 ad (r t ) ∼ r g (0.5α c) −1 (r t /r g ) 3/2 , with α = 0.25 the viscosity coefficient and r g = GM 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 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, 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 ∼ 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, sligthly changing footpoint radii could possibly compensate for this. The lighthouse model was originally designed to account for observed QPOs with P obs ≤ 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 ∼ 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., ∼ 2 − 3 yr) is expected to significantly improve the situation and to help to clarify their putative presence and physical implication.