Preamble

In the 34 years that have elapsed since the discovery [97] 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 include the confirmation of the existence of gravitational radiation [241] as predicted by general relativity [239, 240] and the first detection of an extra-solar planetary system [273, 186].

The diverse zoo of radio pulsars currently known is summarized graphically by the Venn diagram in Fig. 1. 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 the results and spin-offs from these surveys which are of particular interest to the relativity community.

Figure 1
figure 1

The numbers and locations of the various types of radio pulsars known as of December 2000. The large and small Magellanic clouds are denoted by LMC and SMC.

What’s new in this review?

Since the first version of this article was written back in 1997/8 [128] a number of pulsar surveys using the Parkes radio telescope [157] have discovered almost 700 pulsars. As a result, the sample size is now double what it was in 1997. Many of the exciting new discoveries from these searches are discussed in this review. Up-to-date tables of parameters of binary and millisecond pulsars are included as an appendix. Several new sections/figures have been added and existing sections reworked and modularized to make the review more self-contained and (hopefully!) easier to read in an html setting. We begin in § 2 with an overview of the pulsar phenomenon, the key observed population properties, the origin and evolution of pulsars and an introduction to pulsar search techniques. In § 3, we review present understanding in pulsar demography, discussing selection effects and the techniques used to correct for them in the observed sample. This leads to robust estimates of the total number of normal and millisecond pulsars (§ 3.3) and relativistic binaries (§ 3.4) in the Galaxy and has implications for the detection of gravitational radiation from these systems. We discuss pulsar timing in § 4. One application of these exceptional clocks, a sensitive detector of long-period gravitational waves, is discussed in § 5. We conclude with a brief outlook to the future in § 6.

An Introduction to Pulsar Astronomy

Many of the basic observational facts about radio pulsars were established shortly after their discovery [97] by Bell and 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–10 M) stars. In the following, we discuss the basic observational properties most relevant to this review.

The lighthouse model

An animation showing the rotating neutron star or “lighthouse model” of the basic pulsar phenomenon is shown in Fig. 2.

Figure 2
figure 2

The rotating neutron star (or “lighthouse”) model for pulsar emission. To see the movie in action go to the online version of this article on http://www.livingreviews.org/Articles/Volume4/2001-5lorimer. Animation designed by Michael Kramer.

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 repetition period of the pulses is therefore simply the rotation period of the neutron star. The moving “tracker ball” on the pulse profile in the animation shows the relationship between observed intensity and rotational phase of the neutron star.

Neutron stars are extremely stable rotators. They are essentially large celestial flywheels with moments of inertia ∼ 1045 g cm2. The rotating neutron star model, independently developed by Pacini and Gold in 1968 [184, 85], predicts a gradual increase in the pulse period as the outgoing radiation carries away rotational kinetic energy. This model became universally accepted when a period increase of 36.5 ns per day was measured for the pulsar in the Crab nebula [202], enabling Gold [86] to show that a rotating neutron star with a large magnetic field must be the dominant energy supply for the nebula.

Pulse profiles

Pulsars are weak radio sources. Measured flux densities, usually quoted in the literature for a radio frequency of 400 MHz, vary between 0.1 and 5000 mJy (1 Jy≡10-26 W m-2 Hz-1). This means that, even with a large radio telescope, the coherent addition of many thousands of pulses is required in order to produce an integrated profile. Remarkably, although the individual pulses vary 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 fingerprint of the emission beam.

The animation in Fig. 3 shows a sequence of consecutive single pulses from PSR B0329+54Footnote 1, one of the brightest pulsars. This pulsar is seen in the animation to stabilise into its characteristic 3-component form after the summation of a number of seemingly erratic single pulses. Stabilization time-scales are typically several hundred pulses [94]. This property is of key importance in pulsar timing measurements discussed in detail in § 4.

Figure 3
figure 3

Single pulses from PSR B0329+54. To see the movie in action go to the online version of this article on http://www.livingreviews.org/Articles/Volume4/2001-5lorimer.

Fig. 4 shows the rich diversity in morphology from simple single-component profiles to examples in which emission is observed over the entire pulse. The astute reader will notice two examples of “interpulses” — a secondary pulse separated by about 180 degrees from the main pulse. The most natural interpretation for this phenomenon is that the two pulses originate from opposite magnetic poles of the neutron star (see however [158]). Since this is an unlikely viewing angle we would expect interpulses to be a rare phenomenon. Indeed this is the case: the fraction of known pulsars in which interpulses are observed in their pulse profiles is only a few percent.

Figure 4
figure 4

A variety of integrated pulse profiles taken from the available literature. References: (a, b, d, f: [ 87]); (c: [ 20]); (e, g, i: [ 122]); (h: [ 26]). Each profile represents 360 degrees of rotational phase. These profiles are part of a database of over 2600 multi-frequency pulse profiles for over 600 pulsars that is available on-line [ 164].

Two contrasting phenomenological models to explain the observed pulse shapes are shown in Fig. 5. The “core and cone” model, proposed by Rankin [198], depicts the beam as a core surrounded by a series of nested cones. Alternatively, the “patchy beam” model, championed by Lyne and Manchester [149, 89], has the beam populated by a series of randomly-distributed emitting regions. Further work in this area, particularly in trying to quantify the variety of pulse shapes (number of distinct components and the relative fraction that they occur) is necessary to improve our understanding of the fraction of sky covered by the radio pulsar emission beam. We return to this topic in the context of pulsar demography later on in § 3.2.

Figure 5
figure 5

Phenomenological models for pulse shape morphology produced by different line-of-sight cuts of the beam (Figure designed by M. Kramer and A. von Hoensbroech).

The pulsar distance scale

From the sky distribution shown in Fig. 6 it is immediately apparent that pulsars are strongly concentrated along the Galactic plane. This indicates that pulsars populate the disk of our Galaxy. Unlike most other classes of astrophysical objects, quantitative estimates of the distances to each pulsar can be made from an effect known as pulse dispersion, the delay in pulse arrival times across a finite bandwidth. Dispersion occurs because the group velocity of the pulsed radiation through the ionised component of the 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 and a low frequency νlo pulse is given [156] by

$$\Delta t = 4150\;{\rm{s}} \times (\nu _{{\rm{lo}}}^{ - 2} - \nu _{{\rm{hi}}}^{ - 2}) \cdot {\rm{DM}},$$
((1))

where the frequencies are in MHz and the dispersion measure DM (cm-3 pc) is the integrated column density of free electrons along the line of sight:

$$ {\rm{DM}} = \int_0^d {{n_{\rm{e}}}dl.} $$
((2))

Here, d is the distance to the pulsar (pc) and ne is the free electron density (cm-3). From Equation (2) it is obvious that a measurement of the delay across a finite bandwidth yields the DM. 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. Hence, given the DM, the distance can be estimated from a model of the Galactic distribution of free electrons.

Figure 6
figure 6

The sky distribution of 1026 pulsars in Galactic coordinates. The plane of the Galaxy is the central horizontal line. The Galactic centre is the midpoint of this line.

The electron density model is calibrated from the pulsars with independent distance estimates and measurements of scattering for lines of sight towards various Galactic and extragalactic sources. Independent distance estimates now exist for over 100 pulsars based on 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). Based on these data, Taylor & Cordes [235] have developed an electron density model which is free from large systematic trends and can be used to provide distance estimates with an uncertainty of ∼30%. However, use of this model to estimate distances to individual pulsars may result in uncertainties by as much as a factor of two. This model is currently being refined following recent pulsar discoveries and independent distance and scattering measurements (J. Cordes, private communication).

Normal and millisecond pulsars

Spin parameters

The present public-domain catalogue, available on-line at Princeton [196], contains up-to-date parameters for 706 pulsars. Parameters for many of the new Parkes multibeam pulsar surveys are also available on-line [28, 197]. Most of these are “normal” in the sense that their pulse periods P are of order one second and 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 rates of slowdown ≲10-19 s/s. The first millisecond pulsar discovered, B1937+21 [17], with P=1.5578 ms, remains the most rapidly rotating neutron star presently known.

A very useful means of demonstrating the distinction between these two classes is the “P- diagram” — a logarithmic scatter plot of the observed pulse period versus the period derivative. As shown in Fig. 7, normal pulsars occupy the majority of the upper right hand part of the diagram, while the millisecond pulsars reside in the lower left hand part of the diagram. The differences in P and 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 (PṖ)1/2 [163]. Lines of constant magnetic field strength are drawn on Fig. 7, together with lines of constant characteristic agec=P/(2Ṗ)). Typical inferred magnetic fields and ages are 1012 G and 107 yr for the normal pulsars and 108 G and 109 yr for the millisecond pulsars.

Figure 7
figure 7

The ubiquitous P-Ṗ diagram shown for a sample of radio pulsars. Those objects known to be members of binary systems are highlighted by a circle (for low-eccentricity orbits) or an ellipse (for elliptical orbits). Pulsars thought to be associated with supernova remnants are highlighted by the starred symbols.

