Taming outliers in pulsar-timing datasets with hierarchical likelihoods and Hamiltonian sampling

Pulsar-timing datasets have been analyzed with great success using probabilistic treatments based on Gaussian distributions, with applications ranging from studies of neutron-star structure to tests of general relativity and searches for nanosecond gravitational waves. As for other applications of Gaussian distributions, outliers in timing measurements pose a significant challenge to statistical inference, since they can bias the estimation of timing and noise parameters, and affect reported parameter uncertainties. We describe and demonstrate a practical end-to-end approach to perform Bayesian inference of timing and noise parameters robustly in the presence of outliers, and to identify these probabilistically. The method is fully consistent (i.e., outlier-ness probabilities vary in tune with the posterior distributions of the timing and noise parameters), and it relies on the efficient sampling of the hierarchical form of the pulsar-timing likelihood. Such sampling has recently become possible with a "no-U-turn" Hamiltonian sampler coupled to a highly customized reparametrization of the likelihood; this code is described elsewhere, but it is already available online. We recommend our method as a standard step in the preparation of pulsar-timing-array datasets: even if statistical inference is not affected, follow-up studies of outlier candidates can reveal unseen problems in radio observations and timing measurements; furthermore, confidence in the results of gravitational-wave searches will only benefit from stringent statistical evidence that datasets are clean and outlier-free.Comment: 6 pages, 2 figures, RevTeX 4.

An Efficient Approximation to the Likelihood for Gravitational Wave Stochastic Background Detection Using Pulsar Timing Data

Direct detection of gravitational waves by pulsar timing arrays will become feasible over the next few years. In the low frequency regime ($10^{-7}$ Hz -- $10^{-9}$ Hz), we expect that a superposition of gravitational waves from many sources will manifest itself as an isotropic stochastic gravitational wave background. Currently, a number of techniques exist to detect such a signal; however, many detection methods are computationally challenging. Here we introduce an approximation to the full likelihood function for a pulsar timing array that results in computational savings proportional to the square of the number of pulsars in the array. Through a series of simulations we show that the approximate likelihood function reproduces results obtained from the full likelihood function. We further show, both analytically and through simulations, that, on average, this approximate likelihood function gives unbiased parameter estimates for astrophysically realistic stochastic background amplitudes.Comment: 10 pages, 3 figure

Analysis of the first IPTA Mock Data Challenge by the EPTA timing data analysis working group

This is a summary of the methods we used to analyse the first IPTA Mock Data Challenge (MDC), and the obtained results. We have used a Bayesian analysis in the time domain, accelerated using the recently developed ABC-method which consists of a form of lossy linear data compression. The TOAs were first processed with Tempo2, where the design matrix was extracted for use in a subsequent Bayesian analysis. We used different noise models to analyse the datasets: no red noise, red noise the same for all pulsars, and individual red noise per pulsar. We sampled from the likelihood with four different samplers: "emcee", "t-walk", "Metropolis-Hastings", and "pyMultiNest". All but emcee agreed on the final result, with emcee failing due to artefacts of the high-dimensionality of the problem. An interesting issue we ran into was that the prior of all the 36 (red) noise amplitudes strongly affects the results. A flat prior in the noise amplitude biases the inferred GWB amplitude, whereas a flat prior in log-amplitude seems to work well. This issue is only apparent when using a noise model with individually modelled red noise for all pulsars. Our results for the blind challenges are in good agreement with the injected values. For the GWB amplitudes we found h_c = 1.03 +/- 0.11 [10^{-14}], h_c = 5.70 +/- 0.35 [10^{-14}], and h_c = 6.91 +/- 1.72 [10^{-15}], and for the GWB spectral index we found gamma = 4.28 +/- 0.20, gamma = 4.35 +/- 0.09, and gamma = 3.75 +/- 0.40. We note that for closed challenge 3 there was quite some covariance between the signal and the red noise: if we constrain the GWB spectral index to the usual choice of gamma = 13/3, we obtain the estimates: h_c = 10.0 +/- 0.64 [10^{-15}], h_c = 56.3 +/- 2.42 [10^{-15}], and h_c = 4.83 +/- 0.50 [10^{-15}], with one-sided 2 sigma upper-limits of: h_c <= 10.98 [10^{-15}], h_c <= 60.29 [10^{-15}], and h_c <= 5.65 [10^{-15}].Comment: 10 pages, 5 figure

