Timing gamma-ray pulsars using Gibbs sampling

Timing analyses of gamma-ray pulsars in the Fermi Large Area Telescope (LAT) dataset 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 stochastic noise processes to be taken carefully into account, 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 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 enables an efficient estimation of timing and noise model parameters, while taking into account uncertainties in the pulse profile shape. We simulated Fermi -LAT datasets for gamma-ray pulsars with power-law timing noise processes, and show 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.