Multiwavelength Analysis of Fermi-LAT Blazars with High-Significance Periodicity: Detection of a Long-Term Rising Emission in PG 1553+113

Blazars display variable emission across the entire electromagnetic spectrum, with timescales that can range from a few minutes to several years. Our recent work has shown that a sample of five blazars exhibit hints of periodicity with a global significance $\gtrsim2\,\sigma$ at $\gamma$-ray energies, in the range of 0.1~GeV$<$E$<$800~GeV. In this work, we study their multiwavelength (MWL) emission, covering the X-ray, ultraviolet, optical, and radio bands. We show that three of these blazars present similar periodic patterns in the optical and radio bands. Additionally, fluxes in the different bands of the five blazars are correlated, suggesting a co-spatial origin. Moreover, we detect a long-term ($\approx$10 year) rising trend in the light curves of PG~1553+113, and we use it to infer possible constraints on the binary black hole hypothesis.


INTRODUCTION
A blazar is an active galactic nucleus (AGN) with a supermassive black hole (SMBH) launching a relativistic jet towards our line of sight (e.g.Ulrich et al. 1997).Blazar emission is highly variable, spanning the entire electromagnetic spectrum (Urry 2011) and a wide range of timescales, from minutes to years (Urry 1996).Some blazars have periodic multiwavelength (MWL) emission (e.g.Ackermann et al. 2015).However, this temporal behavior is usually interpreted as stochastic (e.g., Covino et al. 2019).
Complementary to periodicity analysis are temporal correlations among the MWL light curves (LCs), which are used to measure time lags between bands.This information helps to constrain the location of the emission region(s), disentangle which jet models are favored, and understand the variability mechanisms (e.g.Cohen et al. 2014;Liodakis et al. 2018).
In Peñil et al. (2022) (P22, hereafter), we reanalyzed 24 blazars exhibiting possible periodicity in Fermi-Large Area Telescope (LAT, Atwood et al. 2009) observations, which were originally presented, as part of a larger sample, in Peñil et al. (2020) (P20, hereafter).From this sample, we selected the five blazars that exhibited the ★ E-mail: ppenil@clemson.edu† E-mail:majello@clemson.edu ‡ NHFP Einstein Fellow most significant periodicity in -ray (≥4.0 of local significance, P22).In this study, we aim to characterize the MWL behavior of these five sources, from radio up to -rays, and to infer its cause.
Employing the same analysis pipeline presented in P22, we search for similar periodic behavior at longer wavelengths.We further perform MWL LC correlation studies and use complementary tools such as fractional variability and structure function to characterize this MWL variability.This paper is the second in a pair of publications.The first publication (hereafter Paper I) focuses on the other 18 blazars from P22 cataloged as low-significance.In Paper I, we employ the same databases and perform the same variability analysis used in this paper.
Regarding periodicity, we used the same methodology as this paper, obtaining 2 sources with a period with ≥3.0 (pre-trial).Moreover, the correlation analysis reveals a high correlation between the -ray, optical, and IR bands with delays < 28 days and the radio band with typical delays of a few hundred days.The structure function shows variability timescales compatible with the periods associated with the emission and shorter timescales, supporting that such variability could originate from instabilities or fluctuations in the accretion disc.
The paper is structured as follows.In Section 2, we introduce our sample of sources and discuss the data collected from Fermi-LAT and MWL sources.Section 3 provides a description of the methodology and analysis methods.In Section 4, we present our results.Section 5 is on the discovery and characterization of new types of variability in PG 1553+113 and PKS 2155−304.Section 6 offers an interpretation of the origin of the variability observed in PG 1553+113.In Section 7, we summarize the main results and conclusions drawn from our study.Additionally, there is an Appendix that presents all the blazars MWL LCs and provides supplementary information on our variability studies.

Source Selection
In this work, we consider five blazars with periods, T, detected in their high-energy (HE, 0.1 GeV<E<800 GeV) -ray emission with global significance (≥2.0,P22).This sample was obtained from the periodicity analysis of 12 years of Fermi-LAT observations, from August 2008 to December 2020.These blazars are listed in Table 1, and their MWL emission are shown in Appendix D.