Binary companions

As can be inferred from Fig. 1, just under 4% of all known pulsars in the Galactic disk 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. Two notable hybrid systems are the “planet pulsars” PSR B1257+12 and B1620–26. B1257+12 is a 6.2-ms pulsar accompanied by at least three Earth-mass bodies [273, 186, 272] while B1620—26, an 11-ms pulsar in the globular cluster M4, is part of a triple system with a white dwarf and a high-mass planet [247, 15, 245].

A very important additional difference between normal and millisecond pulsars is the presence of an orbiting companion. Orbital companions are much more commonly observed around millisecond pulsars (∼ 80% of the observed sample) than around the normal pulsars (≲ 1%). Fig. 8 is a scatter plot of orbital eccentricity versus mass of the companion. The dashed line serves merely to guide the eye in this figure. Binary systems lying below the line are those with low-mass companions (≲0.7M — predominantly white dwarfs) and essentially circular orbits: 10-5≲e≲0.1. Binary pulsars with high-mass companions (≳1M — neutron stars or main sequence stars) are in eccentric orbits, 0.15≳e≳0.9, and lie above the line.

Figure 8
figure 8

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

Evolutionary scenarios

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 [34, 77, 218, 2]. The model is sketched in Fig. 9 and is now qualitatively summarised.

Figure 9
figure 9

Cartoon showing various evolutionary scenarios involving binary pulsars.

Starting with a binary star system, a 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 [98, 31]. In addition, the fraction of surviving binaries is affected by the magnitude and direction of any impulsive “kick” velocity the neutron star receives at birth [98, 18]. Those binary systems that disrupt produce a high-velocity isolated neutron star and an OB runaway star [35]. The high binary disruption probability during the explosion explains, qualitatively at least, why so few normal pulsars have companions. Over the next 107 yr or so after the explosion, the neutron star may be observable as a normal radio pulsar spinning down to a period ≳ several seconds. After this time, the energy output of the star diminuishes to a point where it no longer produces significant radio emission.

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 of life as a pulsar by accreting matter and therefore angular momentum at the expense of the orbital angular momentum of the binary system [2]. The term “recycled pulsar” is often used to describe such objects. During this accretion phase, the X-rays produced by the liberation of gravitational energy of the infalling matter 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 [31].

The high-mass companions are massive enough to explode 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 [102], a 59-ms radio pulsar with a characteristic age of ∼108 yr which orbits its companion every 7.75 hr [239, 240]. In this formation scenario, PSR B1913+16 is an example of the older, first-born, neutron star that has subsequently accreted matter from its companion. So far there are no clear examples of systems where the second-born neutron star is observed as a radio pulsar. In the case of PSR B1820-11 [155], which may be an example, the mass of the companion is not well determined, so either a main-sequence [193] or a white dwarf companion [194] are plausible alternatives. This lack of observation of second-born neutron stars as radio pulsars is probably reasonable when one realises that the observable lifetimes of recycled pulsars are much larger than those of normal pulsars. As discussed in § 3.4.1, double neutron star binary systems are very rare in the Galaxy — another indication that the majority of binary systems get disrupted when one of the components explodes as a supernova. Systems disrupted after the supernova of the secondary form a mildly-recycled isolated pulsar and a young pulsar formed during the explosion of the secondary.

Although no system has so far been found in which both neutron stars are visible as radio pulsars, timing measurements of three systems show that the companion masses are 1.4M⊙ — as expected for neutron stars [215]. In addition, no optical counterparts are seen. Thus, we conclude that these unseen companions are neutron stars that are either too weak to be detected, no longer active as radio pulsars, or their emission beams do not intersect our line of sight. The two known young radio pulsars with main sequence companions massive enough to explode as a supernova probably represent the intermediate phase between high-mass X-ray binaries and double neutron star systems [110, 115].

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 [2]. This model has gained strong support in recent years from the discoveries of quasi-periodic kHz oscillations in a number of low-mass X-ray binaries [266], as well as Doppler-shifted 2.49-ms X-ray pulsations from the transient X-ray burster SAX J1808.4-3658 [267, 53]. 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 [25, 27, 139, 138]. Perhaps the best example is the white dwarf companion to the 5.25-ms pulsar J1012+5307 [177, 133]. This 19th magnitude white dwarf is bright enough to allow measurements of its surface gravity and orbital velocity [257].

The range of white dwarf masses observed is becoming broader. Since this article originally appeared in 1998, the number of “intermediate-mass binary pulsars” [43] has grown significantly [49]. These systems are distinct to the “classical” millisecond pulsar-white dwarf binaries like PSR J1012+5307 in several ways: (1) the spin period of the radio pulsar is generally longer (9–200 ms); (2) the mass of the white dwarf is larger (typically close to 1M); (3) the orbital eccentricity, while still essentially circular, is often significantly larger (∼10-3). It is not presently clear whether these systems originated from either low- or high-mass X-ray binaries. It was suggested by van den Heuvel [254] that they have more in common with high-mass systems, the difference being that the secondary star was not sufficiently massive to explode as a supernova. Instead it formed a white dwarf. Detailed studies of this sub-population of binary pulsars are required for further understanding in this area.

Another relatively poorly understood area is the existence of solitary millisecond pulsars in the Galactic disk (which comprise just under 20% of all Galactic millisecond pulsars). Although it has been proposed that the millisecond pulsars have got rid of their companion by ablation, as appears to be happening in the PSR B1957+20 system [83], it is not clear whether the time-scales for this process are feasible. There is some observational evidence that suggests that solitary millisecond pulsars are less luminous than binary millisecond pulsars [20, 122]. If confirmed by future discoveries, this would need to be explained by any viable evolutionary model.

Space 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 large velocities came from optical observations of the Crab pulsar in 1968 [251], 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 [142, 21, 78, 92]. These data imply a broad velocity spectrum ranging from 0 to over 1000 km s-1 [148].

Such large velocities are perhaps not surprising, given the violent conditions under which neutron stars are formed. Shklovskii [217] 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. In addition, if the neutron star progenitor was a member of a binary system prior to the explosion, the pre-supernova orbital velocity will also contribute to the resulting speed of the newly-formed pulsar. High-velocity pulsars born close to the Galactic plane quickly migrate to higher Galactic latitudes. This migration is seen in Fig. 10, a dynamical simulation of the orbits of 100 neutron stars in a model of the Galactic gravitational potential. Given such a broad velocity spectrum, as much as half of all pulsars will eventually escape the gravitational potential of the Galaxy and end up in intergalactic space [148, 58].

Figure 10
figure 10

A simulation following the motion of 100 pulsars in a model gravitational potential of our Galaxy for 200 Myr. The view is edge-on, i.e. the horizontal axis represents the Galactic plane (30 kpc across) while the vertical axis represents ±10 kpc from the plane. This snapshot shows the initial configuration of young neutron stars. To see the movie in action go to the online version of this article on http://www.livingreviews.org/Articles/Volume4/2001-5lorimer.

Based on the proper motion data, recent studies have demonstrated that the mean birth velocity of normal pulsars is ∼450 km s-1 ([148, 132, 58, 84]; see, however, also [93, 91]). This is significantly larger than the velocities of millisecond and binary pulsars. Recent studies suggest that their mean birth velocity is likely to be in the range ∼80–140 km s-1 [129, 57, 152]. 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 could not have survived had the neutron star received a substantial kick velocity.

Searching for pulsars

Pulsar searching is, conceptually at least, a rather simple process — the detection of a dispersed, periodic signal hidden in a noisy time series collected using a large radio telescope. The search strategy is to form a large number of time series for different trial DM values. These data are then analysed for periodic signals. We give here a brief description of the basic search techniques. More detailed discussions can be found elsewhere [141, 178, 130].

A schematic pulsar search is shown in Fig. 11. The finite bandwidth is split up into a number of channels, typically using a filterbank or a correlator (see e.g. [13]), either of which usually provides a much finer frequency channelisation than the eight channels shown for illustrative purposes in Fig. 11. The channels are then de-dispersed (see § 3.1.2) to form a single noisy time series. An efficient way to find a periodic signal in these data is to take the Fast Fourier Transform (FFT) and plot the resulting amplitude spectrum. For a narrow pulse the spectrum will show a family of harmonics. To detect weaker signals still, a harmonic summing technique is usually implemented [141]. The best candidates are saved and the whole process is repeated for another trial DM.

Figure 11
figure 11

Schematic summarising the essential steps in a “standard” pulsar search.

After the data have been processed for a suitable range of DM, a list of pulsar candidates is compiled and the raw time series data are folded modulo each candidate period. In practise the analysis is often hampered by the presence of periodic interference sources which can often look very “pulsar-like”. Although interference excision schemes (usually based on coincidence analyses of data taken from different points on the sky) work fairly well, interference is an ever-increasing problem in radio astronomy and considerable efforts are required to carry out sensitive searches.

