Pulsars and Gravity

R. N. Manchester

Pulsars are wonderful gravitational probes. Their tiny size and stellar mass give their rotation periods a stability comparable to that of atomic frequency standards. This is especially true of the rapidly rotating “millisecond pulsars” (MSPs). Many of these rapidly rotating pulsars are in orbit with another star, allowing pulsar timing to probe relativistic perturbations to the orbital motion. Pulsars have provided the most stringent tests of theories of relativistic gravitation, especially in the strong-field regime, and have shown that Einstein’s general theory of relativity is an accurate description of the observed motions. Many other gravitational theories are effectively ruled out or at least severely constrained by these results. MSPs can also be used to form a “Pulsar Timing Array” (PTA). PTAs are Galactic-scale interferometers that have the potential to directly detect nanohertz gravitational waves from astrophysical sources. Orbiting super-massive black holes in the cores of distant galaxies are the sources most likely to be detectable. Although no evidence for gravitational waves has yet been found in PTA data sets, the latest limits are seriously constraining current ideas on galaxy and black-hole evolution in the early Universe.

Keywords: gravity — gravitational waves — pulsars — pulsar timing — general relativity

1 Introduction

Pulsars are rotating neutron stars that emit beams of radiation which sweep across the sky as the star rotates. A beam sweeping across the Earth may be detected as a pulse that repeats with a periodicity equal to the rotation period of the star. Because of the large mass of neutron stars,  M, and their tiny size, radii  km, (see Ref. ? for a review of neutron-star properties) the rotation period of neutron stars is incredibly stable, with a stability comparable to that of the best atomic clocks on Earth. This great period stability, combined with the fact that pulsars are often in a binary orbit with another star, makes them wonderful probes of relativistic gravity. Tiny perturbations to their period resulting from, for example, relativistic effects in a binary orbit, may be detected and compared with the predictions of a gravitational theory. Pulsars may also be used as detectors for gravitational waves passing through the Galaxy. To separate the effects of gravitational waves from other perturbations, signals from pulsars in different directions on the sky must be compared – exactly analogous to the way laser-interferometer gravitational-wave detectors compare laser phases in perpendicular arms.

More than 2400 pulsars are now known.aaaSee the ATNF Pulsar Catalogue: http://www.atnf.csiro.au/research/pulsar/psrcat. The vast majority of them reside in our Galaxy, typically at distances of a few thousand light-years from the Sun. Their pulse periods, , range between 1.4 milliseconds and 12 seconds and fall into two main groups. The so-called “normal” pulsars have periods longer than about 30 milliseconds and the “millisecond” pulsars (MSPs) have shorter periods. MSPs comprise about 15% of the known pulsar population.

Although pulsar periods are very stable, they are not constant. In their own reference frame, all pulsars are slowing down, albeit very slowly. Pulsars are powered by their rotational energy. They have extremely strong magnetic fields and, as they spin, they emit streams of relativistic particles and so-called “dipole radiation”, electromagnetic waves with period equal to the rotation period of the star. These carry energy away from the star resulting in a steady increase in the pulse period. Typical rates of period increase, , are a part in for normal pulsars and much less for MSPs. Assuming that the surrounding magnetic field is predominantly dipole, the characteristic age of a pulsar is given by


and their surface dipole magnetic field strength (in gaussbbb1 gauss (G) =  tesla) is


Figure 1 shows the distribution of pulsars on the plane, with several different types of pulsars indicated. For most normal pulsars, is between and years and is between and  G. For MSPs, the corresponding ranges are to  yrs and to  G.

Fig. 1: Distribution of pulsars on the plane with radio-loud and radio-quiet pulsars indicated. Binary pulsars are indicated by a circle around the point and, for those with a neutron-star companion, the circle is filled in. Lines of constant characteristic age () and surface-dipole magnetic field () are shown.

Normal pulsars and MSPs have quite different evolutionary histories. Most if not all normal pulsars are formed in supernova explosions at the death of a massive star. They age with relatively constant until the pulse emission mechanism begins to fail when reaches about  yrs. Many young pulsars ( yrs) are located within supernova remnants, with the most famous example being the Crab pulsar, PSR B0531+21, located near the centre of the Crab Nebula. Most young pulsars lie relatively close to the Galactic Plane, consistent with the idea that they are formed from massive stars. MSPs on the other hand are much more widely distributed in the Galactic halo. They are believed to originate from old, slowly rotating and probably dead neutron stars that accrete matter and angular momentum from an evolving binary companion. This “recycling” process increases their spin rate so that they have periods in the millisecond range and re-energises the beamed emission (see, e.g., Ref. ?).

Figure 1 shows that the majority of MSPs are members of a binary system, consistent with this formation scenario. The accretion also suppresses their apparent magnetic field, so that MSPs have values about five orders of magnitude less than normal pulsars. Since the level of intrinsic period irregularities is related to , it is this low that makes MSPs extremely stable clocks and suitable tools for the study of relativistic gravitation. Figure 1 also indicates the class of radio-quiet pulsars. Essentially all of the 72 known radio-quiet pulsars are young and solitary. They include the so-called “magnetars” which lie in the upper right side of the diagram where dipole magnetic fields are strongest.

Most double-neutron-star systems lie in the zone between the MSPs and the normal pulsars. These systems are believed to have been partially recycled prior to the formation of the second-born neutron star. In almost every case, the observed pulsar is the recycled one since it has a much longer active lifetime than the newly formed young pulsar. A famous exception to this is the Double Pulsar (PSR J07373039A/B) in which the second-born star (B) is (or was) still visible . On Figure 1 the B pulsar is the solitary double-neutron-star system on the right side of the plot. The double-neutron-star system identified near the middle of the plot contains the relatively young pulsar, PSR J1906+0746. There is some doubt about the nature of the companion in this system – it could be a heavy white dwarf . These double-neutron-star systems and their use as probes of relativistic gravity are discussed in some detail in Section 2.1 below.

1.1 Pulsar Timing

The most important contributions of pulsars to investigations of gravitational theories and gravitational waves rely on precision pulsar timing observations. These allow both the relativistic perturbations to binary orbits to be studied in detail and the potential detection of the tiny period fluctuations generated by gravitational waves passing through our Galaxy. Because of the importance of pulsar timing to these studies, we give here a brief review of its basic principles.

The basic observable in pulsar timing is the time of arrival (ToA) of a pulse at an observatory. In fact, because of signal/noise limitations and the intrinsic fluctuation of individual pulse shapes, ToA measurements are based on mean pulse profiles formed by synchronously averaging the data, typically for times between several minutes and an hour. The time at which a fiducial pulse phase (usually near the pulse peak) arrives at the telescope is determined by cross-correlating the observed mean pulse profile with a standard pulse template. A series of these ToAs is measured over many days, months, years and even decades for the pulsar of interest.

These observatory ToAs are affected by the rotational and orbital motion of the Earth (and for satellite observatories, the orbital motion of the satellite). To remove these effects, the observed ToAs are referred to the solar-system barycentre (centre of mass) which is assumed to be inertial (unaccelerated) with respect to the distant Universe.cccThis neglects any acceleration of the solar-system barycentre resulting from, for example, Galactic rotation. For some precision timing experiments, such effects are taken into account at a later stage of the analysis. This correction makes use of a solar-system ephemeris giving predictions of the position of the centre of the Earth with respect to the solar-system barycentre. Such ephemerides, for example, the Jet Propulsion Laboratory ephemeris DE 421, are generated by fitting relativistic models to planetary and spacecraft data. The correction also takes into account the relativistic variations in terrestrial time resulting from the Earth’s motion.

The resulting barycentric ToAs are then compared with predicted pulse ToAs based on a model for the pulsar. The pulsar model can have 20 or more parameters; generally included are the pulse frequency (), frequency time derivative (), the pulsar position and the five Keplerian binary parameters if the pulsar is a member of a binary system. If the data are recorded at different radio frequencies, it is also necessary to include the dispersion measure (DM) which quantifies the frequency-dependent delay suffered by the pulses as they propagate through the interstellar medium.

The differences between the observed and predicted ToAs are known as “timing residuals”. Errors in any of the model parameters result in systematic variations in the timing residuals as a function of time. For example, if the model pulse frequency is too small, the residuals will grow linearly as illustrated in Figure 2. The required correction to the pulse frequency is given by the slope of this variation. An error in the pulsar position results in an annual sine curve which arises from the barycentric correction. The phase and amplitude of this curve give the corrections to the two position coordinates. Similarly, a pulsar proper motion results in a linearly growing sine curve (away from the reference epoch). For a sufficiently close pulsar, the curvature of the wavefront results in a biannual term in the residuals, and offsets in binary parameters result in terms which vary with the orbital periodicity.

