Merger Rate Distribution of Primordial Black Hole Binaries

It is believed that primordial black holes (PBHs) could have formed in the early Universe from the collapse of large density fluctuations (Hawking 1971; Carr & Hawking 1974; Carr 1975). On the other hand, one of the challenges for fundamental physics is the understanding of the nature of dark matter (DM). Among a large variety of models, the speculation that the DM is composed totally or partially by PBHs has attracted much attention, especially since the discovery of black hole coalescence by LIGO (Abbott et al. 2016a), because the coalescence of PBH binaries can be candidates for the observed gravitational-wave events (Bird et al. 2016; Sasaki et al. 2016).

Up to now, there are several gravitational-wave events from binary black hole (BBH) mergers reported by LIGO and VIRGO collaborations: GW150914 (${36}_{-4}^{+5}\,{M}_{\odot }$, ${29}_{-4}^{+4}\,{M}_{\odot }$) (Abbott et al. 2016a), GW151226 (${14.2}_{-3.7}^{+8.3}\,{M}_{\odot }$, ${7.5}_{-2.3}^{+2.3}\,{M}_{\odot }$) (Abbott et al. 2016d), GW170104 (${31.2}_{-6.0}^{+8.4}\,{M}_{\odot }$, ${19.4}_{-5.9}^{+5.3}\,{M}_{\odot }$) (Abbott et al. 2017d), GW170608 (${12}_{-2}^{+7}\,{M}_{\odot }$, ${7}_{-2}^{+2}\,{M}_{\odot }$) (Abbott et al. 2017a), GW170814 (${30.5}_{-3.0}^{+5.7}\,{M}_{\odot }$, ${25.3}_{-4.2}^{+2.8}\,{M}_{\odot }$) (Abbott et al. 2017c), as well as a less significant candidate, LVT151012 (${23}_{-6}^{+18}\,{M}_{\odot }$, ${13}_{-5}^{+4}\,{M}_{\odot }$) (Abbott et al. 2016e, 2016f). These events indicate that the black holes should have an extended mass function. In fact, the generic initial conditions of PBH formation also suggest that the PBH mass should also extend over a wide range.

In the literature there are two main paths for the PBH binary formation. One is formed in the early Universe (Nakamura et al. 1997; Sasaki et al. 2016; Ali-Hamoud et al. 2017) and another is formed in the late Universe (Bird et al. 2016; Ali-Hamoud et al. 2017; Nishikawa et al. 2017), respectively, and the former generically makes the dominant contribution to the PBH merger rate. However, the mass function of PBHs is usually assumed to be monochromatic in Sasaki et al. (2016), Nakamura et al. (1997), Ali-Hamoud et al. (2017), Bird et al. (2016), and Nishikawa et al. (2017). Recently, the merger rate of PBHs with an extended mass function is investigated in Raidal et al. (2017) and Kocsis et al. (2018), but only the tidal force from the PBH closest to the center of mass of the PBH binary is taken into account in Raidal et al. (2017) and a flat mass function of PBHs over a small mass range is considered in Kocsis et al. (2018).

In this paper we consider the torques of all PBHs and linear density perturbations, and calculate the merger rate distribution for the PBH binaries with a general mass function. In the upcoming decades, many more BBH mergers are expected to be detected, and they will provide better information about the black hole mass function; this may finally help us answer the question: What is the origin of black holes, and how are black hole binaries formed in our Universe?

The probability distribution function (PDF) of PBHs, P(m), is normalized to be

$$\begin{eqnarray}&&{\int }_{0}^{\infty }P(m){dm}=1,\end{eqnarray} \tag{ 1 }$$

and the abundance of PBHs in the mass interval, $(m,m+{dm})$, is given by

$$\begin{eqnarray}&&{fP}(m){dm},\end{eqnarray} \tag{ 2 }$$

where f is the total abundance of PBHs in nonrelativistic matter. In this paper, for convenience, the PBH mass is in units of ${M}_{\odot }$. The fraction of PBHs in cold DM is related to f by ${f}_{\mathrm{pbh}}\equiv {{\rm{\Omega }}}_{\mathrm{pbh}}/{{\rm{\Omega }}}_{\mathrm{cdm}}\approx f/0.85$. We introduce a cross-grained discrete PDF, namely

