© The Authors 2025

1. Introduction

Active galactic nuclei (AGN) are among the most luminous and distant persistent objects in the Universe. These are supermassive black holes (SMBH) at the center of some galaxies that accrete matter from surrounding disks of gas and dust (e.g., Urry 1996; Combes 2021). These SMBHs can produce highly collimated jets of particles, accelerated to relativistic velocities (e.g., Blandford et al. 2019). When the jet is closely aligned with our line of sight, these objects are referred to as blazars.

Blazars exhibit variability across the entire electromagnetic spectrum, with timescales ranging from seconds to years (e.g., Urry 2011). While this variability is typically considered stochastic, it can also reveal characteristic patterns (e.g., Raiteri et al. 2001; Sandrinelli et al. 2016; Peñil et al. 2020, 2022). These patterns provide insight into the internal structure of blazars and their emission processes (e.g., Bhatta 2019). Periodicity in blazar emissions can result from several physical mechanisms, including (1) jet precession (e.g., Camenzind & Krockenberger 1992), (2) modulations in the accretion flow (e.g., Gracia et al. 2003), or (3) the presence of binary SMBH systems (e.g., Sobacchi et al. 2016; Celoria et al. 2018). Binary SMBHs, which typically result from galaxy mergers and play a fundamental role in galaxy evolution (e.g., Jiang et al. 2022), are a key focus of this work. Theoretical studies exploring periodic variability linked to binary SMBHs (e.g., Sobacchi et al. 2016; Celoria et al. 2018) highlight the relevance of searching for year-scale periodicities to investigate their dynamic evolution in the gas-driven regime and role in galaxy evolution (e.g., Peñil et al. 2024a, b).

In addition to periodic signals, long-term trends can also be observed in blazars (e.g., Marscher et al. 2011; Rani 2018; Valverde et al. 2020). The variability timescales in SMBH systems reflect their physical properties and dynamic processes. For single SMBH systems, variability spans from minutes to years, with short-term changes often caused by instabilities in the inner accretion disk or jet phenomena, while long-term changes may result from variations in the accretion rate or disk precession (e.g., Peterson 2001; McHardy et al. 2006). In binary SMBH systems, additional periodicities from the orbital motion of the black holes are superimposed on the variability from individual accretion disks (e.g., Graham et al. 2015; Sobacchi et al. 2017). During SMBH mergers, tidal interactions and dynamic evolution lead to complex variability patterns, including quasi-periodic oscillations (QPOs) with timescales from weeks to years (e.g., Merritt & Milosavljević 2005; Valtonen et al. 2008). While linking periodic behavior to specific physical scenarios is an essential goal, achieving this will require more detailed analyses, such as multiwavelength modeling and theoretical simulations, which are beyond the scope of this work.

Identifying QPOs and long-term trends poses an observational challenge due to the need for continuous monitoring over extended periods. In the γ-ray domain, the Large Area Telescope (LAT) aboard the Fermi Gamma-ray Space Telescope, launched in June 2008, monitors the entire sky with high sensitivity. Using this data, Peñil et al. (2020) conducted the first systematic periodicity study based on the third release of the Fermi-LAT catalog (3FGL, Acero et al. 2015), analyzing 9 years of LAT data and 10 different techniques. Peñil et al. (2020) identified 11 blazars with periodic signals showing local significance ≥2.0σ.

However, these studies primarily relied on Fourier-domain techniques, which face inherent limitations in addressing nonstationary signals typical in astrophysical observations (e.g., Prokhorov & Moraghan 2017; Zhang et al. 2017; Sandrinelli et al. 2018; Gupta et al. 2019; Ait Benkhali et al. 2020; Peñil et al. 2020, 2022; Wang et al. 2022; Ren et al. 2023; Li et al. 2023). The complexity of blazar emissions, including noise structures (e.g., Vaughan et al. 2016), long-term trends, and flares, can be better addressed using time-domain methods that can effectively isolate periodic components from the overall signal (see also, e.g., Rueda et al. 2022).

Singular spectrum analysis (SSA) emerges as a promising solution (e.g., Hassani 2007; Nina Golyandina 2020). This statistical method, designed for time series data, decomposes a signal into its trend, oscillatory components, and noise structure. By separating these components, SSA provides a clearer picture of periodic patterns and helps in developing more accurate predictive models. SSA has been successfully applied across various scientific domains, including economics (e.g., Hassani et al. 2009), climate science (e.g., Dabbakuti 2019), and biomedicine (e.g., Sanei & Hassani 2015).

In this study, we performed a systematic search for periodic emissions in a sample of AGN, including blazars, radio galaxies, and narrow-line Seyfert 1 galaxies (see Sect. 2), detected by Fermi-LAT in the second release of the fourth Fermi-LAT source catalog (4FGL-DR2, Abdollahi 2020; Ballet et al. 2020).

The paper is organized as follows: Sect. 2 presents the blazar sample and data reduction methodology. Sect. 3 describes the SSA method, including its application to periodicity, trends, and forecasts. The results are discussed in Sect. 4, and a summary and conclusion are provided in Sect. 5.

2. Data reduction and sample selection

In this section, we describe the γ-ray data reduction and sample selection.

2.1. Analysis of Gamma-ray Data

Fermi-LAT is a γ-ray telescope and pair-conversion detector, capable of measuring energies from below 20 MeV to over 300 GeV (Atwood et al. 2009). Thanks to its wide field of view (2.4 sr), it continuously scans the entire sky with the possibility of observing the entire sky every 3 hours. Data collected by LAT are made publicly available through the Fermi Science Support Center¹ (FSSC). The Fermi-LAT data used in this study are processed with Fermipy, a Python package, following the procedure outlined by Wood et al. (2017). We specifically use the Pass 8 SOURCE class, from the P8R3_SOURCE_V2 instrument response function (IRF), to select the photons (Atwood et al. 2013). Our region of interest (RoI) is a 15°×15° square, centered at the target. To mitigate contamination from γ-rays produced in the Earth's upper atmosphere, we apply a zenith angle cut of θ<90°. Standard quality cuts are applied (DATA QUAL >  0), (LAT CONFIG == 1), and time periods coinciding with solar flares and γ-ray bursts are excluded. The standard LAT analysis of a RoI includes sources located within 20° from the RoI center, using the 4FGL catalog (Abdollahi 2020), which includes the Galactic and isotropic diffuse emissions (gll_iem_v07.fits and iso_P8R3_SOURCE_V6_v06.txt).