Since offsets in different parameters in general result in different systematic residual variations, it is possible to do a least-squares fit to the observed residuals for the correction to any desired parameter.dddThe data span must be sufficiently long to avoid excessive covariance between the variations for different fitted parameters. The precision of the fitted parameters is often amazingly high. For example, the relative precision of the pulse frequency determination is , where is the typical uncertainty in the ToAs and is the data span. Since, for MSPs at least, is often  s and is often many years, the relative precision of the measured can easily exceed 1:10. Similarly pulsar positions can be measured to micro-arcseconds and binary eccentricities measured to . These very high precisions often allow higher-order terms, for example, resulting from relativistic perturbations to the binary orbit, to be measured.

Fig. 2: Variations in timing residuals for offsets in several parameters. The model pulsar has a pulse frequency of 300 Hz (3.3 ms period) and is in an eccentric () binary orbit of period 190 days. The reference epochs for period, position and binary phase are at the middle of the plotted range.

2 Tests of Relativistic Gravity

2.1 Tests of General Relativity with Double-Neutron-Star Systems

2.1.1 The Hulse-Taylor binary, PSR B1913+16

The discovery at Arecibo Observatory in 1974 of the first-known binary pulsar, PSR B1913+16, by Hulse and Taylor was remarkable in a number of respects. Firstly, it showed that pulsars with short pulse periods but large characteristic ages (59 ms and  yr, respectively, for PSR B1913+16) could exist. The period of PSR B1913+16 was second only to the Crab pulsar, but its age was enormously greater than that of other short-period pulsars known at the time. Secondly, it was in a binary orbit with a relatively massive star, very likely another neutron star, showing that an evolutionary pathway to such systems existed. Thirdly, its orbital period was extraordinarily short, only 7.75 h, its eccentricity large, , and, as shown in Figure 3, its maximum orbital velocity very high,  km s or 0.1% of the velocity of light. As was immediately recognised by Hulse and Taylor, these properties opened up the possibility that relativistic perturbations to the orbit were potentially measurable. Lowest-order relativistic effects go as , and so the variations are of order 1:10, enormous by the standards of pulsar measurements.

Fig. 3: Velocity curve for the Hulse-Taylor binary pulsar, PSR B1913+16. (Ref. ?)

Relativistic effects in binary motion can be expressed in terms of “post-Keplerian” parameters that describe departures from Keplerian motion (see, e.g., Ref. ?). The first such parameter to be observed was periastron precession . In Einstein’s general theory of relativity (GR) the rate of periastron precession (averaged over the binary orbit) is given by:


where is the longitude of periastron measured from the ascending node, is the orbital period, s, is the Gravitational Constant, is the mass of the Sun, is the orbital eccentricity and , the sum of the pulsar mass and the companion mass in solar units. It is worth noting that this is the same effect as the excess perihelion advance of Mercury that was used by Einstein as an observational verification of GR. The relativistic effect for Mercury is just 43 arcsec per century, minuscule compared to the per year predicted and observed for PSR B1913+16.

The next most significant parameter, normally labelled , describes the combination of gravitational redshift and 2nd-order (or transverse) Doppler shift, both of which have the same dependence on orbital phase. In GR


Since the Keplerian parameters are very well determined, measurement of and gives two equations in two unknowns, and , and so the two stellar masses can be determined. For PSR B1913+16, these are both close to 1.4 , confirming the double-neutron-star nature of the system. An important consequence of this was that the two stars could safely be treated as point masses in GR, thereby allowing precise tests of the theory.

Given the Keplerian parameters and the two masses, the next post-Keplerian parameter, orbital decay due to the emission of gravitational waves from the system, given in GR by


could be predicted and compared with observation. This therefore constituted a unique test of GR in strong-field gravity. After just four years, the term was well measured and shown to be consistent with the GR prediction . The orbital decay, including more recent results, is illustrated in Figure 4; the ratio of observed to predicted orbit decay is . It should be mentioned that this near-perfect agreement relies on correcting the observed for the differential acceleration of the solar system and the pulsar system in the gravitational field of the Galaxy. This correction is or about 1% of the observed value, and the uncertainty in the final result is dominated by the uncertainty in this correction. Unfortunately, it is unlikely that this uncertainty will be significantly reduced in the near future since it mainly depends on the poorly known distance to the binary system.

Fig. 4: Comparison of the observed and predicted orbital period decay in PSR B1913+16. The orbital decay is quantified by the shift in the time of periastron passage with respect to a non-decaying orbit. The parabolic curve is the predicted decay from GR. (Ref. ?)

Two more post-Keplerian parameters, denoted by and , relate to the Shapiro delay suffered by the pulsar signal while passing through the curved spacetime surrounding the companion star. The relations for them in GR are as follows:


where is the projected semi-major axis of the pulsar orbit (in time units), a Keplerian parameter. The Shapiro delay is generally only observable when the orbital inclination is relatively close to 90, that is, the orbit is seen close to edge-on. For the Hulse-Taylor binary pulsar the orbital inclination is about 47and the Shapiro delay is small and covariant with Keplerian parameters. However it has been observed in a number of other neutron-star binary systems as will be discussed below.

Another relativistic effect, geodetic precession, has observable consequences for PSR B1913+16 and several other binary pulsars. A neutron star formed in a supernova explosion receives a “kick” during or shortly after the explosion which typically gives the star a velocity of several hundred km s . If the pulsar is a member of a binary system which is not disrupted by the kick, the orbital axis is changed so that pulsar spin axis, which was most likely aligned with the orbital axis prior to the explosion, is no longer aligned. Since kick velocities are often comparable to or even larger than the orbital velocities, the misalignment angle can be large. Precession of the spin axis will therefore alter the aspect of the radio beam seen by an observer and may even move the beam out of the line of sight.

In GR, the precessional angular frequency is given by


For PSR B1913+16, the corresponding precessional period is about 297 years, so the aspect changes over observational data spans are small. Never-the-less changes in the relative amplitude of the two peaks in the PSR B1913+16 profile were reported by Weisberg et al. and attributed to the effects of geodetic precession with a “patchy” beam. For a basically conal beam geometry, the separation of the two profile components would be expected to change and evidence for this was first found by Kramer leading to an estimate of the misalignment angle of about 22. A data set extending to 2001 was analysed by Weisberg and Taylor, obtaining results similar to those of Kramer. Their best-fitting model gives the “peanut” shaped beam shown in Figure 5. Clifton and Weisberg have shown that a set of circular nested emission cones can also give apparent pulse-width variations similar to those observed. Over the 20-year interval covered by the data set, the “impact parameter” (minimum angle between the beam symmetry axis and the observer’s line of sight) has changed by a rather small amount, about 3, from to . Extrapolation of this model suggests that the pulsar will become unobservable in about 2025. While this result is compatible with relativistic precession, it is not possible to derive an independent measure of the precessional rate from these data.

Fig. 5: Contours of the symmetric part of the radio emission beam for PSR B1913+16, obtained by fitting to the variation in pulse width over the interval 1981 to 2001. (Ref. ?)

2.1.2 Psr B1534+12

PSR B1534+12, discovered by Wolszczan in 1990, is a binary system with parameters quite similar to those of B1913+16, notably a short orbital period ( hours), relatively high eccentricity () and a compact orbit about 60% larger than that of PSR B1913+16. Analysis of less than a year’s data already allowed measurement of two post-Keplerian parameters, and , thereby determining the masses of the two stars and confirming that they are both neutron stars. These results also showed that the orbit was more edge-on than that of PSR B1913+16, with an implied inclination angle of about 77. Analysis of timing data extending over 22 years by Fonseca et al. built on earlier results by Stairs et al., with significant detections of and , the Shapiro-delay parameters, and the orbital decay, , giving five post-Keplerian parameters and three independent tests of GR. A fourth test of GR, albeit less precise, was provided by an analysis of the evolution of the profile shape and polarisation, yielding a rate for the geodetic precession of the pulsar spin axis  yr which is consistent with the value predicted by GR. Figure 6 shows the so-called “mass-mass” diagram for PSR B1534+12, illustrating these constraints.

If GR gives a correct description of the post-Keplerian parameters, all of these constraints should be consistent with an allowed range of and . For PSR B1534+12, the masses are most accurately constrained by and but the predicted constraint on appears to be inconsistent. As for PSR B1913+16, the measured value of must be corrected for kinematic effects resulting from the differential acceleration of the binary and solar system in the Galaxy. PSR B1534+12 is much closer to the Sun that PSR B1913+16 and the so-called “Shklovskii” effect, an apparent radial acceleration resulting from transverse motion of the binary system (, where is the pulsar proper motion and is the pulsar distance), is also important. The distance estimate is based on the pulsar DM and a model for the free electron density in the Galaxy and is quite uncertain. Stairs and her colleagues inverted the problem, assuming that the GR prediction for is correct, thereby deriving an improved value for the pulsar distance, a technique first suggested by Bell and Bailes .

Fig. 6: Plot of companion mass () versus pulsar mass () for PSR B1534+12. Constraints derived from the measured post-Keplerian parameters assuming GR are shown pairs of lines with separation indicating the uncertainty range. For and , the uncertainty ranges are too small to see on the plot. In addition, a constraint from the measured spin precession rate is shown. (Ref. ?)

2.1.3 The Double Pulsar, PSR J07373039a/b

