1 Introduction and Overview

Pulsars — rapidly rotating highly magnetised neutron stars — have resulted in many applications in physics and astronomy. Striking examples include the confirmation of the existence of gravitational radiation [314] as predicted by general relativity [312, 313], the first detection of an extra-solar planetary system [346, 244] and the discovery of the first double-pulsar binary system [44, 198]. The diverse zoo of radio pulsars currently known is summarized in Figure 1.

Figure 1
figure 1

Venn diagram showing the numbers and locations of the various types of radio pulsars known as of January 2005. The large and small Magellanic clouds are denoted by LMC and SMC.

Pulsar research has proceeded at a rapid pace since the first two versions of this article [179, 180]. Surveys mostly with the Parkes radio telescope [259], but also at Green Bank [232], Arecibo [228] and the Giant Metre Wave Radio Telescope [231] have more than doubled the number of pulsars known back in 1997. The most exciting new results and discoveries from these searches are discussed in this updated review.

We begin in Section 2 with an overview of the pulsar phenomenon, the key observed population properties, the origin and evolution of pulsars and the main search strategies. In Section 3, we review present understanding in pulsar demography, discussing selection effects and their correction techniques. This leads to empirical estimates of the total number of normal and millisecond pulsars (see Section 3.3) and relativistic binaries (see Section 3.4) in the Galaxy and has implications for the detection of gravitational radiation from coalescing neutron star binaries which these systems are the progenitors of. Our review of pulsar timing in Section 4 covers the basic techniques (see Section 4.2), timing stability (see Section 4.3), binary pulsars (see Section 4.4), and using pulsars as sensitive detectors of long-period gravitational waves (see Section 4.5). We conclude with a brief outlook to the future in Section 5. Up-to-date tables of parameters of binary and millisecond pulsars are included in Appendix A.

2 Pulsar Phenomenology

Many of the basic observational facts about radio pulsars were established shortly after their discovery [119] in 1967. Although there are still many open questions, the basic model has long been established beyond all reasonable doubt, i.e. 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 observational properties that are most relevant to this review.

2.1 The lighthouse model

Figure 2 shows an animation depicting the rotating neutron star or “lighthouse” model. As the neutron star spins, charged particles are accelerated out along magnetic field lines in the magnetosphere (depicted by the light blue cones). The accelerating particles 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.

Figure 2
figure 2

gif-Movie (861 KB) Still from a GIF movie showing the rotating neutron star (or “lighthouse”) model for pulsar emission. Animation designed by Michael Kramer. (For video see appendix)

Neutron stars are essentially large celestial flywheels with moments of inertia ∼ 1038 kg m2. The rotating neutron star model [242, 105] predicts a gradual slowdown and hence an 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 [270], which implied that a rotating neutron star with a large magnetic field must be the dominant energy supply for the nebula [106].

2.2 Pulse periods and slowdown rates

The present public-domain catalogue [208], available on-line [13], contains up-to-date parameters for over 1500 pulsars. Most of these are “normal” with pulse periods P ∼ 0.5 s which increase secularly at rates ∼ 10−15 s/s. A growing fraction are “millisecond pulsars”, with 1.5 ms ≲ P ≲≲ 30 ms and ≲ 10−19 s/s. For many years, the most rapidly rotating neutron star known was the original millisecond pulsar B1937+21 [19], with P = 1.5578 ms. However, one of the recent discoveries in the globular cluster Terzan 5 appears to be rotating even more rapidly, with P = 1.4 ms [233]. Confirmation of this exciting result should be published in due course. While the hunt for “sub-millisecond pulsars” continues, and most neutron star equations of state allow shorter periods, it has been suggested [36, 60] that the lack of pulsars with P < 1.5 ms is caused by gravitational wave emission from R-mode instabilities [7].

As can be seen from the “P diagram” in Figure 3, normal and millisecond pulsars are distinct populations. The differences in P and imply fundamentally different magnetic field strengths and ages. Treating the pulsar as a rotating magnetic dipole, one may show [185] that the surface magnetic field strength B ∝ (PṖ)1/2 and the characteristic age τc = P/(2).

Figure 3
figure 3

The P−Ṗ diagram showing the current sample of radio pulsars. Binary pulsars are highlighted by open circles. Theoretical models [64] do not predict radio emission outside the dark blue region. Figure provided by Michael Kramer.

Lines of constant B and τc are drawn on Figure 3, from which we infer typical values of 1012 G and 107 yr for the normal pulsars and 108 G and 109 yr for the millisecond pulsars. For the rate of loss of kinetic energy, sometimes called the spin-down luminosity, we have \(\dot E \propto \dot P/{P^3}\). The lines of constant \({\dot E}\) shown on Figure 3 show that the most energetic objects are the very young normal pulsars and the most rapidly spinning millisecond pulsars.

2.3 Pulse profiles

Pulsars are weak radio sources. Measured intensities, usually quoted in the literature for a radio frequency of 400 MHz, vary between 0.1 mJy and 5 Jy (1 Jy ≡ 10−26 W m−2 Hz−1). As a result, even with a large radio telescope, the coherent addition of many hundreds or even thousands of pulses is usually required in order to produce a detectable “integrated profile”. Remarkably, although the individual pulses vary dramatically, the integrated profile at any particular observing frequency is very stable and can be thought of as a fingerprint of the neutron star’s emission beam. Profile stability is of key importance in pulsar timing measurements discussed in Section 4.

The selection of integrated profiles in Figure 4 shows a rich diversity in morphology including 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 [209]). Since this is an unlikely viewing angle we would expect interpulses to be a rare phenomenon. Indeed, the fraction of known pulsars in which interpulses are observed in their pulse profiles is only a few percent [159].

Figure 4
figure 4

A variety of integrated pulse profiles taken from the available literature. References: Panels a, b, d, f [107], Panel c [23], Panels e, g, i [165], Panel h [29]. Each profile represents 360 degrees of rotational phase. These profiles are freely available from an on-line database [327].

Two contrasting phenomenological models to explain the observed pulse shapes are shown in Figure 5. The “core and cone” model [261] depicts the beam as a core surrounded by a series of nested cones. Alternatively, the “patchy beam” model [202, 109] has the beam populated by a series of randomly-distributed emitting regions. Further work in this area is necessary to improve our understanding of the shape and evolution of pulsar beams and fraction of sky they cover. This is of key importance to the results of population studies reviewed in Section 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 Michael Kramer.

2.4 Interstellar dispersion and the pulsar distance scale

From the sky distribution shown in Figure 6 it is immediately apparent that pulsars are strongly concentrated along the Galactic plane. This indicates that pulsars populate the disk of our Galaxy. Quantitative estimates of the distance to each pulsar can be made from the measurement of 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 lower radio frequencies travel slower through the interstellar medium, arriving later than those emitted at higher frequencies. The delay Δt in arrival times between a high frequency νhi and a low frequency νlo pulse is given [185] by

$$\Delta t = 4.15\;{\rm{ms}} \times \left[ {{{\left({{{{\nu _{{\rm{lo}}}}} \over {{\rm{GHz}}}}} \right)}^{- 2}} - {{\left({{{{\nu _{{\rm{hi}}}}} \over {{\rm{GHz}}}}} \right)}^{- 2}}} \right] \times \left({{{{\rm{DM}}} \over {{\rm{c}}{{\rm{m}}^{- 3}}{\rm{pc}}}}} \right),$$
(1)

where the dispersion measure is

$${\rm{DM =}}\int\nolimits_0^d {{n_{\rm{e}}}dl}$$
(2)

is the integrated free electron column density ne out to the pulsar at a distance d. This equation may be solved for d given a measurement of DM and a model of the free electron distribution.

Figure 6
figure 6

The sky distribution of pulsars in Galactic coordinates. The plane of the Galaxy is the central horizontal line. The Galactic centre is the midpoint of this line. Millisecond pulsars are indicated in red. Binary pulsars are highlighted by the open circles. The more isotropic distribution of the millisecond and binary pulsars reflects the differences in detecting them out to large distances cf. the normal population (see Section 3 ).

In practice, the electron density model is calibrated from the 100 or so pulsars with independent distance estimates and measurements of scattering for lines of sight towards various Galactic and extragalactic sources [310, 335]. The most recent model of this kind [73, 74] provides distance estimates with an average uncertainty of ∼ 30%.

2.5 Pulsars in binary systems

