Bayesian search for gravitational wave bursts in pulsar timing array data

The detection of nanohertz gravitational waves (GWs) is the main objective of pulsar timing arrays (PTAs). These large scale experiments regularly monitor a collection of millisecond pulsars to achieve this goal (see e.g. [1]). The three major PTAs currently operating are the North American Nanohertz Observatory for Gravitational Waves (NANOGrav, [2]), the European Pulsar Timing Array (EPTA, [3]), and the Parkes Pulsar Timing Array (PPTA, [4]). In addition, there are emerging PTA efforts in India (InPTA, [5]), China (CPTA, [6]), and South Africa (SAPTA, [7]). These project are collaborating under the International Pulsar Timing Array (IPTA, [8–10]) consortium.

The most promising GW sources in the nHz regime are supermassive black hole binaries (SMBHBs). These can be detected through observing their collective effect as a red noise process with quadrupolar correlations between pulsars (see e.g. [11]). The most massive and nearby of these SMBHBs might also be detectable individually (see e.g. [12]). It has also been proposed that these two searches can be carried out simultaneously, thus accounting for any interaction between the two types of signatures [13]. Searches have also been carried out for bursts with memory [14], i.e. a permanent deformation of spacetime after a violent astrophysical event, like the merger of two black holes.

In this paper, we focus on searching for generic GW transients (a.k.a. GW bursts) in PTA data. The subject has been studied extensively in the context of ground-based interferometric GW detectors, where many different algorithms are in use [15–20] and several searches have been carried out throughout different observing runs of the LIGO [21] and Virgo [22] detectors [23, 24]. The detection problem of GW bursts has also been considered for space-based detectors [25].

In the PTA context, no analysis of real data has been carried out so far. However, several methods have been suggested, including an analytical Bayesian framework [26], a Bayesian nonparametric approach [27], and a frequentist search working in the time–frequency domain [28]. Reference [29] suggests a possible improvement by considering the coherence between the pulsar terms, which appears in some cases. In this paper we present BayesHopperBurst ¹ , a Bayesian search for GW bursts in PTA data, which is based on a trans-dimensional reversible jump Markov chain Monte Carlo (RJMCMC) sampler [30, 31] akin to that used in the BayesWave algorithm [15] to search for and reconstruct GW bursts in the data of ground-based GW detectors. The work described in this paper is the continuation of that presented in [32], where the authors used similar methods to model noise transients in PTA data. One promising source of nHz GW bursts are parabolic (or highly eccentric) encounters of two SMBHs, which can occur in a hierarchical triple system of SMBHs [33]. Cosmic string kinks and cusps can also produce GW busts in this frequency range [34].

The paper is organized as follows. In section 2 we describe the methods employed by BayesHopperBurst including the model used to describe transient signals and the sampling techniques that enables it to effectively explore the parameter space. In section 3 we perform various injection tests, demonstrating that BayesHopperBurst recovers a wide variety of signals and noise transients. In section 4 we analyze a simulated dataset similar to the latest NANOGrav data release [2] to make a prediction of what we can expect from analyzing the real dataset. In section 5 we analyze B1855 + 09 from the NANOGrav 9 year dataset [35] to demonstrate our algorithm's noise transient modeling capabilities on real data. Finally, we offer concluding remarks in section 6.

In this section, we introduce the model used to describe signals and noise transients with arbitrary morphology, and we provide some details of the sampling techniques employed to efficiently explore the resulting high-dimensional parameter space. We use the enterprise ² [36] software package for handling PTA data and calculating the likelihood used in our Bayesian analysis. Timing models were produced by libstempo ³ , which is a python wrapper for the tempo2 ⁴ timing package. We also make use of the la_forge ⁵ [37] package for some of our figures.

2.1. Model

The model we use to describe our dataset has many similarities with the one used in the BayesHopper algorithm. The only fundamental difference is that the collection of sinusoids are replaced by a set of 'signal' wavelets coherent across detectors and a set of 'noise transient' wavelets which only appear in a given pulsar and describe transient noise features in the data. In this section we give a brief description of the model and focus on differences between BayesHopper and BayesHopperBurst. More details can be found in [13].

Consider a PTA consisting of N pulsars. The ith residual in the kth pulsar of the array is modeled as:

$$\begin{equation}\delta {t}_{ki}=\sum\limits _{l}{M}_{kil}\delta {\xi }_{kl}+{n}_{ki}+{g}_{ki}+{w}_{ki}\left({\boldsymbol{\theta }}_{\text{s}}\right)+{v}_{ki}\left({\boldsymbol{\theta }}_{\text{n}}\right),\end{equation} \tag{ 1 }$$

