Pulsars Probe the Low-Frequency Gravitational Sky: Pulsar Timing Arrays Basics and Recent Results

Abstract

Pulsar Timing Array experiments exploit the clock-like behaviour of an array of millisecond pulsars, with the goal of detecting low-frequency gravitational waves. Pulsar Timing Array experiments have been in operation over the last decade, led by groups in Europe, Australia, and North America. These experiments use the most sensitive radio telescopes in the world, extremely precise pulsar timing models and sophisticated detection algorithms to increase the sensitivity of Pulsar Timing Arrays. No detection of gravitational waves has been made to date with this technique, but Pulsar Timing Array upper limits already contributed to rule out some models of galaxy formation. Moreover, a new generation of radio telescopes, such as the Five hundred metre Aperture Spherical Telescope and, in particular, the Square Kilometre Array, will offer a significant improvement to the Pulsar Timing Array sensitivity. In this article, we review the basic concepts of Pulsar Timing Array experiments, and discuss the latest results from the established Pulsar Timing Array collaborations.

1 INTRODUCTION

Pulsars are highly magnetised, rapidly rotating neutron stars that convert their rotational kinetic energy into magnetic dipole radiation. Although their emission can extend to the entire electromagnetic spectrum, they are typically observed at radio wavelengths.

In a simplified picture, pulsar radio emission is generated in proximity of the magnetic poles, and forms radiation beams. If the magnetic and the spin axes are misaligned, then the two beams rotate with the neutron star and sweep through space. An observer whose line of sight crosses one or both of the beams, will observe a pulsed emission, the period of which corresponds to the pulsar’s spin period (a general introduction to pulsar astronomy can be found, e.g., in Lorimer & Kramer 2005). Pulsars are often referred to as ‘cosmic clocks’, because it is possible to predict the arrival time of each pulse at a telescope, sometimes with sub-microsecond precision (Desvignes et al. 2016; Arzoumanian et al. 2015b; Reardon et al. 2016), through the pulsar timing technique (see Section 3).

Due to their high densities, rapid rotations, strong magnetic fields, and high surface gravity, neutron stars are ideal laboratories for tests of nuclear physics (Lattimer & Prakash 2004), general relativity (GR) (Kramer et al. 2006) and alternative theories of gravity (Shao et al. 2013) in extreme conditions not feasible in Earth-based laboratories (see also Stairs 2003; Chamel & Haensel 2008). In the context of GR tests, Pulsar Timing Array (PTA, Foster & Backer 1990) experiments are among the most exciting projects of the last decade. The primary aim of PTAs is the direct detection of low-frequency gravitational waves (GWs) (Rajagopal & Romani 1995; Wyithe & Loeb 2003; Sesana et al. 2004).

In this article, we describe the basic concepts and approaches of PTA experiments, and we review the recent results from the established PTA experiments. In Section 2, we review the efforts to detect GWs in different parts of the spectrum. In Section 3, we give an overview of the technique of ‘pulsar timing’, used to interpret pulsar data for the purpose of PTAs. In Section 4, we outline the potential sources of GWs at low frequencies, and in Section 5, we describe the basic concepts of PTA experiments. In Section 6, we summarise the latest results from the existing PTA experiments, and in Section 7, we discuss the future prospects of PTAs, also considering the new radio astronomical facilities.

2 THE QUEST FOR GRAVITATIONAL WAVES

2.1. Gravitational waves

GWs were an early prediction of GR (Einstein 1916), and are a consequence of a small (i.e., linearisable) perturbation h _{μ, ν} to an otherwise flat (or Minkowskian) metric η_μν of space–time, produced by asymmetric and accelerated mass distributions:

(1) $$\begin{equation} g_{\mu \nu } = \eta _{\mu \nu } + h_{\mu \nu }. \end{equation}$$

It is possible to demonstrate that the perturbation h _μν propagates in the metric as a transverse wave at the speed of light, and in its propagation, it induces quadrupolar perturbations of space–time. Given a mass distribution (we recall that, in GR, the presence of a mass distribution curves the space–time), it is also possible to demonstrate that such perturbations h _μν are generated if the second time derivative of the quadrupole mass moment Q is not zero (Maggiore 2007).

The amplitude of GWs is typically expressed in terms of the dimensionless strain h, i.e., the fractional change δL induced by GWs over a distance L:

(2) $$\begin{equation} h = \frac{\delta L}{L}. \end{equation}$$

GWs can have two polarisations, commonly referred to as ‘plus’ and ‘cross’. If a ‘plus’-polarised GW propagates along the z axis, then it will alternatively stretch and compress space–time along the y and x axes in the orthogonal direction. A ‘cross’-polarised GW will have the same effect, although rotated by 45°.

Pulsar astronomy brought the first indirect confirmation of the existence of GWs, through observations of PSR B1913+16 (Hulse & Taylor 1974). This object is a pulsar in a ~7.7-h orbit with another neutron star. By assuming the existence of GWs, GR can predict the rate of orbital decay that can be attributed to GW emission due to the orbital motion of the two neutron stars. The orbital decay of this binary system was found in agreement with the predictions, today to a precision greater than 99.5% (Weisberg & Taylor 1981; Weisberg, Nice, & Taylor 2010; Weisberg & Huang 2016, and see Figure 1).

Figure 1. Figure taken from Weisberg & Huang (2016) ¹ ; orbital phase shift of the double neutron star system that includes PSR B1913+16 versus time. The plot shows the perfect agreement between the observed orbital decay of the (black dots) and the prediction by GR (solid line).

2.2. Searching for GWs in the Cosmic Microwave Background

The cosmological inflation is an epoch in the early history of the Universe, that is conjectured to seed structure formation and primordial GWs. The quasi-exponential expansion of the Universe during this phase is thought to have generated a stochastic background of GWs (Starobinskiǐ 1979), that cannot be detected directly with current instrumentation. However, indirect detections may be possible. The inflationary GW background is predicted to have excited both of the polarisation patterns of the Cosmic Microwave Background (CMB); the E-mode pattern (curl-free) and the B-mode pattern (curl). Although the GW-induced E-mode it is not expected to be detectable, the signature in the otherwise quiescent B-mode should be measurable (Polnarev 1985) as an excess power at large angular scales (the recombination bump at l ~ 100, where l is the multipole moment). A detection of the B-mode would provide crucial information in support of the inflationary model. However, this achievement is challenging because of B-mode contaminations given by the gravitational lensing of the E-mode on small angular scales (l ~ 1000), and the polarised foreground emission (such as from dust and synchrotron radiation) from our Galaxy (Tucci, Martínez-González, Vielva & Delabrouille 2005) on spatial scales that are searched for the inflationary signature.