As can be inferred from Figure 1, only 4% of all known pulsars in the Galactic disk are members of binary systems. Timing measurements (see Section 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” B1257+12 and B1620−26. PSR B1257+12 is a 6.2-ms pulsar accompanied by at least three terrestrial-mass bodies [346, 244, 345] while B1620−26, an 11-ms pulsar in the globular cluster M4, is part of a triple system with a 1–2 MJupiter planet [317, 17, 315, 284] orbiting a neutron star-white dwarf stellar binary system. The current limits from pulsar timing favour a roughly 45-yr low-eccentricity (e ∼ 0.16) orbit with semi-major axis ∼ 25 AU. Despite several tentative claims over the years, no other convincing cases for planetary companions to pulsars exist. Orbiting companions are much more common around millisecond pulsars (∼ 80% of the observed sample) than around the normal pulsars (≲ 1%). Binary systems below the line have low-mass companions (≲ 0.7 M — predominantly white dwarfs) and essentially circular orbits: 10−5e ≲ 0.01. Binary pulsars with high-mass companions (≳ 1 M — massive white dwarfs, other neutron stars or main sequence stars) tend to have more eccentric orbits, 0.15 ≳ e ≳ 0.9.

2.6 Evolution of normal and millisecond pulsars

A simplified version of the presently favoured model [37, 95, 285, 2] to explain the formation of the various types of systems observed is shown in Figure 7. Starting with a binary star system, a neutron star is formed during the supernova explosion of the initially more massive star. From the virial theorem, in the absence of any other factors, the binary will be disrupted if more than half the total pre-supernova mass is ejected from the system during the (assumed symmetric) explosion [120, 33]. In practice, the fraction of surviving binaries is also affected by the magnitude and direction of any impulsive “kick” velocity the neutron star receives at birth from a slightly asymmetric explosion [120, 20]. Binaries that disrupt produce a high-velocity isolated neutron star and an OB runaway star [38]. The high probability of disruption explains qualitatively why so few normal pulsars have companions. There are currently four known normal radio pulsars with massive main sequence companions in eccentric orbits which are examples of binary systems which survived the supernova explosion [140, 153, 295, 296, 193]. Over the next 107–8 yr 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 diminishes to a point where it no longer produces significant amounts of radio emission.

Figure 7
figure 7

Cartoon showing various evolutionary scenarios involving binary pulsars.

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 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. Two classes of X-ray binaries relevant to binary and millisecond pulsars exist: neutron stars with high-mass or low-mass companions. Detailed reviews of the X-ray binary population, including systems likely to contain black holes, can be found elsewhere [33].

In a high-mass X-ray binary, the companion is massive enough that it might also explode as a supernova, producing a second neutron star. If the binary system is lucky enough to survive the explosion, the result is a double neutron star binary. At least five and probably eight such systems are now known, the original example being PSR B1913+16 [129] — a 59-ms radio pulsar which orbits its companion every 7.75 hr [312, 313]. 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.

For many years, no clear example was known where the second-born neutron star was observed as a pulsar. The discovery of the double pulsar J0737-3039 [44, 198], where a 22.7-ms recycled pulsar “A” orbits a 2.77-s normal pulsar “B” every 2.4 hr, has now provided a dramatic confirmation of this evolutionary model in which we identify A and B as the first and second-born neutron stars respectively. Just how many more observable double pulsar systems exist in our Galaxy is not clear. Although the population of double neutron star systems in general is reasonably well understood (see Section 3.4.1), given that the lifetime of the second born pulsar is less than one tenth that of the recycled pulsar, and that its radio beam is likely to be much smaller (see Section 3.2.3), the prospects of ever finding more than a few 0737-like systems are rather low.

The companion in a low-mass X-ray binary evolves and transfers matter onto the neutron star on a much longer time-scale, spinning it up to periods as short as a few ms [2]. Tidal forces during the accretion process serve to circularize the orbit. 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. 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 [338], as well as Doppler-shifted 2.49-ms X-ray pulsations from the transient X-ray burster SAX J1808.4-3658 [340, 59]. Six other “X-ray millisecond pulsars” are now known with spin rates and orbital periods ranging between 185–600 Hz and 40 min–4.3 hr respectively [339, 220, 144]. Despite intensive searches [43], no radio pulsations have so far been detected in these binaries. This could be a result of free-free absorption of any radio waves by the thick accretion disk, or perhaps quenching the accelerating potential in the neutron star magnetosphere by infalling matter. More sensitive radio observations are ultimately required to place stringent limits on the physical processes responsible for the lack of emission.

Numerous examples of these systems in their post X-ray phase are now seen as the millisecond pulsar-white dwarf binary systems. Presently, 20 of these systems have compelling optical identifications of the white dwarf companion, and upper limits or tentative detections have been found in about 30 others [330]. Comparisons between the cooling ages of the white dwarfs and the millisecond pulsars confirm the age of these systems and suggest that the accretion rate during the spin-up phase was well below the Eddington limit [114].

Further support for the above evolutionary scenarios comes from two correlations in the observed sample of low-mass binary pulsars. Firstly, as seen in Figure 8, there is a strong correlation between orbital period and eccentricity. The data are very good agreement with a theoretical relationship which predicts a relic orbital eccentricity due to convective eddy currents in the accretion process [249]. Secondly, as shown in Panel b of Figure 9, where companion masses have been measured accurately, through radio timing (see Section 4.4) and/or through optical observations [330], they are in good agreement with a relation between companion mass and orbital period predicted by binary evolution theory [307]. A word of caution is required in using these models to make predictions, however. When confronted with a larger ensemble of binary pulsars using statistical arguments to constrain the companion masses (see Panel a of Figure 9), current models have problems in explaining the full range of orbital periods on this diagram [294].

Figure 8
figure 8

Eccentricity versus orbital period for a sample of 21 low-mass binary pulsars which are not in globular clusters, with the triangles denoting three recently discovered systems [294]. The solid line shows the median of the predicted relationship between orbital period and eccentricity [249]. Dashed lines show 95% the confidence limit about this relationship. The dotted line shows Pbe2. Figure provided by Ingrid Stairs [294] using an adaptation of the orbital period-eccentricity relationship tabulated by Fernando Camilo.

Figure 9
figure 9

Companion mass versus orbital period for binary pulsars showing the whole sample where, in the absence of mass determinations, statistical arguments based on a random distribution of orbital inclination angles (see Section 4.4 ) have been used to constrain the masses as shown (Panel a), and only those with well determined companion masses (Panel b). The dashed lines show the uncertainties in the predicted relation [307]. This relationship indicates that as these systems finished a period of stable mass transfer due to Roche-lobe overflow, the size and hence period of the orbit was determined by the mass of the evolved secondary star. Figure provided by Marten van Kerkwijk [330].

The range of white dwarf masses observed is becoming broader. Since this article originally appeared [179], the number of “intermediate-mass binary pulsars” [49] has grown significantly [54]. These systems are distinct to the millisecond pulsar-white dwarf binaries in several ways:

  1. 1.

    The spin period of the radio pulsar is generally longer (9–200 ms).

  2. 2.

    The mass of the white dwarf is larger (typically ≳ 0.5 M).

  3. 3.

    The orbit, while still essentially circular, is often significantly more eccentric (e ≳ 10−3).

  4. 4.

    They do not necessarily follow the mass-period or eccentricity-period relationships.

It is not presently clear whether these systems originated from low- or high-mass X-ray binaries. It was suggested by van den Heuvel [329] that they have more in common with high-mass systems. More recently, it has been proposed [172] that a thermal-viscous instability in the accretion disk of a low-mass X-ray binary could truncate the accretion phase and produce a more slowly spinning neutron star.

2.7 Isolated recycled pulsars

The scenarios outlined qualitatively above represent a reasonable understanding of binary evolution. There are, however, a number of pulsars with spin properties that suggest a phase of recycling took place but have no orbiting companions. While the existence of such systems in globular clusters are more readily explained by the high probability of stellar interactions compared to the disk [283], it is somewhat surprising to find them in the Galactic disk. Out of a total of 66 millisecond pulsars in the Galactic disk, 15 are isolated (see Table 2). Although it has been proposed that these millisecond pulsars have ablated their companion via their strong relativistic winds [158] as may be happening in the PSR B1957+20 system [103], it is not clear whether the energetics or time-scales for this process are feasible [170]. There is some observational evidence that suggests that solitary millisecond pulsars are less luminous than binary millisecond pulsars [23, 165]. If confirmed by future discoveries, this would need to be explained by any viable evolutionary model.

Table 2 Parameters for the 15 isolated millisecond pulsars currently known in the Galactic disk. Two of the pulsars (J0609+2130 and J2235+1506) are thought to be the descendants of high-mass binary systems. Listed are the spin period P, the base-10 logarithms of the characteristic age τc and surface magnetic field strength B (see Section 2.2 ), the distance d derived from the Cordes & Lasio electron density model [73, 74] 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.

There are two further “anomalous” isolated pulsars with periods in the range 55–60 ms [55, 188]. When placed on the P diagram, these objects populate the region occupied by the double neutron star binaries. The most natural explanation for their existence, therefore, is that they are “failed double neutron star binaries” which disrupted during the supernova explosion of the secondary [55]. A simple calculation [188], however, suggests that for every double neutron star we should see of order ten such isolated objects. Although some other examples of such pulsars are known [150, 91, 42], exactly why so few are observed is currently not clear.

2.8 Pulsar velocities

Pulsars have long been known to have space velocities at least an order of magnitude larger than those of their main sequence progenitors, which have typical values between 10 and 50 km s−1. The first direct evidence for large velocities came from optical observations of the Crab pulsar [321], showing that the neutron star has a velocity in excess of 100 km s−1. Proper motions for 233 pulsars have subsequently been measured largely by radio timing and interferometric techniques [194, 24, 96, 115, 122]. These data imply a broad velocity spectrum ranging from 0 to over 1000 km s−1 [201].

Such large velocities are perhaps not surprising, given the violent conditions under which neutron stars are formed. Shklovskii [282] demonstrated that, if the explosion is only slightly asymmetric, an impulsive “kick” velocity of up to 1000 km s−1 can be 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. As Figure 10 illustrates, high-velocity pulsars born close to the Galactic plane quickly migrate to higher Galactic latitudes. Given such a broad velocity spectrum, as many as half of all pulsars will eventually escape the Galactic gravitational potential [201, 71].

Figure 10
figure 10

gif-Movie (213 KB) Still from a GIF movie showing 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. (For video see appendix)

The distribution of pulsar velocities remains an issue of contention. For example, Monte Carlo simulations by Arzoumanian et al. [8] strongly favour a bimodal distribution with low and high velocity peaks. On the other hand, a recent study [122] found the mean birth velocity of normal pulsars to be consistent with a Maxwellian distribution with a mean of ∼ 400 km s−1. Regardless of the form of the distribution, however, we can say that the mean velocity of young pulsars is significantly larger than the 80–140 km s−1 range for millisecond and binary pulsars [178, 70, 205]. The main reasons for the lower velocities are the fact that the kick must have been small to avoid disruption, and the surviving neutron star has to pull the companion along with it, thus slowing the system down.

2.9 Current and future binary and millisecond pulsar search strategies

2.9.1 All-sky searches

The oldest radio pulsars form a relaxed population of stars oscillating in the Galactic gravitational potential [113]. 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 began with the discovery of two recycled pulsars at high latitudes with the 305-m Arecibo telescope: the double neutron star binary B1534+12 [343] and the “planets pulsar” B1257+12 [346]. Surveys carried out at Arecibo, Parkes, Jodrell Bank and Green Bank (using the 140 ft telescope) by others in the 1990s (summarised in three review papers [47, 50, 51]) have found many other millisecond and recycled pulsars in this way.

2.9.2 Searches close to the plane of our Galaxy

Young pulsars are most likely to be found near to their places of birth, before they have had time for their velocity to move them away, and hence they lie close to the Galactic plane (see Figure 6). This was the target region of the main Parkes multibeam survey and has already resulted in the discovery of over 700 new pulsars [210, 221, 163, 123, 94], 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.8-hr orbit around a white dwarf [154, 241, 25, 127], PSR J1740-3052, a young pulsar orbiting an ∼ 11 M star (probably a main sequence B-star [295, 296]), PSR 1638-4715, a young pulsar in a ∼ 5 yr-eccentric orbit (e ∼ 0.9) around around a 10–20 M companion [193], several intermediate-mass binary pulsars [54], and two double neutron star binaries [199, 93]. Although the main survey has now been completed, extensions of the survey region, and re-analyses of existing data [94] will ensure further discoveries in the near future.

2.9.3 Searches at intermediate Galactic latitudes

To probe more deeply into the population of millisecond and recycled pulsars than possible at high Galactic latitudes, the Parkes multibeam system was also used to survey intermediate latitudes [92, 90]. Among the 69 new pulsars found in the survey, 8 are relatively distant recycled objects. Two of the new recycled pulsars from this survey [90] 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. Arecibo surveys at intermediate latitudes also continue to find new pulsars, such as the long-period binaries J2016+1948 and J0407+1607 [234, 190], and the likely double neutron star system J1829+2456 [62].

2.9.4 Targeted searches of globular clusters

Globular clusters have long been known to be breeding grounds for millisecond and binary pulsars [58]. The main reason for this is the high stellar density and consequently high rate of stellar interaction in globular clusters relative to most of the rest of the Galaxy. As a result, low-mass 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 [283]. To date, searches have revealed 103 pulsars in 24 globular clusters (see Table 5 and [230, 58]). Early highlights include the double neutron star binary in M15 [257] and a low-mass binary system with a 95-min orbital period in 47 Tucanae [53], one of 22 millisecond pulsars currently known in this cluster alone [53, 183].

On-going surveys of clusters continue to yield new surprises [264, 77] including the discovery of the most eccentric binary pulsar so far — J0514-4002 is a 4.99 ms pulsar in a highly eccentric (e = 0.89) binary system in the globular cluster NGC 1851 [100]. The latest sensation, however, is the discovery [266, 233] of 28 millisecond pulsars in Terzan 5 with the Green Bank Telescope [232]. This brings the total known in this cluster to 26. The spin periods and orbital parameters of the new pulsars reveal that, as a population, they are significantly different to the pulsars of 47 Tucanae which have periods in the range 2–8 ms [183]. The spin periods of the new pulsars span a much broader range (1.67–80 ms) including the second and third shortest spin periods of all pulsars currently known. The binary pulsars include two systems with eccentric orbits and likely white dwarf companions. No such systems are known in 47 Tucanae. The difference between the two pulsar populations may reflect the different evolutionary states and physical conditions of the two clusters. In particular, the central stellar density of Terzan 5 is about twice that of 47 Tucanae, suggesting that the increased rate of stellar interactions might disrupt the recycling process for the neutron stars in some binary systems and induce larger eccentricities in others.

2.9.5 Current and future surveys

Motivated by the successes at Parkes, a multibeam receiver has recently been installed at the Arecibo telescope [227]. A preliminary survey using this instrument began in 2004 and has so far discovered 10 pulsars [229]. The main virtue of this survey is the shorter integration times employed and resulting higher sensitivity to accelerated systems than the Parkes survey. Depending on the stamina of the observers, up to 1000 pulsars could be found in this survey over the next decade. Other surveys with the Green Bank Telescope [232] and with the Giant Metre Wave Radio Telescope [231] are also expected to make major contributions.

All surveys that have so far been conducted, or will be carried out in the next few years, are likely to be surpassed by the Square Kilometre Array [132] — the next generation radio telescope, planned to come online around 2020. Simulations suggest [162] that at least 104 pulsars, including 103 millisecond pulsars, could be detected in our Galaxy.

2.10 Going further

In this brief review of pulsar astronomy, we have only skimmed the surface of many topics which are covered more thoroughly elsewhere. Two recent textbooks, Pulsar Astronomy [200] and Handbook of Pulsar Astronomy [185], cover the literature and techniques and provide excellent further reading. The morphological properties of pulsars have recently been comprehensively discussed in two recent reviews [286, 280]. Those wishing to approach the subject from a more theoretical viewpoint are advised to read The Theory of Neutron Star Magnetospheres [218] and The Physics of the Pulsar Magnetosphere [32]. Our summary of evolutionary aspects serves merely as a primer to the vast body of literature available. Further insights can be found from more detailed reviews [33, 251, 292].

Pulsar resources available on the Internet are continually becoming more extensive and useful. A good starting point for pulsar-surfers is the Handbook of Pulsar Astronomy website [186], as well as the pages maintained at Arecibo [226], Berkeley [325], Bonn [215], Cagliari [324], Jodrell Bank [326], Princeton [258], Swinburne [305], UBC [323] and Sydney [259].

3 Pulsar Statistics and Demography

The observed pulsar sample is heavily biased towards the brighter objects that are the easiest to detect. What we observe represents only the tip of the iceberg of a much larger underlying population [108]. 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 Figure 11. The clustering of sources around the Sun seen in the left panel of Figure 11 is clearly at variance with the distribution of other stellar populations which show a radial distribution about the Galactic centre.

Figure 11
figure 11

Left panel: The current sample of all known radio pulsars projected onto the Galactic plane. The Galactic centre is at the origin and the Sun is at (0, 8.5) kpc. Note the spiral-arm structure seen in the distribution which is now required by the electron density model [73, 74]. Right panel: Cumulative number of pulsars as a function of projected distance from the Sun. The solid line shows the observed sample while the dotted line shows a model population free from selection effects.

The extent to which the pulsar sample is incomplete can be seen from the cumulative number of pulsars as a function of the projected distance from the Sun. In Figure 11 the observed distribution is compared to the expected distribution for a simple model population with no selection effects. The observed number distribution becomes strongly deficient beyond a few kpc.

3.1 Selection effects in pulsar searches

3.1.1 The inverse square law and survey thresholds

The most prominent selection effect is the inverse square law, i.e. for a given intrinsic luminosityFootnote 1, 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 high luminosity objects. Beyond distances of a few kpc from the Sun, the apparent flux density falls below the detection thresholds Smin of most surveys. Following [88], we express this threshold as follows:

$${S_{\min}} = {{{\rm{S}}/{{\rm{N}}_{\min}}} \over {\eta \sqrt {{n_{{\rm{pol}}}}}}}\left({{{{T_{{\rm{rec}}}} + {T_{{\rm{sky}}}}} \over {\rm{K}}}} \right){\left({{G \over {{\rm{K}}\;{\rm{J}}{{\rm{y}}^{- 1}}}}} \right)^{- 1}}{\left({{{\Delta \nu} \over {{\rm{MHz}}}}} \right)^{- 1/2}}{\left({{{{t_{{\mathop{\rm int}}}}} \over s}} \right)^{- 1/2}}{\left({{W \over {P - W}}} \right)^{1/2}}{\rm{mJy}},$$
(3)

where S/Nmin is the threshold signal-to-noise ratio, η is a generic fudge factor (≲ 1) which accounts for losses in sensitivity (e.g., due to sampling and digitization noise), npol is the number of polarizations recorded (either 1 or 2), Trec and Tsky are the receiver and sky noise temperatures, G is the gain of the antenna, Av is the observing bandwidth, tint is the integration time, W is the detected pulse width and P is the pulse period.

3.1.2 Interstellar pulse dispersion and multipath scattering

It follows from Equation (3) that the sensitivity decreases as W/(PW) and hence 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 may be broader than the intrinsic value largely as a result of pulse dispersion and multipath scattering by free electrons in the interstellar medium. 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. Multipath scattering from electron density irregularities results in a one-sided broadening due to the delay in arrival times. A simple scattering model is shown in Figure 12 in which the scattering electrons are assumed to lie in a thin screen between the pulsar and the observer [278]. The timescale of this effect varies roughly as ν−4, which can not be removed by instrumental means.

Figure 12
figure 12

Left panel: Pulse scattering caused by irregularities in the interstellar medium. The different path lengths and travel times of the scattered rays result in a “scattering tail” in the observed pulse profile which lowers its signal-to-noise ratio. Right panel: A simulation showing the percentage of Galactic pulsars that are likely to be undetectable due to scattering as a function of observing frequency. Low-frequency (≲ 1 GHz) surveys clearly miss a large percentage of the population due to this effect.

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 [67, 139] compared to the 400-MHz search frequency used in early surveys. 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 [169]. Pulsars themselves have typical spectral indices of −1.6 [191], so that flux densities are roughly an order of magnitude lower at 1400 MHz compared to 400 MHz. Fortunately, this can be at least partially 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 [199] compared to the 430-MHz system, where nominally 32 MHz is available [211].

3.1.3 Orbital acceleration

Standard pulsar searches use Fourier techniques [185] 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 binary systems with orbital periods longer than about a day. For shorter-period binary systems, 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 S/N [136]. 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) pulsar period in its rest frame and c is the speed of light. Given that the width of a frequency bin in the Fourier domain 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, Figure 13 shows a 22.5-min search mode observation of the binary pulsar B1913+16 [129, 312, 313]. Although this observation covers only about 5% of the orbit (7.75 hr), the severe effects of the Doppler smearing on the pulse signal are very apparent. While the standard search code nominally detects the pulsar with S/N = 9.5 for this observation, it is clear that this value is significantly reduced due to the Doppler shifting of the pulse period seen in the individual sub-integrations.

