The variability patterns of the TeV blazar PG 1553+113 from a decade of MAGIC and multi-band observations

PG 1553+113 is one of the few blazars with a convincing quasi-periodic emission in the gamma-ray band. The source is also a very high-energy (VHE;>100 GeV) gamma-ray emitter. To better understand its properties and identify the underlying physical processes driving its variability, the MAGIC Collaboration initiated a multiyear, multiwavelength monitoring campaign in 2015 involving the OVRO 40-m and Medicina radio telescopes, REM, KVA, and the MAGIC telescopes, Swift and Fermi satellites, and the WEBT network. The analysis presented in this paper uses data until 2017 and focuses on the characterization of the variability. The gamma-ray data show a (hint of a) periodic signal compatible with literature, but the X-ray and VHE gamma-ray data do not show statistical evidence for a periodic signal. In other bands, the data are compatible with the gamma-ray period, but with a relatively high p-value. The complex connection between the low and high-energy emission and the non-monochromatic modulation and changes in flux suggests that a simple one-zone model is unable to explain all the variability. Instead, a model including a periodic component along with multiple emission zones is required.


INTRODUCTION
Variability on a wide range of time scales is a common trait of active galactic nuclei (AGN, see Hovatta & Lindfors 2019, for a recent review).Among AGN, the blazar subclass is dominated by the emission from the relativistic outflow emerging from the supermassive black hole (SMBH) and structured in beamed jets, namely, forming a small angle with the observer.The relativistic beaming amplifies the observed radiation and results in a broad spectral energy distribution (SED) dominated by the non-thermal continuum, showing two main humps.The first peaks in the infrared (IR) to X-ray region, presumably originated by synchrotron radiation, and the second hump spans from the UV up to the TeV.Blazars are distinguishable by their flux variability, with a large amplitude -up to two orders of magnitude-and fast variations, down to a few minutes.Depending on the equivalent width (EW) of the absoption line in the optical band, blazars are further divided into Flat Spectrum Radio Quasars (FSRQs, EW > 5 Å) and BL Lac objects (BL Lacs, EW < 5 Å).The location of the synchrotron peak in blazars determines a further division into low-(ν S peak < 10 14 Hz), intermediate-(10 14 < ν S peak < 10 15 Hz), and high-synchrotron peak (ν S peak > 10 15 Hz) BL Lacs.The mechanisms driving the variability in blazars are still debated, and several different interpretations were suggested, related to the physical processes in the jet or in the accretion mechanism (e.g., Raiteri et al. 2017b; Marscher 2016, and references therein).
Jet emission models can be tested by looking at the variability pattern and correlations between the low-energy and high-energy humps of the SED (Falomo et al. 2014;Liodakis et al. 2019), built through contemporaneous multiwavelength (MWL) monitoring campaigns (Ulrich et al. 1997).
In this framework, the discovery of regular and periodic patterns in otherwise apparently stochastic variability can provide a deeper insight into the underlying processes.
PG 1553+113 is one of the brightest HBL emitting at gamma-ray energies.The search of intervening absorption in its far-UV spectrum suggest a redshift in the range 0.413 < z < 0.56 (Danforth et al. 2016(Danforth et al. , 2010)).Recent optical/UV obser-vations seems to converge on a value of 0.433 for the redshift of this source (Dorigo Jones et al. 2022;Johnson et al. 2019).
Its comparably large redshift entails a strong attenuation of the gamma-ray flux above energies E > 250 GeV due to pair production with photons of the extragalactic background light (EBL).Albeit this attenuation, PG 1553+113 is a well known TeV gamma-ray emitter observed by all major imaging air Cherenkov telescopes (IACTs) (Abramowski et al. 2015;Aleksić et al. 2012;Aliu et al. 2015).
The continuous gamma-ray lightcurve was measured with the Large Area Telescope (LAT) on board the Fermi Gammaray Space Telescope since 2008 and shows a clear signature of a periodic modulation of 2.18 ± 0.08 years at E > 100 MeV and E > 1 GeV, covering 3.5 cycles (Ackermann et al. 2015).The periodicity has < 1% chance of being due to random variability.This was the first time such a periodicity has been found convincingly in a gamma-ray blazar.The signature has been confirmed by several other works even in recent times (e.g.Covino et al. 2020;Peñil et al. 2020).Optical fluxes correlate with gamma-ray emission at 99% confidence and the optical lightcurve shows evidence for modulation of 2.06±0.05years over 4.5 cycles.
The mechanisms underlying the daily/weekly variability typical of blazars and the claimed periodicity are not fully understood.The long-term monitoring of the variable lightcurves of such objects is a powerful tool for discovering such processes, and observations with a complete MWL coverage are needed to understand them.
The last years have witnessed the growth of new key observations of possible periodic behaviours in AGN and blazars, evaluated with different methods Bhatta (e.g. 2017); Ait Benkhali et al. (e.g. 2020), and also Covino et al. (2019) for a cautious approach.The source PG 1553+113 has been used extensively to this end (e.g., Sobacchi et al. 2017).The periodicity should be compared to the lifetime and duty cycle of AGN activity, generally assumed as 10 7 − 10 8 years, or to episodic nuclear activity and jet formation ∼ 10 5 years (Haehnelt & Rees 1993;Konar et al. 2013).For this reason, any claim of periodicity must establish compelling statistical evidence emerging from stochastic fluctuations that can mimic a periodic pattern.
Because the emission of HBLs across the electromagnetic spectrum is dominated by the jet, the quasi-periodic modu-lation is most probably associated with the jet itself or with the processes at its base.In the latter case, disk perturbations or instabilities can induce a variation of the accretion rate advected to the jet, with quasi-periodic behaviour (e.g., Tchekhovskoy et al. 2011a).Among the possible interpretations, the modulation can be driven by the interaction of two SMBHs (see e.g., Ackermann et al. 2015, and references therein).Therefore, PG 1553+113 is a candidate for harboring a close binary SMBH system with milli-pc separation in the early inspiral gravitational-wave driven regime prior to coalescence (see e.g., Tavani et al. 2018).However, different viable solutions are possible.The observed modulation can be related to the jet itself, such as jet precession or intrinsically rotating flow or helical jet, or to the process feeding the jet.Each interpretation may lead to different expectations about the evolution of the modulation at different wavelengths.Geometrical models, e.g., due to jet precession (Romero et al. 2000), rotation (Camenzind & Krockenberger 1992;Rieger & Mannheim 2000), or helical structure (Conway & Murphy 1993), would produce a quasi-periodic variation of the beaming factor due to the change of the viewing angle (Rieger 2004).Therefore, almost achromatic variability is expected at all wavelengths, or with clear correlations between different wavelengths (Villata & Raiteri 1999).To test this hypothesis, PG 1553+113 was monitored with the very long baseline array (VLBA) at 15, 24, and 43 GHz in the period 2015-2017, a full cycle of gamma-ray activity (Lico et al. 2020).VLBA data provided evidence of jet angle variations, indicating that geometric effects could play a role in the observed emission variability through Doppler boosting modulation.However, no clear connection was found between the observed variations and the quasi-periodic variability patterns reported in optical and gamma rays.Therefore, additional mechanisms are necessary to explain the variable broadband emission.
Pulsating instabilities occurring in the disk can explain the quasi-periodicity, also in the case of a single SMBH (Honma et al. 1992), as in magnetically dominated and magnetically arrested accretion flows in the inner disc part (Tchekhovskoy et al. 2011b).
In case the observed periodicity is interpreted as a periodic perturbation due to the interaction with a secondary black hole, a double or multiple peak substructure is expected at different bands, resulting from the interaction of the secondary black hole with the accretion disk near the periastron (Lehto & Valtonen 1996).Such a structure is apparent in the optical lightcurve, and the detection of similar double peaks in the X-ray and gamma-ray lightcurves would confirm this interpretation.
PG 1553+113 is also known to be variable on a weekly time scale (Aleksić et al. 2015;Abramowski et al. 2015;Aliu et al. 2015).The study of flaring episodes provided an important input for the modelling of intrinsic emission from the source and was used as a probe of fundamental physics (Guo et al. 2020) and of the EBL (Acciari et al. 2019;Korochkin et al. 2020).Finally, the reported periodicity could be accidental due to stochastic background fluctuations, typically found in lightcurves of AGN and blazars.For this reason, proper statistical approaches such as those discussed in Vaughan et al. (2003) are needed when studying long-term light curves (e.g., Covino et al. 2019).Multiwavelength observations are needed to support the physical interpretation of the emission from the innermost regions of the blazar and its jet.
With the purpose of characterizing PG 1553+113 broadband variability and testing physical scenarios, the MAGIC collaboration initiated a multi-year, MWL monitoring campaign in 2015.In the campaign, various instruments observing in the radio, infrared, optical (both photometry and polarimetry), UV, soft X-ray, and gamma-ray bands were involved.In this article we report the results of this campaign, including data from previous observations.The paper is structured as follows: in §2, we present the details of the MAGIC and MWL data analyses.§ 3 is dedicated to the characterization of the source variability, while § 4 is focused on periodicity study of the data presented in the paper.In § 5 the conclusions of this work are outlined.

