Gravitational wave research using pulsar timing arrays

A pulsar timing array (PTA) refers to a program of regular, high-precision timing observations of a widely distributed array of millisecond pulsars. Here we review the status of the three primary PTA projects and the joint International Pulsar Timing Array project. We discuss current results related to ultra-low-frequency gravitational wave searches and highlight opportunities for the near future.


Introduction
Pulsar observations have been used for numerous astrophysical applications. Not long after the discovery of pulsars, Counselman & Shapiro (1968) [1] described how observations of pulsars could be used to "test general relativity, to study the solar corona, and to determine the earth's orbit and ephemeris time . . . [and to determine] the average interstellar electron density". Most studies to date have concentrated on analysing observations of specific pulsars. For instance, observations of one pulsar may provide excellent tests of general relativity, whereas another pulsar will be observed to probe the solar corona.
During 1982 the first millisecond pulsar was discovered [2]. A few hundred millisecond pulsars are now known. Their rotation is significantly more stable than the normal pulsars and their pulse arrival times can both be measured and also predicted with high accuracy. Foster & Backer (1990) [3] showed how a comparison of timing observations from multiple millisecond pulsars (a spatial array of pulsars) could be used to provide a time standard, to detect perturbations of the Earth's orbit and to search for gravitational waves (GWs). They initiated an observing program (which they termed a "pulsar timing array program") to observe three pulsars using the National Radio Astronomy Observatory 43 m telescope. During 2004, a much larger pulsar timing array (PTA) project began with the Parkes 64 m telescope and is known as the Parkes Pulsar Timing Array (PPTA). This project is ongoing (an overview and the first data release was described by Manchester et al. (2013) [4]) and now the team undertakes regular observations of 25 pulsars in three observing bands. The North American Nanohertz Observatory for Gravitational Waves (NANOGrav) in North America [5] was founded in 2007 and uses the Arecibo and Green Bank telescopes to observe 36 pulsars. Observations are also carried out for 42 pulsars with the Sardinian, Effelsberg, Nancay, Westerbork and Jodrell Bank telescopes by the European Pulsar Timing Array (EPTA) project team [6,7], which was also founded in 2007. The three project teams combine their expertise and data sets as part of the International Pulsar Timing Array (IPTA) [8,9].
In this review article, we will concentrate on one aspect of pulsar timing array research: searching for GW signals. The first observational evidence for GWs came from observations of a binary pulsar system (PSR B1913+16). That system was discovered by Hulse & Taylor (1975) [10] and was shown to be losing energy at exactly the rate predicted by the theory of general relativity for GW emission. The first direct detection of GWs was recently made by the LIGO/Virgo collaboration. During 2015 they detected two bursts of GW emission coming from the coalescence of stellar mass black holes [11]. This exciting result has opened the field of observational GW astronomy, and pulsar timing array projects provide a complementary view of the gravitational wave sky. Whereas the LIGO/Virgo detectors allow us to detect high-frequency GWs from stellar mass systems, the pulsar observations will allow us to detect ultra-low-frequency GWs from supermassive binary black holes. In contrast to the Hulse & Taylor (1975) work that provided evidence for GWs, the PTAs will enable a direct detection of GWs. For completeness we note that spacebased detectors (such as the Laser Interferometer Space Antenna; LISA; Amaro-Seoane et al. 2017 [12]) will be sensitive to GWs in a frequency range between the PTA experiments and LIGO.
In this paper we describe how GWs affect pulsar observations ( § 2). We give a brief summary of current data sets ( § 3) and the techniques being applied to hunt for the signals ( § 4). In § 5 we summarise the results obtained to date and in § 6 we highlight some of the current limitations of the data sets and techniques. In § 7 we consider the future of pulsar timing arrays and conclude in § 8.