$$\begin{eqnarray}&&\int P(m){dm}=1\ \to \ \displaystyle \sum _{{m}_{\min }\leqslant {m}_{i}\leqslant {m}_{\max }}{P}_{i}{\rm{\Delta }}\simeq 1,\end{eqnarray} \tag{ 3 }$$

where $P({m}_{i})\to {P}_{i}$ is the binned PDF and ${{dm}}_{i}\to {\rm{\Delta }}$ denotes the resolution of the PBH mass. Roughly speaking, ${{fP}}_{i}{\rm{\Delta }}\equiv {f}_{i}{\rm{\Delta }}$ is taken as the abundance of PBHs with mass m_i. At matter-radiation equality the total energy density of matter is

$$\begin{eqnarray}&&{\rho }_{\mathrm{eq}}={{\rm{\Omega }}}_{m}{\rho }_{\mathrm{crit}}{(1+{z}_{\mathrm{eq}})}^{3},\end{eqnarray} \tag{ 4 }$$

and then the average distance, ${\bar{x}}_{i}$, between two PBHs with mass m_i is

$$\begin{eqnarray}&&{\bar{x}}_{i}={\left(\displaystyle \frac{3}{4\pi }\displaystyle \frac{{m}_{i}}{{\rho }_{\mathrm{eq}}{f}_{i}{\rm{\Delta }}}\right)}^{1/3},\end{eqnarray} \tag{ 5 }$$

which depends on both mass m_i and its abundance ${f}_{i}{\rm{\Delta }}$. We need to stress that the number densities for PBHs with different masses can be quite different from each other, and in principle, there is no well-defined number density for the overall PBHs. The average distance, $\langle {x}_{{ij}}\rangle$, between two neighboring PBHs with different masses m_i and m_j is estimated as

$$\begin{eqnarray}&&\langle {x}_{{ij}}\rangle ={({\bar{x}}_{i}^{-3}+{\bar{x}}_{j}^{-3})}^{-1/3}={\mu }_{{ij}}^{1/3}{\bar{x}}_{{ij}},\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&{\mu }_{{ij}}=\displaystyle \frac{2{m}_{i}{m}_{j}{f}_{b}}{{m}_{b}({f}_{j}{m}_{i}+{f}_{i}{m}_{j})},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\bar{x}}_{{ij}}^{3}=\displaystyle \frac{3}{8\pi }\displaystyle \frac{{m}_{b}}{{\rho }_{\mathrm{eq}}{f}_{b}{\rm{\Delta }}},\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&{f}_{b}={f}_{i}+{f}_{j},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{m}_{b}={m}_{i}+{m}_{j}.\end{eqnarray} \tag{ 10 }$$

The above formula are valid for ${m}_{i}\ne {m}_{j}$, and can be generalized to cover ${m}_{i}={m}_{j}=m$ if we take $P({m}_{i})\,=P({m}_{j})=P(m)/2$. From now on, for simplicity, we omit the subscript _ij unless it is necessary.

In order for the formation of a PBH binary, the two neighboring PBHs necessarily decouple from the background expansion and form a bound system. In Newtonian approximation, the equation governing the evolution of the proper separation, r, of the black hole binary with masses m_i and m_j along the axis of motion takes the form

$$\begin{eqnarray}&&\ddot{r}-(\dot{H}+{H}^{2})r+\displaystyle \frac{{m}_{b}}{{r}^{2}}\displaystyle \frac{r}{| r| }=0,\end{eqnarray} \tag{ 11 }$$

where the dot denotes the derivative with respect to the proper time. In this paper we work in the geometric units $G=c=1$. Defining $\chi \equiv r/x$, we rewrite Equation (11) as

$$\begin{eqnarray}&&\chi ^{\prime\prime} +\displaystyle \frac{{sh}^{\prime} +h}{{s}^{2}h}(s\chi ^{\prime} -\chi )+\displaystyle \frac{1}{\lambda }\displaystyle \frac{1}{{({sh})}^{2}}\displaystyle \frac{1}{{\chi }^{2}}\displaystyle \frac{\chi }{| \chi | }=0,\end{eqnarray} \tag{ 12 }$$

where x is the comoving separation between these two PBHs, primes denote the derivative with respect to scale factor s, which is normalized to be unity at equality, and $h(s)\equiv H(s)/{\left(\tfrac{8\pi }{3}{\rho }_{\mathrm{eq}}\right)}^{1/2}=\sqrt{{s}^{-3}+{s}^{-4}}$. Here the dimensionless parameter $\lambda$ is