MULTIWAVELENTH DATA AND ANALYSIS
The MWL data used in this study span several orders of magnitude, from radio to VHE gamma rays.In this section, we give a brief description of the collected data and their analyses (from higher to lower frequencies), while in the following sections we report on the scientific interpretation of the data.

VHE gamma-ray data
MAGIC is a system of two IACTs located in La Palma, Canary Islands, at 2200 m asl.It observes VHE gamma rays from 50 GeV up to tens of TeV.The angular resolution at around 200 GeV energies at low zenith angles (0 • -30 • ), is < 0.07 • , while the energy resolution is 16% (Aleksić et al. 2016b).From 2004 to 2009, MAGIC was composed of a single IACT, MAGIC-I (Baixeras et al. 2004).
The MAGIC data are a collection of images registered by the camera of each telescope and are processed with a standard analysis chain (Zanin et al. 2013).
Since its detection in 2004, MAGIC observed PG 1553+113 every year.In the first two years of observation (2005 and 2006) the data were affected by large systematic and statistical errors, therefore, in this work we consider only data taken from 2007 on, where an upgrade of the telescope readout sensibly increased its performances (Aleksić et al. 2016a).For the study, we use the dark night data from 2007 to 2017, including already published data from 2007, 2008, 2009, and 2012campaigns Aleksić et al. (2012, 2015)).
Table 1 summarises the results of the MAGIC data analysis of the samples collected from 2007 to 2017.Data from 2010 on were analysed with monthly binning, while for previous data, the overall yearly sample is considered as reported in Aleksić et al. (2012).The total observation time is ∼100 hours, non uniformly distributed across the years (third column).Average fluxes above 150 GeV in physical units and in Crab Nebula units 1 are reported in columns four and five.The flux varies from 4 to 18% C.U. across time.Finally, the last three columns list the results of the differential flux analysis, reporting the results of a fit with a simple power law function in the form where f0 is the flux at 200 GeV (column six) and Γ is the power law index (column seven in Table 1).
The overall emission above 150 GeV of PG 1553+113 observed with MAGIC from 2007 to 2017 is reported in Figure 1.Data from 2007Data from , 2008Data from , and 2009 were collected with a single telescope and are from Aleksić et al. (2012).For the more recent data, a daily binning was adopted.The 2010-2017 monthly-averaged values are listed in Table 1, along with 2007-2009 yearly values from Aleksić et al. (2012).MAGIC started a regular monitoring of the source for seven months per year in 2015 (MJD ∼ 57000).This explains the irregular and scarce sampling of the curve before 2015.
The daily flux above 150 GeV shows variations within a factor of ∼ 10, and reached its maximum in 2012 during a historical flare reported in Aleksić et al. (2015).The average flux is (2.74±0.04)• 10 −11 cm −2 s −1 .The hypothesis of constant flux can be discarded, based on the χ 2 test for goodness of fit (χ 2 /(degrees of freedom)=1339/157).

