Bayesian Implications for the Primordial Black Holes from NANOGrav’s Pulsar-Timing Data Using the Scalar-Induced Gravitational Waves

Abstract

Assuming that the common-spectrum process in the NANOGrav 12.5-year dataset has an origin of scalar-induced gravitational waves, we study the enhancement of primordial curvature perturbations and the mass function of primordial black holes, by performing the Bayesian parameter inference for the first time. We obtain lower limits on the spectral amplitude, i.e., $\mathcal{A}\gtrsim {10}^{-2}$ at 95% confidence level, when assuming the power spectrum of primordial curvature perturbations to follow a log-normal distribution function with width $\sigma$. In the case of $\sigma \to 0$, we find that the primordial black holes with $2\times {10}^{-4}-{10}^{-2}$ solar mass are allowed to compose at least a fraction ${10}^{-6}$ of dark matter. Such a mass range is shifted to more massive regimes for larger values of $\sigma$, e.g., to a regime of $4\times {10}^{-3}-0.2$ solar mass in the case of $\sigma =1$. We expect the planned gravitational-wave experiments to have their best sensitivity to $\mathcal{A}$ in the range of ${10}^{-4}$ to ${10}^{-7}$, depending on the experimental setups. With this level of sensitivity, we can search for primordial black holes throughout the entire parameter space, especially in the mass range of ${10}^{-16}$ to ${10}^{-11}$ solar masses, where they could account for all dark matter. In addition, the importance of multi-band detector networks is emphasized to accomplish our theoretical expectation.

1. Introduction

It is well-known that the primordial black holes (PBHs) could be produced by gravitational collapses of the enhanced primordial curvature perturbations on small scales [1]. The mass function of PBHs is determined by the power spectrum of primordial curvature perturbations [2,3,4]. To produce a sizeable quantity of PBHs, the scalar spectral amplitude should conquer a critical magnitude ∼$0.1$ [5], which is about eight orders of magnitudes larger than the spectral amplitude on large scales. On the other hand, the dynamics of cosmic inflation [6,7,8] is one of the greatest puzzles of cosmology to be uncovered. In particular, the small-scale primordial curvature perturbations were generated by quantum fluctuations at the late-time stage of inflation. Their power spectrum encodes characteristic information about the late-time dynamics of inflation, which could be different from that of the early-time one measured precisely with experiments on cosmic microwave background (CMB) [9,10] and large-scale structure (LSS) [11,12]. Therefore, we can gain more information about cosmic inflation by investigating the PBHs.

The scenario of PBHs can also be a reasonable candidate of dark matter, or compose a part of the latter [13,14]. The mass function of PBHs is defined as a ratio between the energy density fractions of PBHs and dark matter. Extensive studies have shown observational constraints on the mass function of PBHs (e.g., see reviews in Refs. [15,16]). In addition, the scenario of PBHs was suggested to interpret the formation mechanism of supermassive black holes (SMBHs) [17] and the events of compact binary coalescences [18,19,20,21,22,23,24,25], which emitted the gravitational waves (GWs) observed by the Advanced LIGO, Virgo and KAGRA detectors [26,27,28].

Accompanied by the formation process of PBHs, the scalar induced gravitational waves (SIGWs) [29,30,31,32] were inevitably produced by the enhanced primordial curvature perturbations. Specifically, they were produced through nonlinear couplings of tensor-scalar-scalar modes in the very early cosmos. Given the power spectrum of primordial curvature perturbations, a semi-analytic formula was provided to compute the energy density fraction spectrum of SIGWs [33,34]. It has been widely used to investigate the mass function of PBHs in the literature.