The discovery of the Double Pulsar system heralded a remarkable era for investigation of relativistic effects in double-neutron-star systems. In this system, the A pulsar was formed first and subsequently spun up to approximately its current period of 23 ms by accretion from its evolving binary companion. The companion then imploded to form the B pulsar which has since spun down to its current period of about 2.8 s. The orbital period is only 2.5 hours, less than a third of that for PSR B1913+16, and the projected semi-major axis is about 60% that of PSR B1913+16. These parameters imply relativistic effects much larger than those seen for PSR B1913+16; for example, the predicted relativistic periastron advance is 169 yr, more than four times the value for PSR B1913+16. Added to that, the orbit is seen within a degree or so of edge-on, not only allowing detailed measurement of the Shapiro delay, but also resulting in eclipses of the A pulsar emission by the magnetosphere of the B pulsar. Finally, the still-unique detection of the second neutron star as a pulsar allows a direct measurement of the mass ratio of the two stars from the ratio of the two non-relativistic Roemer delays ().eeeThe B pulsar became undetectable as a radio pulsar in 2008, most probably because the beam precessed out of our line of sight. Because of uncertainty about the beam shape, the date of its return to visibility is very uncertain, but it should be before 2035.

Four post-Keplerian parameters (, , and ) were detected in just seven months of timing data from the Parkes 64-m and Lovell 76-m Telescopes . Further observations, including data from the Green Bank 100-m Telescope, with a 2.5-year data span give the currently most stringent test of GR in strong-field conditions . Figure 8 shows the mass-mass diagram based on these results together with the measurement of geodetic spin precession for the B pulsar described below. A total of six post-Keplerian parameters together with the mass ratio gives five independent tests of GR. As well as the three post-Keplerian parameters, , and observed for PSRs B1913+16 and PSR B1534+12 (Figure 6), for the Double Pulsar we have the mass ratio , the Shapiro delay parameters and and a measurement of the rate of geodetic precession . The observed Shapiro delay (Figure 7) shows that the J07373039A/B orbit inclination angle is .

Fig. 7: Observed Shapiro delay as a function of orbital phase for PSR J07373039A. Upper panel: timing residuals after fitting for all parameters except the Shapiro-delay terms, and , with these set to zero. Lower panel: the full Shapiro delay obtained by taking the best-fit values from a full solution, but setting and to zero. The red line is the prediction based on GR. (Ref. ?)

As mentioned above, this nearly edge-on view of the orbit results in eclipses of of the A-pulsar emission by the magnetosphere of the B pulsar. These eclipses last only about 30 s, showing that the B-pulsar magnetosphere is highly modified by the relativistic wind from pulsar A . Remarkably, observations with high time resolution made with the Green Bank Telescope showed that the eclipse is modulated at the spin period of pulsar B . Modelling of the detailed eclipse profile by absorption in the doughnut-shaped closed-field-line region of the magnetosphere by Lyutikov and Thompson allowed determination of the geometry of the binary system including the offset between the B-pulsar spin axis and the orbital angular momentum axis, the so-called “misalignment angle” which they estimate to be about 60. Even more remarkably, detailed measurements of the eclipse profile over about four years enabled Breton et al. to directly estimate the rate of geodetic precession as yr, consistent with the predicted precessional period of 70.95 years based on GR (Equation 2.1.1). It is notable that no secular profile evolution has been observed for PSR J07373039A, implying that, unlike for the B pulsar, the misalignment angle for pulsar A is very small .

The most precisely measured post-Keplerian parameter is . This constraint is nearly orthogonal to that from the mass ratio (Figure 8) giving values for the two neutron-star masses of  M and  M . Together with the accurately measured Keplerian parameters, these two masses can be used to predict the remaining five post-Keplerian parameters using Equations As Figure 8 shows, GR gives a self-consistent description of the orbital motions, with all five independent constraints consistent with the masses derived from and .

The most stringent test comes from the derived value of , . The ratio of the observed and predicted values of is , by far the most constraining test of GR in the strong-field regime. Furthermore, this test is qualitatively different to that for PSR B1913+16 in that it is based on a non-radiative prediction. With just 2.5 years of data analysed, the orbital decay term for the Double Pulsar is not measured as precisely as that for PSR B1913+16. However, since the phase offset grows quadratically (Figure 4) and the effect of receiver noise decreases as where is the data span (assuming approximately uniform sampling), the precision of the measurement should improve as . Furthermore, J07373039A/B is much closer to the Sun than PSR B1913+16, thus reducing the magnitude and uncertainty of the correction due to differential Galactic acceleration. Also, using very long-baseline interferometry (VLBI), Deller et al. have shown that the transverse space velocity of the binary system is small, just  km s, and so the Shklovskii correction to the observed is also small and accurately known. For these reasons, the orbit-decay test with the Double Pulsar will not be limited by these corrections, at least for another decade or so. Beyond that, improved models for the Galactic gravitational potential can be expected from analysis of GAIA data and improved parallax and proper motion measurements may come from further VLBI observations, both reducing the uncertainty in these kinematic corrections.

Fig. 8: Plot of companion mass () versus pulsar mass () for PSR J07373039 with observed constraints interpreted in the framework of GR. The inset shows the central region at an expanded scale, illustrating that GR is consistent with all constraints. (M. Kramer, private communication)

Several other post-Keplerian parameters can in principle be observed in binary systems , testing other aspects of relativity. Kramer and Wex show that only one of these, , which quantifies relativistic deformations of the binary orbit, is potentially measurable with the Double Pulsar system. Even this will take more than a decade to reach an interesting level of precision given reasonable expectations about future observations.

A potentially exciting prospect is to use measurement of higher-order terms for relativistic periastron precession to put constraints on the moment of inertia of a neutron star . At the second post-Newtonian level, the periastron precession has three components:


where is given by Equation 2.1.1, is the 2nd post-Newtonian contribution and gives the contribution from spin-orbit coupling . This latter term depends on the moment of inertia and spin of the A pulsar (the spin angular momentum of the B pulsar is negligible) and can in principle be measured if the two leading components in Equation 2.1.3 can be measured to sufficient precision. Kramer and Wex show that, given expected advances in precision timing and astrometry, a significant result could be obtained in 20 years or so.

Spin-orbit coupling can also lead to a nonlinear variation in as a function of time and a time variation in the projected semi-major axis of pulsar A, . However, the nonlinear terms depend on , where is the misalignment angle. As mentioned above, it seems as though is very small, so unfortunately these terms may be difficult to detect.

2.1.4 Measured post-Keplerian Parameters

Table 1 summarises the measured post-Keplerian parameters for the eleven pulsar binary systems where three or more post-Keplerian parameters have been measured. Since only two measurements are required to determine the stellar masses, these systems provide the opportunity for tests of theories of relativistic gravity. Of the eleven systems, in six cases the companion star is believed to be a neutron star, in a further four cases a white dwarf companion is more likely and in one case the companion is probably a main sequence star. This latter system, PSR J1903+0327, evidently has had an unusual evolutionary history as a triple system in which one component was ejected . It is currently not clear if there is a significant contribution to the observed from kinematic effects and/or spin-orbit coupling, so the utility of this system for tests of relativistic gravity is limited. There a further ten pulsar systems in the literature where just two post-Keplerian parameters have been measured, enabling estimates of the stellar masses, but no tests of relativistic gravity.

Although an independent measurement of relativistic spin precession has been possible for just two systems so far, as shown in Table 1, secular changes in observed pulse profiles have been observed for four other pulsars and attributed to geodetic precession of the pulsar spin axis.

  Pulsar/Parameter J04374715 J07373039A/B J0751+1807 J1141-6545 B1534+12 J17562251 Peri. advance (yr) 0.016(8) 16.8995(7) 5.3096(4) 1.7557950(19) 2.58240(4) Time dilation (ms) 0.3856(26) 0.773(11) 2.0708(5) 0.001148(9) Orb.P deriv. () 3.73(6) 0.6746(28) 0.99974(39) 0.90(5) 0.9772(16) 0.93(4) Comp. mass () 0.254(18) 1.2489(7) 0.191(15) 1.35(5) 1.6(6) Geod. prec. (yr) 4.77(66) Note d 0.59(10) Binary companion He WD NS He WD CO WD NS NS References ? ?, ? ?, ? ? ? ? Pulsar/Parameter J1807-2459B J1903+0327 J1906+0746 B1913+16 B2127+11C Peri. advance (yr) 0.018339(4) 0.0002400(2) 7.5841(5) 4.226598(5) 4.4644(1) Time dilation (ms) 26(14) 0.470(5) 4.2992(8) 4.78(4) Orb.P deriv. () 0.99715(20) 0.9759(16) Comp. mass () 1.02(17) 1.03(3) Geod. prec. (yr) Note d Note d Note d Binary companion CO WD(?) MS NS(?) NS NS References ? ? ?, ? ?, ? ?  

