Effects of the surrounding primordial black holes on the merger rate of primordial black hole binaries

TL;DR

This paper presents an analytic model for calculating the merger rate of primordial black hole binaries, considering effects of surrounding black holes and density perturbations, revealing mass function independence in key parameters.

Contribution

It introduces a formalism that accounts for torques from surrounding black holes and density perturbations, providing insights into merger rate calculations independent of the mass function.

Findings

The parameter α is found to be 36/37, independent of the mass function.

The ratio of merger rates considering surrounding torques to nearest black hole effects is mass-independent.

The formalism enables more accurate predictions of primordial black hole binary mergers.

Abstract

We develop an analytic formalism for computing the merger rate of primordial black hole binaries with a general mass function by taking into account the torques by the surrounding primordial black holes and linear density perturbations. We find that $\alpha=-(m_{i}+m_{j})^2\partial^{2} \ln {\mathcal R}(m_{i},m_{j})/\partial m_{i} \partial m_{j}=36/37$ is independent of the mass function. Moreover, the ratio of the merger rate density of primordial black hole binaries by taking into account the torques by the surrounding primordial black holes to by the nearest primordial black hole is independent of the masses of binaries.

Equations

$${\int }_0^{\infty }P(m)dm=1.$$

$$fP(m)dm,$$

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

$$n_T\equiv f{\rho }_{\mathrm{eq}}{\int }_0^{\infty }\frac{P(m)}{m}dm.$$

$$\frac{1}{m_{\mathrm{pbh}}}={\int }_0^{\infty }\frac{P(m)}{m}dm,$$

$$F\left(m\right)\equiv \frac{n\left(m\right)}{n_T}=P\left(m\right)\frac{m_{\mathrm{pbh}}}{m},$$

$${\int }_0^{\infty }F(m)dm=1,$$

$$\sum _{l=1}^{L}F(m_l)d{m}_l=1,$$

$$\ddot{r}-\left(\dot{H}+H^2\right)r+\frac{m_b}{r^2}\frac{r}{|r|}=0,$$

$${\chi }^{\prime \prime }+\frac{s{h}^{\prime }+h}{s^{2}h}\left(s{\chi }^{\prime }-\chi \right)+\frac{1}{\lambda }\frac{1}{{\left(sh\right)}^2}\frac{1}{{\chi }^2}\frac{\chi }{|\chi |}=0,$$

$$\lambda =\frac{8\pi {\rho }_{\mathrm{eq}}{x}^3}{3m_b},$$

$$a\approx 0.1\lambda x.$$

$$dP$$

$${\hat{x}}_{ij}$$

$${\bar{x}}_{ij}={\left(\frac{3}{4\pi }\right)}^{1/3}{n}_T^{-1/3}\approx 0.620n_T^{-1/3},$$

$$X\equiv (x/\overline{x}{)}^3.$$

$$j\equiv \sqrt{1-e^2},$$

$$\mathbf{j}\approx 3\frac{m}{m_b}\frac{x^3}{y^3}\left(\hat{\mathbf{x}}\cdot \hat{\mathbf{y}}\right)\left(\hat{\mathbf{x}}\times \hat{\mathbf{y}}\right).$$

$${\mathbf{j}}_l\approx 3\sum _{p=1}^{N_l}\frac{m_l}{m_b}\frac{x^3}{y_{l,p}^3}(\hat{\mathbf{x}}\cdot {\hat{\mathbf{y}}}_{l,p})(\hat{\mathbf{x}}\times {\hat{\mathbf{y}}}_{l,p}).$$

$$\mathbf{j}=\sum _{l=1}^L{\mathbf{j}}_l\approx 3\sum _{l=1}^L\sum _{p=1}^{N_l}\frac{m_l}{m_b}\frac{x^3}{y_{l,p}^3}(\hat{\mathbf{x}}\cdot {\hat{\mathbf{y}}}_{l,p})(\hat{\mathbf{x}}\times {\hat{\mathbf{y}}}_{l,p}).$$

