Constraints on ultra-low-frequency gravitational waves with statistics of pulsar spin-down rates

TL;DR

This study uses pulsar spin-down rate statistics to set new upper limits on ultra-low-frequency gravitational waves in a frequency range inaccessible to traditional methods, providing constraints on super-massive black hole binaries.

Contribution

It introduces a novel statistical method analyzing pulsar spin-down rates to constrain ultra-low-frequency GWs, extending the accessible frequency range beyond conventional pulsar timing arrays.

Findings

Upper bounds on GW amplitude derivatives in the Galactic Center and M87 directions.

Derived constraints on GW strain amplitude at frequencies around 1/(100 years).

Implications for the existence of super-massive black hole binaries.

Abstract

We probe ultra-low-frequency gravitational waves (GWs) with statistics of spin-down rates of milli-second pulsars (MSPs) by a method proposed in our prevous work (Yonemaru et al. 2016). The considered frequency range is $10^{-12}{\rm Hz} \lesssim f_{\rm GW} \lesssim 10^{-10}$Hz, which cannot be accessed by the conventional pulsar timing array. The effect of such low-frequency GWs appears as a bias to spin-down rates which has a quadrupole pattern in the sky. We use the skewness of the spin-down rate distribution and the number of MSPs with negative spin-down rates to search for the bias induced by GWs. Applying this method to 149 MSPs selected from the ATNF pulsar catalog, we derive upper bounds on the time derivative of the GW amplitudes of $\dot{h} < 6.2 \times 10^{-18}~{\rm sec}^{-1}$ and $\dot{h} < 8.1 \times 10^{-18}~{\rm sec}^{-1}$ in the directions of the Galactic Center and M87, respectively. Approximating the GW amplitude as $\dot{h} \sim 2 \pi f_{\rm GW} h$, the bounds translate into $h < 3 \times 10^{-9}$ and $h < 4 \times 10^{-9}$, respectively, for $f_{\rm GW} = 1/(100~{\rm yr})$. Finally, we give the implications to possible super-massive black hole binaries at these sites.

Equations

$$r_{GW}(t)$$

$$F^A(\hat{\Omega },\hat{p})$$

$$e_{ij}^{+}(\hat{\Omega })$$

$$e_{ij}^{\times }(\hat{\Omega })$$

$$\Delta h_A(t^{\prime },\hat{\Omega },\theta )$$

$$r_{GW}(t)$$

$$r_{\dot{p}}(t)$$

$$\frac{{\dot{p}}_{\mathrm{obs}}}{p}$$

$$\alpha (\hat{\Omega },\hat{p},\theta )$$

$${\dot{h}}_{+}(\hat{\Omega },\theta )$$

$${\dot{h}}_{\times }(\hat{\Omega },\theta )$$

$$S_{\alpha +(-)}$$

$${\mu }_{+(-)}$$

$${\sigma }_{+(-)}^2$$

$$\Delta S$$

$$f(x)$$

$$y$$

$$\mu$$

$$\sigma$$

$$S$$

$$h_{\mathrm{GC}}$$

$$h_{M87}$$

Full text

††pubyear: 2019††pagerange: Constraints on ultra-low-frequency gravitational waves with statistics of pulsar spin-down rates–References

1 Introduction

Laser Interferometer Gravitational Wave Observatory (LIGO) has succeeded in detecting gravitational waves (GWs) with frequencies $\sim 100~{}{\rm Hz}$ radiated from black hole binaries with masses of about $30M_{\odot}$ (Abbott et al., 2016). With KAGRA (KAGRA Collaboration, 2019) and VIRGO (Acernese et al., 2015), ground-based interferometers will pave the way for gravitational-wave astronomy. In fact, this is the first step toward multi-wavelength gravitational-wave astronomy and lower-frequency GWs are to be probed: space interferometers such as Laser Interferometer Space Antenna (LISA) (Amaro-Seoane et al., 2017) and DECIGO (Kawamura et al., 2011), Pulsar Timing Arrays (PTAs) such as Parkes Pulsar Timing Array (PPTA) (Manchester et al., 2012), European Pulsar Timing Array (EPTA) (Kramer & Champion, 2013), North American Nanohertz Observatory for Gravitational Waves (NANOGrav) (Jenet et al., 2009) and International Pulsar Timing Array (IPTA) (Verbiest et al., 2016), and observations of B-mode polarization of the cosmic microwave background such as Simons Observatory (Ade et al., 2018) and LiteBIRD (Ishino et al., 2016).