Multiwavelength Archival Data
We use data from several observatories covering a broad swath of the electromagnetic spectrum.For the X-ray and UV bands, we employ data from the Swift/XRT)1 .The raw count-rate data were taken from the automatic processing of Stroh & Falcone (2013) 2 .Hardness ratio estimations (HR, ratio between soft (0.3−2.0 keV) and hard (2.0−10.0keV) X-rays), also by Stroh & Falcone (2013), were used in conjunction with the Swift/XRT detector response to estimate a photon index, , under the assumption that blazar emission in X-rays can be represented by a simple power-law (PL) A(E)=KE −  (Ghisellini & Tavecchio 2009; Middei et al. 2022).The count rate and HR were then used to estimate a flux at each given epoch.We use data from Swift-UVOT (Ultraviolet and Optical Telescope)3 for the filters 'uvw2' (1928 Å), 'uvm2' (2246 Å) and 'uvw1' (2600 Å)4 .We perform a data reduction of all available Swift-UVOT archival observations to provide the optical-to-UV LCs via the standard pipeline, detailed in Poole et al. (2008).For all UVOT filters, the source regions are selected as circles of 5 ′′ , centered on the source.The background is defined as a circle of 30 ′′ away from the source to avoid contamination.The task uvotsource is employed to extract the magnitudes, which are then corrected for Galactic extinction according to the recommendations in Roming et al. (2008).Finally, the fluxes are derived using the standard zero points listed in Breeveld et al. (2011).
For the optical bands, data from KAIT (Katzman Automatic Imaging Telescope, R-band)5 , CSS (Catalina Sky Survey, V-band)6 , ASAS-SN (All-Sky Automated Survey for Supernovae, V-band)7 , and Tuorla (R-band, Nilsson et al. 2018) 8 are used.We also employ the optical V-and R-bands and near-infrared (IR) I-band data of the American Association of Variable Star Observers (AAVSO) 9 .Data from SMARTS (Small and Moderate Aperture Research Telescope System, Bonning et al. 2012) 10 in the optical B-, V-, and R-bands, and the near-IR J-and K-bands are also used.The Steward Observatory11 provides public optical data in V-and R-bands, with photometric and polarimetric observations obtained approximately simultaneously.
We have combined the V-and R-band data from the optical observatories, and the combinations are herein denoted as "V-band" and "R-band."In some cases, non-calibrated V-and R-band data from the Steward Observatory are used, which are not combined with data from other observatories.We denote this data as "Steward-V" and "Steward-R" in the figures of Appendix D.
Finally, 15 GHz data from OVRO (the Owens Valley Radio Observatory 40-m radio telescope) are used.OVRO has been engaged in a blazar monitoring program supporting the Fermi satellite (Richards  et al. 2011) 12 .These data extend for longer than 12 years, to June 2020, when monitoring ceased, and the data were made publicly available.
For all the MWL data, we use the same -ray LCs and 28-day binning as in P22, where the Fermi-LAT data analysis is described.We search for long-term periodicity (∼years) in the range of [1-6] years.

METHODOLOGY
We apply the pipeline described in P22 to all data sets for the periodicity search.The pipeline includes methods that can manage the gaps: Lomb-Scargle (LSP, Lomb 1976;Scargle 1982), Generalized LSP (GLSP, Lomb 1976;Scargle 1982), Phase Dispersion Minimization (PDM, Stellingwerf 1978), Weighted Wavelet Z-transform (WWZ, Foster 1996), Enhanced Discrete Fourier Transform with Welch's method (DFT-Welch, Welch 1967), and Markov Chain Monte Carlo Sinusoidal Fitting (MCMC Sine, see P20).These methods can provide different results depending on different factors, independently if they are based on the same algorithm.The performance of such methods depends on the sensitivity to the gaps in the LC (see P20), the use or not of the uncertainties of the data to search for the periodicity (e.g., GLSP), and the sensitivity to noise (see P22).
Analyzing non-stationary time series can be challenging (Feigelson et al. 2022).Specifically, a non-stationary time series is characterized in the frequency domain through its time-varying power spectrum (Vaughan et al. 2003).This implies that estimating the power spectral density (PSD) cannot be applied uniformly to the entire time series.In such instances, the underlying physical process driving the variability changes over time, resulting in variations in properties such as the PSD (Vaughan 2013).These random fluctuations, typical of a red noise process, do not yield meaningful physical insights (Vaughan et al. 2003).
One method to ensure stationarity in a time series is detrending.Detrending is a recommended preprocessing step (e.g., Welsh 1999) since a trend can introduce contamination in the low-frequency components, potentially leading to the detection of false periodicities in the time series (McQuillan et al. 2013).This helps in avoiding erroneous results associated with any underlying trends in the data.Nonetheless, detrending can introduce correlations in the data, particularly when data transformations are involved, such as the differencing method, which can introduce autocorrelations if the differenced series still retains some underlying patterns.In our case, we employ linear detrending.Linear detrending can have the effect of seemingly increasing noise correlation in the detrended data, even though it does significantly reduce autocorrelation in the time series.This noise correlation can happen if the linear component is not adequately removed or if the noise in the original data contains some non-random systematic patterns or fluctuations.However, in our specific case, we can rule out the latter scenario, as blazar LCs are known for exhibiting red noise characteristics (e.g., Vaughan et al. 2003).
In our analyses, the methods of our methodology are effective in managing non-stationary LCs.Moreover, we take steps to verify the stationary nature of the data, post-detrending, to ensure that any previous spurious effects associated with non-stationary time series are eliminated.To accomplish this, we utilize the augmented Dickey-Fuller test (Dickey & Fuller 1979) for validation of the stationarity of our data (Feigelson et al. 2018).
We compute the cross-correlations of all the MWL data with the corresponding Fermi-LAT LC for each source.The correlations are performed by the z-transformed discrete correlation function (z-DCF, Alexander 2013).We compute the cross-correlation of contemporary high-flux emission states selected with a Bayesian block algorithm (Scargle et al. 2013).
The PSD, typically described as a power law, quantifies the variation of a time series as a function of frequency, which relates to the physical origin of the emission (Abdo et al. 2010b).We estimate its index with the Power Spectrum Response Method (PSRESP, Uttley et al. 2002).Based on the PSD model, we employ the method based on simulating LCs (Emmanoulopoulos et al. 2013) to estimate the local significance (pre-trial correction, see §3.1).These simulated LCs exhibit identical PSD and probability density function (PDF) characteristics as the original LC, implying that they share the same noise properties as the original LCs.Additionally, these LCs undergo the same detrending process as the original ones.Consequently, the significance is appropriately calibrated for detrending across all the methods we employ.We use the implementation from Connolly (2015) for simulating 20,000 non-periodic LCs.

