<figure><img src="image_1.png" /></figure>

Periodic Gamma-Ray Modulation of the Blazar PG 1553+113 Confirmed by Fermi-LAT and Multiwavelength Observations

S. Abdollahi¹, L. Baldini², G. Barbiellini³,⁴, R. Bellazzini⁵, B. Berenji⁶, E. Bissaldi⁷,⁸, R. D. Blandford⁹, R. Bonino¹⁰,¹¹, P. Bruel¹², S. Buson¹³, R. A. Cameron⁹, P. A. Caraveo¹⁴, F. Casaburo¹⁵,¹⁶, E. Cavazzuti¹⁷, C. C. Cheung¹⁸, G. Chiaro¹⁴, S. Ciprini¹⁵,¹⁶, G. Cozzolongo¹⁹,²⁰, P. Cristarella Orestano²¹,²², S. Cutini²², F. D'Ammando²³, N. Di Lalla⁹, F. Dirirsa²⁴, L. Di Venere⁸, A. Domínguez²⁵, S. J. Fegan¹², E. C. Ferrara²⁶,²⁷,²⁸, A. Fiori², Y. Fukazawa²⁹, S. Funk¹⁹, P. Fusco⁷,⁸, F. Gargano⁸, S. Garrappa³⁰, D. Gasparrini¹⁵,¹⁶, S. Germani²²,³¹, N. Giglietto⁷,⁸, F. Giordano⁷,⁸, M. Giroletti²³, D. Green³², I. A. Grenier³³, S. Guiriec²⁸,³⁴, E. Hays²⁸, D. Horan¹², M. Kuss⁵, S. Larsson³⁵,³⁶, M. Laurenti¹⁵,¹⁶, J. Li³⁷,³⁸, I. Liodakis³⁹, F. Longo³,⁴, F. Loparco⁷,⁸, B. Lott⁴⁰, M. N. Lovellette⁴¹, P. Lubrano²², S. Maldera¹⁰, D. Malyshev¹⁹, A. Manfreda², L. Marcotulli⁴²,⁴³, G. Martí-Devesa⁴, M. N. Mazziotta⁸, I. Mereu²¹,²², P. F. Michelson⁹, W. Mitthumsiri⁴⁴, T. Mizuno⁴⁵, M. E. Monzani⁹,⁴⁶, A. Morselli¹⁵, I. V. Moskalenko⁹, M. Negro⁴⁷, N. Omodei⁹, M. Orienti²³, E. Orlando⁹,⁴⁸, J. F. Ormes⁴⁹, D. Paneque³², M. Perri¹⁶,⁵⁰, M. Persic³,⁵¹, M. Pesce-Rollins⁵, T. A. Porter⁹, G. Principe³,⁴,²³, S. Rainò⁷,⁸, R. Rando⁵²,⁵³,⁵⁴, B. Rani²⁸,⁵⁵, M. Razzano², A. Reimer⁵⁶, O. Reimer⁵⁶, P. M. Saz Parkinson⁵⁷, L. Scotton⁵⁸, D. Serini⁸, A. Sesana⁵⁹,⁶⁰, C. Sgrò⁵, E. J. Siskind⁶¹, G. Spandre⁵, P. Spinelli⁷,⁸, D. J. Suson⁶², H. Tajima⁶³,⁶⁴, M. N. Takahashi³²,⁶⁵, D. Tak⁶⁶, J. B. Thayer⁹, D. J. Thompson²⁸, D. F. Torres⁶⁷,⁶⁸, J. Valverde²⁸,⁵⁵, F. Verrecchia¹⁶,⁵⁰, and G. Zaharijas⁶⁹

¹IRAP, Université de Toulouse, CNRS, UPS, CNES, F-31028 Toulouse, France
²Università di Pisa and Istituto Nazionale di Fisica Nucleare, Sezione di Pisa I-56127 Pisa, Italy
³Istituto Nazionale di Fisica Nucleare, Sezione di Trieste, I-34127 Trieste, Italy
⁴Dipartimento di Fisica, Università di Trieste, I-34127 Trieste, Italy
⁵Istituto Nazionale di Fisica Nucleare, Sezione di Pisa, I-56127 Pisa, Italy
⁶California State University Los Angeles, Department of Physics and Astronomy, Los Angeles, CA 90032, USA
⁷Dipartimento di Fisica “M. Merlin” dell’Università e del Politecnico di Bari, via Amendola 173, I-70126 Bari, Italy
⁸Istituto Nazionale di Fisica Nucleare, Sezione di Bari, I-70126 Bari, Italy
⁹W. W. Hansen Experimental Physics Laboratory, Kavli Institute for Particle Astrophysics and Cosmology, Department of Physics and SLAC National Accelerator Laboratory, Stanford University, Stanford, CA 94305, USA
¹⁰Istituto Nazionale di Fisica Nucleare, Sezione di Torino, I-10125 Torino, Italy
¹¹Dipartimento di Fisica, Università degli Studi di Torino, I-10125 Torino, Italy
¹²Laboratoire Leprince-Ringuet, CNRS/IN2P3, École polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France
¹³Institut für Theoretische Physik und Astrophysik, Universität Würzburg, D-97074 Würzburg, Germany
¹⁴INAF-Istituto di Astrofisica Spaziale e Fisica Cosmica Milano, via E. Bassini 15, I-20133 Milano, Italy
¹⁵Istituto Nazionale di Fisica Nucleare, Sezione di Roma “Tor Vergata”, I-00133 Roma, Italy; stefano.ciprini@roma2.infn.it
¹⁶Space Science Data Center - Agenzia Spaziale Italiana, Via del Politecnico, snc, I-00133, Roma, Italy
¹⁷Italian Space Agency, Via del Politecnico snc, I-00133 Roma, Italy
¹⁸Space Science Division, Naval Research Laboratory, Washington, DC 20375-5352, USA
¹⁹Friedrich-Alexander Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics, Erwin-Rommel-Str. 1, 91058 Erlangen, Germany
²⁰Friedrich-Alexander-Universität, Erlangen-Nürnberg, Schlossplatz 4, 91054 Erlangen, Germany
²¹Dipartimento di Fisica, Università degli Studi di Perugia, I-06123 Perugia, Italy; paolo.cristarellaostano@studenti.unipg.it
²²Istituto Nazionale di Fisica Nucleare, Sezione di Perugia, I-06123 Perugia, Italy; sara.cutini@pg.infn.it
²³INAF Istituto di Radioastronomia, I-40129 Bologna, Italy
²⁴Astronomy and Astrophysics Research Development Department, Entoto Observatory and Research Center, Space Science and Geospatial Institute, Addis Ababa, Ethiopia
²⁵Grupo de Altas Energías, Universidad Complutense de Madrid, E-28040 Madrid, Spain
²⁶Department of Astronomy, University of Maryland, College Park, MD 20742, USA
²⁷Center for Research and Exploration in Space Science and Technology (CRESST) and NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA
²⁸Astrophysics Science Division, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA
²⁹Department of Physical Sciences, Hiroshima University, Higashi-Hiroshima 739-8526, Japan
³⁰Department of Particle Physics and Astrophysics, Weizmann Institute of Science, 76100 Rehovot, Israel
³¹Dipartimento di Fisica e Geologia, Università degli Studi di Perugia, via Pascoli snc, I-06123 Perugia, Italy
³²Max-Planck-Institut für Physik, D-80805 München, Germany
³³Université Paris Cité, Université Paris-Saclay, CEA, CNRS, AIM, F-91191 Gif-sur-Yvette, France
³⁴The George Washington University, Department of Physics, 725 21st St., NW, Washington, DC 20052, USA
³⁵Department of Physics, KTH Royal Institute of Technology, AlbaNova, SE-106 91 Stockholm, Sweden; stefan@astro.su.se
³⁶The Oskar Klein Centre for Cosmoparticle Physics, AlbaNova, SE-106 91 Stockholm, Sweden
³⁷CAS Key Laboratory for Research in Galaxies and Cosmology, Department of Astronomy, University of Science and Technology of China, Hefei 230026, People's Republic of China
³⁸School of Astronomy and Space Science, University of Science and Technology of China, Hefei 230026, People's Republic of China
³⁹NASA Marshall Space Flight Center, Huntsville, AL 35812, USA
⁴⁰Université Bordeaux, CNRS, LP2I Bordeaux, UMR 5797, F-33170 Gradignan, France
⁴¹The Aerospace Corporation, 14745 Lee Rd., Chantilly, VA 20151, USA
⁴²Department of Astronomy, Department of Physics and Yale Center for Astronomy and Astrophysics, Yale University, New Haven, CT 06520-8120, USA
⁴³Department of Physics and Astronomy, Clemson University, Kinard Lab of Physics, Clemson, SC 29634-0978, USA
⁴⁴Department of Physics, Faculty of Science, Mahidol University, Bangkok 10400, Thailand
⁴⁵Hiroshima Astrophysical Science Center, Hiroshima University, Higashi-Hiroshima, Hiroshima 739-8526, Japan
⁴⁶Vatican Observatory, Castel Gandolfo, V-00120, Vatican City State, Italy
⁴⁷Department of physics and Astronomy, Louisiana State University, Baton Rouge, LA 70803, USA
⁴⁸Istituto Nazionale di Fisica Nucleare Sezione di Trieste and Università di Trieste, I-34127 Trieste, Italy
<sup>49</sup>Department of Physics and Astronomy, University of Denver, Denver, CO 80208, USA

<sup>50</sup>INAF Astronomical Observatory of Rome, via Frascati 33, I-00078, Monte Porzio Catone, Roma, Italy

<sup>51</sup>INAF Astronomical Observatory of Padova, Vicolo dell'Osservatorio 5, I-35122 Padova, Italy

<sup>52</sup>Dipartimento di Fisica e Astronomia “G. Galilei,” Università di Padova, Via F. Marzolo, 8, I-35131 Padova, Italy

<sup>53</sup>Istituto Nazionale di Fisica Nucleare, Sezione di Padova, I-35131 Padova, Italy

<sup>54</sup>Center for Space Studies and Activities “G. Colombo,” University of Padova, Via Venezia 15, I-35131 Padova, Italy

<sup>55</sup>Center for Space Science and Technology, University of Maryland Baltimore County, 1000 Hilltop Cir., Baltimore, MD 21250, USA

<sup>56</sup>Institut für Astro- und Teilchenphysik, Leopold-Franzens-Universität Innsbruck, A-6020 Innsbruck, Austria

<sup>57</sup>Santa Cruz Institute for Particle Physics, Department of Physics and Department of Astronomy and Astrophysics, University of California at Santa Cruz, Santa Cruz, CA 95064, USA

<sup>58</sup>Center for Space Plasma and Aeronomic Research (CSPAR), University of Alabama in Huntsville, Huntsville, AL 35899, USA

<sup>59</sup>INFN Sezione di Milano-Bicocca, Piazza della Scienza 3, 20126 Milano, Italy

<sup>60</sup>Dipartimento di Fisica, Università degli Studi di Milano-Bicocca, I-20126 Milano, Italy

<sup>61</sup>NYCB Real-Time Computing Inc., Lattingtown, NY 11560-1025, USA

<sup>62</sup>Purdue University Northwest, Hammond, IN 46323, USA

<sup>63</sup>Nagoya University, Institute for Space-Earth Environmental Research, Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan

<sup>64</sup>Kobayashi-Maskawa Institute for the Origin of Particles and the Universe, Nagoya University, Furo-cho, Chikusa-ku, Nagoya, Japan

<sup>65</sup>Institute for Cosmic-Ray Research, University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba, 277-8582, Japan

<sup>66</sup>SNU Astronomy Research Center, Seoul National University, Gwanak-rho, Gwanak-gu, Seoul, Republic of Korea

<sup>67</sup>Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Magrans s/n, E-08193 Barcelona, Spain; and Institut d'Estudis Espacials de Catalunya (IEEC), E-08034 Barcelona, Spain

<sup>68</sup>Institució Catalana de Recerca i Estudis Avançats (ICREA), E-08010 Barcelona, Spain

<sup>69</sup>Center for Astrophysics and Cosmology, University of Nova Gorica, Nova Gorica, Slovenia

Received 2024 April 18; revised 2024 July 6; accepted 2024 July 16; published 2024 November 25

## Abstract

A 2.1 yr periodic oscillation of the gamma-ray flux from the blazar PG 1553+113 has previously been tentatively identified in ~7 yr of data from the Fermi Large Area Telescope. After 15 yr of Fermi sky-survey observations, doubling the total time range, we report >7 cycle gamma-ray modulation with an estimated significance of $4\sigma$ against stochastic red noise. Independent determinations of oscillation period and phase in the earlier and the new data are in close agreement (chance probability <0.01). Pulse timing over the full light curve is also consistent with a coherent periodicity. Multiwavelength new data from Swift X-Ray Telescope, Burst Alert Telescope, and UVOT, and from KAIT, Catalina Sky Survey, All-Sky Automated Survey for Supernovae, and Owens Valley Radio Observatory ground-based observatories as well as archival Rossi X-Ray Timing Explorer satellite-All Sky Monitor data, published optical data of Tuorla, and optical historical Harvard plates data are included in our work. Optical and radio light curves show clear correlations with the gamma-ray modulation, possibly with a nonconstant time lag for the radio flux. We interpret the gamma-ray periodicity as possibly arising from a pulsational accretion flow in a sub-parsec binary supermassive black hole system of elevated mass ratio, with orbital modulation of the supplied material and energy in the jet. Other astrophysical scenarios introduced include instabilities, disk and jet precession, rotation or nutation, and perturbations by massive stars or intermediate-mass black holes in polar orbit.

