A $\gamma$-ray Quasi-Periodic modulation in the Blazar PKS 0301$-$243?

We report a nominally high-confidence $\gamma$-ray quasi-periodic modulation in the blazar PKS 0301$-$243. For this target, we analyze its \emph{Fermi}-LAT Pass 8 data covering from 2008 August to 2017 May. Two techniques, i.e., the maximum likelihood optimization and the exposure-weighted aperture photometry, are used to build the $\gamma$-ray light curves. Then both the Lomb-Scargle Periodogram and the Weighted Wavelet Z-transform are applied to the light curves to search for period signals. A quasi-periodicity with a period of $2.1\pm0.3$ yr appears at the significance level of $\sim5\sigma$, although it should be noted that this putative quasi-period variability is seen in a data set barely four times longer. We speculate that this $\gamma$-ray quasi-periodic modulation might be evidence of a binary supermassive black hole.


INTRODUCTION
Blazars are a subclass of radio-loud active galactic nuclei (AGNs) whose relativistic jets almost point at observers (Urry & Padovani 1995). It is generally thought that a central supermassive black hole (SMBH) provides the energy that powers the relativistic jet through BH spin or rotating accretion disk. The emissions from blazar are dominated by the nonthermal emission from the relativistic jet, extending from MHz radio frequencies to TeV γ-rays energies, also exhibiting variabilities at all energies on a wide range of timescales. The typical mutiwavelength spectral energy distribution (SED) is distinguished by two broad peaks: a synchrotron component peaking at infrared to X-ray bands and a Compton component peaking in MeV to GeV energies.
We present the results of searching for QPO in the γray light curves of blazar PKS 0301−243. A clear quasiperiodic signal with a period cycle of ∼ 2.1-yr is found at the significance level of ∼5σ, though given that the full data set is only 8.78 years long this signal can easily have arisen randomly (e.g. Press 1978). The paper is organized as follows: the detailed LAT data analysis and the main results are reported in Section 2. In Section 3 we summary the results and present a brief discussion on the findings.
The events are collected between 2008 August 4 and 2017 May 19 (Modified Julian Date, MJD: 54,892.66) in the energy range from 100 MeV to 500 GeV, and in a square region of interest (ROI) of 20 • × 20 • centered at the position of PKS 0301−243. The position of the target is located at right ascension (R.A.) = 03 h 03.442 m , declination (decl.) = −24 h 07.192 m (J2000; l = 214.621, b = −60.177). The analysis is performed with the Fermi Science Tools version v10r0p5 package which is provided in the Fermi Science Support Center (FSSC). 7 The Pass 8 LAT data (Atwood et al. 2013) are used with keeping only the SOURCE class photonlike events (with options evclass = 128 and evtype = 3 in the Tool gtselect). To minimize the contamination due to the gamma-rays bright Earth limb, we exclude the events with zenith angles 90 • . By running Tool gtmktime, we obtain the good time intervals with high-quality photons. All the data reductions follow the data analysis thread provided by FSSC 8 . We adopt the instrumental response function (IRF) 'P8R2 SOURCE V6' in the analysis. Two diffuse model files 9 , namely gll iem v06.fit and iso P8R2 SOURCE V6 v06.txt, are used to model the Galactic and extragalactic diffuse γ-rays . A binned maximum likelihood is adopted to fit the events in the whole time range with the model file generated with the script make3FGLxml.py. This file contains the information on the spectral parameters of all known 3FGL sources (Acero et al. 2015) in the ROI. The γ-ray spectrum of the target is power-law in the Fermi 3FGL. The best-fitting results are derived with Fermi Tool gtlike, and are saved as a new model file. We also use the spectra in Fermi 3FGL model file to fit the events in the square ROI. The integrated photon flux of the best-fitting results above 100 MeV is F 0.1−500 GeV = (4.2 ± 0.1) × 10 −8 photons cm −2 s −1 , and the index of power law is 1.90 ± 0.01 with the TS value of 9391.7 (the results in this paper with statistical errors only). We construct the light curves based on this new model file.