Note: a: Binary companion types: CO WD: Carbon-Oxygen White Dwarf; He WD: Helium White Dwarf; MS: Main-sequence star; NS: Neutron star
b: Dominated or significantly biased by kinematic effects
c: For PSR J07373039B
d: Effects of precession observed, but no independent determination of

Table 1: Binary pulsars with three or more significant measured post-Keplerian parameters

2.2 Tests of Equivalence Principles and Alternative Theories of Gravitation

Pulsars and especially binary pulsars have unique advantages in testing theories of relativistic gravitation as a result of their often rapid spin, short orbital periods and the ultra-high density of the under-lying neutron stars. As we have shown above, GR has been amazingly successful in describing all measurements to date. Never-the-less, investigations of quantum gravity and cosmology suggest that, in some regimes, extensions or modifications of GR may be required. This strongly motivates a search for departures from GR within existing experimental capabilities.

Gravitational theories have Equivalence Principles at their heart. The Weak Equivalence Principle (WEP) is basic to Newtonian gravity, stating that acceleration in a gravitational field is independent of mass or composition. The Strong Equivalence Principle (SEP) adds Lorentz invariance (no preferred reference frame) and position invariance (no preferred location) for both gravitational and non-gravitational interactions. GR satisfies the SEP whereas other theories may violate the SEP or even the WEP in one or more respects.

Comparison of different gravitational theories has been greatly facilitated by the “parametrized post-Newtonian” (PPN) formalism which describes observable or potentially observable phenomena in a theory-independent way. This formalism was first developed by Will and Nordtvedt for “weak-field” situations, that is, where , where is the Newtonian gravitational constant, and characterise the size and mass of the object and is the velocity of light. Many gravitational tests are performed within the solar system where , firmly in the weak-field regime. However in the vicinity of a neutron star and so strong-field effects are potentially important. A number of authors have considered the generalisation of the PPN formalism to strong-field situations (see, e.g., Refs ?, ?, ?) allowing investigation of these effects.

Some Lorentz-violating theories predict a dependence of the velocity of light on photon energy or polarisation . Pulsar observations can be used to limit these theories, but potentially stronger limits come from observations of gamma-ray bursts, polarised extra-galactic radio sources and the cosmic microwave background .

A recent comprehensive review of observational limits on theories of gravitation by C. M. Will can be found in Ref. ?. Pulsar tests of relativistic gravitation have been reviewed by I. H. Stairs and, more recently, by N. Wex. Further details on many of the topics discussed here may be found in these reviews.

2.2.1 Limits on PPN parameters

The standard PPN formalism has ten parameters: , , , , , , , , and . In GR , describing space curvature per unit mass, and , describing superposition of gravitational fields, are unity and all others are zero. describes preferred location effects, the preferred frame effects and the others describe violations of momentum conservation. ( also may be non-zero in this case.) Pulsar observations do not directly constrain or but are important in constraining many of the remaining PPN parameters and in fact currently place the strongest constraints on several parameters.

Damour and Schäfer recognised that wide-orbit low-eccentricity binary pulsars, which generally have a white-dwarf companion, could provide a valuable test of the SEP through a strong-field extension of the solar-system tests pioneered by Nordtvedt. A violation of the SEP would cause bodies with different gravitational self-energy to fall at different rates in an external gravitational field. In a pulsar – white-dwarf binary system, this results in a forced eccentricity in the direction of the gravitational field, that of the Galaxy in this case. This eccentricity is given by


where is the projection of the Galactic gravitational field on to the orbital plane, is the relativistic periastron advance, is the orbit semi-major axis and the orbital period. The ratio of the gravitational mass and the inertial mass (which is exactly one if the SEP is obeyed) for body is described by


where is the Nordtvedt parameter, a function of several PPN parameters (see, e.g., Ref. ?) and . Violations of the SEP will result in non-zero .

Since Damour and Schäfer first proposed this method of testing SEP violations, the number of suitable pulsar – white-dwarf binary systems has increased greatly. Gonzalez et al. combined data for 27 systems to place a 95% confidence upper limit on of .fffThis limit may be slightly under-estimated – see Wex. Since for a neutron star, this limit is not as strong as the weak-field limit on from lunar-laser ranging . However, it does enter the strong-field regime and test possible violations of the SEP that solar-system tests cannot reach.

The recent discovery by Ransom et al. of a remarkable triple system containing a pulsar, PSR J0337+1715, and two white-dwarf stars in essentially coplanar orbits, one in a relatively tight 1.6-day orbit with the pulsar and the other in a much wider 327-day orbit around the inner system, opens up the possibility of a much more sensitive test of the SEP. Precise timing observations of the pulsar have already shown that the motion of the inner system is strongly affected by the gravitational field of the outer white dwarf. The gravitational field of this star, which has a mass of about  M, at the inner system is at least six orders of magnitude larger than the Galactic gravitational field at the position of a typical pulsar and so the effect of SEP violations may be expected to be correspondingly larger. Observations over several orbital periods of the outer star will almost certainly be necessary to isolate any SEP-related effects from other orbit perturbations.

The wide-orbit low-eccentricity binary pulsars can also be used test for violations of local Lorentz invariance (LLI) of the gravitational interaction and momentum conservation that are described by the parameter and its strong-field generalisation . Bell and Damour show that such violations produce a forced eccentricity analogous to that produced by the Nordtvedt effect given by


where is a dimensionless “compactness” parameter, for neutron stars about 0.2 , and is the angle between , the velocity of the system with respect to a reference frame defined (for example) by the cosmic microwave background, and the pulsar’s spin axis. A similar analysis to that for the generalised Nordtvedt effect resulted in a 95% confidence limit of , some 13 orders of magnitude lower than the best solar-system test . It is worth noting that the observed small scatter in the period derivatives of MSPs had already been used by Bell to limit to .

Strong-field limits on the other two PPN parameters describing preferred frame effects, specifically LLI, and , can also be obtained from low-eccentricity binary pulsar observations . Non-zero induces a forced eccentricity in the direction of motion of the binary system velocity , analogous to Nordtvedt tests discussed above, whereas non-zero induces a precession of the orbital angular momentum about . The best current limits come from observations of the binary pulsar systems PSRs J1012+5307 and J1738+0333, both of which have short orbital periods, and  days respectively, and extremely low orbital eccentricities . Furthermore, both of these pulsars have optically identified binary companions which, together with proper motion measurements from the timing observations, allow the three-dimensional space velocity of the binary system to be determined. For these pulsars, the observational data span is sufficiently long ( years) that relativistic periastron advance has significantly changed the orientation of the intrinsic eccentricity vector relative to the direction of any forced eccentricity, resulting in a potentially detectable change in the total eccentricity. The strongest limit on comes from observations of PSR J1738+0333 as shown in Figure 9, conservatively . This is not only better than the best weak-field limit of from lunar laser ranging but also constrains strong-field effects as well.

Fig. 9: Constraints on the strong-field PPN parameter from timing observations of the low-eccentricity binary pulsar PSR J1738+0333. The constraints are a function of the unknown orientation of the binary system on the sky (described by the longitude of the ascending node ) and the unknown sign of . If the PPN parameter is assumed to be zero, then only certain ranges of are permitted. The shading corresponds to 68%, 95% and 99.7% confidence limits on . (Ref. ?)

A precession of the orbital axis about the system velocity vector induced by a non-zero value of is potentially observable as a time variation in the projected semi-major axis of the pulsar orbit . is one of the possible post-Keplerian parameters in standard timing solutions and significant values have been measured for both PSRs J1012+5307 and J1738+0333 . There are several possible contributions to the observed but in Ref. ? it is shown that all of these are negligible in these systems except that due to the changing orbit inclination resulting from proper motion of the system. This is a function of the unknown values and for certain values there is no constraint. Consequently only a probabilistic limit on can be obtained. By combining results for the two pulsars, a 95% confidence limit of is obtained.

This is not as constraining as a solar-system limit obtained from the present deviation of the Sun’s spin axis from the solar-system orbit normal , a limit that rests on the assumption that the two axes were aligned at the time of formation of the solar system.

However, pulsars provide an even stronger constraint based on the stability of the spin axis of isolated pulsars. Any precession of the spin axis of a pulsar is likely to result in changes in the observed pulse profile (as observed with geodetic precession in binary pulsars as discussed in §2 above). Shao et al. compared mean pulse profiles for the isolated MSPs B1937+21 and J17441134 taken 10 – 12 years apart with the same observing system and found no perceptible change in the pulse width at 50% of the peak amplitude. They interpreted these results by assuming a circular beam profile for the main pulse in PSR B1937+21 and for PSR J17441134. All of the angles in the problem can be estimated from modelling of radio and gamma-ray observations of the pulsar, taken together with known direction of the system velocity with respect to the cosmic microwave background, except the angle of the projected spin axis on the sky. Probability histograms for allowing for this unknown angle are shown in Figure 10. The final result for the 95% confidence upper limit is , about four orders of magnitude better than the limit described above based on orbital precession in pulsar binary systems and two orders of magnitude better than the limit based on solar spin precession. The assumption of circular beams in these MSPs is problematic, since there is good evidence for caustic enhancement in MSP pulse profiles, both radio and gamma-ray, which would tend to elongate the beam in the latitude direction, reducing the effect of precession on the observed pulse profiles. However, taking this into account would probably increase the limit by a factor less than ten, and so this limit would remain the best available.

