MULTIWAVELENGTH EVIDENCE FOR QUASI-PERIODIC MODULATION IN THE GAMMA-RAY BLAZAR PG 1553+113

Among active galactic nuclei (AGNs), blazars are distinguished by erratic variability at all energies on a wide range of timescales. They are generally thought to be powered by supermassive black holes (SMBHs; 10⁸–10⁹ M_⊙). PG 1553+113 (1ES 1553+113, z ∼ 0.49; Danforth et al. 2010; Abramowski et al. 2015; Aliu et al. 2015) is an optically/X-ray selected BL Lac object (Falomo & Treves 1990) emitting variable GeV/TeV γ radiation (Abramowski et al. 2015; Aleksic et al. 2015). As is typical in very high energy (VHE) BL Lacs, the energetic non-thermal emission of PG 1553+113 originates in a relativistic jet and has a spectral energy distribution (SED) with two humps, overwhelming any other component from either the nucleus or the host galaxy.

The Large Area Telescope (LAT) on the Fermi Gamma-ray Space Telescope is providing continuous monitoring of the high-energy γ-ray sky. The apparent modulation noted in the γ-ray flux of PG 1553+113 initiated the multifrequency and long-term variability study described in this paper.

In Section 2 we describe the Fermi-LAT data analysis and the sources of multiwavelength data; Section 3 details the multiple approaches used for light curves and the cross-correlation analysis; Section 4 outlines preliminary scenarios to interpret these results.

The LAT is a pair conversion telescope with a 2.4 sr field of view, sensitive to γ rays from ∼20 MeV to >300 GeV (Atwood et al. 2009). The present work uses the new Pass 8 LAT database (Atwood et al. 2013). The LAT operating mode allows it to cover the entire sky every two ∼1.6 hr spacecraft orbits, providing a regular and uniform view of γ-ray sources and sampling timescales from hours to years. This work uses observations of PG 1553+113 covering ∼6.9 years (2008 August 4 to 2015 July 19; Modified Julian Day, MJD, 54682.65–57222.65). The LAT data analysis employed the standard ScienceTools v10r0p5⁶⁶ package, selecting events from 100 MeV to 300 GeV with P8R2_SOURCE_V6 instrument response functions, in a circular region of interest of 10^◦ radius centered on the position of PG 1553+113. It used files gll_iem_v06 and iso_P8R2_SOURCE_V6_v06 to model the Galactic and isotropic diffuse emission. Contamination due to the γ-ray-bright Earth limb is avoided by excluding events with zenith angle >90°. An unbinned maximum likelihood model fit technique is applied to each time bin with a power-law spectral model and photon index fixed to the 3FGL Catalog average value (1.604 ± 0.025; Acero et al. 2015) for PG 1553+113. The resulting light curves are shown in Figure 1.

Zoom In Zoom Out Reset image size

Optical R-band data covering an interval of ∼9.9 years (2005 April 19 to 2015 March 29; MJD 53479-57110) are reported in Figure 2. Most unpublished observations were performed as part of the Tuorla blazar monitoring program (Takalo et al. 2008).⁶⁷ The data are reduced using a semi-automatic pipeline (K. Nilsson et al. 2015, in preparation). Public data from the Katzman Automatic Imaging Telescope (KAIT) and the Catalina Sky Survey (CSS) programs are also added. V-band magnitudes are scaled to the R-band values.

Zoom In Zoom Out Reset image size

As part of an ongoing blazar monitoring program supporting Fermi (Richards et al. 2011), the Owens Valley Radio Observatory (OVRO) 40 m radio telescope has been observing PG 1553+113 continually (about every 1–23 days) since 2008 August. Figure 2 reports published 15 GHz light curves for the period from 2008 August 19 to 2014 May 18 (MJD 54697-56795). OVRO instrumentation, data calibration, and reduction are described in Richards et al. (2011).