Figure 13
figure 13

Left panel: 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 quadratic drifting in phase of the pulse during the observation (linear grey scale plot). Right panel: 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” [219], assumes the pulsar has a constant acceleration during the observation. 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 Figure 13). This technique was first used to find PSR B2127+11C [5], a double neutron star binary in M15 which has parameters similar to B1913+16. More recently, its application to 47 Tucanae [53] resulted in the discovery of nine binary millisecond pulsars, including one in a 96-min orbit around a low-mass (0.15 M) companion. This is currently the shortest binary period for any known radio pulsar.

For intermediate orbital periods, in the range 30 min-several hours, another promising technique is the dynamic power spectrum search. Here the time series is split into a number of smaller contiguous segments which are Fourier-transformed separately. The individual spectra are displayed as a two-dimensional (frequency versus time) image. Orbitally modulated pulsar signals appear as sinusoidal signals in this plane as shown in Figure 14.

Figure 14
figure 14

Dynamic power spectra showing two recent pulsar discoveries in the globular cluster M62 showing fluctuation frequency as a function of time. Figure provided by Adam Chandler.

This technique has been used by various groups where spectra are inspected visually [206]. Much of the human intervention can be removed using a hierarchical scheme for selecting significant events [63]. This approach was recently applied to a search of the globular cluster M62 resulting in the discovery of three new pulsars. One of the new discoveries — M62F, a faint 2.3-ms pulsar in a 4.8-hr orbit — was detectable only using the dynamic power spectrum technique.