In this work, we will investigate the mass function of PBHs, as well as the power spectrum of primordial curvature perturbations, by analyzing a 12.5-year dataset from North American Nanohertz Observatory for Gravitational Waves (NANOGrav)1 [35]. Strong evidence of a stochastic common-spectrum process was reported by NANOGrav Collaboration. If this signal has an origin of SIGWs, we can straightforwardly obtain constraints on the power spectrum of primordial curvature perturbations. Since the latter can collapse to form the PBHs due to gravity, we therefore obtain corresponding constraints on the mass function of PBHs. In fact, the scenario of PBHs with different masses has been suggested to interpret such a signal (see, e.g., Refs. [36,37,38,39,40,41]). However, there is still lack of systematic Bayesian parameter inferences, which will be performed for the first time in this paper.

The rest of the paper is arranged as follows. We briefly review the theory of SIGWs in Section 2. By performing Bayesian analysis, we obtain the NANOGrav constraints on the enhancement of primordial curvature perturbations in Section 3. Correspondingly, we demonstrate the physical implications for the mass function of PBHs in Section 4. We forecast the sensitivity curves of ongoing and planned gravitational-wave detectors in Section 5. The conclusion and discussion are shown in Section 6.

2. Theory of Scalar-Induced Gravitational Waves

As previously shown [33,34], there is a lengthy and tedious derivation process to obtain a semi-analytic formula for the energy density fraction spectrum of SIGWs in the early cosmos. Therefore, we summarize some theoretical results that would be essential to our data analysis. The energy density fraction spectrum of GWs is defined as ${\mathsf{\Omega }}_{\mathrm{GW}}\left(k\right)={\rho }_{\mathrm{GW}}\left(k\right)/{\rho }_c$, where k denotes the GW wavenumber, ${\rho }_c$ is the critical energy density of the cosmos, and $\int {\rho }_{\mathrm{GW}}\left(k\right)dlnk$ denotes the total energy density of GWs. The SIGWs in the subhorizon have the energy density spectrum as ${\rho }_{\mathrm{GW}}\sim \langle \overline{h_{ij,l}{h}_{ij,l}}\rangle$ [42], where the overbar stands for the oscillation average and the angle brackets defines the power spectrum. The strain of SIGWs, as the second-order tensor perturbations produced by the linear scalar perturbations, can be roughly expressed as $h\sim {\varphi }^2$, where $\varphi$ denotes the linear scalar perturbation. If considering the Gaussian scalar perturbations, we obtain ${\rho }_{\mathrm{GW}}\sim \langle {\varphi }^4\rangle \sim \langle {\varphi }^2\rangle \langle {\varphi }^2\rangle$, where $\langle {\varphi }^2\rangle$ is determined by the initial value and transfer function of $\varphi$. Therefore, following Equation (14) of Ref. [33], the energy density fraction spectrum of SIGWs in the radiation-dominated cosmos is given as

$$\begin{array}{c}\hfill {\mathsf{\Omega }}_{\mathrm{GW}}\left(k\right)={\int }_0^{\infty }dv{\int }_{|1-v|}^{|1+v|}duF(v,u){\mathcal{P}}_{\mathcal{R}}\left(vk\right){\mathcal{P}}_{\mathcal{R}}\left(uk\right),\end{array}$$

(1)

where ${\mathcal{P}}_{\mathcal{R}}\left(uk\right)$ (and ${\mathcal{P}}_{\mathcal{R}}\left(vk\right)$) is power spectrum of primordial curvature perturbations at wavenumber $uk$ (and $vk$) with u (and v) being a dimensionless variable, and $F(v,u)$ is a functional of the scalar transfer function describing the post-inflationary evolution of $\varphi$ and can be simplified as

$$\begin{array}{ccc}\hfill F(v,u)& =& \frac{3}{2^{10}{v}^8{u}^8}{\left[4v^2-{(v^2-u^2+1)}^2\right]}^2{\left(v^2+u^2-3\right)}^2\hfill \\ & & \times \{{\left[\left(v^2+u^2-3\right)ln\left(\left|\frac{3-{(v+u)}^2}{3-{(v-u)}^2}\right|\right)-4vu\right]}^2\hfill \\ & & +{\pi }^2{\left(v^2+u^2-3\right)}^2\mathsf{\Theta }\left(v+u-\sqrt{3}\right)\}.\hfill \end{array}$$

