EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 707 (March 2026) A&A, 707 (2026) L18 Full HTML
Open Access
Issue
A&A
Volume 707, March 2026
Article Number L18
Number of page(s) 4
Section Letters to the Editor
DOI https://doi.org/10.1051/0004-6361/202659239
Published online 20 March 2026

© The Authors 2026

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1. Introduction

Blazars are radio-loud active galactic nuclei with relativistic jets aligned with the line of sight and that show extreme variability across the electromagnetic spectrum (e.g., Urry & Padovani 1995; Abdo et al. 2010). Although this variability is mainly stochastic, some sources display quasi-periodic oscillations (QPOs) in their light curves, whose physical origin remains widely debated in high-energy astrophysics (e.g., Sandrinelli et al. 2016; Jorstad et al. 2022).

The BL Lacertae object PG 1553+113, classified as a high-synchrotron peaked (HSP) blazar (e.g., Nievas Rosillo et al. 2022; Láinez et al. 2025; Dinesh et al. 2025), represents a benchmark case for such studies. Ackermann et al. (2015) reported a compelling QPO in its gamma-ray flux with a period of ≈2.2 years, a result that has been subsequently confirmed and extended by later analyses (Tavani et al. 2018; Peñil et al. 2020, 2024; Abdollahi et al. 2024; Rico et al. 2025). We note that continuous wavelet transforms may yield lower QPO significance due to high trial factors (Ren et al. 2023). The persistence of this oscillation has fueled the hypothesis that PG 1553+113 hosts a supermassive binary black hole (SMBBH) system (e.g., Sobacchi et al. 2017; Agazie et al. 2024), where the periodicity arises from orbital interaction or hydromagnetic instabilities. Alternative scenarios propose geometric origins, such as jet precession (e.g., Rieger 2004; Liska et al. 2018), which could itself be dynamically driven by a binary companion (Begelman et al. 1980).

Distinguishing whether the QPO is driven by intrinsic plasma mechanisms (e.g., periodic accretion changes or shocks) or geometric effects (e.g., a binary-driven precessing jet) requires a robust investigation of the spectral variability. Standard correlations between flux and photon index, such as the softer-when-brighter or harder-when-brighter trends, can trace the underlying particle acceleration processes (e.g., Fossati et al. 1998; Bhatta et al. 2016). However, blazar light curves are dominated by red noise, which can introduce spurious correlations if not properly treated (e.g., Vaughan et al. 2003). Recent works using singular spectrum analysis (SSA) have demonstrated the importance of decomposing blazar light curves to separate intrinsic variability and noise from baseline trends (e.g., Peñil et al. 2025; Rico et al. 2025).

In this letter we constrain the physical driver of the QPO in PG 1553+113 by analyzing 17 years of Fermi-LAT data. We employ a robust statistical framework, combining SSA detrending with block bootstrap resampling (e.g., Künsch 1989), to test for coupling between the photon index and the flux, as well as the QPO phase.

2. Theoretical background

2.1. Emission models, photon index, and variability

The broadband spectral energy distribution (SED) of HSP blazars such as PG 1553+113 is typically described by the one-zone synchrotron self-compton (SSC) model (e.g., Katarzyński et al. 2001; Aleksić et al. 2015). In this framework, the emission originates from a spherical region of plasma within the relativistic jet, characterized by a magnetic field B and moving with a Doppler factor δ. The low-energy peak is produced by synchrotron radiation from ultra-relativistic electrons, while the high-energy gamma-ray peak arises from inverse Compton scattering of these same synchrotron photons by the parent electron population (e.g., Rybicki & Lightman 1979).

The observed gamma-ray photon index, Γ, derived from the power-law approximation dN/dE ∝ E−Γ, serves as a proxy for the underlying electron energy distribution index, p. Under standard radiative cooling conditions, these quantities are related by Γ = (p + 1)/2 (e.g., Abdollahi et al. 2020). Variations in Γ therefore probe the microphysics of the emitting region, reflecting changes in the balance between particle acceleration and radiative cooling rates.