For the shortest orbital periods, the assumption of a constant acceleration during the observation clearly breaks down. In this case, a particularly efficient algorithm has been developed [262, 143, 263] which is optimised to finding binaries with periods so short that many orbits can take place during an observation. This “phase modulation” technique exploits the fact that the Fourier components are modulated by the orbit to create a family of periodic sidebands around the nominal spin frequency of the pulsar. While this technique has so far not resulted in any new discoveries, the existence of short period binaries in 47 Tucanae [53], Terzan 5 [266] and the 11-min X-ray binary X1820-303 in NGC 6624 [301], suggests that there may be more ultra-compact binary pulsars that await discovery.

3.2 Correcting the observed pulsar sample

In the following, we review common techniques to account for the various selection effects and form a less biased picture of the true pulsar population.

3.2.1 Scale factor determination

A very useful technique [250, 334], 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 is detectable:

$$\xi (P,L) = {{\int {\int\nolimits_{\rm{G}} {\Sigma (R,\,z)} RdRdz}} \over {\int {\int\nolimits_{P,L} {\Sigma (R,\,z)} RdRdz}}}.$$
(4)

Here, Σ(R, z) is the assumed pulsar space density 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.

In practice, ξ is calculated for each pulsar separately using a Monte Carlo simulation to model the volume of the Galaxy probed by the major surveys [222]. For a sample of Nobs observed pulsars above a minimum luminosity Lmin, the total number of pulsars in the Galaxy is

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

where f is the model-dependent “beaming fraction” discussed below in Section 3.2.3. Note that this estimate applies to those pulsars with luminosities ≳ Lmin. Monte Carlo simulations have shown this method to be reliable, as long as Nobs is reasonably large [182].

3.2.2 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 (5), therefore, will tend to underestimate the true size of the population. This “small-number bias” was first pointed out [145, 151] for the sample of double neutron star binaries where we know of only four systems relevant for calculations of the merging rate (see Section 3.4.1). Only when Nobs ≳ 10 does the sum of the scale factors become a good indicator of the true population size.

Despite a limited sample size, recent work [156] has demonstrated that rigorous confidence intervals of NG can be derived. Monte Carlo simulations verify that the simulated number of detected objects Ndetected closely follows a Poisson distribution and that Ndetected = αNG, where α is a constant. By varying the value of NG in the simulations, the mean of this Poisson distribution can be measured. Using a Bayesian analysis it was shown [156] that, for a single object, the probability density function of the total population is

$$P({N_{\rm{G}}}) = {\alpha^2}{N_{\rm{G}}}\exp (- \alpha {N_{\rm{G}}}).$$
(6)

Adopting the necessary assumptions required in the Monte Carlo population about the underlying pulsar distribution, this technique can be used to place interesting constraints on the size and, as we shall see later, birth rate of the underlying population.

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 panel: An edge-on view of a model Galactic source population. Right panel: The thick line shows Ng, the true number of objects in the model Galaxy, plotted against Nobserved, 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.

3.2.3 The beaming correction

The “beaming fraction” f in Equation (5) is the fraction of 4π steradians swept out by a pulsar’s radio beam during one rotation. Thus f is the probability that the beam cuts the line-of-sight of an arbitrarily positioned observer. A naïve estimate for f of roughly 20% assumes a circular beam of width ∼ 10° and a randomly distributed inclination angle between the spin and magnetic axes [311]. Observational evidence suggests that shorter period pulsars have wider beams and therefore larger beaming fractions than their long-period counterparts [223, 202, 34, 306]. As can be seen in Figure 16, however, a consensus on the beaming fraction-period relation has yet to be reached.

Figure 16
figure 16

Beaming fraction plotted against pulse period for four different beaming models: Narayan & Vivekanand 1983 (NV83) [223], Lyne & Manchester 1988 (LM88) [202], Biggs 1990 (JDB90) [34] and Tauris & Manchester 1998 (TM88) [306].

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. An analysis of a large sample of millisecond pulsar profiles [165] suggests that their beaming fraction lies between 50 and 100%. The large beaming fraction and narrow pulses often observed strongly suggests a fan beam model for millisecond pulsars [218].

3.3 The population of normal and millisecond pulsars

3.3.1 Luminosity distributions and local number estimates

Based on a number of all-sky surveys carried out in the 1990s, the scale factor approach has been used to derive the characteristics of the true normal and millisecond pulsar populations and is based on the sample of pulsars within 1.5 kpc of the Sun [205]. Within this region, the selection effects are well understood and easier to quantify than in the rest of the Galaxy. These calculations should therefore give a reliable local pulsar population estimate.

The luminosity distributions obtained from this analysis are shown in Figure 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 a beaming model [34], of 156±31pulsars kpc−2 for luminosities above 1 mJy kpc2. The same analysis for the millisecond pulsars, assuming a mean beaming fraction of 75% [165], leads to a local surface density of 38 ± 16 pulsars kpc−2 also for luminosities above 1 mJy kpc2.

Figure 17
figure 17

Left panel: 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 panel: 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.

