Constraints on primordial curvature perturbations from primordial black hole dark matter and secondary gravitational waves

TL;DR

This paper derives analytical formulas for secondary gravitational waves and uses observational constraints on primordial black holes to set limits on primordial curvature perturbations, aiding early universe physics studies.

Contribution

It introduces a piecewise power law parametrization for the primordial curvature spectrum and links PBH constraints to gravitational wave observations for the first time.

Findings

Derived analytical formulas for gravitational wave energy density evolution.

Established upper bounds on primordial curvature spectrum from PBH constraints.

Discussed future gravitational wave detection prospects for testing the model.

Abstract

Primordial black holes and secondary gravitational waves can be used to probe the small scale physics at very early time. For secondary gravitational waves produced after the horizon reentry, we derive an analytical formula for the time integral of the source and analytical behavior of the time dependence of the energy density of induced gravitational waves is obtained. By proposing a piecewise power law parametrization for the power spectrum of primordial curvature perturbations, we use the observational constraints on primordial black hole dark matter to obtain an upper bound on the power spectrum, and discuss the test of the model with future space based gravitational wave antenna.

Equations

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

$$h_{ij}(\mathbf{x},\eta )=\frac{1}{(2\pi {)}^{3/2}}\int d^{3}k{e}^{i\mathbf{k}\cdot \mathbf{x}}[h_{\mathbf{k}}(\eta )e_{ij}(\mathbf{k})+{\tilde{h}}_{\mathbf{k}}(\eta ){\tilde{e}}_{ij}(\mathbf{k})],$$

$$\begin{array}{rl}{e}_{ij}(\mathbf{k})& =\frac{1}{\sqrt{2}}[e_i(\mathbf{k})e_j(\mathbf{k})-{\tilde{e}}_i(\mathbf{k}){\tilde{e}}_j(\mathbf{k})],\\ {\tilde{e}}_{ij}(\mathbf{k})& =\frac{1}{\sqrt{2}}[e_i(\mathbf{k}){\tilde{e}}_j(\mathbf{k})+{\tilde{e}}_i(\mathbf{k})e_j(\mathbf{k})],\end{array}$$

$${\Phi }_{\mathbf{k}}(\eta )={\varphi }_{\mathbf{k}}\Psi (k\eta ).$$

$$⟨{\varphi }_{\mathbf{k}}{\varphi }_{\tilde{\mathbf{k}}}⟩={\delta }^{(3)}(\mathbf{k}+\tilde{\mathbf{k}})\frac{2{\pi }^2}{k^3}{\left(\frac{3+3w}{5+3w}\right)}^2{\mathcal{P}}_{\zeta }(k),$$

$$\Psi (x)=\frac{9}{x^2}\left(\frac{\sin (x/\sqrt{3})}{x/\sqrt{3}}-\cos (x/\sqrt{3})\right).$$

$$h_{\mathbf{k}}^{\prime \prime }+2\mathcal{H}{h}_{\mathbf{k}}^{\prime }+k^2{h}_{\mathbf{k}}=4S_{\mathbf{k}},$$

$$S_{\mathbf{k}}=\int \frac{d^3\tilde{k}}{(2\pi {)}^{3/2}}{e}_{ij}(\mathbf{k}){\tilde{k}}^i{\tilde{k}}^j\left(2{\Phi }_{\tilde{\mathbf{k}}}{\Phi }_{\mathbf{k}-\tilde{\mathbf{k}}}+\frac{4}{3(1+w){\mathcal{H}}^2}\left({\Phi }_{\tilde{\mathbf{k}}}^{\prime }+\mathcal{H}{\Phi }_{\tilde{\mathbf{k}}}\right)\left({\Phi }_{\mathbf{k}-\tilde{\mathbf{k}}}^{\prime }+\mathcal{H}{\Phi }_{\mathbf{k}-\tilde{\mathbf{k}}}\right)\right).$$

$$⟨h_{\mathbf{k}}(\eta )h_{\tilde{\mathbf{k}}}(\eta )⟩=\frac{2{\pi }^2}{k^3}{\delta }^{(3)}(\mathbf{k}+\tilde{\mathbf{k}}){\mathcal{P}}_h(k,\eta ),$$

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

$$G_{\mathbf{k}}^{\prime \prime }(\eta ,\tilde{\eta })+\left(k^2-\frac{a^{\prime \prime }(\eta )}{a(\eta )}\right)G_{\mathbf{k}}(\eta ,\tilde{\eta })=\delta (\eta -\tilde{\eta }),$$

