A 34.5 day quasi-periodic oscillation in γ-ray emission from the blazar PKS 2247–131

Since 2016 October, the active galaxy PKS 2247−131 has undergone a γ-ray outburst, which we studied using data obtained with the Fermi Gamma-ray Space Telescope. The emission arises from a relativistic jet in PKS 2247−131, as an optical spectrum only shows a few weak absorption lines, typical of the BL Lacertae sub-class of the blazar class of active galactic nuclei. Here we report a ≃34.5 day quasi-periodic oscillation (QPO) in the emission after the initial flux peak of the outburst. Compared to one-year time-scale QPOs, previously identified in blazars in Fermi energies, PKS 2247−131 exhibits the first clear case of a relatively short, month-like oscillation. We show that this QPO can be explained in terms of a helical structure in the jet, where the viewing angle to the dominant emission region in the jet undergoes periodic changes. The time scale of the QPO suggests the presence of binary supermassive black holes in PKS 2247−131.

G alaxies all appear to contain a supermassive black hole (SMBH) at their centers, with a mass usually in the range of~10 6 -10 10 M ⊙ (where M ⊙ is the solar mass) 1 . The SMBH is sometimes fed by enough matter, accreted through a surrounding disc, to produce strong emissions from the central region of the galaxy that can sometimes outshine the emission from all the stars in the galaxy. These Active Galactic Nuclei (AGN) 2 naturally are objects of great interest, and relativistic jets, launched from the immediate vicinity of the SMBHs, are found to be associated with~10% of AGN 3 . If a relativistic jet happens to be pointing close to our line of sight, the emission observed from the jet dominates other contributions over nearly the entire electromagnetic spectrum, because of the special relativistic (Doppler) beaming effect. For the same reasons, such AGN, classified as blazars, show rapid and large-amplitude variability 4,5 .
Studies of the variability of jets allow us to explore their structures, physical properties, and radiation processes 4,5 . Although rare, one particularly intriguing phenomenon revealed by some of these studies are quasi-periodic oscillations (QPOs). QPOs may be interpreted as an evidence of the binary nature of the central SMBH system, such as in the most well-known case OJ 287 6,7 , or reflect the innermost stable orbit of a black hole or oscillation modes of the surrounding accretion disc 8,9 . Therefore, searches for QPOs from blazars have been conducted at different wavelengths. At γ-rays, the Large Area Telescope (LAT) on board the Fermi Gamma-Ray Space Telescope (Fermi) has provided a powerful tool for monitoring blazars at γ-ray energies 10 . Candidate γ-ray QPOs have been discovered; [11][12][13][14][15][16][17][18] these QPOs all have yearly time scales. Here we report our discovery of a possible month-long QPO in γ-ray emission of the blazar PKS 2247−131; the only previous claim of a QPO on this time-scale was reported to have been detected at γ-ray TeV energies during a flare of the blazar Markarian 501 in 1997 19,20 , but fewer cycles were observed then and the strength of that oscillatory signal is weak.

