Multiwavelength Variability Analysis of Fermi-LAT Blazars

Blazars present highly variable $\gamma$-ray emission. This variability, which can range from a few minutes to several years, is also observed at other wavelengths across the entire electromagnetic spectrum. We make use of the first 12 years of data from the Fermi Large Area Telescope (LAT), complemented with multiwavelength (MWL) archival data from different observatories and facilities in radio, infrared and optical bands, to study the possible periodic emission from 19 blazars previously claimed as periodic candidates. A periodicity analysis is performed with a pipeline for periodicity searches. Moreover, we study the cross-correlations between the $\gamma$-ray and MWL light curves. Additionally, we use the fractional variability and the structure function to evaluate the variability timescales. We find five blazars showing hints of periodic modulation with $\geq$3.0$\sigma$ ($\approx$0$\sigma$ post-trials), with periods ranging from 1.2 to 4 years, both in their $\gamma$-ray and MWL emission. The results provide clues for understanding the physical mechanisms generating the observed periodicity.


INTRODUCTION
Supermassive black holes (SMBHs) are found in the centers of almost all galaxies.Accretion of matter onto SMBHs powers some of the most luminous sources in the Universe known as active galactic nuclei (AGNs, Urry & Padovani 1995).Some of these AGNs produce powerful, highly collimated relativistic jets.When these jets are aligned with the line of view of the observer, the jets' emission is relativistically boosted from radio to  rays (Ghisellini et al. 1998), and these AGNs are classified as blazars.One of the main features of blazars is their strong variability at all wavelengths and at different timescales (Fan 2005).
This variability is generally interpreted as arising due to stochastic and unpredictable processes (e.g.Ruan et al. 2012).However, several studies have claimed the detection of quasi-periodic signals coming from blazars (see, e.g.Ackermann et al. 2015;Sandrinelli et al. 2016;Peñil et al. 2020).This latter scenario has grown in interest over the last decades due to the substantial implications derived from periodic emissions in these objects.It can provide information about the nature of the source and the physical processes involved in the most violent environments.Quasi-periodic oscillations in blazars ★ E-mail: ppenil@clemson.edu† E-mail: joteros@iaa.es‡ E-mail: majello@clemson.edu§ E-mail: sara.buson@uni-wuerzburg.de have been interpreted, for instance, in the framework of a binary SMBH system (e.g.Sobacchi et al. 2017).This scenario has been proposed especially for the two best candidates: PG 1553+113 (Tavani et al. 2018) and OJ 287 (Sillanpaa et al. 1988;Valtonen et al. 2011).Nevertheless, periodicity studies and detections of periodic oscillations in blazars are still a controversial topic (see, e.g.Covino et al. 2019;Yang et al. 2021).The impact of the noise in data can generate stochastic uncertainty in the periodicity search, which can provoke the detection of fake periodicity (e.g., Vaughan et al. 2016;Liu et al. 2019).
The recent study from Peñil et al. (2020) (P20 hereafter) presents a systematic search for periodicity in blazars observed by the Fermi Large Area Telescope (LAT, Atwood et al. 2009).As a result, a sample of 24 (out of 351 analyzed, more details in P20) blazars present hints of periodic emission with confidence between 2 and ∼4 (pre-trials).In Peñil et al. (2022, hereafter P22), these 24 blazars were reanalyzed with three more years of Fermi-LAT observations, for a total of 12 years.Here, we analyze the low-significance 19 candidates with a period in -ray of <4 (pre-trials) reported in the latter study as the first publication of two papers, extending the analysis to the MWL data collected for the 24 periodic candidates.We report the results for PG 1553+113, PKS 2155-304, S5 0716+714, OJ 014, and PKS 0454-304 in a separate study (Peñil et al. 2023) since these sources stand out from the sample in showing higher significance in comparison with the rest of the sample for the quasiperiodicity (with pre-trials significance of ≈4 in the period observed in  rays).
The paper is structured as follows.In Section 2, we introduce the sample of sources analyzed in this work and the Fermi-LAT and MWL data collected.In Section 3, we describe the methodology and analysis tools.In Section 4, the results of the analysis are presented.Section 5 provides a possible interpretation of the results.Section 6 summarizes the main results and conclusions.Finally, we also include an appendix in Section 7 where all the MWL LCs of the blazars are presented.

𝛾-ray Sample
In P20, the analysis of 351 blazars detected by Fermi-LAT in 9 years of data (0.1 GeV<E<800 GeV) was carried out.This analysis was performed with a pipeline consisting of several methods (see §3), leading to the discovery of 24 objects with evidence of periodic emission with confidence between 2 and ∼4 (computing the median significance across the analysis methods to sort the blazars).These values refer to the "local significance" since no trial-factor corrections were included in P20.Therefore, no "global significance" was reported.These 24 blazars are reanalyzed in P22, updating the sample with 12 years of Fermi-LAT data at energies >0.1 GeV, where the evidence of periodicity is confirmed for 5 of them with a significance (pre-trials) of >3 (again, using the same pipeline and the median significance across the analysis methods).The "global significance" of these blazars is ≳2.0.
The remaining 19 objects studied here reported a local significance (pre-trials) between 1.0 and ∼3.0.Their global significance is ≈ 0 (see details in § 4.1.1).This global significance does not allow us to claim any periodicity.However, we will search for the same periods in the MWL bands observed in  rays.These results allow us to select the most promising in order to follow them up to confirm or rule out the periods inferred.
The list of the selected blazars can be found in Table 1.We use the Fermi-LAT light curves (LCs) from P22, complemented with MWL data from several observatories listed in Section 2.2.An example of the LCs used in this work is shown in Figure 1 for S4 0814+42.