where δξ_kl is the offset from the best fit value of the lth timing model parameter, and M_kil is the design matrix, which represents a timing model linearized around the best fit parameter values. n_ki represents all noise processes unique to each pulsar, while g_ki is the contribution of an isotropic stochastic GW Background (GWB) which is correlated between different pulsars. The contribution of the transient GW signal is represented by w_ki ( θ _s), while that of incoherent transient noise features is described by v_ki ( θ _n), where θ _s and θ _n represent the parameters describing the GW signal and the transient noise respectively.

2.1.1. Transient noise model

The contribution of noise transients to the timing residual (v_ki ( θ _n)) is modeled as a sum of sine-Gaussian (Morlet–Gabor) wavelets:

$$\begin{equation}{v}_{ki}\left({\boldsymbol{\theta }}_{\text{n}}\right)=\sum\limits _{j=1}^{{M}_{k}}{\Psi}\left({t}_{ki};{\boldsymbol{\lambda }}_{kj}\right),\end{equation} \tag{ 2 }$$

where M_k is the number of wavelets used in the kth pulsar, t_ki is the ith observing time of the kth pulsar, and λ _kj is the parameter vector describing the jth wavelet in the kth pulsar, which is related to the full parameter vector as ${\boldsymbol{\theta }}_{\text{n}}=\left({M}_{1},{\boldsymbol{\lambda }}_{11},\dots ,{\boldsymbol{\lambda }}_{1{M}_{1}};\dots ;{M}_{N},{\boldsymbol{\lambda }}_{N1},\dots ,{\boldsymbol{\lambda }}_{N{M}_{N}}\right)$, and the wavelets are defined as:

$$\begin{equation}{\Psi}\left({t}_{i};{\boldsymbol{\lambda }}_{kj}\right)=A\enspace {\mathrm{e}}^{{\left(t-{t}_{0}\right)}^{2}/{\tau }^{2}}\enspace \mathrm{cos}\left(2\pi {f}_{0}\left(t-{t}_{0}\right)+{\phi }_{0}\right),\end{equation} \tag{ 3 }$$

where A is the amplitude, t₀ is the central time, τ is the characteristic duration, f₀ is the central time, and ϕ₀ is the initial phase. We can see that a single wavelet can be described by five parameters, i.e. λ _kj = (A, t₀, τ, f₀, ϕ₀). We use Morlet–Gabor wavelets as they have the compelling property of having the smallest time–frequency area allowed by the Heisenberg–Gabor limit [38]. Other functions, such as shapelets [39], might be better suited to certain types of signals, so we plan to investigate using them in the future.

2.1.2. Signal model

In order to impose the proper coherence between pulsars, in the signal model we model the GW waveform itself as a sum of wavelets, and then project that onto the pulsar lines of sight taking into account the corresponding antenna factors. Since PTAs are sensitive to the time integral of the metric perturbation, we introduce:

$$\begin{equation}{H}_{+,{\times}}\left(t\right)=\int {h}_{+,{\times}}\left(t\right)\mathrm{d}t,\end{equation} \tag{ 4 }$$

where h_+,× are the usual components of the metric perturbation corresponding to plus and cross polarized GWs. We model H_+,×(t) as:

$$\begin{equation}{H}_{+}\left(t\right)=\sum\limits _{j=1}^{N}{\Psi}\left(t;{t}_{0,j},{f}_{0,j},{\tau }_{j},{A}_{+,j},{\phi }_{0,+,j}\right),\end{equation} \tag{ 5 }$$

$$\begin{equation}{H}_{{\times}}\left(t\right)=\sum\limits _{j=1}^{N}{\Psi}\left(t;{t}_{0,j},{f}_{0,j},{\tau }_{j},{A}_{{\times},j},{\phi }_{0,{\times},j}\right),\end{equation} \tag{ 6 }$$

so that the two polarizations have independent amplitude and initial phase, but t₀, f₀, and τ are common parameters. This simplification will not affect the generality of our model, but helps in efficiency by reducing the number of parameters. Then a rotation around the propagation direction is allowed, which introduces the transformation:

$$\begin{equation}{\bar{H}}_{+}\left(t\right)={H}_{+}\left(t\right)\mathrm{cos}\left(2\psi \right)-{H}_{{\times}}\left(t\right)\mathrm{sin}\left(2\psi \right),\end{equation} \tag{ 7 }$$