3.3.2 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 [137, 181]. One approach is to assume that pulsars have a radial distribution similar to that of other stellar populations and to scale the local number density with this distribution in order to estimate the total Galactic population. The corresponding local-to-Galactic scaling is 1000 ± 250 kpc2 [268]. This implies a population of ∼ 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 required to sustain the observed population. From the P−Ṗ diagram in Figure 3, 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 [328]. As noted in Section 2.2, the millisecond pulsars are much older, with ages close to that of the Universe τu (we assume here τu = 13.8 Gyr [347]). Taking the maximum age of the millisecond pulsars to be τu, we infer a mean birth rate of at least 1 per 345, 000 yr. This is consistent, within the uncertainties, with the birth-rate of low-mass X-ray binaries [189].

3.4 The population of relativistic binaries

Although no radio pulsar has so far been observed in orbit around a black hole companion, 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 timescale. The current sample of objects is shown as a function of orbital period and eccentricity in Figure 18. Isochrones showing various coalescence times τg are calculated using the expression

$${\tau _{\rm{g}}} \simeq 9.83 \times {10^6}{\rm{yr}}{\left({{{{P_b}} \over {{\rm{hr}}}}} \right)^{8/3}}{\left({{{{m_1} + {m_2}} \over {{M_ \odot}}}} \right)^{- 2/3}}{\left({{\mu \over {{M_ \odot}}}} \right)^{- 1}}{(1 - {e^2})^{7/2}},$$
(7)

where m1 and m2 are the masses of the two stars, μ = m1m2/(m1 + m2) is the so-called “reduced mass”, Pb is the binary period and e is the eccentricity. This formula is a good analytic approximation (within a few percent) to the numerical solution of the exact equations for τg [245, 246].

Figure 18
figure 18

The relativistic binary merging plane. Top: Orbital eccentricity versus period for eccentric binary systems involving neutron stars. Bottom: Orbital period distribution for the massive white dwarf-pulsar binaries.

In addition to testing general relativity through observations of these systems (see Section 4.4), estimates of their Galactic population and merger rate are of great interest as one of the prime sources for current gravitational wave detectors such as GEO600 [112], LIGO [46], VIRGO [131] and TAMA [225]. In the following, we review recent empirical determination of the population sizes and merging rates of binaries where at least one component is visible as a radio pulsar.

3.4.1 Double neutron star binaries

As discussed in Section 2.2, double neutron star (DNS) binaries are expected to be rare. This is certainly the case; despite extensive searches, only five certain DNS binaries are currently known: PSRs J0737-3039 [44], B1534+12 [343], J1756-2251 [93], B1913+16 [129] and B2127+11C [257]. Although we only see both neutron stars as pulsars in J0737-3039 [198], we are “certain” of the identification in the other four systems from precise component mass measurements from pulsar timing observations (see Section 4.4). The spin and orbital parameters, timescales and mass constraints for these systems are listed in Table 1, along with three other binaries with eccentric orbits, mass functions and periastron advance measurements that are consistent with a DNS identification, but for which there is presently not sufficient component mass information to confirm their nature.

Table 1 Known and likely DNS binaries. Listed are the pulse period P, orbital period Pb, orbital eccentricity e, characteristic age τc, and expected binary coalescence time-scale τG due to gravitational wave emission calculated from Equation (7 ). To distinguish between definite and candidate DNS systems, we also list whether the masses of both components have been determined via the measurement of two or more post Keplerian parameters as described in Section 4.4 .

Despite the uncertainties in identifying DNS binaries, for the purposes of determining the Galactic merger rate, the systemsFootnote 2 for which τg is less than τu (i.e. PSRs J0737-3039, B1534+12, J1756-2251, 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 [257, 252] and is thought to make only a negligible contribution to the merger rate [248]. The general approach with the remaining systems is to derive scale factors for each object, construct the probability density function of their total population (as outlined in Section 3.2.1) and then divide these by a reasonable estimate for the lifetime. It is generally accepted [151] that the observable lifetimes for these systems are determined by the timescale on which the current orbital period is reduced by a factor of two [9]. Below this point, the orbital smearing selection effect discussed in Section 3.1.3 will render the binary undetectable by current surveys.

The results of the most recent study of this kind [147, 148], which take into account the discovery of PSR J0737-3039, are summarised in the left panel of Figure 19. The combined Galactic merger rate, dominated by the double pulsar, is found to be \(80_{- 70}^{+ 210}\) Myr−1, where the uncertainties reflect the 95% confidence level using the techniques summarised in Section 3.2.2. Extrapolating this number to include DNS binaries detectable by LIGO in other galaxies [248], the expected event rate is \(35_{- 30}^{+ 90} \times {10^{- 3}}{\rm{y}}{{\rm{r}}^{- 1}}\) for initial LIGO and \(190_{- 150}^{+ 470}{\rm{y}}{{\rm{r}}^{- 1}}\) for advanced LIGO. Future prospects for detecting gravitational wave emission from binary neutron star inspirals are therefore very encouraging. Since much of the uncertainty in the rate estimates is due to our ignorance of the underlying distribution of double neutron star systems, future detection rates could ultimately constrain the properties of this exciting binary species.

Figure 19
figure 19

The current best empirical estimates of the coalescence rates of relativistic binaries involving neutron stars. The individual contributions from each known binary system are shown as dashed lines, while the solid line shows the total probability density function on a logarithmic and (inset) linear scale. The left panel shows DNS binaries [147, 148], while the right panel shows the equivalent results for NS-WD binaries [157].

Although the double pulsar system J0737-3039 will not be important for ground-based detectors until its final coalescence in another 85 Myr, it may be a useful calibration source for the future space-based detector LISA [224]. Recent calculations [146] show that a 1-yr observation with LISA would detect (albeit with S/N ∼ 2) the continuous emission at a frequency of 0.2 mHz based on the current orbital parameters. Although there is the prospect of using LISA to detect similar systems systems through their continuous emission, current calculations [146] suggest that significant (S/N > 5) detections are not likely. Despite these limitations, it is likely that LISA observations will be able to place independent constraints on the Galactic DNS binary population after several years of operation.

3.4.2 White dwarf-neutron star binaries

Although the population of white dwarf-neutron star (WDNS) binaries in general is substantial, the fraction which will merge due to gravitational wave emission is small. Like the DNS binaries, the observed WDNS sample suffers from small-number statistics. From Figure 18, we note that only three WDNS systems are currently known that will merge within τu, PSRs J0751+1807 [192], J1757-5322 [90] and J1141-6545 [154]. Applying the same techniques as used for the DNS population, the merging rate contributions of the three systems can be calculated [157] and are shown in Figure 19. The combined Galactic coalescence rate is \(4_{- 3}^{+ 5}\) Myr−1 (at 68% confidence interval). Although the orbital frequencies of these objects at coalescence are too low to be detected by LIGO, they do fall within the band planned for LISA [224]. Unfortunately, an extrapolation of the Galactic event rate out to distances at which such events would be detectable by LISA does not suggest that these systems will be a major source of detection [157]. Similar conclusions were reached by considering the statistics of low-mass X-ray binaries [69].

3.5 Going further

Studies of pulsar population statistics represent a large proportion of the pulsar literature. Good starting points for further reading can be found in other review articles [218, 33]. Our coverage of compact object coalescence rates has concentrated on empirical methods. An alternative approach is to undertake a full-blown Monte Carlo simulation of the most likely evolutionary scenarios described in Section 2.6. In this 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. Although selection effects are not always taken into account in this approach, the final census is usually normalized to the star formation rate.

Numerous population syntheses (most often to populations of binaries where one or both members are NSs) can be found in the literature [86, 272, 322, 254, 173, 27]. A group in Moscow has made a web interface to their code [302]. 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 [148, 157, 149]. 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.

4 Principles and Applications of Pulsar Timing

It became clear soon after their discovery that pulsars are excellent celestial clocks. The period of the first pulsar [119] was found to be stable to one part in 107 over a few months. Following the discovery of the millisecond pulsar B1937+21 [19] it was demonstrated that its period could be measured to one part in 1013 or better [84]. This unrivaled stability leads to a host of applications including uses as time keepers, probes of relativistic gravity and natural gravitational wave detectors.

4.1 Observing basics

Each pulsar is typically observed at least once or twice per month over the course of a year to establish its basic properties. Figure 20 summarises the essential steps involved in a “time-of-arrival” (TOA) measurement. Pulses from the neutron star traverse the interstellar medium before being received at the radio telescope where they are de-dispersed and added to form a mean pulse profile.

Figure 20
figure 20

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; see [83]). The TOA is defined as the arrival time of some fiducial point on the integrated profile with respect to either the start or the midpoint of the observation. Since the profile has a stable form at any given observing frequency (see Section 2.3), the TOA can be accurately determined by cross-correlation of the observed profile with a high S/N “template” profile obtained from the addition of many observations at the particular observing frequency.

Successful pulsar timing requires optimal TOA precision which largely depends on the signal-to-noise ratio (S/N) of the pulse profile. Since the TOA uncertainty ϵTOA is roughly the pulse width divided by the S/N, using Equation (3) we may write the fractional error as

