Modelling variability power spectra of active galaxies from irregular time series

TL;DR

PIORAN presents a Bayesian Gaussian-process framework to infer broad-band power spectra of active galaxies from irregular time series. The PSD is modeled as a bending power-law and approximated with a small set of fast, basis-function components, enabling likelihood computations in $O(N J^2)$ via the celerite framework. The method is validated with simulations and applied to Ark 564 using 22 years of Swift-XRT and XMM-Newton data, revealing strong evidence for two bends and high-frequency slopes $\alpha_3$ consistent with prior X-ray timing studies, while achieving good inter-instrument calibration. This approach provides a robust, scalable tool for characterizing AGN variability in the time domain and is suitable for large, irregularly sampled datasets from future surveys, with potential extensions to multivariate time series and quasi-periodic features.

Abstract

A common feature of Active Galactic Nuclei (AGN) is their random variations in brightness across the whole emission spectrum, from radio to $γ$-rays. Studying the nature and origin of these fluctuations is critical to characterising the underlying variability process of the accretion flow that powers AGN. Random timing fluctuations are often studied with the power spectrum; this quantifies how the amplitude of variations is distributed over temporal frequencies. Red noise variability -- when the power spectrum increases smoothly towards low frequencies -- is ubiquitous in AGN. The commonly used Fourier analysis methods, have significant challenges when applied to arbitrarily sampled light curves of red noise variability. Several time-domain methods exist to infer the power spectral shape in the case of irregular sampling but they suffer from biases which can be difficult to mitigate, or are computationally expensive. In this paper, we demonstrate a method infer the shape of broad-band power spectra for irregular time series, using a Gaussian process regression method scalable to large datasets. The power spectrum is modelled as a power-law model with one or two bends with flexible slopes. The method is fully Bayesian and we demonstrate its utility using simulated light curves. Finally, Ark 564, a well-known variable Seyfert 1 galaxy, is used as a test case and we find consistent results with the literature using independent X-ray data from XMM-Newton and Swift. We provide publicly available, documented and tested implementations in Python and Julia.

Figures (16)

• Figure 1: True (solid line) and approximated (dashed-dot line) models with the basis functions (dashed lines) for a single bending (blue) and double bending (orange) power spectrum model.
• Figure 2: Residuals (top) and ratios (bottom) between the intended power spectrum model and its approximation as a function of frequency. The median, the $68^\mathrm{th}$ and $95^\mathrm{th}$ percentiles are also shown. The minimum and maximum frequencies of the observed data are shown respectively as black dashed-dot and dotted lines.
• Figure 3: Distributions of the mean, median, minimum and maximum of the frequency-averaged residuals (top) and ratios (bottom) between the true power spectrum model. The extreme values of each distribution are shown as whiskers.
• Figure 4: Likelihood evaluation time in seconds for the approximation method and the direct method as a function of the number points and number of basis functions $J$ for the approximation. The black squares show regression with direct method scaling as $\mathcal{O}(N^3)$ (see \ref{['sec:fftmethod']}). The circles and triangles respectively represent the Python tinygp and Julia Pioran.jl codes where colours encode the number of SHO basis functions.
• Figure 5: Distributions of the posterior samples for one simulated time series. The median of the distributions is shown by a blue vertical line while the true value is shown by a dotted magenta vertical lines. The prior distributions are shown with dashed red lines.