$$\begin{equation}{\bar{H}}_{{\times}}\left(t\right)={H}_{+}\left(t\right)\mathrm{sin}\left(2\psi \right)+{H}_{{\times}}\left(t\right)\mathrm{cos}\left(2\psi \right),\end{equation} \tag{ 8 }$$

where ψ is the polarization angle. Then the timing residuals are formed with the use of the antenna pattern functions:

$$\begin{equation}{w}_{ki}\left({\boldsymbol{\theta }}_{\text{s}}\right)=-{F}_{+}\left({{\Omega}}_{k},{{\Omega}}_{\text{GW}}\right){\bar{H}}_{+}\left({t}_{ki};{\boldsymbol{\lambda }}^{\prime }\right)-{F}_{{\times}}\left({{\Omega}}_{k},{{\Omega}}_{\text{GW}}\right){\bar{H}}_{{\times}}\left({t}_{ki};{\boldsymbol{\lambda }}^{\prime }\right),\end{equation} \tag{ 9 }$$

where Ω_k = (θ_k , ϕ_k ) and Ω_GW = (θ_GW, ϕ_GW) are the sky location of the kth pulsar and the GW source, respectively. F_+,×(Ω_k , Ω_GW) are the antenna pattern functions, and λ' contains the internal parameters of the GW signal.

Note that equation (9) assumes that only the Earth term (see e.g. [40]) contributes to the signal. This simplification can be justified as follows: the time difference between the Earth and pulsar terms is given by

$$\begin{equation}{\Delta}t=L\left(1-\mathrm{cos}\enspace \beta \right)\approx 50\enspace \mathrm{y}\mathrm{r}{\left(\frac{\beta }{10{\circ}}\right)}^{2}\left(\frac{L}{1\enspace \mathrm{k}\mathrm{p}\mathrm{c}}\right),\end{equation} \tag{ 10 }$$

where L is the distance to the pulsar, β is the angle between the pulsar and the GW source sky location, and for the second equality we made an approximation valid for small β (which is where we can get the smallest delays). We can see that for typical millisecond pulsar distances of 1 kpc and observing campaigns not longer than a few decades, we will not see the pulsar term signals unless the sky location of the GW source is within a few degrees of a pulsar. Even if that happens, it will only appear in that single pulsar, so it can be modeled as a noise transient.

Note that pulsar term signals from GW bursts with no Earth term observed will appear as noise transients in our dataset. For an array with N_p pulsars, these correspond to an effective observation time N_p times larger than the actual observation time. However, the signal-to-noise ratio (SNR) of a pulsar term signal is a factor of ${N}_{\text{p}}^{1/2}$ lower compared to the Earth term signal, since it appears in just a single pulsar. This smaller SNR corresponds to a decrease in the observable volume by a factor of ${N}_{\text{p}}^{3/2}$. As a result, we expect that the rate of pulsar term bursts will be a factor of ${N}_{\text{p}}^{1/2}$ lower than Earth term bursts, corresponding to a factor of 7 for the 45 pulsars in the NANOGrav 12.5 year dataset. Pulsar term signals can also be coherent between two (or a few) pulsars, but for an array with a few tens of pulsars, the expected rate of such events are even lower [29].

2.2. Sampler

As described above, the dimensionality of our model depends on whether the GWB is turned on or not, and also on how many wavelets are used in the signal and the transient noise model. A GWB can carry one or two parameters depending on whether the spectral slope of the power-law model is fixed or varied. Each noise transient wavelet carries six parameters including the five internal parameters plus an 'external' parameter indexing which pulsar the noise transient belongs to. The signal model has three external parameters common to all wavelets (sky location and polarization angle) plus seven internal parameters for each wavelet.

To sample this variable-dimension parameter space we use an RJMCMC sampler called BayesHopperBurst, which not only gives posterior distributions on the model parameters, but also on the number of different components included in the model. BayesHopperBurst uses many of the proposals developed for BayesHopper such as Fisher matrix proposals, parallel tempering and trans-dimensional proposals. The description of these can be found in [13]. We also need a global proposal to help convergence, however the one used in BayesHopper is specific to continuous wave sources, so instead we use a global proposal based on the τ-scan, which we describe below.

The τ-scan is defined as:

$$\begin{equation}\mathcal{T}\left(\left\{{x}_{i}\right\};{t}_{0}^{\prime },{f}_{0}^{\prime },{\tau }^{\prime }\right)=\sum\limits _{i=1}^{{N}_{\text{p}}}\vert {\left({x}_{i}\vert {{\Psi}}_{\mathrm{c},i}\right)}_{i}{\vert }^{2}+\vert {\left({x}_{i}\vert {{\Psi}}_{\mathrm{s},i}\right)}_{i}{\vert }^{2},\end{equation} \tag{ 11 }$$

where x_i is the time series in the ith pulsar, N_p is the number of pulsars in the array, (.|.)_i is the noise-weighted inner product in the ith pulsar, ${{\Psi}}_{\mathrm{c},i}={\Psi}\left(A={\bar{A}}_{i},{t}_{0}={t}_{0}^{\prime },{f}_{0}={f}_{0}^{\prime },\tau ={\tau }^{\prime },{\phi }_{0}=0\right)$ is a 'cosine wavelet' and ${{\Psi}}_{\mathrm{s},i}={\Psi}\left(A={\bar{A}}_{i}^{\prime },{t}_{0}={t}_{0}^{\prime },{f}_{0}={f}_{0}^{\prime },\tau ={\tau }^{\prime },{\phi }_{0}=\pi /2\right)$ is a 'sine wavelet', where ${\bar{A}}_{i}$ and ${\bar{A}}_{i}^{\prime }$ are normalization constants chosen for each pulsar so that (Ψ_c,i |Ψ_c,i ) = (Ψ_s,i |Ψ_s,i ) = 1.

By calculating $\mathcal{T}$ for various t'₀, f'₀, and τ' values on a grid, we can get a 3D map showing where the data has some excess power. To determine the necessary spacing between grid points so that they properly cover the parameter space, we calculate the match (M) between two neighboring wavelets and require that it is above some threshold, e.g. M > 0.9. The match between two waveforms (h and h') is defined as:

$$\begin{equation}M=\frac{\left(h\vert {h}^{\prime }\right)}{\sqrt{\left(h\vert h\right)\left({h}^{\prime }\vert {h}^{\prime }\right)}},\end{equation} \tag{ 12 }$$

where $\left(.\vert .\right)={\sum }_{i=1}^{{N}_{\text{p}}}{\left(.\vert .\right)}_{i}$ is the array-wide noise-weighted inner product. M is normalized so that −1 ⩽ M ⩽ 1, where M = 1 indicates that h and h' are exactly the same, while M = −1 means there is a perfect anticorrelation between them. We can approximate M for two wavelets close to each other as [15]:

$$\begin{equation}M=1-\frac{1}{4{\bar{\tau }}^{2}}{\left({\Delta}\tau \right)}^{2}-\frac{1}{2{\bar{\tau }}^{2}}{\left({\Delta}{t}_{0}\right)}^{2}-\frac{{\pi }^{2}{\bar{\tau }}^{2}}{2}{\left({\Delta}{f}_{0}\right)}^{2},\end{equation} \tag{ 13 }$$

where Δ indicates the difference in the given parameter between the two wavelets and $\bar{\tau }=\left({\tau }_{i}+{\tau }_{j}\right)/2$ is the average τ between the two wavelets.

We can then require that the mismatch (1 − M) when moving only one of the parameters be less than 0.1. This implies:

$$\begin{equation}{\Delta}\mathrm{ln}\enspace \tau {\leqslant}\sqrt{\frac{2}{5}},\end{equation} \tag{ 14 }$$

$$\begin{equation}{\Delta}{t}_{0}{\leqslant}\frac{\tau }{\sqrt{5}},\end{equation} \tag{ 15 }$$

$$\begin{equation}{\Delta}{f}_{0}{\leqslant}\frac{1}{\sqrt{5}\pi \tau }.\end{equation} \tag{ 16 }$$

We see that the τ layers should be logarithmically spaced with each layer having a τ value maximum $\mathrm{exp}\left(\sqrt{2/5}\right)\approx 1.9$ times larger than the previous one. Time and frequency should be linearly spaced, with the spacing depending on the τ of the given layer. One can always choose to make the spacing finer, but that results in increasing computational costs.

Figure 1 shows the τ-scan calculated for a white noise burst (WNB) signal (see section 3.1) with central time and frequency marked with the white crosses. Here we use five layers with different τ values. We can see that as a result of the τ-dependent time and frequency spacing, pixels at higher τ values are more elongated in time but narrower in frequency (resulting in the same time–frequency area). We can see that the τ-scan clearly lights up around the injected signal indicating that, after proper normalization, this provides an excellent proposal. Indeed, we find that including this proposal speeds up convergence by multiple orders of magnitude, even though the remaining parameters of the wavelet are drawn from their priors. This proposal is also used when proposing to add a new signal wavelet, which makes it possible to have a reasonable acceptance rate for these complicated jumps.