Global significance
We correct the local significance inferred from the different methods since we did not have prior knowledge of the periods of the potential signal in our periodicity analysis.This correction is implemented using the look-elsewhere effect to estimate the "global significance" (Gross & Vitells 2010).We evaluate the "global significance" by applying the trial factor to the local significance of each periodicity.The trial factor combines the number of independent periods we search for periodicity, 35, and the size of our sample of blazars, 351 (more details in P22).This yields: (i) ≈3.5 for a local significance of 5.5 (ii) ≈2.8 for local significance of ≈5 (iii) ≈1.8 for a local significance of ≈4.5 (iv) <1 for a local significance <4.5.
For PG 1553+113 (see Figure A4), we find a periodicity of ≈2.213 yr in its UV, optical, and radio emission with local significance >2, ≥5, and ≥5, respectively.Periods with local significance <3 are considered non-significant in this work.These results are compatible with the observed periodicity in -rays and consistent with Ackermann et al. (2015), who reported a ≈2.2-yr period in ray, optical and radio bands.In the X-ray band, we find a period of ≈1.5 yr with local significance ≈2.A secondary period compatible with the finding by Huang et al. (2021), of ≈2.2 yr, is inferred by LSP and PDM.
Finally, for PKS 2155−304 (see Figure A5), we find a period of ≈1.7 yr for the X-ray band with local significance ≈2.5.Additionally, we find a period of ≈1.7 yr in the UV and optical bands, a period compatible with that observed in -rays, but not significant (≤2.5 of local significance).However, our result is incompatible with the period reported by Sandrinelli et al. (2016), who reported a period of ≈0.9 yr.Bhatta & Dhital (2020) reported a period of ≈1.7 yr for the optical band.
The local significance of the results typically ranges from 2 to 5, depending on the method.Even for the same dataset, the results in terms of period and significance can differ.The disparity in outcomes can be due to several factors.The shorter time coverage and uneven sampling of the MWL data sets also lead to much larger errors of the derived periods than those from the Fermi-LAT LCs.Specifically, the methods are affected differently by the gaps in the LCs (see P20 for the complete study).For instance, LSP is the most robust against the missing data, and DFT is the most affected method (with differences of 50% in detecting period and significance).Another factor is the impact of red noise.Specifically, each method has a different bias to red noise.For instance, the most robust is the GLSP, and the most sensitive is the DFT, with differences in the detection capacity of 50% (see P22).For this reason, some methods tend to have incompatible periodicities and significance for the same dataset (e.g., S5 0716+714 in the UV band, see §9).

Correlation
Correlation results of PG 1553+113 and PKS 2155-304 are shown in Table 3, including the periods inferred with the z-DCF.The correlation results of PKS 0454−234, S5 0716+71, and OJ 014 are included in Appendix A.2.All delays < ±28 days are compatible with zero lag due to the 28-day binning we use for the Fermi-LAT LCs.Negative lags imply that the -rays are leading the MWL emission (i.e. the -ray flare precedes the other wavelengths), and vice versa.
Several studies have looked at MWL correlations for PG 1553+113.Liodakis et al. (2018) find time lags of ≈10 days with ≈2 (significance before applying any trial correction, as for the other significance reported in this section) between Fermi-LAT and KAIT, and ≈100 days (≈2) between Fermi-LAT and OVRO.In Ackermann et al. (2015), a zero time lag is observed for the optical band and ≈ −100 days for radio.A negative time lag denotes that the -ray emission is leading the radio.A time lag in the radio band of −530 days is claimed by Max-Moerbeck et al. (2014).In this work, we do not obtain any time lag between the -ray, X-ray, and optical bands (with >3), with a delay of radio of ≈ −200 days (>2).
The autocorrelation shows a period in optical and radio of ≈2.2 yr, compatible with -rays, but not significant (see Figure 1).Finally, for PKS 2155−304, Patiño-Álvarez et al. ( 2013) report zero lag in the optical band, in agreement with our results.We also find zero lag correlation with the X-ray band.In the autocorrelation, we obtain a compatible period in X-ray band with respect to the one inferred in -rays, but it is not significant.Similar results are obtained in the optical band.

