Principles of Gravitational-Wave Detection with Pulsar Timing Arrays

Pulsar timing uses the highly stable pulsar spin period to investigate many astrophysical topics. In particular, pulsar timing arrays make use of a set of extremely well-timed pulsars and their time correlations as a challenging detector of gravitational waves. It turns out that pulsar timing arrays are particularly sensitive to ultra-low-frequency gravitational waves, which makes them complementary to other gravitational-wave detectors. Here, we summarize the basics, focusing especially on supermassive black-hole binaries and cosmic strings, which have the potential to form a stochastic gravitational-wave background in the pulsar timing array detection band, and the scientific goals on this challenging topic. We also briefly outline the recent interesting results of the main pulsar timing array collaborations, which have found strong evidence of a common-spectrum process compatible with a stochastic gravitational-wave background and mention some new perspectives that are particularly interesting in view of the forthcoming radio observatories such as the Five hundred-meter Aperture Spherical Telescope, the MeerKAT telescope, and the Square Kilometer Array.


Introduction
Gravitational waves (GWs), predicted by Albert Einstein's theory of General Relativity (GR) [1], are ripples in space-time propagating at the speed of light in a vacuum, emitted by compact object systems characterized by a quadrupole moment with a non-null second time derivative. Since the strain associated with GWs is proportional to their amplitude h, which is in most cases less than 10 −16 , even Einstein himself wondered if they could ever be discovered [2]. Over the years, several different techniques have been developed, up to the first successful event detection, named GW150914 [3], emitted during a black-hole (BH) merger, which was achieved in 2015 by the LIGO (Laser Interferometer Gravitational-Wave Observatory) collaboration, thanks to the Livingston and Hanford gravitational interferometers. Then, a new window for observational astronomy opened up, indicating, especially with the discovery of the first binary neutron star merging event (GW170817) [4] observed alongside its electromagnetic counterpart, i.e., a Short Gamma-Ray Burst (GRB170817) [5] detected by the Fermi satellite, the start of the multi-messenger era.
Ground-based laser interferometers, such as those of the LIGO and VIRGO collaborations, are sensitive only to high-frequency (i.e., in the frequency range [10, 10 4 ] Hz) GWs, so the

