Searching for gravitational wave memory bursts with the Parkes Pulsar Timing Array Free

Abstract

1 INTRODUCTION

As supermassive black hole binary (SMBHB) systems coalesce, they are expected to produce gravitational wave (GW) emission. At the time of coalescence, a permanent change in the space–time metric will propagate away from the source (Payne 1983; Christodoulou 1991; Blanchet & Damour 1992; Thorne 1992; Favata 2009). The permanent change is known as the ‘gravitational wave memory’ (GWM) effect. Throughout this paper, we mainly consider GWM events caused by SMBHB coalescences. However, GWM events can also come from other sources such as cosmic strings, supernovae or during the flyby of massive objects (Pshirkov, Baskaran & Postnov 2010).

The passage of a GWM past the Earth or a pulsar will cause a change in the observed frequency of that pulsar's rotation. By observing a sufficiently large number of stable millisecond pulsars, it is expected that an unambiguous detection of the GWM effect could be made. In this paper, we describe a GWM search algorithm and apply it to the recent Parkes Pulsar Timing Array (PPTA) data set (Manchester et al. 2013).

As part of the pulsar timing technique (for details, see Edwards, Hobbs & Manchester 2006), a ‘pulsar timing model’ is developed that describes the pulsar, its companions and the propagation of the pulses from the pulsar to the Earth. The model predicts the pulse times of arrival (ToAs) in an inertial reference frame. In all current pulsar timing experiments, this reference frame is taken to be the Solar system barycentre. The barycentric ToAs are then compared with the predictions of the timing model and the differences are identified and termed the ‘pulsar timing residuals’. As GWs are not, by default, included in the timing model, they will induce timing residuals.

It is thought that observations of millisecond pulsars will lead to the direct detection of GWs (Sazhin 1978; Detweiler 1979; Jenet et al. 2005, 2006). A promising class of GW source potentially detectable through pulsar timing is a stochastic background, which could be generated by an ensemble of individually unresolvable inspiralling SMBHBs scattered throughout the Universe (e.g. Sesana, Vecchio & Colacino 2008). Such GWs are inherently stochastic and will induce correlated noise-like structure in the timing residuals. Other classes of GW sources are continuous waves generated by relatively nearby and massive inspiralling SMBHBs (Sesana, Vecchio & Volonteri 2009; Ravi et al. 2012) and the GWM events (e.g. Seto 2009; Pshirkov et al. 2010; van Haasteren & Levin 2010; Cordes & Jenet 2012) that are the focus of this paper. Such GW signals are deterministic and can be included in pulsar timing models.

Several observing programme have now been started with the goal of observing a large number of pulsars with sufficient precision to detect GW signals (Jenet et al. 2009; Ferdman et al. 2010; Manchester et al. 2013). Such projects are known as pulsar timing arrays (PTAs; Romani 1989; Foster & Backer 1990) and currently three exist. The North American PTA (NANOGrav; McLaughin 2013) was formed in 2007 and carries out observations with the Arecibo and Green Bank Telescopes. The European Pulsar Timing Array (Kramer & Champion 2013) was established in 2004 and includes telescopes in England, France, Germany, the Netherlands and Italy. For this paper, we make use of data from the PPTA project (Manchester et al. 2013) which commenced in 2004 and uses the 64-m diameter Parkes radio telescope. Parkes observations have been used to place an upper bound on a stochastic background of GWs (Shannon et al. 2013) and to search for GW signals from individual, non-evolving, SMBHBs (Zhu et al. 2014).

In this paper, we focus on the GWM phenomenon. Seto (2009), van Haasteren & Levin (2010) and Pshirkov et al. (2010) have independently shown that PTAs would be sensitive to sufficiently strong GWM events. GWM events passing a pulsar will lead to a glitch event in the timing residuals of only that pulsar (Cordes & Jenet, 2012) and may be indistinguishable from a rotational glitch.¹ GWM events passing the Earth will lead to simultaneous glitch events that are potentially detectable in the timing residuals of multiple pulsars in the array. The size and sign of the glitches will depend upon the angle between the source, Earth and pulsar in a quadrupolar fashion. Thus, such events can be separated from rotational glitches in individual pulsars.

Along with GWs, many other physical phenomena are not included in the timing model and will also induce timing residuals that may mask the signals of interest. These include errors in the terrestrial time standard (see e.g. Hobbs et al. 2012), errors in the Solar system ephemeris (e.g. Champion et al. 2010) and uncorrected dispersion measure variations (see e.g. Keith et al. 2013). Pulsars are also known to exhibit intrinsic variations in the timing residuals which include stochastic spin noise (e.g. Hobbs et al. 2010; Shannon & Cordes 2010) and glitch events (e.g. Cognard & Backer 2004; Yu et al. 2013).

In this paper, we try to answer the following four questions.

• How can we detect GWM signals in pulsar data sets?

• Do GWM signals exist in the PPTA data set?

• If no signal is detected, then what is the maximum rate estimate of such GWM events?

• What are the astrophysical implications of this bound on the GWM amplitude?

In Section 2, we describe the observations used in this analysis. In Section 3, we describe the GWM signal. In Section 4, we present our detection algorithm for searching and limiting the GWM signal and answer question (i). In Section 5, we apply our algorithm to the PPTA data and present the results, answering question (ii). In Section 6, we discuss our results and their astrophysical implications. This leads to answers to questions (iii) and (iv). In the appendix, we describe updates made to the software package tempo2 needed for this work.