Where to look for binary and millisecond pulsars

Following the above discussions on demography and evolution, it is instructive to briefly summarize the rationale behind the major searches being carried out at the present time.

All-sky searches

The oldest radio pulsars form a virialised population of stars oscillating in the Galactic gravitational potential. The scale height for such a population is at least 500 pc, about 10 times that of the massive stars which populate the Galactic plane. Since the typical ages of millisecond pulsars are several Gyr or more, we expect, from our vantage point in the Galaxy, to be in the middle of an essentially isotropic population of nearby sources. All-sky searches for millisecond pulsars at high Galactic latitudes have been very effective in probing this population. Much of the initial interest and excitement in this area was started at Arecibo when Wolszczan discovered two classic recycled pulsars at high latitudes: the neutron star binary B1534+12 [270] and the planets pulsar B1257+12 [273]. Surveys carried out at Arecibo, Parkes, Jodrell Bank and Green Bank by others in the 1990s found many other millisecond and recycled pulsars in this way. Camilo has written several excellent reviews of these surveys [41, 44, 45]. See also Tables 2, 3 and 4 in the appendix.

Searches close to the plane of our Galaxy

Young pulsars are most likely to be found near to their place of birth, close to the Galactic plane. This is the target region of one of the Parkes multibeam surveys and has already resulted in the discovery of around 600 new pulsars [48, 159], almost half the number currently known! Such a large haul inevitably results in a number of interesting individual objects such as: PSR J1141-6545, a young pulsar in a relativistic 4-hr orbit around a white dwarf [116]; PSR J1740-3052, a young pulsar orbiting an ∼11M star (probably a giant [159]); several intermediate-mass binary pulsars [49] and a likely double neutron star system, PSR J1811-1736 [147].

Searches at intermediate Galactic latitudes

In order to probe more deeply into the population of millisecond and recycled pulsars than possible at high Galactic latitudes, Edwards et al. have recently completed a survey of intermediate latitudes with the Parkes multibeam system [74, 75]. The results of this survey are extremely exciting — 58 new pulsars including 8 relatively distant recycled objects. Two of the new recycled pulsars from this survey recently announced by Edwards & Bailes [75] are mildly relativistic neutron star-white dwarf binaries. An analysis of the full results from this survey should significantly improve our knowledge on the Galaxy-wide population and birth-rate of millisecond pulsars.

Targeted searches of globular clusters

Globular clusters have long been known to be breeding grounds for millisecond and binary pulsars. The main reason for this is the high stellar density in globular clusters relative to most of the rest of the Galaxy. As a result, lowmass X-ray binaries are almost 10 times more abundant in clusters than in the Galactic disk. In addition, exchange interactions between binary and multiple systems in the cluster can result in the formation of exotic binary systems. Since a single globular cluster usually fits well within a single telescope beam, deep targeted searches can be made. Once the DM of a pulsar is known in a globular cluster, the DM parameter space for subsequent searches is essentially fixed. This allows computation power to be invested in so-called acceleration searches for short-period binary systems (see § 3.1.3). To date, searches have revealed 47 pulsars in globular clusters (see Table 5 in the appendix for a list and the review by Kulkarni & Anderson [123]). Highlights include the double neutron star binary in M15 [195] and a low-mass binary system with a 95-min orbital period in 47 Tucanae [47], one of 20 millisecond pulsars currently known in this cluster alone [81]. On-going surveys of clusters continue to yield new discoveries [199, 62].

Going further

Two excellent graduate-level monographs are available: the classic 1970s text Pulsars by Manchester & Taylor [163] and the more up-to-date Pulsar Astronomy by Lyne & Smith [156]. Those wishing to approach the subject from a more theoretical viewpoint are advised to read Michel’s The Theory of Neutron Star Magnetospheres [167] and The Physics of the Pulsar Magnetosphere by Beskin, Gurevich & Istomin [30]. Our summary of evolutionary aspects serves merely as a primer to the vast body of literature available. The reader is referred to the excellent review by Bhattacharya & van den Heuvel [31] for further insights. For an excellent overview of pulsar distance measurements and their implications, see the review by Weisberg [263].

Pulsar resources available on the Internet are continually becoming more extensive and useful. Good starting points for pulsar-surfers are the pages maintained at Arecibo [173], Berkeley [10], Bonn [165], Jodrell Bank [105], Princeton [196], Swinburne [230] and Sydney [197].

The Galactic Pulsar Population

Soon after the discovery of pulsars, it was realised that the observed sample is heavily biased towards the brighter objects that are the easiest to detect. What we observe therefore most likely represents only the tip of the iceberg of a much larger underlying population [88]. The extent to which the sample is incomplete is well demonstrated by the projection of pulsars onto the Galactic plane and their cumulative number distribution as a function of distance shown in Fig. 12. Although the clustering of sources around the Sun seen in the left panel of Fig. 12 would be consistent with Ptolemy’s geocentric picture of the heavens, it is clearly at variance with what we now know about the Galaxy, where the massive stars show a radial distribution about the Galactic center.

Figure 12
figure 12

Left: The sample of radio pulsars from the Princeton catalog [ 196] projected onto the Galactic plane. The Galactic center is at (0, 0) and the Sun is at (-8.5, 0). Right: Cumulative number of observed pulsars (solid line) as a function of projected distance, d. The dashed line shows the expected distribution for a model population (see text).

The extent to which the pulsar sample is incomplete is shown in the right panel of Fig. 12 where the cumulative number of pulsars is plotted as a function of the projected distance from the Sun. The observed distribution is compared to the expected distribution for a simple model population in which there are errors in the distance scale, but no selection effects. We see that the observed sample becomes strongly deficient in terms of the number of sources for distances beyond a few kpc. We now discuss the main selection effects that hamper the detection of pulsars in some detail.

Selection effects in pulsar searches

The inverse square law and survey thresholds

The most prominent selection effect at play in the observed pulsar sample is the inverse square law, i.e. for a given intrinsic luminosityFootnote 2, 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 falls below the flux thresholds Smin of most surveys. Following [71], we write:

$$\begin{array}{*{20}c} {S_{\min } = SNR_{\min } \left( {\frac{{T_{rec} + T_{sky} }} {K}} \right)\left( {\frac{G} {{K\,Jy^{ - 1} }}} \right)^{ - 1} \left( {\frac{{\Delta \nu }} {{MHz}}} \right)^{ - 1/2} \left( {\frac{{t_{\operatorname{int} } }} {s}} \right)^{ - 1/2} \cdot} \\ {\left( {\frac{W} {{P - W}}} \right)^{1/2} mJy,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;} \\ \end{array} $$
((3))

where SNRmin is the threshold signal-to-noise ratio, Trec and Tsky are the receiver and sky noise temperatures, G is the gain of the antenna, Δν is the observing bandwidth, tint is the integration time, W is the detected pulse width and P is the pulse period.

Pulse dispersion and scattering

It follows from Equation (3) that the sensitivity decreases as W/(P−W) increases. Also note that if WP, the pulsed signal is smeared into the background emission and is no longer detectable, regardless of how luminous the source may be. The detected pulse width W will be broader than the intrinsic value largely as a result of pulse dispersion and scattering by free electrons in the interstellar medium. As discussed above, the dispersive smearing scales as Δν/ν3, 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 profile. This process is known as incoherent dedispersion.

The smearing across the individual frequency 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. A simple scattering model is shown in Fig. 13 in which the scattering electrons are assumed to lie in a thin screen between the pulsar and the observer [210].

Figure 13
figure 13

Pulse scattering caused by irregularities in the interstellar medium. The difference in path lengths and therefore in arrival times of the scattered rays results in a “scattering tail” in the observed pulse profile which lowers its signal-to-noise ratio.

Dispersion and scattering are most severe for distant pulsars in the inner Galaxy where 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 [56, 109] compared to the usual 400 MHz search frequency. An added bonus for such observations is the reduction in Tsky, since the spectral index of the non-thermal Galactic emission is about -2.8 [126]. Pulsars themselves have steep radio spectra. Typical spectral indices are -1.6 [136], so that flux densities are roughly an order of magnitude lower at 1400 MHz compared to 400 MHz. Fortunately, this can usually be compensated for by the use of larger receiver bandwidths at higher radio frequencies. For example, the 1370-MHz system at Parkes has a bandwidth of 288 MHz [147] compared to the 430-MHz system, where 32 MHz is available [160].

Orbital acceleration