$$\begin{eqnarray}&&\lambda =\displaystyle \frac{8\pi {\rho }_{\mathrm{eq}}{x}^{3}}{3{m}_{b}}=\displaystyle \frac{X}{{f}_{b}{\rm{\Delta }}},\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}&&X\equiv {x}^{3}/{\bar{x}}^{3},\end{eqnarray} \tag{ 14 }$$

and $\bar{x}$ is given in Equation (8). The solution of Equation (12) in Ali-Hamoud et al. (2017) implies that the decoupling occurs before the equality if $\lambda \lt 1$ and the semimajor axis a of the formed binary is given by

$$\begin{eqnarray}&&a\approx 0.1\lambda x=\displaystyle \frac{0.1}{{f}_{b}{\rm{\Delta }}}\displaystyle \frac{{x}^{4}}{{\bar{x}}^{3}}=\displaystyle \frac{0.1\bar{x}}{{f}_{b}{\rm{\Delta }}}{X}^{\tfrac{4}{3}}.\end{eqnarray} \tag{ 15 }$$

Without the tidal force from other PBHs and density perturbations, these two PBHs will just collide head-on with each other. However, the tidal force will provide an angular momentum to prevent this system from direct coalescence. For simplicity, we introduce a dimensionless angular momentum, j, defined by

$$\begin{eqnarray}&&j\equiv {\ell }/\sqrt{{m}_{b}a}=\sqrt{1-{e}^{2}},\end{eqnarray} \tag{ 16 }$$

where ℓ is the angular momentum per unit reduced mass and $e\in [0,1]$ is the eccentricity. In order to estimate the initial orbital parameters for the PBH binary in which two PBHs can have different masses, we generalize the method in Ali-Hamoud et al. (2017). Here we do not plan to repeat all of the calculations in Ali-Hamoud et al. (2017), but only highlight the key results. The local tidal field is related to the Newtonian potential ϕ by ${T}_{{ij}}=-{\partial }_{i}{\partial }_{j}\phi$ which exerts a perturbative force per unit mass ${\boldsymbol{F}}={\boldsymbol{T}}\cdot {\boldsymbol{r}}$. Supposing that the initial comoving separation of the binary is small relative to the the mean separation, this tidal force does not significantly affect the orbit of the binary, but it produces a torque which yields

$$\begin{eqnarray}&&{\boldsymbol{\ell }}=\int {dt}\ {\boldsymbol{r}}\times [{\boldsymbol{T}}\cdot {\boldsymbol{r}}].\end{eqnarray} \tag{ 17 }$$

Since the tidal field generated by other PBHs and density perturbation in the radiation-domination era is ${s}^{-3}$, ${\boldsymbol{T}}\simeq {s}^{-3}{{\boldsymbol{T}}}_{\mathrm{eq}}$ and then

$$\begin{eqnarray}&&{\boldsymbol{j}}\approx {x}^{3}\hat{x}\times \left[\displaystyle \frac{{{\boldsymbol{T}}}_{\mathrm{eq}}}{{m}_{b}}\cdot \hat{x}\right],\end{eqnarray} \tag{ 18 }$$

where $\hat{x}$ is the unit vector along ${\boldsymbol{x}}$ and ${{\boldsymbol{T}}}_{\mathrm{eq}}$ is the local tidal field at equality. The tidal field generated by a PBH with mass m_l at a comoving separation $y\gg x$ is given by

$$\begin{eqnarray}&&{T}_{\mathrm{eq}}^{{ij}}={m}_{l}\displaystyle \frac{3{\hat{y}}^{i}{\hat{y}}^{j}-{\delta }^{{ij}}}{{y}^{3}},\end{eqnarray} \tag{ 19 }$$

and then

$$\begin{eqnarray}&&{\boldsymbol{j}}\approx 3\displaystyle \frac{{m}_{l}}{{m}_{b}}\displaystyle \frac{{x}^{3}}{{y}^{3}}(\hat{x}\cdot \hat{y})(\hat{x}\times \hat{y}).\end{eqnarray} \tag{ 20 }$$

Similar to Ali-Hamoud et al. (2017), following Chandrasekhar (1943), the two-dimensional PDF of j from torques by all other PBHs is given by

