Timing analysis for 20 millisecond pulsars in the Parkes Pulsar Timing Array

Abstract

1 INTRODUCTION

The Parkes Pulsar Timing Array (PPTA; Manchester et al. 2013), like its North American (Demorest et al. 2013) and European (Kramer & Champion 2013) counterparts, is a programme in which an array of millisecond pulsars (MSPs) is observed regularly over many years. The times of arrival (ToAs) of pulses from MSPs are highly predictable using timing models that describe the spin evolution and astrometric properties of the pulsar and any companions, as well as taking into account the motion of Earth and pulse propagation through curved space–time and the interstellar medium (ISM). The parameters of the timing model are determined, or improved, by a least-squares fit of the model to the ToAs. The differences between the measured and predicted ToAs (the ‘timing residuals’) after this fit contain the measurement error, stochastic fluctuations in the apparent pulsar rotation rate (known as timing noise; Shannon & Cordes 2010), and other unmodelled effects such as those of gravitational waves (GWs), errors in the assumed time standard, and errors in the Solar system ephemeris. Most of these effects are stronger at lower frequencies, i.e. they have a ‘red’ power spectrum.

The main goal of a pulsar timing array (PTA) is to search for and eventually detect and study nanohertz frequency GWs (e.g. van Haasteren et al. 2011; Demorest et al. 2013; Zhu et al. 2014; Wang et al. 2015), but there are many secondary objectives such as testing general relativity (GR; e.g. Freire et al. 2012; Zhu et al. 2015), constraining common models of supermassive black hole and galaxy formation (e.g. Shannon et al. 2015), measuring planetary masses (e.g. Champion et al. 2010), studying the ISM (e.g. You et al. 2007; Keith et al. 2013, hereafter K13), developing pulsar-based time standards (e.g. Hobbs et al. 2012), and precise measurements of properties of the pulsars themselves (e.g. Verbiest et al. 2008, hereafter V08). The latter includes for example the much improved distance and mass measurements for PSR J0437−4715 presented in this paper, which will be important for future limits to changes of Newton's gravitational constant, and the Neutron Star Interior Composition Explorer (NICER) mission that will attempt to measure its radius (Gendreau et al. 2012). The new and improved distance measurements for pulsars presented in this paper are also useful for future galactic electron density models (cf. Cordes & Lazio 2002).

It has long been known that least-squares fitting for the parameters of the timing model can be biased and can underestimate the uncertainties on the parameters when the residuals contain significant red noise. To account for this effect, V08 used Monte Carlo simulations for PSR J0437−4715 to determine parameter uncertainties. Verbiest et al. (2009, hereafter V09), however, prewhitened the residuals for the three pulsars with the most red noise in their sample by fitting harmonically related sine/cosine pairs (Hobbs, Edwards & Manchester 2006). Coles et al. (2011, hereafter C11) demonstrated that this method can result in biased parameter measurements and underestimated uncertainties. For all pulsars, V09 then doubled the formal uncertainties obtained from the fit.

For our work, we use an extension of the ‘Cholesky’ algorithm developed by C11 and implemented in the timing software package tempo2 (Edwards, Hobbs & Manchester 2006; Hobbs et al. 2006). For the results presented in C11, the red noise was modelled as wide-sense stationary, i.e. having a single power spectrum. However, the algorithm only requires that the red noise be described by a covariance matrix. In our data set, the red noise is not stationary because the earlier data contain uncorrected fluctuations in the dispersion measure (DM; the column density of electrons along the line of sight to the pulsar) of the ISM, while the later data do not. Accordingly, we use, and describe below, a modification of the algorithm in C11 that we refer to as the ‘split-Cholesky’ algorithm, which allows for two different red-noise models in the data set.

In Section 2, we describe the observations and methodology for determining white-noise parameters. In Section 3, we describe the parameters of the timing model. In Section 4, we describe the new split-Cholesky algorithm, modelling of the DM variations, and present parameters describing the red-noise and DM-noise models. In Section 5, for each pulsar, we present the new timing model parameter values. In Section 6, we present simulations of the split-Cholesky algorithm and compare the method to alternate Bayesian pulsar timing analysis algorithms (e.g. van Haasteren et al. 2011; van Haasteren & Levin 2013; Lentati et al. 2014a), and derive a precise pulsar distance for PSR J0437−4715.

2 OBSERVATIONS

The observations used here were published as the extended first PPTA data release (DR1E) by Manchester et al. (2013). All observations were taken using the Parkes 64 m radio telescope. The data set includes observations at three observing bands (with approximate centre wavelengths of 10, 20, and 50 cm) from the PPTA project that commenced in 2005 (these observations alone are referred to as data release one, DR1), along with observations prior to 2005 (in the 20 cm band only) from previous observing programmes. The earliest data were obtained from a timing programme that commenced during the Parkes 70 cm survey (Bailes et al. 1994) and were published by Bell et al. (1997) and Toscano et al. (1999). The sample of pulsars was increased by MSPs discovered during the Swinburne intermediate-latitude survey (Edwards et al. 2001) and elsewhere. Updated timing solutions were published by Hotan, Bailes & Ord (2006) and Ord et al. (2006). Throughout this paper, we refer to the archival 20 cm observations taken prior to the PPTA as the ‘early’ data and the multi-band PPTA observations as ‘recent’ data.

An intensive observing campaign was used to study PSR J0437−4715 in detail. Results were published in van Straten et al. (2001). V08 (for PSR J0437−4715) and V09 (for the other 19 pulsars in the PPTA) combined the earlier data with the initial PPTA data to determine timing ephemerides. Here we use the extra ∼3 years of data provided by Manchester et al. (2013) to improve on the results of V08 and V09. Along with the extra data span, our new data set provides significantly improved observing cadence and, for the recent data, the ability to remove the effects of DM variations (ΔDM) more precisely than previously possible. Most of the raw observation files used in this analysis are available from the Parkes pulsar data archive (Hobbs et al. 2011).

• Estimate the red noise by fitting a smooth model to the residuals. The default smooth model is a linear interpolation through a set of samples on 100 d intervals. The smooth model (red-noise estimate) is then subtracted from the residuals, leaving only the white noise.

• Divide the data into groups based on observing systems that are expected to have the same efac and equad values. For example, our data set includes data taken using different ‘backend’ instruments, some of which have identical firmware, bandwidths, etc., and are therefore expected to share the same noise properties.

