Binary and Millisecond Pulsars

Abstract

Our knowledge of binary and millisecond pulsars has greatly increased in recent years. This is largely due to the success of large-area surveys which have brought the known population of such systems in the Galactic disk to around 50. As well as being interesting as a population of astronomical sources, many pulsars turn out to be superb celestial clocks. In this review we summarise the main properties of binary and millisecond pulsars and highlight some of their applications to relativistic astrophysics.

1 Introduction

In the 30 years since the discovery of pulsars, rapidly rotating highly magnetised neutron stars, the study of these fascinating objects has resulted in many applications in Physics and Astronomy. Striking examples have been the confirmation of the existence of gravitational radiation [150] as predicted by Einstein’s general theory of relativity [146, 147] and the first detection of an extra-solar planetary system [166, 122]. Currently, many new binary systems containing neutron stars are being discovered as a result of the latest generation of pulsar surveys. This review is concerned primarily with some of the results and spin-offs from these surveys which will be of particular interest to the Relativity community. The surveys themselves have been extensively reviewed by several authors [92, 27, 28] to which the interested reader is referred for further details.

By way of introduction, and to make the review fairly self-contained, we begin with an overview of the pulsar phenomenon (§2). This includes a brief review of the key observed population properties as well as a summary of the main theories concerned with the origin and evolution of pulsars. In §3, we review present understanding of the Galactic population of pulsars, discussing selection effects in the major surveys (§3.1), and techniques used to correct the observed sample (§3.2). These studies lead to robust calculations of the total number of normal and millisecond pulsars (§3.3) and neutron star binaries (§3.4) in the Galaxy and have implications for the detection of gravitational radiation from these systems. Perhaps the most important application of pulsars to relativity is the high precision science known as “pulsar timing” discussed in §4. Here it is seen that some pulsars are exceptional celestial clocks (§4.4). One application, a sensitive detector of long-period gravitational waves, is discussed in §5. Finally, in §6, we outline likely areas where progress may be made in the future.

2 The Pulsar Phenomenon

Many of the basic observational facts about radio pulsars were established shortly after their discovery [65] by Jocelyn Bell and Anthony Hewish in 1967. In the intervening years, theoretical and observational progress has flourished. Although there are many remaining questions, particularly about the emission mechanism, the basic model has long been established beyond all reasonable doubt viz: Pulsars are rapidly rotating, highly magnetised neutron stars formed during the supernova explosions of massive (∼ 5–10M_⊙) stars. In the following, we discuss the basic observational properties most relevant to this review.

2.1 The lighthouse model

The so-called “lighthouse model” used to explain the basic pulsar phenomenon is demonstrated in the animated image shown in Fig. 1. As the neutron star spins, charged particles are accelerated out along magnetic field lines in the magnetosphere (depicted by the light blue cones). This acceleration causes the particles to emit electromagnetic radiation, most readily detected at radio frequencies as a sequence of observed pulses produced as the magnetic axis (and hence the radiation beam) crosses the observer’s line of sight each rotation. The period of the pulses is simply the rotation period of the neutron star. The position of the tracker ball on the idealised pulse profile in the animation shows precisely the relationship between observed intensity and rotational phase of the neutron star.

Figure 1

mpg-Movie (146 KB) The rotating neutron star (or “lighthouse”) model for pulsar emission(For video see appendix).

Neutron stars are extremely stable rotators. They are essentially equivalent to large celestial flywheels with moments of inertia ∼ 10⁴⁵ g cm². The rotating neutron star model, which was independently developed by Pacini and Gold in 1968 [112, 55] correctly predicts that the pulse period should gradually increase as the outgoing radiation carries away rotational kinetic energy. When a period increase of 36.5 ns per day was measured for the pulsar in the Crab nebula [124], Gold showed that a rotating neutron star with a large magnetic field must be the dominant energy supply to the nebula and the model became universally accepted.

2.2 Pulse Profiles

Pulsars are weak radio sources. Mean flux densities, usually quoted in the literature at a radio frequency of 400 MHz, vary between 1 and 100 mJy (1 Jy = 10^-26 Wm^-2 Hz^-1). This means that the addition of many thousands of pulses is required in order to produce a discernible profile. A remarkable fact is that, although the individual pulses vary quite dramatically from pulse to pulse, at any particular observing frequency the integrated profile is very stable. The pulse profile can thus be thought of as a finger print of the emission beam.

This process is demonstrated in the animation shown in Fig. 2 — a sequence of consecutive single pulses for one of the brightest pulsars PSR B0329+54¹. This pulsar is seen in the animation to stabilise into its characteristic 3-component form after the summation of about 10 seemingly erratic single pulses. This property is of key importance in pulsar timing measurements (§4).

Figure 2

mpg-Movie (13.1 KB)Single pulses from PSR B0329+54.(For video see appendix).

A sample of pulse profiles is presented in Fig. 3 in which a great diversity in profile morphology can be seen.

Figure 3

A variety of integrated pulse profiles taken from the available literature. References: (a,b,d,f: [56]); (c: [10]); (e,g,i: [80]); (h: [14]). Each profile represents 360 degrees of rotational phase. These, and over 2600 other pulse profiles are available on-line [148].

2.3 The Pulsar Distance Scale

From the sky distribution shown in Fig. 4 we see straight away that pulsars are strongly concentrated along the Galactic plane, indicating that they populate the disk of our Galaxy.

Figure 4

The sky distribution of pulsars in Galactic coordinates. The plane of the Galaxy is depicted by the central horizontal line with the centre of the Galaxy in the middle. Bold points denote pulsars associated with supernova remnants.

Direct measurements of the distances to pulsars are notoriously difficult to obtain. There are three basic techniques: neutral hydrogen absorption, trigonometric parallax (measured either with an interferometer or through pulse time-of-arrival techniques) and from associations with objects of known distance (i.e. supernova remnants, globular clusters and the Magellanic Clouds). Together, these provide measurements of (or limits on) the distances to over 100 pulsars. For an excellent review of these measurements and their implications, see [161].

In the absence of a direct measurement, the distances to most pulsars can be estimated from an effect known as pulse dispersion which arises from the fact that the group velocity of the pulsed radiation through the ionised interstellar medium is frequency dependent: Pulses emitted at higher radio frequencies travel faster through the interstellar medium, arriving earlier than those emitted at lower frequencies. The delay Δt in arrival times between a high frequency ν_hi (MHz) and a low one ν_lo (MHz), can be shown [95] to be

$$\Delta t = 4.15 \times {10^6}\;{\rm{ms}}\; \times \;\left( {\nu _{{\rm{lo}}}^{ - 2} - \nu _{{\rm{hi}}}^{ - 2}} \right)\; \times \;{\rm{DM,}}$$

((1))

where the dispersion measure DM (cm^-3 pc) is the integrated column density of free electrons along the line of sight:

$${\rm{DM}}\;{\rm{ = }}\;\int_0^d {{n_{\rm{e}}}} \;dl\;{\rm{pc}}\;{\rm{c}}{{\rm{m}}^{ - 3}}.$$

((2))

Here, d is the distance to the pulsar (pc) and n_e is the free electron density (cm^-3). Pulsars at large distances have higher column densities and therefore larger DMs than those pulsars closer to Earth so that, from equation 1, the dispersive delay across the bandwidth is greater. Thus, from a measurement of the delay across a finite bandwidth we infer the DM. Hence, the distance can be estimated from a model of the Galactic distribution of free electrons, n_e. Taylor & Cordes [144] have developed such a model calibrated from the sample of pulsars with independently known distances plus measurements of interstellar scattering for various Galactic and extragalactic sources. The model appears to be free of large systematic trends and can be used to provide distance estimates with an uncertainty of ∼ 30%.