$$\frac{dP}{d^{2}j}=\underset{V\to \infty }{\lim }\prod _{l=1}^L\prod _{p=1}^{N_l}{\int }_V\frac{d^3{y}_{l,p}}{V}{\delta }_{\mathrm{D}}\left[\mathbf{j}-3\sum _{n=1}^L\sum _{q=1}^{N_n}\frac{m_n}{m_b}\frac{x^3}{y_{n,q}^5}{y}_{n,q||}{\mathbf{y}}_{n,q\perp }\right],$$

$${\delta }_{\mathrm{D}}(\mathbf{X})={\int }_{\mathbf{k}\perp \hat{x}}\frac{d^{2}k}{(2\pi {)}^2}{\mathrm{e}}^{i\mathbf{k}\cdot \mathbf{X}}.$$

$$\frac{dP}{d^{2}j}$$

$${\mathcal{I}}_l$$

$$\underset{V\to \infty }{\lim }$$

$${\mathcal{J}}_l$$

$$A_l(v)=\int \frac{d^2\hat{y}}{4\pi }\left(1-{\mathrm{e}}^{2i\frac{m_l}{m_b}v(\hat{\mathbf{y}}\cdot \hat{\mathbf{x}})(\hat{\mathbf{y}}\cdot \hat{\mathbf{k}})}\right).$$

$$\hat{\mathbf{y}}\cdot \hat{\mathbf{x}}$$

$$\hat{\mathbf{y}}\cdot \hat{\mathbf{k}}$$

$$(\hat{\mathbf{y}}\cdot \hat{\mathbf{x}})(\hat{\mathbf{y}}\cdot \hat{\mathbf{k}})=\frac{\sin 2\varphi }{2}(1-{\mu }^2),$$

Full text

Primordial black holes (PBHs) produced in the early Universe due to the collapse of large energy density fluctuations, as a promising candidate for dark matter (DM), have recently attracted much attention Hawking (1971); Carr and Hawking (1974); Carr (1975); Gao and Guo (2018); Cai et al. (2018); Chen et al. (2018); Khlopov (2010); Belotsky et al. (2018); Saito and Yokoyama (2009).

On the other hand, the gravitational wave events observed by the LIGO detectors Abbott et al. (2016a) can be explained by the coalescence of PBH binaries Kashlinsky (2016); Bird et al. (2016); Clesse and García-Bellido (2017); Sasaki et al. (2016); Nakamura et al. (1997); Ali-Haïmoud et al. (2017); Bird et al. (2016); Nishikawa et al. (2017); Raidal et al. (2017); Kocsis et al. (2018); Sasaki et al. (2018); Kocsis et al. (2018); Raidal et al. (2017). At present, we do not know how to discriminate the PBH scenario from other astrophysical scenarios. In Ref. Kocsis et al. (2018) a new method is proposed to test the PBH scenario. A quantity $\alpha$ can be constructed from the mass distribution of the merger rate density per unit cosmic time and comoving volume ${\cal R}(m_{i},m_{j})$. It is found that the quantity, $\alpha=-(m_{i}+m_{j})^{2}\partial^{2}\ln{\cal R}(m_{i},m_{j})/\partial m_{i}\partial m_{j}$, is independent of the black hole mass function. In the PBH scenario, $\alpha$ is closed to unity. For black holes from dynamical formation in a dense stellar system, $\alpha\approx 4$ O’Leary et al. (2016). Therefore, it provides us a possibility to discriminate PBHs from astrophysical black holes. However, in Ref. Chen and Huang (2018) it is pointed out that if $P(m)/m$ is a constant, $\alpha=36/37$ which is consistent with the result derived in Ref. Kocsis et al. (2018) while for a general mass function, the value of $\alpha$ deviates from $36/37$. In this letter, we re-analyze the merger rate distribution of PBH binaries with a general mass function by taking into account the torques by the surrounding PBHs and linear density perturbations. We find that $\alpha=36/37$ for a general mass function.

In this letter, we use units of $c=G=1$. Whenever relevant, we adopt the values of cosmological parameters from the Planck measurements Ade et al. (2016). The scale factor $s(t)$ is normalized to unity at the matter radiation equality.

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

[TABLE]

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

[TABLE]

