The Pulsar Timing Array Signal from Infrared Regions of Scalar-Induced Gravitational Waves

Abstract

The common-spectrum process, characterized by the Hellings–Downs angular correlation and observed by pulsar timing array collaborations, such as NANOGrav, PPTA, EPTA, and CPTA, can be explained by the scalar-induced gravitational waves (SIGWs). The energy density of SIGWs exhibits universal behavior in the infrared regions. Utilizing a broken power law parameterization for the primordial curvature power spectrum, we clarify the PTA signal through the infrared characteristics of the SIGWs, using Bayesian analysis to provide posterior distributions. Bayesian factors emphasize the statistical preference for the SIGW model over explanations involving supermassive black hole binaries.

1. Introduction

Linear scalar perturbations, seeded by the primordial curvature perturbations generated during inflation, can construct a second-order perturbation source to induce gravitational waves. These scalar-induced gravitational waves (SIGWs), rich in early Universe information, such as the primordial curvature perturbation at small scales [1,2] and the equation of state of the early Universe [3,4], span a broad frequency range and can be detected by pulsar timing arrays (PTAs) or space-based GW detectors, such as the Laser Interferometer Space Antenna (LISA) [5,6], Taiji [7], or TianQin [8,9]. To induce detectable GWs, the amplitude of the primordial curvature power spectrum should be enhanced by about seven orders of magnitude at small scales compared with large-scale observations, which can be achieved with ultra-slow-roll inflation models [10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34]. Alternatively, accompanied by the SIGWs, the large primordial curvature perturbations can give rise to primordial black holes (PBHs) [1,2,16,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69] (for the review, see [3,70]), offering explanations for LIGO-Virgo GW events [71,72], dark matter [73,74,75,76,77,78,79,80,81], and the anomalous orbits of trans-Neptunian objects [82].

Typically, scalar-induced gravitational waves exhibit a peak in the amplitude of the energy density at a specific scale, and this specific scale is denoted as $k_p$. In the infrared regions $k\ll k_p$, SIGWs resulting from primordial curvature perturbations with a narrow peak follow a universal behavior, specifically

$${\mathsf{\Omega }}_{\mathrm{GW}}\propto {\left(f/f_{\mathrm{year}}\right)}^3{log}_{10}^2(f_p/f),$$

(1)

where the term “universal” means that a broad range of models gives the same log-dependent index of the energy density spectrum of SIGWs in the infrared region [54,83,84,85]. A more effective approach to examining the energy density behavior of SIGWs in the infrared regions involves parameterizing the primordial curvature power spectrum using a broken power law form, ${\mathcal{P}}_{\zeta }=A_{\zeta }{(k/k_{p_0})}^{n_1}$ for $k<k_{p_0}$ and ${\mathcal{P}}_{\zeta }=A_{\zeta }{(k/k_{p_0})}^{n_2}$ for $k\ge k_{p_0}$. Within this broken power law framework, a narrow peak in the power spectrum is achieved for sufficiently large absolute values of $n_1$ and $n_2$, while a broad peak results for sufficiently small absolute values of these parameters. For the case of a narrow peak with $n_1>3/2$, the primordial curvature power spectra generate SIGWs with energy densities in the infrared regions as in Equation (1), which is independent of $n_1$ and consistent with the conclusion displayed in Ref. [54,83,84,85]. For the situation with $n_1=3/2$, the energy density in the infrared regions is ${\mathsf{\Omega }}_{\mathrm{GW}}\propto {\left(f/f_{\mathrm{year}}\right)}^3{log}_{10}^3(f_p/f)$, exhibiting a slight variation from the situation with a narrower peak, $n_1>3/2$. Finally, for the case with $n_1<3/2$, the broken power law form yields SIGWs with an exact power-law form in the infrared regions.