$$h_{\mathbf{k}}(\eta )=\frac{4}{a(\eta )}{\int }_{{\eta }_k}^{\eta }d\tilde{\eta }{G}_{\mathbf{k}}(\eta ,\tilde{\eta })a(\tilde{\eta })S_{\mathbf{k}}(\tilde{\eta }).$$

$$G_{\mathbf{k}}(\eta ,\tilde{\eta })=\frac{1}{k}\sin [k(\eta -\tilde{\eta })].$$

$${\mathcal{P}}_h(k,\eta )=4{\int }_0^{\infty }dv{\int }_{|1-v|}^{1+v}du{\left[\frac{4v^2-{\left(1-u^2+v^2\right)}^2}{4uv}\right]}^2{I}_{\mathrm{RD}}^2(u,v,x){\mathcal{P}}_{\zeta }(kv){\mathcal{P}}_{\zeta }(ku),$$

$${\Omega }_{\mathrm{GW}}(k,\eta )=\frac{1}{6}{\left(\frac{k}{aH}\right)}^2{\int }_0^{\infty }dv{\int }_{|1-v|}^{1+v}du{\left[\frac{4v^2-{\left(1-u^2+v^2\right)}^2}{4uv}\right]}^2\overline{I_{\mathrm{RD}}^2(u,v,x)}{\mathcal{P}}_{\zeta }(kv){\mathcal{P}}_{\zeta }(ku).$$

$$I_{\text{RD}}(u,v,x)=\frac{1}{9x}\left(I_s\sin x+I_c\cos x\right),$$

$$I_c(u,v,x)=-4{\int }_1^{x}y\sin (y)f(y)dy=T_c(u,v,x)-T_c(u,v,1),$$

$$I_s(u,v,x)=4{\int }_1^{x}y\cos (y)f(y)dy=T_s(u,v,x)-T_s(u,v,1),$$

$$T_c(u,v,x)=-4{\int }_0^{x}y\sin (y)f(u,v,y)dy,$$

$$T_s(u,v,x)=4{\int }_0^{x}y\cos (y)f(u,v,y)dy,$$

$$f(u,v,x)=2\Psi (vx)\Psi (ux)+\left[\Psi (vx)+vx{\Psi }^{\prime }(vx)\right]\left[\Psi (ux)+ux{\Psi }^{\prime }(ux)\right].$$

$$\begin{array}{rl}{T}_c=& \frac{-27}{8u^3{v}^3{x}^4}[-48uv{x}^2(x\cos x+3\sin x)\cos \frac{ux}{\sqrt{3}}\cos \frac{vx}{\sqrt{3}}\\ & +48\sqrt{3}{x}^2\cos x\left(v\cos \frac{vx}{\sqrt{3}}\sin \frac{ux}{\sqrt{3}}+u\cos \frac{ux}{\sqrt{3}}\sin \frac{vx}{\sqrt{3}}\right)\\ & +8\sqrt{3}x\sin x([18-x^2(u^2+3-v^2)]v\cos \frac{vx}{\sqrt{3}}\sin \frac{ux}{\sqrt{3}}\\ & +[18-x^2(v^2+3-u^2)]u\cos \frac{ux}{\sqrt{3}}\sin \frac{vx}{\sqrt{3}})\\ & +24x[-6+x^2(3-u^2-v^2)]\cos x\sin \frac{ux}{\sqrt{3}}\sin \frac{vx}{\sqrt{3}}\\ & +24[-18+x^2(3+u^2+v^2)]\sin x\sin \frac{ux}{\sqrt{3}}\sin \frac{vx}{\sqrt{3}}]\\ & -\frac{27(u^2+v^2-3{)}^2}{4u^3{v}^3}(\text{Si}\left[\left(1-\frac{u-v}{\sqrt{3}}\right)x\right]+\text{Si}\left[\left(1+\frac{u-v}{\sqrt{3}}\right)x\right]\\ & -\text{Si}\left[\left(1-\frac{u+v}{\sqrt{3}}\right)x\right]-\text{Si}\left[\left(1+\frac{u+v}{\sqrt{3}}\right)x\right]),\end{array}$$