Standard pulsar searches use Fourier techniques to search for a-priori unknown periodic signals and usually assume that the apparent pulse period remains constant throughout the observation. For searches with integration times much greater than a few minutes this assumption is only valid for solitary pulsars, or those in binary systems where the orbital periods are longer than about a day. For shorter-period binary systems, as noted by Johnston & Kulkarni [106], the Doppler-shifting of the period results in a spreading of the signal power over a number of frequency bins in the Fourier domain, leading to a reduction in signal-to-noise ratio. An observer will perceive the frequency of a pulsar to shift by an amount aT/(Pc), where a is the (assumed constant) line-of-sight acceleration during the observation of length T, P is the (constant) pulse period in its rest frame and c is the speed of light. Given that the width of a frequency bin is 1/T, we see that the signal will drift into more than one spectral bin if aT2/(Pc)>1. Survey sensitivities to rapidly-spinning pulsars in tight orbits are therefore significantly compromised when the integration times are large.

As an example of this effect, as seen in the time domain, Fig. 14 shows a 22.5-min search mode observation of Hulse & Taylor’s famous binary pulsar B1913+16 [102, 239, 240]. Although this observation covers only about 5% of the orbit (7.75 hr), the effects of the Doppler smearing on the pulse signal are very apparent. While the standard search code (seeking constant periodicity) nominally detects the pulsar with a signal-to-noise ratio of 9.5 for this observation, it is clear that the Doppler shifting of the pulse period seen in the individual sub-integrations results in a significant reduction in signal-to-noise.

Figure 14
figure 14

Left: A 22.5-min Arecibo observation of the binary pulsar B1913+16. The assumption that the pulsar has a constant period during this time is clearly inappropriate given the drifting in phase of the pulse during the integration (grey scale plot). Right: The same observation after applying an acceleration search. This shows the effective recovery of the pulse shape and a significant improvement in the signal-to-noise ratio.

It is clearly desirable to employ a technique to recover the loss in sensitivity due to Doppler smearing. One such technique, the so-called “acceleration search” [168], assumes the pulsar has a constant acceleration during the integration. Each time series can then be re-sampled to refer it to the frame of an inertial observer using the Doppler formula to relate a time interval τ in the pulsar frame to that in the observed frame at time t, as τ(t)∝(1+at/c). Searching over a range of accelerations is desirable to find the time series for which the trial acceleration most closely matches the true value. In the ideal case, a time series is produced with a signal of constant period for which full sensitivity is recovered (see right panel of Fig. 14). Anderson et al. [5] used this technique to find PSR B2127+11C, a double neutron star binary in M15 which has parameters similar to B1913+16. Camilo et al. [47] have recently applied the same technique to 47 Tucanae to discover 9 binary pulsars, including one in a 96-min orbit around a low-mass (0.15M) companion. This is currently the shortest binary period for any known radio pulsar.

For the shortest orbital periods, the assumption of a constant acceleration during the observation clearly breaks down. Ransom et al. [199] have developed a particularly efficient algorithm for finding binaries whose orbits are so short that many orbits can take place during an integration. This phase modulation technique exploits the fact that the pulses are modulated by the orbit to create a family of periodic sidebands around the nominal spin period of the pulsar. This technique has already been used to discover a 1.7-hr binary pulsar in NGC 6544 [199]. The existence of these short-period radio pulsar binaries, as well as the 11-min X-ray binary X1820–303 in NGC 6624 [226], implies that there must be many more short-period binaries containing radio or X-ray pulsars in globular clusters that are waiting to be discovered by more sensitive searches.

Correcting the observed pulsar sample

Scale factor determination