• For a given group selected from the whitened data set, the reduced-χ² value (⁠|$\chi ^2_{\rm r}$|⁠) is calculated. If |$\chi ^2_{\rm r} < 1,$| then |${\rm \small {EFAC}} = \sqrt{\chi ^2_{\rm r}}$| and |${\rm \small {EQUAD}} = 0$|⁠. Note that this is the only means to reduce a ToA uncertainty as an equad will always increase the uncertainty.

• If |$\chi ^2_{\rm r} > 1,$| then the normalized residuals (r_i/σ_i) are determined for a particular grid of efac and equad values.

• For each grid position, we determine the probability that the normalized residuals are drawn from a Gaussian distribution (using a Kolmogorov–Smirnov test) to determine optimal efac and/or equad.

• By default, we do not include the efacs and equads for a group that consists of less than 10 ToAs. They are added if it is necessary to produce normally distributed residuals.

3 THE TIMING MODEL

The timing model describes the spin, astrometric, and orbital properties of a pulsar, the ISM along the line of sight, and requires the use of a terrestrial time standard, Solar system ephemeris, and solar wind model. For this work, we use the DE421 Jet Propulsion Laboratory (JPL) Solar system ephemeris and as a time reference use TT(BIPM2013) with units of time in Barycentric Coordinate Time (TCB). We use the default model within tempo2 to account for DM variations caused by the solar wind (wind density at 1 au: 4 cm⁻³; Edwards et al. 2006). The pulsars have been observed over many years with various different ‘backend’ instruments. These instruments have different time offsets which we also measure as part of the usual timing fit (with the split-Cholesky method).

The 13 pulsars in our sample that are in binary systems each have a white dwarf companion. For such systems, we do not expect any time dependences for the orbital parameters caused by mass-loss or spin–orbit coupling. Because of the relatively low mass of the companion stars and relatively long orbital periods, relativistic effects are small in such systems. However, V08 do report a detection of the advance of the longitude of periastron, |$\dot{\omega }$|⁠, for PSR J0437−4715 that is consistent with that predicted from GR, where the component masses were derived from the Shapiro delay measurement and mass function (Thorsett & Chakrabarty 1999). van Straten (2013) also measures |$\dot{\omega }$| for PSR J1022−1001 that is consistent with GR.

For binary pulsars in orbits with small eccentricities, the longitude and epoch of periastron are not well defined and are highly correlated. The ELL1 model (Lange et al. 2001) is used to describe such systems, since it uses a small-eccentricity approximation to avoid this high correlation. For PSRs J1022+1001, J1600−3053, and J1643−1224, this approximation is not valid and so we instead use the DD (Damour & Deruelle 1986) model to describe the binary orbit.

Observed changes in the orbital parameters can be caused by kinematic effects. Changes in the apparent viewing geometry of the orbit caused by proper motion can lead to an apparent time derivative of the projected semi-major axis of the orbit, |$\dot{x}$|⁠, and/or |$\dot{\omega }$| (Kopeikin 1996). For some pulsars, this kinematic |$\dot{x}$| or |$\dot{\omega }$| may be detected individually. However, if both are well determined, we instead parametrize the effect with Ω and i, which are the longitude of the ascending node and inclination angle of the orbit, respectively (van Straten & Bailes 2003). This parametrization of the orbit describes the annual orbital parallax and is implemented through the T2 model (Edwards et al. 2006), which we use for PSRs J0437−4715, J1713+0747, and J1909−3744. For PSRs J1022+1001, J1600−3053, J1603−7202, J1643−1224, and J2145−0750, we use the |$\dot{x}$| measurement to place an upper limit on the inclination angle of the orbit (using |$\tan {i}\le x\mu /\dot{x}$|⁠, where μ is the total proper motion), as was done by Sandhu et al. (1997).

3.1 Choosing parameters to include in the model

For all pulsars, we start with the timing models presented by V09 (and, for PSR J0437−4715, the model from V08). For solitary pulsars, model parameters include the spin (ν, |$\dot{\nu }$|⁠), astrometric (position, proper motion, and parallax), and ISM parameters (Section 4.1). In all cases, we fit for the spin, position (right ascension α and declination δ), and DM variation parameters (from Manchester et al. 2013). These parameters are all measured at a reference epoch of MJD 54500. In a few cases (in particular where the pulsar is either close to or almost perpendicular to the ecliptic plane, respectively), we cannot obtain a significant measurement of the proper motion in declination (μ_δ) or the parallax (π).

• Remove non-essential parameters (proper motion in right ascension μ_α and declination μ_δ, and parallax π) from the timing model (if present) and note the χ² value of residuals.

• Fit for each parameter separately and note the new χ² value in each case.

• Include the parameter that results in the lowest χ² value permanently into the timing model if this parameter also satisfies the AIC.

• Repeat with remaining parameters until either all parameters are included in the timing model, or all remaining parameters fail to improve the timing solution, as determined by the AIC.

This procedure is also applied to pulsars in binary systems. All non-essential binary parameters (post-Keplerian) are initially removed and the AIC is used to find which should be included. If a Shapiro delay (Shapiro 1964) is detectable for the binary system, we parametrize this with the companion mass, M_c, and sine of inclination angle of the orbit, sin i, except in the cases of pulsars described with the T2 binary model, for which we link sin i with the measured inclination angle from the Kopeikin terms (Kopeikin 1996).

4 PARAMETER MEASUREMENT IN THE PRESENCE OF NON-STATIONARY RED NOISE

In the original implementation of the C11 algorithm, the red noise was characterized by a power spectrum, from which a covariance function was estimated and finally a covariance matrix was constructed. This was satisfactory for analysis of the DR1 data set. However, DR1E data contain uncorrected DM variations in the 20 cm residuals prior to multi-band observations and thus have additional red noise. We cannot produce a single power-law model for the entire data set because of this extra red noise which may dominate the total red noise in a subset of the data. Instead, we produce two separate models to describe the two sources of red noise in the residuals. One model describes the frequency-independent noise present throughout the data set using a power law, while the other describes the additional DM noise present only in the early data (Section 4.1). We have updated the implementation of the algorithm to synthesize a covariance matrix for DR1E observations using these two red-noise models, which we refer to as the split-Cholesky algorithm.