Correlations between the gamma-ray flux (Fγ) and photon index (Γ) provide additional information on the physical drivers of variability (e.g., Abdo et al. 2010). A harder-when-brighter behavior, characterized by a decrease in Γ as the flux increases, is commonly interpreted as evidence for the injection of fresh high-energy electrons through efficient acceleration processes, such as shocks, and is often associated with flaring states (e.g., Sambruna et al. 2000). In contrast, a softer-when-brighter trend, in which Γ increases with flux, may indicate variability driven primarily by changes in particle density or enhanced radiative cooling rather than by spectral hardening (e.g., Abdo et al. 2010).

2.2. Geometric versus plasma-driven modulation

Physical models for QPOs can be broadly categorized by their prediction of chromatic versus achromatic variability. Intrinsic plasma-driven models, such as regular accretion changes induced by a binary companion (Valtonen et al. 2009) or gravitational stresses (Tavani et al. 2018), predict phase-dependent spectral changes in the electron distribution, leading to chromatic signatures like periodic spectral hardening (e.g., Sobacchi et al. 2017). Conversely, geometric models, such as jet precession, attribute the variability to periodic changes in the Doppler factor δ. Since the bolometric flux scales as Fγ ∝ δ4, while the spectral shape remains largely achromatic in the comoving frame (e.g., Rieger 2004; Malik et al. 2025), geometric origins are expected to produce flux modulation with little to no dependence of the photon index on the QPO phase. Finding an achromatic QPO would thus disfavor accretion- or shock-driven periodicities, while remaining fully compatible with a binary-driven precessing jet.

3. Methodology

3.1. Data and source selection

The data for this analysis were obtained from the Fermi Large Area Telescope (LAT), the primary instrument on board the Fermi Gamma-ray Space Telescope (Atwood et al. 2009). We accessed the data through the Fermi-LAT Light Curve Repository1 (LCR; Abdollahi et al. 2023), a public database providing flux and spectral information for sources in the 4FGL catalog (Abdollahi et al. 2020), whose automated products for PG 1553+113 have been validated against targeted manual reductions (e.g., Peñil et al. 2024; Abdollahi et al. 2024; Rico et al. 2025). The dataset for PG 1553+113 spans a 17-year interval, covering the period from August 2008 to September 2025 (approximately MJD 54 600 to MJD 61 000).

To characterize this multi-year variability and ensure sufficient statistics for spectral fitting, we selected a time binning of 30 days. This integration time is significantly shorter than the QPO period, thus allowing a detailed phase-resolved analysis, while effectively smoothing out short-term stochastic fluctuations that could mask the periodic modulation. The analysis covers the standard Fermi-LAT energy range from 100 MeV to 100 GeV, following the standard LCR methodology.

We utilized the Free Index fitting option provided by the LCR, which allows the photon index Γ to vary in each time bin. To ensure data quality, we applied a test statistic (TS) threshold of TS ≥ 4 (approximately equivalent to a 2σ detection) for each bin. Time intervals failing this criterion, which often correspond to high background noise or low source flux states, were excluded to avoid introducing spurious outliers in the correlation analysis. The final dataset consists of simultaneous measurements of Fγ and Γ, along with their respective uncertainties, serving as the input for the SSA and correlation routines described below.

3.2. Correlation analysis and treatment of red noise

We quantified the relationship between flux and photon index using two independent metrics: the Pearson correlation coefficient (r), which evaluates linear dependence, and the Spearman rank correlation coefficient (ρ), which assesses monotonic relationships and is less sensitive to outliers (Walpole et al. 2012).

To account for experimental uncertainties in the Fermi-LAT data, we employed a Monte Carlo bootstrap approach. We generated N = 2000 synthetic datasets by perturbing both flux and photon index values within their Gaussian observational uncertainties. Accordingly, the analysis is anchored to the results obtained from the unperturbed data, which are adopted as the nominal values, while the 1σ and 2σ uncertainties are inferred from the spread of the corresponding bootstrap distributions in order to mitigate attenuation bias associated with the randomization (e.g., Peterson et al. 1998).