Correlation results between X-rays and 𝛾-rays
A lack of correlation between X-rays and -rays would signify that different regions are involved in the origin of both emissions.However, the results of Table 3 show correlations approximately compatible with zero lag between -rays and X-rays (all lags < ±28 days are compatible with 0 lag due to the 28-day binning of the Fermi-LAT).These results may imply a common origin of both emissions, considering the uncertainties due to the 28-day sampling.Sikora et al. (2013) discuss two emission scenarios that could help in understanding the disparity in the X-ray and -ray correlations in FSRQs.If the X-ray emission produced by the synchrotron self Compton (SSC) and generated at  > 10 pc dominates over the X-ray emission produced in the BLR, no correlation is expected.On the other hand, a correlation between X-rays and -rays (and, thus, a co-spatial origin) is possible when the Comptonization of hot dust radiation produces the external Compton (EC), dominating the X-ray emission.This takes place at a parsec scale.Regarding the BL Lac objects, the 0-lag correlation of our results may support the model that the -ray and X-ray emissions are likely produced by the same population of relativistic electrons through synchrotron and SSC processes, respectively (e.g., Abdo et al. 2010a).

Long-Term Trend in the MWL Emissions of PG 1553+113
Observing the MWL bands of Figure A4, we can see an increase of the flux (trend) with time in PG 1553+113 in the -ray (Rueda et al. 2022), UV, optical, and radio bands.We estimate each trend by fitting the corresponding LC to a first-degree polynomial function14 .The R-squared (R 2 ) criterion is used to measure the goodness of the fit.We consider a fit reliable when R 2 ⩾75% (Hair et al. 2013).
Examples of this analysis are shown in Figure 2 and Figure 3.The different slopes inferred from the trend analysis for each waveband have compatible values (with a slope ≈2×10 −4 ).We also estimate the average amplitudes of the oscillations in the different bands.Firstly, we normalize the LC with the maximum flux of each LC.Then, we measure amplitudes from the consecutive peak-minimum flux in the different bands.Then, we obtain the averaged amplitude.The normalized amplitudes present compatible values within uncertainties: 1.4±0.2 for  rays, 1.4±0.4 for X rays, 1.5±0.2 in the optical band, and 1.5±0.1 for radio.Finally, our results indicate that the -ray flux (≈2.5×10 −8 ph cm −2 s −1 ) and radio (≈0.15 Jy) begin to increase around the same time, between the years 2011-2013, but in X-ray, UV, and optical bands, the minimum is observed around 2014-2015.

Decreasing Long-Term Trend in the MWL emissions of PKS 2155−304
A trend is also observed in PKS 2155−304, but, in this case, it is a decrease in the flux (observed in the optical band in Zhang et al. 2014).We perform the same trend fit presented in §5.1.The slopes of the -ray and X-ray bands have comparable values of ≈−3×10 −3 .Regarding the optical band, the slope inferred is ≈−2×10 −4 .The normalized amplitudes are comparable values: 1.3±0.1 for -ray, 0.8±0.3 for X-ray, and 1.2±0.1 in the optical band.Finally, the results show that the end of the flux decrease occurs at approximately the same time, around 2014, in all the wavebands.