$$\begin{array}{rl}{T}_s=& \frac{27}{8u^3{v}^3{x}^4}[48uv{x}^2(x\sin x-3\cos x)\cos \frac{ux}{\sqrt{3}}\cos \frac{vx}{\sqrt{3}}\\ & -48\sqrt{3}{x}^2\sin x\left(v\cos \frac{vx}{\sqrt{3}}\sin \frac{ux}{\sqrt{3}}+u\cos \frac{ux}{\sqrt{3}}\sin \frac{vx}{\sqrt{3}}\right)\\ & +8\sqrt{3}x\cos x([18-x^2(u^2+3-v^2)]v\cos \frac{vx}{\sqrt{3}}\sin \frac{ux}{\sqrt{3}}\\ & +[18-x^2(v^2+3-u^2)]u\cos \frac{ux}{\sqrt{3}}\sin \frac{vx}{\sqrt{3}})\\ & +24x[6-x^2(3-u^2-v^2)]\sin x\sin \frac{ux}{\sqrt{3}}\sin \frac{vx}{\sqrt{3}}\\ & +24[-18+x^2(3+u^2+v^2)]\cos x\sin \frac{ux}{\sqrt{3}}\sin \frac{vx}{\sqrt{3}}]-\frac{27(u^2+v^2-3)}{u^2{v}^2}\\ & +\frac{27(u^2+v^2-3{)}^2}{4u^3{v}^3}(\text{Ci}\left[\left(1-\frac{u-v}{\sqrt{3}}\right)x\right]+\text{Ci}\left[\left(1+\frac{u-v}{\sqrt{3}}\right)x\right]\\ & -\text{Ci}\left[\left|1-\frac{u+v}{\sqrt{3}}\right|x\right]-\text{Ci}\left[\left(1+\frac{u+v}{\sqrt{3}}\right)x\right]+\ln \left|\frac{3-(u+v{)}^2}{3-(u-v{)}^2}\right|).\end{array}$$

$$\text{Si}(x)={\int }_0^{x}dy\frac{\sin y}{y},\text{Ci}(x)=-{\int }_x^{\infty }dy\frac{\cos y}{y}.$$

$$\begin{array}{rl}{I}_{\text{RD}}(u,v,x\to \infty )=& -\frac{3\pi (u^2+v^2-3{)}^2\Theta (u+v-\sqrt{3})}{4u^3{v}^{3}x}\cos x\\ & -\frac{1}{9x}\left(T_c(u,v,1)\cos x+{\tilde{T}}_s(u,v,1)\sin x\right),\end{array}$$

$${\tilde{T}}_s(u,v,1)=T_s(u,v,1)+\frac{27(u^2+v^2-3)}{u^2{v}^2}-\frac{27(u^2+v^2-3{)}^2}{4u^3{v}^3}\ln \left|\frac{3-(u+v{)}^2}{3-(u-v{)}^2}\right|.$$

$$\begin{array}{rl}\overline{I_{\text{RD}}^2(u,v,x\to \infty )}=& \frac{1}{2x^2}[{\left(\frac{3\pi (u^2+v^2-3{)}^2\Theta (u+v-\sqrt{3})}{4u^3{v}^3}+\frac{T_c(u,v,1)}{9}\right)}^2\\ & +{\left(\frac{{\tilde{T}}_s(u,v,1)}{9}\right)}^2].\end{array}$$

$${\Omega }_{\text{GW}}\left(k,{\eta }_0\right)={\Omega }_{\text{GW}}(k,\eta )\frac{{\Omega }_{r0}}{{\Omega }_r(\eta )},$$

$$\Omega (k,\eta )\approx 0.7859A_{\zeta }^2.$$

$${\mathcal{P}}_{\zeta }(k)=A_{\zeta }{\left(\frac{k}{k_p}\right)}^{n_s-1},$$

Full text

1 Introduction

As a result of gravitational collapse, primordial black holes (PBHs) [1, 2, 3, 4] form in a region with its density contrast at horizon reentry during radiation domination exceeding the threshold value. Since the temperature and polarization measurements on cosmic microwave background anisotropy only constrain the primordial perturbations to be very small at large scales, large perturbations at small scales that cause the formation of PBHs [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20] are not constrained and they may produce observable secondary gravitational waves (induced GWs) [21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38]. Therefore, both PBHs and secondary GWs can be used to probe the small scale physics at very early time [39, 40, 41, 42, 43].

