Multiwavelength Evidence for Quasi-periodic Modulation in the Gamma-ray Blazar PG 1553+113
Fermi LAT collaboration, M. Ackermann, M. Ajello, A. Albert, W. B. Atwood, L. Baldini, J. Ballet, G. Barbiellini, D. Bastieri, J. Becerra Gonzalez, R. Bellazzini, E. Bissaldi, R. D. Blandford, E. D. Bloom, R. Bonino, E. Bottacini, J. Bregeon, P. Bruel, R. Buehler, S. Buson, G. A. Caliandro, R. A. Cameron, R. Caputo, M. Caragiulo, P. A. Caraveo, E. Cavazzuti, C. Cecchi, A. Chekhtman, J. Chiang, G. Chiaro, S. Ciprini, J. Cohen-Tanugi, J. Conrad, S. Cutini, F. D'Ammando, A. De Angelis, F. De Palma, R. Desiante, L. Di Venere, A. Dominguez, P. S. Drell, C. Favuzzi, S. J. Fegan, E. C. Ferrara, W. B. Focke, L. Fuhrmann, Y. Fukazawa, P. Fusco, F. Gargano, D. Gasparrini, N. Giglietto, P. Giommi, F. Giordano, M. Giroletti, G. Godfrey, D. Green, I. A. Grenier, J. E. Grove, S. Guiriec, A. K. Harding, E. Hays, J.W. Hewitt, A. B. Hill, D. Horan, T. Jogler, G. Johannesson, A. S. Johnson, T. Kamae, M. Kuss, S. Larsson, L. Latronico, J. Li, L. Li, F. Longo, F. Loparco, B. Lott, M. N. Lovellette, P. Lubrano, J. Magill, S. Maldera, A. Manfreda, W. Max-Moerbeck, M. Mayer, M. N. Mazziotta, J. E. Mcenery, P. F. Michelson, T. Mizuno, M. E. Monzani, A. Morselli, I. V. Moskalenko, S. Murgia, E. Nuss, M. Ohno, T. Ohsugi, R. Ojha, N. Omodei, E. Orlando, J. F. Ormes, D. Paneque, T. J. Pearson, et al. (42 additional authors not shown)
aa r X i v : . [ a s t r o - ph . H E ] O c t A CCEPTED TO T HE A STROPHYSICAL J OURNAL L ETTERS
Preprint typeset using L A TEX style emulateapj v. 5/2/11
MULTIWAVELENGTH EVIDENCE FOR QUASI-PERIODIC MODULATIONIN THE GAMMA-RAY BLAZAR PG 1553+113
M. A
CKERMANN , M. A JELLO , A. A LBERT , W. B. A TWOOD , L. B ALDINI , J. B
ALLET , G. B ARBIELLINI , D. B
ASTIERI ,J. B
ECERRA G ONZALEZ , R. B
ELLAZZINI , E. B ISSALDI , R. D. B LANDFORD , E. D. B LOOM , R. B ONINO , E. B
OTTACINI ,J. B REGEON , P. B RUEL , R. B UEHLER , S. B USON , G. A. C
ALIANDRO , R. A. C
AMERON , R. C APUTO , M. C ARAGIULO ,P. A. C ARAVEO , E. C AVAZZUTI , C. C ECCHI , A. C
HEKHTMAN , J. C HIANG , G. C HIARO , S. C IPRINI ,J. C
OHEN -T ANUGI , J. C ONRAD , S. C
UTINI , F. D’A
MMANDO , A. DE A NGELIS , F. DE P ALMA ,R. D
ESIANTE , L. D I V ENERE , A. D OM ´ INGUEZ , P. S. D RELL , C. F AVUZZI , S. J. F
EGAN , E. C. F ERRARA ,W. B. F OCKE , L. F UHRMANN , Y. F UKAZAWA , P. F USCO , F. G
ARGANO , D. G ASPARRINI , N. G
IGLIETTO ,P. G
IOMMI , F. G IORDANO , M. G
IROLETTI , G. G ODFREY , D. G REEN , I. A. G
RENIER , J. E. G ROVE , S. G UIRIEC ,A. K. H
ARDING , E. H AYS , J.W. H EWITT , A. B. H ILL , D. H
ORAN , T. J OGLER , G. J ´ OHANNESSON , A. S. J OHNSON ,T. K AMAE , M. K USS , S. L ARSSON , L. L
ATRONICO , J. L I , L. L I , F. L ONGO , F. L
OPARCO , B. L
OTT ,M. N. L OVELLETTE , P. L UBRANO , J. M
AGILL , S. M ALDERA , A. M ANFREDA , W. M AX -M OERBECK , M. M
AYER ,M. N. M AZZIOTTA , J. E. M C E NERY , P. F. M
ICHELSON , T. M IZUNO , M. E. M ONZANI , A. M ORSELLI , I. V. M OSKALENKO ,S. M URGIA , E. N USS , M. O HNO , T. O HSUGI , R. O JHA , N. O MODEI , E. O RLANDO , J. F. O RMES , D. P ANEQUE ,T. J. P
EARSON , J. S. P ERKINS , M. P ERRI , M. P ESCE -R OLLINS , V. P
ETROSIAN , F. P IRON , G. P IVATO , T. A. P ORTER ,S. R AIN ` O , R. R ANDO , M. R
AZZANO , A. R
EADHEAD , A. R EIMER , O. R
EIMER , A. S
CHULZ , C. S GR ` O ,E. J. S ISKIND , F. S PADA , G. S PANDRE , P. S PINELLI , D. J. S
USON , H. T AKAHASHI , J. B. T HAYER ,D. J. T HOMPSON , L. T
IBALDO , D. F. T ORRES , G. T
OSTI , E. T
ROJA , Y. U
CHIYAMA , G. V IANELLO ,K. S. W OOD , M. W OOD , S. Z IMMER , A. B
ERDYUGIN , R. H. D. C ORBET , T. H
OVATTA , E. L INDFORS , K. N ILSSON ,R. R EINTHAL , A. S ILLANP ¨ A ¨ A , A. S TAMERRA , L. O. T
AKALO , M. J. V ALTONEN Accepted to The Astrophysical Journal Letters
ABSTRACTWe report for the first time a γ -ray and multiwavelength nearly-periodic oscillation in an active galacticnucleus. Using the Fermi
Large Area Telescope (LAT) we have discovered an apparent quasi-periodicity inthe γ -ray flux ( E >
MeV) from the GeV/TeV BL Lac object PG 1553+113. The marginal significance ofthe . ± . year-period γ -ray cycle is strengthened by correlated oscillations observed in radio and opticalfluxes, through data collected in the OVRO, Tuorla, KAIT, and CSS monitoring programs and Swift
UVOT. Theoptical cycle appearing in ∼ years of data has a similar period, while the 15 GHz oscillation is less regularthan seen in the other bands. Further long-term multi-wavelength monitoring of this blazar may discriminateamong the possible explanations for this quasi-periodicity. Subject headings: gamma rays: galaxies — gamma rays: general — BL Lacertae objects: general — BLLacertae objects: individual (PG 1553+113) — galaxies: jets — accretion, accretion disks Deutsches Elektronen Synchrotron DESY, D-15738 Zeuthen, Ger-many Department of Physics and Astronomy, Clemson University, KinardLab of Physics, Clemson, SC 29634-0978, USA W. W. Hansen Experimental Physics Laboratory, Kavli Institutefor Particle Astrophysics and Cosmology, Department of Physics andSLAC National Accelerator Laboratory, Stanford University, Stanford, CA94305, USA Santa Cruz Institute for Particle Physics, Department of Physics andDepartment of Astronomy and Astrophysics, University of California atSanta Cruz, Santa Cruz, CA 95064, USA Universit`a di Pisa and Istituto Nazionale di Fisica Nucleare, Sezionedi Pisa I-56127 Pisa, Italy Laboratoire AIM, CEA-IRFU/CNRS/Universit´e Paris Diderot, Ser-vice d’Astrophysique, CEA Saclay, F-91191 Gif sur Yvette, France Istituto Nazionale di Fisica Nucleare, Sezione di Trieste, I-34127 Tri-este, Italy Dipartimento di Fisica, Universit`a di Trieste, I-34127 Trieste, Italy Istituto Nazionale di Fisica Nucleare, Sezione di Padova, I-35131Padova, Italy Dipartimento di Fisica e Astronomia “G. Galilei”, Universit`a diPadova, I-35131 Padova, Italy NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA Department of Physics and Department of Astronomy, University ofMaryland, College Park, MD 20742, USA Istituto Nazionale di Fisica Nucleare, Sezione di Pisa, I-56127 Pisa,Italy Istituto Nazionale di Fisica Nucleare, Sezione di Bari, I-70126 Bari, Italy Istituto Nazionale di Fisica Nucleare, Sezione di Torino, I-10125Torino, Italy Dipartimento di Fisica Generale “Amadeo Avogadro” , Universit`adegli Studi di Torino, I-10125 Torino, Italy Laboratoire Univers et Particules de Montpellier, Universit´e Mont-pellier, CNRS/IN2P3, Montpellier, France Laboratoire Leprince-Ringuet, ´Ecole polytechnique, CNRS/IN2P3,Palaiseau, France Consorzio Interuniversitario per la Fisica Spaziale (CIFS), I-10133Torino, Italy INAF-Istituto di Astrofisica Spaziale e Fisica Cosmica, I-20133 Mi-lano, Italy Agenzia Spaziale Italiana (ASI) Science Data Center, I-00133 Roma,Italy Istituto Nazionale di Fisica Nucleare, Sezione di Perugia, I-06123 Pe-rugia, Italy Dipartimento di Fisica, Universit`a degli Studi di Perugia, I-06123 Pe-rugia, Italy College of Science, George Mason University, Fairfax, VA 22030,resident at Naval Research Laboratory, Washington, DC 20375, USA INAF Osservatorio Astronomico di Roma, I-00040 Monte PorzioCatone (Roma), Italy email: [email protected] Department of Physics, Stockholm University, AlbaNova, SE-106 91Stockholm, Sweden The Oskar Klein Centre for Cosmoparticle Physics, AlbaNova, SE-106 91 Stockholm, Sweden
Ackermann et al. INTRODUCTION
Among active galactic nuclei (AGN), blazars are distin-guished by erratic variability at all energies on a wide rangeof timescales. They are generally thought to be pow-ered by supermassive black holes (SMBHs, 10 –10 M ⊙ ).PG 1553+113 (1ES 1553+113, z ∼ . , Danforth etal. 2010; Aliu et al. 2015; Abramowski et al. 2015) is anoptically/X-ray selected BL Lac object (Falomo & Treves The Royal Swedish Academy of Sciences, Box 50005, SE-104 05Stockholm, Sweden email: [email protected] INAF Istituto di Radioastronomia, I-40129 Bologna, Italy Dipartimento di Astronomia, Universit`a di Bologna, I-40127Bologna, Italy Dipartimento di Fisica, Universit`a di Udine and Istituto Nazionale diFisica Nucleare, Sezione di Trieste, Gruppo Collegato di Udine, I-33100Udine Universit`a Telematica Pegaso, Piazza Trieste e Trento, 48, I-80132Napoli, Italy Universit`a di Udine, I-33100 Udine, Italy Dipartimento di Fisica “M. Merlin” dell’Universit`a e del Politecnicodi Bari, I-70126 Bari, Italy Max-Planck-Institut f¨ur Radioastronomie, Auf dem H¨ugel 69, D-53121 Bonn, Germany Department of Physical Sciences, Hiroshima University, Higashi-Hiroshima, Hiroshima 739-8526, Japan Space Science Division, Naval Research Laboratory, Washington, DC20375-5352, USA NASA Postdoctoral Program Fellow, USA University of North Florida, Department of Physics, 1 UNF Drive,Jacksonville, FL 32224 , USA School of Physics and Astronomy, University of Southampton, High-field, Southampton, SO17 1BJ, UK Science Institute, University of Iceland, IS-107 Reykjavik, Iceland Department of Physics, Graduate School of Science, University ofTokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan Department of Physics, KTH Royal Institute of Technology, Al-baNova, SE-106 91 Stockholm, Sweden email: [email protected] Institute of Space Sciences (IEEC-CSIC), Campus UAB, E-08193Barcelona, Spain Centre d’ ´Etudes Nucl´eaires de Bordeaux Gradignan, IN2P3/CNRS,Universit´e Bordeaux 1, BP120, F-33175 Gradignan Cedex, France National Radio Astronomy Observatory (NRAO), Socorro, NM87801, USA Cahill Center for Astronomy and Astrophysics, California Institute ofTechnology, Pasadena, CA 91125, USA Hiroshima Astrophysical Science Center, Hiroshima University,Higashi-Hiroshima, Hiroshima 739-8526, Japan Istituto Nazionale di Fisica Nucleare, Sezione di Roma “Tor Vergata”,I-00133 Roma, Italy Center for Cosmology, Physics and Astronomy Department, Univer-sity of California, Irvine, CA 92697-2575, USA Department of Physics and Astronomy, University of Denver, Den-ver, CO 80208, USA Max-Planck-Institut f¨ur Physik, D-80805 M¨unchen, Germany Funded by contract FIRB-2012-RBFR12PM1F from the Italian Min-istry of Education, University and Research (MIUR) Institut f¨ur Astro- und Teilchenphysik and Institut f¨ur TheoretischePhysik, Leopold-Franzens-Universit¨at Innsbruck, A-6020 Innsbruck, Aus-tria NYCB Real-Time Computing Inc., Lattingtown, NY 11560-1025,USA Department of Chemistry and Physics, Purdue University Calumet,Hammond, IN 46323-2094, USA email: [email protected] Instituci´o Catalana de Recerca i Estudis Avanc¸ats (ICREA),Barcelona, Spain Tuorla Observatory, University of Turku, FI-21500 Piikki¨o, Finland Center for Research and Exploration in Space Science and Technol-ogy (CRESST) and NASA Goddard Space Flight Center, Greenbelt, MD20771, USA γ radiation (Aleksic et al.2015; Abramowski et al. 2015). As typical in very-high en-ergy (VHE) BL Lacs, the energetic non-thermal emission ofPG 1553+113 originates in a relativistic jet and has a spec-tral energy distribution (SED) with two humps, overwhelm-ing any other component from either the nucleus or the hostgalaxy.The Large Area Telescope (LAT) on the Fermi Gamma-raySpace Telescope is providing continuous monitoring of thehigh-energy γ -ray sky. The apparent modulation noted in the γ -ray flux of PG 1553+113 stimulated the multi-frequencyand long-term variability study described in this paper.In § Fermi
LAT data analysis and thesources of multiwavelength data; § § FERMI LAT AND RADIO, OPTICAL, X-RAY DATA
The LAT is a pair conversion detector with a 2.4 sr field ofview, sensitive to γ rays from ∼ MeV to > GeV (At-wood et al. 2009). The present work uses the new Pass 8 LATdatabase (Atwood et al. 2013). The LAT operating mode al-lows it to cover the entire sky every two ∼ γ -ray sources,sampling timescales from hours to years. This work uses ob-servations of PG 1553+113covering ∼ . years (2008 Au-gust 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 from100 MeV −
300 GeV with
P8R2 SOURCE V6 instrument re-sponse functions, in a circular Region of Interest of 10 ◦ ra-dius centered on the position of PG 1553+113. It used files gll iem v06 and iso P8R2 SOURCE V6 v06 to modelthe Galactic and isotropic diffuse emission. Contaminationdue to the γ -ray-bright Earth limb is avoided by excludingevents with zenith angle > ◦ . An unbinned maximum like-lihood model fit technique is applied to each time bin with apower-law spectral model and photon index fixed to the 3FGLCatalog average value ( . ± . , Acero et al. 2015) forPG 1553+113. The resulting lightcurves are shown in Fig. 1.Optical R-band data covering an interval of ∼ . years(2005 April 19 to 2015 March 29, MJD 53479-57110) arereported in Fig. 2. Most unpublished observations wereperformed as part of the Tuorla blazar monitoring program(Takalo et al. 2008) . The data are reduced using a semi-automatic pipeline (Nilsson et al. in prep.). Public data fromthe Katzman Automatic Imaging Telescope (KAIT) and theCatalina Sky Survey (CSS) programs are also added. V-bandmagnitudes are scaled to the R-band values. Department of Physics and Center for Space Sciences and Technol-ogy, University of Maryland Baltimore County, Baltimore, MD 21250,USA Aalto University, Mets¨ahovi Radio Observatory, Mets¨ahovintie 114,Kylmala, Finland Finnish Centre for Astronomy with ESO (FINCA), University ofTurku, FI-21500 Piikii¨o, Finland INAF, Osservatorio Astronomico di Torino, I-10025 Pino Torinese(TO), Italy Scuola Normale Superiore, Piazza dei Cavalieri, 7, I-56126 Pisa, Italy email: [email protected] * Correspondence: [email protected]; [email protected];[email protected]; [email protected]; [email protected] http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/ users.utu.fi/kani/1m uasi-periodic modulation in PG 1553+113 3 F IG . 1.— Fermi
LAT γ -ray lightcurves of PG 1553+113 over ∼ As part of an ongoing blazar monitoring program support-ing
Fermi (Richards et al. 2011), the Owens Valley RadioObservatory (OVRO) 40-m radio telescope has been observ-ing PG 1553+113 continually (about every 1 to 23 days) since2008 August. Figure 2 reports published 15 GHz lightcurvesfor the period from 2008 August 19 to 2014 May 18 (MJD54697-56795). OVRO instrumentation, data calibration andreduction are described in Richards et al. (2011).
Swift observed PG 1553+113 110 times between 2005 April20 and 2015 July 18 (unabsorbed 0.3–2 keV flux lightcurve inFigure 2). X-Ray Telescope (XRT) data were first calibratedand cleaned ( xrtpipeline , XRTDAS v.3.0.0) and energyspectra extracted from a region of 20 pixel ( ∼ ′′ ) radius,with a nearby 20 pixel radius region for background. Individ-ual XRT spectra are well fitted with a log-parabolic model,with column density fixed to the Galactic value of . × cm − (Kalberla et al. 2011). Aperture photometry (5 ′′ radius)for the UVOT V-band filter was performed. TEMPORAL VARIABILITY ANALYSIS AND CROSSCORRELATION ANALYSIS
We performed continuous wavelet transform (CWT)and Lomb-Scargle Periodogram (LSP) analyses on thelightcurves. Fig. 3 shows clear peaks at ∼ years for γ -ray and optical power spectra. We also made an epoch folding(pulse shape) analysis used to extract the period, shape, ampli-tude and phase, with uncertainties (Larsson 1996). The χ forthe folded pulse as a function of trial periods was fitted witha model containing 4 Fourier components, giving a period of ± days ( . ± . years), consistent with the CWTand LSP findings (Fig. 3). The value of the signal power peakdoes not change using regular 20-day and 45-day bins or anadaptive-bin technique (Lott et al. 2012) for construction ofthe LAT lightcurve.A direct Power Density Spectrum (PDS) constructed from aLAT count-rate lightcurve using exposure-weighted aperturephotometry (Corbet et al. 2007; Kerr 2011) above 100 MeVfor a region with 3 ◦ radius with 600 second time bins (Fig. Ackermann et al. F IG . 2.— Multifrequency lightcurves of PG 1553+113 at X-ray, optical and radio bands. Top panel: Swift -XRT integrated flux (0.3-2.0 keV). Central panel:optical flux density from Tuorla (R filter, black filled circle points), Catalina CSS (V filter rescaled, blue filled squared points), KAIT (V filter rescaled, red filleddiamond points),
Swift -UVOT (V filter rescaled, green filled circle points). Dotted line: LAT flux (
E >
MeV) with time bins of 20 days, scaled and y-shiftedfor comparison. Bottom panel: 15 GHz flux density from OVRO 40m (black filled circle points) and parsec-scale 15 GHz flux density by VLBA (MOJAVEprogram, yellow diamond filled points). . ± . years,at × the mean power level. The low-frequency modulationprevents an easy fit subtraction to the PDS continuum. Thepeak is ∼ times the mean level using a 4th order polynomialfit.The significance of any apparent periodic variation dependson what assumption is made about spurious stochastic vari-ability mimicking a periodic variation. The significance ofthe ∼ -year γ -ray periodicity is difficult to assess given thelimited length of the γ -ray lightcurve. Red-noise, i.e. randomand relatively enhanced low-frequency fluctuations over inter-vals comparable to the sample length, hinders the evaluationof 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 am-plitude flares (as seen in PG 1553+113 and some other LATblazars), and the probability for these to line up in a regularpattern is estimated. The coherence of the periodic modula-tion was investigated by studying phase variations along thelightcurve. The local phase at each minimum and maximum was estimated by correlating a one-period long data segmentwith the Fourier template of the full lightcurve. The rms vari-ations relative to a perfectly coherent modulation was 27.4days. The chance probabilities for 3, 4 and 5 random eventsto be distributed with at least this coherence, as estimated byMonte Carlo simulations, are 0.0535, 0.0105 and 0.0027 re-spectively, implying a chance probability of a few percent forthe 3.5-peak γ -ray lightcurve of PG 1553+113.2) We modeled the red-noise using Monte Carlo simu-lations with a first-order autoregressive process as the nullhypothesis to assess whether the signal is consistent with astochastic origin. Non-linear influence on the PDS is mini-mal thanks to the evenly spaced γ -ray lightcurve. The powerpeak in Fig. 3 is above the 99% confidence contour level, i.e.has < chance of being a statistical fluctuation. The opticalpower peak has < chance of being a statistical fluctuation.Although the γ -ray periodicity signal alone is not com-pelling, the 9.9-years of optical data support the finding ofa periodic oscillation in PG 1553+113. The optical data, al-though affected by seasonal gaps, were analyzed using theuasi-periodic modulation in PG 1553+113 5 F IG . 3.— Left top panel: pulse shape (epoch-folded) γ -ray ( E >
MeV) flux lightcurve at the 2.18 year period (two cycles shown). Left bottom panels:2D plane contour plot of the CWT power spectrum (scalogram) of the γ -ray lightcurve, using a Morlet mother function (filled color contour). The side panel tothis is the 1D smoothed, all-epoch averaged, spectrum of the CWT scalogram showing a signal power peak in agreement with the 2.18-year value, also showingthe LSP. Dashed lines depict increasing levels of confidence against red-noise calculated with Monte Carlo simulation. The γ -ray signal peak is above the 99%confidence contour level ( < % chance probability of being spurious). Right top panel: pulse shape from epoch folding of the optical flux lightcurve at the 2.18year period (two cycles shown). Right bottom panels: the same CWT and LSP diagrams for the optical lightcurve. The optical signal peak is above the 95%confidence contour level. same techniques as for the γ -ray data. This analysis gives aperiod of ± days ( . ± . years), consistent withinuncertainties with the γ -ray results (Fig. 3). −1 R e l a ti v e P o w e r F IG . 4.— Power Density Spectrum (PDS) of the LAT . − GeVcount rate lightcurve of PG 1553+113 from a 3 ◦ exposure-weighted aperturephotometry technique with 600-second time bins. The less coherent 15 GHz lightcurve (5.7-years OVROdata) shows a signal power peak at . ± . year, with an additional power component at a 1.2-year timescale. Swift
XRT data show a factor of 5 variation linearly correlated withthe γ -ray flux, while the synchrotron peak frequency shows afactor ∼ increase during high X-ray states, as suggested byReimer et al. (2008).The long-term X-ray count rate lightcurve from the Rossi -XTE ASM instrument (1996 February 20 to 2010 September11) and the
Swift -BAT (from 2005 May 29) were also ana-lyzed but do not show any signal above the low-frequencynoise, because of insufficient statistics.An important diagnostic for multi-frequency periodicityanalysis is the discrete cross-correlation function (DCCF)used with two independent and complementary approaches.In the first procedure, flux variations are modeled assum-ing a simple power law ∝ /f α (with f = 1 /t ) in the PDSas measured directly from the lightcurve data, allowing us toestimate the cross-correlations significance avoiding the as-sumption of equal variability in all sources at the cost of amodel assumption (Max-Moerbeck et al. 2014). For the γ -raylightcurve with 20-day binning we obtain a best fit α = 0 . ,but the error is unconstrained, indicating that the length ofthe data set is too short (i.e. below five cycles), relative tothe suspected periodic modulation, to enable a reliable datacharacterization. The 45-day bin lightcurve yields a best fit α = 0 . with unconstrained error. The optical PSD is con-strained: the best fit value is α = 1 . , with 1 σ limits at [1 . , . . The 15 GHz flux light curve a slope of α = 1 . ,with unconstrained limits on the α values as for the γ -ray data. Ackermann et al. −400 −200 0 200 400 τ [days] −1.0−0.50.00.51.0 D CC F −400 −200 0 200 400 τ [days] −1.0−0.50.00.51.0 D CC F −400 −200 0 200 400 τ [days] −1.0−0.50.00.51.0 D CC F F IG . 5.— Discrete cross-correlation plots from the approach with PDSmodel measured from the lightcurve data (Max-Moerbeck et al. 2014). Ineach plot the black dots are the DCCF estimates, and the red, orange andgreen lines are the 1 σ , 2 σ and 3 σ significance levels respectively. Top panel:DCCF between the radio 15GHz and γ -ray (20-day time bins) lightcurve.Central panel: DCCF between the unbinned optical lightcurve and γ -ray (20-day time bins) lightcurve. Bottom panel: DCCF between the 20-day rebinnedoptical lightcurve and γ -ray (20-day time bins) lightcurve. The oscillatingshape of the significance contours for this case is due to the number of sam-ples in each bin. The DCCF between the unbinned radio lightcurve and the 20-day bin γ -ray lightcurve results in a most probable time lagfor radio-flux lagging the γ -ray flux by ± days, witha 98.14% significance for the best PSD fit with a range of[89.56%-99.99%] when fit errors are taken into account (Fig.5), using the fitting procedure of Max-Moerbeck et al. (2014).The DCCF between the unbinned optical lightcurve and the20-day bin γ -ray lightcurve results in a most probable timelag for γ -ray flux lagging the optical flux by ± days,with a 99.14% significance for the best PSD fit and [96.09%-99.97%] when fit errors are taken into account (Fig. 5). TheDCCF 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 mostprobable lag is ± days.In the second procedure, the significance of the γ -ray– radio correlation was estimated to be 95% by a mixed source correlation procedure (Fuhrmann et al. 2014), cross-correlating the PG 1553+113 lightcurve with those of 132comparison sources in that work, and evaluating the averageDCCF level for time lags − to +100 days. The γ -ray –optical correlation is significant at the 99% level, even thoughpartly limited by the number of comparison sources and opti-cal lightcurve gaps. With only 132 comparison lightcurves wecan measure a minimum probability-value of 0.0075, there-fore in principle a 99% level of significance, but in this ap-proach the error in that estimate is hard to determine. Withthe mixed source methods there are two limitations: 1) the as-sumption that all the sources can be described with the samemodel for the variability, and 2) the sample variance due to thelimited number of lightcurves must be assessed. The opticalflux is found to lead the γ -ray variations by ± days andthe radio by ± days ( γ -ray variations lead the 15 GHzflux variations by ± days). The possible reverse γ -ray-optical time lag decreases to ± days when the opticallightcurve is binned.The possible optical- γ -ray lag was already pointed out byCohen et al. (2014), using KAIT unbinned optical lightcurvesand LAT data. The high degree of γ -ray-radio correlationin PG 1553+113 is not typically found in other individ-ual blazars/AGN (see Max-Moerbeck et al. 2014). Signifi-cant cross-correlations are, nevertheless, found when stackingblazar samples (radio lagging γ rays; Fuhrmann et al. 2014). DISCUSSION AND CONCLUSIONS
Factors that led to the indication of a possible ∼ -yearperiodic modulation in PG 1553+113 are: the continuousall-sky survey of Fermi ; the increased capability of the new
Fermi
LAT Pass 8 data; and the long-term radio/optical mon-itoring of γ -ray blazars. Although the statistical significanceof periodicity is marginal in each band, the consistent positivecross-correlation between bands strengthens the case, mak-ing PG 1553+113 the first possible quasi-periodic GeV γ -rayblazar and a prime candidate for further studies. Hints ofpossible γ -ray periodicities are rare in literature (for exampleSandrinelli et al. 2014). The similarity of the low- and high-energy modulation in PG 1553+113 is also a novel behaviorfor AGN (Rieger 2004, 2007). Any periodic driving scenarioshould be related to the relativistic jet itself or to the processfeeding the jet for this VHE BL Lac object. We outline, asexamples, four possibilities:1. Pulsational accretion flow instabilities, approximatingperiodic behavior, are able to explain modulations inthe energy outflow efficiency. Magnetically-arrestedand magnetically-dominated accretion flows (MDAF)could be suitable regimes for radiatively inefficientBL Lacs (Fragile & Meier 2009), characterized byadvection-dominated accretion flows and subluminal,turbulent and peculiar radio kinematics (Karouzos etal. 2012; Piner & Edwards 2014). Such kinematicsare sometimes explained as a precessing or helical jet(Conway & Murphy 1993). MDAF in a inner disc por-tion can be able to efficiently impart energy to parti-cles in the jets of VHE BL Lacs (Tchekhovskoy et al.2011). Periodic instabilities are believed to have shortperiods, ∼ s · ( M SMBH / M ⊙ ) (Honma et al.1992), but MHD simulations of magnetically chokedaccretion flows are seen to produce longer periods forslow-spinning SMBH (McKinney et al. 2012).2. Jet precession (e.g., Romero et al. 2000; Stirling etuasi-periodic modulation in PG 1553+113 7al. 2003; Caproni et al. 2013), rotation (Camenzind& Krockenberger 1992; Vlahakis & Tsinganos 1998;Hardee & Rosen 1999) or helical structure (e.g., Con-way & Murphy 1993; Roland et al. 1994; Villata &Raiteri 1999; Nakamura & Meier 2004; Ostorero et al.2004), i.e. geometrical models (Rieger 2004), in thepresence of a jet wrapped by a sufficiently strong mag-netic field, could have a net apparent periodicity fromthe change of the viewing angle. Correspondingly theresulting Doppler magnification factor changes periodi-cally without the need for intrinsic variation in outflowsand efficiency. Non-ballistic hydrodynamical jet pre-cession may explain variations with periods > year(Rieger 2004). A differential Doppler factor ∆ D ( t ) =Γ − (1 − β ( t ) cos θ ( t )) − . variation (precessionangle ∼ ◦ ) might be sufficient to support the ∼ . amplitude flux modulation seen in γ rays. A homoge-neous curved helical jet scenario for PG 1553+113 wasproposed in Raiteri et al. (2015).3. A mechanism analogous to low-frequency QPO fromGalactic high-mass binaries/microquasars could pro-duce an accretion-outflow coupling mechanism as thebasis of the periodicity (Fender & Belloni 2004). Kinget al. (2013) ascribed the radio QPO in the FSRQCGRaBS J1359+4011 to this mechanism. However BLLac objects like PG 1553+113 are thought to possessa lower accretion rate. The microquasar QPO mecha-nism of Lense-Thirring precession (Wilkins 1972) re-quires that the inner accretion flow forms a geometri-cally thick torus rather than a standard thin disc as thelatter warps (Bardeen-Petterson effect, Bardeen & Pet-terson 1975) rather than precesses (Ingram et al. 2009).A low mass accretion rate means that the accretion pro-cess probably forms an Advection-Dominated Accre-tion Flow (ADAF), so it can precess (Fragile & Meier2009). The X-ray emission in PG 1553+113 is prob-ably from the jet rather than from the flow, making itunlikely that the changing inclination of the hot flowcauses the QPO. However, Lense-Thirring precessionof the flow could affect the jet direction, giving the QPOas in (2) above.4. The presence of a gravitationally bound binary SMBHsystem (Begelman et al. 1980; Barnes & Hernquist1992) with a total mass ∼ M ⊙ , and a milli-pc sep-aration in the early inspiral gravitational-wave drivenregime, might be another hypothesis. Keplerian binaryorbital motion, would induce periodic accretion per-turbations (Valtonen et al. 2008; Pihajoki et al. 2013;Liu et al. 2015) or jet nutation expected from the mis-alignment of the rotating SMBH spins, or the gravi-tational torque on the disc exerted by the companion(Katz 1997; Romero et al. 2000; Caproni et al. 2013;Graham et al. 2015). Significant acceleration of the discevolution and accretion onto a binary SMBH system isdepicted by modeling (Nixon et al. 2013; Do˘gan et al.2015) .Binary SMBH induced periodicities have timescales ranging from ∼ to ∼ years (Komossa 2006;Rieger 2007). The SMBH total mass in PG 1553+113,estimated utilizing the putative link between in-flow/accretion (disc luminosity) and outflow/jet (jetpower) in blazars (Ghisellini et al. 2014), is ≃ . × M ⊙ , using a 0.1 ˙ M Edd rate and Doppler factor D = 30 ,in agreement with estimates for VHE BL Lacs (Woo etal. 2005).The observed 2.18-year period is equivalent to an in-trinsic orbital time T ′ Kep ≤ T obs / (1 + z ) ≃ . years,and the binary system size would be . pc ( ∼ Schwarzschild radii). The probability to be observingsuch milli-pc system, estimated from the binary massratios ∼ . − . and the GW-driven regime lifetime(Peters 1964), t GW ≃ − years might be toosmall.Periodicities claimed for AGN 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 asso-ciated PeV neutrino emission (Padovani & Resconi 2014).VLBI structure observations, radio/optical polarization data,and a prolonged multifrequency monitoring campaign willshed light on the situation. If the periodic modulation is realand coherent, as would be expected for a binary scenario, thensubsequent 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, andProf. R. W. Romani of Stanford University, USA, for useful comments dur-ing the course of this work. The
Fermi -LAT Collaboration acknowledgessupport for LAT development, operation and data analysis from NASA andDOE (United States), CEA/Irfu and IN2P3/CNRS (France), ASI and INFN(Italy), MEXT, KEK, and JAXA (Japan), and the K.A. Wallenberg Foun-dation, the Swedish Research Council and the National Space Board (Swe-den). Science analysis support in the operations phase from INAF (Italy) andCNES (France) is also gratefully acknowledged. Tuorla blazar monitoringprogram has been partially supported by Academy of Finland grant 127740.KAIT telescope program is supported by Katzman Foundation and the Na-tional Science Foundation. The CSS survey is funded by the National Aero-nautics and Space Administration under Grant No. NNG05GF22G issuedthrough the Science Mission Directorate Near-Earth Objects ObservationsProgram. The CRTS survey is supported by the U.S. National Science Foun-dation under grants AST-0909182. The OVRO 40-m program is supportedin part by NASA grants NNX08AW31G and NNX11A043G and NSF grantsAST-0808050 and AST-1109911. The MOJAVE program is supported un-der NASA-Fermi grant NNX12A087G. The National Radio Astronomy Ob-servatory (NRAO) is a facility of the National Science Foundation operatedunder cooperative agreement by Associated Universities, Inc. The NASA
Swift γ -ray burst explorer is a MIDEX Gamma Ray Burst mission led byNASA with participation of Italy and the UK. This research has made useof the Smithsonian/NASA’s ADS bibliographic database. This research hasmade 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 ScienceData Center (ASDC), a facility of the Italian Space Agency (ASI Headquar-ter, Rome, Italy). This research has made use of the XRT Data AnalysisSoftware (XRTDAS) developed under the responsibility of the ASDC. Thiswork is a product of the ASDC Fermi team developed in the frame of theINAF Senior Scientists project and the foreign visiting scientists program ofASDC.Facilities: Fermi Gamma-ray Space Telescope — Swift — OVRO — Tuorla — KVA — KAIT — CSS — CRTS
REFERENCESAbramowski, A., Aharonian, F., Ait Benkhali, F., et al. 2015, ApJ, 802, 65Acero, F., Ackermann, M., Ajello, M., et al. 2015, ApJS, 218, 23Aleksi´c, J., Ansoldi, S., Antonelli, L. A., et al. 2015, MNRAS, 450, 4399 Aliu, E., Archambault, S., Archer, A., et al. 2015, ApJ, 800, 61Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071