Zoom In Zoom Out Reset image size

A variant of the τ-scan proposal is used for noise transient wavelets. For these we create a τ-scan for each pulsar only using the data for the pulsar in question. These are not only used to draw f₀, t₀ and τ values for the wavelet, but also to draw which pulsar to put in the wavelet based on how much excess power is in each pulsar.

To test BayesHopperBurst we analyze a few different simulated datasets. All the datasets in this section correspond to a simulated PTA with 20 pulsars at random sky locations. Each has been observed for 10 years every 30 days, has the same constant white noise level of 0.5 μs and no red noise. We add different kinds of GW signals and noise transients to this base dataset. In section 3.1 we show results for datasets only containing GW signals and not noise transients; in section 3.2 we look at a scenario with no GW signals but noise transient in many pulsars of the array; while in section 3.3 we look at a few examples of how BayesHopperBurst can deal with datasets containing both GW signals and noise transients.

3.1. Signal without noise transients

First we look at a dataset with a GW signal that consists of two sine-Gaussian wavelets that are close to each other in the time–frequency plane. This can be considered as one of the easiest signals to reconstruct, since it matches our model, i.e. we can perfectly reconstruct such a signal with two signal wavelets. Figure 2 shows the reconstruction of this signal in one of the pulsars in the array. We can see that BayesHopperBurst uses no noise transient wavelets and two signal wavelets most of the time, as expected. We can also see that the reconstructed signal matches well the injected GW signal. Note that while figure 2 only shows the reconstruction in a particular pulsar, the projection of the signal in the datastream of other pulsars can be significantly different. The reconstructions for all 20 pulsars in the array are shown in figure A1, indicating similarly accurate waveforms in all pulsars. To quantitatively characterize the goodness of fit throughout the whole array, we calculate the match between the median reconstructed signal and the injected signal, which we find to be 0.96. The SNR of the simulated signal is 17.9, which we define as $\mathrm{S}\mathrm{N}\mathrm{R}=\sqrt{\left({h}_{\text{inj}}\vert {h}_{\text{inj}}\right)}$, where h_inj is the waveform of the simulated (or injected) signal.

Zoom In Zoom Out Reset image size

The next dataset contains a WNB, which is produced from a white noise time series by applying a Gaussian window both in time and frequency. The variance of the white noise sets the amplitude of our WNB, while the windowing determines its time–frequency location and spread. While a WNB is similar to a sine-Gaussian wavelet in that it is also well localized on the time–frequency plane, it has a more complicated time dependence. This also means that a WNB cannot be perfectly reproduced by a finite number of sine-Gaussian wavelets, so the number of wavelets used will be determined by the trade-off between getting an increasingly better fit and paying an Occam penalty for each additional wavelet included.

Figure 3 shows the reconstruction of a WNB signal. We can see that BayesHopperBurst uses a single wavelet to fit the majority of the signal, while low-amplitude parts of the signal are missed. This also leads to a lower match between the median reconstructed signal and the injected signal (0.812) compared to the previous example, even though it has a slightly higher SNR of 18.3. This suggests that BayesHopperBurst performs slightly worse for a WNB signal compared to sine-Gaussians. Note however that if we take into account the uncertainty in the linearized timing model (orange band on figure 3), the reconstruction is consistent with the injected GW signal at the 90% confidence level.

Zoom In Zoom Out Reset image size

One of the potential sources of GW bursts in the PTA band are SMBHBs on parabolic or highly eccentric orbits. Here we analyze the simulated signal from a parabolic encounter of two SMBHs. We calculate the waveform from equations (3) and (4) of [26]. Beyond being a potential astrophysical source, parabolic SMBHBs are also a good test of the algorithm, because they are not elliptically polarized, so BayesHopperBurst must use its flexibility to fit two very different waveforms for the plus and cross polarizations.

Figure 4 shows the waveform reconstructed by BayesHopperBurst for an equal mass SMBHB on a parabolic orbit with a total mass of M = 2 × 10⁹ M_⊙, an impact parameter of b = 60M, at a distance of 15 Mpc, with an SNR of 16.8. Unlike wavelets or WNBs, these parabolic encounter waveforms have a significant low-frequency component. We can see that the linearized timing model can fit this component with a quadratic function by changing the period and period derivative parameters of the timing model. The high-frequency (and high-amplitude) part of the signal is reconstructed by a single wavelet. The match between the median reconstructed signal and the injected signal is 0.956.