GWs radiated from possible super-massive black hole (SMBH) binaries in the Galactic Center (GC) will greatly improve our understanding of gravity theories, astrophysics of SMBHs and environment of the GC region. In the late stage of SMBH binary evolution, GWs with frequencies of $10^{-9}~{}{\rm Hz}\lesssim f_{\rm GW}\lesssim 10^{-6}~{}{\rm Hz}$ are radiated and the range is main target of PTAs. The frequency range of PTA is determined by the observational time span and cadence, and lower frequencies $f_{\rm GW}\lesssim 10^{-10}~{}{\rm Hz}$ are difficult to probe. In fact, the sensitivety is expected to scale as $f_{\rm GW}^{-2}$ toward lower frequencies (Blandford et al., 1984; Moore et al., 2015).

In Yonemaru et al. (2016), we proposed a new detection method for GWs with ultra-low-frequencies of $f_{\rm GW}\lesssim 10^{-10}~{}{\rm Hz}$ (for other methods, see Bertotti et al. (1983); Kopeikin (1997); Potapov et al. (2003)). The method utilizes the fact that the spin-down rate of milli-second pulsars (MSPs) is biased by such GWs, since both give the same quadratic time dependence to the time of arrival of pulses. This effect depends on the relative direction of the GW source and a pulsar. Thus, statistics of the spin-down rate distributions can probe such GWs as we describe later. In Yonemaru et al. (2018), by using simulated 3,000 MSPs, which are expected to be discovered by the Square Kilometre Array (SKA) survey, we estimated the sensitivity of this method in a simple situation, where we assume that MSPs are located uniformly in the sky and the “pulsar term” is neglected. We concluded that GWs with the derivative of amplitude as small as $3\times 10^{-19}$ s*-1* could be detected. Then, in Hisano et al. (2019), we considered a more realistic model of MSP distribution in the Galaxy and took the pulsar term into account in order to obtain more accurate estimates of the sensitivity by extending the analysis of Yonemaru et al. (2018). We found that the sensitivity depends on the direction, polarization and frequency of GWs and becomes worse at low frequencies ($f_{\rm GW}\lesssim 10^{-12}~{}{\rm Hz}$) because of the pulsar term.

This work is the first attempt to apply the above method to real data. We use MSPs selected from the current ATNF pulsar catalog (Manchester et al., 2005) and probe GWs with a frequency of $10^{-12}~{}{\rm Hz}\lesssim f_{\rm GW}\lesssim 10^{-10}~{}{\rm Hz}$. In Section 2, we give a brief summary on our method proposed in Yonemaru et al. (2016). Then, after describing our data set, upper bounds on the derivative of GW amplitude are derived in Section 3. In Section 4, we discuss the implication of the upper bounds to possible SMBH binaries at the Galactic Center and M87. Finally, our results are summarized in Section 5.

2 Detection Principle

Let us first describe the detection method of ultra-low-frequency GWs following Yonemaru et al. (2016). Timing residuals induced by GW are given by (Detweiler, 2010),

[TABLE]

where we denote the direction of pulsar as $\hat{p}$, the propagation direction of GW as $\hat{\Omega}$ and the GW polarization angle as $\theta$. Here, antenna beam pattern $F^{A}(\hat{\Omega},\hat{p},\theta)$ is the geometric factor written by (Anholm et al., 2009),

[TABLE]

where $e^{A}_{ij}(\hat{\Omega})$ $(A=+,\times)$ are the GW polarization tensor given by

[TABLE]

with $\hat{m}$ and $\hat{n}$ being the polarization basis vectors. In Eq.(1), $\Delta h_{A}(t^{\prime},\hat{\Omega},\theta)$ is the difference of geometric perturbation between the earth and the pulsar. This is given by,

[TABLE]

where $t_{p}=t-\tau$ with $\tau=L/c(1+\hat{\Omega}\cdot\hat{p})$ being the pulse propagation time from the pulsar at the distance $L$ to the earth.