where $f\equiv\Omega_{\mathrm{pbh}}/\Omega_{\mathrm{m}}$ is the total abundance of PBHs in non-relativistic matter. The fraction of PBHs in DM is related to $f$ by $f_{\rm{pbh}}\equiv\Omega_{\rm{pbh}}/\Omega_{\rm{dm}}\approx f/0.85$. The average number density of PBHs in mass interval $(m,m+dm)$ at equality is given by

[TABLE]

where $\rho_{\rm{eq}}$ is the energy density of the Universe at the matter radiation equality, and the comoving total average number density of PBHs $n_{T}$ is given by

[TABLE]

For simplicity, we define $m_{pbh}$ as

[TABLE]

and $F(m)$ as

[TABLE]

where $F\left(m\right)$ is the fraction of the average number density of PBHs with mass $m$ in the total average number density of PBHs. It is easily obtained

[TABLE]

which can be rewritten as

[TABLE]

where $m_{1}$ represents the minimum mass of PBH and $m_{L}$ represents the maximum mass of PBH.

Let us consider the condition two nearest PBHs with masses $m_{i}$ and $m_{j}$ decouple from the expanding Universe, assuming negligible initial peculiar velocities here and throughout. Considering the different scaling with time of the two competing effects (their gravitational attraction versus the expansion of the Universe), the equation of motion for their proper separation $r$ in Newtonian approximation is given by

[TABLE]

where the dot denotes differentiation with respect to the proper time and $m_{b}=m_{i}+m_{j}$ is total mass of the PBH binary. By defining $\chi\equiv r/x$, we can rewrite Eq. (9) as

[TABLE]

where $x$ is the comoving distance between these two nearest PBHs and $h(s)\equiv H(s)/\left({8\pi\over 3}\rho_{\rm{eq}}\right)^{1/2}=\sqrt{s^{-3}+s^{-4}}$. Primes denote differentiation with respect to scale factor $s$ and the dimensionless parameter $\lambda$ is

[TABLE]

The solution of Eq. (10) in Ali-Haïmoud et al. (2017) implies the semi-major axis $a$ of the formed binary is given by

[TABLE]

Assuming that the distribution of PBHs is random distribution, the probability distribution of the separation $x$ between two nearest PBHs with masses $m_{i}$ and $m_{j}$ and without other PBHs in the volume of ${4\pi\over 3}x^{3}$ is given by

[TABLE]

The average distance ${\hat{x}}_{ij}$ between two nearest PBHs with masses $m_{i}$ and $m_{j}$ is

[TABLE]

Following Ref. Ali-Haïmoud et al. (2017), we denoted by ${\bar{x}}_{ij}$ characteristic comoving separation between nearest PBHs with mass $m_{i}$ and $m_{j}$,

[TABLE]

which is independent of the PBH binary masses. Therefore, from now on, we omit the subscript ‘ij’ unless it is necessary. Given a comoving separation $x$, we define the dimensionless variable $X$ as

[TABLE]

When the two PBHs come closer and closer, the surrounding PBHs, especially the nearest PBH, will exert torques on the PBH binaries. As the result, the two PBHs avoid a head-on collision with each other. The tidal force will provide a angular momentum to prevent this system from direct coalescence. Here, we introduce a dimensionless angular momentum $j$ defined by

[TABLE]

where $e$ is the eccentricity of the binary at the formation time. From Ali-Haïmoud et al. (2017), the angular momentum generated by a PBH with mass $m$ at a comoving separation $y\geq x$ is given by

[TABLE]

where $\hat{\bm{x}}$ is the unit vector along $\bm{x}$ and where $\hat{\bm{y}}$ is the unit vector along $\bm{y}$ .

We consider $N_{l}$ PBHs with mass $m_{l}$ uniformly distributed within a volume $V=\frac{4\pi}{3}R^{3}$ and take the limit $N,V\rightarrow\infty$ at constant density $n(m_{l})dm_{l}=N_{l}/V$. Therefore, the reduced angular momentum $\bm{j}_{l}$ resulting from $N_{l}$ PBHs with mass $m_{l}$ is given by

[TABLE]

where $y_{l,p}$ is the comoving distance from the binary to the $p$-th PBH with mass $m_{l}$. The total reduced angular momentum $\bm{j}$ resulting from the surrounding PBHs with masses from $m_{1}$ to $m_{L}$ is given by