2.4 Normal and Millisecond Pulsars

The present public-domain catalogue, available on-line at Princeton [121], contains up-to-date parameters for some 706 pulsars. Most of these are “normal” in the sense that their pulse periods are of order one second and, with few exceptions, are observed to increase secularly at rates of typically 10^-15 s/s. A growing fraction of the observed sample are the so-called “millisecond pulsars”, which have spin periods primarily in the range 1.5 and 30 ms and have rates of slowdown < 10^-19 s/s. The first millisecond pulsar discovered, PSR B1937+21 [6], with a period of 1.5578 ms, remains the most rapidly rotating neutron star presently known to man.

A very useful means of demonstrating the distinction between these two classes is the “P - P” diagram — a logarithmic scatter plot of the observed pulse period versus the period derivative. This is shown in Fig. 5, in which normal pulsars occupy the majority of the upper right hand part of the diagram, whilst the millisecond pulsars reside in the lower left hand part of the diagram.

Figure 5

Pulse periods plotted against period derivatives for a sample of 639 pulsars. Filled red stars denote those pulsars known to be associated with supernova remnants. The blinking blue circles highlight those pulsars which are members of binary systems.

The differences in P and ̇P imply different ages and surface magnetic field strengths. By treating the pulsar as a rotating magnetic dipole, one may show that the surface magnetic field strength is proportional to $\sqrt {P\dot P}$[102]. Lines of constant magnetic field strength are drawn on Fig. 5, together with lines of constant characteristic age (τ = P/(2̇P)). From these, one infers typical magnetic fields and ages of 10¹² G and 10⁷ yr for the normal pulsars and 10⁸ G and 10⁹ yr for the millisecond pulsars respectively.

A very important additional difference between normal and millisecond pulsars is the presence of an orbiting companion. Presently, about 7% of all known pulsars are members of binary systems. Timing measurements (§4) place useful constraints on the masses of the companions which, supplemented by observations at other wavelengths, tell us a great deal about their nature. The present sample of orbiting companions are either white dwarfs, main sequence stars, or other neutron stars. A notable additional hybrid system is the so-called “planets pulsar” PSR B1257+12 — a 6.2 ms pulsar accompanied by at least two ∼ Earth-mass bodies [166, 122]. Orbital companions are much more commonly observed around millisecond pulsars (∼ 80% of the observed sample) than the normal pulsars (< 1%). The sample of binary pulsars is delineated in Fig. 6 as orbital eccentricity against mass of the companion. The dashed line serves merely to guide the eye in this figure. From this we note that the binary systems below the line are those with low-mass companions (< 0.7M_⊙ — predominantly white dwarfs) and essentially circular orbits: 10^-5 < e < 0.1. The binary pulsars with high-mass companions (> 1M_⊙ — neutron stars or main sequence stars) are in much more eccentric orbits: 0.2 > e > 0.9 and lie above the line.

Figure 6

Companion mass versus orbital eccentricity for the sample of binary pulsars.

The presently favoured model to explain the formation of the various types of systems has been developed over the years by a number of authors [21, 49, 138, 1] and can be qualitatively summarised as follows: Starting with a binary star system, the neutron star is formed during the supernova explosion of the initially more massive star which has an inherently shorter main sequence lifetime. From the virial theorem, it follows that the binary system gets disrupted if more than half the total pre-supernova mass is ejected from the system during the explosion. In addition, the fraction of surviving binaries is affected by the magnitude and direction of the impulsive “kick” velocity the neutron star receives at birth [66, 7]. Those binary systems that disrupt produce a high velocity isolated neutron star and an OB runaway star [22]. Over the next 10⁷ yr or so after the explosion, the neutron star may be observable as a normal radio pulsar spinning down to periods of several seconds or longer. The high disruption probability explains why so few normal pulsars have companions.

For those few binaries that remain bound, and in which the companion is sufficiently massive to evolve into a giant and overflow its Roche lobe, the old spun-down neutron star can gain a new lease on life as a pulsar by accreting matter and therefore angular momentum from its companion [1]. The term “recycled pulsar” is often used to describe such objects. During this accretion phase, the X-rays liberated by heating the infalling material onto the neutron star mean that such a system is expected to be visible as an X-ray binary system. Two classes of X-ray binaries relevant to binary and millisecond pulsars exist, viz: neutron stars with high-mass or low-mass companions. For a detailed review of the X-ray binary population, including systems likely to contain black holes rather than neutron stars, the interested reader is referred to [19].

The high-mass companions are massive enough to explode themselves as a supernova, producing a second neutron star. If the binary system is lucky enough to survive this explosion, it ends up as a double neutron star binary. The classic example is PSR B1913+16 [68], a 59 ms radio pulsar with a characteristic age of ∼ 10⁸ yr which orbits its companion every 7.75 hr [146, 147]. Other examples include PSR B2303+46, a 1 s pulsar with a characteristic age of ∼ 107.5 yr in a 12 day orbit [141]. Within the framework of this formation scenario, PSR B1913+16 is an example of the older, first-born, neutron star that has subsequently accreted matter from its companion. PSR B2303+46 is, on the other hand, likely to be the younger, second-born neutron star in its binary system. As we shall see (§3.4), double neutron star binary systems are very rare in the Galaxy, a direct indication that the majority of binary systems get disrupted by the exploding star.

Although no such system has yet been found in which both neutron stars are visible as radio pulsars, timing measurements (§4.3) show that the companion masses are ∼ 1.4 M_⊙ — exactly that expected for a neutron star [135]. In addition, no optical counterparts are seen. Thus, we conclude that these unseen companions are neutron stars that are either too weak to be seen/no longer active as radio pulsars or their emission beams do not intersect our line of sight. The two known radio pulsars with main sequence companions may well represent the “missing link” between high-mass X-ray binaries and double neutron star systems [73, 76].

By definition, the companions in the low-mass X-ray binaries evolve and transfer matter onto the neutron star on a much longer time-scale, spinning it up to periods as short as a few ms [1]. This model has recently gained strong support from the discoveries of kHz quasi-periodic oscillations in a number of low-mass X-ray binaries [162], as well as Doppler-shifted 2.49 ms X-ray pulsations from the transient X-ray burster SAX J1808.4-3658 [163, 30].

At the end of the spin-up phase, the secondary sheds its outer layers to become a white dwarf in orbit around a rapidly spinning millisecond pulsar. Presently ∼ 10 of these systems have compelling optical identifications of the white dwarf companion [13, 15, 91, 90]. Perhaps the best example is the white dwarf companion to the 5.25 ms pulsar J1012+5307 [110, 87]. This 19^th magnitude white dwarf is bright enough to allow precise measurements of its surface gravity, as well as the Doppler shifts of its spectral lines as it moves in its orbit [159].

2.5 Pulsar Velocities

Pulsars have long been known to have space velocities at least an order of magnitude larger than those of their main sequence progenitors which have typical values between 10 and 50 km s^-1. The first direct evidence for this came from optical observations of the Crab pulsar 1968 [155], showing that the neutron star has a velocity in excess of 100 km s^-1. Proper motions for about 100 pulsars have subsequently been measured largely by radio interferometric techniques [94, 9, 50, 60]. These data imply a broad velocity spectrum ranging from 0 to 1000 km s^-1 [96].

Such large velocities are perhaps not surprising, given the violent conditions which form neutron stars. Shklovskii [137] demonstrated that, if the explosion is only slightly asymmetric, an impulsive “kick” velocity of up to 1000 km s^-1 is imparted to the neutron star. As a result, a newly formed high-velocity pulsar quickly leaves its birth site close to the Galactic plane and on average it migrates to higher Galactic latitudes. This effect is seen most dramatically in Fig. 7, a dynamical simulation of the orbits of 100 neutron stars in a model of the Galactic gravitational potential.

