Effects of the merger history on the merger rate density of primordial black hole binaries

TL;DR

This paper introduces a formalism to compute the merger rate density of primordial black hole binaries considering their merger history, highlighting how different mass functions influence the dominance of single or multiple merger events.

Contribution

It presents a novel formalism for calculating merger rate density of primordial black holes with various mass functions, incorporating merger history effects.

Findings

Monochromatic mass function dominated by single-merger events.

Power-law and log-normal mass functions show significant multiple-merger contributions.

Merger history effects vary depending on the black hole mass function.

Abstract

We develop a formalism to calculate the merger rate density of primordial black hole binaries with a general mass function, by taking into account the merger history of primordial black holes. We apply the formalism to three specific mass functions, monochromatic, power-law and log-normal cases. In the former case, the merger rate density is dominated by the single-merger events, while in the latter two cases, the contribution of the multiple-merger events on the merger rate density can not be ignored. The effects of the merger history on the merger rate density depend on the mass function.

Equations

$$x<x_{\max }\equiv (f_{\mathrm{pbh}}/n_{\mathrm{pbh}}{)}^{1/3}={\left(M/{\rho }_{\mathrm{dm}}\right)}^{1/3},$$

$$1+z_{\mathrm{dec}}=(1+z_{\mathrm{eq}}){\left(x_{\max }/x\right)}^3,$$

$$a=A\frac{x}{1+z_{\mathrm{dec}}}=A\frac{{\rho }_{\mathrm{dm}}{x}^4}{(1+z_{\mathrm{eq}})M},$$

$$b=B{\left(\frac{x}{y}\right)}^{3}a,$$

$$j\equiv \sqrt{1-e^2}=B{\left(\frac{x}{y}\right)}^3,$$

$$\tau \simeq \frac{3a^4{j}^7}{170M^3}.$$

$$dP=\frac{4\pi x^{2}dx}{n_{\mathrm{pbh}}^{-1}}\frac{4\pi y^{2}dy}{n_{\mathrm{pbh}}^{-1}}\exp \left(-\frac{4\pi y^3}{3n_{\mathrm{pbh}}^{-1}}\right)\Theta (y-x).$$

$$dP=\frac{4\pi x^{2}dx}{n_{\mathrm{pbh}}^{-1}}\frac{4\pi y^{2}dy}{n_{\mathrm{pbh}}^{-1}}\Theta (y-x)\Theta (y_{\max }-y),$$

$$G(t)=\int dxdy\frac{dP}{dxdy}\Theta (t-\tau (x,y))\Theta (x_{\max }-x).$$

$$x=(\frac{t}{k}{)}^{\frac{1}{37}}{y}^{\frac{21}{37}},$$

$$k=\frac{3}{170}\frac{1}{M^3}(\frac{{\rho }_{\mathrm{dm}}}{(1+z_{\mathrm{eq}})M}{)}^4{A}^4{B}^7.$$

$$x_{\max }=(\frac{t}{k}{)}^{\frac{1}{37}}{y}_{\max }^{\frac{21}{37}}$$

$$f_c$$

$$y_1=x,x=(\frac{t}{k}{)}^{\frac{1}{37}}{y}_1^{\frac{21}{37}},$$

$$y_1=(\frac{t}{k}{)}^{\frac{1}{16}}.$$

$$x_{\max }=(\frac{t}{k}{)}^{\frac{1}{37}}{y}_2^{\frac{21}{37}}$$

$$y_2=(\frac{k}{t}{)}^{\frac{1}{21}}{x}_{\max }^{\frac{37}{21}}$$

$$G(t)$$

$$R\left(t\right)=\frac{1}{2}{n}_{\mathrm{pbh}}\underset{dt\to 0}{\lim }\frac{G\left(t+dt\right)-G\left(t\right)}{dt},$$