Pulse Time of Arrivals
The pulsar ToAs are not directly measurable but must be inferred by the observed pulse profile, which shows how the pulsar luminosity changes with the pulsar phase. However, single pulses are usually too weak to be detected one by one, and there is a large variability between each observed pulse profile. A common procedure used to increase the signal-tonoise ratio [22] and to improve the detection stability is to produce the integrated pulse profile, taking the average of a series of pulses. This operation reduces the uncorrelated Gaussian noise, giving a more clean and stable pulse profile. The corresponding ToA is then obtained by cross-correlating it with a high signal-to-noise template, which is typically the sum of pulse profiles over many epochs [23].
The observed ToAs carry a lot of useful information, which can be extracted by comparing them to the predicted ToAs. To obtain a timing model, one must express the N-th rotation of the pulsar as Taylor series, with respect to time [24]: Equation (1) can be inverted in order to obtain the expression of the N-th ToA: Note that Equation (2) only holds at the position of the pulsar. However, the ToAs are measured with respect to the Earth, which is not a good reference frame because of its motion around its axis and around the Sun. Therefore, it is more convenient to use a reference frame centered in the Solar System barycenter (SSB) that, within this context, can be considered as inertial. The reference frame currently used for pulsar timing is the International Celestial Reference Frame (ICRF). This is a realization of the International Celestial Reference System (ICRS), which is the standard reference system adopted by the International Astronomical Union (IAU), centered in the SSB. Once the position of the SSB [25] is determined with sufficient precision, the timing model, in the general case of a binary pulsar, must include all the effects that can cause delays or advances in the ToAs [26]: Here, the term ∆ includes the effects pertinent to the transformation of the pulse ToAs at the Earth (t ⊕ ) to the equivalent at the SSB, namely atmospheric delays, vacuum retardation due to observatory motion (Roemer delay and parallax), the dispersion delay due to the free electrons content of the solar wind, the effects of relativistic frame transformation (Einstein delay) and the delay experienced by photons as they go through the gravitational potential generated by Solar System bodies (Shapiro delay). The ∆ ism term includes the effects due to the propagation time of the signal from the pulsar (or binary system barycenter) to the SSB, namely the vacuum propagation delay (the path length divided by the light speed in vacuum), the dispersion delay due to the free electron content of the interstellar medium and other frequency-dependent delays, and ∆ bin includes the effects caused by the presence of a binary companion, namely the excess in the vacuum light travel time due to the Euclidean displacement of the pulsar (binary Roemer delay), a pseudo-delay that accounts for the aberration of the radio beam by binary motion, the gravitational redshift and special relativistic time dilation in the pulsar frame (binary Einstein delay) and the gravitational time dilation in the proximity of the companion (binary Shapiro delay). Obviously, in the case of an isolated pulsar, all the terms associated with the binary companion can be set to zero. For an exhaustive compendium of all these effects and their mathematical expressions, the reader is referred to Ref. [27].
The majority of the effects affecting ToAs are, in general, well modeled and, therefore, they can be removed with the available pulsar timing analysis software such as Tempo2 [27,28], which is one among those most commonly used. It follows that the most significant effects for pulsar timing analysis are those that cannot be modeled due to their stochastic nature. This is the case, for example, of the dispersion delay ∆ dis , which is given by: where e is the electron charge, m e is the electron mass, ν rad is the MSP radio-emission frequency, d p is the MSP distance from the observer, and n e is the electron number density [22,29]. For convenience, geometrical units c = G = 1 have been adopted unless otherwise specified. As can be seen from Equation (4), the dispersion delay is subject to stochastic variations caused by the random interstellar medium density fluctuations.

Timing Residuals
The values of the timing model's parameters can be obtained by fitting the function that predicts the ToAs to the observed ones, including the corrections expressed by Equation (3). The timing model has to be continuously updated to include newly discovered effects, as long as there are clock corrections, if any. If all known effects that can influence the ToAs have been taken into account, the timing residuals, calculated by the difference between the fitting function and the ToAs data-set, are expected to be uniformly distributed around zero. The presence of structures in the timing residuals can be, therefore, the signature of unmodeled phenomena. In particular, as will be discussed in Section 3, GWs can induce timing residuals in the pulsar ToAs, with a shape dependent on the GW source.

