Pulsar Timing

For a detailed look at pulsar timing and other pulsar observing techniques, see the Handbook of Pulsar Astronomy by Duncan Lorimer and Michael Kramer.

Pulsars are intrinsically interesting and exotic objects, but much of the best science based on pulsar observations has come from their use as tools via pulsar timing.  Pulsar timing is the regular monitoring of the rotation of the neutron star by tracking (nearly exactly) the times of arrival of the radio pulses. The key point to remember is that pulsar timing unambiguously accounts for every single rotation of the neutron star over long periods (years to decades) of time. This unambiguous and very precise tracking of rotational phase allows pulsar astronomers to probe the interior physics of neutron stars, make extremely accurate astrometric measurements, and test gravitational theories in the strong-field regime in unique ways.

For pulsar timing, astronomers "fold" radio data modulo the instantaneous pulse period $P$ or pulse frequency $f = 1 / P = d\phi / dt$.  Averaging over many pulses yields a high signal-to-noise average pulse profile. Although individual pulse shapes vary considerably, the shape of the average profile is quite stable. Typically, the average profile is correlated with a template or model profile so that a phase offset can be determined.  When multiplied by the instantaneous pulse period, that phase yields a time offset that can be added to a high-precision reference point on the profile (for example, the left edge of the profile based on the recorded time of the first sample of the observation) to create the time-of-arrival or TOA.  The precision with which a TOA can be determined is approximately equal to the duration of a sharp pulse feature (e.g., the leading edge) divided by the signal-to-noise ratio of the average profile.  It is usually expressed in terms of the width of the pulse features $W_f$ in units of the period $P$, the pulse period $P$, and the signal-to-noise ratio $SNR$ such that $\sigma_{\rm TOA} \propto W_fP / SNR$.  Therefore strong, fast pulsars with narrow pulse profiles provide the best arrival times.

In the nearly inertial frame of the Solar-system barycenter, the rotational period of a pulsar is nearly constant, so the time-dependent phase $\phi(t)$ of a pulsar can be approximated by a Taylor expansion
$$\phi(t) = \phi_0 + f(t-t_0) + \frac{1}{2}\dot f(t-t_0)^2 + \dots,$$ where $\phi_0$ and $t_0$ are arbitrary reference phases and times for each pulsar. The important thing about pulsar timing, though, is that the observed rotational phase difference between each of the TOAs must contain an integer number of rotations.  Since each TOA corresponds to a different time $t$, the parameters that we are fitting for, such as $f$ and $\dot f$, must result in a phase change between each pair of TOAs $i$ and $j$ that is an integer number of turns, or $\Delta\phi_{\rm ij} = n$ turns (1 turn = $2 \pi$ radians). Since all measurements are made with regard to the integrated pulse phase rather than the instantaneous pulse period, the precision with which astronomers can make long-term timing measurements can be quite extraordinary.
Example: With what precision can we determine the spin frequency $f$ of a pulsar using pulsar timing?

Since $f = d\phi / dt$ when $\phi$ is measured in turns, the precision is based on how precisely we can measure a change in phase $\Delta\phi$ over some time interval $\Delta T$.  Typically, $\Delta T$ is a long period of time (up to several tens of years for many pulsars now) over which a pulsar's phase has been tracked through regular monitoring.  $\Delta\phi$ is determined principally by the individual TOA precisions, although for some types of measurments a statistical component is important as well since precision improves as the number of measurements $N^{-1/2}$ if the random errors are larger than the systematic errors.