2 OBSERVATIONS

We use the initial PPTA data set (Manchester et al. 2013). This data set is available from the Commonwealth Scientific and Industrial Research Organization (CSIRO) data archive² and has a digital object identifier of 10.4225/08/534CC21379C12.³ The raw observations that made up that data release are also available from the same website as part of the Parkes pulsar data archive (Hobbs et al. 2011).

The data set includes regular observations of 20 millisecond pulsars at intervals of 2–3 weeks between the years 2005 and 2011.⁴ For each pulsar, ToAs for the band that has the lowest overall rms timing residuals after the data have been corrected for dispersion measure variations (Keith et al. 2013) have been selected. All the observations were performed with the Parkes 64-m radio telescope with typical integration times of 1 h.

Timing residuals were formed using the tempo2 software package (Hobbs, Edwards & Manchester 2006) making use of the JPL DE421 Solar system ephemeris (Folkner, Williams & Boggs 2008) and referred to terrestrial time as realized by the Bureau International des Poids et Mesures⁵ (BIPM2011). The post-fit timing residuals for the 20 pulsars are shown in Fig. 1.

Figure 1.

The post-fit timing residuals for the PPTA data set. The dashed, horizontal lines indicate zero residual. The pulsar name and the range of the timing residuals are labelled on each subplot.

Open in new tabDownload slide

The first three columns in Table 1 provide, for each pulsar, its name, pulse period and dispersion measure. Pulsar timing data sets vary significantly, and we currently do not have a simple way to quantify the quality of different data sets. Usually, the weighted rms of the timing residuals is used. We present this value, σ_w, in column 4. However, we note that pulsars scintillate and therefore some ToAs can have much larger uncertainties than other ToAs for the same pulsar. We therefore also present the unweighted rms timing residuals, σ_uw, in column 5. Both of these statistics are affected by any non-white noise process in the data and so the uncertainties on individual ToAs are often significantly lower than the rms values. To quantify this, we give, in column 6, the median ToA uncertainty for each pulsar, m_orig, as measured during the ToA determination procedure. Below, we show that additional white noise, which is not described by the ToA uncertainties, may also be present. After correction for such noise, the median uncertainty will be increased. This corrected median ToA uncertainty, m, is listed in column 7 of the table. The remaining three columns give the number of ToAs (N_obs), the dates of the first and last observations (as Modified Julian Dates, MJDs) and the data span, respectively.

All pulsar data sets have a number of properties that make searching for the GWM effect challenging. Although the data span is similar for each pulsar in the PPTA data set, the data sampling is irregular and is not the same for each pulsar. The ToA uncertainties are time variable and tend to decrease with time as new receivers are commissioned and/or wider bandwidths become available. The variability of the residuals is quite different between pulsars (depending on, for instance, the scintillation properties of that pulsar) and low-frequency fluctuations in the residual time series are evident in many pulsars. Such low-frequency variations are probably dominated by stochastic spin noise, but can also be caused by imperfect correction for dispersion measure variations or slight errors that arise when combining data sets from different observing instruments.

Figure 2.

Top panel: power spectrum density of the timing residuals for PSR J1643−1224 (solid line). A model of the red noise is indicated as the dotted line. The dashed line is the mean power spectrum of 100 simulations of the white and red noise. Lower panel: power spectrum for the whitened residuals of PSR J1643−1224 (solid line). The expected mean and ±2σ confidence intervals are shown as horizontal lines.

Open in new tabDownload slide

The objective of this noise modelling is to find a linear transformation that whitens and normalizes the residuals. The efacs and equads used to model the white noise are not unique, nor is the red noise spectral model. The test we use for a satisfactory linear transformation is that the power spectrum of the whitened residuals fits within the ±2σ error bars for a random variable with 2 degrees of freedom. We use the Lomb–Scargle algorithm to estimate this spectrum because it is unbiased for a truly white time series. The power spectrum for PSR J1643−1224 after whitening is shown in the lower panel of Fig. 2. We have included the parameters of the red noise models we used in Table 2 so that others can duplicate our analysis. The efac and equads are tabulated in Appendix. We note that the development of these noise models is subjective and discuss the implications of this in Section 6.2.

3 THE GWM SIGNAL IN PULSAR TIMING

In the example of two coalescing equal mass black holes, the amplitude of the GWM signal grows rapidly. The growth time-scale of the metric change is ∼10⁴ s (M/10⁸ M_⊙)(1 + z) (van Haasteren & Levin 2010), where z is the redshift of the source and M is the mass of each black hole (the black holes are assumed to have equal mass). This is short compared with typical observation intervals for existing PTAs (which is normally one observation every 2–3 weeks). We assume such growth time-scales for all GWM sources and treat the signal as a discrete jump of the metric propagating through space.

Figure 3.

The geometry used to describe the GWM emission. The top panel has the Earth at the centre and the GWM source indicated by the star symbol. The lower panel represents the GWM emission coming out of the page.

Open in new tabDownload slide

