Timing Gamma-ray Pulsars using Gibbs Sampling

Technical Summary: Timing Gamma-ray Pulsars using Gibbs Sampling

Problem Statement
Timing analyses of gamma-ray pulsars in the Fermi Large Area Telescope (LAT) data set offer unique astrophysical insights, including the study of timing noise (TN) in young pulsars, orbital period variations (OPVs) in binary systems, and the stochastic gravitational wave background (GWB). However, existing radio pulsar timing methods, which rely on Gaussian approximations of time-of-arrival (ToA) uncertainties, are not directly applicable to gamma-ray data. Gamma-ray data consists of discrete photon arrival times with uncertain origins (signal vs. background) and non-Gaussian likelihoods due to the multi-peaked, broad nature of gamma-ray pulse profiles.

Previous approaches have struggled with these limitations:

1. ToA-based methods: Integrating data to create pulse profiles and extracting ToAs requires long integration times, often longer than the timescales of binary Doppler shifts or solar-system delays, forcing binary and astrometric parameters to be fixed rather than fitted.
2. Unbinned likelihood methods: While these use every photon and allow all parameters to be fitted, they typically rely on Markov Chain Monte Carlo (MCMC) sampling of the full posterior. This becomes computationally inefficient as the number of parameters increases (e.g., for complex noise models) and struggles with "Neal's Funnel" problems when fitting low-amplitude stochastic noise. Furthermore, these methods often assume a fixed pulse profile template, ignoring uncertainties in the pulse shape.

Methodology
The authors present a new method that transforms the gamma-ray timing problem into a weighted least-squares problem solvable via Gibbs sampling. The core innovation involves introducing latent variables to handle the mixture-model nature of the gamma-ray likelihood.

1. Latent Variable Formulation: The likelihood function for gamma-ray data is a weighted sum of a constant background term and a mixture-model pulse profile. The authors introduce two latent vectors, $z$ and $m$, for each photon:

   • $z$ assigns each photon to a specific component of the pulse profile (one of $K$ Gaussian peaks, the unpulsed pulsar emission, or the background).
   • $m$ assigns integer "wraps" (rotations) to photons to resolve phase ambiguities relative to the timing model.
     By conditioning on these assignments, the complex mixture likelihood reduces to a product of simple Gaussian distributions.

2. Gibbs Sampling Scheme: Instead of sampling the high-dimensional joint posterior directly, the method iteratively samples from conditional distributions:

   • Template and Latent Sampling: Given the current timing model, the method samples the pulse profile parameters ($\tau$) and the latent assignments ($z, m$). The pulse profile parameters are sampled using a Metropolis-Hastings (MH) step, while assignments are drawn based on the probability weights and the current template.
   • Timing Model and Hyperparameter Sampling: Given the fixed template and assignments, the problem becomes a standard Gaussian process timing problem. The timing model parameters ($\theta$) and noise hyperparameters ($\lambda$) are sampled. The authors use a collapsed Gibbs approach, sampling $\lambda$ from the marginal likelihood first to avoid funnel-shaped posteriors, and then sampling $\theta$ from the resulting conditional Gaussian distribution.

3. Orbital Period Variations (OPVs): The method is specifically applied to model OPVs in "spider" binary systems (black widows and redbacks). These variations are modeled as a Gaussian process using a Fourier basis expansion (the "ORBWAVES" model in the PINT software), allowing for the fitting of complex, stochastic orbital phase deviations without the strong parameter correlations found in Taylor expansion models.

4. Implementation: The method is implemented in a Python package named shoogle, which interfaces with the PINT pulsar timing package to evaluate timing models and their derivatives.

Key Contributions

• Efficient Marginalization: The method enables the efficient marginalization over uncertain pulse profile shapes and latent photon assignments, which was previously computationally prohibitive in unbinned gamma-ray timing.
• Robust Noise Estimation: It allows for the robust estimation of stochastic noise processes (like TN and OPVs) within the unbinned framework, overcoming the inefficiencies of standard MCMC in high-dimensional spaces.
• Unified Framework: It provides a unified approach to fit timing models, noise hyperparameters, and pulse profile templates simultaneously, accounting for the non-Gaussian nature of gamma-ray data.

Results
The authors validate the method through three demonstrations:

1. Basic Timing (PSR J1526−2744): The Gibbs sampling results for a binary millisecond pulsar without complex noise were compared to an existing emcee-based method sampling the full posterior. The posterior distributions were found to be indistinguishable, confirming the method's accuracy in recovering the full posterior.
2. Timing Noise Simulations: The authors simulated 200 gamma-ray pulsars with injected power-law timing noise. Probability-Probability (PP) plots showed that the method recovers the amplitude and spectral index of the noise processes without intrinsic bias, even for low-amplitude signals.
3. PSR B1957+20 (Black Widow): The method was applied to the 12-year LAT data set of PSR B1957+20, a system with significant OPVs.
   • The model successfully fit the orbital period variations as a Gaussian process, finding a fractional amplitude of $\Delta P_{orb}/P_{orb} \sim 10^{-7}$ and a spectral index $\gamma \approx 5$.
   • The analysis yielded a 95% confidence upper limit on the GWB amplitude of $A_{GWB} < 4.6 \times 10^{-14}$ yr$^{3/2}$. This limit is consistent with previous results but obtained while simultaneously marginalizing over the pulse profile and fitting complex orbital variations, rather than fixing these parameters.

Significance and Claims
The paper claims that this new method significantly advances gamma-ray pulsar timing by enabling the analysis of faint pulsars and complex systems (like spider binaries) that were previously difficult to time accurately. By marginalizing over the pulse profile shape and using Gibbs sampling, the method provides more robust estimates of timing noise and orbital variations compared to previous unbinned approaches that assumed fixed templates.

The authors highlight that this method facilitates the construction of a Gamma-ray Pulsar Timing Array (GPTA) with improved sensitivity. Specifically, it allows for the inclusion of spider binaries (often excluded from radio PTAs due to eclipses) and enables joint radio-gamma-ray timing analyses. The paper notes that when combined with Fourier-domain likelihoods for inter-pulsar correlations, this approach provides a more robust estimate of the GWB amplitude, avoiding the underestimation tendencies of ToA-based methods. The method is presented as a practical tool for analyzing current and future Fermi-LAT data releases, particularly for characterizing the stochastic nature of orbital period variations in binary systems.