Now that we have a flavour for the variety and severity of the selection effects that plague the observed sample of pulsars, how do we decouple these effects to form a less biased picture of the true population of objects? A very useful technique, first employed by Phinney & Blandford and Vivekanand & Narayan [191, 259], is to define a scaling factor ξ 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{{\iint_{Galaxy} {\sum ( R,z)RdRdz}}}{{\iint_{P,L} {\sum ( R,z)RdRdz}}}.$$
((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 approach is similar to the classic V/Vmax technique first used to correct observationally-biased samples of quasars [211].

This technique 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 [170]. For a sample of Nobs observed pulsars above a minimum luminosity Lmin, the total number of pulsars in the Galaxy with luminosities above this value is simply

$$\sum\limits_{i = 1}^{{N_{{\rm{obs}}}}} {\frac{{{\xi _i}}}{{{f_i}}},} $$
((5))

where f is the model-dependent “beaming fraction” discussed below in § 3.2.3. Monte Carlo simulations of the pulsar population incorporating the aforementioned selection effects have shown this method to be reliable, as long as Nobs is reasonably large [131].

The small-number bias

For small samples of observationally-selected objects, the detected sources are likely to be those with larger-than-average luminosities. The sum of the scale factors (Equation (5)), therefore, will tend to underestimate the true size of the population. This “small-number bias” was first pointed out by Kalogera et al. [112, 113] for the sample of double neutron star binaries where we know of only three clear-cut examples (§ 3.4.1). Only when the number of sources in the sample gets past 10 or so does the sum of the scale factors become a good indicator of the true population size.

The beaming correction

The “beaming fraction” f in Equation (5) is simply 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 [238]. Observational evidence suggests that shorter period pulsars have wider beams and therefore larger beaming fractions than their long-period counterparts [171, 149, 32, 231]. It must be said, however, that a consensus on the beaming fraction-period relation has yet to be reached. This is shown in Fig. 16 where we compare the period dependence of f as given by a number of models. 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.

Figure 16
figure 16

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

When most of these beaming models were originally proposed, the sample of millisecond pulsars was ≲5 and 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 by Kramer et al. [122] suggests that the beaming fraction of millisecond pulsars lies between 50 and 100%.

The population of normal and millisecond pulsars

Luminosity distributions and local number estimates

The most recent use of the scale factor approach to derive the characteristics of the true normal and millisecond pulsar populations is based on the sample of pulsars within 1.5 kpc of the Sun [152]. 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 distributions obtained from this analysis are shown in Fig. 17. For the normal pulsars, integrating the corrected distribution above 1 mJy kpc2 and dividing by π×(1.5)2 kpc2 yields a local surface density, assuming Biggs’ beaming model [32] of 156±31 pulsars kpc-2. The same analysis for the millisecond pulsars, assuming a mean beaming fraction of 75% [122], leads to a local surface density of 38±16 pulsars kpc-2 for luminosities above 1 mJy kpc2.

Figure 17
figure 17

Left: The corrected luminosity distribution (solid histogram with error bars) for normal pulsars. The corrected distribution before the beaming model has been applied is shown by the dot-dashed line. Right: The corresponding distribution for millisecond pulsars. In both cases, the observed distribution is shown by the dashed line and the thick solid line is a power law with a slope of -1. The difference between the observed and corrected distributions highlights the severe under-sampling of low-luminosity pulsars.

Galactic population and birth-rates

Integrating the local surface densities of pulsars over the whole Galaxy requires a knowledge of the presently rather uncertain Galactocentric radial distribution [153, 107]. One approach is to assume that pulsars have a radial distribution similar to that of other stellar populations and scale the local number density with this distribution to estimate the total Galactic population. The corresponding local-to-Galactic scaling is 1000±250 kpc2 [200]. With this approach 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. From the P diagram in Fig. 7, we infer a typical lifetime for normal pulsars of ∼107 yr, corresponding to a Galactic birth rate of ∼1 per 60 yr — consistent with the rate of supernovae [253]. As noted in § 2.4, the millisecond pulsars are much older, with ages close to that of the Universe τu (we assume here τu=10 Gyr [104]). Taking the maximum age of the millisecond pulsars to be τu, we infer a mean birth rate of at least 4×10-6 yr-1. This is consistent, within the uncertainties, with the birth-rate of low-mass X-ray binaries [135].

Implications for gravitational wave detectors

The estimates of the local surface density of active pulsars allow us to deduce the likely distance of the nearest neutron star to 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. 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 [213]. According to Thorne [244], currently planned detectors will be able to detect neutron stars with ellipticities greater than 7.5×10-11(Pd)2, where P is the rotation period in ms and d is the distance in kpc. The recent probable detection of free precession in the radio pulsar B1828-11 [221] does indicate that ellipticities exist in neutron stars so that nearby objects may be continuous sources of gravitational radiation.

Thus, in order to detect small ellipticities, nearby sources with short spin periods are required. One of the best candidates is the nearby 5.75-ms pulsar J0437-4715 [108]. At a distance of 178±26 pc [208] this is currently the closest known millisecond pulsar to the Earth. The closest known neutron star is RX J185635-3754 discovered in the ROSAT all-sky survey [261]. Multi-epoch HST observations show that this isolated neutron star is located at a distance of 61±9 pc [260]. In keeping with other radio-quiet isolated neutron stars, the period of this pulsar is likely to be several seconds [176].

The population of relativistic binaries

Of particular interest to the astronomical community at large are the numbers of relativistic binary systems in the Galaxy. Systems involving neutron stars are: double neutron star binaries, neutron star-white dwarf binaries and neutron star-black hole binaries. Interest in these systems is two-fold: (1) to test the predictions of general relativity against alternative theories of strong-field gravity using the radio pulsar as a highly stable clock moving in the strong gravitational field; (2) to detect strong gravitational wave emission from coalescing binaries with upcoming gravitational-wave observatories like GEO600 and LIGO, VIRGO and TAMA [244, 212]. For a Living Reviews article on this new generation of interferometric gravitational wave detectors see [100].

Although no radio pulsar in orbit around a black hole companion has so far been observed, we now know of several double neutron star and neutron-star white dwarf binaries which will merge due to gravitational wave emission within a reasonable time-scale. The merging time τg of a binary system containing two compact objects due to the emission of gravitational radiation can be calculated from the following formula which requires only the component masses and current orbital period Pb and eccentricity e:

$${\tau _g} \simeq {10^7}{\rm{yr}}{\left( {\frac{{{P_{\rm{b}}}}}{{{\rm{hr}}}}} \right)^{8/3}}{\left( {{{\frac{{{m_1} + {m_2}}}{M}}_ \odot }} \right)^{1/3}}{\left( {{{\frac{\mu }{M}}_ \odot }} \right)^{ - 1}}{(1 - {e^2})^{7/2}}.$$
((6))

Here m1,2 are the masses of the two stars and μ=m1m2/(m1+m2). This formula is a good analytic approximation (within a few percent) to the numerical solution of the exact equations for τd in the original papers by Peters & Mathews [187, 188]. In the following subsections we review current knowledge on the population sizes and merging rates of such binaries where one component is visible as a radio pulsar.

Double neutron star binaries

As noted in § 2.4, double neutron star (DNS) binaries are expected to be rare since the binary system has to survive two supernova explosions. This expectation is certainly borne out by radio pulsar searches which have revealed only three certain DNS binaries so far: PSRs B1534+12 [270], B1913+16 [102] and B2127+11C [195]. Although we cannot see the companion neutron star in any of these systems, we are “certain” of the identification from the precise measurements of the component masses via relativistic effects measured in pulsar timing observations (see § 4.1). The spin and orbital parameters of these pulsars are listed in Table 1. Also listed in Table 1 are three further DNS candidates with eccentric orbits and large mass functions but for which there is presently not sufficient component mass information to confirm their nature.

Table 1 Known DNS binaries and candidates. Listed are the pulse period P, the orbital period Pb, the orbital eccentricity e, the pulsar characteristic age τc, and the expected binary coalescence time-scale τg due to gravitational wave emission calculated from Equation (6). Cases for which τg is a factor of 100 or more greater than the age of the Universe are listed asτu. To distinguish between definite and candidate DNS systems, we also list whether the masses of both components have been determined.

Despite the uncertainties in identifying DNS binaries, for the purposes of determining the Galactic merger rate, the systems for which τg is less than τu (i.e. PSRs B1534+12, B1913+16 and B2127+11C) are primarily of interest. Of these PSR B2127+11C is in the process of being ejected from the globular cluster M15 [195, 192] and is thought to make only a negligible contribution to the merger rate [190]. The general approach with the remaining two systems is to derive scale factors for each object (as outlined in § 3.2.1) and then divide these by a reasonable estimate for the lifetime. In what follows we summarize the main studies of this kind. The most comprehensive investigation of the DNS binary population to date is the recent study by Kalogera et al. (hereafter KNST; [113]).

As discussed in § 3.2.1, scale factors are dependent on the assumed pulsar distribution. The key parameter here is the scale height of the population with respect to the Galactic plane which itself is a function of the velocity distribution of the population. KNST examined this dependence in detail and found scale heights in the range 0.8–1.7 kpc. Based on this range, KNST revised earlier scale factor estimates [60] to 145–200 for B1534+12 and 45–60 for PSR B1913+16. As mentioned in § 3.2.2 scale factors calculated from a small sample of objects are subject to a significant bias. KNST find the bias in their sample to be anywhere between 2 and 200. This boosts the scale factors to the range 190–40000 for B1534+12 and 90–12000 for B1913+16.

The above scale factors also require a beaming correction. As noted in § 3.2.3, current radio pulsar beaming models vary considerably. Fortunately, for the two pulsars under consideration, detailed studies of the beam sizes [9, 119, 262] lead KNST to conclude that both pulsars beam to only about a sixth of the entire sky. The beaming-corrected numbers suggest a total of between 1680 and 312,000 active DNS binaries in our Galaxy. Many of these systems will be extremely faint objects. These estimates are dominated by the small-number bias factor. KNST’s study highlights the importance of this effect.

Some debate exists about what is the most reasonable estimate of the lifetime. Phinney [190] defines this as the sum of the pulsar’s spin-down age plus τg defined above. A few years later, van den Heuvel and myself argued [255] that a more likely estimate can be obtained by appealing to steady-state arguments where we expect sources to be created at the same rate at which they are merging. The mean lifetime was then found to be about three times the current spin-down age. This argument does, however, depend on the luminosity evolution of radio pulsars which is currently only poorly understood. Arzoumanian, Cordes & Wasserman [7] used kinematic data to constrain the most likely ages of the DNS binaries. They note that the remaining detectable lifetime should also take account of the reduced detectability at later epochs due to acceleration smearing as the DNS binary becomes more compact due to gravitational wave emission. KNST concluded that the lifetimes are dominated by the latter time-scale which, following Arzoumanian et al., they took to be the time for the orbital period to halve. The resulting lifetimes are 2.5×109 yr for B1534+12 and 2.5×108 yr for B1913+16.

Taking these number and lifetime estimates, KNST find the Galactic merger rate of DNS binaries to range between 3×10-6 and 4×10-4 yr-1. Extrapolating this number out to include DNS binaries detectable by LIGO in other galaxies á la Phinney [190], KNST find the expected event rate to be <0.25 yr-1 for LIGO-I and 2–1300 yr-1 for LIGO-II. Thus, despite the uncertainties, it seems that the prospects for detecting gravitational-wave emission from DNS inspirals in the near future are most promising.

White dwarf-neutron star binaries

The population of white dwarf-neutron star (WDNS) sources containing a young radio pulsar has only recently been confirmed by observers following the identification of a white dwarf companion to the binary pulsar B2303+46 by van Kerkwijk et al. [256]. Previously, this eccentric binary pulsar was thought to be an example of a DNS binary in which the visible pulsar is the second-born neutron star [229, 143]. The optical identification rules this out and now strongly suggests a scenario in which the white dwarf was formed first. In this case, material was transfered onto the secondary during the giant phase of the primary so that the secondary became massive enough to form a neutron star [256, 232].

PSR B2303+46 has a long orbital period and does not contribute significantly to the overall merger rate of WDNS binaries. The new discovery of PSR J1141-6545 [116], which will merge due to gravitational-wave emission within 1.3 Gyr, is suggestive of a large population of similar binaries. This is particularly compelling when one considers that the radio lifetime of the visible pulsar is only a fraction of total lifetime of the binary before coalescence due to gravitational-wave emission. Edwards & Bailes [75] estimate there to be 850 WDNS binaries within 3 kpc of the Sun which will merge within τu.

Population syntheses by Tauris & Sennels [232] suggest that the formation rate of WDNS binaries is between 10–20 times that of DNS binaries. Based on the merging rate estimates for DNS binaries discussed in the previous section, this translates to a merging rate of WDNS binaries of between 3×10-5 yr-1 and 8×10-3 yr-1. In summary, although statistics are necessarily poor at this stage, coalescing WDNS binaries look to be very promising sources for gravitational wave detectors.

Going further

Studies of pulsar population statistics represent a large proportion of the pulsar literature. During this section we have tried to cite many of the key papers in this field. Good starting points for further reading can be found in other review articles [167, 31, 156]. Our coverage of compact object coalescence rates has concentrated on empirical methods and we have hopefully convinced the reader that these are fair and straightforward.

An alternative approach is to undertake a full-blown Monte Carlo simulation of the most likely evolutionary scenarios described in § 2.4.3. In this “scenario-machine” approach, a population of primordial binaries is synthesized with a number of underlying distribution functions: primary mass, binary mass ratio, orbital period distribution etc. The evolution of both stars is then followed to give a predicted sample of binary systems of all the various types. Since the full range of binary parameters is known, the merger rates of each type of binary are then automatically predicted by this model without the need to debate what the likely coalescence times will be. Selection effects are not normally taken into account in this approach. The final census is usually normalized to the star formation rate. Numerous examples of the scenario-machine approach (most often to populations of binaries where one or both members are NSs) can be found in the literature [69, 204, 252]. These include the widely-cited code, developed by Lipunov and collaborators to perform population syntheses of binary stars [227]. Although extremely instructive, the uncertain assumptions about initial conditions, the physics of mass transfer and the kicks applied to the compact object at birth result in a wide range of predicted event rates which are currently broader than the empirical methods [112]. Ultimately, the detection statistics from the gravitational wave detectors could provide far tighter constraints on the DNS merging rate than the pulsar surveys from which these predictions are made.

An excellent overview of gravitational-wave astronomy and the detection of gravitational waves from inspiraling binaries is presented by Thorne in his presentation at the centennial meeting of the American Physical Society which is available on-line [243].

Pulsar Timing

It became clear soon after their discovery that pulsars are excellent celestial clocks. In the original discovery paper [97], the period of the first pulsar to be discovered, PSR B1919+21, was found to be stable to one part in 107 over a time-scale of a few months. Following the discovery of the millisecond pulsar B1937+21 in 1982 [17] it was demonstrated that its period could be measured to one part in 1013 or better [65]. This unrivaled stability leads to a host of applications including time keeping, probes of relativistic gravity and natural gravitational wave detectors.

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. 18 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 and added to form a mean pulse profile.

Figure 18
figure 18

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 Global Positioning System of satellites (GPS; [64]). The TOA of this mean pulse is then defined as the arrival time of some fiducal 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 the 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. Using Equation (3), we can express this as a fraction of the pulse period:

$$\frac{{ \in _{TOA} }} {P} \simeq \left( {\frac{{T_{rec} + T_{sky} }} {K}} \right)\left( {\frac{G} {{KJy^{ - 1} }}} \right)^{ - 1} \left( {\frac{{\Delta \upsilon }} {{MHz}}} \right)^{ - 1/2} \left( {\frac{{t_{\operatorname{int} } }} {s}} \right)^{ - 1/2} \left( {\frac{W} {P}} \right)^{3/2} .$$
((7))

In this expression Trec and Tsky are the receiver and sky noise temperatures, G is the gain of the antenna, Δν is the observing bandwidth, tint is the integration time, W is the pulse width and P is the pulse period (we assume WP). 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 receivers 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 [90] 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 signal processing can be done either on-line using finite impulse response filter devices [14] or off-line in software [219, 222]. The on-line approach allows for large bandwidths to be employed and real-time viewing of the data. Off-line reduction, while slow and computationally expensive, allows for more flexible data reduction schemes as well as periodicity searches to be carried out.

The maximum time resolution obtainable via coherent dedispersion is the inverse of the receiver bandwidth. For bandwidths of 10 MHz, this technique makes it possible to resolve features on time-scales as short as 100 ns. This corresponds to probing regions in the neutron star magnetosphere as small as 30 m!

The timing model

Ideally, in order to model the rotational behaviour of the neutron star, we require TOAs measured by an inertial observer. An observatory located on Earth experiences accelerations with respect to the neutron star due to the Earth’s rotation and orbital motion around the Sun and 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 [103] to transform the observed TOAs to this frame using a planetary ephemeris such as the JPL DE200 [224]. The transformation is summarised as the difference between barycentric (T) and observed (t) TOAs:

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

Here is the position of the Earth with respect to the barycentre, ŝ̱ 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 presently only measurable for four nearby millisecond pulsars [117, 46, 208]. The term Δtrel represents the Einstein and Shapiro corrections due to general relativistic effects within the solar system [16]. Since measurements are often carried out at different observing frequencies with different dispersive delays, the TOAs are generally referred to the equivalent time that would be observed at infinite frequency. This transformation corresponds to the term τtDM 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 to. The model pulse phase φ as a function of barycentric time is thus given by:

$$\varphi (\mathcal{T}) = \varphi _{\text{o}} + (\mathcal{T} - \mathcal{T}_o )\Omega _{\text{o}} + \frac{1} {2}(\mathcal{T} - \mathcal{T}_{\text{o}} )^2 \dot \Omega _{\text{o}} + \cdots ,$$
((9))

where φo is the pulse phase at To. 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.

A set of timing residuals for the nearby pulsar B1133+16 spanning almost 10 years is shown for illustrative purposes in Fig. 19. Ideally, the residuals should have a zero mean and be free from any systematic trends (Fig. 19a). Inevitably, however, due to our a-priori ignorance of the rotational parameters, the model needs to be refined in a bootstrap fashion. 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.

Figure 19
figure 19

Timing model residuals versus date for PSR B1133+16. Case (a) shows the residuals obtained from the best fitting model which includes period, the period derivative, position and proper motion. Case (b) is the result of setting the period derivative term to zero in this model. Case (c) shows the effect of a 1 arcmin error in the assumed declination. Case (d) shows the residuals obtained assuming zero proper motion. The lines in (b)–(d) show the expected behaviour in the timing residuals for each effect (see text).

From Equation (9), an error in the assumed Ωo results in a linear slope with time. A parabolic trend results from an error in Ωo (Fig. 19b). Additional effects will arise if the assumed position of the pulsar (the unit vector ŝ̱ in equation (8)) used in the barycentric time calculation is incorrect. A position error of just one arcsecond results in an annual sinusoid (Fig. 19c) 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. A proper motion produces an annual sinusoid of linearly increasing magnitude (Fig. 19d).

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. Timing observations of the original millisecond pulsar B1937+21, spanning almost 9 years (exactly 165,711,423,279 rotations!), measure a period of 1.5578064688197945±0.0000000000000004 ms [117, 114] 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.

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. 20. These are taken from the Jodrell Bank timing program [216]. Such “timing noise” is most prominent in the youngest of the normal pulsars [162, 59] and virtually absent in the much older millisecond pulsars [117]. While the physical processes of this phenomena are not well understood, it seems likely that it may be connected to superfluid processes and temperature changes in the interior of the neutron star [3], or processes in the magnetosphere [55, 54].

Figure 20
figure 20

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”.

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 dedispersion outlined in § 4.1. Much of the pioneering work in this area has been made by Joseph Taylor and collaborators at Princeton University [196]. From high quality observations made using the Arecibo radio telescope spanning almost a decade [206, 207, 117], 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. 21. This figure shows σz, a parameter closely resembling the Allan variance used by the clock community to estimate the stability of atomic clocks [233, 1]. Atomic clocks are known to have σ∼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 [117]. No such noise component is observed for 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. Timing observations of millisecond pulsars are discussed further in the context of the pulsar timing array in § 5.2.

Figure 21
figure 21

Fractional timing instabilities for PSRs B1855+09 and B1937+21 as a function of time. (After Kaspi, Taylor & Ryba 1994 [ 117].)

Binary pulsars and Kepler’s laws

For binary pulsars, the simple timing model introduced in § 4.2 needs to be extended to incorporate the additional radial acceleration 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 (Pb), projected semi-major orbital axis (ap sin i, see below), orbital eccentricity (e), longitude of periastron (ω) and the epoch of periastron passage (T0). This description, using five “Keplerian parameters”, is identical to that used for spectroscopic binary stars.

For spectroscopic binaries the orbital velocity curve shows the radial component of the star’s velocity as a function of time. The analogous plot for pulsars is the apparent pulse period against time. Two examples are given in Fig. 22.

Figure 22
figure 22

Orbital velocity curves for two binary pulsars. Left: PSR J1012+5307, a 5.25-ms pulsar in a 14.5-hour circular orbit around a low-mass white dwarf companion [ 177, 256, 125]. Right: PSR J1811-1736, a 104-ms pulsar in a highly eccentric 18.8-day orbit around a massive companion (probably another neutron star) [ 147].

Constraints on the mass of the orbiting companion can be placed by combining the projected semi-major axis ap sin i and the orbital period to obtain the mass function:

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

where G is the universal gravitational constant. Assuming a pulsar mass mp of 1.35Mȩ (see below), the mass of the orbiting companion mc can be estimated as a function of the (initially unknown) angle i between the orbital plane and the plane of the sky. The minimum companion mass mmin 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 i0 is p(<i0)=1−cos(i0). 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 m90 constrains the mass range between mmin and m90 at the 90% confidence level.

Post-Keplerian parameters

Although many of the presently known binary pulsar systems can be adequately described by a Keplerian orbit, there are several systems, including the original binary pulsar B1913+16, which exhibit relativistic effects that require an additional set of up to five “post-Keplerian” parameters. Within the framework of general relativity, these can be written [37] as:

$$\dot \omega = 3{\left( {\frac{{{P_{\rm{b}}}}}{{2\pi }}} \right)^{ - 5/3}}{({T_ \odot }M)^{2/3}}{(1 - {e^2})^{ - 1}},$$
((11))
$$\gamma = e\left( {\frac{{P_{\text{b}} }} {{2\pi }}} \right)^{1/3} T_ \odot ^{2/3} M^{ - 4/3} m_{\text{c}} (m_{\text{p}} + 2m_{\text{c}} ),$$
((12))
$${\dot P_{\rm{b}}} = - \left( {\frac{{192\pi }}{5}} \right){\left( {\frac{{{P_{\rm{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_{\rm{p}}}{m_{\rm{c}}}{M^{ - 1/3}},$$
((13))
$$r = {T_ \odot }{m_{\rm{c}}},$$
((14))
$$s = x\left( {\frac{{P_{\text{b}} }} {{2\pi }}} \right)^{ - 2/3} T_ \odot ^{ - 1/3} M^{2/3} m_{\text{c}} ^{ - 1} .$$
((15))

In addition to the symbols defined above for Equation (10), Mmp+mc, xapsin i/c, s≡sin i and TGM/c3≃4.925 μs. All masses are in solar units.

Measurements of post-Keplerian parameters for PSR B1913+16 have been carried out by Taylor and a number of collaborators over the years with steadily improving precision. The first of these parameters to be measured was the advance of the longitude of periastron ω̇. This measurement is analogous to the perihelion advance of Mercury [169]. For PSR B1913+16 this amounts to about 4.2 degrees per year [237], some 4.6 orders of magnitude larger than for Mercury. A measurement of ω̇ alone yields the total mass for this system, M=2.83M⊙, assuming this advance is due to general relativity.

Measurement of a second post-Keplerian parameter for B1913+16, γ (gravitational redshift and transverse Doppler shifts in the orbit), permits an unambiguous determination of mp, mc and i when combined with ω̇ and the five Keplerian parameters. The original measurements [236] have since been substantially refined [239, 240] and the mass of the pulsar and its unseen companion have been determined to be 1.442±0.003M and 1.386±0.003M respectively. Such phenomenal precision is a testament to the timing stability of radio pulsars as clocks, and the diligence of Taylor and collaborators in carrying out these long-term measurements. Similar mass measurements now exist for the two other double neutron star binary systems discussed in § 3.4.1: B1534+12 [220] and B2127+11C [66]. For a number of other systems, ω̇ measurements allow interesting constraints to be placed on the component masses [246, 179, 248].

An important general relativistic prediction for eccentric double neutron star systems is the orbital decay due to the emission of gravitational radiation (Ṗb in Equation (13)). Taylor et al. [236, 239, 240] 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 3.2 mm per orbit, is seen most dramatically as the gradually increasing shift in orbital phase for periastron passages with respect to a non-decaying orbit shown in Fig. 23. This figure includes recent Arecibo data taken in 1998 and 1999 following the upgrade of the telescope in the mid 1990s. The observations of the orbital decay, now spanning a 25-year baseline, are in agreement with general relativity at the level of about 0.5% and provide the first (indirect) evidence for the existence of gravitational radiation. Hulse and Taylor were awarded the Nobel prize in Physics in 1993 [241, 101, 234] in recognition of their discovery of this remarkable laboratory for testing general relativity.

Figure 23
figure 23

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.

For those binary systems which are oriented nearly edge-on to the line-ofsight, 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. The so-called “range” and “shape” of the Shapiro delay effect are parameterized by the last two post-Keplerian parameters r and s≡sin i that were introduced in Equations (14) and (15). 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 [206, 117, 46], and for the double neutron star binaries B1534+12 [220] and B1913+16 [240].

Geodetic precession

Shortly after the discovery of PSR B1913+16 it was realized that, if the spin axis of the visible pulsar was misaligned with the angular momentum axis of the binary system, the perturbing effect of the companion on the space-time around the radio pulsar would cause it to precess around the angular momentum axis [63, 76]. Within the framework of general relativity, the rate of precession Ωp was shown [22] to be

$${\Omega _{\rm{p}}} = \frac{{{{(2\pi )}^{5/3}}T_ \odot ^{2/3}{m_{\rm{c}}}(4{m_{\rm{p}}} + 3{m_{\rm{c}}})}}{{P_{\rm{b}}^{5/3}{{({m_{\rm{p}}} + {m_{\rm{c}}})}^{4/3}}(1 - {e^2})}},$$
((16))

where we assume the same notation used for the discussion in § 4.4 and § 4.5. Inserting the parameters of PSR B1913+16 yields Ωp=1.21 deg yr-1. The period of the precession is 297.5 yr. The observational consequence of geodetic precession is a secular change in the pulse profile as the line-of-sight cut through the emission beam changes (recall Fig. 5).

Early qualitative evidence for profile evolution due to this effect [236] was substantiated with long-term Arecibo measurements of component changes by Weisberg et al. [264]. Further changes were seen by Kramer with new Effelsberg data acquired in the 1990s [119]. In addition to relative amplitude variations, the expected changes in component separation for a hollow-cone beam model were also seen in the Effelsberg data. These observations are summarized in Fig. 24.

Figure 24
figure 24

Geodetic precession in the binary pulsar B1913+16 system: (a) Changes in the observed pulse shapes between 1981–1988 seen as a decrease in amplitude between the left and right components in the pulse profile [ 240, 264]; (b) relative heights of the two amplitudes plotted between 1975–1998 [ 119]; (c) component separation change [ 119].

In addition to the above results, there is now evidence for geodetic precession in the other classic neutron star binary, PSR B1534+12 [9, 223]. Although geodetic precession in binary pulsars is another successful test of general relativity (albeit at a lower precision than e.g. orbital decay measurements), what is perhaps more interesting are the various consequences it has. Geodetic precession only occurs when the spin and orbital axes are misaligned [22]. This is most likely to occur if the neutron star received an impulsive “kick” velocity at birth (§ 2.4.4). Wex et al. [265] have investigated the B1913+16 observations and find that the kick magnitude was at least 250 km s-1 and was directed almost perpendicular to the spin axis of the neutron star progenitor. This places stringent constraints on any kick mechanism. Detailed monitoring of the pulse profile and polarization properties now underway [262, 120] will allow the first map of the emission beam of a neutron star to be made. This has important implications for the various beaming models described in § 3.2.3. There are already indications that the beam is circular [120].

The current results predict that B1913+16 will completely precess out of the line of sight by around 2025 and re-appear some 240 years later [119]. Although we shall lose a most treasured pulsar, we can take comfort from the fact that other pulsars will precess into our field of view. Perhaps one example is the newly-discovered relativistic binary J1141-6545 [116] discussed in § 2.6.2 and § 3.4.2. This relatively bright object was apparently missed by two previous searches during the early 1990s [109, 160, 152].

Going further

This chapter has outlined past and present progress in a number of areas related to pulsar timing. For further details on the technical details and prospects of pulsar timing, the interested reader is referred to a number of excellent review articles [269, 16, 233, 23, 24]. Two freely available software packages which are routinely used for time-of-arrival analyses by the pulsar community are available, viz. TEMPO [240, 196] and TIMAPR [72, 166]. These packages are based on more detailed versions of the timing model outlined in § 4.2. An up-to-date list summarising the various timing programmes is kept by Don Backer [10]. An audio file and slides from a lecture on pulsar timing presented by Backer at the centennial meeting of the American Physical society is also available online [11]. Another relevant lecture from that meeting is Will’s presentation [268] on tests of Einstein’s relativity which includes an excellent overview of Taylor and Weisberg’s measurements of PSR B1913+16. Kramer [120] has written a lucid review article discussing measurements of geodetic precession in binary pulsars and their implications.

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 three Earth-mass bodies [273, 186, 272]. Subsequent measurements of B1257+12 were even able to measure resonance interactions between two of the planets [271], 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 [247, 15, 245]. For detailed reviews of these systems, and their implications for planetary formation scenarios, the interested reader is referred to [189].

Pulsars as Gravitational 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 [185]. The idea to use pulsars as natural detectors of gravitational waves was first explored independently by both Sazhin and Detweiler in the late 1970s [209, 68]. 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 [29, 36]. This method, which already yields interesting upper limits on the GWB, is reviewed in § 5.1. The idea of a more sensitive detector based on an array of pulsar clocks distributed over the sky is discussed in § 5.2.

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 [68] 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 _{\rm{g}}} = \frac{{243{\pi ^3}{f^4}\sigma _{\rm{g}}^2}}{{208G}}.$$
((17))

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 _{\rm{c}}} = \frac{{3H_0^2}}{{8\pi G}} \simeq 2 \times {10^{ - 29}}{h^2}{\rm{g}}\;{\rm{c}}{{\rm{m}}^{ - 3}},$$
((18))

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

Romani & Taylor [205] applied this technique to a set of TOAs for PSR B1237+12 obtained from regular observations over a period of 11 years as part of the JPL pulsar timing programme [73]. 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.3). 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.3). In the intervening period, more elaborate techniques had been devised [29, 36, 228] 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 [29] 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, P(f), can be written [99, 36] as

$$\mathcal{P}(f) = \frac{{G_{\rho _g } }} {{3\pi ^3 f^4 }} = \frac{{H_0^2 \Omega _g }} {{8\pi ^4 f^5 }} = 1.34 \times 10^4 \Omega _g h^2 f_{{\text{yr}}^{ - 1} }^{ - 5} \mu {\text{s}}^2 {\text{ yr}},$$
((19))

where \(f_{{\rm{yr}}}^{ - 1}\) is frequency in cycles per year. The timing residuals for B1937+21 shown in Fig. 25 are clearly non-white and, as we saw in § 4.3, 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 [249] place a 95% confidence upper limit of Ωgh2<10-8 for f=4.4×10-9 Hz. This limit is diffult to reconcile with most cosmic string models for galaxy formation [40, 249].

Figure 25
figure 25

Timing residuals for PSRs B1855+09 (left panel) and B1937+21 (right panel) obtained from almost a decade of timing at Arecibo (Kaspi, Taylor & Ryba 1994 [ 117]).

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 [29]. The most stringent results presently available are based on the B1855+09 limit Ωgh2<2:7×10-4 in the frequency range 10-11<f<4.4×10-9 Hz. Kopeikin [118] has recently presented this limit and discusses the methods in detail.

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 [96]. 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 ith pulsar in the array:

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

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 gravitational wave of dimensionless amplitude A. The timing noise intrinsic to the pulsar is characterised by the function Ni. The result of a cross-correlation between pulsars i and j is then

$$\alpha _i \alpha _j \left\langle {\mathcal{A}^2 } \right\rangle + \alpha _i \left\langle {\mathcal{A}\mathcal{N}_j } \right\rangle + \alpha _j \left\langle {\mathcal{A}\mathcal{N}_i } \right\rangle + \left\langle {\mathcal{N}_i \mathcal{N}_j } \right\rangle ,$$
((21))

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 αiαjA2〉 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 [95]. This allows the separation of the effects of terrestrial clock and solar system ephemeris errors from the GWB [79].

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. 26, 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 has resulted in essentially such a distribution, shown in the right panel of Fig. 26.

Figure 26
figure 26

Hammer-Aitoff projections showing the known Galactic disk millisecond pulsar population in 1990 and 2001. The impact of the new discoveries is seen by comparing the sample circa 1990, where all the known sources had been discovered at Arecibo, with the present sample where the sources are much more uniformly distributed on the sky.

A number of long-term timing projects are now underway to monitor these millisecond pulsars with a goal of detecting low-frequency gravitational radiation. At Arecibo, regular timing of a dozen or more millisecond pulsars has been carried out following the completion of the upgrade to the telescope in 1997. A summary of these observations is shown in Fig. 27. The rms timing residuals for several of the pulsars are now approaching the 100 ns level. This degree of precision demands a high level of commitment to investigate possible causes of systematic errors in the signal path through the telescope. Combining datasets from several observatories is also challenging. The Berkeley pulsar group lead by Don Backer [10] are among the most active observers in this area. Backer and collaborators have now installed identical sets of datataking equipment at a number of radio telescopes around the world in an attempt to ensure a homogeneous set of residuals. Continued timing of these millisecond pulsars in the coming years should greatly improve the sensitivity and will perhaps allow the detection of gravitational waves, as opposed to upper limits, in the not-too-distant future.

Figure 27
figure 27

Arecibo timing residuals of millisecond pulsars monitored by groups at Berkeley and Princeton on a regular basis as part of a long-term project to detect low-frequency gravitational waves. The level of precision being achieved is well below the vertical bars surrounding each set of residuals (±5 µs).

Going further

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

Summary and Future Prospects

The main aim of this article was to review some of the many astrophysical applications provided by the present sample of binary and millisecond radio pulsars. The topics covered here, along with the bibliography and associated tables of observational parameters, should be useful to those wishing to delve deeper into the vast body of literature that exists. We now briefly recap on the main issues.

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 relativistic binaries at cosmic distances (§ 3.4). Continued improvements in gravitational wave detector sensitivities 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.5) 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 in these areas and perhaps detect the gravitational wave background (if it exists) within the next decade.

These applications will benefit greatly from the continued discovery of new systems by the present generation of radio pulsar searches which 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. It is therefore likely that we have not seen all of the pulsar zoo. More sensitive surveys are being planned both in the short term (a multibeam system on the Arecibo telescope [174]) and in the longer term (the Square Kilometer Array [242]). These should provide a far more complete census of the Galactic pulsar population. Possible discoveries in the future include:

  • A dual-line binary pulsar, i.e. a double neutron star system in which both components are observable as radio pulsars. The additional clock in such a binary system would be most valuable in further tests of strong-field gravity.

  • A radio pulsar with a black-hole companion would undoubtably also be a fantastic laboratory for studying gravity in the strong-field regime.

  • A sub-millisecond pulsar. The original millisecond pulsar, B1937+21, rotating at 642 Hz is still the most rapidly rotating neutron star known. Do kHz neutron stars exist? Searches now have sensitivity to such objects [39] and a discovery of even one would constrain the equation of state of matter at high densities.

  • A binary system in which the neutron star is in the process of transforming from an X-ray-emitting neutron star to a millisecond radio pulsar.

Spurred on by recent discoveries [47, 154, 199, 62], a number of high-sensitivity searches for pulsars in globular clusters are being conducted. These have tremendous potential for discovering new and exotic binary systems like a millisecond pulsar-black hole binary.

Acknowledgments

Many thanks to Maura McLaughlin, Jiannis Seiradakis and Michael Kramer who read and commented on earlier incantations of this revised review, as well as a number of other colleagues who gave me useful feedback on the original article. Jiannis Seiradakis urged me to include the tables of parameters given in the appendix after I promised to put them in the original article, but didn’t. The tables and references should be useful to both observers and theorists. Thanks also to Fernando Camilo for allowing me to include details on a number of pulsars in these tables prior to publication.

I am indebted to a number of colleagues who kindly gave permission to use a selection of figures in this article. Michael Kramer provided the cute animation of the rotating neutron star presented in Fig. 2. Vicky Kalogera provided the graph used in Fig. 15. Joe Taylor and Joel Weisberg provided the updated orbital decay curve of PSR B1913+16 shown in Fig. 23. Fig. 20, based on unpublished timing observations carried out at Jodrell Bank, was supplied by Andrew Lyne. Vicky Kaspi provided the millisecond pulsar timing residuals and comparison of their timing stabilities shown in Figs. 25 and 21 respectively. Andrea Lommen provided the pulsar timing residuals from the Arecibo timing program used to produce Fig. 27. Frequent use was made of NASA’s magnificent Astrophysics Data System [172] and LANL’s preprint archives [137] during the literature searches. Finally, I’d like to thank the Living Reviews editor, Theresa Velden, for being extraordinarily patient, and for putting up with many a feeble excuse from me during the writing of this long-overdue update.

Figure 15
figure 15

Small-number bias of the scale factor estimates derived from a synthetic population of sources where the true number of sources is known. Left: An edge-on view of a model Galactic source population. Right: The thick line shows NG, the true number of objects in the model Galaxy, plotted against Nobs, the number detected by a flux-limited survey. The thin solid line shows Nest, the median sum of the scale factors, as a function of Nobs from a large number of Monte-Carlo trials. Dashed lines show 25 and 75% percentiles of the Nest distribution.

Appendix: Tables of Binary and Millisecond Pulsars

Table 2 Parameters for the 9 isolated millisecond pulsars currently known in the Galactic disk. Listed are the spin period P, the base-10 logarithms of the characteric age τc and surface magnetic field strength B (§ 2.4.1), the distance d derived from the Taylor & Cordes electron density model [235] or independently (when available), and the transverse speed vt inferred from d and a proper motion measurement (when available). Key publications for each pulsar are referenced to the bibliography.
Table 3 Parameters for the 9 high-eccentricity (e>0.15) binary pulsars currently known in the Galactic disk. Listed are the spin period P, the base-10 logarithms of the characteric age τc and surface magnetic field strength B (§ 2.4.1), the distance d derived from the Taylor & Cordes electron density model [235] or independently (when available), the transverse speed vt inferred from d and a proper motion measurement (when available), the binary period Pb, the projected semi-major axis of the orbit x in units of light seconds, the orbital eccentricity e, and the companion mass m2 evaluated from the mass function assuming a pulsar mass of 1.4M and an inclination angle of 60 degrees (§ 4.4) or (when known) from independent measurements. Key publications for each pulsar are referenced to the bibliography.
Table 4 Parameters for 38 low-eccentricity binary pulsars currently known in the Galactic disk. Listed are the spin period P, the base-10 logarithms of the characteric age τc and surface magnetic field strength B (§ 2.4.1), the distance d derived from the Taylor & Cordes electron density model [235] or independently (when available), the transverse speed vt inferred from d and a proper motion measurement (when available), the binary period Pb, the projected semi-major axis of the orbit x in units of light seconds, the orbital eccentricity e (when measured), and the companion mass m2 evaluated from the mass function assuming a pulsar mass of 1.4M and an inclination angle of 60 degrees (§ 4.4) or (when known) from independent measurements. Key publications for each pulsar are referenced to the bibliography.
Table 5 Parameters for the 47 pulsars currently known in globular clusters. Listed are the spin period P, cluster name and the distance d to the cluster. For binary pulsars we also list the binary period Pb, the projected semi-major axis of the orbit x in units of light seconds, the orbital eccentricity e, and the companion mass m2 evaluated from the mass function assuming a pulsar mass of 1.4M and an inclination angle of 60 degrees (§ 4.4) or (when known) from independent measurements. Key publications for each pulsar are referenced to the bibliography.