4.1 Modelling the DM variations

With the advent of the PPTA in 2005, regular observations occurred in multiple observing bands (10, 20, and 50 cm). For these data, we are, in principle, able to obtain a measurement of the DM because the group delay is ∝λ². However, as the pulse profiles of MSPs in our sample evolve significantly with frequency (Dai et al. 2015), the effect of an absolute DM on the residuals is coupled with the frequency evolution of the pulse profile and is therefore difficult to determine. We are however able to measure changes in DM, which we refer to as ΔDM(t), using interband measurements with the required accuracy (≈1: 10⁵). Residual errors in these corrections are an important source of noise in the PPTA (and all other PTAs). Such time series (obtained from the same observations we analyse here but with different efacs and equads) were analysed by K13.

The annual variation, linear trend, and covariance function parameters that we have used to construct a DM model for each pulsar are given in Table 1. In some cases, no dDM/dt is apparent in the data, but the ΔDM(t) noise is nevertheless well modelled by the covariance function since it is small. For all pulsars with the exception of PSR J0437−4715 (which is very well modelled by a Kolmogorov power law), we find that these covariance function models successfully whiten the residuals and we therefore include them in our combined red-noise models. For PSR J0437−4715, we use the Kolmogorov power law presented in K13 to model the DM noise present in the early data.

For PSR J1603−7202, there is an extreme scattering event (ESE; Fiedler et al. 1987), lasting ∼250 d, which was reported by K13 and is described in detail in Coles et al. (2015). The ESE dominates the shape of the non-DM-corrected 20 cm residuals in that region. We examined the early 20 cm data searching for ESEs comparable with the one reported by K13 and found none. We obtained the covariance model for ΔDM(t) by linearly interpolating across the ESE before computing the covariance function.

Maitia, Lestrade & Cognard (2003) reported on a 3-year-long ESE (centred on the year 1998) detected in the direction of PSR J1643−1224 by studying flux variability of the pulsar using observations undertaken at the Nancay observatory. Unfortunately, we have poor data during this time and consequently do not find evidence for such an event. For PSR J1713+0747, we see a peculiar ‘dropout’ in the ΔDM(t) that is probably related to an ESE-like structure. We therefore computed the covariance of ΔDM(t) for PSR J1713+0747 by interpolating across this dropout.

The deterministic terms in the model for each pulsar include an absolute DM measurement (which we hold constant because it is covariant with pulse profile evolution; Dai et al. 2015), and if required, the additional terms dDM/dt, A and B, or ΔDM(t) measurements taken at a spacing of Δt_DM days to remove any residual red noise. The sampling interval Δt_DM is ideally the widest spacing required to remove the red noise since the measurements add white noise to the residuals. We start with the Δt_DM published with the data set in Manchester et al. (2013). We sample more frequently if there is residual red noise because the published Δt_DM were selected to minimize the rms residual in the best band for each pulsar, and not to absorb the most red noise. However, adding the deterministic DM parameters often allows us to reduce the ΔDM(t) measurement cadence relative to that published in Manchester et al. (2013).

4.2 Modelling the red timing noise

The analytical model (equation 4) for the frequency-independent red noise (e.g. timing noise) must be estimated with data that have the frequency-dependent red-noise contributions from DM variations already removed. We therefore use the DM-corrected DR1 data and fit the analytical model (equation 4) to the power spectrum. In most cases, we find that f_c = 1/T₁ fits the data reasonably well, where T₁ is the length of the DR1 data. For these cases, we assume that the model will fit the entire DR1E with f_c = 1/T_E, where T_E is the length of the DR1E data set, and we scale P₀ to this new f_c accordingly.

In many cases, extrapolating the timing noise model by adjusting f_c will underestimate the noise at f ≈ 1/T_E because fitting for the spin frequency and its derivative removes much of the low-frequency power at f ≈ 1/T₁. For pulsars where the timing noise dominates the DM noise, we create a single spectral model for the entire DR1E span directly. If this model ultimately whitens the residuals adequately, we do not need to estimate the covariance matrix of the DM variations. If however the timing noise model does not extrapolate well and the DM noise is too significant to ignore, we construct a red-noise model using the DR1E data, but reduce the amplitude, P₀, to account for the known DM noise in the data. The amplitude needs to be reduced such that the residuals are sufficiently whitened by the red-noise model.

For PSRs J0613−0200, J1600−3053, and J2145−0750, we find that the red-noise model from the DR1 data alone does not extrapolate well over the entire data set since the power spectrum of the DR1E data is observed to turn-over at f_c > 1/T_E. This was not obvious when analysing the DR1 data alone. We observe a similar turn-over for PSR J1909−3744; however, the data set used for this pulsar includes multi-band observations across the entire span and so is largely DM corrected (see Section 5.16). For this pulsar, we therefore model the power spectrum of the entire data set and use this with the original Cholesky algorithm.

The parameters describing the red timing noise models for each pulsar are given in Table 2. The DM models that are used with these red-noise models are described in Section 4.1, with parameters in Table 1.

4.3 How do we know when our models are optimal?

The split-Cholesky algorithm produces the best linear unbiased estimators if, and only if, the residuals after whitening are white and normally distributed. If this is not the case, one or more of the red-noise model, efacs, or equads are incorrect. We utilize two tests to check that the final residuals are what we require to have confidence in the parameter measurements and uncertainties.

The first of these is an Anderson–Darling (AD) test for normality (Anderson & Darling 1954). We apply this to our whitened, normalized residuals to determine whether they are consistent with a normal distribution with μ = 0 and σ = 1. The result is the modified AD statistic, A*². This is used to test the hypothesis that the residuals obey the described normal distribution. The hypothesis is rejected if A*² > 2.492 with 5 per cent significance (Stephens 1974), since the expected distribution function for the normalized residuals is known. For example, panel (a) of Fig. 1 shows the post-fit whitened and normalized residuals for PSR J1713+0747 and panel (c) shows the cumulative distribution of these residuals. The modified AD statistic, A*² = 1.03, indicates that the normalized residuals are consistent with a standard normal distribution.

Figure 1.

Open in new tabDownload slide

