TEMPO 2, a new pulsar timing package – II. The timing model and precision estimates

TEMPO 2 is a new software package for the analysis of pulsar pulse times of arrival. In this paper, we describe in detail the timing model used by TEMPO 2, and discuss limitations on the attainable precision. In addition to the intrinsic slow-down behaviour of the pulsar, TEMPO 2 accounts for the effects of a binary orbital motion, the secular motion of the pulsar or binary system, interstellar, Solar system and ionospheric dispersion, observatory motion (including Earth rotation, precession, nutation, polar motion and orbital motion), tropospheric propagation delay, and gravitational time-dilation due to binary companions and Solar system bodies. We believe that the timing model is accurate in its description of predictable systematic timing effects to better than 1 ns, except in the case of relativistic binary systems where further theoretical development is needed. The largest remaining sources of potential error are measurement error, interstellar scattering, Solar system ephemeris errors, atomic clock instability and gravitational waves.


I N T RO D U C T I O N
The core of any pulsar timing software package is the mathematical model for the times of arrival (TOAs) of pulses at the observatory. As the pulsar population grows more diverse and observing systems increase in sensitivity, the accuracy demanded of timing models grows accordingly. For many years, the standard choice for pulsar timing was the TEMPO 1 software package. Accurate at the 100-ns level, the timing model of TEMPO was adequate for numerous important discoveries (e.g. Taylor, Fowler & McCulloch 1979;Taylor & Weisberg 1989;Wolszczan & Frail 1992;Kaspi, Taylor & Ryba 1994;van Straten et al. 2001). However, motivated primarily by the possibility of directly detecting gravitational radiation through its effect on pulse TOAs, current observing campaigns are setting highly ambitious goals in terms of data precision and volume. These projects aim to obtain of the order of ∼10 4 independent TOA measurements with a typical root-mean-square (rms) uncertainty of 100 ns (e.g. Jenet et al. 2005), dictating a maximum allowable systematic error of 1 ns.
Motivated by the need for improved accuracy as well as greater flexibility in fitting and analysis procedures, we have developed a new pulsar timing package called TEMPO2. The architecture, capabilities and usage of TEMPO2 are described by Hobbs, Edwards & Manchester (2006) (hereafter Paper I). The purpose of this work is to describe the timing model in detail, and provide estimates of the accuracy of each component of the model. Section 2 begins with an overview of the model and definitions of the various relativistic coordinate frames and proper times (Section 2.2). This is followed by an expansion of the complete geometric propagation delay accurate to better than 1 ns (Section 2.3). The remainder of Section 2 discusses all oWther delay terms in the timing formula, and provides a detailed description of how each contribution is evaluated according to the parameters of the model and other information. Section 3 gives a detailed account of the accuracy with which each term is computed, and discusses possible avenues for further improvement of the timing model. Finally, for the convenience of comparison with the commonly used TEMPO software package, Section 4 details the differences between TEMPO and TEMPO2. Tabular summaries of timing model parameters and variables are provided in Appendix A.

Coordinate frames and time variables
The relativistic frame transformation between observatory proper time and pulsar proper time is computed in several steps. This section outlines the steps, gives the names of the various reference frames and introduces the notations used to refer to time variables in each frame. The time variable for each frame may appear subscripted with the letter 'a' or 'e', where 'a' denotes a TOA and 'e' a time of emission, at the origin of the frame denoted by the superscript. A time variable may also appear without a subscript, meaning that it is an arbitrary coordinate time variable. As such, time variables without subscripts are related by coordinate transformations, whereas subscripted time variables additionally include propagation delays. Spatial three-dimensional vectors are defined in terms of a right-handed Cartesian system.
The pulsar timing procedure begins with a measured TOA, which after correcting for imperfections in the observatory clock yields a TOA t obs a . The clock-correction process in TEMPO2 operates through an interpolation of tables of measured offsets between pairs of clocks and their variation with time (see Paper I), and includes offsets relating to the fact that observatory clocks usually approximate Coordinated Universal Time (UTC 2 ) rather than Terrestrial Time (TT). The corrected TOA is measured in the TT time-scale, which is defined to differ from coordinate time at the geocentre by a constant rate. It should be noted then that times referred to TT are only approximately proper times of the observatory. This transformation to Geocentric Coordinate time-scale (TCG) is done in the same step as transformation to Barycentric Coordinate Time (TCB), so no variable is assigned to this time scale. The spatial part of the terrestrial frame is the International Terrestrial Reference System (ITRS), while the geocentric counterpart is the Geocentric Celestial Reference System (GCRS).
Coordinate times t SSB in the frame of the Solar system barycentre (SSB) are obtained by application of the gravitational redshift and special relativistic time-dilation integral ('Einstein delay') described in Section 2.5.1. Coordinate time (TCB) in this frame represents the proper time that would be measured by an observer located at the SSB were the gravitational field of the Sun and planets not present. The spatial part of the barycentric frame is the Barycentric Celestial Reference System (BCRS), while directions in this frame define the International Celestial Reference System (ICRS). The origin of the BCRS is the SSB, hence t SSB a refers to the TOA of the pulse at the barycentre. Tracing the propagation backwards towards the source, the next reference frame is that which is comoving with the pulsar, or if it is a member of a binary system, the binary barycentre (BB). Coordinates in this frame are defined such that the coordinate time t BB would be equal to the proper time of an observer at the origin of the frame, were the gravitational field of the pulsar and its companion not present. The transformation of the TOA to this frame is discussed in Section 2.6.3. The orientation of this frame is chosen such that the Lorentz transformation relating it to the SSB frame is rotation-free, so that to zeroth order the corresponding basis vectors are parallel. The spatial origin of this frame is the BB. The time origin is chosen such that t BB has zero offset from t SSB at t BB = t SSB = t pos , where the latter is an epoch of position set by the user.
The time-scale for measurements is the proper time t psr measured at the pulsar, or more precisely at its centre, were its gravitational field not present. For isolated pulsars t psr = t BB . For binary pulsars, t psr − t BB is the binary Einstein delay (Section 2.7.4). As with the SSB-BB case, the BB-pulsar transformation is rotation-free. The spatial origin of this frame is the centre of gravitational field of the pulsar, so that t psr e is the equivalent time of emission from the centre of the star. This differs from the actual time of emission by (to first order) the light traveltime to the true point of emission, which is assumed to be stable on the time-scale of interest.

Geometric propagation delay
Although it is convenient to break the propagation delay into separate contributions and treat them in the order in which they occur along the photon path, some terms in the model depend jointly on the motion of the observatory and the motion of the pulsar (relative to the SSB). For this reason, in this section, we derive the complete geometric propagation delay, that is the time of propagation neglecting dispersion, refraction and relativistic delays. The resultant terms are then assigned to different parts of the propagation path and treated in the relevant sections below.
The geometric propagation delay is simply the Euclidean distance from the observatory to the pulsar, divided by the vacuum speed of light, c. The displacement vector from the observatory to the pulsar is the sum of the barycentric position of the observatory (r), the barycentric position of the pulsar (or BB) (R 0 ) at a given epoch (t pos ), the position of the pulsar with respect to the BB (b, zero for isolated pulsars), and the displacement (k) of the BB (or the pulsar itself, if isolated) in the time elapsed since the epoch t pos , owing to the initial velocity of the system (relative to the SSB) and acceleration in the Galactic or cluster gravitational field: (1) Squaring equation (1) and expanding the square root to third order, where Neglecting terms of the order of |R 0 | −3 and denoting the initially radial and transverse components of each vector using subscripts, that is, a = a · R 0 , a ⊥ = a − a R 0 /|R 0 | and hence, a ⊥ · b ⊥ = (a × R 0 ) · (b × R 0 ), the following relation is found: The first four terms in equation (5) are the initial SSB-BB distance (Section 2.6.2), the secular displacement in the initially radial direction (Section 2.6.2), the projection of the observatory-SSB vector on the initial line of sight (Section 2.5.3), and the projection of the binary motion on the initial line of sight (Section 2.7.1). The terms in the first pair of parentheses correspond to the Shklovskii effect (Section 2.6.2), annual proper motion (Section 2.5.3), secular changes in the apparent orbital viewing geometry (Section 2.7.1), annual parallax (Section 2.5.3), annual-orbital parallax (Section 2.7.1) and orbital parallax (Section 2.7.1). The second, third and fourth terms in the second pair of parentheses represent variations in the effects just discussed due to radial motion of the pulsar or binary system, the motion of the Earth around the SSB, and the motion of the pulsar about the BB. Although these terms give rise to an additional 18 contributions to the timing formula, to reach the stated accuracy goal of 1 ns, for a 20-yr observing campaign on any of the presently known millisecond pulsars only three of them need to be considered: the secular change (Section 2.6.2) and annual modulation (Section 2.5.3) of the Shklovskii effect and the secular change of the annual proper motion (Section 2.5.3). The displacement due to secular motion, k may be broken into first and second derivatives: where μ is the velocity divided by the distance, or in essence a three-dimensional proper motion, and a is an acceleration vector accounting for Galactic differential rotation and gravitational acceleration. The acceleration is of the order of ∼ 10 −11 m s −2 (Bell & Bailes 1996), contributing only ∼500 km to k at either end of a 20-yr observation campaign centred on t pos . Equation (5) can be expanded in terms of this parametrization of k, giving rise to an additional 27 terms. We have estimated the magnitude of these terms for all pulsars in the Australia Telescope National Facility (ATNF) catalogue 3 , including an allowance for greater acceleration for pulsars in globular clusters (a ∼ 10 −8 m s −2 ; Freire et al. 2001). Only two of the extra terms involving a exceed 1 ns for a 20-yr observing campaign, the first contributing to k and the second contributing to the Shklovskii term of equation (5). For all other terms of equation (5) involving k, the acceleration is neglected.

