Prospective constraints on the primordial black hole abundance from the stochastic gravitational-wave backgrounds produced by coalescing events and curvature perturbations

TL;DR

This paper evaluates how current and future gravitational-wave experiments can constrain the abundance of primordial black holes in dark matter by analyzing stochastic GW backgrounds from binary coalescences and curvature perturbations, focusing on sub-solar mass ranges.

Contribution

It provides the first detailed analysis of constraints on primordial black hole abundance from both coalescence and curvature-induced gravitational waves across a wide mass range.

Findings

Experiments can constrain PBH fraction to between 10^{-5} and 1 from coalescence GWs.

Constraints from curvature perturbations are effective for PBH fractions between 10^{-13} and 1.

Narrow mass range around 10^{-7} M_sun shows limited constraining power for very low PBH fractions.

Abstract

For a variety of on-going and planned gravitational-wave (GW) experiments, we study expected constraints on the fraction ($f_{\rm PBH}$) of primordial black holes (PBHs) in dark matter by evaluating the energy-density spectra of two kinds of stochastic GW backgrounds. The first one is produced from an incoherent superposition of GWs emitted from coalescences of all the binary PBHs. The second one is induced through non-linear mode couplings of large primordial curvature perturbations inevitably associated with the generation of PBHs in the early Universe. In this paper, we focus on the PBHs with their masses of $10^{-8}M_{\odot}\leq M_{\mathrm{PBH}} < 1M_{\odot}$, since they are not expected to be of a stellar origin. In almost all ranges of the masses, we show that the experiments are sensitive to constrain the fraction for $10^{-5} \lesssim f_{\rm PBH} \lesssim 1$ by considering the GWs from coalescing events and $10^{-13} \lesssim f_{\rm PBH} \lesssim 1$ by considering the GWs from curvature perturbations. Exceptionally, only in a narrow range of masses for $M_{\mathrm{PBH}} \simeq 10^{-7} M_{\odot}$, the fraction cannot be constrained for $f_{\rm PBH} \lesssim 10^{-13}$ by those two GW backgrounds.

Equations

$$\frac{k}{k_{\ast }}=7.49\times 10^7{\left(\frac{M_{\odot }}{M_H}\right)}^{1/2}{\left(\frac{g_{\ast ,\rho }(T(M_H))}{106.75}\right)}^{1/4}{\left(\frac{g_{\ast ,s}(T(M_H))}{106.75}\right)}^{-1/3},$$

$$M=K{M}_H{\left(\delta -{\delta }_c\right)}^{\gamma },$$

$${\sigma }^2(k)={\int }_{-\infty }^{+\infty }d\ln q{w}^2(q/k){\left(\frac{4}{9}\right)}^2{\left(\frac{q}{k}\right)}^4{T}^2(q,\tau =1/k)P_{\zeta }(q),$$

$$P_{\zeta }(k)=A\delta (\ln k-\ln k_0),$$

$${\sigma }^2(k)=16A{e}^{-1/x^2}\left[{\cos }^2\left(\frac{1}{\sqrt{3}x}\right)+x\left(3x{\sin }^2\left(\frac{1}{\sqrt{3}x}\right)-\sqrt{3}\sin \left(\frac{2}{\sqrt{3}x}\right)\right)\right],$$

$${\beta }_{M_H}={\int }_{{\delta }_c}^{\infty }\frac{M}{M_H}{\mathcal{P}}_{M_H}(\delta (M))d\delta (M)={\int }_{-\infty }^{\infty }\frac{M}{M_H}{\mathcal{P}}_{M_H}(\delta (M))\frac{d\delta (M)}{d\ln M}d\ln M\equiv {\int }_{-\infty }^{\infty }{\tilde{\beta }}_{M_H}(M)d\ln M,$$

$${\mathcal{P}}_{M_H}(\delta (M))=\frac{1}{\sqrt{2\pi {\sigma }^2(k(M_H))}}\exp \left(-\frac{{\delta }^2(M)}{2{\sigma }^2(k(M_H))}\right),$$

$${\tilde{\beta }}_{M_H}(M)=\frac{K}{\sqrt{2\pi }\gamma \sigma (k(M_H))}{\left(\frac{M}{K{M}_H}\right)}^{1+\frac{1}{\gamma }}\exp \left(-\frac{1}{2{\sigma }^2(k(M_H))}{\left({\delta }_c+{\left(\frac{M}{K{M}_H}\right)}^{\frac{1}{\gamma }}\right)}^2\right).$$

$$f(M)=\frac{{\Omega }_{\mathrm{m}}}{{\Omega }_{\mathrm{CDM}}}{\int }_{-\infty }^{\infty }\left(\frac{g_{\ast ,\rho }(T(M_H))}{g_{\ast ,\rho }(T_{\mathrm{eq}})}\frac{g_{\ast ,s}(T_{\mathrm{eq}})}{g_{\ast ,s}(T(M_H))}\frac{T(M_H)}{T_{\mathrm{eq}}}\right){\tilde{\beta }}_{M_H}(M)d\ln M_H,$$

$${\Omega }_{\text{GW}}(\nu )=\frac{\nu }{{\rho }_{\text{c}}}{\int }_0^{\frac{{\nu }_{\mathrm{cut}}}{\nu }-1}\frac{R_{\text{PBH}}(z)}{(1+z)H(z)}\frac{d{E}_{\text{GW}}}{d{\nu }_{\text{s}}}({\nu }_{\text{s}})dz,$$

$${\Omega }_{\mathrm{IGW}}\left(\nu =\frac{k}{2\pi }\right)$$

$$k=a(M_H)H(M_H).$$

$$\frac{a(M_H)}{a_0(=1)}={\left(\frac{g_{\ast ,s}(T_0)}{g_{\ast ,s}(T(M_H))}\right)}^{1/3}\frac{T_0}{T(M_H)},$$