Results
Quasi-periodic oscillation signal. In 2016 October, a γ-ray outburst event in the direction of PKS 2247−131 was detected with Fermi LAT 21 , and based on multi-wavelength observations, we confirmed the outburst to have originated from PKS 2247 −131 (see Methods). Moreover, this source is a blazar-type active galactic nucleus (AGN); more specifically a BL Lacertae sub-type. The γ-ray light curve of PKS 2247−131, constructed by binning 5-day Fermi LAT data in the energy range 0.1-300 GeV, shows a clear periodic modulation after the initial main flux peak (around MJD 57675; see Fig. 1). Our initial search for periodic modulation in the light curve indicated a period of~34 days, the temporal duration of which ranges from MJD 57693 to 57903 (Fig. 1). These results were obtained from the Weighted Wavelet Ztransform 22 (WWZ) of the entire light curve (see the left panel of Fig. 2). The WWZ method on the light curve during MJD 57693-57903 provided clearer results (see the right panel of Fig. 2). The WWZ power, as well as its time-averaged power, strongly indicates the presence of a~34 day QPO.
In addition, the Lomb-Scargle periodogram 23,24 (LSP) of the light curve during MJD 57693-57903 provided similar results (see the right panel of Fig. 2). From the LSP power density curve, the oscillation period P obs was determined to be P obs =34.5 ± 1.5 days.  We used the light curve simulation method 25 to estimate the significance of the oscillation signal, in which the noise in the LSP power density curve was modeled with a smoothly bending power law (see Supplementary Methods and Supplementary Figures 1-3 for details). The 4σ and 5σ significance curves were obtained from simulating the light curves 10 7 times (see the right panel of Fig. 2). The results indicate a significance of 5.2σ for this signal, or 4.6σ after considering a trial number of 20 (see Supplementary Methods). We thus concluded that a~4.6σ QPO exists in the light curve during MJD 57693-57903. This is the first clear case that a comparatively short, month-long QPO has been seen in γ-ray emission emanating from a blazar.
Within the time span of the periodic modulation, six full cycles are visible. We modeled the data in each 5-day time bin with a power law NðEÞ $ E ÀΓ p (where E and Γ p are the photon energy and photon index, respectively), and no significant spectral variations were found. The photon indices during the six cycles have an average of 2.26 with a standard deviation of 0.14, while none of them are significantly away from the 2.26 ± 0.14 range (see Fig. 1). We constructed a folded light curve by performing the likelihood analysis on the data of the six cycles in each of 16 phase ranges (phase zero was set at MJD 57692.66; see Fig. 3). The light curve is not symmetric, dropping gradually from the peak to the lowest point of the curve, and a second component appears to exist, as hints of the second component may be seen in the individual modulations. More detailed analysis shows that the second component actually has a flux peak nearly as high as that of the main one in the energy range 1-300 GeV, and the different heights seen between the two components in the folded light curve in the whole energy range are because of somewhat lower fluxes of the second component in the low energy portion, <1 GeV, although their spectral shapes are similar (see Supplementary Figures 4,5). We note that two Lorentz functions provide a good fit to the light curve profile ( Fig. 3; the details of the Lorentz functions are given in Supplementary Table 1).