$$\begin{eqnarray}&&\displaystyle \frac{{dP}}{{d}^{2}j}=\mathop{\mathrm{lim}}\limits_{V\to \infty }\int \displaystyle \frac{{d}^{2}k}{{(2\pi )}^{2}}{e}^{i{\boldsymbol{k}}\cdot {\boldsymbol{j}}}\displaystyle \prod _{l}{{ \mathcal I }}_{l}^{{N}_{l}},\end{eqnarray} \tag{ 21 }$$

where ${N}_{l}={n}_{l}V$ is total number of PBHs with mass m_l,

$$\begin{eqnarray}&&{{ \mathcal I }}_{l}={\int }_{V}\displaystyle \frac{{d}^{3}y}{V}\exp \left[-3\displaystyle \frac{{m}_{l}}{{m}_{b}}i\displaystyle \frac{{x}^{3}}{{y}^{5}}{y}_{| | }\,{\boldsymbol{k}}\cdot {{\boldsymbol{y}}}_{\perp }\right],\end{eqnarray} \tag{ 22 }$$

and ${y}_{| | }\equiv {\boldsymbol{y}}\cdot \hat{x}$ and ${{\boldsymbol{y}}}_{\perp }\equiv \hat{x}\times {\boldsymbol{y}}$. After a tedious computation, we find

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{V\to \infty }{{ \mathcal I }}_{l}^{{N}_{l}}={e}^{-\displaystyle \frac{4\pi }{3}\displaystyle \frac{{m}_{l}}{{m}_{b}}{n}_{l}{x}^{3}k}.\end{eqnarray} \tag{ 23 }$$

Since ${m}_{l}{n}_{l}={\rho }_{l}$ is the energy density of PBHs with mass m_l, ${\sum }_{l}{\rho }_{l}={\rho }_{\mathrm{pbh}}=f{\rho }_{\mathrm{eq}}$ and then

$$\begin{eqnarray}&&\displaystyle \frac{{dP}}{{dj}}=j\int {{kdkJ}}_{0}({kj}){e}^{-{j}_{X}k},\end{eqnarray} \tag{ 24 }$$

where

$$\begin{eqnarray}&&{j}_{X}=0.5\displaystyle \frac{f}{{f}_{b}{\rm{\Delta }}}X,\end{eqnarray} \tag{ 25 }$$

which encodes the torques by all other PBHs. Integrating over Equation (24) gives

$$\begin{eqnarray}&&{\left.j\displaystyle \frac{{dP}}{{dj}}\right|}_{X}={ \mathcal P }(j/{j}_{X}),\quad { \mathcal P }(\gamma )=\displaystyle \frac{{\gamma }^{2}}{{(1+{\gamma }^{2})}^{3/2}},\end{eqnarray} \tag{ 26 }$$

where $\gamma =j/{j}_{X}$. In addition, the variance of ${\boldsymbol{j}}$ due to the torques by density perturbations is

$$\begin{eqnarray}&&\langle {j}^{2}{\rangle }^{1/2}\approx 0.5\displaystyle \frac{8\pi }{3}\displaystyle \frac{{\sigma }_{\mathrm{eq}}{\rho }_{\mathrm{eq}}}{{m}_{b}}{x}^{3}\,=\,0.5\displaystyle \frac{{\sigma }_{\mathrm{eq}}}{{f}_{b}{\rm{\Delta }}}X,\end{eqnarray} \tag{ 27 }$$

where ${\sigma }_{\mathrm{eq}}\equiv \langle {\delta }_{\mathrm{eq}}^{2}{\rangle }^{1/2}$ is the variance of density perturbations of the rest of DM on a scale of the order of ${ \mathcal O }({10}^{0}\sim {10}^{3}){M}_{\odot }$ at equality. Taking into account both the torques by all of other PBHs and density perturbations, the characteristic value of j_X in Equation (26) reads as

$$\begin{eqnarray}&&{j}_{X}\approx 0.5{({f}^{2}+{\sigma }_{\mathrm{eq}}^{2})}^{1/2}\displaystyle \frac{X}{{f}_{b}{\rm{\Delta }}}.\end{eqnarray} \tag{ 28 }$$

After the formation of a PBH binary, the orbit of these two binary PBHs shrinks due to the gravitational waves, and the coalescence time is given by Peters (1964) as

$$\begin{eqnarray}&&t=\displaystyle \frac{3}{85}\displaystyle \frac{{a}^{4}}{{m}_{i}{m}_{j}{m}_{b}}{j}^{7}.\end{eqnarray} \tag{ 29 }$$

