Towards robust gravitational wave detection with pulsar timing arrays

Precision timing of highly stable millisecond pulsars is a promising technique for the detection of very low frequency sources of gravitational waves. In any single pulsar, a stochastic gravitational wave signal appears as an additional source of timing noise that can be absorbed by the noise model, and so it is only by considering the coherent response across a network of pulsars that the signal can be distinguished from other sources of noise. In the limit where there are many gravitational wave sources in the sky, or many pulsars in the array, the signals produce a unique tensor correlation pattern that depends only on the angular separation between each pulsar pair. It is this distinct fingerprint that is used to search for gravitational waves using pulsar timing arrays. Here we consider how the prospects for detection are diminished when the statistical isotropy of the timing array or the gravitational wave signal is broken by having a finite number of pulsars and a finite number of sources. We find the standard tensor-correlation analysis to be remarkably robust, with a mild impact on detectability compared to the isotropic limit. Only when there are very few sources and very few pulsars does the standard analysis begin to fail. Having established that the tensor correlations are a robust signature for detection, we study the use of “sky scrambles” to break the correlations as a way to increase confidence in a detection. This approach is analogous to the use of “time slides” in the analysis of data from ground-based interferometric detectors.

Figure 1

The upper two maps show the distribution of the gravitational wave power across the sky, scaled by the all-sky average for two signal model. The map on the upper left is for a simulated black hole population, smoothed with a two degree Gaussian blur and clipped at a contrast of 100 to enhance weaker features. The peak intensity in the un-clipped map exceeds 800. The map on the upper right is for a statistically isotropic signal with the same power spectrum as the black hole simulation. The map has been smoothed with a two degree Gaussian blur. Clearly, the power distribution from a realistic black hole population is far from isotropic. The lower two panels show the detected power in the Earth term for pulsars at different sky locations, in other words, the raw signals convolved with the antenna patterns, summed and squared. Despite the large differences in the underlying power distribution, the response to the anisotropic BH background (lower left) is qualitatively identical to the response to the statistically isotropic signal (lower right).

Figure 2

The detected power for pulsars at different sky locations for a single BH located at the center of the “petals”. The orientation of the petals rotates depending on the polarization of the signal. The broad spread in the power distribution for PTAs explains why the response to isotropic and anisotropic signals is so similar.

Figure 3

The average power spectrum in a 20-pulsar array for a black hole population (in red), a stochastic, isotropic model with the same average power level as the black hole population (in blue), and a power law, stochastic, isotropic model with amplitude $A={10}^{-15}$.

Figure 4

The detectability of the full simulated black hole population (red) and a statistically isotropic stochastic signal with the same power spectrum (blue) as a function of the white timing noise level in a simulated pulsar timing array with 20 equally sensitive pulsars, as measured by the signal-to-noise Bayes factor. The spread in the Bayes factors is computed by considering ten realizations of each signal type. The shaded bands cover a one standard deviation spread about the mean.

Figure 5

The detectability of down-sampled black hole populations (red) and a statistically isotropic stochastic signal with the same power spectrum (blue) as a function of the white timing noise level in a simulated pulsar timing array with 20 equally sensitive pulsars, as measured by the signal-to-noise Bayes factor. The upper panel is for the full population down-sampled by a factor of 10, the middle panel is down-sampled by 100, and the lower panel is for the ten brightest binaries. The spread in the Bayes factors is computed by considering ten realizations of each signal type. The shaded bands covers a one standard deviation spread about the mean.

Figure 6

The detectability of signals from a population of black holes (red) and a statistically isotropic stochastic signal with the same power spectrum (blue) for a small pulsar timing array made up of five equally sensitive pulsars, as measured by the signal-to-noise Bayes factor. The upper panel is for the full black hole population, while the lower panel is for the brightest ten black holes. The spread in the Bayes factors is computed by considering ten realizations of each signal type. The shaded bands covers a one standard deviation spread about the mean.

Figure 7

The detectability of signals from a population of black holes (red) and a statistically isotropic stochastic signal with the same power spectrum (blue) for the IPTA pulsar array, again for the full spectrum in the upper panel and the brightest 10 black holes in the lower panel. The horizontal axis shows the percentage of the actual timing noise in each pulsar that is included in the injected data. The results are more similar to the 5-pulsar array than the 20-pulsar case, despite the fact that there are 36 pulsars in the IPTA. This is because a few of the pulsars are much better timed than the others.

Figure 8

The detectability of simulated black hole populations as a function of the noise level for a simulated pulsar timing array with 20 equally sensitive pulsars. The red bands show the spread in log Bayes factors for the full tensor correlation model computed from multiple realizations of a black hole population model. The green bands show the spread in log Bayes factors for the diagonal correlation model applied to the same set of simulated signals. The gravitational wave signal model is favored when the log Bayes factor for the tensor correlation model versus the noise model is positive and that it exceeds the log Bayes factor for the diagonal correlation model versus the noise model. From top to bottom the simulations are for the full black hole background, the populations down-sampled by 10, then 100, and finally for just the ten brightest systems.

Figure 9

The log Bayes factors between the full and diagonal correlation model. Values greater than zero indicate that the Hellings-Downs correlation pattern has been detected. The red bands are for simulated black hole populations, and the blue bands are for the isotropic equivalents. For the highly down-sampled cases (lower panels), the detectability of the Hellings-Downs pattern is much more realization-dependent for the black hole populations than for the isotropic backgrounds, indicating a greater difficulty in differentiating between a GW signal and a different correlated noise source.

Figure 10

Histograms of the full match $M$ (in blue), and the off-diagonal match $\overline{M}$, for randomly drawn sky scrambles of an equal-sensitivity 20-pulsar timing array.

Figure 11

Log Bayes factors in favor of a detection as a function of the match $\overline{M}$ for an equal-sensitivity 20-pulsar timing array. The green points compare the full signal model to the noise model, while the purple points compare the full signal model to a signal model with a diagonal correlation matrix.

Figure 12

A histogram of the log Bayes factors between the full correlation model and a diagonal correlation model for the 18 independent sky scrambles for an equal-sensitivity 20-pulsar timing array. Here independent is defined as matches with the true correlation pattern and each other of $\overline{M}<0$. For reference, the analysis with unscrambled sky locations gave a log Bayes factor of 12.