Zoom In Zoom Out Reset image size

3.2. Multiple noise transients

We have seen that BayesHopperBurst is able to reconstruct GW bursts in the presence of Gaussian noise. In this section, we test if it can correctly identify incoherent noise transients. We consider four different datasets: (i) with a single pulsar containing a sine-Gaussian noise transient; (ii) with three pulsars containing the same noise transient; (iii) with ten pulsars (half the array) containing the same noise transient; and (iv) all 20 pulsars containing the same noise transient. Having the same noise transient in many pulsars can be considered a worst case scenario for discriminating signals and noise transients. It is virtually impossible to have such a situation by chance alignment of incoherent noise transients, but it provides a good test case. Note that even if all the pulsars have the exact same noise transient, in principle, that can be distinguished from a coherent GW burst, because there is no sky location and polarization that would result in the exact same waveform (with the same amplitude) in each pulsar.

Figure 5 shows the histogram of the number of wavelets used in the signal model and the noise transient model for each of these datasets. We can see that for the first three datasets we are recovering the correct model of having 1/3/10 noise transient wavelets and no signal wavelets. For the dataset with all the pulsars containing the same noise transient we instead see that BayesHopperBurst is using a single signal wavelet and 12 to 15 noise transient wavelets. This can be understood as a result of the inbuilt parsimony of Bayesian statistical analysis. As noted above, there are no such external parameters that could result in the exact same amplitude in the datastream of all the pulsars. However, one can use the signal model to account for noise transients in a few pulsars, and use noise transient wavelets to model the rest. A single-wavelet signal model has ten parameters, while each noise transient wavelet has five parameters. Thus such a solution could be statistically favored if the signal model can account for noise transients in more than two pulsars. Looking at the number of noise transient wavelets used, we can see that in this example the signal model is able to mimic the residuals from up to eight noise transients.

Zoom In Zoom Out Reset image size

3.3. Signal and noise transients

We have seen how BayesHopperBurst can reconstruct signals or noise transients. In this section we test it on datasets containing both signals and noise transients to see how it is able to distinguish between the two. First we analyze a dataset containing a signal and a noise transient in one of the pulsars, both of which have sine-Gaussian waveforms but with significantly different t₀, f₀ and τ parameters. We produce two versions of this dataset with different amplitudes: (i) with signal SNR of 26.8, and noise transient SNR of 10; (ii) with signal SNR of 8.9, and noise transient SNR of 5. The reconstructions for these datasets are shown in figure 6. We can see that both the signal and the noise transients are correctly identified and reconstructed for the higher amplitude case, and we find M = 0.994 between the injected and recovered waveforms. For the lower amplitude case only the signal has been found with a moderate match of 0.87. Note however that both the relatively lower match value and the fact that no noise transient has been found are consistent with the significantly lower SNRs of both signal and noise transient in this dataset.

Zoom In Zoom Out Reset image size

In the next example, we create a dataset with a sine-Gaussian signal (SNR = 16.1) and a noise transient (SNR = 10.8) in one of the pulsars with the exact same waveform as the signal. This can be considered a worst case scenario, since BayesHopperBurst can only rely on the coherence of the signal when trying to distinguish GW signal and noise transient. We make two versions of this dataset with the noise transient being in two different pulsars. Figure 7(a) shows the reconstruction for the case where the noise transient is in a pulsar far away from the GW source on the sky as shown by the sky location posterior map (see figure 7(b)). We can see that both the waveforms and the sky location is accurately reconstructed in this case (M = 0.98). However, having a noise transient in a pulsar that lies close to the true GW source sky location results in no significant noise transient reconstructed (see figure 7(c)). Instead, BayesHopperBurst fits the data by changing the sky location slightly to account for the effect of the noise transient (see corresponding skymap in figure 7(d)). Since the noise transient is in one of the pulsars closest to the GW source on the sky, the location of the GW source can be changed slightly in a way that it can account for the noise transient in that pulsar but not change significantly the projection of the signal in other pulsars. This reconstruction results in a match of 0.96, which is only slightly lower than the previous example, even though it is not using the true model.

Zoom In Zoom Out Reset image size