Searches for the B-mode of the CMB polarisation are currently ongoing, through experiments such as POLARBEAR (Kermish et al. 2012), the ongoing observations with the South Pole Telescope (see, e.g., Benson et al. 2014), and the Background Imaging of Cosmic Extragalactic Polarization (see, e.g., Keating et al. 2003). These experiments, that are focusing on smaller l than the satellites, are expected to detect effects of the inflationary GW background at ultra-low frequencies (below 10⁻¹⁶ Hz, Lasky et al. 2016). No detection of the inflationary-induced B-mode has been made to date, while the E-mode was observed for the first time with the Degree Angular Scale Interferometer (Kovac 2002) and followed-up in more detail by the Wilkinson Microwave Anisotropy Probe and the Planck satellites (see, e.g., Komatsu et al. 2011; Planck Collaboration, et al. 2014).

2.3. Searching for GWs with interferometers

In addition to indirect detections based on pulsars or the CMB, it is possible to make direct detections of GWs. This has been achieved by the Advanced Laser Interferometer Gravitational-wave Observatory (aLIGO), and in the future, direct detections will be possible with the Evolved Laser Interferometer Space Antenna (eLISA), and PTAs.

aLIGO comprises two laser interferometers (in Hanford and Livingston) capable of detecting changes in the length of the interferometer arms induced by GWs. As other ground-based laser interferometers such as Virgo in Italy (Acernese et al. 2015), GEO 600 (Willke et al. 2002) in Germany, or KAGRA (Araya et al. 2017, not online yet) in Japan, aLIGO explores GW frequencies from approximatively 1 to 10³ Hz. The lower limit is set by gravity gradients in the Earth gravitational potential, while the upper limit is given by the shot noise of the laser photons (Aasi et al. 2013). GWs emitted in this frequency range are predicted to be generated by coalescing binary systems of neutron stars or stellar-sized black holes. In 2015 September, aLIGO achieved the first direct detection of GWs (Abbott et al. 2016a) from a coalescing binary of stellar-mass black holes. This detection (followed by other six events and one candidate since 2015, Abbott et al. 2016b, 2017) signed the beginning of the era of GW astronomy.

However, more is required to explore this branch of science—more detections, and a wider range of GW frequencies.

eLISA (eLISA Consortium et al. 2013; Amaro-Seoane et al. 2017) is a project of the European Space Agency to deploy a three-body, space-based interferometer with arms 2.5 million km long, that will probe the GW spectrum in a frequency range from 10⁻¹ down to 10⁻⁵ Hz. An eLISA pathfinder, a 40-cm one-armed miniature of the future device, was launched in 2015 and reported significantly lower noise levels than expected (Armano et al. 2017). This success grew confidence and expectations for the mission, planned for launch in 2034. The frequency range limit of eLISA is given by the measurement accuracies of the free-falling test mass accelerations (Amaro-Seoane et al. 2017). The predicted sources of GWs in this frequency range are inspiralling binary systems of white dwarfs and super-massive black holes (SMBH).

The only experiment that can currently provide longer interferometric baselines, on a parsec scale, are PTAs (see Section 5). PTAs explore the frequency range from about 10⁻⁶ to 10⁻⁹ Hz, where the most likely source of GW emission are coalescing SMBH binaries (SMBHB). Other sources might be cosmic strings and relic GWs from inflation. PTAs are experiments based on the monitoring of an ensemble of selected pulsars, in order to search for spatially correlated deviations in the arrival times of their pulses. A number of phenomena can induce such correlations, included GWs.

The GW frequency bands explored by these three kinds of experiments are complementary, as shown in Figure 2.

Figure 2. Figure taken from Janssen et al. 2015 ² ; GW amplitude versus GW frequency, and frequency ranges explored by the interferometric experiments searching for GWs, aLIGO, eLISA, and PTAs. In the ‘PTA band’, the nominal sensitivities for the International Pulsar Timing Array are shown and the Square Kilometer Array, together with a representation of the expected emission from the SMBHB population (solid blue line) in the universe, the emission from ‘GW-loud’ SMBHBs (blue triangles) and from the unresolvable SMBHBs (light purple squares). In the ‘eLISA band’, the nominal eLISA sensitivity curve is shown, together with the expected GW signals from different masses of merging SMBHBs (cyan), a binary with a very high mass ratio (aquamarine), and from the Galactic population of inspiralling white dwarf binaries (yellow). In the ‘aLIGO band’, the sensitivity curve of aLIGO (as of 2015) is shown, together with the expected signals from different inspiralling compact-object binaries (purple). In brown, orange, and red are the GW background expected from inflation and two models of cosmic strings.

3 PULSAR TIMING

As mentioned in Section 1, pulsars are often referred to as ‘cosmic clocks’, as it is possible to predict their times of arrival (ToAs) to high accuracy. Essentially, there are three requirements to enable accurate predictions of the arrival times:

1. 1. The pulsar is a stable rotator. As mentioned in Section 1, pulsars lose rotational energy via magnetic dipole radiation, and therefore they spin down. However, the spin-down might suffer from irregularities, in the form of abrupt (‘glitches’, Downs 1981) or long-term (‘timing noise’, see Section 3) variations in the spin frequency. For high-precision experiments, pulsars should have predictable spin-evolutions.

2. 2. The shape of the integrated pulse profile is stable in time. Pulsars are intrinsically weak sources, with fluxes in the order of a few milliJansky (Lorimer et al. 1995; Kramer et al. 1998). Typically, the individual pulses do not exceed the radiometer noise of the telescope. Thus, many pulsar studies use the integrated pulse profile, i.e., the coherent sum of many thousands of individual pulses. While individual pulses differ often from each other (both in flux distribution and phase), the integrated profile is statistically stable. Known sources of variations in the integrated profile are pulse jitter (Liu et al. 2016a), plasma propagation effects (Geyer et al. 2017), or magnetospheric instabilities (Lyne et al. 2010). The long-term temporal stability of the integrated profile is an assumption for high-precision experiments, and efforts are ongoing to mitigate the impact of integrated profile variations in PTA experiments (Lentati et al. 2017b).

3. 3. The timing model of the pulsar is well known. The timing model of a pulsar (or ephemeris) is a set of parameters that describes the pulsar spin and spin-down, its orbital parameters (if any), its astrometry, and the dispersive influence of the ionised interstellar medium (IISM) along the line of sight to the pulsar. The frequency-dependent dispersive effect of the IISM on radio pulses is quantified by the dispersion measure (DM). The DM is defined as the integrated column density of free electrons along the line of sight:

   (3) $$\begin{equation} \mathrm{DM} = \int _{0}^{{\it d}} n_e \text{d}l, \end{equation}$$
   where d is the distance to the pulsar (pc), and n _e is the free electron number density (cm⁻³).