PBHs are also dark matter candidate. Observations from extragalactic gamma ray background (EG$\gamma$) [44], femtolensing of gamma-ray bursts [45, 46], millilensing of compact radio sources [47], microlensing of quasars [48], the Milky way and Magellanic Cloud stars [49, 50, 51] constrained the abundance of PBH dark matter [52, 53]. For a recent summary of the constraints, please see Ref. [53]. These constraints can be used to probe the primordial curvature perturbations at small scales. The large curvature perturbations may behave like a narrow peak [5] which can be parameterized as a delta or Gaussian function, or an extended dome shape [10, 28, 54, 37]. Near the inflection point of the potential of the inflaton, ultra-slow-roll inflation may enhance the curvature perturbations by several orders of magnitude at small scales, and the power spectrum for curvature perturbations remains nearly scale invariant [55, 28, 56]. Furthermore, an extremely blue power-law power spectrum may be generated from axionlike curvaton inflationary model at small scales [57]. Therefore, the power spectrum may be parameterized as power-law form at both large and small scales. In this paper, we propose a piecewise power-law parametrization for the power spectrum of primordial curvature perturbations, and use the constraints on the abundance of PBH dark matter to obtain an upper limit on the power spectrum at small scales. With the power spectrum, we calculate the secondary GWs induced by the large density perturbations at small scales. The induced GWs can be tested by space based GW observatory like Laser Interferometer Space Antenna (LISA) [58, 59], TianQin [60] and TaiJi [61], and the Pulsar Timing Array (PTA) [62, 63, 64, 65] including the Square Kilometer Array (SKA) [66] in the future. For simple test, we compare the strength of induced GWs with the sensitivity curves of the space based detectors [67, 68, 69]. On the other hand, the observations of induced GWs can also be used to constrain the power spectrum.

This paper is organized as follows. In section 2, we review the computation of the energy density of induced GWs and derive the formula for the induced GWs produced after the horizon reentry. We propose a piecewise power-law parametrization for the power spectrum of primordial curvature perturbation in section 3, and use current observations on PBH dark matter to obtain an upper bound on the power spectrum. Then we use the formula derived in section 2 and the upper bound to calculate the induced GWs and discuss the possible detection of the induced GWs by future GW observations. The conclusions are drawn in section 4.

2 The induced GWs

Working in the Newtonian gauge, we write the perturbed metric as

[TABLE]

where the scalar perturbation $\Phi$ is the Bardeen potential. The Fourier component of the tensor perturbation $h_{ij}$ is

[TABLE]

where the plus and cross polarization tensors $e_{ij}(\bm{k})$ and $\tilde{e}_{ij}(\bm{k})$ are

[TABLE]

the orthonormal basis vectors $\bm{e}$ and $\tilde{\bm{e}}$ are orthogonal to $\bm{k}$, $\bm{e}\cdot\tilde{\bm{e}}=\bm{e}\cdot\bm{k}=\tilde{\bm{e}}\cdot\bm{k}=0$. The Fourier component of the Bardeen potential $\Phi_{\bm{k}}$ is related with the primordial value $\phi_{\bm{k}}$ by the transfer function $\Psi(k\eta)$

[TABLE]

The primordial value $\phi_{\bm{k}}$ is determined by the primordial curvature perturbation $\mathcal{P}_{\zeta}(k)$ as

[TABLE]

where $w$ is determined by the time when the perturbations reenter the horizon. In this paper, we are interested in those scales that reenter the horizon during radiation domination, so we take $w=1/3$. During radiation domination, the transfer function is

[TABLE]

To the first order, the scalar perturbation decouples from tensor perturbations $h_{ij}$, and the cosmological equation for $h_{ij}$ is homogeneous. But to the second order, they are coupled. The equation for induced GWs with either polarization in Fourier space with $\Phi_{\bm{k}}$ being the source is given by

[TABLE]

where $\mathcal{H}=a^{\prime}/a$ is the conformal Hubble parameter and the prime denotes the derivative with respect to conformal time. The source $S_{\bm{k}}$ is given by

[TABLE]

The power spectrum of the induced GWs is defined as

[TABLE]

and the fractional energy density is

[TABLE]

where the Hubble parameter $H=\mathcal{H}/a$. Before presenting the detailed derivation of the induced GWs, we discuss its qualitative behavior first. Following [23], we assume that the induced GWs are produced instantaneously when the relevant scales reenter the horizon. At the horizon reentry, $h_{\bm{k}}\sim S_{\bm{k}}/k^{2}$ and it gets contributions from all scalar modes $\Phi_{\tilde{\bm{k}}}$. However, combining Eqs. (2.8) and (2.9), it is easy to see that $k^{3}\tilde{k}^{3}/|\bm{k}-\tilde{\bm{k}}|^{3}$ appears in the integrand in $\mathcal{P}_{h}$, so the main contributions to $\mathcal{P}_{h}$ are from $\tilde{\bm{k}}$ that are close to $\bm{k}$. Since the source $S_{\bm{k}}$ decays as $a^{-\gamma}$ with $3\leq\gamma\leq 4$ [23], soon after the horizon reentry GWs propagate freely and $h_{\bm{k}}\propto a^{-1}$, so $\Omega_{\text{GW}}(k,\eta)$ is a constant well within the horizon.