For the original millisecond pulsar B1937+21, the TOA precision is approximately 1 $\mu$s (which is a fractional error in phase of about $6\times10^{-4}$ turns, and it has been timed for over 25 years: $$\Delta f \sim \frac{\sigma_{\rm TOA}}{\Delta T} = \frac{6\times10^{-4}}{25\,{\rm yrs} \times 3.15\times10^7 {\rm s\,yr}^{-1}} = 8\times10^{-13}\,{\rm Hz}.$$

In order to measure $\phi(t)$ in this form, though, many corrections have to be applied to the observed TOAs first.  If we measure a pulse at our observatory on Earth at topocentric (topocentric means measured from a fixed point on the Earth's surface) time $t_{\rm t}$, we can correct this time to the time $t$ in the nearly inertial Solar-system center of mass or barycentric frame, which we assume to be the nearly the same as the time in the frame comoving with the pulsar.  Note that the measured pulse rates will differ from the actual pulse rates in the pulsar frame by the unknown Doppler factor resulting from the unknown line-of-sight pulsar velocity.
$$t = t_{\rm t} - t_{t_0} + \Delta_{\rm clock} - \Delta_{\rm DM} + \Delta_{{\rm R\odot}} + \Delta_{{\rm E\odot}} + \Delta_{{\rm S\odot}} + \Delta_{\rm R} + \Delta_{\rm E} + \Delta_{\rm S} .$$ As before, $t_{t_0}$ is a reference epoch,  $\Delta_{\rm clock}$ represents clock correction that accounts for differences between the observatory clocks and terrestrial time standards, $\Delta_{\rm DM}$ is the dispersion delay caused by the ISM, and the other $\Delta$ terms are delays from within the Solar System and, if the pulsar is in a binary, from within its orbit. The Roemer delay $\Delta_{\rm R\odot}$ is the classical light travel time across the Earth's orbit, with a magnitude of $\sim500\cos\beta$ s, where $\beta$ is the ecliptic latitude of the pulsar, and $\Delta_{\rm R}$ is the corresponding delay across the orbit of a pulsar in a binary or multiple system.  The Einstein delay $\Delta_{\rm E}$ accounts for the time dilation from the moving pulsar (and observatory) and the gravitational redshift caused by the Sun and planets or the binary companion, and the Shapiro delay $\Delta_{\rm S}$ is the extra time required by the pulses to travel through the curved space-time containing the Sun/planets/companions.  Errors in any of these parameters, as well as other parameters such as $f$, $\dot f$, and proper motion, give very specific systematic signatures in plots of the timing residuals, which are simply the phase differences between the observed TOAs and the predicted TOA times based on the current model parameters.

Example:  How (and how accurately) can we measure positions using pulsar timing?

Pulsar positions on the sky are determined by timing a pulsar over the course of a year as the Earth orbits the Sun and tracking the changing time delays (i.e. the Roemer delay) of pulses as the Earth moves.

The Roemer delay $\tau$ across the Solar System from a pulsar at ecliptic coordinates $\lambda$ (longitude) and $\beta$ (latitude) is: $$\tau \simeq 500\,{\rm s}\ \cos(\beta)\cos\left(\theta(t) + \lambda\right),$$ where $\theta(t)$ is the orbital phase of the Earth with respect to the vernal equinox.  This is an approximate time delay since we are assuming that the Earth's orbit is circular.

If there is an error in our position estimate, the individual position errors components $\Delta\lambda$ and $\Delta\beta$ cause a differential time delay $\Delta\tau$ to be present in the timing residuals with respect to the correct Roemer delay:
$$\Delta\tau \simeq 500\,{\rm s}\ \left[\cos(\beta+\Delta\beta) \cos\left(\theta(t) + \lambda + \Delta\lambda\right) - \cos(\beta)\cos(\theta(t) + \lambda)\right].$$
If the positional errors are small, such that we can use $\sin x\sim x$, $\cos x\sim 1$, and $\Delta\beta\,\Delta\lambda\sim0$, we can use trigonometric angle-sum identities and then simplify to get:
$$\Delta\tau \simeq -500\,{\rm s}\ \left[\Delta\lambda\cos(\beta)\sin(\theta(t) + \lambda) + \Delta\beta\sin(\beta)\cos\left(\theta(t) + \lambda\right)\right].$$
Comparing the trig identity $A\sin\left(\theta(t)+\phi\right) = A\cos\phi\sin\theta(t) + A\sin\phi\cos\theta(t)$ to the equation for $\Delta\tau$, we see that: $$A\cos\phi = -500\,{\rm s}\ \Delta\lambda\cos\beta$$ $$A\sin\phi = -500\,{\rm s}\ \Delta\beta\sin\beta,$$
and therefore: $$\Delta\lambda = -\frac{A\cos\phi}{500\,{\rm s}\ \cos\beta}$$ $$\Delta\beta = -\frac{A\sin\phi}{500\,{\rm s}\ \sin\beta},$$ where $A$ and $\phi$ are the amplitude and phase of the error sinusoid in the timing residuals.

When the pulsar is located near the ecliptic plane (with $\beta \sim 0$), $\cos\beta \sim 1$ and there is maximum timing leverage (and therefore minimum error) to determine $\lambda$. However, $\sin\beta \sim 0$ and so the errors on $\beta$ are huge.  If astrometric accuracy for pulsars near the ecliptic is necessary, VLBI positions are the best way to go.

For the timing fits themselves, the amplitude of the sinusoid $A$ is in time units (i.e. light travel time) and a timing fit will determine $A$ to an absolute precision $\Delta A$ approximately equal to the TOA uncertainty.  If that uncertainty is small, say 2 $\mu$s for a millisecond pulsar, and there are a large number of measurements (say $N = 16$) over the course of a year, the averaged phase errors (and therefore the errors on $A$) will be approximately $\sim 2\,\mu s / \sqrt{16} \sim 5\times10^{-7}$ sec.  The overall position errors for an MSP 30$^\circ$ off the ecliptic plane are approximately $$\Delta\lambda \sim \frac{5\times10^{-7}}{500\,\cos\beta} = 1\times10^{-9}\,{\rm radians}$$ and
$$\Delta\beta \sim \frac{5\times10^{-7}}{500\,\sin\beta} = 2\times10^{-9}\,{\rm radians}.$$ These correspond to errors in both directions of only a few hundred micro-arcsec!  Even normal pulsars with slow spin periods provide astrometric precisions typically of 0.1 arcsec or better.

pulsar timing
Figure 1:  Establishing a timing solution for an isolated pulsar.  In panel (e), you identify closely spaced days with unambiguous phase connection and fit for spin frequency.  In panel (f), you extend that phase connection until either RA or Dec errors dominate and then fit for it.  In panel (g), you fit for the other position component.  Finally, in panel (h), you fit for frequency derivative, which completes the timing solution.
timing examples
Figure 2: Pulsar timing examples.  Panel (a) shows a "good" timing solution with no unmodeled effects.  The sinusoidal ripple in Panel (c) indicates an error in position.  Panel (b) shows an error in the frequency derivative ($f = d\phi/dt$ so $\dot f = d^2\phi/dt^2$).  Panel (d) shows unmodeled pulsar proper motion.  From Lorimer and Kramer, 2005.

For binary pulsars, the pulsar Roemer delays comprise up to five Keplerian parameters: the projected semi-major axis $x\equiv a_1\sin i/c$, the longitude of periastron $\omega$, the time of periastron passage $T_0$, the orbital period $P_b$, and the orbital eccentricity $e$. Relativistic binaries may allow the measurement of up to 5 post-Keplerian (PK) parameters: the rate of periastron advance $\dot{\omega}$, the orbital period decay $\dot{P_b}$, the so-called relativistic $\gamma$ (i.e. the Einstein term corresponding to time dilation and gravitational redshift), and the Shapiro delay terms $r$ (range) and $s$ (shape).

PSR J0437-4715 timing
Table 1:  Millisecond pulsar timing example.  A timing ephemeris for the nearby MSP J0437$-$4715 by van Straten et al. 2001.  This is one of the best "timing" pulsars known (post-fit RMS timing residuals of $\sim$100 ns), and this measurement is one of the most accurate astrometric measurements ever made.  In addition, the timing accuracy allowed a fundamentally new test of general relativity.

In any theory of gravity, the five PK parameters are functions only of the pulsar mass $m_1$,  the companion mass $m_2$, and the standard five Keplerian orbital parameters. For general relativity, the formulas are: $$\dot\omega = 3 \left(\frac{P_b}{2\pi}\right)^{-5/3} (T_\odot M)^{2/3}\,(1-e^2)^{-1}$$ $$\gamma = e \left(\frac{P_b}{2\pi}\right)^{1/3} T_\odot^{2/3}\,M^{-4/3}\,m_2\,(m_1+2m_2)$$ $$\dot P_b = -\,\frac{192\pi}{5} \left(\frac{P_b}{2\pi}\right)^{-5/3} \left(1 + \frac{73}{24} e^2 + \frac{37}{96} e^4 \right) (1-e^2)^{-7/2}\,T_\odot^{5/3}\, m_1\, m_2\, M^{-1/3}$$ $$r = T_\odot\, m_2$$ $$s = x \left(\frac{P_b}{2\pi}\right)^{-2/3} T_\odot^{-1/3}\,M^{2/3}\,m_2^{-1}.$$
In these equations, $T_\odot\equiv GM_\odot/c^3 = 4.925490947\,\mu$s is the solar mass in time units, $m_1$, $m_2$, and $M\equiv m_1+m_2$ are in solar masses, and $s\equiv\sin i$ (where $i$ is the orbital inclination). If any two of these PK parameters are measured, the masses of the pulsar and its companion can be determined. If more than two are measured, each additional PK parameter yields a different test of a gravitational theory.

For the famous case of the Hulse-Taylor binary pulsar B1913+16, high-precision measurements of $\dot\omega$ and $\gamma$ were first made to determine the masses of the two neutron stars accurately. The Nobel-prize-winning measurement came with the eventual detection of $\dot P_b$, which implied that the orbit was decaying in accordance with general relativity's predictions for the the emission of gravitational radiation. The recently discovered double-pulsar system J0737$-$3039 is in a more compact orbit (2.4 hrs compared to 7.7 hrs for PSR B1913+16), which allows the measurement of all five PK parameters as well as the mass ratio $R$, giving a total of four tests of general relativity.  Kramer et al. (2006) showed that GR is correct at the 0.05% level and measured the masses of the two neutron stars to better than 1 part in $10^4$.

PSR B1913+16 MvMPSR B1913+16 orbital decay
Figure 3. Timing results for the Hulse-Taylor binary pulsar B1913+16.  The left panel shows the mass vs. mass plot for the pulsar and its companion neutron star.  The three lines correspond to the three measured post-Keplerian parameters.  The right panel shows the periastron shift caused by the decay of the orbit via emission of gravitational radiation.  The detection of gravitational radiation resulted in a Nobel prize for Hulse and Taylor.  (Figure provided by J. Weisberg).
PSR J0737-3039 MvM
Figure 4:  PSR J0737$-$3039 mass vs. mass diagram.  As in Figure 3, the diagram shows lines corresponding to the post-Keplerian parameters measured for the system.  In this case, though, six parameters were measured, including the mass ratio R since both neutron stars are pulsar clocks.  These measurements have tested GR to ~0.05% (Kramer et al. 2006).