Taking into account Equation (15), the dimensionless angular momentum is

$$\begin{eqnarray}&&j(t;X)={\left(\displaystyle \frac{85}{3}\displaystyle \frac{{{tm}}_{i}{m}_{j}{m}_{b}{({f}_{b}{\rm{\Delta }})}^{4}}{{(0.1\bar{x})}^{4}{X}^{16/3}}\right)}^{1/7}.\end{eqnarray} \tag{ 30 }$$

Assuming that PBHs possess a random distribution, the probability distribution of the separation, x, between two nearest PBHs with mass m_i and m_j and without other PBHs in the volume of $\tfrac{4\pi }{3}{x}^{3}$ becomes

$$\begin{eqnarray}&&\displaystyle \frac{{dP}}{d\tilde{X}}={e}^{-\displaystyle \frac{4\pi }{3}{x}^{3}{n}_{T}}={e}^{-\tilde{X}\cdot \displaystyle \frac{4\pi }{3}\langle {x}_{{ij}}{\rangle }^{3}{n}_{T}},\end{eqnarray} \tag{ 31 }$$

where $\tilde{X}\equiv {x}^{3}/\langle {x}_{{ij}}{\rangle }^{3}=X/\mu$, μ is given in Equation (7), and ${n}_{T}\equiv f{\rho }_{\mathrm{eq}}{\int }_{0}^{\infty }\tfrac{P(m)}{m}{dm}$.⁵ Therefore we have

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}P}{d\tilde{X}{dt}}=\displaystyle \frac{1}{7t}{e}^{-\tilde{X}\cdot \displaystyle \frac{4\pi }{3}\langle {x}_{{ij}}{\rangle }^{3}{n}_{T}}\,{ \mathcal P }({\gamma }_{X}),\quad {\gamma }_{X}\equiv \displaystyle \frac{j(t;X)}{{j}_{X}},\end{eqnarray} \tag{ 32 }$$

where j_X is given in Equation (28). The probability distribution of the time of merger becomes

$$\begin{eqnarray}&&\displaystyle \frac{{dP}}{{dt}}=\displaystyle \frac{{\mu }^{-1}}{7t}\int {{dXe}}^{-\displaystyle \frac{X}{\mu }\cdot \displaystyle \frac{4\pi }{3}\langle {x}_{{ij}}{\rangle }^{3}{n}_{T}}{ \mathcal P }({\gamma }_{X}),\end{eqnarray} \tag{ 33 }$$

and the comoving merger rate at time t reads as

$$\begin{eqnarray}&&{R}_{{ij}}(t)\equiv \displaystyle \frac{{{dN}}_{\mathrm{merger}}}{{dtdV}}={\rho }_{m}^{0}\min \left(\displaystyle \frac{{f}_{i}{\rm{\Delta }}}{{m}_{i}},\displaystyle \frac{{f}_{j}{\rm{\Delta }}}{{m}_{j}}\right)\displaystyle \frac{{dP}}{{dt}},\end{eqnarray} \tag{ 34 }$$

where ${\rho }_{m}^{0}\simeq 4\times {10}^{19}\,{M}_{\odot }{\mathrm{Gpc}}^{-3}$ is the matter density at present. Since ${ \mathcal P }({\gamma }_{X})$ has a sharp peak at

$$\begin{eqnarray}&&{X}_{* }(t)\approx 0.032{\left(\displaystyle \frac{t}{{t}_{0}}\right)}^{\displaystyle \frac{3}{37}}{f}_{b}{\rm{\Delta }}{({f}^{2}+{\sigma }_{\mathrm{eq}}^{2})}^{-\displaystyle \frac{21}{74}}{({m}_{i}{m}_{j})}^{\displaystyle \frac{3}{37}}{m}_{b}^{-\displaystyle \frac{1}{37}},\end{eqnarray} \tag{ 35 }$$

if $\displaystyle \frac{{X}_{* }}{\mu }\cdot \tfrac{4\pi }{3}\langle {x}_{{ij}}{\rangle }^{3}{n}_{T}\ll 1$,⁶ the comoving merger rate at time t becomes

$$\begin{eqnarray}&&{R}_{{ij}}(t)={{ \mathcal R }}_{{ij}}(t){{\rm{\Delta }}}^{2},\end{eqnarray} \tag{ 36 }$$