(a) Whitened and normalized post-fit residuals for PSR J1713+0747. (b) Power spectra of whitened and normalized post-fit residuals. Dotted line marks f = 1 yr⁻¹ and dashed lines indicate expected mean and variance for the spectrum. The pseudo-frequency is determined by converting the whitened components to a pseudo-time series using the ToAs of the unwhitened residuals. (c) Cumulative distribution of whitened and normalized post-fit residuals (solid line) with expected distribution based on normal distribution with zero mean and unit variance (dashed line). Modified AD statistic for this distribution with the expected distribution is labelled. (d) Power spectra of post-fit residuals. Dotted line marks f = 1 yr⁻¹ and solid flat line is an estimate of the white-noise level. Dashed line is the power-law model of the frequency-independent timing noise.

We then test for ‘whiteness’ by inspecting the power spectrum of the whitened, normalized residuals (see also C11). Panel (b) of Fig. 1 shows the Lomb–Scargle periodogram of these residuals for PSR J1713+0747. We compute this power spectrum by first converting the whitened components in panel (a) to a time series using the ToAs of the unwhitened residuals. The frequency axis of the resulting power spectrum is not well defined, and we term it the pseudo-frequency. However, the whitening process provides a diagonally dominant whitening matrix such that low frequencies in the unwhitened residuals translate to low frequencies in the pseudo-time series of the whitened components. The pseudo-frequency power spectrum can therefore be a useful test of the red-noise model used in the whitening process, since it must be flat. The power spectrum of the unwhitened residuals is shown in panel (d) of Fig. 1 with the timing noise model used for this pulsar. If the test fails for any pulsar, we update the red-noise model and re-fit the timing model until the residuals are successfully whitened as required.

5 RESULTS

In this section, we present our final timing solutions for each pulsar, with the post-fit residuals given in Fig. 2. We compare our measured parameters with those in the literature, including very long baseline interferometry (VLBI) measurements and distances derived from the pulsar DMs using the Taylor & Cordes (1993, hereafter TC93) and NE2001 (Cordes & Lazio 2002) galactic free electron distribution models (giving distances accurate to approximately 25 and 20 per cent, respectively).

Figure 2.

Open in new tabDownload slide

Final post-fit residuals for each of the pulsars in our sample. The vertical range of each subplot is given below the pulsar name.

Much of the comparison in this section will be with V08 and V09, which used a subset of our data set and is therefore not independent. However, V08 and V09 used a different Solar system ephemeris, DE405, and time standard, TT(TAI), to our analysis, which results in apparent changes to some parameters. Use of the DE405 ephemeris in particular induces significantly different position parameters compared to the newer DE421 ephemeris, while the use of TT(TAI) changes the apparent spin frequency and its derivative. For each pulsar, we find that inconsistencies in the position parameters are explained by the use of this different ephemeris.

Although K13 measured DM variations for all PPTA pulsars and determined an optimal sampling time, Manchester et al. (2013) did not publish any ΔDM measurements for seven pulsars. This is because ΔDM values were only published by Manchester et al. (2013) if their inclusion improved the rms residuals. For our work, we must account for all red noise in the residuals since this is required for the split-Cholesky algorithm (this process is described in Section 4). The testing of red- and white-noise models for each pulsar is described in Section 4.3, where we present the analysis of PSR J1713+0747 as an example. We do not present the analysis of the other pulsars, with the exception of PSR J0437−4715, which shows some variation that we discuss below in Section 5.1.

For each pulsar, we also compare the various distance measurements available in the literature to check that our new distance measurements from π or |$\dot{P_b}$| are consistent with these. For each parallax measurement, we calculate the Lutz–Kelker (Lutz & Kelker 1973) corrected parallax value and the corresponding corrected distance using the method of Verbiest et al. (2012) and with the mean flux of each pulsar given in Manchester et al. (2013). The parallax measurements for PSRs J1024−0719, J1045−4509, J1603−7202, and J1857+0943 have high fractional uncertainties (parallax measured with less than 3σ confidence) and therefore the bias-corrected values for these closely resemble the prior distributions used. The measurements however are important for placing an upper bound on the parallax value, or similarly, a lower bound on the distance measurement. A table of these values, and the distances from |$\dot{P_b}$| and DM are given in Table 3.

For each pulsar, we derive the characteristic age, surface magnetic field strength, and energy loss rate from the spin-down. Where available, we use the pulsar mass function (Thorsett & Chakrabarty 1999) and a precise Shapiro delay measurement to calculate a pulsar mass, and use the kinematic |$\dot{x}$| measurements to place a limit on the inclination angle of the orbit. These derived parameters are given in Appendix A with the tables of parameters.

The red-noise models, ToA files, and parameter files for each pulsar are available from the Commonwealth Scientific and Industrial Research Organisation (CSIRO) Data Access Portal.²

The final parameters for the seven solitary pulsars J0711−6830, J1024−0719, J1730−2304, J1744−1134, J1824−2452A, J1939+2134, and J2124−3358 are given in Table A1. The binary pulsars are separated by the binary model used to describe their orbit. Parameters for the small-eccentricity pulsars described by the ELL1 model, PSRs J0613−0200, J1045−4509, J1603−7202, J1732−5049, J1857+0943, J2129−5721, and J2145−0750, are presented in Table A2; DD model pulsars J1022+1001, J1600−3053, and J1643−1224 are presented in Table A3; and T2 model pulsars J0437−4715, J1713+0747, and J1909−3744 are presented in Table A4.

5.1 PSR J0437−4715

PSR J0437−4715 is the closest MSP currently known and the brightest at radio wavelengths. Van Straten et al. (2001) presented a timing solution including the full three-dimensional geometry of the binary orbit. An updated model was presented by V08. This included a precise distance estimate derived from an orbital period-derivative measurement.

In Fig. 3, we show the two components of the red-noise model with the power spectrum of the post-fit residuals (panel d). As shown in panel (b), the model successfully whitens the residuals. Significant uncorrected DM noise is present in the early data. Therefore, as expected, the timing noise model underestimates the total noise, whereas the DM-noise model alone overestimates the total noise (since it does not apply to the entire data set). There is excess noise at all frequencies (the mean of the power spectrum of normalized residuals is >1). This could result, for example, by additional uncorrected short-time-scale correlated noise in the residuals. The normalized residuals do not pass the AD test (A*² = 6.75). This may be because of instrumental effects or because of pulse jitter (Shannon et al. 2014). Non-Gaussianity has been detected previously for PSR J0437−4715 by Lentati, Hobson & Alexander (2014b). However, it was shown that, at this level, the non-Gaussianity and high mean spectral power do not significantly affect the parameter measurements or uncertainties.