(2)

Due to a lack of direct measurements, the nature of primordial curvature perturbations on small scales remains unknown, compared with that on the largest scales measured by the cosmic microwave background. For computational simplicity, we take two concrete expressions of ${\mathcal{P}}_{\mathcal{R}}\left(k\right)$, which are frequently adopted in the literature (e.g., Refs. [33,43,44,45,46]), as follows

$$\begin{array}{ccc}& & {\mathcal{P}}_{\mathcal{R}}\left(k\right)=\mathcal{A}\delta \left(ln\frac{k}{k_{*}}\right),\hfill \end{array}$$

(3)

$$\begin{array}{ccc}& & {\mathcal{P}}_{\mathcal{R}}\left(k\right)=\frac{\mathcal{A}}{\sqrt{2\pi }\sigma }exp\left(-\frac{{ln}^2\left(k/k_{*}\right)}{2{\sigma }^2}\right),\hfill \end{array}$$

(4)

where $\mathcal{A}$ is an amplitude, $\sigma$ is a standard variance, and $k_{*}$ is a pivot wavenumber. We can relate the wavenumber with the frequency, i.e., $k=2\pi f$ and $k_{*}=2\pi f_{*}$. During the process of parameter inferences, both $\mathcal{A}$ and $f_{*}$ are independent parameters. However, we will not let $\sigma$ be independent, even though different values of $\sigma$ would alter our results. In fact, it is challenging to explore the full parameter space of $\sigma$ based on the current dataset. Therefore, we let $\sigma =1$ for simplicity in this work. We can test the robustness of our conclusions by comparing the results obtained when taking two different values of $\sigma$. In particular, we will compare the results from Equations (3) and (4). We leave this point to be demonstrated further in the following sections.

In Figure 1, we show the energy density fraction spectrum of SIGWs, which is normalized with ${\mathcal{A}}^2$. For comparison, we show the results when taking different versions of $\sigma$. Equation (3) can be viewed as a limit of Equation (4) when $\sigma \to 0$. In the regime of $f\ll f_{*}$, the choice of $\sigma \to 0$ displays the highest spectral amplitude than non-vanishing choices of $\sigma$. This implies that in such a case, we can obtain the best sensitivity on measurements of $\mathcal{A}$ for a given GW experiment. Correspondingly, we can obtain the strictest constraints on the mass function of PBHs. We will demonstrate these points in the following sections.

Once the early-universe quantity ${\mathsf{\Omega }}_{\mathrm{GW}}\left(k\right)$ is given in Equation (1), we have the energy density fraction spectrum of SIGWs in the present cosmos to be [43]

$${\mathsf{\Omega }}_{\mathrm{GW},0}\left(f\right)={\mathsf{\Omega }}_{\mathrm{r},0}\left[\frac{g_{*,\rho }\left(T\right)}{g_{*,\rho }\left(T_{\mathrm{eq}}\right)}\right]{\left[\frac{g_{*,s}\left(T_{\mathrm{eq}}\right)}{g_{*,s}\left(T\right)}\right]}^{\frac{4}{3}}{\mathsf{\Omega }}_{\mathrm{GW}}\left(k\right),$$

(5)

where ${\mathsf{\Omega }}_{\mathrm{r},0}$ is the present energy density fraction of radiations, and the subscript ${}_{\mathrm{eq}}$ stands for cosmological quantities at the epoch of matter-radiation equality. Here, we have taken into account contributions from the effective relativistic degrees of the cosmos, i.e., $g_{*,\rho }$ and $g_{*,s}$, as functions of f by interpolating their tabulated data in terms of cosmic temperature T in Ref. [47]2 as well as by considering a relation between T and f in Wang et al. [43], i.e.,