High energy gamma-ray data
Further gamma-ray data considered in the study are those collected with the Fermi-LAT and analysed above 100 MeV in Tavani et al. (2018).In that work the authors added 2016 and 2017 data to the sample used in the original paper claiming the quasi-periodic behaviour (Ackermann et al. 2015).The Fermi-LAT data therefore cover almost continuously the 2008-2017 period and represent the only continuous monitoring considered.These data are analysed with 20 days binning.

X-ray data
The Neil Gehrels Swift observatory (Swift) (Gehrels et al. 2004) observed PG 1553+113 since 2005 during outbursts and almost regularly since 2013.We have collected all snapshots in the period from 2005 up to the end of 2017.PG 1553+113 was observed simultaneously with the the Xray Telescope (XRT; Burrows et al. 2005, 0.2-10.0keV), and with all six filters of the Ultraviolet/Optical Telescope (UVOT; Roming et al. 2005, 170-600 nm).The Swift-XRT data, reported in Fig. 2, were collected in photon counting mode (PC) and windowed timing mode (WT).In both cases, the data were processed using the FTOOLS task xrtpipeline (version 0.13.5), which is distributed by HEASARC within the HEASoft package (v6.28).Events with grades 0−12 were selected for the data in PC mode, and with grades 0 − 2 for the data in WT mode.The corresponding response matrices available in the Swift CALDB version were used.When the source count rate in PC mode was higher than 0.6 counts s −1 the pile-up was evaluated following the standard procedure 2 .Observations affected by pile-up were corrected by masking the central 7 arcsec region.For each observation the extraction region was checked visually on the image and slightly centred to the peak of the signal, when needed.The signal was extracted within an annulus with an inner radius of 3 pixels (7 arcsec) and an outer radius of 30 pixels (70 arcsec).Events in different channels were grouped with the corresponding redistribution matrix (rmf), and ancillary (arf) files with the task grppha, setting a binning of at least 25 counts for each spectral channel in order to use the chi-squared statistics.The resulting spectra were analyzed with Xspec version 12.11.1.We fitted the spectrum with an absorbed power-law using the photoelectric absorption model tbabs (Wilms et al. 2000), with a neutral hydrogen column density fixed to its Galactic value (NH = 3.67 × 10 20 cm −2 ; Kalberla et al. 2005).The results are shown in Fig 2.

UV data
Swift/UVOT data in the v, b, u, w1, m2, and w2 filters are reported in Fig. 3 and were reduced with the HEAsoft package v6.28 using the uvotsource task.We extracted the source counts from a circle with 5 arcsec radius centred on the source nominal position, corresponding to the optimal aperture on the source count rate (Poole et al. 2008).The background counts were extracted from a circle with 60 arcsec radius in a near, source-free region.Conversion of magnitudes into dereddened flux densities was obtained by adopting the extinction value E(B-V) = 0.054 as in Raiteri et al. (2015a), the mean galactic extinction curve in Fitzpatrick (1999) and the magnitude-flux calibrations by Poole & Breeveld (2005).Statistical uncertainty on magnitudes of the order of 0.03 mag, on the zero-point UVOT calibration 0.02-0.06mag (Poole et al. 2005) and the count ratio to flux correction (Poole & Breeveld 2005) have been propagated to estimate the error on the flux, resulting in a 4% to 6% uncertainty.

Optical data
The optical R-band data were obtained as part of Tuorla blazar monitoring program3 .The observations are described in Nilsson et al. (2018).The data have been analysed using the semiautomatic pipeline for differential photometry developed at the Tuorla Observatory (Nilsson et al. 2018) using the comparison and control star magnitudes from Raiteri et al. (2015b).The observed fluxes have been corrected for galactic extinction using a value of 0.113 from Schlafly & Finkbeiner (2011).
Other optical data were provided by the Whole Earth Blazar Telescope (WEBT) Collaboration4 .WEBT observations up to 2015 October were analysed in Raiteri et al. (2015bRaiteri et al. ( , 2017a)).New data in the 2016 and 2017 optical observing seasons were acquired at the following observatories: Abastumani (Georgia), Aoyama Gakuin (Japan), Crimean (Crimea 5 ), Hans Haffner (Germany), Mt.Maidanak (Uzbekistan), Perkins (US), Rozhen (Bulgaria), Siding Spring (Australia), Siena (Italy), Sirio (Italy), St. Petersburg (Russia), 54000 54500 55000 55500 56000 56500 57000 57500 58000  Teide (Spain), Tijarafe (Spain), and at the Astronomical Station Vidojevica (Serbia).Additional WEBT observations were carried out with telescopes belonging to the Las Cumbres Observatory global telescope network at the Haleakala, Siding Spring, and Teide observing sites.Calibration was performed using the same photometric sequence as in the case of KVA data.The R-band lightcurve obtained by assembling all datasets is shown in Fig. 3 and was carefully inspected and, when necessary, processed to obtain a homogeneous and reliable time series.Indeed, even if WEBT observers use the same photometric sequence, differences in equipment may lead to some offset between various datasets.These offsets clearly appear when datasets overlap in time and can consequently be corrected for.Moreover, we removed a few data points strongly deviating from the main trend traced by the bulk of the datasets, and mostly affected by large uncertainties.Finally, noisy intra-night sequences from the same telescope were binned.The above processing is a necessary step to undertake if one wants to deal with light curves that can be used for robust analysis and modelling.The details of the data reduction of the AZT-8+ST7 data are described in Larionov et al. (2008), and the Perkins telescope observations in Jorstad et al. (2010).The Steward Observatory data are publicly available and the polarimetric data are described in detail in Smith et al. (2009).
The data were obtained using the R-band filter except for the Steward Observatory where data were obtained using a filter between 5000 − 7000 Å.All data were checked for consistency and the polarisation degree was corrected for positive bias using the formula in Wardle & Kronberg (1974).We removed six data points, which had a signal-to-noise ratio less than two in fractional polarisation.In 2016, RoboPol observed the source with a faster cadence of multiple observa-tions per night, and we averaged these to a single observation per night to avoid biasing our analysis with more densely sampled curves during that time.
In Fig. 3, we show the EVPAs starting from a range between 0 • and 180 • .The difference between the EVPAs of consecutive points is minimised to be less than 90 • by adding or subtracting 180 • from the following points.If the time gap between the points is longer than 50 days, we set the EVPA to the original range of 0 • − 180 • as we do not know the evolution of the EVPA over such long gaps.
In all other cases, except for the Steward Observatory observations, photometry is also determined from the observations.The R-band magnitudes from these observations are also included in Fig. 3 along with the KVA and WEBT data.The polarisation data are used in this paper for charaterising the general variability patterns, while more detailed physical modelling of the polarisation is the subject of a separate paper (Nilsson et al. in prep.).