The first draft of ephemeris for a certain pulsar can be obtained from its discovery, and provides an approximate estimate of the pulsar spin, position, and DM. A precise knowledge of the timing model can be achieved through the technique of pulsar timing (Lorimer & Kramer 2005).

Let us assume that an observing campaign is performed on a pulsar with a given radio telescope. For each observation, we can obtain an integrated pulse profile P, so that the pulse profile is statistically stable and has a suitably high signal-to-noise ratio (S/N). In pulsar timing, average ToAs for each observation are computed via a cross-correlation of each integrated pulse profile with a high-S/N reference template S (Taylor 1992), which is typically a noise-free representation of the pulse profile. This yields a phase shift τ between P and S, if we consider P to be described as

(4) $$\begin{equation} P(t) = a + bS(t - \tau ) + n(t), \end{equation}$$

where a, b, and n(t) represent, respectively, an intensity baseline, an intensity normalisation, and the noise level. The topocentric ToA, ToA_topo is the sum of τ to a time stamp associated with the observation.

The topocentric ToAs are then transformed to the (at first-order) inertial reference frame of the Solar System barycentre (SSB). This conversion is based on the parameters included in the timing models, a reference for the time standard, and for the planetary ephemeris (Edwards et al. 2006):

(5) $$\begin{equation} {\rm ToA}_{{\rm SSB}} = {\rm ToA}_{{\rm topo}} + t_{{\rm clk}} - \frac{D}{f^2} + \Delta _{\rm R} + \Delta _{\rm E} + \Delta _{\rm S}. \end{equation}$$

In this equation, t _clk transforms the reference time standard from the (typically) maser-based clock at the observatory to a world-wide recognised time standard such as Terrestrial Time. The third term removes the effects of observing at non-infinite frequency:

(6) $$\begin{equation} D = \frac{e^2}{2 \pi m_e c} {\rm DM}, \end{equation}$$

where e and m _e are the charge and the mass of an electron, and c is the light speed. Δ_R is the Roemer delay, that corrects for the difference in travel time between the observatory and the SSB. The Roemer delay is purely based on geometrical considerations, and uses the astrometric parameters of the studied pulsar and the planetary ephemeris. Δ_E, the Einstein delay, is based on the planetary ephemeris and corrects for the effects of the gravitational redshift induced by the bodies of the Solar system. Δ_S, the Shapiro delay, accounts for the additional time travel required to the light waves for travelling across the gravitational field of the Solar system. Additional corrective parameters are required if the pulsar is part of a binary system.

Once the barycentric ToAs t has been derived, we compute the pulse number N that represents a ‘counter’ for the number of pulsar rotations:

(7) $$\begin{equation} N (t) = N_0 + \nu _0 (t - t_0) + \frac{1}{2}\dot{\nu }(t - t_0)^2 + ..., \end{equation}$$

where N ₀ is the pulse number at the reference time t ₀, and ν₀ and $\dot{\nu }$ are the spin period at t ₀ and the spin down rate, respectively. The right-hand side of equation (7) is the Taylor expansion of the pulsar spin.

In the last steps of a timing analysis, the parameters included in the timing model can be varied so that the ToAs are spaced of an integer number of pulsar rotations. The refinement of the timing model can be achieved through dedicated software such as tempo2 (Hobbs et al. 2006) and the inspection of the ‘timing residuals’, i.e., the difference between the closest integer number of pulsar rotations and actual number of pulsar rotation among the ToAs. If parameters in the timing model are imprecisely estimated or missing, then we expect to see structures in the timing residuals. For example, we see from equation (7) that an incorrect spin frequency or spin-down rate will show as a linear and a parabolic trend in the timing residuals (see Lorimer & Kramer 2005). If the timing model is sufficiently accurate, then the timing residuals will look ‘white’, i.e., with no correlations (see Figure 3). For reasons that will be explained in Section 5, PTAs are particularly interested in the study of the ‘red noise’. A time series affected by red noise shows long-term correlated structures in the time domain, and an excess in the low-frequency bins of its power spectrum (see Figure 4). Several phenomena can induce red noise, for example, DM variations (You et al. 2007), or intrinsic instabilities in the pulsar spin (better known as ‘spin noise’ or simply ‘timing noise’, Caballero et al. 2016), instrumental imperfections, or GWs.

Figure 3. Timing residuals versus time for PSR J1012+5307. The used observations were obtained at L-band with the Effelsberg radio telescope, and the Effelsberg–Berkeley Pulsar Processor (EBPP) backend. The used ephemeris, with no additional fitting applied, were obtained by Verbiest et al. (2016), based on all the available IPTA datasets, including the EBPP one (i.e., the IPTA data release, see Section 6).

Figure 4. Timing residuals versus time for PSR J1939+2134 (upper panel), and corresponding power spectrum versus frequency (lower panel). The used observations were obtained at L-band with the Effelsberg radio telescope, and the Effelsberg–Berkeley Pulsar Processor (EBPP) backend. The used ephemeris, with no additional fitting applied, were obtained by Verbiest et al. (2016), based on all the available IPTA datasets, including the EBPP one (i.e., the IPTA data release, see Section 6). The timing residuals of PSR J1939+2134 are clearly affected by red noise, most likely spin noise.

4 SOURCES OF GRAVITATIONAL WAVES AT LOW FREQUENCIES

As mentioned in Section 1, GWs are produced by the second time derivative of the quadrupole moment of the mass distribution that distorts the space–time.

Following a dimensional analysis, the amplitude h of a GW is given by Hughes (2003):

(8) $$\begin{equation} h \propto \frac{G}{r c^4} \frac{d^2 Q}{{\rm d}t^2}, \end{equation}$$

where G is the gravitational constant, r is the distance to the GW source, and Q is the quadrupole moment. Due to the factor G/c ⁴, the peak GW amplitude is expected to be small. GWs are therefore more easily detectable when the second derivative of the quadrupole moment is large (Thorne 1987), as in the case of massive, fast-moving objects. The most likely source of GWs at low frequencies are coalescing SMBHBs, although other potential sources have been identified: GWs from inflation (Grishchuk 1974; Starobinskiǐ 1979) or cosmic strings (Kibble 1976).

4.1. Super-massive black hole binaries

