A Bayesian parameter estimation approach to pulsar time-of-arrival analysis

The increasing sensitivities of pulsar timing arrays to ultra-low frequency (nHz) gravitational waves promises to achieve direct gravitational wave detection within the next 5-10 years. While there are many parallel efforts being made in the improvement of telescope sensitivity, the detection of stable millisecond pulsars and the improvement of the timing software, there are reasons to believe that the methods used to accurately determine the time-of-arrival (TOA) of pulses from radio pulsars can be improved upon. More specifically, the determination of the uncertainties on these TOAs, which strongly affect the ability to detect GWs through pulsar timing, may be unreliable. We propose two Bayesian methods for the generation of pulsar TOAs starting from pulsar “search-mode” data and pre-folded data. These methods are applied to simulated toy-model examples and in this initial work we focus on the issue of uncertainties in the folding period. The final results of our analysis are expressed in the form of posterior probability distributions on the signal parameters (including the TOA) from a single observation.

💡 Research Summary

The paper presents a Bayesian framework for estimating pulsar times‑of‑arrival (TOAs) and their uncertainties directly from raw observational data, addressing a key limitation in current pulsar timing array (PTA) analyses. Two distinct data products are considered: (1) “search‑mode” data, which consist of high‑time‑resolution intensity measurements across many frequency channels, and (2) pre‑folded data, which have already been summed over many pulse periods using an assumed period.

The authors construct a simple yet illustrative signal model. Each pulse is represented by a Gaussian profile with amplitude A, width w, and a fixed period P. Random pulse‑to‑pulse jitter ξα with variance σ²ξ is added, and dispersion is modeled by a frequency‑dependent delay Δt(fk) proportional to the dispersion measure DM. The data are assumed to be corrupted by independent Gaussian noise of known variance.

For search‑mode data the time series in each frequency channel is Fourier transformed. In the frequency domain the signal appears as a series of narrow harmonics whose complex amplitudes factor into (i) a real envelope proportional to A·w·√(2π) multiplied by an exponential decay exp