*Unified Astronomy Thesaurus concepts: Active galactic nuclei (16); Blazars (164); Gamma-rays (637); Gamma-ray sources (633); Gamma-ray telescopes (634); Gamma-ray astronomy (628); BL Lacertae objects (158); Jets (870); High energy astrophysics (739); Relativistic jets (1390); Galaxy accretion disks (562); Supermassive black holes (1663)*

## 1. Introduction

Blazars are active galactic nuclei (AGNs), harboring accreting supermassive black holes (SMBHs) with powerful relativistic jets almost aligned with our line of sight (Urry & Padovani 1995). One of the more remarkable such objects is PG 1553+113, a high-energy-peaked BL Lac object, which is observed over the full electromagnetic spectrum, from radio to very high energy (VHE; >100 GeV) gamma rays. PG 1553+113 was discovered as an optical and X-ray source (1ES 1553+113, $z \simeq 0.43$, Falomo & Treves 1990; Abramowski et al. 2015; Nicastro et al. 2018; Dorigo Jones et al. 2022) but has been most systematically monitored in gamma rays by the Large Area Telescope (LAT) on board the Fermi Gamma-ray Space Telescope, which was launched in 2008.

Early LAT observations raised suspicions of an almost sinusoidal modulation in the light curve of PG 1553+113, and a tentative identification of a 2.1 yr periodic or possibly quasiperiodic oscillation (QPO) was described in Ackermann et al. (2015a, hereafter A15a).

A15a reported results based on LAT data from 2008 August 5 to 2015 July 19, MJD 54682.65 to MJD 57222.65 (a 6.9 yr range), with radio and optical flux modulations in association with the gamma-ray oscillation.

Several papers have followed up on the results in A15a, mainly corroborating the QPO discovery in gamma-ray and optical bands with more data and more analyses (for example, Cavaliere et al. 2017; Sobacchi et al. 2017; Sandrinelli et al. 2018; Tavani et al. 2018; Yan et al. 2018; Cavaliere et al. 2019; Ait Benkhali et al. 2020; Covino et al. 2020; Peñil et al. 2020; Agarwal et al. 2021, 2022; Ren et al. 2023; Chen et al. 2024; Peñil et al. 2024a). Further blazars are also tentatively identified as producing correlated GeV gamma-ray and optical QPOs behavior (for example, Otero-Santos et al. 2020; Sarkar et al. 2021; Peñil et al. 2024a, 2024b).

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
The fact that Fermi LAT is producing long and essentially continuous light curves makes these data ideal to use for searches for periodic or quasiperiodic variability. This, together with the indication of a periodic oscillation in PG 1553+113, has stimulated systematic studies of sources sampled from the Fermi-LAT Third and Fourth Catalogs of Active Galactic Nuclei (3LAC; Ackermann et al. 2015b and 4LAC; Ajello et al. 2020; for example, Finke & Becker 2014; Itoh et al. 2016; Sandrinelli et al. 2016; Prokhorov & Moraghan 2017; Sandrinelli et al. 2018; Meyer et al. 2019; Ryan et al. 2019; Ait Benkhali et al. 2020; Peñil et al. 2020; Tarnopolski et al. 2020; Yang et al. 2021; Peñil et al. 2024a, 2024b).

The literature contains many studies, based on LAT and other data, with claims of QPO detections in variable blazar light curves. Many of these claims are controversial, primarily because blazars and other AGN exhibit erratic red noise variability that can easily mimic periodicity, in short time sequences containing few oscillation cycles (Corbet et al. 2007; Do et al. 2009; Lasky et al. 2015; Vaughan et al. 2016). A common approach to avoid, or at least limit, this risk is to estimate period significance relative to simulated red noise. The present consensus is that the most sound simulation method available is the one described in Emmanoulopoulos et al. (2013).

A claim of binary SMBH, with characteristics similar to PG 1553+113, was recently presented by O'Neill et al. (2022), for the BL Lac object PKS 2131-021, with a 4.8 yr radio-band flux periodicity (2.1 yr in the rest frame). The three brightest epochs in the LAT gamma-ray light curve<sup>70</sup> appear to precede the radio maxima by about a year.

In A15a, we introduced a close (0.005–0.01 pc separation) gravitationally bound binary system of SMBHs (total mass ~5 × 10<sup>8</sup> M<sub>☉</sub>) as a possible explanation for the observed periodicity. This would place the system in a state prior to the early inspiral, low-frequency, gravitational-wave (GW) driven regime, in agreement with Dhiman et al. (2021). This state, with weak GW emission, would last t ~ 10<sup>5</sup>–10<sup>6</sup> yr, followed by a rapid inspiral and merger of the two BHs (Peters 1964; Shapiro & Teukolsky 1983).

The scenario with a sub-parsec binary has been adopted or discussed in follow-up papers, for example, Prokhorov & Moraghan (2017), Ciprini (2018), Ait Benkhali et al. (2020), Covino et al. (2020), and Peñil et al. (2020, 2024a), generally with the gamma-ray QPO confirmed at an ~3σ–4σ level.

Evaluating the different proposed models for the QPO requires additional information, which might be provided by other types of observations. These include polarization and radio (parsec scale) flux structure long-term data as important probes of any regular wobbling (Caproni et al. 2017; Lico et al. 2020), helical patterns, and polarimetric rotations in the jet, and potential identification of double-component spectroscopic line profiles (Graham et al. 2015a).

In multi-epoch GHz milliarcsecond resolution maps by VLBA contained in the MOJAVE/2 cm Survey Data Archive (Lister et al. 2009), the radio core-dominated and magnetic-field-dominated BL Lac object PG 1553+113 shows a limb-brightened and diffuse jet structure. These features might be the first indications of a jet with a precessing baseline, wobbling motion by torques of the inner disk, or alternatively, an intrinsic pulsational accretion flow (Lister et al. 2013; Caproni et al.

2017; Lico et al. 2020). Such geometric dynamics would cause periodic changes in the relativistic Doppler boosting, magnetic field stresses, magnetic reconnections, and jet MHD instabilities (Cavaliere et al. 2017; Sobacchi et al. 2017; Tavani et al. 2018; Cavaliere et al. 2019; Rieger 2019; Lico et al. 2020).

In this paper we describe in Section 2 the Fermi-LAT gamma-ray data, together with multiwavelength X-ray (Swift-XRT and Rossi X-Ray Timing Explorer (RXTE) All Sky Monitor (ASM)) data, optical (Tuorla, KAIT, Catalina Sky Survey (CSS), Swift-UVOT, All-Sky Automated Survey for Supernovae (ASAS-SN)) data, and 15 GHz data from the Owens Valley Radio Observatory (OVRO) radio telescope. In Section 3, we report multiple approaches in time-domain variability/periodicity analysis, cross correlation, and estimation of the gamma-ray QPO coherence, with the findings suggesting a potential energy dependence. In Section 4, we introduce some possible astrophysical interpretations, and in Section 5, we summarize our findings. Appendix B contains optical historical data from plates (years 1912–1988), and RXTE ASM data (from 1996 to 2011) for reference.

## 2. Gamma-Ray, X-Ray, Optical, and Radio-band Light-curve Data

### 2.1. Fermi-LAT Data

The Fermi-LAT is a pair conversion detector with a 2.4 sr field of view, sensitive to gamma rays from ~20 MeV to >300 GeV (Atwood et al. 2009). The present work uses the Pass 8 LAT database (Atwood et al. 2013). The Fermi-LAT operating mode allows it to cover the entire sky every two spacecraft orbits (~1.6 hr orbital period), providing a regular and uniform view of gamma-ray sources, sampling timescales from hours to years. This work uses observations of PG 1553+113 covering ~15.3 yr (2008 August 4 to 2023 November 13, Modified Julian Day, MJD, 54682.65–60261.65). The LAT data analysis employed the standard *ScienceTools v11r01p01*<sup>71</sup> package, selecting events from 100 MeV to 300 GeV with *P8R3_SOURCE_V2* instrument response functions, in a circular region of interest (ROI) of 20° radius centered on the position of PG 1553+113 (15<sup>b</sup>55<sup>m</sup>43<sup>s</sup>0440 +11<sup>d</sup>11<sup>m</sup>24<sup>s</sup>365, J2000, Beasley et al. 2002). We applied a zenith angle cut of >90° in order to reduce the contamination from the Earth limb and the standard data quality cuts (*DATA QUAL* > 0 && *LAT_CONFIG* == 1) for the extraction of good time intervals. We selected only “SOURCE” event class (*LAT EVENT_CLASS* = 128) and event type “FRONT+BACK” (*LAT EVENT_TYPE* = 3). We used files *gll_iem_v07* and *iso_P8R3_SOURCE_V2_v02* to model the Galactic and isotropic diffuse emission, respectively.

A binned maximum likelihood model fit technique was applied to each time bin with a power-law spectral model and photon index left free to vary for PG 1553+113. The background model included sources from the 4FGL catalog (Abdollahi et al. 2020) within the ROI of 20°; the fluxes and photon indices of background sources within a radius of 10° from the source of interest were left free to vary.

The resulting light curves, in different time and energy bins are shown in Figure 1.

The gamma-ray flux modulation oscillation is also found using different time bin sizes. In particular in Appendix A we

<footnote n="70"><sup>70</sup> https://fermi.gsfc.nasa.gov/ssc/data/access/lat/LightCurveRepository/</footnote>

<footnote n="71"><sup>71</sup> http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Fermi-LAT gamma-ray integrated flux light curves and spectral index time profile of PG 1553+113 over more than 15 yr, from 2008 August to 2023 November. The light curves above 100 MeV and above 1 GeV are shown with a constant 45 day time bin (first and middle panels); in the last panel we show the Fermi LAT spectral index ($E > 100$ MeV) as a function of time for the 45 day time bins, with a scaled version of the ($E > 100$ MeV) gamma-ray light curve superimposed as a dashed line.</figcaption></figure>

report results from the public data from the Fermi-LAT light curve repository at the FSSC (Fermi Large Area Telescope Collaboration 2021; Abdollahi et al. 2023).

## 2.2. Neil Gehrels Swift X-Ray Telescope and Burst Alert Telescope Data

Four hundred thirty five, 435, Swift observations of PG 1553+113 between 2005 April 20 and 2023 May 31 (MJD 53480.2–MJD 60095.8) are analyzed. The Swift X-Ray Telescope (XRT) unabsorbed 0.3–2 keV flux light curve is presented in Figure 2. While all the PC-mode data were used, WT-mode data are omitted for the last 1.5 yr of this data range.

XRT data were first calibrated and cleaned (xrtpipeline, XRTDAS v.3.7.0), and energy spectra were extracted from a region of 20 pixel (~47′′) radius, with a nearby 20 pixel radius region for background. Individual XRT spectra were well fitted with a log-parabolic model, with neutral hydrogen equivalent column density fixed to the Galactic value of $3.6 \times 10^{20} \text{ cm}^{-2}$ (Kalberla et al. 2005). For completeness we have also reported the historical public light curve from the Swift-Burst Alert Telescope (BAT). We depict only data points that have a relative uncertainty flux/Δflux > 0.5. The data are publicly

provided by the Swift-BAT team (Krimm et al. 2013) through a dedicated webpage at NASA Goddard Space Flight Center.<fnref n="72" />

## 2.3. RXTE ASM Data

X-ray flux densities in the 1.5–12 keV band were calculated from the public high-level archive (Levine et al. 1996) of the ASM on board RXTE for PG 1553+113 in the 15.7 yr long data set from 1996 January 6 through 2011 September 23 (MJD 50088–55827), using daily binning (see Appendix A for the complete 1 day bin count rate time series). Data after MJD 55450 (last year of the light curve) are of poor quality because of the RXTE cathode loss.

Since ASM source X-ray counts were given as negative on many days (fluxes below the detection threshold), timing analysis via fast Fourier transform and Lomb–Scargle periodogram (LSP) power spectra, was applied to both the original data and to data where we omitted negative counts and upper limits, using 1 day bins and other binnings as checks. The equivalence 1 Crab unit = 75 ASM count s⁻¹ and the conversion 1 Crab [2–10 keV] = $2.4 \times 10^{-8} \text{ erg cm}^{-2} \text{ s}^{-1}$ are adopted