Figure 3.

Open in new tabDownload slide

(a) Whitened and normalized post-fit residuals for PSR J0437−4715. (b) Power spectra of whitened and normalized post-fit residuals. Dotted line marks f = 1 yr⁻¹ and dashed lines indicate expected mean and variance for the spectrum. The pseudo-frequency is determined by converting the whitened components to a pseudo-time series using the ToAs of the unwhitened residuals. (c) Cumulative distribution of whitened and normalized post-fit residuals (solid line) with expected distribution based on normal distribution with zero mean and unit variance (dashed line). Modified AD statistic for this distribution with the expected distribution is labelled. (d) Power spectra of post-fit residuals. Dotted line marks f = 1 yr⁻¹ and solid flat line is an estimate of the white-noise level. The timing noise model applies over entire data set in combined noise model, while the Kolmogorov DM-noise model only applies to residuals prior to MJD 53430.

Using our new, precise measurement of |$\dot{P_b}$|⁠, we can calculate an improved distance to PSR J0437−4715 of D = 156.79 ± 0.25 pc. We discuss this measurement in detail in Section 6.2. As shown in Table 3, this measurement is consistent with independent distances, including our parallax distance measurement of D = 157 ± 2 pc. V08 measured a high pulsar mass for this pulsar, of M_p = 1.76 ± 0.2 M_⊙. Our improved measurement of the Shapiro delay reduces the uncertainty on the pulsar mass and we find M_p = 1.44 ± 0.07 M_⊙, significantly smaller at the 1.5σ level. This improved mass measurement will be important for the NICER mission which will attempt to measure the neutron star radius, probing the neutron star equation of state.

Kinematic contributions to the measured |$\dot{\omega }$| are included in the timing model through the measurement of Kopeikin terms. We therefore expect that our |$\dot{\omega }$| measurement is solely due to the effects of GR, as was reported in V08. Under this assumption, the reported |$\dot{\omega }$| corresponds to a combined pulsar and companion mass of M_p + M_c = 1.44 ± 0.2 M_⊙, which is just consistent with the measured masses at the 1σ level. In analysis of future data sets for this pulsar, the improving |$\dot{\omega }$| measurement will be able to be used in combination with the Shapiro delay measurement to further constrain the pulsar mass.

5.2 PSR J0613−0200

The only new parameters included in our PSR J0613−0200 timing model are those describing the sinusoidal annual variations in DM that were presented in Section 4.1.

5.3 PSR J0711−6830

The pulsar timing model for this pulsar contains the same parameters as in V09. As expected, our measurements are more precise than previous work. We have not been able to obtain a parallax measurement for this pulsar because of its proximity to the ecliptic pole.

5.4 PSR J1022+1001

This pulsar has the smallest ecliptic latitude (−0.06 deg) of all pulsars in the sample. Observations when the line of sight to the pulsar passes close to the Sun have been used by You et al. (2007) to study the solar corona. However, these observations are generally removed for high-precision timing applications. For our timing solutions, we removed the ToAs that were obtained when the line of sight to the pulsar passed within 5° of the Sun. The removal of these residuals produces normalized residuals that are consistent with zero mean and unit variance.

van Straten (2013) analysed 7.2 years of data for this pulsar from the Parkes radio telescope using a new method of polarimetric calibration to improve timing precision. The DR1 and DR1E data sets did not include this new calibration procedure. With our data set, we have detected a significant Shapiro delay using the traditional parameters M_c and sin i. However, while we found sin i = 0.69 ± 0.18, we did not measure a significant companion mass, M_c = 2.2 ± 2.4 M_⊙. Our parameters are consistent with those of van Straten (2013).

