JAVELIN: A Probabilistic AGN Lag Estimator

• The paper presents a fully probabilistic framework that replaces traditional cross-correlation methods with a DRW model and convolution transfer functions for AGN lag estimation.
• It details a mathematical formulation using Gaussian processes and explicit transfer function kernels to yield posterior distributions and error covariances essential for AGN structural inference.
• The approach enables accurate emission-line and continuum lag mapping while addressing challenges like irregular sampling, seasonal gaps, and cadence-induced bias.

JAVELIN (Just Another Vehicle for Estimating Lags in Nuclei) is a fully probabilistic framework for measuring reverberation lags in active galactic nuclei (AGN), with applications ranging from emission-line echo mapping to continuum inter-band lag analyses in quasars. It replaces traditional cross-correlation techniques with a model-driven approach, using damped random walk (DRW) processes for stochastic continuum variability and explicit convolution transfer functions for emission-line or multi-band responses. JAVELIN robustly accommodates irregular sampling, correlated errors, seasonal gaps, and multiple responding light curves by fitting the data directly within a Gaussian-process likelihood framework, thereby yielding parameter posterior distributions and error covariances essential for AGN structural inference (Zu et al., 2010). The method is widely implemented for both emission-line and continuum lag studies via the open-source JAVELIN Python package.

1. Mathematical Formulation of Quasar Variability

JAVELIN models AGN light curve variability as a stationary Gaussian process, specifically a Damped Random Walk (DRW) or Ornstein–Uhlenbeck process. For a continuum light curve $s_c(t)$, the covariance between two epochs $t_i$ and $t_j$ is

$$\langle s_c(t_i)\,s_c(t_j)\rangle = \sigma^2 \exp\left(-\frac{|t_i-t_j|}{\tau_d}\right)$$

where $\sigma$ is the long-term process variance and $\tau_d$ is the relaxation timescale. This parameterization ensures that short-timescale structure ($\Delta t \ll \tau_d$) grows as $\hat{\sigma} \sqrt{\Delta t}$, while long-timescale variability saturates at $\sigma$ (Zu et al., 2010, Jiang et al., 2016, Secunda et al., 2023, Jha et al., 2021). The statistical likelihood for a vector of flux measurements $x$ with mean subtraction is

$$\ln \mathcal{L} = -\frac{1}{2} N \ln(2\pi) - \frac{1}{2} \ln|C| - \frac{1}{2} x^T C^{-1}x$$

where $C$ is the DRW covariance matrix plus measurement noise (Secunda et al., 2023).

2. Transfer Function Formalism and Lag Recovery

Emission-line and multi-band continuum responses are modeled as scaled, smoothed, and displaced copies of the driving continuum via a linear convolution ansatz:

$$s_\ell(t) = \int_{-\infty}^\infty \Psi(\tau) s_c(t-\tau) d\tau$$

JAVELIN typically adopts an explicit top-hat kernel for the transfer function,

$$\Psi(\tau) = \begin{cases} A/(t_2-t_1), & t_1 \leq \tau \leq t_2 \ 0, & \text{otherwise} \end{cases}$$

where $A$ is the responsivity, $t_1$ and $t_2$ define the leading and trailing edges of the top-hat, giving mean lag $\langle \tau \rangle = (t_1 + t_2)/2$ and width $\Delta \tau = t_2 - t_1$ (Zu et al., 2010, Jha et al., 2021, Jiang et al., 2016). In practice, extensions to the transfer function include Gaussian, two-component top-hat, or Gaussian + skew kernels to accommodate more complex or multi-peaked responses (Secunda et al., 2023).

3. Likelihood Construction and Bayesian Inference

The full covariance matrix for continuum and multiple response light curves is constructed by combining signal covariance $S$ with measurement error covariance $N$:

$C = S + N$

where $S_{ij}$ includes DRW autocovariance and convolution with the transfer function; $N_{ij}$ may accommodate non-diagonal, epoch-wise correlated measurement errors. The cross-covariance between a line and continuum is

$$\langle s_\ell(t_i) s_c(t_j) \rangle = \int \Psi(t_i - t') \langle s_c(t') s_c(t_j) \rangle dt'$$

and the line autocovariance is

$$\langle s_\ell(t_i) s_\ell(t_j) \rangle = \iint \Psi(t_i - t') \Psi(t_j - t'') \langle s_c(t') s_c(t'') \rangle dt' dt''$$

Parameter inference is performed by maximizing the marginalized log-likelihood using the downhill simplex ("amoeba") algorithm, followed by Markov Chain Monte Carlo (Metropolis–Hastings or Differential Evolution MCMC) to sample the full joint posterior. Priors are chosen log-uniform on DRW parameters and uniform on lag kernel parameters. Posterior credible intervals for recovered lag(s), responsivities, and error covariances are extracted from marginalized posteriors (Zu et al., 2010, Jha et al., 2021, Secunda et al., 2023, Jiang et al., 2016).