$${{{\epsilon _{{\rm{TOA}}}}} \over P} \simeq {\left({{{{S_{{\rm{psr}}}}} \over {{\rm{mJy}}}}} \right)^{- 1}}\left({{{{T_{{\rm{rec}}}} + {T_{{\rm{sky}}}}} \over {\rm{K}}}} \right){\left({{G \over {{\rm{K}}\;{\rm{J}}{{\rm{y}}^{- 1}}}}} \right)^{- 1}}{\left({{{\Delta \nu} \over {{\rm{MHz}}}}} \right)^{- 1/2}}{\left({{{{t_{{\mathop{\rm int}}}}} \over {\rm{s}}}} \right)^{- 1/2}}{\left({{W \over P}} \right)^{3/2}}.$$
(8)

Here, Spsr is the flux density of the pulsar, Trec and Tsky are the receiver and sky noise temperatures, G is the antenna gain, Δν is the observing bandwidth, tint is the integration time, W is the pulse width and P is the pulse period (we assume WP). Optimal results are thus obtained for observations of short period pulsars with large flux densities and small 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 Section 2.4, pulses emitted at lower radio frequencies travel slower and arrive later than those emitted at higher frequencies. This process has the effect of “stretching” the pulse across a finite receiver bandwidth, increasing W and therefore increasing ϕTOA. For normal pulsars, dispersion can largely be compensated for by the incoherent de-dispersion process outlined in Section 3.1.

The short periods of millisecond pulsars offer the ultimate in timing precision. In order to fully exploit this, a better 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 [110, 185] 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 on-line either using finite impulse response filter devices [16] or completely in software [290, 297]. Off-line data reduction, while disk-space limited, allows for more flexible analysis schemes [21].

The maximum time resolution obtainable via coherent dedispersion is the inverse of the total receiver bandwidth. The current state of the art is the detection [111] of features on nanosecond timescales in pulses from the 33-ms pulsar B0531+21 in the Crab nebula shown in Figure 21. Simple light travel-time arguments can be made to show that, in the absence of relativistic beaming effects [104], these incredibly bright bursts originate from regions less than 1 m in size.

Figure 21
figure 21

A 120 µs window centred on a coherently-dedispersed giant pulse from the Crab pulsar showing high-intensity nanosecond bursts. Figure provided by Tim Hankins [111].

4.2 The timing model

Ideally, in order to model the rotational behaviour of the neutron star, we require TOAs measured by an inertial observer. 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 solar system centre-of-mass (barycentre) can be regarded as an inertial frame. It is now standard practice [130] to transform the observed TOAs to this frame using a planetary ephemeris such as the JPL DE405 [299]. The transformation is summarised as the difference between barycentric (\({\mathcal T}\)) and observed (t) TOAs:

$${\mathcal T} - t = {{\underline r \cdot\underline {\hat s}} \over c} + {{{{(\underline r \cdot\underline {\hat s})}^2} - \vert\underline r {\vert ^2}} \over {2cd}} + \Delta {t_{{\rm{rel}}}} - \Delta {t_{{\rm{DM}}}},$$
(9)

where r is the position of the observatory with respect to the barycentre, ŝ is a unit vector in the direction of the pulsar at a distance d and c is the speed of light. The first term on the right hand side of Equation (9) is the light travel time from the observatory to the solar system barycentre. Incoming pulses from all but the nearest pulsars can be approximated by plane wavefronts. The second term, which represents the delay due to spherical wavefronts, yields the parallax and hence d. This has so far only been measured for five nearby millisecond pulsars [276, 320, 174, 287]. The term Δtrel represents the Einstein and Shapiro corrections due to general relativistic time delays in the solar system [18]. Since measurements can be carried out at different observing frequencies with different dispersive delays, TOAs are generally referred to the equivalent time that would be observed at infinite frequency. This transformation is the term ΔtDM (see also Equation (1)).

Following the accumulation of a number of TOAs, a surprisingly simple model is usually sufficient to account for the TOAs during the time span of the observations and to predict the arrival times of subsequent pulses. The model is a Taylor expansion of the rotational frequency Ω = 2π/P about a model value Ω0 at some reference epoch \({\mathcal T_0}\). The model pulse phase is

$$\phi ({\mathcal T}) = {\phi _0} + ({\mathcal T} - {{\mathcal T}_0}){\Omega _0} + {1 \over 2}{({\mathcal T} - {{\mathcal T}_0})^2}{\dot \Omega _0} + \ldots,$$
(10)

where \({\mathcal T}\) is the barycentric time and φ0 is the pulse phase at \({{\mathcal T}_0}\). 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 in Figure 22. Ideally, the residuals should have a zero mean and be free from any systematic trends (see Panel a of Figure 22). To reach this point, however, 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 incorporated into the model.

Figure 22
figure 22

Timing model residuals for PSR B1133+16. Panel a: Residuals obtained from the best-fitting model which includes period, period derivative, position and proper motion. Panel b: Residuals obtained when the period derivative term is set to zero. Panel c: Residuals showing the effect of a 1-arcmin position error. Panel d: Residuals obtained neglecting the proper motion. The lines in Panels b−d show the expected behaviour in the timing residuals for each effect. Data provided by Andrew Lyne.

From Equation (10), an error in the assumed Ω0 results in a linear slope with time. A parabolic trend results from an error in \({{\dot \Omega}_0}\) (see Panel b of Figure 22). Additional effects will arise if the assumed position of the pulsar (the unit vector ŝ in Equation (9)) used in the barycentric time calculation is incorrect. A position error results in an annual sinusoid (see Panel c of Figure 22). A proper motion produces an annual sinusoid of linearly increasing magnitude (see Panel d of Figure 22).

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 with 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 [155, 152] defined at midnight UT on December 5 1988! Astrometric measurements are no less impressive, with position errors of ∼ 20 μarcsec being presently possible [333].

4.3 Timing stability

Ideally, after correctly applying a timing model, we would expect a set of uncorrelated timing residuals with a zero mean and a Gaussian scatter with a standard deviation consistent with the measurement uncertainties. As can be seen in Figure 23, this is not always the case; the residuals of many pulsars exhibit a quasi-periodic wandering with time.

Figure 23
figure 23

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”. Data taken from the Jodrell Bank timing program [281, 124]. Figure provided by Andrew Lyne.

Such “timing noise” is most prominent in the youngest of the normal pulsars [213, 72] and present at a lower level in the much older millisecond pulsars [155, 11]. While the physical processes of this phenomenon 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 to processes in the magnetosphere [66, 65].

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 Section 4.1. Much progress in this area has been made by groups at Princeton [258], Berkeley [325], Jodrell Bank [326], UBC [323], Swinburne [305] and ATNF [12]. From high quality observations spanning over a decade [274, 275, 155], these groups have demonstrated that the timing stability of millisecond pulsars over such time-scales is comparable to terrestrial atomic clocks.

This phenomenal stability is demonstrated in Figure 24 which shows σz [214], a parameter closely resembling the Allan variance used by the clock community to estimate the stability of atomic clocks [308, 1]. Both PSRs B1937+21 and B1855+09 seem to be limited by a power law component which produces a minimum in σz after 2 yr and 5 yr respectively. This is most likely a result of a small amount of intrinsic timing noise [155]. Although the baseline for the bright millisecond pulsar J0437-4715 is shorter, its σz is already an order of magnitude smaller than the other two pulsars or the atomic clocks. Timing observations of an array of millisecond pulsars in the context of detecting gravitational waves from the Big Bang are discussed further in Section 4.5.3.

Figure 24
figure 24

The fractional stability of three millisecond pulsars compared to an atomic clock. Both PSRs B1855+09 and B1937+21 are comparable, or just slightly worse than, the atomic clock behaviour over timescales of a few years [214]. More recent timing of the millisecond pulsar J0437-4715 [126] indicates that it is more stable than the atomic clock.

4.4 Timing binary pulsars

For binary pulsars, the simple timing model introduced in Section 4.2 needs to be extended to incorporate the additional motion 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 x, 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. Analogous to the radial velocity curve in a spectroscopic binary, for binary pulsars the orbit is described by the apparent pulse period against time. An example of this is shown in Panel a of Figure 25. Alternatively, when radial accelerations can be measured, the orbit can also be visualised in a plot of acceleration versus period as shown in Panel b of Figure 25. This method is particularly useful for determining binary pulsar orbits from sparsely sampled data [102].

Figure 25
figure 25

Panel a: Keplerian orbital fit to the 669-day binary pulsar J0407+1607 [190]. Panel b: Orbital fit in the period-acceleration plane for the globular cluster pulsar 47 Tuc S [102].

Constraints on the masses of the pulsar mp and the orbiting companion mc can be placed by combining x and Pb to obtain the mass function

$${f_{{\rm{mass}}}} = {{4{\pi ^2}} \over G}{{{x^3}} \over {P_{\rm{b}}^2}} = {{{{({m_{\rm{c}}}\sin i)}^3}} \over {{{({m_{\rm{p}}} + {m_{\rm{c}}})}^2}}},$$
(11)

where G is Newton’s gravitational constant and i is the (initially unknown) angle between the orbital plane and the plane of the sky (i.e. an orbit viewed edge-on corresponds to i = 90°). In the absence of further information, the standard practice is to consider a random distribution of inclination angles. Since the probability that i is less than some value i0 is p(< i0) = 1 − cos(i0), the 90% confidence interval for i is 26° < i < 90°. For an assumed pulsar mass, the 90% confidence interval for mc can be obtained by solving Equation (11) for i = 26° and 90°. If the sum of the masses M = mp + mc can be determined (e.g., through a measurement of relativistic periastron advance described below), then the condition sin i < 1 sets a lower limit on the companion mass mc > (fmassM2)1/3 and a corresponding upper limit on the pulsar mass.