Gravitational-Wave Detection
The detection of GWs is essentially based on a GR result according to which their passage can "stretch" and "compress" the space-time. This implies that the time of flight of an electromagnetic signal propagating from one position to another is not a constant, but can be longer or shorter than it would have been in the absence of GWs. This principle, which has already been successfully applied to the detection of high-frequency GWs through gravitational interferometry, would also allow revealing ultra-low-frequency GWs by pulsar timing. The frequency range to which PTAs are sensitive is determined by the duration of the observation campaign and the cadence of observations. The current observation campaign is thirteen years long, and observations are performed mainly once per week. This leads to sensitivity in the frequency range [10 −9 , 10 −6 ] Hz.
To explain how a GW can influence the ToA of a pulse, let us first consider an ideal situation where an observer is located in the SSB, which is the origin of a right-handed threedimensional reference frame, and a MSP located at a distance d p from the observer. As next, let us consider a GW emitted by a source at distance d s (with d s d p ), which goes through this two-body system, ignoring for the moment the presence of other effects which will be discussed in Section 2. This allows us to use the plane-wave approximation, thereby writing the GW as a superposition of its polarization states. Since the GW amplitude is small, the space-time metric tensor g µν can be written as the sum of the Minkowski tensor η µν and the perturbation tensor h µν , with |h µν | |η µν |. Writing down the squared distance ds 2 between two space-time events and considering that for light-type events ds 2 = 0 one obtains a differential equation that, if integrated along the unperturbed light path for two consecutive MSP pulses, gives the period variation that can be measured by the observer. If the pulsar period is much less than the time of flight of the GW period, a condition that especially holds for ultra-low-frequency GWs, one has: where p i is the versor in the direction of the pulsar, t e is the time of the emission of the first pulse, t p is the light travel time between the pulsar and the observer and the perturbation tensor depends on the time t and the position x i . From Equation (5), writing the tensor h ij as the superposition of the "plus" (+) and "cross" (×) polarization states of the GW, and making explicit the dependencies of the amplitude h A of the A-th polarization state as is the wavenumber four-vector of the GW, one finds [30]: where z = −∆P/P is the observed relative period variation, the superscript A indicates the A-th polarization, and F A is the antenna pattern function, given by: where e A ij is the polarization tensor, and Ω i is the versor in the direction of the GW emission. In Equation (6), the two terms in the square brackets are often called, respectively, the Earth term and pulsar term because the former indicates the metric perturbation in the proximity of the Earth (at SSB), while the latter refers to the pulsar.
The quantity z will be here referred to as pulse redshift to distinguish it from the photon redshift [31], which is the observed frequency variation of an electromagnetic signal. Actually, it can be shown that a GW can also induce a photon redshift described by the same equation, but in this case, the detection technique is not based on pulsar timing [32]. The GW-induced timing residuals can be defined by integrating Equation (6) with respect to the time: Equation (8) describes the observed variation of the ToAs of the pulses induced by the GW [33]. The antenna pattern function, defined in Equation (7), describes how the gravitational radiation is irradiated into space. This appears more clearly when the antenna pattern function is written explicitly. In the adopted reference frame, the z-axis is oriented in the direction of the GW source and the coordinates θ and φ indicate, respectively, the angle between the direction of the GW source and that of the MSP and the angle between the x-axis and the projection on the xy-plane of the direction of MSP. With these choices, one obtains Ω i = (0, 0, 1) and p i = (sin θ cos φ, sin θ sin φ, cos θ). Considering that in GR the polarization tensors are: and using the Weinberg sign convention for the Minkowski tensor, one has [34]: As can be noticed from Equation (10) and Figure 1, the maximum of the antenna pattern function occurs for θ = π, which is when the GW source and the MSP are aligned with the observer. Here, there is a removable discontinuity, that can be fixed by replacing F A (π) with lim θ→π F A . Figure 1. Plot of the antenna pattern function for the + polarization. The function is taken as the absolute value and is plotted in Cartesian coordinates x, y, and z. The z axis is oriented along the direction of the GW source. The color indicates the magnitude of the function, as described by the color bar. The antenna pattern function for the × polarization has the same shape, but is rotated by ∆φ = π/4.