How GWs affect pulsar observations?
Various authors [13,14] determined how GWs affect an electromagnetic signal propagating from an emitting object to a detector. Such calculations were applied to pulse arrival times from pulsars by Sazhin (1978) [15] and Detweiler (1979) [16]. A GW induces a fluctuation in the observed pulse frequency, δν/ν, of: where H i j is a geometrical term that depends on the position of the GW source, the Earth and the pulsar (at distance D from the Earth). The GW strain, h i j (t, x), is evaluated at the Earth (at time t e and position x e ) and at the pulsar (at the time the GW signal passed the pulsar, t p , and at position x p ). The shift of the pulse frequency is not directly measured. Instead pulse times-of-arrival (ToAs) are determined. These ToAs are then compared with predictions for the arrival times based on a pulsar timing model. The differences between the predictions and the measurements are known as the pulsar "timing residuals". A GW signal will induce timing residuals at time t from the initial observation of The theory of general relativity predicts two polarisation states, A + and A × , for GWs (see Lee, Jenet & Price 2008 [17] for non-GR predictions). We can therefore re-write the Earth term as (note that the pulsar term is the same, but with an extra phase): in which P + and P × are geometrical terms and γ is the GW-Earth-pulsar angle. For a non-evolving, continuous wave source (i.e., from a non-evolving, supermassive, binary black hole system), the A +,× will oscillate with an angular frequency of the GWs being ω g . For a supermassive, circular, binary, black hole system the GWs will be emitted at twice the orbital frequency. Eccentric binaries radiate GWs over a spectrum of harmonics of the orbital frequency. An estimation of the amplitude of the induced timing residuals caused by a binary system can be determined from [18]: where d is the luminosity distance to the system which has a total mass of M/(1 + z) (where z is the redshift) and the GW frequency is f . More details are provided in Rosado et al. (2016) [19] who considered the detectability of binary systems at high redshift. They showed that very high mass (> 10 10 M ⊙ ) binary systems could be detected by current PTAs at arbitrarily high redshifts. Of course, our universe will contain a large number of supermassive, binary black hole systems. To determine the total GW signal from these systems we therefore need to sum Equation 1 over all the sources. This results in a background of GW signals. For an isotropic, stochastic, unpolarised background signal, Hellings & Downs (1983) [20] showed that the timing residuals for each pulsar pair will be correlated as: where x = [1 − cos θ]/2 for an angle θ on the sky between two pulsars and δ(x) is the Dirac delta function 1 . This analytic expression is plotted in Figure 1 and is commonly referred to as the Hellings-and-Downs curve 2 . When searching for a background of GWs, the PTA teams therefore determine how correlated the timing residuals are for each pulsar pair. A convincing detection of the GW background will be made if those correlations are shown to follow the Hellings-and-Downs curve. Typically pulsars are observed every few weeks and the longest data spans are now a few decades (millisecond pulsars were discovered in 1982). This implies that PTA data sets are sensitive to GWs with wavelengths from weeks to years. These correspond to ultra-low-frequency (10 −9 -10 −8 Hz) GWs 3 .