<footnote n="72">⁷² https://swift.gsfc.nasa.gov/results/transients/weak/QSOB1553p113/</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Multiwavelength light curves of PG 1553+113 at X-ray and soft gamma-ray, optical, and radio bands with the LAT gamma-ray light curve superimposed as a dashed gray line (*E* > 100 MeV). Top panel: Swift-XRT integrated flux (0.3–2.0 keV) shown as blue-filled circles, and RXTE data points averaged into 30 day time bins extrapolated to the 0.3–2.0 keV energy flux (Swift-XRT energy band) shown as blue empty diamonds. Swift-BAT count rate at 15–50 keV, multiplied by an arbitrary constant to scale it to the X-ray light curves, shown in green. Central panel: optical flux density from Tuorla (*R* filter, black-filled circle points), CSS (V filter rescaled, blue-filled squared points), KAIT (V filter rescaled, red-filled diamond points), Swift-UVOT (V filter rescaled, green-filled circles) and ASAS (V filter rescaled, light-blue-filled points). Bottom panel: 15 GHz flux density from OVRO 40 m (black-filled circles) and parsec-scale 15 GHz flux density from VLBA (MOJAVE program, yellow-filled diamonds).</figcaption></figure>

(factor 2.79 × 10⁻⁸ erg cm⁻² s⁻¹ to compare with the XRT 0.3–2 keV fluxes in Figure 2, evaluated from Kirsch et al. 2005).

## 2.4. Optical, Swift-UVOT, KAIT, Tuorla, Catalina Sky Survey, and ASAS-SN

Aperture photometry for the UVOT V-band filter was performed. The task is included in the HEASoft package (v6.23 by NASA HEASARC; 5" radius source aperture and two 18" apertures for the background evaluation). Counts were converted to fluxes using the standard zero-points (Breeveld et al. 2010), and dereddened using the appropriate values of *E*(B–V) (Schlafly & Finkbeiner 2011) for the effective wavelengths of the UVOT filters. Tuorla blazar monitoring program (Takalo et al. 2008) data are those published in A15a, obtained with a 35 cm telescope attached to the Kungliga Vetenskapsakademien (KVA) telescope on La Palma and the Tuorla 1 m telescope in Finland (A15; Nilsson et al. 2018), with optical observations typically performed two to three times per week.

Public data from the Katzman Automatic Imaging Telescope (KAIT), CSS, and ASAS programs were also added. V-band

magnitudes were scaled with an inter-calibrated fixed offset, (from an average of daily V–R) to the R-band values.

Figure 2 shows multifrequency light curves. The light curves in different wavelengths indicate a coherent periodic behavior from radio to gamma ray. This is true especially for the optical and gamma-ray fluxes, also apparent in Figure 3.

## 2.5. Radio OVRO Data