Swift observed PG 1553+113 110 times between 2005 April 20 and 2015 July 18 (unabsorbed 0.3–2 keV flux light curve in Figure 2). X-Ray Telescope (XRT) data were first calibrated and cleaned (xrtpipeline, XRTDAS v.3.0.0) and then energy spectra were extracted from a region of 20 pixel (∼47'') radius, with a nearby 20 pixel radius region as a background. Individual XRT spectra are well fitted with a log-parabolic model, with column density fixed to the Galactic value of 3.6 × 10²⁰ cm⁻² (Kalberla et al. 2005). Aperture photometry (5'' radius) for the UVOT V-band filter was performed.

We performed continuous wavelet transform (CWT) and Lomb–Scargle Periodogram (LSP) analyses on the light curves. Figure 3 shows clear peaks at ∼2 years for γ-ray and optical power spectra. We also made an epoch folding (pulse shape) analysis to extract the period, shape, amplitude, and phase with uncertainties (Larsson 1996). The χ² for the folded pulse as a function of the trial periods was fitted with a model containing four Fourier components, giving a period of 798 ± 30 days (2.18 ± 0.08 years), consistent with the CWT and LSP findings (Figure 3). The value of the signal power peak does not change using regular 20 day and 45 day bins or an adaptive-bin technique (Lott et al. 2012) for construction of the LAT light curve.

Zoom In Zoom Out Reset image size

A direct power density spectrum (PDS) constructed from a LAT count rate light curve using exposure-weighted aperture photometry (Corbet et al. 2007; Kerr 2011) above 100 MeV for a region with 3° radius with 600 s time bins (Figure 4) confirms previous results with a peak at 2.16 ± 0.08 years, at 82× the mean power level. The low-frequency modulation prevents an easy fit subtraction to the PDS continuum. The peak is ∼5 times the mean level using a 4th order polynomial fit.

Zoom In Zoom Out Reset image size

The significance of any apparent periodic variation depends on what assumption is made about spurious stochastic variability mimicking a periodic variation. The significance of the ∼2 year γ-ray periodicity is difficult to assess given the limited length of the γ-ray light curve. Red noise, i.e., random and relatively enhanced low-frequency fluctuations over intervals comparable to the sample length, hinders the evaluation of the periodicity significance (e.g., Hsieh et al. 2005; Lasky et al. 2015). We have approached the problem with two procedures.

(1) The red noise is assumed to be produced by similar amplitude flares (as seen in PG 1553+113 and some other LAT blazars), and the probability for these to line up in a regular pattern is estimated. The coherence of the periodic modulation was investigated by studying phase variations along the light curve. The local phase at each minimum and maximum was estimated by correlating a one-period-long data segment with the Fourier template of the full light curve. The rms variations relative to a perfectly coherent modulation were 27.4 days. The chance probabilities for three, four, and five random events to be distributed with at least this coherence, as estimated by Monte Carlo simulations, are 0.0535, 0.0105, and 0.0027, respectively, implying a chance probability of a few percent for the 3.5 peak γ-ray light curve of PG 1553+113.

(2) We modeled the red noise using Monte Carlo simulations with a first-order autoregressive process as the null hypothesis to assess whether the signal is consistent with a stochastic origin. The nonlinear influence on the PDS is minimal thanks to the evenly spaced γ-ray light curve. The power peak in Figure 3 is above the 99% confidence contour level, i.e., has <1% chance of being a statistical fluctuation. The optical power peak has <5% chance of being a statistical fluctuation.

Although the γ-ray periodicity signal alone is not compelling, the 9.9 years of optical data support the finding of a periodic oscillation in PG 1553+113. The optical data, although affected by seasonal gaps, were analyzed using the same techniques as for the γ-ray data. This analysis gives a period of 754 ± 20 days (2.06 ± 0.05 years), consistent within uncertainties with the γ-ray results (Figure 3).

The less coherent 15 GHz light curve (5.7 years OVRO data) shows a signal power peak at 1.9 ± 0.1 year, with an additional power component at a 1.2 year timescale. Swift-XRT data show a factor of 5 variation linearly correlated with the γ-ray flux, while the synchrotron peak frequency shows a factor of ∼6 increase during high X-ray states, as suggested by Reimer et al. (2008).