Figure 7

mpg-Movie (213 KB) A simulation following the motion of 100 pulsars in a model gravitational potential of our Galaxy for 200 Myr. The view is side-on i.e. the horizontal axis represents the galactic plane (30 kpc across) whilst the vertical axis represents ±10 kpc from the plane. This snapshot shows the initial configuration of young neutron stars.(For video see appendix)

Using the proper motion data, recent studies have demonstrated that the mean birth velocity of normal pulsars is ∼ 450 km s^-1 ([96, 86, 36, 54]; see, however, also [61, 59]). This is significantly larger than the millisecond pulsars — recent studies suggest that their mean birth velocity is likely in the range ∼ 80 – 140 km s^-1 [83, 35, 99]. The main reason for this difference surely lies in the fact that about 80% of the millisecond pulsars are members of binary systems (§2.4) which couldn’t have survived if the neutron star had received a substantial kick velocities.

2.6 Going further

This section has served as an introduction to the basic aspects of pulsar astronomy relevant to subsequent sections of this review. For those wishing to delve deeper, two excellent graduate-level monographs are available: Pulsars by Manchester & Taylor [102] and Pulsar Astronomy by Lyne & Smith [95]. Those wishing to approach the subject from a more theoretical viewpoint are advised to read Michel’s The Theory of Neutron Star Magnetospheres [105] and The Physics of the Pulsar Magnetosphere by Beskin, Gurevich & Istomin [18]. Our summary of evolutionary scenarios is meant primarily as an introduction to the vast body of literature on this subject; the interested reader is referred to the excellent review by Bhattacharya & van den Heuvel [19] for further insights, as well as to recent observationally-based reviews concerning pulsar velocities [8, 84]. A useful up-to-date list of binary and millisecond pulsar parameters is available in a recent review by Camilo [27].

The list of pulsar-related resources available on the Internet is continually becoming more extensive and useful. Some of the more relevant links known to me are contained throughout this review. Good starting points for pulsar-surfers are the pages maintained by research groups at Bonn [106], Jodrell Bank [70] and Princeton [121].

3 The Galactic Pulsar Population

Soon after the discovery of pulsars, it was realised that the observed is heavily biased towards the brighter, nearby objects and probably represents only the tip of the iceberg of a much larger underlying population [57]. The extent to which the sample is incomplete is well demonstrated by the cumulative number distribution of pulsars as a function of distance projected onto the Galactic plane shown in Fig. 8. Here the observed distribution (solid line) is compared to the expected distribution for a uniform disk population in which there are no such selection effects (dashed line). We see that the observed sample becomes strongly deficient in terms of the number of sources for distances beyond a few kpc. This incompletion is a result of a number of different selection effects present in pulsar searches which we now discuss in some detail.

Figure 8

Cumulative number of observed pulsars (solid line) as a function of projected distance from the Sun. The dashed line shows the expected distribution for a hypothetical uniform disk Galaxy with a radius of 15 kpc in which all pulsars can be seen.

3.1 Selection Effects in Pulsar Searches

The most prominent selection effect at play in the observed pulsar sample is the inverse square law, i.e. for a given luminosity², the observed flux density falls off as the inverse square of the distance. This results in the observed sample being dominated by nearby and/or bright objects. Beyond distances of a few kpc from the Sun, the apparent flux density of most pulsars falls below the flux thresholds S_min of the a survey. Following [45], we may parameterise the survey threshold by:

$${S_{{\rm{min}}}} = \frac{{\beta \;{\rm{SN}}{{\rm{R}}_{{\rm{min}}}}\;\left( {{T_{{\rm{rec}}}} + {T_{{\rm{sky}}}}} \right)}}{{G\sqrt {{N_p}\Delta \nu \;{t_{{\rm{int}}}}} }}\sqrt {\frac{W}{{P - W}}} .$$

((3))

In this expression, β is a correction factor∼ 1.3 which reflects losses to hardware limitations, SNR_min is the threshold signal-to-noise ratio (typically 7–10), T_rec and T_sky are the receiver and sky noise temperatures, G is the antenna gain, N_p is the number of polarisations observed, Δν is the observing bandwidth, t_int is the integration time, W is the observed pulse width and P is the pulse period.

It follows from this equation that the sensitivity decreases as W/(P - W) increases. Also note that if W > P, the pulsed signal is smeared into the background emission and is no longer detectable. The observed pulse width W is in fact broader than the intrinsic value for a number of reasons: finite sampling effects; pulse dispersion, as well as scattering due to the presence of free electrons in the interstellar medium. As discussed above, the dispersive smearing scales as Δν/ν³, where ν is the observing frequency. This can largely be removed by dividing the pass-band into a number of channels and applying successively longer time delays to higher frequency channels before summing over all channels to produce a sharp “de-dispersed” profile³. The smearing across the individual channels, however, still remains and becomes significant at high dispersions when searching for short-period pulsars. Multi-path scattering results in a one-sided broadening due to the delay in arrival times which scales roughly as ν^-4, which can not be removed by instrumental means.

Dispersion and scattering become most severe for distant pulsars in the inner Galaxy as the number of free electrons along the line of sight becomes large. The strong frequency dependence of both effects means that they are considerably less of a problem for surveys at observing frequencies > 1400 MHz [34, 72] compared to the usual 400 MHz search frequency. An added bonus for such observations is the reduction in T_sky, since the spectral index of the non-thermal emission is about -2.8 [81]. Pulsars themselves have steep radio spectra. Typical spectral indices are -1.6 [100, 88], so that flux densities are an order of magnitude lower at 1400 MHz compared to 400 MHz. Fortunately, this can usually be compensated somewhat by the use of larger receiver bandwidths at higher radio frequencies. For example, the 1380 MHz system at Parkes has a bandwidth of 270 MHz compared to their 430 MHz system, where 32 MHz is available.

In the past, the main disadvantage in surveying at high frequencies has been the sky coverage rate which scales with the solid angle of the telescope beam. The current generation of high-frequency pulsar searches at Parkes and Jodrell Bank tackles this problem by installing multi-beam receivers in these telescopes. At Parkes, a 13 beam system [113] has been installed for use in neutral hydrogen surveys. The system is also being used for pulsar searches and can cover the sky at the same rate as the recent Parkes low-frequency 430 MHz survey [103, 99]. Together with a 4-beam system being installed at Jodrell Bank, the sensitivity of these systems is about 7 times better than previous surveys at 1400 and 1520 MHz [34, 72] and should thus discover several hundred new pulsars. Indeed, the Parkes survey has recently begun and has already discovered over 100 new pulsars. For an update, and more information, see Fernando Camilo’s multibeam page [26].

3.2 Correcting the observed pulsar sample

In order to decouple the selection effects from the observed sample, after Phinney & Blandford and Vivekanand & Narayan [118, 160], we define a scaling function ξ as the ratio of the total Galactic volume weighted by pulsar density to the volume in which a pulsar could be detected by the surveys:

$$\xi (P,L) = \frac{{\int {\int_{{\rm{Galaxy}}} {\;\sum {(R,\;z)\;R} } } \;dR\;dz}}{{\int {\int_{P,L} {\;\sum {(R,\;z)\;R\;dR\;dz} } } }}.$$

((4))

In this expression, ∑(R, z) is the assumed pulsar distribution in terms of galactocentric radius R and height above the Galactic plane z. Note that ξ is primarily a function of period P and luminosity L such that short period/low-luminosity pulsars have smaller detectable volumes and therefore higher ξ values than their long period/high-luminosity counterparts.