$$3H^2(M_H)M_{\text{G}}^2=\rho \approx {\rho }_{\text{rad}}=\frac{{\pi }^2{g}_{\ast ,\rho }(T(M_H))}{30}{T}^4(M_H),$$

$$M_H=$$

$$H(M_H)=$$

$$M_H=12{\left(\frac{10}{g_{\ast ,\rho }(T)}\right)}^{1/2}\frac{M_G^3}{T^2}.$$

$$\frac{k}{k_{\ast }}=7.49\times 10^7{\left(\frac{M_{\odot }}{M_H}\right)}^{1/2}{\left(\frac{g_{\ast ,\rho }(T(M_H))}{106.75}\right)}^{1/4}{\left(\frac{g_{\ast ,s}(T(M_H))}{106.75}\right)}^{-1/3},$$

$$\frac{k_{(\mathrm{peak}\mathrm{of}\mathrm{PBHs})}}{k_{(\mathrm{peak}\mathrm{of}\mathrm{primordial}\mathrm{curvature}\mathrm{perturbations})}}=0.730715.$$

$$d{P}_t=\{\begin{array}{l}\frac{3}{58}\left[-{\left(\frac{t}{t_0}\right)}^{\frac{3}{8}}+{\left(\frac{t}{t_0}\right)}^{\frac{3}{37}}\right]\frac{dt}{t}\text{for}t<t_c\\ \frac{3}{58}{\left(\frac{t}{t_0}\right)}^{\frac{3}{8}}\left[-1+{\left(\frac{t}{t_c}\right)}^{-\frac{29}{56}}{\left(\frac{4\pi }{3}{f}_{\mathrm{PBH}}\right)}^{-\frac{29}{8}}\right]\frac{dt}{t}\text{for}t\ge t_c,\end{array}$$

$$R_{\text{PBH}}(z)=\left(\frac{f_{\text{PBH}}{\Omega }_{\text{CDM}}{\rho }_{\mathrm{c}}}{M}\right)\frac{d{P}_t}{dt}.$$

$$\frac{d{E}_{\mathrm{GW}}}{d{\nu }_{\mathrm{s}}}\left({\nu }_{\mathrm{s}}\right)=\frac{{\left(G\pi \right)}^{2/3}{M}_c^{5/3}}{3}\{\begin{array}{l}{\nu }_{\mathrm{s}}^{-1/3}\mathrm{for}{\nu }_{\mathrm{s}}<{\nu }_1\\ w_1{\nu }_{\mathrm{s}}^{2/3}\mathrm{for}{\nu }_1\le {\nu }_{\mathrm{s}}<{\nu }_2\\ w_2\frac{{\sigma }^4{\nu }_{\mathrm{s}}^2}{{\left({\sigma }^2+4{\left({\nu }_{\mathrm{s}}-{\nu }_2\right)}^2\right)}^2}\mathrm{for}{\nu }_2\le {\nu }_{\mathrm{s}}\le {\nu }_3\\ 0\mathrm{for}{\nu }_3<{\nu }_{\mathrm{s}}\end{array}$$

$$\pi M_t{\nu }_1=(1-4.455+3.521)+0.6437\eta -0.05822{\eta }^2-7.092{\eta }^3$$

$$\pi M_t{\nu }_2=(1-0.63)/2+0.1469\eta -0.0249{\eta }^2+2.325{\eta }^3$$

$$\pi M_t\sigma =(1-0.63)/4-0.4098\eta +1.829{\eta }^2-2.87{\eta }^3$$

$$\pi M_t{\nu }_3=0.3236-0.1331\eta -0.2714{\eta }^2+4.922{\eta }^3$$

$${{\Omega }_{\mathrm{GW}}(k)|}_{\mathrm{today}}=\frac{{\Omega }_{\mathrm{r},0}}{24}\left(\frac{g_{\ast ,\rho }(T)}{g_{\ast ,\rho }(T_{\mathrm{eq}})}\right){\left(\frac{g_{\ast ,s}(T)}{g_{\ast ,s}(T_{\mathrm{eq}})}\right)}^{-4/3}{\left(\frac{k}{aH}\right)}^2\overline{{\mathrm{P}}_h(\tau ,k)},$$

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

$$(k\tau {)}^2\overline{I^2(v,u,k\tau \gg 1)}$$

Full text

\acresetall

I Introduction

The first detection of \acpGW from a \acBBH merger by the first \acaLIGO observing run Abbott et al. (2016a) has revived extensive interests in \acpPBH Hawking (1971); Carr and Hawking (1974), which are produced directly from the gravitational collapses of the enhanced inhomogeneities in the primordial Universe. Until now, the origin of these \acpBH and the formation mechanism of \acpBBH are still under debate. Besides an astrophysical origin Abbott et al. (2016b); Belczynski et al. (2016); Coleman Miller (2016), the possibility that these \acpBH are of a primordial origin is also considered Bird et al. (2016); Sasaki et al. (2016); Clesse and Garcia-Bellido (2016a); Eroshenko (2018); Carr et al. (2016); Kashlinsky (2016); Bartolo et al. (2016); Cholis (2016); Harada et al. (2016); Georg and Watson (2017); Nakamura et al. (2016); Raccanelli et al. (2016); Carr et al. (2017); Sasaki et al. (2018). Recently, it has been proposed that the \acpPBH are capable of accounting for the event rate of \acBBH mergers observed by \acaLIGO Bird et al. (2016); Sasaki et al. (2016), although the formation mechanisms of \acPBH binaries bring about uncertainties of a couple of orders of magnitude (see e.g. Ref. Sasaki et al. (2018) and references therein). The \acpPBH can be one of the most promising candidates for the \acCDM Carr et al. (2016). Currently, the nature of \acCDM is still uncertain Bertone and Hooper (2018). There is not definitive evidence for the \acpWIMP which is a prime candidate for \acCDM Tan et al. (2016); Akerib et al. (2016); ATLAS Collaboration and CMS Collaboration (2016); Accardo et al. (2014); Ackermann et al. (2012). Conventionally, one defines the abundance of \acpPBH in \acCDM as a dimensionless fraction of the form $f_{\mathrm{PBH}}=\Omega_{\mathrm{PBH}}/\Omega_{\mathrm{CDM}}$, where $\Omega_{\mathrm{PBH}}$ and $\Omega_{\mathrm{CDM}}$ denote the present energy-density fractions of \acpPBH and \acCDM, respectively. This quantity has been constrained in a variety of mass ranges by a variety of observations (see e.g. Refs. Carr et al. (2010); Sasaki et al. (2018) and references therein), for example, the microlensing events caused by massive astrophysical compact halo objects Calchi Novati et al. (2013); Mediavilla et al. (2009); Green (2016); Chapline and Frampton (2016), the gas accretion effect of \acpPBH on the \acCMB Chen et al. (2016); Ali-Haimoud and Kamionkowski (2017); Poulin et al. (2017), the null detection of a third-order Shapiro time delay using a pulsar timing array Schutz and Liu (2016), and the claimed event rate of \acBBH mergers from \acaLIGO Bird et al. (2016); Sasaki et al. (2016); Abbott et al. (2018), etc.

The \acpPBH can be also a useful probe to the primordial curvature perturbations Mukhanov et al. (1992), since the former are formed via directly gravitational collapses of the latter Hawking (1971); Carr and Hawking (1974). Contrary to the astrophysical processes for which only \acpBH heavier than $\mathcal{O}(1)$ solar mass can be produced Rhoades and Ruffini (1974), the small-mass \acpBH could be also produced by the strong gravity inside the highly compressed overdensities in the early Universe Carr (2005). The \acPBH mass depends on the \acPBH formation redshift $z_{f}$, namely $M\simeq 30M_{\odot}[4\times 10^{11}/(1+z_{f})]^{2}$ Sasaki et al. (2016), where $M_{\odot}$ is the solar mass ($=2\times 10^{33}{\rm g}$). Since inflation models Starobinsky (1979, 1980); Guth (1981); Linde (1982); Albrecht and Steinhardt (1982); Sato (1981); Linde (1983) predict the properties of the primordial curvature perturbations, which determine the mass and abundance of \acpPBH (see e.g. Refs. Kohri et al. (2008); Carr et al. (2010); Sasaki et al. (2018); Byrnes et al. (2018); Yoo et al. (2018) and references therein), our observational knowledge of the \acpPBH is important to learn about the physics of the inflationary Universe.

Recently, it has been proposed that the energy-density fraction of \acpPBH can be constrained by measuring the energy-density spectrum of the \acSGWB. First, the \acSGWB can be produced from an incoherent superposition of \acpGW emitted from all the coalescing \acPBH binaries. The null detection of such a \acSGWB by the first \acaLIGO observing run Abbott et al. (2017) has been used to independently constrain $f_{\mathrm{PBH}}$ Wang et al. (2018); Raidal et al. (2017); Mandic et al. (2016); Clesse and Garcia-Bellido (2016b). For example, Ref. Wang et al. (2018) obtained the tightest observational constraint on $f_{\mathrm{PBH}}$ in the mass range $1-10^{2}M_{\odot}$, pushing the existing observational constraints tighter by one order of magnitude. The possibility to detect the \acSGWB from \acpPBH, in particular from subsolar-mass \acpPBH, by upcoming \acaLIGO observing runs was also predicted Wang et al. (2018). Second, the \acSGWB is induced from the enhanced primordial curvature perturbations Mollerach et al. (2004); Ananda et al. (2007); Baumann et al. (2007); Assadullahi and Wands (2010) 111Based on the inflation model, primordial \acpGW Grishchuk (1975); Starobinsky (1979) are decoupled with primordial curvature perturbations at the first order. However, the induced \acpGW can be generated from primordial curvature perturbations at the second order. Whether or not the primordial \acpGW are detected in the future Creminelli et al. (2015); Huang et al. (2015); Huang and Wang (2019); Cabass et al. (2016); Escudero et al. (2016); Errard et al. (2016); Kamionkowski and Kovetz (2016); Santos et al. (2016); Guzzetti et al. (2016); Lasky et al. (2016), the induced \acpGW could be sizable and even be larger than the primordial \acpGW, if the primordial curvature perturbations are significantly enhanced Ananda et al. (2007); Bugaev and Klimai (2010); Alabidi et al. (2012); Mollerach et al. (2004); Baumann et al. (2007); Assadullahi and Wands (2009); Alabidi et al. (2013); Cai et al. (2019, 2018); Unal (2019).. By making use of the semi-analytic calculation of the induced \acGW spectrum Kohri and Terada (2018); Espinosa et al. (2018), the null detection of such a \acSGWB by a variety of \acGW detectors has been used to obtain constraints on the spectral amplitude of primordial curvature perturbations Inomata et al. (2017); Inomata and Nakama (2018); Orlofsky et al. (2017). The constraints on the induced \acSGWB can be recast as the constraints on the abundance of \acpPBH, and vice versa Bugaev and Klimai (2010, 2011); Saito and Yokoyama (2009, 2010); Nakama and Suyama (2016).

In this paper, we focus on the small-mass \acpPBH for $10^{-8}M_{\odot}\leq M_{\mathrm{PBH}}\leq 1M_{\odot}$. Correspondingly, we calculate the energy-density fraction of the above two kinds of \acpSGWB, and report the expected constraints on the energy-density fraction of \acpPBH from the null detection of the \acpSGWB by several on-going and planned \acGW experiments (see details in Ref. Pitkin et al. (2011)), which include \acSKA Moore et al. (2015), \acLISA Audley et al. (2017); Robson et al. (2019), \acDECIGO Sato et al. (2017) and B-\acDECIGO Isoyama et al. (2018), \acBBO Harry et al. (2006), \acET Punturo et al. (2010), and \acaLIGO design sensitivity Abbott et al. (2016c). Although we focus on the mass range $10^{-8}M_{\odot}\leq M_{\text{PBH}}\leq 1M_{\odot}$, the method of our analysis is equally applicable to the PBS masses outside of this range. In this context, the authors of Refs. Bartolo et al. (2018a, b) studied the SGWB induced by the curvature perturbations associated with the PBHs of masses around $10^{-12}M_{\odot}$ as it may explain the whole abundance of the dark matter. The authors of Ref. Inomata and Nakama (2018) obtained the constraints on the primordial curvature perturbations by studying the detectability of the curvature-induced SGWB in a wide frequency range corresponding to a wide PBH mass range.

First, following Ref. Wang et al. (2018), we evaluate the energy-density fraction of the \acSGWB from binary \acPBH coalescence, by assuming a monochromatic mass distribution of \acpPBH. This choice of the delta function is reasonable since the mass distribution of \acpPBH is insensitive to the details of the spectral shape of primordial curvature perturbations especially after taking into account coarse graining within the Hubble horizon and the effects of critical collapse Byrnes et al. (2018). In addition, the inflation scenario does not favor a significantly extended \acPBH mass distribution Carr et al. (2016). Second, following Ref. Kohri and Terada (2018), we evaluate the energy-density fraction of the induced \acSGWB, by assuming a delta function for the power spectrum of primordial curvature perturbations. In principle, the spectrum of the induced \acSGWB depends on the details of the spectral shape of primordial curvature perturbations. Recently, Ref. Inomata and Nakama (2018) found a spread of the \acSGWB spectrum by studying a log-normal distribution for the power spectrum of primordial curvature perturbations. So the results obtained by this work can be regarded as the conservative one.222This is not that simple because the spectral index of the tails of the SGWB is also relevant as well as the width around the peak, and the spectral index depends on the shape of the curvature perturbations. See the discussion around Eq. (58) of Ref. Espinosa et al. (2018). Anyway, we focus on the delta function case for definiteness. See also Ref. Clesse et al. (2018) which discusses the effects of a broad spectrum.

The rest of the paper is arranged as follows. In Sec. II, we briefly review the formation of \acpPBH in the early Universe, given the power spectrum of primordial curvature perturbations. In Sec. III, we evaluate the energy-density fraction of the \acSGWB from binary \acPBH coalescence, and use it to obtain expected constraints on $f_{\mathrm{PBH}}$ from a variety of on-going and planned \acGW detectors. In Sec. IV, we evaluate the induced \acSGWB from the enhanced primordial curvature perturbations, and obtain the expected constraints on the energy-density fraction of \acpPBH from \acSKA and \acLISA. The conclusions and discussions are given in Sec. V.

II Formation of primordial black holes

Given the power spectrum of the primordial curvature perturbations, we can evaluate a probability of the \acPBH production, the mass function of \acPBHs and the \acPBH abundance Carr et al. (2010); Sasaki et al. (2018). In this work we assume that the \acpPBH are formed in the early Universe which is \acRD. First of all, we need to estimate the wavenumber scale $k$ which is related with a given mass scale $M_{H}$ within the Hubble horizon at the time of horizon re-entry. According to Appendix A, it is represented by

[TABLE]

where $k_{\ast}=0.05\mathrm{Mpc}^{-1}$. Here we can numerically obtain the temperature at the formation $T(M_{H})$ by using Eq. (17). The effective degrees of freedom of relativistic particles, i.e. $g_{\ast\rho}$ and $g_{\ast s}$, are precisely calculated for the Standard Model in Ref. Saikawa and Shirai (2018). Here we interpolate the tabulated data provided by the associated website 333http://member.ipmu.jp/satoshi.shirai/EOS2018 to this reference.

The phenomena of critical collapse Yokoyama (1998); Carr et al. (2016) could describe the formation of \acpPBH with mass $M$ in the early Universe, depending on the horizon mass $M_{H}$ and the amplitude of density fluctuation $\delta$. We have the following relation

[TABLE]

where $K=3.3$, $\gamma=0.36$ and $\delta_{c}={0.45}$ are numerical constants 444Analytically it is estimated to be $\delta_{c}={0.41}$ Harada et al. (2013). In fact, however, they depend on the radial profile of the density perturbations. Harada et al. (2015). The above equation can be inverted to express $\delta$ in terms of $M/M_{H}$, namely, $\delta=(M/(KM_{H}))^{1/\gamma}+\delta_{c}$, which is useful in the following calculations.

In the \acRD Universe, the coarse grained density perturbation is given by

[TABLE]

where $w(q/k)=\mathrm{exp}(-q^{2}/(2k^{2}))$ is a Gaussian window function, and $T(q,\tau)=3(\sin y-y\cos y)/y^{3}$ ($y\equiv q\tau/\sqrt{3}$) is a transfer function (see e.g. Refs. Young et al. (2014); Ando et al. (2018) for details). We consider the power spectrum of primordial curvature perturbations $P_{\zeta}(k)$ to be a delta function of $\ln k$, i.e.,

[TABLE]

where $k_{0}$ is a given constant wavenumber, and $A$ is a dimensionless amplitude. By substituting Eq. (4) into Eq. (3), we obtain

[TABLE]

where $x\equiv k/k_{0}$ is a dimensionless quantity. We show $\sigma^{2}(k)/A$ versus $k/k_{0}$ in a figure at the end of Appendix A.

To convert $\sigma(k)$ to the mass function of \acpPBH, by making use of the Press-Schechter formalism Press and Schechter (1974), we calculate the probability of the \acPBH production, i.e.,

[TABLE]

which accounts for the fraction of the Hubble volumes that collapse into \acpPBH when the horizon mass is $M_{H}$. $\widetilde{\beta}_{M_{H}}(M)$ is the distribution of the (logarithmic) masses of \acpPBH resulting after the critical collapse. Here $\mathcal{P}_{M_{H}}(\delta)$ denotes a Gaussian probability distribution of primordial density perturbations at the given horizon scale corresponding to $M_{H}$. It is represented by

[TABLE]

where $\sigma(k(M_{H}))$ is computed by making use of Eq. (5), and $k(M_{H})$ is given by Eq. (1). The explicit form of $\widetilde{\beta}_{M_{H}}(M)$ is written to be Niemeyer and Jedamzik (1998)

[TABLE]

The mass function of \acpPBH is defined as $f(M)=\frac{1}{\Omega_{\mathrm{CDM}}}\frac{d\Omega_{\mathrm{PBH}}}{d\ln M}$, and the abundance of \acpPBH in \acCDM is given by $f_{\mathrm{PBH}}=\int f(M)\text{d}\ln(M/M_{\odot})$. We obtain the mass function of \acpPBH as follows (see e.g. Ref. Byrnes et al. (2018))

[TABLE]

where $T_{\mathrm{eq}}$ is the temperature of the Universe at the epoch of matter-radiation equality.

In Fig. 1, we depict several examples (colored solid) for the mass function $f(M)$ which is originated from $P_{\zeta}(k)$ in Eq. (4). To be specific, we choose a horizon mass to be $M_{H}=10^{-i}M_{\odot}~{}(i=0,1,2,...,8)$, each of which determines the value of its own $k_{0}=k_{0}(M_{H})$. Here $(10^{2}A)$ is 5.8898, 5.4658, 5.1116, 4.8011, 4.6799, 4.6490, 4.3002, 3.8740 and 3.6813, respectively, so that $f_{\mathrm{PBH}}$ equals the upper limit on $f(M)$ with the aforementioned values of $M_{H}$. For comparison, we plot the existing observational constraint (red dashed) on the \acPBH mass function. The constraint used here arises from the microlensing observations of Subaru/HSC Niikura et al. (2017), OGLE Niikura et al. (2019), EROS-2 Tisserand et al. (2007), MACHO Allsman et al. (2001), and the caustic crossing Oguri et al. (2018).

III Stochastic gravitational-wave background due to binary primordial black hole mergers

Two different mechanisms have been proposed to form binaries from the \acpPBH. One scenario assumes that two \acpPBH could form a binary due to the energy loss via gravitational radiation when they pass by each other accidentally in the late Universe Bird et al. (2016); Clesse and Garcia-Bellido (2016a). The other one assumes that two nearby \acpPBH form a binary due to the tidal force from a third neighboring \acPBH in the early Universe Nakamura et al. (1997); Ioka et al. (1998); Sasaki et al. (2016). Both scenarios are capable of explaining the merger rates of \acpBBH reported by \acaLIGO. However, the first one requires the \acpPBH to contribute most of the \acCDM, which is disfavored by various observational constraints in the relevant mass range. On the other hand, the second one is still allowed. In this work, we thus adopt the formation scenario 555In this section, we use the revised formalism in Ref. Sasaki et al. (2018), instead of the original one in Ref. Sasaki et al. (2016). of \acPBH binaries proposed in Ref. Nakamura et al. (1997) and revisited by Refs. Ioka et al. (1998); Sasaki et al. (2016); Ali-Haimoud et al. (2017); Eroshenko (2018); Hayasaki et al. (2016); Ballesteros et al. (2018); Chen and Huang (2018); Raidal et al. (2018); Liu et al. (2018); Abbott et al. (2018); Magee et al. (2018). In Appendix B we show a brief summary of the formalism for such a scenario.

We calculate the \acSGWB spectrum produced from the coalescing \acPBH binaries. In general, the dimensionless energy-density spectrum of the \acSGWB is defined as $\Omega_{\mathrm{GW}}=\rho_{\mathrm{c}}^{-1}d\rho_{\mathrm{GW}}/d\ln\nu$, where $\rho_{\mathrm{GW}}$ is the \acGW energy density, and $\nu$ is the \acGW frequency Allen and Romano (1999). Knowing the merger rate of \acPBH binaries in Eq. (21), according to Ref. Wang et al. (2018), we can compute the \acSGWB energy-density spectrum within the frequency interval $(\nu,\nu+d\nu)$. It is given by

[TABLE]

where $\frac{dE_{\textrm{GW}}}{d\nu_{\textrm{s}}}(\nu_{\textrm{s}})$ is the \acGW energy spectrum of a \acBBH coalescence (see details in Refs. Ajith et al. (2008, 2011), or a brief summary in Appendix C), $\nu_{\textrm{s}}$ is the frequency in the source frame and is related to the observed frequency $\nu$ through $\nu_{\textrm{s}}=(1+z)\nu$, and $\nu_{\mathrm{cut}}$ is the cutoff frequency for a given \acBBH system.

For the \acPBH binaries for the component masses $10^{-i}M_{\odot}~{}(i=0,1,2,...,8)$, which correspond to the examples of the \acPBH mass function in Fig. 1, we plot the corresponding energy-density fractions of the \acSGWB due to binary \acPBH coalescence at the existing observational constraints on the \acPBH abundance in Fig. 2. The color coding is the same as that in Fig. 1. For comparison, we depict the sensitivity curves 666Sometimes, only the amplitude spectral density $S_{n}(f)$ is shown for a given gravitational-wave detector. We have $\sqrt{S_{n}(f)}=h_{n}(f)f^{-1/2}$, which has the unit of $\mathrm{Hz}^{-1/2}$, and $h_{n}$ is the noise amplitude. The sensitivity to the \acSGWB energy density is related with $S_{n}(f)$ by $\Omega_{\mathrm{GW},n}(f)=3.132\times 10^{35}h^{-2}(f/\mathrm{Hz})^{3}(\sqrt{S_{n}(f)}/\mathrm{Hz}^{-1/2})^{2}$ Moore et al. (2015); Kikuta et al. (2014). The reduced Hubble constant is $h=0.678$ in this paper. of several \acGW experiments (colored dashed/dotted curves), which include pulsar timing array (\acSKA Moore et al. (2015)), space-based \acGW interferometers (\acLISA Robson et al. (2019), \acDECIGO Sato et al. (2017) and B-\acDECIGO Isoyama et al. (2018), \acBBO Harry et al. (2006)), third-generation ground-based \acGW interferometer (\acET Punturo et al. (2010)) and second-generation ground-based \acGW interferometer (\acaLIGO Abbott et al. (2016c)). If the spectrum predicted in a model intersects the sensitivity curve of a given experiment, the expected \acSNR is equal to or greater than unity, which means a possible detection of such a spectrum by this experiment.

Null detection of the \acSGWB by the given future or on-going \acGW experiment can place an upper bound on the magnitude of the energy-density fraction of the \acSGWB at a given frequency band, and can be further recast to constrain the maximum \acPBH abundance. From Fig. 2, all the \acGW experiments have possible contributions to improve the existing observational constraints on the \acPBH abundance, since their sensitivity curves intersect some spectra. Therefore, by regarding the sensitivity curves of all these experiments as upper bounds on the \acSGWB spectrum, we evaluate the expected upper limits on the \acPBH abundance from these experiments. We depict our results in Fig. 3.

Our results are as follows. \acSKA, \acLISA and \acaLIGO will give us relatively weak constraints in future. It is notable that this expected limit from aLIGO is surely stronger than the current one which was reported recently by Ref. Abbott et al. (2018). Both \acET and B-\acDECIGO also have similar constraints on the abundance. All the above four experiments are expected to improve the existing observational constraints on the subsolar-mass \acpPBH. However, both \acDECIGO and \acBBO are expected to significantly improve the existing constraints over the mass range $\mathcal{O}(10^{-6})\leq M/M_{\odot}\leq\mathcal{O}(10^{0})$.

IV Stochastic gravitational-wave background induced by primordial curvature perturbations

The \acSGWB can be also induced by the enhanced primordial curvature perturbations via the scalar-tensor mode coupling in the second-order perturbation theory Ananda et al. (2007). In Appendix D, we show a brief summary of the evaluations of the induced \acSGWB spectrum. For details, see Ref. Kohri and Terada (2018) and references therein. In the following, we will use Appendix D to calculate the energy-density spectrum of the induced \acSGWB, given the form of $P_{\zeta}(k)$ in Eq. (4). We consider the minimal case in which the statistics of the curvature perturbations is Gaussian 777The statistical properties of curvature perturbations can modify the relation between the amount of induced \acSGWB and the \acPBH abundance. Even the curvature perturbations are completely Gaussian, the density contrasts are non-Gaussian due to the nonlinear nature of the gravity, as shown recently by Refs. Yoo et al. (2018); Kawasaki and Nakatsuka (2019); De Luca et al. (2019); Young et al. (2019). and neglect the time evolution of the mass function of PBHs due to accretion, but generalizations can be found in Ref. Garcia-Bellido et al. (2017).

According to Eqs. (27)–(29), we obtain the dimensionless energy-density spectrum of the induced \acSGWB as

[TABLE]

where $\nu=k/2\pi$ denotes the frequency of \acGW, and the dimensionless wavenumber $\tilde{k}=k/k_{0}$ is introduced for simplicity. Based on Appendix A, cosmic temperature $T$ can be numerically related with $M_{H}$ and $k$, and then with $\nu$. Here $\Theta(x)$ denotes the Heaviside theta function with variable $x$.

Similarly to Fig. 2, we plot the energy-density fractions of the induced \acSGWB due to the enhanced primordial curvature perturbations in Fig. 4. Both $A$ and $k_{0}$ are chosen as those in Sec. II. The same color coding is used as in Fig. 1. The double-peak structures arise from the property of the delta function for $P_{\zeta}(k)$ in Eq. (4). For a broader distribution for $P_{\zeta}(k)$, e.g., a log-normal distribution in Ref. Inomata and Nakama (2018), one could find a spread of the \acSGWB spectrum. Therefore, our discussions in the next two paragraphs could be regarded as conservative.

In Fig. 5, besides the energy-density spectra of the induced \acSGWB (colored solid, same as Fig. 4), we depict the sensitivity curves of \acSKA Moore et al. (2015) (red dashed) and \acLISA Robson et al. (2019) (orange dashed) for comparison. For a given spectrum of the induced \acSGWB, we conservatively drop the right-handed peak since such a spiky structure exists only for source spectra with a tiny width. Similarly to the discussions in last section, if a model-predicting spectrum intersects the sensitivity curve of a given experiment, it is possible to measure such a spectrum by this experiment. In such a case, both \acSKA and \acLISA are expected to exclude most of the parameter space, or equivalently improve the existing observational constraints on the \acPBH abundance significantly.

Assuming the null detection of the induced \acSGWB from the enhanced primordial curvature perturbations, similarly to Fig. 3, we plot the expected constraints on the \acPBH abundance from \acSKA (red shaded) and \acLISA (orange shaded) in Fig. 6. The shaded regions mean the excluded parts of the parameter space by these experiments. In fact, here we first obtain the constraints on $A$ from the induced \acSGWB, and then recast them as the upper limits on $f_{\text{PBH}}$ according to the formulae in Sec. II.

Finally, we can combine the results in Fig. 3 with Fig. 6 to obtain Fig. 7. Generally speaking, the slopes of the upper bounds (i.e., boundaries of shaded regions) from the induced \acSGWB are significantly sharper than those (i.e., colored curves) from the \acSGWB due to the coalescing \acPBH binaries. This property can be easily understood as follows. On the one hand, we directly constrained the magnitude of $f_{\mathrm{PBH}}$ by calculating the \acSGWB from the coalescing \acPBH binaries. The detection of such a \acSGWB requires a significant amount of \acPBH binaries in the Universe. This implies that the \acGW detectors can probe the enhanced primordial curvature perturbations only if $A\sim\mathcal{O}(0.1)$. When $A\ll\mathcal{O}(0.1)$, there would be so few \acpPBH in the Universe that the thresholds of \acGW detectors are not triggered. On the other hand, by detecting the induced \acSGWB, we directly obtained the constraints on $A$, which were recast as the indirect constraints on $f_{\mathrm{PBH}}$. In fact, by detecting the induced \acSGWB, the \acGW detectors can probe the primordial curvature perturbations of arbitrary amplitudes within their sensitivities. Since the induced \acSGWB spectrum is proportional to $A^{2}$ while $f_{\mathrm{PBH}}$ is exponentially sensitive to $A$, we obtained the sharper slopes for the upper bounds from the induced \acSGWB than those from the \acSGWB due to coalescing \acPBH binaries in Fig. 7. Thus, the constraints on $f_{\text{PBH}}$ from the \acSGWB induced by the curvature perturbations are stronger than those from the \acSGWB whose origin is merger events except for a narrow gap around $10^{-7}M_{\odot}$ corresponding to the relatively weak observational sensitivity around about $10^{-6}$Hz. Nevertheless, both types of the \acSGWB are complementary and useful for the consistency check of the PBH hypothesis since those two types of the \acSGWB have their own individual features in the spectra and are probed by different observations which are supposed to measure GWs at different frequency bands.

V Conclusions

It has been known that \acpPBH can form binaries in the early Universe, and a \acPBH binary can merge to a new heavier \acBH due to the energy loss via gravitational radiation. Based on Ref. Sasaki et al. (2016), the merger rate of \acPBH binaries depends on the abundance and mass of \acpPBH. Given the existing constraints on the mass function of \acpPBH, following Ref. Wang et al. (2018), we have evaluated the energy-density spectra of \acpSGWB which arise from coalescences of \acPBH binaries with component masses $10^{-i}M_{\odot}~{}(i=0,1,2,...,8)$. From Fig. 2, we found that some of them intersect the sensitivity curves of several future \acGW experiments. This means that the existing limits can be improved by these experiments in the future if these experiments do not detect the \acSGWB. By making use of these sensitivity curves as upper limits on the \acSGWB energy-density fraction, we have evaluated the expected upper limits on the abundance of \acpPBH, and shown our results in Fig. 3. In particular, both DECIGO and BBO are expected to significantly improve the existing limits over the mass range $10^{-6}M_{\odot}-10^{0}M_{\odot}$.

The generation of \acpPBH in the early Universe requires large amplitudes of the primordial curvature perturbations, which can always induce the \acSGWB. By taking into account the existing constraints on the mass function of \acpPBH and making use of the semi-analytic formula in Ref. Kohri and Terada (2018), we have calculated the energy-density spectrum of the induced \acSGWB, and shown our results in Fig. 4. We find several intersections between the induced \acSGWB spectra and the sensitivity curves of \acSKA and \acLISA in Fig. 5. This implies that these experiments can improve the existing upper limits on the mass function of \acpPBH in the future if they claim the null detection of the induced \acSGWB energy-density fraction. In this case, the shaded regions in Fig. 6 will be excluded by \acSKA and \acLISA, respectively.

Finally, by combining Fig. 3 with Fig. 6 to obtain Fig. 7, we found stronger constraints on $f_{\mathrm{PBH}}$ from the \acSGWB induced by curvature perturbations than those from the \acSGWB due to coalescing events, except for a narrow gap around $10^{-7}M_{\odot}$. However, both types of the \acSGWB are complementary and useful for the consistency check of the \acPBH hypothesis.

Acknowledgements.

This work is supported in part by the JSPS Research Fellowship for Young Scientists (TT) and JSPS KAKENHI Grants No. JP17H01131 (SW and KK) and No. JP17J00731 (TT), and MEXT KAKENHI Grants No. JP15H05889 (KK), and No. JP18H04594 (KK).

Appendix A Relation between $k$ and $M_{H}$ in the radiation dominated Universe

During the radiation dominated (\acRD) era of the Universe, the relation between the wavenumber $k$ and the horizon mass $M_{H}$ is obtained as follows. By definition, we have

[TABLE]

The value of the scale factor $a(M_{H})$, when the mode corresponding to $M_{H}$ re-enters the Hubble horizon, is obtained by using the entropy conservation to be

[TABLE]

where $T_{0}=2.725\mathrm{K}$ denotes the present temperature of the \acCMB, and the temperature $T(M_{H})$ is given by the Friedmann equation, i.e.,

[TABLE]

where $M_{\text{G}}=M_{\text{P}}/\sqrt{8\pi}$ is the reduced Planck mass. The relation between the horizon mass $M_{H}$ and the Hubble radius $H^{-1}$ is given by

[TABLE]

Combining Eq. (15) with the left equality of Eq. (14), we have the following formula between $H$ and $M_{H}$, i.e.,

[TABLE]

By combining Eq. (16) with the right equality of Eq. (14), we thus obtain a relation between $M_{H}$ and $T$, i.e.,

[TABLE]

Combining Eqs. (12), (13), (16) and (17) together, we obtain

[TABLE]

where $M_{\odot}$ denotes the solar mass and $k_{*}=0.05\text{Mpc}^{-1}$.

In the following, we depict Fig. 8 to show $\sigma^{2}(k)/A$ versus $k/k_{0}$. In the wavenumber space, the peak of the coarse grained perturbations shifts from the original peak $k_{0}$. Numerically, the shifted peak is obtained as $k=0.730715k_{0}$.

[TABLE]

In our example by assuming $P_{\zeta}(k)$ to be the delta function, the wavenumber in the denominator is nothing but $k_{0}$. When a \acPBH forms, the shorter scales had already experienced the radiation pressure and have been smoothened. Therefore, the \acPBH mass scale corresponds to the coarse grained perturbation scale. In other words, $k$ in Eqs. (1) and (18) should be the one appearing in the numerator of the left-hand side of Eq. (19).

Appendix B Formalism for the merger rate of \acPBH binaries

Given the fraction of \acpPBH in CDM, namely $f_{\mathrm{PBH}}$ 888As was discussed in the Introduction, we use a monochromatic mass distribution of \acpPBH as a reasonable approximation here. Therefore, we use the \acPBH abundance $f_{\mathrm{PBH}}$ instead of the \acPBH mass function $f(M)$. , for a fixed \acPBH mass $M$, the probability that a \acPBH binary coalesces within the cosmic time interval $(t,t+dt)$ is given by (see e.g. Refs. Sasaki et al. (2016, 2018) for details)

[TABLE]

where we define $t_{0}=({3}/{170})\{{\bar{x}}^{4}/[(GM)^{3}\left(4\pi f_{\mathrm{PBH}}/3\right)^{4}]\}$ and $t_{c}=t_{0}\left(4\pi f_{\mathrm{PBH}}/3\right)^{37/3}$, and $\bar{x}=[3M/\left(4\pi\rho_{\mathrm{PBH,eq}}\right)]^{1/3}$ is the physical mean separation of \acpPBH at the epoch of matter-radiation equality. Here $\rho_{\mathrm{PBH,eq}}$ is the energy density of the \acpPBH at the epoch of matter-radiation equality. Multiplying $dP_{t}/dt$ by the present average number density of \acpPBH, one can obtain the merger rate of \acPBH binaries as

[TABLE]

The redshift $z$ is related to the cosmic time $t$ through $t=\int_{z}^{\infty}dz^{\prime}/[(1+z^{\prime})H(z^{\prime})]$, where $H(z)=H_{0}[\Omega_{r,0}(1+z)^{4}+\Omega_{m,0}(1+z)^{3}+\Omega_{\Lambda}]^{1/2}$ is Hubble parameter at redshift $z$. The quantities $\Omega_{r,0}$ and $\Omega_{m,0}$ denote the present energy-density fractions of radiations and non-relativistic matter, respectively. The present energy-density fraction of dark energy is derived as $\Omega_{\Lambda}=1-\Omega_{r,0}-\Omega_{m,0}$. Here $\rho_{\mathrm{c}}=3H_{0}^{2}M_{G}^{2}$ is the critical energy density of the Universe. Throughout this paper, we adopt the $\Lambda$CDM model with cosmological parameters measured by Planck satellite Ade et al. (2016).

Appendix C Energy spectrum of gravitational waves

In the non-spinning limit, the inspiral-merger-ringdown energy spectrum for a \acBBH coalescence takes the following form Ajith et al. (2008, 2011)

[TABLE]

where $\nu_{s}$ is a \acGW frequency in the source frame, $w_{1}$ and $w_{2}$ are two normalization constants that make the spectrum to be continuous. The parameters $\nu_{i}$ ($i=1,2,3$) and $\sigma$ can be expressed in terms of $M_{t}$ and $\eta$ as follows

[TABLE]

which can be found in Table 1 of Ref. Ajith et al. (2011). Here $M_{c}$ is the chirp mass, i.e., $M_{c}^{5/3}=m_{1}m_{2}(m_{1}+m_{2})^{-1/3}$, and $M_{t}=m_{1}+m_{2}$ is the total mass. The symmetric mass ratio is defined by $\eta=m_{1}m_{2}(m_{1}+m_{2})^{-2}$, which gives $0.25$ in this work, since we assume the monochromatic mass of \acpPBH. The cutoff frequency is given to be $\nu_{\mathrm{cut}}=\nu_{3}$.

Appendix D Curvature-induced gravitational waves in a nutshell

We briefly summarize the semi-analytic calculation of the \acSGWB spectrum induced in the \acRD era from the non-linear (tensor-scalar-scalar) mode coupling. The details are described in Ref. Kohri and Terada (2018) and references therein. The energy-density fraction of the induced \acSGWB is given by

[TABLE]

where cosmic temperature $T=T(M_{H}(k))$ with the horizon mass $M_{H}(k)$, $aH$ and conformal time $\tau$ are to be evaluated at (a time somewhat after) the horizon entry of the relevant mode (when the $\Omega_{\mathrm{GW}}$ has reached a temporary asymptotic value). In fact, $T(M_{H}(k))$ can be numerically evaluated by combining Eq. (17) with Eq. (18). The last two factors in the above formula are given by

[TABLE]

In the above equation, we have

[TABLE]

Then, we combine the above three equations to compute the energy-density spectrum of the induced \acSGWB in this work. In case of the delta function source (Eq. (4)), the integral is easily calculated to obtain Eq. (11).