Recently, a common-spectrum process, which appears compatible with the Hellings–Downs angular correlation has been announced by PTA collaborations, including the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) [86,87], Parkers Pulsar Timing Array (PPTA) [88,89], European Pulsar Timing Array (EPTA), Indian Pulsar Timing Array (InPTA) [90,91], and Chinese Pulsar Timing Array (CPTA) [92]. Assuming that the signals come from an ensemble of binary supermassive black hole binary (SMBHB) inspirals and the characteristic-strain spectrum is proportional to $f^{-2/3}$, the strain amplitude of the signals is around $2.4_{-0.6}^{+0.7}\times {10}^{-15}$ at a reference frequency of $1{\mathrm{yr}}^{-1}$ [86]. While supermassive black hole binaries are considered as an explanation [93,94,95,96,97,98,99], alternative interpretations include cosmic inflation [94,95], first-order phase transitions [100,101,102,103,104,105,106,107], cosmic strings [108,109,110,111,112,113,114,115], domain walls [116,117,118], or scalar-induced gravitational waves [4,119,120,121,122,123,124,125,126,127,128,129,130,131]. For more about the physical explanation of signals in the PTA band, please refer to the references [132,133,134,135,136,137,138,139,140,141,142]. The feature of the energy density of SIGWs in the infrared regions aligns with the PTA signals, prompting our investigation using the NANOGrav 15-year dataset and EPTA Data Release 2 (DR2).

The organization of this paper is as follows. In Section 2, we give a brief review of the energy density of SIGWs in the infrared regions. The constraint on the SIGWs in the infrared regions from the PTAs data is presented in Section 3, and conclusions are given in Section 4.

2. The Infrared Behavior of SIGWs

In this section, we will give a brief review of the calculation of the energy density of SIGWs and their universal behavior in the infrared regions with the power spectrum of the primordial curvature being chosen as a broken power law form. The metric with perturbations in the Newtonian gauge and cosmological background is [143]

$$\mathrm{d}{s}^2=-a^2\left(\eta \right)(1+2\mathsf{\Phi })\mathrm{d}{\eta }^2+a^2\left(\eta \right)\left[(1-2\mathsf{\Phi }){\delta }_{ij}+\frac{1}{2}{h}_{ij}\right]\mathrm{d}{x}^{i}d{x}^j,$$

(2)

where $\mathsf{\Phi }$ is the Bardeen potential, a is the scale factor of the Universe, $\eta$ is the conformal time, and $h_{ij}$ is the tensor perturbation. For the second order of tensor perturbations, the linear order scalar perturbations can construct a second-order source to induce the gravitational waves. The power spectrum of these scalar-induced gravitational waves is [1,2,38,40,144,145]

$${\mathcal{P}}_h\left(k,\eta \right)=4{\int }_0^{\infty }\mathrm{d}v{\int }_{\left|1-v\right|}^{1+v}\mathrm{d}u{\left(\frac{4v^2-{\left(1+v^2-u^2\right)}^2}{4vu}\right)}^2\times I_{\mathrm{RD}}^2\left(u,v,x\right){\mathcal{P}}_{\zeta }\left(vk\right){\mathcal{P}}_{\zeta }\left(uk\right),$$

(3)

where $x=k\eta$, and ${\mathcal{P}}_{\zeta }$ is the power spectrum of primordial curvature perturbations generated during inflation. After the horizon re-entry, and at the later time $k\eta \gg 1$, the integrated kernel $I_{\mathrm{RD}}^2$ can be expressed as [38]

$$\begin{array}{cc}\hfill \overline{I_{\mathrm{RD}}^2(v,u,x\to \infty )}& =\frac{1}{2x^2}{\left(\frac{3(u^2+v^2-3)}{4u^3{v}^3}\right)}^2\times \{{\pi }^2{(u^2+v^2-3)}^2\mathsf{\Theta }\left(v+u-\sqrt{3}\right)\hfill \\ & -{\left(4uv-(u^2+v^2-3)log\left|\frac{3-{(u+v)}^2}{3-{(u-v)}^2}\right|\right)}^2\},\hfill \end{array}$$