We conduct a binned analysis, using 10 energy bins per decade from 0.1 to 800 GeV with 0.1° spatial bins, which produces a maximum likelihood analysis over the entire observation period (2008 Aug 04 15:43:36 UTC to 2020 Dec 10 00:01:26 UTC, approximately 12 years). Using the parameters and spectral shapes from the 4FGL catalog, we model the sources within the RoI. The normalization and spectral index of the target source, as well as sources within a 3° radius centered at the target position, are allowed to vary, along with the two diffuse backgrounds. We construct a test statistic (TS) map for searching for newly detected sources relative to the 4FGL catalog. The TS is defined as 2log(L/L₀) (Mattox et al. 1996), where L is the likelihood of the model with a point source at a given position and L₀ is the likelihood without the point source. This statistic quantifies how significantly a source emerges from the background (Abdollahi 2020). A power-law spectrum is used to model new test sources, where only the normalization is allowed to vary during the fitting process, while the photon index is fixed at 2. We look for significant peaks, retaining only sources with TS >25 and separated by at least 0.5° from existing sources in the model. We then add a new point source at the position of the most significant peak. The model for the RoI is readjusted, and a new TS map is produced. This iterative process continues until no further significant excesses are found, generally leading to the addition of two point sources. Any γ-ray excess not listed in the DR2 catalog but with TS >  25 in the full-time interval analysis is included in the RoI as part of our analysis strategy.

To generate light curves (LCs), we split the data for each source into 28-day bins and perform a full likelihood fit for each bin. The best-fit RoI model from the time-interval analysis is adopted. We first focus on the target source and all nearby sources within the RoI (i.e., within 3°), including the diffuse components. During the fitting process, only the normalization parameters are allowed to vary. If the fitting does not converge, we gradually reduce the number of free parameters until convergence is achieved.

An iterative approach is used, starting by fixing sources with low detection significance (TS <  4) within the RoI. Next, we address sources with TS <  9, and then fix sources within a 1° radius of the RoI center if their TS <  25. Finally, we leave only the normalization of the target source free, fixing all other parameters. Note these LCs are the same as those used by Peñil et al. (2022, 2024a, b).

2.2. Sample selection

Our sample is drawn from the 4FGL-DR2 catalog, which includes data from approximately 3500 blazars (Ajello et al. 2021). We filtered the 4FGL-DR2 sample by selecting only sources that show variability in their emissions. Specifically, we used a variability index threshold, requiring the variability index to be greater than 18.48 (Abdollahi 2020). This threshold indicates that the probability of the source remaining steady is less than 1% (see Table 12 of Abdollahi 2020). We used 28-day binned LCs to be compatible with our data reduction, which is optimized for detecting periodicities on the order of 1 year, as this provides a practical balance between computational efficiency and sensitivity to long-term variations.

In cases of non-detection, we established an upper limit for the flux in all time bins where the TS is below 4, making our sample more conservative. The choice of using upper limits for <50% of the data was tested by Peñil et al. (2020), who found that LCs with more upper limits were too incomplete for reliable analysis. We followed the same approach as Peñil et al. (2022) for handling upper limits. To address potential inaccuracies from weak signals or non-detections (TS <  1), we replaced flux values with those that maximize the likelihood function for each time bin. This approach was applied to about 8% of the data. For the sources listed in Table A.1, the average percentage of time bins with upper limits (TS <  4) is approximately 16%.

3. Methodology

3.1. Theoretical background on singular spectrum analysis

SSA is a statistical technique for time series analysis that integrates elements of classical time series analysis, dynamical systems, multivariate statistics, multivariate geometry, and signal processing (e.g., Hassani 2007; Nina Golyandina 2020). The core idea of SSA is to decompose the original signal into a sum of components, such as a smooth trend, an oscillatory component, and a noise structure.

The SSA technique consists of two main stages: decomposition and reconstruction. Both stages include two separate steps. In the first stage, the series is decomposed to enable signal extraction and noise reduction. In the second stage, the original series is reconstructed, and the less noisy reconstructed series is utilized for forecasting new data points. Below, we provide a concise explanation of the SSA methodology.

1. Embedding. The starting point of SSA is to embed a one-dimensional time series Y_N=(y₁, …, y_N) of length N into a multidimensional series X₁, …, X_K with vectors $X_i=\left(y_i,\dots ,y_{i+L-1}\right)\in {ℝ}^L$, where K=N−L+1. The single parameter for embedding is the window length L, an integer such that 2≤L≤N/2. The result of this step is the trajectory matrix X = [X₁,…,X_K]:

   $$\text{X}=\left(\begin{array}{ccccc}\hfill y_1\hfill & \hfill y_2\hfill & \hfill y_3\hfill & \hfill \cdots \hfill & \hfill y_K\hfill \\ \hfill y_2\hfill & \hfill y_3\hfill & \hfill y_4\hfill & \hfill \cdots \hfill & \hfill y_{K+1}\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill y_L\hfill & \hfill y_{L+1}\hfill & \hfill y_{L+2}\hfill & \hfill \cdots \hfill & \hfill y_N\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \end{array}\right),$$(1)

   where the trajectory matrix X is a Hankel matrix, meaning that all elements along the diagonal i+j=const are equal. This step breaks down the LC into overlapping time windows, similar to a sliding window approach, to analyze trends and oscillations in smaller sections.

2. Decomposition. The trajectory matrix X is then decomposed using singular value decomposition (SVD):

   $$\text{X}=\text{U}\mathrm{\Sigma }{\text{V}}^T,$$(2)

   where $\mathrm{\Sigma }=\text{diag}(\sqrt{{\lambda }_1}\ge \sqrt{{\lambda }_2}\ge \dots \ge \sqrt{{\lambda }_L})\in {ℝ}^{L\times L}$, $\text{V}=[V_1,V_2,\dots ,V_L]\in {ℝ}^{K\times L}$ and $\text{U}=[U_1,U_2,\dots ,U_L]\in {ℝ}^{L\times L}$. Here, λ₁, λ₂, …, λ_L are the eigenvalues of XX^T in decreasing order (λ₁≥λ₂≥…≥λ_L≥0), and U₁, U₂, …, U_L are the corresponding eigenvectors. Set d as the maximum i, such that λ_i>0, representing the rank of X. If we denote $V_i={\text{X}}^T{U}_i/\sqrt{{\lambda }_i}$, the SVD of the trajectory matrix can be written as:

   $$\text{X}={\text{X}}_1+{\text{X}}_2+\dots +{\text{X}}_d,$$(3)

   where ${\text{X}}_i=\sqrt{{\lambda }_i}{U}_i{V}_i^T(i=1,2,\dots ,d)$. The matrices X_i have rank 1. SVD (3) is optimal in the sense that, among all matrices of rank r<d, the matrix ${\text{X}}^{(r)}={\sum }_{i=1}^r{\text{X}}_i$ provides the best approximation to the trajectory matrix X, minimizing ∥X−X^(r)∥_F, where ∥.∥_F is the Frobenius norm. Note that $\parallel \text{X}{\parallel }^2={\sum }_{i=1}^d{\lambda }_i$ and ∥X_i∥²=λ_i for i = 1, …, d. The ratio ${\lambda }_i/{\sum }_{i=1}^d{\lambda }_i$ quantifies the contribution of matrix X_i to the expansion (3). Consequently, ${\sum }_{i=1}^r{\lambda }_i/{\sum }_{i=1}^d{\lambda }_i$ is a measure of the approximation of the trajectory matrix by matrices of rank r. This step separates the LC into its main components, such as long-term trends, oscillatory behaviors, and random noise, allowing us to analyze each aspect separately.