Continuous Waves
By now, it is well known that almost all bright galaxies harbor a supermassive black hole (SMBH) [35]. Observations on several stars close to the center of the Milky Way, for example, confirmed the presence of a SMBH of about 4 × 10 6 M [36,37]. So, massive black holes cannot be the remnant of stellar evolution; their formation likely started in the early Universe from smaller BH seeds that evolve through accretion and successive merging events. It is thought that BH seeds may be the remnant of Population III stars that, due to their low metallicity, could grow up to one hundred solar masses [38], or else are primordial BHs formed during the early stages of the Universe [39]. It is also possible that BH seeds are formed directly by the collapse of large gas clouds [40]. In any case, it is expected that the formation of SMBHBs occurs when two galaxies merge. In this scenario, the whole system loses energy due to dynamical friction, and this process reduces the separation between the two SMBHs at the center of the galaxy. Dynamical friction energy loss becomes less effective as the SMBHB orbital period shortens and eventually stops when it becomes smaller than the interaction time of the galaxy stars. At this stage, the GW emission is negligible, so in the absence of other energy subtraction processes, the SMBHB would stall indefinitely. This is known as the last parsec problem because it is still not clear which mechanism is responsible for the shrinking of the SMBHB down to sub-parsec scales, where the GW emission starts to dominate. Once this phase is reached, the SMBHB emits continuous GWs with a frequency of the order of 10 −9 Hz. Such GWs induce a space-time perturbation that, within the context of GR, can be described by: where i = √ −1 denotes the imaginary unit, h A cgw and α A are, respectively, the amplitude and the initial phase of the A-th polarization. Notice that it has been specified that one must take the real part (indicated by Re) of the right member since h A is a physical quantity. Here, it is important to remark that the GW frequency is not a constant. It can be shown that for an SMBHB that has been circularized by GW emission, the emitted GW frequency is about two times the size of the orbital frequency f orb of the system. Hence, accounting for the cosmological redshift, the observed GW frequency f is given by: The energy loss due to GW emission makes the SMBHB shrink so that the frequency of the emitted GWs grows with time. It can be shown that the variation of the orbital frequency is given by [41]: where µ is the reduced mass of the SMBHB and M is the total mass of the binary. Equation (13) implies that an SMBHB in the ultra-low-frequency GWs emission regime evolves very slowly. For this reason, the frequency of a GW that passes through the MSPs and the Earth can be considered, in good approximation, the same. Combining Equations (6), (8) and (11), one obtains: The function in Equation (14) is periodic with respect to both variables θ and t and, as can be seen from Figure 2, it rapidly oscillates with respect to θ. Moreover, due to the GW transverse nature, it is null when the GW source and the MSP are aligned with the observer, and assumes its global maximum for θ π. For a given value of θ, as can be seen from Figure  3, the GW-induced timing residuals show a sinusoidal behavior. The GW angular frequency ω determines its period, while the factor h A cgw /ω gives an upper limit on timing residuals that can be induced by a GW emitted by an SMBHB. For ω 10 −8 Hz, which is in the frequency range of GWs detectable by current PTAs, the typical value of the GW amplitude turns out to be about 10 −16 . Therefore, one has h A cgw /ω 10 ns. Unfortunately, there are very few MSPs that can be timed with such high precision, so at the moment it is unlikely to detect single SMBHBs by pulsar timing.

Gravitational-Wave Background
The ensemble of GWs emitted by a large number of independent sources is expected to form the so-called stochastic gravitational-wave background (GWB). Like the Cosmic Microwave Background (CMB), GWB carries a lot of information about the early Universe and possibly allows us to extend our observations beyond the last-scattering surface, from where the electromagnetic radiation could not escape. One of the largest contributions to GWB should be due to the population of SMBHBs. Since each of these BH binaries is in a different evolution stage and a different position, the pulse redshift associated with the ensemble can be obtained by integrating Equation (6) with respect to frequency and direction. Thus, using the GW amplitude expressed in Equation (11), this leads to: A useful way to characterize statistically any stochastic background is by computing the ensemble average of its components, which are random variables. This quantity is equivalent to the time average only if the ergodic hypothesis holds. In the case of the GWB, some well-motivated assumptions can be made to simplify the calculation [42]. For a stationary, Gaussian, isotropic and unpolarized GWB, one has: where S h ( f ) is the spectral density of the GWB [30]. Equations (15) and (16) allow us to determine the ensemble average of the pulse redshift of a pair of MSPs, labeled with a and b respectively, which is given by: where the asterisk indicates the operation of complex conjugation. Notice that the term in the square bracket in Equation (15), does not appear in Equation (17). Indeed, for MSPs typically used in PTAs, t p is larger with respect to the inverse of ω and, for this reason, the complex exponentials oscillate very rapidly with the frequency. Therefore, in this limit, their overall contribution to the integral can be neglected. The second integral in Equation (17) is the cross correlation between the antenna pattern function of the MSPs. The estimation of this integral can be performed by adopting a more convenient reference frame. If the z axis is chosen to be aligned to one MSP, it results in: where ξ is the angular separation between the MSPs and m i and n i are orthogonal versors lying in the plane perpendicular to Ω i . Since the polarization tensors are given by: the antenna pattern function can be written explicitly, in the adopted reference frame, using Equations (18) and (19). Then, from Equation (17) one obtains: where the term: is known as the Hellings and Downs function [43], which is plotted in Figure 4. The ensemble average between the timing residuals of a pair of MSPs can be determined following the same approach. Since timing residuals are directly related to the ToAs of the pulses, which are actually observable, this quantity is much more useful in searching for GWB evidence. From Equation (14), one can derive: The Hellings and Downs function is considered the "smoking gun" of the GWB because it may allow the distinguishing of quadrupolar correlations, due to the GWs, from monopolar or dipolar correlations [29], which can be due to the choice of reference clocks and to misplacement of the position of SSB, respectively.
It can be shown that the characteristic strain spectrum of the GWB is well approximated by a power law: where h c ( f ) is the characteristic strain amplitude, which is linked to the spectral density by [44], and α is the spectral index which, in this case, is α = 2/3.

