Logo
ScienceStack
Table of Contents
Fetching ...
Sign In

Timing Gamma-ray Pulsars using Gibbs Sampling

Colin J. Clark, Serena Valtolina, Lars Nieder, Rutger van Haasteren

TL;DR

The paper tackles gamma-ray pulsar timing with LAT data, where photon-by-photon measurements yield a non-Gaussian, multi-modal likelihood. It introduces a latent-variable Gibbs-sampling framework that integrates a Gaussian-process timing approach with marginalization over uncertain pulse-profile shapes, enabling robust estimation of timing parameters, timing-noise properties, and orbital-period variations. Implemented in the shoogle package, the method is validated on simulations and real data, producing unbiased TN parameter recovery and competitive GWB upper limits for PSR B1957+20. This approach enhances gamma-ray pulsar timing analyses and supports future joint radio-gamma timing efforts to improve stochastic background constraints and astrophysical inferences.

Abstract

Timing analyses of gamma-ray pulsars in the Fermi Large Area Telescope data set can provide sensitive probes of many astrophysical processes, including timing noise in young pulsars, orbital period variations in redback binaries, and the stochastic gravitational wave background (GWB). These goals can require careful accounting of stochastic noise processes, but existing methods developed to achieve this in radio pulsar timing analyses cannot be immediately applied to the discrete gamma-ray arrival time data. To address this, we have developed a new method for timing gamma-ray pulsars, in which the timing model fit is transformed into a weighted least squares problem by randomly assigning each photon to an individual Gaussian component of a template pulse profile. These random assignments are then numerically marginalised over through a Gibbs sampling scheme. This method allows for efficient estimation of timing and noise model parameters, while taking into account uncertainties in the pulse profile shape. We simulated Fermi-LAT data sets for gamma-ray pulsars with power-law timing noise processes, showing that this method provides robust estimates of timing noise parameters. We also describe a Gaussian-process model for orbital period variations in black-widow and redback binary systems that can be fit using this new timing method. We demonstrate this method on the black-widow binary millisecond pulsar B1957+20, where the orbital period varies significantly over the LAT data, but which provides one of the most stringent gamma-ray upper limits on the GWB.

Timing Gamma-ray Pulsars using Gibbs Sampling

TL;DR

The paper tackles gamma-ray pulsar timing with LAT data, where photon-by-photon measurements yield a non-Gaussian, multi-modal likelihood. It introduces a latent-variable Gibbs-sampling framework that integrates a Gaussian-process timing approach with marginalization over uncertain pulse-profile shapes, enabling robust estimation of timing parameters, timing-noise properties, and orbital-period variations. Implemented in the shoogle package, the method is validated on simulations and real data, producing unbiased TN parameter recovery and competitive GWB upper limits for PSR B1957+20. This approach enhances gamma-ray pulsar timing analyses and supports future joint radio-gamma timing efforts to improve stochastic background constraints and astrophysical inferences.

Abstract

Timing analyses of gamma-ray pulsars in the Fermi Large Area Telescope data set can provide sensitive probes of many astrophysical processes, including timing noise in young pulsars, orbital period variations in redback binaries, and the stochastic gravitational wave background (GWB). These goals can require careful accounting of stochastic noise processes, but existing methods developed to achieve this in radio pulsar timing analyses cannot be immediately applied to the discrete gamma-ray arrival time data. To address this, we have developed a new method for timing gamma-ray pulsars, in which the timing model fit is transformed into a weighted least squares problem by randomly assigning each photon to an individual Gaussian component of a template pulse profile. These random assignments are then numerically marginalised over through a Gibbs sampling scheme. This method allows for efficient estimation of timing and noise model parameters, while taking into account uncertainties in the pulse profile shape. We simulated Fermi-LAT data sets for gamma-ray pulsars with power-law timing noise processes, showing that this method provides robust estimates of timing noise parameters. We also describe a Gaussian-process model for orbital period variations in black-widow and redback binary systems that can be fit using this new timing method. We demonstrate this method on the black-widow binary millisecond pulsar B1957+20, where the orbital period varies significantly over the LAT data, but which provides one of the most stringent gamma-ray upper limits on the GWB.
Paper Structure (16 sections, 38 equations, 5 figures, 1 table)

This paper contains 16 sections, 38 equations, 5 figures, 1 table.