3. Grouping. The grouping step involves splitting the elementary matrices X_i into several groups and adding the matrices within each group. The grouping procedure partitions the set of indices {1, …, L} into m disjoint subsets {I₁, …, I_m} (eigentriple grouping), where each group I_c contains a set of principal components {i₁, …, i_p} representing specific components of the signal (trend, oscillations, harmonic components, etc.). For example, if the periodic component can be obtained using I_period={1, 2}, then ${\text{X}}_1+{\text{X}}_2={\text{X}}_{I_{\mathrm{period}}}$. For more details, see Section 2.2.2 in Sanei & Hassani (2015). Let I_c={i₁, …, i_p} be a group of indices. Then the matrix ${\text{X}}_{I_c}$ corresponding to group I_c is defined as ${\text{X}}_{I_c}={\text{X}}_{i_1}+\dots +{\text{X}}_{i_p}$. These matrices are computed for I_c=I₁, …, I_m, and the SVD expansion leads to the decomposition $\text{X}={\text{X}}_{I_1}+\dots +{\text{X}}_{I_m}$. Each ${\text{X}}_{I_c}$ represents a set of eigentriples describing a specific component of the original time series. In short, after separating the different components of the LC, we group those that represent similar behaviors. For example, if two components both show a periodic pattern, we group them together to extract the full periodic signal. This helps us isolate different types of behavior in the data, such as long-term trends or periodic oscillations.

4. Reconstruction. The final step in basic SSA is to determine the components corresponding to the original series. This involves deriving m new time series ${\tilde{Y}}^{(c)}$ of length N, where each time series corresponds to a specific component of the signal (e.g., trend, oscillations). Diagonal averaging of each ${\text{X}}_{I_c}$ provides the c-th component of the series Y_N. The n-th sample is obtained by averaging over the cross-diagonal i+j=n+1=const of ${\text{X}}_{I_c}$, since each ${\text{X}}_{I_c}$ is expected to be a Hankel matrix for the corresponding embedded component series. This procedure is called diagonal averaging or Hankelization of a matrix. Applying Hankelization to all matrix components of ${\text{X}}_{I_1}+\dots +{\text{X}}_{I_m}$ yields another expansion:

   $$\text{X}={\tilde{\text{X}}}_{I_1}+\dots +{\tilde{\text{X}}}_{I_m}.$$(4)

   This is equivalent to decomposing the initial series Y_N=(y₁, …, y_N) into a sum of m series:

   $$y_t={\displaystyle \sum _{c=1}^m}{\tilde{y}}_t^{(c)}$$(5)

   where ${\tilde{Y}}_N^{(c)}=({\tilde{y}}_1^{(c)},\dots ,{\tilde{y}}_N^{(c)})$ corresponds to the matrix $X_{I_c}$.

   In summary, after decomposition and grouping, we combine the components we are interested in, such as the periodic signal, while filtering out noise.

3.2. Further notes on SSA: Separability and grouping selection

As mentioned earlier, SSA decomposes a time series into various components, which can be categorized as trend, oscillatory components, and noise. The window length L is the only parameter in the decomposition stage. If the time series contains a periodic component with an integer period, it is recommended to set the window length proportional to that period to enhance the separability of the periodic component. A periodic component generates two eigentriples with closely matched singular values, making it beneficial to check the breaks obtained from the eigenvalue spectra. Additionally, a pure noise series typically produces a slowly decreasing sequence of singular values.

To assess the quality of the decomposition of the observed data into signal and noise components, we determine separability. To measure this separability, we use the w-correlation matrix, which is useful for finding correlations and identifying groups. This matrix consists of weighted cosines of angles between the reconstructed time series components. Specifically, we analyze the matrix of absolute values of the w-correlations. If the absolute value of the w-correlations is small, then the corresponding series (here, signal and noise) are almost w-orthogonal. However, if the w-correlation is large, then the two series are far from being w-orthogonal and are thus poorly separable. If two reconstructed components have zero w-correlations, they are separable. Large w-correlations between reconstructed components suggest they should be grouped together and correspond to the same component in SSA decomposition.

Let us consider a time series Y_N, where $Y_N=Y_N^{(1)}+Y_N^{(2)}$, with $Y_N^{(1)}$ representing the main signal (trend + periodic components), and $Y_N^{(2)}$ representing the noise. If $Y_N^{(1)}$ can be separated from Y_N, it means there exists a partition into groups at the grouping stage such that ${\tilde{Y}}_N^{(m)}=Y_N^{(m)}$.

We apply the separability formalism defined in Section 2.1.3 of Golyandina et al. (2018), summarized as follows. Let X^(m) be the trajectory matrix of the series, with their SVDs, ${\text{X}}^{(m)}={\sum }_{i=1}^{d_m}\sqrt{{\lambda }_{m,i}}{U}_{m,i}{V}_{m,i}^{\text{T}}$. The column and row spaces of the trajectory matrix are called column and row trajectory spaces, respectively. The separability conditions are:

1. Two series $Y_N^{(1)}$ and $Y_N^{(2)}$ are weakly separable if their column and row trajectory spaces are orthogonal.

2. Two series are strongly separable if they are weakly separable and the sets of singular values of their L-trajectory matrices are disjoint, λ_{1, i}≠λ_{2, j}, for any i and j.

Let w_n (n = 1, …, N) be the number of times the series element x_n appears in the trajectory matrix. Consider two series of length N, $Y_N^{(1)}$ and $Y_N^{(2)}$. Set $(Y_N^{(1)},Y_N^{(2)}{)}_w={\sum }_{n=1}^N{w}_n{y}_n^{(1)}{y}_n^{(2)}$, the w-scalar product of the time series. The w-correlation matrix is defined as:

$${\rho }_w(Y_N^{(1)},Y_N^{(2)})=\frac{\left(Y_N^{(1)},Y_N^{(2)}\right)}{(\parallel Y_N^{(1)}{\parallel }_w\parallel Y_N^{(2)}{\parallel }_w)}·$$(6)

The well-separated components in Equation (5) are weakly correlated, while poorly separated components have high correlation. By examining the matrix of w-correlations between elementary reconstructed series ${\tilde{Y}}_N^{(i)}$ and ${\tilde{Y}}_N^{(j)}$, it is possible to identify groups of correlated series components. It is important to avoid grouping highly correlated components separately. For more details, see Hassani (2007).