Observational evidence shows that SMBHs are hosted at the centre of the most or all galaxies (Kormendy & Richstone 1995; Magorrian et al. 1998). The hierarchical or ‘bottom-up’ scenario (White & Rees 1978) predicts that larger galaxies are generated via merging of smaller galaxies at high redshifts (z). When two galaxies merge, we then expect that the two SMBHs at their centres form a binary system (Begelman et al. 1980; Volonteri, Haardt, & Madau 2003). A binary system is characterised by a non-zero value of the second time derivative of its mass quadrupole moment, thus it is a continuous source of GWs. Defining the masses of the individual SMBHs as m ₁ and m ₂, and assuming a circular orbit for simplicity, with total mass M = m ₁ + m ₂, reduced mass μ = m ₁ m ₂/M, and $\mathcal {M}=(m_1m_2)^{3/5}/M^{1/5}$ (the chirp mass), and using geometricised units such that G = c = 1, the luminosity emitted in GWs (L _gw) by SMBHB is given by (Thorne 1987; Sesana 2013a)

(9) $$\begin{equation} L_{{\rm gw}} = \frac{32}{5} (\pi \mathcal {M} f)^{10/3}, \end{equation}$$

where f is the observed frequency of the emitted GWs, equal to twice the orbital frequency f _B.

The inclination-polarisation averaged amplitude, h, of the radiated GWs is given by (Sesana 2013a)

(10) $$\begin{equation} h = \sqrt{\frac{32}{5}} \frac{G \mathcal {M}}{r c^4} (\pi f)^{2/3}. \end{equation}$$

Equations (9) and (10) describe the simple case of a SMBHB in the local universe (i.e., with zero redshift). As pointed out by Vecchio (2004), it is possible to ‘move’ the GW source at a different redshift by substituting m _x with m _x (1 + z), r with r(1 + z) and f with f(1 + z).

Note that both of the expressions for the GW luminosity and strain contain the GW frequency f. During the binary inspiral, f _B (and hence f) changes in time. This means that the strain h of the propagated GW might not be the same in two different points in space–time. However, during short time scales over which the change in orbital separation is negligible, the SMBHB can be considered a monochromatic source of GWs (Sesana & Vecchio 2010).

We can identify three stages in the evolution of an SMBHB (Flanagan & Hughes 1998):

1. 1. Inspiral, The two SMBHs orbit each other and such orbital separation shrinks due to environmental effects and GW emission.

2. 2. Merger, The two SMBHs coalesce, emitting a GW burst that permanently modifies space-time, called a ‘memory’ event.

3. 3. Ringdown, the SMBH and the nearby space–time undergo to a relaxation that leads to a spherical configuration.

f _gw = 2f _B , the low-frequency GW emission is supposed to happen during the inspiral and merger phases, while the ringdown stage is predicted to generate GWs at higher frequencies.

The ‘memory’ phenomenon (Braginskii & Thorne 1987) is a non-oscillatory GW emission that should occur before the actual coalescence of the two BHs. The final stages of the merger generate a net non-zero contribution to the ‘plus’ polarisation mode of the GW emission. Such a DC offset induces a permanent deformation in the metric (Favata 2009). Memory events can only be detected during their passage through the detectors, that, in the moments when they operate the metric deformation. The strain of a memory event, h _m , is predicted to be (Madison, Cordes, & Chatterjee 2014)

(11) $$\begin{equation} h_m \approx \frac{1-\sqrt{8}/3}{24} \frac{G\mu }{c^2 r} \sin {\mathcal {I}}(17 + \cos ^2{\mathcal {I}})\times [1 + \mathcal {O}(\mu ^2/M^2)], \end{equation}$$

where $\mathcal {I}$ is the inclination angle of the binary before the merger. For r = 1 Gpc, and M ₁ = M ₂ = 10⁹ solar masses, h _m ≈ 10⁻¹⁵ (Madison et al. 2014).

The amplitude of a memory event should rapidly increase in the very final stages of the coalescence before the merging, in a time scale τ ≈ 2πR _s /c, where R _s is the Schwarzschild radius (Cordes & Jenet 2012; Madison et al. 2014).

5 PULSAR TIMING ARRAYS

Here we outline the expected signatures in pulsar timing data given by the various types of GW emission from SMBHBs.

5.1. Modelling a GW signal from single sources

Two processes are predicted to lead to the orbital shrinking of the SMBHB. First, the shrinking is led by environmental effects until the GW emission reaches approximatively nHz frequencies (see Sesana 2013a for a review). Let us assume a circular, evolving (i.e., shrinking) SMBHB, and a pulsar p, characterised by an angle i between the orbital plane of the SMBHB and the line-of-sight towards the pulsar. Because the GW emission from the SMBHB permeates the entire sky in a quadrupolar fashion, we can expect that ToAs from p will be affected, arriving slightly advanced or slightly delayed than the timing model prediction (Sazhin 1978). It can be demonstrated that the effect on the timing residuals R(t) of a single pulsar is independent from the travel path of the radiation (Detweiler 1979). In particular (Babak et al. 2016),

(12) $$\begin{equation} R(t) = R_{\rm E}(t) - R_{\rm p}(t), \end{equation}$$

R _E(t) and R _p(t), called the Earth term and the pulsar term, describe the residuals induced by GWs passing over the Earth and the pulsar, respectively. In the quadrupolar approximation (Lommen 2015), we have that (Babak et al. 2016)

(13) $$\begin{eqnarray} R_{\rm E}(t) &=& \frac{h}{\omega } \lbrace (1+{\rm cos}^2{i})F^+[{\rm sin}(\omega t + \Phi ) - {\rm sin}{\Phi }] \nonumber\\[-3pt] &&+\, 2{\rm cos}{i}F^{\times } [{\rm cos}(\omega t + \Phi ) -{\rm cos}{\Phi }] \rbrace , \nonumber\\[-3pt] R_{\rm p}(t) &=& \frac{h_{\rm p}}{\omega _{\rm p}} \lbrace (1+{\rm cos}^2{i})F^+[{\rm sin}(\omega t + \Phi + \Phi _{\rm p}) - {\rm sin}(\Phi + \Phi _{\rm p})]\nonumber\\[-3pt] &&+ \,2{\rm cos}{i}F^{\times } [{\rm cos}(\omega t + \Phi + \Phi _{\rm p}) -{\rm cos}{\Phi + \Phi _{\rm p}}] \rbrace , \end{eqnarray}$$