In the following, we will discuss GWs with periods much longer than the observational time span, and in this case, the GW amplitude changes linearly with time. At the same time, we use the assumption that the second term ("pulsar term") is a random noise with zero average, which is reasonable when the GW wavelength is much shorter than the typical distance to pulsars. Thus, the GW frequency range we consider here is $10^{-12}{\rm Hz}\lesssim f_{\rm GW}\lesssim 10^{-10}$Hz. For such GWs, we can rewrite Eq.(5) as

[TABLE]

Then, substituting Eq.(5) into Eq.(1), we find the timing residual induced by ultra-low frequency GWs is described by

[TABLE]

This time dependence is the same as timing residual induced by pulsar spin down, which is given by

[TABLE]

where $p$ and $\dot{p}$ are the pulse period and spin-down rate, respectively. Therefore, the influence of ultra-low-frequency GWs is absorbed into the spin-down rate of the pulsar and cannot be identified in the standard analysis of PTA. On one hand, in the presence of ultra-low-frequency GWs, spin-down rate is biased as,

[TABLE]

where $\dot{p}_{\rm obs}$ and $\dot{p}_{0}$ are observed and intrinsic spin-down rates, respectively, and the bias factor $\alpha(\hat{\Omega},\hat{p})$ is given by,

[TABLE]

The bias factor depends on the relative direction between the GW propagation and each pulsar and the spatial pattern is plotted in Fig. 1. Here, $\dot{h}_{+}(\Omega,\theta)$ and $\dot{h}_{\times}(\Omega,\theta)$ are depend on GW polarization $\theta$ and given by

[TABLE]

In Yonemaru et al. (2016), we proposed a method to utilize the spatial pattern of $\alpha(\hat{\Omega},\hat{p},\theta)$ to probe ultra-low-frequency GWs. First, by assuming GW source position and polarization, we divide pulsars into two groups according to the sign of the bias factor, depending on the location of each pulsar. Although GW signals cannot be extracted from individual pulsars, since spin-down rates are biased to positive and negative values in the two groups respectively, it is possible to detect GWs by measuring the systematic difference in the spin-down rate distribution between the two groups. We use the skewness of the spin-down rate distribution to characterize the bias induced by GWs, and convert the difference in the skewness of the two groups to the value of $\dot{h}$. Below, we apply this method to real data and derive constraints on $\dot{h}$ by analyzing the skewness difference. In addition, if the amplitude of GW is too strong, some of pulsars in the region with negative $\alpha(\hat{\Omega},\hat{p},\theta)$ get negative spin-down rate. We derive constraints on the GW amplitude by using the number of pulsars with negative spin-down rates in the real pulsar observation.

3 Results

3.1 Pulsar data

We use data of observed MSPs in ATNF pulsar catalog ver. 1.59. The data set includes 181 MSPs with the measured periods shorter than 30 msec and the time derivatives. Among 181 MSPs, we exclude 30 MSPs in globular clusters since they would be biased significantly by the gravitational potential and complicated dynamics inside the cluster. In addition, two MSPs are removed as outliers: one with a negative spin-down rate $(\dot{p}_{\rm obs}/p=-10^{-20.2}~{}[{\rm sec}^{-1}]$, J1801-3210) and one with an exceptionally large spin-down rate $(\dot{p}_{\rm obs}/p=10^{-11.5}~{}[{\rm sec}^{-1}]$, J0537-6910). Thus, 149 MSPs are used for our analysis below.

Fig. 2 shows the histogram of $\dot{p}_{\rm obs}/p$ of 149 MSPs. Mean, standard deviation, skewness and kurtosis of the distribution are $-17.4$, $0.36$, $1.2$ and $5.9$, respectively, and the deviation from Gaussian distribution was shown to be statistically significant by the Jarque-Bera test (Yonemaru et al., 2018).

3.2 GW search

Given the propagation direction and polarization of GWs, the sky is divided into two areas according to the sign of the bias factor $\alpha(\hat{\Omega},\hat{p},\theta)$. Then MSPs are classified into two groups and skewness of the $\dot{p}_{\rm obs}/p$ distribution is calculated for each group. The skewness in the positive (negative) $\alpha(\hat{\Omega},\hat{p},\theta)$ region is given by

[TABLE]

where $i=1,\cdots,N_{+(-)}$ is the number of MSP in the positive (negative) $\alpha(\hat{\Omega},\hat{p},\theta)$ region, and $\mu_{+(-)}$ and $\sigma^{2}_{+(-)}$ are the mean value and variance of the $\log_{10}\dot{p}/p$ distribution, respectively,

