A γ-ray Quasi-periodic Modulation in the Blazar PKS 0301–243?

Blazars are a subclass of radio-loud active galactic nuclei 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 through a rotating accretion disk. The emissions from blazars are dominated by the non-thermal emission from the relativistic jet, extending from MHz radio frequencies to TeV γ-ray energies, and also exhibit variabilities at all energies on a wide range of timescales (e.g., Liao et al. 2014; Liao & Bai 2015). 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 at MeV to GeV energies.

The periodic variabilities of blazars have been investigated extensively in the optical band (e.g., Bai et al. 1998, 1999; Fan & Lin 2000; Xie et al. 2008; Li et al. 2009; Urry 2011; King et al. 2013; Zhang et al. 2014; Bhatta et al. 2016; Fan et al. 2016). An interesting case is that of OJ 287, which has a ∼12 yr period cycle (Kidger et al. 1992; Valtonen et al. 2006). Searches for γ-ray quasi-periodic oscillations (QPOs) became possible after the launch of the Fermi Gamma-ray Space Telescope (Fermi) in 2008. So far, the Large Area Telescope (LAT; Abdo et al. 2009; Atwood et al. 2009) on board Fermi has collected gamma-rays for over 8 years. The possible quasi-periodic variabilities of blazars with periods of ∼2–3 yr have been reported in γ-ray fluxes of several blazars (e.g., Sandrinelli et al. 2014; Ackermann et al. 2015; Sandrinelli et al. 2016a, 2016b, 2017; Zhang et al. 2017a, 2017b). In particular, PKS 2155–304 has been reported to have a significance of ∼4σ (Zhang et al. 2017b). The quasi-periodic modulations in the blazars carry important information on the BH—jet system.

We present the results of searching for QPOs in the γ-ray light curves of blazar PKS 0301–243. A clear quasi-periodic signal with a period cycle of ∼2.1 yr is found at a 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. A detailed LAT data analysis and the main results are reported in Section 2. In Section 3 we summarize the results and provide a brief discussion of the findings.

PKS 0301–243 is a high-synchrotron-peaked blazar (HSP) with a synchrotron peak frequency $\approx {10}^{15}$ Hz (Abramowski et al. 2013), and its redshift is 0.266 (Pita et al. 2012). The High Energy Stereoscopic System (H.E.S.S.) has detected TeV photons from this source (Abramowski et al. 2013).

The events are collected between 2008 August 4 and 2017 May 19 (Modified Julian Date, MJD: 54,682.66—57,892.66) in the energy range from 100 MeV to 500 GeV, and in a square region of interest (ROI) of $20^\circ \times 20^\circ$ centered at the position of PKS 0301–243. The position of the target is located at right ascension (R.A.) = ${03}^{{\rm{h}}}03\buildrel{\rm{m}}\over{.} 442$, declination (decl.) = $-{24}^{{\rm{h}}}\,07\buildrel{\rm{m}}\over{.} 192$ (J2000; $l=214.621,\,b=-60.177$). The analysis is performed with the Fermi Science Tools version v10r0p5 package, which is provided at the Fermi Science Support Center (FSSC).⁷ The Pass 8 LAT data (Atwood et al. 2013) are used, but keeping only the SOURCE class photon-like events (with options evclass = 128 and evtype = 3 in the tool gtselect). To minimize the contamination due to the gamma-ray-bright Earth limb, we exclude the events with zenith angles $\geqslant 90^\circ$. By running the tool gtmktime, we obtain good time intervals with high-quality photons. All the data reductions follow the data analysis thread provided by FSSC.⁸ We adopt the instrumental response function (IRF) "P8R2_SOURCE_V6" in the analysis. Two diffuse model files,⁹ 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 time range, with the model file generated with the script make3FGLxml.py. This file contains information on the spectral parameters of all the known 3FGL sources (Acero et al. 2015) in the ROI. The γ-ray spectrum of the target is a power law in the Fermi 3FGL. The best-fitting results are derived with the Fermi tool gtlike, and are saved as a new model file. We also use the spectra in the 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\mathrm{GeV}}=(4.2\pm 0.1)$ × ${10}^{-8}\,\mathrm{photons}\ {\mathrm{cm}}^{-2}\ {{\rm{s}}}^{-1}$, and the index of the power law is 1.90 ± 0.01, with a TS value of 9391.7 (the results in this paper feature statistical errors only). We construct the light curves based on this new model file.

