Probing Primordial-Black-Hole Dark Matter with Scalar Induced Gravitational Waves

TL;DR

This paper investigates how third-order corrections in scalar-induced gravitational waves can enhance detection prospects and constrain primordial black hole dark matter models, especially with upcoming LISA and PTA observations.

Contribution

It calculates third-order corrections to scalar-induced GWs, showing increased amplitude and extended frequency cutoff, impacting PBH dark matter constraints and detection strategies.

Findings

Third-order corrections increase GW signal amplitude by about 10%.

Extended cutoff frequency from 2k* to 3k* broadens detectable PBH mass range.

Null GW detections can exclude PBHs as all dark matter and link GW events to PBHs.

Abstract

The possibility that primordial black holes (PBHs) represent all of the dark matter (DM) in the Universe and explain the coalescences of binary black holes detected by LIGO/Virgo has attracted a lot of attention. PBHs are generated by the enhancement of scalar perturbations which inevitably produce the induced gravitational waves (GWs). We calculate the induced GWs up to the third-order correction which not only enhances the amplitude of induced GWs, but also extends the cutoff frequency from $2k_*$ to $3k_*$. Such effects of the third-order correction lead to an around $10\%$ increase of the signal-to-noise ratio (SNR) for both LISA and pulsar timing array (PTA) observations, and significantly widen the mass range of PBHs in the stellar mass window accompanying detectable induced GWs for PTA observations including IPTA, FAST and SKA. On the other hand, the null detections of the induced GWs by LISA and PTA experiments will exclude the possibility that all of the DM is comprised of PBHs and the GW events detected by LIGO/Virgo are generated by PBHs.

Figures

Equations

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

$${\varphi }_{\mathbf{k}}(\eta )\equiv {\varphi }_{\mathbf{k}}T(k\eta ),$$

$$T(k\eta )=\frac{9}{(k\eta {)}^2}\left[\frac{\sin (k\eta /\sqrt{3})}{k\eta /\sqrt{3}}-\cos (k\eta /\sqrt{3})\right],$$

$$h_{ij}^{\prime \prime }+2\mathcal{H}{h}_{ij}^{\prime }-{\nabla }^2{h}_{ij}=-4{\mathcal{T}}_{ij}^{\ell m}{S}_{\ell m},$$

$$S_{ij}^{(2)}=4\varphi {\partial }_i{\partial }_j\varphi +2{\partial }_i\varphi {\partial }_j\varphi -{\partial }_i\left(\varphi +\frac{{\varphi }^{\prime }}{\mathcal{H}}\right){\partial }_j\left(\varphi +\frac{{\varphi }^{\prime }}{\mathcal{H}}\right),$$

$$\begin{array}{rl}{S}_{ij}^{(3)}=& \frac{1}{\mathcal{H}}\left(12\mathcal{H}\varphi -{\varphi }^{\prime }\right){\partial }_i\varphi {\partial }_j\varphi -\frac{1}{{\mathcal{H}}^3}\left(4\mathcal{H}\varphi -{\varphi }^{\prime }\right){\partial }_i{\varphi }^{\prime }{\partial }_j{\varphi }^{\prime }\\ & +\frac{1}{3{\mathcal{H}}^4}\left(2{\partial }^2\varphi -9\mathcal{H}{\varphi }^{\prime }\right){\partial }_i\left(\mathcal{H}\varphi +{\varphi }^{\prime }\right){\partial }_j\left(\mathcal{H}\varphi +{\varphi }^{\prime }\right),\end{array}$$