[TABLE]

Using Eq. (20), the two-dimensional probability distribution for $\bm{j}$ is given by

[TABLE]

where $y_{||}\equiv\bm{\bm{y}}\cdot\hat{\bm{x}}$, $\bm{y}_{\bot}\equiv\hat{\bm{x}}\times\bm{y}$, and $\delta_{\rm D}$ is the two-dimensional Dirac function, which we rewrite as

[TABLE]

We hence get

[TABLE]

where

[TABLE]

When $V\rightarrow\infty$ the latter integral is convergent, then we arrive

[TABLE]

where

[TABLE]

By rescaling $y\rightarrow(1.5k)^{1/3}xy$ and defining $v=1/y^{3}$, the integral $\mathcal{J}_{l}$ becomes

[TABLE]

where

[TABLE]

By using

[TABLE]

we can get

[TABLE]

where $\mu=\cos{\theta}$. From (28) and (30) and using $d^{2}\hat{y}=\sin{\theta}d\theta d\phi$, $A_{l}(v)$ is given by

[TABLE]

where $J_{0}$ is the zeroth-order Bessel function. Since $J_{0}(x)=1+\mathcal{O}(x^{2})$ for $x\rightarrow 0$, we could compute the integral over $v$ first. From (27) and (31), $\mathcal{J}_{l}$ is given by

[TABLE]

where $u=\frac{vm_{l}}{m_{b}}(1-\mu^{2})$. The last two integrals are analytic. By using

[TABLE]

we can get a simple expression

[TABLE]

So, (25) becomes

[TABLE]

Since $m_{l}n(m_{l})=\rho_{l}$ is the energy density of PBHs with mass $m_{l}$, we can get

[TABLE]

So we can arrive

[TABLE]

where

[TABLE]

We hence arrive at the probability distribution

[TABLE]

For a given $X$,

[TABLE]

Similar to Ali-Haïmoud et al. (2017), taking into account both the torques by the surrounding PBHs and density perturbations, we can rewrite the characteristic value of $j_{X}$ as

[TABLE]

where $\sigma_{\rm{eq}}\equiv\langle\delta_{\rm{eq}}^{2}\rangle^{1/2}$ is the variance of density perturbations of the rest of DM at the matter radiation equality. Once the PBHs decouple from the expanding Universe and form a binary, they gradually shrink by gravitational radiation and finally merge. The coalescence time can be estimated as Peters (1964)

[TABLE]

We can rewrite (42) as

[TABLE]

The differential probability distribution of $X$ and $t$ is given by

[TABLE]

From (Effects of the surrounding primordial black holes on the merger rate of primordial black hole binaries) and (16), we can get

[TABLE]

Given that $j\propto t^{1/7}$, $\partial j/\partial t=j/(7t)$. By using (40), we can get

[TABLE]

where

[TABLE]

From Bayes’ theorem, we get the probability distribution of $X$ for binaries merging after a time $t$

[TABLE]

We now find the value of $X_{*}$ to satisfy the probability is maximized. Since $X_{*}\ll 1$, we approximate ${\rm e}^{-X}\approx 1$. The equation we need to solve is

[TABLE]

Since $\gamma_{X}$ is strictly monotonic, this implies $\mathcal{P}^{\prime}(\gamma_{X_{*}})=0$, then we get

[TABLE]

By solving (41), (43) and (50), the most probable value of $X$ for binaries merging at time $t$ is given by

[TABLE]

The probability distribution of the time of PBH merger with mass $m_{i}$ and $m_{j}$ is given by

[TABLE]

Since the integrand peaks at $X_{*}\ll 1$, we get ${\rm e}^{-X}=1$ By using $\gamma_{X}\propto X^{-37/21}$, and $\gamma_{X_{*}}=\sqrt{2}$, we can find

[TABLE]

The total probability distribution of the time of merger is given by

[TABLE]

The merger rate per unit volume at time $t$ is obtained from

[TABLE]

where the factor $1/2$ account for that each merger event involves two PBHs and $z_{\rm eq}\simeq 3400$ is the redshift at the matter radiation equality. Finally, we arrive at

[TABLE]

where

[TABLE]