In a follow-up study, we plan to use BayesHopperBurst to search for nHz GW bursts in the NANOGrav 12.5 year dataset [2]. To assess the sensitivity we might expect from such an analysis we run BayesHopperBurst on a simulated dataset made to resemble the NANOGrav 12.5 year dataset and used in [41]. The dataset has been made using the properties (timing parameters, white and red noise parameters) of pulsars and observational time stamps from the NANOGrav 12.5 year dataset. Note that the red noise parameters are derived by filtering out a common red noise process, so that they truly characterize pulsar intrinsic red noise. To reduce computational costs the residuals are epoch averaged. More details on this simulated dataset can be found in [41].

This dataset does not contain any GW bursts, so by running BayesHopperBurst on it we can set upper limits on the amplitude of generic GW transients. We consider these upper limits to be indicative of the upper limits we might expect to get from running on the real NANOGrav 12.5 year dataset assuming we do not find any significant GW candidates. Figure 8 shows the upper limits we get from this simulated dataset as a function of central time (t₀) and central frequency (f₀) defined e.g. by equations (8a) and (8b) in [42], which coincide with the parameters used to describe the wavelets if only one wavelet is used. We quote these upper limits in terms of the root-sum-squared strain amplitude defined as:

$$\begin{equation}{h}_{\mathrm{r}\mathrm{s}\mathrm{s}}=\sqrt{{\int }_{-\infty }^{\infty }\left[{h}_{+}^{2}+{h}_{{\times}}^{2}\right]\mathrm{d}\mathrm{t}}=\sqrt{{\int }_{-\infty }^{\infty }\left[{\left(\frac{\mathrm{d}{\mathrm{H}}_{+}}{\mathrm{d}\mathrm{t}}\right)}^{2}+{\left(\frac{\mathrm{d}{\mathrm{H}}_{{\times}}}{\mathrm{d}\mathrm{t}}\right)}^{2}\right]\mathrm{d}t}.\end{equation} \tag{ 17 }$$

To produce this plot we calculate t₀, f₀ and h_rss for each sample from our MCMC. Then we bin the samples based on a grid on the t₀–f₀ plane, and in each bin we calculate the 95th percentile of the h_rss samples from that bin, which we quote as the 95% upper limit.

Zoom In Zoom Out Reset image size

A clear trend from figure 8 is that the upper limit decreases as we move to later times. This is because the array contains more and more pulsars as time goes on, which greatly improves the array's sensitivity. In addition, we see drastic improvements in the second half of the dataset, which can be explained by the improved backends used after ∼7.5 year and a high-cadence observing campaign, which started at ∼8.5 year (see e.g. [2]). We can also see that the sensitivity decreases at the temporal edges of the dataset, which is due to the fact that waveforms with a central time near the edges will have a significant portion that lies outside of the observational time span of the dataset. The worst sensitivity in the analyzed time–frequency region is achieved at very high frequencies and very early times. This is due to the fact that there are less than ten pulsars observed at those times. There is also a slight increase in the upper limit around the frequency corresponding to a period of one year (log₁₀  f₀ ≃ −7.5), where PTAs are known to be less sensitive due to the sky location fitting in the timing model.

The best sensitivity is achieved around t₀ ≃ 10 years and f₀ ≃ 8 nHz at an amplitude of h_rss ≃ 5 × 10⁻¹¹ Hz^−1/2. This can be converted into the minimum GW energy that can be detected at a fiducial distance using e.g. equation (4.2) of [43]. We get that at a fiducial distance of 100 Mpc we are sensitive to sources emitting an energy of more than ∼40M_⊙ c² in GWs.

Note that we use a uniform prior on the wavelet amplitude which pushes up the amplitude values resulting in a conservative upper limit. However, as we also allow the sky location to change, and the distribution of the pulsars is highly anisotropic, this also pushes the sampler to the least sensitive sky locations where the highest amplitude signals can remain hidden. This means that the sensitivity at the best sky location could be significantly better than the value quoted above. On the other hand, in this run we fix the red noise parameters in each pulsar to the true values due to computational constraints. This in principle could make our upper limit estimates overly optimistic as it does not take into account any covariances between the red noise and GW burst signal models. However, since the red noise only dominates at the lowest frequencies, we do not expect this to have a drastic effect at the most sensitive time–frequency location reported above. In addition, even while fixing the red noise parameters to their true values, we still allow the realization of the red noise to vary, which can accommodate some covariance between the signal and red noise models.