$$\begin{split}S_{ij}^{(4)}=&16\phi^{3}\partial_{i}\partial_{j}\phi+\frac{1}{3\mathcal{H}^{3}}\Big{[}2\phi^{\prime}\partial^{2}\phi-9\mathcal{H}\phi^{\prime 2}-8\mathcal{H}\phi\partial^{2}\phi\\ &\qquad\qquad\qquad\quad+18\mathcal{H}^{2}\phi\phi^{\prime}+96\mathcal{H}^{3}\phi^{2}\Big{]}\partial_{i}\phi\partial_{j}\phi\\ &+\frac{2}{3\mathcal{H}^{5}}\Big{[}-\phi^{\prime}\partial^{2}\phi+3\mathcal{H}\phi^{\prime 2}+4\mathcal{H}\phi\partial^{2}\phi\\ &\qquad\qquad\qquad\quad\quad+3\mathcal{H}^{2}\phi\phi^{\prime}-12\mathcal{H}^{3}\phi^{2}\Big{]}\partial_{i}\phi^{\prime}\partial_{j}\phi^{\prime}\\ &+\frac{1}{36\mathcal{H}^{6}}\Big{[}-16(\partial^{2}\phi)^{2}-3\partial_{k}\phi^{\prime}\partial^{k}\phi^{\prime}+120\mathcal{H}\phi^{\prime}\partial^{2}\phi\\ &\qquad\quad-6\mathcal{H}\partial_{k}\phi\partial^{k}\phi^{\prime}+144\mathcal{H}^{2}\phi\partial^{2}\phi-180\mathcal{H}^{2}\phi^{\prime 2}\\ &\qquad\qquad+33\mathcal{H}^{2}\partial_{k}\phi\partial^{k}\phi-504\mathcal{H}^{3}\phi\phi^{\prime}-144\mathcal{H}^{4}\phi^{2}\Big{]}\\ &\qquad\qquad\times\partial_{i}\left(\mathcal{H}\phi+\phi^{\prime}\right)\partial_{j}\left(\mathcal{H}\phi+\phi^{\prime}\right).\end{split}$$

$$⟨h(\mathbf{k},\eta )h({\mathbf{k}}^{\prime },\eta )⟩\equiv \frac{2{\pi }^2}{k^3}{\mathcal{P}}_h(k,\eta )\delta (\mathbf{k}+{\mathbf{k}}^{\prime }).$$

$${\rho }_{\mathrm{GW}}=\frac{M_p^2}{16a^2}⟨\overline{{\partial }_k{h}_{ij}{\partial }^k{h}^{ij}}⟩,$$

$${\Omega }_{\mathrm{GW}}(\eta ,k)\equiv \frac{{\rho }_{\mathrm{GW}}(k,\eta )}{{\rho }_{\mathrm{cr}}}=\frac{1}{12}{\left(\frac{k}{\mathcal{H}}\right)}^2\overline{{\mathcal{P}}_h(k,\eta )},$$

$$h(\mathbf{k},\eta )=\frac{1}{ka(\eta )}\int \mathrm{d}\tilde{\eta }\sin (k\eta -k\tilde{\eta })a(\tilde{\eta }){\mathcal{S}}_{\mathbf{k}}(\tilde{\eta }),$$

$$⟨{\zeta }_{\mathbf{k}}{\zeta }_{{\mathbf{k}}^{\prime }}⟩\equiv \frac{2{\pi }^2}{k^3}{\mathcal{P}}_{\zeta }(k)\delta (\mathbf{k}+{\mathbf{k}}^{\prime }).$$

$${\mathcal{P}}_{\zeta }(k)=A{k}_{\ast }\delta \left(k-k_{\ast }\right),$$

$$\beta ={\int }_{{\zeta }_c}^{+\infty }\frac{\mathrm{d}\zeta }{\sqrt{2\pi }\sigma }{e}^{-{\zeta }^2/2{\sigma }^2}=\frac{1}{2}\mathrm{erfc}\left(\frac{{\zeta }_c}{\sqrt{2A}}\right),$$

$$f_{\mathrm{pbh}}\simeq 2.5\times 10^8\beta {\left(\frac{g_{\ast }^{\mathrm{form}}}{10.75}\right)}^{-\frac{1}{4}}{\left(\frac{m_{\mathrm{pbh}}}{M_{\odot }}\right)}^{-\frac{1}{2}},$$

$$\frac{m_{\mathrm{pbh}}}{M_{\odot }}\approx 2.3\times 10^{18}{\left(\frac{3.91}{g_{\ast }^{\mathrm{form}}}\right)}^{1/6}{\left(\frac{H_0}{f_{\ast }}\right)}^2,$$

