Pulsars and neutron stars/Overview of the pulsar timing method

Introduction

The assumption that a pulsar is a regular rotator that follows a simple slow-down model forms the basis of a powerful technique that is used for many applications in pulsar research. The technique is known as the Pulsar Timing technique. The pulsar timing method relies on having an initial physical model of the pulsar and a set of pulse arrival times.

Pulsar timing can give information on many pulsar properties depending upon the complexity of the model and the amount (and quality) of the observational data available. For a pulsar which has been observed for a few years it is normally possible to calculate its:

• period: an estimate of the period is usually obtained during the discovery of the pulsar. However, the pulsar timing method provides period determinations with exquisite precision.

• period derivative: this is usually determined after a few days of observations, although for some millisecond pulsars period derivatives are only measurable after many months.

• position: at least a year's worth of data is required for accurate position determinations. Positions are naturally measured in ecliptic coordinates, but usually equatorial coordinates are used.

• dispersion measure: if multiple frequency observations are available then the pulsar timing method can provide a determination of the pulsar's dispersion measure (and, in many cases, its variation with time).

• orbital motion parameters: the set of Keplerian (orbital period, epoch of periastron, eccentricity, longitude of periastron and the projected semimajor axis of the orbit) parameters can usually be determined for pulsars in binary systems.

With sufficient data it is also possible to determine post-Keplerian (relativistic) binary parameters, the pulsar's parallax and proper motion.

Various software packages (described below) have implemented the pulsar timing method. These packages require a set of "pulse arrival times" and an initial "timing model" for each pulsar to be processed. The pulse arrival times are often called "Times of Arrival" (ToAs) or "site times of arrival" (SATs) and are usually presented in an "arrival time file", which traditionally have the file extension ".tim". The timing model is sometimes called the "timing ephemeris", the "pulsar model" or the "pulsar parameters" and is usually given in a file with extension ".par".

Measuring pulsar arrival times

A folded pulse observation has a known start time. The pulse arrival time is usually determined by cross-correlating a known template of the pulse shape with the actual observation to determine the actual arrival time. This cross-correlation is usually carried out in the Fourier domain. The Taylor (1990) paper provides the basic overview, but is difficult now to obtain. Here we re-describe his work.

We assume that the profile has been binned into N bins and is defined by $p_{j}$. The standard template is the same $s_{j}$. The discrete Fourier transform of the profile and the template can be written as:

$P_{k}e^{i\theta _{k}}=\sum _{j=0}^{N-1}p_{j}e^{i2\pi jk/N}$

$S_{k}e^{i\phi _{k}}=\sum _{j=0}^{N-1}s_{j}e^{i2\pi jk/N}$

The template and the profile are related by an unknown phase offset, a scaling factor and noise:

$P_{k}e^{i\theta _{k}}=aN+bS_{k}e^{i(\phi _{k}+k\tau )}+G_{k}$

The best match is obtained when the difference between the template and the profile is minimised and so we want to minimise:

$\chi ^{2}(b,\tau )=\sum _{k=1}^{N/2}\left[{\frac {P_{k}-bS_{k}e^{i(\phi _{k}-\theta _{k}+k\tau )}}{\sigma _{k}}}\right]^{2}$

These can be calculated analytically:

${\frac {\delta \chi ^{2}}{\delta b}}={\frac {2b}{\sigma ^{2}}}\sum S_{k}^{2}-{\frac {2}{\sigma ^{2}}}\sum P_{k}S_{k}\cos(\phi _{k}-\theta _{k}+k\tau )=0$

${\frac {\delta \chi ^{2}}{\delta b}}={\frac {2b}{\sigma ^{2}}}\sum kP_{k}S_{k}\sin(\phi _{k}-\theta _{k}+k\tau )=0$

These equations allow the routine to determine $\tau$ and then $b$. Of course, it is important to determine uncertainties on the pulse arrival times. These can be estimated using

$\sigma _{b}^{2}=\left({\frac {\delta ^{2}\chi ^{2}}{\delta b^{2}}}\right)^{-1}={\frac {\sigma ^{2}}{2\sum S_{k}^{2}}}$