Fig. 10: Constraints on the strong-field PPN parameter based on the long-term stability of the mean pulse profiles for the isolated MSPs B1937+21 and J17441134. (Ref. ?)

The stable pulse profiles of PSRs B1937+21 and J17441134 have also been used to limit the PPN parameter describing local position invariance (LPI), also known as the Whitehead parameter, and its strong-field counterpart (Ref. ?). The centripetal acceleration of Galactic rotation results in an anisotropy in the local gravitational field at a pulsar resulting in a precession of the pulsar spin vector around the direction of the Galactic acceleration with period.


where is the velocity of Galactic rotation at the pulsar and is the (unknown) angle between the pulsar spin vector and the Galactic acceleration. Combining the results for the two pulsars in a way analogous to that for the test described above, Shao and Wex obtain a 95% confidence limit . Even with the qualification about beam shapes mentioned above, this limit is at least two orders of magnitude better than the next best limit obtained from considering of the evolution of the solar spin-misalignment angle over the lifetime of the solar system .

The PPN parameter is one of a number of parameters that may be non-zero in gravitational theories that violate conservation of total momentum . A non-zero would result in an acceleration of the centre of mass of a binary system that changes direction with periastron precession. This is best measured by looking for a change in the rate of orbital decay in a binary system with large periastron precession and a long data span. Will used observations of PSR B1913+16 with a 15-year data span to obtain a limit . Obviously, this test could be made more stringent using the long data spans now available for PSR B1913+16 and other double-neutron-star systems with high .

2.2.2 Dipolar gravitational waves and the constancy of G

As described above, because of its tiny eccentricity, the pulsar – white-dwarf binary system PSR J1738+0333 has played a significant role in limiting the PPN parameters describing preferred frame effects. This system is composed of two very different stars, the neutron star we see as a pulsar and the companion which we know to be a white dwarf of mass 0.181 M through its optical identification . This and its short orbital period ( h) make it an ideal system for testing theories of gravity that predict a dipolar component to gravitational-wave damping as well as general scalar-tensor theories .

From analysis of about 10 years of timing data obtained at Parkes and Arecibo observatories, Freire et al. find an observed rate of orbital decay for PSR J1738+0333 of . To obtain the intrinsic rate of orbital decay, kinematic contributions from the differential acceleration of the binary system and the solar system in the Galactic gravitational field and the Shklovskii effect due to transverse motion must be subtracted, giving . Since the orbital parameters and the masses of the two stars are well known, the orbit decay due to GR can be accurately determined: , leaving a residual orbit decay of .

This residual orbit decay is consistent with zero, which can be interpreted as a further confirmation of the accuracy of GR. However, because of the very different nature of the two stars in this binary system, this result also places strong constraints on theories of gravity that predict a dipolar component to gravitational-wave emission. Besides a dipolar component , there are several other possible contributions to :


where is due to mass loss from the binary system, is a term resulting from tidal effects on the white dwarf (tidal effects on the neutron star are negligible) and is decay resulting from a possible variation in the gravitational ‘constant’ . Freire et al. show that the likely and tidal terms are small for this system, , and so the limit on is effectively a limit on .

Within certain restrictions on strong-field effects , the dipole term is given by:


where is the mass ratio (with subscript 1 refering to the neutron star), is the difference in “sensitivity” of the mass of each body to a scalar field , where


and is a body-independent constant that describes the dipole self-gravity contribution in a given theory of gravity (see, e.g., Ref. ?). The sensitivity depends on the stellar equation of state and, for neutron stars, is typically about 0.15, whereas for a white dwarf it is . Therefore if is non-zero, dipole radiation will contribute to the orbit decay.

The remaining term in the residual orbit decay is that due to possible variations in . In weak gravity, has been constrained to be less than  yr from lunar laser ranging experiments , giving


and hence . However, it is also possible to obtain independent estimates of the effects of dipole radiation and by combining the J1738+0333 results with those for PSR J04374715 . This southern pulsar is in a wider orbit than PSR J1738+0333 and hence has a different mix of the dipole and components, allowing them to be separated. After accounting for the fact that a changing will also change the stellar masses , the formal results are  yr and , both effectively upper limits. While the limit on is the best available, the derived limit on is about an order of magnitude weaker than the result (actually a limit on the variation of ) from the Mars Reconnaissance Orbiter and from lunar laser ranging .

Interestingly, pulsars have provided two other independent limits on . Thorsett used determinations of neutron star masses from timing observations of double-neutron-star systems that formed many gigayears ago. In standard formation scenarios, the mass of a neutron star depends on the Chandrasekhar mass, the maximum possible mass of a white dwarf star, just prior to the collapse to a neutron star. The Chandrasekhar mass is proportional to and so the observed small range of neutron star masses implies that . This limit has been somewhat weakened by recent discoveries of both less massive and more massive neutron stars in pulsar binary systems (see Ref. ? for a recent review).

The very small observed rate of change of pulsar period observed in some pulsars (after correction for kinematic effects) may be used to set a further independent limit . A variation in will result in an inverse variation in the stellar moment of inertia, with the exact relation depending on the neutron-star structure. If the observed (intrinsic) is entirely attributed to this effect, a limit of is obtained.

2.2.3 General scalar-tensor and scalar-vector-tensor theories

Many alternate theories of gravity can be expressed in a “tensor-scalar” framework in which a scalar field contributes to the “physical metric” through a coupling function :


where is the usual tensor metric. The coupling constant may be expressed in different ways , one of which is as an expansion around the asymptotic value of the scalar field :


(Ref. ?). For GR, . In the well-known example of a tensor-scalar theory, the Brans-Dicke theory, the scalar coupling is described by a single parameter . For , the Brans-Dicke theory approaches GR. For this theory and . In other theories both the linear and quadratic terms in Equation 2.2.3 (and higher-order terms) may be non-zero.

The various observational constraints on PPN and post-Keplerian parameters can be expressed as limits in the (,) space for tensor-scalar theories as shown in Figure 11 . The Cassini experiment placed a strong limit on the PPN parameter which translates to a limit on of about 0.003. Only the limits on dipole gravitational radiation from the asymmetric binary systems PSR J11416545 and PSR J1738+0333 rival the Cassini limit over most of the space. Since the precision of these measurements will increase with time, it seems likely that ultimately binary systems such as these will provide the strongest constraints on tensor-scalar theories.

PSR J0348+0432 and its white-dwarf companion form another asymmetric binary system, one that is distinguished by its very short orbital period (2.46 hours) and massive neutron star ( M). The high neutron-star mass is of interest since some scalar-tensor theories predict a strongly non-linear relationship between the strength of dipole GW emission and neutron-star mass or self-gravity (cf., Equation 2.2.2). For most of the parameter space of the class of theories discussed by Freire et al., the limits on effective scalar coupling from PSR J0348+0432 are currently not as strong as those from PSR J1738+0333 but, because of the high neutron-star mass, they place stronger limits on some other theories with a greater degree of non-linear coupling.

Fig. 11: Constraints on the scalar-field parameters in the (,) plane from various observational tests. Only the region below each line is allowed by the corresponding test. “SEP” refers to the test of the strong equivalence principle based on low-eccentricity pulsar – white-dwarf binary systems , “LLR” refers to lunar laser ranging results and Cassini to the “Shapiro delay” experienced by signals to and from the Cassini spacecraft (on its way to Saturn) as its line of sight passed close to the Sun . Other binary-system tests are labelled according to the pulsar concerned. In GR, both and are zero. (Ref. ?)

Bekenstein has proposed a relativistic generalisation of the so-called “MOND” theory of gravity that seeks to avoid the need for dark matter in galactic dynamics. The generalisation invokes an additional vector field and hence is known as a tensor-vector-scalar theory. Such theories relax some of the constraints on tensor-scalar theories, in particular, the dipole radiation constraints, and allow significantly larger values of . However, as shown by Freire et al., the binary pulsar results still significantly constrain theories of this type and in fact are more constraining than solar-system tests. With future observations, binary pulsar tests have the potential to make this class of theories untenable.

In another example of the use of pulsar observations, especially the limits on dipolar-GW radiation, to place limits on gravitational theories, Yagi et al. have strongly limited the allowed parameter space for the LLI-violating “Einstein-Æther” and the “Khronometric” theories.

2.3 Future Prospects

As described above, GR has provided an accurate description of all pulsar timing results obtained so far. However, continued refinement of existing methods and development of new tests is highly desirable. Continued pulsar timing measurements, especially with the advent of new and more sensitive observing facilities such as the 500-m Arecibo-type FAST radio telescope in China and the Square Kilometre Array (SKA) in South Africa and Australia will certainly improve on existing limits and enable new tests of gravitational theories. They may even demonstrate a failure of GR and hence a need for a modified or conceptually different theory of gravity. Conversely, if GR is assumed to be valid, results of astrophysical significance can be deduced from the observations. For example, as described in Section 2.1.3, observations of higher-order terms in the relativistic perturbations may enable a measurement of the moment of inertia of a neutron star.