From bright binaries to bumpy backgrounds: Mapping realistic gravitational wave skies with pulsar-timing arrays

Within the next several years, pulsar-timing array programs will likely usher in the next era of gravitational-wave astronomy through the detection of a stochastic background of nanohertz-frequency gravitational waves, originating from a cosmological population of inspiraling supermassive binary black holes. While the source positions will likely be isotropic to a good approximation, the gravitational-wave angular power distribution will be anisotropic, with the most massive and/or nearby binaries producing signals that may resound above the background. We study such a realistic angular power distribution, developing fast and accurate sky-mapping strategies to localize pixels and extended regions of excess power while simultaneously modeling the background signal from the less massive and more distant ensemble. We find that power anisotropy will be challenging to discriminate from isotropy for realistic gravitational-wave skies, requiring SNR >10 in order to favor anisotropy with 10:1 posterior odds in our case study. Amongst our techniques, modeling the population signal with multiple point sources in addition to an isotropic background provides the most physically motivated and easily interpreted maps, while spherical-harmonic modeling of the square-root power distribution, P(Ω)^(1/2), performs best in discriminating from overall isotropy. Our techniques are modular and easily incorporated into existing pulsar-timing array analysis pipelines

The Need For Speed: Rapid Refitting Techniques for Bayesian Spectral Characterization of the Gravitational Wave Background Using PTAs

Current pulsar timing array (PTA) techniques for characterizing the spectrum of a nanohertz-frequency stochastic gravitational-wave background (SGWB) begin at the stage of timing data. This can be slow and memory intensive with computational scaling that will worsen PTA analysis times as more pulsars and observations are added. Given recent evidence for a common-spectrum process in PTA data sets and the need to understand present and future PTA capabilities to characterize the SGWB through large-scale simulations, we have developed efficient and rapid approaches that operate on intermediate SGWB analysis products. These methods refit SGWB spectral models to previously-computed Bayesian posterior estimations of the timing power spectra. We test our new methods on simulated PTA data sets and the NANOGrav $12.5$-year data set, where in the latter our refit posterior achieves a Hellinger distance from the current full production-level pipeline that is $\lesssim 0.1$. Our methods are $\sim10^2$--$10^4$ times faster than the production-level likelihood and scale sub-linearly as a PTA is expanded with new pulsars or observations. Our methods also demonstrate that SGWB spectral characterization in PTA data sets is driven by the longest-timed pulsars with the best-measured power spectral densities which is not necessarily the case for SGWB detection that is predicated on correlating many pulsars. Indeed, the common-process spectral properties found in the NANOGrav $12.5$-year data set are given by analyzing only the $\sim10$ longest-timed pulsars out of the full $45$ pulsar array, and we find that the ``shallowing'' of the common-process power-law model occurs when gravitational-wave frequencies higher than $\sim 50$~nanohertz are included. The implementation of our methods is openly available as a software suite to allow fast and flexible PTA SGWB spectral characterization and model selection.Comment: 19 pages, 12 figures. Submitting to Physical Review

Hyper-efficient model-independent Bayesian method for the analysis of pulsar timing data

A new model independent method is presented for the analysis of pulsar timing data and the estimation of the spectral properties of an isotropic gravitational wave background (GWB). We show that by rephrasing the likelihood we are able to eliminate the most costly aspects of computation normally associated with this type of data analysis. When applied to the International Pulsar Timing Array Mock Data Challenge data sets this results in speedups of approximately 2 to 3 orders of magnitude compared to established methods. We present three applications of the new likelihood. In the low signal to noise regime we sample directly from the power spectrum coefficients of the GWB signal realization. In the high signal to noise regime, where the data can support a large number of coefficients, we sample from the joint probability density of the power spectrum coefficients for the individual pulsars and the GWB signal realization. Critically in both these cases we need make no assumptions about the form of the power spectrum of the GWB, or the individual pulsars. Finally we present a method for characterizing the spatial correlation between pulsars on the sky, making no assumptions about the form of that correlation, and therefore providing the only truly general Bayesian method of confirming a GWB detection from pulsar timing data.Comment: 9 pages, 4 figure