Implication of Pulsar Timing Array Experiments on Cosmological Gravitational Wave Detection

Gravitational waves provide a new probe of the Universe which can reveal a number of cosmological and astrophysical phenomena that cannot be observed by electromagnetic waves. Different frequencies of gravitational waves are detected by different means. Among them, precision measurements of pulsar timing provides a natural detector for gravitational waves with light-year scale wavelengths. In this review, first a basic framework to detect a stochastic gravitational wave background using pulsar timing array is introduced, and then possible interpretations of the latest observational result of 12.5 year NANOGrav data are described.


INTRODUCTION
Gravitational waves are ripples of spacetime first predicted by Einstein in his theory of general relativity in 1916. In Newtonian theory of gravitation and mechanics, both space and time are rigid in the sense that they are not affected by any material content existing in the Universe, and gravity is a nonlocal force. In contrast, Einstein advocated that matter and energy make the spacetime curved and the curvature itself is the essence of gravitation.
In his theory a massive body curves the surrounding space and another particle moves in the curved space approaching the former as if it is directly attracted by the body. This is in close analogy of the electromagnetism of Faraday and Maxwell, where a moving charged particle produces electric and magnetic fields in the surrounding space, and another charged particle feels a force from these fields. In the case of gravity, the spacetime curvature is characterized by a metric tensor, g µν (x), with which the spacetime interval is expressed as ds 2 = g µν (x)dx µ dx ν where repeated Greek indices are summed over spacetime coordinates 0 to 3. Without any gravitational effect, g µν (x) coincides with the Minkowski metric η µν = diag(−1, 1, 1, 1), taking the speed of light unity.
Gravitational waves can be expressed by a small perturbation around the Minkowski metric in terms of the transverse-traceless gauge as which satisfies ∂ i h T T i j = 0, h T T i i = 0. Here repeated Latin indices are summed over the spatial coordinates 1 to 3. In the linearized weak-gravity limit of general relativity, the tensor perturbation h T T ij (t, x) satisfies the same wave equation as the electromagnetic field, which is also equivalent to the Klein-Gordon equation for a massless scalar field, so that the gravitational wave propagates with the speed of light.
One can easily see that propagating gravitational waves do not change the position of whatever body initially at rest whose motion obeys the geodesic equation: where τ is the proper time. Since Γ i 00 = 0 for the metric (1), a star or a particle initially at rest does not feel any acceleration by the tensor perturbation and so remains at rest. Thus the only effect gravitational wave causes is the change of the light path or distance between two objects at rest.
In the case of a ground-based detector such as the advanced LIGO (aLIGO) [1] in USA, advanced Virgo (aVirgo) [2] in Europe, and KAGRA [3] in Japan, two reflective laser lines are set with a rectangular configuration whose interference measures modulations of their relative path length to detect gravitational waves. Their baseline is 3-4 kilometers long and these detectors are sensitive to gravitational waves in deca-to hecto-Hertz band.
While they were originally designed to detect gravitational wave event from binary neutron star coalescence, the first direct detection of gravitational waves, which was achieved by aLIGO on September 17, 2015, was from coalescence of two massive black holes around 30 solar mass [4]. Up to then, even the existence of such a massive black hole had not been appreciated in the community. This manifests another example of common practice of the history of astronomy that every time the mankind obtained a new means of observation, a new unexpected event had been discovered.
The first discovery of gravitational waves from binary neutron star coalescence was made on August 17, 2017 [5] which led to joint observation of multiband photons ranging from gamma-ray to infrared [6]. We note that this epoch making discovery was successfully made only by the GstLAL pipeline jointly developed by Kipp Cannon at RESCEU and Chad Hanna at Penn State University among the six pipelines of gravitational wave burst detection at work.
In order to detect lower frequency gravitational waves by a laser interferometer, we must extend the baseline, which may be achieved by launching a set of three satellite to form a laser interferometer of triangular configuration. This idea was first put into serious investigation by the Laser Interferometer Space Antenna (LISA) project which aims at launching three spacecrafts separated by 2.5 million km in a triangular formation orbiting around the Sun following the Earth in the same orbit [7]. It will have the best sensitivity at milli-Hertz frequency range and can probe black holes with much larger masses than those discovered by aLIGO and aVirgo in the range 10 2−7 solar mass. The path-finder satellite was launched in December 3, 2015 which was one day after the centennial of the publication of the Einstein equation [8]. The satellite exhibited much better performance than originally planned providing a hope to launch the full mission in 2034 [9].
DECIGO in Japan, on the other hand, sets a goal of the direct detection of the stochastic gravitational wave background which was produced quantum mechanically during inflation in the early Universe, observing a decihertz range of gravitational waves [10]. In fact, DECIGO will ultimately be able to measure the reheating temperature after inflation which indicates when the Big Bang happened [11]. While it has been shown that the baseline should be as long as 1500km with a Fabry-Perot cavity [12], which induces a resonance of laser by reflecting the light many times by mirrors at both ends of the cavity, current plan under discussion is to realize its degraded version called B-DECIGO whose baseline is only 150km first which is still in a phase of conceptual study.
In contrast, Chinese space-based projects are advancing much more rapidly with two independent projects of space-based gravitational wave detection, namely, Taiji [13] and TianQin [14]. Both projects have successfully launched the initial path finder satellites already.
The scope of this article is to focus on a set of gravitational wave detectors with by far the longer baseline prepared by nature, namely, timing observation of millisecond pulsars, or pulsar timing arrays.
Pulsars are rotating neutron stars with a strong magnetic field whose axis is not identical to the rotation axis, so that pulse-like periodic electromagnetic radiation is observed with the rotational period. So far about 2.5 thousand pulsars have been discovered since its first discovery in 1967 by Bell and Hewish [15], and their rotation period spans from a millisecond to about ten seconds. Among them, the rotation period of millisecond pulsars with their period less than 30 msec is especially stable, most of which exist in a binary system and have a relatively weak magnetic field of order of 10 8 G [16]. They are so stable that by observing possible modulation of the arrival times of their pulses, we can in principle detect stochastic gravitational wave background using observational timing data of many pulsars.
Let us first summarize how the time of arrival (TOA) of the pulsar's signal is modulated due to the propagating gravitational waves. For more detailed introduction, one may consult e.g. [17].