We measured the orbital period derivative, |$\dot{P}_b = (5.5 \pm 2.3)\times 10^{-13}$|⁠, for the first time [this parameter was not measured by van Straten (2013). Assuming a pulsar mass of M_p = 1.4 M_⊙ and companion mass of M_c = 0.2 M_⊙, the expected |$\dot{P}_b$| contribution from quadrupolar GW emission is |$\dot{P}^{\rm GR}_{b}=-1.7\times 10^{-16}$|⁠, three orders of magnitude smaller than this measurement. We therefore expect that our measurement is an apparent orbital period increase caused by the Shklovskii effect. Given the proper motion of the pulsar, we can derive the distance to the pulsar, D = 1.5 ± 0.5 kpc, which, as shown in Table 3, is consistent with other distance measurements.

The measured |$\dot{x}$| is expected to be from the Kopeikin kinematic effects discussed in Section 3. From this, we place an upper limit on the inclination angle of the orbit of i ≤ 84 deg using the total proper motion presented in van Straten (2013), since we do not measure the proper motion in declination. The |$\dot{\omega }$| measurement can have contributions from the same kinematic effects, but may also be consistent with the periastron advance expected from GR. Assuming that the observed |$\dot{\omega }$| comes entirely from GR, we derive the combined mass of the system to be M_p + M_c = 2.5 ± 1.3 M_⊙, which is consistent with a neutron star–white dwarf binary system.

5.5 PSR J1024−0719

The red timing noise for PSR J1024−0719 has a large spectral exponent (i.e. it is very steep). Our usual procedure, as described earlier, requires that we extrapolate the red-noise model obtained from the recent, multi-wavelength data into the earlier data. For PSR J1024−0719, we found that the red-noise model obtained from the recent data did not extrapolate well. Since DM noise is not detectable in the detrended data, we simply modelled the red noise in the entire data set and did not use an additional DM-noise model.

Our measurements are consistent with those of V09, and the measurement precision is improved in all cases. We measure a parallax of π = 0.5 ± 0.3 mas (prior to Lutz–Kelker bias correction), which was undetected by V09. A parallax was previously measured by Hotan et al. (2006) to be π = 1.9 ± 0.8 mas, which is also consistent with our measurements. However, since a parallax was undetected by V09, and PSR J1024−0719 has steep red noise that was unaccounted for by Hotan et al. (2006), the uncertainty for their original parallax measurement is likely to be severely underestimated.

5.6 PSR J1045−4509

A parallax measurement was presented by V09 of π = 3 ± 4 mas, but this measurement may have been affected by uncorrected red noise. We now make the first significant parallax measurement of π = 2.2 ± 1.1 mas.

5.7 PSR J1600−3053

We present the first significant measurement of |$\dot{x}=(-4.2\pm 0.7)\times 10^{-15}$|⁠, which is attributed to the proper motion of the system (Kopeikin 1996), and gives us an upper limit on the inclination angle of i ≤ 67 deg. We detect the first significant parallax of π = 0.48 ± 0.11 mas (compared with π = 0.2 ± 0.3 mas given by V09). From the Shapiro delay companion mass measurement M_c = 0.34 ± 0.15 M_⊙ and the mass function, we can provide a constraint on the pulsar's mass: M_p = 2.4 ± 1.7 M_⊙.

5.8 PSR J1603−7202

For this pulsar, we have the first measurement of a parallax π = 1.1 ± 0.8 mas (prior to Lutz–Kelker bias correction). We also present the measured first derivatives of the orbital period, |$\dot{P_b}=(3.1\pm 1.5)\times 10^{-13}$|⁠, and the projected semi-major axis, |$\dot{x}=(1.36\pm 0.16)\times 10^{-14}$|⁠. Since the GR contribution (for both |$\dot{P_b}$| and |$\dot{x}$|⁠) is negligible for this system, we use the proper motion of the pulsar and our |$\dot{P}_b$| measurement to derive the distance to the pulsar, D = 3.9 ± 1.8 kpc, which is marginally consistent with other distance measurements (Table 3). However, all distance measurements for this pulsar are poor. Using the |$\dot{x}$| measurement, we place an upper limit on the orbital inclination angle of i ≤ 31 deg.

5.9 PSR J1643−1224

We present the first measurement of |$\dot{\omega } = -0.007\pm 0.004$| deg. An improved measurement of |$\dot{x} = (-5.25\pm 0.16) \times 10^{-14}$| allows us to derive an upper limit on the inclination angle of i ≤ 28 deg. For the measured |$\dot{\omega }$| to be the result of GR effects, rather than the assumed kinematic effects, the combined mass of the system would need to be M_p + M_c = 54 M_⊙, an order of magnitude larger than expected. We therefore expect that this measurement is not contaminated by GR contributions and instead results from kinematic effects. However, since these measurements are not well determined and we do not detect a Shapiro delay, we are unable to find a unique solution for the Kopeikin terms i and Ω, and we therefore do not re-parametrize the orbit.

5.10 PSR J1713+0747

Splaver et al. (2005) reported on 12 years of timing observations of this pulsar from the Arecibo observatory. Their analysis was carried out with the JPL DE405 Solar system ephemeris and the time reference was TT(BIPM03). They obtained π = 0.89 ± 0.08 mas, μ_α cos δ = 4.917 ± 0.004 mas yr⁻¹, and μ_δ = −3.933 ± 0.01 mas yr⁻¹. The orbital projection effects caused by this proper motion allowed them to determine Ω = (87 ± 6)°. Their analysis was carried out by ‘whitening’ the residuals using eight time derivatives of the pulse frequency. C11 showed that such whitening can lead to underestimated parameter uncertainties.

V09 obtained a less precise parallax determination of π = 0.94 ± 0.10 mas and proper motion components of μ_α cos δ = 4.924 ± 0.10 mas yr⁻¹ and μ_δ = −3.85 ± 0.02 mas yr⁻¹. They obtained i = (78.6 ± 1.7)° and Ω = (67 ± 17)°. V09 also included |$\dot{P}_b = (41\pm 20)\times 10^{-13}$|⁠.

The most recent VLBI observations of this pulsar (Chatterjee et al. 2009) give proper motion components of |$\mu _\alpha = 4.75^{+0.16}_{-0.07}$| mas yr⁻¹ and |$\mu _\delta = -3.67^{+0.06}_{-0.15}$| mas yr⁻¹ and parallax of |$\pi = 0.95^{+0.06}_{-0.05}$| mas. These values are in fair agreement with our values in Table A4.

Zhu et al. (2015) have analysed 21 years of timing data from the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) for this pulsar to conduct tests of theories of gravity. In their analysis, they measured parameters using a number of different noise models. Using tempo2 with a jitter-based white-noise model and a red-noise model, they measured π = 0.87 ± 0.03 mas, μ_α cos δ = 4.915 ± 0.003 mas yr⁻¹, and μ_δ = −3.914 ± 0.005 mas yr⁻¹. For the binary model, they measured the Kopeikin terms i = (71.9 ± 0.7)° and Ω = (88 ± 2)°, as well as a companion mass of M_c = 0.286 ± 0.012 M_⊙ and |$\dot{P_b}=(0.36\pm 0.17)\times 10^{-12}$|⁠. These parameters are consistent with, and more precise than, our measurements below because of the longer data span and higher timing precision of the NANOGrav data set for this pulsar.

With our analysis we obtain |$\dot{P_b}=(1.7\pm 0.7)\times 10^{-12}$|⁠. The intrinsic |$\dot{P}^{\rm GR}_{b}$| from GW emission is negligible and so we expect that this result comes from the Shklovskii effect. This provides a pulsar distance of D = 3.1 ± 1.2 kpc, which is marginally consistent with other distance measurements (Table 3). We also use the Shapiro delay companion mass measurement M_c = 0.32 ± 0.05 M_⊙ and the mass function, to calculate the pulsar mass, M_p = 1.7 ± 0.4 M_⊙.

5.11 PSR J1730−2304

We present the first measurement of a parallax for this pulsar of π = 1.5 ± 0.3 mas, prior to Lutz–Kelker bias correction. All other parameters are consistent with the previous values from V09.

5.12 PSR J1732−5049

V09 was only able to determine μ_δ = −9.3 ± 0.7 mas yr⁻¹. We now also present a measurement of the proper motion in right ascension, μ_α cos δ = −0.41 ± 0.09 mas yr⁻¹, but parallax was not detected.

5.13 PSR J1744−1134

All parameters are consistent with V09 after accounting for apparent changes resulting from the different Solar system ephemeris and time standard used in the analysis.

5.14 PSR J1824−2452A