An example of the w-correlation matrices is shown in Figure 1. The matrix plots different components of the decomposed LC along both the x-axis and y-axis. Each point $\left({\tilde{Y}}_N^{(i)},{\tilde{Y}}_N^{(j)}\right)$ represents the degree of correlation between the i-th and j-th component. The firsts components are usually the trend and periodicity, whereas the rest of the components can be characterized as noise. The color represents the level of correlation, where yellow indicates a high correlation (close to 1) and dark purple shows near-zero correlation. A correlation of 1 means the components are almost identical, while a correlation close to 0 suggests the components are well-separated (orthogonal). If the components are well-separated (low correlation), we can confidently treat them as distinct parts of the signal, such as separating periodic behavior from noise. The diagonal elements (always 1) show the self-correlation of each component. Typically, a value greater than 0.8 means that the components are significantly correlated, and we consider them from the same group, that is, either trend, periodicity, or noise. For instance, in this example plotted in Figure 1, the components 1 $\left({\tilde{Y}}_N^{(1)}\right)$ and 2 $\left({\tilde{Y}}_N^{(2)}\right)$ are considered correlated and they both store the periodicity information, while the component 0 $\left({\tilde{Y}}_N^{(0)}\right)$ is characterized as trend. The dashed red lines separate the trend and periodicity components from the noise. Figure 1 (right) shows a zoom of the first ten components.

Fig. 1. w-correlation matrices from SSA decomposition of S5 1044+71 LC. Dashed red lines separate the top-left region, containing the trend and periodic components, from the rest of the matrix, where noise components are more dominant. (Left panel) General w-correlation matrix, where each point $\left({\tilde{Y}}_N^{(i)},{\tilde{Y}}_N^{(j)}\right)$ represents the correlation between the i-th and j-th component. The color scale indicates the strength of the correlation, with yellow representing high correlation and dark purple representing low correlation. The diagonal represents self-correlation. (Right panel) Zoomed-in view focusing on the first 10 components, providing a detailed examination of the correlations between the initial components that typically capture key trends and periodic signals.

3.3. Challenges and limitations in the search for periodicity

Several limitations are inherent in any periodicity analysis, including the presence of different noise types. Noise can be characterized by an exponent β, which represents the power-law index of the power spectral density (PSD). The PSD is crucial for characterizing variability emissions and quantifies the variability power as a function of the frequency ν, with PSD(ν)∝ν^−β, where β = 0 for white noise (spurious background independent of frequency), β = 1 for pink noise, and β = 2 for red or Brownian noise (Rieger 2019). These noise types can distort the signal and make the detection of periodic patterns more challenging (VanderPlas 2018), potentially leading to false detections (Vaughan et al. 2016).

Underlying trends and/or unpredictable stochastic emissions, such as flares, can also distort the periodicity analysis. SSA effectively mitigates these challenges by decomposing the signal, allowing us to search for inherent periodic patterns by isolating the oscillatory component from the underlying trend and noise structure.

In Section 3.2, we discussed the significance of the window length as a key parameter in our analysis. There are no strict guidelines for selecting a specific window length in SSA, as long as it falls within the range 2≤L≤N/2, where N is the data length. However, larger window lengths (up to 30–45% of the length of the time series) are sometimes necessary to adequately separate periodicities from the overall trend. For more details on window length selection, see Nina Golyandina (2020).

We find that a window length of 40% of the data length is optimal for our analysis. When the window length is lower, stochastic events tend to dominate the signal, whereas when the window length is larger, there are too many components for adequate grouping.

We search for long-term periodicities in the range of 1–6 years, as our 12-year LCs allow for the detection of at most two cycles within this timeframe. This range is particularly relevant for binary SMBH candidates in the gas-driven regime (see, e.g., Peñil et al. 2024a), where systems with parsec to sub-parsec separations between the black holes, which cannot generally be spatially resolved, are expected to show periodicities. Focusing on this range also minimizes the impact of red noise contamination, which primarily affects longer periodicities, while shorter periodicities are less relevant to this study.

3.4. Periodicity search pipeline

We applied SSA to the LCs of the blazars selected in our sample, as detailed in Section 2.2.

To calculate the period, we employed a widely recognized technique for identifying periodicity in time series data, the Lomb-Scargle periodogram (LSP, Lomb 1976; Scargle 1982). We applied the LSP to the periodicity component obtained after applying the SSA code to the original signal. The period corresponded to the highest peak in the periodogram.

3.4.1. Systematic search for periodicity candidates

We performed a systematic search for periodicity in our sample to select the most promising blazars. The LSP was computed for simulated LCs to determine their power values. These simulated LCs were generated using the best-fitting parameters derived from the PSD and the Probability Density Function (PDF) of the original LC. We applied the Emmanoulopoulos technique² to generate artificial LCs and evaluated the significance of the periods (Emmanoulopoulos et al. 2013). The mean power value within each period bin was multiplied by various percentiles corresponding to confidence significance levels. The significance of the detected period was assessed by comparing the periodicity peak derived from the LSP for the periodicity component of the studied blazar with the computed significance levels. As an example, Figure 2 shows the significance levels for PG 1553+113. We selected all blazars where the significance of the peak in the LSP is ≥2.0σ.

Fig. 2. Lomb-Scargle periodogram of PG 1553+113. A significant period peak is detected at ∼2.2 yr. The different dashed lines represent the significance levels computed by simulating artificial LCs (see Section 3.4.1).

3.4.2. Local significance

For this subsample, we refined the analysis and calculated a more robust significance. Specifically, we generated simulated LCs² to estimate the significance of the identified periods by computing the local significance, as explained below. We used the same procedure outlined in Appendix A of O’Neill et al. (2022). Note that this procedure already accounts for the uncertainty arising from the fact that the signal frequency is initially unknown, necessitating a systematic search over a range of possible periods.

First, we calculate the period significance, defined as the probability of identifying a periodic signal with the same period, p_peak, and equal or greater power, P_peak, at the LSP. This involves evaluating how often the observed period of the AGN is replicated in the simulated LC. Consequently, we derive a p-value, denoted as p-value_peak. The local significance is obtained by applying the same procedure with the periods obtained in the different simulated LCs.

In summary, the steps of the significance procedure are as follows:

1. For each simulated LC, we apply the SSA decomposition and obtain the periodicity component.

2. We compute the LSP for each periodicity component, obtaining the period, p_sim, and the power of the period peak, P_sim.

3. For each simulated LC, the period significance is calculated, p-value_sim. We count all simulations where the period p_sim is detected with power ≥P_sim.

4. To estimate the local significance, we compute the probability that p-value_sim≤p-value_peak.

3.4.3. Global Significance

In periodicity analyses, we must correct the local significance for the look-elsewhere effect (e.g., Ait Benkhali et al. 2020). This occurs in situations where a signal is searched for only in a specific region of parameter space (Vitells 2011). In our case, we are searching for periodicity signals within a sample of sources. The look-elsewhere effect describes the ratio between the probability of observing the excess at a fixed value and the probability of observing it anywhere in the range. It is quantified in terms of a specific number of trials (Gross & Vitells 2010; Peñil et al. 2022) as follows,

$$p_{\text{global}}=1-(1-p_{\text{local}}{)}^N,$$(7)

where p_local is the p-value associated with the local significance, and N is the trial factor, computed from the total number of blazars in our sample (N = 494).

Since our local significance was limited by the total number of simulations, we performed a large number of artificial LC simulations to obtain a meaningful significance. One million simulations were conducted for each AGN in our sample, with the highest significance obtained being 4.8σ, 1 event as significant as the one observed in 1 million simulated LCs assuming no real periodic signal. After applying the trial factor corrections, we obtained the global significances listed in Table 1.