Discussion
The recently identified candidate γ-ray QPOs [11][12][13][14][15][16][17][18] , having yearlong periods, were often discussed in scenarios involving a binary SMBH AGN system. These are plausible, as if the binary black holes have a total mass of $ 10 8 M and the separation, R, between them is at a milli-parsec scale (~100 times the Schwarzschild radius R s = 2GM/c 2 , where G is the gravitational constant and c is the speed of light, for a mass M ¼ 10 8 M ), the binary orbital period is in the range of several years. The companion black hole could perturb the accretion periodically at the orbital period 7 . Second, the companion black hole could also induce long-term periodicities by exerting a torque on the accretion disc 26 . However, the driven preccesion of the accretion disc (and the jet) has a time scale of hundreds of years. On the other hand, single SMBH systems can give rise to short-term, sub-day QPOs, resembling some QPO cases seen in stellar-mass black hole binaries in the Milky Way 8 . Pulsational accretion flow instabilities may also drive jet-disc QPOs, which are predicted to have intra-day time scales 27,28 .
The 34.5 day QPO we found in PKS 2247−131 is obviously different from the time scales considered in the above scenarios. A geometrical origin instead is the likely explanation. At radio wavelengths, Very Long Baseline Interferometry (VLBI) has revealed that in some blazars, the parsec-scale cores appear to be misaligned with the large, kiloparsec-scale structures of jets. Such misalignment can be naturally explained by a helical structure in these inner jets 29 . Significant flux variations can arise due to the changing relativistic beaming effect when different parts of such a helical jet pass closest to the line of sight, even if no intrinsic variations of the emission from the jet are present 30 . Furthermore, our viewing angle to the helical motion changes essentially periodically as the emission blob in a jet moves towards us, thus giving rise to QPO-like flux modulations (see a schematic illustration in Fig. 4).
In this geometrical scenario, the viewing angle θ of an emitting blob's motion as a function of phase t (normalized to the minimum value of θ = θ min at t = 0) 31 is cosθ = sin ϕ sin ψ cos(2πt) + cosϕcosψ, where ϕ is the pitch angle between the emitting blob's motion and the jet's axis, and Ψ is the inclination angle of the jet to the line of sight (Fig. 4). We then have the relativistic beaming is the bulk Lorentz factor of the moving blob, which has a bulk speed of βc. The flux intensity I(ν) in the observer's frame is transformed from I′(ν′) in the blob's comoving frame by 32 νIðνÞ ¼ δ 4 ν′I′ðν′Þ. We may set ϕ = 2°, as it usually is in the range of ≤1°−4°3 1 . The other parameters that provide a decent fit are ψ = 5°and Γ = 8.5, both of which are rather typical for blazars. This simple model can produce a light curve that well describes the main part of the observed modulation, but the model light curve is symmetric and does not provide a fit to the right wing of the modulation (Fig. 3). If we consider that there are two components represented by the Lorentz functions used to fit, the model can be adjusted (with only the normalization factor changed) to provide a fit to the observed main modulation (Fig. 3). The origin of the two components is not clear, their difference being the second one has lower fluxes at <1 GeV energies, but we suspect that there are two significant emission regions, which may correspond to a forward shock and a reverse shock 33 . Evidence supporting this geometrical origin over, for example, one involving emission changes resulting from evolution of the emitting particle energies, is provided by the essential constancy of the photon power-law indices (Fig. 1), even when the flux varied periodically by a factor of~6. We note that this geometrical origin might also be revealed from multi-wavelength observations in optical and X-ray bands; however, the data we have collected are not only too sparse but unfortunately do not cover the long modulation cycles (see Supplementary Table 2).
Due to the effect of special relativity, the observed period P obs is much shorter than the physical period P at the host galaxy: 34 P obs = (1−β cosϕ cosψ)P. Using the above parameters, P ≈ 7.1 years. The distance that the blob moves in one cycle would be approximately D = cβP cosϕ ≈ 2.2 pc, and the total projected distance in six cycles would be D p ¼ 6Dsinψ % 1:1 pc. Such helical motions of newly born radio emitting knots, along with the parsec-scale jets, have been detected in several blazars [35][36][37] . For PKS 2247−131, because its distance is large,~1100 Megaparsec (estimated from the redshift z ≃ 0.22; see Methods), the helical motion would not be resolved by VLBI.
A helical structure could be intrinsic in a jet, driven by a coiled magnetic field 38,39 . This structure is supported by optical polarization changes seen in jets 40 . In our case, we note that the putative physical period P is actually on the order of orbital periods of close binary SMBHs, but boosted to appear much shorter, and thus the orbital motion of such a binary could be the cause of the helical motion seen in PKS 2247−131 41 . Due to the orbital velocity v of the jet-launching from the vicinity of the primary black hole, the direction of the ejected blob (at a velocity close to c) would have a small angle Δα from the axis of the accretion disc (or the black hole spin) 31 Δα ≃ v/c ≃ 1.3°q/(q + 1) 1/2 (R/1000R s ) −1/2 , where q is the mass ratio of the secondary black hole relative to the primary. As a result, the jet appears to be rotating with respect to an observer, and Δα corresponds to the pitch angle ϕ in the above helical scenario.
Our analysis of the light curve of PKS 2247−131 may appear to be the first time that a jet has been found to exhibit a helical structure through clear flux modulation of its high-energy flux. However, the complexity of the QPOs in blazars or AGN in general must be noted. Our explanation is based only on the monthly modulation time scale, but other QPO driving mechanisms 11,14 , such as those mentioned above, may not be absolutely excluded. The binary SMBH possibility is often considered as the cause for yearly QPOs claimed to be found in different searches, particularly at optical bands, but the number of these binaries should be limited, and are constrained by the gravitational wave background 42 . Even in the binary SMBH scenario, the cause of flux modulation in a blazar may be the jet  Fig. 4 Schematic illustration of a helical jet that produces periodically modulated emission. The emitting blob's motion has a pitch angle ϕ from the jet's axis, which has an inclination angle ψ from the line of sight. As the emitting blob moves towards the observer, the viewing angle to the blob changes periodically instabilities, dynamically perturbed by the orbiting companion 43,44 . The helical structure explanation could be confirmed by nearfuture multi-wavelength monitoring of PKS 2247−131 since the QPO should appear again, although the modulation is extremely sensitive to our viewing angle to the helical motion. Searches for similar QPOs in the γ-ray band can also help as their results could constrain the frequency of the appearance of such QPOs. Previous systematic searches in Fermi LAT data have found several candidates with yearly QPOs 45 . Thus far, more than 1500 blazars have been detected with Fermi LAT 10 . We have analyzed most of the LAT data for these known blazars and only 19 of them (see Supplementary Table 3) had γ-ray fluxes as high as that seen during the outburst of PKS 2247−131. However, none of the 19 sources have been found to exhibit similar QPOs. These results may suggest that either PKS 2247−131 is a rather special case because of its putative binary SMBH nature, or similar QPOs are hardly detectable with current high-energy facilities.