Gravitational-Wave Bursts with and without Memory
A sudden and temporary variation of the quadrupole moment of a system may lead to impulsive space-time perturbations, known as GW bursts. Possible sources of GW bursts are then SMBHs in hyperbolic or highly eccentric elliptic orbits [45]. While in the former case there are multiple GW emissions, which mainly occur close to each periastron passage, in the latter case there is a single GW emission concentrated mostly at the distance of the closest approach; therefore, these events are much harder to detect. Even if PTAs are sensitive to GW bursts, they must last at least some weeks to be detectable, because of the limit imposed by the cadence of observations. As a result of extremely energetic events, GW bursts may lead to a permanent space-time modification, and for this reason, these are referred to as GW bursts with memory (BWM). Possible sources of BWM are then core-collapse supernovae, where the BWM is produced by the asymmetric ejection of a large amount of matter and radiation [46], or SMBHBs, where the BWM is associated with the strong and anisotropic GW emission in the final stage of the merger [47]. The upper limit on the BWM strain, which can be estimated by the ratio between its Schwarzschild radius and its distance from the observer [48], turns out to be of the order of 10 −22 for core-collapse supernovae and 10 −15 for SMBHBs.
Analogously to the case of continuous GWs, BWM can also induce timing residuals. The space-time perturbation can be modeled as a step function: where h bwm is the BWM amplitude, t b is the observation time of the BWM and Θ(t − t b ) is the Heaviside function, which is null before the burst (i.e., for t ≤ t b ) and it is unitary after it (i.e., for t > t b ). Taking into account Equations (6) and (24), by integrating Equation (8) results [49]: where t 0 = t b and t 1 = t b − ωt p (1 + cos θ), and therefore, the two terms in the right member are, respectively, the Earth term and the pulsar term. The BWM-induced timing residuals, expressed in Equation (25), can be rewritten in a simpler form by neglecting the pulsar term because it gives rise to an uncorrelated contribution in the timing residuals of each MSP, and is then treated as a source of stochastic noise in the whole PTA analysis [49,50]. Moreover, it has to be considered that uncertainties in the MSP spin-period and spin-period first-derivative lead, respectively, to linear and quadratic contributions in the timing residuals [51]. Therefore, after removing them, the GW-induced timing residuals for a BWM become: where I(t) = (t − t b )Θ(t − t b ) − I qdr (t) is the Earth term subtracted by the quadratic term I qdr (t) = a(t − t b ) 2 + b(t − t b ) + c, plotted in Figure 5. The coefficients a, b and c are chosen to minimize the root mean square (RMS) of I(t).

Cosmological Origin of the Cosmic Strings
According to the Cosmological Standard Model, the Universe began with an immediate expansion, referred to as the Big Bang, from an initial singularity characterized by extreme values of density and temperature. Right after the Big Bang, an accelerated expansion epoch started called inflation, which lasted for about 10 −32 s, during which the Universe's size increased due to the inflation field, and nowadays it is still increasing due to the dark energy, which became relevant about 10 Gyr after the Big Bang. Since the Universe can be considered as an isolated system, its expansion is adiabatic and led to a decrease in temperature. Therefore, the Universe, at some point in its evolution, may have undergone a phase transition, during which some of the existing symmetries may have been broken and, as a consequence, topological defects might have arisen. In the case of the U(1) symmetry, these structures are one dimensional and, for this reason, take the name of cosmic strings.