$\sigma _{\tau }^{2}=\left({\frac {\delta ^{2}\chi ^{2}}{\delta \tau ^{2}}}\right)^{-1}={\frac {\sigma ^{2}}{2b\sum k^{2}P_{k}S_{k}\cos(\phi _{k}-\theta _{k}+k\tau )}}$

These algorithms have been updated to determine arrival times in the time domain (instead of the frequency domain as shown above) and new algorithms are being developed to measure pulse arrival times over wide-bandwidths.

For the tempo or tempo2 software packages, pulse arrival times are presented in files similar to the following (note that there are slight differences in the formatting between the software packages):

FORMAT 1
t110120_155047.czFTp 3094.50000000 55581.68263874327820062 2.33600 PKS
s110120_155046.czFTp 732.27500000 55581.68276815691150006 2.64800 PKS

Each observation is represented on a single line in the file. The first column gives an identifier for the observation (usually a filename). The second column provides the observation frequency (in MHz), the third column gives the actual pulse arrival time (in MJD) and the fourth column provides the uncertainty on the arrival time measurement (in microseconds). The final column indicates which telescope was used for the observation.

Obtaining an initial timing model

When a pulsar is discovered the following parameters are usually known:

1. Position (RAJ, DECJ)
2. Period (P0) or frequency (F0)
3. The dispersion measure (DM)
4. Whether the pulsar is in a fast binary

A parameter file consists of the following required parameters:

PSRJ <pulsar name>
RAJ <right ascension hh:mm:ss.ss>
DECJ <declination dd:mm:ss.ss>
P0 <pulse period in seconds> or F0 <pulse frequency in Hz>
PEPOCH <MJD epoch of period determination>
DM <dispersion measure (cm^-3 pc)>

As more parameters are measured they can be included in the parameter file. For a known pulsar, an initial timing model can be obtained from the ATNF pulsar catalogue by typing the pulsar name into the "Pulsar name" box and then selecting "Get ephemeris". For instance, the current ephemeris for PSR J1939+2134 is:

PSRJ            J1939+2134
RAJ             19:39:38.561297          2.000e-06
DECJ            +21:34:59.12950          4.000e-05
DM              71.0227                  9.000e-04
PEPOCH          52601
F0              641.928233642            1.200e-08
F1              -42.91E-15               6.000e-17
PMRA            0.072                    2.000e-03
PMDEC           -0.415                   3.000e-03
RM              +6.7                     6.000e-01
PX              0.13                     1.300e-01
EPHVER          2
UNITS           TDB

The values in the third column represent the uncertainties on the parameters. However, such uncertainties are currently not accounted for in the existing pulsar timing software packages. Note that the pulsar catalogue is a catalogue of the "best" available pulsar parameters. In a few cases such parameters do not form the best timing model, but usually the catalogue can be used to get a sufficient model for initial processing.

Calculating barycentric arrival times

Pulse arrival times are measured at an observatory that is (usually) on the surface of the Earth. Of course, the Earth is rotating and moving in its orbit. It is therefore necessary to convert the measured pulse arrival times to the times that would have been measured if the observatory were situated at the centre of mass of the solar system (known as the solar system barycentre or SSB). After the site arrival times have been modified to represent an observatory at the barycentre they are referred to as barycentric arrival times (BATs).

Throughout this section we use the example of an observation of PSR J1744+1134 from the Parkes observatory. The arrival time file contains:

s110205_233319.zCFTp 1367.89400000 55598.00404304428645119 0.07600 PKS

The measured pulse arrival time is therefore 55598.00404304428645119. We will use this to demonstrate the processes followed during the pulsar timing method.

Converting to Terrestrial Time

The measured site-arrival-time was determined using an observatory time-standard. No time standard is perfect and so it is necessary to convert the measured time to a realisation of terrestrial time, TT. It is not possible to convert directly from an observatory clock to Terrestrial Time. Instead a set of clock corrections are applied. In some cases multiple different methods for reaching Terrestrial Time are possible. For our example, the Parkes clock is first converted to the global positioning satellites (GPS) and from there to UTC and finally to Terrestrial Time as realised by International Atomic Time, TT(TAI). Within the tempo/tempo2 software packages these clock corrections are updated manually and provided as ascii files that contain the difference between two time standards at a given time. An example is given in the Figure in which we show the corrections (in seconds) between the Parkes observatory and GPS.

The clock corrections between the Parkes observatory clock and GPS

