Power density spectrum of nonstationary short-lived light curves

The power density spectrum of a light curve is often calculated as the average of a number of spectra derived on individual time intervals the light curve is divided into. This procedure implicitly assumes that each time interval is a different sample function of the same stochastic ergodic process. While this assumption can be applied to many astrophysical sources, there remains a class of transient, highly nonstationary and short-lived events, such as gamma-ray bursts, for which this approach is often inadequate. The power spectrum statistics of a constant signal affected by statistical (Poisson) noise is known to be a chi2(2) in the Leahy normalisation. However, this is no more the case when a nonstationary signal is also present. As a consequence, the uncertainties on the power spectrum cannot be calculated based on the chi2(2) properties, as assumed by tools such as XRONOS powspec. We generalise the result in the case of a nonstationary signal affected by uncorrelated white noise and show that the new distribution is a non-central chi2(2,lambda), whose non-central value lambda is the power spectrum of the deterministic function describing the nonstationary signal. Finally, we test these results in the case of synthetic curves of gamma-ray bursts. We end up with a new formula for calculating the power spectrum uncertainties. This is crucial in the case of nonstationary short-lived processes affected by uncorrelated statistical noise, for which ensemble averaging does not make any physical sense.

💡 Research Summary

The paper addresses a fundamental limitation in the standard practice of computing power density spectra (PDS) for astrophysical time series. Traditionally, a light curve is divided into several contiguous intervals, each interval’s Fourier spectrum is calculated, and the resulting spectra are averaged. This approach implicitly assumes that each interval is an independent realization of the same stationary, ergodic stochastic process. While this assumption holds for many persistent sources, it breaks down for highly non‑stationary, short‑lived transients such as gamma‑ray bursts (GRBs). In such cases, averaging over intervals is physically meaningless; only a single PDS computed over the entire observation is appropriate.

The authors start from the Leahy‑normalised power definition
(P_j = \frac{2}{N_{\rm ph}}|a_j|^2 = \frac{2}{N_{\rm ph}}\sum_{k,l} x_k x_l e^{2\pi i (k-l)j/N}),
where (x_k) are the observed counts, (N_{\rm ph} = \sum_k \sigma_k^2) is the expected total variance (or total counts for Poisson statistics), and (a_j) are the discrete Fourier transform (DFT) amplitudes. They treat the observed series as a deterministic signal (\eta_k) contaminated by uncorrelated white noise with variance (\sigma_k^2).

Gaussian noise case
When each (x_k) follows a normal distribution (N(\eta_k,\sigma_k^2)), the power can be written as a quadratic form (P_j = X^{\rm T} A X), where (A) is a symmetric matrix built from cosine terms. By analysing the matrix product (A\Sigma A) (with (\Sigma) the covariance matrix), they show that, for equal variances or for large total counts (N_{\rm ph}), the condition (A\Sigma A \approx A) holds. Consequently, (P_j) follows a non‑central chi‑square distribution with 2 degrees of freedom (or 1 degree at the Nyquist frequency):
(P_j \sim \chi^2_{r}(\lambda)), with (r=2) (or (r=1) at Nyquist) and non‑centrality parameter (\lambda = P(\eta)_j), i.e. the power spectrum of the underlying deterministic signal. The mean and variance are then
(\mathbb{E}