γ-ray light-curve
We use the two techniques, the maximum likelihood optimization (ML) and the exposure-weighted aperture photometry (AP; Corbet et al. 2007;Kerr 2011), to construct the γ-ray light curves of PKS 0301−243. The 30-day-bin ML light curve is generated by employing the unbinned maximum likelihood fitting technique. In this step, the Tool gtlike is employed for each time-bin, and the events are selected in a circle ROI of 15 • centered at the coordinates of the target. We use the same parameter value as that in the new model file for all the sources in the ROI and freeze the spectral parameters except for the target. The 30-day-bin ML light curve is shown in the left upper panel of Fig. 1. For testing whether the power peaks vary with different length of time-bin, we also produce the ML light curve with the 10-day-bin, which is shown in the left upper panel of Fig. 2. The light curve can also be produced by the method of the exposure-weighted aperture photometry. In this method, we calculate the probabilities for each photon with the Fermi -Tool gtsrcprob, and then sum the probabilities of each photon within 1 • radius centered on the position of target for each 2.5-day-bin, in which the counts are weighted by its relative exposure for each time-bin. The AP light curve is shown in the left upper panel of Fig. 3.
We notice that there is an isolated large flare around MJD 55320. To avoid its impact on searching for quasiperiodic variability, we remove this flare in the following quasi-periodicity analyses.

Analyses on γ-ray data
We apply the two widely used methods, Lomb-Scargle Periodogram (LSP;Lomb 1976;Scargle 1982) and Weighted Wavelet Z-transform (WWZ ;Foster 1996), to the γ-ray light curves. For the 30-day-bin ML lightcurve, three power spectra, LSP power, WWZ power and time-averaged WWZ power, are shown in Fig. 1. A strong peak near a period cycle of 2.1 ± 0.3 yr appears, in which the maximum power is > 18.6 times of the mean power value. The probability (Prob) for obtaining a power larger than the maximum power from the noise is < 1.58 × 10 −9 (corresponding to a > 6.0σ significance level). The Prob(P > Pn) is assessed through the formula: Prob(P > Pn) = (1 − 2 × Pn N−1 ) (N−3)/2 with the normalization from Horne & Baliunas (1986), where N=105 is the number of time-bin in the month-bin light curve. We correct the probability in the range of 1/3000 day −1 -1/60 day −1 with the "trial factor = 50" (the number of sampled independent frequencies) (Zechmeister & Kürster 2009), and find that the false-alarm probability (FAP) is less than 7.8 × 10 −8 , corresponding to > 5.4σ. By fitting the power peak with Gaussianfunction, we derive the period cycle of 763.3±114.9 days. The uncertainty of the period is evaluated based on the half width at half maximum (HWHM) of the Gaussian fitting.
In order to evaluate the precise significance of the signal, we use the method in Emmanoulopoulos et al. (2013) (also see Ackermann et al. 2015;Bhatta et al. 2016) to simulate light curves 3 × 10 6 times based on the obtained best-fitting result of power spectral density (with the form of P (f ) ∼ 1/f α +c, where c represents the Poisson noise level) and the probability density function of observed variation. We then derive the significance curves of 5-σ and 4-σ based on the simulations, which are shown in the lower right panel of Fig. 1. The significance of the signal is 5.4 σ. We also calculate the power spectra of the 10-day ML light curve and 2.5-daybin AP light-curve, which are shown in Fig. 2 and Fig.  3, respectively. In these two power spectra, we also find strong signals at ∼2.1 yr.
In order to further check the reliability of the quasiperiodic signal, we fit γ-ray light curve with autoregressive integrated moving average (ARIMA) models (Box & Jenkins 1976;Hamilton 1994;Chatfield 2003) to assess whether the signal is consistent with a stochastic origin of autoregressive noise. We use the Akaike Information Criterion (AIC; Akaike 1973) to select the best-fit model. In Table 1, we show the AIC values for 72 ARIMA models fitting the 10-day-bin γ-ray light curve. One can see that the ARIMA (3,0,2) model, with the minimum AIC value of 1355 10 , is the best-fit one. In Fig. 4, we show the standard residuals and the auto-correlation function (ACF) of the residuals for the best-fit model. It can be seen that there is a spike at the lag of 660 days that exceeds the 95% confidence limit. This marginal evidence indicates that the γ-ray quasi-periodic variability may not be produced by such type of stochastic processes.
We fold the events within a square region of interest (ROI) of 20 • × 20 • centered at the position of PKS 0301−243 into 15 uniform bins based on orbital phase with the phase zero corresponding to MJD 54,682.66. We then fit the data in each phase bin by using the above best-fitting model-file to obtain the phase-resolved likelihood results. In Fig. 5, one can see that this folded light curve varies with the phase, indicating substantial variability in the source brightness (see the upper panel of Fig. 5); but no variability appears in its spectral shape (see the lower panel of Fig. 5).

Analyses on optical and X-ray data
We also search for quasi-periodic signal in the optical and X-ray data from this source. The long-term optical data from the Catalina Sky Surveys covering from 2005 October to 2013 October and daily averaged X-ray data from Swift-BAT covering from 2005 February to 2017 January are shown in the upper panels of Fig. 6 and Fig. 7, respectively. The LSP powers of the optical data and X-ray data are respectively shown in the lower panels of Fig. 6 and Fig. 7. No obvious peak is found in the corresponding powers. It is noted that the X-ray data are weakly variable.