Effects of gravitational waves on the timing of a single pulsar
Suppose that an observer is at rest at the origin of the spacetime (1) and receive a signal emitted at t em from the a-th pulsar at the distance d a and directional cosine n a . The observer receives the signal at Repeating the same calculation one cycle (T a ) later and subtracting (2), the modulation of pulsar period observed is given by where the integral is calculated along the unperturbed photon path For a plane gravitational wave propagating along the direction n, with τ a ≡ t obs − t em . In general, it can be expressed as The first term in the right-hand side is called the earth term and the second term the pulsar term. In the actual data analysis, the earth term is calculated at the barycenter of the solar system to remove seasonal and other modulations associated with the dynamics of the solar system. Choosing the coordinate frame so that the gravitational wave propagates along the z direction, n = (0, 0, 1), h T T ij reads If the a-th pulsar is located in the angular direction in terms of the polar coordinate, we find The dependence on 2φ a manifests the spin 2 nature of the graviton. We now apply it to the case of isotropic stochastic background of gravitational wave radiation where h T T ij (x, t)is expressed in terms of the stochastic mode functionh A (f, n) as where e A ij (n) is the symmetric polarization tensor which satisfies n i e A ij (n) = 0 and e A ij (n)e A ij (n) = 2δ AA . We will return to their explicit form shortly. The stochastic background is characterized by the power spectrum S h (f ) which satisfies The modulation of pulsar's period is then expressed as with F A a (n) = n i a n j a e A ij (n) 2(1 + n · n a ) .
The first term in the square bracket represents the earth term of (5). The primary quantity of observational importance is the time integral of (5) and (10), r a (t) ≡ t tini z a (t )dt , which is called the timing residual.