Standard bootstrap methods assume data independence, an assumption violated by the temporal autocorrelation intrinsic to blazar variability. To preserve the time-dependent structure of the light curves and mitigate red-noise effects, we adopted a block bootstrap approach (e.g., Künsch 1989). The 17-year light curve is divided into consecutive blocks of size k = 4 bins (corresponding to 120 days), a scale chosen to retain local temporal correlations, while still allowing adequate randomization (e.g., Politis & Romano 1994). Synthetic light curves are then constructed by resampling these blocks with replacements, thus ensuring that the resulting confidence intervals for r and ρ remain robust against red-noise artifacts. To verify that our results are not sensitive to this specific choice, we performed a consistency check by varying the block size over the range k ∈ [2, 6].

3.3. Detrending via singular spectrum analysis

Baseline trends in blazar emission can introduce artificial correlations that mask shorter-term variability. To isolate the QPO and flaring components from this slow evolution, we applied SSA, a nonparametric spectral estimation technique well suited for nonstationary astrophysical time series; we refer to Rico et al. (2025) for complete SSA configuration details.

The SSA decomposition separates the light curve into a slowly varying trend, oscillatory components (including the QPO), and noise. We subtracted the reconstructed trend component from the original light curve to obtain detrended flux values following the method described by Rico et al. (2025). This trend may be part of a longer ∼22 yr QPO (Adhikari et al. 2024). We then repeated the correlation analysis on these detrended data to test whether the spectral trends persist in the absence of secular baseline variations.

3.4. Phase-dependent spectral analysis

Finally, to test the hypothesis of geometric precession versus binary interaction, we folded the photon index data over the reported QPO period of PG 1553+113 (PQPO ≈ 2.2 years; e.g., Ackermann et al. 2015; Abdollahi et al. 2024; Peñil et al. 2025). We assigned a phase ϕ ∈ [0, 1] to each 30-day bin and analyzed the dependence of Γ on ϕ. A significant correlation here would imply that the spectral hardness is modulated by the mechanism driving the periodicity, as predicted by binary shock models.

4. Results and discussion

4.1. Flux and spectral variability

The 17-year light curve of PG 1553+113 (MJD 54 600–61 000) is shown in Fig. 1. The source shows significant variability on multiple timescales, including short-term flares and the well-reported multi-year modulation (e.g., Ackermann et al. 2015). SSA successfully isolates the baseline trend (Fig. 1), thus capturing the slow evolution of the baseline flux without distorting the short-term flaring activity (Rico et al. 2025).

Thumbnail: Fig. 1. Refer to the following caption and surrounding text. Fig. 1.

Gamma-ray light curve of PG 1553+113 (30-day bins, 2008–2025): Fermi-LAT flux (blue circles, 1σ errors) and the SSA baseline (red line) used for detrending.

We investigated the relationship between Fγ and Γ. As shown in Fig. 2, we observe a clear positive correlation, indicating a chromatic softer-when-brighter behavior, contrasting with the harder-when-brighter trend typical of HSPs (Abdo et al. 2010). The Pearson (r) and Spearman (ρ) correlation coefficients for the original dataset are r = 0.57 ± 0.06 and ρ = 0.59 ± 0.06.

Thumbnail: Fig. 2. Refer to the following caption and surrounding text. Fig. 2.

Left: Photon index vs. Fermi-LAT flux, showing the correlation reported in Table 1 and the weighted linear fit Γ = (0.28 ± 0.03)(Fγ × 107)+(1.43 ± 0.02). Right: Photon index vs. QPO phase, showing no correlation (ρ ≈ 0.04, p = 0.57).