[TABLE]

Then the skewness difference is given by

[TABLE]

Fig. 3 shows the skewness difference as a function of GW polarization angle for the cases where we assume that the GW source is located in the direction of the Galactic Center and M87. It should be noted that polarization angles of [math] and $90$ degree correspond to the same polarization, but the sign of $\alpha$ is inverted so that the sign of skewness difference is also inverted. Maximum values for the case of the Galactic Center and M87 are 0.672 and 0.676, respectively.

In the same way, we calculate the skewness difference for all directions of GW source and all angles of GW polarization. In Fig. 4, at each point of the sky where GW source is assumed to be, we have searched for the largest skewness difference by changing the GW polarization. The skewness difference is mostly smaller than unity and the largest value in the whole sky is $1.07$. The two red regions are antipodes.

In order to study the statistical significance of this value, we perform a series of simulations. First, we make mock data of 149 MSPs in absence of GWs. MSPs are located at the same position as that of the real data. The value of $\dot{p}_{\rm obs}/p$ is randomly allocated to each MSP according to a generalized Gaussian distribution with the same values of mean, standard deviation and skewness as the ones obtained from the catalog. The generalized Gaussian distribution is given by

[TABLE]

where $\phi(y)$ is the standard Gaussian distribution and $y$ is given by

[TABLE]

Here, $\xi$, $\alpha$ and $\kappa$ are the location, scale and shape parameters, respectively, and the mean $\mu$, standard deviation $\sigma$ and skewness $S$ are expressed by these parameters.

[TABLE]

For each realization of mock data, we obtain the skewness difference in the same way as above and search for the maximum varying the position and polarization.

Fig. 5 shows the probability distribution function of the maximum skewness difference in the sky obtained through 10,000 realizations of mock MSP data without GW injection. We find that the distribution extends from $0.6$ to $2.0$ and, as a result, the obtained value $1.07$ is fairly consistent with the statistical fluctuations without GWs.

3.3 GW limits from skewness difference

In the previous subsection, we have shown that the current pulsar data is consistent with the non-existence of GWs within the statistical error. In this subsection, we derive upper bounds on the derivative of the GW amplitude $\dot{h}$. Due to the limited computational power, we focus on two astrophysically important directions of Galactic Center and M87 where the existence of supermassive black hole binaries has been suggested (Yu & Tremaine, 2003; Oka et al., 2016; Yonemaru et al., 2016).

As we saw in the previous subsection, the maximum skewness differences in the direction of the Galactic Center and M87 are 0.672 (at polarization angle of 25 deg) and 0.676 (at polarization angle of 108 deg), respectively. In the presence of GWs with large enough value of $\dot{h}$, the probability of obtaining such small values is low. We place the upper bound on $\dot{h}$ by using the threshold where the probability of having skewness difference less than 0.672 (0.676) is $2\%$ for the direction of the Galactic Center (M87).

In order to evaluate the upper bounds, we make mock data of 149 MSPs in the same way as the previous subsection and inject the bias due to GWs to the simulated intrinsic spin-down rates $\dot{p}_{0}/p$. The polarization angle of GWs is set to be the same as the above: 25 deg for the Galactic Center and 108 deg for M87. In this way, a simulated value of $\dot{p}_{\rm obs}/p$ is allocated to each MSP and then we compute the skewness difference for the data set. It should be noted that MSPs with very small values of intrinsic spin-down rate ($\sim 10^{-18}~{}{\rm sec}$) can have negative values of observed spin-down rate if they are located at an area with negative bias. They are removed from the computation of skewness difference and the total number of used MSPs is slightly smaller than 149.

Fig. 6 shows the probability distribution function of skewness difference in the directions of the Galactic Center and M87 estimated from 10,000 realizations of simulations with and without GWs. Here, we fix GW polarization angles at $25~{}{\rm deg}$ (GC) and $108~{}{\rm deg}$ (M87), which give the largest skewness difference in Fig. 3. We see that, as the value of $\dot{h}$ increases, the probability of having a larger value of skewness difference becomes higher. Comparing them with the observed values (0.672 and 0.676), we obtain an upper bound of $\dot{h}<6.2\times 10^{-18}~{}{\rm sec}^{-1}$ for the Galactic Center and $\dot{h}<8.1\times 10^{-18}~{}{\rm sec}^{-1}$ for M87. The implication of the upper bounds will be discussed in Section. 4.