which can be interpreted as the merger rate density in unit of Gpc*-3* yr*-1*$M_{\odot}^{-2}$. From (57), we confirm the result derived in Kocsis et al. (2018), without considering the merger history of PBHs, $\alpha=-(m_{i}+m_{j})^{2}\partial^{2}\ln\mathcal{R}(t,m_{i},m_{j})/\partial m_{i}\partial m_{j}=36/37$, which is independent of the PBH mass function, by taking into account the torques by the surrounding PBHs.

Only accounting for tidal torquing by the nearest PBH, the merger rate per unit volume at time $t$ is given byLiu et al. (2019)

[TABLE]

where

[TABLE]

which can be interpreted as the merger rate density in unit of Gpc*-3* yr*-1*$M_{\odot}^{-2}$. As a comparation, we define $\lambda$ as

[TABLE]

By using $f_{\rm pbh}=f/0.85$, $\lambda$ is given by

[TABLE]

which is independent of $m_{i},m_{j}$and $t$.

Now let us consider two typical PBH mass functions. One is the monochromatic case Sasaki et al. (2016); Bird et al. (2016); Nishikawa et al. (2017),

[TABLE]

The other takes the power-law form as follows Carr (1975)

[TABLE]

with $m\geq M$ and $\beta>1$. In the monochromatic case, we arrive at

[TABLE]

Neglecting the density perturbation, $\lambda=0.61$ is independent of the mass of PBH. In the power-law case, $\lambda$ is given by

[TABLE]

Neglecting the density perturbation, in the case $\beta>1.05$, $\lambda\sim{\cal O}(1)$ is almost independent on the mass function of PBH.

Fig. 1 and Fig. 2 show that the merger rate constrained by LIGO/VIRGO could be explained by mergers of PBH binaries. Here we take $M=30M_{\odot}$ for the monochromatic case and $M=5M_{\odot},\beta=2.3$ for the power-law case. In the monochromatic case, LIGO/VIRGO implies that $0.0019\lesssim f_{\rm{pbh}}\lesssim 0.014$ for taking into account the torque by the nearest PBH and $0.0013\lesssim f_{\rm{pbh}}\lesssim 0.0098$ for taking into account the torques by all PBHs. In the power-law case, LIGO/VIRGO implies that $0.0014\lesssim f_{\rm{pbh}}\lesssim 0.010$ for taking into account the torque by the nearest PBH and $0.0010\lesssim f_{\rm{pbh}}\lesssim 0.0077$ for taking into account the torques by all PBHs. Such an abundance of PBHs satisfies current constraints from other observations Chen et al. (2016); Green (2016); Schutz and Liu (2017); Wang et al. (2018); Gaggero et al. (2017); Ali-Haïmoud and Kamionkowski (2017); Aloni et al. (2017); Horowitz (2016); Kühnel and Freese (2017); Inoue and Kusenko (2017); Carr et al. (2017); Green (2017); Guo et al. (2017); Poulin et al. (2017).

So far, several gravitational wave events from black hole binary mergers have been detected by LIGO/VIRGO collaboration, such as GW150914 Abbott et al. (2016a), GW151226 Abbott et al. (2016b), GW170104 Abbott et al. (2017a), GW170608 Abbott et al. (2017b) and GW170814 Abbott et al. (2017c). One of the most important question is how to discriminate PBHs and astrophysical black holes. In this letter, we develop a formalism for calculating the merger rate distribution of PBH binaries with a general mass function by taking into account the torques by the surrounding PBHs and linear density perturbations. We find $\alpha=36/37$ is independent of the mass function of PBHs. The ratio of the merger rate density of PBHs by taking into account the torques by the surrounding PBHs to the nearest PBH is independent of $m_{i}$ and $m_{j}$. We apply our formalism to two specific examples, the monochromatic and power-law cases. In these cases, three body approximation is a good approximation. In the future, more and more coalescence events of black hole binaries will be detected by LIGO/VIRGO. This will provide more information of the merger rate distribution of black hole binaries to test the PBH scenario.

Note added

In finishing this letter, we find a parallel independent work Raidal et al. (2018) that has some overlap with our calculation of the merger rate of primordial black hole binaries.