We have incorporated the effect of a GWM event into the tempo2 timing model. This allows us to include the GWM in a fit and to simulate residuals (or ToAs) that include such an event. The new timing model parameters are (h^mem, t₀, α_g, δ_g, ζ). Tempo2 uses a linear least-squares-fitting algorithm. If the GWM epoch, polarization angle and source position were known, it would be possible to fit for the amplitude of the GWM source as part of the standard tempo2 timing fit. However, if these parameters are not known, then a non-linear fitting routine is needed to determine their values.

We have found this parametrization convenient for simulating timing residuals affected by a GWM signal. However, in order to search for the GWM events, we have found it useful to provide a second parametrization of the GWM effect within tempo2. In this parametrization, we describe the GWM using two orthogonal components, A₁ and A₂ where A₁ = h^memcos (2ζ) and A₂ = h^memsin (2ζ). This formulation has the advantage that A₁ and A₂ enter the timing model linearly and can be fitted with linear least squares. We emphasize that, even with this parametrization, the position of the source and epoch of the event cannot be obtained using a linear-fitting routine. We therefore fit for A₁ and A₂ at a grid of points for every possible sky direction and epoch.

As an example, we show in Fig. 4 simulated timing residuals for five pulsars. These timing residuals were formed by determining ToAs that are exactly predicted (i.e. yield zero residual; see Hobbs et al. 2009) by a timing model that included a GWM event (with h^mem = 3 × 10⁻¹³, α_g = 21^h51^m, δ_g = −30| $_{.}^{\circ}$|3) for each pulsar. We subsequently added 100 ns of Gaussian white noise to each observation (an observing cadence of 14 d was assumed). The resulting residuals, shown in the figure, were obtained from the initial pulsar timing models that did not include the GWM event. The GWM event is clearly seen in the centre of the data span. As expected, the size of the induced timing residuals depends upon the pulsar position. The same data set is shown in Fig. 5 after fitting for the pulsars’ pulse frequency and frequency derivative parameters.

Figure 4.

Simulated timing residuals for five pulsars. The pulsars are affected by white noise and a GWM event that occurred at the centre of the data span. No pulsar parameters have been fitted to the timing residuals. The value underneath the pulsar's name gives the range of the timing residuals for each pulsar.

Open in new tabDownload slide

Figure 5.

As in Fig. 4, but each pulsar's pulse frequency and its derivative have been fitted and post-fit timing residuals are shown. The value underneath the pulsar's name gives the range of the timing residuals for each pulsar.

Open in new tabDownload slide

4 THE DETECTION ALGORITHM

The response of each pulsar to a single GWM burst is completely determined by specifying the source parameters: A₁, A₂, position and epoch. For a given position and epoch, we jointly fit A₁ and A₂ with the pulsar parameters (spin-down, astrometry, orbital configuration, etc.) by adding the GWM response to the tempo2 timing model and minimizing the whitened timing residuals. We use the algorithm described in Coles et al. (2011) to account for the correlations in the pre-fit timing residuals caused by unmodelled red noise. To determine the sky position and epoch, we search over a regular three-dimensional grid whose spacing we describe below.

At each position and epoch, the tempo2 fit returns the parameter vector |$\boldsymbol {A} = [A_1; A_2]$| and their covariance matrix, |${\bf C}_{\rm o}$|⁠. From these, we require a detection statistic which provides an optimal estimate of the amplitude of the GWM. We can then use that statistic to locate the GWM in the grid of possible positions and epoch. While this approach is not as computationally efficient as a non-linear fit, it provides an opportunity to study the statistics of the noise by examining the response over the entire three-dimensional grid.

4.1 A detection statistic

A₁ and A₂ can be viewed as the GWM amplitude modified by the response of our ‘detector’. If the pulsars in the array were distributed uniformly in position and brightness, A₁ and A₂ would be orthogonal (independent) and equally sensitive. For this ideal case, |$D\equiv A_1^2 + A_2^2$| is an optimal detection statistic.

If |${\bf C}_{\rm o}$| is exact then, in the absence of a GWM signal, D is the sum of the squares of two unit variance Gaussians, and thus follows a χ² distribution with 2 degrees of freedom. We note that the mean and standard deviation of such a distribution both equal two, its probability density is exponential p(D) = (1/2)exp ( − D/2) and its cumulative probability is c(D) = 1 − exp ( − D/2). If one chooses a detection threshold, D^★, then the false-alarm probability is 1 − c(D^★) = exp ( − D^★/2). For a 5 per cent false-alarm probability, D^★ = 6.

Two factors prevent the actual statistic from following the ideal distribution. First, we estimate |${\bf C}_{\rm o}$| from the data, and the statistical uncertainty in the correction of |$\boldsymbol {A}$| to |$\boldsymbol {A}_{\rm w}$| increases the variance and biases D. As the observation span increases and |${\bf C}_{\rm o}$| becomes better characterized, the distribution of D approaches χ². Secondly, as discussed earlier, we have not actually fitted optimized noise models to each pulsar. Each pulsar noise model requires five parameters and we are reluctant to fit an additional 100 parameters to the data set and reduce the degrees of freedom by a corresponding amount. An error in the noise model will not bias the parameters A₁ or A₂, but it will alter the estimate of |${\bf C}_{\rm o}$|⁠, biasing D. In both cases, as long as the shape of the distribution is unchanged, the expected false-alarm probability can be recovered by renormalizing D such that its mean is 2.