These more sensitive radio telescopes can also be used to search for previously unknown pulsars and binary systems that are suitable for tests of gravitational theories. Past experience has shown that pulsar searches repeatedly turn up new classes of object. This potential is wonderfully illustrated by the recent discovery of the pulsar triple system PSR J0337+1715 which promises to provide a strong limit on violations of the Strong Equivalence Principle. A dream for such searches is the discovery of a pulsar in a close orbit around a black hole as this would offer much more stringent tests of gravity in the strong-field regime. Discovery of a pulsar orbiting the black hole at the centre of our Galaxy with an orbital period of a few months or less could even allow a test of the so-called “no-hair” theorem for black holes.

3 The Quest for Gravitational-Wave Detection

The direct detection of the gravitational waves (GWs) predicted by Einstein’s general theory of relativity and other relativistic theories of gravity is one of the major goals of current astrophysics. As described in §2 above, we have excellent evidence from the orbital decay of binary systems for the existence of GWs at the level predicted by GR, but up to now there has been no direct detection of the changing curvature of spacetime induced by a passing GW. This changing curvature induces a change in the proper distance between two test masses, described by the gravitational strain . The problem is that, for any likely source, is tiny. For example, the LIGO gravitational-wave detector hopes to detect the merger of two neutron stars at a distance of 100 Mpc for which , a change in the length of its 4-km arms of  m or of the diameter of a proton!

Pulsar timing can measure a change in the proper distance between the pulsar and the telescope. Systematic changes in timing residuals for a given pulsar reflect unmodelled changes in the effective time of emission, the pulsar position, the propagation path or the position of the telescope. With care, and with observations of multiple pulsars, residual delays due to changing proper distances can be isolated, effectively giving us a set of interferometers with baselines of  km! However, even in the best cases, we can only measure the interferometer “phase” to about 100 ns, so that the limiting strain is about . Unlike ground-based laser-interferometer systems, which are most sensitive to GW signals with frequencies around 100 Hz, pulsar timing systems are most sensitive to signals with frequencies around the inverse of the data span, typically a few nanohertz. Potential sources of detectable GWs in this low-frequency band include super-massive black-hole (SMBH) binary systems in distant galaxies and cosmic strings in the early Universe.

3.1 Pulsar Timing Arrays

The effect of GWs on pulsar timing signals was first considered by Sazhin and Detweiler, with the latter being the first to consider the effect of a GW from a distant source passing over a pulsar and the Earth. For this case, it can be shown that the net effect on the observed pulsar arrival times is simply the difference between the effect of the GW passing over the pulsar and the effect of the GW passing over the Earth. For a GW travelling in the direction, the redshift of the pulse frequency for a pulsar at distance with direction cosines (, , ) is given by:


where , with representing the two possible wave polarisation states (,) in GR, and and the gravitational strain at the pulsar and the Earth, respectively.gggSee Ref. ? for a rederivation of the Detweiler result. The observed timing residuals are then given by the integral of the redshift,


The coefficients multiplying the and terms are the “antenna patterns” for the two polarisations as illustrated in Figure 12. A pulsar located in the GW propagation direction has zero sensitivity to the GW since its direction cosines and are zero. Furthermore, despite the term in the denominator of Equation 3.1, the response for a pulsar exactly in the direction (the same direction as the GW source) is also zero. This comes about because of a cancellation of the term by the expansion of when (cf. Ref. ?).

Fig. 12: Effective “antenna pattern” for detection of a GW with pulsar timing. The wave is propagating in the direction and is assumed to have the + polarisation. The pattern for the polarisation is the same but rotated by about the axis. (Ref. ?)

A pulsar timing array (PTA) consists of a set of pulsars spread across the sky which have precise timing measurements over a long data span. The detection of GWs by PTAs depends on the correlated timing residuals for different pulsars given by the Earth term in Equation 3.1. GWs passing over the pulsars produce uncorrelated residuals because of both the retarded time and the different GW environment for each pulsar. Also the pulsars themselves have uncorrelated timing noise at some level, either intrinsic or resulting from uncorrected variations in interstellar delays. Because the expected GW strain is so weak, only MSPs have sufficient timing precision to make GW detection with PTAs feasible.

For an isolated source of continuous GWs, say an SMBH binary system in a nearby galaxy, in principle, both the Earth term and and the pulsar term in Equation 3.1 could be detected. For a rapidly evolving source, the pulsar term and the Earth term may be at different frequencies because of the retarded time of the pulsar term (see, e.g., Ref. ?). They can then be added incoherently to increase the detection sensitivity. For a non-evolving source, i.e., where is the change in binary orbital frequency over the pulsar timing data span , in principle the pulsar term and the Earth term could be summed coherently for optimal sensitivity. As is discussed further in Section 3.4 below, unfortunately we currently don’t know enough pulsar distances to sufficient accuracy to make this coherent addition possible. In a PTA, the pulsar terms therefore add with random phase, washing out the fringes in the antenna pattern (see also Ref. ?) and adding “self-noise” to the signal from the Earth term. Since the antenna pattern (Figure 12) has a maximum for pulsars roughly in the same direction as the GW source, the maximum response of a PTA is toward the greatest concentration of pulsars in the array.

A stochastic background of nanohertz GW from many SMBH binary systems in distant galaxies is likely to be the signal first detected by PTAs. To a first approximation, this background is also likely to be statistically isotropic, i.e., the expectation value is independent of direction when averaged over typical data spans. Hellings and Downs were the first to show that, in this case, the correlation between GW-induced timing residuals for two pulsars separated by an angle on the sky is dependent only on and not on the sky positions of the two pulsars. The zero-lag correlation function, commonly known as the Hellings & Downs curve and obtained by integrating the product of the antenna patterns (Figure 12) for the two pulsars over all possible GW propagation directions, is given by:


where , and is plotted in Figure 13. goes negative for angular separations around and then positive again for pulsars that are more-or-less opposite on the sky – this is a direct consequence of the quadrupolar nature of GWs. It is also important to note that the limiting value as is 0.5, not 1.0. This is a consequence of the fact that the pulsar terms in Equation 3.1 are uncorrelated and, on average, of equal amplitude to the Earth term. The scatter in the simulated correlations results from the random phases of the pulsar terms and illustrates the “self-noise” that limits the sensitivity of PTA experiments in the strong-signal limit.

Fig. 13: The Hellings & Downs correlation function, i.e., the correlation between timing residuals for pairs of pulsars as a function of their angular separation for an isotropic stochastic background of GWs. Also shown are simulated correlations between the 20 pulsars of the Parkes PTA for a single realisation of a strong GW signal that dominates all other noise contributions. (Ref. ?)

3.2 Nanohertz Gravitational-Wave Sources

3.2.1 Massive black-hole binary systems

There is good evidence that massive black holes form in the centre of galaxies at very early times (see, e.g., Ref. ?) and also that merger events play a major role in galaxy growth (see, e.g., Ref. ?). When two galaxies, each containing a central massive black hole, merge, dynamical friction will result in the two black holes migrating to the centre of the merged galaxy to form a binary system, with an estimated timescale for the migration of the order of giga-years (see, e.g., Ref. ?). When the binary separation is less than about 1 pc, loss of energy to GWs becomes the dominant orbital decay mechanism and the binary system will ultimately coalesce, becoming a strong GW source as it spirals in. There remains much controversy about the efficiency of orbital decay mechanisms as the binary separation approaches a parsec – known as the “last parsec problem”. Some (see, e.g., Ref. ?, ?) argue that dissipation mechanisms will quickly move the binary system through this phase, whereas others (see, e.g., Ref. ?) argue that the binary is likely to stall at separations where the gravitational decay is ineffective. Detection of a stochastic GW background (GWB) would resolve this issue.

Since large numbers of binary systems with different orbital periods contribute to the GWB, it is a broadband signal which is best described in the spectral domain. It is convenient to express the amplitude of the GW signal in terms of the dimensionless “characteristic strain”, defined by


where is the Fourier transform of and is the GW frequency. (Note that, for a circular binary system, the frequency of the emitted GW is twice the binary orbital frequency .) Two other quantities that are often used to parameterise GW spectra and detector spectral sensitivities are the square root of the one-sided strain power spectral density


and the GW energy density as a fraction of the closure energy density of the Universe


where is the Hubble constant.

In order to understand the astrophysical implications of results obtained from PTA experiments, it is necessary to have estimates of the likely strength of signals from potential sources of nanohertz GW. For a cosmological population of SMBH binary systems at luminosity distance and redshift , the local energy density in GW at frequency is given by:


where is the binary chirp mass and and are the masses of the binary components, is the comoving density of binary systems with redshift and chirp mass between and and and , respectively, and is the total energy emitted by a single binary system in the logarithmic frequency interval , where . The local GW energy density is related to the local characteristic strain by