Observations and timing data sets
The three PTAs carry out regular timing observations of their sample of millisecond pulsars. Details of the observing systems have been presented in the various data release papers (see Desvignes et al. 2016 [7], Arzoumanian et al. 2015 [25] and Manchester et al. 2013). In brief, the data from a given telescope is generally folded online using the known timing model for the pulsar being observed. The resulting data files are processed to remove radio-frequencyinterference and to apply various calibration procedures (such as removing instrumental delays and calibrating the polarisation and flux density of the signal). Pulse ToAs are determined for each observation by cross-correlating the observed pulse profile with a template providing a high S/N representation of the expected pulse shape. One of the primary noise sources that affect searches for GWs are variations in electron densities in the interstellar medium (see, for example, Keith et al. 2013 [26] and Lee et al. 2014 [27]). Such changes can be monitored and (to some extent) removed or modelled by observing the pulsars over a wide range of frequencies. The PPTA currently uses a dual band receiver providing simultaneous observations in the 10 cm (3 GHz) and 40 cm (700 MHz) observing bands along with a 20 cm receiver (1400 MHz). The EPTA uses their large number of telescopes to obtain observations of each pulsar at frequencies between ∼300 MHz with the Westerbork Synthesis Radio telescope and ∼2.6 GHz with the Effelsberg radio telescope. Data in the 20 cm observing band from five of the European telescopes are also combined as part of the Large European Array for Pulsars (LEAP) to form a tied-array telescope with an effective aperture equivalent to a 195 m diameter telescope [28]. NANOGrav carries out observations between ∼300 MHz and 2.4 GHz.
Some of the pulsars in the Southern hemisphere can only be observed by the PPTA. The Northern hemisphere pulsars are generally observed by a large number of telescopes in both Europe and North America. A few pulsars are observed by all three PTAs. This has led to some observing campaigns in which a large number of IPTA telescopes observe the same source. For instance, PSR J1713+0747 was observed non-stop for 24 hours (Dolch et al. 2016). A list of publically-accessible timing array data sets is given in Table 1. It is expected that new data sets will continue to be released from the IPTA and the three PTA projects on a regular basis.
Even though they are not official members of the IPTA, other telescopes are used to observe millisecond pulsars and are likely to contribute to PTA research. In China, the Nanshan, Yunnan, Shanghai and Jiamusi telescopes observe pulsars at a wide range of observing frequencies. The GMRT in India, LOFAR in Europe and the MWA in Australia observe pulsars at low frequencies (10 to 240 MHz and 80 to 300 MHz respectively). Papers related to PTA research have also been published using observations from Kalyazin in Russia (e.g., Ilyasov et al. 2004[30]).  [29] 4 Techniques GWs will induce timing residuals. The form of those timing residuals will depend on the nature of the GWs (single, non-evolving sources will induce sinusoidal residuals, a background will induce timing residuals that have a powerlaw spectrum). The statistical challenge is therefore to first search for statistically significant residuals and then to prove that they arise because of GWs. If no GWs are detected then upper bounds on the GW amplitude can be determined.
All existing techniques are based on the "pulsar timing method". Timing software packages (such as tempo, tempo2 or pint 4 ) are used to compare the measured pulse ToAs with predictions for those ToAs based on a model for the astrometric, pulse and, interstellar medium and orbital parameters of each pulsar. Both "frequentist" and "Bayesian" methods exist for searching GWs. We summarise these methods in Table 2. The various algorithms described in the table can be split into routines for specific types of GW signals (such as the findCW and detectGWB plugins to tempo2) or more general codes that can search for various GW types. In all these cases, the user provides a set of high precision pulsar observations and defines the type of GW signal to be searched for along with information on what other noise processes are likely to be present in the data. For instance, the measured timing residuals are not only induced by GWs as pulse ToAs are affected by many phenomena. Along the line-of-sight to the pulsar the interstellar medium and the Solar wind can contribute significant delays to the measured ToAs. The measurement of a ToA is also sensitive to instrumental errors and incomplete polarisation calibration. A GW background produces low-frequency timing residuals, but so do errors in terrestrial time standards, intrinsic pulsar instabilities and much more. A detection of a GW signal therefore requires the confirmation of the expected spatial angular signature (the Hellings-Downs curve for a background or a quadrupolar spatial signature for a single GW source). When searching for GWs it is therefore necessary either to first determine and then remove the non-GW noise processes or to search for the GWs whilst simultaneously modelling these other phenomena that affect the pulse arrival times. Bayesian algorithms can be significantly slower than frequentist algorithms and GW searches using large data sets often require high-performance-computing facilities 5 . Frequentist-based methods often require that each noise process is dealt with in turn and, without care, this can lead to correlations between the different processes being unaccounted for in the final results. All the available algorithms should be used with care and tested on simulated data sets that have similar properties to the real data (i.e., different noise properties, irregular sampling, different data spans for different pulsars, etc.). The IPTA produced a set of simulated data sets 6 with the primary goal of testing and comparing different GW detection algorithms.

Results to date
Ultra-low-frequency GWs have not yet been detected. Work has therefore been split between 1) predicting the expected signal and time to detection, 2) making more-and-more sensitive searches for the GWs and 3) understanding the astrophysical implications of our non-detections. A summary is given in Figure 2 where we show the current upper-bounds from the three PTAs as dotted lines. A theoretical prediction (from Sesana et al. 2016) for the likely GW background signal from coalescing supermassive, binary black holes is shown in the shaded region. One possible realisation of such a background (made up from numerous individual black hole binaries) is shown as the jagged, solid line. We note that the current PTAs are starting to constrain some models of the GW signal. However, it is likely  that we will require the sensitivity of the entire IPTA project to finally detect the GWs. A bound that potentially could be reached by the IPTA within a decade or so is shown as the dotted, black line. A more in-depth analysis of the likely time to detection has recently been published by Kelley et al. (2017).
In the following subsections we describe the results to date for different types of GW searches. We note (and see e.g., Rosado et al. 2015 [39]) that it is not yet clear whether we will first detect an individual supermassive binary black hole system or a background made up from a large number of GW sources. Lommen & Backer (2001) [40] unsuccessfully searched for GW emission from Sagittarius A * that had been proposed to be a binary system. They also searched for continuous GW emission from other nearby galaxies. Sudou et al. (2003) [41] presented evidence for a supermassive black hole binary system in the radio galaxy 3C66B. It was shown by  [44] have presented recent bounds on individual GW sources over the entire sky or in particular sky directions with a representative value being that the GW strain amplitude:

Individual supermassive binary black holes
at a frequency of 10 nHz. The exact bound depends upon the sky direction and GW frequency. Such bounds constrain the local merger rate density of supermassive binary black holes, but the current limits are higher than the theoretical expectations (e.g., Sesana 2013 [45] and Ravi et al. 2012 [23]) for such GWs.

Burst events
Even though the first GW detection was a burst signal in the audio band (Abbott et al. 2016), the expectations for burst GW signals are not well developed in the pulsar timing band. To date, the PTA teams have considered burst emission from the formation of supermassive black holes, highly eccentric black hole binaries, close encounters of massive objects, cosmic string cusps and memory events. Seto (2009) [46] showed how permanent distortions can occur in spacetime during mergers of supermassive binary black holes. Such a "memory event" will lead to a step change in the pulse frequency of all pulsars. Cordes & Jenet (2012) [47] gave an order-of-magnitude estimation of the signal-strength of: where µ is the reduced mass of the system and D is the distance. Searches for bursts with memory have been carried out by Arzoumanian et al. (2015) and Wang et al. (2015). No detection was made and the current conclusion is that GW memory events are unlikely to be detected in the near future.

Background of GWs
Most PTA research has been related to detecting or bounding the GW background. It is common to define the background spectrum as: where f 1yr = 1/(1yr). The spectral exponent is thought to be close to α = −2/3 for a background caused by coalescing black holes, −1 for cosmic strings and −7/6 for relic GWs (see Jenet [8]). The amplitude of the astrophysical GW background predicted by different models is based on various physical assumptions (see, e.g., Sesana 2013 for an overview). Firstly, models for the binary supermassive black hole (SMBH) population rely on measurements of the galaxy merger rate. Secondly, it has often been assumed that all galaxy mergers form binary SMBHs that coalesce well before a subsequent galaxy merger. Thirdly, the binary orbital decay has been assumed to be driven only by losses of energy to GWs when radiating in the pulsar timing frequency band.
The current GW background limits are inconsistent with some of the early theoretical models for the expected GW background signal, which suggests that at least one of the physical assumptions underlying such GWB models is likely to be incorrect. For instance, the galaxy-merger timescales might be longer than we currently expect. This would result in a lower merger rate and, hence, fewer binary SMBHs. It is also possible that not all large galaxies host SMBHs. SMBH binaries may not efficiently reach the GW emitting stage in our assumed time scale (they "stall"; e.g., Simon & Burke-Spolaor 2016 [52]) and other mechanisms in addition to GW emission drive binaries to coalescence, such as the coupling of binary SMBHs to their environments (e.g., Kocsis & Sesana 2011 [53]) or orbital eccentricity [54,55]). Binary SMBHs could lose energy and momentum because of three-body scattering of stars and viscous friction against circumbinary gaseous disk. Such coupling mechanisms also bend the GW background spectrum at low frequencies (i.e., large orbital separations). Arzoumanian et al. (2016) checked the consistency of their limit with previously reported scaling relations between SMBH mass and galactic bulge mass, using fiducial estimates for galaxy merger rates and the stellar mass function. Under the assumption of circular GW-driven binaries, they found that the scaling relations of Kormendy & Ho (2013) [56] and McConnell & Ma (2013) [57] to be inconsistent with their data at the 95% and 90% level respectively. They also placed constrains on the strength of environmental coupling effects via parameterisation of the GW background spectrum that allows for a turn-over frequency (see also Sampson, Cornish & McWilliams 2015 [58]). More recently Kelley et al. (2017) [59] used the iIllustris simulation to make a new prediction for the GW background signal in the pulsar timing band and their models are not ruled out by the current upper bounds.
SMBHB mergers are not the only possible source of a GW background in the pulsar timing band. Quantum fluctuations of the gravitational field in the early universe, amplified by an inflationary phase could produce a stochastic relic GW background (e.g., Grishchuk 1977 [60]

Current limitations
Pulsar timing data processing packages are sufficiently advanced to produce pulsar timing residuals at the level needed for GW detection (the tempo2 software package already accounts for all known phenomena that affects the propagation of the pulse from the pulsar to the observatory at the 1 ns level; Hobbs et al. 2006 [64]). For instance, the PTA teams are now achieving sub-100 ns rms timing residuals on a handful of pulsars (this implies that, for at least a few pulsars, calibration errors, radio interference mitigation, interstellar medium propagation effects, orbital motion and many more phenomena can be accounted for). GW detection codes have also been well tested on both real and simulated data sets. The primary issue is that we are simply not achieving the necessary precision for a large sample of pulsars and/or we are not observing enough pulsars. Jenet et al. (2005) [65] demonstrated that at least 20 pulsars are needed to be observed over 5-10 years with an rms timing residual of 100 ns in order to make a low-sigma detection of a GW background. Siemens et al. (2013) [66] considered whether a small number of well timed pulsars or a larger number of poorly-timed pulsars would be preferable and showed that adding in a sufficient number of pulsars with rms residuals of ∼ 1µs can significantly increase the sensitivity of the array.
In many cases the timing precision is not limited simply by the sensitivity of the observing systems. For the very brightest pulsars (such as PSR J0437−4715) we are limited by intrinsic pulse variations (a noise process known as "jitter"; e.g., Shannon et al. 2014 [67]). Jitter currently only affects a relatively small number of the PTA pulsars, but the era of much larger and more sensitive telescopes jitter will become a major limitation for PTA research. On long time scales many of the pulsars are affected by long-term timing variations known as timing noise (see, e.g., Lentati et al. 2016 [68], Caballero et al. 2016 [69], Lam et al. 2017 [70]) and two of the PTA pulsars (PSRs J0613−0200 and J1824−2452A) have been observed to glitch (McKee et al. 2016 [71] and Cognard & Backer 2005 [72], respectively). Even though mitigation routines exist and wide observing bandwidths are used, the interstellar medium is still a limiting factor for PTAs. Levin (2015) [73] reviewed the effects of dispersion and scattering for PTAs.
Ensuring an un-ambiguous detection of a GW signal is also non-trivial. Taylor et al. (2016) [75] presented two methods for assessing the significance of a stochastic GW background signal in a given data set. Zhu et al. (2015) [74] discuss the construction of null streams as a means to provide a consistency check when a single source GW signal is detected. Tiburzi et al. (2016) showed that, without care, false detections of a GW background could be made even if no GW signal were present. These false detections can come from incomplete removal (or modelling) of effects such as errors in terrestrial clocks, the solar wind or in the solar system planetary ephemerides. In particular the planetary ephemerides are currently being actively investigated as the more recent ephemerides seem to induce low-frequency noise into the timing residuals (that can mimic the signature of a gravitational wave background; Shannon et al., private communication). The reason for this is not currently understood, but is being explored by all three PTAs and it is clear that the choice of planetary ephemeris has a significant effect on the resulting data sets and sensitivity to GWs.
7 What next for pulsar timing arrays?
The three main PTA projects have now been ongoing for more than a decade. The data sets produced are the most hightiming-precision data sets yet made. Processing such data has required huge improvements in understanding noise processes and calibration methods. The observations have been used to place stringent constraints on the amplitude of GWs. The bounds have been used to constrain cosmological and cosmic string models and to rule out models of GW emission from supermassive binary black holes. The data have also been used to study atomic clocks and the solar system ephemeris. Unfortunately GWs have not been detected by the PTAs.
Determining when we should expect our first detection is non-trivial. The most recent GW models from supermassive binary black holes predict a background with an amplitude A ∼ 5 × 10 −16 (e.g., Sesana et al. 2016 [76] and Kelley et al. 2017), but the uncertainty is still large. Predictions from cosmic strings or the early Universe are even more uncertain. Groups have attempted to calculate the time to detection based on estimates of the GW amplitude (see, e.g., Siemens et al. 2013 [66] and Taylor et al. 2016 [77]). The estimated times range from a few years from now to a few decades. These calculations are made even more uncertain through the assumptions made about how many new pulsars will be discovered (and added into PTAs) over the coming years, what the timing precision will be of those pulsars, whether we will be able to correct for interstellar medium effects at the necessary level, how intrinsically stable the pulsars are and how will new telescopes (that are currently not part of the PTA experiments) contribute.
Some pulsars (such as PSRs J0437−4715, J1713+0747 and J1909−3744) have significantly smaller timing residuals than many of the other IPTA pulsars. It is therefore probable that the first evidence of GWs will be in the timing residuals of just a few pulsars. The timing residuals for those pulsars will exhibit low-frequency noise or sinusoidal signatures (caused by the GWs) well before correlations can be measured between a large number of pulsar pairs. Such noise will have two effects: 1) GW bounds will stop decreasing as they will be dominated by the noise in these pulsars and 2) even though pulsar astronomers will not be able to state unambiguously that these residuals have been induced by GWs, statements such as "if these residuals are induced by GWs then the amplitude and spectral properties of those GWs are . . . " will become possible.