Our statistical analysis must also account for the requirement to search over a grid of possible positions and epochs. We expect detectable GWM events to be extremely rare (see Section 6.3), and so we adopt D_max, the maximum D in a search over epoch and position, as our final detection statistic. Because the angular and time resolution of our array is modest, only a handful of epochs and positions are independent, and we must carry out simulations to determine the false-alarm probability of D_max.

4.2 Search over possible sky positions and GWM arrival times

The array sensitivity is non-uniform in both position and epoch. In particular, sensitivity to GWM events drops to zero at the edges of the observational span, as shown in Fig. 6, which depicts the sensitivity for equally spaced observations with regular white noise. This is discussed by van Haasteren & Levin (2010), and our results concur with theirs. Accordingly, we restrict our search over epoch to the central 80 per cent of the observing span, shown in the figure with a dotted line representing the region with roughly constant sensitivity.

Figure 6.

The normalized S/N ratio for detection of a GWM signal in white noise as a function of the epoch of the event (calculated by fitting the amplitude of a ramp function that occurs at the specified epoch and then subsequently dividing the value by its uncertainty). The S/N ratio is normalized to the average over the range 10–90 per cent of the observing interval.

Open in new tabDownload slide

To determine N_tot, we performed 1000 simulations of the full PPTA data set using the red and white noise models discussed earlier. Each of the 1000 realizations was normalized, as discussed above, so the sample mean D = 2 for each realization. We then find D_max for each realization and compute the cumulative distribution of D_max, shown in Fig. 7. We determine N_tot by matching this distribution to equation (11), yielding N_tot = 80 and |$D_{\rm max}^* = 14.7$| for a 5 per cent false-alarm probability. In these simulations, we used N_e = 16 and N_s = 34, so the simulations were oversampled by approximately a factor of 7.

Figure 7.

False-alarm probability (solid blue line) 1 − c(D_max) obtained from 1000 simulations of the PPTA data set with both red and white noise. The blue dashed lines are the ≈2σ statistical errors on the measurement. The red line is the theoretical expression for N_tot = 80. For a false-alarm probability of 5 per cent (horizontal dotted line), D_max = 14.7 (vertical dotted line).

Open in new tabDownload slide

To verify that our rescaling of D maintains the correct false-alarm probability, we examined the 544 000 simulated raw D values. While the sample mean was |$\bar{D} = 2.35 \pm 0.01$|⁠, the normalized D values matched the expected |$\chi _2^2$| distribution within the error bars.

5 APPLICATION TO THE PPTA OBSERVATIONS

We applied our algorithm to the PPTA data set by calculating D over a fine grid with 1034 sky positions and 150 epochs. Our grid is therefore significantly oversampled compared with the 80 independent degrees of freedom determined above. The additional computational burden in this oversampling is insignificant for a single ‘realization’ and allows us to study the correlations in detail. We demonstrate in Section 6.1 that statistically identical results are obtained using a smaller grid.

For each grid position, we perform a global fit for A₁ and A₂ whilst simultaneously fitting for the parameters specific to each pulsar. The mean, reduced χ² of the timing model fits is 0.9 which indicates that the noise models were reasonably accurate, but slightly conservative, i.e. about 10 per cent higher in variance than the observations. The mean value of the unnormalized statistic, |$\bar{D} = 2.1$|⁠, is lower than the mean value of the simulations (2.35), but within 1σ of it. We therefore rescaled our values of D by a small factor of 2.0/2.1 = 0.95 for the remainder of the analysis.

For each sky position, we calculated D_e, the maximum statistic determined over all epochs. These values are shown as a function of sky position in Fig. 8. The pulsar positions are indicated on this figure using star symbols. The peak D_e over the sky is D_max = 12.4 (corresponding to a false-alarm probability of 15 per cent) and occurred at MJD 54986 (corresponding to 2009 June 04) at (α_g = 02^h24^m, δ_g = −15| $_{.}^{\circ}$|8).

Figure 8.

Normalized detection statistic, D_e, for the PPTA observations measured at each of 1034 sky positions. The expected maximum value of D_e (see the text) is D^★ = 14.7. The locations of the pulsars are denoted by star symbols. The white box indicates the position corresponding to the maximum value of D_e.

Open in new tabDownload slide

Our value of D_max = 12.4 is lower than the 5 per cent false-alarm probability threshold of 14.7. We therefore do not claim a detection of any GWM event in the current PPTA data set.

5.1 Sensitivity over the sky

Using our value of D_max = 12.4, we can determine the sensitivity of the data set to GWMs as a function of sky position, i.e. we determine the amplitude, h^mem(α_g, δ_g), of a GWM event that would give D = D_max. We assume, as discussed earlier, that the sensitivity is independent of epoch provided the epoch is restricted to the central 80 per cent of the observing interval. The sensitivity h^mem(α_g, δ_g) is expected to vary significantly over the sky because the pulsars are not distributed uniformly, do not have the same ToA precision nor do they have the same red noise properties.

Figure 9.

Sensitivity of the PPTA data set to GWM events. Note that the scale is in 10⁻¹³ leading to a maximum sensitivity of ∼10⁻¹⁴ and a worst case sensitivity of ∼3 × 10⁻¹³.

Open in new tabDownload slide

5.2 Bounding

Fig. 9 shows that we could detect events at a level of h^mem = 2.4 × 10⁻¹³ anywhere on the sky (with 95 per cent confidence) if they occurred during our five-year effective observing interval. Our data set is sensitive to events with amplitude h^mem ∼ 2 × 10⁻¹⁴, but only over a smaller area of the sky.