Real data can bring complications compared to simulated data like uneven sampling, temporal gaps in the datastream, chromatic effects, etc. To test how BayesHopperBurst can handle these complications we run it on a single pulsar (B1855 + 09) from the NANOGrav 9 year dataset [35]. We show the epoch-averaged residuals from B1855 + 09 in figure 9(a) along with the reconstruction of the linearized timing model and the red noise model. Figure 9(b) shows the epoch averaged residuals corrected with the median timing and red noise model. These indicate a transient feature around 55 000 Modified Julian Date (MJD), which is reconstructed by the noise transient model also shown on this figure. Note that although we show the epoch-averaged residuals for plot clarity, our analysis used all individual residuals.

Zoom In Zoom Out Reset image size

We can see in figure 9 that BayesHopperBurst adjusts the timing parameters such that some features of the data were smoothed out while others were amplified. Then it used a few wavelets to correct those amplified features. In particular, we can see that the observations made at different frequencies at the same epoch are sometimes inconsistent in the original dataset. One can correct that by changing the value of the dispersion measure parameter (DMX, see e.g. [35]) for that epoch, but that will result in a net residual in both observations. The legacy analysis thus cannot do that, since it cannot account for transient noise features. On the other hand, BayesHopperBurst can change the dispersion measure values and model the resulting noise transient with a collection of wavelets. Thus, our analysis suggests that there is a noise transient in this dataset. Such a noise transient was not found before, since the canonical analysis does not model noise transients and their presence can be masked by inflating the white noise parameters. In fact, BayesHopperBurst only finds the noise transient when we allow the white noise parameters to vary, so that it can lower their values compared to the canonical analysis.

Based on figure 9, it might look like that the feature around 55 000 MJD fitted by the noise transient model was introduced by the similar feature with opposite sign in the timing and red noise model. Note however that the linearized timing model introduces perturbations around timing model parameters previously determined by an analysis which does not allow for unmodeled noise transients. Such an analysis will try to adjust its free parameters to model out any noise transients present in the dataset. Our results indicate that if we introduce a model for noise transients, the Bayesian analysis prefers to undo that timing model adjustment, and instead model the feature with the dedicated noise transient model.

We have presented BayesHopperBurst, a Bayesian algorithm to search for nHz GW burst using PTA data, by modeling them as a collection of sine-Gaussian wavelets, where the number of wavelets is found by a trans-dimensional RJMCMC sampler. We demonstrated how it can reconstruct signals with a wide range of morphologies, and also how it can distinguish between coherent GW burst and incoherent noise transients that may occur in the datastream of some pulsars.

In the near future, we plan to use BayesHoperBurst to carry out the first ever search for generic GW bursts in the nHz frequency regime using the NANOGrav 12.5 year dataset. In preparation for that we analyzed a simulated dataset similar to the NANOGrav 12.5 year dataset and made predictions on the sensitivity we might expect from such a search. To demonstrate that BayesHopperBurst is capable of dealing with the additional complications arising when analyzing real data, we analyzed a single pulsar (B1855 + 09) from the NANOGrav 9 year dataset.

In addition to searching for generic GW bursts, an algorithm like BayesHopperBurst could also be a useful tool for analyzing transient noise features in the datastream of single pulsars as illustrated in [32]. We are planning to further explore the potential of such an approach in a follow-up study. Systematically analyzing a large collection of pulsars will also help us better understand the covariances we have seen in section 5 between transient noise, DMX, and red noise models.

Important areas to investigate in future studies are further development and optimization of BayesHopperBurst. One particularly interesting direction would be experimenting with different types of functions in our signal and noise transient models (e.g. shapelets [39]), which might be more suited for certain signals or noise transients.

The authors thank David Nice and Justin Ellis for informative discussions about B1855 + 09, and Sarah Vigeland for feedback on the manuscript. We appreciate the support of the NSF Physics Frontiers Center Award PFC-1430284. Some computations were performed on the Hyalite High Performance Computing System, operated and supported by University Information Technology Research Cyberinfrastructure at Montana State University.

The data that support the findings of this study are available upon reasonable request from the authors.

In section 3, we only show the waveforms reconstructed in a single pulsar (see figures 2–4, 6, and 7), but it is also interesting to look at how the signal appears in different pulsars. Figure A1 shows the simulated data and the reconstructed signal in all 20 pulsars of the array for the injection with two sine-Gaussian wavelets (cf figure 2). While the amplitude and waveform of the signal appearing in datastreams of different pulsars are significantly different, BayesHopperBurst consistently reconstructs the signal in each pulsar. This is not surprising given that we employ a coherent signal model which uses information from all the pulsars to constrain the waveform of the signal.

Zoom In Zoom Out Reset image size