Effects of gravitational waves on the timing of a pair of pulsars
In order to detect stochastic gravitational wave background in terms of pulsar timing array, it is important to study the correlation of the timing residual of two or more pulsars. The primary quantity is the equal-time correlation function of the timing residual of pulsars a and b, which is well approximated by unity as the other terms are highly oscillatory. Let us express the directional vector of gravitational waves as n = (sin θ cos φ, sin θ sin φ, cos θ).
Then introducing orthogonal vectors u = (sin φ, − cos φ, 0) and v = (cos θ cos φ, cos θ sin φ, − sin θ), we can express the polarization tensors as with which we find Thanks to the assumed isotropy of the stochastic background, (11) depends only on the relative angle between pulsars a and b, θ ab . Then we may choose the z axis along the direction to the a-th pulsar, so that n a = (0, 0, 1) and take n b = (sin θ ab , 0, cos θ ab ). Inserting them to (11) to calculate (10) we find where The most important feature of (11) is the fact that the form of the angular part is independent of the spectrum of gravitational radiation, S h (f ),which was first found by Hellings and Downs [18] and the functional form of C (θ ab ) is often called the Hellings-Downs curve.
The correlation function of the timing residual reads We can extend the above to the case of unequal time but the result is always proportional to the Hellings-Downs function C (θ ab ) which reflects the quadrupole nature of gravitational radiation.

Characteristic strain and density parameter of gravitational wave background
The characteristic strain, h c (f ), is defined by It is often modeled by a single power law: where it is convenient to take the characteristic frequency at f yr = 1yr −1 when we analyze observation using pulsar timing.
To discuss cosmological implication of stochastic gravitational wave background, it is convenient to express the energy density of gravitational waves, which is naturally defined by in terms of its contribution to the density parameter per logarithmic frequency interval, Ω gw (f ). Since (20) is expressed as Ω gw (f ) reads where ρ cr = 3H 2 0 (8πG) is the critical density with H 0 being the current Hubble parameter.

Observables and data analysis
The observed quantity is a time sequence of timing residuals of all the observed pulsars denoted as r a (t i ), a = 1, 2, ...N p ; i = 1, 2, 3, ...with N p being the total number of pulsars regularly monitored. The observables may be decomposed into the deterministic part and noise part as where collectively denotes the parameters characterizing the nature of the a-th pulsar and other deterministic parameters modeling the observational setup [19].
Assuming that noises are Gaussian distributed, we can perform statistical analysis once the two-point correlation function of noises is obtained. The auto-correlation function of noises of the same pulsar, N aa;ij = n a (t i )n a (t j ) , consists of a frequency-independent white noise part such as instrumental errors and frequency dependent red noise term with excess power at lower frequencies such as spin noise, pulse profile changes, and imperfectly modeled dispersion measure variations. The stochastic gravitational wave background also contributes to the latter. Noises acting on different pulsars at different time, N ab;ij = n a (t i )n b (t j ) , is subdominant compared with N aa;ij but plays an important role to detect the stochastic background of gravitational waves, as their effect appear as its quadrupole term, while their monopole and dipole terms are mainly due to the errors in the reference clock on the earth and solar system modeling, respectively.
In order to measure the amplitude of stochastic gravitational wave background, we first marginalize over the pulsar parameters ξ a assuming that gravitational waves are absent and then write down the likelihood function incorporating it. The log-likelihood function has a form wherer a andN collectively denotes the sequence of observed data with the deterministic part subtracted, r a (t i ) − r det a (t i , ξ a ), and the covariance matrix of noises both after marginalization of deterministic parameters, respectively.
We assumer a andr b with a = bhas a vanishing correlation apart from the quadrupolar one due to the gravitational waves. We may express the marginalized noise covariance matrix asN = P + A 2 GW S where the matrix S is estimated as S ab = r a tr b with temporal coefficients being suppressed [20]. Since A GW is a tiny quantity if not zero, the inverse matrix is easily obtained as We therefore find the log-likelihood ratio as We thus find the optimal statistic to calculate the amplitude of gravitational waves asÂ which yields Â 2 = A 2 GW . If the signal is weak, the standard deviation is given by , so that the signal-to-noise ratio reads S/N =Â 2 /σ 0 [21] [22].
In order to judge if nonvanishing stochastic gravitational wave background is present, one may adopt the Baysian analysis to obtainP (M i |d), the probability of the model M i being correct for a given set of observational data d. The Bayes theorem tells us where B ij is the Bayes factor and the last equality applies for the flat prior case. If it is larger than unity the data prefers model i to j, and its degree of preference has been classified by Jeffreys as follows. For log 10 B ij = 0 − 0.5, the preference is not worth more than a bare mention, for log 10 B ij = 0.5 − 1, it is substantial, for log 10 B ij = 1 − 2, strong, and finally for log 10 B ij > 2, it is decisive [23].
As the observation continues, sensitivity to the gravitational wave background improves and we obtain more stringent upper bound on its amplitude with the frequency range extending toward lower frequencies. If there exist a nonvanishing background, we would first encounter stagnation of improvement in the upper limit, followed by emergence of spatially uncorrelated commonspectrum red noises to all the pulsars being monitored. Finally, we would find the quadrupolar signature in the spatial correlation described by the Hellings-Downs curve, which would serve as the final confirmation of the very existence of the gravitational wave background.