This scaling function can be used to estimate the total number of active pulsars in the Galaxy. In practice, this is achieved by calculating ξ for each pulsar separately using a Monte Carlo simulation to model the volume of the Galaxy probed by the major surveys [107]. For a sample of N_obs observed pulsars above a minimum luminosity L_min, the total number of pulsars in the Galaxy with luminosities above this value is

$${N_{{\rm{gal}}}} \approx \sum\limits_{i = 1}^{{N_{{\rm{obs}}}}} {\frac{{{\xi _i}}}{{{f_i}}}} .$$

((5))

Monte Carlo simulations of the pulsar population incorporating the aforementioned selection effects have shown this approximation to be reliable [85]. The factor f in this expression, known as the “beaming factor”, is the fraction of 4π steradians swept out by the radio beam during one rotation. Thus f gives the probability that the beam cuts the line-of-sight of an arbitrarily positioned observer. A naïve estimate of f is 20%; this assumes a beam width of ∼ 10° and a randomly distributed inclination angle between the spin and magnetic axes [145]. Observational evidence suggests that shorter period pulsars have wider beams and therefore larger beaming fractions than their long-period counterparts [109, 97, 20, 142]. It must be said, however, that a consensus on the beaming fraction-period relation has yet to be reached. This is shown below where we compare the period dependence of f as given by a number of models.

Figure 9

Beaming fraction plotted against pulse period for four different beaming models: Tauris & Manchester 1998 (TM88; [142]), Lyne & Manchester 1988 (LM88; [97]), Biggs 1990 (JDB90; [20]) and Narayan & Vivekanand 1983 (NV83 [109]).

Adopting the Lyne & Manchester model, pulsars with periods ∼ 100 ms beam to about 30% of the sky compared to the Narayan & Vivekanand model in which pulsars with periods below 100 ms beam to the entire sky. Note that, when many of these models were proposed, the sample of millisecond pulsars was < 5; hence their predictions about the beaming fractions of short-period pulsars relied largely on extrapolations from the normal pulsars. A recent analysis of a large sample of millisecond pulsar profiles [80] suggests that the beaming fraction of millisecond pulsars lies between 50% and unity.

3.3 The Population of Normal and Millisecond Pulsars

The most recent analysis to use the scaling function approach to derive the characteristics of the true normal and millisecond pulsar populations is based on the sample of pulsars within 1.5 kpc from the Sun [99]. The rationale for this cut-off is that, within this region, the selection effects are well understood and easier to quantify by comparison with the rest of the Galaxy. These calculations should give reliable estimates for the local pulsar population.

The luminosity distribution obtained from this analysis is shown in Fig. 10. These calculations lead to a local surface density of 30 ± 6 pulsars kpc^-2 for luminosities greater than 1 mJy kpc². Using Biggs’ [20] beaming model, the mean surface density of active pulsars with luminosities above 1 mJy kpc² is 156 ± 31 pulsars kpc^-2. Applying the same techniques to the sample of millisecond pulsars, and assuming a mean beaming fraction of 75% [80], the local surface density of millisecond pulsars with luminosities above 1 mJy kpc² is 37 ± 16 kpc^-2.

Figure 10

The corrected luminosity distribution (solid histogram with error bars) derived from a sample of 84 normal pulsars within 1.5 kpc of the Sun [99]. The observed distribution is shown by the dashed line. The difference between the observed and derived distributions below 10 mJy kpc² highlights the severe under-sampling of low-luminosity pulsars in the observed population.

These estimates of the local surface density of active pulsars allow us to deduce the likely distance to the nearest neutron star to Earth. This number is of interest to those building gravitational wave detectors, since it determines the likely amplitude of gravitational waves emitted from nearby rotating neutron stars. According to Thorne [151], currently planned detectors will be able to detect neutron stars with ellipticities greater than

$${ \in _{{\rm{min}}}}\; = 3 \times {10^{ - 10}} \times {\left( {\frac{P}{{200\;{\rm{ms}}}}} \right)^2} \times {\left( {\frac{d}{{1\;{\rm{kpc}}}}} \right)^2},$$

((6))

where P is the rotation period of a neutron star at a distance d from the Earth. For the combined millisecond and normal pulsar populations, with a surface density of 193 ± 35 pulsars kpc^-2, the nearest neutron star is thus likely to be < 40 pc. Future detections of such sources would be able to determine whether neutron stars have such ellipticities. One of the best known candidates is the nearby 5.75 ms pulsar J0437-4715 [74] which, at a distance of 178 ± 26 pc [130], is the closest known millisecond pulsar to the Earth.

Integrating the local surface densities of pulsars over the whole Galaxy requires a knowledge of the presently rather uncertain Galactocentric radial distribution [98, 71]. One approach is to assume that pulsars have a radial distribution similar to that of other stellar populations. The corresponding scale factor is then 1000±250 kpc² [123]. With this factor, we estimate there to be ∼ 160,000 active normal pulsars and ∼ 40, 000 millisecond pulsars in the Galaxy. Based on these estimates, we are in a position to deduce the corresponding rate of formation or birth-rate ${\mathcal B}$. From the P-̃P diagram in Fig. 5, we infer a typical lifetime for normal pulsars of ∼ 10⁷ yr, corresponding to a Galactic birth rate of ${\mathcal B}$ ∼ 1 per 60 yr — consistent with the rate of supernovae [157]. Different techniques yield consistent results [99]. As noted in §2.4, the ages of the millisecond pulsars are much older — close to that of the Universe ∼ 10¹⁰ yr. Taking the maximum age of the millisecond pulsars to be 10¹⁰ yr, we infer a mean birth rate of at least 4 × 10^-6 yr^-1. Given the uncertainties involved, this agrees satisfactorily with the birth-rate of low-mass X-ray binaries [89].

3.4 The Population of Neutron Star Binaries

Neutron star-neutron star (NS-NS) binaries similar to the classic Hulse-Taylor pulsar B1913+16 are of great interest, since they are expected to coalesce due to the emission of gravitational radiation. Only during the final few seconds of coalescence are gravitational waves strong enough to detect [133]. As noted in §2.4, NS-NS systems are expected to be rare since the binary system has to survive two supernova explosions. However, by including galaxies out to distances of ∼ 200 Mpc, it is hoped that merging NS-NS binaries will prove numerous enough to make a detection likely by the present and future generation of gravitational wave detectors [151, 134].

Of key importance here is the rate ${\mathcal R}$ of NS-NS mergers in the Galaxy which can then be extrapolated to more distant galaxies to provide event rate predictions for gravity wave detectors [117]. Statistical studies of the Galactic population of NS-NS systems are, however, hampered by small number statistics. The present sample summarised in Table 1 includes only three systems which have merging times significantly smaller than a Hubble time-scale viz: B1913+16 [68], B1534+12 [164], which are in the disk of our Galaxy, and B2127+11C which is located in the Globular Cluster M15 [120].

Table 1 Parameters of the known NS-NS binaries. Listed are the pulse period P, the orbital period P_orb, the pulsar characteristic age τ_c, the expected binary coalescence time-scale due to gravitational wave emission τ_g calculated according to [115].

