LITMUS: Bayesian Lag Recovery in Reverberation Mapping with Fast Differentiable Models

TL;DR

LITMUS addresses the challenge of robust lag recovery in reverberation mapping by casting lag inference as Bayesian forward modelling with a damped random walk description of AGN variability. Its core innovation, the Laplace Quadrature, maps the posterior over lags via Gaussian slices and computes model evidence to enable principled Bayes factors for lag significance, while grid smoothing and preconditioning enhance speed and reliability. The framework is modular, supports alternative statistical models, and includes additional algorithms like SVI Quadrature and Nested Sampling for cross-validation. Validation on OzDES-like mocks shows superior discrimination of true lags from aliases and strong suppression of false positives compared to JAVELIN and ICCF, enabling more reliable RM in large surveys. The approach promises to extend RM constraints across redshift and luminosity space and to other single-lag timing problems, with potential insights into BLR structure, MgII/Hβ relationships, and beyond.

Abstract

Reverberation mapping is a technique in which the mass of a Seyfert I galaxy's central supermassive black hole is estimated, along with the system's physical scale, from the timescale at which variations in brightness propagate through the galactic nucleus. This mapping allows for a long baseline of time measurements to extract spatial information beyond the angular resolution of our telescopes, and is the main means of constraining supermassive black hole masses at high redshift. The most recent generation of multi-year reverberation mapping campaigns for large numbers of active galactic nuclei (e.g. OzDES) have had to deal with persistent complications of identifying false positives, such as those arising from aliasing due to seasonal gaps in time-series data. We introduce LITMUS (Lag Inference Through the Mixed Use of Samplers), a modern lag recovery tool built on the "damped random walk" model of quasar variability, built in the autodiff framework JAX. LITMUS is purpose built to handle the multimodal aliasing of seasonal observation windows and provides evidence integrals for model comparison, a more quantified alternative to existing methods of lag validation. LITMUS also offers a flexible modular framework for extending modelling of AGN variability, and includes JAX-enabled implementations of other popular lag recovery methods like nested sampling and the interpolated cross correlation function. We test LITMUS on a number of mock light curves modelled after the OzDES sample and find that it recovers their lags with high precision and a successfully identifies spurious lag recoveries, reducing its false positive rate to drastically outperform the state of the art program JAVELIN. LITMUS's high performance is accomplished by an algorithm for mapping the Bayesian posterior density which both constrains the lag and offers a Bayesian framework for model null hypothesis testing.

Figures (10)

• Figure 1: A demonstration of the sort of light curves that GP modelling can reconstruct from observations. For some time-series observations (error bars) a particular GP models the entire family of underlying light curves that exhibit the power spectral density of the GP, conditioned on how well they fit the observations. In this example the light curves is fit as a DRW with $\tau=200 \mathrm{d}$ and $\sigma=1$, both in arbitrary units for this demonstrative example. The shaded regions represent the $1$ and $2 \sigma$ contours of the distribution of all such walks.
• Figure 1: An exaggerated demonstration of the grid smoothing algorithm for a simple multimodal function using $\alpha = 0.8$ up to $j=5$ iterations with $32$ points. The top panel shows the true distribution (black) with its estimate from the first evenly spaced grid (red) and the final smoothed grid (blue). The bottom panel shows how the spacing of the grid updates over each iteration, progressing from top to bottom, with the first and last iterations coloured for emphasis, and gray dots representing samples from previous iterations. The initial spacing is so coarse that it misses much of the detail of the left mode, and cuts off the right-mode entirely. By final iteration, the estimate of the mode is significantly more accurate.
• Figure 1: Histogram of the run-times for the five fitting methods in LITMUS over all mocks, using the fitting parameters described in Table \ref{['tab: fitting_params']}. The ICCF method, which requires no matrix inversion owing to its absence of GP fitting, is consistently the fastest. The Laplace Quadrature can perform very fast except for cases where it get stuck optimising at a new test lag when the local optimum changes quickly over the lag axis. The SVI Quadrature has a similar issue, but runs overall somewhat slower.
• Figure 2: A demonstration of the source of the aliasing problem, specifically in the context of a parametric GP model. Top shows mock data with cadence, measurement uncertainty and baseline similar to OzDES with a DRW timescale of $\tau=200 \mathrm{d}$ and a true lag of $\Delta t=360 \mathrm{d}$. From left to right the sub-panels show lags being tested at $\Delta t=0 \mathrm{d}$, $180 \mathrm{d}$ and $360 \mathrm{d}$. The left panel is clearly a bad fit as near simultaneous observations are in clear tension, and the right panel is a clear good fit as we see very little tension. The middle panel, corresponding to the first aliasing peak, is an ambiguously good fit; the lack of overlap means we cannot observe clear tensions between the light curves. The bottom panel shows the (un-normalised) log-natural of the posterior distribution, with all non-lag parameters fixed at their true values. At "on-season" lags (un-shaded) we can easily reject bad fits, and so the posterior is extremely low. During the off-season lags (blue shading) there are local optima arising from the ambiguity. The mode associated with the true lag (red dot) is clearly defined and dominates over aliasing modes, with the rest of the posterior being $<1 \%$ of the maximum posterior density in this well behaved, high SNR example. Even so, the posterior still suffers from the rough geometry and multimodality that introduces numerical challenges in navigating it.
• Figure 3: A demonstration of the failure mode of the Affine-Invariant Ensemble Sampler (AIES), the MCMC proposal algorithm used by emcee, in multi-modal distributions. Both top and bottom panels are posterior distributions generated from the same mock data with a true lag at $\Delta t = 854 \mathrm{d}$ (dashed line), with the bottom panel being the result from the AIES, the same MCMC sampler as JAVELIN, while the top is found from exhaustive sampling of the prior range. The AIES estimate for the posterior has produced an aliasing peak at $\Delta t = 540 \mathrm{d}$ where none truly exists due to its ensemble of live sampling points becoming pinned at this minor mode.