Current telescopes
Most of the existing telescopes being used for PTA observations will continue to be upgraded. For instance, the Parkes telescope will soon be commissioning a new ultra-wide-bandwidth receiver covering from 700 MHz to ∼4 GHz. This receiver will produce the same observations as the existing receiver suite, but 1) in half the time, 2) with continuous frequency coverage and 3) with cutting-edge radio-frequency-interference mitigation and backend processing. Similar receiver upgrades are being made at many of the other IPTA telescopes and it is clear that future instrumentation for PTA research will be based around wide-band receivers. In the near future it is likely that low-frequency observations from e.g., MWA, LOFAR and the GMRT will play an essential role in monitoring the interstellar medium in the direction of the IPTA pulsars (although note the cautions in e.g., Cordes, Shannon & Stinebring 2016 [78]).

New telescopes
A large number of new telescopes will soon to be significantly contributing to PTA projects (see Figure 3 for a World-map showing current and future PTA telescopes). The most exciting new developments are the Five Hundred Metre Aperture Spherical Telescope (FAST; e.g., Li et al. 2013 [85]) that recently opened in China and MeerKAT in South Africa (e.g., Foley et al.. 2016 [86]). Both telescopes will contribute to 1) finding new millisecond pulsars, 2) high-precision pulsar timing and 3) studies of pulsar noise processes that can only be probed by a detailed study of individual pulses (for instance, jitter noise). On the longer term the Chinese are also planning a large, single-dish, steerable telescope known as the QiTai (or QTT) telescope (see Xu & Wang 2016 [87]). This telescope will operate over a wide range of frequencies and will be able to observe a much larger number of pulsars than is available from FAST. An analysis of how the Qitai and FAST telescopes will contribute to PTAs was provided in Hobbs et al. (2014) [88]. In the Southern hemisphere, the Square Kilometre Array will be built in Southern Africa and Australia. One of the key science goals for the mid-frequency part of this telescope is to detect and study ultra-low-frequency GWs using the methods described in this review (see, e.g., Janssen et al. 2015 [89]).
Over time more and more telescopes will contribute observations to the global IPTA effort. This means that it will become important to optimise which telescopes observe which pulsars. Lee et al. (2012) [90] and Burt, Lommen & Finn (2011) [91] have considered methods for such optimisation. Of course, such optimisation is non-trivial as pulsar observations are used for so much more than GW detection.

Conclusion
Pulsar timing arrays are ongoing with existing telescopes and form key science goals for many of the next generation radio telescopes. PTAs combine observations of pulsars to probe the population of supermassive black holes, the early universe and cosmic strings. These are all exciting objects and topics for the general public. Projects such as PULSE@Parkes (Hobbs et al. 2009) [92] in which high school students observe pulsars using the Parkes telescope, the Arecibo Remote Command Center (ARCC) project in which students use the Arecibo telescope to search for new pulsars and the Pulsar Search Collaboratory (Rosen et al. 2010) [93] with Green Bank ensure that the enthusiasm that the IPTA members have for this project gets passed down to the next generation of pulsar astronomers.
Within decades it is likely that ground-based detectors will be analysing high-frequency GWs in detail, spacebased detectors will study mid-frequency GWs and PTAs will be probing ultra-low frequency GWs. Astronomy over the entire spectrum of GW frequencies will soon become a reality.