$${\mathcal{S}}_{\mathbf{k}}^{(2)}(\eta )$$

$${\mathcal{S}}_{\mathbf{k}}^{(3)}(\eta )$$

$${\mathcal{S}}_{\mathbf{k}}^{(4)}(\eta )$$

$$f_2(x)$$

$$f_3(x)$$

$$f_4(x)$$

$$+4$$

$$\Omega_{\mathrm{GW}}(\eta,k)=\frac{A^{2}}{192\tilde{k}^{2}}\Big{[}\overline{I_{2}^{2}}M_{1}+A\left(M_{2}\overline{I_{3}^{2}}+M_{1}\overline{I_{2}I_{4}}\right)\Big{]},$$

$$I_i=\underset{x\to \infty }{\lim }{\int }_0^x\mathrm{d}\tilde{x}{f}_i\left(\frac{\tilde{x}}{\tilde{k}}\right)\frac{\tilde{x}}{x}\sin (x-\tilde{x}),$$

$$M_1(k)$$

$$M_2(k)$$

$$\Delta =$$

$$M_0=$$

$$\displaystyle\Big{[}(1-\nu_{i}^{2}$$

$$\mu$$

$${\nu }_{1,2}$$

$${\Omega }_{\mathrm{GW}}({\eta }_0,f)\simeq {\Omega }_r\times {\Omega }_{\mathrm{GW}}(\eta ,f),$$

$$\rho =\sqrt{T}{\left[\int \mathrm{d}f{\left(\frac{{\Omega }_{\mathrm{GW}}(f)}{{\Omega }_n(f)}\right)}^2\right]}^{1/2},$$

$$\rho =\sqrt{2T}{\left(\sum _{I,J}^M{\chi }_{IJ}^2\left[\int \mathrm{d}f{\left(\frac{{\Omega }_{\mathrm{GW}}(f)}{{\Omega }_n(f)+{\Omega }_{\mathrm{GW}}(f)}\right)}^2\right]\right)}^{1/2},$$

Full text

Various independent cosmological observations indicate the existence of dark matter (DM) in our Universe. The nature of DM remains highly elusive despite decades of dedicated searches. However, cosmological observations are only sensitive to the macroscopic properties of DM. Even though a new elementary particle is postulated in standard DM scenarios, primordial black hole (PBH) DM has attracted a lot of attention Bird:2016dcv ; Garcia-Bellido:2017fdg ; Sasaki:2018dmp ; Barack:2018yly , ever since the first direct detection of gravitational waves (GWs) from a binary black hole (BBH) coalescence Abbott:2016blz . The primordial-origin BBHs are appealing candidates of LIGO/Virgo BBHs if the abundance of stellar mass PBHs in DM is a few part in thousand Chen:2018czv ; Chen:2019irf ; Chen:2018rzo . There are various powerful observational constraints in literature Carr:2009jm ; Barnacka:2012bm ; Graham:2015apa ; Niikura:2017zjd ; Griest:2013esa ; Niikura:2019kqi ; Tisserand:2006zx ; Brandt:2016aco ; Gaggero:2016dpq ; Ali-Haimoud:2016mbv ; Blum:2016cjs ; Horowitz:2016lib ; Chen:2016pud ; Wang:2016ana ; Abbott:2018oah ; Magee:2018opb ; Wang:2019kaf ; Chen:2019irf ; Montero-Camacho:2019jte ; Laha:2019ssq , but a substantial window remains open for PBHs as all of DM in the approximate range $[10^{-16},10^{-14}]\cup[10^{-13},10^{-12}]M_{\odot}$. See a recent summary in Chen:2019irf .

PBHs are supposed to form from the enhancement of the scalar perturbations Hawking:1971ei ; Carr:1974nx . The process during which the PBHs are formed would be inevitably accompanied by GWs Tomita1967 ; Matarrese:1992rp ; Matarrese:1993zf ; Matarrese:1997ay ; Noh:2004bc ; Carbone:2004iv ; Nakamura:2004rm ; Ananda:2006af . These so-called induced GWs are driven by scalar perturbations during radiation-dominated (RD) era and could leave detectable signal at present for testing the hypothesis of PBH DM Saito:2008jc ; Bugaev:2010bb ; Sasaki:2018dmp ; Inomata:2018epa ; Baumann:2007zm ; Clesse:2018ogk ; Nakama:2016enz ; Saito:2009jt ; Bugaev:2009zh ; Assadullahi:2009jc . Recently, a semi-analytical expression for the GWs induced by second-order scalar perturbations has been derived in Espinosa:2018eve ; Kohri:2018awv , and the discussions have also been extended to detect the primordial non-Gaussianity Garcia-Bellido:2017aan ; Unal:2018yaa ; Cai:2018dig ; Cai:2019amo through induced GWs and the bispectrum Bartolo:2018evs ; Bartolo:2018rku ; Espinosa:2018eve of the induced GWs. PBHs could have been produced during RD era when relatively large scalar perturbations with amplitudes $\mathcal{O}(0.01-0.1)$ re-entered the Hubble horizon Ivanov:1994pa ; GarciaBellido:1996qt ; Ivanov:1997ia ; Yokoyama:1995ex ; Kawasaki:2006zv . The non-linearity is an intrinsic property for gravity in general relativity, and therefore the higher-order corrections to the scalar induced GWs are expected to be important. However, previous studies only focused on the GWs induced by second-order scalar perturbations, and as far as we know, no higher-order effects have been taken into account in literature. In this letter, we give the first result of the induced GWs up to the third-order, and find that the third-order correction would be detectable by some future GW detectors, such as LISA Audley:2017drz , IPTA Hobbs:2009yy , FAST Nan:2011um and SKA Kramer:2015jsa , and has significant observational implications.

The perturbed metric in the FRW spacetime with Newtonian gauge takes the form, Ananda:2006af ,

[TABLE]

where $\phi$ and $h_{ij}$ are the scalar and tensor perturbations respectively. The scalar perturbation in Fourier space has a solution Baumann:2007zm ; Kohri:2018awv ; Sasaki:2018dmp ; tomita1971

[TABLE]

where $\phi_{\bm{k}}$ is the primordial perturbation and $T(k\eta)$ is the transfer function

[TABLE]

which oscillates and decays as $\sim 1/\eta^{2}$ in the RD era. In order to extract the induced GWs up to the third-order, it is necessary to expand Einstein equations up to the fourth-order (perturbations are performed utilizing the xPand Pitrou:2013hga package). After a tedious but straightforwards computation, the evolution of the GWs up to the fourth-order is given by

[TABLE]

where the prime denotes the derivative with respect to the conformal time $\eta$, $\mathcal{H}\equiv a^{\prime}/a$ is the conformal Hubble parameter, and $\mathcal{T}_{ij}^{\ell m}$ is the projection operator Ando:2017veq onto the transverse and traceless tensor. Although Eq. (4) has the same form as the evolution of GWs at second-order (see e.g. Ando:2017veq ; Baumann:2007zm ), the source term $S_{ij}=S_{ij}^{(2)}+S_{ij}^{(3)}+S_{ij}^{(4)}$ has been calculated up to fourth-order as follows

[TABLE]

[TABLE]

[TABLE]

After solving Eq. (4) in Fourier space by Green’s function, we can use the two-point correlation of the GWs to calculate their power spectrum, namely

[TABLE]

The angle brackets stand for ensemble average which can be calculated using timing average instead. The GW energy density, $\rho_{\mathrm{GW}}=\int\rho_{\mathrm{GW}}(k,\eta)\,\mathrm{d}\ln k$, can be evaluated as Maggiore:1999vm

[TABLE]

where the overline stands for the average of the oscillating effect of the time-varying phase and $M_{p}$ is the Planck mass. The dimensionless GW energy density parameter $\Omega_{\mathrm{GW}}$ is defined as the energy density of GWs per logarithmic frequency normalized by the critical density $\rho_{\rm{cr}}$,

[TABLE]

where we have summed over the two polarization modes of $+$ and $\times$. Solving Eq. (4) by the Green’s function method, one obtains Baumann:2007zm

[TABLE]

where $\mathcal{S}_{\bm{k}}(\eta)\equiv-4e^{ij}(\bm{k})\tilde{S}_{ij}(\bm{k},\eta)$ with $\tilde{S}_{ij}(\bm{k},\eta)$ being the source term transformed into Fourier space. The polarization tensors $e_{ij}(\bm{k})$ are defined as $(e_{i}e_{j}-\bar{e}_{i}\bar{e}_{j})/\sqrt{2}$ and $(e_{i}\bar{e}_{j}+\bar{e}_{i}e_{j})/\sqrt{2}$ for $+$ and $\times$ polarizations respectively, where $e_{i}(\bm{k})$ and $\bar{e}_{i}(\bm{k})$ are two independent unit vectors orthogonal to $\bm{k}$. Since $\phi$ is related to the comoving curvature perturbation $\zeta$ by $\phi=(2/3)\zeta$ on superhorizon scales, $\Omega_{\mathrm{GW}}(\eta,k)$ can be calculated using the power spectrum of the comoving curvature perturbation, $\mathcal{P}_{\zeta}(k)$, defined by

[TABLE]

From now on we will dedicated to the following monochromatic power spectrum

[TABLE]

to illustrate the effects of third-order correction. Here $A$ is an overall normalization coefficient and $k_{*}$ is the wavenumber at which the power spectrum has a $\delta$-function peak. Note that it is robust to neglect the contribution of the long wavelength modes on the CMB scales, because not only the time delay effect of the long wavelength will not affect the GWs power spectrum Bartolo:2018rku , but also the long wavelength mode is well outside the horizon during the formation of PBHs and thus should not affect any local physical processes. On the other hand, the $\delta$-spectrum also corresponds to a monochromatic PBH formation. The formation of PBHs is a threshold process where the comoving curvature perturbation $\zeta{\bm{(x)}}$ exceeds a threshold value $\zeta_{c}$ and causes an overdensed region. The possibility of forming a PBH can be evaluated statistically by integrating the probability density function over the threshold region Carr:2016drx

[TABLE]

where $\zeta_{c}\simeq 1$ is the threshold value Musco:2008hv ; Musco:2004ak ; Musco:2012au ; Harada:2013epa to form a PBH and $\sigma^{2}\equiv\left\langle\zeta^{2}\right\rangle=\int\mathcal{P}_{\zeta}(k)\mathrm{d}\ln k=A$ is the variance of the curvature perturbation. For monochromatic PBHs, the possibility to form a single PBH, $\beta$, is equivalent to the abundance of PBHs which is related to the fraction of PBHs by Nakama:2016gzw

[TABLE]

where $g_{*}^{\mathrm{form}}$ is the effective degrees of freedom when PBHs are formed, and the mass of the PBH roughly equals to the horizon mass, namely

[TABLE]

where we have used $k_{*}=aH_{*}=2\pi f_{*}$ and $H_{0}$ is the Hubble constant at present.

For the $\delta$-spectrum in Eq. (13), the source terms in Fourier space are given by

[TABLE]

where we have defined $\mathbf{e}(\bm{p},\bm{q})\equiv e^{ij}(\bm{k})p_{i}q_{j}$, and $f_{i}(x)$ ($i=2,3,4$) have the following functional forms

[TABLE]

The explicit expression for transfer function $T=T(x)$ can be found in Eq. (3). From Eqs. (8) and (11), we see that only $\langle\mathcal{S}^{(2)}_{\bm{k}}\mathcal{S}^{(2)}_{\bm{k^{\prime}}}\rangle$ contributes to second-order induced GWs. Meanwhile, both $\langle\mathcal{S}^{(3)}_{\bm{k}}\mathcal{S}^{(3)}_{\bm{k^{\prime}}}\rangle$ and $\langle\mathcal{S}^{(2)}_{\bm{k}}\mathcal{S}^{(4)}_{\bm{k^{\prime}}}\rangle$ contribute to the third-order correction. Following the pioneering work of Espinosa:2018eve ; Kohri:2018awv , for the $\delta$-spectrum in Eq. (13), we obtain

[TABLE]

where an overbar denotes the oscillation average Kohri:2018awv and $I_{i}$ ($i=2,3,4$) are defined as

[TABLE]

which reflects the phase oscillation of the GWs. For convenience, we have defined some dimensionless parameters, i.e. $\tilde{k}\equiv k/k_{*}$, $x\equiv k\eta$ and $\tilde{x}\equiv k\tilde{\eta}$. Similar to Kohri:2018awv , Eq. (24) can be analytically integrated by multiple usages of the trigonometric addition theorem and the properties of sine integral $\rm{Si}(\theta)$ and cosine integral $\rm{Ci}(\theta)$. The angle integrals $M_{1}$ and $M_{2}$ in Eq. (23) are defined as

[TABLE]

where $\Theta$ is the Heaviside step function, $p_{\rm{min}}=|1-\tilde{k}|$, and $p_{\rm{max}}=\min(2,1+\tilde{k})$. Due to its complexity, the expression of $M_{2}$ is evaluated numerically and the definition of $M_{0}$ and $\Delta$ are given by

[TABLE]

where $\mu$ and $\nu_{i}$ ($i=1,2$) are defined as

[TABLE]

Note that $M_{1}$ is originated from the $\langle\mathcal{S}^{(2)}_{\bm{k}}\mathcal{S}^{(2)}_{\bm{k^{\prime}}}\rangle$ and $\langle\mathcal{S}^{(2)}_{\bm{k}}\mathcal{S}^{(4)}_{\bm{k^{\prime}}}\rangle$ terms, and indicates a cutoff frequency at $k=2k_{*}$; while $M_{2}$ is originated from the $\langle\mathcal{S}^{(3)}_{\bm{k}}\mathcal{S}^{(3)}_{\bm{k^{\prime}}}\rangle$, and indicates a cutoff at $k=3k_{*}$. Therefore the third-order correction not only enhances the amplitude of GW energy density, but also extends the cutoff frequency from $2k_{*}$ to $3k_{*}$.

Note that Eq. (23) is only valid from the horizon re-entry to matter-radiation equality. Because the energy density of GWs decays as radiation, current density parameter of GWs can be approximated by Espinosa:2017sgp

[TABLE]

where $\Omega_{r}$ is the density parameter of radiation at present. From Eq. (16), heavier masses of PBHs corresponds to lower peak frequencies of induced GWs. Fig. 1 shows the GWs induced by scalar perturbation up to third order compared with the sensitivity curves of LISA, IPTA, FAST and SKA. Because a substantial windows for PBH as all of DM in the approximate range $[10^{-16},10^{-14}]\cup[10^{-13},10^{-12}]M_{\odot}$ is still available, for example, we choose $m_{\rm{pbh}}=10^{-12}M_{\odot}$ and $f_{\mathrm{pbh}}=1$ and then the peak frequency of the induced GW is $f\sim 10^{-3}$ Hz which is within the LISA frequency band. In addition, roughly speaking, since the abundance of stellar mass PBHs has been constrained to be less than $0.01$ Chen:2019irf , we choose $m_{\rm{pbh}}=0.2M_{\odot}$ and $f_{\mathrm{pbh}}=10^{-3}$, and then the peak frequency of induced GWs is just located at the PTA frequency band. In fact, since Eqs. (14) and (15) imply $f_{\rm{pbh}}\propto\sqrt{A}\exp(-A^{-1})$ for a relatively small value of $A$, even decreasing $f_{\mathrm{pbh}}$ by multiple orders, the consequent value of $A$ will almost not change. Therefore the amplitude of $\Omega_{\mathrm{GW}}$ will almost be irrespective with the choice of $f_{\mathrm{pbh}}$ unless $f_{\mathrm{pbh}}$ is dramatically much smaller than our choice in Fig. 1. From Fig. 1, we also see that the third-order correction not only enhances the magnitude of $\Omega_{\mathrm{GW}}$ (especially smooth the deep valley at $f=\sqrt{2/3}f_{*}$ from the second-order effect), but also extends cutoff frequency from $2f_{*}$ to $3f_{*}$. Besides, for a narrow spectrum, it would also generate two new resonant peaks at $f=(1/\sqrt{3})f_{*}$ and $f=\sqrt{3}f_{*}$. Moreover, the typical $\Omega_{\mathrm{GW}}(\eta_{0},f)$ in Fig. 1 shows that the third-order correction of the induced GWs is also far beyond the sensitivity curves for both LISA and PTA observations, rendering third-order correction being detectable and thus, indispensable in testing the hypothesis of PBHs.

In order to quantitatively evaluate the effects of the third-order correction in observations, we need to estimate the SNR, $\rho$, for different GW experiments. For LISA, it is given by Thrane:2013oya

[TABLE]

where $\Omega_{n}(f)=2\pi^{2}f^{3}S_{n}/(3H_{0}^{2})$ and $S_{n}$ is the strain noise power spectral density Cornish:2018dyw . For PTA observations, if we assume pulsars are distributed homogeneously and all pulsars have the same noise characteristics, the SNR is given by Siemens:2013zla

[TABLE]

where $\chi_{IJ}$ is the Hellings and Downs coefficient for pulsars $I$ and $J$ Hellings:1983fr . Fig. 2 shows the expected SNR obtained for LISA and PTA experiments, and indicates that the third-order correction would raise the SNR as expected. For LISA, we expect around a $15\%$ increase in the relative SNR; while for PTA observations, we expect the increase is more than $5\%$ and could even reach $20\%$ for stellar mass PBHs. Most importantly, since the third-order correction will extend the cutoff frequency from $2k_{*}$ to $3k_{*}$, the induced GWs accompanying heavier PBHs are supposed to be detected by PTA observations. For instance, as shown in Fig. 2, the third-order correction will extend the largest detectable mass of PBHs from around $30M_{\odot}$ to $65M_{\odot}$, which will be an invaluable complimentary tool to test the PBH scenario in addition to the analysis LIGO/Virgo.

In this letter, we compute the third-order correction to the induced GWs generated by the scalar perturbations accompanying the formation of PBHs during the RD era. After deriving a general expression for the GW energy density, we investigate an infinite narrow power spectrum of the scalar perturbation and obtain a semi-analytical expression for $\Omega_{\mathrm{GW}}$. Our result implies that the third-order correction to the induced GWs accompanying formation of PBH DM will generate detectable effects on the waveform and the amplitude in observation data. The third-order correction not only enhances the magnitude of $\Omega_{\mathrm{GW}}$, but also extends the cutoff frequency from $2f_{*}$ to $3f_{*}$. We also forecast the SNR for LISA and PTA observations, including IPTA, FAST and SKA. These planned GW projects cover a wide frequency band from $10^{-9}$ Hz to $10^{-3}$ Hz corresponding to PBHs with mass range $10^{-18}M_{\odot}-10^{2}M_{\odot}$. Our results indicate that the third-order correction could not only lead to an increase in the relative SNR around $10\%$ for these planned GW projects, but also extend the maximum detectable PBH mass. On the other hand, if all of these projects are unable to detect such induced GWs, we could rule out PBH DM in a wide mass ranges of $\sim 10^{-18}M_{\odot}-10M_{\odot}$.

Acknowledgments. We would like to thank Lu Chen, Yun Fang, Fan Huang, Jun Li, Lang Liu, Shi Pi, You Wu, Yu Sang, Sai Wang, Hao Wei and Xue Zhang for useful conversations. CY is indebted to Yu-Jia Zhai for her continuous support and encouragement. We acknowledge the use of HPC Cluster of ITP-CAS. This work is supported by grants from NSFC (grant No. 11690021, 11575271, 11747601), the Strategic Priority Research Program of Chinese Academy of Sciences (Grant No. XDB23000000, XDA15020701), and Top-Notch Young Talents Program of China. This research has made use of GWSC.jl gwsc package to calculate the SNR for various gravitational-wave detectors.