IR data
We observed PG 1553+113 with the Rapid Eye Mounting Telescope (REM, Zerbi et al. 2004), a robotic telescope located at La Silla Observatory (Chile).It performed photometric observations using NIR filters in the period from April 08, 2005 (MJD 53468) to June 20, 2017 (MJD 57802).REM data are shown in Fig. 3.The telescope is able to operate in a fully autonomous way (Covino et al. 2004), and data are reduced and analysed following standard procedures.Aperture photometry was derived using custom tools, and the calibration was based on the scheme described by Sandrinelli et al. (2014).
We used reference stars from the Two Micron All Sky Survey (2MASS) Catalog 6 (Skrutskie et al. 2006).All images have been visually checked, eliminating those where the targets or the reference stars are close to the borders of the frame, and where obvious biases were present.

Radio data
Regular 15 GHz observations of PG 1553+113 were carried out as part of a high-cadence gamma-ray blazar monitoring programme using the Owens Valley Radio Observatory (OVRO) 40 m telescope (Richards et al. 2011).PG 1553+113 was observed with a nominal twice-per-week cadence.
The OVRO 40 m uses off-axis dual-beam optics and a cryogenic receiver with a 3 GHz bandwidth centered at 15 GHz.The two sky beams are Dicke switched using the off-source beam as a reference, and the source is alternated between the two beams in an ON-ON fashion to remove atmospheric and ground contamination.In May 2014, a new dual-beam correlation receiver was installed on the 40 m telescope and the fast gain variations are corrected using a 180 degree phase switch instead of a Dicke switch.The performance of the new receiver is very similar to the old one and no discontinuity is seen in the lightcurves (see Fig. 3).Flux density calibration is achieved using a temperature-stable diode noise source to remove receiver gain drifts and the flux density scale is derived 6 http://www.ipac.caltech.edu/2mass/from observations of 3C 286 assuming the Baars et al. (1977) value of 3.44 Jy at 15.0 GHz.The systematic uncertainty of about 5% in the flux density scale is not included in the error bars.Complete details of the reduction and calibration procedure are found in Richards et al. (2011).
Radio observations were also carried out with the 32-m dishes located in Medicina (at 8 and 24 GHz) and Noto (at 5 GHz).Continuum acquisitions were performed exploiting On-The-Fly cross-scans in Equatorial coordinates.Flux density calibration was carried out observing 3C 286, 3C 123, NGC 7027 and, prior to January 2018, also 3C 48.Reference flux densities for the calibrator sources were computed for the observed band central frequency, according to Perley & Butler (2013).For 24-GHz observations, the atmospheric contribution was also taken into account in the calibration procedure; the zenith opacity was estimated by means of skydip acquisitions.The data reduction was performed using the Cross-Scan Analysis Pipeline described in Giroletti & Righini (2020).A large part of the HE gamma-ray data dataset, as well as radio and optical datasets, are published in Abramowski et al. (2015) and were used for the periodicity analysis.

CHARACTERIZATION OF THE VARIABILITY
The only instrument considered in this work that performed a continuous monitoring is Fermi-LAT.Also radio data have very good coverage, followed by optical data that suffered only from a few months break per year related to the visibility of the source.The coverage is more scattered for IR, UV, X-rays, and VHE gamma-ray data that are strongly affected by sparse sampling and often the observations are driven by a high state alert trigger, and therefore may be biased towards high states.In these bands, the coverage had a clear improvement starting from late 2014 (MJD ∼57000).This is the result of an intense MWL and multi-year campaign aimed at a precise monitoring of the source state for the detailed modelling of the source emission and periodicity.
From a visual inspection of Figure 3 we can conclude that the source shows high variability over the years in all bands, with moderate variations in radio and more pronounced variations in the other bands.This behaviour is quite common in HBLs (Acciari et al. 2021).
A detailed characterisation of the variability is the key to investigate the physical phenomenology responsible for the broad-band emission as detailed in Rieger (2019) and references therein.In the following subsections, several variability studies are presented.The aim is twofold: first, the characterisation of the variability (and periodicity) at different bands.Secondly, the identification of inter-band connections.These connections are a powerful tool to unveil single/multiple regions responsible for the observed emission.

Flux-spectral index correlation
A clear correlation between the integral flux and the slope of the power law approximating the differential energy flux in X-rays and gamma rays characterises many flaring events of BL Lacs.In the case of negative correlation, the effect is often referred as harder-when-brighter behaviour (see, e.g., Albert et al. 2007).This behaviour is associated to a shift of the synchrotron peak during flares, meaning that more energetic electrons are responsible for the bulk of the emission (see e.g.Acciari et al. 2021).
A study of the flux-slope correlation was performed for the VHE gamma-ray and X-ray data.The results are shown in Fig. 4. The monthly averaged MAGIC data from 2010 to 2017 were considered for the study.No correlation appears in the MAGIC data, but the large time interval considered and the relatively low statistics involved in the study may have diluted this correlation.X-ray data show instead a hint of anticorrelation between the spectral index and the flux state, indicating a harder-when-brighter trend.To evaluate the level of correlation we adopted the Spearman correlation coefficient, a value close to (-)1 pointing to a strong (anti-)correlation.The Spearman correlation coefficient in X-ray data is -0.39, and the p-value of the null-hypothesis (i.e., nocorrelation) is ∼ 10 −10 .
It is important to emphasise that, contrary to the vast majority of studies available in literature claiming a harderwhen-brighter trend in the considered bands (e.g., Pian et al. 1998;Albert et al. 2007;Acciari et al. 2020), the results presented here are not from a single campaign/flare, but they come from a large time interval spanning more than 10 years.