Multiwavelength Data
For the MWL analysis, we employ archival data from different databases and observatories, covering most of the electromagnetic spectrum.In the X-ray bands, we use data from the Neil Gehrels Swift Observatory, specifically Swift-XRT (X-Ray Telescope) 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 opticalto-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).
The optical data are obtained from KAIT (Katzman Automatic Imaging Telescope, R band, Li et al. 2003) 5 , CSS (Catalina Sky Survey, V band, Drake et al. 2009) 6 , and ASAS-SN (All-Sky Automated Survey for Supernovae, V band, Shappee et al. 2014;Kochanek et al. 2017) 7 .Data from the Tuorla blazar monitoring program8 (R band, Takalo et al. 2008) are also included, using the data from Nilsson et al. (2018).We also employ the optical V-and R-band data from the American Association of Variable Star Observers (AAVSO) 9 International Database.We use data from SMARTS (Small and Moderate Aperture Research Telescope System, Bonning et al. 2012) 10 in the optical B, V, and R bands, and near-infrared (IR) J and K bands.Finally, we also retrieve data from the Astronomy & Steward Observatory11 (Smith et al. 2009), with optical V-and R-band observations.For the analysis, we combine all the data from different observatories in the V and R bands, denoted in the tables and figures as "V-band" and "R-band".Different calibrations from different observatories can introduce offsets that affect the data.To evaluate and correct for the presence of such an effect, we compare all the simultaneous data from the different databases.If the simultaneous data show a systematic offset, the mean difference between each simultaneous pair of measurements is calculated.This estimation is used as the mean offset to correct the data from different observatories and ensure their compatibility.
Since we are interested in studying the long-term evolution of these data series, this scaling is irrelevant to the analysis performed in this work.Some blazars use non-calibrated V-and R-band data from the Steward Observatory.These data can be used for analyzing long-term time series.However, they are not combined with the rest of the LCs of the corresponding band.These data are denoted in the tables as "Steward-V" and "Steward-R", respectively.
Moreover, we also include radio data in the 15 GHz band for those sources observed by the 40m radio telescope from Owens Valley Radio Observatory (OVRO) as part of their blazar monitoring program (Richards et al. 2011a) 12 .These data extend for longer than 12 years, to June 2020, when monitoring ceased, and the data were made publicly available.Radio data are available for 5 of the 19 sources studied in this work.Data with a signal-to-noise ratio lower than 3 are considered non-detection and excluded from the analysis to avoid such data.To perform reliable periodicity and correlation analyses, we need to balance reducing the uncertainties introduced by the noisiest observations and maintaining sufficient sampling and time coverage in the LCs.Otherwise, the analysis could become unreliable or result in a poor estimation of the period/correlation (Otero-Santos et al. 2023).
As we are interested in yearly-scale variations, we use 28-day bins for the MWL LC periodicity analysis, matching the bin size from the Fermi-LAT -ray LCs.The data is organized into 28-day bins, a crucial step that allows us to explore variations beyond the shortterm fluctuations, including those occurring within a single day or week.Our binning method relies on the median value within each bin, which has proven to be effective in managing data with erratic noise, as demonstrated in previous studies (e.g., Rani et al. 2013;Negi et al. 2023).This approach strikes a balance between making the analysis computationally feasible while still maintaining sensitivity to longer-term variations, which may extend over the course of a year.On the other hand, the correlation with the -ray LCs (to have an accurate estimation of the time lags) and the variability analyses (the exception is the fractional variability, see §3.3 for details) are performed using the original sampling of the LCs.