where h and h _p are the amplitudes of the GW at the Earth and at the pulsar, ω and ω_p are the GW angular frequency at the Earth and at the pulsar, Φ and Φ_p are the GW phases at the Earth and at the pulsar, and F ⁺ and F ^× are the antenna response functions for the two GW polarisation at the pulsar (i.e., how space–time around the pulsar is affected by the GW). If the SMBHB does not evolve, the power spectrum of the signature described in equation (12) is described by two Dirac delta functions. This signal would be indistinguishable from the signature given by an error in the orbital period and thus not detectable with pulsar timing.

As mentioned in Section 4.1, the last stage of an SMBHB merger lead to a permanent deformation of the metric, similar to a DC offset in space–time. This non-oscillatory phenomenon is propagated, and affects the timing residuals of a pulsar in a way that is equivalent to a change in its rotational spin frequency (vaan Haasteren & Levin 2010). The timing residuals induced by a memory event will be given by (vaan Haasteren & Levin 2010)

(14) $$\begin{equation} R(t) = h_m B(\theta ,\phi )\times [\Theta (t-t_0) - \Theta (t-t_1)], \end{equation}$$

where B(θ, ϕ) = 1/2cos (2ϕ)(1 − cos θ), θ is the angular separation between the pulsar and the SMBHB, ϕ is the angular separation between the principal polarisation of the GW signal and the projected line-of-sight to the pulsar onto the plane perpendicular to the GW propagation direction. Θ is the Heaviside function, while t ₀ and t ₁ are the instants in which the memory event passes the Earth and the pulsar, respectively. In equation (14), it is thus possible to identify an Earth term and a pulsar term, as for the oscillatory contribution to the GW emission shown in equation (12). In particular, the Earth term sensitivity to a memory event is found to increase with the square root of the number of pulsars included in a PTA (Cordes & Jenet 2012).

5.2. Modelling a GW background signal

The expected number of SMBHB systems is extremely large, up to 10⁶ depending on the redshift and the mass range of the involved BHs (Sesana, Vecchio, & Colacino 2008). The choral GW signal coming from such a population of SMBHBs gives rise to an incoherent superposition of the individual GW signals, that effectively generates a stochastic background of GWs (GWB, Sesana et al. 2008; Ravi et al. 2015), usually considered isotropic.

The GWB is predicted to induce a red-noise signal in pulsar timing residuals, with a power spectrum P(f) that can be described by a steep power-law (Phinney 2001):

(15) $$\begin{equation} P(f) = \frac{h^2}{12\pi ^2} \left( \frac{f}{f_{{\rm 1yr}}} \right) ^{2\alpha -3}, \end{equation}$$

where f _1yr normalises the GW frequency at 1/1_yr, h is now the amplitude of the GWB, and α is a coefficient whose value is 2/3 in the case of an isotropic and stochastic GWB, thus the spectral index for a GWB is expected to be −13/3. In the case of a GWB, Ω_gw(f), the ratio between the energy density ρ_gw of the GWs (per unit logarithmic frequency) and the critical energy density of the Universe ρ _c , is related to the strain h as (Allen & Romano 1999):

(16) $$\begin{equation} \Omega _{{\rm gw}} (f) = \frac{2 \pi ^2}{3 H_0^2} f^{2} h^2(f), \end{equation}$$

where H ₀ is is the Hubble expansion rate (100 h _H km s⁻¹ Mpc⁻¹, with h _H being the dimensionless Hubble parameter). In the case of relic GW from inflation and cosmic strings, the spectral index of equation (15) is expected to be −5 (Grishchuk 2005) and −16/3 (Damour & Vilenkin 2005), respectively.

5.3. The Hellings and Downs curve

Although SMBHBs are considered the loudest sources of GWs in the universe, the amplitude of such emission is predicted to be extremely tiny (Arzoumanian et al. 2014; Babak et al. 2016; Zhu et al. 2016). The GWB amplitude is expected to exceed the amplitude of the signals from the vast majority of individual SMBHBs. This implies that the detection of a GWB is much more likely than GWs from an individual SMBHB (Rosado, Sesana, & Gair 2015).

As already discussed, the GWB is a stochastic signal that can be described as a red noise process [equation (15)]. Given an individual pulsar, a GWB signal cannot be distinguished unequivocally from other red noise processes such as timing noise, IISM effects, clock noise, or ephemeris errors. Additionally, the pulsar timing procedure absorbs all of the power present in the first two bins of the power spectrum of the timing residuals (corresponding to the spin period and spin period derivative) and in the 1/1yr frequency bin (due to the orbital period of the Earth), effectively decreasing the detectability of a GW signal at this frequency.

However, a GW emission (both from SMBHBs and a background) would affect different pulsars in a correlated fashion depending on their respective sky position. Therefore, detection statistics are based on the correlation among the timing residuals of an array of pulsars (Romani 1989; Foster & Backer 1990).

The correlation C among the timing residuals of pairs of pulsars perturbed by an isotropic and stochastic GWB was studied by Hellings & Downs (1983). Given a pair of pulsars (i, j), separated by an angle θ _{i, j} in the sky, they demonstrated that C takes a specific functional form, known as ‘Hellings and Downs curve’:

(17) $$\begin{equation} C(\theta _{ij}) = \left[ \frac{3}{2} x\, {\rm log}(x) - \frac{x}{4} + \frac{1}{2} \right](1+\delta _{{\rm {\it i,j}}}), \end{equation}$$

where x = (1 − cosθ_{i, j})/2. The Hellings and Downs curve is the sky- and polarisation-averaged angular correlation between pairs of pulsars. In the computation, the GWB is assumed to be isotropic (i.e., the power spectrum of the GWB does not have an angular dependence), and the short-wavelength approximation to be valid (i.e., f _gw r > >1, that is, the distance between the Earth and the pulsar, and between the pulsars in the array, is large if compared to the wavelength of the GWs). It should be noted that the Hellings and Downs curve is computed using the Earth term only. The pulsar term is estimated to bring a significant contribution only at angular distances close to zero, and only if the pulsar pair is effectively close in space. In this case, the contribution of the pulsar term brings the angular correlation to 1. For pulsar pairs at even smaller angular distances, the contribution of the pulsar term becomes rapidly negligible (Mingarelli & Sidery 2014). Figure 5 shows the Hellings and Downs curve, without taking into account the additional correlation that would occur at θ_{i, j} = 0 when considering the pulsar term.

Figure 5. Angular correlation C given by the Hellings & Downs curve as described by equation (17) (minus the contribution of the pulsar term) versus angular distance θ.

5.4. Aims and characteristics of Pulsar Timing Arrays