Although most of the presently known binary pulsar systems can be adequately timed using Kepler’s laws, there are a number which require an additional set of “post-Keplerian” (PK) parameters which have a distinct functional form for a given relativistic theory of gravity. In general relativity (GR) the PK formalism gives the relativistic advance of periastron

$$\dot \omega = 3{\left({{{{P_{\rm{b}}}} \over {2\pi}}} \right)^{- 5/3}}{({T_ \odot}M)^{2/3}}{(1 - {e^2})^{- 1}},$$
(12)

the time dilation and gravitational redshift parameter

$$\gamma = e{\left({{{{P_{\rm{b}}}} \over {2\pi}}} \right)^{1/3}}T_ \odot ^{2/3}{M^{- 4/3}}{m_{\rm{c}}}({m_{\rm{p}}} + 2{m_{\rm{c}}}),$$
(13)

the rate of orbital decay due to gravitational radiation

$${\dot{P}_{\rm{b}}} = - {{192\pi} \over 5}{\left({{{{P_{\rm{b}}}} \over {2\pi}}} \right)^{- 5/3}}\left({1 + {{73} \over {24}}{e^2} + {{37} \over {96}}{e^4}} \right){(1 - {e^2})^{- 7/2}}T_ \odot ^{5/3}{m_{\rm{p}}}{m_{\rm{c}}}{M^{- 1/3}}$$
(14)

and the two Shapiro delay parameters

$$r = {T_ \odot}{m_{\rm{c}}}$$
(15)

and

$$s = x{\left({{{{P_{\rm{b}}}} \over {2\pi}}} \right)^{- 2/3}}T_ \odot ^{- 1/3}{M^{2/3}}m_{\rm{c}}^{- 1}$$
(16)

which describe the delay in the pulses around superior conjunction where the pulsar radiation traverses the gravitational well of its companion. In the above expressions, all masses are in solar units, M = mp + mc, x = ap sin i/c, s ≡ sin i and TGM/c3 = 4.925490947 µs. Some combinations, or all, of the PK parameters have now been measured for a number of binary pulsar systems. Further PK parameters due to aberration and relativistic deformation [80] are not listed here but may soon be important for the double pulsar [185].

The key point in the PK definitions is that, given the precisely measured Keplerian parameters, the only two unknowns are the masses of the pulsar and its companion, mp and mc. Hence, from a measurement of just two PK parameters (e.g., \({\dot \omega}\) and γ) one can solve for the two masses and, using Equation (11), find the orbital inclination angle i. If three (or more) PK parameters are measured, the system is “overdetermined” and can be used to test GR (or, more generally, any other theory of gravity) by comparing the third PK parameter with the predicted value based on the masses determined from the other two.

The first binary pulsar used to test GR in this way was PSR B1913+16 discovered by Hulse & Taylor in 1974 [129]. Measurements of three PK parameters (\({\dot \omega}\), γ and b) were obtained from long-term timing observations at Arecibo [312, 313]. The measurement of 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 Figure 26. This figure includes recent Arecibo data taken since the upgrade of the telescope in the mid 1990s. The measurement of orbital decay, now spanning a 30-yr baseline [336], is within 0.2% of the GR prediction and provided the first indirect evidence for the existence of gravitational waves. Hulse and Taylor were awarded the 1993 Nobel Physics prize [314, 128, 309] in recognition of their discovery of this remarkable laboratory for testing GR.

Figure 26
figure 26

Orbital decay in the binary pulsar B1913+16 system demonstrated as an increasing orbital phase shift for periastron passages with time. The GR prediction due entirely to the emission of gravitational radiation is shown by the parabola. Figure provided by Joel Weisberg [336].

More recently, five PK parameters have been measured for PSRs B1534+12 [298] and J0737-3039A [161]. For PSR B1534+12, the test of GR comes from measurements of \({\dot \omega}\), γ and s, where the agreement between theory and observation is within 0.7% [298]. This test will improve in the future as the timing baseline extends and a more significant measurement of r can be made. Although a significant measurement of exists, it is known to be contaminated by kinematic effects which depend on the assumed distance to the pulsar [293]. Assuming GR to be correct, the observed and theoretical values can be reconciled to provide a “relativistic measurement” of the distance d = 1.04 ± 0.03 kpc [288]. Prospects for independent parallax measurements of the distance to this pulsar using radio interferometry await more sensitive telescopes [289].

For PSR J0737-3039, where two independent pulsar clocks can be timed, five PK parameters of the 22.7-ms pulsar “A” have been measured as well as two additional constraints from the measured mass function and projected semi-major axis of the 2.7-s pulsar “B”. In terms of a laboratory for GR, then, J0737-3039 promises to go well beyond the results possible from PSRs B1913+16 and B1534+12. A useful means of summarising the limits so far is Figure 27 which shows the allowed regions of parameter space in terms of the masses of the two pulsars. The shaded regions are excluded by the requirement that sin i < 1. Further constraints are shown as pairs of lines enclosing permitted regions as predicted by GR. The measurement of \(\dot \omega = 16.899 \pm 0.001\) deg yr−1 gives the total system mass M = 2.587±0.001 M. The measurement of the projected semi-major axes of both orbits gives the mass ratio R = 1.071 ± 0.001. The mass ratio measurement is unique to the double pulsar system and rests on the basic assumption that momentum is conserved. This constraint should apply to any reasonable theory of gravity. The intersection between the lines for \({\dot \omega}\) and R yield the masses of A and B as mA = 1.338±0.001 M and mB = 1.249±0.001 M. From these values, using Equations (13–16) the expected values of γ, b, r and s may be calculated and compared with the observed values. These four tests of GR all agree with the theory to within the uncertainties. Currently the tightest constraint is the Shapiro delay parameter s where the observed value is in agreement with GR at 0.1% level.

Figure 27
figure 27

‘Mass-mass’ diagram showing the observational constraints on the masses of the neutron stars in the double pulsar system J0737-3039. Inset is an enlarged view of the small square encompassing the intersection of the tightest constraints. Figure provided by Michael Kramer [161].

Less than two years after its discovery, the double pulsar system has already surpassed the three decades of monitoring PSR B1913+16 and over a decade of timing PSR B1534+12 as a precision test of GR. On-going precision timing measurements of the double pulsar system should soon provide even more stringent and new tests of GR. Crucial to these measurements will be the timing of the 2.7-s pulsar B, where the observed profile is significantly affected by A’s relativistic wind [198, 216]. A careful decoupling of these profile variations is required to accurately measure TOAs for this pulsar and determine the extent to which the theory-independent mass ratio R can be measured.

The relativistic effects observed in the double pulsar system are so large that corrections to higher post-Newtonian order may soon need to be considered. For example, \({\dot \omega}\) may be measured precisely enough to require terms of second post-Newtonian order to be included in the computations [81]. In addition, in contrast to Newtonian physics, GR predicts that the spins of the neutron stars affect their orbital motion via spin-orbit coupling. This effect would most clearly be visible as a contribution to the observed \({\dot \omega}\) in a secular [26] and periodic fashion [337]. For the J0737-3039 system, the expected contribution is about an order of magnitude larger than for PSR B1913+16 [198]. As the exact value depends on the pulsars’ moment of inertia, a potential measurement of this effect allows the moment of inertia of a neutron star to be determined for the first time [81]. Such a measurement would be invaluable for studies of the neutron star equation of state and our understanding of matter at extreme pressure and densities [168].

The systems discussed above are all double neutron star binaries. A further self-consistency test of GR has recently been made in the 4.7-hr relativistic binary J1141-6545, where the measurement [25] of \({\dot \omega}\), γ and b yield a pulsar mass of 1.30 ± 0.02 and a companion mass of 0.99 ± 0.02 M. Since the mass of the companion is some seven standard deviations below the mean neutron star mass (see Figure 28), it is most likely a white dwarf. The observed b = −(4±1) × 10−13 is consistent, albeit with limited precision, with the predicted value from GR (−3.8 × 10−13). Continued timing should reduce the relative error in b down to 1% by 2010 [25].

Figure 28
figure 28

Distribution of neutron star masses as inferred from timing observations of binary pulsars [292]. The vertical dotted line shows the canonical neutron star mass of 1.4 M. Figure provided by Ingrid Stairs.

PK parameters have now been measured for a number of other binary pulsars which provide interesting constraints on neutron star masses [318, 292]. Figure 28 shows the distribution taken from a recent compilation [292]. While the young pulsars and the double neutron star binaries are consistent with, or just below, the canonical 1.4 M, we note that the millisecond pulsars in binary systems have, on average, significantly larger masses. This provides strong support for their formation through an extended period of accretion in the past, as discussed in Section 2.6.

4.5 Pulsar timing and gravitational wave detection