A critical concern in blazar studies is that red noise can mimic apparent correlations. To assess the robustness of our results, we compared the correlation coefficients derived from the original data with those obtained after detrending the light curve through removal of the SSA trend, using both standard bootstrap and block bootstrap methods. The corresponding results are summarized in Table 1. The correlation persists after detrending, with coefficients of r ≈ 0.51 and ρ ≈ 0.50. In all cases, the 2σ confidence intervals overlap and exclude zero, providing a statistically robust confirmation that the observed relationship is not a product of chance or red-noise fluctuations. Even when varying the block size over the range k ∈ [2, 6], the resulting 95% confidence interval remains [0.42, 0.73], strictly excluding zero and confirming the robustness of the correlation. Importantly, the softer-when-brighter spectral trend is already present in the original light curve and remains unchanged after removal of the secular baseline component via SSA, indicating that the correlation is not driven by slow baseline evolution, but reflects an intrinsic property of the stochastic variability.

Table 1.

Correlation coefficients and 95% confidence intervals.

4.2. Phase-resolved spectral behavior

Figure 2 shows the photon index folded over the ∼2.2-year QPO period. No phase-dependent spectral modulation is detected. The Spearman rank coefficient between phase and index is negligible (ρ = 0.04) with a p-value of 0.57, which indicates statistical independence.

This result clarifies the flux–spectral behavior discussed in Sect. 4.1. Although a clear chromatic trend is present in the flux variations, it disappears when the data are folded on the QPO phase. The softer-when-brighter behavior therefore arises from fast variability processes, on timescales shorter than the 30-day binning, rather than from the periodic modulation. Rapid spectral evolution is governed by plasma-related mechanisms, while the QPO modulates the flux without introducing additional spectral changes, thus implying an effectively achromatic origin.

In binary scenarios involving disk interactions or periodic shocks, spectral hardening is expected to coincide with the flux maximum (e.g., Sobacchi et al. 2017). The absence of phase-locked spectral evolution disfavors plasma-driven binary shock models and instead supports a geometric origin for the QPO variability. In a jet precession scenario, the modulation is produced by periodic changes in the Doppler factor δ. Since the observed flux scales as F ∝ δp (with p > 3), while the photon index Γ remains invariant in the comoving frame, precession naturally produces strong flux variations with minimal spectral impact (e.g., Rieger 2004). The data are therefore consistent with plasma-driven spectral variability superimposed on a geometrically modulated baseline.

For HSP blazars, the measured correlation coefficients (r ≈ 0.57, ρ ≈ 0.59) indicate a genuine physical coupling given the limited dynamical range of the photon index in these sources. This correlation persists after removing secular trends via SSA detrending (ρ ≈ 0.50 ± 0.07) and remains significant at the 2σ level when temporal autocorrelation is accounted for using block bootstrap resampling. This behavior contrasts with the near-zero correlation between the photon index and the QPO phase, reinforcing the conclusion that fast variability is driven by intrinsic plasma processes, while the QPO is effectively achromatic and consistent with a geometric origin such as jet precession. This interpretation does not exclude the presence of a SMBBH, which could naturally induce precession without producing phase-locked spectral signatures.

The coincidence of high-flux soft-index states reveals a robust softer-when-flaring behavior, indicating that individual flares are intrinsically softer rather than simply reflecting a global flux–index correlation. These episodes likely arise from stochastic plasma processes within the precessing jet, such as plasmoid motion (e.g., Sobacchi et al. 2017) or magnetic reconnection, superimposed on geometric precession. Doppler boosting at QPO maxima amplifies concurrent flares, naturally producing the observed softer-when-brighter trend as the macroscopic imprint of these intrinsically softer flaring events over an otherwise achromatic QPO modulation.

5. Summary and conclusions

In this letter we presented a rigorous spectral analysis of the blazar PG 1553+113 using 17 years of Fermi-LAT data (2008–2025). Our primary goal was to discriminate between plasma-driven and geometric physical scenarios for the reported ∼2.2-year QPO by investigating the coupling between the photon index and the flux modulation. To ensure the validity of our results against the characteristic red noise of blazar light curves, we employed SSA for detrending and a block bootstrap method for uncertainty estimation.