In terms of Green’s function $G_{\bm{k}}(\eta,\tilde{\eta})$ satisfying the equation

[TABLE]

the solution to Eq. (2.7) is

[TABLE]

Note that we assume that the induced GWs are produced after the horizon reentry, so we take $k\eta_{k}=1$. During radiation domination, the Green’s function is

[TABLE]

Combining Eqs. (2.4), (2.6), (2.8), (2.9) and (2.12), after a straightforward and tedious calculation, we obtain the power spectrum of the induced GWs [22, 23, 27, 32]

[TABLE]

where $u=|\bm{k}-\tilde{\bm{k}}|/k$, $v=\tilde{k}/k$, $x=k\eta$, the power spectrum $\mathcal{P}_{\zeta}(k)$ for the primordial curvature perturbation is evaluated at horizon exit during inflation. Combining Eqs. (2.10), (2.14) and (2.27), we get induced GWs in radiation dominated era,

[TABLE]

For the convenience of taking the time average, we split the source term $I_{\text{RD}}$ in the radiation era into the combinations of two oscillations [13],

[TABLE]

where $I_{c}$ and $I_{s}$ are given by 222We learned from Davide Racco that similar results were obtained in [37].

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Note that we take the point of view that induced GWs are produced after the relevant modes reenter the horizon, the lower limit of the integrals (2.17) and (2.18) should be 1 [23, 37], so we need to subtract the terms $T_{c}(u,v,1)$ and $T_{s}(u,v,1)$ in Eqs. (2.17) and (2.18). In [27, 32], the lower limit of the integrals (2.17) and (2.18) was chosen to be zero, i.e., it was assumed that the production of induced GWs begins long before the horizon reentry. We expect the choice of the lower limit of the integral will affect small $k$ modes more because they stay outside the horizon longer. If we take $I_{c}(u,v,x)=T_{c}(u,v,x)$ and $I_{s}(u,v,x)=T_{s}(u,v,x)$, then we recover the result for $I_{\text{RD}}(u,v,x)$ in [32].

Substituting the transfer function (2.6) into Eqs. (2.19) and (2.20), we get

[TABLE]

and

[TABLE]

The sine-integral function Si$(x)$ and cosine-integral function Ci$(x)$ are defined as

[TABLE]

They have the asymptotic behavior $\text{Si}(x)\rightarrow x$ and $\text{Ci}(x)\rightarrow\ln(x)+\gamma_{E}$ as $x\rightarrow 0$, here $\gamma_{E}$ is the Euler number.

At late times, $\eta\gg\eta_{k}$ and $x\rightarrow\infty$,

[TABLE]

where

[TABLE]

So the time average is

[TABLE]

Substituting (2.27) into (2.14), we find that $\overline{\mathcal{P}_{h}(k,\eta)}\sim 1/\eta^{2}$ for the modes well inside the horizon in the radiation dominated era. During radiation domination, $\mathcal{H}=aH\sim 1/\eta$, so $\Omega_{\text{GW}}$ is time independent late in the radiation dominated era as discussed above. Since GWs behave like radiation, the current energy densities of GWs are related to their values well after the horizon reentry in the radiation dominated era

[TABLE]

where $\Omega_{r}$ is the fractional energy density of radiation, $\eta\gg\eta_{k}$ is chosen to be earlier than the matter radiation equality and late enough so that $\Omega_{\text{GW}}(k,\eta)$ is a constant, and the subscript [math] denotes for quantities evaluated at today.

Once we are given the power spectrum $\mathcal{P}_{\zeta}(k)$ for the primordial curvature perturbation, we combine Eqs. (2.15) and (2.27) to calculate induced GWs in radiation dominated era, and obtain $\Omega_{\text{GW}}(k,\eta_{0})$ from Eq. (2.28). In the following, we use several examples to calculate $\Omega_{\text{GW}}$.

2.1 The scale invariant power spectrum

For the scale invariant power spectrum, $\mathcal{P}_{\zeta}(k)=A_{\zeta}$, the numerical integration gives

[TABLE]

Comparing with the result $\Omega(k,\eta)\approx 0.8222A_{\zeta}^{2}$ obtained in [32] by assuming that the production of induced GWs starts well before the horizon reentry, this value is about 4.6% smaller, so the contribution by the induced GWs produced before the horizon reentry is small. In this case, all the modes contribute equally in the integration (2.15) and the contribution to the total integral by $T_{c}(u,v,1)$ and $T_{s}(u,v,1)$ in Eq. (2.16) is small.