$$\begin{array}{rl}& R\left(t\right)\approx \\ & \{\begin{array}{l}1.61\times 10^{12}{\left(\frac{M}{M_{\odot }}\right)}^{-\frac{26}{21}}{\left(\frac{t}{t_0}\right)}^{-\frac{8}{7}}{f}_{\mathrm{pbh}}^3,\mathrm{for}{f}_{\mathrm{pbh}}<f_c,\\ 1.86\times 10^6{\left(\frac{M}{M_{\odot }}\right)}^{-\frac{32}{37}}{\left(\frac{t}{t_0}\right)}^{-\frac{34}{37}}{f}_{\mathrm{pbh}}^{\frac{53}{37}},\mathrm{for}{f}_{\mathrm{pbh}}\ge f_c,\end{array}\end{array}$$

$$x<x_{\max }={\left(\frac{m_b}{2{\rho }_{\mathrm{dm}}}\right)}^{1/3},$$

$$a\approx A\frac{2{\rho }_{\mathrm{dm}}{x}^4}{(1+z_{\mathrm{eq}})m_b},$$

$$j\approx B\frac{2m_l}{m_b}{\left(\frac{x}{y}\right)}^3,$$

$$\tau =\frac{3}{85}\frac{a^4}{m_i{m}_j{m}_b}{j}^7.$$

$$M_2=m_b-E_{\mathrm{GW}}\approx \gamma m_b,$$

$$\int dmP(m)=1.$$

$$fP(m)dm,$$

$$n\left(m\right)dm=\frac{fP(m)dm{\rho }_{\mathrm{m}}}{m}=\frac{f_{\mathrm{pbh}}P(m)dm{\rho }_{\mathrm{dm}}}{m},$$

$$n_T\equiv f_{\mathrm{pbh}}{\rho }_{\mathrm{dm}}\int dm\frac{P(m)}{m}.$$

Full text

I Introduction

Various astrophysical and cosmological observations provide substantial evidences firmly establishing the existence of dark matter (DM) in our Universe. However, the nature of DM remains one of the major unsolved problems in fundamental physics. Primordial black holes (PBHs) produced in the radiation-dominated era of the early universe due to the collapse of large energy density fluctuations, as a promising candidate for dark matter, have attracted much attention Hawking (1971); Carr and Hawking (1974); Carr (1975); Khlopov (2010); Carr et al. (2010, 2016); Gao and Guo (2018); Cai et al. (2018); Sasaki et al. (2018); Saito and Yokoyama (2009); Cai et al. (2019); Chen et al. (2018); Carr et al. (2017); Kannike et al. (2017); Kuhnel and Freese (2019); Kühnel et al. (2016).

Two neighboring PBHs can form a binary in the early Universe and coalesce within the age of the Universe. The merge rate of PBH binaries was first estimated through the three-body interaction for the case where all PBHs have the same mass Nakamura et al. (1997); Ioka et al. (1998). In the PBH binary formation scenario, the gravitational wave event GW150914 detected by LIGO Abbott et al. (2016a) and the merger rate estimated by the LIGO-Virgo Collaboration can be explained by the coalescence of PBH binaries if PBHs have the mass about $30M_{\odot}$ and constitute a tiny fraction of DM Sasaki et al. (2016). The binary formation was extended to account for an arbitrary PBH mass function based on the three-body approximation Raidal et al. (2017) or to account for the torque from the surrounding PBHs as well as standard large-scale adiabatic perturbations assuming a monochromatic mass function Ali-Haïmoud et al. (2017). The mechanism has recently been developed for a general mass function by taking into account the torque from the surrounding PBHs Kocsis et al. (2018); Chen and Huang (2018); Raidal et al. (2018); Liu et al. (2019).

However, these studies ignore the possibility that a PBH binary merges into a new black hole which together with another PBH form a new PBH binary. Such a second-merge event can in principle be detected by LIGO-Virgo at the present time. In this paper, we develop an analytic formalism to work out the merger rate density of PBH binaries with a general mass function, by taking into account the merger history of PBHs.

The paper is organized as follows. In the next section, we summarize the basic equation for the primordial input parameters of PBHs and revisit the merger rate for a monochromatic mass function as the first-merger process. In Sec. III, we develop a formalism to calculate the merger rate density of PBH binaries with a general mass function, by taking into account the merger history of PBHs. In Sec. IV, we consider three specific examples, monochromatic, power-law mass and log-normal functions, to investigate the effects of the merger history on the merger rate density of PBH binaries. The final section is devoted to conclusions.

In this paper, we use units of $c=G=1$. Whenever relevant, we adopt the values of cosmological parameters consistent with the Planck measurements Ade et al. (2016). The scale factor is normalized to unity at the present time.

II Single-merger events

Let us start with deriving the basic equation of the merger rate of PBH binaries. It could be easily checked that the gravitational attraction between two approximately isolated PBHs dominates their dynamics if their average mass is bigger than the background mass contained in a comoving sphere whose radius equals to their conformal distance. Considering the different scaling with time of the two competing effects (their gravitational attraction versus the expansion of the Universe) in the equation of motion for their separation Ali-Haïmoud et al. (2017). Following Ref. Sasaki et al. (2016), in this section, we assume that all PBHs have the same mass, $M$, and PBH binaries decouple from the expansion of the Universe during radiation domination provided that their comoving separation, $x$, approximately satisfies

[TABLE]

where $f_{\rm pbh}$ is the fraction of PBHs in DM, $n_{\rm pbh}$ denotes the comoving average number density of PBHs and $\rho_{\rm dm}$ denotes the present energy density of DM. The redshift $z_{\rm dec}$ at which the binary decoupling occurs is given by

[TABLE]

where $z_{\rm eq}\simeq 3400$ is the redshift at matter-radiation equality, assuming negligible initial peculiar velocities here and throughout. Therefore, given PBH mass $M$ and the initial comoving distance of PBHs $x$, the decoupling time is determined by PBHs. In this work, we assume that accretion and evaporation are negligible before the epoch of binary formation. When two PBHs come closer, the nearest PBH exert torque on the bound system. As a result, the two PBHs avoid a head-on collision and form a highly eccentric binary. The major and minor axes are given by (denoted by $a$ and $b$, respectively)

[TABLE]

[TABLE]

where $y$ is the comoving distance to the third PBH, $A$ and $B$ are numerical factors of ${\cal O}(1)$. A detailed investigation of the dynamics of the binary formation suggests $A=0.4$ and $B=0.8$ Ioka et al. (1998). To be exact, in the following calculation, we adopt $A=0.4$ and $B=0.8$. The dimensionless angular momentum of PBH binaries is given by

[TABLE]

where $e$ is the eccentricity of the binary at the formation time. Once two PBHs form a binary, they gradually shrink through the emission of gravitational radiation and eventually merge at the time $\tau$ after its formation, which can be estimated as Peters (1964)

[TABLE]

To calculate the merger rate of PBH binaries, we have to know the spatial distribution of PBHs. Assuming that the spatial distribution of PBHs is random one, the probability that the comoving distances, $x$ and $y$, are in the intervals $(x,x+dx)$ and $(y,y+dy)$ is given by

[TABLE]

To deal with this probability distribution, we can rewrite Eq. (7) as follows

[TABLE]

where $y_{\rm max}={\left(4\pi n_{\rm pbh}/3\right)}^{-1/3}$, which is adopted in Nakamura et al. (1997). In Fig. 1 we show the merger rate estimated by using the initial distribution (7) and the simplified distribution (8), which indicates that the difference between the two cases is insignificant compared to the uncertainty of the merger rate estimated by the LIGO-Virgo Collaboration.

The fraction of PBHs which have merged before the time $t$ is given by

[TABLE]

In Fig. 2 is a schematic illustration on calculating $G(t)$.

From Eqs. (3), (5) and (6), we can get

[TABLE]

where

[TABLE]

When $f=f_{c}$, there is

[TABLE]

By solving Eq. (12), we can get

[TABLE]

By solving

[TABLE]

we arrive

[TABLE]

By solving

[TABLE]

we arrive

[TABLE]

For $f>f_{c}$, $G(t)$ is given by

[TABLE]

For $f<f_{c}$, $G(t)$ is given by

[TABLE]

Therefore, the merger rate of PBH binaries per unit volume per unit time (at the time $t$) can be easily obtained by

[TABLE]

where the factor $1/2$ accounts for that each merger event involves two PBHs. From Eq. (20), the final result is given by

[TABLE]

which can be interpreted as the merger rate in Gpc*-3* yr*-1*. We show the single-merger rate of PBH binaries as a function of the PBH abundance in Fig. 3. For $f_{\rm{pbh}}>f_{c}$ it scales as $f_{\rm{pbh}}^{53/37}$ and for $f_{\rm{pbh}}<f_{c}$ it scales as $f_{\rm{pbh}}^{3}$.

Now we have to emphasize what is the difference between our formalism and the one developed in Refs Sasaki et al. (2016). For the single-merger case, the merger rate of PBH binaries is usually calculated by converting the probability distribution function $P$ of $x$ and $y$ into the one of $a$ and $e$. However, for the multiple-merger case, the probability distribution is a function of more variables than $x$ and $y$. It becomes hard to convert the probability distribution function of $(x,y,z,...)$ into the one of $a$ and $e$. Therefore, the known formalism does not work in the multiple-merger case. To get the merger rate of PBH binaries, we directly deal with the probability distribution in the $x-y$ plane to find which PBHs have been merged. It becomes easy to extend our formalism to the second and third merger events. In this section, to warm up we consider the merger rate of PBH binaries in the single-merger case. In the next section, we shall extend the formalism to the second- and third-merger cases.

III Multiple-merger events

So far, several gravitational wave events from black hole binary mergers have been detected by the LIGO-Virgo collaboration, such as 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. (2016b), GW170104 ($31.2_{-6.0}^{+8.4}M_{\odot}$, $19.4_{-5.9}^{+5.3}M_{\odot}$) Abbott et al. (2017a), GW170608 ($12_{-2}^{+7}M_{\odot}$, $7_{-2}^{+2}M_{\odot}$) Abbott et al. (2017b) and GW170814 ($30.5_{-3.0}^{+5.7}M_{\odot}$, $25.3_{-4.2}^{+2.8}M_{\odot}$) Abbott et al. (2017c). These events detected by LIGO-Virgo suggest that the black holes should have an extended mass function. In this section, we calculate the merger rate distribution for PBH binaries with a general mass function by taking into account the effect of merger history on the merger rate density of PBH binaries.

First of all, we consider the condition that two neighboring PBHs with the masses $m_{i}$ and $m_{j}$ decouple from the expansion of the Universe and form a bound system. Their comoving separation, $x$, approximately satisfies

[TABLE]

where $m_{b}=m_{i}+m_{j}$ is the total mass of the PBH binary. When two PBHs come closer, the nearest PBH with the mass $m_{l}$, exert torque on the bound system. As a result, the two PBHs avoid a head-on collision and form a highly eccentric binary. The major axis $a$ of the binary orbit and the dimensionless angular momentum are given by

[TABLE]

[TABLE]

where $y$ is the comoving distance to the third PBH with the mass $m_{l}$. Once two PBHs form a binary, they gradually shrink through the emission of gravitational radiation and eventually merge at the time $\tau$ after its formation, which can be estimated as Peters (1964)

[TABLE]

The two neighboring PBHs with the masses $m_{i}$ and $m_{j}$ merge into a bigger black hole. The mass is given by

[TABLE]

where $E_{\rm GW}$ is the energy of gravitational wave and $\gamma$ is a factor of ${\cal O}(1)$. In the monochromatic case, $\gamma=0.95$ is adopted in Bringmann et al. (2018). For simplicity, in this paper, we take $\gamma=1$, which means we assume that the energy of gravitational wave is zero.

In this paper, the probability distribution function of PBHs $P(m)$ is normalized to be

[TABLE]

Therefore, the abundance of PBHs in the mass interval $(m,m+dm)$ can be easily obtained by

[TABLE]

where $f$ is a fraction of PBHs in non-relativistic matter including DM and baryons. The fraction of PBHs in DM $f_{\rm{pbh}}$ is given by $f_{\rm{pbh}}\equiv\Omega_{\rm{pbh}}/\Omega_{\rm{dm}}\approx f/0.85$. At the present time, the average number density of PBHs in the mass interval $(m,m+dm)$ is given by

[TABLE]

where $\rho_{\rm{m}}$ is the total energy density of matter and the present total average number density of PBHs, $n_{T}$, is obtained by

[TABLE]

For simplicity, here we define $m_{\rm pbh}$ as

[TABLE]

We define $F(m)$ as

[TABLE]

which is the fraction of the present average number density of PBHs with the mass $m$ in the present total average number density of PBHs.

The result in Ali-Haïmoud et al. (2017) indicates that in the case of $f_{\rm{pbh}}<f_{c}$, the effects of the linear density perturbations on the merger rate of PBH binaries is significant. Here, we only consider the the case of $f_{\rm{pbh}}>f_{c}$ which is shown to be relevant to the LIGO observations Sasaki et al. (2016). In other words, we ignore the bound (22).

The only essential ingredient that we need is the spatial distribution of PBHs. We firstly consider the spatial distribution of two PBHs. The probability distribution of the comoving separation $x$ between two nearest PBHs with the masses $(m_{i},m_{i}+dm_{i})$ and $(m_{j},m_{j}+dm_{j})$ and without other PBHs in the comoving volume of $4\pi x^{3}/3$ is given by

[TABLE]

Clearly, in the non-monochromatic case, to calculate the merger rate in the first-merger process, the differential probability distribution is given by

[TABLE]

where $x$ is the comoving separation between two nearest PBHs with the masses $m_{i}$ and $m_{j}$ and $y$ is the comoving distance to the third PBH with the mass $m_{l}$ which provides the angular momentum for the bound system. The fraction of PBHs that have merged before the time $t$ is given by

[TABLE]

So, we can arrive

[TABLE]

${\cal R}_{1}(t,m_{i},m_{j},m_{l})$ is given by

[TABLE]

where the factor $1/2$ accounts for that each merger event involves two PBHs. From Eq. (III), one has

[TABLE]

The merger rate density of PBH binaries with the masses $m_{i}$ and $m_{j}$ in the first-merger process is

[TABLE]

Let us estimate the merger rate density in the second-merger process. In the first-merger process, two neighboring PBHs decouple from the expansion of the Universe and then merge into a new black hole with the mass $m_{i}+m_{j}$. In the second-merger process, the new black hole and the nearest PBH with mass $m_{k}$ form a new binary. The merge event of the new binary is detected by LIGO-Virgo at the time $t$. Statistically, the second coalescence time is larger than the first one, therefore, we can ignore the first coalescence time. The differential probability distribution is given by

[TABLE]

So, the fraction of PBHs that have merged in the second-merger process is given by

[TABLE]

Then, we can arrive

[TABLE]

${\cal R}_{2}(t,m_{i},m_{j},m_{l})$ is given by

[TABLE]

where the factor $1/3$ accounts for that each merger event in second-merger process involves three PBHs. From Eq. (III), the final result is given by

[TABLE]

The merger rate density of PBH binaries with the masses $m_{i}$ and $m_{j}$ in the second-merger process is given by

[TABLE]

Similarly, $G_{3}(t,m_{i},m_{j},m_{k},m_{f},m_{l})$ and ${\cal R}_{3}(t,m_{i},m_{j},m_{k},m_{f},m_{l})$ are given by

[TABLE]

[TABLE]

The merger rate density of PBH binaries with the masses $m_{i}$ and $m_{j}$ in the third-merger process is given by

[TABLE]

The total merger rate density of PBH binaries with the masses $m_{i}$ and $m_{j}$ detected by LIGO-Virgo is given by

[TABLE]

In the single-merger case, we have $\alpha=-(m_{i}+m_{j})^{2}\partial^{2}\ln{\cal R}(t,m_{i},m_{j})/\partial m_{i}\partial m_{j}=36/37$ which is independent of the PBH mass function. It is consistent with the result obtained in Kocsis et al. (2018). However, by taking account into the merger history of PBHs, $\alpha$ depends on the PBH mass function, which could help us reconstruct the mass function of PBHs.

IV Applications

The total fraction of PBH binaries that have merged before the time $t$ in single-merger events is given by

[TABLE]

The merger rate of PBH binaries in single-merger events at the time $t$ is given by

[TABLE]

$G_{N}(t)$ is the total fraction of PBH binaries that have merged before the time $t$ in $N$-th merger process and $R_{N}(t)$ is merger rate of PBH binaries at time $t$ in $N$-th merger process.

Let us consider three typical PBH mass functions: monochromatic, power-law and log-normal function.

IV.1 Monochromatic mass function

In this subsection, we consider the following monochromatic mass function Sasaki et al. (2016); Bird et al. (2016); Nishikawa et al. (2017)

[TABLE]

In this case, we can rewrite (31) and (32) as

[TABLE]

[TABLE]

From Eqs. (36), (38), (39), (50), (IV), (53) and (54), the total fraction of PBH binaries that have merged before the time $t$ and the merger rate of PBH binaries at the time $t$ in the first-merger process are given by

[TABLE]

[TABLE]

which is consistent with (20). Similarly, the total fraction of PBH binaries that have merged before the time $t$ and the merger rate of PBH binaries at the time $t$ in the second-merger process are given by

[TABLE]

[TABLE]

In Fig. 4, we show the merger rate of PBH binaries in the second-merger process as a function of $f_{\rm{pbh}}$, which scales as $f_{\rm{pbh}}^{69/37}$. The total fraction of PBH binaries that have merged before the time $t$ and the merger rate of PBH binaries at the time $t$ in the third-merger process are given by

[TABLE]

[TABLE]

In Fig. 5, we show the merger rate of PBH binaries in the third-merger process as a function of $f_{\rm{pbh}}$, which scales as $f_{\rm{pbh}}^{69/37}$. In the case of $M=30M_{\odot}$ and $f_{\rm{pbh}}=0.01$, we can find $R_{1}(t_{0})=187$ Gpc*-3* yr*-1*, $R_{2}(t_{0})=2.35$ Gpc*-3* yr*-1* and $R_{3}(t_{0})=3.29\times 10^{-2}$ Gpc*-3* yr*-1*, as shown in Fig. 6. It indicates that, in the monochromatic case, although the merger events of both $30~{}M_{\odot}-30~{}M_{\odot}$ PBH binaries and $60~{}M_{\odot}-30~{}M_{\odot}$ PBH binaries could occur at the same time, the major of merger events detected by LIGO-Vigo is the merger event of $30~{}M_{\odot}-30~{}M_{\odot}$ PBH binaries. Therefore, in the monochromatic case, the effect of the merger history on the merger rate of PBH binaries is negligible.

IV.2 Power-law mass function

In this subsection, we take the PBH mass function as a power-law form Carr (1975):

[TABLE]

with $500M\geq m\geq M$ and $q>1.5$. In the power-law case, we can rewrite (31) and (32) as

[TABLE]

[TABLE]

Choosing $f_{\rm pbh}=0.01$, $q=2.3$, $M=0.2M_{\odot}$, we can get $R_{1}(t_{0})=9.66\times 10^{3}$ Gpc*-3* yr*-1*, $R_{2}(t_{0})=1.15\times 10^{2}$ Gpc*-3* yr*-1*, $R_{3}(t_{0})=5.00$ Gpc*-3* yr*-1*. In power-law case, the effect of the merger history on the merger rate of PBH binaries is small. However, the effect of the merger history on the merger rate density is significant in some region of the parameter space. For example, ${\cal R}_{1}(t_{0},30M_{\odot},30M_{\odot})=8.55\times 10^{-7}$Gpc*-3* yr*-1*$M_{\odot}^{-2}$, ${\cal R}_{2}(t_{0},30M_{\odot},30M_{\odot})=8.90\times 10^{-7}$Gpc*-3* yr*-1*$M_{\odot}^{-2}$, ${\cal R}_{3}(t_{0},30M_{\odot},30M_{\odot})=4.88\times 10^{-8}$Gpc*-3* yr*-1*$M_{\odot}^{-2}$. In Fig. 7, we show the ratio of the total merger rate density to the single-merger one in the PBH mass plane. There are several gravitational wave events detected by LIGO-Virgo. Masses of black hole all are in $(5M_{\odot},50M_{\odot})$. In such region, in the future, more and more coalescence events of black hole binaries will be detected by LIGO-Virgo Wang et al. (2019); Chen and Huang (2019). When we use the merger rate distribution to fit the mass function of PBH, the effect of merger history on the merger rate density of PBH binaries can not be ignored.

IV.3 Log-normal mass function

In this subsection, we take the PBH mass function as a log-normal form Dolgov and Silk (1993); Green (2016); Kühnel and Freese (2017):

[TABLE]

In the power-law case, we can rewrite (31) and (32) as

[TABLE]

[TABLE]

Choosing $f_{\rm pbh}=0.01$, $m_{c}=15M_{\odot}$, $\sigma=0.5$, we can get $R_{1}(t_{0})=423$ Gpc*-3* yr*-1*, $R_{2}(t_{0})=6.5$ Gpc*-3* yr*-1*, $R_{3}(t_{0})=0.1$ Gpc*-3* yr*-1*. In log-normal case, the effect of the merger history on the merger rate of PBH binaries is also small. According to ${\cal R}_{1}(t_{0},30M_{\odot},30M_{\odot})=2.16\times 10^{-2}$Gpc*-3* yr*-1*$M_{\odot}^{-2}$, ${\cal R}_{2}(t_{0},30M_{\odot},30M_{\odot})=2.14\times 10^{-3}$Gpc*-3* yr*-1*$M_{\odot}^{-2}$, ${\cal R}_{3}(t_{0},30M_{\odot},30M_{\odot})=2.31\times 10^{-5}$Gpc*-3* yr*-1*$M_{\odot}^{-2}$, the effect of the merger history on the merger rate density of PBH binaries could not be negligible in some region of the parameter space. In Fig. 8, we show the ratio of the total merger rate density to the single-merger one in the PBH mass plane in the case of $f_{\rm pbh}=0.01$, $m_{c}=15M_{\odot}$ and $\sigma=0.5$. In Fig. 9, we also plot the contour of $({\cal R}(t_{0},30M_{\odot},30M_{\odot})/{\cal R}_{1}(t_{0},30M_{\odot},30M_{\odot})-1)$ in the parameter space of PBH mass function to show that the effect of the merger history on the merger rate density depend on the mass function.

V Conclusions

We have developed the formalism to calculate the merger rate density of PBH binaries with a general mass function, by taking into account the merger history of PBHs. In the monochromatic case, we find that $R_{1}\gg R_{2}\gg R_{3}$, which is independent on $f_{\rm pbh}$. Therefore, the effect of the merger history on the merger rate of PBH binaries is negligible. However, the multiple-merger events may play an important role in the merger rate density of PBH binaries in the non-monochromatic case. For example, for the power-law and log-normal mass function, the effect of the merger history on the merger rate density of PBH binaries could not be negligible. In the future, more and more coalescence events of black hole binaries will be detected by LIGO-Virgo. This will provide more rich information on the merger rate distribution of black hole binaries to test the PBH scenario.

We calculate the merger rate density of PBH binaries up to three mergers. In principle, one can directly calculate it at more than three mergers by using the formalism developed in the present paper. Since the contribution of the merger history on the merger rate density of PBH binaries depends on the mass function and the mass region, it is hard to judge whether mergers of higher order should be computed for a generic mass function.

The effects of the tidal field from the smooth halo, the encountering with other PBHs, the baryon accretion and present-day halos, are carefully investigated in Ali-Haïmoud et al. (2017). It is found in Ali-Haïmoud et al. (2017) that these effects make no significant contributions to the overall merger rate. We therefore neglected these subdominant effects throughout our computation.

Acknowledgements.

This work is supported in part by the National Natural Science Foundation of China Grants No.11575272, No.11435006, No.11690021, No.11690022, No.11851302 and No.11821505, in part by the Strategic Priority Research Program of the Chinese Academy of Sciences Grant No. XDB23030100, No. XDA15020701 and by Key Research Program of Frontier Sciences, CAS.