Our main findings are summarized as follows: (1) We confirmed a statistically significant positive correlation between the photon flux and the photon index (ρ = 0.59 ± 0.06). This softer-when-brighter behavior, often observed in flat-spectrum radio quasars, but rare in HSP blazars, suggests that the short-term variability is driven by mechanisms where radiative cooling or particle density variations dominate over particle acceleration during high states. Crucially, this correlation remains stable (ρ ≈ 0.5) even after removing baseline trends and accounting for temporal autocorrelation via block bootstrapping. (2) We found no significant correlation between the photon index and the phase of the QPO (p-value = 0.57, ρ = 0.04). While the total flux is strongly modulated over the ∼2.2-year cycle, the spectral shape of the gamma-ray emission does not vary systematically with the oscillation phase. (3) The absence of phase-dependent spectral modulation provides a strong constraint on the physical origin of the QPO. Models involving intrinsic periodic instabilities or binary black hole disk interactions are expected to induce spectral changes (e.g., hardening due to shocks) coupled to the flux periodicity. Conversely, our results are fully consistent with a geometric origin, such as jet precession. In this scenario, the periodic variation in the Doppler factor modulates the observed gamma-ray flux achromatically without altering the intrinsic electron distribution of the emitting region.

We conclude that PG 1553+113 shows a dual variability nature. Fast stochastic softer-when-brighter chromatic variability is driven by internal plasma processes, while an achromatic QPO arises from a separate large-scale geometric precession that modulates the observed flux on multi-year timescales. Future work should extend this methodology to the broader population of QPO blazar candidates, and to analogous multi-wavelength studies, to assess whether this geometric signature is a common feature of quasi-periodic gamma-ray sources.

Data availability

Data are available in the Fermi-LAT Light Curve Repository.

Acknowledgments

We acknowledge the use of data from the Fermi Large Area Telescope (LAT) Light Curve Repository, supported by the Fermi LAT Collaboration.