The long-term X-ray count rate light curve from the RXTE ASM instrument (1996 February 20 to 2010 September 11) and the Swift-BAT (from 2005 May 29) were also analyzed but do not show any signal above the low-frequency noise because of insufficient statistics.

An important diagnostic for multifrequency periodicity analysis is the discrete cross-correlation function (DCCF) used with two independent and complementary approaches.

In the first procedure, flux variations are modeled assuming a simple power law ∝ 1/f^α (with f = 1/t) in the PDS as measured directly from the light curve data, allowing us to estimate the cross-correlation significance, and avoiding the assumption of equal variability in all sources at the cost of a model assumption (Max-Moerbeck et al. 2014). For the γ-ray light curve with 20 day binning we obtain a best fit α = 0.8, but the error is unconstrained, indicating that the length of the data set is too short (i.e., below five cycles) relative to the suspected periodic modulation to enable a reliable data characterization. The 45 day bin light curve yields a best fit α = 0.1 with unconstrained error. The optical PSD is constrained: the best fit value is α = 1.85, with 1σ limits at [1.75, 2.00]. The 15 GHz flux light curve has a slope of α = 1.4, with unconstrained limits on the α values as for the γ-ray data. The DCCF between the unbinned radio light curve and the 20 day bin γ-ray light curve results in a most probable time lag for the radio flux lagging the γ-ray flux by 50 ± 20 days, with 98.14% significance for the best PSD fit with a range of [89.56%–99.99%] when fit errors are taken into account (Figure 5) using the fitting procedure of Max-Moerbeck et al. (2014). The DCCF between the unbinned optical light curve and the 20 day bin γ-ray light curve results in a most probable time lag for the γ-ray flux lagging the optical flux by 130 ± 14 days, with 99.14% significance for the best PSD fit and [96.09%–99.97%] when fit errors are taken into account (Figure 5). The DCCF peak is broad, however, and consistent with no lag. This is also seen when the optical data are rebinned into 20 day intervals, as shown in the bottom panel, where the most probable lag is 10 ± 51 days.

Zoom In Zoom Out Reset image size

In the second procedure, the significance of the γ-ray–radio correlation was estimated to be 95% using a mixed source correlation procedure (Fuhrmann et al. 2014), cross-correlating the PG 1553+113 light curve with those of 132 comparison sources in that work, and evaluating the average DCCF level for time lags of −100 to +100 days. The γ-ray–optical correlation is significant at the 99% level, even though it is partly limited by the number of comparison sources and optical light curve gaps. With only 132 comparison light curves we can measure a minimum probability value of 0.0075, therefore in principle a 99% level of significance, but in this approach the error in that estimate is hard to determine. With the mixed source methods there are two limitations: (1) the assumption that all sources can be described with the same model for the variability, and (2) the sample variance due to the limited number of light curves must be assessed. The optical flux is found to lead the γ-ray variations by 75 ± 27 days and the radio by 158 ± 10 days (γ-ray variations lead the 15 GHz flux variations by 83 ± 27 days). The possible reverse γ-ray–optical time lag decreases to 28 ± 27 days when the optical light curve is binned.

The possible optical–γ-ray lag was already pointed out by Cohen et al. (2014) using KAIT unbinned optical light curves and LAT data. The high degree of γ-ray–radio correlation in PG 1553+113 is not typically found in other individual blazars/AGNs (see Max-Moerbeck et al. 2014). Significant cross-correlations are, nevertheless, found when stacking blazar samples (radio lagging γ-rays; Fuhrmann et al. 2014).

Factors that led to the indication of a possible ∼2 year periodic modulation in PG 1553+113 are the continuous all-sky survey of Fermi, the increased capability of the new Fermi-LAT Pass 8 data, and the long-term radio/optical monitoring of γ-ray blazars. Although the statistical significance of periodicity is marginal in each band, the consistent positive cross-correlation between bands strengthens the case, making PG 1553+113 the first possible quasi-periodic GeV γ-ray blazar and a prime candidate for further studies. Hints of possible γ-ray periodicities are rare in the literature (for example Sandrinelli et al. 2014). The similarity of the low- and high-energy modulation in PG 1553+113 is also a novel behavior for AGNs (Rieger 2004, 2007). Any periodic driving scenario should be related to the relativistic jet itself or to the process feeding the jet for this VHE BL Lac object. We outline, as examples, four possibilities.