Figure 2 shows flux density, long-term monitoring of PG 1553+113 at 15 GHz obtained by the OVRO 40 m radio telescope,<fnref n="73" /> MOJAVE program, with the observation cadence of this blazar between 1 and 23 days since August 2008 (Richards et al. 2011). A cooled receiver (3.0 GHz bandwidth, centered on 15.0 and 2.5 GHz reception bandwidth) is placed at the prime focus of the 40 m radio telescope, with dual off-axis corrugated horn feed, projecting to Gaussian beams (157" FWHM) on the sky separated in azimuth by 12.95, and with flux densities measured using azimuth double switching after peaking up on-source (Readhead et al. 1989). The OVRO flux

<footnote n="73">⁷³ http://www.astro.caltech.edu/ovroblazars/</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Left panel: Fermi-LAT epoch-folded ($E > 100$ MeV) flux light curve with superimposed Fourier pulse fit, showing a slight bimodal peak. Right panel: superposed Fermi-LAT ($E > 1$ GeV) and optical, rescaled, light curves for visual comparison (Tuorla, CSS, KAIT, Swift-UVOT). The periodic gamma-ray Fourier oscillating pulse precursor/maximum (P/M) epoch zones are also highlighted. These zones are taken to be 120 days long with a start at MJD 57440 and period of 796.5 days.</figcaption></figure>

density data have a minimum uncertainty of 4 mJy in 32 s of on-source integration, and a typical rms relative error of 3%.

### 3. Variability Analysis

The ~2.1 yr period in PG 1553+113 tentatively identified in A15a was based on data covering only approximately three oscillation cycles. The light curves used in the present study are more than twice as long, which allows us to make several analysis improvements:

1. Use a pulsation model built on the A15a data and compare it with the phase and frequency of a similar model built on the new data. Since the two data sets are independent we can estimate the probability that two oscillations by chance will have a particular phase and frequency.

2. Make a more strict significance estimate for the periodic oscillation in comparison with simulated stochastic red noise.

3. Perform a wavelet time-frequency study of pulsation evolution along time, and significance in relation to a theoretical stochastic (random walk) first-order autoregressive, AR(1), model.

4. Improve on earlier estimates of the correlation between the optical, radio, and gamma-ray variations in PG 1553+113.

5. Set limits on the coherence of the oscillation.

6. Study the energy dependence of the gamma-ray QPO.

Below we describe each of these topics and present the associated results.

#### 3.1. Testing the Periodicity Hypothesis with New Data and Power Density Spectra Pulse Analysis

The analysis in A15a showed that the Fermi-LAT light curve for the first 6.9 yr could be described as a periodic modulation. By assuming a strictly periodic modulation and estimating its

parameters from the data in A15a it is now possible to test that hypothesis by comparison with an independent fit of period and phase to the new data. In the present analysis, we have done so by fitting a Fourier pulse (fundamental period plus three overtones) to the light curve ($E > 1$ GeV) of the first ~7 yr and then extrapolating that modulation to later times. The energy selection and fitting procedure are identical to the corresponding analysis in A15a. To avoid biasing the hypothesis test we tried no other procedures. Before the analysis, the full light curve was detrended by subtracting a least square fitted line. The light curve was then normalized to mean 1 and the same fractional variability (rms/flux) as the original light curve. The normalized light curve is shown in the upper part of Figure 4. The resulting Fourier fit is also shown in the figure. The dashed vertical line at MJD 57203 indicates the separation between the data that were used for the Fourier fit (corresponding to A15a) and the new data. The phase difference between the two oscillations is estimated by a cross correlation between the extrapolated Fourier template and the new data after MJD 57203. The estimated time lag in phase units of the template with respect to the new data is $0.117 \pm 0.033$, where the error value only takes the statistical uncertainties of the new data into account. The total uncertainty of the phase difference is dominated by the model extrapolation from the first observing epoch. The uncertainty range for the extrapolated phase at the midpoint of the new data is between +0.094 and -0.062. We conclude that the observed phase difference between the predicted and observed phase in the new data is marginally consistent when phase uncertainties are taken into account.

The oscillation periods for the two data segments, before and after MJD 57203, can be obtained from the power density spectra (PDS). These are shown in the lower part of Figure 4. Each PDS is oversampled by a factor of 4, which matches the frequency uncertainty in our Fourier fit analysis. The uncertainty range is approximately two-thirds of the width of the frequency bins. The difference in peak frequency for the two segments is approximately $5 \times 10^{-5}$ day$^{-1}$, which is less
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Top panel: a Fourier fit to the linearly detrended 1 GeV light curve up to MJD 57203, the data used in A15a. The fitted oscillation is extrapolated and shown as a dashed line after MJD 57203 for comparison with the new data. Lower panel: comparison of PDSs for the two independent data segments (blue: *t* < MJD 57203, purple: *t* > MJD 57203).</figcaption></figure>

than the widths of the frequency bin (1.01 × 10⁻⁴ and 0.83 × 10⁻⁴ day⁻¹, respectively). Assuming that the period in the new data could fall in any of the 118 frequency bins (neglecting the first four due to the linear detrending of the light curve), the chance probability for the two oscillations to be so close in frequency is 0.0085. The similarity in the pulsation phase described above gives additional weight to the confirmation of the pulsation model for the previous data. A conservative conclusion is that the chance probability for the pulsation period and phase of the new independent data to be close to those of the earlier data is less than 0.01.

## 3.2. Estimating QPO Significance with Pure-noise Light-curve Simulations

Periodicity studies of astronomical time series deal not only with the estimate of the period and its uncertainty, but with whether there is a (quasi)-regular periodicity or not. The statistical significance for such periodicity must then be defined. This calculation is not trivial: significance can be estimated using the false alarm probability (FAP; VanderPlas 2018), which measures the probability that the data are dominated by white-colored noise would lead to a signal peak in the periodogram of higher power than the peak due to a truly periodic signal. The FAP can be estimated through light curve simulations and the LSP (Lomb 1976; Scargle 1982), as the

<table><caption>**Table 1**
Summary of the Gamma-Ray QPO Significance in σ for the Various LAT Light Curves of PG 1553+113 Presented or Discussed in This Work</caption><thead><tr><th rowspan="2">Gamma-Ray Band</th><th colspan="4">E > 100 MeV</th><th>E > 1 GeV</th></tr><tr><th colspan="2">Photon Flux</th><th colspan="2">Energy Flux</th><th colspan="2">Photon Flux</th></tr></thead><tbody><tr><td>Light Curve Type</td><td colspan="2"></td><td colspan="2"></td><td colspan="2"></td></tr><tr><td>Bin size</td><td>20 days</td><td>45 days</td><td>20 days</td><td>45 days</td><td>45 days</td><td></td></tr><tr><td>Light curve original</td><td>1.5σ</td><td>1.7σ</td><td>3σ</td><td>>4σ</td><td>3.5σ</td><td></td></tr><tr><td>Light curve detrended</td><td>2σ</td><td>2.3σ</td><td>3σ</td><td>4σ</td><td>4σ</td><td></td></tr></tbody></table>

ratio *n*/N: the number *n* of simulated time series with power peaks higher than the original light curve peak found, over the total number of simulations (*N* = 10⁶ in our case). The FAP is not a direct estimate of the periodicity significance, but rather the probability of a false alarm, which is substantially different. Simulations are randomizations of the original, true, light-curve data. However, the brown/red/pink, colored, noise power spectra of true light curves of highly variable blazars (and many other astronomical source types) tend to result in over-estimated, not reliable, periodicity significances (Emmanoulopoulos et al. 2013; Vaughan et al. 2016).

In order to simulate pure-noise data with the same characteristics as the true gamma-ray light curves of PG 1553+113, we used the Emmanoulopoulos algorithm (Emmanoulopoulos et al. 2013), to generate surrogate simulated light curves with PDF and PDS similar to the true data light curve. The resulting *n*/N FAP value, evaluated as the greatest LSP power peak of each simulated light curve, is difficult to interpret; therefore, a corresponding significance, equivalent to the number of standard deviations σ calculated from the cumulative of a normal distribution up to 1 − *n*/N, is introduced. If none of the simulations has a greater peak than the original peak (*n* = 0), a lower limit on the FAP must be set. For PG 1553+113 we choose a lower limit of 1/*N*, and a corresponding upper limit on the value of σ is set too.

We apply this approach for the different gamma-ray light curves for PG 1553+113 presented in this work (Table 1): photon/energy flux, fluxes for (E > 1 GeV) and (E > 100 MeV), 45 and 20 day time bin light curves (also 3 day bin data reported in Appendix A). For each, we calculate 10⁶ simulations. The gamma-ray light curves have a hint of a secular increase in flux, and detrended light curves (and their proper simulations) are also evaluated for comparison. The results are reported in Table 1. The distributions for both of the simulated data for the true 45 day photon and energy flux light curves, are shown as an example in Figure 5.

## 3.3. Wavelet Time-frequency Study of the QPO Signal's Evolution

In a generalized Fourier space, the continuous wavelet transform (CWT) is able to identify both the dominant modes of variability and how those modes vary in time. CWT can also be quite sensitive to signals that would go undetected by other methodologies (e.g., Mallat 2008). In Figure 6, we report the diffuse two-dimensional, time-frequency *t* versus *f* = 1/*T*, “scalogram” plots, with the normalized 2D modulus of the CWT energy density function ||*W*<sub>*n*</sub>(*s*)||²/σ², where 1/σ² measures the power relative to white noise, corresponding to
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Signal power with the LSP peak distributions of $10^6$ simulated light curves, for the 45 day bin ($E > 100$ MeV) energy flux light curve in the inset and for the 45 day bin ($E > 100$ MeV) photon flux light curve in the main panel. Dashed lines are the LSP power for the two, true data, light curves. The significance of the period here is >4σ and 1.7σ for the energy flux and photon flux light curves, respectively.</figcaption></figure>

the LAT $E > 1$ GeV and $E > 100$ MeV gamma-ray, XRT X-ray, and optical light curves of Figures 1 and 2. The Morlet mother waveform, a wave packet of a plane wave modulated by a Gaussian, is used here as the best trade-off between epoch localization in *t* and resolution in frequency *f* (or equally period/periodicity timescale *T*). Lighter colors describe higher peaks, bumps, and regions, and hence stronger values of the 2D CWT power spectrum.

The orange-filled color scale 2D image contour plots in the top panel of Figure 6, result from the Morlet CWT gamma-ray power spectrum for the $E > 100$ MeV and $E > 1$ GeV, 45 day bin, light curves reported in Figure 1. Side plots show the global CWT power spectrum (i.e., averaged along the light curve MJD epochs on the *x*-axis) and the LSP power spectrum, in the bias-corrected implementation given by the REDFIT method and software (Schulz & Mudelsee 2002), which also provides a mathematical first-order autoregressive AR(1) model fitting, adequate for stochastic natural and astronomical processes. The method works directly on unevenly binned/gapped data, avoiding the introduction of bias during interpolation. Also, we did not manipulate the light curve data, avoiding any linear detrending prior to analysis, even if this could have given us slightly higher significance values. With REDFIT, we applied a Welch window overlapped segment averaging, choosing to split the time series into two segments that overlap by 50% to reduce noise, while avoiding too great a reduction of the spectral resolution, and then averaging their power spectra. We obtained the same results with rectangular and triangle window functions. The REDFIT software also calculates the significance levels on the LSP spectrogram based on parametric approximation with the $\chi^2$ calculated with respect to the computed AR(1)-noise spectrum (the null hypothesis), where the degrees of freedom depend on the number of data points. It also calculates the significance level with 1000 Monte Carlo random realizations of the AR(1) random walk, a.k.a. red/brown/Brownian, noise model.

In general, any peak of the one-dimensional LSP spectro-
gram power above the 99% confidence curve level, has less
than a 1% chance of being a false alarm product of stochastic
red/Brownian AR(1) process realizations. In this case, the
parametric $\chi^2$ smooth and Monte Carlo AR(1) irregular
confidence curve levels reported in the four side panels in
Figure 6 have similar trends and intensity.

Both the >100 MeV and >1 GeV panels in Figure 6 show a coherent gamma-ray global CWT peak centered at ≈780 days (2.1 yr) with a significance of >3σ, which approaches 4σ for the 1 GeV flux light curve, as shown by the prominence of both of the bumps with respect to the 99% AR(1) noise false alarm confidence levels. Most of the signal power during this 15.3 yr gamma-ray light-curve epoch is manifestly placed around this timescale, well represented by the clear/white horizontal strip in both of the CWT 2D orange-scale image plots at the top, with such strips sufficiently outside the, spurious, “cone of influence” region. The CWT indicates therefore a dominant 2.1 yr QPO gamma-ray flux modulation, along the light curve time range.

In Figure 6, the 2D Morlet CWT (blue color scale) plot of the Swift-XRT integrated (0.3–2.0 keV) flux from MJD 55973, shows erratic variability and no significant QPO, i.e., a behavior substantially different from the GeV gamma-ray, optical, and radio periodic modulation.

The LSP of X-ray data indicates a <2σ potential peak around ~498 day timescale (~1.4 yr), not corroborated by the wavelet. The CWT has some X-ray signal in a region between ~1.4 and ~1.8 yr, extended at least from 2013 June to 2020 October. Two high-power white ovals (~900–1000 day timescale) are both placed at the edges of the X-ray light curve, denoting when relevant X-ray outbursts from 2012 and 2023 occurred. In such a cross-hatched region, spurious power by edge effects on the CWT scalogram is important. A stronger hint is recently reported in Aniello et al. (2024).

The 2D CWT (green color scale) plot of the 19.7 yr optical light curve has well-defined global CWT and LSP peaks at ≈790 days (≈2.16 yr), quite similar to the gamma-ray QPO, and a rather stable evolution in time along the *x*-axis of the 2D scalogram image. The significance is lower (~2σ) as shown by the 99% and 95% false alarm confidence levels. This supports a multiband QPO, also because the optical-gamma-ray cross correlation, even if such wavelet and LSP of the optical data, when taken standalone, can still be consistent with erratic red/brown-noise variability. We remember that, already in A15a, we reported an apparent 754 ± 20 day periodicity in the 9.9 yr range (MJD 53479-57110) for the optical data.

In this work, the wavelet analysis for the 15 GHz radio light
curve is not presented since the date range is similar to that
considered in A15a, and the results are similar.

The results of the wavelet analysis and the REDFIT spectrogram (>3σ significance gamma-ray QPO, 2σ significance optical QPO), although insufficient to claim a deterministic periodicity, raise a more solid possibility of there being a true QPO, when supported, as in this work, by the other methods that we illustrated in the previous sections.

3.4. Gamma Ray—Optical and Radio Correlation

The relationship between the gamma-ray, optical, and radio
variability in PG 1553+113 was investigated in A15a, using
the first 6.9 yr of Fermi-LAT observations. The time lag
between the optical and the Fermi-LAT ($E > 1$ GeV) light
curves, both 20 days binned, was found to be $10 \pm 51$ days,
with the optical leading but still consistent with zero lag. The
relation between the gamma-ray and radio was found to be
more complex with the 15 GHz flux lagging the gamma-ray
flux by some 2–3 months.

Using the updated light curves in Figure 2, we have
calculated the discrete cross-correlation function (DCCF) for
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Top panels: 2D image filled (orange color scale) contour plots of the Morlet CWT gamma-ray power spectrum, the “scalogram,” respectively for the Fermi-LAT $E > 100$ MeV and $E > 1$ GeV, 45 day bin, photon flux light curves of Figure 1. Side plots are the global CWT power spectrum averaged along the MJD epochs, and the one-dimension LSP power spectrum in the bias-corrected REDFIT implementation. The LSP power axis is rescaled to the scale of the global CWT comparing the relative main bumps. The 99% false alarm confidence level lines (parametric $\chi^2$ and Monte Carlo) against the null hypothesis of an origin by pure stochastic AR(1)-noise process realization, are shown. Lower panels: 2D scalogram-filled color scale contour plot of the Morlet CWT (blue color scale) of the Swift-XRT integrated (0.3–2.0 keV) flux density (omitting the few, too sparse, XRT data prior to the year 2012), and 2D scalogram-filled color scale contour plot of the Morlet CWT (green color scale) of the optical flux density, corresponding to light curves presented in Figure 2. Cross-hatched regions in all the panels are the “cone of influence,” where spurious edge effects caused by finite time series boundaries become important. Side plots are the corresponding global CWT and the LSP power with their false alarm confidence levels lines, scaled to the global CWT peak.</figcaption></figure>

Fermi-LAT ($E > 1$ GeV) with respect to optical and OVRO (15 GHz) flux. The optical R- and V-band observations shown in Figure 2 were normalized and merged into a single optical light curve. The optical measurements are very irregularly distributed in time and were therefore averaged into 20 day bins before the analysis. Since early OVRO measurements had lower signal to noise than later observations, the data up to MJD 57200 were sampled in 15 day bins (with typically three original points per bin). The cross correlations are shown in Figure 7. Time lags were estimated by fitting a Gaussian function to the DCCF peak. Lag uncertainties were estimated by a Monte Carlo method (bootstrap resampling and adding white noise equal to total standard deviation (see Larsson 2012, for details)). The estimated time lags relative to the gamma ray are $6 \pm 18$ days for the optical and $188 \pm 28$ days for radio. While the gamma-ray-optical correlation is consistent with zero lag, the large 15 GHz time lag is not unusual for blazars (Fuhrmann et al. 2014; Max-Moerbeck et al. 2014). It is also clear that the radio light curve is less regular than the one in gamma rays. This is true for the first part of the OVRO light

<figure><img src="image_9.png" /><figcaption>**Figure 7.** DCCFs between Fermi-LAT ($E > 1$ GeV) and the rebinned optical and radio (15 GHz) light curves. Blue (dashed) curve: gamma vs. optical. Purple (solid) curve: gamma vs. radio. Estimated time lags relative to gamma rays are $6 \pm 18$ days for the optical and $188 \pm 28$ days for the radio.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** OVRO 15 GHz light curve (in black with error bars) together with LAT ($E > 1$ GeV) flux (purple). Inset panels show cross correlations (DCCF) between the light curves of segments (a)–(c). While the overall time lag of the 15 GHz flux relative to gamma ray is on the order of 200 days there are substantial cycle-to-cycle variations. Lag 0 in the DCCF plots is indicated by vertical lines.</figcaption></figure>

curve, segments (a) and (b) in Figure 8. By contrast, the light curve after approximately MJD 57200 (segment (c)) shows a much more clean, quasi-sinusoidal modulation with a visually apparent time lag relative to the gamma-ray flux. This is corroborated by the cross correlation for that part of the data. A fit to the DCCF in the lower right of Figure 8 gives an estimated lag = 226 ± 25 days. The division into time intervals (a)–(c) here is only meant to illustrate what seems like a variation in correlation behavior between the radio and gamma-ray bands.

### 3.5. Coherency

The question of whether the periodicity in PG 1553+113 is quasiperiodic or strictly coherent has important implications for the interpretation and theoretical modeling of this periodicity. The comparison of pulse phases for the two segments made in Section 3.1 implies an absolute phase drift of ≲0.12 over a few cycles. It is also consistent with a strictly coherent periodicity. Since the presently available gamma-ray data cover almost seven pulse cycles it is also possible to look for any significant phase jitters over this epoch. The times of individual pulse peaks were estimated by cross correlating one cycle of Fermi-LAT data centered on the peak with a Fourier pulse. Our timings of the seven peaks were based on the Fermi-LAT ($E > 1$ GeV) light curve and one peak timing from the optical observations preceding the Fermi-LAT data. Combining the optical and gamma-ray timings is motivated by the strong correlation with a small or negligible time lag between the two bands as shown in Section 3.4. A linear fit to the timing of the gamma-ray peaks results in the following ephemeris, where N is the cycle number:

$$T_{\text{peak}} = \text{MJD} 54659(\pm 31) + 771.5(\pm 7) \times N \text{ days.} \quad (1)$$

The value 771.5 ± 7 (i.e., 2.11 yr) is the modulation period. Deviations relative to this ephemeris are shown in Figure 9. The observed phase jitter is slightly larger than expected for a strictly coherent periodicity ($\chi^2 = 11$ for 6 dof, probability = 0.09). We cannot exclude that this is an effect of a slight underestimation of the data errors.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Pulse phase variation relative to a strictly coherent period. The pulse phase for each cycle is estimated by cross correlation with a Fourier model fitted to the full light curve. The first point is from the optical data in Figure 2. All other points are from the $E > 1$ GeV light curve. The cycle number corresponds to N in Equation (1). An observed-computed ($O-C$) diagram compares the observed time of maximum flux O with the calculated time C, assuming a known constant period T, and is a classical method for studies of period changes in variable stars.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** SED of PG 1553+113 in the precursor phase shown in red (gray) circles as highlighted in the right panel of Figure 3 and the main phase shown in blue (dark) circles. A cutoff power-law function was used, shown by the green dashed lines.</figcaption></figure>

### 3.6. Energy Dependence of the Gamma-Ray QPO

The spectral energy distributions (SEDs) in different phase ranges of the periodic variations, were extracted with the binned maximum likelihood method with the same analysis parameters as were adopted for light curve extraction. Each of the spectral bins was obtained using a power-law model for the source of interest with a spectral index fixed to 2. We selected data in two different phase ranges, selecting the substructure in the pulse shape with four Fourier components as shown in Figure 3 (first panel). We estimated the SED for the precursor phase interval of about 0.5–0.7 in Figure 3. Note that these phases are based on a slightly different period than that given in Equation (1). The precursor time intervals are highlighted in Figure 3 (second panel). For comparison we also extracted the SED in the main part of the peak with a complementary time selection (phase interval of about 0.7–1.5). The two SEDs are shown in Figure 10. We fitted the SEDs with a cutoff power law with best-fit spectral index values of 1.62 ± 0.03
($\chi^2_{\text{red}}$ of 1.2) and $1.51 \pm 0.02$ ($\chi^2_{\text{red}}$ of 1.2) and an energy cutoff equal to $84 \pm 23$ GeV and $140 \pm 30$ GeV for the precursor and main peak phases, respectively. The precursor phase shows a softening in the spectral index with respect to the main part of the periodic modulation with a significance of about $3\sigma$.

This behavior is corroborated by the spectral index evolution as a function of time using the 45 day time bins of gamma-ray light curves as reported in Figure 1. A few stochastic flaring activities show a softening of the spectrum during the precursor time frames. Moreover, we also performed periodicity searches on the energy flux light curve, and using this different set of data we found a higher significance as reported in Section 3.2, confirming the hypothesis that the periodic modulation is more prominent at higher energies.

# 4. Discussion

## 4.1. Candidate Astrophysical Scenarios for Periodicity

The fairly long-lived ~2.1 yr period and correlated gamma-ray and optical cyclic modulation in PG 1553+113 have been discovered thanks to the ~15 yr continuous all-sky survey performed by Fermi-LAT, with no degradation of its performance.

Advection-dominated accretion regimes (magnetically dominated/arrested accretion flows: MDAF and MAAF) related to QPOs, could explain turbulent, peculiar radio kinematics related to radiatively inefficient, high-energy, TeV, BL Lac objects like PG 1553+113 (Fragile & Meier 2009; Karouzos et al. 2012; Piner & Edwards 2014). Contrary to this, if PG 1553+113 is assumed to be an intermediate blazar, the stronger accretion power enabling QPOs scenarios, might be explained by a pulsational perturbational and/or helical structure, also precessional (i.e., geometrical) mechanisms. Precession can naturally be produced by a binary SMBH companion, with jet wiggling providing a QPO contribution.

On the other hand, an intrinsic internal helical structure with a single central SMBH (e.g., Camenzind & Krockenberger 1992; Villata & Raiteri 1999; Nakamura & Meier 2004; Rieger 2004) predicts that magnetic fields play a major role in QPOs of blazars (e.g., McKinney et al. 2012). VLBI parsec-scale structures, showing wiggling radio knots (e.g., Piner et al. 2010; Piner & Edwards 2014; Lico et al. 2020), are believed to originate from shocks propagating through cylindrical or conical jets, interacting with a preexisting helical structure. Magnetic field turbulence can lead to recurring Kelvin-Helmholtz instabilities, capable of generating weak shocks (e.g., Lister et al. 2013; Hughes et al. 2015), that are advected into the jet base. The toroidal component of the jet magnetic field is wound up by the accretion flow from the central rotating SMBH. Lense-Thirring (rotational dragging in general relativity) precession (e.g., Wilkins 1972; Bardeen & Petterson 1975; Ingram & Done 2009) of the inner portions of the accretion disk, also can generate QPOs.

The contribution by bulk jet precession, via differential Doppler boosting, is expected to be modest. Changes in direction at the jet nozzle can occur by accretion disk Lense-Thirring precession, or orbital Keplerian motion of the accretion system of a close, <10¹⁸ cm separation, binary SMBHs. Parsec-scale jet curvature, wobbling, precession, rotation, nutation (rocking and nodding), or non-ballistic helical motion of in-jet emission components are all able to produce flux modulations, via cyclical magnification provided

by recurrent Doppler beaming variations. Nozzle and jet base precession and rotation can then produce particle-accelerating, MHD stresses and turbulence. Works following up on A15a have reported some of these mechanisms for PG 1553+113 (Caproni et al. 2017; Cavaliere et al. 2017; Sobacchi et al. 2017; Sandrinelli et al. 2018; Tavani et al. 2018; Cavaliere et al. 2019; Lico et al. 2020).

Among many open astrophysical interpretations, we list five scenarios for the gamma-ray QPO of PG 1553+113:

(a) Pulsational accretion instabilities produce efficiency modulations in the energy outflow. MDAF with subluminal, turbulent, and peculiar radio kinematics (Fragile & Meier 2009; Karouzos et al. 2012; Piner & Edwards 2014) could be explained as a precessing or helical jet (Conway & Murphy 1993), and is also able, periodically, to efficiently impart energy to particles in the jet (Tchekhovskoy et al. 2011). Periodicity here could be too short (~$10^5$ s $\cdot M_{\text{SMBH}}/10^8 M_\odot$) (Honma et al. 1992), but magnetically choked accretion flows can produce longer periods for slow-spinning SMBHs (McKinney et al. 2012).

(b) Geometrical jet precession (e.g., Romero et al. 2000; Stirling et al. 2003; Rieger 2004; Caproni et al. 2013), rotation and nutation (Camenzind & Krockenberger 1992; Vlahakis & Tsinganos 1998; Hardee & Rosen 1999; Valtonen et al. 2008), or an intrinsic helical structure where strong winds and a magnetic field wraps the jet (e.g., Conway & Murphy 1993; Roland et al. 1994; Villata & Raiteri 1999; Nakamura & Meier 2004; Ostorero et al. 2004; Raiteri et al. 2015), all are able to produce a QPO, through the cyclical change of the viewing angle $\theta$ and change of the bulk Doppler boosting factor $D = 1/(\Gamma(1 - \beta \cos \theta))$. Here, $\beta = v/c$, $\Gamma = 1/\sqrt{1 - \beta^2}$. The observed flux $F_\gamma \propto E^{-\alpha}$ is proportional to $D^{3-\alpha}$, for a boosted high-energy emission blob in the jet. In PG 1553+113 a variation $\Delta D(t) = \Gamma^{-1}(1 - \beta(t)\cos\theta(t))^{-1}$ (Rieger 2004) greater than 40% with precession angle ~1° is required to explain the ~2.8 flux modulation factor seen in the gamma-ray light curves obtained by the LAT. On the other hand, $D$ asymptotically diverges for small $\beta \cos\theta$, and $\Delta D(t)$ changes substantially if the wobbling spans very small angles. In addition, concurrent intrinsic outflow and efficiency variations (point (a)) could help to relax this constraint.

(c) A mechanism analogous to a low-frequency QPO from high-mass X-ray binaries is supported by the many BH timescales that depend inversely on the mass (Fender & Belloni 2004; King et al. 2013), where the accretion-outflow coupling is the basis of the periodicity. The microquasar mechanism of Lense-Thirring precession of the jet's nozzle (Wilkins 1972) requires that the inner accretion flow forms a geometrically thick and viscous torus. A standard thin disk instead warps (Bardeen-Petterson effect, Bardeen & Petterson 1975) rather than precesses (Ingram & Done 2009). In the case of a binary SMBH system, the primary minidisk precesses with nutation by the gravitational influence of the secondary, changing jet direction, viewing angle, and Doppler boosting factor $D$. It also changes the variability timescales ($\Delta t/\Delta t' = 1/D$, but these usually act on much longer timescales than the 2 yr period).

(d) Eccentric orbiting massive stars (or an intermediate-mass black hole (IMBH)) with inclined/polar revolutions around a single SMBH. This should be very common and could provide periodic modulations of the accretion flow as introduced in Section 4.4.
(e) A gravitationally bound binary SMBH system (Begelman et al. 1980; Barnes & Hernquist 1992) with a total mass $5.8 \times 10^8 M_\odot$ (Dhiman et al. 2021), and a milliparsec separation in the early inspiral GW-driven regime, as introduced in A15a. Most binary SMBHs will spend the majority of their lifetime between 0.01 and 1 pc separations, in an intermediate phase of evolution between scattering any stars in the nuclear region and gravitational radiation dominance (Sesana & Khan 2015). Keplerian binary orbital motion would induce periodic accretion perturbations (Bogdanović et al. 2008; Valtonen et al. 2008; Pihajoki et al. 2013; Liu et al. 2015; Ciprini 2018) or jet nutation due to 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. 2015a).

Significant acceleration of the disk evolution and accretion onto a binary SMBH system, where there is excited eccentricity in the inner region of the circumbinary disk, can create an overdense lump giving rise to enhanced periodicity in the accretion rate (Roedig et al. 2012; Nixon et al. 2013; Farris et al. 2014; Gold et al. 2014; Doğan et al. 2015).

Binary SMBH-induced periodicities have timescales ranging from ~1 to ~25 yr (Komossa 2006; Rieger 2007). The total mass of the SMBH in PG 1553+113 can also be estimated by the putative link between inflow/accretion (disk luminosity) and outflow (jet power; Ghisellini et al. 2014) with a 0.1 $dM_{Edd}/dt$ rate and $D = 30$. Comparing luminosities $L_{disk}$, $L_{rad}$, and jet radiative power $P_{rad}$, we obtain $M_{SMBH} \approx 1.6 \times 10^8 M_\odot$, in agreement with Woo et al. (2005), but a factor of 3.5 lower than that in Dhiman et al. (2021). The observed (2.1 ± 0.2) yr period is equivalent to an intrinsic orbital time $T'_{Kep} \le T_{obs}/(1+z) \approx 1.5$ yr. Considering a Keplerian orbit for the bound binary SMBH gives a binary size of 0.005 pc $\approx 100 R_s$, with $R_s = 2G(M_1)/c^2$ the Schwarzschild radius and $M_1$ the mass of the main SMBH. 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} = (5c^5a^4)/(256qG^3M^3) \approx 10^5-10^6$ yr (with $q = M_1M_2/(M_1+M_2)^2 \sim 0.1-0.01$) might be too small. In Section 4.2, we briefly evaluate the gravitational lifetime and reliability of the binary SMBH hypothesis, while in Section 4.3, we consider the consistency with a binary SMBH population and the gravitational radiation background.

## 4.2. Binary SMBH Model Reliability

In principle, a binary system of two accreting and jetted SMBHs with $M_1 = 1.6 \times 10^8 M_\odot$ and mass ratio in the range $q = M_2/M_1 = 0.01-0.1$, can induce the coherent, high-energy/multifrequency, short (~2.1 yr) periodicity, with $T'_{Kep} \approx 1.5$ yr in the source rest frame, as observed in PG 1553+113. The lifetime of the system, which is directly proportional to its probability of being observed, can be evaluated as the standard GW merger timescale $t_{GW}$, as defined in the previous section, where circular orbits are here considered, and separation $a$ depending on the characteristic orbital period $T$ of each of the three hypotheses below.

(1) The periodicity of the binary SMBH causes a periodic accretion supply of gas in the minidisk of the primary SMBH, which is the direct source of the variability of the relativistic jet emission through quadrupolar torque on the large circumbinary disk (e.g., Rieger 2004; Fragile et al. 2007; Roedig et al. 2012; Farris et al. 2014; Tang et al. 2017; Kelley et al. 2019). This

binary Keplerian orbital period $T_k$ dictates the observed periodicity and $a = (T_k^2 GM/(4\pi^2))^{1/3}$ yielding the lifetime $t_{GW} \sim 2 \times 10^6$ yr for mass ratio $q = M_2/M_1 = 0.01$ ($2 \times 10^5$ yr for $q = 0.1$), in agreement with our estimates.

(2) Considering the jet aligned with the primary BH spin, the secondary jet precesses around the total angular momentum of the binary spin-orbit (SO) coupling, with the same periodicity. This SO precession has $a = (T_{SO}q(4+3q)GM_1\sqrt{G(M_1+M_2)}/(4\pi(1+q)c^2))^{2/5}$ (Apostolatos et al. 1994), yielding a too-short lifetime of $t_{GW} \sim 0.8$ yr for $q = 0.1$ (0.2 yr for $q = 0.01$); such a system would certainly not persist 15 yr.

(3) If the primary jet, perpendicular to the SMBH accretion minidisk, is misaligned with the binary orbital plane, the minidisk precesses due to the torque induced by the secondary SMBH with a timescale depending on its distance from the primary SMBH. Torques that tend to align the minidisk with the orbital plane, compete with Bardeen-Petterson torques that tend to align the minidisk with the SMBH spin. With the jet launched by a magnetic funnel anchored to the inner edge of the disk, the relevant distance is the innermost stable circular orbit (ISCO). With the precession of the primary accretion disk induced by the secondary SMBH, similar to the second case with SO precession periodicity observed $T_{SO}$, now for the inner ring of the accretion disk, around the ISCO and assumed as the jet launcher, we have the disk semimajor axis: $a_d = ((3T_dq(6GM_1/c^2))^{3/2}\sqrt{G(M_1+M_2)} \cos i)/(8\pi\sqrt{1+q})^{1/3}$ (Miller & Krolik 2013), where $i$ is the inclination of the disk at the ISCO (for $i \sim \pi/4$ or smaller yields similar results depending only weakly on $i$), and $T_d$ the precession period of the inner ring of the accretion disk. Again, such a system would certainly not live 15 yr ($t_{GW} \sim 0.6$ yr for $q = 0.1$ and 0.3 yr for $q = 0.01$).

The 1.5 yr rest-frame periodicity cannot therefore be ascribed to precession timescales induced by the binary companion, either on the primary SMBH spin or on its accretion disk because of the too-short binary merging timescale ($\lesssim 1$ yr). If we postulate a binary SMBH system in PG 1553+113, the only reasonable mechanism is a Keplerian periodic orbital modulation of the accretion rate onto the primary SMBH induced by the companion, which dominates the hydrodynamic periodic variability with a recurrent fueling of the primary SMBH accretion disk and the activity of its relativistic jet. A tentative graphical toy model illustration of accretion geometries for such a binary SMBH system is shown in Figure 11.

A periodic perturbation of the primary ($\#1, M_1 \approx 10^8 M_\odot$) SMBH's minidisk and accretion power, feeding the primary relativistic jet responsible for gamma-ray emission, is provided through an intrinsically oscillating, gravitational quadrupolar torque of the large and very dense circumbinary disk that couples with the two SMBHs' dynamics. The cavity gap of the circumbinary disk has a radius roughly twice the binary separation. Periodic accretion streams ($\dot{M}_1, \dot{M}_2$), extending from the circumbinary disk, orbiting at longer periods than the circum-single minidisks, feed the two BHs. The relatively high mass ratio (e.g., $q = M_2/M_1 \approx 0.01$ or $q = 0.1$) also allows the secondary, \#2 MBH ($\#2, M_2 \approx 10^6 M_\odot$, or $10^7 M_\odot$) to reach the vicinity of the circumbinary annulus, with a consequent larger accretion flow and greater gas received ($\dot{M}_2 > \dot{M}_1$), causing an additional periodic perturbation in accretion power and disk tearing of \#1, during orbit of \#2 around the center of mass (the + cross in the drawing center). The gravitational
<figure><img src="image_13.png" /></figure>

**Figure 11.** Cartoon of the hydrodynamical accretion geometries of a binary SMBH theoretical scenario for PG 1553+113. Around each SMBH a “mini circum-single” disk is fed by variable accretion streams ($M_1, M_2$) extending from a large and dense circumbinary disk, having a gap radius of about twice the binary separation. The relatively high mass ratio (e.g., $q = M_2/M_1 \simeq 0.01$ or $q = 0.1$) also allows the secondary, #2 SMBH ($\#2, M_2 \simeq 10^5 M_\odot$, or $10^6 M_\odot$) to reach the vicinity of the circumbinary annulus, tending to receive a disproportionate share of the accretion rate ($\dot{M}_2 > \dot{M}_1$), causing an additional periodic perturbation. The circumbinary disk also has intrinsic oscillations by gravitational quadrupolar torque and orbits at longer periods than the two minidisks feed, producing further cyclic variations in the accretion rate, jet matter, and energy feeding.

lifetime for GW radiation losses of such a binary system geometry would be $t_{GW} \sim 2 \times 10^6$ yr for $q = 0.01$, and $2 \times 10^5$ yr for $q = 0.1$.

### 4.3. Pulsar Timing Array Constraints on Close Binary SMBHs

It is interesting to frame our findings for PG 1553+113 within the recent results presented by the pulsar timing array (PTA) experiments ongoing around the world. In fact, evidence for a common red signal, with properties consistent with a GW origin has been found in the latest data releases of the European PTA (Antoniadis et al. 2023), NANOGrav (Agazie et al. 2023), the Parkes PTA (Reardon et al. 2023) and the Chinese PTA (Xu et al. 2023). The significance of the detection is between $2\sigma$ and $4\sigma$ depending on the data set, and the amplitudes of the signals are consistent with a stochastic gravitational-wave background (GWB) produced by a cosmic population of SMBH binary (SMBHB) systems, although an early-Universe origin cannot be excluded (Afzal et al. 2023; Antoniadis et al. 2024). Assuming an SMBHB origin of the signal, Holgado et al. (2018) constructed a theoretical framework to translate its amplitude into the fraction of expected periodic sources observed by Fermi-LAT. The details of the model can be found in their paper, but the main result is shown in their Figure 4. This figure shows that, under the assumption of a GWB signal amplitude of $A_{1\,yr} = 1.45 \times 10^{-15}$ (the published limit at the time), the expected fraction of Fermi-LAT BL Lacs showing flux oscillations with a period less than 5 yr is about $3 \times 10^{-4}$. The signal now observed in the PTA data has an

amplitude of about $A_{1\,yr} = 2.5 \times 10^{-15}$. Since the amplitude scales with the square root of the number of SMBHBs contributing to the signal, this implies an expected fraction of periodic BL Lacs of about $10^{-3}$.

This is consistent with finding one confirmed periodic source (e.g., PG 1553+113) out of thousands of sources in the Fermi-LAT catalog. It should be noted, however, that PG 1553+113 is among the brightest sources in the Fermi-LAT sample. Moreover, periodicity (albeit at a lower significance) has been suggested for PKS 2155-304 and BL Lac (Sandrinelli et al. 2018) and for a few other systems (Prokhorov & Mora-ghan 2017; Peñil et al. 2020, 2024a, 2024b), all of which are in the bright flux tail of the Fermi-LAT source distribution (see Figure 4 in Sandrinelli et al. 2018). So it might be that the detection of periodicity in Fermi-LAT data is flux limited, and the identification of PG 1553+113 (plus perhaps a few other systems) as the only periodic Fermi-LAT sources, is due to selection effects. Therefore, no firm indication either in favor or against PG 1553+113 being a binary SMBH system can be drawn from the GW signal observed in the PTA data.

### 4.4. Perturbation by Inclined Orbit Star-like Objects

Eccentric orbiting star-like objects, or IMBH around a single SMBH should be common. With an orbital binding energy $-3GM^2/(5R) \simeq$ some MeV per nucleon (stellar interiors kinetic energies of nuclei have a Maxwellian internal energy distribution with $\sim$ keV averaged value), the Roche density of the star $\rho_R = 3.5M/r^3$ required to resist tidal heating and disruption is very small for the intrinsic gamma-ray 1.5 yr periodicity of PG 1553+113 in the case of nonpolar orbits. Highly inclined stellar orbits of $\sim 10^{16}$ cm radius can be perturbed more easily than an MBH/SMBH orbit, with a gravitational energy loss and orbit coalignment with the disk timescale of $> 10^9$ yr ($10^3 \times$ the stellar interaction time; for example, Dai & Blandford 2013; Linial & Metzger 2023). Local fire bubbles ($\sim 10^5$ K thermal bremsstrahlung) or magnetic reconnection flare avalanches in the disk and corona are expected from disk piercing, twice, by the star (or white dwarf, neutron star, or star-sized BH) with energy release $\Delta E = \frac{1}{2}(\Sigma\sigma^2v_r^2)$; where $\Sigma$ is the disk surface density, $\sigma^2$ the disk-star cross section, and $v_r$ the hypersonic relative velocity of the star to the disk matter. The disk serves not only as the source of the outflow material, but it also illuminates it and reflects the outflow radiation back. Accretion flow perturbations propagating to the inner disk are the result of double disk-star impact fireballs, modulating the accretion rate, nozzle, and wind-confining relativistic jet of a blazar like PG 1553+113 (Ivanov et al. 1998). General relativity MHD simulations, with an orbiting IMBH ($10^2-10^5 M_\odot$) companion, indicate mildly relativistic particle acceleration in the stable poloidal magnetic field of the jet funnel, and ADAF ejecting outflow clumps (Kejriwal et al. 2024; Pasham et al. 2024). High eccentricity of the IMBH's orbit around the SMBH, can also explain variations from a strict periodicity.

## 5. Conclusions

The primary aim of this paper is to follow up on the tentative identification of a $2.1 \pm 0.2$ yr periodic oscillation in the blazar PG 1553+113. The additional data used in this study provide further and independent support for the existence of the
periodic flux modulation, which was described in the initial study in A15a.

The Fermi-LAT gamma-ray data show a persistent 2.1 yr oscillation throughout the 15 yr of observations. Compared to the A15a analysis, the QPO significance has increased to $4\sigma$ with respect to simulations using two different types of red noise models.

An essential benefit of the longer time series is that it allows us to compare oscillation properties in two independent parts of the data. Both sections of the light curve are found to contain an oscillation with a well-determined period and phase. The chance probability for the two oscillations to be so close in period and phase is less than 0.01.

A multiwavelength analysis is used to study correlations between the gamma-ray flux and the X-ray, optical, and radio measurements. No significant gamma-ray flux correlation is found with X-rays, while both optical and OVRO 15 GHz radio show strong correlations with time lags of $6 \pm 18$ and $188 \pm 28$ days, respectively, relative to the Fermi-LAT ($E > 1$ GeV) gamma rays. Pulse timing over eight gamma-ray periodic cycles shows arrival time variations that are consistent with a strictly coherent period. Additional years of long-term and regular flux monitoring of PG 1553+113 are needed, especially at X-ray and VHE gamma-ray ($E > 100$ GeV) energy bands, possibly with a regular radio/optical polarization monitor.

A spectral gamma-ray analysis reveals an energy dependence of the QPO such that its amplitude is greater above 1 GeV than at lower energies. We also find hints for softening during a precursor feature in the rising part of the oscillation pulse.

The binary SMBH interpretation, with $1.6 \times 10^8 M_{\odot}$ total mass, $\sim 100 R_s$ separation and lifetime (residence time $\sim |a/(da/dt)| \sim 10^5$, $10^6$ yr, where $a$ is the semimajor axis of the binary) has been discussed in Section 4, although at least four additional astrophysical scenarios, including precessing, nutation, or helical and wiggling jet structures, MHD instabilities, and perturbations by inclined massive stellar objects or IMBHs, are possible. The $\sim 1.5$ yr rest-frame QPO of PG 1553+113 cannot be ascribed to precession timescales (binary SO coupling or spin/accretion disk precession of the primary SMBH induced by the lower-mass companion), because of the too-short ($\lesssim 1$ yr) merging scale with these mechanisms. The optical gamma-ray QPO can be caused by Keplerian periodic orbital modulation of the supplied gas and matter in the main jet by the primary minidisk of the SMBH, with the contribution also of a substantial quadrupolar torque of the large and dense circumbinary disk. The QPO might be a hydrodynamic periodic variability of the fueling of the primary SMBH accretion power, where the toroidal component of the jet magnetic field is wound up by the accretion flow, and magnetorotational stresses are induced in the jet, energizing electron populations and producing an underlying periodic modulation in gamma-ray flux enhancements. This would survive for the rather long gravitational losses lifetime, prior to the disruptive coalescence of $t_{GW} \sim 2 \times 10^6$ yr for $q=0.01$, and $2 \times 10^5$ yr for $q=0.1$. The case of PG 1553+113 (chirp mass $M_8 \lesssim 1$), representing the strongest case of significant QPO out of a sample of $\sim 1000$ Fermi-LAT bright gamma-ray BL Lac objects, could suggest a possible contribution to the GW stochastic background probed by international PTAs, from the entire object class of BL Lacs and their misaligned AGN parent population. Instabilities and perturbations are easier in the radiatively inefficient regimes of BL Lac objects, and their

supposed longer lifetimes $t_{GW}$ in close orbits, with respect to those of heavier and more distant quasars and FSRQs. This also supports a GW stochastic contribution, which could be determined in the future by the Square Kilometre Array. There is already a candidate detection of 4–5 nHz continuous GW emission generated by binary SMBHs in the local Universe (Antoniadis et al. 2024).

Minor mergers, i.e., $M_2 \lesssim M_1/4$, are more likely to be observed electromagnetically because the time to merge is $t_{GR} \propto M_2^{-1}M_1^{-2/3}$ and therefore these should be more frequent events. In this view, perturbations by a highly eccentric and inclined polar-orbiting massive star or an IMBH are interesting to consider for the case of PG 1553+113.

QPO and related astrophysical models for PG 1553+113 could be evaluated in the future by observations of helical or wiggling and wobbling patterns in the parsec-scale radio jet, by the identification of regular radio components and polarization patterns and by double-peaked spectroscopic emission lines. The recent X-ray observations of PG 1553+113 by IXPE and simultaneous multifrequency data in a range of $\sim 8$ days in 2023 February (Middei et al. 2023) point out turbulence as a contributor to the jet magnetic field that, eventually, may produce instabilities generating weak shocks advected into the jet base. TESS optical observations over a range of 24 days in 2022 April–May could add more information about jet turbulence.

The unique Fermi-LAT all-sky and time-domain survey will continue, providing essential data for PG 1553+113 and triggering the scientific interest in this potential multimessenger science case, while also promoting efforts in multifrequency follow-up observations.

## Acknowledgments

We thank the anonymous referee for the useful comments. The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration (NASA) and the Department of Energy (DoE) in the United States, the Commissariat à l'Energie Atomique (CEA) and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules (CNRS IN2P3) in France, the Agenzia Spaziale Italiana (ASI) and the Istituto Nazionale di Fisica Nucleare (INFN) in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council (Vetenskapsrådet, VR) and the Swedish National Space Board (SNSB) in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica (INAF) in Italy and the Centre National d'Études Spatiales (CNES) in France. This work was performed in part under DOE Contract DE-AC02-76SF00515.

The NASA satellite Neil Gehrels Swift Gamma-Ray Burst Explorer is a MIDEX Gamma-Ray Burst mission led by NASA with the participation of Italy and the UK. The NASA satellite RXTE was a mission managed and controlled by NASA's Goddard Space Flight Center (GSFC) in Greenbelt, Maryland. The RXTE ASM light curves project was possible thanks to
RXTE teams at MIT and GSFC, and supported by NASA contract No.NAS5-30612.

This research has made use of data from the OVRO 40 m monitoring program (Richards et al. 2011), which is supported in part by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G 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. This work has made use of data obtained by the Tuorla optical blazar monitoring program (mainly the KVA Observatory and the 1.03 m telescope at the Tuorla Observatory), in this case all included in paper A15a. This work has made use of data obtained with the 0.76 m Katzman Automatic Imaging Telescope (KAIT), supported by donations from corporations, foundations, the Lick Observatory, and the US NSF. This work has made use of data obtained with the Catalina Real-Time Transient Survey (CRTS), supported by the US NSF under grant Nos. AST-0909182 and AST-1313422, and the CSS funded by NASA grant No. NNG05GF22G. This work has made use of data obtained from the ASAS-SN project, supported by a foundation grant GBMF5490 to the Ohio State University and NSF grants AST-1515927 and AST-1908570. The DASCH project at Harvard is grateful for partial support from US NSF grants AST-0407380, AST-0909073, and AST-1313370.

This research has made use of data and/or software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC. This research has made use of the Smithsonian/NASA's ADS bibliographic database. This research has made use of the NASA/IPAC Extragalactic Database (NED) database, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA. This research has made use of archival data, software, and/or online services provided by the

Space Science Data Center (SSDC) operated by ASI. This research has made use of the XRT Data Analysis Software (XRTDAS) developed under the responsibility of the SSDC. This research has made use of the SIMBAD database, CDS, Strasbourg Astronomical Observatory, France.

This work was supported by a grant from the Simons Foundation (00001470, R.B.). S.Ci. and INFN personnel in SSDC also performed under contract agreement ASI-INFN 2021-43-HH.0. S.G. acknowledges the Italian MUR, program "Dipartimenti di Eccellenza 2018-2022" (Grant SUPER-C).

*Facilities:* Fermi, Swift, OVRO:40m, RXTE, KAIT, SO:1.5m, SO:1m, SO:Schmidt, ASAS-SN.

## Appendix A

### The 3 Day and Weekly Bin Public LAT Light Curves

The evidence of a potential ~2.1 yr periodic gamma-ray flux modulation is independent of the bin value adopted to calculate the Fermi-LAT light curve of PG 1553+113 as an example here of the analysis of the public Fermi-LAT light-curve data of PG 1553+113 from 2008 August to 2023 November, with a fine, 3 day binning and publicly available at the Fermi-LAT light-curve repository (LCR) of the FSSC<sup>74</sup> (Fermi Large Area Telescope Collaboration 2021; Abdollahi et al. 2023) is shown in Figure 12. In addition, other techniques are able to highlight such a signal of periodic oscillation, including the phase dispersion minimization, the discrete autocorrelation function (DACF), and the structure function (Abdo et al. 2010). The LSP spectrogram and DACF are included in the right panel of Figure 12. The LSP power spectrum shows a break and flattening at about 10 × the light-curve bin size (34 days), and a signal power-law index $\alpha = 0.8$ ($1/f^\alpha$ spectrum, $f = 1/t$) for longer timescales. The power peak at 2.13 ± 0.22 yr is also evident. *T* in the main plot and inset (DACF) represents the ~2.1 yr timescale. The distribution of the LSP signal peaks for $10^6$ simulated light curves corresponding to properties of the 3 day bin light curve of Figure 12, is calculated and reported in Figure 13.

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Left panel: 1 week bin and 3 day bin, science-ready and public Fermi-LAT light curve of PG 1553+113 from 2008 August 8 to 2023 November 16, i.e., from 15 yr of Fermi mission survey operations, extracted from the LAT light curves repository (LCR) at the Fermi Science Support Center (GSFC-NASA). Right panel: corresponding periodogram power density spectrum (main panel) from the 3 day bin light curve, and DACF in the inset, with the periodic component *T* = 2.13 yr highlighted. This period is consistent with the findings reported in the main text of this paper.</figcaption></figure>

<footnote n="74"><sup>74</sup> https://fermi.gsfc.nasa.gov/ssc/data/access/lat/LightCurveRepository/source.html?source_name=4FGL_J1555.7+1111</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Distribution of the LSP signal peaks for $10^6$ simulated light curves, respectively for the public 3 day bin LCR light curve shown in Figure 12 in the main panel, and for the detrended version of it in the inset panel. Dashed lines are the LSP power for the two, true data, light curves. Significance of the period here is $4\sigma$ and $>4\sigma$ for the 3 day bin LCR light curve and the detrended version of it, respectively.</figcaption></figure>

## Appendix B
1996–2011 RXTE ASM Light Curve and 1912–1988 Historical Optical Light Curve

A 1 day bin averaged count rate X-ray (1.5–12 keV) light curve was calculated from the ASM archive$^{75}$ (Levine et al.

1996) of RXTE, from 1996 January 6 to 2011 September 23 (Figure 14). Only spurious signal power peaks (~1 yr, ~1/2 yr, and ~50 days) are found (Figure 14).

PG 1553+113 is an optically bright blazar (range $15 \lesssim R \lesssim 13$ mag), with non-negligible chances for detection in optical observatory plates, despite the fact that it is not very close to relevant, historically intensively observed, M/NGC/IC catalog objects in the field. In Figure 14, we present a first version of the PG 1553+113 secular optical light curve (1912–1988) extracted from the “Digital Access to a Sky Century at Harvard” (DASCH) project for plate scanning and digitization, and based on 450000 plates (years from 1890 to 1990) of the Harvard College Observatory (Grindlay et al. 2012). R-band flux densities (1912–1988) from DASH DR5 and magnitudes calibrated with APASS AAVSO Photometry All-Sky Survey Rel. 8 (B-band mag) and ATLAS All-Sky Stellar Reference Catalog v2 ($G \simeq V$ band mag), are extrapolated using fixed color indexes $B-R=0.55$ and $B-V=0.25$. The various periodicity analysis methods are not successful. Future searches for plates and scanning projects around the world would likely help to improve this large-gapped and sparse historical time series.

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Top panel: the Rossi XTE ASM 1.5–12 keV X-ray count rate light curve with 1 day bin averaged data points from 1996 to 2011 (a range of 15.7 yr) with negative rates and upper limits removed. In the inset panel the corresponding LSP showing only spurious 1 yr periodicity owing to data gaps due to the Sun vicinity with PG 1553+113, its 1/2 yr harmonic and the 50 day scale of the typical Sun gap duration and the spacecraft precession. Bottom panel: secular historical optical flux light curve of PG 1553+113 with data points from 1912 to 2020 obtained using DASCH DR5 with all data scaled with constant color indexes to approximately the R-band flux densities. A strict 2.1 yr sinusoidal line is added as background. In the inset panel the LSP for the 1910–1990 historical and highly gapped portion of the light curve is shown, highlighting only low-frequency noise and the usual 1 yr timescale (mild peaks at 371 and 310 days), common in ground-based optical light curves. The data after 2005 were already presented in Sections 2.4 and 2.3.</figcaption></figure>

<footnote>$^{75}$ https://heasarc.gsfc.nasa.gov/docs/xte/</footnote>
ORCID iDs

S. Abdollahi https://orcid.org/0000-0002-6803-3605

L. Baldini https://orcid.org/0000-0002-9785-7726

R. Bellazzini https://orcid.org/0000-0002-2469-7063

B. Berenji https://orcid.org/0000-0002-4551-772X

E. Bissaldi https://orcid.org/0000-0001-9935-8106

R. D. Blandford https://orcid.org/0000-0002-1854-5506

R. Bonino https://orcid.org/0000-0002-4264-1215

P. Bruel https://orcid.org/0000-0002-9032-7941

S. Buson https://orcid.org/0000-0002-3308-324X

R. A. Cameron https://orcid.org/0000-0003-0942-2747

P. A. Caraveo https://orcid.org/0000-0003-2478-8018

F. Casaburo https://orcid.org/0000-0002-2260-9322

E. Cavazzuti https://orcid.org/0000-0001-7150-9638

C. C. Cheung https://orcid.org/0000-0002-4377-0174

G. Chiaro https://orcid.org/0000-0001-9328-6439

S. Ciprini https://orcid.org/0000-0002-0712-2479

G. Cozzolongo https://orcid.org/0009-0001-3324-0292

P. Cristarella Orestano https://orcid.org/0000-0003-3219-608X

S. Cutini https://orcid.org/0000-0002-1271-2924

F. D'Ammando https://orcid.org/0000-0001-7618-7527

N. Di Lalla https://orcid.org/0000-0002-7574-1298

F. Dirirsa https://orcid.org/0000-0002-3909-6711

L. Di Venere https://orcid.org/0000-0003-0703-824X

A. Domínguez https://orcid.org/0000-0002-3433-4610

S. J. Fegan https://orcid.org/0000-0002-9978-2510

E. C. Ferrara https://orcid.org/0000-0001-7828-7708

A. Fiori https://orcid.org/0000-0003-3174-0688

Y. Fukazawa https://orcid.org/0000-0002-0921-8837

S. Funk https://orcid.org/0000-0002-2012-0080

P. Fusco https://orcid.org/0000-0002-9383-2425

F. Gargano https://orcid.org/0000-0002-5055-6395

S. Garrappa https://orcid.org/0000-0003-2403-4582

D. Gasparrini https://orcid.org/0000-0002-5064-9495

S. Germani https://orcid.org/0000-0002-2233-6811

N. Giglietto https://orcid.org/0000-0002-9021-2888

F. Giordano https://orcid.org/0000-0002-8651-2394

M. Giroletti https://orcid.org/0000-0002-8657-8852

D. Green https://orcid.org/0000-0003-0768-2203

I. A. Grenier https://orcid.org/0000-0003-3274-674X

S. Guiriec https://orcid.org/0000-0001-5780-8770

E. Hays https://orcid.org/0000-0002-8172-593X

D. Horan https://orcid.org/0000-0001-5574-2579

M. Kuss https://orcid.org/0000-0003-1212-9998

S. Larsson https://orcid.org/0000-0003-0716-107X

M. Laurenti https://orcid.org/0000-0001-5762-6360

J. Li https://orcid.org/0000-0003-1720-9727

I. Liodakis https://orcid.org/0000-0001-9200-4006

F. Longo https://orcid.org/0000-0003-2501-2270

F. Loparco https://orcid.org/0000-0002-1173-5673

B. Lott https://orcid.org/0000-0003-2186-9242

M. N. Lovellette https://orcid.org/0000-0002-0332-5113

P. Lubrano https://orcid.org/0000-0003-0221-4806

S. Maldera https://orcid.org/0000-0002-0698-4421

D. Malyshev https://orcid.org/0000-0002-9102-4854

A. Manfreda https://orcid.org/0000-0002-0998-4953

L. Marcotulli https://orcid.org/0000-0002-8472-3649

G. Martí-Devesa https://orcid.org/0000-0003-0766-6473

M. N. Mazziotta https://orcid.org/0000-0001-9325-4672

I. Mereu https://orcid.org/0000-0003-0219-4534

P. F. Michelson https://orcid.org/0000-0002-1321-5620

W. Mitthumsiri https://orcid.org/0000-0002-3776-072X

T. Mizuno https://orcid.org/0000-0001-7263-0296

M. E. Monzani https://orcid.org/0000-0002-8254-5308

A. Morselli https://orcid.org/0000-0002-7704-9553

I. V. Moskalenko https://orcid.org/0000-0001-6141-458X

M. Negro https://orcid.org/0000-0002-6548-5622

N. Omodei https://orcid.org/0000-0002-5448-7577

M. Orienti https://orcid.org/0000-0003-4470-7094

E. Orlando https://orcid.org/0000-0001-6406-9910

J. F. Ormes https://orcid.org/0000-0002-7220-6409

D. Paneque https://orcid.org/0000-0002-2830-0502

M. Perri https://orcid.org/0000-0003-3613-4409

M. Persic https://orcid.org/0000-0003-1853-4900

M. Pesce-Rollins https://orcid.org/0000-0003-1790-8018

T. A. Porter https://orcid.org/0000-0002-2621-4440

G. Principe https://orcid.org/0000-0003-0406-7387

S. Rainò https://orcid.org/0000-0002-9181-0345

R. Rando https://orcid.org/0000-0001-6992-818X

B. Rani https://orcid.org/0000-0001-5711-084X

M. Razzano https://orcid.org/0000-0003-4825-1629

A. Reimer https://orcid.org/0000-0001-8604-7077

O. Reimer https://orcid.org/0000-0001-6953-1385

P. M. Saz Parkinson https://orcid.org/0000-0001-6566-1246

L. Scotton https://orcid.org/0000-0002-0602-0235

D. Serini https://orcid.org/0000-0002-9754-6530

A. Sesana https://orcid.org/0000-0003-4961-1606

C. Sgrò https://orcid.org/0000-0001-5676-6214

E. J. Siskind https://orcid.org/0000-0002-2872-2553

G. Spandre https://orcid.org/0000-0003-0802-3453

P. Spinelli https://orcid.org/0000-0001-6688-8864

D. J. Suson https://orcid.org/0000-0003-2911-2025

H. Tajima https://orcid.org/0000-0002-1721-7252

M. N. Takahashi https://orcid.org/0000-0002-0574-6018

D. Tak https://orcid.org/0000-0002-9852-2469

D. J. Thompson https://orcid.org/0000-0001-5217-9135

D. F. Torres https://orcid.org/0000-0002-1522-9065

J. Valverde https://orcid.org/0000-0002-8090-6528

F. Verrecchia https://orcid.org/0000-0003-3455-5082

G. Zaharijas https://orcid.org/0000-0001-8484-7791

References

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 722, 520

Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33

Abdollahi, S., Ajello, M., Baldini, L., et al. 2023, ApJS, 265, 31

Abramowski, A., Aharonian, F., Ait Benkhali, F., et al. 2015, ApJ, 802, 65

Ackermann, M., Ajello, M., Albert, A., et al. 2015a, ApJL, 813, L41

Ackermann, M., Ajello, M., Atwood, W. B., et al. 2015b, ApJ, 810, 14

Afzal, A., Agazie, G., Anumarlapudi, A., et al. 2023, ApJL, 951, L11

Agarwal, A., Mihov, B., Andruchow, I., et al. 2021, A&A, 645, A137

Agarwal, A., Mihov, B., Andruchow, I., et al. 2022, JApA, 43, 9

Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023, ApJL, 951, L8

Ait Benkhali, F., Hofmann, W., Rieger, F. M., et al. 2020, A&A, 634, A120

Ajello, M., Angioni, R., Axelsson, M., et al. 2020, ApJ, 892, 105

Aniello, T., Antonelli, L. A., Tombesi, F., et al. 2024, A&A, 686, A300

Antoniadis, J., Arumugam, P., Arumugam, S., et al. 2023, A&A, 678, A50

Antoniadis, J., Arumugam, P., Arumugam, S., et al. 2024, A&A, in press

Apostolatos, T. A., Cutler, C., Sussman, G. J., et al. 1994, PhRvD, 49, 6274

Atwood, W., Albert, A., Baldini, L., et al. 2013, arXiv:1303.3514

Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071

Bardeen, J. M., & Petterson, J. A. 1975, ApJL, 195, L65

Barnes, J. E., & Hernquist, L. 1992, ARA&A, 30, 705

Beasley, A. J., Gordon, D., Peck, A. B., et al. 2002, ApJS, 141, 13

Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Natur, 287, 307

Bogdanović, T., Smith, B. D., Sigurdsson, S., et al. 2008, ApJS, 174, 455
Breeveld, A. A., Curran, P. A., Hoversten, E. A., et al. 2010, MNRAS, 406, 1687

Camenzind, M., & Krockenberger, M. 1992, A&A, 255, 59

Caproni, A., Abraham, Z., & Monteiro, H. 2013, MNRAS, 428, 280

Caproni, A., Abraham, Z., Motter, J. C., et al. 2017, ApJL, 851, L39

Cavaliere, A., Tavani, M., Munar-Adrover, P., et al. 2019, ApJL, 875, L22

Cavaliere, A., Tavani, M., & Vittorini, V. 2017, ApJ, 836, 220

Chen, Y., Yi, T., Chen, J., et al. 2024, NewA, 108, 102186

Cipriani, S. 2018, in Accretion Processes in Cosmic Sources - II (Trieste: Proc. of Science), 41

Conway, J. E., & Murphy, D. W. 1993, ApJ, 411, 89

Corbet, R. H. D., Markwardt, C. B., & Tueller, J. 2007, ApJ, 655, 458

Covino, S., Landoni, M., Sandrinelli, A., et al. 2020, ApJ, 895, 122

Dai, L., & Blandford, R. 2013, MNRAS, 434, 2948

Dhiman, V., Gupta, A. C., Gaur, H., et al. 2021, MNRAS, 506, 1198

Do, T., Ghez, A. M., Morris, M. R., et al. 2009, ApJ, 691, 1021

Dorigo Jones, J., Johnson, S. D., Muzahid, S., et al. 2022, MNRAS, 509, 4330

Doğan, S., Nixon, C., King, A., & Price, D. J. 2015, MNRAS, 449, 1251

Emmanoulopoulos, D., McHardy, I. M., & Papadakis, I. E. 2013, MNRAS, 433, 907

Falomo, R., & Treves, A. 1990, PASP, 102, 1120

Farris, B. D., Duffell, P., MacFadyen, A. I., et al. 2014, ApJ, 783, 134

Fender, R., & Belloni, T. 2004, ARA&A, 42, 317

Fermi Large Area Telescope Collaboration 2021, ATel, 15110

Finke, J. D., & Becker, P. A. 2014, ApJ, 791, 21

Fragile, P. C., Blaes, O. M., Anninos, P., et al. 2007, ApJ, 668, 417

Fragile, P. C., & Meier, D. L. 2009, ApJ, 693, 771

Fuhrmann, L., Larsson, S., Chiang, J., et al. 2014, MNRAS, 441, 1899

Ghisellini, G., Tavecchio, F., Maraschi, L., Celotti, A., & Sbarrato, T. 2014, Natur, 515, 376

Gold, R., Paschalidis, V., Etienne, Z. B., et al. 2014, PhRvD, 89, 064060

Graham, M. J., Djorgovski, S. G., Stern, D., et al. 2015a, Natur, 518, 74

Grindlay, J., Tang, S., Los, E., et al. 2012, in IAU Symp. 285, New Horizons in Time Domain Astronomy, ed. R. E. Griffin, R. J. Hanisch, & R. L. Seaman (Cambridge: Cambridge Univ. Press), 29

Hardee, P. E., & Rosen, A. 1999, ApJ, 524, 650

Holgado, A. M., Sesana, A., Sandrinelli, A., et al. 2018, MNRAS, 481, L74

Honma, F., Matsumoto, R., & Kato, S. 1992, PASJ, 44, 529

Hughes, P. C., Aller, M. F., & Aller, H. D. 2015, ApJ, 799, 207

Ingram, A., & Done, C. 2009, MNRAS, 397, L101

Itoh, R., Nalewajko, K., Fukazawa, Y., et al. 2016, ApJ, 833, 77

Ivanov, P. B., Igumenshchev, I. V., & Novikov, I. D. 1998, ApJ, 507, 131

Kalberla, P. M. W., Burton, W. B., Hartmann, D., et al. 2005, A&A, 440, 775

Karouzos, M., Britzen, S., Witzel, A., Zensus, J. A., & Eckart, A. 2012, A&A, 537, AA112

Katz, J. I. 1997, ApJ, 478, 527

Kejriwal, S., Witzany, V., Zajacek, M., et al. 2024, MNRAS, 532, 2143

Kelley, L. Z., Haiman, Z., Sesana, A., et al. 2019, MNRAS, 485, 1579

King, O. G., Hovatta, T., Max-Moerbeck, W., et al. 2013, MNRAS, 436, L114

Kirsch, M. G., Briel, U. G., Burrows, D., et al. 2005, Proc. SPIE, 5898, 22

Komossa, S. 2006, MmSAI, 77, 733

Krimm, H. A., Holland, S. T., Corbet, R. H. D., et al. 2013, ApJS, 209, 14

Larsson, S. 2012, arXiv:1207.1459

Lasky, P. D., Melatos, A., Ravi, V., & Hobbs, G. 2015, MNRAS, 449, 3293

Levine, A. M., Bradt, H., Cui, W., et al. 1996, ApJL, 469, L33

Lico, R., Liu, J., Giroletti, M., et al. 2020, A&A, 634, A87

Linial, I., & Metzger, B. D. 2023, ApJ, 957, 34

Lister, M. L., Aller, H. D., Aller, M. F., et al. 2009, AJ, 137, 3718

Lister, M. L., Aller, M. F., Aller, H. D., et al. 2013, AJ, 146, 120

Liu, T., Gezari, S., Heinis, S., et al. 2015, ApJL, 803, L16

Lomb, N. R. 1976, Ap&SS, 39, 447

Mallat, S. G., & Peyre, G. 2008, A Wavelet Tour of Signal Processing: The Sparse Way (3rd Edition) (New York: Academic)

Max-Moerbeck, W., Hovatta, T., Richards, J. L., et al. 2014, MNRAS, 445, 428

McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083

Meyer, M., Scargle, J. D., & Blandford, R. D. 2019, ApJ, 877, 39

Middei, R., Perri, M., Puccetti, S., et al. 2023, ApJL, 953, L28

Miller, M. C., & Krolik, J. H. 2013, ApJ, 774, 43

Nakamura, M., & Meier, D. L. 2004, ApJ, 617, 123

Nicalastro, F., Kaastra, J., Krongold, Y., et al. 2018, Natur, 558, 406

Nilsson, K., Lindfors, E., Takalo, L. O., et al. 2018, A&A, 620, A185

Nixon, C., King, A., & Price, D. 2013, MNRAS, 434, 1946

O'Neill, S., Kiehlmann, S., Readhead, A. C. S., et al. 2022, ApJL, 926, L35

Ostorero, L., Villata, M., & Raiteri, C. M. 2004, A&A, 419, 913

Otero-Santos, J., Acosta-Pulido, J. A., Becerra González, J., et al. 2020, MNRAS, 492, 5524

Pasham, D. R., Tombesi, F., Sukova, P., et al. 2024, SciA, 10, 13

Peñil, P., Domínguez, A., Buson, S., et al. 2020, ApJ, 896, 134

Peñil, P., Otero-Santos, J., Ajello, M., et al. 2024b, MNRAS, 529, 1365

Peñil, P., Westernacher-Schneider, J. R., Ajello, M., et al. 2024a, MNRAS, 527, 10168

Peters, P. C. 1964, PhRv, 136, 1224

Pihajoki, P., Valtonen, M., & Ciprini, S. 2013, MNRAS, 434, 3122

Piner, B. G., & Edwards, P. G. 2014, ApJ, 797, 25

Piner, B. G., Pant, N., & Edwards, P. G. 2010, ApJ, 723, 1150

Prokhorov, D. A., & Moraghan, A. 2017, MNRAS, 471, 3036

Raiteri, C. M., Stamerra, A., Villata, M., et al. 2015, MNRAS, 454, 353

Readhead, A. C. S., Lawrence, C. R., Myers, S. T., et al. 1989, ApJ, 346, 566

Reardon, D. J., Zic, A., Shannon, R. M., et al. 2023, ApJL, 951, L6

Ren, H. X., Cerruti, M., & Sahakyan, N. 2023, A&A, 672, A86

Richards, J. L., Max-Moerbeck, W., Pavlidou, V., et al. 2011, ApJS, 194, 29

Rieger, F. 2019, Galax, 7, 28

Rieger, F. M. 2004, ApJL, 615, L5

Rieger, F. M. 2007, Ap&SS, 309, 271

Roedig, C., Sesana, A., Dotti, M., et al. 2012, A&A, 545, A127

Roland, J., Teyssier, R., & Roos, N. 1994, A&A, 290, 357

Romero, G. E., Chajet, L., Abraham, Z., & Fan, J. H. 2000, A&A, 360, 57

Ryan, J. L., Siemiginowska, A., Sobolewska, M. A., et al. 2019, ApJ, 885, 12

Sandrinelli, A., Covino, S., Dotti, M., et al. 2016, AJ, 151, 54

Sandrinelli, A., Covino, S., Treves, A., et al. 2018, A&A, 615, A118

Sarkar, A., Gupta, A. C., Chitnis, V. R., et al. 2021, MNRAS, 501, 50

Scargle, J. D. 1982, ApJ, 263, 835

Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103

Schulz, M., & Mudelsee, M. 2002, CG, 28, 421

Sesana, A., & Khan, F. M. 2015, MNRAS, 454, L66

Shapiro, S. L., & Teukolsky, S. A. 1983, Black Holes, White Dwarfs, and Neutron Stars: the Physics of Compact Objects (New York: Wiley), 633

Sobacchi, E., Sormani, M. C., & Stamerra, A. 2017, MNRAS, 465, 161

Stirling, A. M., Cawthorne, T. V., Stevens, J. A., et al. 2003, MNRAS, 341, 405

Takalo, L. O., Nilsson, K., Lindfors, E., et al. 2009, in AIP Conf. Proc. 1085, HIGH ENERGY GAMMA-RAY ASTRONOMY, ed. F. A. Aharonian, W. Hofmann, & F. Rieger (Melville, NY: AIP), 705

Tang, Y., MacFadyen, A., & Haiman, Z. 2017, MNRAS, 469, 4258

Tarnopolski, M., Żywucka, N., Marchenko, V., et al. 2020, ApJS, 250, 1

Tavani, M., Cavaliere, A., Munar-Adrover, P., et al. 2018, ApJ, 854, 11

Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2011, MNRAS, 418, L79

Urry, C. M., & Padovani, P. 1995, PASP, 107, 803

Valtonen, M. J., Lehto, H. J., Nilsson, K., et al. 2008, Natur, 452, 851

VanderPlas, J. T. 2018, ApJS, 236, 16

Vaughan, S., Uttley, P., Markowitz, A. G., et al. 2016, MNRAS, 461, 3145

Villata, M., & Raiteri, C. M. 1999, A&A, 347, 30

Vlahakis, N., & Tsinganos, K. 1998, MNRAS, 298, 777

Wilkins, D. C. 1972, PhRvD, 5, 814

Woo, J.-H., Urry, C. M., van der Marel, R. P., Lira, P., & Maza, J. 2005, ApJ, 631, 762

Xu, H., Chen, S., Guo, Y., et al. 2023, RAA, 23, 075024

Yan, D., Zhou, J., Zhang, P., et al. 2018, ApJ, 867, 53

Yang, S., Yan, D., Zhang, P., et al. 2021, ApJ, 907, 105