3.5. Evaluation of the method

In this section, we perform tests to evaluate the performance of SSA under different conditions.

3.5.1. Tests against noise

Our evaluation involves a purely periodic signal, resembling a sinusoidal waveform with a predefined period, using the same time sampling as the original LCs. To simulate realistic conditions, we introduced random noise to the signal, specifically white, pink, and red noise, following Timmer & Koenig (1995). This process was repeated multiple times to observe how often SSA successfully detects the predetermined period in the presence of noise.

In this study, we used the Timmer & Koenig (1995) method via the colorednoise Python package to simulate red, pink, and white noise components. This approach evaluates the robustness of SSA in detecting periodic signals under realistic conditions for blazar emissions.

We note that this methodology is consistent with Monte Carlo SSA techniques (Allen & Smith 1996), which combine SSA with Monte Carlo simulations using AR(1) noise models (autoregressive processes where each value depends on the previous one) to assess the statistical significance of detected components.

The periodic signal is computed as:

$$\varphi (t)=Asin\left(\frac{2\pi t}{T}+\theta \right)+O,$$(8)

where A is the amplitude, T the period, θ the phase, and O the offset of the signal. These parameters are randomly selected within a range of values. The ranges were computed using the Markov Chain Monte Carlo sinusoidal fitting method from the blazars in Table A.1, as explained in Section 3.1.6 of Peñil et al. (2020). The period range is selected from the periods found in this work,

$$\begin{array}{cc}\hfill A& =[1--20]\times {10}^{-6}\mathrm{MeV}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1},\hfill \\ \hfill O& =[0-150]\times {10}^{-6}\mathrm{MeV}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1},\hfill \\ \hfill T& =[1.0-4.5]\text{yr},\hfill \\ \hfill \theta & =[0,2\pi ]\text{rad}.\hfill \end{array}$$

We generated the noise structure using the Python package colorednoise³, based on Timmer & Koenig (1995), and considered red, pink, and white noise.

We simulated 100 000 LCs by adding random noise with the same data length and amplitude to the sinusoidal time series for each period. The detection rate depends on the period, as shown in Figure 3.

Fig. 3. Detection rates as a function of period. Rates are computed using simulated LCs with added noise (see Section 3.5.1).

As expected, the detection rate decreases as the signal period increases. This is a logical outcome since longer periods result in fewer cycles within the analyzed time interval. The methods utilized in Peñil et al. (2022) report a detection range spanning approximately 12–65%, depending on the period range under analysis. For periods in the [1–2.5] year range, SSA exhibits a detection rate of 78%, improving the methods in Peñil et al. (2022) by 18%. For longer periods in the [3.5–4.5] year range, SSA shows a 50% enhancement in detection compared to the methods in Peñil et al. (2022). These findings highlight the advantages of SSA in detecting periods in signals contaminated by stochastic noise.

3.5.2. False positive rate

We also assessed whether the algorithm used to obtain the period component produces genuine results. For this, we generated pure red, pink, and white noise LCs (using Emmanoulopoulos et al. 2013) with the same time sampling as the original LCs. We applied our algorithm to determine how often we obtained significant results. This allows us to quantify the probability of detecting a period with high significance that has occurred stochastically. We performed 500 000 simulations and find that 0.13% have a significance ≥3σ, and 2.27% have a significance ≥2σ. A false positive rate of 2.27% implies that, out of 494 stochastic LCs, you would expect 13 to show periodic signals at 2.0σ.

3.6. Long-term trends

Identifying periodic patterns in time series data can be challenging when there is a prominent trend. Applying the LSP to analyze signals with both periodicity and trend components often amplifies the periodogram's peaks. This complicates the detection of underlying periodicity and increases uncertainty in the derived period. It is recommended to detrend the signal as an initial step in periodicity analysis (e.g., Welsh 1999), since trends can impact the low-frequency components of the PSD. Such trends may falsely suggest periodic patterns within the time series, leading to erroneous detection of false periodicity (McQuillan et al. 2013). Therefore, separating the periodic and trend components is essential for accurate periodicity estimates. Additionally, understanding the general condition of the blazar in which periodicity is detected is crucial, as it can significantly impact emissions and the observed periodicity.

In some cases, the trend can dominate the periodic structure, such as when the periodicity and amplitude of the signal increase or decrease as the trend rises or falls. Since the first components obtained through SSA decomposition represent the trend of the time series, we analyze these components to characterize the trend.

It is common to apply statistical time series tests such as the F-test to assess trends. The null hypothesis of the F-test assumes that the trend can be represented as a constant. Applying the F-test to our analysis indicated that all blazars exhibited a trend. However, note that this test assesses the entire time series signal, while we focus solely on the trend reported by SSA. The F-test concludes that the trend is not constant because SSA does not report a perfectly consistent line; instead, it reveals small-scale variations around the overall blazar emission state. Since we are interested in long-term trends, the F-test is unsuitable for this analysis, as it consistently indicates a trend due to small-scale fluctuations. To address this, we propose an ad-hoc conservative criterion to select blazars with trends. A linear regression of the trend component was performed to find the best-fit slope. We considered the trend as increasing when the slope was >0.02, decreasing when the slope was <−0.02, and constant when the slope was in the range [−0.02, 0.02].

In Section 4.2, we characterize the trends of our periodicity candidates.

3.7. Forecasts

Predicting future states of γ-ray blazars is important, especially for planning observational proposals. In the context of periodicity, understanding and characterizing underlying trends and periodic patterns enables the development of predictive models for future source behavior. To achieve this, we applied SSA, specifically recurrent SSA forecasting (SSA-R⁴, Kalantari et al. 2020).

Our method forecasts under the assumption that the time series satisfies a Linear Recurrent Relation (LRR). The time series (y₁, …, y_N) satisfies an LRR of order d if there are coefficients a₁, …, a_d such that:

$$y_{i+d}={\displaystyle \sum _{k=1}^d}{a}_k{y}_{i+d-k},1\le i\le N-d,a_d\ne 0,d<N.$$(9)

The coefficients (a₁, …, a_d) are the LRR coefficients.

Let I be the set of eigentriples determined at the grouping step of SSA, as described in Section 3.1. The first r eigentriples are selected to obtain forecasts of the reconstructed series. In our case, we considered only the group containing the trend and the periodicity component. If i∈I, and the eigenvectors of the chosen eigentriples are U_i∈R^L, denoted by $\underset{\_}{U_i}\in {\text{R}}^{L-1}$ (the vector consisting of the first L−1 components of U_i, where L is the window length referenced earlier), and π_i the last component of U_i, then $v^2={\sum }_{i\in I}{\pi }_i^2$. Let ℒ⊂R^L be the linear space spanned by the U_i vectors, where ℒ=span(U_i, i∈I). Assuming ℒ is a non-vertical space (i.e., e_L∉ℒ, with e_L=(0, 0, …, 1)^T∈R^L), then v²<1.

We define the vector R=(a_L−1, …, a₁)^T as:

$$R=\frac{1}{1-v^2}{\displaystyle \sum _{i\in I}}{\pi }_i\underset{\_}{U_i},$$(10)

and the SSA-R forecasting algorithm can be summarized as follows:

The i-th element of the time series Y_N+h={y₁, …, y_N+h} is defined as:

$$y_i=\{\begin{array}{cc}\widetilde{y_i}& \text{for}\hspace{0.5em}i=1,\dots ,N\\ \sum _{j=1}^{L-1}{a}_j{y}_{i-j}& \text{for}\hspace{0.5em}i=N+1,\dots ,N+h,\end{array}$$(11)

where $\widetilde{y_i}$ (i = 1, …, N) are the reconstructed series. The values y_N+1, …, y_N+h are the h-step-ahead recurrent forecasts (for further details, see Section 3 in Golyandina et al. 2018).

The confidence intervals are obtained by considering the residuals, that is, the components not used for forecasting (in particular, noise). The procedure is as follows.

Set the Y_N time series, which can be described as a sum of components $Y_N^{(1)}+Y_N^{(2)}$, where $Y_N^{(1)}$ is governed by an LRR of relatively small dimension, and the residual series $Y_N^{(2)}$, $Y_N^{(2)}=Y_N-Y_N^{(1)}$, is approximately strongly separable from $Y_N^{(1)}$ at some window length L.

To establish confidence bounds for the forecast, we assume that Y⁽²⁾ is a finite random noise series perturbing the signal Y⁽¹⁾.

Consider a bootstrap method to construct the confidence bounds for the signal $Y_N^{(1)}$ with rank r<L. By applying SSA and reconstructing the signal, we obtain $Y_N={\tilde{Y}}_N^{(1)}+{\tilde{Y}}_N^{(2)}$. In our case, ${\tilde{Y}}_N^{(1)}$ is the reconstructed component related to the trend and periodicity of the studied source. The next step is to calculate the empirical distribution of the residuals ${\tilde{Y}}_N^{(2)}$ and simulate Q independent copies ${\tilde{Y}}_{N,q}^{(2)}$ of the $Y_N^{(2)}$ series. We apply the forecasting procedure to Q independent time series defined as ${\tilde{Y}}_{N,q}={\tilde{Y}}_N^{(1)}+{\tilde{Y}}_{N,q}^{(2)}(q=1,\dots ,Q)$. With an M step-ahead forecast, we obtain the elements ${\tilde{y}}_{N+M,q}^{(1)}(1\le q\le Q)$. Afterward, we calculate the (empirical) lower and upper quantiles at some level (in our case, 95%) and obtain the corresponding confidence interval for the forecast, known as bootstrap confidence intervals.

3.8. Software

We used Python to perform our analysis, documented in a Jupyter Notebook that is made public online⁵ for clarity and reproducibility. The underlying SSA analysis relies on commonly used libraries such as numpy for matrix operations (e.g., SVD using numpy.linalg.svd).

Established Python libraries such as pyrssa⁶ and pyts.decomposition⁷ can be used to reproduce the results presented in this paper. A tutorial available on Kaggle⁸ offers a helpful introduction to SSA.

To compute LSP, the scipy.signal.lombscargle function was utilized. For forecasting, SSA-R forecasting techniques were employed to predict the future behavior of the sources. Python packages such as pyrssa, MSSA⁹, and the implementation available on GitHub¹⁰ provide robust tools for reproducing forecasting results.

4. Results and discussion

In this section, we provide a comprehensive discussion of the results derived in this work, listing periodicity candidates and conducting an extensive literature search to compare our findings with prior research. In Section 3.5, we evaluated SSA's detection capabilities in the presence of random noise and assessed the likelihood of obtaining high-significance candidates that may correspond to purely stochastic emission sources.

4.1. Periodicity candidates

From our systematic periodicity search pipeline, we identified 46 AGN with evidence of periodic γ-ray emissions, including 24 flat spectrum radio quasars (FSRQs), 18 BL Lacertae objects (BL Lacs), 3 blazars of uncertain type (BCUs), and 1 Narrow-Line Seyfert 1 (NLSy1) galaxy. Table A.1 lists these candidates, along with their periods and local significance calculated with 10⁶ simulations.

Figure 4 shows the results obtained using SSA for the blazar PG 1553+113. This blazar has gained significant attention due to the detection of a 2.2-year period with notable significance (Ackermann et al. 2015; Prokhorov & Moraghan 2017; Peñil et al. 2022). Our analysis yielded a 2.2±0.2 yr period, which is consistent with previous studies, supported by a local significance of ∼4σ (global significance ∼2.2σ), as reported in Table A.1. Detailed plots corresponding to the candidates presented in this paper can be found in Appendix B.

Fig. 4. SSA applied to PG 1553+113. Blue dots represent the original signal with error bars. Orange, green, and red lines correspond to the periodicity component, the trend, and the sum of the trend and the periodicity reported by the SSA algorithm, respectively. The spacing between the blue bars indicates the period computed through the LSP (2.2 yr). The width of the blue bars is the period uncertainty (0.2 yr) computed from the FWHM of the period peak in the LSP.

Figure 5 shows the relationship between the photon index and integrated photon flux of all the blazars detected by Fermi-LAT and our periodicity candidates. Our periodicity candidates are typically detected for bright sources with integrated flux ≳10⁻⁹ph cm⁻² s⁻¹.

Fig. 5. Photon spectral index as a function of integrated photon flux. We show all the blazars detected by Fermi-LAT (gray circles) and blazars selected for periodicity analysis with a variability index above 18.48 (dark gray circles). Our periodicity candidates are shown according to their AGN type: BL Lacs (blue circles), FSRQs (red triangles), NLSy1 galaxies (orange square), and BCUs (green crosses).

Furthermore, Figure 6 shows the distribution of redshifts for our candidates according to the detected periodicities. This distribution aligns with expectations, as BL Lac objects tend to be found at lower redshifts and FSRQs at higher redshifts (Ajello et al. 2020). Candidates with longer periods tend to exhibit greater uncertainty, as red noise has more impact at lower frequencies, leading to broader peaks in the LSP and fewer cycles in the data. Figure 7 shows the sky map of the periodicity candidates.

Fig. 6. Period in years as a function of redshift, with sources classified according to AGN type. PG 1553+113 and other candidates discussed in Sections 4.1.1 and 4.1.2 are labeled.

Fig. 7. Sky map of periodicity candidates classified by AGN type.

We categorized the selected sources into two groups: those for which previous studies have reported periodicity in γ-ray emissions and those that are newly identified periodicity candidates presented in this work.

4.1.1. New γ-ray periodicity candidates

Twenty-five sources are reported as new γ-ray periodicity candidates in our analysis with a local significance ≳2σ. Of these, 14 are classified as FSRQs, 8 as BL Lac-type blazars, and 3 as BCUs. BCUs have a broad-band SED resembling that of blazars but lack spectroscopic confirmation. Below are examples of new candidates, with their periods and local significance. Note that a local significance of 4.8σ is the maximum significance achievable given the simulations we conducted.

• OC 457. We find a 1.8±0.2 yr (4.8σ) period. The LC is shown in Appendix B. This blazar is a case where the trend decreases, causing the periodicity to disappear toward the end. Consecutive upper limits and low emission also impact the analysis, but this is the first time OC 457 has been identified as a periodicity candidate.

• PKS 0524−485. A 2.1±0.2 yr (4.8σ) period was detected. Unlike OC 457, this blazar shows a clear periodic pattern over time, with a consistent and steady increase in the long-term trend, suggesting a strong connection between the trend and the observed periodicity.

• TXS 1530−131. A 1.4±0.1 yr (4.8σ) period was detected. The trend is highly consistent, with the periodic cycles maintaining both their frequency and amplitude, indicating regularity in their behavior.

• PKS 2155−83. A 4.7±0.9 yr (4.8σ) period was found. This is the longest-period candidate, but the uncertainty is higher due to red noise contamination at lower frequencies, leading to a broader LSP peak.

• 3C 66A. Previous studies report quasi-periodic behavior in the optical band but not in γ-rays, with a period of ∼2.3±0.3 yr in the V-band (Cheng et al. 2022; Otero-Santos et al. 2020). We find a 2.3±0.2 yr (∼3.8σ) period in γ-rays, indicating a connection between the optical and γ-ray emission mechanisms.

4.1.2. Periodicity candidates in literature

Twenty-one of the 46 blazars with significant periodic emission in this work have a period reported in the literature across multiple bands. Table A.2 lists these blazars, ordered by highest to lowest local significance according to our results. Many of the results align with those in previous works. Below, we discuss some specific cases.

• S3 0458−02. Peñil et al. (2022) find a ∼3.8±0.9 yr period, while we find a ∼1.00±0.04 yr (4.8σ) period. Although SSA can detect both periods, the significance of the 3.8-year period is lower than 2σ. The 1-year period is more distinct in the LC, demonstrating SSA's ability to uncover underlying periodicities more effectively.

• B2 1520+31. We detect a 3.4±0.6 yr (4.8σ) period. The literature reports a periodicity of a few days (Gupta et al. 2019; Ren et al. 2023), but no long-term period has been reported.

• S5 0716+71. Peñil et al. (2022) find a ∼2.7±0.4 yr period, compatible with our result of 2.9±0.3 yr (∼4.5σ). Ren et al. (2023) and Peñil et al. (2022) also report a secondary period of ∼0.9 yr, which we also detect with SSA.

• PKS 0208−512. Peñil et al. (2020) reported a ∼2.7±0.1 yr period, while Peñil et al. (2022) reported a ∼3.8±0.5 yr period. We detect a 2.5±0.3 yr (∼4.5σ) period. The discrepancy arises from a flare in 2020 that distorts the periodicity, leading to a detected 4-year period. SSA more accurately detects the underlying ∼2.5 yr period.

• PKS 1510−089. Li et al. (2023) find a period of ∼3.6±0.2 yr, while we detect a period of 1.7±0.1 yr (∼2.5σ). Although no other studies report long-term periodicity, there have been studies reporting transient periodicities of a few days (Roy et al. 2022).

Figure 8 presents a comparison between the periods detected using SSA and those reported in the literature, as summarized in Table A.2.

Fig. 8. Periodicity results using SSA compared with those reported in the literature, as summarized in Table A.2. Blue bars represent SSA results from this work, while orange and green bars represent periodicity values reported in the literature. Orange bars indicate periodicity values from individual studies, while green bars are added when other studies report additional periodicity values for the same blazar. Error bars indicate the uncertainties associated with each period.

4.2. Long-term trends

We conducted a comprehensive trend analysis, following the methods outlined in Section 3.6, for all the periodicity candidates in this study. The results are presented in Table A.3, showing that 21 candidates exhibit a constant long-term trend, 13 show increasing trends, and 12 show decreasing trends.

An example is the NLSy1 PKS 0250−225, where the trend is decreasing, and the periodicity disappears over time, following the trend (see LC in Appendix B). The opposite case is the blazar 4FGL MH 2136−428, where the periodicity becomes more evident and the amplitude increases as the trend increases. SSA is particularly effective in such cases for detecting significant periods.

4.3. Forecast

For all candidates, we used the algorithm described in Section 3.7 to predict the ∼10-year upcoming gamma-ray emission peaks from 2021. Table A.3 details the next four emission cycles. The uncertainty in peak prediction was based on the period uncertainty calculated with the LSP. In some cases, there were no values for the second peak prediction because the next emission occurs beyond the time interval used for forecasting. For the NLSy1 galaxy PMN J0948+0022, predictions were not possible because the amplitude trends toward zero, leading to the disappearance of periodicity.

In cases where periodicity disappears, the computed confidence interval becomes a rectangle, where the periodic time structure fades, and successive cycles cannot be predicted reliably. To determine significant forecasts and calculate the dates of future emission peaks, we applied a methodology that fits the upper bound of the confidence interval to a linear regression. We considered candidates with an R² parameter of less than 50%, a moderate value according to Hair et al. (2011). Forecasting is inadequate for blazars where periodicity was detected in more advanced SSA components, such as those associated with flares (e.g., PKS 0208-512). Our forecasting models are designed to work with AGN both with and without long-term trends. When a source lacks a long-term trend, the model adapts to the data without introducing an artificial trend. As a result, the absence of long-term trends does not negatively impact the overall results, as the models are flexible and reflect the true characteristics of each source.

We successfully predicted the next four γ-ray emission cycles for 28 candidates. Figure 9 shows the forecasting model for PG 1553+113. We used the 28-day time-bin data from 2008 to 2021 to predict the next two cycles. Appendix C contains forecasting plots for all candidates, except those where the methodology was inadequate.

Fig. 9. SSA-R forecasting models for PG 1553+113. (Top panel) SSA-R forecasting model for PG 1553+113 using 28-day time-bin data. (Bottom panel) SSA-R forecasting model for PG 1553+113 using 56-day time-bin data until 2023.

As a proof of concept, we evaluated the model's efficacy by forecasting the 56-day time-bin data of PG 1553+113 up to early 2023. We initiated the forecast around the end of 2020 and predicted the next two cycles. The 56-day time-bin data is shown as blue points in Figure 9. The results demonstrate a good fit between the forecasted and observed data, as can be visually observed in Figure 9 (bottom).

When comparing our predicted dates with updated data from the Fermi-LAT Light Curve Repository¹¹, we observe that the next cycle occurred in February 2021, compared to our initial prediction of March 2021. The uncertainties arise from the period uncertainty. The subsequent cycle occurred on May 3, 2023, while our prediction was May 8, 2023. These results for PG 1553+113 show a good agreement with Fermi-LAT data, considering the uncertainties in cycle period determination.

5. Summary and conclusions

In this study, we employed SSA to conduct a systematic search for periodic patterns in AGN within the time domain. This marks the first application of this powerful statistical technique in this field. Our analysis includes 494 gamma-ray sources with well-sampled LCs. SSA effectively isolates underlying oscillatory components within the signal, detecting periodicities that may remain undetected by conventional methods. We employed a robust methodology to determine the significance of the detected periodicities by generating artificial LCs using Monte Carlo methods. SSA also facilitates long-term trend characterization and the development of forecasting models capable of predicting future emission peaks.

The main results of our study are as follows:

1. We identified 46 γ-ray periodicity candidates with a local significance of ≥2σ, including 21 known candidates and 25 newly discovered in this work, effectively doubling the number of known candidates.

2. Long-term trends were characterized for 25 of these AGN.

3. Forecasting models were provided for 28 of the AGN.

While this study successfully identifies γ-ray periodicities in a sample of blazars, the task of linking these periodicities to specific physical mechanisms remains unresolved. The periodicities observed on year-long timescales could plausibly be attributed to processes such as jet precession (e.g., Camenzind & Krockenberger 1992; Britzen et al. 2019), accretion flow instabilities such as warped disks or quasi-periodic oscillations (e.g., Gracia et al. 2003; Mishra et al. 2022), or SMBBH orbital motion (e.g., Graham et al. 2015; Sobacchi et al. 2017). However, distinguishing between these scenarios necessitates detailed multiwavelength observational campaigns to investigate correlations across different parts of the electromagnetic spectrum (e.g., Peñil et al. 2024a), the study of the time-dependence of the QPOs to examine their stability or evolution (e.g., Abdollahi et al. 2024), and theoretical simulations (e.g., Westernacher-Schneider et al. 2022). Such analyses would enable a more comprehensive understanding of the physical origins of the observed periodic behavior but fall outside the scope of this work. This study provides a foundation for future investigations by identifying promising candidates and highlighting their potential significance within the context of blazar variability.

Acknowledgments

The authors thank Sarah Wagner and the anonymous referee for their valuable comments, which improved the manuscript. A.R. acknowledges the support of an Investigo grant funded by the European Union, Next Generation EU. A.D. is thankful for the support of Proyecto PID2021-126536OA-I00 funded by MCIN / AEI / 10.13039/501100011033. A.R. and M.A. are grateful for the support of the Fermi Guest Investigator Cycle 17 program (proposal number 171063). The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration 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 was performed in part under DOE Contract DE-AC02-76SF00515.

References

Appendix A: Tables

Appendix B: Singular spectrum analysis plots

This section presents the LCs of all our candidates. The blue dots represent the γ-ray data from Fermi-LAT. Black arrows indicate the upper limits, which are considered when the flux value has a test statistic lower than 1. The upper-limit flux value is the flux that maximizes the likelihood function for the corresponding time-bin. The orange, green, and red solid lines represent the periodicity, trend, and the sum of the trend and periodicity from SSA, respectively. The distance between the blue rectangles represents the source period, and their width indicates the period uncertainty.

Fig. B.1. Light curves of quasi-periodic candidates presented in Table A.1.

Fig. B.1. (Continued).

Fig. B.1. (Continued).

Fig. B.1. (Continued).

Fig. B.1. (Continued).

Appendix C: Forecasting plots

This section presents the forecasting plots for our periodicity candidates. The blue dots represent the γ-ray data from Fermi-LAT. Black arrows indicate the upper limits. The orange and green solid lines represent the periodicity + trend and its forecast, respectively. The gray area shows the 95% confidence interval of the forecast.

Fig. C.1. Forecasting LCs of quasi-periodic candidates. Predicted cycles are reported in Table A.3.

Fig. C.1. (Continued).

Fig. C.1. (Continued).

Fig. C.1. (Continued).

Fig. C.1. (Continued).

All Tables

All Figures

Fig. 1. w-correlation matrices from SSA decomposition of S5 1044+71 LC. Dashed red lines separate the top-left region, containing the trend and periodic components, from the rest of the matrix, where noise components are more dominant. (Left panel) General w-correlation matrix, where each point $\left({\tilde{Y}}_N^{(i)},{\tilde{Y}}_N^{(j)}\right)$ represents the correlation between the i-th and j-th component. The color scale indicates the strength of the correlation, with yellow representing high correlation and dark purple representing low correlation. The diagonal represents self-correlation. (Right panel) Zoomed-in view focusing on the first 10 components, providing a detailed examination of the correlations between the initial components that typically capture key trends and periodic signals.

In the text

	Fig. 2. Lomb-Scargle periodogram of PG 1553+113. A significant period peak is detected at ∼2.2 yr. The different dashed lines represent the significance levels computed by simulating artificial LCs (see Section 3.4.1).
In the text

	Fig. 3. Detection rates as a function of period. Rates are computed using simulated LCs with added noise (see Section 3.5.1).
In the text

Fig. 4. SSA applied to PG 1553+113. Blue dots represent the original signal with error bars. Orange, green, and red lines correspond to the periodicity component, the trend, and the sum of the trend and the periodicity reported by the SSA algorithm, respectively. The spacing between the blue bars indicates the period computed through the LSP (2.2 yr). The width of the blue bars is the period uncertainty (0.2 yr) computed from the FWHM of the period peak in the LSP.

In the text

Fig. 5. Photon spectral index as a function of integrated photon flux. We show all the blazars detected by Fermi-LAT (gray circles) and blazars selected for periodicity analysis with a variability index above 18.48 (dark gray circles). Our periodicity candidates are shown according to their AGN type: BL Lacs (blue circles), FSRQs (red triangles), NLSy1 galaxies (orange square), and BCUs (green crosses).

In the text

	Fig. 6. Period in years as a function of redshift, with sources classified according to AGN type. PG 1553+113 and other candidates discussed in Sections 4.1.1 and 4.1.2 are labeled.
In the text

	Fig. 7. Sky map of periodicity candidates classified by AGN type.
In the text

Fig. 8. Periodicity results using SSA compared with those reported in the literature, as summarized in Table A.2. Blue bars represent SSA results from this work, while orange and green bars represent periodicity values reported in the literature. Orange bars indicate periodicity values from individual studies, while green bars are added when other studies report additional periodicity values for the same blazar. Error bars indicate the uncertainties associated with each period.

In the text

	Fig. 9. SSA-R forecasting models for PG 1553+113. (Top panel) SSA-R forecasting model for PG 1553+113 using 28-day time-bin data. (Bottom panel) SSA-R forecasting model for PG 1553+113 using 56-day time-bin data until 2023.
In the text

	Fig. B.1. Light curves of quasi-periodic candidates presented in Table A.1.
In the text

	Fig. B.1. (Continued).
In the text

	Fig. B.1. (Continued).
In the text

	Fig. B.1. (Continued).
In the text

	Fig. B.1. (Continued).
In the text

	Fig. C.1. Forecasting LCs of quasi-periodic candidates. Predicted cycles are reported in Table A.3.
In the text

	Fig. C.1. (Continued).
In the text

	Fig. C.1. (Continued).
In the text

	Fig. C.1. (Continued).
In the text

	Fig. C.1. (Continued).
In the text