Time delays and stationarity in quasar light curves

TL;DR

Time-delay cosmography with gravitationally lensed quasars requires robust inference of intrinsic variability and inter-image delays. The authors introduce a fully Bayesian marginalised-GP framework that jointly models a non-stationary drift via a flexknot mean and stationary stochasticity via Matérn and Spectral Mixture kernels, marginalising over kernel and mean choices using Bayesian evidence $p(\mathcal{D}|\mathcal{M})$ to obtain posterior model probabilities. Applied to WFI J2033-4723, B1608+656, and HE 0435-1223, they report system-specific stationarity properties and kernel preferences, with time delays broadly consistent with prior measurements within the uncertainties. The study demonstrates that time-delay inferences are robust to model assumptions, emphasizes convergence diagnostics for nested sampling, and outlines directions for permutation-invariant mean models and multi-wavelength extensions that could improve physically interpretable constraints on quasar variability.

Abstract

We present a fully Bayesian framework for time delay inference and stationarity tests in quasar light curves using marginalised Gaussian processes. The model separates a deterministic, non-stationary drift (piecewise linear mean) from stationary stochastic variability (Matérn and Spectral Mixture kernels), and jointly models multiple images with per-image microlensing. Bayesian evidence and parameter posteriors are obtained via nested sampling and marginalised over model choices. Applied to the quasars WFI J2033 - 4723, B 1608 + 656, and HE 0435 - 1223, we find strong evidence for non-stationarity in B 1608 + 656 and HE 0435 - 1223, while WFI J2033 - 4723 is consistent with stationarity. The stochastic component favours an Markovian exponential kernel for B 1608 + 656 and a non-Markovian Matérn-$\frac32$ kernel for WFI J2033 - 4723 and HE 0435 - 1223. Multi-length-scale Spectral Mixture kernels are disfavoured. Time delays are shown to be robust to model assumptions and consistent with prior work within the error. We further identify and mitigate a likelihood pathology which biases toward large delays, providing a practical nested sampling convergence protocol.

Figures (4)

• Figure 1: Evidence $\log Z$, posterior-averaged log-likelihood $\langle\log L\rangle_P$ and -divergence $D_\mathrm{KL}$ for the three quasar light curve data sets. The three quantities are related by Occam's razor equation, $\log Z=\langle\log L\rangle_P- D_\mathrm{KL}$. Each coloured square in a subplot corresponds to a fit, performed with a given number of knots $n_\mathrm{fk}$ in the flexknot mean function and a given kernel. Zero and one knots correspond to a fixed and constant mean (with the constant a parameter), respectively. Note that the numbers in brackets denotes the standard deviation on the last digits. As shown by the -divergence, model complexity increases roughly from the top left to the bottom right corner. The goodness of fit, $\langle\log L\rangle_P$, roughly shows the opposite trend, with more complex models producing better fits. As the balance of these two, the evidence $\log Z$ shows an overall preference to the and kernels. Non-stationarity (more than one knot) is preferred for B 1608 + 656 and HE 0435 -- 1223, whereas there is no clear indication of non-stationarity for WFI J2033 -- 4723.
• Figure 2: (left column) and mean function (right column) posterior predictive distributions for the three quasar light curve data sets, each consisting of a number of light curves: A, B, C and (if existing) D. The is the sum of the deterministic mean function and the stochastic contribution from the kernel, which are fit jointly to the data. Plotting just the mean function therefore extracts the non-stationarity of the . The fits are marginalised over all choices in Figure \ref{['fig:real-data-evidences']}. Visually, the mean function of WFI J2033 -- 4723 is consistent with a constant, whereas for B 1608 + 656 and HE 0435 -- 1223, the mean function is clearly non-stationary. Note that the full on the right reverts to the mean function outside the observation window, i.e. extrapolation is dictated by the mean function. The time axis is set to zero at the first observation.
• Figure 3: Corner plots of the time delay posteriors for the three data sets. The posteriors are marginalised over all model choices in Figure \ref{['fig:gp-posterior-plots']}. For comparison, the inferred values (mean and standard deviation) from hojjati2013robust are overplotted. All posteriors are unimodal and well-constrained, with the exception of $\Delta t_\mathrm{AB}$ of B 1608 + 656, which displays a weakly protruding secondary peak. The numerical means and standard deviations are also shown in Table \ref{['tab:timedelays']}.
• Figure 4: Scatter of the light curve data points around the mean function of the fit in units of magnitude. The distribution of the residuals is expected to follow a normal distribution. For comparison, a Gaussian distribution with the same mean and standard deviation is shown as well. For WFI J2033 -- 4723 and B 1608 + 656, the agreement in the tails is better than in the bulk. For HE 0435 -- 1223, the observed deviation from the exact Gaussian on the left half is attributed to a mis-fit resulting from lack of permutation-invariance of the model.