$$\frac{f}{\mathrm{nHz}}=26.5\left(\frac{T}{\mathrm{GeV}}\right){\left[\frac{g_{*,\rho }\left(T\right)}{106.75}\right]}^{\frac{1}{2}}{\left[\frac{g_{*,s}\left(T\right)}{106.75}\right]}^{-\frac{1}{3}}.$$

(6)

Throughout this paper, we use the measured value of cosmological parameters in the data release 2018 of Planck satellite [10]. The publicly available Astropy [48,49,50] software3 is adopted to evaluate all cosmological quantities. Please note that Equation (5) is one of the leading formulas that will be used during Bayesian analysis in the next section.

3. NANOGrav Constraints on Primordial Curvature Perturbations

When a pulsar timing array (PTA) experiment measures the stochastic gravitaitonal wave background, the timing residual cross-power spectral density is given by $S_{ab}\left(f\right)={\mathsf{\Gamma }}_{ab}{h}_c^2/\left(12{\pi }^2{f}^3\right)$, which can be shown by combining Equations (1) and (2) in Arzoumanian et al. [51]. Here, ${\mathsf{\Gamma }}_{ab}$ stands for the overlap reduction function (ORF) that describes the correlation between two pulsars a and b as a function of line-of-sight separation angle between them. Following Equation (17) in Maggiore [52], $h_c\left(f\right)$ is the characteristic strain defined as $h_c^2=3H_0^2{\mathsf{\Omega }}_{\mathrm{GW},0}/\left(2{\pi }^2{f}^2\right)$, where $H_0=100h_0\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ is the Hubble constant with $h_0$ being the reduced Hubble constant. Therefore, the timing residual cross-power spectral density becomes [51]

$$S_{ab}={\mathsf{\Gamma }}_{ab}\frac{H_0^2{f}_{\mathrm{yr}}^{-5}}{8{\pi }^4}{\left(\frac{f}{f_{\mathrm{yr}}}\right)}^{-5}{\mathsf{\Omega }}_{\mathrm{GW},0}\left(f\right),$$

(7)

where $f_{\mathrm{yr}}$ is a pivot frequency corresponding to a duration time of 1 year, and the formula for ${\mathsf{\Omega }}_{\mathrm{GW},0}\left(f\right)$ is shown in Equation (5). For an isotropic background of GW, e.g., the SIGWs considered in this work, we take ORF to be the Hellings and Downs coefficients [53].

The timing residual data for a single pulsar is decomposed in its individual constituents, i.e., [54]

$$\delta \mathit{t}=M\mathit{ϵ}+F\mathit{a}+U\mathit{j}+\mathit{n}.$$

(8)

The term $M\mathit{ϵ}$ stands for the inaccuracies in the subtraction of timing model, where M denotes the timing model design matrix, and $\mathit{ϵ}$ is a vector describing small offsets for the timing model parameters. The matrix M is computed through libstempo [55]4, which is a python interface for TEMPO2 [56,57]5 timing package. We use the latest Jet Propulsion Laboratory (JPL) solar system ephemeris (SSE) model, DE438 [58], in the timing model fits. The term $F\mathit{a}$ accounts for all low-frequency signals, including the pulsar-intrinsic red noise. The Fourier design matrix F has alternating sine and cosine functions, and $\mathit{a}$ is a vector comprised of Fourier coefficients at the integer ($1,2,\dots ,N_{\mathrm{mode}}$) multiples of the harmonic base frequency $1/T_s$, where $T_s$ denotes the span between the minimum and maximum time of arrivals (TOAs) in the array [59]. Described by a per-epoch variance (ECORR) for each receiver and backend system [54], the term $U\mathit{j}$ denotes the white noise which is fully correlated for simultaneous observations at different frequencies, but fully uncorrelated in time. The matrix U maps all TOAs observed simultaneously at different frequencies to a total TOA, and $\mathit{j}$ is the per-epoch white noise which is fully correlated across all observing frequencies. The term $\mathit{n}$ is the timing residual introduced by Gaussian white noise, which is described by the parameters of the TOA uncertainties (EFAC) and an additive variance (EQUAD) for each receiver and backend system [54].