PSR J1824−2452A is a solitary pulsar located in the globular cluster M28. The timing residuals for this pulsar exhibit red noise, which may be caused by acceleration within the cluster potential, or timing noise. Our timing model includes dDM/dt but, as we do not include any single-frequency data for this pulsar, we use the original Cholesky routines with a single noise model. We did not measure a significant proper motion in declination nor a parallax even though V09 did publish a proper motion in declination.

5.15 PSR J1857+0943

This pulsar is in an orbit that is highly inclined to our line of sight, allowing for a precise measurement of the Shapiro delay. We do not improve parameter uncertainties for every parameter since V09 made use of publicly available data from the Arecibo observatory to extend the data set, while we chose to use only the PPTA DR1E data set. Using our Shapiro delay companion mass measurement M_c = 0.25 ± 0.03 M_⊙ and the mass function, we calculate the pulsar mass to be M_p = 1.5 ± 0.2 M_⊙.

5.16 PSR J1909−3744

The narrow pulse width, particularly at 10 cm, allows us to achieve very low ToA uncertainties. Recent PPTA data for this pulsar, timed to sub-100 ns precision over more than 10 years, have led to the most stringent limit on the stochastic GW background to date (Shannon et al. 2015). For this reason, this pulsar is an important tool for testing models of galaxy and supermassive black hole formation.

For this pulsar, we use a corrected version of the DR1E data set that is described in Shannon et al. (2013). Previously undetected instrumental offsets were found and corrected, and additional archival 50 cm observations were included to allow measurement of ΔDM over the entire data set. We also include an additional jump corresponding to a software upgrade at MJD 55319.8 that was identified by Shannon et al. (2015). We include dDM/dt in the timing model, which removes the majority of the DM noise. While there is no evidence for red noise in the 10 cm residuals (Shannon et al. 2015), we identify some slight red noise originating from the 20 and 50 cm residuals. This could be the result of instrumental noise or residual interstellar dispersion noise. Because of this noise, we included a red-noise model that sufficiently whitens the residuals.

V09 did not include Kopeikin terms in the timing model, but instead fitted for |$\dot{x}$|⁠; |$\dot{\omega }$| was not measured. We now include the Kopeikin terms, giving the inclination angle, i = 93.52 ± 0.09 deg, and the longitude of ascending node, Ω = 39 ± 10 deg. We measure the orbital period derivative to be |$\dot{P}_b = (5.03 \pm 0.06)\times 10^{-13}$|⁠. The expected contribution from quadrupolar GW emission to this measurement is |$\dot{P}^{\rm GR}_{b}=-2.7\times 10^{-15}$|⁠, two orders of magnitude smaller than this measurement. This expected value was calculated from the measured companion mass M_c = 0.2067 ± 0.0019 M_⊙ and the calculated pulsar mass from the Shapiro delay and mass function of M_p = 1.47 ± 0.03 M_⊙. From this |$\dot{P}_b$| measurement and the proper motion, we derive a distance of D = 1.140 ± 0.012 kpc, which is consistent with the parallax distance. The distances derived from the DM and galactic electron density models are evidently underestimated.

5.17 PSR J1939+2134

PSR J1939+2134 was the first MSP discovered (Backer et al. 1982), and it is currently the second fastest spinning pulsar known. The timing residuals for this pulsar are dominated by red noise.

5.18 PSR J2124−3358

All measured parameters for this pulsar are consistent with V09 and the uncertainties have decreased in all cases.

5.19 PSR J2129−5721

We have the first measurement of the orbital period derivative, |$\dot{P_b}=(7.9\pm 3.6)\times 10^{-13}$|⁠. Using this measurement and the proper motion, we derive a distance of D = 3.2 ± 1.5 kpc, which is consistent with the DM distances (Table 3); however, all distance measurements for this pulsar are poor. We do not yet detect a parallax for this pulsar because of its proximity to the ecliptic pole.

5.20 PSR J2145−0750

We have the first measurement of |$\dot{x}=(8.0\pm 0.8)\times 10^{-15}$|⁠, resulting from the proper motion of the pulsar. Using this value, we place an upper limit on the inclination angle of the orbit of i ≤ 69 deg. All parameters are consistent with the previous results. V09 published |$\dot{x}$|⁠, |$\dot{P_b}$|⁠, and |$\dot{\omega }$| values, but the measurements were not significant. Furthermore, the uncertainties were likely to be underestimated because of the red noise present in these observations.

6 DISCUSSION

6.1 Advantages of using the split-Cholesky algorithm

The methodology that we have used in this paper is based on traditional, frequentist analysis of pulsar timing residuals. An alternate approach is through Bayesian algorithms, such as those described in van Haasteren et al. (2009), van Haasteren & Levin (2013), or the temponest algorithm developed by Lentati et al. (2014a). These algorithms have successfully been used by the NANOGrav and European Pulsar Timing Array groups (e.g. van Haasteren et al. 2011; Arzoumanian et al. 2014). When the same algorithms are applied to PPTA data sets, uncertainties arise in the noise models because of the covariance between DM variations and timing noise processes in the single-wavelength (20 cm) early data. At present, there is no way for current implementations of the Bayesian algorithms to model such non-stationary red noise in the way that the split-Cholesky algorithm allows. Instead, the Bayesian algorithms assume that the noise is wide-sense stationary. Constructing separate red-noise models for the timing noise and DM variations allows us to better understand our noise model, by avoiding large uncertainties in our early data. We therefore chose to use a frequentist approach to analyse our data set since it is less computationally expensive than the Bayesian alternatives and gives us greater control of our noise models.

To demonstrate the necessity for the split-Cholesky algorithm with our data set, we created 500 realizations of PSR J0437−4715 data with red noise and DM noise at the level presented in Section 4. The parameters in the timing model were fitted (including jumps and DM variations) with three different noise treatments: no red-noise model, extrapolated DR1 red-noise model in the Cholesky algorithm only, and a two-component red-noise model with Kolmogorov DM noise in the split-Cholesky algorithm. In Fig. 4, we show the distribution of post-fit parallax values represented by the number of standard deviations from the true value. When no noise model is used (panel a), the parameter uncertainties are clearly underestimated. When the single-component red-noise model is used (panel b), there is a significant improvement but the parameter uncertainties remain underestimated. Finally, when a two-component red-noise model is used with the split-Cholesky algorithm (panel c), we can accurately model the total red noise for the pulsar, and as a result we avoid underestimation of uncertainties. This is true for all parameters with the exception of ν and |$\dot{\nu }$| (as was the case in the original Cholesky algorithm; see C11). The distributions of post-fit values for each parameter in the PSR J0437−4715 timing model (excluding ν and |$\dot{\nu }$|⁠), using the split-Cholesky method, are given in Fig. 5 with the AD statistic used to test the distribution. We see that the distributions are consistent with the expected zero mean, unit variance distribution for all parameters except for declination and proper motion, which have slightly overestimated uncertainties.