Top-level timing formula
Having described the coordinate frames used and derived the basic Euclidean propagation length, it remains to specify the details of the calculation of the various delay terms associated with frame transformations, the vacuum propagation delay, and corrections due to non-unity refractive indices and space-time curvature along the propagation path. For convenience, the timing formula is divided into three main contributions, which will be further broken down in the sections that follow: (Section 2.5) includes the coordinate transformation to the SSB frame, vacuum propagation delays associated with the Earth's orbital motion and the spin, precession and nutation of the Earth (on which the observatory is located) and excess delays owing to the passage of the signal through the Earth's atmosphere and the Solar system. This term relates the measured TOA to the TCB TOA at the SSB: t SSB a = t obs a − . IS (Section 2.6) includes the transformation to the BB frame, vacuum propagation delays due to the secular motion of the system, and excess propagation delays due to passage of the signal through the interstellar medium (ISM). This term relates the TCB TOA at the SSB to the BB coordinate TOA at the BB: t BB a = t SSB a − IS .
B includes the transformation to the pulsar frame, vacuum delays due to the binary orbital motion, and excess delays due to passage of the signal through the gravitational field of the companion. This term relates the BB coordinate TOA at the BB to the pulsar proper time of emission: t psr e = t BB a − B . Note that since some of the geometric delay terms involve two of the three components (Solar system, secular or binary) of the relative motion of the observatory and the pulsar, the choice of which of these terms are assigned to each of the three top-level terms is somewhat arbitrary. See Section 2.3 for a detailed listing of where each term is discussed.
Once the time of emission is obtained, it remains to compute the phase of the pulse train at that time. This relationship, φ(t psr e ), is discussed in Section 2.8. For the correct set of model parameters and in the absence of measurement noise, the phase for every computed time of emission is an integer number of cycles (turns). The fractional part of φ(t) is the timing residual, which is more often expressed in time units through division by the pulse frequency.

Forming barycentric TOAs
The term denoted by in equation (7) includes all effects pertinent to the transformation of the observed TOA (t obs a ) to an equivalent TOA for the same pulse wavefront at the SSB. This, in turn, can be broken into several steps. These are: atmospheric delays, vacuum retardation due to observatory motion ('Roemer delay' and parallax), excess propagation delay due to dispersion, the effects of relativistic frame transformation ('Einstein delay'), and the excess delay experienced by rays as they pass through the gravitational potential of Solar system bodies ('Shapiro delay'). These terms are denoted as follows:

Einstein delay
The space-time coordinates of pulse arrival events are specified in the coordinate frame of the observatory, in which pulse timing anomalies may manifest due to spin and orbital acceleration and variations in gravitational potential. To avoid these effects, pulse arrival events are transformed to the quasi-inertial frame of the SSB. The two frames are related by a relativistic four-dimensional space-time transformation. The spatial part of the event is the displacement of the observatory from the geocentre at the instant of reception, which is altered by a negligible amount due to special and general relativistic length contraction. The effects of relativistic time-dilation, on the other hand, cannot be neglected. TEMPO2 uses the numerical time-dilation results of Irwin & Fukushima (1999), who used the DE405 Solar system ephemeris ) to compute the time-dilation integral: where U ⊕ is the gravitational potential at the geocentre due to all Solar system bodies except the Earth, and v ⊕ is the velocity of the geocentre relative to the SSB. The first two terms describe, to order 1/c 2 , the gravitational redshift and special relativistic time-dilation, respectively. These contribute in roughly 2:1 ratio ( 1 2 U ⊕ v 2 ⊕ 10 −8 c 2 ), for a mean drift of ∼ 1.5 × 10 −8 , or ∼0.5 s yr −1 . The main source of timevariation in the integrand is due to the orbital acceleration of the Earth, which via conservation of energy causes approximately equal changes (of the same sign) in the first two terms, amounting to ∼3 × 10 −10 in rate amplitude, or ∼2 ms in integrated delay amplitude. The remaining terms apply a correction for higher-order relativistic terms [ L (PN) C = 1.097 × 10 −16 ; Fukushima 1995] and asteroids [ L (A) C = 5 × 10 −18 ; Fukushima 1995], made in mean rate only.
The time-dilation integral of equation (9) relates the coordinate times of clocks at the geocentre and the SSB. For observers located on the surface of the Earth, corrections must be made for the differential time-dilation and gravitational redshift between the observatory and the geocentre, yielding where s is a vector from the geocentre to the observatory,ṙ ⊕ is the velocity of the geocentre with respect to the barycentre, and W 0 = 6.969 290 134 × 10 −10 c 2 approximates the gravitational plus spin potential of the Earth at the geoid. (See Section 2.5.3 for the calculation of s andṙ ⊕ .) Note that, while the value of the potential varies geographically and temporally, the most recent definition of TT (Rickman 2001), the terrestrial time-scale to which TEMPO2 refers pulse TOA measurements, takes the exact value of W 0 /c 2 quoted above as definitive of the rate difference between TT and SI coordinate time at the geocentre. Given a common epoch of t 0 = 43 144.000 3725 (Modified Julian Date), equations (9) and (10) relate measured TT to the IAU recommended coordinate time-scales, TCG (geocentric) and TCB (barycentric).
It should be noted that the time ephemeris of Irwin & Fukushima (1999) in fact takes its time argument, T eph , in a reference frame that is a linear transformation of TCB, making the computation of the Einstein delay strictly an inverse problem. In practice, since the argument time-scale has zero mean rate difference to TT, the latter can be substituted with negligible error in the result (Irwin & Fukushima 1999). In addition, the evaluation of equation (10) requires consulting the Solar system ephemeris (Section 2.5.3) which also takes T eph times as its argument. Since T eph is computed from TCB (equation 10), a circular dependency problem arises. This is solved to sufficient accuracy by two passes through an iterative refinement process.

Atmospheric delays
The group velocity of radio waves in the atmosphere differs from the vacuum speed of light. Refractivity is induced both by the ionized fraction of the atmosphere (mainly in the ionosphere) and the neutral fraction (mainly in the troposphere).
The tropospheric propagation delay be separated into the so-called 'hydrostatic' and 'wet' components. In Very Long Baseline Interferometry (VLBI) and Global Positioning System (GPS) applications, these components are typically modelled as the product of the delay induced at the zenith, and a so-called 'mapping function', which specifies the ratio between the zenith delay and the delay in a given direction. For a planar atmosphere, the mapping function is simply given by csc where is the source elevation angle. More advanced mapping functions take into account the curvature of the atmosphere, assuming azimuthal symmetry. The effect of azimuthal asymmetry is typically less than a nanosecond even at very low elevation angles (e.g. MacMillan 1995), and need not be considered for the present purpose.
The hydrostatic component contributes approximately 90 per cent of the total delay, and may be computed a priori, assuming a given mixture of gases in hydrostatic equilibrium. TEMPO2 uses the formula of Davis et al. (1985), re-written as a time-delay hz = P/43.921 kPa c(1 − 0.002 66 cos 2ϕ − 3.6 × 10 6 m/H ) , where hz is the hydrostatic zenith delay, P is the surface atmospheric pressure, ϕ is the geodetic latitude of the site, and H is its height above the geoid. In combination with the mapping function, this leads to a timing term of amplitude ∼7.7 ns × csc , mostly on a diurnal time-scale. If atmospheric pressure data are unavailable, TEMPO2 uses a canonical value of one standard atmosphere (101.325 kPa). TEMPO2 uses the Niell mapping function (NMF; Niell 1996), which is of comparable accuracy to other published mapping functions but does not require meteorological data: it depends only on the source elevation. The wet component of the tropospheric propagation delay is highly variable and cannot be predicted accurately. Fortunately, it is small. The zenith wet delay (ZWD; wz ) may be measured using appropriate analyses of radiosonde, water vapour radiometer, GPS or VLBI observations (e.g. Niell et al. 2001), of which only GPS is likely to be obtainable on a routine basis at most radio observatories. If such measurements are available, TEMPO2 can make use of them in conjunction with the NMF to account for the effects of tropospheric water vapour on pulse TOAs. The ZWD could conceivably be included as a free parameter in the pulsar timing model, but because the effect is small and varies from day to day, this would not be possible with the sensitivity of presently attainable observing systems. For this reason, if no tabulated ZWD information is available, the effect is neglected.
Neglecting dispersion (below), the total atmospheric delay is written as where m h ( ) and m w ( ) are the NMFs for the hydrostatic and wet components. Values for P and wz are to be provided as input data to TEMPO2, or default values of P = 101.325 kPa and wz = 0 are used. The passage of the signal through the ionosphere induces a dispersive delay that varies strongly with the solar activity cycle and also on seasonal, diurnal and shorter time-scales. The integrated column density of electrons ('total electron content', TEC) in the ionosphere typically lies in the range 5-100 TECU (1 TECU =10 16 m −2 3.2 × 10 −7 cm −3 pc), corresponding to variable propagation delays of the order of 7-130 ns ×(f/1 GHz) −2 , where f is the observing frequency. Models of the ionosphere are available (e.g. Schaer 1999;Bilitza 2001); however, the predicted TEC is accurate to only ∼3-20 TECU (Mamoru, Kondo & Kawai 2003). In any case, ionospheric dispersion is inseparable from and smaller than uncertainties in interplanetary (Section 2.5.4) and interstellar (Section 2.6.1) dispersion, necessitating the inclusion of a fittable, fully time-variable dispersion parameter (Section 2.6.1).