To estimate the allowed parameter space, we perform Bayesian parameter inferences by analyzing a dataset of 45 pulsars in the 12.5-year data release of NANOGrav Collaboration [35]6. We list all independent parameters as well as their priors in

Table 1. Priors used in all analyses performed in this paper.

Parameter	Description	Prior	Comment
White Noise
$E_k$	EFAC per backend/receiver system	Uniform $[0,10]$	single-pulsar analysis only
$Q_k$ [s]	EQUAD per backend/receiver system	log-Uniform $[-8.5,-5]$	single-pulsar analysis only
$J_k$ [s]	ECORR per backend/receiver system	log-Uniform $[-8.5,-5]$	single-pulsar analysis only
Red Noise
$A_{\mathrm{red}}$	power-law spectral amplitude	log-Uniform $[-20,-11]$	one parameter per pulsar
${\gamma }_{\mathrm{red}}$	power-law spectral index	Uniform $[0,7]$	one parameter per pulsar
Primordial curvature perturbations
$log\mathcal{A}$	spectral amplitude	Uniform $[-3,0]$	one parameter for PTA
$log{f}_{*}$ [Hz]	spectral pivot frequency	Uniform $[-10,0]$	one parameter for PTA

. We adopt $N_{\mathrm{mode}}=30$ frequency bands to the power-law spectrum of pulsar-intrinsic red noise and the common-spectrum process. In practice, we will use the publicly available enterprise [60]7 to compute the likelihoods and PTMCMCSampler [61]8 to perform Markov-Chain Monte-Carlo sampling.

Our results are shown in Figure 2, which displays the 95% CL contour plots of $log\mathcal{A}$ and $log{f}_{*}$ for the power spectra of primordial curvature perturbations with $\sigma \to 0$ (red solid curve) and $\sigma =1$ (blue solid curve), respectively. We find that $\mathcal{A}$ has lower limits around ${10}^{-2}$ for the two different choices of $\sigma$. In contrast, the NANOGrav dataset prefers different frequency bands. For a given value of $\mathcal{A}$, the frequency range in the case of $\sigma \to 0$ is almost always larger than that in the case of $\sigma =1$ by a few times. This result is consistent with the expectations of Figure 1. However, we find strong positive correlations between $log\mathcal{A}$ and $log{f}_{*}$ in either choice of $\sigma$. Furthermore, in Figure 2, we also label the spectral amplitudes that produce the PBHs with total abundance $\overline{f}{}_{\mathrm{pbh}}=1$ (dashed curves) and $\overline{f}{}_{\mathrm{pbh}}={10}^{-10}$ (dotted curves) that will be interpreted in the next section.

4. Implications for Primordial Black Holes

We define the mass function of PBHs as $f_{\mathrm{pbh}}\left(m\right)={\mathsf{\Omega }}_{\mathrm{dm}}^{-1}d{\mathsf{\Omega }}_{\mathrm{pbh}}/dlnm$, in which ${\mathsf{\Omega }}_{\mathrm{dm}}$ and ${\mathsf{\Omega }}_{\mathrm{pbh}}$ stand for the present energy density fractions of dark matter and PBHs with mass m, respectively. We adopt the concept of critical collapse [62,63] and Press-Schechter formalism [64]. Therefore, following Equation (9) in Wang et al. [43], we have

$$\begin{array}{c}\hfill f_{\mathrm{pbh}}\left(m\right)=\frac{{\mathsf{\Omega }}_{\mathrm{m}}}{{\mathsf{\Omega }}_{\mathrm{dm}}}\int g\left(T\right)\tilde{\beta }(m,m_H)dln{m}_H,\end{array}$$

(9)