(4)

where a bar on the top denotes the mean value among several periods. The energy density in gravitational waves is

$${\mathsf{\Omega }}_{\mathrm{GW}}\left(\eta ,k\right)=\frac{1}{24}{\left(\frac{k}{\mathcal{H}\left(\eta \right)}\right)}^2\overline{{\mathcal{P}}_h\left(k,\eta \right)}.$$

(5)

To describe the power spectrum of the primordial curvature perturbations, we use a broken power-law parameterization, which takes the following form

$${\mathcal{P}}_{\zeta }\left(k\right)=\left\{\begin{array}{cc}{A}_{\zeta }{\left(\frac{k}{k_{p_0}}\right)}^{n_1},\hfill & k\le k_{p_0},\hfill \\ A_{\zeta }{\left(\frac{k}{k_{p_0}}\right)}^{n_2},\hfill & k>k_{p_0},\hfill \end{array}\right.$$

(6)

where $k_{p_0}$ is the peak scale and $A_{\zeta }$ is the amplitude at the peak. To ensure a single peak in the power spectrum at the scale $k_{p_0}$, the indices of the parameterization (6) should satisfy the condition $n_1>0$ and $n_2<0$. The magnitude of $n_1$ and $n_2$ determines the nature of the peak: small absolute values result in a broad peak, while large absolute values produce a narrow peak. Consequently, the broken power-law parameterization is versatile and capable of representing both narrow and broad peak power spectra. By substituting the parameterization (6) into the above equations, we can derive the energy density of scalar-induced gravitational waves in the infrared region,

$$h^2{\mathsf{\Omega }}_{\mathrm{GW}}\left(k\right)=\left\{\begin{array}{cc}{A}_{GW}{\left(\frac{f}{f_{\mathrm{year}}}\right)}^{2n_1},\hfill & n_1<3/2\hfill \\ A_{GW}{\left(\frac{f}{f_{\mathrm{year}}}\right)}^3{log}_{10}^3(f_p/f),\hfill & n_1=3/2,\hfill \\ A_{GW}{\left(\frac{f}{f_{\mathrm{year}}}\right)}^3{log}_{10}^2(f_p/f),\hfill & n_1>3/2,\hfill \end{array}\right.$$

(7)

where the relation of the wave number k and the frequency f is

$$f=\frac{k}{2\pi }\simeq 1.6\mathrm{nHz}\left(\frac{k}{{10}^6{\mathrm{Mpc}}^{-1}}\right),$$

(8)

and $f_{\mathrm{year}}$ represents the frequency corresponding to one year. For the scenario of a narrow peak with $n_1>3/2$, the behavior of the energy density of SIGWs in the infrared regions is independent of the model $n_1$ and consistent with the conclusion in Ref. [54,83,84,85]. In contrast, for the broad peak with $n_1<3/2$, the behavior of the energy density of SIGWs in the infrared regions follows an exact power law form.

3. Results

The increasing feature of the energy density of SIGWs in infrared regions aligns with the recent PTA signals. To explain the PTA signal and obtain the constraints on the SIGWs described by Equation (7), we independently perform a Bayesian analysis on the NANOGrav 15-year dataset and the EPTA DR2, respectively. In the Bayesian analysis, we employ the Bilby gravitational wave inference library [146] using nested sampling with the dynesty algorithm [147]. The NANOGrav dataset consists of 14 frequency bins from the NANOGrav 15-year dataset [86,94] and the EPTA dataset consists of 9 bins from the EPTA DR2 [91,95]. To compute the log-likelihood, we first calculate the SIGW energy density at all the frequency bins for each dataset. We then evaluate the log probabilities from independent kernel density estimates of the energy density at each bin. Finally, we obtain the log-likelihood by taking the sum of all these log probability terms for each dataset [148]. Consequently, the likelihood can be expressed as

$$ln{\mathcal{L}}_{\mathcal{D}}\left(\mathsf{\Theta }\right)=\sum _{i=1}^{n_{\mathcal{D}}}ln{\mathcal{L}}_i\left({\mathsf{\Omega }}_{\mathrm{GW}}\left(f_i,\mathsf{\Theta }\right)\right),$$

(9)

where $\mathsf{\Theta }=\{{\mathcal{A}}_{\mathrm{GW}},n_1,f_p\}$ is the collection of parameters presented in Equation (7), $f_i$ is the frequency of each bin, the label $\mathcal{D}=\{\mathrm{nano},\mathrm{epta}\}$ denotes the NANOGrav dataset and EPTA dataset, and $\{n_{\mathrm{nano}}=14$, $n_{\mathrm{epta}}=9\}$. The natural logarithm of the Bayes factors is defined as

$$ln{\mathcal{B}}_{10}\left(\mathcal{D}\right)=ln\left(\frac{p\left(\mathcal{D}\right|{\mathcal{M}}_1)}{p\left(\mathcal{D}\right|{\mathcal{M}}_0)}\right),$$

(10)

where $p\left(\mathcal{D}\right|\mathcal{M})$ is the model evidence. In this work, we choose the fiducial model ${\mathcal{M}}_0$ as the SMBHB model, and the model ${\mathcal{M}}_1$ is given by Equation (7).

In

Table 1. The priors, maximum posterior values, 1-$\sigma$ credible interval bounds of posteriors, and Bayes factor for the model are described by Equation (7). The labels “nano” and “epta” denote the Bayesian analysis conducted using the NANOGrav 15-year dataset and EPTA DR2, respectively. The fiducial model ${\mathcal{M}}_0$ is SMBHB model. The symbol $\mathcal{U}$ denotes the uniform distribution.

Model	Parameters	Prior	Posterior (Nano)	Posterior (Epta)	$ln\mathcal{B}$ (Nano)	$ln\mathcal{B}$ (Epta)
$n_1<3/2$	${\log }_{10}{A}_{\mathrm{GW}}$	$U[-10,-5]$	$-7.{17}_{-0.25}^{+0.24}$	$-7.{06}_{-0.32}^{+0.23}$	$4.4$	$4.8$
	$n_1$	$U[0,1.5]$	$+0.{95}_{-0.17}^{+0.16}$	$+1.{17}_{-0.24}^{+0.19}$
$n_1=3/2$	${\log }_{10}{A}_{\mathrm{GW}}$	$U[-10,-5]$	$-8.{02}_{-0.68}^{+0.42}$	$-9.{27}_{-0.43}^{+0.68}$	$4.5$	$5.3$
	${\log }_{10}(f_p/\mathrm{Hz})$	$U[-9,-3]$	$-6.{67}_{-0.37}^{+0.96}$	$-4.{94}_{-1.25}^{+1.27}$
$n_1>3/2$	${\log }_{10}{A}_{\mathrm{GW}}$	$U[-10,-5]$	$-7.{20}_{-0.53}^{+0.31}$	$-8.{31}_{-0.40}^{+0.54}$	$4.0$	$5.2$
	${\log }_{10}(f_p/\mathrm{Hz})$	$U[-9,-3]$	$-7.{09}_{-0.26}^{+0.80}$	$-5.{20}_{-1.31}^{+1.49}$
SMBHB	${\log }_{10}{\mathrm{A}}_{\mathrm{GW}}$	$U[-10,-5]$	$-8.{14}_{-0.16}^{+0.13}$	$-8.{21}_{-0.08}^{+0.18}$	−	−

, we display the priors of the parameters $\{{\mathcal{A}}_{\mathrm{GW}},n_1,f_p\}$. The priors of the parameters are uniform or log-uniform distributions.

The posteriors for the parameters are shown in Figure 1, Figure 2 and Figure 3, where the left panel and right panel use the NANOGrav 15-year dataset and EPTA DR2, respectively. For the priors listed in Table 1, increasing the interval of the bounds results in almost the same posterior, with the Bayes factor changing slowly as the interval becomes larger. Therefore, priors with larger interval bounds do not significantly affect the posterior and the final conclusions. The corresponding maximum posterior values, 1-$\sigma$ credible interval bounds, and Bayes factors are presented in Table 1, where the labels “nano” and “epta” denote the result by using the NANOGrav dataset and EPTA DR2, respectively. For the broad peak with $n_1<3/2$, the energy density of SIGWs obeys the power law form, and the posterior, displayed in Figure 1, is consistent with the results given in the NANOGrav collaboration [94]. By using the NANOGrav 15-year dataset, the maximum posterior values and 1-$\sigma$ credible interval bounds of posteriors are $n_1=0.{95}_{-0.17}^{+0.16}$ and ${log}_{10}{A}_{\mathrm{GW}}=-7.{17}_{-0.25}^{+0.24}$; by using the EPTA DR2; they are $n_1=1.{17}_{-0.24}^{+0.19}$ and ${log}_{10}{A}_{\mathrm{GW}}=-7.{06}_{-0.32}^{+0.23}$. For the intermediate case with $n_1=3/2$, the posterior of the parameters is illustrated in Figure 2. By using the NANOGrav 15-year dataset, the maximum posterior values and 1-$\sigma$ credible interval bounds of posteriors are ${log}_{10}{f}_p/\mathrm{Hz}=-6.{67}_{-0.37}^{+0.96}$ and ${log}_{10}{A}_{\mathrm{GW}}=-8.{02}_{-0.68}^{+0.42}$, and by using the EPTA DR2, they are ${log}_{10}{f}_p/\mathrm{Hz}=-4.{94}_{-1.25}^{+1.27}$ and ${log}_{10}{A}_{\mathrm{GW}}=-9.{27}_{-0.43}^{+0.68}$. For the narrow peak with $n_1>3/2$, the energy density of SIGWs in infrared regions has a universal log-dependent form, which is independent of the index $n_1$. The posterior for the parameters is shown in Figure 3; by using the NANOGrav 15-year dataset, we have ${log}_{10}{f}_p/\mathrm{Hz}=-7.{09}_{-0.26}^{+0.80}$ and ${log}_{10}{A}_{\mathrm{GW}}=-7.{20}_{-0.53}^{+0.31}$; and by using the EPTA DR2, we have ${log}_{10}{f}_p/\mathrm{Hz}=-5.{20}_{-1.31}^{+1.49}$ and ${log}_{10}{A}_{\mathrm{GW}}=-8.{31}_{-0.40}^{+0.54}$.

Utilizing the maximum posterior values from Table 1, and taking them into Equation (7), we visualize the energy density of the SIGWs in the infrared regions in Figure 4. The left panel and right panel are the fittings of the NANOGrav 15-year dataset and EPTA DR2, respectively. The blue and gray violins represent the observational bins from the NANOGrav 15-year dataset and EPTA DR2 dataset. Meanwhile, the black, green, and red curves correspond to the energy density of SIGWs from broad ($n_1<3/2$), intermediate ($n_1=3/2$), and narrow peaks ($n_1>3/2$), respectively. The three energy density curves are closely packed both for the NANOGrav dataset and EPTA dataset, making it challenging to distinguish between the three models. This result is also supported by the natural logarithm of the Bayes factors listed in Table 1.

4. Discussion and Conclusions

The common-spectrum process, characterized by the Hellings–Downs angular correlation detected by NANOGrav, PPTA, EPTA, and CPTA collaborations can be accounted for by scalar-induced gravitational waves. The energy density of SIGWs arising from the primordial curvature power spectrum with a narrow peak has a universal log-dependent form in the infrared regions. An effective parameterization of the primordial curvature power spectrum is the broken power law form, given by ${\mathcal{P}}_{\zeta }=A_{\zeta }{(k/k_{p_0})}^{n_1}$ for $k<k_{p_0}$ and ${\mathcal{P}}_{\zeta }=A_{\zeta }{(k/k_{p_0})}^{n_2}$ for $k\ge k_{p_0}$. For $n_1>3/2$, this form gives a power spectrum with a narrow peak, and the corresponding SIGWs exhibit the universal form $h^2{\mathsf{\Omega }}_{\mathrm{GW}}=A_{GW}{\left(f/f_{\mathrm{year}}\right)}^3{log}_{10}^2(f_p/f)$; for $n_1=3/2$, this intermediate case induces $h^2{\mathsf{\Omega }}_{\mathrm{GW}}=A_{GW}{\left(f/f_{\mathrm{year}}\right)}^3{log}_{10}^3(f_p/f)$; and for $n_1<3/2$, the primordial curvature power spectrum with a broad peak produces SIGWs with an energy density in the infrared regions exhibiting an exact power-law form.

Utilizing Bayesian analysis on the dataset of NANOGrav 15-year data and EPTA DR2 independently, we determine the posterior distribution of parameters in the SIGW energy density in Equation (7). For the scenario where the primordial curvature power spectrum features a broad peak with $n_1<3/2$, the maximum posterior values and 1-$\sigma$ credible interval bounds of posteriors by using NANOGrav dataset are $n_1=0.{95}_{-0.17}^{+0.16}$ and ${log}_{10}{A}_{\mathrm{GW}}=-7.{17}_{-0.25}^{+0.24}$, and they are $n_1=1.{17}_{-0.24}^{+0.19}$ and ${log}_{10}{A}_{\mathrm{GW}}=-7.{06}_{-0.32}^{+0.23}$ by using EPTA DR2. For the intermediate case with $n_1=3/2$, we have ${log}_{10}{f}_p/\mathrm{Hz}=-6.{67}_{-0.37}^{+0.96}$ and ${log}_{10}{A}_{\mathrm{GW}}=-8.{02}_{-0.68}^{+0.42}$ by using NANOGrav dataset, and ${log}_{10}{f}_p/\mathrm{Hz}=-4.{94}_{-1.25}^{+1.27}$ and ${log}_{10}{A}_{\mathrm{GW}}=-9.{27}_{-0.43}^{+0.68}$ by using EPTA DR2. In the case that the primordial curvature power spectrum has a narrow peak with $n_1>3/2$, we obtain ${log}_{10}{f}_p/\mathrm{Hz}=-7.{09}_{-0.26}^{+0.80}$ and ${log}_{10}{A}_{\mathrm{GW}}=-7.{20}_{-0.53}^{+0.31}$ by using the NANOGrav dataset, and ${log}_{10}{f}_p/\mathrm{Hz}=-5.{20}_{-1.31}^{+1.49}$ and ${log}_{10}{A}_{\mathrm{GW}}=-8.{31}_{-0.40}^{+0.54}$ by using the EPTA DR2. Compared to interpreting the signal as a stochastic background from supermassive black hole binaries, the natural logarithm of the Bayes factors for these three models are around $ln\mathcal{B}\approx 4$ by using the NANOGrav dataset, and $ln\mathcal{B}\approx 5$ by using the EPTA DR2. Consequently, the PTA signal favors the SIGW explanation over SMBHBs, yet cannot distinguish between broad and narrow peak scenarios. Besides SIGWs, other cosmological processes, such as domain walls and phase transitions, are also favored by the PTA dataset in terms of Bayesian analysis. To confirm the source of this PTA signal, we require information from other frequency bands. Future space-based gravitational wave detectors, such as LISA, TianQin, and Taiji, may have the opportunity to provide that information.