2.2 The power-law power spectrum

For a nearly scale invariant power spectrum with the power-law form,

[TABLE]

we get

[TABLE]

where the factor $Q(n_{s})$ needs to be calculated numerically. We show the numerical results for $Q(n_{s})$ in Fig. 1. Again, the results are about 5% smaller than those in [32]. Similar to the scale invariant case, all the modes contribute in the integration (2.15) and the contribution to the total integral by $T_{c}(u,v,1)$ and $T_{s}(u,v,1)$ in Eq. (2.16) is small. In [23], it was estimated that $Q(n_{s})\approx 10$, so that estimate is an order of magnitude larger than the more accurate result $Q(n_{s})\approx 0.8$.

2.3 The monochromatic power spectrum

Now we consider induced GWs generated by the monochromatic curvature perturbations with the delta-function-type power spectrum which can be taken as the idealized limit of a peak in the power spectrum [22, 27, 32]

[TABLE]

where $A_{\zeta}$ is the amplitude and $k_{p}$ is the wave number at which the delta-function peak occurs. Plugging Eq. (2.32) into Eq. (2.15), we get the corresponding induced GWs

[TABLE]

where $\bar{k}\equiv k/k_{p}$,

[TABLE]

and

[TABLE]

The result for Eq. (2.33) is shown as the red solid line in Fig. 2. In the integral (2.15), only the mode $u=v=\bar{k}^{-1}$ contributes to the integration, so $\Omega_{GW}$ is determined by Eq. (2.27) with $u=v=\bar{k}^{-1}$. Around $\bar{k}=2/\sqrt{3}$, the last Ci term in (2.33) is logarithmic divergent, so there is a sharp peak in $\Omega_{\text{GW}}$ at $\bar{k}=2/\sqrt{3}$ as shown in Fig. 2. Before the sharp peak, some terms in (2.33) may cancel each other, so there is a dip in $\Omega_{\text{GW}}$ before the sharp peak. The sharp peak at $k=2k_{p}/\sqrt{3}$ is due to the resonant amplification [22, 32]. The factor 2 is from the second order effect (the source is $\Phi^{2}$) and the factor $1/\sqrt{3}$ is due to the sound speed of radiation background. For comparison, the blue solid line in Fig. 2 shows the induced GWs from the monochromatic power spectrum obtained in [32]. The difference comes from the generation of GWs before the horizon reentry. In [32], they assume that the production of induced GWs starts when the Universe begins radiation domination. In deriving Eq. (2.33), we assume that induced GWs are produced after the relevant scales reenter the horizon. For the monochromatic power spectrum, only single mode $u=v=\bar{k}^{-1}$ contributes to the integration, and large scales reenter the horizon at later time, so the contribution by the production of induced GWs well outside the horizon will be larger. As shown in Fig. 2, the difference can be several orders of magnitude at large scales ($k/k_{p}\lesssim 0.1$), but the difference at small scales ($k/k_{p}\gtrsim 1$) is small.

2.4 The Gaussian power spectrum

More generally, the peak in the power spectrum may be parameterized as Gaussian form [29, 35]

[TABLE]

In the limit $\sigma\rightarrow 0$, the sharp peak with delta function (2.32) is recovered. The smaller width $\sigma$ corresponds to a sharper peak on power spectrum. Combining Eqs. (2.15), (2.27) and (2.36), we calculate induced GWs produced by Gaussian power spectrum with $\sigma=0.2$ and $\sigma=0.5$, and the results are shown as the dashed green and black lines in Fig. 2. Since more modes contribute to the integration (2.15), the resonance peak is broadened and the amplitude is reduced as $\sigma$ becomes larger. The dip before the peak disappears when $\sigma$ is large enough. We also show the induced GWs by considering the production of GWs before the horizon reentry with the black dot dashed line in Fig. 2. The difference due to the production of GWs before the horizon reentry is large at large scales.

3 PBH and the observational constraints

PBHs form in the region with its density contrast at horizon reentry exceeding the threshold $\delta_{c}$. Suppose the density perturbations are Gaussian, the probability distribution of the smoothed density contrast $\delta(R)$ over a sphere with comoving radius $R$ is [70]

[TABLE]

where the smoothing scale $R$ is the horizon size, $R=\mathcal{H}^{-1}$ and the mass variance $\sigma(R)$ associated with the PBH mass $M_{\text{PBH}}$ is