PTA experiments aim to detect signals that are angularly correlated across the sky, using the clock-like behaviour of an array of hyper-stable pulsars. The primary goal of PTAs is the detection of low-frequency GWs, and the most likely GW source to emit in the PTA band is coalescing SMBHBs. In this sense, PTAs can be considered as interferometer on Galactic scales, although instead of lasers, PTAs exploit the pulsed radio emission from the pulsars in the array.

To aid our detection prospects, we select pulsars with high rotational stability for PTA analysis (Shannon & Cordes 2010). The most rotationally stable pulsars are millisecond pulsars (MSPs; Alpar et al. 1982). MSPs are pulsars that have been spun up via a transfer of mass and angular momentum by a companion star, which accelerates the neutron star to spin periods in the order of milliseconds (Bhattacharya & van den Heuvel 1991). Following the mass transfer, both the magnetic field intensity and the spin-down rate are remarkably lowered. MSPs are much more stable (Verbiest et al. 2009) than normal pulsars, and are characterised by timing residuals that are typically lower and whiter (Hobbs et al. 2004). As such, they are the only class of pulsars that are included in PTA monitoring campaigns.

The sensitivity of PTA experiments lies in the frequency range from approximatively 10⁻⁶ to 10⁻⁹ Hz. The two boundaries are due to the limits imposed by the Nyquist theorem, and are set by the observing cadence at the higher frequency (assumed to be once per month) and the total timespan at the lower frequency (assumed to be around 20 yrs.

No GW detection has yet been made by PTA experiments. However, the upper limits on the GW amplitudes estimated by PTAs have already given powerful insights in the models for Galaxy formation, aiding to exclude a fraction of them (Sesana 2013b; Shannon et al. 2015).

With the current sensitivities, and amplitude and rate predictions, it is unlikely that PTAs will detect individual SMBHBs in the near future, either in the form of continuous wave or in the form of a memory event (Babak et al. 2016; Ravi et al. 2015; Wang et al. 2015). Concerning a GWB, as mentioned, although its amplitude is thought to be higher than that of individual SMBHB, the peak strain is expected to be very low. Sesana et al. (2016), an update of Sesana (2013b), estimated the spectrum of the GWB amplitude generated by a population of GW-driven, adiabatically inspiralling SMBHBs in quasi-circular orbits, and demonstrated that the majority of the GWB is due to major mergers, where the mass ratio between the two SMBHs is >0.25 within z = 1.5, and for black hole masses larger than 10⁸ M_⊙ (Sesana et al. 2008). The study also assumes values from different studies available in literature to account for the SMBHB merger rates and masses. Combining the values from observational constraints, the authors generated more than 2 500 realisations of a GWB and computed a distribution for its amplitude. At 3σ, they predict 1.4 × 10⁻¹⁶ < A < 1.1 × 10⁻¹⁵ at 95% confidence. Such uncertainty mainly stems from poorly constrained estimates for the galaxy merger rate and the relation between the mass of the SMBH and the mass of the host galaxy. The influence of the SMBH-host relations is shown in Sesana et al. (2016), where the authors compare the GWB predictions obtained by using two of these relations (Kormendy & Ho 2013 and Shankar et al. 2016). Kormendy & Ho (2013) is claimed to be biased high, due to a number of overestimated SMBH masses obtained through dynamic measurements, while Shankar et al. (2016) claims to have corrected the bias. The results, shown in Figure 6, indicate that the GWB predictions based on the two BH-host relations differ by a factor 3. This highlights the importance of the relations in these studies, and the necessity of refining them.

Figure 6. Figure taken from Sesana et al. (2016) ³ ; GWB amplitude versus f _gw. The plot shows a comparison between the GWB predictions as based on the BH-host relations from Kormendy & Ho (2013) (left panel) and Shankar et al. (2016) (right panel). The sensitivity curves for EPTA, PPTA, and NANOGrav are shown in green, blue, and red, respectively, and the differently shaded area represents 99.7, 95, and 65% of probability.

GWs, both in the form of emission from single sources and a background, are not the only signals that can be angularly correlated among pulsars in a PTA. For example, imperfections in the reference time standards and in the planetary ephemeris used to identify the SSB would induce, respectively, a monopolar and a dipolar angular correlation in the timing residuals (Foster & Backer 1990), and the creation of a pulsar-based time reference and the improvement of planetary ephemerides are ongoing projects within the framework of PTAs (Hobbs et al. 2012; Champion et al. 2010). Tiburzi et al. (2016) studied the impact of correlated noise processes other than GWs (such as errors in time standards, planetary ephemeris, and unmodelled effects of instrumentation and the Solar wind) on PTA sensitivity to the GWB. The study demonstrated that, without including mitigation techniques in the detection pipelines, such signals can induce false detections (see also Taylor et al. 2017). Mitigation is feasible, especially for the monopolar signal, while the dipolar signal is more difficult to subtract without compromising the GWB search.

5.5. Current PTA collaborations

There are currently three well-established collaborations in the world that are leading PTA experiments: the European Pulsar Timing Array (EPTA, Desvignes et al. 2016) in Europe, the Parkes Pulsar Timing Array (PPTA, Reardon et al. 2016; Manchester et al. 2013) in Australia, and the North American NanoHertz Observatory for Gravitational waves (NANOGrav, Arzoumanian et al. 2015b) in the North America. EPTA, PPTA, and NANOGrav, all based on MSP observations with 100-m class radio telescopes, collaborate as the International Pulsar Timing Array (IPTA, Verbiest et al. 2016).

5.5.1. EPTA

The EPTA was officially established in 2005, and currently monitors 42 MSPs (Desvignes et al. 2016) at an approximatively monthly cadence with each of the five largest radio telescopes in Europe: the Effelsberg Radio Telescope (Eff, Germany), the Nançay Radio Telescope (NRT, France), the Westerbork Synthesis Radio Telescope (WSRT, the Netherlands), the Lovell Telescope at Jodrell Bank Observatory (JBO, UK), and the Sardinia Radio Telescope (SRT, Italy). In addition, a special programme within the EPTA, the Large European Array for Pulsars (LEAP), effectively acts as a sixth EPTA telescope (Bassa et al. 2016a).

Besides the dataset collected since its establishment, the EPTA uses archival data of MSPs dating back to the 1990s, and collected under different timing proposals with ‘historical’ backends and receivers (i.e., pulsar instruments now decommissioned).

The current observing setup of the EPTA telescopes is as follows:

• • EFF. Performs coherently dedispersed observations of MSPs at three different frequencies with the PSRIX backend (Lazarus, Karuppusamy, Graikou, Caballero, Champion, Lee, Verbie 2016): 1360 MHz, 2639 MHz, and 4800 MHz.

• • JBO. Observes MSPs with two backends in parallel, the DFB (incoherent dedispersion, see Manchester et al. 2013) and the ROACH (coherent dedispersion, see Bassa et al. 2016a), at 1532 MHz.

• • NRT. Performs coherently dedispersed observations of MSPs in two frequency ranges, between 1100 and 1800 MHz, and between 1700 and 3500 MHz with the NUPPI backend (Liu et al. 2014).

• • WSRT. WSRT is currently unavailable for EPTA observations, as a new backend and frontend (ARTS and APERTIF) are being commissioned. The previous setup performed coherently dedispersed observations at 345, 1380, and 2273 MHz with the PuMa II backend (Karuppusamy, Stappers, & van Straten 2008). The receiver at 345 MHz has been officially decommissioned.

• • SRT. The first official EPTA observing run commenced in 2016, performing observations between 305 and 410 MHz and between 1300 and 1800 MHz (sometimes simultaneously) with a DFB and a ROACH backend.

• • LEAP: Performs coherently dedispersed, interferometric observations of MSPs with the five EPTA telescopes at 1396 MHz, using ROACH backends, and collects dual-polarisation baseband data, that are then correlated offline.

5.5.2. PPTA

The PPTA project commenced in 2005, and currently monitors 24 MSPs (Reardon et al. 2016) with the Parkes Radio Telescope (NSW, Australia) every two to three weeks. In addition to observations obtained within the PTA programme, the PPTA uses data sets collected since the 1990s for other timing campaigns. Currently, the observations are carried out at three different frequencies: 3000, 1500, and 600 MHz, using a DFB (incoherent dedispersion) and the CASPSR ⁴ (coherent) backends.

5.5.3. NANOGrav

NANOGrav was officially established in 2007, and currently monitors 59 MSPs (a selection that has been expanded after Arzoumanian et al. 2016) with the Arecibo Observatory (AO, Puerto Rico), the Green Bank Telescope (GBT, West Virginia, USA) and the Very Large Array (VLA, New Mexico, USA), every three or four weeks (Arecibo and Green Bank are also carrying out weekly observations of a subset of the monitored MSPs).

The current observing setup of the NANOGrav telescopes is as follows:

• • AO. Performs coherently dedispersed observations of MSPs at 430, 1410, and 2030 MHz with the PUPPI backend (DuPlain et al. 2008).

• • GBT. Performs coherently dedispersed observations of MSPs at 820 and 1500 MHz with the GUPPI backend (DuPlain et al. 2008).

• • VLA. The newest addition to the NANOGrav programme, it observes MSPs between 1000 and 2000 MHz and between 2000 and 4000 MHz since 2017.

5.5.4. New PTA collaborations

Efforts to establish PTA experiment are ongoing in India, China, and South Africa.

The Indian PTA observes MSPs with the Ooty Radio Telescope (ORT) and the Giant Metrewave Radio Telescope (GMRT, both conventional and upgraded). In particular, conventional GMRT (that is timing nine MSPs together with ORT) has 32 MHz of bandwidth available in coherent dedispersion, while updated GMRT (that is timing 18 MSPs) has 200 MHz of bandwidth available in incoherent dedispersion (coherently dedispersion will be available in the near future). ORT has a central frequency of 334 MHz and can observe with coherent dedispersion within 16 MHz of bandwidth (M. Bagchi, private communication).

The Chinese PTA had an inaugural meeting in 2017 May. The Chinese PTA operates several 100-m class telescopes (e.g., NSRT, Kunming, Tianma), but the two most important facilities will be the Five hundred meter Aperture Spherical Telescope (FAST, Peng et al. 2001) and the QiTai Radio Telescope (QTT). Once combined, FAST and QTT will be sensitive to a GWB amplitude of 2 × 10⁻¹⁶ within a few years of observations (Lee 2016).

MeerTIME is an approved proposal dedicated to pulsar timing, that will use the MeerKAT telescope (South Africa, Booth et al. 2009). MeerKAT is one of the numerous pathfinders for the Square Kilometer Array (SKA) (see Section 7), and is currently under deployment. Among the planned pulsars that will be observed with an increasingly larger number of antennas, are several PTA-relevant sources.

6 RECENT RESULTS FROM PTAS

6.1. EPTA

The current EPTA data set is comprised of 42 MSPs, and is presented in Desvignes et al. (2016). It includes updated timing solutions and ToAs spanning more than 15 yrs of data for many of the presented MSPs, besides deepening the astrometric properties of the sources. Caballero et al. (2016) studied the red noise properties of the EPTA dataset, and found a significant level of red noise in 25 MSPs. Errors in the time standards were estimated to affect for at most 1% of the total noise budget, reducing the sensitivity to the GWB and resolvable SMBHBs.

The six most stable EPTA MSPs were used to derive the upper limits on the GWB amplitude and GWs from individual SMBHB, and to search for anisotropies in the GWB. Lentati et al. (2015) computed a robust upper limit on the GWB amplitude of A < 3.0 × 10⁻¹⁵, taking into account the presence of other spatially correlated noise (Tiburzi et al. 2016). Babak et al. (2016) shows that the highest sensitivity to resolvable sources is reached by EPTA between 5 and 7 × 10⁻⁹ Hz, with a strain amplitude limit at 95% between 6 and 14 × 10⁻¹⁵.

Taylor et al. (2015) assessed that the current EPTA dataset cannot constrain the angular distribution of the anisotropies yet, but their amplitude is 40% of the effect given by the isotropic GWB.

The dataset presented in Desvignes et al. (2016) and complemented with historical data was used to carry out individual-pulsar studies of MSPs J1024−0719 (Bassa et al. 2016b), J0613−0200 (McKee et al. 2016), and J2051−0827 (Shaifullah et al. 2016). The EPTA project LEAP (Bassa et al. 2016a; Smits et al. 2017) presented a single pulse analysis of MSP J1713+0747 (Liu et al. 2016b), important to assess the impact of pulse jitter on timing precision.

6.2. PPTA

Reardon et al. (2016) presented an extension to the first PPTA data release (Manchester et al. 2013), which included new timing solutions for 20 MSPs and their red noise analysis based on a new version of the Cholesky method (Coles et al. 2011). This study includes the first distance to a pulsar, MSP J0437−4715, measured to sub-parsec precision. A multi-frequency polarisation and spectral analysis of the PPTA MSPs was presented in Dai et al. (2015), finding deviations from the models commonly applied to study pulsar spectra and Faraday rotation in some of the pulsars.

Additional studies of pulse jitter (Shannon et al. 2014), extreme scattering events (Coles et al. 2015), differences in measured positions between VLBI, and pulsar timing studies (Wang et al. 2017), and variations in the pulse profiles (Shannon et al. 2016) were presented between 2014 and 2016.

The now-established technique of profile-domain pulsar timing (Lentati et al. 2014) has been expanded to include the frequency evolution of pulse profiles (Lentati et al. 2017a) and the impact on the pulse profile due to the variable scattering effects of the IISM (Lentati et al. 2017c).

Shannon et al. (2015) used the four most stable PPTA sources and placed the most constraining upper limit on the GWB amplitude to date, 1 × 10⁻¹⁵. Madison et al. (2016) developed a new technique to search for individual GW sources, without constrains on the waveform, and the PPTA dataset was searched for individual SMBHBs and memory bursts by, respectively, Zhu et al. (2014) and Wang et al. (2015). No evidence for GWs was found, and the two studies placed upper limits on the amplitude of the two events. In the case of resolvable SMBHB, an upper limit of A < 1.7 × 10⁻¹⁴ was found at 10⁻⁸ Hz, while no burst events with an amplitude lower than 2 × 10⁻¹⁴ could have been detected in the PPTA dataset studied in Wang et al. (2015).

6.3. NANOGrav

The NANOGrav collaboration published its latest data release in 2015 (Arzoumanian et al. 2015b), which included the timing solutions for 37 MSPs obtained from datasets spanning up to 9 yrs. New methods were developed to account for variable DM and profile evolution with frequency, and 10 pulsars were found to be affected by red noise. Fonseca et al. (2016) measured the Shapiro delay and masses for 14 MSPs in binary systems in the NANOGrav dataset, while Matthews et al. (2016) studied the astrometry of the 37 sources, finding the velocity dispersions to be much smaller than for the general pulsar population.

Detailed analyses of the effects induced by the turbulent IISM were also carried out. Levin et al. (2016) analysed the scattering contribution, and concluded that the effect on the ToA errors due to variable multi-path propagation effects is negligible. Jones et al. (2017) studied the DM variations, finding incompatibility with a Kolmogorov spectrum (Armstrong, Rickett, & Spangler 1995) in four of the pulsars, but the discrepancies can be explained by the presence of unaccounted trends in the data.

Noise analyses were conducted both on short (Lam et al. 2016) and long timescales (Lam et al. 2017), and five more pulsars in addition to those identified by Arzoumanian et al. (2015b) were found to be affected by red noise.

The 9-yr NANOGrav dataset was searched for GWs (Arzoumanian et al. 2016). No evidence of a GWB was found, and an upper limit on the GWB amplitude was set at 1.5 × 10⁻¹⁵. The previous NANOGrav dataset (Demorest et al. 2013) was also searched for GWs from individual sources, in the form of continuous GW emission (Arzoumanian et al. 2014) and memory bursts (Arzoumanian et al. 2015a). No evidence for GWs was found, but upper limits were placed on the amplitude of continuous waves (A < 3.0 × 10⁻¹⁴ at 10⁻⁸ Hz) and for the occurrence rate of memory bursts depending on their amplitude (e.g., memory bursts with an amplitude larger than 4 × 10⁻¹⁴ at 6.2 yr⁻¹).

6.4. IPTA

Verbiest et al. (2016) and Lentati et al. (2016) presented the first IPTA data release, based on the combination of the EPTA, PPTA, and NANOGrav datasets for 49 MSPs (see Figure 7). The IPTA dataset consists of the ToAs time series, timing solutions, and noise models for the 49 sources. The noise analysis carried out by Lentati et al. (2016) showed that the two main sources of red noise are variable DM, and intrinsic timing noise. However, these two sources of noise are often indistinguishable, due to a lack of multifrequency data.

Figure 7. Figure taken from Verbiest et al. (2016) ⁵ ; Aitoff projection of the IPTA MSPs. The two axes represent Galactic longitude (l) and latitude (b), while the stars represent the position of the IPTA MSPs. The crosses represent the MSPs that have been detected at radio wavelengths which are not part of a globular cluster present in the ATNF Pulsar Catalogue at the time of writing.

A basic search for a GWB was carried out in Verbiest et al. (2016), using all of the pulsars in the array. No evidence for GWs was found, and the IPTA placed an upper limit on the GWB amplitude of 1.7 × 10⁻¹⁵. This value, higher than the most stringent upper limit from PTA experiments (1 × 10⁻¹⁵, Shannon et al. 2015), is more constraining than that obtained by the individual PTAs. This indicates that the IPTA as a whole is more sensitive than the individual PTAs by at least a factor of two (Verbiest et al. 2016).

7 FUTURE PROSPECTS

The main challenge of PTAs progresses lies in identifying and correcting for corrupting effects on the ToAs, many of which were of no importance until recently. As far as current efforts go, the bulk of the ongoing research is dedicated to study several long-period processes affecting the residuals—effects such as inaccuracies in the Solar System ephemerides, the IISM, intrinsic pulsar-timing noise, and instrumental instabilities. In addition to those, there are continuous efforts to increase the number of highly precise MSPs in the arrays, decreasing the levels of white noise, increasing observing baselines, cadence, and frequency coverage, improving analysis methods for multifrequency data, tackling previously intractable issues visible in unprecedented high-S/N data.

The new generation of radio telescopes that are now coming online will greatly increase our sensitivity to low-frequency GWs. The most sensitive instrument will be the SKA (Braun et al. 2015). The SKA will be built in Western Australia (low-band antennae) and in South Africa (mid- and high-band antennae), and will boost the sensitivity beyond the limits currently set by radio telescopes. The predicted probability of a GW detection after 5 yrs of observations with a SKA-based PTA (even without taking into account the current IPTA dataset, and assuming the original SKA design) is 50% (Janssen et al. 2015).

In preparation for the SKA, several pathfinders have been deployed, such as MeerKAT and the MWA (Western Australia, Tingay et al. 2013), LWA (New Mexico, USA, Ellingson, Clarke, Cohen, Craig, Kassim, Pihlstrom, Rickard 2009), and the LOw Frequency ARray (Europe, van Haarlem 2013), and they are proving to be vital to tackle some of the mentioned main challenges. For example, the low-frequency facilities such as LOFAR, LWA, and MWA, are fundamental instruments to monitor the turbulent IISM and its effects on pulsar timing, due to the frequency dependence of IISM effects on the propagation of radio waves. In particular, DM variations are one of the main sources of red noise in the ToA time series. IISM studies at low frequencies will be able to provide invaluable insights to improve the red noise models, and to disentangle the IISM contribution from intrinsic timing noise generated from instabilities in the pulsar spin.