and for a circular binary system


Therefore we have


As Phinney has emphasised, the result that for a cosmological population of circular binary systems decaying through GW emission is quite general and independent of any particular cosmology, black-hole mass function or galaxy merger scenario. Consequently the spectrum of the GWB is often parameterised as follows:


where , is the characteristic strain at and for the case described above. For pulsar timing experiments, the one-sided power spectrum of the timing residuals is given by


Consequently, a GWB produces a very “red” modulation of the timing residuals with a spectral index of for .

In order to estimate the likely strength of this modulation, the factor in Equation 3.2.1 must be evaluated. This requires a prescription for the cosmological evolution of massive black-hole binary systems in galaxies. Different approaches to this problem have been taken by different authors. An early paper by Jaffe and Backer used observational constraints on close galaxy pairs coupled with a black-hole mass function, whereas another early paper by Wyithe and Loeb used a prescription for merger of dark-matter halos coupled with different scenarios for growth of massive black holes in galaxies. The latter approach was developed further by Sesana et al. who showed that the GWB spectrum steepens at frequencies above about  Hz since the number of binary systems contributing to the background at these frequencies becomes small. This is illustrated in the left panel of Figure 14 which shows that binary systems at contribute most of the strain to the GWB. The right panel shows that massive binary systems with  M contribute most of the low-frequency GWB but that only systems with  M contribute to the high-frequency end. Importantly, at the high-mass end, only a few systems contribute to the GWB. As shown in Figure 15, this leads to a steepening and large uncertainties in the expected GWB spectrum at frequencies  Hz. These results were confirmed and extended by Sesana et al. and Ravi et al. using the Millenium simulation of the cosmological evolution of dark matter structures to define a merger history together with various prescriptions for galaxy and black-hole formation and growth.

Fig. 14: Left panel: Number of binary black-hole systems and their contribution to the characteristic strain of the GWB as function of source redshift . The solid lines are results from a Monte Carlo approach and the dotted lines are from a semi-analytic analysis. In both panels, the upper histograms are for a GW frequency  Hz and the lower histograms for  Hz. Right panel: The lower two panels are as for the left panel but as a function of source chirp mass . The upper panel shows the number of frequency bins spanned by the chirp over the 5-year span of the simulation. (Ref. ?)

While there is some consensus about the form of the GWB spectrum, there remain significant uncertainties. For example, Ravi et al. consider the effects of the stellar environment on the late evolution of massive black-hole binary systems in the cores of galaxies and conclude that the effect of dynamical friction is important, both in extracting energy from the binary system and inducing an eccentricity. Both of these have the effect of reducing the strength of the predicted GWB, especially at the low-frequency end, consequently making its detection by PTAs more difficult. On the other hand, McWilliams et al. argue for a model in which all black-hole growth is by merger rather than by accretion after coalescence, which is the main contributor to black-hole growth in the models discussed above. This leads to predictions of a significantly larger GWB characteristic strain compared to previous predictions and, hence, the imminent detection of the GWB by PTAs.

Fig. 15: Characteristic strain spectrum for the GWB based on different prescriptions for black-hole growth by accretion between mergers (solid, dashed and dotted lines) and different black-hole mass functions (different panels). (Ref. ?)

As Figure 15 indicates, there is a possibility that the nearby universe could contain a massive black-hole binary system with an orbital period of the order of a few years that actually dominates the nanohertz GW spectrum. This opens up the exciting prospect of the GW detection and study of an isolated supermassive black-hole binary system using pulsar timing and even the possibility of detection and study of the system in the electromagnetic bands – so-called “multi-messenger” astronomy. Searches for binary GW sources will be described in §3.3 and the prospects for their detailed study will be discussed in §3.4.

3.2.2 Cosmic strings and the early Universe

Cosmic strings and the related cosmic super-strings are one-dimensional topological defects which may have formed in phase transitions in the early Universe. Cosmic strings occur in standard field-theory inflation models, whereas superstrings are found in brane inflationary models. The idea that such strings will oscillate and hence emit GWs was first proposed by Vilenkin. Such oscillations may contribute to the stochastic GWB (see, e.g., Ref. ?) or generate bursts of GW radiation from string cusps and kinks (see, e.g., Ref. ?). The amplitude of GWs from cosmic strings is dependent on a large number of poorly known (or unknown) parameters and hence is very uncertain (see, e.g., Ref. ?). Key parameters are the string tension , usually parameterised by the dimensionless quantity , and the size of string loops relative to the horizon radius at the time of birth. Other significant parameters for the GWB are the intrinsic spectral index of the GW emission, a characteristic node number for the high-frequency cutoff in the emission spectrum and the probability of “intercommutation”, that is, intersecting strings dividing and the two parts exchanging. Such intercommutation can, for example, form two smaller loops from an intersecting twist in a larger loop. For standard strings, but it may be less for superstrings.

Vibrating cosmic strings are likely to decay by emission of GWs in a series of harmonics with fundamental frequency , where is the length of the loop, with an initial value , where is the horizon distance at the time of loop creation. The rate of energy loss for vibration mode of a given loop is:


where is factor depending on the shape of the loop, typically about 50 . As the loop loses energy, it shrinks and eventually disappears. The creation of loops through intercommutation and their decay through GW emission sets up an equilibrium distribution of loop sizes. Sanidas et al. compute the number density of loops as a function of loop length and time and hence, using Equation 3.2.2, the predicted spectrum of the GWB from string loops as a function of the various parameters. As Figure 16 shows, the spectrum is very broad, extending all the way from nanohertz to Megahertz. Since the lowest () frequency is proportional to the string length, the weaker and smaller loops do not contribute to the nanohertz background.

Fig. 16: Energy density spectrum for GWs emitted by cosmic strings as a function of the dimensionless string tension . Other parameters are held fixed at values , , and . The dashed spectrum is for the critical point where ; spectra above this are for large loops and spectra below are for small loops. (Ref. ?)

3.2.3 Transient or Burst GW sources

The prime target of ground-based laser-interferometer GW detectors is the burst emitted at the coalescence of a double-neutron-star system. Such a burst is intense for just a few milliseconds and clearly cannot be detected by PTA experiments which typically are sensitive to signals of duration between a few weeks and a few years. Possible sources of longer-duration bursts are coalescence of SMBH binary systems, highly eccentric massive black-hole binary systems, and formation and decay of cusps and kinks in cosmic strings.

A major difference between detection of GW bursts and continuous GW sources is that there is generally no interference between the Earth term and the pulsar terms in the signal detected by PTAs (Equation 3.1). This is because the duration of the burst (by definition, less than the observational data span) is much less than the light-time to the pulsars and so, except for a source in the same direction as a pulsar, the burst will occur at very different times in the Earth and pulsar terms. The pulsar term reflects the effect of the burst on the pulsar at a time before the burst arrives at the Earth, but the pulsar term is detected at a time later than the Earth term, where is the pulsar distance and is the angle between the pulsar and the GW propagation direction as seen from the Earth.

Figure 17 shows the GW waveform produced by the coalescence of two black holes in a coordinate system where all the signal is in . The maximum amplitude of the waveform is about 0.1 in the time units of Figure 17, or about or where is the total system mass in solar units and is the (comoving) distance in Gpc. Although this is comparable to the strain sensitivities achieved by current PTAs for continuous GW signals (see §3.3 below), even for the largest SMBHs, the timescale of the burst, is only of order 10 days and the period of the oscillation is about an order of magnitude less than that. Not only is this too short to be resolved by any existing PTA, but the sensitivity over this short interval would be much less than that achieved for CW signals integrated over the entire data span. Space-based laser interferometer systems such as the planned eLISA have the potential to directly detect these bursts.

Fig. 17: Gravitational waveform resulting from the coalescence of two equal-mass black holes (reduced mass ratio ) with total mass at distance . Both and are expressed in time units; the conversions to conventional units are and . and are assumed source directions. The dashed line is the predicted waveform if the gravitational memory effect is ignored. (Ref. ?)

However, Figure 17 also shows another effect known as GW “memory” which is potentially detectable with pulsar timing. During the coalescence event, a non-oscillatory component to the gravitational strain builds up, so that at the end of the “ring-down” phase, the strain has a permanent offset from the pre-coalescence value. The amplitude of the memory effect is


where is the mass fraction contributing to the memory effect and is the reduced mass fraction, 0.25 for equal-mass binary components .

This step change in produces a step change in the observed pulse frequency , i.e., a “glitch”. This glitch persists until it is reversed by the pulsar term at a time later. In a PTA, these reversals will occur at different times for different pulsars. Of course, it is also possible that a GW memory jump could be detected in the pulsar terms, but there it may be confused with a real glitch in the intrinsic pulse frequency, whereas in the Earth term there is a correlation in the effect on different pulsars. However, glitches in MSPs are rare (only one very small glitch detected so far: Ref. ?). Also, real pulsar glitches are generally spin-ups, whereas a GW-memory jump may be of either sign, depending on pulsar-source angle. Therefore, as Cordes and Jenet have discussed, the pulsar terms may give an improved probability of detection.