The pks2gps.clk file tells us that the following clock corrections straddle our requested ToA:

55597.16493 -0.000000967864
55598.16840 -0.000000963384

The first column is the MJD of the clock correction measurement and the second column is the value of the clock correction (in seconds). We usually use a simple linear interpolation to determine the clock correction for the actual measured ToA. In this case we have:

$\Delta T_{\rm {pks2gps}}=-9.6412\times 10^{-7}s$

and using a similar procedure:

$\Delta T_{\rm {gps2utc}}=-4.49\times 10^{-9}s$
$\Delta T_{\rm {utc2tai}}=34s$
$\Delta T_{\rm {tai2tt}}=32.184s$

(Note that since MJD 41317 $\Delta T_{\rm {utc2tai}}$ simply represents the number of leap seconds.)

For our example, the total time correction needed to convert from the measured arrival time to TT(TAI) is therefore the summation of the above values:

$\Delta T_{\rm {clk}}=\Delta T_{\rm {pks2gps}}+\Delta T_{\rm {gps2utc}}+\Delta T_{\rm {utc2tai}}+\Delta T_{\rm {tai2tt}}=66.18399903139s$

Converting to Barycentric Time

It is now necessary to determine the extra delays required to convert the clock-corrected arrival time to the barycentric arrival time.

The Roemer Delay

The most obvious delay is simply the extra travel time for the pulse to reach the SSB after it was detected at the observatory. This requires knowledge of the observatory position with respect to the centre of the Earth and the Earth's position with respect to the SSB. It is also necessary to know the position of the pulsar.

Determining the Earth's position with respect to the SSB

The Earth's position with respect to the SSB is determined using a planetary ephemeris. Current ephemerides are described in Section XX. For this example we will make use of the Jet Propulsion Laboratory (JPL) DE421 ephemeris.

Shapiro delays

The Shapiro delay, $\Delta S_{i}$, represents the Shapiro delay caused by the i'th mass in the Solar System. Here we must include the Sun (as it has the largest effect and the major planets). The Shapiro delay corrects for the curvature of space-time caused by the presence of masses (Shapiro 1964):

$\Delta S_{i}=-2{\frac {GM_{i}}{c^{3}}}\ln \left[R(1+\cos \theta \right]$

where $G$ is the gravitational constant, $M_{i}$ is the mass of the i'th object, $R$ is the magnitude of the vector between the observatory and the mass and $\theta$ is the pulsar-observatory-mass angle. CHECK THIS

The total Shapiro delay from all masses is simply the addition of the individual delays:

$\Delta T_{\rm {shapiro}}=\sum _{i}\Delta S_{i}$

Give table of Shapiro delay sizes

Calculating pre-fit pulsar timing residuals

The derived pulsar emission time and the given pulsar timing model is used to form the timing residuals, $R_{i}$:

$R_{i}={\frac {\phi _{i}-N_{i}}{\nu }}$

where $\phi _{i}$ describes the time evolution of the pulse phase based on the model pulse frequency ($\nu$) and its derivatives (and glitch parameters). $N_{i}$ is the nearest integer to $\phi _{i}$.

Improving the timing model

Terms corresponding to offsets in model parameters are fitted to the residuals in order to improve the measurements of the parameters.

Studying post-fit timing residuals

Pulsar timing residuals that statistically deviate from zero are caused by 1) incorrect parameters in the timing model, 2) an inadequate timing model, 3) errors in the timing software package or 4) unmodelled, or incorrectly modelled, physical phenomena that affect the pulse arrival times. In contrast if the timing residuals are consistent with zero then the software is accurately accounting for all the physical affects within the uncertainties of the timing residuals. It is now possible to produce timing residuals with an rms timing residual of ~100ns over a decade or more.

Pulsar timing residuals are analysed using various statistical tests:

How good is the data set?

The most common measure of the "quality" of a pulsar timing residual data set is the root-mean-square (rms) residuals. As ToA uncertainties significantly vary, it is common to present a weighted rms value.

should mention chisq

Pulsar timing software

The following software packages exist for pulsar timing:

• PSRTIME (Jodrell Bank)
• TIMAPR (Bonn)
• ANTIOPE (Nancay)
• CPHAS (Hartebeesthoek)
• TEMPO
• TEMPO2
• PINT