NANOGrav Observation
There are three major pulsar timing array [24] (PTA) experiments worldwide, the European Pulsar Timing Array (EPTA) [25], Parkes Pulsar Timing Array (PPTA) [26], and the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) [27], which constitute the International Pulsar Timing Array (IPTA) [28]. Accumulation of observation for more than a decade now, the upper bound on the amplitude of the stochastic gravitational wave background has been improving to the frequency reaching . Recently, however, the 12.5 year NANOGrav data [29] has found a strong evidence of a common-spectrum stochastic process [30].
The data consists of observation of 47 millisecond pulsars made from July 2004 to June 2017 by 305m Arecibo Observatory in Puerto Rico and 100m Green Bank Telescope in West Virginia [29], and 45 pulsars observed for more than three years among them were used in the gravitational wave analysis.

Latest Result of NANOGrav Observation
The 12.5year NANOGrav data set shows possible existence of common-spectrum process as depicted in Fig. 3 and fitted by a power law model in analogy with (15): Figures 3 and 4 indicate evidence for common red process at lower frequencies. Indeed the Bayes factor (BF ) of the model with and without it reads log 10 BF = 4.5 for the case the solar system ephemeris is fixed to DE436 provided by Jet Propulsion Laboratory and log 10 BF = 2.4 for the case uncertainties in solar system ephemeris are marginalized using an algorithm called BayesEphem. The latter provides a more conservative result but in both cases, the existence of red noise is decisive. The next question is if the correlation among different pulsars satisfies the Hellings-Downs curve, or has monopolar or dipolar features.  The monopolar correlation is strongly disfavored as log 10 BF = −2.3 with DE438 and −1.3 with BayesEphem. The dipolar correlation is also strongly disfavored as log 10 BF = −2.4 with DE438 and −2.3 with BayesEphem. On the other hand, the Bayes factor of the Hellings-Downs curve has been calculated as log 10 BF = 0.64 with DE438 and 0.37 with BayesEphem, and mixtures with either common, monopole, or dipole process do not yield any higher BF. Hence it is premature to conclude that they have observed stochastic gravitational wave background yet.  (16) for three types of correlations as well as the common-spectrum amplitudeA 2 CP . Although the monopolar distribution shows somewhat higher S/N ratio, the central value is lower than the case of Hellings-Downs correlation, and so the result is consistent with the Bays factor analysis.
Thus their result as a whole indicates a possible hint of the existence of stochastic gravitational wave background but by no means confirms it at this stage. We should still wait for the good news patiently. It is remarkable, however, a number of hasty theorists have already written many papers on many kinds of exotic interpretations of the result assuming that the result actually originates from gravitational wave background besides the more conventional possibility that it is from supermassive black hole binaries. In the next section we overview some of them in turn.

Generation of gravitational waves
Before discussing specific sources of gravitational wave background, let us introduce the generic quadrupole formula to provide the order-of-magnitude relation between the parameters of the source and the resultant amplitude of gravitational radiation. In this subsection we recover the dependence on the speed of light c. The linearized Einstein equation reads with being the energy-momentum tensor. The solution can easily be found with the retarded Green function as with the electromagnetism as Thanks to the conservation law, ∂ µ T µν = 0, one can express it in terms of the quadrupole moment as Thus in terms of the mass, M , size, R, time scale, T , of the source, we can express it as where r g is the Schwarzschild radius of the source.
To discuss gravitational wave background of cosmological origin, the D'Alembert operator of (29) must be replaced by the one in the Lemaitre-Robertson-Walker metric, and we must take relativistic effects in the energy momentum tensor into account to solve the Einstein equation instead of adopting slow-motion approximation used in (31).