where

$$\begin{eqnarray}{{ \mathcal R }}_{{ij}}(t) & \approx & 3.9\cdot {10}^{6}\times {\left(\displaystyle \frac{t}{{t}_{0}}\right)}^{-\displaystyle \frac{34}{37}}{f}^{2}{({f}^{2}+{\sigma }_{\mathrm{eq}}^{2})}^{-\displaystyle \frac{21}{74}}\\ & & \times \,\min \left(\displaystyle \frac{P({m}_{i})}{{m}_{i}},\displaystyle \frac{P({m}_{j})}{{m}_{j}}\right)\left(\displaystyle \frac{P({m}_{i})}{{m}_{i}}+\displaystyle \frac{P({m}_{j})}{{m}_{j}}\right)\\ & & \times \,{({m}_{i}{m}_{j})}^{\displaystyle \frac{3}{37}}{({m}_{i}+{m}_{j})}^{\displaystyle \frac{36}{37}},\end{eqnarray} \tag{ 37 }$$

which can be interpreted as the comoving merger rate density in units of Gpc⁻³ yr⁻¹, and PBH masses are in units of ${M}_{\odot }$. Again, we want to remind readers that $P({m}_{i})=P({m}_{j})=P(m)/2$ for ${m}_{i}={m}_{j}=m$. If $P(m)/m\,=$ constant, $\tilde{\alpha }\equiv -{({m}_{i}+{m}_{j})}^{2}{\partial }^{2}\mathrm{ln}{{ \mathcal R }}_{{ij}}/\partial {m}_{i}\partial {m}_{j}=36/37$, which is consistent with Kocsis et al. (2018). However, for a general mass function, $\tilde{\alpha }$ can be quite different from 36/37.

Let's consider two typical PBH mass functions in the literature. One takes the power-law form (Carr 1975)

$$\begin{eqnarray}&&P(m)\approx \displaystyle \frac{\alpha -1}{M}{\left(\displaystyle \frac{m}{M}\right)}^{-\alpha }\end{eqnarray} \tag{ 38 }$$

for $m\geqslant M$ and $\alpha \gt 1$, and the other has a lognormal distribution (Dolgov & Silk 1993) of

$$\begin{eqnarray}&&P(m)=\displaystyle \frac{1}{\sqrt{2\pi }\sigma m}\exp \left(-\displaystyle \frac{{\mathrm{log}}^{2}(m/{m}_{c})}{2{\sigma }^{2}}\right).\end{eqnarray} \tag{ 39 }$$

LIGO/VIRGO gives the merger rate for BBHs with ${m}_{1},{m}_{2}\,\geqslant 5\,{M}_{\odot }$ and ${m}_{1}+{m}_{2}\leqslant 100\,{M}_{\odot }$ as ${R}_{T}=12\sim 213$ Gpc⁻³ yr⁻¹ in Abbott et al. (2017d). Similar to Ali-Hamoud et al. (2017), we take ${\sigma }_{\mathrm{eq}}\approx 0.005$. Figure 1 indicates that the merger rate constrained by LIGO/VIRGO can be explained by mergers of PBH binaries. Here, for the purpose of illustration, we simply take $M=5\,{M}_{\odot }$ and $\alpha =1.6$ for the power-law PDF, and ${m}_{c}=15\,{M}_{\odot }$ and $\sigma =0.6$ for the lognormal PDF. Therefore LIGO/VIRGO implies that $1.4\times {10}^{-3}\lesssim {f}_{\mathrm{pbh}}\lesssim 6.6\times {10}^{-3}$ for the power-law PDF and $1.2\times {10}^{-3}\lesssim {f}_{\mathrm{pbh}}\lesssim 5.7\times {10}^{-3}$ for the lognormal PDF. Such an abundance of PBHs is consistent with current constraints from other observations (Chen et al. 2016; Green 2016; Horowitz 2016; Ali-Hamoud & Kamionkowski 2017; Aloni et al. 2017; Carr et al. 2017; Clesse & García-Bellido 2017; Gaggero et al. 2017; Green 2017; Guo et al. 2017; Inoue & Kusenko 2017; Khnel & Freese 2017; Schutz & Liu 2017; Zumalacarregui & Seljak 2017; Wang et al. 2018). In order to break the degeneracy for different PBH mass functions, we need more information. Keeping the total merger rate of R_T = 100 Gpc⁻³ yr⁻¹ fixed, we obtain ${f}_{\mathrm{pbh}}=4.3\times {10}^{-3}$ for the power-law PDF and ${f}_{\mathrm{pbh}}=3.7\times {10}^{-3}$ for the lognormal PDF, respectively, and then plot the one-dimensional merger rate distribution (where we integrate the mass of the lighter black hole in the binary from $5\,{M}_{\odot }$ over the mass of the heavier black hole) in Figure 2. We see that one-dimensional merger rate distributions for the lognormal and power-law PDFs are quite different from each other even though both PDFs give the same total merger rate, R_T. Furthermore, more information will be obtained in the two-dimensional merger rate distributions in Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to compare to the events detected by LIGO, we need to take the sensitivity of LIGO into account. Since LIGO probes mergers approximately in the redshift range of $z\in [0,1]$, the expected number of triggers, Λ, is therefore estimated as (Abbott et al. 2016c, 2016b, 2017b; Kavanagh et al. 2018)

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{{ij}}={\int }_{0}^{1}{R}_{{ij}}(z)\displaystyle \frac{d\langle {VT}\rangle }{{dz}}{dz},\end{eqnarray} \tag{ 40 }$$

