Publications

Fermi LAT observations of gamma-ray pulsars can be used to build a pulsar timing array (PTA) experiment to search for gravitational wave (GW) signals at nanohertz frequencies. At those frequencies, the dominant signal is expected to be a stochastic gravitational wave background (GWB) produced by the incoherent superposition of the quasi-monochromatic GW emissions from a population of supermassive black hole binaries. While the radio PTAs have recently announced compelling evidence for a GWB signal, in 2022 an analysis of 12.5 years of Fermi data for 35 pulsars led to an upper limit of 1e-14 for the GWB amplitude. Here, we reanalyse the same dataset using a regularized likelihood method that correctly models cross-pulsar correlations directly from the photons, while additionally marginalising over the uncertain pulse profile shape. We prove through simulations that the photon-by-photon method for GWB recoveries is, statistically, more robust. The resulting upper limit obtained for the GWB strain amplitude is 1.2e-14, consistent with previous analyses.

Pulsar Timing Array (PTA) projects have reported various lines of evidence suggesting the presence of a stochastic gravitational wave (GW) background in their data. One key line of evidence involves a detection statistic sensitive to interpulsar correlations, such as those induced by GWs. A p-value is then calculated to assess how unlikely it is for the observed signal to arise under the null hypothesis H_0, purely by chance. However, PTAs cannot empirically draw samples from H_0. As a workaround, various techniques are used in the literature to approximate p-values under H_0. One such technique, which has been characterized as a model-independent method, is the use of “scrambling” transformations that modify the data to cancel out pulsar correlations, thereby simulating realizations from H_0. In this work, scrambling methods and the detection statistic are investigated from first principles. The p-value methodology that is discussed is general, but the discussions regarding a specific detection statistic apply to the detection of a stochastic background of gravitational waves with PTAs. All methods in the literature to calculate p-values for such a detection statistic are rigorously analyzed, and analytical expressions are derived for the distribution of the detection statistic and the corresponding p-values. All this leads to the conclusion that scrambling methods are not model independent and thus not completely empirical. With a single realization of data our results are necessarily always model dependent, which any analysis will need to accept. Instead of scrambling approaches, rigorous Bayesian and Frequentist p-value calculation methods are advocated, the evaluation of which depend on the generalized chi^2 distribution. This view is consistent with the posterior predictive p-value approach that is already in the literature. Efficient expressions are derived to evaluate the generalized chi^2 distribution of the detection statistic on real data. It is highlighted that no Frequentist p-values have been calculated correctly in the PTA literature to date.

Pulsar timing arrays (PTAs) seek to detect a nano-Hz stochastic gravitational-wave background (GWB) by searching for the characteristic Hellings and Downs angular pattern of timing residual correlations. So far, the evidence remains below the conventional 5-sigma threshold, as assessed using the literature-standard “optimal cross-correlation detection statistic”. While this quadratic combination of cross-correlated data maximizes the deflection (signal-to-noise ratio), it does not maximize the detection probability at fixed false-alarm probability (FAP), and therefore is not Neyman-Pearson (NP) optimal for the assumed noise and signal models. The NP-optimal detection statistic is a different quadratic form, but is not used because it also incorporates autocorrelations, making it more susceptible to uncertainties in the modeling of pulsar timing noise. Here, we derive the best compromise: a quadratic detection statistic which is as close as possible to the NP-optimal detection statistic (minimizing the variance of its difference with the NP statistic) subject to the constraint that it only uses cross-correlations, so that it is less affected by pulsar noise modeling errors.

The recent announcement of evidence for a stochastic background of gravitational waves (GWB) in pulsar timing array (PTA) data has piqued interest across the scientific community. A combined analysis of all currently available data holds the promise of confirming the announced evidence as a solid detection of a GWB. However, the complexity of individual pulsar noise models and the variety of modeling tools used for different types of pulsars present significant challenges for a truly unified analysis. In this work we propose a novel approach to the analysis of PTA data: first a posterior distribution over Fourier modes is produced for each pulsar individually. Then, in a global analysis of all pulsars these posterior distributions can be re-used for a GWB search, which retains all information regarding the signals of interest without the added complexity of the underlying noise models or implementation differences. This approach facilitates combining radio and gamma-ray pulsar data, while reducing the complexity of the model and of its implementations when carrying out a GWB search with PTA data.