Variability timescale
Apart from few exceptions (e.g., Acciari et al. 2020), blazars are highly variable objects in almost all bands.In HBLs, the variability timescale may range from months down to a minute timescale.According to special relativity, the comoving size of the emitting region, ∆r ′ , can be constrained by the variability timescale, (as detailed in Rieger 2019).The variability timescale pinpoints the properties of the region responsible for the observed radiation.Interestingly, subday variations represent a challenge for the simplest emission models in blazars (e.g.Tavecchio et al. 1998), and for the subclass of Flat Spectrum Radio Quasars (e.g., fast variability-10 minutes doubling time-observed in VHE gamma rays in PKS 1222+21, Aleksić et al. 2011).
To constrain the variability timescale of PG 1553+113, a search for intraday variations was performed on MAGIC and Swift-XRT data.For statistical reasons, the study was limited to 10 snapshots with the highest flux recorded in both bands.The average duration of MAGIC observations was 1 hour, while the average duration of Swift-XRT observations was 1.2 ks.The analyses revealed no hint of intraday variability in Swift-XRT and MAGIC data.

Interband correlation
Short or long-term correlation between the fluxes emitted at different bands allows us to track down the connection between photons emitted at different possible regions in the jet or with different, but correlated, mechanisms.This is the case of synchrotron self-Compton emission (SSC), where lowenergy synchrotron photons are emitted together with inverse Compton, high-energy photons produced by the same electrons upscattering the synchrotron photons.
In our study, we focus on the interband correlation search on the IR, optical, UV, X-ray, HE and VHE gamma-ray bands.Radio data are excluded from this study in consideration of the well known lag due to the different location of the emission zone, that will be further discussed in the text.
The results of the correlation analysis performed with the Spearman test as implemented in the SciPy python package are listed in Tab. 2, where the Spearman coefficient and the p-value appear in the third and fourth columns, respectively.
For this study, only data within a 1.5-day window have been considered simultaneous apart for the correlation studies involving Fermi-LAT data where the simultaneity window has been extended to ± 10 days, to have sufficient statistics and in agreement with the Fermi-LAT data binning of 20 days.The Table is ordered by decreasing the correlation levels according to this indicator.Similar results are obtained with the weighted Pearsons coefficient (that has the advantage of taking into account the errors of the flux, but the disadvantage of assuming a Gaussian distribution of values).Fig. 5 and 6 show some selected scatter plots, also indicated in the last column of Table 2.
The main results of the correlation analysis are as follows: -A strong correlation (Spearman coefficient ⩾ 0.9) is observed both between optical and UV (UW2 band) and optical and IR (H band) data; -Optical and UV data show a net correlation (Spearman coefficient ⩾ 0.6) also with the HE gamma ray data.A similar relation is also observed between X-rays and VHE gammarays; -X-ray and UV data, with strictly simultaneous sampling, show a milder correlation (Spearman coefficient = 0.55); -Only a hint of correlation (Spearman coefficient ⩽ 0.4) with close to zero time lag is observed between the other bands (HE and VHE gamma rays, optical/X-ray, and optical VHE gamma rays).
-In the optical/X-ray case, a careful analysis of the scatter plot allows us to identify episodes with different correlation behaviours: from a clear correlation, corresponding to the Xray and optical flare at MJD ∼ 57000 to anti-correlated events (X-ray enhanced state at MJD ∼ 57500), as highlighted in Fig. 7.These results suggest a common origin for the spectral features observed in IR, optical, UV, and Fermi-LAT bands.In particular, IR, optical, and UV photons are likely synchrotron photons from the same emitting region.Single-zone SSC process is the cogent mechanism connecting optical and HE gamma-ray photons.The same process may be responsible for the X-ray to VHE gamma-ray connection, even if the radiation should come at least in part from a different (or additional) region with respect to the low-energy counterpart, to explain the weaker correlation with the other bands.
The possibility of a delay in the PG 1553+113 correlation between bands has recently been investigated in Liodakis et al. (2018) for the radio, optical, and gamma-ray bands.Although gamma-ray and optical data are consistent with no time-lag correlation, a delay of ∼3-4 months appears between radio and both optical and Fermi-LAT data.We have inves- Table 2. Results of the correlation study between integral flux in different bands ordered by decreasing Spearman coefficient (third column).The simultaneity window assumed is ±1.5-day apart for the correlation studies involving Fermi-LAT data, where it has been extended to ± 10 days.The last two columns report the p-value (null hypothesis: no correlation with close to zero time lag) and, if available, the panel with the scatter plot in Fig. 5, respectively.

Band-1
Band- tigated the possibility of a delay between radio and X-ray and VHE gamma-ray data, between Fermi-LAT and X-ray and VHE gamma-ray data, and between X-rays and VHE gamma rays with the same method presented in Liodakis et al. (2018).In all cases, from the discrete correlation function study no significant time lag emerged.The analysis of radio and optical/gamma-ray data, instead, are fully consistent with those reported in Liodakis et al. (2018).A delayed correlation in the radio band is well known effect in blazars and is due to the self-absorption of radio emission, which becomes visible when the density of the region drops, inducing a delay with respect to high-energy emission.

Bayesian block analysis
In order to determine variations and flares in the MWL light curves, we model flux variations in a model-independent manner using Bayesian blocks (Scargle et al. 2013).Figure 8 shows the Bayesian block representation of the MWL lightcurves, including the optical polarisation degree.We have optimised the prior of the slope on the number of bins for individual lightcurves to better match their sampling and variability.We find interesting flaring behaviour across all wavelengths.For example, we find contemporaneous flares in Xrays, UV, optical, and elevated activity in VHE gamma rays and polarisation at ∼56060 MJD that does not seem to have a counterpart in the Fermi-LAT band.On the other hand, we find contemporaneous flares in optical, UV, X-rays and gamma-rays at ∼57020 MJD.The 15 GHz radio also seems to be in an elevated state during that time, although the radio lightcurve does not follow the variability patterns in other wavelengths.Unfortunately, apart from ∼57020 MJD, the prominent Fermi flares fall into optical and/or VHE gammaray gaps, although in all of them the optical seems to be in an elevated state.