2.1.  $\gamma$-ray Light Curve

We use two techniques, maximum likelihood optimization (ML) and 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 circular ROI of 15° centered at the coordinates of the target. We use the same parameter value as that used 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 upper left panel of Figure 1. To test whether the power peaks vary with different lengths of time bins, we also produce the ML light curve with a 10-day bin, which is shown in the upper left panel of Figure 2. The light curve can also be produced by 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 a 1° radius centered on the position of the target for each 2.5-day bin, in which the counts are weighted by their relative exposure for each time bin. The AP light curve is shown in the upper left panel of Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We note that there is an isolated large flare around MJD 55320. To avoid its impact on searching for quasi-periodic variability, we remove this flare in the following quasi-periodicity analyses.

2.2. Searching for Quasi-periodic Variability

2.2.1. Analyses of $\gamma$-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 light curve, three power spectra, LSP power, WWZ power, and time-averaged WWZ power, are shown in Figure 1. A strong peak near a period cycle of 2.1 ± 0.3 yr appears, in which the maximum power is >18.6 times that of the mean power value. The probability (Prob) for obtaining a power larger than the maximum power from the noise is $\lt 1.58\times {10}^{-9}$ (corresponding to a $\gt 6.0\sigma$ significance level). The Prob (P > Pn) is assessed through the formula: $\mathrm{Prob}{\left(P,\gt ,\mathrm{Pn},),=,(,1,-,2,\times ,\tfrac{\mathrm{Pn}}{N-1}\right)}^{({\rm{N}}-3)/2}$, with the normalization from Horne & Baliunas (1986), where N = 105 is the number of time bins in the month-bin light curve. We correct the probability in the range of $1/3000\,{\mathrm{day}}^{-1}-1/60\,{\mathrm{day}}^{-1}$ with the "trial factor = 50" (the number of sampled independent frequencies) (Zechmeister & Kürster 2009), and find that the false-alarm probability is less than $7.8\times {10}^{-8}$, corresponding to $\gt 5.4\sigma$. By fitting the power peak with a Gaussian function, we derive a period cycle of 763.3 ± 114.9 days. The uncertainty of the period is evaluated based on the half width at half maximum 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⁶ times, based on the obtained best-fitting result of the power spectral density (with the form of $P(f)\sim 1/{f}^{\alpha }+c$, where c represents the Poisson noise level) and the probability density function of observed variation. We then derive significance curves of 5σ and 4σ based on the simulations, which are shown in the lower right panel of Figure 1. The significance of the signal is $\simeq 5.4\sigma$. We also calculate the power spectra of the 10-day ML light curve and 2.5-day bin AP light curve, which are shown in Figures 2 and 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 quasi-periodic signal, we fit the γ-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,¹⁰ is the best-fit model. In Figure 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 a 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 stochastic processes.

Zoom In Zoom Out Reset image size

Table 1.  AIC Values for the ARIMA Models Fitting the 10-day Bin Light Curve

		MA (q)
ARIMA(p, d, q)		d = 0						d = 1
		MA(0)	MA(1)	MA(2)	MA(3)	MA(4)	MA(5)	MA(0)	MA(1)	MA(2)	MA(3)	MA(4)	MA(5)
	AR(0)	1475	1388	1370	1368	1366	1368	1425	1391	1365	1363	1365	1367
	AR(1)	1359	1361	1359	1359	1361	1363	1410	1363	1364	1365	1366	1368
