Rapid inference for individual binaries and a stochastic background with pulsar timing array data

TL;DR

The paper tackles the problem of jointly inferring a stochastic gravitational wave background and deterministic continuous-wave signals in pulsar timing array data. It extends a Fourier-basis framework to include deterministic CW signals while modeling inter-pulsar correlations and red noise with a fixed white-noise model, enabling per-pulsar likelihoods and efficient, scalable likelihood evaluations. In simulations, the approach recovers injected GWB and CW parameters, matches standard analyses for background parameters, and delivers substantial computational speed-ups; this makes joint detection and characterization feasible for larger PTA datasets. The method thus provides a practical path toward resolving anisotropy and individual binaries in growing PTA data, with potential extensions to trans-dimensional and more complete noise modeling.

Abstract

The analysis of pulsar timing array data has provided evidence for a gravitational wave background in the nanohertz band. This raises the question of what is the source of the signal, is it astrophysical or cosmological in origin? If the signal originates from a population of supermassive black hole binaries, as is generally assumed, we can expect to see evidence for both anisotropy and to be able to resolve signals from individual binaries as more data are collected. The anisotropy and resolvable systems are caused by a small number of loud signals that stand out from the crowd. Here we focus on the joint detection of individual signals and a stochastic background. While methods have previously been developed to perform such an analysis, they are currently held back by the cost of computing the joint likelihood function. Each individual source is described by $N=8+2N_p$ parameters, where $N_p$ are the number of pulsars in the array. With the latest combined datasets having over one hundred pulsars, the parameter space is very large, and consequently, it takes a large number of likelihood evaluations to explore these models. Here we present a new approach that extends the Fourier basis method, previously introduced to accelerate analyses for stochastic signals, to also include deterministic signals. Key elements of the method are that the likelihood evaluations are per-pulsar, avoiding expensive operations on large matrices, and the templates for individual binaries can be computed analytically or using fast Fourier methods on a sparsely sampled grid of time samples. This analysis method scales better than quadratically with the size of the dataset, while the approach currently being used in most analyses scales quartically or worse with the number of data points. As datasets grow with more observations, this analysis will be orders of magnitude faster than previous approaches.

Figures (11)

• Figure 1: Samples of selected parameters for the stochastic processes obtained via HMC. The first two columns correspond to the power law hyper-parameters of the GWB, and subsequent columns are power law hyper-parameters for intrinsic pulsar RN and the Fourier coefficients. The "$a$" Fourier coefficients scale sine functions and "$b$" scale cosine. Parameter subscripts index the pulsars in the array, superscripts on the Fourier coefficients index the corresponding frequency bins. The dashed lines are the parameter values injected into the simulated data. The summary statistics and shading of one-dimensional marginal distributions is $1\sigma$ either side the median computed with the CDF. Contours of the two-dimensional distributions enclose the $1\sigma$ and $2\sigma$ credible regions.
• Figure 2: Samples of continuous wave model parameters obtained via HMC. The CW SNR is 13.7. The dashed lines are the parameter values injected into the simulated dataset. The summary statistics and shading of one-dimensional marginal distributions is $1\sigma$ either side the median computed with the CDF. Contours of the two-dimensional distributions enclose the $1\sigma$ and $2\sigma$ credible regions.
• Figure 3: Samples of selected pulsar distance and phase parameters obtained via HMC. Parameter subscripts index pulsars in the array. The dashed lines are the parameter values injected into the simulated dataset. The summary statistics and shading of one-dimensional marginal distributions is $1\sigma$ either side the median computed with the CDF. Contours of the two-dimensional distributions enclose the $1\sigma$ and $2\sigma$ credible regions.
• Figure 4: Sky map illustrating resolvability of the CW source's sky location. The red stars are pulsars observed in the PTA. The blue circle is the CW source. The color of the sky map is the log-posterior density (numerically) marginalized over all parameters except sky location.
• Figure 5: Violin plot illustrating the recovered power for the stochastic GWB from a free spectral run. The injected GWB power law is shown in orange. The dark (light) green shading corresponds to a $1\sigma$ ($2\sigma$) uncertainty bound in the hyper-parameters from a hierarchical power law analysis on the same dataset. The vertical red line is the frequency of the injected CW signal.