• 1.  
  Pulsational accretion flow instabilities, approximating periodic behavior, are able to explain modulations in the energy outflow efficiency. Magnetically arrested and magnetically dominated accretion flows (MDAFs) could be suitable regimes for radiatively inefficient BL Lacs (Fragile & Meier 2009), characterized by advection-dominated accretion flows and subluminal, turbulent, and peculiar radio kinematics (Kharb et al. 2008; Karouzos et al. 2012; Piner & Edwards 2014). Such kinematics are sometimes explained as a precessing or helical jet (Conway & Murphy 1993). MDAFs in an inner disk portion are able to efficiently impart energy to particles in the jets of VHE BL Lacs (Tchekhovskoy et al. 2011). Periodic instabilities are believed to have short periods, ∼10⁵ s · (M_SMBH/10⁸ M_⊙) (Honma et al. 1992), but MHD simulations of magnetically choked accretion flows are seen to produce longer periods for slow-spinning SMBHs (McKinney et al. 2012).
• 2.  
  Jet precession (e.g., Romero et al. 2000; Stirling et al. 2003; Caproni et al. 2013), rotation (Camenzind & Krockenberger 1992; Vlahakis & Tsinganos 1998; Hardee & Rosen 1999), or helical structure (e.g., Conway & Murphy 1993; Roland et al. 1994; Villata & Raiteri 1999; Nakamura & Meier 2004; Ostorero et al. 2004; Mohan & Mangalam 2015), i.e., geometrical models (Rieger 2004), in the presence of a jet wrapped by a sufficiently strong magnetic field, could have a net apparent periodicity from the change of the viewing angle. Correspondingly, the resulting Doppler magnification factor changes periodically without the need for intrinsic variation in outflows and efficiency. Non-ballistic hydrodynamical jet precession may explain variations with periods >1 year (Rieger 2004). A differential Doppler factor ${\rm{\Delta }}{\mathcal{D}}(t)={{\rm{\Gamma }}}^{-1}{(1-\beta (t)\mathrm{cos}\theta (t))}^{-1}\lesssim 40\%$ variation (precession angle ∼1°) might be sufficient to support the ∼2.8 amplitude flux modulation seen in γ rays. A homogeneous curved helical jet scenario for PG 1553+113 was proposed in Raiteri et al. (2015).
• 3.  
  A mechanism analogous to low-frequency QPO from Galactic high-mass binaries/microquasars could produce an accretion–outflow coupling mechanism as the basis of the periodicity (Fender & Belloni 2004). King et al. (2013) ascribed the radio QPO in the FSRQ CGRaBS J1359+4011 to this mechanism. However BL Lac objects like PG 1553+113 are thought to possess a lower accretion rate. The microquasar QPO mechanism of Lense–Thirring precession (Wilkins 1972) requires that the inner accretion flow forms a geometrically thick torus rather than a standard thin disk as the latter warps (Bardeen–Petterson effect, Bardeen & Petterson 1975) rather than precesses (Ingram et al. 2009). A low mass accretion rate means that the accretion process probably forms an advection-dominated accretion flow (ADAF) so it can precess (Fragile & Meier 2009). The X-ray emission in PG 1553+113 is probably from the jet rather than from the flow, making it unlikely that the changing inclination of the hot flow causes the QPO. However, Lense–Thirring precession of the flow could affect the jet direction, giving the QPO as in (2) above.