Early estimates indicated a Galactic rate as high as ${\mathcal R} \sim 3 \times {10^{ - 4}}{\rm{y}}{{\rm{r}}^{ - 1}}$ [33]. These rather optimistic estimates had to be reduced by a factor of 30 [108, 117] following revised scale factor calculations based on the surveys that discovered PSRs B1534+12 and B2127+11C. These estimates implied that many new NS-NS systems should be found by the all-sky millisecond pulsar surveys [38]. This was surprisingly not the case; the only discovery being PSR J1518+4904 — a mildly relativistic binary system that will merge on a time-scale of ∼ 10¹² yr [111]. The most recent estimates of the NS-NS population, which take into account the lack of detections by the all-sky surveys, place a lower limit on the Galactic NS-NS population of ∼ 250 and a merging rate of 8 × 10^-6 yr^-1 [38, 158].

It is important to stress that these estimates are insensitive to systems fainter than some defined luminosity limit and therefore only provide a lower limit to ${\mathcal R}$. Bailes [8] proposed a means of estimating an upper bound on ${\mathcal R}$. Bailes postulates that the birth rates of normal pulsars in NS-NS systems must be equal to those of recycled pulsars in NS-NS systems. We know of only one NS-NS system in which the normal pulsar is visible: B2303+46 (see Table 1). This system will, however, not merge within a Hubble time-scale. In this case, ${\mathcal R}$ cannot be larger than ${\mathcal B}{\rm{/}}{N_{{\rm{obs}}}}$, where N_obs is the number of observed normal pulsars and ${\mathcal B}$ is their birth rate. Inserting the most recent numbers into this equation (N_obs ≃ 700; ${\mathcal B} \simeq {\rm{0}}{\rm{.015}}\;{\rm{y}}{{\rm{r}}^{{\rm{ - 1}}}}$) yields an upper limit of the formation rate of merging NS-NS binaries of ${\mathcal R} < 2 \times {10^{ - 5}}{\rm{y}}{{\rm{r}}^{ - 1}}$.

An independent method of estimating ${\mathcal R}$ can be made via population syntheses of binary stars [79, 44, 126, 156]. The essence of this approach is to follow the orbital and stellar evolution of a large number (> 10⁵) of binary star systems of varying mass and orbital separation. Based on a number of plausible physical arguments it is possible to predict the relative fractions of the various types of binary systems containing neutron stars. The most recent estimates using this method [82, 119] find ${\mathcal R} \sim 2 - 5 \times {10^{ - 5}}{\rm{y}}{{\rm{r}}^{ - 1}}$. Given the uncertainties involved, we conclude that both methods yield consistent results i.e.:

$${\mathcal R} \sim {10^{ - 5.0 \pm 0.5}}{\rm{y}}{{\rm{r}}^{ - 1}}.$$

Extrapolations including Galaxies out to ∼ 200 Mpc suggest that the rate of NS-NS mergers will be high enough to yield detection rates of ∼ several sources per year. Ultimately, the detection statistics from the gravity wave detectors should provide far tighter constraints on the NS-NS merging rate.

3.5 Going further

Studies of population statistics represents a large proportion of the yearly publications in pulsar astronomy. Good starting points for further reading can be found in [145, 95, 105]. A demonstration version of the widely-cited “Scenario Machine” code, developed by Lipunov and collaborators to perform population syntheses of binary stars, is available on-line [132].

4 Pulsar Timing

It became clear soon after their discovery that pulsars are excellent celestial clocks. In the original discovery paper [65], the period of the first pulsar to be discovered, PSR B1919+21, was found to be stable to one part in 10⁷ over a time-scale of a few months. Following the discovery of the millisecond pulsar B1937+21 in 1982 [6] it was demonstrated that its period could be measured to one part in 10¹³ or better [41]. This unrivaled stability leads to a host of applications including time keeping, probes of relativistic gravity and natural gravitational wave detectors. Subsequently a whole science has developed into measuring the pulse time-of-arrival as accurately as possible.

4.1 Observing basics

As each new pulsar is discovered, the standard practice is to add it to a list of pulsars which are regularly observed at least once or twice per month by large radio telescopes throughout the world. The schematic diagram shown in Fig. 11 summarises the essential steps involved in such a “time-of-arrival” (TOA) measurement. Incoming pulses emitted by the rotating neutron star traverse the interstellar medium before being received by the radio telescope. After amplification by high sensitivity receivers, the pulses are de-dispersed (see below) and added to form a mean pulse profile.

Figure 11

Schematic showing the main stages involved in pulsar timing observations.

During the observation, the data regularly receive a time stamp, usually based on a caesium time standard or hydrogen maser at the observatory plus a signal from the GPS (Global Positioning System of satellites) time system [40]. The TOA of this mean pulse is then defined as the arrival time of some fiducial point on the profile. Since the mean profile has a stable form at any given observing frequency (§2.2) the TOA can be accurately determined by a simple cross-correlation of the observed profile with a high signal-to-noise “template” profile — obtained from the addition of many observations of the pulse profile at the particular observing frequency.

Success in pulsar timing hinges on how precisely this fiducial point can be determined. This is largely dependent on the signal-to-noise ratio (SNR) of the mean pulse profile. The uncertainty in a TOA measurement ∈_TOA is given roughly by the pulse width divided by the SNR. Based on equation 3, we can express this as follows

$${ \in _{{\rm{TOA}}}} = \frac{W}{{{\rm{SNR}}}} \simeq \frac{{{S_{{\rm{sys}}}}}}{{{S_{{\rm{psr}}}}}}\frac{P}{{\sqrt {\Delta \nu {t_{{\rm{int}}}}} }}{\left( {\frac{W}{P}} \right)^{3/2}},$$

((7))

where W is the equivalent width of a top hat function having the same mean flux density S_psr as the observed pulse with a period P. The observing system is parameterised by the available bandwidth Δν, the integration time t_int, and the system equivalent flux density S_sys. Optimum results are thus obtained for observations of short period pulsars with large flux densities and narrow duty cycles (W/P) using large telescopes with low-noise receiver systems (both effects reduce S_sys) and large observing bandwidths.

One of the main problems of employing large bandwidths is pulse dispersion. As discussed in §2.3, the velocity of the pulsed radiation through the ionised interstellar medium is frequency-dependent: Pulses emitted at higher radio frequencies travel faster and arrive earlier than those emitted at lower frequencies. This process has the effect of “stretching” the pulse across a finite receiver bandwidth, reducing the apparent signal-to-noise ratio and therefore increasing ∈_TOA. For most normal pulsars, this process can largely be compensated for by the incoherent de-dispersion process outlined in §3.1.

To exploit the precision offered by millisecond pulsars, a more precise method of dispersion removal is required. Technical difficulties in building devices with very narrow channel bandwidths require another dispersion removal technique. In the process of coherent de-dispersion [58] the incoming signals are de-dispersed over the whole bandwidth using a filter which has the inverse transfer function to that of the interstellar medium. The maximum resolution obtainable is then the inverse of the receiver bandwidth. For typical bandwidths of 10 MHz, this technique makes it possible to resolve features on time-scales as short as 100 ns. This corresponds to regions in the neutron star magnetosphere as small as 30 m!

4.2 The Timing Model

Ideally, in order to model the rotational behaviour of the neutron star, we require TOAs measured by an inertial observer. Due to the Earth’s orbit around the Sun, an observatory located on Earth experiences accelerations with respect to the neutron star. The observatory is therefore not in an inertial frame. To a very good approximation, the centre-of-mass of the solar system, the solar system barycentre, can be regarded as an inertial frame. It is standard practice [69] to transform the observed TOAs to this frame using a planetary ephemeris such as the JPL DE200 [139]. The transformation is summarised in the following equation as the difference between barycentric (${\mathcal T}$) and observed (t) TOAs:

$$\mathcal{T} - t = \frac{{\underline r .\hat {\underline {s}}}}{c} + \frac{{{{(\underline r.\hat {\underline {s}})}^2} - |\underline r{|^2}}}{{2cd}} + \Delta {t_{{\text{rel}}}} - \Delta {t_{{\text{DM}}}}.$$