Roemer delay and parallax
The Roemer delay is the simple vacuum delay between the arrival of the pulse at the observatory and the SSB, not including effects related to the binary motion of, or finite distance to, the pulsar: whereR BB ≡ R BB /|R BB | is a unit vector in the direction of the BB at the time of observation. TEMPO2 performs this calculation in the BCRS. The vectorR BB is constructed from the spherical coordinate angles of the ICRS source direction, right ascension (α) and declination (δ), at time t pos , and the Cartesian components of the proper motion of the source in the plane of the sky (μ α and μ δ ) and along the line of sight (μ ), all of which are fittable parameters in TEMPO2: wherê and Alternatively, ecliptic coordinates (λ, β, μ λ , μ β ) may be used, in which case: and and Here 0 = 84 381.405 78 arcsec (Harada & Fukushima 2004) is the mean obliquity of the ecliptic at J2000.0. The four terms that result from the expansion of equation (14) derive from corresponding terms in the expansion of equation (5), such that Here we have assumed that k =μ |R 0 |(t BB a − t pos ), safely neglecting any acceleration or higher-order secular motion of the pulsar or BB (Section 2.3). In principle, t BB a cannot be computed without first knowing R ; this is resolved by initially setting t BB a to the TCB TOA at the observatory, and iteratively refining and IS until convergence. Equation (24) consists of four terms corresponding to the expression in equation (14) of the instantaneous source direction as the sum of four contributions. These are the initial direction (yielding the first term of equation 24), the proper motion (second term), a correction to maintain unit length after addition of the proper motion (third term) and the acceleration or deceleration of the proper motion owing to the fact that over time, the proper motion moves the source vector, causing part of the initially radial velocity to acquire a transverse component (fourth term). The third term can also be interpreted as the proper-motion-induced second-order reduction in the projection of the source direction on the radial component of r that accompanies the first-order increase in its projection on the transverse component of r. Another alternative interpretation of the third term is that it corrects the Shklovskii term (Section 2.6.2) for the distance variations induced by the Earth's orbital motion. The fourth term, first introduced by N. Wex (1997, unpublished contribution to TEMPO), provides one of two possible ways of accessing the elusive radial velocity, the integrated first-order contribution of which (k ) is otherwise inseparable from a Doppler shift of the pulse and orbital frequencies and their derivatives. The second method is described in Section 2.6.2.
The first term of equation (13) describes the difference in TOA of a pulse at the SSB compared to that at the point obtained by projecting the SSB-observatory vector on the SSB-BB vector. The TOA at that point differs again from the TOA at the observatory due to the curvature of spherical 'wavefronts' connecting photons emitted simultaneously from the pulsar. This term is referred to as the 'parallax' term, although it differs somewhat from stellar parallax: positional astrometry measures the three-dimensional orientation of the wavefront normal, whereas pulsar timing measures the position of the wavefront in its intersection with the ecliptic. To first order, the curvature is proportional to the square of the lateral displacement of the observatory from the SSB-BB vector, yielding an excess delay corresponding to one of the terms of equation (5): Here the 'parallax distance', d p is substituted for |R 0 | to allow this effect to be separated from others involving the source distance. The conventional 'parallax' ≡ 1 au/d p , being the angle subtended from the source by a line of 1 au in length at the SSB, is a fittable parameter in TEMPO2. Its value is specified in units of milliarcseconds, corresponding numerically to the reciprocal of the distance in kiloparsecs.
In computing the Roemer and parallax delays, the vector r is constructed in two steps, from the SSB to the geocentre and from the geocentre to the observatory, that is, r = r ⊕ + s. The first part (r ⊕ ) is accomplished with the aid of a numerical Solar system ephemeris. The current default ephemeris is DE405 , which is aligned within a known uncertainty to the ICRS, but uses a time-scale and length-scale slightly different from the SI second and metre . TEMPO2 uses the equations and constants of Irwin & Fukushima (1999) to appropriately transform the TCB site TOA (Section 2.5.1) for input to the ephemeris (using an offset and scale), and appropriately scales the output vectors to SI units. It should be noted that the use of alternative ephemerides that require a different transformation procedure will introduce errors (see Section 3.5).
The second part of the observatory position vector, s, is obtained by transformation of the terrestrial observatory coordinates into the GCRS (which to sufficient accuracy is equivalent to the BCRS). This is accomplished using an ITRS site position vector (s ITRS ) supplied to TEMPO2, transformed to the GCRS using algorithms consistent with the IAU 2000 Resolutions. The transformation proceeds as follows: where the matrices W, R and Q account for polar motion, Earth rotation, and the motion of the Earth spin axis in the ICRS. The latter term, in turn, consists of a product of a frame bias and the IAU2000B precession-nutation matrix (McCarthy & Luzum 2003). Details of the construction of these matrices are available in the IERS Conventions (McCarthy & Petit 2004) and in the documentation of the IAU Standards of Fundamental Astronomy (SOFA) software library 4 used by TEMPO2 to compute them. Polar motion and Earth rotation both exhibit unpredictable variations that must be corrected post facto using measured parameters. TEMPO2 uses the 'C04' series of Earth orientation parameters, provided by the IERS, as input to the SOFA routines. For the purposes of transforming the observing frequency to the barycentric frame (Section 2.6.1), the time-derivative of the Roemer delay is needed: where dots denote the time-derivative and the effects of proper motion and parallax are negligible. The barycentric velocity of the geocentre, r ⊕ , is provided by the Solar system ephemeris. The barycentric velocity of the site is dominated by the velocity imparted by Earth rotation. This is calculated usinġ Here a unit North vector (ẑ) is transformed using the inverse polar motion matrix and multiplied by the instantaneous Earth angular rotation rate, ω ⊕ , to give the angular velocity of the Earth in the ITRS. The cross-product of this with the site radius vector gives the tangential velocity, which is transformed to the GCRS in the same manner as s. Of the neglected terms in this approximation, the largest is the precessional velocity, at a negligible ∼10 −13 c.