4. Application to Emission-Line Reverberation Mapping and Multi-Band Continuum Lags

JAVELIN is employed to measure BLR reverberation lags (e.g. Hβ, Hα), multi-line, and velocity-resolved delay maps. It robustly addresses irregular sampling, correlated errors, and seasonal gaps, outperforming conventional interpolated cross-correlation (ICCF) and discrete correlation function (DCF) methods, which are susceptible to interpolation error, binning artifacts, and aliasing. In joint fits, JAVELIN recovers more sharply peaked posteriors and suppresses false lag solutions due to cadence aliases (Zu et al., 2010, Secunda et al., 2023).

In continuum reverberation mapping, JAVELIN measures interband continuum lags (e.g. $\tau_{g-r}$, $\tau_{g-i}$, $\tau_{g-z}$), elucidating accretion disk size and structure. These lags are consistently found to be several times larger than standard Shakura–Sunyaev thin disk model predictions, revealing the necessity for revised radial temperature profiles or extra reprocessing physics (Jiang et al., 2016, Jha et al., 2021). Simultaneous "disk-mode" fitting constrains global disk size ($R_0$) as a function of wavelength ($$\tau(\lambda) = \tau_0[(\lambda/\lambda_0)^{4/3}-1]$$), with empirical sizes systematically larger than analytic predictions—average disk size being $\sim 3.9$ times Shakura–Sunyaev analytic values (Jha et al., 2021).

5. Recovery of Viscous (Negative) Lags in Quasar Photometry

JAVELIN is particularly suited for detection of long, negative lags in quasar photometry associated with the viscous timescale in accretion disks. Simulating LSST-like light curves with injected long lags ($\tau_{\text{in}} = -50, -130, -400$ days), JAVELIN recovers the median lag values within statistical uncertainties (e.g. $-50^{+1}_{-1}$ d, $-130^{+0}_{-0}$ d, $-410^{+10}_{-20}$ d for various injection cases) (Secunda et al., 2023). Robustness is noted even at low signal-to-noise (SNR ∼ 30) and with non-DRW driving processes. When both short (light-travel) and long (viscous) lags are present, JAVELIN reliably recovers the long lag component as long as its amplitude is $\gtrsim 0.1$ of the short-lag signal; for lower amplitude ratios, bias towards zero is introduced, making only upper limits feasible. A plausible implication is that long-separation LSST cadences or gap-filling external visits are required for unambiguous recovery of long negative lags.

6. Comparison with Alternative Lag Estimation Techniques

JAVELIN outperforms ICCF and DCF for irregularly sampled, seasonally gapped light curves. In alias-prone cases, traditional methods yield multiple peaks, whereas JAVELIN identifies the correct lag with lower uncertainty, suppressing cadence-induced false solutions (Zu et al., 2010, Secunda et al., 2023). For dual-lag (short + long) scenarios, ICCF reliably detects both only when long/short amplitude ratio $\gtrsim 0.2$, and Fourier maximum-likelihood methods are more effective for isolating short-lags. The Von-Neumann estimator can recover long lags on baseline cadences if one lag dominates but is otherwise single-lag constrained (Secunda et al., 2023).

7. Practical Implementation, Extensions, and Computational Considerations

JAVELIN’s inference pipeline involves fitting the continuum alone for DRW parameters (broad log-priors), followed by joint fits to continuum + response curves with Gaussian priors on DRW parameters and uniform priors on transfer function parameters ($A, t_1, t_2$). Computation scales as $\mathcal{O}(K^3)$ in the number of epochs due to matrix inversion in the likelihood evaluation; sparse or tridiagonal algorithms can accelerate the continuum-only stage. JAVELIN supports multi-line and multi-band fits, providing covariance matrices between lag estimates—critical for interpreting velocity–delay maps and for simultaneous modeling of luminosity-dependent lags ($\tau \propto L^\alpha$) (Zu et al., 2010). Priors on response width and amplitude, multi-component kernels, and direct "disk-mode" fitting can be invoked to address broad or multi-peaked transfer functions (Secunda et al., 2023, Jha et al., 2021).

Summary tables of methods and performance metrics demonstrate JAVELIN’s robust applicability to contemporary photometric survey data, enabling detailed AGN structure inference, echo mapping, and accretion disk modeling in regimes previously inaccessible to interpolation-based lag estimators (Secunda et al., 2023).

---

References

• "An Alternative Approach To Measuring Reverberation Lags in Active Galactic Nuclei" (Zu et al., 2010)
• "Negative Lags on the Viscous Timescale in Quasar Photometry and Prospects for Detecting More with LSST" (Secunda et al., 2023)
• "Accretion Disk Sizes from Continuum Reverberation Mapping of AGN Selected from the ZTF Survey" (Jha et al., 2021)
• "Detection of Time Lags Between Quasar Continuum Emission Bands based on Pan-STARRS Light-curves" (Jiang et al., 2016)