[TABLE]

$\mathcal{P}_{\delta}$ is the power spectrum of the matter perturbation and the window function is ${W}(kR)=\exp\left(-k^{2}R^{2}/2\right)$. During radiation domination, the matter perturbation relates to the primordial curvature perturbation as

[TABLE]

Using Press-Schechter theory [71], we get the fraction of the energy density in the Universe going to PBHs 333There should be a factor $\gamma$ in (3.4) [26]. However, it has very little effect on the result, so we ignore this factor here.

[TABLE]

where $\delta_{c}=0.42$ [72]. Combining Eqs. (3.2) and (3.3), we see that the dominant contribution to the mass variance $\sigma^{2}(R)$ comes from the scale $k=1/R$, so $\sigma^{2}(R)\propto\mathcal{P}_{\zeta}(1/R)$. Following Ref. [53], at each $k$, we calculate $\sigma^{2}(R)$ with scale invariant $\mathcal{P}_{\zeta}$, so we have

[TABLE]

Since PBH forms in the radiation dominated era, the mass of PBH is of the order of the horizon mass $M_{H}=4\pi\rho/(3H^{3})=(2GH)^{-1}$ [42]

[TABLE]

where the order one ratio $\gamma$ is chosen as $\gamma=3^{-3/2}\approx 0.2$ [3], $\Omega_{\text{r}0}=9.17\times 10^{-5}$, $M_{0}=(2GH_{0})^{-1}\approx 4.63\times 10^{22}~{}M_{\odot}$, $H_{0}=67.27\text{ km/s/Mpc}$ [73], $g_{*}^{0}\approx 3.36$ and $g_{*}^{i}$ denote the effective degrees of freedom for energy density at present and at the formation of PBH respectively. In this paper, we don’t distinguish the difference between the effective degrees of freedom for the entropy and energy density. For the mass scale of PBHs we are interested in, we take $g_{*}^{i}\approx 10.75$. After their formation, PBHs behave like matter, so the energy fraction of PBHs increases until the matter radiation equality. Ignoring the mass accretion and evaporation, the energy fraction of PBHs at their formation is

[TABLE]

where $f_{\text{PBH}}=\Omega_{\text{PBH}}/\Omega_{\text{DM}}$ is the current energy fraction of PBHs $\Omega_{\text{PBH}}$ to dark matter $\Omega_{\text{DM}}$.

Combining Eqs. (3.5) and (3.7), we can obtain $\mathcal{P}_{\zeta}$ for a given $f_{\text{PBH}}$ and vice versa. This allows us to use the observational constraints on PBH abundance, namely $f_{\text{PBH}}$, to constrain the power spectrum for primordial curvature perturbations at small scales. Alternatively, it allows us to use $f_{\text{PBH}}$ to constrain some inflationary models. The current observational constraints on $f_{\text{PBH}}$ and $\mathcal{P}_{\zeta}$ at small scales were summarized in Ref. [53] and we show them in Fig. 3. On observable scales $10^{-4}\text{ Mpc}^{-1}\lesssim k\lesssim 10^{-1}\text{ Mpc}^{-1}$, the temperature and polarization measurements on the cosmic microwave background anisotropy constrain the nearly scale invariant power spectrum for the primordial curvature perturbation as [74]

[TABLE]

where $k_{*}=0.05\text{ Mpc}^{-1}$, $A_{s}=2.1\times 10^{-9}$ and $n_{s}=0.9649\pm 0.0044$.

To get large enhancement on the power spectrum from a single field inflation, an ultra-slow-roll inflation near the inflection point may be used, and the power spectrum is also nearly scale invariant [55, 28, 56]. Furthermore, an extremely blue power-law power spectrum may be generated from axionlike curvaton model at small scales [57]. For simplicity, we use the piecewise power-law parametrization for the power spectrum to fit the observational bounds in Fig. 3 and the fitting results are

[TABLE]

We show this piecewise power-law parametrization of the power spectrum in Fig. 3 by the solid black line. The power spectrum well fits the upper bound from PBHs and the CMB constraints.