SUMMARY AND DISCUSSION
Possible γ-ray QPOs have been reported in several blazars (e.g., Sandrinelli et al. 2014;Ackermann et al. 2015;Sandrinelli et al. 2016aSandrinelli et al. ,b, 2017Zhang et al. 2017a,b). However, the significance of the claimed QPOs is not very high. In this paper, we report the first detection of γ-ray quasi-periodic modulation at a nominal confidence level of ∼ 5σ in PKS 0301−243. No quasiperiodic modulation is found in its optical and X-ray data.
In PG 1553+113, the quasi-periodic variabilities in optical and γ-ray data have the same period cycle . In PKS 2155−304, the periods of optical and γ-ray quasi-periodic variabilities are different (Sandrinelli et al. 2014). In PKS 0426−380, no optical quasi-periodic variability is found (Zhang et al. 2017a). The lack of optical and X-ray quasi-periodic variabilities may be because of the optical and X-ray originating from the different region that does not contribute γ-rays. If the lack of optical and X-ray quasi-periodic variabilities is confirmed by futuer long-term monitoring, it would challenge the most popular one-zone blazar emission model in which optical, X-ray and γ-ray emissions are assumed to be produced in the same region (Zhang et al. 2017a).
The mechanism causing the γ-ray quasi-periodic modulation in blazars is poorly understood. Given that the γ-rays are produced in the jet, two possibilities may account for the γ-ray quasi-periodic variabilities in blazars (e.g., Ackermann et al. 2015): (i) pulsational accretion flow instabilities may induce a quasi-periodic injection of plasma into the jet, hence a quasi-periodic modulation appears in the γ-ray flux from the jet; and (ii) the Doppler magnification factor changes periodically caused by jet precession/rotation.
Note that in our case the γ-ray photon index does not vary with the phase (see the lower panel of Fig. 5). The gamma-ray photon index is mainly determined by the high-energy electrons distribution. For HSP, the electron cooling is inefficient (e.g., Ghisellini & Tavecchio 2008;Yan et al. 2014), and the electron distribution is mainly governed by the acceleration mechanism in the jet. This result indicates that the process yielding the QPO in the γ-ray flux would not have an impact on the acceleration process. The first possible origin for the QPO outlined above would have an impact on the energy outflow efficiency which is relative to the acceleration process in the jet (e.g., Ackermann et al. 2015). Therefore, our results may prefer to the second origin, i.e., jet precession. The jet precession could be the result of a helical jet (e.g, Rieger 2004;Komossa & Zensus 2016). Furthermore, a binary SMBH system would be involved in the formation of a helical jet (e.g, Komossa & Zensus 2016;Sobacchi et al. 2017). Within such a scenario, the observed 2.1 yr period is the orbital time, and the equivalent intrinsic orbital time P int =P obs /(1 + z). The central SMBH of PKS 0301−243 is ∼ 8 × 10 8 M (Ghisellini et al. 2010). Assuming the total mass of the binary SMBH of 10 9 M , the binary system size would be ∼ 0.006 pc. At this stage, gravitational wave emission may be non-negligible in carrying away the energy.
In the jet precession model, the issue of the lack of optical and X-ray quasi-periodic variabilities could be resolved if the optical and X-ray radiations originate from a large region where the Doppler boosting is weak. Systematic sample study on QPOs at different electromagnetic frequencies in blazars may reveal deep physics of the jet (e.g., Sandrinelli et al. 2016a).
The γ-ray QPO in PKS 0301−243 is the first detection of such a kind of signal in blazars at a confidence level of ∼ 5σ. Since there were barely four nominal quasi-periods in the currently Fermi-LAT data, this result certainly requires confirmation. Fortunately, our claim for a QPO should be tested rather soon, as the next flux maximum would be expected in 2018.
Finally, we would like to stress these claimed γ-ray QPOs in blazars are different from the X-ray QPO in BH X-ray binaries and narrow-line Seyfert 1 galaxy (Zhang et al. 2017b). For the X-ray QPOs, there is an inverse linear relation between QPO frequency and BH mass (e.g., Abramowicz et al. 2004;Török 2005;Remillard & McClintock 2006;Pan et al. 2016). This relation spans from stellar-mass to SMBH. No such relation is found in γ-ray QPO in blazars (Fig. 8). It seems that the intrinsic period of γ-ray QPO in blazars is independent on the SMBH mass. Moreover, the relation of the γ-ray QPO frequency-BH mass significantly deviates from the inverse relation found in the X-ray QPOs. The X-ray and γ-ray QPOs provide us different insights into the BH -jet system.