References

  1. Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 710, 1271 [NASA ADS] [CrossRef] [Google Scholar]
  2. Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33 [Google Scholar]
  3. Abdollahi, S., Ajello, M., Baldini, L., et al. 2023, ApJS, 265, 31 [NASA ADS] [CrossRef] [Google Scholar]
  4. Abdollahi, S., Baldini, L., Barbiellini, G., et al. 2024, ApJ, 976, 203 [Google Scholar]
  5. Ackermann, M., Ajello, M., Albert, A., et al. 2015, ApJ, 813, L41 [NASA ADS] [CrossRef] [Google Scholar]
  6. Adhikari, S., Peñil, P., Westernacher-Schneider, J. R., et al. 2024, ApJ, 965, 124 [CrossRef] [Google Scholar]
  7. Agazie, G., Arzoumanian, Z., Baker, P. T., et al. 2024, ApJ, 963, 144 [Google Scholar]
  8. Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2015, MNRAS, 450, 4399 [Google Scholar]
  9. Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071 [CrossRef] [Google Scholar]
  10. Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307 [Google Scholar]
  11. Bhatta, G., Zola, S., Stawarz, Ł., et al. 2016, ApJ, 832, 47 [CrossRef] [Google Scholar]
  12. Dinesh, A., Domínguez, A., Paliya, V., et al. 2025, A&A, 703, A162 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Fossati, G., Maraschi, L., Celotti, A., et al. 1998, MNRAS, 299, 433 [Google Scholar]
  14. Jorstad, S. G., Marscher, A. P., Raiteri, C. M., et al. 2022, Nature, 609, 265 [NASA ADS] [CrossRef] [Google Scholar]
  15. Katarzyński, K., Sol, H., & Kus, A. 2001, A&A, 367, 809 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Künsch, H. R. 1989, Ann. Stat., 17, 1217 [Google Scholar]
  17. Láinez, M., Nievas-Rosillo, M., Domínguez, A., et al. 2025, A&A, 700, A229 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Liska, M., Hesp, C., Tchekhovskoy, A., et al. 2018, MNRAS, 474, L81 [Google Scholar]
  19. Malik, Z., Akbar, S., Shah, Z., et al. 2025, MNRAS, 539, 2185 [Google Scholar]
  20. Nievas Rosillo, M., Domínguez, A., Chiaro, G., et al. 2022, MNRAS, 512, 137 [NASA ADS] [CrossRef] [Google Scholar]
  21. Peñil, P., Domínguez, A., Buson, S., et al. 2020, ApJ, 896, 134 [CrossRef] [Google Scholar]
  22. Peñil, P., Westernacher-Schneider, J. R., Ajello, M., et al. 2024, MNRAS, 527, 10168 [Google Scholar]
  23. Peñil, P., Domínguez, A., Buson, S., et al. 2025, ApJ, 980, 38 [Google Scholar]
  24. Peterson, B. M., Wanders, I., Horne, K., et al. 1998, PASP, 110, 660 [Google Scholar]
  25. Politis, D. N., & Romano, J. P. 1994, J. Am. Stat. Assoc., 89, 1303 [Google Scholar]
  26. Ren, H. X., Cerruti, M., & Sahakyan, N. 2023, A&A, 672, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Rico, A., Domínguez, A., Peñil, P., et al. 2025, A&A, 697, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Rieger, F. M. 2004, ApJ, 615, L5 [NASA ADS] [CrossRef] [Google Scholar]
  29. Rybicki, G. B., & Lightman, A. P. 1979, Radiative Processes in Astrophysics (New York: Wiley) [Google Scholar]
  30. Sambruna, R. M., Chou, L. L., & Urry, C. M. 2000, ApJ, 533, 650 [Google Scholar]
  31. Sandrinelli, A., Covino, S., Dotti, M., & Treves, A. 2016, AJ, 151, 54 [NASA ADS] [CrossRef] [Google Scholar]
  32. Sobacchi, E., Sormani, M. C., & Stamerra, A. 2017, MNRAS, 465, 161 [NASA ADS] [CrossRef] [Google Scholar]
  33. Tavani, M., Cavaliere, A., Munar-Adrover, P., & Argan, A. 2018, ApJ, 854, 11 [NASA ADS] [CrossRef] [Google Scholar]
  34. Urry, C. M., & Padovani, P. 1995, PASP, 107, 803 [NASA ADS] [CrossRef] [Google Scholar]
  35. Valtonen, M. J., Nilsson, K., Villforth, C., et al. 2009, ApJ, 698, 781 [NASA ADS] [CrossRef] [Google Scholar]
  36. Vaughan, S., Edelson, R., Warwick, R. S., & Uttley, P. 2003, MNRAS, 345, 1271 [Google Scholar]
  37. Walpole, R. E., Myers, R. H., Myers, S. L., & Ye, K. 2012, Probability& Statistics for Engineers& Scientists, 9th edn. (México: Pearson) [Google Scholar]

All Tables

Table 1.

Correlation coefficients and 95% confidence intervals.

In the text

All Figures

Thumbnail: Fig. 1. Refer to the following caption and surrounding text. Fig. 1.

Gamma-ray light curve of PG 1553+113 (30-day bins, 2008–2025): Fermi-LAT flux (blue circles, 1σ errors) and the SSA baseline (red line) used for detrending.
In the text
Thumbnail: Fig. 2. Refer to the following caption and surrounding text. Fig. 2.

Left: Photon index vs. Fermi-LAT flux, showing the correlation reported in Table 1 and the weighted linear fit Γ = (0.28 ± 0.03)(Fγ × 107)+(1.43 ± 0.02). Right: Photon index vs. QPO phase, showing no correlation (ρ ≈ 0.04, p = 0.57).
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.