where $\langle {VT}\rangle$ is the averaged sensitive spacetime volume of LIGO, which depends on the masses of the merging BBH. We adopt the semi-analytical approximation presented in Veitch et al. (2015), Abbott et al. (2016b, 2016c), and Usman et al. (2016) to calculate $\langle {VT}\rangle$. Here we assume LIGO's O1 and O2 runs share a common sensitive volume, with an observing time of 48.6 days for O1 (Abbott et al. 2016e) and 117 days for O2 (Abbott et al. 2017e). Note that Λ is not the mean number of confidently detected BBH events, but the mean number of signals above the chosen threshold (Abbott et al. 2016c). The two-dimensional distributions of Λ, along with the six events detected by LIGO/VIRGO, are then shown in Figure 4. Since there are only a few events available so far, we cannot make a decisive conclusion as to which PDF fits the data better. As the data accumulating, we may finally pin down the shape of the black hole mass function and nail down whether lower and/or upper mass cutoffs (Fishbach & Holz 2017; Kovetz et al. 2017) of merging black hole binaries exist or not.

Zoom In Zoom Out Reset image size

In this paper we work out the merger rate distribution of PBH binaries with a general mass function, taking into account the torques by all other PBHs and the linear density perturbations. In Ali-Hamoud et al. (2017), the effects of the tidal field from the smooth halo, the PBH binary encountering with other PBHs, the baryon accretion, and present-day halos were carefully investigated, and they concluded that all of these effects make no significant contribution to the overall merger rate. So it is reasonable for us to neglect these subdominant effects throughout our estimation as well. In addition, we find that the evolution of the merger rate of PBH binaries is ${t}^{-34/37}$, which is quite different from that for the astrophysical black hole binaries.

For the power-law and lognormal cases, the abundance of PBH is roughly constrained to the range of ${10}^{-3}\lesssim {f}_{\mathrm{pbh}}\lesssim {10}^{-2}$, which is consistent with other observations (Chen et al. 2016; Green 2016; Horowitz 2016; Ali-Hamoud & Kamionkowski 2017; Aloni et al. 2017; Carr et al. 2017; Clesse & García-Bellido 2017; Gaggero et al. 2017; Green 2017; Guo et al. 2017; Inoue & Kusenko 2017; Khnel & Freese 2017; Schutz & Liu 2017; Zumalacarregui & Seljak 2017; Wang et al. 2018). Our results hence confirm that the dominant fraction of DM should not originate from the stellar mass PBHs (Sasaki et al. 2016; Ali-Hamoud et al. 2017; Raidal et al. 2017; Kocsis et al. 2018).

A possible method to discriminate a PBH scenario from the others is to measure the merger rate distribution of black hole binaries. In particular, we show that the one-dimensional and two-dimensional merger rate distributions are quite sensitive to the mass function of PBHs. In the near future, many more coalescences of black hole binaries will be detected and provide the distribution of the binary parameters, and hence the mass function of a PBH (or black hole) could be reconstructed from the one-dimensional and two-dimensional merger rate distributions. This may finally help us answer the questions: What is the origin of black holes detected by LIGO and VIRGO collaborations and how are black hole binaries formed?