Due to the uncertainties in the value of $\delta_{c}$ and the effect of non spherical collapse, the upper limit on the power spectrum by the non detection of PBH dark matter can be much different [75, 53]. However, the method discussed here can be easily applied to those cases. Using the power spectrum (3.9) and the method of calculating induced GWs presented in the previous section, we obtain the energy density of secondary GWs and the result is shown in Fig. 4. In Fig. 4, we also plot the sensitivity curves for the ground based detector advanced Laser Interferometer Gravitational-Wave Observatory (aLIGO) [76, 77], future space based GW detectors LISA [58, 59] and TianQin [60], and PTA [62, 63, 64, 65] including the European PTA (EPTA) and SKA [66]. It is obvious that the secondary GWs can be detected by EPTA, SKA, LISA and TianQin although there is no detection of PBH dark matter. In other words, the observation of induced GWs puts stronger constraint on the primordial curvature perturbation at small scales. Since the current PTA observations don’t find stochastic GWs yet, so the upper limit (3.9) for $k\gtrsim 10^{4}\text{ Mpc}^{-1}$ is overestimated. Using the power-law power spectrum (3.9), we calculate the $\mu$ distortion [78, 79]

[TABLE]

where

[TABLE]

$k_{\text{min}}\approx 1~{}$Mpc*-1*, $A\approx 0.9$ and $\hat{k}=k/[1\ \text{Mpc}^{-1}]$, and we get $\mu_{\text{ac}}=0.03$. Again this result shows that the upper limit (3.9) for $k\gtrsim 10^{4}\text{ Mpc}^{-1}$ is too large. In other words, the upper limit set by PTA and $\mu$ distortion is more stringent.

For the power-law power spectrum, if there is no detection of induced GWs by LISA, then the constraint is

[TABLE]

If we choose $\delta_{c}=0.42$, plugging the constraint (3.12) into Eqs. (3.5) and (3.7), we get $f_{\text{PBH}}<10^{-400}$. This means if LISA does not observe induced GWs, then the contribution from PBHs with the mass around $10^{-14}M_{\odot}$ to dark matter is negligible. In Fig. 4, we also show the secondary GWs produced by the inflationary model with the polynomial potential [28]. For convenience, we call the model as D-G model. From Fig. 4, we find that the D-G model can be tested by SKA, LISA and TianQin in the future.

4 Conclusion

In the case that the production of secondary GWs starts long before the horizon reentry, there was an analytical formula for the time integral of the source $I_{\text{RD}}(u,v,\eta)$. For secondary GWs produced only after the horizon reentry, we derive similar analytical formula for $I_{\text{RD}}(u,v,\eta)$ by splitting $I_{\text{RD}}(u,v,\eta)$ into the combinations of two oscillations $\sin(k\eta)$ and $\cos(k\eta)$. With this analytical formula, it is easy to obtain the $1/\eta^{2}$ behavior of the power spectrum of induced GWs and hence it helps to understand why induced GWs evolve as radiation at late time. For nearly scale invariant primordial curvature perturbations, we find that the GWs produced before the horizon reentry contribute about 5% to the total energy density of induced GWs because all the modes accounts for the production. For the power spectrum of curvature perturbations with a sharp peak which is parameterized as a delta function, there exists resonant amplification because only one single mode contributes to the integration, and the production of GWs before the horizon reentry becomes dominant at large scales. Since the amplitude of induces GWs is proportional to the square of the peak amplitude of the power spectrum of curvature perturbations, $A_{\zeta}^{2}$, the upper bound on $A_{\zeta}$ from the observational constraints on PBH dark matter can be used to discuss the possible detection of induced GWs.

Using the piecewise power-law parametrization for the power spectrum of primordial curvature perturbations and the observational constraints on PBH dark matter, the best fit upper bound on primordial curvature perturbations was obtained. We find that at small scales $k\gtrsim 10^{4}$ Mpc*-1*, the upper limit on the power spectrum is $\mathcal{P}_{\zeta}\lesssim 0.05$. However, this upper limit gives large stochastic GW background which is inconsistent with the observations of EPTA and the $\mu$ distortion caused by this upper limit is also too large. This means that in nanohertz bands, PTA observations set more stringent upper bound. On the other hand, if the power spectrum peaks at some particular small scales, it evades the constraint by EPTA. Therefore, the detection of induced GWs in the future puts more stringent constraint on the power spectrum and the abundance of PBH dark matter. The non-detection of induced GWs by LISA constrains the power spectrum in the LISA band to be $\mathcal{P}_{\zeta}\lesssim 4\times 10^{-4}$, so the contribution from PBHs with the mass around $10^{-14}M_{\odot}$ to dark matter is negligible if induced GWs are not observed by LISA in the future.

5 acknowledgments

This research was supported in part by the National Natural Science Foundation of China under Grant No. 11875136 and the Major Program of the National Natural Science Foundation of China under Grant No. 11690021. The authors would like to thank Davide Racco for pointing out their work which derived similar formulae as ours in Section 2.