Cosmic Strings from U(1) Spontaneous Symmetry Breaking
A possible origin of the cosmic strings is the spontaneous symmetry breaking of a theory characterized by a U(1) global symmetry and Lagrangian density given by: where φ is the only scalar field of the theory and V(φ) is a Mexican hat potential: where λ is a coupling constant and η is the vacuum state of the field φ, corresponding to the minimum of its potential V(φ). For convenience, in this section, natural units c =h = 1 have been adopted. The equations of motion for the density Lagrangian in Equation (27), which can be derived from the Euler-Lagrange equations, results: where = ∂ µ ∂ µ denotes the d'Alembert operator and m s = √ λη. Equation (29) can be solved by expressing the scalar field φ as: where the cylindrical coordinate system (ρ, θ, z) has been adopted, f is a function and n is an integer. From Equation (29), both a massless scalar field, associated with the spontaneous symmetry breaking (i.e., is the Goldstone boson) and a massive scalar field with mass m s , arise. Since the function depends on the product between the mass of the scalar field φ and the radial distance, it is convenient to define ξ = m s ρ. By substituting Equation (30) in Equation (29), one has: Equation (31) implies that the function f , and hence the scalar field φ, has two different limits depending on ξ. For ξ → 0, one must have f → 0 to guarantee the validity of Equation (31). For ξ 1, instead, from Equation (31) results f → 1. The energy density associated with the scalar field φ in Equation (30) is: The function f determines the behavior of the energy density in Equation (32). Indeed, for ρ m −1 s , one has E ∝ ρ −2 . The energy density is hence confined within a distance ρ m −1 s , and since this is true ∀z ∈ R, that defines a string [42]. The integer n can be, therefore, interpreted as the winding number of the string. The energy density can be integrated over the element of surface dσ = ρdρdθ, to obtain the string tension µ. Since for ρ → ∞ the integral diverges, it must be performed within a cutoff radius r cut : µ = πn 2 η 2 log(m s r cut ) (33) A different origin of the cosmic strings is the spontaneous symmetry breaking of a theory characterized by a U(1) local symmetry and an Abelian-Higgs Lagrangian density, given by: where F µν = ∂ µ A ν − ∂ ν A µ is the curl of the vector field A µ , D µ = ∂ µ + igA µ is the covariant derivative operator, g is a coupling constant, and V(φ) is the same Mexican hat potential introduced above (see Equation (28)). Therefore, in addition to a massive scalar field with mass m s = √ λη, which is present also in the U(1) global symmetry, there is a massive vector field with mass m v = gη, while is absent the massless scalar field, which is absorbed by the vector field and it is responsible for its mass. As for the U(1) global symmetry, the equation of motion can be derived from the Lagrangian density in Equation (34), and it can be eventually solved by expressing the scalar and the vector fields as: where f and a are two functions andθ i is the versor in the direction of the vector field. The solutions of the equation of motion are not trivial, but it can be found that asymptotically f → 1 and a → 1. Again, it can be shown that the energy density is localized in a string [52].
The main difference with the U(1) global symmetry is that, for ρ → ∞, the integral is finite, and the string tension, for n = 1, is: where β = m v /m s , and g(β) is a slowly varying function that satisfies the condition g(1) → 1.

