Recurrence plot analysis of blazar gamma-ray light curves: Exploiting the time-domain capabilities of Fermi-LAT

Variability studies of jetted AGN, in particular blazars, have been used to gain a better understanding of the particle acceleration mechanisms in jets. However, statistical methods used for the characterization of variability often rely on stationary time series data, which is not fulfilled for most blazar light curves. We introduce the recurrence plot method for long-term $γ$-ray light curves sampled by Fermi-LAT and present our results for the BL Lac object Mkn 421 and the FSRQ PKS 1424-41. Using surrogates to determine the significance of our findings, we conclude that Mkn 421 exhibits more determinism than PKS 1424-41, and that both sources potentially show nonlinearity. However, the latter has to be tested against more advanced surrogates that are able to replicate the nonstationarity of the original light curves. In future work, we will extend our recurrence analysis to a sample of $\sim50$ $γ$-ray bright sources to probe the jet dynamics of different blazar classes.

💡 Research Summary

The paper introduces recurrence plot (RP) analysis as a novel, non‑linear time‑series tool for studying the γ‑ray variability of blazars, addressing the limitation of traditional methods that assume stationarity. Blazars, with jets aligned close to our line of sight, display rapid and large‑amplitude flux changes that encode information about particle acceleration, jet dynamics, and central engine processes. However, most variability analyses—fractional variability, periodograms, power spectral densities (PSDs), wavelet transforms—require the mean and variance of the light curve to remain constant over time, an assumption violated by the flaring, non‑stationary nature of blazar emission.

The authors apply RP analysis to long‑term, high‑cadence Fermi‑LAT γ‑ray light curves of two bright sources: the BL Lac object Markarian 421 (Mkn 421) and the flat‑spectrum radio quasar PKS 1424‑41. Light curves are binned in 7‑day intervals, and only bins with a test statistic (TS) ≥ 1 are retained (discarding 1.9 % of Mkn 421 data and 3.9 % of PKS 1424‑41 data). The RP construction follows the classic time‑delay embedding method (Sauer et al. 1991). For Mkn 421 the embedding delay is 8 bins with an embedding dimension of 4; for PKS 1424‑41 the delay is 21 bins and the dimension is 6. The recurrence matrix is built by marking a ‘1’ whenever two points in the reconstructed phase space lie within a distance ε, and ‘0’ otherwise. The resulting square matrices are visualised as recurrence plots.

Two quantitative RP descriptors are examined: Determinism (DET) and Lmax. DET is the fraction of recurrence points that form diagonal lines of length ≥ 2, serving as a proxy for predictability and deterministic structure. Lmax is the length of the longest diagonal line, related to the largest positive Lyapunov exponent, and is often used as an indicator of non‑linearity or chaos.

To assess statistical significance, the authors generate ensembles of Iterative Amplitude Adjusted Fourier Transform (IAAFT) surrogate light curves (100 realizations per source). IAAFT surrogates preserve the original power spectrum and flux distribution but destroy any deterministic dynamics, thereby representing the null hypothesis of a stationary linear Gaussian process. By comparing DET and Lmax of the real data against the surrogate distributions across a range of recurrence rates (0.05–0.20, corresponding to different ε thresholds), the authors avoid bias from selecting a single ε value.

Results show that Mkn 421 exhibits a robust DET excess over its surrogates across a broad range of recurrence rates, indicating a higher degree of deterministic, repeatable behavior. PKS 1424‑41 only shows a modest DET excess at the lowest recurrence rates, suggesting weaker determinism and stronger non‑stationarity. Both sources display Lmax values significantly larger than those of their surrogates, hinting at underlying non‑linear dynamics. However, the authors caution that Lmax can be inflated by non‑stationarity (Kugiumtzis 2001), and note that standard IAAFT surrogates do not preserve non‑stationary features (Borgnat & Flandrin 2009). Consequently, they propose using more sophisticated surrogates such as Pinned Wavelet IAAFT (PWIAAFT) in future work to disentangle true non‑linearity from mere non‑stationarity.

The paper concludes that RP analysis, together with DET and Lmax, provides a valuable complementary perspective to traditional frequency‑domain methods for blazar variability studies. Mkn 421 appears more deterministic than PKS 1424‑41, while both may harbor non‑linear processes that merit deeper investigation. The authors outline a future program to extend this methodology to a sample of ~40–50 γ‑ray bright blazars, spanning both BL Lac and FSRQ classes, with the goal of probing jet dynamics, searching for quasi‑periodic oscillations on week‑to‑year timescales, and statistically comparing RP metrics across blazar subclasses.

Overall, the study demonstrates that recurrence plot techniques can capture dynamical information hidden in non‑stationary γ‑ray light curves, opening a new avenue for exploring the complex physics of relativistic jets in active galactic nuclei.