3.4 GW limits from the number of MSPs with negative spin-down rates

As we mentioned in the previous subsection, when MSPs with very small intrinsic spin-down rates are biased negatively, they can have negative values of observed spin-down rate. The number of such MSPs will increase for stronger GWs. Therefore, the number of MSPs with negative observed spin-down rates could be used for another measure to probe ultra-low-frequency GWs. In the current data set, there is only one MSP with a negative $\dot{p}_{\rm obs}/p$ except ones in globular clusters. Here, we set upper bounds on $\dot{h}$ using the threshold where the probability of having two or more MSPs have negative $\dot{p}_{\rm obs}/p$ is $98\%$.

Fig. 7 shows upper bounds on $\dot{h}$ as a function of GW polarization angle for the Galactic Center and M87. The upper bounds are of order $10^{-17}-10^{-17.4}~{}{\rm sec}^{-1}$ and comparable to those from skewness difference obtained in the previous subsection. GWs from the Galactic Center are slightly well constrained than those from M87 and the dependence on the polarization angle is very weak.

4 Discussion

In the previous section, we have obtained upper bounds on the time derivative of GW amplitudes, rather than the amplitudes themselves. Using an approximation $\dot{h}\sim 2\pi f_{\rm GW}h$ where $f_{\rm GW}$ is the frequency of GW, which is reasonable for most of the periods, our constraints for the Galactic Center and M87 can be approximated as,

[TABLE]

respectively. On the other hand, the recent PTA analyses put upper bounds of $h_{GW}\sim 7.3\times 10^{-15}$ at the frequency of 8 nHz and the bound scales as $f^{-2}_{\rm GW}$ (Aggarwal et al, 2018; Babak et al., 2015). Thus, our constraints are comparable to those of standard PTAs at $f_{\rm GW}\sim 1/30000~{}{\rm years}\sim 10^{-12}~{}{\rm Hz}$ and are better at even lower frequencies. Although our constraints are weaker for $f_{\rm GW}\gtrsim 10^{-12}~{}{\rm Hz}$, they are still valuable as independent constraints. It should be noted that at frequencies lower than $10^{-13}~{}{\rm Hz}$, the pulsar term cannot be treated as random noise and pulsar distances are necessary to account for it (Yonemaru et al., 2018; Hisano et al., 2019). Currently, the distance is not available for most pulsars and we do not consider this frequency range here. Thus, the constraints Eqs. (24) and (25) are applicable for $f_{\rm GW}\gtrsim 10^{-13}~{}{\rm Hz}$.

Let us consider possible SMBH binaries at these cites. A SMBH with mass of $4.0\times 10^{6}M_{\odot}$ is known to reside in the Galactic Center and the possibility of the existence of another SMBH has been discussed (e.g. Oka et al. (2016)). If there exists a SMBH orbiting around the known SMBH, it could be a source of GWs. Assuming the period of the binary motion to be $100~{}{\rm years}$, the upper bound of Eq.(24) translates into a upper bound on the companion mass of $2\times 10^{16}~{}M_{\odot}$.

Concerning the M87, the mass of a SMBH is estimated to be $6.6\times 10^{9}M_{\odot}$ at the center of M87. it has been indicated to have secondary SMBH and would be binary and GWs from such a potential pc-scale SMBH binary has been discussed (Betcheldor et al., 2010; Yonemaru et al., 2016), while constraints on the amplitude of GWs emitted by a milli-pc scale SMBH binary in the PTA frequency bands has been already studied (Schutz & Ma, 2016). Assuming the orbital period to be $100~{}{\rm years}$, the upper bound of Eq.(25) results in a upper bound on the companion mass of $4\times 10^{16}M_{\odot}$.

In the above considerations, binaries are assumed to have circular orbits and inclination is zero degree (face-on). In our previous works (Yonemaru et al., 2018; Hisano et al., 2019), we estimated future constraints on ultra-low-frequency GWs with 3,000 MSPs, which are expected to be found by the SKA2. There, we found that GWs with $\dot{h}$ as small as $\sim 3\times 10^{-19}~{}{\rm sec}^{-1}$ can be detected and the second SMBH mass as small as $3\times 10^{14}~{}M_{\odot}$ could be probed in the case of circular orbits and zero inclination. In fact, the GW amplitude is sensitive to the eccentricity and the phase of the binary motion.