Gravitational-Wave Emission from Cosmic Strings
Frequently, cosmic strings interact with each other or with themselves. From these processes, new structures, such as loops and kinks, might arise. Loops are cosmic strings in a ring-like configuration rapidly oscillating due to the high string tension and radiate energy via gravitational waves thereby shrinking and, eventually, decaying. During loop oscillations, the string can double back on itself, producing a cusp [53]. Kinks are cosmic strings in an open configuration characterized by a vertex. GW emission from cosmic strings might form a GWB. For a string network, the characteristic strain spectrum of the GWB can be described by a power law with a spectral index α = −7/6 [54].
Cusps and kinks might emit GW bursts strong enough to stand above the GWB. In particular, cusps are expected to be powerful GW burst sources [8]. The correspondent space-time perturbation can be written as: where w is the duration of the GW burst emitted by the cusp, t 0 is the time relative to its peak, and h csb is its amplitude, given by: where Γ is the Euler function, l is the string length, and r(z) is the comoving distance at redshift z [8] (see Figure 6). The string length places a limit on GW frequency f , indeed, since the wavelength of the emitted GW cannot be greater than the string length, f l −1 . Therefore, PTAs can be used to detect GW burst from cosmic string with a length in the range [10 14 , 10 17 ] m. Within the context of Grand Unified Theory, cosmic strings should be characterized by a string tension µ 10 −6 G −1 (i.e., µ 10 22 kg m −1 ) [55]. This is a huge value, just think that a cosmic string with a length of 10 pc should have a mass of 10 9 M concentrated in a narrow filament, with a thickness, determined by m −1 s and m −1 v (see Section 5.2), that can be smaller than the proton radius [42]. The GW-induced timing residuals for a GW burst due to a cusp can be found using Equations (7), (8), (37) and (38) [56]. For simplicity, the pulsar term has been neglected, as was the case in in Section 4.3. Due to the presence of the Heaviside function in Equation (37), the result is a piecewise function expressed by: where sign(t − t 0 ) is the sign function, which is positive if the argument is greater than zero and negative otherwise. Moreover, as explained in Section 4.3, the linear contribution in the timing residuals can be subtracted to remove the uncertainties in the MSP period (see Figure  7).

Figure 7.
Plot of the GW-induced timing residual function for a GW burst due to a cusp. The dotted line describes the timing residuals due to the linear term, the dashed and the solid lines describe, respectively, the GW-induced timing residuals for a GW burst due to a cusp before and after the linear term subtraction, represented by the dotted line. A cosmic string at z = 1000, with length l = 10 kpc and tension Gµ 10 −6 , emitting a GW burst in the direction (θ, φ) = (π/2, 0) with duration w = 15 years has been considered, reaching peak at t 0 = 15 years. The horizontal axis indicates the observation time, expressed in years. The vertical axis indicates the timing residuals, expressed in µs.