PERIODIC MODULATION IN THE MWL LIGHTCURVE
As first proposed in Ackermann et al. (2015) and then confirmed in several other studies (Covino et al. 2020;Sandrinelli et al. 2018;Prokhorov & Moraghan 2017;Gupta et al. 2016), a periodic modulation of the HE gamma-ray emission is firmly established in Fermi-LAT data.Although the optical curve appears to be much more complex than Fermi-LAT curve, the modulation with a period similar to that observed at higher frequencies is found, with a smaller significance (Tavani et al. 2018).
As a first step in the periodicity study, we visually inspected if the lightcurves at different bands are in agreement with the hypothesis of a periodic modulation of period P f ermi = 798 days.To this purpose, we built a normalised MWL folded lightcurve assuming a period of 798 days (Ackermann et al. 2015).The final result is shown in Fig. 9  cal, UV, X-ray, HE gamma-ray, and VHE gamma-ray bands.In each bin, the average value and its error are reported with a continuous line.The bars instead represent the square root of the variance of the folded data of each phase bin.We underline that, while the error on the average is strongly dependent on the number of points in that bin, the variance is not.Therefore, the latter is an indication of the dispersion of the sample, except for a few cases with bins with a single point (hence with a low uncertainty that is by no means representative of the real dispersion).
As expected, the HE gamma-ray folded curve displays a clear peak at phase ∼ 0. The minimum emission is instead located at phase 0.3-0.5.The same trend is observed in the optical and UV curves, while the radio curve behaviour seems shifted, in agreement with the delayed emission resulting from the DCF analysis in Liodakis et al. (2018), also confirmed in our study.
It is interesting to note that, despite the poor sampling of the lightcurves, X-ray, and VHE gamma ray-folded curves are almost fully characterised, meaning that the observations carried out in the 11 years considered were diluted in different phases of the assumed periodic modulation.This allows us a first qualitative study of the trend of these curves.The main finding of this study is that while the folded curves, including the optical polarisation, show a minimum in the phase interval 0.3-0.5, in agreement with the observations in the HE gamma-ray band, no clear maximum is observed in the optical polarisation, X-ray, and VHE gamma-ray bands.
In parallel with this visual inspection, we have also car- ried out a study of periodicity, which is reported in the next section.

Systematic search for periodicity
In all bands we search for a sinusoidal periodic signal using the Generalized Lomb Scargle (GLS) periodogram (Zechmeister & Kürster 2009) as implemented in the PyAstronomy python package 7 (Czesla et al. 2019) with frequencies ranging from 1/T to N/(2T ) sampled in steps of 1/(10T ), where T is the total time of the light curve and N is the number of data points.
We identify the period and power of the strongest peak.In the periodogram of the radio data the GLS power at the lowest frequency is higher than the next strongest peak; however, we do not take the low-frequency peak into consideration because the peak period lies outside of the covered frequency range and less than one full cycle would be covered by the full data.
We need to assess whether a detected peak provides significant evidence for an intrinsic periodicity or whether it is a sporadic result of the typical flaring behaviour.We follow the procedure described in Appendix A of O'Neill et al. (2022).Our Null hypothesis is that the light curves follow a stochastic red-noise process, with the same statistical properties of the original data -namely the Power Spectral Density (PSD) and Probability Density Function (PDF).We assume a power-law PSD ∼ ν −β , where ν is the frequency, and estimate the index, β, with an implementation of the method introduced by Uttley et al. (2002).In the radio band the index was found to be 2.0 ± 0.51 and for the MAGIC data 1.0 ± 0.42, for the remaining bands ranging from 1.4-1.5 with uncertainties ranging from 0.1-0.26.
We estimated the PDFs through the Empirical Cumulative Distribution Functions (ECDFs) of the light curves.For the Fermi-LAT analysis we created 50 000 simulations, for all other bands 10 000 simulations each, that implement the Null hypothesis.We calculate the GLS for all simulations and count the simulations that have a power equal to or higher than the peak power at the peak frequency of the data GLS; this gives us the local p-value, i.e. the probability that a rednoise process results in an apparent periodicity as strong or stronger than observed at specifically the frequency where it was detected.The local p-value does not take into account that we test many frequencies.To take the look-elsewhereeffect into account we identify the strongest peak (discarding peaks at the edge of the frequency range) in the GLS of each simulation and calculate its local p-value, then we count all simulations with peaks that have a local p-value equal to or lower than the local p-value of the data GLS peak to estimate the global p-value.
Our results are reported in Table 3, where the PSD index, the peak period in days, the peak power and the local and global p-values are listed as a function of the considered band (first column).The GLS of the Fermi-LAT light curve shows a prominent peak at a period of 786 days, with a global pvalue of 1.0 × 10 −3 , which corresponds to 3.1σ in a one-sided test if the statistic were normal distributed.At an a priori chosen significance level of 3σ we reject the Null hypothesis that the detected peak in the periodogram is a likely result of a red-noise process and conclude that light curve very likely contains a truly periodic signal.The Swift-XRT and MAGIC data do not show evidence for a significant periodicity, as the p-values suggest that the data behaviour is fully consistent with the null hypothesis.The Swift-UVOT and radio light curves show the strongest GLS peaks at periods at 806 and 865 days, but cannot be considered significant with global p-values as high as 10% and 3%.
Before the analysis of the optical light curve, including all the optical data collected in the campaign, we averaged data points within time ranges of one day weighted by the corresponding uncertainties.In the GLS of the optical light curve we identify the strongest peak at a period of 957 days with a local p-value of ∼10% and a global p-value of 47%.This result appears to be in conflict with those of Sandrinelli et al. (2018), who claim a period of 810 ± 52 days with p-values of 1% or 5% depending on the method used, and the results of Agarwal et al. (2021), who found periods in the range of 801-812 days with uncertainties ranging from 20-70 days, depending on the analysis method and specific band (V or R) with a p-value < 1%.However, we find that the strongest GLS peak, shown in Fig. 10 (upper panel) shows a broad, flat plateau covering the period range from ∼800-1000 days.Therefore, the period is poorly constrained.Furthermore, the local p-value strongly depends on the period in that period range.Towards shorter periods the local p-value decreases and is ∼ 5% at a period of ∼ 800 days, comparable to one of the results of Sandrinelli et al. (2018).Regarding the optical data we currently do not find convincing evidence for  a true periodicity and we conclude that if there truly is a periodicity in the light curve more cycles need to be covered by observations to get a better constraint on the period and the significance.
To complete our analysis, we have estimated the local pvalue in the above mentioned GLS at the literature period of 798 days (Ackermann et al. 2015).This approach, also in line with the folded lightcurve presented in the previous Section, has the advantage of minimizing the trial factors, since only a single period is tested.The resulting p-values are listed in the last column of Table 3.The Fermi-LAT p-value reflects the local p-value, which is expected as it is the period determined by a subsample of the same data.The radio, optical, Swift-UVOT p-values decrease to 2%, 7.8%, and 0.56%, respectively, while the p-values in case of MAGIC and Swift-XRT data are well compatible with the null hypothesis with values 70% and 86%, respectively.The optical and MAGIC GLS are diplayed in Figure 10.