THE BINARY HYPOTHESIS FOR PG 1553+113
Trends can be associated with red noise, which can mimic periodicity (e.g.Vaughan et al. 2016) or lead to the detection of false periodicity in a time series (McQuillan et al. 2013).Nonetheless, it is instructive to consider interpretations based on the binary hypothesis since there are well-established binary accretion dynamics that might naturally explain such trends.
As the basis for an illustrative discussion on long-term trends in AGN LCs, we consider the binary hypothesis for PG 1553+113.In addition to near-orbital modulations in the jet luminosity, for certain binary parameters (D'Orazio et al. 2013;Miranda et al. 2017;Duffell et al. 2020;Zrake et al. 2021), a longer period is observed in circumbinary accretion simulations (MacFadyen & Milosavljević 2008).This longer period is set by the orbital period associated with an overdensity (the "lump") in the circumbinary cavity wall (see the description in, e.g.Miranda et al. 2017).For year-like orbital periods appropriate for PG 1553+113, the lump period is measured to be ≈ 5 − 10 orbital periods in simulations with baroclinic thermodynamics and radiative cooling (e.g.Farris et al. 2015;Westernacher-Schneider et al. 2022).If the long-term MWL trends for PG 1553+113 are due to a lump, then visual inspection of the V-band currently indicates that the lump period is at least roughly eight orbits -see Figure A4.Whether the observed long-term trends are due to a lump implies different constraints on the binary parameters described in this section.
Firstly, we can derive a constraint on the binary eccentricity independently of the interpretation of the long-term MWL trends by showing that the binary is likely in the gravitational wave-driven (GW-driven) regime.A mass estimate of (0.4 − 8) × 10 8  ⊙ for the central black hole in PG 1553+113 was made based on LC variability (Dhiman et al. 2021).Holding a binary SMBH's semi-major axis fixed, the binary is further from the GW-driven regime for lower total mass, higher accretion rate, lower eccentricity, and lower mass ratio (Westernacher-Schneider et al. 2022).Thus, we can attempt to place the binary outside of the GW-driven regime by assuming its mass is the lower end of the above estimate 0.4 × 10 8  ⊙ , it accretes at the Eddington limit (with a fiducial radiative efficiency of  = 0.1), and it is circular eccentricity ( = 0).We must also assume a fiducial accretion "eigenvalue"  (we choose  = 1) (Paczynski 1991;Popham & Narayan 1991), and we assume that the periodicity obtained in the LCs is roughly the redshifted orbital period ≈ 1.5 × (1 + ) yr for  = 0.433.Doing so (following Westernacher-Schneider et al. 2022), we find the binary is in the GW-driven regime for all mass ratios  ≳ 0.02. 15Simulations show that mass ratios of accreting binaries are driven rapidly upward from such low values as  = 0.02 (e.g.Duffell et al. 2020).Thus, the hypothetical binary in PG 1553+113 likely has  > 0.02, and is, therefore, likely in the GW-driven regime.
The fact that the hypothetical binary in PG 1553+113 is likely in the GW-driven regime, combined with simulation results, implies a constraint on the binary eccentricity.Simulations suggest an equilibrium value of binary eccentricity in the gas-driven regime of  ≈ 0.4-0.6 (Roedig et al. 2011;Roedig & Sesana 2014;Zrake et al. 2021;D'Orazio & Duffell 2021), with the most recent results pointing to  ≈ 0.4.Since GW emission circularizes binaries, we expect the putative binary in PG 1553+113 to therefore have eccentricity  ≲ 0.4.This is our first constraint on the supposed binary.Note that simulations also predict that lump periodicity manifests when binaries have mass ratios  ≳ 0.2 (e.g.D 'Orazio et al. 2013;Duffell et al. 2020) and eccentricities  ≲ 0.1 (e.g.Miranda et al. 2017;Zrake et al. 2021).
We consider the following three scenarios for a hypothetical binary in PG 1553+113: (i) The long-term MWL trends are due to a lump.(ii) The long-term MWL trends are not due to a lump because the lump is not present in the system.
(iii) The long-term MWL trends are not due to a lump, even though the lump is present in the system.
In the first scenario, a lump must exist; thus, the constraints on the binary from simulations are  ≳ 0.2 and  ≲ 0.1.In the second scenario, the lump does not exist; thus, we either have a constraint on the mass ratio ( ≲ 0.2) or a constraint on the eccentricity ( ≳ 0.1), or perhaps both.In either case, since the binary is likely in the GWdriven regime, the eccentricity is also constrained by  ≲ 0.4.All of the constraints in this second scenario are consistent with the binary model proposed by Cavaliere et al. (2017) ( = 0.1,  = 0.2), which was arrived at on the basis of different considerations.In this case, the long-term MWL trends are left unexplained and are presumably stochastic in the absence of another plausible mechanism.In the third scenario, lump periodicity does not transmit to jet variability in sufficient measure.
Two mechanisms by which lump periodicity could imprint on blazar emission are via transmittance of lump periodicity to mass accretion rates (a phenomenon which depends upon disk parameters; compare the simulations from, e.g., Farris et al. 2015;Westernacher-Schneider et al. 2022), and periodic modulations in the supply of seed photons from the disk.In the accretion rate mechanism, one expects a systematic modulation of jet processes, manifesting in the SED via both the low-energy synchrotron component (e.g. via modulation in the supply of electrons to the jet) and, consequently, the high-energy SSC component.In the seed photon mechanism, seed photons from the disk, modulated on the lump period, undergo inverse-Compton scattering in the jet(s) and/or disk corona ("external Compton" (EC, e.g.Band & Grindlay 1986;Levinson & Blandford 1996), thereby resulting in a modulation of the EC component of the SED on the lump period.The third scenario enumerated above would require that both of these imprint mechanisms are inefficient relative to other mechanisms of long-term variability.It is important to note that the SED of a BL Lac blazar like PG 1553+113 is often fit without an EC component (see, e.g., Aleksić et al. 2012), so if an EC component is non-negligible in the high-energy SED, then it must appear sufficiently degenerate with the SSC bump.
On the other hand, the first scenario enumerated above instead seems to require that both of these imprint mechanisms are efficient and comparably so.This is because the minimum of the long-term trend in optical, apparently occurring near 2015, does not appear to correspond to a minimum in -rays (see Figure A4), although a longer temporal baseline is needed for higher confidence.Simulations have shown that there is a lag of order ≈ 20 − 30% of a lump period between lump modulations in thermal disk emission and accretion rates (see Figure 5 in Farris et al. 2015).Supposing that high-energy emission receives non-negligible contributions from both the SSC and EC components and that lower-energy emission does not receive a significant contribution from the EC component, one may therefore expect the lag between accretion and the external supply of seed photons could manifest as a multi-year shift of the minimum in the long-term MWL trend in -rays with respect to e.g.optical & UV.
Future theoretical work should address the efficiency of the mechanisms by which lump periodicity can imprint on blazar SEDs.If the long-term MWL trends in PG 1553+113 are due predominantly to lump periodicity, then we expect the recent upward trends reverse over the next few years, which is a strong motivation to continue monitoring PG 1553+113 across the entire electromagnetic spectrum.

SUMMARY
In this work, we have implemented a variability study of the MWL (radio, IR, optical, UV, and X ray) emission of the five blazars identified as having evidence of periodicity in Peñil et al. (2022).We find that two of them, PG 1553+113 and PKS 2155−304 show similar periodic behavior in their MWL emissions.The correlation  rayoptical does not indicate any lags with the limitation of 28 days due to the sampling of the Fermi-LAT LCs.This result suggests a common spatial origin of both emissions.In the  ray-radio correlation, 0 and ≈200 days of delay have been observed in two blazars.No significant correlation above 3 is observed between -rays and X-rays.Regarding PG 1553+113, we made the first analysis of a long-term trend of increasing/decreasing flux in all bands.We explored an interpretation in terms of a hypothetical supermassive black hole binary central engine.

ACKNOWLEDGEMENTS
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 We also acknowledge support from NSF grant AST-1715661.
We also want to thank all the observatories from which we used data.We thank the Las Cumbres Observatory and its staff for their continuing support of the ASAS-SN project.ASAS-SN is supported by the Gordon and Betty Moore Foundation through grant GBMF5490 to the Ohio State University, and NSF grants AST-1515927 and AST-1908570.Development of ASAS-SN has been supported by NSF grant AST-0908816, the Mt.Cuba Astronomical Foundation, the Center for Cosmology and AstroParticle Physics at the Ohio State University, the Chinese Academy of Sciences South America Center for Astronomy (CAS-SACA), the Villum Foundation, and George Skestos.The AAVSO database: Kafka, S., 2021, Observations from the AAVSO International Database, https://www.aavso.org.The CSS survey is funded by the National Aeronautics and Space Administration under Grant No. NNG05GF22G was issued through the Science Mission Directorate Near-Earth Objects Observations Program.The Catalina Real-Time Transient Survey is supported by the U.S. National Science Foundation under grants AST-0909182 and AST-1313422.This paper has made use of upto-date SMARTS optical/near-infrared light curves that are available at www.astro.yale.edu/smarts/glast/home.php.Data from the Steward Observatory spectropolarimetric monitoring project were used.This program is supported by Fermi Guest Investigator grants NNX08AW56G, NNX09AU10G, NNX12AO93G, and NNX15AU81G.This research has made use of data from the OVRO 40-m monitoring program (Richards et al. 2011), supported by private funding from the California Institute of Technology and the Max Planck Institute for Radio Astronomy, and by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G and NSF grants AST-0808050 and AST-1109911.We also acknowledge the use of public data from the Swift data archive.
Table 2. List of periods and their associated uncertainties (superscript), local significances (left subscript), and global significances (right subscript) of PG 1553+113 and PKS 2155−304.The symbol "★" is utilized to represent two periods that have been derived from the same energy band and exhibit similar significance.These periods are organized based on their peak amplitudes.The symbol † is used to denote the periods where the PDM results exhibit the harmonic effect as described in P20.The symbol ‡ is used to indicate periods that are consistently presented in the WWZ for all the LC.The MCMC sine fitting results provide information solely about the inferred period and its associated uncertainties.The period values are expressed in years.1.9 ±0.5

Association
1.6/0.0 1.9 ±0.5 1.4/0.0 1.9 Table 3. Results of the z-DCF cross-correlation analysis for PG 1553+113 and PKS 2155−304.The table also shows the period inferred from the cross-correlation with the -ray LC and the auto-correlation.These periods are associated with the significance levels of the oscillations (high-level, low-level, see Figure 1).We report one significance when the high and low significance levels are the same.The correlation with the  rays is limited by the 28 days sampling of the Fermi-LAT LCs.

A.1 Periodicity
The results of the periodicity analysis are shown in Table A1.PKS 0454−234 has no significant evidence of periodicity inferred from the analysis in the available bands.However, the optical data have insufficient temporal coverage (≈7.0 yr, see Figure A6) to reliably rule out the period observed in -rays (≈3.6 yr).
For S5 0716+71 (see Figure A7), we find a period compatible with that observed in -rays and in the V-band, ≈2.7 yr (similar to Bhatta 2021).A period of ≈1 yr is also inferred by some of the employed methods.This is compatible with the results from P22, where both periods were obtained in -rays.However, this V-band period of ≈2.7 yr is not significant (<3 of local significance).Sandrinelli et al. (2017) report no significant evidence of periodicity in the Rband.We do not observe significant periodicity in X-ray and radio bands.
We do not find any evidence of periodicity for OJ 014 in the radio band.The periodicity analysis for the optical band was not performed due to insufficient temporal coverage of the data (Figure A8).

A.2 Correlation
Correlation results are shown in Table A2, including the periods inferred with the z-DCF.
For PKS 0454−234, Cohen et al. (2014) claim a delay between -ray and the optical bands of ≈150 days; instead, we find a zero lag with local significance of >4.No period is inferred from the autocorrelation, in agreement with the results in §A.1.
Regarding S5 0716+71, we find a lag of ≈0 days in the X-ray, optical, and radio LCs (≥2).Additionally, the cross-correlation analysis suggests a ≈2.7 yr period for both optical and radio bands, compatible with the one from P22, but not significant.The autocorrelation also suggests a compatible period for the V-band, but not at a significant level.
OJ 014 presents no time lag between its -ray and optical emission.For radio, the lag is −40 days (>4).Both correlations show a compatible period with -rays but it is not significant.
In some cases, the results of Table A1 and the autocorrelations are significantly different (e.g., S5 0716+71).There are several factors that contribute to the divergent results.The shorter time coverage and uneven sampling of the MWL data can lead to much larger errors in the derived periods than those derived from the Fermi-LAT LCs, including resulting in the absence of a period (note that all of the autocorrelations for -rays are compatible with those reported in P22).Additionally, the methods handle LC gaps differently (see P20).Finally, each method is impacted differently by the choice of binning in each analysis (periodicity and autocorrelation Otero-Santos et al.

B Variability Study
LC variability is studied using the fractional variability (  ), the structure function (SF), as well as the PSD and PDF.  and the SF are respectively used to quantify the variability (Vaughan et al. 2003) and measure the characteristic variability timescales. var is affected by the time coverage, sampling, and binning of the data (see Vaughan et al. 2003;Schleicher et al. 2019).Consequently, employing a similar time window and binning to analyze MWL data.Additionally, using nonsimultaneous data can also lead to inconsistencies in  var (Schleicher et al. 2019).
Finally, We also analyze the polarization of blazar emission using polarimetric data from the Steward Observatory.

B.1 Fractional Variability
We have evaluated the amount of variation displayed by our blazar sample through the evaluation of  var .The results are presented in Table A4.The  var values, typically >0.25, prove the variable nature of the sources studied here.We also observe that the most variable source of the sample (that with the highest  var in the different bands) is the only FSRQ included in this study, PKS 0454−234.On the other hand, BL Lac objects tend to have smaller values of  var (e.g., Bhatta & Dhital 2020).In fact, for the four BL Lac objects included here, the two showing the synchrotron peak at lowest frequencies (OJ 014 and S5 0716+714) are those with the highest  var among the BL Lacs, while PG 1553+113 and PKS 2155−304, with higher synchrotron peak frequencies, are less variable in our study.Despite the small number of sources studied here, this trend of increasing  var with the decreasing value of   is in line with the results reported in other studies (e.g., Bhatta & Dhital 2020).Hence, FSRQs and low-synchrotron-peak BL Lacs are typically more variable than high-synchrotron-peak BL Lacs.
We note that the calculation and comparison of the MWL  var might be highly affected by the non-simultaneity of the data between the different bands (Schleicher et al. 2019).Therefore, as a crosscheck, we have calculated  var using only the simultaneous data in all the bands. var derived for simultaneous data (i.e., taken within the same 24-hour period) yields values with no significant differences in comparison to  var obtained with all the data of the LCs (see figures of §D), except for the case of the X-ray  var values of PG 1553+113 and S5 0716+714, which  var of the simultaneous data are ≈15% higher than  var obtained with all the X-ray data (see Figure A1).Nevertheless, the evolution with the frequency of  var remains unchanged with both approaches.
By studying the evolution with the frequency of  var in the different bands, we can also retrieve information on the processes and particle populations dominating the variability.However, for the five sources analyzed here, the MWL coverage is not always optimal for evaluating the evolution of  var with the frequency.Nevertheless, we observe the minimum variability for PG 1553+113 at radio wavelengths (see Figure A1), increasing up to its maximum in the X-ray band.Then,  var decreases in the HE -ray band.We do not have data in the very-high-energy (VHE, E>100 GeV) -ray band.We note that Aleksić et al. (2015a) report for VHE an  var close to that in the X-ray band.This could point towards a double-peaked shape, observed in the past for several sources (Aleksić et al. 2015b).Since this source is a high-synchrotron-peak BL Lac object with its X-ray emission corresponding to the end of the synchrotron bump, a higher variability in this regime and in the VHE -ray band could indicate a higher variability of the high-energy electron population responsible for the emission.In comparison, we can see the case of S5 0716+714 (see Figure A1), where the maximum, ≈1.2, is found at UV wavelengths, in comparison with its minimum in the X-ray band.In this case, the source is an intermediate-synchrotron-peak BL Lac, and its X-ray emission is found in the transition between the low-and highenergy SED peaks.Hence, depending on  var structure and SED characteristics, we can observe differences between sources that may reveal information regarding the particle populations responsible for the emission and variability observed.

B.2 Characteristic timescales
The results of Table A4 reveal different variability timescales, ranging from ≈100 days to ≈4 yr.The emission is most stable in the radio.The most variable is in the optical, as inferred by the set of timescales derived from the structure functions (SF, see Figure A2).
The average slope of the SFs analyzed in this work is 0.5±0.1.This result is compatible with a slope of 0.44±0.03associated with a model based on instabilities in the accretion disk, as derived by Hawkins (2002).Regarding the periodicity, most of the periods inferred with the SF (see Table A4) are compatible with the ones shown in Table 2 and Table A1.

B.3 Power Spectrum Index Estimation
The results of the power spectrum index estimation are shown in Table A3.The slopes obtained for the optical and radio bands range between 1.2 and 1.5, compatible with slopes of the -ray emission presented in P22.This index range can be associated with an admixture of flickering (pink) noise, consisting of short-scale variability and red noise associated with long-term variability.According to Abdo et al. (2010b), these oscillations are associated with instability and turbulence in the accretion flow through the disk or the jet.

B.4 Polarized Light
Blazars are characterized by high levels of radio and optical polarization.This polarized emission is also extremely variable, and it depends on the structure of the magnetic field of the emitting region.Thus, it can provide valuable information about the origin of the emission in blazars.Raiteri et al. (2013) interpret the long-term flux variability of the polarized optical emission according to two different models: helical magnetic fields and transverse shock waves.Both models predict two different behaviors of the polarization degree w.r.t. the viewing angle  (see Figure A3).In both models, when  <  max (with  max ≈ 1/Γ the angle of maximum polarization and Γ the bulk Lorentz factor of the jet), the polarization degree increases with the viewing angle of the jet.The opposite is expected when  >  max (distinguishing two regions, see Figure A3).Additionally, the observed flux increases for decreasing the viewing angle (Raiteri et al. 2013).Therefore, in the first (second) region, the polarization degree is anti-correlated (correlated) with the observed flux.Consequently, the correlation between the polarized and total emission can constrain the scenarios responsible for the variability.Specifically, inferring a correlation for viewing angles <  max between the polarized and total flux could indicate that the variability is not associated with internal-jet processes but produced by external phenomena affecting the jet (e.g., the accretion disc, or the single/binary black hole system, Raiteri et al. 2013).In Figure A3, the polarization degree as a function of the viewing angle is represented for the helical magnetic field model (the transverse shock model has a similar behavior, see Figure 17 in Raiteri et al. 2013).We use Γ = 12.2 and  obs = 3.0 • for PKS 0454−234 (Ghisellini et al. 2009), and Γ = 10.3 and  obs = 5.2 • for S5 0716+714 (Hovatta et al. 2009).Observing the plot for both blazars, polarization increases as the viewing angle increases from =0 • to  max , and then slowly decreases for higher viewing angles.
In this context, we evaluate the correlation of V-band emission of PKS 0454−234 and S5 0716+714 between the polarized and total flux (since only these sources have available data in the Steward database).For PKS 0454−234, both emissions are correlated, with a time lag of 4.9 ± 11.9 days (>2, significance before applying any trial correction, as in the other significance reported in this section).A similar result is observed for S5 0716+714, where a Table A1.List of periods and their associated uncertainties (superscript), local significances (left subscript), and global significances (right subscript) of PKS 0454−234, S5 0716+714, OJ 014.The symbol "★" is utilized to represent two periods that have been derived from the same energy band and exhibit similar significance.These periods are organized based on their peak amplitudes.The symbol † is used to denote the periods where the PDM results exhibit the harmonic effect as described in P20.The symbol ‡ is used to indicate periods that are consistently presented in the WWZ for all the LC.The MCMC sine fitting results provide information solely about the inferred period and its associated uncertainties.The period values are expressed in years.

Association
Table A3.List of power spectral indices inferred with the PSRESP method for the MWL LCs of the blazars.The uncertainty and the probability of each index are also shown.The 'X' denotes that the PSRESP method does not converge to a specific value.The "success fraction" indicates the goodness of fit of PSRESP by giving the probability of a model being accepted.The power spectral indices of the  rays are obtained from P22.

Figure 1 .
Figure 1.z-DCF autocorrelation of PG 1553+113 (V band).The figure shows a period of ≈2.2 yr. denotes the time lag in days.

Figure 2 .
Figure 2. Trend decomposition of the -ray emission of PG 1553+113.R 2 is the R-squared criterion to measure the goodness of the fit.

Figure 3 .
Figure 3.Long-term trend decomposition of the V-band emission of PG 1553+113.
Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare 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 and the Swedish National Space Board in Sweden.Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France.This work performed in part under DOE Contract DE-AC02-76SF00515.P.P. and M.A. acknowledge funding under NASA contract 80NSSC20K1562.S.B. acknowledges financial support by the European Research Council for the ERC Starting grant MessMapp, under contract no.949555.A.D. is thankful for the support of the Ramón y Cajal program from the Spanish MINECO, Proyecto PID2021-126536OA-I00 funded by MCIN / AEI / 10.13039/501100011033, and Proyecto PR44/21-29915 funded by the Santander Bank and Universidad Complutense de Madrid.L.M. acknowledges that support for this work was provided by NASA through the NASA Hubble Fellowship grant #HST-HF2-51486.001-Aawarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555.J.O.S thanks the support from grant FPI-SO from the Spanish Ministry of Economy and Competitiveness (MINECO) (research project SEV-2015-0548-17-3 and predoctoral contract BES-2017-082171) and financial support from the Spanish Ministry of Science and Innovation (MICINN) through the Spanish State Research Agency, under Severo Ochoa Programme 2020-2023 (CEX2019-000920-S) and the project PID2019-107988GB-C22. J.O.S. acknowledges financial support from the Severo Ochoa grant CEX2021-001131-S funded by MCIN/AEI/ 10.13039/501100011033.

Figure A2 .
Figure A2.Structure function of the UVOT LC of PKS 2155−304.The vertical line shows the dip that denotes the presence of a period in the signal.

Figure A8 .
Figure A8.MWL light curves of OJ 014.From top to bottom: Fermi-LAT (E > 0.1 GeV), and V-band (combination of CSS and ASAS-SN) light curves.

Table 1 .
Our sample of blazars given by their association name, AGN type, redshift, period (in years), local significance, and global significance obtained by P22.The AGN types are flat-spectrum radio quasar (fsrq) and BL Lacertae (bll).

Table A4 .
Results of the fractional variability and the periods inferred from the structure functional (see §B).