Current Results
Recently, significant steps towards the detection of ultra-low-frequency GWs by PTAs have been made. All the main PTA collaborations might have observed a common-spectrum process compatible with a GWB, described by a power law characterized by a spectral index α = −2/3 and a characteristic strain amplitude h c 10 −15 at a frequency f = 1 year −1 , produced either by a population of SMBHBs or by a cosmic string network [57][58][59]. However, it is important to remark that, according to the NANOGrav analysis, the no-correlations hypothesis is rejected only mildly, with a 2σ significance, which is lower than the 3σ and 5σ significance standards of particle physics, required to claim for evidence and detection, respectively. Note that this result does not refer to time correlation, but instead to spatial correlation. The latter, as illustrated in Section 4.2, must show the Hellings and Downs signature (i.e., it must be a quadrupolar correlation) in order to claim a GWB detection consistent with GR. For detailed information about the observed common-spectrum process and complete data analysis of the most recent data sets released by EPTA, NANOGrav, and PPTA collaborations, the reader is referred to Refs. [21,59,60].
In addition to SMBHBs and cosmic strings, there might be other sources of a commonspectrum process, which are worth considering. Some of the most suggestive and physically impactful, if confirmed, are inflation and reheating [61], BBHs from the previous aeon merged during the Big Crunch [62,63] in cyclic conformal cosmological models of the Universe [64], massive gravity [65], and axion-like particles [66]. In particular, a GWB produced by quantum fluctuations of the gravitational field during the inflationary epoch is actually predicted within the context of the standard cosmological model Λ-Cold Dark Matter (ΛCDM), so its discovery would be crucial for probing the evolution of the Universe. For a complete and detailed description of the inflationary epoch characteristic imprint expected on the observations and the upper limits that can be placed by PTAs, the reader is referred to Ref. [61].
The GW nature of the common-spectrum process can be proved only by observing the Hellings and Downs signature in the cross-correlated timing residuals. While the GW detection cannot be claimed with the data currently available, significant constraints have been placed. For example, assuming a GWB produced by a population of SMBHBs, the merger rate has been constrained in the range [1.2, 45] × 10 −5 Mpc −3 Gyr −1 and has been placed a lower limit on the merger timescale, which is 2.7 Gyr. According to these results, the formation and the merging of a significant number of SMBHBs would be allowed [67]. Assuming instead a GWB produced by a cosmic string network, the string tension has been constrained in the range [4,10] × 10 −11 G −1 [68]. This improved by almost four orders of magnitude the upper limit placed by the Planck data on the CMB [69]. Cosmic strings with such a low string tension would not emit GW bursts detectable by ground-based interferometers or PTAs [70]. It is worth mentioning that some very recent analyzes have highlighted the existence of some "tension" between the estimates on the string tension resulting from the data collected respectively by the NANOgrav and PPTA collaboration. In fact, while both could be compatible with a stochastic gravitational-wave background produced by cosmic strings, the string tension is constrained to be in the range [2,30] × 10 −11 G −1 by the NANOgrav data and 4 × 10 −11 G −1 by the PPTA Collaboration. However, it has been recently shown that this tension can be somewhat relaxed by invoking primordial Coleman-Weinberg inflation which partially inflates the string network [71].
PTA data can also be used in multi-messenger astronomy to confirm the presence of an SMBHB (in galaxies centers), showing a periodicity in the electromagnetic observations. Although an unambiguous GW detection has not yet been accomplished, PTAs can be used to constrain the physical parameters of a possible SMBHB hosted by a galaxy. For instance, one of the most intriguing galaxies suspected of hosting an SMBHB within its core is 3C 66B [72]. The early constraint on its chirp mass, obtained using long-baseline interferometry, has been recently improved by almost one order of magnitude using PTA data, placing an upper limit M (1.65 ± 0.02) × 10 9 M [73]. This is still compatible with the SMBHB hypothesis, and hopefully, soon, PTAs will be able to either confirm or rule it out.

Conclusions and Further Perspectives
PTAs are very promising tools for ultra-low-frequency GW detection. They can be used to observe sources, such as SMBHBs and cosmic strings, by detecting continuous GW, GW bursts, and the GWB. This would allow addressing some extremely relevant astrophysical and cosmological questions about the Universe. Although PTAs already give interesting constraints on the GW sources above (see Section 6), no GW detection has been claimed yet. However, the latest results of the PTA major collaborations suggest that the first GW detection is not that far from present times. Next-generation radio telescopes, as the Five hundred-meter Aperture Spherical Telescope (FAST) [74], the MeerKAT radio telescope [75], and the Square Kilometre Array (SKA) [76], will improve the timing precision. Therefore, in the next decade, the PTAs will probably be used as ultra-low-frequency GW detectors.
Eventually, it has to be said that even with better detectors and larger data sets, the possibility of having no GW detection cannot be excluded. There may be several explanations behind that case. One possibility is that there is not any signal to detect. For example, the last parsec problem (see Section 4.1) might prevent the merging of the SMBHBs, and cosmic strings (see Section 5) might not exist. Another possibility is that the GW detection method currently adopted by PTA, mainly based on the Hellings and Downs cross correlation, needs to be revised or supplemented.
One way to do this could be, for example, by including MSPs inside the core of GCs in PTAs. Indeed, since these MSPs are confined in such a tiny angular region, which radius can be of the order of some arcseconds, the angular separation between them is so small that they should have strongly correlated GW-induced timing residuals (see Figure 4). This property, once all the possible ToAs variations due to the GC have been taken into account in the timing model, makes them very useful for GW detection. For a detailed discussion on this novel and important issue, the reader is referred to Ref. [77].