Time delays and stationarity in quasar light curves

1. Problem Statement

The accurate measurement of time delays between multiple images of gravitationally lensed quasars is critical for time-delay cosmography, a method used to measure the Hubble constant ($H_0$) and resolve the current "Hubble tension." However, this measurement faces significant challenges:

• Data Characteristics: Quasar light curves are irregularly sampled, noisy, and often contain seasonal gaps.
• Complex Variability: Quasar variability is largely stochastic (often modeled as red noise) but exhibits complex behaviors, including long-term trends and non-stationarity (changes in statistical properties over time).
• Model Limitations: Traditional methods often assume stationarity or use simple mean functions (e.g., constant or linear). These assumptions can lead to biased time delay estimates, particularly when extrapolating data. Furthermore, standard Gaussian Process (GP) models with fixed kernels may overfit or fail to capture the underlying physical drift, leading to pathological likelihood structures where the posterior favors unphysical large time delays.
• Microlensing: Extrinsic variability caused by microlensing must be disentangled from the intrinsic quasar variability.

2. Methodology

The authors propose a fully Bayesian framework using Marginalised Gaussian Processes (GPs) to jointly model time delays, stationarity, and microlensing.

A. The Flux Variability Model

The model decomposes the observed light curve $y(t)$ into three components:

1. Deterministic Non-Stationary Drift (Mean Function):
   • Instead of a fixed mean, the authors use a flexknot model: a piecewise linear spline where the knot locations ($t$) and values ($y$) are parameters to be sampled.
   • This allows the model to adaptively learn long-term trends without assuming a specific functional form (e.g., polynomial).
   • The number of knots ($n_{fk}$) is treated as a hyperparameter.
2. Stationary Stochastic Variability (Kernel):
   • The residuals around the mean are modeled using stationary kernels.
   • The authors marginalize over a wide family of kernels: Matérn (Exponential/E, M3/2, M5/2, M7/2), Squared Exponential (SE), Rational Quadratic (RQ), and Spectral Mixture (SM) kernels.
   • This approach avoids the bias of pre-selecting a single kernel (e.g., the Damped Random Walk/Exponential kernel).
3. Microlensing:
   • Modeled as an additional GP with a Squared Exponential kernel for each image, added to the diagonal of the covariance matrix.

B. Multi-Image Joint Modeling

• The model treats the light curves of the multiple lensed images (e.g., A, B, C, D) as shifted instances of a single latent GP.
• Time Delays ($\Delta t$): The time shifts between images are parameters.
• Magnitude Offsets ($\Delta m$): Constant magnitude differences between images are parameters.
• The covariance matrix is constructed as a block matrix, where off-diagonal blocks represent the correlation between shifted time series.

C. Inference and Marginalization

• Nested Sampling (NS): The authors use the PolyChord algorithm to perform Bayesian inference. NS is chosen for its ability to handle multimodal posteriors, hard constraints (like sorted knot positions), and singular covariance matrices.
• Marginalization: Instead of selecting the "best" model, the framework marginalizes over:
  • Kernel choice.
  • Mean function complexity (number of knots).
  • Noise parameters.
  • This yields a robust posterior for time delays that accounts for model uncertainty.
• Convergence Protocol: The authors identify a likelihood pathology where NS runs can get stuck in local modes favoring large time delays. They implement a convergence protocol involving increasing the number of "live points" and, for specific datasets, narrowing the prior on time delays based on previous literature to ensure convergence.

3. Key Contributions

1. Novel GP Framework: Introduction of a hierarchical GP model that separates non-stationary drift (flexknots) from stationary stochastic variability, allowing for data-driven detection of non-stationarity.
2. Robust Time Delay Estimation: Demonstration that marginalizing over kernel and mean function choices yields time delays consistent with prior work but with a more rigorous accounting of systematic uncertainties.
3. Stationarity Testing: A Bayesian method to test for stationarity in quasar light curves by comparing the evidence of models with different mean function complexities.
4. Pathology Mitigation: Identification and mitigation of a specific likelihood pathology in time delay inference that biases results toward large delays, proposing a practical convergence protocol for Nested Sampling.
5. Model Selection Insights: Systematic comparison of various GP kernels (Matérn, SM, etc.) for quasar data, moving beyond the standard Damped Random Walk assumption.

4. Results

The framework was applied to three benchmark lensed quasars: WFI J2033 – 4723, B 1608 + 656, and HE 0435 – 1223.

A. Stationarity and Kernel Preferences

• WFI J2033 – 4723: The data is consistent with stationarity. The mean function is effectively constant. The stochastic component strongly favors the Matérn-3/2 (M3/2) kernel, indicating a non-Markovian process (memory of previous values).
• B 1608 + 656: Strong evidence for non-stationarity. The mean function shows significant long-term drift. The stochastic component favors the Exponential (E) kernel (Markovian/Damped Random Walk).
• HE 0435 – 1223: Strong evidence for non-stationarity. The mean function captures large-amplitude variations. The stochastic component favors the Matérn-3/2 (M3/2) kernel.
• Spectral Mixture (SM) Kernels: Multi-length-scale SM kernels were generally disfavoured, suggesting that the long-term variations are better captured by the deterministic flexknot mean function rather than by adding more length scales to the stochastic kernel.

B. Time Delays

• The inferred time delays are robust to model assumptions (kernel and mean function choices).
• The results are consistent with previous studies (e.g., Hojjati et al. 2013) within error bars, validating the framework.
• Uncertainty: Marginalization over models generally reduces or maintains the precision of time delay estimates compared to fixed-model approaches, while providing a more honest estimate of uncertainty.

C. Residual Analysis

• After subtracting the inferred non-stationary mean function, the residuals for all three quasars were found to be consistent with a Gaussian distribution (in magnitude units).
• This validates the core assumption that the underlying stochastic process is Gaussian, supporting the use of GP models.
• A minor artifact was identified: the model showed a slight preference for fitting the arbitrarily chosen "reference" light curve (Image A) better than others, suggesting a need for permutation-invariant models in future work.

5. Significance

• Cosmological Impact: By providing a more robust and less biased method for measuring time delays, this work strengthens the reliability of time-delay cosmography, which is crucial for resolving the Hubble tension.
• Physical Insight: The distinction between stationary and non-stationary quasars, and the preference for different kernels (Markovian vs. non-Markovian), offers new insights into the physical mechanisms driving accretion disk variability (e.g., thermal fluctuations vs. other processes).
• Methodological Advancement: The paper sets a new standard for handling irregularly sampled astronomical time series by combining flexible mean functions with kernel marginalization. It highlights the importance of rigorous convergence testing in Nested Sampling for complex likelihood surfaces.
• Future Directions: The authors suggest extending this to multi-wavelength data (to constrain physical models of the accretion disk) and developing permutation-invariant models to eliminate reference-image bias.