AR(p)	AR(2)	1361	1363	1360	1360	1362	1364	1397	1364	1363	1366	1360	1369
	AR(3)	1362	1360	1355	1362	1364	1364	1387	1365	1366	1368	1370	1371
	AR(4)	1362	1363	1362	1364	1365	1366	1383	1367	1368	1370	1360	1371
	AR(5)	1363	1362	1364	1356	1360	1368	1373	1368	1368	1370	1371	1373

Download table as:  ASCIITypeset image

We fold the events within a square ROI of $20^\circ \times 20^\circ$ 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 using the above best-fitting model file to obtain the phase-resolved likelihood results. In Figure 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 Figure 5); but no variability appears in its spectral shape (see the lower panel of Figure 5).

Zoom In Zoom Out Reset image size

2.2.2. Analyses of Optical and X-Ray Data

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

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Possible γ-ray QPOs have been reported for several blazars (e.g., Sandrinelli et al. 2014; Ackermann et al. 2015; Sandrinelli et al. 2016a, 2016b, 2017; Zhang et al. 2017a, 2017b). 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 $\sim 5\sigma$ in PKS 0301–243. No quasi-periodic 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 (Ackermann et al. 2015). 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. 2017b). The lack of optical and X-ray quasi-periodic variabilities may be because the optical/X-ray and γ-ray emissions originate from different regions. If the lack of optical and X-ray quasi-periodic variabilities is confirmed by future long-term monitoring, this 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. 2017b).

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 due to jet precession/rotation.

Note that in our case, the γ-ray photon index does not vary with the phase (see the lower panel of Figure 5). The gamma-ray photon index is mainly determined by the high-energy electron 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 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}_{\mathrm{int}}={P}_{\mathrm{obs}}/(1+z)$. The central SMBH of PKS 0301–243 is $\sim 8\times {10}^{8}\ {M}_{\odot }$ (Ghisellini et al. 2010). Assuming a total mass of the binary SMBH of ${10}^{9}\ {M}_{\odot }$, the binary system size would be $\sim 0.006$ pc. At this stage, gravitational wave emission would 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. A systematic sample study on QPOs at different electromagnetic frequencies in blazars could reveal the deep physics of the jet (e.g., Sandrinelli et al. 2016a).

The γ-ray QPO in PKS 0301–243 is the first detection of its kind of a signal in blazars at a confidence level of $\sim 5\sigma$. Since there were barely four nominal quasi-periods in the current Fermi-LAT data, this result certainly requires confirmation. Fortunately, our claim for a QPO should be tested rather soon, as the next flux maximum is expected in 2018.

Finally, we would like to stress that these claimed γ-ray QPOs in blazars are different from the X-ray QPOs in BH X-ray binaries and narrow-line Seyfert 1 galaxies (Zhang et al. 2017a). 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 masses to SMBHs. No such relation is found in γ-ray QPOs in blazars (Figure 8). It seems that the intrinsic period of γ-ray QPOs in blazars is independent of the SMBH mass. Moreover, the relation of the γ-ray QPO frequency to BH mass significantly deviates from the inverse relation found in the X-ray QPOs. The X-ray and γ-ray QPOs provide us with different insights into the BH—jet system.

Zoom In Zoom Out Reset image size

We thank the anonymous referees for useful and constructive comments. Part of this work is based on archival data, software, or online services provided by the ASI Science Data Center (ASDC). We acknowledge the financial support from the 973 Program of China under grant 2013CB837000, the National Natural Science Foundation of China (NSFC-11433004, NSFC-11525313, NSFC-11573060, NSFC-11573026, and NSFC-11661161010), and the Key Laboratory of Astroparticle Physics of Yunnan Province (No. 2016DG006). D. H. Yan is very grateful to Prof. Yefei Yuan (USTC) for helpful discussions. The work of D. H. Yan is supported by the CAS "Light of West China" Program.