• 4.  
  The presence of a gravitationally bound binary SMBH system (Begelman et al. 1980; Barnes & Hernquist 1992) with a total mass of ∼10⁸ M_⊙ and a milliparsec separation in the early inspiral gravitational-wave driven regime might be another hypothesis. Keplerian binary orbital motion would induce periodic accretion perturbations (Valtonen et al. 2008; Pihajoki et al. 2013; Liu et al. 2015) or jet nutation expected from the misalignment of the rotating SMBH spins or the gravitational torque on the disk exerted by the companion (Katz 1997; Romero et al. 2000; Caproni et al. 2013; Graham et al. 2015). Significant acceleration of the disk evolution and accretion onto a binary SMBH system is depicted by modeling (Nixon et al. 2013; Doğan et al. 2015).Binary SMBH induced periodicities have timescales ranging from ∼1 to ∼25 years (Komossa 2006; Rieger 2007). The SMBH total mass in PG 1553+113, estimated utilizing the putative link between inflow/accretion (disk luminosity) and outflow/jet (jet power) in blazars (Ghisellini et al. 2014), is ≃1.6 × 10⁸ ${M}_{\odot }$ using a 0.1 ${\dot{M}}_{\mathrm{Edd}}$ rate and Doppler factor ${\mathcal{D}}=30,$ in agreement with estimates for VHE BL Lacs (Woo et al. 2005). The observed 2.18 year period is equivalent to an intrinsic orbital time ${T}_{\mathrm{Kep}}^{\prime }\leqslant {T}_{\mathrm{obs}}/(1+z)\simeq 1.5$ years, and the binary system size would be 0.005 pc (∼100 Schwarzschild radii). The probability of observing such a milliparsec system, estimated from the binary mass ratios ∼0.1–0.01 and the GW-driven regime lifetime (Peters 1964), t_GW ≃ 10⁵−10⁶ years, might be too small.

Periodicities claimed for AGNs are often controversial; however, PG 1553+113 may potentially represent a key γ-ray/multimessenger laboratory in the hypothesis of low-frequency gravitational-wave emission and may have associated PeV neutrino emission (Padovani & Resconi 2014). Very long baseline interferometry structure observations (Kharb et al. 2008), radio/optical polarization data, and a prolonged multifrequency monitoring campaign will shed light on the situation. If the periodic modulation is real and coherent, as would be expected for a binary scenario, then subsequent maxima would be expected in 2017 and 2019, well within the possible lifetime of the Fermi mission.

We thank the anonymous referee for useful and constructive comments. We extend special thanks to Prof. C. Done of Durham University, UK, and Prof. R. W. Romani of Stanford University, USA, for useful comments during the course of this work. The Fermi-LAT Collaboration acknowledges support for LAT development, operation and data analysis from NASA and DOE (United States); CEA/Irfu and IN2P3/CNRS (France); ASI and INFN (Italy); MEXT, KEK, and JAXA (Japan); and the K.A. Wallenberg Foundation, the Swedish Research Council, and the National Space Board (Sweden). Science analysis support in the operations phase from INAF (Italy) and CNES (France) is also gratefully acknowledged. The Tuorla blazar monitoring program has been partially supported by the Academy of Finland grant 127740. The KAIT telescope program is supported by Katzman Foundation and the National Science Foundation. The CSS survey is funded by the National Aeronautics and Space Administration under grant No. NNG05GF22G issued through the Science Mission Directorate Near-Earth Objects Observations Program. The CRTS survey is supported by the U.S. National Science Foundation under grants AST-0909182. The OVRO 40 m program is supported in part by NASA grants NNX08AW31G and NNX11A043G and NSF grants AST-0808050 and AST-1109911. The MOJAVE program is supported under NASA-Fermi grant NNX12A087G. The National Radio Astronomy Observatory (NRAO) is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The NASA Swift γ-ray burst explorer is a MIDEX Gamma Ray Burst mission led by NASA with participation of Italy and the UK. This research has made use of the Smithsonian/NASA's ADS bibliographic database. This research has made use of the NASA/IPAC NED database (JPL CalTech and NASA, USA). This research has made use of the archives and services of the ASI Science Data Center (ASDC), a facility of the Italian Space Agency (ASI Headquarters, Rome, Italy). This research has made use of the XRT Data Analysis Software (XRTDAS) developed under the responsibility of the ASDC. This work is a product of the ASDC Fermi team developed in the frame of the INAF Senior Scientists project and the foreign visiting scientists program of ASDC.

Facilities: Fermi Gamma-ray Space Telescope, Swift, OVRO:40m, KAIT.