Pulsar Timing Array (PTA) searches for nHz gravitational-wave backgrounds (GWBs) typically model time-correlated noise by assuming a diagonal covariance in Fourier space, neglecting inter-frequency correlations introduced by the finite observation window. We show that this diagonal approximation can lead to biased estimates of spectral parameters, especially for the common red process that represents the GWB. To address these limitations, we present a method that (i) computes the time-domain autocorrelation on a coarse grid using a fast Fourier transform (FFT), (ii) interpolates it accurately to the unevenly sampled observation times, and (iii) incorporates it into a low-rank likelihood via the Sherman–Morrison–Woodbury identity. Using both analytic covariance comparisons and end-to-end simulations inspired by the NANOGrav 15-year dataset, we demonstrate that our method captures frequency correlations faithfully, avoids Gibbs ringing, and recovers unbiased spectral parameters with modest computational cost. As PTA datasets increase in sensitivity and complexity, our approach offers a practical and scalable path to fully accurate covariance modeling for current and future analyses.

Over the past decade and a half, adoption of Bayesian inference in pulsar timing analysis has led to increasingly sophisticated models. The recent announcement of evidence for a stochastic background of gravitational waves by various pulsar timing array (PTA) projects highlighted Bayesian inference as a central tool for parameter estimation and model selection. Despite its success, Bayesian inference is occasionally misused in the pulsar timing community. A common workflow is that the data is analysed in multiple steps: a first analysis of single pulsars individually, and a subsequent analysis of the whole array of pulsars. A mistake that is then sometimes introduced stems from using the posterior distribution to craft the prior for the analysis of the same data in a second step, a practice referred to in the statistics literature as ‘circular analysis’. This is done to prune the model for computational efficiency. Multiple recent high-profile searches for gravitational waves by PTA projects have this workflow. This letter highlights this error and suggests that Spike and Slab priors can be used to carry out model averaging instead of model selection in a single pass. Spike and Slab priors are proved to be equal to log-uniform priors.

Pulsar timing array (PTA) projects have found evidence of a stochastic background of gravitational waves (GWB) using data from an ensemble of pulsars. In the literature, minimal assumptions are made about the signal and noise processes that affect data from these pulsars, such as pulsar spin noise. These assumptions are encoded as uninformative priors in Bayesian searches, though frequentist approaches make similar assumptions. Uninformative priors are not suitable for (noise) properties of pulsars in an ensemble, and they bias estimates of model parameters such as gravitational-wave signal parameters. Both frequentist and Bayesian searches are affected. In this article, more appropriate priors are proposed in the language of hierarchical Bayesian modeling, where the properties of the ensemble of pulsars are jointly described with the properties of the individual components of the ensemble. Results by PTA projects should be reevaluated using hierarchical models.

Pulsar-timing data sets have been analysed 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 data sets: 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 data sets are clean and outlier-free.

Many data-analysis problems involve large dense matrices that describe the covariance of wide-sense stationary noise processes; the computational cost of inverting these matrices, or equivalently of solving linear systems that contain them, is often a practical limit for the analysis. We describe two general, practical, and accurate methods to approximate stationary covariance matrices as low-rank matrix products featuring carefully chosen spectral components. These methods can be used to greatly accelerate data-analysis methods in many contexts, such as the Bayesian and generalized-least-squares analysis of pulsar-timing residuals.

In this work we review the application of the theory of Gaussian processes to the modeling of noise in pulsar-timing data analysis, and we derive various useful and optimized representations for the likelihood expressions that are needed in Bayesian inference on pulsar-timing-array data sets. The resulting viewpoint and formalism lead us to two improved parameter-sampling schemes inspired by Gibbs sampling. The new schemes have vastly lower chain autocorrelation lengths than the Markov-chain Monte Carlo methods currently used in pulsar-timing data analysis, potentially speeding up Bayesian inference by orders of magnitude. The new schemes can be used for a full-noise-model analysis of the large data sets currently being assembled by pulsar-timing-array collaborations, which generally present a serious computational challenge to existing methods.

We describe a new method for extracting gravitational wave signals from pulsar timing data. We show that any gravitational wave signal can be decomposed into an orthogonal set of sky maps, with the number of maps equal to the number of pulsars in the timing array. These maps may be used as a basis to construct gravitational wave templates for any type of source, including collections of point sources. A variant of the standard Hellings-Downs correlation analysis is recovered for statistically isotropic signals. The template based approach allows us to probe potential anisotropies in the signal and produce maps of the gravitational wave sky.