Methods
Multi-wavelength identification. PKS 2247−131 has an accurate sky position of RA = 22 h 49 m 59 s .6125, Dec. = −12°51′16′′.825 (equinox J2000.0), and has been detected by various radio surveys 46 and occasionally optical and infrared surveys, for example the SuperCOSMOS Science Archive 47 and WISE all-sky data 48 . Due to its radio-emission properties and the possibility of being a counterpart to the EGRET γ-ray source 3EG J2251−1341 49 , PKS 2247−131 has been considered to be an AGN.
For the outburst starting from 2016 October, 16 Swift X-ray observations of the field containing PKS 2247−131 were conducted since 31 July 2016 (see Supplementary Table 2). We analysed all the Swift XRT data by using its standard data reduction tool xrtpipeline (version 0.13.4). An X-ray source was detected in 11 of the 16 observations. We used one observation (conducted on 2016 October 24) in which the source had the highest count rate, and determined its position to be RA = 22 h 49 m 59 s .8, Dec. = −12°51′16′′.6 (equinox J2000.0, with a positional uncertainty of 5′′.2 at a 90% confidence level). This X-ray position matches well with those of the radio and optical counterparts (see Supplementary Figure 6).
Using the 6.5-meter Clay Telescope, we took two spectroscopic exposures of PKS 2247−131 on 2017 November 22. The exposure times were 20 min and 15 min. The instrument used was the Low Dispersion Survey Spectrograph 3 (LDSS3), with a wavelength coverage of 4250-9500 Å. The spectrum combined from the two exposures is basically featureless, having no emission lines but a few weak absorption lines (~3-10σ detections; see Supplementary Figure 7). The source thus is a blazar because emission from a jet, arising from non-thermal radiation processes, dominates the host galaxy emission and only weak lines from the host galaxy may be present 50 . We determined the redshift z ≈ 0.22, which implies a source distance of~1100 Mega-parsec, where the cosmological parameters from the Planck mission 51 are used (the Hubble constant H 0 ≃ 67 km s −1 Mpc −1 ).
Fermi LAT data analysis. We used Fermi LAT data between MJD 54682 (4 August 2008) and MJD 58177 (28 February 2018) in the energy range of 0.1-300 GeV. Events with zenith angles ≤90°were selected, in order to avoid contamination from the Earth's limb. The region of interest (ROI) we set was 20°× 20°, centered at the position of PKS 2247−131. The standard binned likelihood analysis was performed on the data in the ROI. Sources in the Fermi LAT Third Source Catalog 52 and the Galactic and isotropic background (files gll_iem_v06.fits and iso_P8R2_SOURCE_V6_v06.txt, respectively) were included in our source model. PKS 2247−131 was not detected by Fermi LAT before April 2016, resulting in a flux upper limit (at a 95% confidence level) of 3.65×10 −9 ph cm −2 s −1 . Starting from 12 April 2016, the source could be detected in most sets of 5-day binned data. We thus constructed its light curve binned in 5-day intervals, in which the spectral parameters of the sources in the ROI were fixed at their best-fit values, except that the normalizations of the two background components and one variable source within 5°(Variability_Index>72.44) 52 were set as free parameters. Detailed flux variations were clearly revealed in the light curve. This choice of 5 day binning provided the shortest intervals for which nearly all bins were long enough to yield detections during the flare. In a few binned data that did not detect PKS 2247−131 (the Test Statistic (TS) 53 values are smaller than 9), we calculated the 95% confidence flux upper limits. We also tried increasing the time bins to 10-40 days for the group of the upper-limit data points before the outburst peak (Fig. 1), but in two 20-day intervals, only upper limits were obtained. In addition, we also constructed a smooth light curve by shifting each 5-day time bin only one day forward (instead of 5 days) and obtaining the flux in such time bins, which helps show fine details of the flux variations. This light curve is over-plotted in Fig. 1.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.