where ${\mathsf{\Omega }}_{\mathrm{m}}$ is the present energy density fraction of non-relativistic matter. Here, we denote

$$\begin{array}{ccc}\hfill g\left(T\right)& =& \frac{g_{*,\rho }\left(T\right)}{g_{*,\rho }\left(T_{\mathrm{eq}}\right)}\frac{g_{*,s}\left(T_{\mathrm{eq}}\right)}{g_{*,s}\left(T\right)}\frac{T}{T_{\mathrm{eq}}},\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill \tilde{\beta }(m,m_H)& =& \frac{\kappa {\mu }^{\gamma +1}}{\sqrt{2\pi }\gamma \mathsf{\Delta }\left(k\right)}\exp \left[-\frac{{\left({\delta }_c+\mu \right)}^2}{2{\mathsf{\Delta }}^2\left(k\right)}\right],\hfill \end{array}$$

(11)

where we have $\mu ={[m/\left(\kappa m_H\right)]}^{1/\gamma }$ with numerical constants $\kappa =3.3$ and $\gamma =0.36$, and the critical overdensity for gravitational collapse is ${\delta }_c=0.45$. Here, we disregard corrections to ${\delta }_c$ from, e.g., QCD equation of state [65]. The quantity $\tilde{\beta }(m,m_H)$ is a mass distribution function of PBHs, which were produced when the horizon mass was $m_H$. The subscript ${}_{\mathrm{eq}}$ stands for cosmological quantities at the epoch of matter-radiation equality. Cosmic temperature T is related to $m_H$ enclosed by the Hubble horizon, i.e., [43]

$$\frac{m_H}{M_{\odot }}=4.76\times {10}^{-2}{\left(\frac{T}{\mathrm{GeV}}\right)}^{-2}{\left[\frac{g_{*,\rho }\left(T\right)}{106.75}\right]}^{-\frac{1}{2}},$$

(12)

where $M_{\odot }$ denotes the solar mass. We can relate f with $m_H$ by combining Equation (6) with Equation (12) and by reducing T from both equations. The coarse-grained fluctuations in the radiation-dominated cosmos are given by [66,67]

$${\mathsf{\Delta }}^2\left(k\right)=\frac{16}{81}\int dlnq{\left[w\left(\frac{q}{k}\right)\right]}^2{\left(\frac{q}{k}\right)}^4{\mathcal{T}}^2\left(q,\frac{1}{k}\right){\mathcal{P}}_{\mathcal{R}}\left(q\right),$$

(13)

where $w\left(y\right)=\exp (-y^2/2)$ is a Gaussian window function and $\mathcal{T}(q,\tau )=3(sinx-xcosx)/x^3$ with $x=q\tau /\sqrt{3}$ is a scalar transfer function. We further define the total abundance of PBHs in dark matter to be $\overline{f}{}_{\mathrm{pbh}}=\int f_{\mathrm{pbh}}\left(m\right)dlnm$, and define the average mass of PBHs to be $\overline{m}=\int f_{\mathrm{pbh}}\left(m\right)dm/\overline{f}{}_{\mathrm{pbh}}$. The latter roughly displays $\overline{m}\simeq m_{*}$, where $m_{*}$ is corresponded to $k_{*}$.