((8))

Here ̱r is the position of the Earth with respect to the barycentre, $\underline {\hat s}$ is a unit vector in the direction towards the pulsar at a distance d, and c is the speed of light. The first term on the right hand side of this expression is the light travel time from the Earth to the solar system barycentre. For all but the nearest pulsars, the incoming pulses can be approximated by plane wavefronts. The second term, which represents the delay due to spherical wavefronts and which yields the trigonometric parallax and hence d, is only presently measurable for four nearby millisecond pulsars [77, 29, 130]. The term Δt_rel represents the Einstein and Shapiro corrections due to general relativistic effects within the solar system [5]. Since measurements are often carried out at different observing frequencies, each with its own dispersive delays, the TOAs are generally shifted back to the equivalent time that would be observed at infinite frequency. This shift corresponds to the term Δ t_DM and may be calculated from equation 1.

Following the accumulation of about ten to twenty barycentric TOAs from observations spaced over at least several months, a surprisingly simple model can be applied to the TOAs and optimised so that it is sufficient to account for the arrival time of any pulse emitted during the time span of the observations and predict the arrival times of subsequent pulses. The model is based on a Taylor expansion of the angular rotational frequency Ω = 2π/P about a model value Ω_o at some reference epoch t_o. The model pulse phase ϕ as a function of barycentric time is thus given by:

$$\phi ({\mathcal T}) = {\phi _{\rm{o}}} + ({\mathcal T} - {{\mathcal T}_{\rm{o}}}){\Omega _{\rm{o}}} + \frac{1}{2}{({\mathcal T} - {{\mathcal T}_{\rm{o}}})^2}{\dot \Omega _{\rm{o}}} + \cdots ,$$

((9))

where ϕ_o is the pulse phase at ${{\mathcal T}_{\rm{o}}}$. Based on this simple model, and using initial estimates of the position, dispersion measure and pulse period, a “timing residual” is calculated for each TOA as the difference between the observed and predicted pulse phases.

Early sets of residuals will exhibit a number of trends indicating a systematic error in one or more of the model parameters, or a parameter not initially incorporated into the model (e.g. ̇Ω_o). From equation 9, an error in the assumed Ω_o results in a linear slope with time. A parabolic trend results from an error ̇Ω_o. Additional effects will arise if the assumed position of the pulsar (the unit vector $\underline {\hat s}$ in equation 8) is incorrect. A position error of just one arcsecond results in an annual sinusoid with a peak-to-peak amplitude of about 5 ms for a pulsar on the ecliptic; this is easily measurable for typical TOA uncertainties of order one milliperiod or better. Similarly, the effect of a proper motion produces an annual sinusoid of linearly increasing magnitude.

After a number of iterations, and with the benefit of a modicum of experience, it is possible to identify and account for each of these various effects to produce a “timing solution” which is phase coherent over the whole data span. The resulting model parameters provide spin and astrometric information about the neutron star to a precision which improves as the length of the data span increases. The latest observations of the original millisecond pulsar, B1937+21, spanning almost 9 yr (exactly 165,711,423,279 rotations!) measure a period of 1.5578064688197945 ± 0.0000000000000004 ms [77, 75] defined at midnight UT on December 5 1988! Astrometric measurements based on these data are no less impressive, with position errors of ∼ 20 μarcsec being presently possible.

4.3 Binary pulsars

For binary pulsars, the timing model described in §4.2 needs to be extended to incorporate the additional radial accelerations of the pulsar as it orbits the common centre-of-mass of the binary system. Treating the binary orbit using Kepler’s laws to refer the TOAs to the binary barycentre requires five additional model parameters: the orbital period, projected semi-major orbital axis, orbital eccentricity, longitude of periastron and the epoch of periastron passage. Light travel times over binary pulsar systems are typically of the order of several seconds. Since this corresponds to a large number of pulse periods, obtaining a phase coherent timing solution thus requires accurate determination of the orbital parameters from a densely sampled initial set of observations.

The Keplerian formularism is analogous to spectroscopic binary stars where constraints on the mass of the orbiting companion can be placed by combining the projected semi-major axis a_p sin i and the orbital period P_o to obtain the familiar mass function:

$$f({m_{\rm{p}}},\;{m_{\rm{c}}}) = \frac{{4{\pi ^2}}}{{\rm{G}}}\frac{{{{({a_{\rm{p}}}\;{\rm{sin}}\;i)}^3}}}{{P_{\rm{o}}^2}} = \frac{{{{({m_{\rm{c}}}\;{\rm{sin}}\;i)}^3}}}{{{{({m_{\rm{p}}} + {m_{\rm{c}}})}^2}}},$$

((10))

where G is the universal gravitational constant. Assuming a pulsar mass m_p of 1.4M_⊙, the mass of the orbiting companion m_c can be estimated as a function of the (unknown) angle i between the orbital plane and the plane of the sky. The minimum companion mass m_min occurs when the orbit is assumed edge-on (i = 90°). For a random distribution of orbital inclination angles, the probability of observing a binary system at an angle less than some value i_o is p(< i_o) = 1 - cos(i_o). This implies that the chances of observing a binary system inclined at an angle < 26° is only 10%; evaluating the companion mass for this inclination angle m₉₀ constrains the mass range between m_min and m₉₀ (in a probabilistic sense) at the 90% confidence level. The mass distribution of the companions in the present sample discussed in §2.4 indicates that the companions are white dwarfs, other neutron stars or main sequence stars.

Figure 12

The variation in apparent pulse period for the binary pulsar J1012+5307 as it moves about the centre of mass every 14.5 hours. Points denote observations of the pulse period, the smooth curve is a Keplerian fit to the data.

An example of a purely Keplerian orbital solution is shown in Fig. 12 by the fit to a set of pulse period measurements for the binary pulsar J1012+5307. The peak-to-peak variations in period correspond to a change in radial velocity during the 14.5 hour binary orbit of about 46 km s^-1. The eccentricity of the orbit dictates the shape of the orbital curve. In this case e < 2 × 10^-5 and the Doppler shifts are indistinguishable from a sinusoid. For the higher-mass binary systems with e > 0.3 the curve becomes increasingly like a saw-tooth in appearance.

Although many of the presently known binary pulsar systems can be adequately described by a Keplerian orbit, there are a number of systems, including the original binary pulsar B1913+16, which exhibit a number of relativistic effects requiring an additional set of “post-Keplerian” parameters [23, 62, 39, 147]. Taylor & Weisberg [146] measured two such parameters from an analysis of timing data for B1913+16 between 1974 and 1981. The first of these parameters is advance of the longitude of periastron, analogous to the general relativistic perihelion advance of Mercury. For PSR B1913+16 this amounts to about 4.2 degrees per year [146], some 4.5 orders of magnitude larger than for Mercury. The second post-Keplerian parameter is related to the gravitational redshift and transverse Doppler shifts in the orbit. Together with the five Keplerian parameters, these two post-Keplerian parameters allow an unambiguous description of the masses of the binary components, together with the inclination angle of the orbit.

The latest published measurements for PSR B1913+16 [147] determine the mass of the pulsar and its (unseen) companion to be 1.442 ± 0.003M_⊙ and 1.386 ± 0.003M_⊙ respectively, demonstrating that the companion is most likely another neutron star (§2.4). Both these post-Keplerian parameters have now been measured for the two other double neutron star binary systems that will merge within a Hubble time: B1534+12 [164] and B2127+11C [42]. These measurements provide similar, but slightly less precise, determinations of neutron star masses. For a number of other systems, the advance of the longitude of periastron has been measured yielding the total system mass and thus allowing interesting constraints to be placed on the component masses [152, 111].

