Pulsars and neutron stars/Overview of the pulsar timing method
- 1 Introduction
- 2 Measuring pulsar arrival times
- 3 Obtaining an initial timing model
- 4 Calculating barycentric arrival times
- 4.1 Converting to Terrestrial Time
- 4.2 Converting to Barycentric Time
- 4.3 The Roemer Delay
- 4.4 Conversion to barycentric time and the Einstein delay
- 4.5 Shapiro delays
- 4.6 Dispersion delays
- 4.7 Atmospheric delays
- 5 Secular motion
- 6 Converting to Pulsar Emission Times
- 7 Calculating pre-fit pulsar timing residuals
- 8 Understanding the timing residuals
- 9 Improving the timing model
- 10 Studying post-fit timing residuals
- 11 Pulsar timing software
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 . The standard template is the same . The discrete Fourier transform of the profile and the template can be written as:
The template and the profile are related by an unknown phase offset, a scaling factor and noise:
The best match is obtained when the difference between the template and the profile is minimised and so we want to minimise:
These can be calculated analytically:
These equations allow the routine to determine and then . Of course, it is important to determine uncertainties on the pulse arrival times. These can be estimated using
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:
- Position (RAJ, DECJ)
- Period (P0) or frequency (F0)
- The dispersion measure (DM)
- 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 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:
and using a similar procedure:
(Note that since MJD 41317 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:
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 telescope position
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.
Conversion to barycentric time and the Einstein delay
The Shapiro delay, , 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):
where is the gravitational constant, is the mass of the i'th object, is the magnitude of the vector between the observatory and the mass and is the pulsar-observatory-mass angle. CHECK THIS
The total Shapiro delay from all masses is simply the addition of the individual delays:
Give table of Shapiro delay sizes
Converting to Pulsar Emission Times
Calculating pre-fit pulsar timing residuals
The derived pulsar emission time and the given pulsar timing model is used to form the timing residuals, :
where describes the time evolution of the pulse phase based on the model pulse frequency () and its derivatives (and glitch parameters). is the nearest integer to .
Understanding the timing residuals
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)