Probing Primordial-Black-Hole Dark Matter with Scalar Induced Gravitational Waves
Chen Yuan Zu-Cheng Chen Qing-Guo Huang






     
   
             
         
     



      


 




 
B_s


 


α
σ
δ
Δ
Λ
∂
ω
Ω
d
[1]𝒪(#1) 
κ


()

    [
]

    [1]Eq. (<ref>)
        [1]Fig. <ref>
        [1]Table <ref>
        [1]Sec. <ref>
        M_⊙
        f_pbh
        f_pbh0
        ℛ
        σ_eqΩ_GWGpc^-3 yr^-1LIGO/Virgo
        SNR
        m_minm_maxM_minf_minVTρ_GWθ⃗d⃗λ⃗N_obs[1]⟨ #1 ⟩
        kmMpcT_obsN_temp
        [add ref] 
        e.g. ≈1/2Add discussion here![1]#1
        ℋc_s^2[1]S_ij^(#1)[1]v_i^(#1)⟨⟩α1 2kqplx𝐞𝒮η̃p̃k̃x̃F̃Ã|--||-||-|[1]^(#1)(, η)k⃗k_*x_*m_pbh
        
        
        	
        
        
        yuanchen@itp.ac.cnCAS Key Laboratory of Theoretical Physics, 
        	Institute of Theoretical Physics, Chinese Academy of Sciences,
        	Beijing 100190, ChinaSchool of Physical Sciences, 
        	University of Chinese Academy of Sciences, 
        	No. 19A Yuquan Road, Beijing 100049, China
        chenzucheng@itp.ac.cn
        CAS Key Laboratory of Theoretical Physics, 
        	Institute of Theoretical Physics, Chinese Academy of Sciences,
        	Beijing 100190, ChinaSchool of Physical Sciences, 
        	University of Chinese Academy of Sciences, 
        	No. 19A Yuquan Road, Beijing 100049, China
        huangqg@itp.ac.cnCAS Key Laboratory of Theoretical Physics, 
        	Institute of Theoretical Physics, Chinese Academy of Sciences,
        	Beijing 100190, ChinaSchool of Physical Sciences, 
        	University of Chinese Academy of Sciences, 
        	No. 19A Yuquan Road, Beijing 100049, ChinaCenter for Gravitation and Cosmology, 
        	College of Physical Science and Technology, 
        	Yangzhou University, Yangzhou 225009, ChinaSynergetic Innovation Center for Quantum Effects and Applications, 
        	Hunan Normal University, Changsha 410081, China
        
        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. 
        		
        
        ???Probing Primordial-Black-Hole Dark Matter with Scalar Induced Gravitational Waves
            Qing-Guo Huang
            December 30, 2023
        =================================================================================
        
        
        
        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 <cit.>,
        ever since the first direct detection of gravitational waves (GWs) from a binary black hole (BBH)
        coalescence <cit.>. 
        The primordial-origin BBHs are appealing candidates of  BBHs if the abundance of stellar mass PBHs in DM is a few part in thousand <cit.>. 
        There are various powerful observational constraints in literature <cit.>, but a substantial window remains open for PBHs as all of DM in the approximate range [10^-16,10^-14] ∪ [10^-13,10^-12] M_⊙. See a recent summary in <cit.>. 
        
        
        
        
        
        
        
        
        
        	
        PBHs are supposed to form from the enhancement of the scalar perturbations <cit.>. The process during which the PBHs are formed would be inevitably accompanied by GWs <cit.>.
        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  <cit.>. 
        Recently, a semi-analytical expression for the GWs induced by second-order scalar perturbations
        has been derived in <cit.>, and the discussions have also been extended to detect the primordial non-Gaussianity 
        <cit.> 
        through induced GWs and the bispectrum
        <cit.> of the induced GWs.
        PBHs could have been produced during RD era when relatively large scalar perturbations with amplitudes 𝒪(0.01-0.1) re-entered the Hubble horizon <cit.>. 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 <cit.>, IPTA <cit.>, FAST <cit.> and SKA <cit.>, and has significant observational implications.
        
        	
        The perturbed metric in the FRW spacetime with Newtonian gauge takes the form, <cit.>, 
         s^2=a^2{-(1+2 ϕ) η^2+[(1-2 ϕ) δ_i j+h_i j/2]  x^i x^j},
        
        where ϕ and h_ij are the scalar and tensor perturbations respectively.
        
        The scalar perturbation in Fourier space has a solution <cit.>ϕ_k⃗(η) ≡ϕ_k⃗T(kη),
        
        where ϕ_k⃗ is the primordial perturbation and T(kη) is the transfer function
        
        	T(kη) =  9/(kη)^2
            sin(kη/√(3))/kη/√(3) - cos(kη/√(3))
        ,
        
        which oscillates and decays as ∼1/η^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 <cit.> package). 
        After a tedious but straightforwards computation, the evolution of the GWs up to the fourth-order is given by
        
        h_i j^''+2 ℋ h_i j^'-∇^2 h_i j=-4 𝒯_i j^ℓ m S_ℓ m,
        
        where the prime denotes the derivative with respect to the conformal time η, 
        ≡ a'/a is the conformal Hubble parameter, and 
        𝒯_i j^ℓ m is the projection operator <cit.> 
        onto the transverse and traceless tensor.
        Although eqh has the same form as the evolution of GWs at second-order 
        (see <cit.>),
        the source term S_ij=2+3+4 has been calculated up to fourth-order as follows 
        
        2= 4 ϕ_i_jϕ + 2_iϕ_jϕ-_i ϕ + ϕ'_jϕ + ϕ',
        
        
        
        
        3 =   1/12ϕ - ϕ'_iϕ_jϕ
        	- 1/^34ϕ - ϕ'_iϕ'_jϕ' 
           + 1/3^42^2ϕ - 9ϕ'_iϕ +ϕ'_jϕ +ϕ',
        
        
        
        
        
        
        4 =    16 ϕ^3 _i_jϕ
        + 1/3^3[2 ϕ' ^2ϕ - 9ϕ'^2 
        - 8ϕ^2ϕ
                         + 18 ^2 ϕϕ'
        + 96^3 ϕ^2] _iϕ_jϕ
           + 2/3^5[- ϕ' ^2ϕ + 3ϕ'^2 
        + 4ϕ^2ϕ
                           + 3 ^2 ϕϕ'- 12^3 ϕ^2] _iϕ' _jϕ'
           + 1/36^6[ 
        -16 (^2ϕ)^2 - 3  ∂_kϕ' ∂^kϕ'
        + 120 ϕ' ^2ϕ
                  - 6 ∂_kϕ∂^kϕ'	+ 144 ^2 ϕ^2ϕ - 180 ^2 ϕ'^2 
                    + 33 ^2 ∂_kϕ∂^kϕ- 504 ^3 ϕϕ'
        - 144 ^4 ϕ^2
        ] 
                   ×_iϕ +ϕ'_jϕ +ϕ'.
        
        
        
        	
        After solving Eq. (<ref>) in Fourier space by Green's function, we can use the two-point correlation of the GWs to calculate their power spectrum, namely 
        ⟨ h(k⃗,η) h(k⃗'⃗,η)⟩≡2 π^2/k^3𝒫_h(k,η) δ(k⃗+k⃗'⃗).
        
        The angle brackets stand for ensemble average which can be calculated using timing average instead. The GW energy density, ρ_GW=∫ρ_GW(k,η) ln k, 
        can be evaluated as <cit.>ρ_GW=M_p^216a^2∂_kh_ij∂^kh^ij,
        
        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  is defined as the energy density of GWs per logarithmic frequency normalized by the critical density ρ_cr,
        (η, k)≡ρ_GW(k,η)/ρ_cr = 1/12k/^2 𝒫_h(k, η),
        
        where we have summed over the two polarization modes of + and ×.
        Solving eqh by the Green's function method, one obtains <cit.> 
        	h(k⃗,η) = 1/k a(η)∫sin(kη - k) a () 𝒮_k⃗(),
         
        where 𝒮_k⃗(η)≡ -4e^ij(k⃗) S̃_ij(k⃗,η) with S̃_ij(k⃗,η) being the source term transformed into Fourier space.
        The polarization tensors e_ij(k⃗) are defined as 
        (e_i e_j - e̅_i e̅_j)/√(2) and (e_i e̅_j + e̅_i e_j)/√(2) 
        for + and × polarizations respectively, 
        where e_i() and e̅_i() are two independent unit vectors orthogonal to .
        Since ϕ is related to the comoving curvature perturbation ζ by ϕ = (2/3) ζ on superhorizon scales, (η, k) can be calculated using the power spectrum of the comoving curvature perturbation, 𝒫_ζ(k), defined by 
        ⟨ζ_ζ_k⃗'⃗⟩≡2 π^2/k^3𝒫_ζ(k) δ(k⃗+k⃗'⃗). 
        
        
        
        From now on we will dedicated to the following monochromatic power spectrum
        𝒫_ζ(k)=A δk-,
        
        to illustrate the effects of third-order correction.
        Here A is an overall normalization coefficient and  is the wavenumber at which the
        power spectrum has a δ-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 <cit.>, 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 δ-spectrum also corresponds to a monochromatic PBH formation. 
        The formation of PBHs is a threshold process where the comoving curvature perturbation ζ(⃗x⃗)⃗ exceeds a threshold value ζ_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 <cit.>β=∫_ζ_c^+∞dζ/√(2 π)σ e^-ζ^2 / 2 σ^2=erfcζ_c√(2A),
        
        where ζ_c≃1 is the threshold value <cit.> to form a PBH and σ^2≡ζ^2=∫𝒫_ζ(k) dln k =A is the variance of the curvature perturbation. 
        For monochromatic PBHs, the possibility to form a single PBH, β, is equivalent to the abundance of PBHs which is related to the fraction of PBHs by <cit.>
        	f_pbh≃ 2.5 × 10^8β(g_*^form/10.75)^-1/4(m_pbh/M_⊙)^-1/2,
        
        where g_*^form is the effective degrees of freedom when PBHs are formed, 
        
        
        and the mass of the PBH roughly equals to the horizon mass, namely 
        m_pbh M_⊙≈2.3×10^18(3.91/g_*^form)^1 / 6(H_0/f_*)^2,
        
        where we have used k_*=aH_*=2π f_* and H_0 is the Hubble constant at present.
        
        For the δ-spectrum in pzeta, the source terms in Fourier space are given by
        𝒮^(2)_k⃗(η)   ≡   ∫^3p/(2π)^3/2(, )
        	f_2(η) ζ_ζ_k-p,
        𝒮^(3)_k⃗(η)   ≡   ∫^3p ^3q/(2π)^3(, )
        	f_3(η) ζ_ζ_ζ_k-p-q,
        𝒮^(4)_k⃗(η)   ≡   ∫^3p ^3q ^3l/(2π)^9/2
            (, )+ (,)
         f_4(η) 
                          ×ζ_ζ_ζ_ζ_k-p-q-l,
        
        where we have defined (, ) ≡ e^ij(k⃗)p_iq_j, and f_i(x) (i = 2, 3, 4) 
        have the following functional forms
        
        	f_2(x)      =893T^2 + 2xTT' + x^2T'^2,
        
        	f_3(x)      = -64/81[ x^2-18 T^3 + 2x 3+x^2 T^2 T'
        	
                     + x^2 15 + x^2 T T'^2 + 3 x^3 T'^3 ], 
        
        	f_4(x)      = 16729[720-29x^2+16x^4T^4
                    +4x144+73x^2+8x^4T^3T'
                    +2x^2864+219x^2+8x^4T^2T'^2
        
        	+4      x^3198 + 31 x^2TT'^3+x^4 108 + 7 x^2 T'^4].
        
        The explicit expression for transfer function T=T(x) can be found in transfer. 
        From Eqs. (<ref>) and (<ref>), we see that only 𝒮^(2)_k⃗𝒮^(2)_k⃗'⃗ contributes to second-order induced GWs.
        Meanwhile, both 𝒮^(3)_k⃗𝒮^(3)_k⃗'⃗ 
        and 𝒮^(2)_k⃗𝒮^(4)_k⃗'⃗ contribute to the third-order correction. 
        Following the pioneering work of <cit.>, for the -spectrum in pzeta, we obtain 
        (η, k)=A^2/192^2[I_2^2M_1+
        	A M_2I_3^2+M_1I_2 I_4],
        
        where an overbar denotes the oscillation average <cit.> and  I_i (i = 2, 3, 4) are defined as
        
        	I_i = lim_x →∞∫_0^x   f_i/ xsin(x-),
        
        which reflects the phase oscillation of the GWs. 
        For convenience, we have defined some dimensionless parameters, i.e. ≡ k/, x≡ kη and ≡ kη̃. 
        Similar to <cit.>, Ii can be analytically integrated by multiple usages of
        the trigonometric addition theorem and the properties of sine integral Si(θ) 
        and cosine integral Ci(θ).
        
        The angle integrals M_1 and M_2 in ogw are defined as
        
        	M_1(k)    =   4 - ^2^2 Θ(2 - ),
        
        	M_2(k)    =   1/π^2∫_p_min^p_max∫_0^2π∫_0^2πϕ M_0 Θ(Δ),
        
        where Θ is the Heaviside step function,
        p_min = |1-|, and p_max = min(2,1+).
        Due to its complexity, the expression of M_2 is evaluated numerically and 
        the definition of M_0 and Δ are given by 
        Δ =       4μ^2+4(1-μ^2)cos(-ϕ)^2-p̃^2,
        
        M_0 =      ∑_i=1^2(1-ν_i^2)/|μ√(1-ν_i^2) 
        	- ν_i √(1-μ^2)cos(-ϕ)|
        [ 
        (1-ν_i^2      )^32cos^22 + 2^3 (1-μ^2)^32cos2ϕcos(+ϕ)
              - 2  (1-ν_i^2) (1-μ^2)^cos2cos(+ϕ)
               - ^2 (1-μ^2) (1-ν_i^2)^cos(+ϕ)^2 ],
        
        where μ and ν_i (i=1,2) are defined as
        μ      = ^2+^2-1  2, 
        ν_1,2      = μ±√(1-μ^2)|cos(-ϕ)|√(Δ)/2(1-μ^2)cos(-ϕ)^2+μ^2.
        
        Note that M_1 is originated from the 
        𝒮^(2)_k⃗𝒮^(2)_k⃗'⃗ and 
        𝒮^(2)_k⃗𝒮^(4)_k⃗'⃗ terms, 
        and indicates a cutoff frequency at k=2;
        while M_2 is originated from the 
        𝒮^(3)_k⃗𝒮^(3)_k⃗'⃗,
        and indicates a cutoff at k=3.
        Therefore the third-order correction not only enhances the amplitude of GW energy density,
        but also extends the cutoff frequency from 2 to 3.
        
        
        
        
        
        
        Note that ogw 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 <cit.>(η_0, f) ≃Ω_r×(η, f),
        
        where Ω_r is the density parameter of radiation at present.
        
        From mkrelation, heavier masses of PBHs corresponds to lower peak frequencies of induced GWs. 
        
        OmegaGW 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]∪ [10^-13,10^-12] M_⊙ is still available, for example, we choose =10^-12 and =1 and then the peak frequency of the induced GW is f∼ 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<cit.>, we choose =0.2 and =10^-3, and then the peak frequency of induced GWs is just located at the PTA frequency band. 
        In fact, since Eqs. (<ref>) and (<ref>) imply f_pbh∝√(A)exp(-A^-1) for a relatively small value of A, 
        even decreasing  by multiple orders, the consequent value of A will almost not change.
        Therefore the amplitude of  will almost be irrespective with the choice of  unless  is dramatically much smaller than our choice in OmegaGW. 
        From OmegaGW, we also see that the third-order correction not only enhances the magnitude of
         (especially smooth the deep valley at f=√(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/√(3))f_* and f=√(3)f_*. 
        
        Moreover, the typical (η_0, f) in OmegaGW 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, ρ, for different GW experiments.
        For LISA, it is given by <cit.>ρ = √(T)
            ∫ f (f)/_n(f)^2
        ^1/2,
        
        where _n(f) = 2π^2 f^3 S_n/(3H_0^2)
        and S_n is the strain noise power spectral density <cit.>. 
        For PTA observations, if we assume pulsars are distributed homogeneously
        and all pulsars have the same noise characteristics, the SNR is given by <cit.>ρ = √(2T)∑_I, J^Mχ_IJ^2
            ∫ f (f)/_n(f) + (f)^2^1/2,
        
        where χ_IJ is the Hellings and Downs coefficient for pulsars I and J<cit.>.
        SNR 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 SNR, the third-order correction will extend the largest detectable mass of PBHs from around 30 to 65, which will be an invaluable complimentary tool to
        test the PBH scenario in addition to the analysis .
        
        
        	
        
        	
        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 . 
        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 , 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^-18M_⊙-10^2M_⊙. 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 ∼10^-18M_⊙-10M_⊙. 
        	
        
        
        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 <cit.> package to calculate
        	the SNR for various gravitational-wave detectors.
        
        	
        	
        
        
        99Bird:2016dcv 
        S. Bird, I. Cholis, J. B. Muoz, Y. Ali-Hamoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli and A. G. Riess,
        
        Phys. Rev. Lett.  116, no. 20, 201301 (2016)
        
        [arXiv:1603.00464 [astro-ph.CO]].
        		
        Garcia-Bellido:2017fdg 
        J. Garca-Bellido,
        
        J. Phys. Conf. Ser.  840, no. 1, 012032 (2017)
        
        [arXiv:1702.08275 [astro-ph.CO]].
        
        Barack:2018yly 
        L. Barack  et al.,
        
        arXiv:1806.05195 [gr-qc].
        
        Sasaki:2018dmp 
        M. Sasaki, T. Suyama, T. Tanaka and S. Yokoyama,
        
        Class. Quant. Grav.  35, no. 6, 063001 (2018)
        
        [arXiv:1801.05235 [astro-ph.CO]].
        
        Abbott:2016blz 
        B. P. Abbott  et al. [LIGO Scientific and Virgo Collaborations],
        
        Phys. Rev. Lett.  116, no. 6, 061102 (2016)
        
        [arXiv:1602.03837 [gr-qc]].
        
        Chen:2018czv 
        Z. C. Chen and Q. G. Huang,
        
        Astrophys. J.  864, no. 1, 61 (2018)
        
        [arXiv:1801.10327 [astro-ph.CO]].
        
        Chen:2018rzo 
        Z. C. Chen, F. Huang and Q. G. Huang,
        
        Astrophys. J.  871, 97 (2019)
        
        [arXiv:1809.10360 [gr-qc]].
        
        Chen:2019irf 
        Z. C. Chen and Q. G. Huang,
        
        arXiv:1904.02396 [astro-ph.CO].
        
        Tisserand:2006zx 
        P. Tisserand  et al. [EROS-2 Collaboration],
        
        Astron. Astrophys.  469, 387 (2007)
        
        [astro-ph/0607207].
        
        Carr:2009jm 
        B. J. Carr, K. Kohri, Y. Sendouda and J. Yokoyama,
        
        Phys. Rev. D  81, 104019 (2010)
        
        [arXiv:0912.5297 [astro-ph.CO]].
        
        Barnacka:2012bm 
        A. Barnacka, J. F. Glicenstein and R. Moderski,
        
        Phys. Rev. D  86, 043001 (2012)
        
        [arXiv:1204.2056 [astro-ph.CO]].
        
        Griest:2013esa 
        K. Griest, A. M. Cieplak and M. J. Lehner,
        
        Phys. Rev. Lett.  111, no. 18, 181302 (2013).
        
        Graham:2015apa 
        P. W. Graham, S. Rajendran and J. Varela,
        
        Phys. Rev. D  92, no. 6, 063007 (2015)
        
        [arXiv:1505.04444 [hep-ph]].
        
        Brandt:2016aco 
        T. D. Brandt,
        
        Astrophys. J.  824, no. 2, L31 (2016)
        
        [arXiv:1605.03665 [astro-ph.GA]].
        
        Chen:2016pud 
        L. Chen, Q. G. Huang and K. Wang,
        
        JCAP  1612, no. 12, 044 (2016)
        
        [arXiv:1608.02174 [astro-ph.CO]].
        
        Wang:2016ana 
        S. Wang, Y. F. Wang, Q. G. Huang and T. G. F. Li,
        
        Phys. Rev. Lett.  120, no. 19, 191102 (2018)
        doi:10.1103/PhysRevLett.120.191102
        [arXiv:1610.08725 [astro-ph.CO]].
        
        Gaggero:2016dpq 
        D. Gaggero, G. Bertone, F. Calore, R. M. T. Connors, M. Lovell, S. Markoff and E. Storm,
        
        Phys. Rev. Lett.  118, no. 24, 241101 (2017)
        
        [arXiv:1612.00457 [astro-ph.HE]].
        
        Ali-Haimoud:2016mbv 
        Y. Ali-Haïmoud and M. Kamionkowski,
        
        Phys. Rev. D  95, no. 4, 043534 (2017)
        
        [arXiv:1612.05644 [astro-ph.CO]].
        
        Blum:2016cjs 
        D. Aloni, K. Blum and R. Flauger,
        
        JCAP  1705, no. 05, 017 (2017)
        
        [arXiv:1612.06811 [astro-ph.CO]].
        
        Horowitz:2016lib 
        B. Horowitz,
        
        arXiv:1612.07264 [astro-ph.CO].
        
        Niikura:2017zjd 
        H. Niikura  et al.,
        
        Nat. Astron.  3, no. 6, 524 (2019)
        
        [arXiv:1701.02151 [astro-ph.CO]].
        
        Abbott:2018oah 
        B. P. Abbott  et al. [LIGO Scientific and Virgo Collaborations],
        
        Phys. Rev. Lett.  121, no. 23, 231103 (2018)
        
        [arXiv:1808.04771 [astro-ph.CO]].
        
        Magee:2018opb 
        R. Magee  et al.,
        
        Phys. Rev. D  98, no. 10, 103024 (2018)
        
        [arXiv:1808.04772 [astro-ph.IM]].
        
        Niikura:2019kqi 
        H. Niikura, M. Takada, S. Yokoyama, T. Sumi and S. Masaki,
        
        Phys. Rev. D  99, no. 8, 083503 (2019)
        
        [arXiv:1901.07120 [astro-ph.CO]].
        
        Wang:2019kaf 
        S. Wang, T. Terada and K. Kohri,
        
        Phys. Rev. D  99, no. 10, 103531 (2019)
        
        [arXiv:1903.05924 [astro-ph.CO]].
        
        Montero-Camacho:2019jte 
        P. Montero-Camacho, X. Fang, G. Vasquez, M. Silva and C. M. Hirata,
        
        arXiv:1906.05950 [astro-ph.CO].
        
        Laha:2019ssq 
        R. Laha,
        
        arXiv:1906.09994 [astro-ph.HE].
        
        Hawking:1971ei 
          S. Hawking,
          
          Mon. Not. Roy. Astron. Soc.  152, 75 (1971).
        
        Carr:1974nx 
          B. J. Carr and S. W. Hawking,
          
          Mon. Not. Roy. Astron. Soc.  168, 399 (1974).
        
        Tomita1967
        K. Tomita, Prog. Theor. Phys.  37, 831 (1967).
        
        Matarrese:1992rp 
        S. Matarrese, O. Pantano and D. Saez,
        
        Phys. Rev. D  47, 1311 (1993).
        
        
        Matarrese:1993zf 
        S. Matarrese, O. Pantano and D. Saez,
        
        Phys. Rev. Lett.  72, 320 (1994)
        [astro-ph/9310036].
        
        Matarrese:1997ay 
        S. Matarrese, S. Mollerach and M. Bruni,
        
        Phys. Rev. D  58, 043504 (1998)
        [astro-ph/9707278].
        
        Noh:2004bc
        H. Noh and J. c. Hwang,
        
        Phys. Rev. D  69 (2004) 104011.
        Carbone:2004iv
        C. Carbone and S. Matarrese,
        
        Phys. Rev. D  71 (2005) 043508
        
        [astro-ph/0407611].
        
        Nakamura:2004rm 
        K. Nakamura,
        
        Prog. Theor. Phys.  117, 17 (2007)
        
        [gr-qc/0605108].
        
        Ananda:2006af 
        K. N. Ananda, C. Clarkson and D. Wands,
        
        Phys. Rev. D  75, 123518 (2007)
        
        [gr-qc/0612013].
        
        Baumann:2007zm 
        D. Baumann, P. J. Steinhardt, K. Takahashi and K. Ichiki,
        
        Phys. Rev. D  76, 084019 (2007)
        
        [hep-th/0703290].
        
        Saito:2008jc 
        R. Saito and J. Yokoyama,
        
        Phys. Rev. Lett.  102, 161101 (2009)
        Erratum: [Phys. Rev. Lett.  107, 069901 (2011)]
        
        [arXiv:0812.4339 [astro-ph]].
        
        Assadullahi:2009jc 
        H. Assadullahi and D. Wands,
        
        Phys. Rev. D  81, 023527 (2010)
        doi:10.1103/PhysRevD.81.023527
        [arXiv:0907.4073 [astro-ph.CO]].
        
        Bugaev:2009zh 
        E. Bugaev and P. Klimai,
        
        Phys. Rev. D  81, 023517 (2010)
        
        [arXiv:0908.0664 [astro-ph.CO]].
        
        Saito:2009jt 
        R. Saito and J. Yokoyama,
        
        Prog. Theor. Phys.  123, 867 (2010)
        Erratum: [Prog. Theor. Phys.  126, 351 (2011)]
        
        [arXiv:0912.5317 [astro-ph.CO]].
        
        Bugaev:2010bb 
        E. Bugaev and P. Klimai,
        
        Phys. Rev. D  83, 083521 (2011)
        
        [arXiv:1012.4697 [astro-ph.CO]].
        
        Nakama:2016enz 
        T. Nakama and T. Suyama,
        
        Phys. Rev. D  94, no. 4, 043507 (2016)
        
        [arXiv:1605.04482 [gr-qc]].
        
        Inomata:2018epa 
        K. Inomata and T. Nakama,
        
        Phys. Rev. D  99, no. 4, 043511 (2019)
        
        [arXiv:1812.00674 [astro-ph.CO]].
        
        Clesse:2018ogk 
        S. Clesse, J. García-Bellido and S. Orani,
        
        arXiv:1812.11011 [astro-ph.CO].
        
        Espinosa:2018eve 
        J. R. Espinosa, D. Racco and A. Riotto,
        
        JCAP  1809, no. 09, 012 (2018)
        
        [arXiv:1804.07732 [hep-ph]].
        
        Kohri:2018awv 
        K. Kohri and T. Terada,
        
        Phys. Rev. D  97, no. 12, 123532 (2018)
        
        [arXiv:1804.08577 [gr-qc]].
        
        Garcia-Bellido:2017aan 
        J. Garcia-Bellido, M. Peloso and C. Unal,
        
        JCAP  1709, no. 09, 013 (2017)
        doi:10.1088/1475-7516/2017/09/013
        [arXiv:1707.02441 [astro-ph.CO]].
        
        Cai:2018dig 
        R. G. Cai, S. Pi and M. Sasaki,
        
        Phys. Rev. Lett.  122, no. 20, 201101 (2019)
        doi:10.1103/PhysRevLett.122.201101
        [arXiv:1810.11000 [astro-ph.CO]].
        
        Unal:2018yaa 
        C. Unal,
        
        Phys. Rev. D  99, no. 4, 041301 (2019)
        
        [arXiv:1811.09151 [astro-ph.CO]].
        
        Cai:2019amo 
        R. G. Cai, S. Pi, S. J. Wang and X. Y. Yang,
        
        JCAP  1905, no. 05, 013 (2019)
        
        [arXiv:1901.10152 [astro-ph.CO]].
        
        Bartolo:2018evs 
        N. Bartolo, V. De Luca, G. Franciolini, A. Lewis, M. Peloso and A. Riotto,
        
        Phys. Rev. Lett.  122, no. 21, 211301 (2019)
        
        [arXiv:1810.12218 [astro-ph.CO]].
        
        Bartolo:2018rku 
        N. Bartolo, V. De Luca, G. Franciolini, M. Peloso, D. Racco and A. Riotto,
        
        Phys. Rev. D  99, no. 10, 103521 (2019)
        
        [arXiv:1810.12224 [astro-ph.CO]].
        
        Ivanov:1994pa 
        P. Ivanov, P. Naselsky and I. Novikov,
        
        Phys. Rev. D  50, 7173 (1994).
        
        Yokoyama:1995ex 
        J. Yokoyama,
        
        Astron. Astrophys.  318, 673 (1997)
        [astro-ph/9509027].
        
        GarciaBellido:1996qt 
        J. Garcia-Bellido, A. D. Linde and D. Wands,
        
        Phys. Rev. D  54, 6040 (1996)
        
        [astro-ph/9605094].
        
        Ivanov:1997ia 
        P. Ivanov,
        
        Phys. Rev. D  57, 7145 (1998)
        
        [astro-ph/9708224].
        
        Kawasaki:2006zv 
        M. Kawasaki, T. Takayama, M. Yamaguchi and J. Yokoyama,
        
        Phys. Rev. D  74, 043525 (2006)
        
        [hep-ph/0605271].
        
        Audley:2017drz 
        H. Audley  et al. [LISA Collaboration],
        
        arXiv:1702.00786 [astro-ph.IM].
        
        Hobbs:2009yy 
        G. Hobbs  et al.,
        
        Class. Quant. Grav.  27, 084013 (2010)
        
        [arXiv:0911.5206 [astro-ph.SR]].
        
        Nan:2011um 
        R. Nan  et al.,
        
        Int. J. Mod. Phys. D  20, 989 (2011)
        
        [arXiv:1105.3794 [astro-ph.IM]].
        
        Kramer:2015jsa 
        M. Kramer and B. Stappers,
        
        PoS AASKA  14, 036 (2015).
        
        tomita1971
        K. Tomita, Prog. Theor. Phys.  45, 6 (1971)
        
        Pitrou:2013hga 
        C. Pitrou, X. Roy and O. Umeh,
        
        Class. Quant. Grav.  30, 165002 (2013)
        
        [arXiv:1302.6174 [astro-ph.CO]].
        
        Ando:2017veq 
        K. Ando, K. Inomata, M. Kawasaki, K. Mukaida and T. T. Yanagida,
        
        Phys. Rev. D  97, no. 12, 123512 (2018)
        
        [arXiv:1711.08956 [astro-ph.CO]].
        
        Maggiore:1999vm 
        M. Maggiore,
        
        Phys. Rept.  331, 283 (2000)
        
        [gr-qc/9909001].
        
        Carr:2016drx 
        B. Carr, F. Kuhnel and M. Sandstad,
        
        Phys. Rev. D  94, no. 8, 083504 (2016)
        doi:10.1103/PhysRevD.94.083504
        [arXiv:1607.06077 [astro-ph.CO]].
        
        Musco:2004ak 
        I. Musco, J. C. Miller and L. Rezzolla,
        
        Class. Quant. Grav.  22, 1405 (2005)
        
        [gr-qc/0412063].
        
        Musco:2008hv 
        I. Musco, J. C. Miller and A. G. Polnarev,
        
        Class. Quant. Grav.  26, 235001 (2009)
        
        [arXiv:0811.1452 [gr-qc]].
        
        Musco:2012au 
        I. Musco and J. C. Miller,
        
        Class. Quant. Grav.  30, 145009 (2013)
        
        [arXiv:1201.2379 [gr-qc]].
        
        Harada:2013epa 
        T. Harada, C. M. Yoo and K. Kohri,
        
        Phys. Rev. D  88, no. 8, 084051 (2013)
        Erratum: [Phys. Rev. D  89, no. 2, 029903 (2014)]
        
        [arXiv:1309.4201 [astro-ph.CO]].
        
        Nakama:2016gzw 
        T. Nakama, J. Silk and M. Kamionkowski,
        
        Phys. Rev. D  95, no. 4, 043511 (2017)
        
        [arXiv:1612.06264 [astro-ph.CO]].
        
        Espinosa:2017sgp 
        J. R. Espinosa, D. Racco and A. Riotto,
        
        Phys. Rev. Lett.  120, no. 12, 121301 (2018)
        
        [arXiv:1710.11196 [hep-ph]].
        
        Thrane:2013oya 
        E. Thrane and J. D. Romano,
        
        Phys. Rev. D  88, no. 12, 124032 (2013)
        [arXiv:1310.5300 [astro-ph.IM]].
        
        Kuroda:2015owv 
        K. Kuroda, W. T. Ni and W. P. Pan,
        
        Int. J. Mod. Phys. D  24, no. 14, 1530031 (2015)
        doi:10.1142/S0218271815300311
        [arXiv:1511.00231 [gr-qc]].
        
        Cornish:2018dyw 
        T. Robson, N. J. Cornish and C. Liug,
        
        Class. Quant. Grav.  36, no. 10, 105011 (2019)
        
        [arXiv:1803.01944 [astro-ph.HE]].
        
        Siemens:2013zla 
        X. Siemens, J. Ellis, F. Jenet and J. D. Romano,
        
        Class. Quant. Grav.  30, 224015 (2013)
        [arXiv:1305.3196 [astro-ph.IM]].
        
        Hellings:1983fr 
        R. w. Hellings and G. s. Downs,
        
        Astrophys. J.  265, L39 (1983).
        
        gwsc 
        https://github.com/bingining/GWSC.jl