An important general relativistic prediction for eccentric double neutron star systems is the orbital decay due to the emission of gravitational radiation. Taylor & Weisberg were able to measure this for B1913+16 and found it to be in excellent agreement with the predicted value. The orbital decay, which corresponds to a shrinkage of about 1.5 cm per orbit, is seen most dramatically as the gradually increasing shift in orbital phase for periastron passages with respect to a non-decaying orbit. This is reproduced from Taylor & Weisberg’s classic paper [146] in Fig. 13. Further observations of this system demonstrate that the agreement with General Relativity is at the level of about 1% [147]. Hulse and Taylor were awarded the Nobel prize in Physics in 1993 [150] in recognition of the discovery of this remarkable laboratory for testing General Relativity.

Figure 13

Orbital decay in the binary pulsar B1913+16 system demonstrated as an increasing orbital phase shift for periastron passages with time. The general relativistic prediction due entirely to the emission of gravitational radiation is shown by the parabola. (After Taylor & Weisberg 1982).

For those binary systems which are oriented nearly edge-on to the line-of-sight, a significant delay is expected for orbital phases around superior conjunction where the pulsar radiation is bent in the gravitational potential well of the companion star. This effect, analogous to the solar system Shapiro delay, has so far been measured for two neutron star-white dwarf binary systems: B1855+09 and J1713+0747 [128, 77, 29].

4.4 Timing Stability

Ideally, after correctly applying a timing model, we would expect a set of uncorrelated timing residuals scattered in a Gaussian fashion about a zero mean with an rms consistent with the measurement uncertainties. This is not always the case; the residuals of many pulsars exhibit a quasi-periodic wandering with time . A number of examples are shown in Fig. 14. These are taken from the Jodrell Bank timing program [136, 104]. Such “timing noise” is most prominent in the youngest of the normal pulsars [101, 37] and virtually absent in the much older millisecond pulsars [77]. Whilst the physical processes of this phenomena are not well understood, it seems likely that this phenomenon may be connected to superfluid processes in the interior of the neutron star and its temperature [2] or processes in the magnetosphere [32, 31].

Figure 14

Examples of timing residuals for a number of normal pulsars. Note the varying scale on the ordinate axis, the pulsars being ranked in increasing order of timing “activity” (Figure provided by Andrew Lyne).

The relative dearth of timing noise for the older pulsars is a very important finding. It implies that, presently, the measurement precision depends primarily on the particular hardware constraints of the observing system. Consequently, a large effort in hardware development is presently being made to improve the precision of these observations using, in particular, coherent de-dispersion outlined in §4.1. Much of the pioneering work in this area has been made by Joseph Taylor and collaborators at Princeton University [121]. From high quality observations made using the Arecibo radio telescope spanning almost a decade [128, 129, 77], the group has demonstrated that the timing stability of millisecond pulsars over such time-scales is comparable to terrestrial atomic clocks.

This phenomenal stability is demonstrated in Fig. 15. This figure shows σ_z, a parameter closely resembling the Allan variance used by the clock community to estimate the stability of atomic clocks [143]. Atomic clocks are known to have σ_z ∼ 5 × 10^-15 on time-scales of order 5 years. The timing stability of PSR B1937+21 seems to be limited by a power law component which produces a minimum in its σ_z after ∼ 2 yr. This is most likely a result of a small amount of intrinsic timing noise [77]. No such noise component is presently seen in the Allan variance of PSR B1855+09. This demonstrates that the timing stability for PSR B1855+09 becomes competitive with the atomic clocks after about 3 yr. The absence of timing noise for B1855+09 is probably related to its characteristic age ∼ 5 Gyr which is about a factor of 20 larger than B1937+21.

Figure 15

Fractional instabilities for PSRs B1855+09 and B1937+21 expressed as the modified Allan variance σ_z as a function of time. (After Kaspi, Taylor & Ryba 1994 [77]).

Recently, a number of millisecond pulsars have been discovered in the all-sky surveys discussed in §3.1. These include PSRs J0437-4715 [74] and J1713+0747 [53] — two bright millisecond pulsars with small duty cycles which allow very high precision measurements (equation 7). Early indications are that their timing stability will be at least as good as B1855+09 [52, 130].

4.5 Going further

For further details on the technical details and prospects of pulsar timing, the interested reader is referred to a number of excellent review articles [5, 143, 93, 11, 12]. An up-to-date list summarising the various timing programmes is kept by Don Backer [16]. Two freely available software packages which are routinely used for such analyses amongst the pulsar community are also available viz: TEMPO [147, 121] and TIMAPR [47, 46]. These packages are based on more detailed versions of the timing model outlined in §4.2.

The remarkable precision of these measurements, particularly for millisecond pulsars, allows the detection of radial accelerations on the pulsar induced by orbiting bodies smaller than the Earth. Alex Wolszczan detected one such “pulsar planetary system” in 1990 following the discovery of a 6.2 ms pulsar B1257+12. In this case, the pulsar is orbited by at least two ∼ Earth-mass bodies [166, 122]. Subsequent measurements of B1257+12 were even able to measure resonance interactions between two of the planets [165], confirming the nature of the system beyond all doubt. Long-term timing measurements of the 11 ms pulsar B1620-26 in the globular cluster M4 indicate that it may also have a planetary companion [153, 4]. For detailed reviews of these systems, and their implications for planetary formation scenarios, the interested reader is referred to [116].

5 Pulsars as Gravity Wave Detectors

Many cosmological models predict that the Universe is presently filled with a stochastic gravitational wave background (GWB) produced during the big bang era [114]. The idea to use pulsars as natural detectors of gravitational waves was first explored independently by Sazhin and Detweiler in the late 1970s [131, 43]. The basic concept is to treat the solar system barycentre and a distant pulsar as opposite ends of an imaginary arm in space. The pulsar acts as the reference clock at one end of the arm sending out regular signals which are monitored by an observer on the Earth over some time-scale T. The effect of a passing gravitational wave would be to cause a change in the observed rotational frequency by an amount proportional to the amplitude of the wave. For regular monitoring observations of a pulsar with typical TOA uncertainties of ∈_TOA, this “detector” would be sensitive to waves with dimensionless amplitudes > ∈_TOA/T and frequencies as low as ∼ 1/T [17, 24]. This method, which already yields interesting upper limits on the GWB, is reviewed in §5.1. The idea of more sensitive detector based on an array of pulsar clocks distributed over the sky is discussed in §5.2.

5.1 Limits from individual pulsars

In the ideal case, the change in the observed frequency caused by the GWB should be detectable in the set of timing residuals after the application of an appropriate model for the rotational, astrometric and, where necessary, binary parameters of the pulsar. As discussed in §4, all other effects being negligible, the rms scatter of these residuals σ would be due to the measurement uncertainties and intrinsic timing noise from the neutron star. Detweiler [43] showed that a GWB with a flat energy spectrum in the frequency band f±f/2 would result in an additional contribution to the timing residuals σ_g. The corresponding wave energy density ͡_g (for fT ≫ 1) is

$${\rho _g} = \frac{{243{\pi ^3}{f^4}{\sigma ^2}g}}{{208\;G}}.$$

((11))

An upper limit to ρ_g can be obtained from a set of timing residuals by assuming the rms scatter is entirely due to this effect (σ = σ_g). These limits are commonly expressed as a fraction of ρ_c the energy density required to close the Universe:

$${\rho _c} = \frac{{3\;H_{\rm{o}}^2}}{{8\pi \;G}} \simeq 2 \times {10^{ - 29}}{h^2}\;g\;{\rm{c}}{{\rm{m}}^{ - 3}},$$