Periodicity
The MWL study presented here is focused on two aspects: (1) periodicity and correlation search and (2) variability.The periodicity search is performed using the pipeline described in P20 and P22, to which we refer the reader for further details.This pipeline searches for periodic variability using multiple methods, which are listed below: (i) Lomb-Scargle periodogram (LSP, Lomb 1976 The selection of these methods for our analysis is based on their strong performance in handling unevenly sampled data, particularly LSP, GLSP, and WWZ, which were specifically developed to address this challenge.The LSP, GLSP, and DFT methods are all Fouriertransform-based, making them well-suited for detecting sinusoidallike signals even in data with irregular sampling.WWZ, while also rooted in the Fourier transform, offers the additional advantage of decomposing the signal into both time and frequency domains, enabling an evaluation of the persistence of potential periods.PDM, in contrast, relies on the dispersion of LC data for various periods, making it a reliable choice for detecting non-sinusoidal periodic patterns, including repeating flares.The MCMC Sine fitting method stands apart from the others as it solely relies on the LC series and searches for sinusoidal signals without encountering the typical limitations associated with Fourier-transform methods, such as restricted frequency resolution, aliasing, and the presence of spurious peaks. As demonstrated in prior studies P20, P22, and Otero-Santos et al. ( 2023), these methods can also exhibit different sensitivities to data gaps or red noise.By employing a combination of these methods, we aim to complement each other's strengths and weaknesses, ultimately providing more robust results when there is a consensus among multiple methods.
Another challenge in periodicity analysis, as highlighted in Feigelson et al. (2022), pertains to non-stationary LCs.This means that estimating the power spectral density (PSD, a measurement of the power in a signal as a function of the frequency) cannot be uniformly applied to the entire time series because the source of the variability changes over time (Vaughan 2013).To address this issue, we employ a detrending step.Detrending is a recommended preprocessing procedure (e.g., Welsh 1999) to mitigate contamination that could otherwise lead to the incorrect inference of false periodicities (Mc-Quillan et al. 2013).In our study, we opt for linear detrending, which can increase noise correlation, especially if the linear component is not effectively removed or if the original data contains non-random systematic patterns or fluctuations.However, we can rule out the latter case, as our blazar LCs are characterized by exhibiting red noise characteristics (e.g., Vaughan et al. 2003).
The methods employed in our methodology are indeed effective in handling non-stationary LCs.Additionally, we take measures to ensure that the data becomes stationary after the detrending process.
To achieve this, we utilize the augmented Dickey-Fuller test (Dickey & Fuller 1979).This test serves as a validation step to confirm the stationarity of the data (Feigelson et al. 2018).
To perform the periodicity analysis, we perform a binning of 28 days in the MWL LC to be consistent with the -rays LCs (see P20 and P22).We search for long-term periodicity (∼years), periods in the range of [1-6] years.
To obtain the local significance of the correlation and periodicity analysis, we employ the method based on simulating LCs (Emmanoulopoulos et al. 2013).These simulated LCs have the same sampling, PSD, and probability distribution function (PDF) as the original LC.Thus, they will have the same statistical properties as the original data set, and accurately modeling the underlying type of noise of the real data.We calculate the PSD of each LC and obtain the spectral index that describes the derived PSD associated with each data set.The PSD denotes the energy variation as a function of Table 1.List of the blazars analyzed here in descending order of the median significance (pre-trials) of the quasi-periodic hypothesis, estimated from the results presented in Table 2.Note that this median significance does not have an actual statistical meaning; it is used as an arbitrary way of combining all of the significance for a given source, sorting the candidates, and comparing with the results obtained in P22 for the -rays.The blazars are characterized by their Fermi-LAT name, coordinates, blazar class, the blazar type according to the frequency of the synchrotron peak (LSP: low synchrotron peaked, ISP: intermediate synchrotron peaked, HSP: high synchrotron peaked), redshift, association name.Additionally, we include the average period (in years) with the uncertainty and the local significance obtained in P22 (as the average of the periodicity) and the corresponding global significance estimated in P22.Note that some sources have two significant periods (organized by the amplitude of the peak), denoted by ★. (1) ( the frequency of the time series and is commonly modeled as a PL function (PSD ∝  −  ).We estimate the PLs describing the different PSDs using the Power Spectrum Response Method (PSRESP, Uttley et al. 2002).PSRESP 13 also provides the uncertainty of the estimated slope and the "success fraction" as a measurement of the goodness of the fit.This "success fraction" estimates the discrepancy between the data and the fit for each scanned slope value, leading to the slope that best reproduces the derived PSD.To estimate the power spectrum index, we simulate 1,000 LCs using the approach from Timmer & Koenig (1995), with the same observational properties of the original LC, i.e., mean, standard deviation, flux PDF distribution (Shah et al. 2020).Ait Benkhali et al. (2020) show that both methods of Timmer & Koenig (1995) and Emmanoulopoulos et al. (2013) yield similar results.We sample the frequency space using a binning of 0.005 and scan slope values between 0.8 and 2.0 using a binning of ∼0.02-0.03.An example of a PSD fit using the PSRESP method can be seen in Figure 2.With these PSD slope and PDF estimations, we employ the implementation from Connolly (2015) for simulating 20,000 non-periodic LCs. 13 We use the implementation of https://github.com/wegenmat-privat/psresp

Correlation
The search for correlations can be performed using either one signal (auto-correlation) or two different signals (cross-correlation). Physically, auto-correlation (cross-correlation) corresponds to a measurement of the similarity of a signal with itself (with respect to a second signal).
The Discrete Correlation Function (DCF) is the traditional approach to search for correlations (Edelson & Krolik 1989).However, as described in Alexander (2013), the DCF presents some inherent problems.First, the DCF adds interpolated points between those from observational data, assuming that the LC varies smoothly, which may be a risky assumption.Another problem of the DCF is the bias in the estimation of the correlation coefficient of each bin, which can produce inconsistencies in the time lag of correlation.To overcome these disadvantages, the z-transformed discrete correlation function was implemented (z-DCF, Alexander 2013).z-DCF avoids the interpolation and performs a correct normalization for the correlation coefficient, ensuring that the absolute value of them is ≤1.This technique estimates the correlation function for sparse, unevenly sampled LCs.To perform the correlation, we identify contemporaneous high-flux states across all wavelengths using the Bayesian block algorithm (Scargle et al. 2013).For this work, we compute the cross-correlations of all the MWL data sets with the corresponding Fermi-LAT LC for each source.As for the periodicity search methods, the z-DCF was also developed for handling unevenly sampled time series with multiple gaps (see Alexander 2013).
No 28-day binning in the MWL LCs (as we do for the periodic analysis) is used for the correlation.The local significance is obtained by applying the same methodology as the previous subsection: simulating 20,000 LCs with the same properties as the original.
The (auto-or cross-) correlation can also be used to search for periodicities.In such a scenario, the correlation function is expected to have a sinusoidal shape with recurrent and approximately equidistant maxima and minima, separated by a delay corresponding to the signal period.It is also a solid method for periodic, non-sinusoidal series, as it is based on the similarity of the LCs at different time lags, regardless of the shape of the data series.To measure the period, we smooth the correlation function with a Savitzky-Golay filter 14 , which reduces the low-frequency variability without distorting the signal tendency (Press & Teukolsky 1990), and identifies the minima and maxima in the resulting output curve.Then, a list of periods can be calculated from the distance between consecutive maxima and minima.The median of the different inferred values is the period of the signal.The uncertainty is obtained by the equation presented by McQuillan et al. (2013) as where  is the number of peaks in the correlation and MAD is the median of the absolute deviations of the periods inferred from the different peaks.For this work, we calculate the periodicity using the cross-correlation between all the MWL data with the corresponding Fermi-LAT LC for each source and the auto-correlation of such MWL data sets.

Variability
To quantify the variability of the studied blazars, we evaluate their fractional variability ( var ) in each wavelength.This parameter was estimated as following the prescription from Vaughan et al. (2003), where its uncertainty can be expressed as where  2 is the variance of the LC, ⟨⟩ is the mean value of the flux in the LC, and ⟨ 2  ⟩ the mean square error. var is affected by the time coverage, sampling, and binning of the data (see Vaughan et al. 2003;Schleicher et al. 2019).Consequently, we employ a similar time window and binning to analyze MWL data.A detailed discussion on the estimation of the  var and caveats such as gaps or uneven sampling can be found in Schleicher et al. (2019).
Additionally, we evaluate the timescales of the observed variability using the structure function (SF).The SF estimates the differences between the squared magnitude as a function of the time separation between measurements as where () and ( + ) are measurements separated by a time interval .Characteristic variability timescales are reflected as flattenings and dips in the SF (see, e.g., Raiteri et al. 2021).Additionally, the 14 Making use of the function "savgol filter" of the Python package "Scipy"  SF can complement the periodicity analysis, as it is expected that periodic variability signatures also appear in the SF as dips at a time interval compatible with the value of the period (Wang et al. 2017).While this method cannot be interpreted as an independent test of periodic variations, it should report the presence of that variability signature, indicating that the reported periodic variability timescale is a real timescale present in the data rather than a red-noise artifact.However, Emmanoulopoulos et al. (2010) presents the problems of the SF, suggesting that the method could provide incorrect timescales depending on the length of the data set and the shape of the associated PSD.Despite this caveat, this tool can complement the periodicity analysis, reveal other variability timescales, and give some information about the origin of the variability (see Section 5).

Periodicity
The periodicity analysis is performed on a total of 9 sources with enough temporal coverage on the MWL LCs to detect periods on the order of a few years (more than 2 cycles using the derived period from P22).The results of Table 2 also reveal compatible periods with those reported in P20 and P22 in 5 of them: PG 1246+586, S4 1144+40, S4 0814+42, PKS 0301−243, and TXS 1902+556.The results of TXS 0518+211 show a slightly longer (∼4 years) period than that reported in P20 and P22 (3.1 ± 0.4 years).We note that the MWL LCs of this source have a shorter coverage w.r.t. the Fermi-LAT curve.Nevertheless, due to the large uncertainties derived by the periodicity analysis for the different inferred values of the period, the two estimates are still compatible.
The period of 2.8 yr 15 obtained for PKS 0447−439 is inconsistent with the 1.9-year period from P20 and P22.Similarly, the period of 1.8 yr for MG1 J021114+1051 is inconsistent with the 2.9-year period from P20 and P22.Finally, a period of 2.6 yr is reported for PKS 0208−512, which is compatible with P20 (2.7 yr) but inconsistent with P22 (3.8 yr).
The local significance of the results typically ranges from 2 to 5, depending on the method.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.The periodicity-search methods respond differently to the gaps of the LCs (see P20).In fact, as stated in Section 3, the LSP, GLSP, WWZ and z-DCF have been specifically developed for having an improved performance under these conditions (see e.g.Lomb 1976;Scargle 1982;Zechmeister & Kürster 2009;Foster 1996;Alexander 2013 for these methods).As demonstrated in P20 and Otero-Santos et al. ( 2023), LSP and WWZ are the most robust against the missing data, and DFT is the most affected method.
The impact of red noise can be seen in Table 2. Specifically, each method has also a different bias to red noise.The most robust ones are again the GLSP and WWZ, as seen in performance tests developed by P22 and Otero-Santos et al. ( 2023), and the most sensitive is again the DFT.This red noise can also result in false peaks and small shifts (see Otero-Santos et al. 2023).For this reason, some methods tend to have incompatible periodicities and local significance (e.g., MG1 J021114+1051 in the radio band).Additionally, when no periodicity is obtained in the original LC, the methods tend to report periods >4 years due to red noise.This is because red noise tends to overproduce periodograms with peaks at long periods (>4 years), which can distort the estimation of this long-period significance.
To be statistically correct, we consider the look-elsewhere effect (Gross & Vitells 2010).We infer the global significance by applying the trial factor to the local significance of each periodicity.This correction is approximated by where  is the trial factor.The trial factor results from the combination of the number of independent frequencies we search for periodicity in this work and the number of blazars in our sample (351, see P20).In P22, we estimated the number of independent frequencies according to the characteristics of the LCs analyzed (e.g., binning, telescope time) and period range of detection ([1-6] years).
To perform this estimation, we consider analyzed LCs that have the same temporal coverage (≈12 years).We also use the same binning for the LCs (28 days).The result is 35 independent frequencies (see P22).Therefore, the trial factor of 351×35.Consequently, the global significance of the periodicity analyses is: (i) ≈2.8 for local significance of ≈5 (13% of the results in Table 2 and Table 3

Correlation
The results of the -ray-MWL cross-correlation analysis are presented in Table 3.Additionally, if the cross-correlation obtained with uncertainties to make the text more readable.However, all reported numbers can be found in their respective tables, with the associated uncertainties.We refer the reader to those for details.
the z-DCF shows a hint/evidence of periodicity, the inferred period is also displayed in this table (i.e., the sinusoidal shape of the derived correlation curve with recurrent and approximately equidistant maxima/minima, see McQuillan et al. 2013).In our notation, positive lags mean that the -ray LC precedes the optical/IR/radio emission.We find that all optical/IR LCs are positively correlated with the -ray emission with time lags compatible with 0 days.Here, all lags < ±28 days are compatible with 0 lag due to the 28-day binning of the Fermi-LAT.We also note that since the optical/IR LCs are expressed in magnitudes, negative values of the correlation function represent a positive correlation between the -ray and optical/IR LCs.This is due to the fact that decreasing magnitudes are translated into increasing fluxes.These correlations show significance ranging from 2 to 5 (pre-trials significance).Our results in the optical and IR bands are compatible with those included in Liodakis et al. (2018): PG 1246+586, S4 0814+42, TXS 0518+211, TXS 1902+556, MG1 J021114+1051, and TXS 0059+581.PKS 0208−512 has a time lag between the -ray and optical bands compatible with 0 days and a local significance of >4.0, in agreement with Chatterjee et al. (2013).
Moreover, radio emission has typically been detected with a delay of a few hundred days (see, e.g.Liodakis et al. 2018).These authors report time lags for the radio--ray correlation of ≈160-240 days for TXS 0059+581, S4 0814+42, TXS 0518+211, and TXS 1902+556, with local significance <3.We find compatible results for S4 0814+42, with a time lag of 207.2±28.1 days.However, we do not have radio data for TXS 0059+581 and TXS 0518+211.We also see a positive correlation for S4 1144+40 with a lag of 51.5±17.3days.Furthermore, no clear correlation is found between the radio and -ray emission for PG 1246+586 and MG1 J021114+1051.For the latter, a hint of correlation (3.0 of local significance) appears with a time lag of −146.5±31.9days, contrary to the typical behavior between these bands.
Especially interesting is the behavior of TXS 1902+556.The crosscorrelations between the optical/radio LCs with the -ray emission display a clear sinusoidal behavior expected from a periodic signal (see McQuillan et al. 2013), with a period of 3.6 ± 0.2 years and a local significance between 2-3 (see Figure 3), which is compatible with the results of Table 2.This is also compatible with the results from P22 and supports the hint of periodic behavior for this blazar.The correlation analysis also reveals that the periodic radio emission is delayed half a period w.r.t. the optical and -ray modulations.Through the correlation analysis between the -ray and MWL data sets, other blazars in our sample also show hints of periodicities in their optical LCs, with the recurrent peaks of the sinusoidal z-DCF functions reaching local significances of 2-3.We find compatible periods w.r.t.those from P22 for PKS 0426-380 (V band), PKS 0447-439, PKS 0250-225, S4 1144+40, S4 0814+42, TXS 0518+211, PKS 2052-47 (V and R bands) and 87GB 164812.2+524023.The blazars S4 0814+42, TXS 0518+211 also have compatible periods with Table 2. Regarding PKS 0208−512, the cross and the autocorrelations show periods compatible of ≈2.7 yr, which again is compatible with P20 but inconsistent with P22.We also observe similar periodicity in X-rays in the autocorrelation, with a lag ≈0 days.The periods are reported in Table 3.Unfortunately, not all the MWL LCs have enough temporal coverage to display a clear sine-like shape in their cross-correlation.MG1 J021114+1051 has lags compatible with 0 days in the optical band, but no periodicity is detected with the previous methods (see Table 2).This result is due to the low sampling of the optical LCs.

Fractional Variability
The results of the variability analysis are shown in Table 3.For 4 out of 5 sources observed by OVRO, we find that the  var is smaller or similar to that in optical/IR wavelengths.Lower  var in the radio band has been the common behavior observed in the past in variability studies.PG 1246+586 is the only source of our sample with a significantly higher  var in radio, with  radio var = 0.71 ± 0.02, compared to  V, R var = 0.27 ± 0.01 found in the V and R optical bands.This indicates that the emission is more variable in radio than in optical wavelengths for this blazar.However, the coverage and sampling of the radio data set are higher than that in the latter wave bands.By computing the  var on the quasi-simultaneous radio and optical data, we obtain again the typically observed behavior for blazars, with a lower  var = 0.14 ± 0.02 in radio.Moreover, the derived values of the  var in the different optical/IR bands are typically compatible with errors.Significant discrepancies (e.g., B band of PKS 2255-282 or B band of PKS 0301-243) are explained as large differences in time coverage of the LCs as for the case of PG 1246+586 already mentioned.In these cases, the  var of the simultaneous data sets are found to show compatible values.
We also evaluate the differences in the  var between BL Lac objects and flat spectrum radio quasars (FSRQs).We find that, while BL Lac objects show rather low  var values (typically between 0.2 to 0.4), FSRQs display a much higher variability ( var ∼ 0.45-0.85),as can be seen in Figure 4, which shows the distribution of the optical  var with respect to that in the -ray regime for the different blazar types.
The estimation of the  var is highly affected by the time coverage, sampling, and binning of the data (see Vaughan et al. 2003;Schleicher et al. 2019).To test the effect of the different sampling in each wave band, we have also estimated the  var with all the simultaneous MWL and -ray data with a matching 28-day binning.This estimation is also represented in Figure 4 with open markers.We observe the same trend as before, where FSRQs tend to be more variable than BL Lac objects.Also, except for a few objects, the optical  var derived for the simultaneous data do not change significantly, pointing towards a dominant long-term variability over shorter timescales.
Table 3. List of periods and the uncertainties (top) with their associated local significance (bottom) for the periodic-emission candidates according to the organization presented in Table 1.The MCMC sine fitting values only show the uncertainties.The -ray period is obtained from the average period (in years) and uncertainty resulting in P22 (as the average of the periodicity analysis).All periods are in years.Some sources have two significant periods (organized by the peak amplitude), denoted by ★.UV1, UV2, UV3 correspond to UVOT 1928 Å, UVOT 2246 Å and UVOT 2600 Å, respectively. corresponds to the data from the Steward Observatory and 15 GHz to the radio band of OVRO. var is the fractional variability, and  SF is the main timescale inferred from the SF.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 peak/valley of the oscillations (see Figure 3).We report one significance when the significance levels of the peak and the valley are the same. (1) (  Finally, as shown in Figure 4, blazars with better sampled MWL LCs and more simultaneous data (e.g., TXS 0518+211) exhibit optical  var values that are consistent between the complete data set and the binned data set.For these blazars, this indicates that long-time scales dominate the variability.For sources with a difference between the  var before and after the binning of the data, there may be significant variability in shorter timescales.Nevertheless, these sources typically correspond to those with the poorest sampling.Therefore, more data are needed to evaluate the presence of significant, faster variability.

Characteristics Timescales
Characteristic variability timescales are also evaluated from the SF.We find timescales compatible within uncertainties with the derived periods by P22 and in this work (see Table 3).However, some discrepant values between the suggested timescales by the SF and the derived periods are observed (for instance, TXS 0059+581, where the derived times scale corresponds to half of the period reported by P22).
Additionally, for some sources, we find more characteristic variability timescales.These timescales show faster variability on the order of a few tens or hundred days.In Table 3, we only show those compatible with the reported periodicities.An example of two SFs in the optical and radio bands is shown in Figure 5.

Power Spectrum Index and Flux Distribution
We estimate the PSDs of each target and the best-fit PL functions that describe these PSDs with the PSRESP package.The resulting indices of the PLs are shown in Table 3.We find that the PL indices range from 1.12 to 1.78.No significant differences in the PSD slope are detected between the BL Lac (with mean =1.37±0.27)and FSRQ (with mean =1.43±0.39)subsamples.Moreover, we see that the mean radio PSDs have slightly steeper slope values (with mean =1.55±0.32)than the optical/IR slopes (with mean =1.40±0.35).However, the sample of blazars with radio data is smaller than that of blazars with optical/IR data (4 vs. 12), which can distort this comparison.
In addition to that, we also evaluate both the flux distributions used for the artificial LC simulation with a normal and a log-normal PDF.We observe that the log-normal distribution provides a better fit to the data.

Periodicity
The periodicity analysis performed on the MWL LCs finds hints of periods compatible with those reported for the -ray LCs in P20 and P22 in some of the analyzed blazars: S4 1144+40, S4 0814+42, TXS 0518+211, TXS 1902+556, PG 1246+586 and PKS 0301−243 (only in the last two blazars the pre-trials median significance is ≳3.5, ≈0 post-trials).This significance is obtained as the median of the results of Table 2 for each method.Note that this median significance does not have an actual statistical meaning; rather, it is used as an arbitrary way of combining all of the significance for a given source to sort the candidates and compare with the results obtained in P22 for the -rays.
This periodicity has been interpreted in the past within several models.The most common interpretations can be summarized in models based on the existence of a binary SMBH system or models based on geometrical effects.
emissions across the electromagnetic spectrum, as seen in OJ 287 and PG 1553+113 (Agudo et al. 2011;Ackermann et al. 2015, respectively).High-resolution simulations confirm that such systems can have yearly orbital periods, influencing jet luminosity through various factors like binary eccentricity, SMBH mass, mass ratio, and the binary's evolutionary stage (Westernacher-Schneider et al. 2022;Zrake et al. 2021).These factors produce a consistent periodic pattern across different wavebands.Another hypothesis involves a helical jet caused by one black hole affecting the other, as suggested by studies on AO 0235+164 (see Ostorero et al. 2004;Raiteri et al. 2006), where the emission's periodic modulation is due to changes in the jet's orientation and the resulting variations in relativistic boosting.
Alternatively, geometrical models are also used to explain MWL periodicity, attributing it to jets influenced by strong magnetic fields.These models include precession-jet (e.g., Villata & Raiteri 1999), rotation-jet (e.g., Hardee & Rosen 1999), and helical structures (e.g., Ostorero et al. 2004), with helical jet models effectively explaining periodicities in blazars like PKS 1830−211 and PKS 2247−131 (Nair et al. 2005;Zhou et al. 2018, respectively).Another possibility sug-gests a twisted jet with periodic orientation changes, impacting the Doppler factor and viewing angle, similar to the effects seen in SMBH binary systems.For example, OJ 287's periodicity has been linked to variations in the Doppler factor due to jet helicity (Butuzova & Pushkarev 2020).

Correlation
The results of the correlation analysis between the -ray and optical/IR bands reveal a high degree of correlation with time lags compatible with <28 days.This behavior implies a co-spatial origin of both emissions, typically expected from leptonic models (Liodakis et al. 2018) rather than from hadronic models (Böttcher 2007).
In leptonic models, the optical emission is produced by synchrotron radiation while the -rays are typically produced by inverse Compton (IC) scattering of either synchrotron photons (Synchrotron Self-Compton, SSC; e.g.Maraschi et al. 1992;Van den Berg et al. 2019;Rajput et al. 2021) or radiation originating from outside of the jet from either the broad-line region or the dusty torus (External Compton, EC; e.g.Sikora et al. 1994;Pacciani et al. 2014).In BL Lacs, particularly HBLs, the process is expected to be SSC, whereas in FSRQs, it is expected to be EC.These strong correlations between the optical and -ray emissions have been observed in the past for other blazars, favoring the leptonic scenario since the same population of electrons is responsible for both emissions (see, e.g.Liao et al. 2014).On the other hand, in hadronic scenarios, the high-and low-energy contributions are caused by different particle populations that are not necessarily easily related.Specifically, protons (and possibly higher-Z nuclei) are responsible for the high-energy emission, and electrons for the low-energy emission (Cerruti 2020).Cohen et al. (2014) claim that FSRQs typically have positive lags, i.e., -rays leading the optical since they are dominated by the EC, while no evident prevailing lag is observed in BL Lacs, with lags ranging between −40 days and 40 days, approximately.The results reported in Table 3 find time lags in our sample of blazars consistent with such short time delays.The optical emission is typically highly correlated with the -ray emission with no delay (<28 days) for all the sources analyzed here.The co-spatial origin can also be applied to the X-rays for PKS 0208−512 since the -ray and the X-ray emissions are correlated without lag with a local significance of >4.
On the contrary, radio emission is typically delayed a few hundred days (see, e.g.Liodakis et al. 2018).We find compatible results with those from Liodakis et al. (2018) for 4 of the 5 sources with radio data available.The most common explanation is that radio emission comes from an outer part of the jet but is triggered by the same physical mechanism (Max-Moerbeck et al. 2014).This is due to the high opacity and self-absorption suffered by radio photons in the inner regions.Thus, radio wavelengths are observable from regions of the jet located further away from the central engine.However, other scenarios, such as differences in the cooling times of the electron populations responsible for the radio and optical emissions, have also been proposed (see for instance Bai & Lee 2003).

Fractional Variability
The fractional variability estimation reveals that FSRQs are typically more variable than BL Lacs.Historically, concerning the BL Lac objects, it has been reported that these sources show lower (higher) variability when the synchrotron peak is located at higher (lower) frequencies (e.g.Otero-Santos et al. 2022).The results found here are in line with this trend observed between the fractional variability and the frequency of the synchrotron peak reported in several works (see, e.g., Figure 26 from Ackermann et al. 2011, where the variability is quantified for all the sources included in the second Fermi-LAT catalog for AGNs, the 2LAC).Here, we observe that ISPs and HSPs show lower  var than FSRQs.Despite the higher variability of the latter, the low number of each type of BL Lacs (especially HSPs) does not allow us to see a clear trend of the variability of the different subtypes with the frequency of the synchrotron peak.Nevertheless, the smaller sample analyzed here is consistent with the results reported for the complete 2LAC catalog.The existence of this trend has been reported by other works like, for instance, Rajput et al. (2020); Bhatta & Dhital (2020).Ackermann et al. (2011) explain these differences in terms of the different cooling timescales of the electrons between the EC and SSC processes of FSRQs and BL Lacs, respectively.Complementary phenomena such as a higher jet power, more efficient accretion disc, and higher electron energies in FSRQs with respect to BL Lac objects lead to a higher variability for the former (Hovatta et al. 2014).Finally, this combination of phenomena also explains the differences between BL Lacs associated with the frequency of their synchrotron peak (Ackermann et al. 2011).
Furthermore, we see that typically, the  var derived for the MWL data set is lower or of the same order as that of the -ray Fermi-LAT LCs.A deep study on the structure of the fractional variability as a function of the energy can also reveal important information regarding the particle population and/or processes producing the broadband emission (see, e.g.Aleksić et al. 2015b).For instance, Mrk 421 shows a  var with its maximum in the X-ray regime, with a second bump in -rays.In this case, this has been interpreted as a higher variability of the high-energy electron population responsible for the high-energy part of the synchrotron emission observed for this source at X-ray wavelengths.
Alternatively, a progressively increasing  var from radio to -rays was observed for Mrk 501 (Aleksić et al. 2015a).This could indicate that for this source, the variability comes from a combination of the low-and high-energy electron populations in the Thomson and Klein-Nishina regimes, respectively.In addition, the  var structure can also change with time, indicating that the processes and particle populations dominating the variability can vary.For example, the double peak structure was also observed for Mrk 501 by Furniss et al. (2015).Therefore, by studying and quantifying the broadband variability of blazars, we can understand the importance of each particle population in the overall variability.However, the lack of ultraviolet and X-ray data for the sources analyzed here does not allow us to draw any reliable conclusion concerning the broadband structure of the  var .More MWL data are needed to extract a firm conclusion in this regard.

Structure Function
The variability analysis performed through the SF reveals the characteristic timescales of the dominant long-term flux variations.For those sources for which the temporal coverage of the data allowed us to perform this analysis, the SF analysis measures variability timescales compatible with the periods previously reported by P20 and P22.Additionally, several shorter timescales are observed in the SFs.The radio is the most stable band with almost no timescales displayed by the SF (see right panel of Figure 5).On the other hand, optical/IR wavelengths are characterized by a variety of variability timescales, ranging from a few tens of days to several hundred days (Figure 5, left panel).
Variability timescales inferred with the SF have been interpreted within several models of different natures.Kawaguchi et al. (1998) explain these timescales based on instabilities in the accretion disc fluctuating in time.Specifically, these timescales would be associated with avalanches of matter on the accretion disc at the same timescales.Three scenarios are presented in Hawkins (2002) to interpret the observed variability.Different timescales are associated with instabilities in the accretion disc in the first models, provoked by the injection of matter with the same timescales.An alternative model associates the variability with supernova events.Finally, in the third scenario, different timescales are induced by microlensing effects caused by MACHOs (dark, compact bodies of planetary sizes), leading to timescales of a few years.The characteristic timescales and slopes of the SFs may indicate which scenario can take place.The average slope of the SFs analyzed in this work is 0.5±0.1,compatible with the 0.44 ± 0.03 slope derived by Hawkins (2002), representative of the accretion disc model.These results are in line with the ones reported by de Vries et al. (2006), where the lensing and starburst models are disfavoured as the origin of the variability in AGNs, supporting the disc instabilities as the cause of the observed variability.This slope also confirms the long-term nature of the variations.Kataoka et al. (2001) claim that short-term variability is modeled in the SF as a PL with index ∼2.0.However, the slope flattens and adopts lower values for longer timescales, as is the case in this analysis.Additionally, the shorter timescales revealed by the SF are explained through models based on instabilities, turbulence, or shocks propagating along the jet (see for instance Marscher & Gear 1985;Marscher 2013).We also compare the mean SF slopes for BL Lacs and FSRQs.No significant differences were found between the slope derived for BL Lacs (0.59±0.10) and the one obtained for the FSRQs (0.50±0.10).
Moreover, we obtain similar values for the optical and radio SF slopes for 3 of the 5 sources with OVRO data: PG 1246+586 (  = 0.61,   = 0.60 and    = 0.63), S4 0814+42 (  = 0.47,   = 0.53 and    = 0.55) and TXS 1902+556 (  = 0.51 and    = 0.56).According to the results presented by Collier & Peterson (2001), this could be evidence of a common mechanism producing both emissions.On the other hand, MG1 J021114+1051 shows a higher slope for the radio SF (   = 1.26) w.r.t. the V-and R-band SFs (  = 0.42 and   = 0.41).This difference between optical and radio can explain the different periods for these wavebands shown in Table 2 (1.8 and 4 yrs, respectively).The low temporal coverage of the V-band LC of S4 1144+40 does not allow us to compare its slope with the one obtained from the radio data set.

PSD characterization
We also evaluate the PSDs of the different LCs.These PSDs are fitted to a PL function to study the type of underlying noise associated with their variability.The range of spectral indices obtained here is placed between the flickering (pink) and red noise regimes, as expected from blazars with rapid flares and short timescale fluctuations, and long-term (periodic) variability (see Abdo et al. 2010).A pure pink noise signal is described as a PL with a dependency of 1/  with the frequency (Hufnagel & Bregman 1992).On the other hand, pure red noise processes are described as 1/  2 and are typically associated with the erratic and stochastic variability displayed by blazars (Vaughan 2005).Periodic and stochastic variability in different timescales can co-exist in these sources, leading to the admixture of pink and red noise nature, with indices ranging between 1.0 and 2.0.This is supported by SFs that reveal variability on different timescales that are not necessarily associated with the periodic modulation of the emission and can be stochastic in nature, which will lead to the observed spectral indices.

Flux Distribution
A log-normal distribution fits the flux distributions analyzed.This result reflects the influence of the accretion disc on the variability.Indeed, fluctuations in the accretion disc at different radii are propagated inwards to produce an aggregate multiplicative effect in the innermost disc.This disturbance is then transmitted to the jet, generating the log-normal distribution of the variability observed in the MWL emission (Shah et al. 2018).The usual normal distribution scenario is associated with additive processes, where independent events are linearly added (see, e.g.Biteau & Giebels 2012 and references therein).However, a log-normal distribution does not necessarily imply multiplicative processes (Scargle 2020).

SUMMARY
We analyze the -ray and MWL (optical, IR, and radio) data of 19 blazars of P22, showing hints of periodic modulation in 5 of them with ≥3.0 (≈0 post-trials).The ≈0 in the post-trials significance does not allow to claim any presence of periodicity in these sources.However, observing the same period could indicate similar regularoscillating behavior in the MWL emission of such blazars.Therefore, these 5 blazars are promising candidates to be monitored in the next years.
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.This can be the result of the radio being emitted from an outer region of the jet due to the absorption suffered at these wavelengths in the inner regions.Additionally, the variability analysis performed with the SF shows variability timescales compatible with the periods associated with the emission.This result supports the existence of a periodic modulation in the emission of 10 blazars.Shorter timescales revealed by the SF support that such variability could originate from instabilities or fluctuations in the accretion disc.
A combination of pink and red noise for the emission of these objects is confirmed through the PSD analysis.The spectral indices derived from this fit range from 1.12 to 1.78.Finally, the flux distribution analysis shows that a log-normal distribution successfully fits the PDF of the MWL LCs of these objects.These PDFs are characteristic of multiplicative variability processes due to fluctuations in the accretion disc transmitted to the jet.

Figure 2 .
Figure 2. Power-law index distribution as obtained using the PSRESP method on the V-band LC of TXS 1425+516.

Figure 3 .
Figure3.Left: Z-DCF cross-correlation results between the R-band and -ray LCs of TXS 1902+556.Right: Z-DCF cross-correlation results between the radio and -ray LCs of TXS 1902+556.Both panels display hints of periodic correlation (approximated by the black line which smooths the correlation with Savitzky-Golay filter), denoting that both signals could present similar regular-oscillating behavior with a significance of 2sigma (pre-trials).The periods are inferred from the distance between the different maxima/minima of the smoothed curves.
have a local significance of ≈5) (ii) ≈1.8 for a local significance of ≈4.5 (11% of the results in Table 2 and Table 3 have a local significance of ≈4.5) (iii) <1 for local significance <4.5 (76% of the results in Table 2 and Table 3 have a local significance <4.5)