For instance, from Fig. 9 we can say that no event of h^mem ≥ 2.4 × 10⁻¹³ occurred during our data set. This gives an event rate of λ < 0.75 yr⁻¹ for events of this size. Events at a level h^mem < 2 × 10⁻¹⁴ could not be detected at all, and no bound at that level can be set. The bounds on λ obtained from the PPTA data are shown in Fig. 10.

Figure 10.

The solid line shows upper bounds on the event rate λ as a function of the GWM amplitude bound using the PPTA data set. The dotted and dot–dashed lines at the bottom of the plot are predictions (detailed in the text) based on models of coalescing SMBHs. The red triangles are event rate predictions from Ravi et al. (2014). The red crosses are the predictions from Cordes & Jenet (2012). Results of simulations are also shown as dotted or dashed green lines with symbols, with the lower set based on simulated 20-year data sets with either optimistic or pessimistic extrapolation of the red noise model (see Section 6.5 for details).

Open in new tabDownload slide

6 DISCUSSION

6.1 Choice of grid sampling

In the previous section, we calculated D over a fine grid with 1034 sky positions and 150 epochs. This number of grid points is significantly greater than the 80 independent degrees of freedom, but does allow an easy way to make sky-maps such as those in Figs 8 and 9. However, in order to confirm that our astrophysical results are not significantly affected by such oversampling, we re-analysed the data in an identical manner, but first using 34 sky positions and 16 epochs, and secondly using 16 sky positions and 34 epochs.

Using the procedure described above, we obtained, for the first case, |$\bar{D} = 1.96$|⁠. Rescaling by 1.96/2.0 gave D_max = 9.9 (compared with D_max = 12.4 from the initial grid). This new value is also lower than the 5 per cent false-alarm probability threshold and we therefore note that this result is consistent with that obtained from the larger grid. To conclude this analysis, we obtain a single sky-averaged bound on the event rate of 0.75 yr⁻¹ for h^mem ≥ 1.7 × 10⁻¹³. For the second case, we obtain |$\bar{D} = 1.95$| and, after rescaling, D_max = 8.4. The sky-averaged bound on the event rate of 0.75 yr⁻¹ for h^mem ≥ 1.6 × 10⁻¹³. We therefore conclude that, as expected, our results do not significantly change with the number of grid points used.

6.2 Sensitivity of D_max to noise models

As we were carrying out the research for this paper, we discovered minor flaws in our initial models for the red and white noise for some pulsars. The models were corrected as the flaws became apparent, but the results we obtained with those initial models allow us to determine the effect of slightly incorrect noise models.

The penultimate noise models were the same as presented in Table 2 except for two pulsars. For PSR J1909−3744, we initially included a red noise model (with α = 4, P₀ = 1.2 × 10⁻²⁹ and f_c = 0.5) and slightly different efac and equad values. For PSR J2129−5721, we originally used only a white noise model. Following our detection procedure using these original noise models leads to an unnormalized D_max = 18 (compared with the unnormalized D_max = 12.9 using the final models). After normalizing the statistic, we obtained D_max = 14 (compared with 12.4 with our final models).

This highlights that even with poor noise models the statistics only change slightly. However, the poor noise models can lead to detection statistics that do slightly exceed the expected false-alarm probability and therefore any ‘detection’ should be treated with caution. The following checks can be carried out on such detections.

• Ensure that the detection does not depend solely on a single pulsar. That pulsar could exhibit a glitch event that would produce a strong response in our detection statistic.

• The detection must exhibit the expected signature in the position and epoch parameter space. A burst will have a deterministic pattern in A₁ and A₂ which depends upon the actual (α_g, δ_g, t₀) of the burst. One could then subtract this pattern from the grid of observations. The residuals should be entirely caused by noise. If the maximum detection statistic is not reduced, then there is either another burst present (which would be unlikely, but could be tested by iteration) or the original burst is not real.

• As long as the GWM event can be distinguished from the low-frequency timing noise, the significance of any such detection would increase with longer data spans and by including extra pulsar data sets.

• It is likely that any such GWM event would be caused by the coalescence of SMBHs. It is possible that such an event may be associated with an observable signal using other telescopes (see e.g. Burke-Spolaor 2013). As the ability to localize the source improves with time, it therefore may be possible to identify the host galaxy of such an event.

6.3 Astrophysical implications

Our results imply that fewer than 0.75 events yr⁻¹ have occurred with h^mem > 10⁻¹³ and fewer than ∼3 events yr⁻¹ for h^mem > 4 × 10⁻¹⁴. In this section, we compare these results with predictions for the event rate from astrophysical models. We note that it is extremely difficult to predict the event rates from existing observational data. Instead, we require the use of models for SMBHB mergers and then predict the event rates and amplitudes. Those predictions also rely on assumptions of the sources, such as their distances, black hole mass ratios, etc. Two papers have already made predictions. Cordes & Jenet (2012) predict an event rate of 0.4 yr⁻¹ for h^mem > 10⁻¹⁶ and 0.02 yr⁻¹ for h^mem > 10⁻¹⁵. Ravi et al. (2014) used the best available observational data on black hole masses and galaxy merger rates to conclude that only ∼10⁻⁵ bursts yr⁻¹ with h^mem > 5 × 10⁻¹⁵ and ∼10⁻³ bursts yr⁻¹ with h^mem > 2 × 10⁻¹⁵ are expected. These predictions are both so much lower than our current bounds that (1) we are unlikely to detect such events in the near future and (2) we can make predictions on the ‘time-to-detection’ using a simple toy model of the GWM event rate. We therefore note that the same black hole binary systems that lead to a GW background (GWB) will also coalesce to form the GWM sources. We can use this to derive a simple relation between the strength of the GWB and the GWM event rate.