Binary supermassive black holes (SMBH)
It is known that supermassive black holes exist ubiquitously in galactic nuclei [31] and in their evolution history it is expected that a number of coalescence took place of binary black holes. Their inspiral motion generates gravitational wave background which may be observed by pulsar timing [32]. The spectrum of stochastic gravitational wave background generated by their superposition was calcualted by Enoki et al [33] based on a semianalytic model of galaxy and quasar formation based on the hierarchical clustering scenario.
Suppose a pair of black holes with masses M 1 and M 2 forming a binary system in a circular orbit whose orbital radius is R. The frequency, f r , of the gravitational waves emitted by this system is twice the orbital frequency, so from the Kepler's law it satisfies f 2 r = π 2 GM R −3 where M is the total mass. From the quadrupole formula (31), we can estimate the amplitude of gravitational waves observed at a distance r as where µ is the reduced mass and M c is the chirp mass defined by M c ≡ µ 3/5 M 2/5 . Then the power emitted by gravitational radiation scales as which also gives the energy loss rate of the binary system. Since the orbital energy is given by , The observed frequency of gravitational waves emitted by a binary system at the redshift z is f = (1 + z) −1 f r , so the total energy density spectrum received by an observer today is  This simple dependence is a nice feature of the binary SMBH scenario to identify the origin of the stochastic background on one hand, but this also means that we can extract information on their mass function only through the overall amplitude , which indicates it is difficult to translate the result to the differential mass spectrum in this simple treatment with circular orbits.
In [36], two models of SMBH evolution have been compared with the NANOGrav 12.5 year data. One is the model studied by Middleton et al (referred to M16), where the mass function of binary black hole system is assumed to be described by a Schechter function in terms of the chirp mass with the redshift evolution also having a Schechter type dependence on with different sets of parameters [37]. Their model has five parameters, namely, power index and exponential scale height for both chirp mass and the redshift to determine the shape of Schechter function as well as the overall amplitude.
The other model studied by Chen et al (referred to C19) is more sophisticated as it incorporates eccentricity of the orbit and interaction with environment. It depends on 18 model parameters and has more complicated frequency dependence than a simple power law [38]. Figure 6 [36] represents the redshift-integrated mass function of the SMBH chirp mass fitted to these two models, showing that the mass function in the simpler model with a single power-law frequency dependence is less constrained than the more sophisticated one.

Cosmic strings
The power index of the gravitational wave background inferred by 12.5 year NANOGrav is also consistent with smaller values than α = −2/3 which may be preferred by cosmological gravitational wave background.
Among them, cosmic strings are line-like topological defects created by a cosmological phase transition when the vacuum state after the symmetry breaking allows a nontrivial mapping to a unit circle such as the case U(1) symmetry is broken [39]. The simplest field theoretic model allowing a string solution is the Abelian Higgs model described by the Lagrangian with a complex scalar field Φ, The cosmic string is characterized by the dimensionless tension Gµ where µis related to the symmetry breaking scale of the theory by µ = bπv 2 with b being a function of g 2 λ [40]. It takes a value between 0.5 and 3 as g 2 λ is varied from 100 to 0.01 [41]. Production of cosmic strings is usually attributed to a thermal phase transition which takes place due to the high-temperature correction to the potential [42]. The critical temperature of the above Abelian Higgs model is which is fairly high compared with typical reheating temperature after cosmic inflation [43]. For models with lower reheating temperature, one can still produce strings after inflation if is nonminimally coupled to scalar curvature R with the term ξRΦΦ † [44]. Then symmetry is restored during inflation if 12ξH 2 > λv 2 [45], or the tensor-to-scalar ratio satisfies It is known that 80% of the energy of strings are in infinitely long ones and they evolve in the expanding universe intersecting with each other to form loops [46]. As a result they evolve according to the scaling solution in which there always exists fixed number of infinitely long strings in the horizon volume at each time and loops are formed with the constant comoving rate with the initial length proportional to the formation time as αt. Numerical simulations show that has an extended spectrum up to α ≈ 0.1 [47]. We take α = 0.1 hereafter as corresponding loops make most important contributions to the gravitational wave background. Loops oscillate and decay by emitting gravitational waves with an energy-loss rate independent of its size,Ė = Γ Gµ 2 , where Γ is a constant around 50 [48]. As a result, a loop formed at t i has a length (t) = αt i − Γ Gµ(t − t i )at t and emits gravitational wave with a frequency f (t) = 2k/ (t), where k is a positive integer. It is known that k-th harmonics has an emission rate proportional to k −4/3 [49].
Noting that the string loop contributing to the current gravitational wave frequency f through k-mode emitted at timet was created at and summing up all the contribution throughout the cosmic evolution, we find the current spectrum of stochastic gravitational wave background where t f is the onset of the scaling regime. C ef f is a factor representing the dependence of the scaling solution on the cosmic expansion law [50]. We find C ef f = 5.4(0.39) in radiation (matter) dominated regime. Figure 7 shows the energy spectrum of stochastic background fitted to the 12.5 year NANOGrav data assuming that they are entirely from cosmic string Fig. 7 Spectrum of stochastic gravitational wave back ground generated by cosmic string networks in the scaling regime fitted to the 12.5 year NANOGrav data assuming that their signal entirely comes from strings. Taken from [51]. signals, where it is found that 68If we allow variation of α, even wider range of tension is consistent [52].
The stochastic background extended to higher frequencies as seen in Fig.  7 will cause big noise sources for the future space-based laser interferometers such as LISA, TianQin, Taiji, and DECIGO as well as to the Einstein Telescope [53]. Before reaching the era of the next and next-to-next detectors, however, once the world-wide network of currently available ground-based detectors, namely, aLIGO, aVirgo, and KAGRA (LVK) are in full operation with their respective design sensitivity, we can probe stochastic background to the level well below Ω gw (f ) < 10 −9 [54] and touch the "NANOGrav curve" in Fig. 7. Thus the cosmic string scenario may be verified even with the current technology.
Even if LVK did not observe the stochastic signal, however, it would not necessarily mean that cosmic string scenario would not work, because Fig.  7 is depicted under the assumption that waves relevant to LVK observation were created under radiation domination. If the reheating temperature after inflation is sufficiently low, they would have been created during the inflaton's oscillation regime when the cosmic expansion law was the same as in the matter domination, and the energy spectrum would acquire suppression factor proportional to f −2 [55]. I estimate the upper bound on the reheating temperature is around 400GeV in the case NANOGrav observation is entirely due to cosmic strings but LVK in full operation does not find any stochastic background.