The analysis of pulsar timing data, especially in pulsar timing array (PTA) projects, has encountered practical difficulties: evaluating the likelihood and/or correlation-based statistics can become prohibitively computationally expensive for large datasets. In situations where a stochastic signal of interest has a power spectral density that dominates the noise in a limited bandwidth of the total frequency domain, a linear transformation exists that transforms the timing residuals to a basis in which virtually all the information about the stochastic signal of interest is contained in a small fraction of basis vectors. By only considering such a small subset of these “generalised residuals”, the dimensionality of the data analysis problem is greatly reduced, which can cause a large speedup in the evaluation of the likelihood: the ABC-method (Acceleration By Compression). Both direct tests on the likelihood, and Bayesian analysis of mock data, show that the signal can be recovered as well as with an analysis of uncompressed data.

Although it is widely understood that pulsar timing observations generally contain time-correlated stochastic signals (TCSSs; red timing noise is of this type), most data analysis techniques that have been developed make an assumption that the stochastic uncertainties in the data are uncorrelated, i.e. “white”. Recent work has pointed out that this can introduce severe bias in determination of timing-model parameters, and that better analysis methods should be used. This paper presents a detailed investigation of timing-model fitting in the presence of TCSSs, and gives closed expressions for the post-fit signals in the data. This results in a Bayesian technique to obtain timing-model parameter estimates in the presence of TCSSs, as well as computationally more efficient expressions of their marginalised posterior distribution. A new method to analyse hundreds of mock dataset realisations simultaneously without significant computational overhead is presented, as well as a statistically rigorous method to check the internal consistency of the results.

Direct detection of low-frequency gravitational waves (GWs, Hz) is the main goal of pulsar timing array (PTA) projects. One of the main targets for the PTAs is to measure the stochastic background of gravitational waves (GWB) whose characteristic strain is expected to approximately follow a power-law of the form , where f is the GW frequency. In this paper we use the current data from the European PTA to determine an upper limit on the GWB amplitude A as a function of the unknown spectral slope alpha with a Bayesian algorithm, by modelling the GWB as a random Gaussian process. For the case alpha=-2/3, which is expected if the GWB is produced by supermassive black hole binaries, we obtain a 95 per cent confidence upper limit on A of 6 × 10-15, which is 1.8 times lower than the 95 per cent confidence GWB limit obtained by the Parkes PTA in 2006. Our approach to the data analysis incorporates the multitelescope nature of the European PTA and thus can serve as a useful template for future intercontinental PTA collaborations.

Pulsar timing arrays (PTAs) are designed to detect gravitational waves with periods from several months to several years, e.g. those produced by wide supermassive black hole binaries in the centres of distant galaxies. Here, we show that PTAs are also sensitive to mergers of supermassive black holes. While these mergers occur on a time-scale too short to be resolvable by a PTA, they generate a change of metric due to non-linear gravitational-wave memory which persists for the duration of the experiment and could be detected. We develop the theory of the single-source detection by PTAs, and derive the sensitivity of PTAs to the gravitational-wave memory jumps. We show that mergers of 108 Modot black holes are 2 -sigma-detectable (in a direction-, polarization- and time-dependent way) out to comoving distances of sim1 billion light-years. Modern prediction for black hole merger rates imply marginal to modest chance of an individual jump detection by currently developed PTAs. The sensitivity is expected to be somewhat higher for futuristic PTA experiments with the Square Kilometre Array.

The long-term precise timing of Galactic millisecond pulsars holds great promise for measuring the long-period (months to years) astrophysical gravitational waves. Several gravitational-wave observational programs, called Pulsar Timing Arrays (PTA), are being pursued around the world.Here, we develop a Bayesian algorithm for measuring the stochastic gravitational-wave background (GWB) from the PTA data. Our algorithm has several strengths: (i) it analyses the data without any loss of information; (ii) it trivially removes systematic errors of known functional form, including quadratic pulsar spin-down, annual modulations and jumps due to a change of equipment; (iii) it measures simultaneously both the amplitude and the slope of the GWB spectrum and (iv) it can deal with unevenly sampled data and coloured pulsar noise spectra. We sample the likelihood function using Markov Chain Monte Carlo simulations. We extensively test our approach on mock PTA data sets and find that the algorithm has significant benefits over currently proposed counterparts. We show the importance of characterizing all red noise components in pulsar timing noise by demonstrating that the presence of a red component would significantly hinder the detection of the GWB.Lastly, we explore the dependence of the signal-to-noise ratio on the duration of the experiment, number of monitored pulsars and the magnitude of the pulsar timing noise. These parameter studies will help formulate observing strategies for the PTA experiments.