Further search for periodicity
A significant indication of periodicity in the gamma-ray lightcurve measured by Fermi-LAT was found in several other studies.Among those, the work by Peñil et al. (2020) implemented different methods to detect periodicity in Fermi-LAT blazars.We used the same pipeline described in Peñil et al. (2020), focussing on the Lomb-Scargle and on the wavelet transform approaches, and extending it to the search for periodicity to the VHE, X-ray, UV and optical bands.The details on the two methods can be found in the corresponding paper.The results on the Fermi-LAT lightcurve is reproduced providing a peak at the period of ∼ 800 days with a p-value < 1 × 10 −4 (pre-trial).The significance of the periodicity was evaluated with respect to the null hypothesis of a stochastic red noise with the same statistical properties of the lightcurve of the original data.Similar period and significance are obtained on the optical lightcurve, in contrast to the analysis described in the previous section, but slightly compatible For these purposes, we have first characterised the variability in the VHE gamma-ray and X-ray bands.In both cases, we have not found evidence of intraday variability.Interestingly, intraday variability in the X-ray band was detected by Raiteri et al. (2017a)  X-ray data in our sample show a hint of harder-whenbrighter behaviour often detected in blazars (even if usually it is observed over a shorter timescale).This indicates that possibly freshly injected, high-energy electrons are responsible for the X-ray variability.Furthermore, the long time span considered ensures that the mechanisms driving the spectral variability did not change substantially over time.The same study applied to the MAGIC data gave inconclusive results, probably due to the 1-month averaging applied to the data.A detailed spectral study is planned in a future publication and is beyond the scope of this work.
Interband correlation studies performed on IR, optical, UV, X-ray, HE and VHE gamma-ray data confirmed the strong IR/optical/UV connection, related to the common synchrotron origin of the emission.An evident correlation between X-rays and VHE gamma rays, and between optical/UV/IR and HE gamma-ray, also emerged from this study, suggesting intertwined emission processes such as that foreseen in the multi-zone, SSC emission scenario.
Another piece of the puzzle is represented by the detected delayed correlation (of about 3-4 months) between the radio and both optical and Fermi-LAT emission (Liodakis et al. 2018).This result is in line with the average behaviour found in gamma-ray detected blazars reported in the same paper.Interestingly, this time-delayed correlation is not present in our radio and X-/VHE gamma-ray data.
After this detailed characterisation of the variability in general, we have focused on the study of the periodicity.
A search for evidence of periodicity in the X-ray and TeV bands, as well as in the optical, UV, and GeV bands was performed with a solid statistical approach.
Our main finding of the periodicity analysis with the GLS method was that the X-ray and VHE gamma-ray data do not show statistical evidence for a periodic signal.Remarkably, a (hint of a) periodic signal compatible with the one published in literature was found only in the gamma-ray data, which is also the only band with a continous coverage.A solid statistical analysis was applied to the data in the other bands.Radio, UV, and optical data show a periodogram with a peak compatible with the one firmly established in Fermi-LAT observations, but with a relatively high p-value (ranging from 4 to ∼50%).This is in agreement with the folded MWL lightcurve built assuming the literature period.The visual inspection of the MWL folded lightcurve suggests a hint of periodic behaviour in the radio, optical, UV, and HE gamma-ray bands.The peak is more pronounced in gamma rays and radio, while it appears broader, and with a similar pattern, in the R and UV2 bands.Polarimetric data, as well as X-ray and VHE gamma-ray data do not show any evident peak in the periodgram.X-ray and VHE gamma-ray folded curves exhibit a similar pattern.Interestingly, a low activity was recorded in all bands at approximately the same phase.
The observed periodicity may be interpreted as a periodic perturbation of the accretion rate on the SMBH and consequently of the fuelling at the base of the jet.The presence of a secondary black hole in a sub-parsec orbit with respect to the primary SMBH originating the jet represents a natural explanation, as previously invoked for OJ 287, despite not unique (e.g., Dey et al. 2018).Different mechanisms such as jet precession, internal jet rotation, or helical jet motion may also be invoked to explain the periodicity.
The most direct way to constrain a simple precessing jet model would be to observe motion of the jet on the sky.This has been studied in the radio band by Caproni et al. (2017) who modelled 15 GHz VLBA data of PG 1553+113 taken be-tween 2009 and 2016 using a precessing jet model.They modelled the jet using individual Gaussian components, which they then connected to episodes of gamma-ray flares.More recently, the radio jet properties of PG 1553+113 were studied by Lico et al. (2020) using VLBA observations taken between 2015-2017and lightcurves from OVRO between 2008-2018.While they found clear enhanced activity periods in the radio data, they were not found to be correlated with the 2.2-year gamma-ray periodicity.Moreover, they concluded that the position angle variations of the jet of PG 1553+113 were not correlated with the gamma-ray periodicity, and a simple geometric model where the variability is caused by changes in the Doppler boosting cannot explain the periodicity, if the gamma-ray and radio variations originate in the same region of the jet.Lico et al. (2020) also studied the radio polarisation of PG 1553+113 using their VLBA data.They found that periods of enhanced polarisation were connected with total intensity flares, indicating that the mechanism producing them is connected.Additionally, they saw a flattening of the radio spectral index at the times of total intensity activity.Such a behaviour could be expected, for example, when shocks compress magnetic field lines, which both increases the fractional polarisation and induces particle acceleration, which flattens the spectral index and increases the total intensity.They also suggested that the low polarisation observed in the core region of PG 1553+113 is due to multiple polarised components blended within the beam.
Multiple emission components are also supported by the lack of clear periodic modulation in the X-ray and VHE data, which is seen in most of the other wavelengths.On the other hand, the short-term variability in all bands is clearly correlated on some occasions (for example, data around MJD 57000 in Fig. 7 connecting optical and X-ray emission), while at other times there can even be an anti-correlation (data around MJD 57500 in the same Figure ).This shows that the situation is very complex.The difficulty to connect the lowenergy part of the synchrotron bump (IR/optical/UV) to the high-energy synchrotron part (X-rays) was extensively studied in Raiteri et al. (2015b).They studied the synchrotron spectrum of the source in multiple activity states and found that the changes in the spectrum can be explained with an inhomogeneous helical jet model, where the high-energy emission originates closer to the black hole than the lowenergy emission.Alternatively, there could be multiple (disconnected) emission components or a more complex electron distribution than typically assumed.
Our analyses on the periodicity show that there clearly must be multiple components contributing to the emission, but that they also cannot be fully disconnected because we (at least sometimes) see simultaneous flaring in all bands.Moreover, the minima in the folded lightcurves seem to be in phase in all bands.Some of the bands (X-ray and VHE) may be more sensitive to the stochastic variations only, while in the other wavelengths, connected with the low-energy part of the SED peaks, we can also see the periodic modulation.This means that any model explaining the periodicity should also be able to explain why it is more prominent in the low-energy part.Confirming this discrepancy would also require longer, densely sampled, lightcurves in X-ray and VHE energies.