For any PTA with a data span of 6 yr, the rate bounds will lie between 0.8 and 3.5 yr⁻¹ and the most constraining bound will be around 1 yr⁻¹. At that rate, the PPTA sensitivity must be increased by a factor of 1700 to begin to constrain the existing theory. The current lack of detections is therefore entirely expected. Note that the observational bound will significantly improve with longer and more precise data sets and will move diagonally down and to the left in Fig. 10. In Section 6.5, we describe the data sets that would be needed in order to reach the required sensitivity to make a detection of a GWM event.

We note also that although SMBH binaries are the most commonly predicted GWM emitters detectable by pulsar timing, our bounds on GWM rates do not only apply to memory events associated with SMBH binaries. In this regard, our results put the strongest limits to date on this ‘gravitational wave discovery space’ as described in Cutler et al. (2014).

6.4 Pulsars that contribute to the bound

The current PPTA data set contains observations of 20 pulsars. It is interesting to consider which pulsars actually contribute to the resulting bounds on GWM events. For instance, our available observing time can be used more efficiently if we know which pulsars do not significantly affect our scientific results. To determine this, we eliminated each of the PPTA pulsars from the existing array in turn, and examined the resulting rate bounds obtained for a GWM amplitude of 5× 10⁻¹⁴. In this way, we ranked each pulsar by the increase in the bound (Δλ) that occurred when it was eliminated from the array. The nine pulsars shown in Table 3 are sufficient to obtain a bound 10 per cent larger than the one we actually obtained.

For this paper, we therefore could have achieved almost the same result by only processing these pulsars. However, we do not recommend that the PPTA significantly reduces the number of pulsars observed even though more observing time could then be spent on the ‘best’ pulsars. In order to confirm a detection of either a GWM event or a GWB, it will be necessary to identify the correlated signals using observations of numerous pulsars.

6.5 Expectations for the future

Our data sets will become more sensitive to GWM events with longer and/or improved data sets. The Parkes data are being combined with observations from Northern hemisphere telescopes as part of the International Pulsar Timing Array (IPTA) project (see e.g. Hobbs et al. 2010; Manchester 2013). The initial IPTA data set is currently not finalized, but should lead to the best available data sets for carrying out this research. Currently, 50 pulsars are being timed as part of the IPTA (see Manchester 2013).

In the more distant future, telescopes such as the Five-hundred-meter Aperture Spherical Telescope (FAST), the South African MeerKAT radio telescope and the Square Kilometre Array (SKA) should produce even more sensitive data sets. Exactly how sensitive these data sets will be depends upon to what degree the pulsars will be affected by jitter noise (e.g. Shannon & Cordes 2012; Osłowski et al. 2013) and/or intrinsic timing noise (e.g. Hobbs et al. 2010; Shannon & Cordes 2010). Jitter noise will limit the achievable ToA precision whereas timing noise will limit the long-term stability of pulsar data sets.

Our event rate bounds are currently orders of magnitude away from the expected event rates. We therefore consider what changes in observations would be necessary to constrain astrophysical models of the GWM rate. Longer data sets and more sensitive arrays are the obvious choices. Longer data sets are more sensitive to the signature of a GWM event in the residuals and the chance of an event occurring within the data set will increase as the data span increases. We therefore estimate the effect of extending the time span of the existing PPTA data set using simulations, and we estimate the effect of more sensitive arrays using an analytical model.

First, we expand the existing PPTA data sets to have a 20-year data span assuming no significant change in the number of pulsars observed, their observing cadence or the timing precision achieved. We also assume that the existing red noise models are adequate for extrapolation to the longer time spans. We try an ‘optimistic (in terms of GW detection) extrapolation which assumes that the red noise plateaus at the current data span (i.e. we leave the corner frequency as is) and a ‘pessimistic’ extrapolation in which we assume that the red noise power law continues without a corner frequency for the extended data span.

We form simulated arrival times for each pulsar. We then calculate the sensitivity of the data set to GWM events that occurred at a time corresponding to 15 per cent of the data span. As before, this provides us with bounds on the event rate for a given GWM amplitude. We overplot in Fig. 10, a green, dashed curve that represents the bound obtained from the pessimistic extrapolation of the red noise (overlaid with triangle symbols) and the green, dotted curve that represents the bound obtained from the optimistic extrapolation (overlaid with cross symbols). For comparison, we also plot the bound obtained using the same data span and red noise models as the actual data (green, dot–dashed line; circle symbols). This should agree with the solid curve. The slight difference occurs as (1) these new curves represent GWM events occurring at 15 per cent of the data span, as opposed to at any time within central 80 per cent of the data span and (2) the real data contains effects that are not perfectly modelled by the red and white noise models.