We have seen in the previous Section 4.4 how pulsar timing can be used to provide indirect evidence for the existence of gravitational waves from coalescing stellar-mass binaries. In this final section, we look at how pulsar timing might soon be used to detect gravitational radiation directly. The idea to use pulsars as natural gravitational wave detectors was first explored in the late 1970s [277, 85]. 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 perturb the local spacetime metric and cause a change in the observed rotational frequency of the pulsar. For regular pulsar timing observations with typical TOA uncertainties of ϵTOA, this “detector” would be sensitive to waves with dimensionless amplitudes hϵTOA/T and frequencies f ∼ 1/T [31, 39].

4.5.1 Probing the gravitational wave background

Many cosmological models predict that the Universe is presently filled with a low-frequency stochastic gravitational wave background (GWB) produced during the big bang era [243]. A significant component [260, 134] is the gravitational radiation from massive black hole mergers. 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 Section 4.2, 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.

For a GWB with a flat energy spectrum in the frequency band f ± f/2 there is an additional contribution to the timing residuals σg [85]. When fT ≫ 1, the corresponding wave energy density

$${\rho _{\rm{g}}} = {{243{\pi ^3}{f^4}\sigma _{\rm{g}}^2} \over {204G}}.$$
(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 the energy density required to close the Universe

$${\rho _{\rm{c}}} = {{3H_0^2} \over {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.

This technique was first applied [273] to a set of TOAs for PSR B1237+25 obtained from regular observations over a period of 11 years as part of the JPL pulsar timing programme [89]. 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 (see Section 4.3). By ascribing the rms scatter in the residuals (σ = 240 ms) to the GWB, the limit ρg/ρc ≲ 4 × 10−3h−2 for a centre frequency f = 7 nHz.

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 (see Section 4.3). In the intervening period, more elaborate techniques had been devised [31, 39, 303] 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 [31] it is convenient to define

$${\Omega _{\rm{g}}} = {1 \over {{\rho _{\rm{c}}}}}{{{\rm{d}}\log {\rho _{\rm{g}}}} \over {{\rm{d}}\log f}},$$
(19)

i.e. the energy density of the GWB per logarithmic frequency interval relative to ρc. With this definition, the power spectrum of the GWB [125, 39] is

$${\mathcal P}(f) = {{G{\rho _{\rm{g}}}} \over {3{\pi ^3}{f^4}}} = {{H_0^2{\Omega _{\rm{g}}}} \over {8{\pi ^4}{f^5}}} = 1.34 \times {10^4}{\Omega _{\rm{g}}}{h^2}f_{{\rm{y}}{{\rm{r}}^{- 1}}}^{- 5}\mu {{\rm{s}}^2}{\rm{yr}},$$
(20)

where \({f_{{\rm{yr}^{- 1}}}}\) is frequency in cycles per year. The long-term timing stability of B1937+21, discussed in Section 4.3, limits its use for periods ≳ 2 yr. Using the more stable residuals for PSR B1855+09, and a rigorous statistical analysis [319], the current 95% confidence upper limit is Ωgh2 < 10−8. This limit is difficult to reconcile with most cosmic string models for galaxy formation [45, 319].

For binary pulsars, the orbital period provides an additional clock for measuring the effects of gravitational waves. 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 [31]. 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 [160].

4.5.2 Constraints on massive black hole binaries

In addition to probing the GWB, pulsar timing is beginning to place interesting constraints on the existence of massive black hole binaries. Arecibo data for PSRs B1937+21 and J1713+0747 already make the existence of an equal-mass black hole binary in Sagittarius A* unlikely [176]. More recently, timing data from B1855+09 have been used to virtually rule out the existence of a proposed supermassive black hole as the explanation for the periodic motion seen at the centre of the radio galaxy 3C66B [304].

A simulation of the expected modulations of the timing residuals for the putative binary system, with a total mass of 5.4 × 1010 M, is shown along with the observed timing residuals in Figure 29. Although the exact signature depends on the orientation and eccentricity of the binary system, Monte Carlo simulations show that the existence of such a massive black hole binary is ruled out with at least 95% confidence [135].

Figure 29
figure 29

Top panel: Observed timing residuals for PSR B1855+09. Bottom panel: Simulated timing residuals induced from a putative black hole binary in 3C66B. Figure provided by Rick Jenet [135].

4.5.3 A millisecond pulsar timing array

A natural extension of the single-arm detector concept discussed above is the idea of using timing data for a number of pulsars distributed over the whole sky to detect gravitational waves [118]. Such a “timing array” 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 common to all pulsars in the array. To illustrate this, consider the fractional frequency shift of the ith pulsar in an array

$${{\delta {\nu _i}} \over {{\nu _i}}} = {\alpha _i}{\mathcal A}(t) + {{\mathcal N}_i}(t),$$
(21)

where α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 \({\mathcal A}\). The timing noise intrinsic to the pulsar is characterised by the function \({{\mathcal N}_i}\). The result of a cross-correlation between pulsars i and j is then

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

where the bracketed terms indicate cross-correlations. Since the wave function and the noise contributions from the two pulsars are independent, the cross correlation tends to \({\alpha _i}{\alpha _j}\left\langle {{\mathcal A^2}} \right\rangle\) 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 [117] and allows the separation of the effects of clock and ephemeris errors from the GWB [97].

A preliminary analysis applying the timing array concept to Arecibo data for three millisecond pulsars (B1937+21, B1855+09 and J1737+0747) now spanning 17 years has reduced the energy density limit to Ωgh2 < 2 × 10−9 [175]Footnote 3. This limit, and the region of the gravitational wave energy density spectrum probed by the current pulsar timing array is shown in Figure 30 where it can be seen that the pulsar timing regime is complementary to the higher frequency bands of LISA and LIGO.

Figure 30
figure 30

Summary of the gravitational wave spectrum showing the location in phase space of the pulsar timing array (PTA) and its extension with the Square Kilometre Array (SKA). Figure updated by Michael Kramer [162] from an original design by Richard Battye.

A number of long-term timing projects are now underway to make a large-scale pulsar timing array a reality. At Arecibo and Green Bank, regular timing of a dozen or more millisecond pulsars are now carried out using on-line data acquisition systems by groups at Berkeley [325], Princeton [258] and UBC [323]. A similar system is also being commissioned at Jodrell Bank [142]. The Berkeley group have installed identical sets of data taking equipment at Effelsberg and Nancay which has enabled these telescopes to make regular high-quality millisecond pulsar timing observations [174, 68]. In the southern hemisphere, a high-precision timing program at Parkes [333] has been in existence since the discovery bright nearby millisecond pusar J0437-4715 [138]. Since February 2004, weekly Parkes observations of ∼20 millisecond pulsars are already achieving a limit of Ωgh2 < 10−4, with the ultimate goal being to reach a limits of Ωgh2 < 10−10 before 2010 [121].

Looking further ahead, the increase in sensitivity provided by the Square Kilometre Array [132, 162] should further improve the limits of the spectrum probed by pulsar timing. As Figure 30 shows, the SKA could provide up to two orders of magnitude improvement over current capabilities.

4.6 Going further

For further details on the techniques and prospects of pulsar timing, a number of excellent reviews are available [18, 308, 28, 291, 185]. Two audio files and slides from lectures at the centennial meeting of the American Physical Society are also available on-line [15, 341]. The tests of general relativity discussed in Section 4.4, along with further effects, are reviewed elsewhere in this journal by Will [342] and Stairs [291]. Further discussions on the realities of using pulsars as gravity wave detectors can be found in other review articles [271, 14, 121].

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

Through an understanding of the Galactic population of radio pulsars summarised in Section 3 it is possible to predict the detection statistics of terrestrial gravitational wave detectors to nearby rapidly spinning neutron stars (see Section 3.3), as well as coalescing relativistic binaries at cosmic distances (see Section 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 (see Section 4.4) and as natural detectors of gravitational radiation (see Section 4.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 continuing 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 Section 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 using a multibeam system on the Arecibo telescope [227] are now being carried out. Future surveys with the Square Kilometre Array [132] will ultimately provide a far more complete census of the Galactic pulsar population. With the double pulsar now checked off the list of pulsar milestones given in the last version of this review [180], we look forward to new discoveries. Two possible “holy grails” of pulsar astronomy which may soon be found are the following ones:

  • A pulsar-black hole binary system: Following the wide variety of science possible from the new double pulsar J0737-3039 (see Sections 2.6, 3.4.1 and 4.4), a radio pulsar with a blackhole companion would without doubt be a fantastic gravitational laboratory. Excellent places to look for such systems are globular clusters and the Galactic centre [247].

  • A sub-millisecond pulsar: With the new discoveries in Terzan 5, it now appears that the original millisecond pulsar, B1937+21 is no longer the most rapidly rotating neutron star known [233]. Do neutron stars with kHz spin rates exist? Searches now have sensitivity to such objects [41] and a discovery of even one would constrain the equation of state of matter at high densities. As mentioned in Section 2.2, the R-mode instability may prevent neutron stars from reaching such high spin rates [36, 60, 7].

Pulsar astronomy remains an extremely active area of modern astrophysics and the next decade will undoubtedly continue to produce new results from currently known objects as well as new surprises. In my opinion, pulsar research is currently limited by a shortage of researchers, and not necessarily resources. Keen graduate students more than ever are needed to help shape the future of this exciting and continually evolving field.