Tensor perturbation from inflation
Cosmic inflation in the early universe is an indispensable ingredient of modern cosmology to realize big and old universe filled with hierarchical structures [43]. During inflation, each Fourier mode of tensor perturbation defined in (6), h A (f, n), satisfies the Klein-Gordon equation for a massless field, and can be quantized in the same way as a massless scalar field minimally coupled to gravity. As a result, it acquires nearly scale-invariant dispersion given by where the pre-factor on the left-hand side is the phase space density with k = 2πf and 64πGoriginates from the normalization factor of the Einstein-Hilbert action through canonical quantization procedure [56]. The spectral index of the tensor perturbation is defined by Here H(t k ) is the Hubble parameter during inflation when the comoving wavenumber k = 2πf left the Hubble horizon during inflation. In the exact de Sitter background we find a scale-invariant spectrum with n t = 0. In vast majority of inflation models n t takes a small negative value as n t > 0 means violation of null energy condition. Since CMB imposes a stringent constraint on the amplitude of gravitational waves or tensor perturbations on large scales we need a positive tensor spectral index if we wish to account for the NANOGrav observation in terms of inflationary tensor perturbation. At the same time we must suppress higher frequency stochastic background as aLIGO has already obtained a constraint more stringent than that imposed by successful big bang nucleosynthesis (BBN), namely, Ω gw < 6 × 10 −8 for a flat spectrum [58]. This can be achieved by prolonged reheating phase after inflation with a low reheating temperature as mentioned in the last subsection.
Parametric analysis blind to specific models has been done in [59] and we can read off range of parameters required to explain the NANOGrav data by focusing on the border line of pulsar constraints. As a result, the possible positive detection can be explained by an inflation model with a large tensor spectral index n t ≈ 0.9 and low reheating temperature Such a blue tensor spectrum may be realized in G-inflation [60] [61] [62] or Gauss-Bonnet inflation [63]. Both of them can be analyzed using the generalized G-inflation framework [64]. A relatively simple model to explain 12.5year NANOGrav data by tensor perturbation of quantum origin has been proposed by Tahara and Kobayashi making use of only G2 and G3 generalized galileon functions [65].