We note that adding 14 extra years does, as expected, make a significant improvement to the bounds (in Fig. 10, the bounds move downwards and to the left). However, the use of the pessimistic red noise extrapolation only slightly reduces the bounds. This is because many of the pulsar data sets that are highly ranked in Table 3 (such as PSRs J1909−3744 and J1744−1134) have no red noise model (they are assumed to have white residuals). It is likely that these pulsars do exhibit timing noise which is currently undetectable and making more realistic predictions of how the bounds will improve with time will require red noise models for those pulsars. Clearly, even with a 20-year observing span, it is unlikely that GWM detection will be made.

Figure 11.

The expected number of detectable sources versus data span for idealized data sets containing: (solid line) 20 pulsars with 1 μs rms timing residuals sampled every 20 d, (dashed line) the same, but for 100 ns rms timing residuals and (dot–dashed line) 50 pulsars are observed with 50 ns white noise and sampled once per week.

Open in new tabDownload slide

These conclusions suggest that it is unlikely that GWM will be detected in the near future. However, it is still important that bursts with memory are searched for as our data sets get longer because of the following.

• As described in this paper, searching for GWM events is reasonably straightforward and the techniques are built into our existing software packages.

• The astrophysical predictions may be incorrect.

• GWM events may occur from sources other than black hole binary systems (for instance, from cosmic strings) and those sources may have a higher event rate.

• As shown in this paper, searching for GWM events leads to an improved understanding of the noise processes within a given data set.

Throughout this paper, we have considered GWM events passing through the solar neighbourhood. A GWM event passing a pulsar will induce a glitch event in that pulsar alone. It would be extremely difficult to prove that any such event in a single pulsar was caused by a GWM, but searching for such events can be used to improve the bounds on the event rates. Whereas in this paper, we can only constrain events within our data span of ∼5 yr, by searching for GWM events in each pulsar independently we would be able to place a constraint over a time-scale of ∼20 × 5 = 100 yr. Such searches have been described by Cordes & Jenet (2012) and Madison et al. (2014).

7 CONCLUSIONS

We have presented straightforward algorithms for detecting the presence of a GWM burst and for bounding the event rate of GWM bursts when no burst is detected. No GWM burst events were detected when applying the algorithm to the first PPTA data release.

The weakest GWM burst we could have detected has an amplitude ∼2 × 10⁻¹⁴. Such a burst would have to occur in a small region of the sky to be detected. By comparison, a burst strong enough that we would have detected it anywhere in the sky would have an amplitude greater than 2.4 × 10⁻¹³. The event rate for such bursts must be less than 0.75 yr⁻¹. These bounds do not significantly constrain models for SMBHB coalescence rates, and, although the bounds will improve with time, the detection of a GWM event may take many years.

The Parkes radio telescope is part of the Australia Telescope, which is funded by the Commonwealth of Australia for operation as a National Facility managed by the CSIRO. This work is supported by National Basic Research Program of China (973 Program 2012CB821800), National Natural Science Foundation of China (nos. 11403086, U1431107, 11173041, 11173042, 11373006 and 11203063) and West Light Foundation of CAS (nos. XBBS201322, XBBS201223 and XBBS201123). GH, LW and YL acknowledge support from the Australian Research Council (ARC) Future Fellowship. VR is a recipient of a John Stocker Postgraduate Scholarship from the Science and Industry Endowment Fund. DRM received support from the NSF/PIRE Grant 0968296. XJZ acknowledges the support of an University Postgraduate Award at UWA.

REFERENCES

APPENDIX A: THE efacEquad PLUGIN

The efac/equad plugin is used to produce values of efacs and equads for a given data set. It relies on the user having identified specific backend and receiver combinations that are expected to have the same efac and equad values. Different combinations are known as different ‘groups’ and are uniquely identified using tempo2 flags. The plugin accepts various command line arguments. The most common usage is

tempo2 -gr efacEquad -f psr.par psr.tim

 -flag <flagID> -plot

where flagID is the flag identifying each group.

The plugin then do the following.

• Switches off all fitting within the parameter file before turning on fits for the pulse frequency (F0) and its first time derivative (F1).

• The red noise is modelled using a constrained, linear interpolation method (known within tempo2 as ifuncs). By default, the linear interpolation is based on a 100 d grid.

• The data are refitted using the linear interpolation, F0 and F1.

• The post-fit residuals are extracted for each group in turn.

• For each group, the reduced-χ² of the fit is determined. If the reduced-χ² < 1, then the uncertainties for that group are decreased by |$\sqrt{\chi ^2_{\rm r}}$|⁠, otherwise the following process is carried out.

• The distribution of normalized residuals (i.e. the residual divided by its uncertainty) is determined and compared (currently using a Kolmogorov–Smirnov test) with a Gaussian distribution with zero mean and unit variance.

• The normalized residuals are re-calculated by modifying the ToA uncertainties with specific efac and equad values. The values are chosen from a grid. For each grid point, the Kolmogorov–Smirnov test probability is recorded and graphically plotted.

• The grid point (efac, equad) that leads to the best match according to the Kolmogorov–Smirnov test is output and used in subsequent processing.

APPENDIX B: THE efac AND equad VALUES

In order to calculate efac and equad values, the observations were first grouped into backend/receiver combinations that we expect to have the same properties. These are listed in Table B1.

The values for each pulsar are as follows:

PSR J0437-4715

T2EQUAD -group 8 0.033

T2EQUAD -group 10a 0.047

T2EFAC -group 10 1.7

T2EQUAD -group 10 0.008

EQUAD 0.080

PSR J0613-0200

T2EQUAD -group 7 0.6

T2EFAC -group 9 1.2

T2EQUAD -group 9 0.2

PSR J0711-6830

T2EFAC -group 4 0.667076

T2EFAC -group 7 0.94107

T2EQUAD -group 9 1

PSR J1022+1001

T2EFAC -group 7 1.3

T2EQUAD -group 7 0.6

T2EQUAD -group 9d1 1.4

T2EFAC -group 9d2 1.7

T2EQUAD -group 9d2 0.5

T2EFAC -group 9 1.9

T2EQUAD -group 9 0.4

EQUAD 0.5

PSR J1024-0719

T2EFAC -group 4 0.956148

T2EFAC -group 7 1.1

T2EQUAD -group 9 0.3

PSR J1045-4509

T2EQUAD -group 4 0.9

T2EFAC -group 7 0.96299

T2EFAC -group 4 1.2

PSR J1600-3053

T2EQUAD -group 7 0.4

T2EFAC -group 4 0.745517

T2EFAC -group 9 1.1

T2EQUAD -group 9 0.3

PSR J1603-7202

T2EQUAD -group 7 1.2

T2EQUAD -group 9 0.7

T2EQUAD -group 9_d2 4

T2EFAC -group 9 1.1

PSR J1643-1224

T2EFAC -group dfb1 0.922936

T2EFAC -group dfbs 0.985818

PSR J1713+0747

T2EQUAD -group 6 0.5

T2EFAC -group 8 1.4

T2EFAC -group 10 1.6

PSR J1730-2304

T2EFAC -group 7 1.2

T2EQUAD -group 9 1

J1732-5049

T2EQUAD -group 9 0.8

T2EFAC -group 7 1.3

T2EQUAD -group 7 0.3

PSR J1744-1134

T2EFAC -group 7 1.3

T2EQUAD -group 7 0.1

T2EFAC -group 4 2.6

T2EFAC -group 9 1.1

T2EQUAD -group 9 0.2

PSR J1857+0943

T2EFAC -group 4 0.741308

T2EFAC -group 7 0.891042

T2EQUAD -group 9 0.1

PSR J1909-3744

EQUAD 0.1

PSR J1939-3744

T2EFAC -group 7 1.3

T2EFAC -group 9 2.3

T2EQUAD -group 9 0.1

T2EFAC -group 13 2.4

T2EQUAD -group 13 0.1

PSR J2124-3358

T2EFAC -group 4 1.1

T2EQUAD -group 4 2.8

T2EFAC -group 3 0.707916

T2EQUAD -group 9 0.6

T2EFAC -group 13 2.6

PSR J2129-5721

T2EFAC -group 7 0.897594

T2EFAC -group 9 0.904696

PSR J2145-0750

T2EQUAD -group 4 0.8

T2EQUAD -group 7 0.5

T2EQUAD -group 9 0.6

EQUAD 0.800000

APPENDIX C: tempo2 AND GWM

As part of this work, various updates have been made to the tempo2 software packages. This updates are available in the current distribution (http://sourceforge.net/projects/tempo2/).

A GWM event can be defined in a pulsar timing model using

GWM_AMP <amp> <fitflag>

GWM_POSITION <ra> <dec>

GWM_EPOCH <mjd>

GWM_PHI <zeta>

where fitflag is usually set to ‘2’ to indicate a global fit. The position (α_g, δ_g) and the position angle (in this paper described as ζ) are given in radians. Instead of using GWM_AMP and GWM_PHI, we have found that the following parameterization is more useful:

GWM_A1 <amp1> <fitflag>

GWM_A2 <amp2> <fitflag>

GWM_POSITION <ra> <dec>

GWM_EPOCH <mjd>

These parameters can be included when fitting or for simulating data sets. In this paper, all fits were carried out using the Cholesky routines using a command similar to

tempo2 -f psr1.par psr1.tim -f psr2.par psr2.tim ...

 -global global.par -dcf model.dat

where the red noise model file was

MODEL T2

PSR J0437-4715

MODEL T2PowerLaw 3 4.34384e-29 0.2

PSR J0613-0200

MODEL T2PowerLaw 5 5.52754e-28 0.4

PSR J1022+1001

MODEL T2PowerLaw 3 1.78292e-27 0.5

PSR J1024-0719

MODEL T2PowerLaw 4 3.96149e-27 0.2

PSR J1045-4509

MODEL T2PowerLaw 2.5 2.0376e-26 0.2

PSR J1600-3053

MODEL T2PowerLaw 4 3.81388e-27 0.2

PSR J1603-7202

MODEL T2PowerLaw 3 2.27845e-27 0.2

PSR J1643-1224

MODEL T2PowerLaw 4 1.5e-26 0.2

PSR J1713+0747

MODEL T2PowerLaw 4 4.13154e-28 0.2

PSR J1824-2452A

MODEL T2PowerLaw 4 1.30232e-25 0.2

PSR J1939+2134

MODEL T2PowerLaw 2.5 2.86815e-27 0.2

PSR J2129-5721

MODEL T2PowerLaw 4 3.56769e-27 0.2