Solar system dispersion
Radio signals encounter significant dispersion in the interplanetary medium, due to the electron content of the solar wind. The electron distribution follows a roughly r −2 form consistent with spherical expansion (Issautier et al. 1998). Integrating along the line of sight yields a dispersion measure (DM): where |r| is the heliocentric radius of the observatory, ρ is the pulsar-Sun-observatory angle, and n 0 is an overall scale parameter, which can be specified in TEMPO2 parameter files as NE1AU.
[It should be noted that equation (30) is an approximation that ignores the effects on the dispersion relation due to the bulk plasma velocity and its variation along the line of sight. However, there are other much more significant sources of error as discussed below.] The default value for n 0 is 4 cm −3 , consistent with recent measurements (Issautier et al. 1998(Issautier et al. , 2001Splaver et al. 2005) and significantly lower than the value of 9.961 cm −3 used by TEMPO. The computed DM is used to calculate the delay due to interplanetary dispersion: where f SSB is the observing frequency transformed to the barycentric frame (see Section 2.6.1). This correction is of limited accuracy. Recent studies with the Ulysses spacecraft have revealed extraordinary complexity in the solar wind. At solar minimum, the electron density (scaled by r 2 ) was found to show rms temporal variations of 10-50 per cent, depending on heliocentric latitude (Issautier et al. 1998), while at solar maximum the mean electron density increased but was subject to very strong modulation (Issautier et al. 2001). The default correction of equation (31) is therefore uncertain within at least a factor of 2, corresponding to a minimum of 130 ns of error at f = 1 GHz for a source located at an ecliptic pole, increasing to many microseconds for sources within a few degrees of the ecliptic.
Errors in the predicted Solar system dispersion delay are of concern not only because they add random noise to the timing residuals and fitted parameters, but also because a portion of the error term will exhibit annual periodicities that could corrupt the fitted pulsar position, proper motion and parallax (Splaver et al. 2005). It is also possible that other periodicities in this term could be misinterpreted as arising due to a planetary companion to the pulsar (Scherer et al. 1997, though see also Wolszczan et al. 2000). Previous attempts to improve the model have focused on adding fittable degrees of freedom (e.g. Cognard et al. 1996;Splaver et al. 2005); however the complex behaviour seen in the Ulysses observations indicates that the only way to adequately remove the effect is to measure it directly using multifrequency observations (see Section 2.6.1).

Shapiro delay
The Shapiro delay accounts for the time-delay caused by the passage of the pulse through curved space-time. The total delay obtained by summing over all the bodies in the Solar system (Backer & Hellings 1986): whereR is a unit vector in the direction of the pulsar, m j is the mass of body j, r j is a vector from body j to the telescope, and S 2 is a second-order correction discussed below. The first term may be re-written in the following more convenient form: where ψ j is the pulsar-telescope-object angle for the jth object. TEMPO2 includes the first-order Shapiro delay for the Sun, Venus, Jupiter, Saturn, Uranus and Neptune. For observations of sources very close to Solar system bodies, higher-order effects may become important. The largest of these is the geometrical excess path length due to gravitational light bending (Richter & Matzner 1983), which amounts to 9.1 ns for a ray with a trajectory grazing the solar limb. 5 From the derivation of Hellings (1986), assuming general relativity and re-parametrizing in terms of ψ for a very distant source, this term is given by For pulsar timing purposes, only the geometric delay due to the Sun needs to be considered.

Propagation through interstellar space
The propagation time of the signal from the pulsar (or binary system barycentre) to the SSB is conveniently broken into three contributions: the vacuum propagation delay, that is, the path length divided by the vacuum speed of light, and the excess delay due to dispersion and other frequency-dependent delays: The fourth term accounts for the special relativistic time-dilation between the SSB reference frame and the binary reference frame. 5 Hellings (1986) quoted a figure of 36 ns, referring to the calculation of Richter & Matzner (1983). However, that calculation incurred an extra pair of factors of two, one due to the fact that the calculation referred to a round trip, and the other due the fact that the reflector was located 1 au past the Sun, rather than being a very distant source.
where D is the so-called dispersion constant, and f SSB is the frequency of the radiation at the SSB. The barycentric frequency differs from the frequency at the observatory, due to a simple Doppler shift (∼10 −4 in magnitude) and a second-order relativistic correction (∼ 10 −8 ), plus the effects of gravitational redshift (∼10 −8 ). The first term is simply the time-derivative of the Roemer delay (the 'Roemer rate'), while the sum of the latter two is derivative of the Einstein delay (the 'Einstein rate'): Although D is the directly measurable quantity of equation (37), it is customary to speak in terms of the inferred column density of electrons along the line of sight, that is, the DM, DM = k D D. The value of k D used by TEMPO2 (and TEMPO) is k D ≡ 2.410 × 10 −16 cm −3 pc; see Paper I for further discussion.
Many pulsars exhibit temporal variation in the interstellar DM. TEMPO2 follows other software packages in modelling these variations by means of a Taylor series: where DM (n) is the (fittable) nth derivative of the DM, and t D is a user-specified epoch. This approach is limited to modelling variations of low to moderate complexity. At the level of accuracy demanded by current highprecision timing campaigns, significant unmodelled variations will appear on short time-scales, due to ionospheric (Section 2.5.2), interplanetary (Section 2.5.4) and interstellar electrons (Foster & Cordes 1990). An alternative and recommended course of action is to use the 'stride fit' feature of TEMPO2 (Paper I), to fit a different DM to every group of simultaneous (or contemporaneous) multifrequency TOA measurements. TEMPO2 also allows for the fitting of an additional frequency-dependent delay term: where k f and ζ define the scale and spectral index of the term. This term may be used in conjunction with the stride fit feature to model, for example, delays due to refraction in the ISM (Foster & Cordes 1990) or deviations from the cold plasma dispersion law (Phillips & Wolszczan 1992).

Vacuum propagation delay
The vacuum propagation delay is affected by the distance to the pulsar and the variation of this with time. The full geometric distance was derived in Section 2.3. The first term from equation (5) to be assigned to IS , |R 0 |/c, is in fact not measurable to even remotely sufficient precision, but is constant and so can be dropped from the timing formula. A side effect of this is that epochs of measured pulse frequencies, times of periastrons, glitches and so on are, in fact, retarded by the initial light propagation time, c|R 0 |. The remaining parts assigned to the interstellar delay are as follows: The first term on the right-hand side of equation (41) simply represents the displacement of the system in the initially radial direction. This receives potentially significant contributions from the radial velocity of the pulsar and the radial acceleration in the gravitational potential of the globular cluster or Galactic environment (Damour & Taylor 1991). The second term describes the so-called Shklovskii effect (Shklovskii 1970), whereby initially transverse motion attains a radial component due to the change in the direction of the line of sight caused by the transverse motion of the pulsar. That is, the pulsar moves tangentially, while a path of constant distance describes a circle centred on the SSB. The distance between the two is the excess propagation length. The third term corresponds to secular change in the size of the Shklovskii effect as the radial motion increases or decreases the distance, and hence changes the curvature of the line of constant distance just mentioned. An alternative interpretation of this term is that it is the cumulative effect on the Shklovskii term caused by the increase or decrease in the apparent proper motion, owing to the increasing transverse component of the initially radial motion (Section 2.5.3). To our knowledge, this effect was first noted by van Straten (2003), in the context of its manifestation as an apparent second derivative of the spin and orbital periods.
These effects are often neglected in timing formulae, because as simple constant or secularly increasing Doppler shifts, they are inseparable from the pulse and orbital periods and their derivatives. A common procedure is to correct the apparent period derivatives by estimating the contributions of transverse motion and gravitational acceleration and subtracting them to obtain the intrinsic value (e.g. Damour & Taylor 1991;Camilo, Thorsett & Kulkarni 1994). Alternatively, in certain binary systems where the true orbital period derivative is expected to be negligible, the measured derivative can be attributed entirely to these effects to obtain, in combination with a measurement of the transverse proper motion, an estimate of |R 0 | (Bell & Bailes 1996). Likewise, a measured second spin frequency derivative could, in principle, yield a measurement of the radial velocity (van Straten 2003) via the third term of equation (41).
As an alternative to modelling these effects via artificial contributions to spin and binary parameters, TEMPO2 can optionally include them directly in the following form: where v Shkdot , a , d Shk and a μ ≡ a ⊥ · μ ⊥ are fittable parameters that specify the radial velocity, radial acceleration, Shklovskii distance, and component of acceleration in the direction of the transverse proper motion. The terms involving acceleration arise due to the expansion of k in terms of proper motion, distance and acceleration (equation 6). The circular dependence of VP and t BB a is resolved by iteration as described in Section 2.5.3.
The parameters a and d Shk are inseparable, and highly covariant with a simultaneous change to the spin frequency derivative and orbital period derivative. It is therefore imperative that at least one of the first two parameters plus any one of the remaining three parameters are held fixed at some independently determined value when fitting. These parameters are included separately because each of them is potentially subject to external constraints. The transverse proper motion, μ ⊥ , is typically immune from this covariance due to its appearance in the annual proper motion term (Section 2.5.3). Likewise, v Shkdot and a μ are inseparable, and highly covariant with the spin frequency second derivative, so only one of these parameters should be allowed to vary in a fit. They are included separately because v may be determined from μ (Section 2.5.3) and one or more of the distance parameters, while a may be estimated from models of the Galactic rotation and gravitational potential. The first-order radial velocity term is discussed in the following section.

Transformation from SSB to BB coordinate time
The final term contributing to IS accounts for the special relativistic time-dilation owing to the relative velocities of the Solar system and BBs: where v is the relative velocity of the frames. This transformation is applied to the TOA of the pulse at the BB, so that neglecting terms of the order of v 4 /c 4 and higher. Here, the notation ES refers to the fact that the term is analogous to the Solar system and binary Einstein delays (Sections 2.5.1 and 2.7.4) and relates to the relative secular motion of the SSB and BB. This transformation combines with the much larger v /c term of equation (42) to effect a Doppler shift between the two frames: An unknown Doppler shift of this kind cannot be distinguished from a re-scaling of various parameters intended to refer to the BB frame, meaning that neither v nor v 2 can be measured via their contribution to . In fact, because the numerical values used for c and G are the same in all frames, the effect of neglecting ES is the same as a change of unit system involving re-scaling of the length mass and time units, and physical tests involving only 'measured' (wrong) parameters and these constants will remain valid (Damour & Deruelle 1986). In the timing model presented here (and in contrast to that of Damour & Deruelle 1986), v can, in principle, be deduced from parameters measured via other terms of equation (5). For this reason, it may be included by specifying non-fittable values for the initial radial velocity, v , and transverse speed, v ⊥ , in the parameter file, which correct for the Doppler effect via equation (42) and the following formulation of the relativistic second-order Doppler shift: However, it must be remembered that, in practice, the uncertainty of many parameters measured in the frame of the BB or pulsar will be dominated by the uncertainties in v and v ⊥ .

The effects of a binary companion
The effects of a binary companion are separated into several components: the geometric Roemer delay, which accounts for excess the vacuum light traveltime due to the Euclidean displacement of the pulsar, a pseudo-delay that accounts for the aberration of the radio beam by binary motion, the Einstein delay (combined gravitational redshift and special relativistic time-dilation in the pulsar frame) and Shapiro delay (gravitational time-dilation in the vicinity of the companion): The binary model used by TEMPO2 is based on those of Blandford & Teukolsky (1976), Damour & Deruelle (1986; hereafter DD), Taylor & Weisberg (1989), Wex (1998 unpublished contribution to TEMPO; see also Lange et al. 2001) and Wex (1998), with extra terms as described by Kopeikin (1995Kopeikin ( , 1996.

Roemer delay and Kopeikin terms
The variation in the distance of the pulsar due to binary motion contributes four significant terms to the full geometric path length of equation (5): The first term of equation (50) is the well-known first-order geometric delay due to the initially radial component of the displacement due to binary motion. The remaining terms are much smaller in magnitude and have only proven measurable in a few pulsars to date. Hereafter, we will refer to these collectively as the 'Kopeikin terms'.
The second term of equation (50) describes changes to the apparent viewing geometry of the orbital motion, due to the proper motion of the system (Arzoumanian et al. 1996;Kopeikin 1996). An approach taken in the past to account for this effect is to allow it to be absorbed in derivatives of the projected semimajor axis, x, and longitude of periastron, ω, with the measured values being converted into constraints on the position angle of ascending node, , and the inclination, i (e.g. Sandhu et al. 1997;Nice, Splaver & Stairs 2001). Alternatively, the binary model can be modified to explicitly include linear changes to x and ω as parametrized by the fittable parameters and i (e.g. van Straten 2004). Here, we take a third approach which we believe is cleaner: rather than defining the orbital parameters in a rotating reference frame, we state explicitly that the orbital parameters refer to the reference frame defined by the 'initial' geometry at time t BB , and hence account for the effect directly via the second term of equation (50).
The third term of equation (50) describes the so-called annual-orbital parallax (Kopeikin 1995). This may be interpreted either as the modulation of the Solar system Roemer delay (Section 2.5.3) due to the transverse orbital motion of the pulsar, or as the modulation of the binary Roemer delay (first term of equation 50) due to the transverse orbital motion of the Earth. As with the secular changes to the apparent viewing geometry, although the annual-orbital parallax has been modelled in the past by perturbing the orbital parameters with terms involving r (van Straten et al. 2001;Splaver et al. 2005), we choose to maintain a consistent reference frame for the orbital parameters and instead model the effect directly via the third term of equation (50).
The fourth term of equation (50) is the orbital parallax (Kopeikin 1995), which accounts for the fact that transverse binary motion alters the line of sight and thereby acquires an apparently radial component: this is the orbital equivalent of the Shklovskii effect (Section 2.6.2). To our knowledge, this is yet to be measured in any pulsar, and naturally cannot be absorbed in time-varying orbital parameters since the changes themselves occur on the orbital time-scale.
The contributions of the four terms of equation (50) are included in the timing model as follows: and where d AOP and d OP are estimates of the initial distance which can be either held fixed at some independently determined value (or omitted from the parameter file, to neglect the effect), tied to other fittable distance parameters (e.g. d p and/or d Shk ), or, given sufficient constraints on the relevant orbital parameters, fitted to obtain independent distance estimates. Likewise, μ ⊥B can be either tied to μ ⊥ , or omitted from the parameter file to neglect the effect. Note that these parameters, in fact, refer to the quantities as observed in a frame that is comoving with the BB, in contrast to corresponding distance and proper motion parameters in preceding sections. However, the transformation can be neglected for a fractional error of O(v 2 /c 2 ) < ∼10 −5 in those parameters, corresponding to changes in the time of emission calculation well below our 1-ns goal.

Post-Newtonian orbital kinematics
Calculation of the various orbital delays depends on a knowledge of the displacement b of the pulsar from the binary system barycentre, or at least its radial component b , depending on the demanded accuracy. TEMPO2 follows the generalized post-Newtonian treatment of DD for the calculation of b, with the addition of secular derivatives for the orbital period, eccentricity and projected semimajor axis after Taylor Weisberg (1989). Specifically, from equations (16) and (17) where |b| = a(1 − e r cos u), where P b0 andṖ b are the initial value of the orbital period and its derivative, ω 0 andω are initial longitude of periastron and the mean of its derivative, T 0 is the proper time of periastron, e 0 andė are the initial 'proper time eccentricity' and its derivative, and δ r and δ θ parametrize relativistic deformations of the orbit. The matrices of equation (54) in the right-to left-hand order account for the inclination of the orbit to the line of sight (i), rotation about the line of sight ( ), and projection from a radial-transverse basis to a frame rotationally aligned to the ICRS (Section 2.2). The matrices are chosen such that for i = 0, the angular momentum vector of the orbit is antiparallel toR 0 , and measures the position angle of the ascending node with respect toê 2 , in the sense of rotation intoê 1 . The definitions of i and therefore correspond to astronomical convention ifê 1 andê 2 are east and north vectors. This is the case by default in TEMPO2:ê 1 =α,ê 2 =δ. Alternatively, if ecliptic coordinates are in use (Section 2.5.3),ê 1 =λ andê 2 =β, so that position angles are measured counterclockwise from an ecliptic meridian. Note that the solution to Kepler's equation (56) involves the proper time of emission, t psr e . This depends on the equivalent coordinate TOA at the BB, which, in turn, is related to the proper time of reception at the observatory by all of the terms in the timing model, including those involving t psr e . This circular dependence is resolved by starting with t psr e = t BB a , computing the orbital delays, making an updated estimate of t psr e , and iterating this process until the change in the orbital delay is less than 100 ps.
Using the above expression for b and writing it in terms of its transverse and radial components, the geometric delay of equation (51) can be found explicitly. Beginning by expanding S ≡ |b|/a sin θ and C ≡ |b|/a cos θ, and where we have followed the appendix of DD in making a minor modification to the definition of e to simplify the expressions. Writing equation (54) in terms of the projections of b onê 1 ,ê 2 andR 0 as a function of S and C, b 1 = a (C sin + S cos cos i) , and b = aS sin i.
The basic Roemer delay may then be written as and x ≡ a/c sin i is the projected semimajor axis as a vacuum light traveltime, defined in terms of its initial value, x 0 (DD) and its derivative, x (Taylor & Weisberg 1989). For the Kopeikin terms, following equation (53), where SR = x t BB a − t pos [(μ 1 sin + μ 2 cos )C csc i + (μ 1 cos − μ 2 sin )S cot i] , and OP = cx 2 2d OP (C 2 csc 2 i + S 2 cot 2 i).
The proper motion components are defined in terms of fittable parameters that can optionally float independently to the annual proper motion parameters (Section 2.5.3), in order to allow the latter to be separated from the secular change of binary viewing geometry. For equatorial coordinates, μ 1 = μ αB and μ 2 = μ αB , while for ecliptic coordinates μ 1 = μ λB and μ 2 = μ βB . In agreement with astronomical convention, the parameter corresponds to the position angle of periastron in terms of eastward rotation from north, while the parameter i is defined such that i = 0 corresponds to a binary orbital angular momentum vector directed towards the observer. As pointed out by Splaver et al. (2005), these differ from the definitions of Kopeikin (1995Kopeikin ( , 1996. We note that some caution is required in this area owing to variations in definitions in the literature: for other examples of the use of the unconventional definition of i, see, for example, DD, Damour & Taylor (1992), Weisberg & Taylor (2002) and Stairs, Thorsett & Arzoumanian (2004).

Aberration
Owing to the relative transverse velocity of the pulsar and the Earth, the direction of the observer as seen from the pulsar differs from that which would be measured in a frame comoving with the observer. The direction is aberrated according to the Lorentz transformation that relates the two frames. Under the assumption that the source of pulses is the periodic rotation of an emission beam, any such change in the direction of the observer alters the rotational phase corresponding to radiation received at a given time. Following DD (equation 27), in TEMPO2 this is converted to an equivalent time-delay: where A and B are parameters related to the orientation of the pulsar spin axis and the size of the orbit (see equations 38 and 39, and section 3.2, of DD). As discussed by DD, these parameters are highly covariant with the other binary parameters, and may only be separated on the time-scale of geodetic precession. Note, here the relative motion of the observatory and the BB is neglected.

Post-Newtonian delays
In addition to the modified Roemer delay of equation (51) and the aberration pseudo-delay, under the post-Keplerian formalism of DD there are two additional delay terms: the Einstein delay and the Shapiro delay. The Einstein delay is the difference between the proper time of emission and the coordinate time in the quasi-inertial frame of the BB: This is due to the combined effect of gravitational redshift and time-dilation. Following DD and in contrast to the TEMPO2 treatment of the Solar system Einstein delay, the binary Einstein delay specifically excludes a linear term by scaling t psr in such a way that the orbital period is numerically the same in either time-scale. It is expressed in a theory-independent manner in terms of a timing model parameter, γ ; from equation (19) of DD: The Shapiro delay is caused by the curvature of space-time in the gravitational field of the binary companion. After equation ( Here the fittable parameters are the theory-independent 'range', r, and 'shape', s. Under general relativity, s = sin i and r = Gm 2 /c 3 , where m 2 is the mass of the binary companion. In TEMPO2, r is expressed in units of T = GM /c 3 , so that its value is numerically equal to the general relativistic (GR) companion mass in units of solar masses. (T is half the light traveltime across the solar Schwarzschild radius.) An alternative parametrization is available, where the parameter s is replaced by the 'SHAPMAX' parameter, z s ≡ −ln (1 − s), This derives from a modification to TEMPO designed to avoid difficulties in fitting highly inclined orbits near |sin i| = 1 (Kramer et al. 2006a). In this case, the above relation applies after substituting s = 1 − exp(−z s ).

Assuming general relativity
The post-Keplerian parameters of Sections 2.7.2 and 2.7.4 (Ṗ b ,ω, δ r , δ θ , γ, s, m 2 ) either can be specified directly, or can be derived from a reduced set of physical parameters based on the assumption that the predictions made by general relativity are correct. These parameters are m 2 and M ≡ m 1 + m 2 , the total system mass (in addition to the Keplerian parameters, P b0 , ω 0 , e, T 0 , and x 0 ). After equations (8.48)-(8.55) of Lorimer & Kramer (2005), anḋ Although these equations account fully for the general relativistic rates of change in orbital period and longitude of periastron, extra variation in these parameters is allowed for via the parametersṖ bx andω x : and Note thatω x is not intended to model the secular change in orbital viewing geometry due to proper motion: this is modelled by a separate term in the timing model, causing ω to be defined in a consistent, non-rotating reference frame (Section 2.7.1). Under the condensed parametrization of this section, the orientation of the spin axis may be separated from the orientation of the orbit normal vector to yield a physical parametrization of the aberration delay (DD). Following equations (8.56) and (8.57) of Lorimer & Kramer (2005): where ν is the pulsar spin frequency (Section 2.8), ξ is the inclination of the pulsar spin axis, and χ is its position angle, 6 relative to the position angle of ascending node, . Here χ and csc ξ are the fittable parameters. Note that because cos i is only ambiguously constrained by the Shapiro delay term (cos i = ± 1 − sin 2 i), χ too is subject to ambiguity. To evaluate equation (90), TEMPO2 assumes that cos i 0. Values of χ determined on this basis are therefore indistinguishable under the transformation χ → π − χ . Only if i is measurable via the Kopeikin terms can this degeneracy be broken. A second, unbreakable degeneracy exists due to the fact that aberration is only sensitive to the projection of the transverse orbital velocity on the spin axis: ξ → π − ξ .

Nearly circular orbits
For low-eccentricity orbits, ω and T 0 are highly covariant, which complicates the analysis of the uncertainty in these parameters. This problem can be avoided by parametrizing the motion in terms of the Laplace-Lagrange parameters, and κ ≡ e cos ω, (92) and the time of ascending node, T asc ≡ T 0 − ω/n (93) (Wex 1998, unpublished contribution to TEMPO as 'ELL1' binary model). Following Lange et al. (2001), the following approximation to equation (70) applies, to order e: where and κ = κ 0 +κ t psr e − T asc .
Here η 0 ,η, κ 0 ,κ and T asc replace e 0 ,ė, ω 0 ,ω and T 0 as fittable parameters, while δ θ and δ r are neglected. The original ELL1 model as presented by Lange et al. (2001) did not include the Kopeikin terms (Section 2.7.1). In order to include them, an expression for C (equation 66) is required. To the same accuracy as equation (95), we find neglecting a constant offset of 3aκ/2 (analogous to an offset of 3xη/2 dropped by Lange et al. 2001). Substitution into equations (73)-(75) yields the Kopeikin terms under this parametrization. Under the Laplace-Lagrange parametrization, the Shapiro delay is computed as follows (Lange et al. 2001): where either of the parametrizations (s, z s ) of sin i of Section 2.7.4 may be used. The aberration delay becomes (Wex 1998, unpublished): The Einstein delay is neglected if this parametrization is chosen.

Main-sequence companions
The orbits of pulsars with main-sequence binary companions show deviations from a Keplerian orbit attributed to spin-orbit coupling (Lai, Bildsten & Kaspi 1995). TEMPO includes two binary models for taking these into account: the BTJ model, which induces steps ('jumps') in several parameters at specified epochs (typically, at each periastron), and the MSS model, which includes orbital and secular variations to the inclination angle (Wex 1998). The insertion of jumps is straightforward and need not be discussed further. The MSS model is included in the TEMPO2 model by changing the manner in whichẋ alters x, and including second period derivatives in x and ω: and

Emulating other models
As noted in Section 2.7.2, the orbital kinematics represent an inverse problem, because the motion is parametrized by the pulsar proper time, which itself depends on the orbital propagation delay. The standard approach in TEMPO2 is to iterate until convergence. However, in order to reproduce results from the implementations of the Blandford & Teukolsky (1976) or DD models, a facility is provided to alter this behaviour. Both models begin with the zeroth-order approximation that t psr e = t BB a to compute the various delay terms, then add correction terms to the zeroth-order approximation to the combined Roemer and Einstein delay, 0 RE , which is expressed (after DD equations 46-48) as follows: and Note that the Kopeikin terms (equation 53) are neglected. The expansion provided by DD is wherê The approximation provided by Blandford & Teukolsky (1976) is equivalent to In addition, Blandford & Teukolsky (1976) neglected BS , BA , δ r and δ θ , and the orbital modulation of ω: that is, in contrast to equation (62), Keeping the DD formulation of ω but the BT formulation of the other terms allows for the emulation of the BT+ model of TEMPO (Damour & Taylor 1992). Additionally, it is possible to emulate the TEMPO implementation of the Laplace-Lagrange parametrization (Section 2.7.6), which made an approximation to the Roemer delay analogous to the DD treatment: where and 0 = n(t BB e − T asc ).

Pulse phase model
The final component of the timing model is the evolution of the phase of the pulse sequence, relative to the proper time t psr of the pulsar centre of mass (equation 7). In most cases, a simple Taylor series expansion is used: The frequency derivative terms [ν (n) ] are the fittable parameters, while the epoch, t P , at whichφ = ν is set by the user. Absolute phase alignment is achieved by means of φ 0 , which is not specified directly. Instead, in order to protect against changes in t P or the choice of barycentric time-scale, φ 0 is defined in terms of a reference TOA for a specified observing site and observing frequency. This TOA is chosen such that it is close to the centre of the time-span of the fitted data, and has a predicted pulse phase that is exactly an integer.
Note that under the assumption that the basic pulse train is a result of a rotating emission beam, the time of zero phase is, in fact, dependent on the position of the observer. Of all the sources of relative motion (Earth orbital, binary orbital and secular motion and acceleration), only the secular motion of the pulsar (or BB) is significant, given our accuracy goal of 1 ns. However, since the transverse speed is nearly constant, so is the rate of change of apparent pulse phase, meaning that the effect cannot be separated from the basic pulse frequency, ν. the effect is neglected in TEMPO2. Likewise, the emitted beam is subject to aberration due to the relative velocity of the pulsar; however, the constant component of this velocity contributes only a constant phase shift that is absorbed in φ 0 . On the other hand, the binary motion induces a variable aberration, which, for simplicity, is included in B as a pseudo-propagation delay A (Section 2.7.1).
Though the pulse frequency of most pulsars shows a continual gradual slow-down, attributed to magnetic braking of a spinning neutron star, certain pulsars exhibit anomalies known as 'glitches'. These events are generally well modelled by permanent increments to the phase, frequency and frequency first derivative, in addition to a frequency increment that later decays exponentially to zero: An arbitrary number of glitches may be modelled in this manner; the corresponding values of φ g being added to the basic Taylor series of equation (120) to give the predicted evolution of pulse phase with time.

Geometric propagation delay
Errors are introduced to the geometric part of the propagation delay in two ways: through the neglect of various terms of equation (5), and through approximations made in the computation of other terms. Cases of the latter are discussed below. The neglected terms are the O(|R 0 | −3 ) terms, and all those involving the expansion of the third, second, and fourth terms in the second set of parentheses, except for the secular and annual changes to the Shklovskii effect and the secular change to the annual proper motion. Based on measured parameters and a timing campaign spanning 20 yr, we have estimated the likely magnitude of the neglected terms for all pulsars in the current ATNF pulsar catalogue having spin periods shorter than 30 ms (slower pulsars being unlikely to yield high-precision timing measurements). The largest of the neglected terms are the change in annual proper motion due to transverse acceleration (∼1 ns for PSR B1620−26), second-order secular changes to the apparent binary orbital viewing geometry (∼600 ps for PSR J0437−4715), the O(|R 0 | −3 ) terms involving secular motion (∼600 ps for PSR B1257+12), the binary modulation of the Shklovskii effect (∼300 ps for PSR J0437−4715) and the second-order annual proper motion (∼100 ps for PSR J0437−4715). It is possible that pulsars yet to be discovered may yield larger errors than the examples listed here, if they are in wider binary orbits, are closer to the Earth, have larger proper motions, and/or experience greater secular acceleration. The current sample is not likely to be biased in relation to proper motion, and any bias in relation to distance will tend to select closer pulsars due to the dependence of detection on flux density. Highly accelerated pulsars are difficult to detect; however, the level of acceleration at which this becomes problematic is many orders of magnitude larger than that may be experienced at the Galactic or intracluster environment. Binary acceleration also hampers detection, but for the known population the lack of any strong dependence of companion mass on orbital separation means that this selection effect biases the sample towards large orbital separations rather than small. However, it is not inconceivable that the present sample has missed an entire class of millisecond pulsars that experience strong binary acceleration due to very massive companions in wide orbits. The discovery of such pulsars may necessitate the re-inclusion of some of the dropped terms of the geometric propagation delay.

Clock corrections
Clock corrections applied by TEMPO2 may suffer from time-transfer errors present in the clock offset files provided by the user, owing either to incorrect tabulated values, or to variations in the clock offset on time-scales shorter than the tabulation interval, over which linear interpolation is performed. In most instances, the dominant source of these errors will be the first step of correction, involving the local observatory clock, since observatory time standards and time-transfer systems are typically less accurate than those involved in GPS and UTC time keeping. As an example of the level at which time-transfer errors may occur, we note that differencing of two independent time-transfer systems at the Parkes Observatory leads to an error term with an rms deviation of 8 ns over the last three years (Fig. 1). Time transfer to UTC from GPS or one of the national time standards contributing to UTC also leads significant error, with a typical rms of ∼5 ns (Lewandowski et al. 2005).
In addition to clock correction errors, the 'corrected' TOA, as measured against a specified realization of TT will differ from that measured against an ideal TT, owing to instability in the contributing atomic clock systems. An indication of the level at which instability contributes can be obtained by means of the difference between the latest revision of the BIPM realization of TT, TT(BIPM2005), and the 2001 version of the same time-scale. As shown in Fig. 2, differences exist at the ∼10 ns level, even after allowing for the absorption of a fraction of the variation by other terms in the timing model.

Solar system Einstein delay
The time-dilation integral of Irwin & Fukushima (1999) (equation 9) includes only terms of the order of 1/c 2 , plus a mean rate correction for higher-order terms. Inaccuracies in the Solar system ephemeris lead to an estimated uncertainty of 100 ps in the time-integral, and 10 −16 in its derivative. The remaining higher-order terms have a maximum amplitude of 33 ps (Fukushima 1995). The effect on the spatial coordinates of the four-dimensional transformation relating the GCRS and BCRS is neglected. The neglected length contraction is of the order of 10 −8 s/c, corresponding to up to 200 ps of excess propagation time.

Atmospheric delays
Equation (11) predicts the hydrostatic zenith delay to an accuracy of 0.5 mm under conditions of hydrostatic equilibrium, with variations of up to 20 mm possible under extreme weather conditions (Davis et al. 1985). These correspond to timing errors of 1.7-67 ps × csc . In addition, the accuracy of NMF degrades at small elevation angles (Niell 1996), leading to ∼30 ps of error at = 5 • . Should the user not provide tabulated surface atmospheric pressure, the use of a constant value leads to errors of the order of 1 ns × csc . Likewise, lack of tabulated values for the ZWD results in error at the level of ∼1 ns× csc . As noted in Section 2.5.2, the dispersive component of the ionosphere causes propagation delays of the order of 7-130 ns ×(f/1 GHz) −2 , where f is the observing frequency. Such effects must be removed by fitting for a time-variable DM with the aid of multifrequency observations.

Solar system Roemer delay
The Roemer delay depends on the position of the observatory in the terrestrial frame, provided as input to TEMPO2. An accuracy of 30 cm is needed to meet the 1-ns criterion. In many cases, the position is known to an accuracy of a few centimetres, by means of VLBI or other geodesy techniques. In addition, most sites are subject to continental drift at a rate of a few centimetres per year, which is neglected by TEMPO2. This may be remedied in future when the sensitivity and time-span of observational data dictate it.
The transformation of the observatory position from the terrestrial to celestial reference frame depends on a knowledge of the Earth orientation parameters. Our requirement of 1-ns precision corresponds to maximum 10 mas of angular displacement on the surface of the Earth, or 600 μs in UT1, the time-scale tied to Earth rotation. 7 TEMPO2 uses the 'C04' Earth orientation parameter series of the IERS, which exceeds the required precision for all dates after 1980. The effects of diurnal and subdiurnal tides are removed from the series, resulting in an additional error of at most 0.7 mas in the location of the pole, and +50 −80 μs in UT1. The IAU2000B precession-nutation model used by TEMPO2 is accurate to 1 mas.
As noted in Section 2.5.3, TEMPO2 applies transformations to the input and output values of the Solar system ephemeris that are specific to the DE405 ephemeris. Should a different ephemeris be specified, it is important that the offset of the input time parameter, and the scale of input time and output spatial vectors match those of DE405, as specified by Irwin & Fukushima (1999). Likewise, it is important that the ephemeris coordinate frame is rotationally aligned with the ICRS, for while the gross effect of a frame rotation on r ⊕ ·R 0 could be absorbed in the fitted pulsar position, the resultant changes inR 0 will induce an error in s ·R 0 with a diurnal signature. For example, the coordinate frame of the earlier DE200 ephemeris is offset from the ICRS by ∼14 mas (Folkner et al. 1994), resulting in up to ± 1.4 ns of error.
Ultimately, the accuracy of the Roemer delay hinges on the chosen Solar system ephemeris. Preliminary investigations indicate that errors in the Solar system ephemeris may be significant and indeed measurable in pulsar timing data.

Solar system Shapiro delays
TEMPO2 includes the first-order Shapiro delay due to all bodies for which the excess delay exceeds 0.5 ns (see Paper I). The next largest effect is due to Mars, with 60 ps of delay. The second-order term is included only for the Sun. The next largest effect is due to Jupiter, with a maximum amplitude of 5 ps. Further terms appear for the Sun at the 500 ps level (Richter & Matzner 1983).

Interstellar propagation delays
The treatment of TEMPO2 is precise in the limit of cold plasma dispersion with negligible refraction. In practice, turbulence in the ISM leads to stochastic refraction and scattering of incident radio waves, giving rise not only to familiar scintillation effects, but also to frequency-dependent variable propagation delays (Foster & Cordes 1990). Such delays arise both due to differences in the propagation path length, and also due to induced errors in the computed Roemer delay, owing to the variations in apparent source direction. According to simulations conducted by Foster & Cordes (1990), the amplitude of these terms in the selected pulsars is of the order of hundreds of nanoseconds, for observations conducted at λ ∼ 20 cm. For future, more sensitive telescopes, these errors will almost certainly dominate over radiometer noise, unless steps are taken to correct for the effect or avoid it by observing at higher frequencies.

Transformation from SSB to BB frame
The transformation from the relativistic coordinate frames of the SSB to BB is modelled as a rotationless Lorentz transformation. This is a pure special relativistic transformation that assumes that the space-time metric tensors are completely determined by bodies in the Solar system and binary system, respectively. Under general relativity, this assumption will be violated by the presence of gravitational waves from external sources, causing unmodelled effects in the pulse TOAs. Indeed, the detection of the gravitational wave background is one of the main aims of contemporary pulsar timing campaigns. Jenet et al. (2005) conducted simulations of a number of different hypothetical data sets to test the detectability of a gravitational wave background at the level predicted by theory. A significant detection was predicted for several simulated data sets of a standard optimistically within the reach of current observing programmes, for which the minimum detectable systematic effect (as indicated by the rms measurement error divided by the square root of the number of independent observations) is one to several nanoseconds.

Binary orbital effects
The principal source of error in the formulation of the timing delay due to the presence of a binary companion is the use of equations of motion derived at only the first post-Newtonian order (Damour & Deruelle 1986). To date no complete timing formula accurate at the second post-Newtonian (2PN) order has been constructed, although the orbital motion has been derived by Wex (1995), including a proper-time representation that thereby provides the 2PN Einstein delay. The largest discrepancies are due to the 2PN contribution to the rate of periastron advance, and the effect of spin-orbit interaction on the periastron advance and orbital orientation (Damour & Taylor 1992), though these effects can be absorbed in changes toω andẋ. The 2PN corrections to the quasi-periodic Einstein delay are of the order 80 ns for PSR B1913+16 (Wex 1995). Additional unmodelled effects include the modification to the Shapiro delay due to the gravitomagnetic field of a spinning companion (Laguna & Wolszczan 1997), and the effect of gravitational light bending on the apparent rotational phase of the pulsar beam (Schneider 1990;Doroshenko & Kopeikin 1995), which cannot be separately measured (Wex & Kopeikin 1999). For double-neutron-star binaries, all of these effects are quite large compared to our 1-ns accuracy goal, especially for nearly edge-on orbits; however, it so happens 7 Errors in UT1 relate, in fact, to a position offset, and should not be confused with errors that appear directly in the timing model (such as clock corrections). This is why the precision demanded of UT1 is some 600 000 times higher than the required timing accuracy. that the binary evolutionary history of such systems rarely results in a pulsar that is spinning rapidly enough to provide timing measurements of high accuracy.
Only since the discovery of PSR J0737−3039A, a 22-ms pulsar in a highly inclined compact orbit with a second pulsar (Lyne et al. 2004), has the measurability of these effects been seriously considered. It is expected that the 2PN contribution toω will be measurable in coming years (Kramer et al. 2006b). Although this effect is easily accommodated in the existing binary model, we understand that efforts are underway to develop a fully consistent 2PN timing formula. An incorporation of this formula in TEMPO2 would be highly desirable for the assurance that the model will continue to be accurate at the level demanded by ever-decreasing measurement errors and potentially more extreme binaries discovered in coming decades.
It should be noted that, owing to approximations made in the orbital kinematics, the Laplace-Lagrange parametrization is only valid for orbits of low eccentricity. The error term is of the order of xe 2 (Lange et al. 2001).

Tropospheric delays
TEMPO makes no attempt to model non-dispersive atmospheric delays. These are dominated by the hydrostatic delay (Section 2.5.2), which is a mainly diurnal term of up to tens of nanoseconds in amplitude, depending on the elevation angle distribution.

Einstein delay
To compute a version of the time-dilation integral of equation (9), TEMPO uses the series approximation of Fairhead & Bretagnon (1990), since shown to be significantly in error (Fukushima 1995;Irwin & Fukushima 1999). An example of the discrepancy is shown in Fig. 3, which plots the difference in barycentric time as computed by TEMPO2 and TEMPO. The comparison is made in terms of the now-obselete Barycentric Dynamical Time (TDB, see below), as used by TEMPO. The difference consists mainly of a constant offset of ∼60 ns, an annual term of ∼2 ns amplitude, and a small noise-like contribution due to truncation of the internal Chebyshev representation used by TEMPO. The constant offset is immaterial since the zero-point of TDB has never been well defined. It arises for the most part due to the use by Irwin & Fukushima (1999) of an erroneous estimate of the offset between TDB and T eph , with an additional ∼6 ns due to software errors in early versions of the Fairhead & Bretagnon (1990) FORTRAN code, which only manifest on certain computer systems (P. Wallace, private communication). The annual term, probably due to inaccuracies in the Fairhead & Bretagnon (1990) series, will act to corrupt the source position at the 0.5-μas level, as fitted by TEMPO. It should also be noted that the time-integral actually computed by Fairhead & Bretagnon (1990) and Irwin & Fukushima (1999) includes a negative term in the integrand designed to give a long-term zero mean value of the integral. Unlike TEMPO2, TEMPO does not add a counteracting linear term to the result, causing an overall difference in the scale of the barycentric coordinate time-scale. In the case of TEMPO, the resultant time refers to the T eph time-scale (Section 2.5.1), which is approximately equivalent to the ill-defined IAU predecessor to TCB, known as TDB (see . The continued use of T eph or TDB as a TCB for pulsar timing is deprecated, since under this condition the scale of units of fitted parameters involving length and/or time differs subtly from the standard SI definition and fails to comply with IAU Resolutions: see Paper I.

Earth orientation
In its calculation of the observatory position in the barycentric frame, TEMPO uses the IAU 1976 precession (Lieske et al. 1977) and IAU 1980 nutation (Seidelmann 1982) theories, which are in error at the 50-mas level. More significantly, it also neglects polar motion, which amounts to up to ±300 mas, corresponding to up to ±35 ns of timing error, with a variable diurnal signature. The polar motion error is dominated by two periodicities, one with a period of one year (due to seasonal influences on oceans and the atmosphere), and the other with a period of 435 d (the Chandler wobble, due to free precession) -see Paper I. Depending on the site position, source direction and hour angle distribution, a part of the annual polar motion may enter the timing model as an annual term at the level of up to ∼10 ns, corrupting the fitted pulsar position at the level of a few microarcseconds.

Shapiro delays
TEMPO includes only the first-order Shapiro delay term due to the Sun. Terms of up to 180 ns for the planetary Shapiro delays, and 9 ns for the second-order term of the Solar Shapiro delay were omitted.

Dispersion
TEMPO includes some small errors in the transformation of the observing frequency to the barycentric reference frame. First, the computation of the site velocity (equation 27) neglected polar motion, precession and nutation. Polar motion changes the direction of the velocity vectoṙ x by up to 300 mas, altering its projection on the line of sight by up to ∼10 −10 c and corrupting the dispersion term IS of equation (37) by a fractional amount ∼10 −9.5 . More significantly, precession moves the vector by ∼20 arcsec yr −1 , leading to a fractional error in IS with an amplitude increasing at a rate of up to ∼2 × 10 −8 yr −1 . This will manifest as an apparent annual variation in the DM, with an amplitude equal to ∼DM · 10 −8 t yr −1 , where t is the time elapsed since Julian year J2000.0. For single-frequency data sets, the effect would be absorbed by the annual proper motion term, corrupting it at the ∼10 μas yr −1 level.
Secondly, TEMPO neglected to account for the Einstein rate (Section 2.6.1). The use by TEMPO of TDB instead of TCG (Section 4.2) removes a large positive mean from the Einstein rate, leaving residual rate errors at the 10 −9.5 level, mostly modulated with an annual periodicity. This corresponds to fractional errors in IS at the 10 −9 level.

Secular and binary motions
The full geometric time-delay of equation (5) contains several terms involving the displacement of the pulsar due to secular motion (k). Of these, TEMPO includes only those terms related both to secular and to annual motions. In doing so, it neglected direct Doppler shift (which includes both classical and relativistic terms) due to radial velocity and its change in time due to radial acceleration, and the Shklovskii effect and its change over time due to the radial motion and transverse acceleration. None of these omissions is necessarily problematic, because all can be absorbed in a redefinition of other measured parameters. In TEMPO2, they are provided directly so that advantage may be taken of additional constraints that may be available: for example, if the intrinsic spin orbital period derivative is known to be small, it can be held fixed at zero and the distance determined directly via the Shklovskii effect.
The public distribution of TEMPO also omits the Kopeikin terms; however, the secular changes to the viewing geometry may be absorbed in other parameters, while modified versions of TEMPO exist that include the annual-orbital parallax. The orbital parallax, not modelled in any known version of TEMPO, contributes a timing term of up to 30 ns in amplitude.

S U M M A RY A N D C O N C L U S I O N S
We have presented a model for pulsar pulse TOAs that exceeds all previous efforts in its accuracy. This timing model is the basis of TEMPO2, a new software package for pulsar timing. The goal of TEMPO2 to provide accuracy at the 1-ns level is largely satisfied by the model presented here, the exception being pulsars that are part of extremely relativistic binary systems. Further theoretical work on second-and higher-order relativistic effects will be required before progress can be made in this area. Fortunately, the pulsars for which the most-precise measurements are possible tend not to be part of such relativistic binaries.
Terrestrial clock instability, Solar system ephemeris errors and the local gravitational wave background stand as the likely sources of significant systematic error. Pulsar timing array campaigns should eventually be able to provide measurements of these effects, allowing for their removal from the residuals and more importantly providing a valuable insight into the causative physical processes.
With the next generation of radio telescopes, such as the Square Kilometre Array, achievable measurement error will be reduced in some cases to below 10 ns. For systematic effects, the timing model presented here should remain applicable in most cases, and will present a useful starting point for extension, if necessary. However, at this level of precision previously negligible stochastic effects will come into play. Possibly, the most problematic of these will be the effect of interstellar scattering.