Once the constraints on the power spectrum of primordial curvature perturbations are obtained, as shown in Figure 2, we can recast them into constraints on the mass function of PBHs, or more precisely, on the average mass and total abundance of PBHs. We show the results in Figure 3. The shaded regions are allowed by the NANOGrav 12.5-year dataset for $\sigma \to 0$ (red region) and $\sigma =1$ (blue region). When at least one fraction (e.g., $f_{\mathrm{pbh}}={10}^{-10}$) of dark matter is composed of PBHs, the mass range $2\times {10}^{-4}-{10}^{-2}{M}_{\odot }$ ($4\times {10}^{-3}-0.2M_{\odot }$) for $\sigma \to 0$ ($\sigma =1$) is preferred by the current dataset. Based on Figure 3 in Wang et al. [43], we find that these mass ranges can be cross-checked with high-precision by observing the GWs emitted from inspiraling stage of PBH binaries. In addition, they might be further tested by measuring the anisotropies in stochastic gravitational-wave background (SGWB) [68]. For comparison, we also depict the existing upper limit (cyan curve) on the mass function of PBHs, as reviewed in Carr et al. [15]. We find the SIGW probe to be more powerful than electromagnetic probes, implying that a larger parameter space can be explored with the SIGW probe. This can also be understood by revisiting Figure 2, in which we depicted the curves labeling $\overline{f}{}_{\mathrm{pbh}}={10}^{-10}$. The latter is corresponded to $\mathcal{A}\simeq \mathrm{few}\times {10}^{-2}$. In contrast, the NANOGrav experiment has reached $\mathcal{A}\simeq {10}^{-2}$ in the most sensitive frequency band ∼$({10}^{-9}-{10}^{-8})$ Hz.

5. Constraints from Ongoing and Planned Gravitational-Wave Detectors

In the future, the power spectrum of primordial curvature perturbations and the mass function of PBHs, which remain unknown until now, can be further explored with ongoing and planned GW experiments, such as Square Kilometer Array (SKA) [69,70,71], $\mu$Ares [72], Laser Interferometer Space Antenna (LISA) [73,74], Big Bang Observer (BBO) [75,76], Decihertz Interferometer Gravitational wave Observatory (DECIGO) [77,78], Einstein Telescope (ET) [79] and Advanced LIGO and Virgo [26,27,28]. Such multi-band observations could explore all possible parameter spaces of primordial curvature perturbations and PBHs. Other experiments are not considered in this work, but our method can be generalized straightforwardly to study them, if needed.

To estimate the sensitivity curve of a given GW experiment consisting of $n_{\det }$ detectors, we define an optimal signal-to-noise ratio (SNR) denoted with $\rho$ as follows [80]

$${\rho }^2=n_{\det }{T}_{\mathrm{obs}}{\int }_{f_{\min }}^{f_{\max }}{\left[\frac{S_{\mathrm{GW}}\left(f\right)}{S_n^{\mathrm{eff}}\left(f\right)}\right]}^{2}df,$$

(14)

where $T_{\mathrm{obs}}$ is a duration time of observing run, the power spectral density (PSD) of SIGWs is defined as

$$S_{\mathrm{GW}}\left(f\right)=\left(\frac{3H_0^2}{2{\pi }^2{f}^3}\right){\mathsf{\Omega }}_{\mathrm{GW},0}\left(f\right),$$

(15)

and the effective noise PSD of the detector network is denoted as $S_n^{\mathrm{eff}}\left(f\right)$ that is a function of f. The concrete setups of aforementioned experiments are summarized in Table 2 of Ref. [81] and references therein. We consider a single detector for LISA, two independent detectors for $\mu$Ares, BBO and DECIGO, three detectors for ET, and 200 pulsars for SKA. For comparison, we consider one detector for Advanced LIGO with an observing duration of four years and 100% duty circle.

Given a desired value of SNR, which is unity in this work, we obtain the minimal detectable ${\mathcal{A}}_{\min }$ for the given experiment by resolving Equation (14). Since ${\mathsf{\Omega }}_{\mathrm{GW},0}\left(f\right)$ is uniquely determined by $\mathcal{A}$ and $f_{*}$, we depict the theoretical expectation of ${\mathcal{A}}_{\min }$ in the $\mathcal{A}-f_{*}$ plane in Figure 4 and Figure 5 for the choices of $\sigma \to 0$ and $\sigma =1$, respectively. For comparison, we plot the exclusion region on $\mathcal{A}$ from the Advanced LIGO and Virgo first three observing runs (red shaded) [82].