Table of Contents

  1. Introduction
  2. Methods
  3. Pulsar timing using Gaussian processes
  4. Gamma-ray pulsar timing
  5. Gibbs sampling
  6. Sampling the template pulse profile and latent parameters
  7. Sampling the timing model and hyperparameters
  8. Summary of Gibbs sampling procedure
  9. Modelling Orbital Period Variations
  10. Demonstrations
  11. PSR J1526$-$1744: basic timing model example
  12. Timing Noise Simulations
  13. PSR B1957$+$20: Orbital period variations and timing noise modelling
  14. Discussion
  15. Conclusions
  16. ...and 1 more sections

Figures (5)

  • Figure 1: Illustration of the Gibbs sampling procedure applied to the binary MSP J1526$-$2744. A: photon phases are determined by the current timing model parameters, with probability weights indicated by the greyscale. B: a short Metropolis-Hastings MCMC chain is produced over the template pulse profiles (faint black curves), given these photon phases, and the final step is chosen as the next sample for $\vec{\tau}$ (solid red curve). C: The relative likelihoods as a function of phase for each Gaussian component in the template pulse profile are then used, alongside the photon weights, to determine the probabilities for assigning each photon to a single peak (or to the background). D: These assignments are performed randomly according to these probabilities, with the outcomes illustrated by the colours of photons in panel corresponding to the relevant peak in the panel above (or the grayscale for photons assigned to the background). E: A short MH-MCMC chain, starting from the red square and following the blue path to the green triangle, is run targeting the marginal likelihood for the hyperparameters, illustrated by the greyscale and red contour lines at the 1$\sigma$, 2$\sigma$ and 3$\sigma$ levels. Here, a power-law timing noise model is assumed, with hyperparameters consisting of the log-amplitude $\log_{10}A_{\rm TN}$ and spectral index $\gamma_{\rm TN}$. The final sample in the chain is chosen as the next sample for $\vec{\lambda}$. F: A weighted least-squares fit is performed using Gaussian likelihoods for all photons assigned to peaks in the template, constrained by the prior distribution defined by the chosen sample of $\vec{\lambda}$. A single sample for the timing model (green curve) is chosen from the posterior distribution that results from this fit (illustrated by the faint black lines). Here, the sinusoidal curve whose amplitude grows with time indicates a significant detection of proper motion. The photon phases in panel A are then updated by these resulting phase shifts, and the process repeats.
  • Figure 2: Comparison between the posterior samples obtained via Gibbs sampling (red) and those obtained using the existing emcee method (black), for the binary millisecond pulsar PSR J1526$-$2744. Template pulse profile parameters are shown in the upper right corner, while timing model parameters are shown in the lower left. There are no strong correlations between pairs of parameters across these two blocks. Contour lines are at the 1$\sigma$, 2$\sigma$ and 3$\sigma$ level.
  • Figure 3: PP-plot for the recovery of the TN amplitude and slope, as defined in Eq. \ref{['e:broken_pwl']}, from 200 simulated pulsars. We also show the case of the amplitude recovery from datasets where the slope was fixed to $13/3$ (predicted value for a GWB signal). The gray areas are the 1$\sigma$ and 3$\sigma$ confidence intervals from the predicted distribution (black diagonal).
  • Figure 4: Results of Gibbs-sampling analysis on PSR B1957$+$20. Left panels show the photon phases (bottom) and integrated pulse profile (top) according to the best-fitting timing model. Faint black curves on the pulse profile plot show individual posterior samples for the template pulse profile, with the best-fitting model shown in orange. Middle panels illustrate the posterior uncertainty on the photon phases (bottom), and on the PSD of the GWB-like timing noise component that we searched for (top). The dashed black line shows the estimated white-noise level, while the red shaded regions show the $1\sigma$ and $2\sigma$ constraints on a putative power-law noise component with spectral index $\gamma = 13/3$. Black lines show the posterior uncertainties on the individual Fourier powers - no significant power is detected above the white noise level at any frequency, and so these posteriors all closely follow the priors. Right panels show the same, but for orbital phase noise, in which a steep spectrum process is detected.
  • Figure 5: Posterior distribution for the GWB and OPV hyperparameters for PSR B1957$+$20. Contours on the 2-D distributions are at the 1$\sigma$ and 2$\sigma$ levels, while dashed vertical lines on the 1-D marginal distributions indicate the 5%, 50% and 95% quantiles.