Generation of Correlated Time Series for X-ray Astronomy Applications

Technical Summary: Generation of Correlated Time Series for X-ray Astronomy Applications

Problem Statement
Fourier-domain techniques, including power spectral densities (PSDs), cross-power spectral densities (CSDs), coherence ($\gamma^2$), and phase lags ($\phi$), are fundamental for studying accretion physics in X-ray binaries (XRBs) and active galactic nuclei (AGN). Recent advances, such as multi-Lorentzian fitting of joint power and cross-spectra, have revealed "hidden" variability components (e.g., quasi-periodic oscillations) that are invisible in power spectra alone. However, validating these complex analysis techniques and understanding the statistical behavior of Fourier-domain estimators requires synthetic time series with realistic, prescribed cross-spectral properties.

Existing methods, such as the Timmer & König (1995) algorithm, can generate time series with a specific power spectrum but produce statistically independent realizations. When averaged, these independent series exhibit zero coherence and zero net phase lag, rendering them unsuitable for testing methods that rely on correlated light curves with non-trivial phase relationships. While extensions exist for non-Gaussian statistics, a general method for generating pairs of time series with arbitrary, frequency-dependent coherence and phase lag profiles was lacking.

Methodology
The authors present an algorithm that generalizes the Timmer & König (1995) approach to generate pairs of time series ($x$ and $y$) with arbitrary power spectra ($P_X, P_Y$), coherence functions ($\gamma^2$), and phase lag profiles ($\phi$).

The core of the method involves decomposing the "dependent" time series ($y$) into two components:

1. An incoherent component ($y_u$): Uncorrelated with the reference series.
2. A coherent component ($y_c$): Linearly correlated with the "reference" time series ($x$).

In the Fourier domain, the dependent series $Y$ is constructed as:
$$Y = Y_u + R(\nu)X$$
where $X$ is the Fourier transform of the reference series, and $R(\nu)$ is a complex linear transfer function.

The authors derive the specific forms for the normalization constant and the transfer function required to match target spectral shapes:

• Transfer Function ($R$): Derived to satisfy the target coherence and phase lag:
  $$R(\nu) = \sqrt{\frac{P_Y(\nu)\gamma^2(\nu)}{P_X(\nu)}} e^{i\phi(\nu)}$$
• Normalization Constant ($K$): Derived to ensure the dependent series preserves its target power spectrum $P_Y$ despite the addition of the coherent term:
  $$K = \sqrt{\frac{P_Y - P_X|R|^2}{2}}$$
  Here, the incoherent part $Y_u$ is generated as $K(H_r + iJ_r)$, where $H_r$ and $J_r$ are independent Gaussian random variables.

The algorithm proceeds by drawing Gaussian random vectors for the reference and incoherent components, applying the derived transfer function and normalization, and performing an inverse Fourier transform to obtain the time-domain realizations.

Key Contributions

1. Synthetic Data Generation: The paper provides a complete algorithm to generate pairs of time series with independently specified power spectra, coherence, and phase lags. This allows for the creation of "ground truth" datasets for testing cross-spectral analysis pipelines.
2. Analytic Variance and Covariance Derivations: Using the delta method, the authors derive approximate analytic expressions for the variance of coherence and phase lag estimators as a function of the number of averaged segments ($n$):
   • $\text{Var}(\hat{\gamma}^2) \approx \frac{2\gamma^2(1-\gamma^2)^2}{n}$
   • $\text{Var}(\hat{\phi}) \approx \frac{1-\gamma^2}{2n\gamma^2}$
     They further construct the full $6 \times 6$ covariance matrix relating the base quantities (real/imaginary cross-spectrum, power spectra) to the derived quantities (coherence, phase lag).
3. Bias Characterization: The work explicitly addresses the known bias in coherence estimators (Nuttall & Carter 1976), confirming that the estimator systematically overestimates coherence, particularly at low coherence values and small segment counts.

Results
The authors validate the method through two primary tests:

• Ensemble Validation: Applying the algorithm to a two-Lorentzian model with distinct phase lags ($\phi_0 = 0.15$ rad, $\phi_1 = -0.8$ rad). The ensemble-averaged power spectra, coherence, and phase lag profiles of the generated data closely matched the input models across four decades in frequency. The coherence profile correctly showed drops at frequencies where multiple components contributed comparable power, mimicking observed behavior in sources like Cyg X-1.
• Monte Carlo Simulations: Using 50,000 realizations, the authors verified that the observed scatter in coherence and phase lag estimates matched the derived analytic variance formulas. The simulations also confirmed the theoretical bias in the coherence estimator and demonstrated that the phase lag estimator is uncorrelated with the power spectra, while the coherence estimator exhibits weak correlations with the base spectral quantities.

Significance and Limitations
The paper positions this method as a foundational tool for the rigorous testing of cross-spectral analysis frameworks, particularly those involving multi-Lorentzian fitting and reverberation mapping. By providing well-characterized statistical properties and analytic expressions for estimator variance, the method aids in experimental design (e.g., determining required data segment counts for a target precision).

The authors explicitly note several limitations and scope boundaries:

• Noise and Detectors: The method generates noiseless time series. Realistic simulations require the post-hoc addition of Poisson noise and the modeling of detector effects (e.g., pile-up, dead time, time-binning), which can alter coherence and phase lag.
• Linearity: The construction captures only linear correlations between bands, neglecting higher-order terms (e.g., non-linear phase coupling).
• Statistics: The generated time series are Gaussian, whereas accreting systems often exhibit log-normal flux distributions.
• Physical Interpretation: The derived transfer function $R(\nu)$ is a mathematical construct to reproduce observables, not necessarily a physical response function. In cases where the dependent series leads the reference (positive phase lag), $R$ encodes a causality-violating relationship; thus, it should not be interpreted as a physical model of the system's response.

The work concludes that robust simulation tools with characterized statistical properties are essential as cross-spectral techniques become more sophisticated in X-ray astronomy.