Based on Figure 4 and Figure 5, we find that the regions allowed by the current dataset of NANOGrav (enclosed by black curves) can be thoroughly tested with SKA and $\mu$Ares. Smaller values of $\mathcal{A}$, i.e., ∼$({10}^{-4}-{10}^{-7})$, are also expected to be reached by the planned GW experiments. The sensitivity of Advanced LIGO to measure the SIGWs is expected to be improved by one order of magnitude in upcoming observing runs. If other values of SNR are desired, we can obtain revised ${\mathcal{A}}_{\min }$ by rescaling the above results by multiplying a factor of ${\left(\mathrm{SNR}\right)}^{1/2}$. After recasting the expected constraints on $\mathcal{A}$ into constraints on the mass function of PBHs, we further find that the full parameter space of PBHs, that account for at least a fraction of dark matter, e.g., $f_{\mathrm{pbh}}\sim {10}^{-10}$, can be thoroughly explored with multi-band GW measurements, e.g., a detector network composed of SKA, $\mu$Ares and ET, or other detector networks. Therefore, we expect that the scenario of PBHs as a candidate of dark matter can be supported or vetoed by future GW observations.

6. Conclusions and Discussion

In this work, we obtained new constraints on the power spectrum of primordial curvature perturbations and the mass function of PBHs by searching for the energy density fraction spectrum of SIGWs in the NANOGrav 12.5-year dataset. We found the lower limits on $\mathcal{A}$, namely, $\mathcal{A}\gtrsim {10}^{-2}$ and showed the 95% CL contours of $\mathcal{A}$ and $f_{*}$ (see Figure 2). Recasting these contours into the PBH mass-abundance plane, we showed the parameter space of PBHs allowed by the NANOGrav (see Figure 3). We found that at least a fraction of dark matter can be interpreted with the scenario of PBHs in the mass range ∼$({10}^{-4}-{10}^{-1})M_{\odot }$. We also studied dependence of the above results on the value of $\sigma$. Furthermore, this mass range is expected to be cross-checked with high-precision by observing the GWs from PBH mergers (see Figure 3 in Wang et al. [43]). We also found that the SIGW probe is much more powerful than other probes to search for the PBHs.

We further forecasted the sensitivity curves of ongoing and planned GW experiments on detection of PBHs and the primordial curvature perturbations by searching for SIGWs. We found that the primordial curvature perturbations with spectral amplitude larger than ∼$({10}^{-4}-{10}^{-7})$ can be measured with planned GW detection programs. The sensitivity of Advanced LIGO-Virgo to detect the SIGWs was expected to be improved by one order of magnitude in the near future. This prediction may promote extensive investigations of cosmic inflation at late-time stages. Meanwhile, the scenario of PBHs within almost whole mass range can be thoroughly tested, since the critical spectral amplitude to form the PBHs is much larger than ${10}^{-4}$ (see Figure 4 and Figure 5). In particular, we can search for the PBHs within mass range $({10}^{-16}-{10}^{-11})M_{\odot }$, which can compose all of the dark matter and otherwise are beyond capabilities of other probes. In addition, we emphasized the importance of multi-band GW detector networks for accomplishing the above theoretical expectations.

In this paper, we made several assumptions to simplify our computations. First, we disregarded possible contributions of primordial non-Gaussianity to the formation of PBHs [46,83,84,85,86] and to the production of SIGWs [45]. It is an interesting topic to study the primordial non-Gaussianity, deserving an independent work. Second, we took into account the median value of the effective relativistic degrees of freedom of the early Universe, but disregarded their uncertainties [47]. In fact, changing the above two assumptions could alter the mass function of PBHs and the energy density fraction spectrum of SIGWs. We would leave a possible study of this question to future works. Third, we disregarded contributions of nonlinear cosmological perturbations to the energy density fraction spectrum of SIGWs, since there is not a complete theory of SIGWs at the third order [87,88,89,90,91]. We might revisit this assumption in future, once the theory is complete.