Figure 4.

Open in new tabDownload slide

Distribution of parallax values from timing model fits to 500 realizations of simulated PSR J0437−4715 data. σ is the number of standard deviations from the true value for each of the N realizations. We used three different noise treatments: (a) no red-noise model, (b) red-noise model from DR1 data, extrapolated to apply over the entire data set, was used with the Cholesky algorithm, and (c) two-component red-noise model with Kolmogorov DM model is used in the split-Cholesky algorithm. For each panel, the black line is a normal distribution fit to the distribution and the red dashed line is a normal distribution with zero mean and standard deviation equal to the average of the standard tempo2 uncertainties for the 500 realizations, scaled to the same area as the black-line distribution.

Figure 5.

Open in new tabDownload slide

Distributions of parameter values from timing model fits with the split-Cholesky algorithm to 500 realizations of simulated PSR J0437−4715 data. σ is the number of standard deviations from the true value for each realization. Data were simulated with statistically identical red-noise properties to those of PSR J0437−4715.

6.2 PSR J0437−4715: kinematic distance measurement from |${\dot{P}_{b}}$|

The galactic acceleration component can be estimated by combining the accelerations from differential rotation and the galactic potential. The differential rotation acceleration is found using the galactic longitude and distance to the pulsar, and the galactocentric distance and circular velocity of the Sun. Acceleration in the galactic potential varies as a function of height and can be computed from a model of the local surface density of the Galaxy, given the galactic latitude of the pulsar and its parallax distance. V08 used the Holmberg & Flynn (2004) and calculate the total galactic contribution to be |$\dot{P}_{b}^{\rm {Gal}}=\left(-1.8 -0.5\right) \times 10^{-14}=-2.3\times 10^{-14}$|⁠. Since the value of the galactic acceleration is larger than the uncertainty of our observed value, the uncertainty in this component may become important. Bovy et al. (2012) find that the uncertainty for the circular velocity of local sources is approximately 3 per cent by using data from the Apache Point Observatory Galactic Evolution Experiment. From Holmberg & Flynn (2004), the uncertainty in the surface density resulting in the vertical component of acceleration is approximately 10 per cent. Using these uncertainties, we have |$\dot{P}_{b}^{\rm {Gal}}=\left(-2.3\pm 0.08\right)\times 10^{-14}$|⁠, which is small compared to the uncertainty in our measurement.

7 CONCLUSION

We have presented new models for red noise caused by DM variations in single-frequency data sets, obtained by first including deterministic components in the timing model, and then modelling the covariance function of the remaining noise. For 14 of the pulsars in our sample, we detect a linear trend in the DM variations and include this in the timing model, and for four pulsars (PSRs J0613−0200, J1045−4509, J1643−1224, and J1939+2134) we also include new parameters that describe annual DM variations. The new DM-noise models apply only to the early data where excess noise is present, and are used in conjunction with a red-noise model for the frequency-independent noise that is present in the entire data set. These two-component models were used in the new ‘split-Cholesky’ algorithm to whiten the residuals to provide unbiased parameter measurements. We have described this algorithm and demonstrated its effectiveness on simulated data. Model parameters were shown to be unbiased and have accurate uncertainties through simulations based on PSR J0437−4715.

Determining new timing models for the 20 PPTA pulsars in the DR1E data set required these new models and algorithm because of the non-stationary red noise for most pulsars. The models we present provide the best description of the noise currently possible with the PPTA data sets, and result in the most accurate and precise parameter measurements to date for most pulsars in our sample, as well as the detection of several new parameters. Most notably, we presented the first significant parallax measurements for PSRs J1024−0719, J1045−4509, J1600−3053, J1603−7202, and J1730−2304, and determined the distance to PSR J0437−4715 at the 0.16 per cent level of precision. We also measured an improved pulsar mass for PSR J0437−4715, which at M_p = 1.44 ± 0.07 M_⊙ is somewhat lower than the previous measurement.

The analysis described here can easily be applied to future PPTA data releases and to any PTA data with non-stationary noise processes. Longer data sets for PSR J0437−4715 will further improve the distance measurements based on both the parallax and the orbital period derivative from Shklovskii acceleration. If the timing parallax distance becomes more precise than the current VLBI distance, these two independent measurements can be used for example to further improve current constraints on the change to Newton's gravitational constant (Freire et al. 2012). It may also be possible to measure, or place limits on, the acceleration of the pulsar system caused by mass in the Galaxy. With longer data sets, improving |$\dot{\omega }$| measurements, particularly for PSR J0437−4715, can complement the Shapiro delay to further improve measurements of the pulsar mass.

We thank J. Verbiest for assistance with the use of the Lutz–Kelker bias correction webpage, and the referee for their useful feedback. The Parkes radio telescope is part of the Australia Telescope National Facility which is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. YL and GH are recipients of ARC Future Fellowships (respectively, FT110100384 and FT120100595). YL, MB, WvS, and PDL are supported by ARC Discovery Project DP140102578. SO is supported by the Alexander von Humboldt Foundation. VR is a recipient of a John Stocker postgraduate scholarship from the Science and Industry Endowment Fund of Australia. LW and XZ acknowledge support from the Australian Research Council. J-BW is supported by NSFC project no. 11403086 and West Light Foundation of CAS XBBS201322. XPY is supported by NSFC project U1231120, FRFCU project XDJK2015B012, and China Scholarship Council (CSC).

REFERENCES

APPENDIX A: TABLES OF PARAMETERS

The pulsar parameters are listed in the tables below, and a discussion for each pulsar is given in Section 5. The model parameter files, ToA files, and red-noise models for each pulsar are available from the CSIRO Data Access Portal: http://dx.doi.org/10.4225/08/561EFD72D0409.