Black-hole binary systems with circular orbits emit GW at the second harmonic of the orbital frequency, i.e., . For eccentric orbits, the GW emission becomes more burst-like as the accelerations and hence GW power are greatest around periastron when the two black holes are closest together . In the spectral domain, power spreads to higher harmonics and also to the fundamental frequency . In the gravitational-decay phase of evolution, when energy loss is dominated by GW emission, the orbit tends to circularise . However, at earlier phases of the orbital decay when three-body stellar interactions or interaction with a gaseous disk surrounding the binary system are important, the eccentricity may grow. Stellar three-body interactions can result in orbital decay through dynamical friction but probably result in a modest increase in the eccentricity of the black-hole binary system (e.g., Ref. ?). However, Roedig et al. find that initially mildly eccentric binary systems decaying through interaction with a gaseous disk evolve toward a limiting eccentricity in the 0.6 – 0.8 range. In a regime of frequent mergers it is even possible that a black-hole triplet could form and, in this case, eccentricities as high as 0.99 could exist.

Finn and Lommen have investigated the GW emission and the resulting timing residuals from a close parabolic encounter of two massive black holes. As Figure 18 shows, an encounter of two black holes located at a distance of 15 Mpc with a minimum separation of 0.02 pc produces a GW burst of duration about 1 year and maximum GW strain . This results in potentially detectable PTA timing residuals of about s amplitude with the same timescale. Unfortunately, the probability of having such a close encounter of two SMBHs in the local universe within PTA observational data spans is not high.

Fig. 18: The upper plot shows the gravitational waveforms in the and polarisations resulting from a parabolic encounter with an impact parameter of 0.02 pc of two black holes located at a distance of 15 Mpc in the direction of the Virgo Cluster. The lower plot shows the resulting timing residuals for several PTA pulsars. (Ref. ?)

Cosmic strings are another potential source of GW burst emission. They can radiate over a wide frequency range with a huge range of possible amplitudes and timescales (cf., Figure 16) depending on the detailed mechanism invoked (loops, cusps, short strings, etc.) and the very wide (virtually unlimited) parameter space . Short bursts radiating in the LIGO and eLISA bands may be frequent and unresolved, producing a GWB at these frequencies. However, bursts with longer timescales, producing radiation in the nanohertz band, are also possible but are likely to be extremely rare . Consequently, while in principle such bursts could be detected, in practice it is unlikely that PTAs will be able to significantly constrain models for GW burst emission from cosmic strings.

3.3 Pulsar Timing Arrays and Current Results

In this section we first describe the three main PTAs currently operating world-wide: the European Pulsar Timing Array (EPTA), the North American pulsar timing array (NANOGrav) and the Parkes Pulsar Timing Array (PPTA), and the collaboration between them, the International Pulsar Timing Array (IPTA). PTAs have many possible applications such as establishing a pulsar-based timescale , investigating the accuracy of solar-system ephemerides , and investigating the properties of the pulsars themselves (e.g., Refs ?, ?) and of the intervening interstellar medium (e.g., Ref. ?). However, here we concentrate on what is undoubtedly their primary scientific goal, the direct detection of gravitational waves. Unfortunately, in common with other GW detection efforts around the world, PTAs have so far only been able to place limits on the strength of signals from potential GW sources. However, these limits are now beginning to seriously constrain the astrophysical source models and the assumptions that go into them and hence have implications that go far beyond the GW studies themselves.

3.3.1 Existing PTAs

The EPTA uses five large radio-telescopes in Europe, the Effelsberg 100-m telescope in Germany, the Nançay Radio Telescope in France (95-m equivalent area), the Westerbork Synthesis Radio Telescope in the Netherlands (similar effective area to the Nançay telescope), the 76-m Lovell Telescope at Jodrell Bank in England and the recently completed 64-m Sardinia Radio Telescope in Italy, to observe about 40 MSPs with a cadence of between a few days and 30 days for different pulsars . Different telescopes observe at different frequencies in the range 0.3 – 2.6 GHz, but all are instrumented at 1.4 GHz. Normally the five telescopes observe independently, but in a project known as the “Large European Array for Pulsars” (LEAP) 1.4 GHz signals over a bandwidth of 128 MHz from the five telescopes can be summed coherently to form a 194-m equivalent diameter radio telescope. The different telescopes use different signal-processing systems, either digital filterbanks or coherent dedisperion systems or both.

NANOGrav makes use of the 300-m Arecibo radio telescope in Puerto Rico and the 100-m Green Bank Telescope (GBT) in West Virginia . A sample of about 36 pulsars is observed, typically at 3-week intervals. At Arecibo, observations are made in bands centred at 430 MHz and 1410 MHz, whereas at the GBT, the observed bands are centred at 820 MHz and 1500 MHz. Currently both radio telescopes use coherent dedispersion systems with bandwidths up to 800 MHz, but in the past a range of filterbank and coherent dedispersion systems with more limited bandwidths have been used.

As the name suggests, the PPTA uses the Parkes 64-m radio telescope located in New South Wales, Australia. A sample of 22 pulsars is currently being observed with regular observations at 2 – 3 week intervals in three bands around 730 MHz, 1400 MHz and 3100 MHz respectively . Coherent dedispersion systems are used at 730 MHz and 1400 MHz with bandwidths up to 310 MHz and digital filterbanks at 1400 MHz and 3100 MHz with bandwidths of 256 MHz and 1024 MHz respectively.

Data sets from all three PTAs have spans ranging from a few years up to about 20 years for different pulsars; three (including the original MSP, PSR B1937+21) have Arecibo data spans of nearly 30 years. The three PTAs together observe about 50 pulsars with some being observed by two or even all three of the PTAs. Their distribution on the sky is shown in Figure 19.

Fig. 19: Distribution on the sky of MSPs being timed by the three PTAs, with different symbols for each PTA. Right ascension increases to the left with at the plot centre. The symbol size is related to the ratio , where and are the pulsar 1400 MHz flux density and pulse period respectively. (Ref. ?)

Given that the combined data set of the three PTAs contains a larger number of pulsars, improved observation cadence and greater frequency diversity than the data set of any one PTA, there is a strong motivation to combine all the available data sets to obtain maximum sensitivity for PTA scientific objectives. The International Pulsar Timing Array (IPTA) consortium was set up to facilitate progress toward this goal . The IPTA also arranges annual science meetings and student workshops and provides a forum for outreach programs and other activities related to PTA research.

3.3.2 Limits on the nanohertz GW background

As discussed in Section 3.2.1, GWs from a cosmological distribution of SMBH binary systems are expected to contribute a very “red” signal to the spectrum of pulsar timing residuals. The expected signal from other GWB sources is similar. Consequently, long-term observations of a single pulsar with little or no detectable intrinsic timing irregularities can be used to place a limit on the strength of the GWB in the Galaxy. Of course, statistical limits can be improved by using data from several such pulsars. An early limit on the GWB at a frequency of 4.5 nHz was set by Kaspi et al. using Arecibo observations of two MSPs, PSR B1855+09 and PSR B1937+21, with a 95% confidence limit on of , where  km s.

Since the advent of the various PTA projects, both the quality and quantity of timing data sets has improved and a variety of analysis techniques have been employed to extract increasingly restrictive limits. Based on early PPTA data on seven pulsars combined with the Kaspi et al. PSR B1855+09 data set, Jenet et al. used a “frequentist” approach with a statistic based on the amplitude of the low-frequency components in the power spectrum of the timing residuals to set a 95% confidence limit of about on at a GW frequency of 1/8 yr or 4 nHz. From Equations 3.2.1 and 3.2.1, this result is equivalent to a characteristic strain at frequency 1/1 yr, .

van Haasteren et al. analysed EPTA 1400-MHz data sets for five MSPs with spans of 5 – 8 years using a Bayesian analysis to place limits on the GWB amplitude as a function of its spectral index . For , the derived limit at the 95% confidence level is , about a factor 1.8 better than the Jenet et al. limit.

NANOGrav multi-band data sets recorded between 2005 and 2010 for 17 MSPs were analysed by Demorest et al.. Timing analyses taking into account time-varying dispersion delays and frequency-dependent pulse profiles were carried out to form sets of post-fit residuals and the corresponding covariance matrices for each pulsar. For most MSPs in the sample, no red noise signal was detectable in the post-fit residuals. Considering just the pulsar with the smallest post-fit residuals, PSR J1713+0747, and taking into account absorption of red noise by the timing fit, Demorest et al. obtained a 95% confidence limit for of . A separate cross-correlation analysis weighted by the expected Hellings & Downs function (Equation 3.1) across all the pulsars in the sample resulted in a somewhat better limit , although this limit was dominated by correlations with the two best pulsars in the sample, PSRs J1713+0747 and J19093744.

Based on PPTA and earlier Parkes timing observations made in three observing bands centred near 700 MHz, 1400 MHz and 3100 MHz respectively with data spans of up to 17 years , together with the PSR B1855+09 archival Arecibo data , Shannon et al. placed a limit of