Tensor perturbation generated by second-order scalar perturbations from inflation
During inflationary expansion in the early universe, curvature perturbations are generated in a similar manner to tensor perturbation [66]. Although the simplest canonical slow-roll models predict nearly scale-invariant spectrum, nontrivial dependence of curvature perturbation on the scalar field model makes it possible to realize much larger amplitudes on smaller scale while preserving the consistency with CMB observations which determine the amplitude of the curvature perturbations of order of on large scale [57]. In particular, if slow-roll parameters satisfy a certain condition, the would-be decaying mode of the curvature perturbation can grow even after horizon exit, and can result in an enhanced spectrum on some specific scales [67]. If this happens, although tensor and scalar perturbations are decoupled at the level of linear perturbation, second-order effect may generate an appreciable amplitude of tensor perturbations [68] [69]. This feature, pioneered in [70] [71], has been extensively discussed recently in relation to the formation of primordial black holes (PBHs) [72] [73] which also requires high amplitude density perturbation of order of 0.1.
PBHs are produced when a region with a large density perturbation of order of unity falls in the horizon in the early universe dominated by radiation [74] [75]. Its mass is of order of the horizon mass at formation. Spectrum of density perturbation with a peak realizing formation of PBHs with mass creates stochastic gravitational wave background at the typical frequency Hence gravitational wave background of nano-hertz range is an important probe of PBHs around the solar mass. In particular, as the observation time gets longer, we will be able to constrain lower frequency background, so that we can eventually probe the mass range of the binary black holes observed by advanced LIGO and advanced Virgo. The possible positive detection of the background by NANOGrav collaboration has somewhat affected this strategy as a number of authors have proposed to attribute it to the positive signature of the existence of PBHs. Among them, Kohri and Terada [76] argue that NANOGrav may suggest the presence of PBHs around the solar mass, while Luca, Franciolini, and Riotto [77] suggest a spectrum of curvature perturbation predicting PBHs with a wide range of masses which as a whole constitute all the dark matter. Vaskonen and Veermäe [78] take into account the effect of critical phenomenon on the mass function of PBHs and obtained a comprehensive contour constraining the average mass and the abundance of PBHs in terms of the amplitude and slope of the gravitational wave background.
All these three works use a formula first obtained by Carr [79], which is often called Press-Schechter type formula even though it is not based on their theory [80], to relate the abundance of PBHs and the amplitude of density perturbations. It is desirable to reanalyze the problem under the improved mass function recently proposed [81] [82] including the dependence on the choice of the window function [83] [84].

First order phase transition
Quantum field theory with spontaneous symmetry breaking may induce a cosmological phase transition in the early universe due to the high temperature or high curvature effects. If a first-order phase transition takes place, it provides several sources of gravitational waves, namely, collisions of true vacuum bubbles, sound waves, and turbulences [85]. In order to produce gravitational waves at nano-hertz frequencies, the phase transition must have taken at a temperature around 10-100MeV. In the standard model of particle physics, only QCD phase transition is relevant in this energy scale, but it is already known to be a cross over without appreciable production of gravitational radiation. Hence new theoretical input such as a dark sector inducing a dark phase transition is required [86] [87].

Turbulence in magnetohydrodynamics
It is known that magnetohydrodynamics (MHD) turbulence in the early universe may induce stochastic gravitational wave background. For example, turbulence in the plasma and magnetic fields produced by a first-order phase transition at the electroweak scale may generate gravitational waves detectable by LISA [88]. In order to create gravitational waves relevant to NANOGrav observation in the same manner, we need a first-order QCD phase transition. Although such a possibility is interesting as we may relate the observation of gravitational waves with primordial magnetic fields, some new physics is required to realize such a phase transition [89].

CONCLUSION
In this article, I have briefly introduced a basic framework to observe gravitational waves by pulsar timing experiments and possible interpretations of the latest 12.5year observation of NANOGrav collaboration which reports a positive common spectrum signature that may be a hint of finite stochastic gravitational wave background in the nano-hertz frequency range. At present it is not guaranteed that the observed hint is due to the isotropic stochastic background of gravitational waves as the characteristic Hellings-Downs curve of angular correlation [18] has not been confirmed. Furthermore, most of the preferred amplitude of NANOGrav observation is in tension with the upper bound which has been obtained by PPTA, although they turn out to be consistent with each other allowing the uncertainties in ephemeris modeling. Hence we should continue careful study both on observational and theoretical analyses [90].