Figure 2 .
Figure 2. Top panel: PG 1553+113 lightcurve in the 0.5-10 keV band as detected with the Swift-XRT satellite.Bottom panel: spectral slope lightcurve resulting from an absorbed power-law fit.

Figure 3
Figure3displays the lightcurves collected from PG 1553+113 at several wavelengths from radio (bottom panel) to VHE gamma-rays (top panel), for 12 years, from 2005 to 2017 as described in the previous section.A large part of the HE gamma-ray data dataset, as well as radio and optical datasets, are published inAbramowski et al. (2015) and were used for the periodicity analysis.The only instrument considered in this work that performed a continuous monitoring is Fermi-LAT.Also radio data have very good coverage, followed by optical data that suffered only from a few months break per year related to the visibility of the source.The coverage is more scattered for IR, UV, X-rays, and VHE gamma-ray data that are strongly affected by sparse sampling and often the observations are driven by a high state alert trigger, and therefore may be biased towards high states.In these bands, the coverage had a clear improvement starting from late 2014 (MJD ∼57000).This is the result of an intense MWL and multi-year campaign aimed at a precise monitoring of the source state for the detailed modelling of the source emission and periodicity.From a visual inspection of Figure3we can conclude that the source shows high variability over the years in all bands, with moderate variations in radio and more pronounced variations in the other bands.This behaviour is quite common in HBLs(Acciari et al. 2021).A detailed characterisation of the variability is the key to investigate the physical phenomenology responsible for the broad-band emission as detailed inRieger (2019) and references therein.In the following subsections, several variability studies are presented.The aim is twofold: first, the characterisation of the variability (and periodicity) at different bands.Secondly, the identification of inter-band connections.These connections are a powerful tool to unveil single/multiple regions responsible for the observed emission.

Figure 5 .
Figure5.Selected scatter plots used to investigate flux correlations.For correlation studies involving Fermi-LAT, the simultaneity window has been set to ± 10 days.In all other cases, the window assumed for simultaneity is ±1.5 days.

Figure 8 .
Figure 8. Bayesian block representation of the MWL lightcurves having excluded the observational gaps.

Figure 9 .
Figure 9. PG 1553+113 MWL folded lightcurves obtained by assuming a period of 798 days.In each panel, the average values are reported with a continuous line.The error bar in each bin is the standard error of the mean values.The colored bar in each bin represents the variance of the data distribution.

Figure 10 .
Figure 10.GLS periodgram of the optical (upper panel) and MAGIC (lower panel) data.The optical data show a broad GLS peak, described in the text.In case of MAGIC data, no significant peak emerged from the study.
in long XMM-Newton observations performed in 2015.From the ∼ 1 h variability time-scale, they inferred a size of the emitting region R ≲ δ × 10 14 cm.Recently, Dhiman et al. (2021) confirmed the intraday variability of PG 1553+113 in the 0.3-10 keV band with XMM-Newton data taken during 2010-2018.The authors found an indication of variability in 16 over 19 observations, where the duration of the observations ranged from 21 to 140 ks.The doubling timescale ranged from 2 to 33 ks, i.e. ∼ 30 min to ∼ 9 h.The short duration of our single pointings prevented us to probe intranight variability in >hour timescales, as the one suggested in Dhiman et al. (2021) study.

Table 1 .
Summary of MAGIC telescopes observations.

Table 3 .
Results of the search for periodicity in the different bands.We have analysed the MWL behaviour of PG 1553+113 using MAGIC and MWL data from 2007 until the end of 2017 covering bands from radio to VHE gamma rays.The main motivation of this work was to study if the 2.2-year periodicity seen in the GeV gamma rays by the Fermi-LAT (Ackermann et al. 2015) can be seen in our MAGIC data, and if the MWL data can be used to constrain the models explaining the periodicity.