Our simulation is based on the assumption that the distribution of spin-down rate in the sky is isotropic in the absense of GWs. However, considering the evolution of MSPs, it may not be the case. Neutron stars are mostly produced in Galactic plane and often have large peculiar velocities due to the kick at supernovae (Hansen & Phinney, 1997). Since MSPs with small values of $\dot{p}/p$ have long characteristic age ($p/\dot{p}$), they may tend to be located far from the birth place outside the Galactic plane. Thus, there is a possibility that MSPs with large (small) $\dot{p}/p$ are populated outside (inside) of the Galactic plane, which induces the anisotropy of $\dot{p}/p$ distribution in the sky and results in systematics in our method.

In the top panel of Fig. 8, the position of 149 MSPs in the galactic coordinates is shown with the indication of spin-down rates. In the bottom panel, we show the histogram of spin-down rates of MSPs within Galactic plane ($|b|<10~{}{\rm deg}$) and outside ($|b|>10~{}{\rm deg}$), separately. The mean and standard deviation of the histograms are -17.4 and 0.5 for $|b|<10~{}{\rm deg}$, and -17.5 and 0.4 for $|b|>10~{}{\rm deg}$, respectively. Thus, significant difference is not found between the two histograms and the systematics is considered to be negligible. Here it should be noted that MSPs with $|b|>10~{}{\rm deg}$ may reside within the Galactic disk but it cannot be known without information on the distance.

Finally, we note that the observed spin-down rates could be biased by other factors than GWs such as the Shklovskii effect (Shklovskii, 1970), the Galactic differential rotation (Damour & Taylor, 1991; Rong & Tan, 1999) and acceleration toward the Galactic disk (Nice & Taylor, 1995). Although the biases from the Galactic differential rotation and acceleration toward the disk would have spatial correlations in the sky, these effects are less significant ($\Delta(\dot{p}/p)$ $\leq 10^{-19}$ for MSPs at $\leq$ 10 kpc) (Nice & Taylor, 1995) and will be removed if the distance to MSPs is measured precisely in the future.

5 Summary

In this paper, we have placed constraints on GWs from a single source with ultra-low frequencies ($10^{-12}~{}{\rm Hz}\lesssim f_{\rm GW}\lesssim 10^{-10}~{}{\rm Hz}$) by applying a method proposed in Yonemaru et al. (2016) to observed milli-second pulsars (MSPs). This method is based on the statistics of spin-down rate distribution, where the skewness difference between two MSP groups divided according to the position in the sky, is used to find a bias induced by GWs. We selected 149 MSPs from the ATNF pulsar catalog and calculated the skewness difference. By comparing with mock MSP data, we have shown that the current MSP data is consistent with no GWs from any direction of the sky. Furthermore, we have derived upper bounds on the time derivative of the GW amplitude $\dot{h}<6.2\times 10^{-18}~{}{\rm sec}^{-1}$ and $\dot{h}<8.1\times 10^{-18}~{}{\rm sec}^{-1}$ for the Galactic Center and M87, respectively. Consistent bounds were derived from the number of MSPs with negative spin-down rates. Approximating the GW amplitude as $\dot{h}\sim 2\pi f_{\rm GW}h$, the bounds respectively translate into $h<3\times 10^{-9}$ and $h<4\times 10^{-9}$ for $f_{\rm GW}=1/(100~{}{\rm yr})$. The constraints will improve by more than one order of magnitude with 3,000 MSPs in the SKA era (Yonemaru et al., 2018; Hisano et al., 2019).

Acknowledgements

We thank George Hobbs for useful discussion. NY was financially supported by the Grant-in-Aid from the Overseas Challenge Program for Young Researchers of JSPS. SK is partially supported by the Grant-in-Aid for Scientific Research from JSPS, Grant Number 17K14282, and by the Career Development Project for Researchers of Allied Universities. KT is partially supported by JSPS KAKENHI Grant Numbers JP15H05896, JP16H05999 and JP17H01110, and Bilateral Joint Research Projects of JSPS. 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 CSIRO. The ATNF pulsar catalogue at http://www.atnf.csiro.au/people/pulsar/psrcat/ was used for this work. We also thank the referee for a careful reading of the manuscript.