((12))

where the Hubble constant H_o = 100 h km s^-1 Mpc.

Romani & Taylor [127] applied this technique to a set of TOAs for PSR B1237+12 obtained from regular observations over a period 11 years as part of the JPL pulsar timing programme [48]. This pulsar was chosen on the basis of its relatively low level of timing activity by comparison with the youngest pulsars, whose residuals are ultimately plagued by timing noise ( §4.4). By ascribing the rms scatter in the residuals (σ = 240 ms) to the GWB, Romani & Taylor placed a limit of ρ_g/ρ_c < 4 × 10^-3h^-2 for a centre frequency f = 7 × 10^-9 Hz.

This limit, already well below the energy density required to close the Universe, was further reduced following the long-term timing measurements of millisecond pulsars at Arecibo by Taylor and collaborators (§4.4). In the intervening period, more elaborate techniques had been devised [17, 24, 140] to look for the likely signature of a GWB in the frequency spectrum of the timing residuals and to address the possibility of “fitting out” the signal in the TOAs. Following [17], it is convenient to define Ω_g, the energy density of the GWB per logarithmic frequency interval relative to ρ_c. With this definition, the power spectrum of the GWB ${\cal P}$(f) can be written [67, 24] as

$${\mathcal P}(f) = \frac{{G\;{\rho _g}}}{{3{\pi ^3}\;{f^4}}} = \frac{{H_{\rm{o}}^2{\Omega _g}}}{{8{\pi ^4}\;{f^5}}} = 1.34\; \times \;{10^4}{\Omega _g}{h^2}\;f_{y{r^{ - 1}}}^{ - 5}\;\mu {{\rm{s}}^2}{\rm{yr,}}$$

((13))

where f_yr-1 is frequency in cycles per year. The timing residuals for B1937+21 shown in Fig. 16 are clearly non-white and, as we saw in §4.4, limit its timing stability for periods > 2 yr. The residuals for PSR B1855+09 clearly show no systematic trends and are in fact consistent with the measurement uncertainties alone. Based on these data, and using a rigorous statistical analysis, Thorsett & Dewey [154] place a 95% confidence upper limit of Ω_gh² < 10^-8 for f = 4.4 × 10^-9 Hz. This limit is difficult to reconcile with most cosmic string models for galaxy formation [25, 154].

Figure 16

Timing residuals for PSRs B1855+09 (top panel) B1937+21 (bottom panel) obtained from almost a decade of timing at Arecibo. (After Kaspi, Taylor & Ryba 1994 [77]).

For those pulsars in binary systems, an additional clock for measuring the effects of gravitational waves is the orbital period. In this case, the range of frequencies is not limited by the time span of the observations, allowing the detection of waves with periods as large as the light travel time to the binary system [17]. The most stringent results presently available are based on B1855+09 limit Ω_gh² < 2.7 × 10^-4 in the frequency range 10^-11 < f < 4.4 × 10^-9 Hz. Kopeikin [78] has recently presented this limit and discusses the methods in detail.

5.2 A pulsar timing array

The idea of using timing data for a number of pulsars distributed on the sky to detect gravitational waves was first proposed by Hellings & Downs [64]. Such a “timing array” of pulsars would have the advantage over a single arm in that, through a cross-correlation analysis of the residuals for pairs of pulsars distributed over the sky, it should be possible to separate the timing noise of each pulsar from the signature of the GWB, which would be common to all pulsars in the array. To quantify this, consider the fractional frequency shift observed for the i^th pulsar in the array:

$$\frac{{\delta {\nu _i}}}{{{\nu _i}}} = {\alpha _i}{\mathcal A}(t) + {{\mathcal N}_i}(t).$$

((14))

In this expression α_i is a geometric factor dependent on the line-of-sight direction to the pulsar and the propagation and polarisation vectors of the gravity wave of dimensionless amplitude ${\mathcal A}$. The timing noise intrinsic to the pulsar is characterised by the function ${{\mathcal N}_i}$. The result of a cross-correlation between pulsars i and j is then

$$\alpha {}_i{\alpha _j} < {{\mathcal A}^2} > + {\alpha _i} < {\mathcal A}{{\mathcal N}_j} > + {\alpha _j} < {\mathcal A}{{\mathcal N}_i} > + < {{\mathcal N}_i}{{\mathcal N}_j} > ,$$

((15))

where the bracketed terms indicate cross-correlations. Since the wave function and the noise contributions from the two pulsars are independent quantities, the cross correlation tends to ${\alpha _i}{\alpha _j} < {{\mathcal A}^2} >$ as the number of residuals becomes large. Summing the cross-correlation functions over a large number of pulsar pairs provides additional information on this term as a function of the angle on the sky [63]. This allows, in principle, the separation of the effects of terrestrial clock and solar system ephemeris errors from the GWB [51].

Applying the timing array concept to the present database of long-term timing observations of millisecond pulsars does not improve on the limits on the GWB discussed above. The sky distribution of these pulsars, seen in the left panel of Fig. 17, shows that their angular separation is rather low. To achieve optimum sensitivity it is desirable to have an array consisting of pulsar clocks distributed isotropically over the whole sky. The flood of recent discoveries of nearby binary and millisecond pulsars by the all-sky searches has resulted in essentially such a distribution, shown in the right panel of Fig. 17. Continued timing of these pulsars in the coming years should greatly improve the sensitivity and will perhaps allow the detection of gravity waves, as opposed to upper limits, in the future.

Figure 17

Hammer-Aitoff projections showing the known Galactic disk millisecond pulsar population in 1990 (left) and in 1997 (right). The impact of the new discoveries is clearly seen comparing the situation circa 1990, in which all of the known sources had been found at Arecibo, with the present situation in which the sources are much more uniformly distributed on the sky.

5.3 Going further

Further discussions on the realities of using pulsars as gravity wave detectors can be found in two excellent review articles by Romani [125] and Backer [3].

6 Summary and Future Prospects

The main aim of this article was to review, from an observational perspective, some of the many astrophysical applications provided by the present sample of binary and millisecond radio pulsars. In writing this article I hope that many of the broad-ranging results in this field were covered and that it will be of use as a starting point for those wishing to explore the vast body of literature that presently exists.

Through an understanding of the Galactic population of radio pulsars summarised in §3, it is possible to predict the detection statistics of terrestrial gravitational wave detectors to nearby rapidly spinning neutron stars (§3.3), as well as coalescing binary neutron stars at cosmic distances (§3.4). I look forward to the continued improvements in gravitational wave detector sensitivities. These should result in a number of interesting developments and contributions in this area. These developments and contributions might include the detection of presently known radio pulsars, as well as a population of coalescing binary systems which have not yet been detected as radio pulsars.

The phenomenal timing stability of radio pulsars leads naturally to a large number of applications, including their use as laboratories for relativistic gravity (§4.3) and as natural detectors of gravitational radiation (§5). Long-term timing experiments of the present sample of millisecond and binary pulsars currently underway appear to have tremendous potential to dig more deeply in these areas and perhaps detect the gravitational wave background (if it exists) in the next 5–10 years.

These and other applications will benefit greatly from the continued discovery of new systems by the present generation of radio pulsar searches. These searches continue to probe new areas of parameter space. Based on the results presented in §3.3, it is clear that we are aware of only about 1% of the total active pulsar population in our Galaxy. Thus, in terms of the exotic systems that exist, we have almost certainly not seen all of the pulsar zoo. I look forward to the discovery of a number of unexpected systems in the coming years….perhaps a dual-line binary pulsar, a neutron star with a black hole companion, as well as the exciting possibility of sub-millisecond pulsars.