Formation and Evolution of Primordial Black Hole Binaries in the Early Universe

TL;DR

This paper analyzes the formation and evolution of primordial black hole binaries in the early universe, deriving merger rates, performing N-body simulations, and constraining PBH abundance based on gravitational wave observations.

Contribution

It provides analytical estimates and N-body simulation results for PBH binary merger rates, highlighting the impact of initial conditions and PBH abundance on observable gravitational wave signals.

Findings

Analytical merger rate overestimates when PBHs are all dark matter.

Most initial binaries are disrupted by surrounding PBHs.

LIGO constraints limit PBH abundance in certain mass ranges.

Abstract

The abundance of primordial black holes (PBHs) in the mass range $0.1 - 10^3 M_\odot$ can potentially be tested by gravitational wave observations due to the large merger rate of PBH binaries formed in the early universe. To put the estimates of the latter on a firmer footing, we first derive analytical PBH merger rate for general PBH mass functions while imposing a minimal initial comoving distance between the binary and the PBH nearest to it, in order to pick only initial configurations where the binary would not get disrupted. We then study the formation and evolution of PBH binaries before recombination by performing N-body simulations. We find that the analytical estimate based on the tidally perturbed 2-body system strongly overestimates the present merger rate when PBHs comprise all dark matter, as most initial binaries are disrupted by the surrounding PBHs. This is mostly due to the formation of compact N-body systems at matter-radiation equality. However, if PBHs make up a small fraction of the dark matter, $f_{\rm PBH} \lesssim 10\%$, these estimates become more reliable. In that case, the merger rate observed by LIGO imposes the strongest constraint on the PBH abundance in the mass range $2 - 160 M_\odot$. Finally, we argue that, even if most initial PBH binaries are perturbed, the present BH-BH merger rate of binaries formed in the early universe is larger than $\mathcal{O}(10)\,{\rm Gpc}^{-3} {\rm yr}^{-1}\, f_{\rm PBH}^3$

Figures

Equations

$$H^2=\frac{8\pi }{3}\left({\rho }_M{a}^{-3}+{\rho }_R{a}^{-4}\right),$$

$$\mathrm{d}{s}^2=-(1+2\varphi (\mathbf{x}))\mathrm{d}{t}^2+(1-2\varphi (\mathbf{x}))a(t{)}^2\mathrm{d}{x}^2,$$

$$a^{-2}\Delta \varphi =4\pi {\rho }_{\mathrm{PBH}}(\mathbf{x})+4\pi {\bar{\rho }}_{\mathrm{M}}{\delta }_{\mathrm{M}}(\mathbf{x}).$$

$$\varphi (\mathbf{x})=-\sum _i\frac{m_i}{a|{\mathbf{x}}_i-\mathbf{x}|}-4\pi {\rho }_{\mathrm{DM}}\int \frac{{\mathrm{d}}^{3}k}{(2\pi {)}^3}\frac{a^2}{k^2}{e}^{-i\mathbf{k}\cdot \mathbf{x}}{\tilde{\delta }}_{\mathrm{DM}}(\mathbf{k}),$$

$$S^{(N)}=\int{\rm d}t\bigg{[}\sum_{i}m_{i}\left(\frac{1}{2}\dot{\bf r}_{i}^{2}+\frac{1}{2}\frac{\ddot{a}}{a}{\bf r}_{i}^{2}-\phi_{\rm ex}({\bf x})\right)+\sum_{i>j}\frac{m_{i}m_{j}}{|{\bf r}_{i}-{\bf r}_{j}|}\bigg{]}\,,$$

$$\displaystyle S^{(2)}\approx\int{\rm d}t\bigg{[}$$

$$\mathbf{F}/\mu =\underset{\text{mbox}Hubbleflow}{\underbrace{\mathbf{r}\ddot{a}/a}}-\underset{\text{mbox}self-gravity}{\underbrace{M\hat{\mathbf{r}}/r^2}}+\underset{\text{mbox}radialtidalforces}{\underbrace{(\hat{\mathbf{r}}\cdot \mathbf{T}\cdot \mathbf{r})\hat{\mathbf{r}}}}+\underset{\text{mbox}tidaltorque}{\underbrace{(\mathbf{r}\times (\mathbf{T}\cdot \mathbf{r}))}}\times (\hat{\mathbf{r}}/r),$$

$${\delta }_{\mathrm{b}}\equiv \frac{M/2}{{\rho }_{\mathrm{M}}V(x_0)},$$

$$a_{\mathrm{dc}}\equiv a_{\mathrm{eq}}/{\delta }_{\mathrm{b}},$$

$$\mathbf{L}=\mu \int \mathrm{d}t\mathbf{r}\times (\mathbf{T}\cdot \mathbf{r}).$$

$$\mathbf{j}\equiv \frac{\mathbf{L}/\mu }{\sqrt{r_{a}M}},$$

$$\tau =\frac{3}{85}\frac{r_a^4}{\eta M^3}{j}^7,$$

$$x^{\prime \prime }+\frac{a}{a_{\mathrm{dc}}}\frac{x_0^3}{x^2}=0,$$

$$r_a=a{x}_a/2\approx 0.1a_{\mathrm{dc}}{x}_0.$$

$$\mathbf{j}$$

$$r_a=c_a\frac{8\pi {\rho }_R}{3}\frac{x_0^4}{M},\mathbf{j}=c_j\frac{x_0^3}{M}\hat{\mathbf{r}}\times (\mathbf{T}{a}^3\cdot \hat{\mathbf{r}}).$$

$$\tau =\frac{4096{\pi }^4{c}_a^4{c}_j^7{\rho }_R^4}{2295}\frac{x_0^{37}}{\eta M^{14}}|\hat{\mathbf{r}}\times (\mathbf{T}{a}^3\cdot \hat{\mathbf{r}}){|}^7.$$

$$\mathrm{d}{n}_{\mathrm{b}}$$

$$\mathrm{d}R$$

$$\displaystyle=\frac{1}{14\tau}{\rm d}n(m_{1}){\rm d}n(m_{2})\int{\rm d}V(x_{0})e^{-\bar{N}(y)}j\frac{{\rm d}P(j|x_{0},y)}{{\rm d}j}\bigg{|}_{j=j(\tau)}\,,$$

$${\mathbf{j}}_1({\mathbf{x}}_i,m_i)=j_0\frac{m_i}{⟨m⟩}\frac{3}{\bar{N}(x_i)}{\hat{\mathbf{x}}}_i\times \hat{\mathbf{r}}({\hat{\mathbf{x}}}_i\cdot \hat{\mathbf{r}}),$$

$$j_0\equiv c_j\bar{N}(x_0)⟨m⟩/M\approx 0.4f_{\mathrm{PBH}}/{\delta }_b,$$

$$\frac{\mathrm{d}P}{{\mathrm{d}}^{3}j}\equiv ⟨\delta (\mathbf{j}-{\mathbf{j}}_{\mathrm{M}}-{\mathbf{j}}_{\mathrm{PBH}})⟩=\int \frac{{\mathrm{d}}^{3}k}{(2\pi {)}^3}{e}^{-i\mathbf{k}\cdot \mathbf{j}+K(\mathbf{k})},$$

$$⟨\mathbf{j}⟩=0,⟨\mathbf{j}\otimes \mathbf{j}⟩=\frac{1}{2}{\sigma }_j^2(1-\hat{\mathbf{r}}\otimes \hat{\mathbf{r}}),$$

$${\sigma }_{j,\mathrm{M}}^2$$

$$=\frac{6}{5}{j}_0^2\frac{{\sigma }_{\mathrm{M}}^2}{f_{\mathrm{PBH}}^2},$$

$${\sigma }_{j,\mathrm{PBH}}^2$$

$$=\frac{6}{5}{j}_0^2\frac{1+{\sigma }_m^2/⟨m{⟩}^2}{\bar{N}(y)},$$

$${\sigma }_j^2={\sigma }_{j,\mathrm{M}}^2+{\sigma }_{j,\mathrm{PBH}}^2=\frac{6}{5}{j}_0^2\left(\frac{1+{\sigma }_m^2/⟨m{⟩}^2}{\bar{N}(y)}+\frac{{\sigma }_{\mathrm{M}}^2}{f_{\mathrm{PBH}}^2}\right).$$

$$K_{\mathrm{M}}(\mathbf{k})=-\frac{1}{2}⟨(\mathbf{k}\cdot \mathbf{j}{)}^2⟩=-\frac{1}{4}{\sigma }_{j,\mathrm{M}}^2{k}_{\perp }^2.$$

$$K_{\mathrm{PBH}}(\mathbf{k})$$

$$-z\int \frac{{\mathrm{d}}^2\Omega }{4\pi }{\int }_0^z\frac{\mathrm{d}u}{u^2}\left[\exp \left(i3u({\hat{\mathbf{k}}}_{\perp }\cdot \hat{\mathbf{x}})(\hat{\mathbf{x}}\cdot \hat{\mathbf{r}})\right)-1\right]$$

Full text

1 Introduction

The Advanced Laser Interferometer Gravitational-Wave Observatory (LIGO) observed ten black hole (BH) binary mergers during its first runs [1, 2, 3, 4, 5, 6]. The constituents of these binaries were relatively heavy, with masses in the range $7-50M_{\odot}$, and it remains unclear whether they have astrophysical or primordial origin [7, 8, 9]. Observables that may probe the origin of these BH binaries include the distribution of their eccentricities [10], masses [11, 12, 13, 14, 15, 16] and spins [17, 18, 16], the redshift dependence of the BH merger rate [16] and the correlation between the gravitational wave (GW) sources and galaxies [19, 20] or dark matter (DM) spikes [21].

A population of primordial black holes (PBHs) may also contribute to the DM abundance, but the possibility that PBHs of any mass make up all of the DM is strongly constrained by a wide range of experimental observations (see e.g. Refs. [22, 23, 24]). However, recent revisions of the femtolensing [25] and the HSC/Subaru [26] surveys have opened new windows for PBH DM for light PBHs with masses below $10^{-11}M_{\odot}$ [27, 28]. Heavier PBH DM, in the mass range $0.1-10^{3}M_{\odot}$, is constrained by microlensing [29, 30], survival of stars in dwarf galaxies [31, 32], the distribution of wide binaries [33] and by modification of the cosmic microwave spectrum due to accreting PBHs [34, 35, 36, 37]. Also the recent 21cm observations by EDGES experiment [38] constrain the energy injection from accreting PBHs putting bounds on the abundance of PBHs in this mass window [39].

The strongest potential bounds on the PBH abundance in a mass range around $10M_{\odot}$ can be derived from LIGO observations. If PBHs made up a significant fraction of DM, they would produce a BH binary merger rate and a gravitational wave (GW) background much larger than what is observed by LIGO [40, 13, 41]. These constraints are, however, subject to large theoretical uncertainties. A reliable estimate of the GW signatures from PBH binary mergers requires a good understanding of both the formation of the PBH binaries and of their subsequent interactions with the surrounding matter that may disrupt the binaries.

In the most common scenario, PBHs are formed when large curvature fluctuations in the early universe directly collapse gravitationally to form BHs [42, 43]. Independent of the formation mechanism, they are thereafter dynamically coupled to cosmic expansion and thus their peculiar velocities are negligible. PBHs become gravitationally bound to each other roughly when their local density becomes equal to the density of surrounding radiation, which generically happens after matter-radiation equality. However, due to large Poisson fluctuations at small scales, some PBHs can decouple much earlier. When this happens, the closest PBH pairs start falling towards each other and their head-on collision may be prevented only by the torque caused by the gravitational field of the surrounding PBHs and other matter inhomogeneities. As a result a population of PBH binaries is formed.

This mechanism was first described about two decades ago in Ref. [44]. It was assumed that all PBHs have the same mass, were uniformly distributed in space, all torque was provided by the PBH closest to the binary, and the early binaries were not disrupted between their formation and merger. This model has since been significantly improved [45, 46, 47, 13, 48, 49, 50]. In this paper we provide a self contained extension of the merger rate derivation of Ref. [48], that included the torques from surrounding PBH and matter inhomogeneities, to a general PBH mass function. We also account for the necessary separation between the initial binary and the surrounding PBHs.

The merger rate may be modified by interactions of the binaries with surrounding matter after their formation. Estimates based on hierarchical 3-body systems show that excluding initial conditions where the pair becomes bound to the nearest PBH cuts the merger rate in half [45]. The required initial distance between the binary and the third PBH was found to be larger than the average distance between the PBHs. In that case, however, it is likely that the 3-body system is also coupled to other surrounding PBHs, so a full $N$-body analysis is needed to determine which binaries are disrupted by the nearest PBH. In Ref. [48] the disruption of PBH binaries after the formation of the first virialised haloes consisting of 10 or more PBHs was estimated, using simple analytic arguments, to have a negligible effect on the merger rate. However, the period between formation of the binary and formation of the first haloes was not considered. In Ref. [49] the interaction with the surrounding CDM was shown to have only a mild effect on the merger rate, because of the high eccentricity of the PBH binaries merging today.

In this paper, we focus on the earliest stages of structure formation before recombination. To check if the binaries are significantly disrupted during this epoch, we perform $70$-body simulations that model the evolution of PBHs for the first 377 kyr. A large number of surrounding PBHs makes it possible to numerically test the analytical predictions for the statistical distribution of orbital characteristics of the initial binaries. Because most binaries with a close third PBH are disrupted, we find that the distribution of eccentricities is better approximated by a Gaussian distribution than the broken power law found in Ref. [48]. If PBHs comprise most of the DM, then they (including any binaries) will rapidly form small $N$-body systems beginning matter-radiation equality. Due to this effect, the simulations show a much larger disruption rate than predicted in Ref. [48]. On the other hand, if PBHs make up only a small fraction of DM, bound structures of PBHs form much later and the early disruption rate becomes negligible.

We re-evaluate the LIGO constraints on the PBH abundance from the observed BH-BH merger rate and from non-observation of the stochastic GW background by accounting for the possible suppression of the merger rate due to the interactions of the binary with the surrounding matter. The BH-BH mergers observed by LIGO require a relatively narrow mass function centred around $20M_{\odot}$ while the merger rate can be reproduced if about $0.2\%$ of DM consists of PBHs. We stress, that these constraints are tentative, as they neglect the late disruption rate, which is likely significant for a large PBH fraction, and they also ignore the contribution from the population of perturbed initial PBH binaries. We argue, however, that, due to the hardness of initial binaries, the stochastic GW background and the merger rate should be observable even if the initial binary population is almost completely perturbed, which is the case if PBHs make up most of DM.

The paper is organised as follows. The analytical estimates of the PBH merger rate are derived in Sec. 2. In Sec. 3 we compare the analytical results to numerical simulations and study the early interactions of PBH binaries with the surrounding PBHs. In Sec. 4 we consider the gravitational wave phenomenology of PBH mergers and revise the constraints on the fraction of DM in PBHs. We summarise our main results in Sec. 5.

Geometric units, $G=c=1$, are used throughout this work.

2 Merger rate from early PBH binary formation

In this section we derive the PBH merger rate under the assumption that the initially formed binary population is not disturbed. The initial condition consists of two PBHs with masses $m_{1}$ and $m_{2}$, proper separation $r$, and a sphere of comoving radius $y$ which contains no other PBHs. The distribution of the initial PBH binaries is set by the distribution of their initial conditions and subsequent emission of GWs. To exclude the possibility of early disruption by nearby PBHs, we will consider a range of possible initial conditions that gives a conservative estimate for the merger rate.

The binary is formed from a pair of close PBH. The surrounding matter contributes with its gravitational potential, which we expand around the centre of mass of the pair, and acts on the pair via a tidal force that generates the angular momentum of the binary. The derivation of the merger rate can be divided into two parts:

• •

The calculation of the orbital parameters of the binaries for a given initial separation and a known distribution of surrounding BHs and other matter inhomogeneities. From the orbital parameters we can estimate the coalescence time of the binary. This will be discussed in Sec. 2.1.

• •

From the probability distribution of possible initial conditions that will lead to binary formation we can then derive the distribution of the orbital parameters and coalescence times. This will be the focus of Sec. 2.2.

We compare these analytic estimates for the formation of the initial binaries with numerical results from $N$-body simulations in Sec. 3.

2.1 Dynamics of early binary formation

Consider first the dynamics of a PBH pair in an expanding background surrounded by other PBHs and matter with small inhomogeneities. At that time the Hubble parameter can be expressed as

[TABLE]

where $\rho_{\rm M}$ and $\rho_{\rm R}$ denote the comoving energy densities of matter and radiation, respectively, and the scale factor is chosen so that today $a=1$. The scale factor at matter-radiation equality is $a_{\rm eq}\equiv\rho_{R}/\rho_{M}\approx 1/3400$.

In the comoving Newtonian gauge and in the absence of anisotropic stress, the line element for spacetimes with small inhomogeneities is [45]

[TABLE]

where the potential $\phi$ is determined by

[TABLE]

Here $\rho_{\rm PBH}$ and $\bar{\rho}_{\rm M}$ are the densities of PBHs and other types of matter respectively, and $\delta_{\rm M}$ denotes the density fluctuations of the latter. Treating the PBHs as point masses, we obtain

[TABLE]

where $m_{i}$ and ${\bf x}_{i}$ denote the masses and comoving positions of the PBHs. The motion of test particles is encoded in the action $m\int{\rm d}s$. A system of $N$ non-relativistic PBHs, $a{\rm d}x/{\rm d}t\ll 1$, obeys the action

[TABLE]

where ${\bf r}_{i}\equiv a{\bf x}_{i}$ is the proper distance. The last term accounts for the pairwise interaction of the PBHs and $\phi_{\rm ex}({\bf x})$ is an external potential describing the effect of the surrounding matter on the $N$-body system. The action of a two PBH system can then be approximated as

[TABLE]

where $\bf r$ denotes the separation of the PBH, ${\bf r}_{c}$ is the centre of mass of the 2-body system, $M\equiv m_{1}+m_{2}$ and $\mu\equiv m_{1}m_{2}/M$ are the total and reduced mass of the binary, and ${\bf T}_{ij}\equiv\partial_{i}\partial_{j}\phi({\bf r}_{c})$ results from the expansion of the external potential around the centre of mass of the system. We assume that the centre of mass is stationary, so the time dependence of ${\bf T}$ can be estimated from the expansion only. The forces acting on the pair are summarised as follows

[TABLE]

where $\hat{\bf r}\equiv{\bf r}/r$ is the unit vector parallel to ${\bf r}$. The first three forces are radial (i.e. parallel to $\hat{\bf r}$), and the last term provides the torque that prevents the head on collision of the two PBHs. For the binary to form, the tidal forces are required to be much weaker than the gravitational attraction of the PBH pair, that is, the two bodies must be more tightly coupled to each other than to any other PBH.

Decoupling from the expansion takes place when the second term in Eq. (2.7) starts to dominate over the first one. We define

[TABLE]

where $V(x)\equiv(4\pi/3)x^{3}$ denotes the comoving volume of comoving radius $x$. The quantity $\delta_{\rm b}-1$ can be interpreted as the effective matter overdensity generated by a PBH pair of a total mass $M$ and a comoving separation $x_{0}$. We are interested in initial PBH pairs with $\delta_{\rm b}\gg 1$, as they produce tightly bound PBH binaries that may merge within the age of the universe. In the radiation dominated epoch such perturbations collapse roughly when $\rho_{\rm R}a^{-4}\approx\delta_{\rm b}\rho_{\rm M}a^{-3}$. So, we define the scale

[TABLE]

as an approximate estimate of decoupling.

After decoupling, the separation of the PBH pair will stop growing and the tidal forces will be damped fast, ${\bf T}\propto a^{-3}$, as the universe expands. This is implied from Eq. (2.4) because, in the first approximation, during radiation domination the surrounding PBHs are assumed to follow the Hubble flow and the matter density perturbations are constant. This means that the potential (2.4) scales roughly as $\phi\propto a^{-1}$, which translates to ${\bf T}\propto a^{-3}$, as claimed.

The angular momentum ${\bf L}$ of the two body system vanishes initially. It is generated by the tidal torque,

[TABLE]

It is more convenient to work with the dimensionless angular momentum,

[TABLE]

where $r_{a}$ is the semi-major axis of the PBH binary. The orbital eccentricity $e$ of the binary is given by $e=\sqrt{1-j^{2}}$, and the coalescence time for eccentric orbits, $j\ll 1$, by [51]

[TABLE]

where $\eta\equiv\mu/M$ is the symmetric mass fraction. For non-eccentric orbits Eq. (2.12) may overestimate the coalescence time by at most a factor of $1.85$ [51].

Following Ref. [48] we estimate the effect of tidal torque perturbatively assuming that the decoupling from expansion takes place in the radiation dominated epoch and that the orbit of the binary remains eccentric, $j\ll 1$. The evolution of ${\bf j}$ can be evaluated perturbatively by first solving for purely radial motion using $\ddot{r}-r\ddot{a}/a+Mr^{-2}=0$ and then plugging the solution into Eq. (2.10). As the relevant binaries are formed in the radiation dominated epoch, we may ignore the contribution of the matter density $\rho_{M}$ to the Hubble parameter (2.1) in the first approximation. The equation of motion for the comoving separation $x\equiv r/a$ is then simply

[TABLE]

where the prime denotes derivation with respect to $\ln(a)$ and the initial conditions are given by $x(a_{0})=x_{0}$, $x^{\prime}(a_{0})=0$. We must assume that $a_{0}\ll a_{\rm dc}$. Eq. (2.13) is solved by $x(a)\approx x_{0}\chi(a/a_{\rm dc})$, where $\chi(y)$ is determined by $(y\partial_{y})^{2}\chi+y\chi^{-2}=0$ and the boundary conditions $\chi(y\to 0)=1$, $\chi^{\prime}(y\to 0)=0$. The function $\chi(y)$ is positive and oscillates with an amplitude that decreases asymptotically as $0.2/y$. So, the semi-major axis of the fully decoupled binary is

[TABLE]

The first root of $\chi(y)$, indicating the first close encounter of the pair, lies at $y\approx 0.54$, which translates to $a\approx 0.54a_{\rm dc}$, so the binary must decouple even earlier than our first estimate $a\approx a_{\rm dc}$ suggested. Although Eq. (2.13) ignores the contribution of matter to cosmic expansion, it yields a relatively good approximation even when $a_{\rm dc}\approx a_{\rm eq}$. Plugging the solution of Eq. (2.13) into Eq. (2.10), and, again, neglecting the contribution of matter to the expansion, gives

[TABLE]

where we used $\int^{\infty}_{0}{\rm d}y\,\chi(y)^{2}\approx 0.3$. Note that ${\bf T}a^{3}$ is constant under the assumptions made above. Because ${\bf T}$ is rapidly damped by the expansion, the binary acquires about 90% of its angular momentum during the first period.111The first period corresponds to the first root of $\chi(y)$ at $y\approx 0.54$, during which the angular momentum generated is proportional to $\int^{0.54}_{0}{\rm d}y\,\chi(y)^{2}\approx 0.26$. The total angular momentum corresponds to $\int^{\infty}_{0}{\rm d}y\,\chi(y)^{2}\approx 0.3$.

The coefficients in Eqs. (2.14) and (2.15) can be determined from the numerical simulations that we discuss in Sec. 3. So, to account for deviations from the simplified case, we will not fix them and define the coefficients $c_{a}$ and $c_{j}$ instead,

[TABLE]

For the above simplified scenario we have $c_{a}=0.1$ and $c_{j}=0.95$.

The orbit of the binary is now explicitly determined by the PBH masses $m_{1}$ and $m_{2}$, the initial comoving separation $x_{0}$ and the tidal torque ${\bf T}$, and the coalescence time (2.12) of the initial binary can be written as

[TABLE]

Of course, this prediction holds only when the orbit is not disturbed between the formation and merger of the binary. Note that the surrounding matter enters these expressions through the radiation density, which, by assumption, determines the Hubble parameter when the binary is formed, and the tidal torque, that implicitly accounts for the surrounding PBH and matter fluctuations. Both $x_{0}$ and ${\bf T}$ should be thought of as random variables that depend on the statistical properties of the initial PBH population. Next we will discuss how the distribution of masses and positions of the PBHs will determine the distribution of the orbital parameters derived above.

2.2 Distribution of initial binaries

Having detailed the properties of the PBH binary given the initial condition with known separation $x_{0}$, masses $m_{1}$, $m_{2}$ and ${\bf T}$, we now turn to discuss the distribution of initial conditions which will eventually determine the distribution of $j$ and the merger rate of the PBH binaries.

Consider a PBH pair with masses $m_{1}$, $m_{2}$ at a comoving separation $x_{0}$ so that they are the only PBHs in spherical volume of comoving radius $y$. This set-up is shown in the interior region of Fig. 1. The reason for forbidding surrounding PBHs closer than $y$ is to exclude initial configurations where the binary gets disrupted by surrounding PBHs shortly after formation. In such cases the perturbative estimate of the coalescence time will inevitably fail. For the sake of generality, we will not fix the value of $y$ in the general discussion. The aim of the following is to estimate the density of viable initial conditions and, from it, the distribution of coalescence times and the merger rate. The spatial PBH distribution is assumed to be Poisson throughout the paper.222Based on general arguments, the spatial distribution at small scales has been shown to be well approximated by the Poisson distribution [41]. It has been argued, however, that accounting for the two-point function of PBHs, $\xi_{\rm PBH}$, may affect the merger rate for wider mass functions [52, 53, 13]. In Ref. [53] this effect was estimated to be irrelevant for PBHs in the LIGO mass range. In addition, a rough estimate yields that for $\xi_{\rm PBH}\leqslant 1$ the contribution of the PBH two-point function to the merger rate is generally subleading to the direct contribution of the width of the mass function [13]. Only the latter will thus be considered in this paper. Formation of initially clustered PBH distributions, enabled by some more exotic PBH formation mechanisms, and the evolution of such clusters has been considered in Ref. [54].

In this case the comoving number density of configurations producing a binary is

[TABLE]

where ${\rm d}n(m)$ is the comoving number density of PBH in the mass range $(m,m+{\rm d}m)$, $\rho_{\rm DM}$ denotes the present DM energy density, $\bar{N}(y)\equiv nV(y)$ is the expected number of PBH in a spherical volume with radius $y$, and $n=\int{\rm d}n(m)$ is the PBH number density. The factor $1/2$ avoids overcounting.

The differential merger rate per unit time and comoving volume is then given by333In Ref. [55] it was argued, that the PBH binary merger rate is enhanced by a cascade of mergers in the early Universe. However, their application of the early binary formation mechanism beyond the first merger step is questionable as the PBHs may have peculiar velocities.

[TABLE]

where $j(\tau)$ is obtained from Eq. (2.12) and ${\rm d}P/{\rm d}j$ is the distribution of the dimensionless angular momentum for a given $y$ and $x_{0}$. We derive ${\rm d}P/{\rm d}j$ in the following section. The choice of $\bar{N}(y)$ should guarantee that the pair is initially a 2-body system and thus not gravitationally coupled to a third PBH. Numerical simulations presented in Sec. 3 indicate that initial configurations with a third PBH in the surrounding volume $V(y)\lesssim M/\rho_{\rm PBH}$ will produce binaries that are likely to be disrupted and will thus not contribute to the present merger rate.

2.2.1 Angular momenta

The distribution of eccentricities or, equivalently, the dimensionless angular momenta $j$, gets contributions from surrounding PBH and matter fluctuations. The dimensionless angular momentum from the surrounding PBHs, for which the tidal tensor reads ${\bf T}_{\rm PBH}a^{3}=\sum_{i}m_{i}({\bf 1}-3\hat{\bf x}_{i}\otimes\hat{\bf x}_{i})/x_{i}^{3}\,$, is given by the sum, ${\bf j}_{\rm PBH}\equiv\sum_{i}{\bf j}_{1}({\bf x}_{i},m_{i})$, of contributions from individual PBHs (see Eq. (2.15))

[TABLE]

where $\langle m\rangle=\rho_{\rm PBH}/n$ is the average PBH mass444In general, the average over masses is defined as $\langle X\rangle\equiv n^{-1}\int X{\rm d}n$. and we defined

[TABLE]

where we used $\Omega_{\rm M}/\Omega_{\rm DM}=1.2$. The quantity $j_{0}$ provides an order of magnitude estimate of the average dimensionless angular momentum, since, as most of the torque is likely generated by the PBH closest to the binary, for which $\bar{N}(x_{i})\approx 1$ and $m_{i}\approx\langle m\rangle$ on average, so the other terms in Eq. (2.20) will on average contribute an order one factor. Binary formation requires $\delta_{b}\gg 1$, thus we must assume $j_{0}\ll 1$. As we will see, this will be true for coalescence times less than the age of the universe.

We assume that the masses $m_{i}$ and positions ${\bf x}_{i}$ of the surrounding PBH and the matter density perturbations are statistically independent. The angular momentum, ${\bf j}={\bf j}_{\rm PBH}+{\bf j}_{\rm M}$, is then composed of a sum of independent variables, so its distribution is most conveniently estimated using the cumulant generating function $K({\bf k})\equiv\ln\left\langle e^{i{\bf k}\cdot{\bf j}}\right\rangle$,

[TABLE]

because of the additive property of cumulants.

We start by considering the first two cumulants of ${\bf j}$. Isotropy implies that

[TABLE]

where $\sigma_{j}^{2}\equiv\langle{\bf j}^{2}\rangle$ and the $\hat{\bf r}\otimes\hat{\bf r}$ term arises because ${\bf j}\perp\hat{\bf r}$ by construction (2.16). In general, the last feature implies that $K$ is a function of ${\bf k}_{\perp}\equiv{\bf k}\times\hat{\bf r}$, while isotropy further constrains it to be a function of $k_{\perp}$. Statistical independence implies that $\sigma_{j}^{2}=\sigma_{j,\rm PBH}^{2}+\sigma_{j,\rm M}^{2}$. The variance due to matter perturbations is obtained by applying Eq. (2.16), averaging over orientations of $\hat{\bf r}$ and using that, by Eq. (2.3), in Fourier space $a^{3}{\bf T}=\hat{\bf q}\otimes\hat{\bf q}\,4\pi\bar{\rho}_{\rm M}\delta_{\rm M}(q)$,555We assumed that $\bar{\rho}_{M}=\rho_{\rm M}$. This introduces a minor error as the contribution of $\sigma_{j,\rm M}$ becomes insignificant for $f_{\rm PBH}=1$, while the distinction is negligible for $f_{\rm PBH}\ll 1$.

[TABLE]

where $\sigma_{\rm M}^{2}\equiv(\Omega_{\rm M}/\Omega_{\rm DM})^{2}\left\langle\delta_{\rm M}^{2}\right\rangle$ is the rescaled variance of matter density perturbations at the time the binary is formed and $f_{\rm PBH}\equiv\rho_{\rm PBH}/\rho_{\rm DM}$ is the fraction of PBHs. Following Ref. [48] we will use the value $\left\langle\delta_{\rm M}^{2}\right\rangle=0.005$ for numerical estimates. We stress, however, that differences from this estimate can arise when PBH formation is accompanied by enhanced density perturbations at small scales, as is often the case.

A similar calculation yields, for the variance due to surrounding PBHs,

[TABLE]

where $\sigma_{m}\equiv\sqrt{\langle m^{2}\rangle-\langle m\rangle^{2}}$ is the width of the mass distribution. The new element in this computation is the average over positions of the surrounding PBHs. We first estimated it for a finite volume $V$ and then took the limit $V\to\infty$ by keeping the PBH number density, $n=N/V$, fixed. The contributions from each PBH are statistically independent and identical, so $\sigma_{j,{\rm PBH}}^{2}$ is just the contribution from a single PBH times the number of surrounding PBHs. The integral over positions has a lower bound because we exclude initial conditions where PBHs can be closer than $y$ to the centre of mass of the binary.

In conclusion, the variance of ${\bf j}$ is

[TABLE]

The limit $\bar{N}(y)\to 0$, where the variance diverges, is never realised as, in order not to disrupt the binary, the distance to the closest PBH must be larger than the separation of the binary, that is $y>x_{0}$. Therefore, since $j_{0}\approx\bar{N}(x_{0})\ll 1$ we have $\sigma_{j}^{2}\lesssim j_{0}$. A comparison of Eqs. (2.24) and (2.25) indicates that the variance due to the surrounding PBHs may be attributed to Gaussian matter perturbations with variance $(1+\sigma_{m}^{2}/\langle m\rangle^{2})f_{\rm PBH}^{2}/\bar{N}(y)$. As $\bar{N}(y)$ is the expected number of PBHs in the empty volume around the binary, this result is not surprising for a Poisson distribution of PBHs. Moreover, when $\bar{N}(y)\gg 1$, it is expected that the distribution of ${\bf j}$ is Gaussian with a width (2.26). This is indeed so, as we will confirm shortly.

The cumulant generating function of ${\bf j}$ decomposes as $K=K_{\rm PBH}+K_{\rm M}$. We will next calculate these separately. As matter density fluctuations are Gaussian only the first two cumulants contribute, thus $K_{\rm M}$ is simply

[TABLE]

To compute $K_{\rm PBH}$ we proceed in a similar manner as in the case of $\sigma_{j,\rm PBH}^{2}$ – starting with a finite volume and then taking the limit $V\to\infty$ by keeping $N/V$ fixed. This gives

[TABLE]

After plugging in Eq. (2.20) and averaging over PBH positions in the region $|{\bf x}|>y$ we get666Defining $z\equiv k_{\perp}j_{0}m/(\rho_{\rm PBH}V(y))$ and changing the integration variable to $u\equiv zV(y)/V(x)$, the spatial integration can be performed as follows

$\displaystyle-z\int\frac{{\rm d}^{2}\Omega}{4\pi}\,\int^{z}_{0}\frac{{\rm d}u}{u^{2}}\left[\exp\left(i3u(\hat{\bf k}_{\perp}\cdot\hat{\bf x})\,(\hat{\bf x}\cdot\hat{\bf r})\right)-1\right]$

$\displaystyle=-z\int^{z}_{0}\frac{{\rm d}u}{u^{2}}\int^{\pi}_{0}\int^{2\pi}_{0}\frac{{\rm d}\cos(\theta)\,{\rm d}\phi}{4\pi}\,\left[\exp\left(i\frac{3u}{2}\sin^{2}(\theta)\sin(2\phi)\right)-1\right]$

$\displaystyle=-z\int^{z}_{0}\frac{{\rm d}u}{u^{2}}\left[\frac{\pi}{2\sqrt{2}}J_{-\frac{1}{4}}\left(\frac{3u}{4}\right)J_{\frac{1}{4}}\left(\frac{3u}{4}\right)-1\right]={}_{1}F{}_{2}\left(-\frac{1}{2};\frac{3}{4},\frac{5}{4};-\frac{9z^{2}}{16}\right)-1.$

[TABLE]

where $F(z)={}_{1}F{}_{2}\left(-1/2;3/4,5/4;-9z^{2}/16\right)-1$ and ${}_{1}F{}_{2}$ is the generalised hypergeometric function. The function $F$ has the following properties

[TABLE]

which we list for later convenience. As a consistency check we find that the first two cumulants obtained from Eq. (2.29) match the earlier direct computation, i.e. $\left.-i\partial_{{\bf k}}K_{\rm PBH}\right|_{{\bf k}=0}=0$, $\left.(-i\partial_{{\bf k}})^{2}K_{\rm PBH}\right|_{{\bf k}=0}=\sigma_{j}^{2}$.

Both $K_{\rm PBH}$ and $K_{\rm M}$ depend only on $k_{\perp}$, so, as expected for an isotropic distribution of matter, it is sufficient to consider the distribution of $j$. From Eq. (2.22) we then obtain

[TABLE]

This distribution peaks at $j\lesssim j_{0}$. Since for PBH binaries merging today $j_{0}\ll 1$, this result is consistent with the assumption of eccentric orbits.

The logarithmic distribution $j{\rm d}P/{\rm d}j$ depends on $j$ only through the ratio $j/j_{0}$, so the characteristics of the binary, which affect only $j_{0}$, do not change the shape of the distribution but only shift it on the logarithmic scale. In particular, for monochromatic mass functions, ${\rm d}n(m)=n\delta(m-m_{c}){\rm d}m$, the shape of this distribution does not explicitly depend on the PBH mass. The distribution (2.31) for a monochromatic mass function is shown in Fig. 2 for $f_{\rm PBH}=1$ by the solid lines for different values of $\bar{N}(y)$.

In the limit $\bar{N}(y)\to 0$, $\sigma_{M}\ll f_{\rm PBH}$ we obtain a power law distribution with a break at $j_{0}$,

[TABLE]

This limiting case matches the result of Refs. [48, 50]. Interestingly, this result does not depend on the PBH mass function as it drops out of the integral in Eq. (2.31). Note, however, that the mass of the binary enters implicitly through $j_{0}$. As shown in Fig. 2, Eq. (2.32) approximates the distribution (2.31) relatively well for $\bar{N}(y)\lesssim 0.2$ if $f_{\rm PBH}\gtrsim\sigma_{\rm M}$. However, if $f_{\rm PBH}\lesssim\sigma_{\rm M}$ matter fluctuations can dominate over the Poisson fluctuations of PBH. In conclusion, the approximation (2.32) holds if the variance of the distribution (2.26) is dominated by PBHs, that is $\bar{N}(y)\ll f_{\rm PBH}^{2}/\sigma_{\rm M}^{2}$.

In the limit $\bar{N}(y)\to\infty$ we obtain a Gaussian distribution,

[TABLE]

where the width $\sigma_{j}$ is given by Eq. (2.26). The dotted lines in Fig. 2 show this limiting case. We see that, for monochromatic mass functions, the Gaussian distribution is a decent approximation already for $\bar{N}(y)=2$. The latter corresponds to the case where the binary lies in an underdense region prompting matter to initially move away from the binary. It also matches the conclusion based on the analysis of 3-body systems which states that the distance to the closest PBH should be at least of the order of the average distance between PBHs [45].

2.2.2 Merger rate

Assuming negligible disruption between formation and merger, the merger rate of the PBH binaries can now be obtained from Eq. (2.19). We take $\bar{N}(y)$ to be independent of $x_{0}$, that is, the binaries are expected to be disrupted if and only if initially there are surrounding PBH closer than $y$. This assumption is in good agreement with the numerical results of Sec. 3 in case $f_{\rm PBH}\ll 1$. The integrals over $x_{0}$ and $u$ can then be factorised by replacing $x_{0}$ with $v=uj_{0}/j(\tau)$. The primordial merger rate reads

[TABLE]

where

[TABLE]

is the rate in the limit $\bar{N}(y)\to 0$ and $\sigma_{M}/f_{\rm PBH}\to 0$. On the second line we have used the usual definition of a PBH mass function,

[TABLE]

that is normalised to unity, $\int\psi(m){\rm d}m=1$, and the numerical values $c_{a}=0.1$, $c_{j}=1$, $t_{0}=13.8\times 10^{9}\,{\rm yr}$. The suppression factor

[TABLE]

quantifies the contribution from matter density fluctuations and modifications due to the size of the empty region assumed around the pair. Note that, since $S$ does not depend on the coalescence time, the merger rate has an universal time dependence given by $\tau^{-\frac{34}{37}}$.

From the asymptotics of $F$ given in Eq. (2.30) we obtain the following limiting cases: First, in the limit $\bar{N}(y)\to 0$ the suppression factor reads

[TABLE]

where $U$ is the confluent hypergeometric function. Note that this is again independent of the shape of the mass function. The approximate merger rate reported in Ref. [48] corresponds to $S=\left(1+\sigma_{\rm M}^{2}/f_{\rm PBH}^{2}\right)^{-21/74}$ and underestimates the suppression factor (2.38) by at most a factor of 1.43 when $f_{\rm PBH}\ll\sigma_{\rm M}$.

Second, in the limit $\bar{N}(y)\to\infty$ we get

[TABLE]

For narrow mass functions, Eq. (2.39) is a good approximation already when $\bar{N}(y)\gtrsim 1$. The inequalities in Eq. (2.30) imply that the suppression factor is bounded between the above asymptotics,

[TABLE]

where in the last inequality can be saturated in the limit, $\sigma_{\rm M}/f_{\rm PBH}\to 0$. In particular, ${\rm d}R_{0}$ constitutes the upper bound on the merger rate of initial PBH binaries.

The dependence of the integrated merger rate on the PBH fraction is shown in Fig. 3. We assumed a lognormal mass function for the PBHs,

[TABLE]

where, following the notation of [23], $m_{c}$ denotes the peak mass of $m\psi(m)$, and $\sigma$ characterises the width of the mass spectrum. This class of mass functions has been found to provide a good representative of a large class of extended mass functions [56, 57, 58, 35, 59]. We stress, however, that the log-normal mass function is not universal as several effects may cause deviations from it [60, 61, 62, 63, 64, 65].

3 Evolution of PBH binaries in the early universe

The PBH binaries are the first gravitationally bound structures to be formed in the early universe. Subsequently, interactions with nearby PBHs may disrupt the binaries. To estimate the effect of the surrounding PBHs on the binary population in the early universe and to determine which initial configurations produce undisrupted PBH binaries, we conducted $N$-body simulations of PBHs in an expanding background, using a custom-written C++ code.

The simulations focus on the early evolution of PBH binaries, from their formation up to $a=3a_{\rm eq}$, which corresponds roughly to the time of recombination. The initial conditions of the central PBH pair are chosen such that a binary with an expected coalescence time equal the age of the universe is formed according to Eq. (2.17). The aim is to study the interaction of this pair with the surrounding PBHs, and to test the validity of the formation mechanism described in Sec. 2.1 and the analytical predictions for the distribution of the PBH binary orbital parameters given in Sec. 2.2.

3.1 Simulation set-up

The simulations are performed in physical coordinates and solve the equations of motion resulting from the action (2.5). The initial state consists of a central PBH pair that will form a binary and $N-2$ randomly distributed PBHs in a sphere around the binary as shown in Fig. 1. The individual particles in the simulation are subject to three different forces:

The gravitational attraction of the other $N-1$ particles is calculated in the Newtonian approximation without any regularisation at small distances. The individual contributions are summed. 2. 2.

The expansion of the universe is simulated by including a Hubble acceleration of the form $\ddot{\vec{r}}=\left(\dot{H}+H^{2}\right)\,\vec{r}$. 3. 3.

The gravitational attraction of the PBHs in the rest of the universe is approximated by the gravitational potential of a spherical underdensity inside otherwise homogeneously distributed matter. Equivalently, this corresponds to a sphere with uniform negative mass density equal to the positive mass density of the simulated black holes. The comoving radius and the total negative mass are kept constant.

We remark that the radius of the simulation is always much smaller than the Hubble radius, so the Newtonian approximation is justified. The simulation neglects the contributions of matter fluctuations as well as from the PBHs outside of the sphere. These are not relevant for $f_{\rm PBH}\gg\sigma_{\rm M}$, as the dominant forces acting on the central binary arise in that case from the closest PBHs surrounding it. For f_{\rm PBH}\mathrel{\raise 1.29167pt\hbox{<\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}\sigma_{\rm M}, in turn, these omissions underestimate the tidal forces. We also neglect the effect of emission of gravitational radiation. A binary that is loosing energy only by gravitational radiation satisfies $r_{a}\propto j^{-2}$ for $j\ll 1$ [51]. Eq. (2.12) implies the time dependence $j(t)=j_{0}/(1-t/\tau)$, where $\tau$ is the coalescence time at formation, so for binaries merging today the impact of gravitational radiation on the evolution of the orbit is of the order $10^{-4}$ within the first millions of years and can thus be safely neglected. We neglect baryons and, when $f_{\rm PBH}<1$, treat the surrounding DM as non-dynamical, accounting it only in the expansion rate.

The simulation uses adaptive, individual time steps for the particles and implements the leapfrog integration algorithm in the drift-kick-drift formulation. The size of the individual time steps is chosen such that the simulation (when ignoring the expansion of the universe) approximately conserves the total energy of the system. To determine the required precision of the energy conservation, we consider a typical binding energy of $E_{\mathrm{ref}}=10^{39}$ J (corresponding e.g. to two lightly bound BHs with masses of 30 $M_{\odot}$ circling each other with a distance of $\sim$4 mpc). The simulation is then required to conserve the total energy with a precision $\Delta E\leq 3\cdot 10^{-5}\,E_{\mathrm{ref}}$ in each global time step. We find that this level of precision allows for sufficiently accurate integration of binary orbits while keeping the total CPU time requirement manageable.

To allow for exact reproduction of the black hole orbits, a large dynamical range is required. The simulation reduces the individual, adaptive time steps until the required precision or the maximal number of substeps is reached. If the minimal substep of $\sim 0.63$ seconds cosmic time is not sufficient to reach the required precision, the code marks an exception in the output and continues with the minimal substep.

To ensure that the simulations finish in finite time, it is furthermore necessary to define a minimal distance of the BHs, which is chosen at $d_{\mathrm{min}}=3\cdot 10^{8}$ m. If two BHs approach each other more closely, they are merged into a single object while conserving the mass and momentum of the previous two objects. In practice, we find that this minimal distance allows for the simulation of binaries with lifetimes corresponding to a small fraction of the age of the universe, with the exact number depending on the other orbital parameters. We also tested that changes in the results are insignificant if $d_{\mathrm{min}}$ is decreased.

The simulations are started with an initial scale factor of $a_{\mathrm{init}}=10^{-3}\,a_{\mathrm{eq}}$. The initial set-up shown in Fig. 1 is achieved in the following way: First, the $N-2$ other particles in the simulation are randomly added to a spherical volume, such that the density of BHs in the spherical volume is equal to $f_{\mathrm{PBH}}$ times the total dark matter density. Then, the projection of the tidal force in a fixed direction from these $N-2$ particles is calculated for the centre of the sphere. Finally, the pair is added. The centre of mass of the pair is set to coincide with the centre of the spherical region and the vector joining the two PBH matches the fixed direction chosen before. The initial separation of the central BH binary is calculated by solving Eq. (2.17) for $x_{0}$ and demanding that $\tau$ equals the current age of the universe. Since the simulation is performed in non-comoving coordinates, each object is introduced with an initial velocity given by the Hubble expansion, i.e. $\vec{v}_{0}=H_{0}\,\vec{r}_{0}$, corresponding to a vanishing peculiar velocity. To guarantee that the initial conditions produce, on average, binaries that are expected to merge within the age of the universe, we used the values $c_{a}=0.129$ and $c_{j}=1.055$ found numerically from $f_{\rm PBH}=1$ simulations. After the central PBH pair has been added, the simulation is run until a scale factor of $a_{\mathrm{final}}=2.97\,a_{\mathrm{eq}}$ is reached, corresponding to a total simulation time of 377 kyr. The simulation therefore finishes after recombination.

We ran 70-body simulations with a monochromatic mass function for $f_{\rm PBH}=1$, $f_{\rm PBH}=0.1$ and $f_{\rm PBH}=0.01$ and with an extended mass function for $f_{\rm PBH}=0.1$. Each set of simulations consisted of 3000 simulations. The data used in the following analysis only includes simulations that reached $a=3a_{\rm eq}$ within a runtime of 10 days for $f_{\rm PBH}=1$ and 6 days for other simulations. We also excluded simulations where the central binary merged due to reaching the minimal distance set in the simulation. The latter makes up about 20% of the finished simulations in the case of monochromatic simulations with $f_{\rm PBH}=1$, 7% for $f_{\rm PBH}=0.1$ and 2% for $f_{\rm PBH}=0.01$, while the corresponding fraction is 2 -8% for simulations with an extended mass function. We checked that increasing the minimal distance or decreasing the runtime of the simulation does not affect our conclusions about undisrupted binaries. However, we see in Figs. 4 and  8, that the total number of used simulations (blue histogram) in the region where the central binary is expected to be disrupted remains below the expectation for Poisson distributed initial conditions (dashed line). This systematic effect can be attributed to the fact that $N$-body collisions are more likely to produce encounters where the distance briefly decreases below $d_{\mathrm{min}}=3\cdot 10^{8}$ m while their accurate simulation requires more computational resources making the simulation less likely to finish.

3.2 Properties of binaries merging within the age of the universe

The binary may be disrupted shortly after formation by a nearby PBH or a small cluster of PBHs, which may appear due to the Poissonian nature of the spatial distribution. We find that a single nearby PBH is the dominant source for disruption at the earliest times, especially when $f_{\rm PBH}\ll 1$. When $f_{\rm PBH}\approx 1$, however, the formation of bound systems of several PBHs can be observed already at matter-radiation equality, and, if the initial binary finds itself close to such clusters, it is very likely to be disrupted.

The initial binaries expected to merge today are highly eccentric, $j\ll 1$. If they interact with other objects, the eccentricity decreases and, since $\tau\propto j^{7}$, the coalescence time will on average be increased by several orders of magnitude. So, in the first approximation, all disrupted binaries are removed form the population of binaries merging today.

Consider collisions between the binary and the PBH initially closest to it. Let $x_{\rm NN}$ denote the initial comoving distance of the nearest neighbour. We would like to choose the size of the empty region around to binary as the smallest $x_{\rm NN}$ that does not lead to the early disruption of the binary. A rough estimate can be obtained by noting that, since the volume $V(x_{\rm NN})$ surrounding the PBH pair does not contain any PBH by construction, the central pair corresponds to an effective matter density fluctuation $\delta_{\rm NN}\approx(M-\rho_{\rm PBH}V(x_{\rm NN}))/(\rho_{\rm M}V(x_{\rm NN}))$. As such configurations are expected to collapse at $a\approx a_{\rm eq}/\delta_{\rm NN}$ we estimate that, at a given $a$, the binary has collided with its neighbour if $x_{\rm NN}<y$, where $y$ is given by

[TABLE]

We stress that this estimate is intended to provide a maximal value for $x_{NN}$ below which most initial binaries will be disrupted by the PBHs initially surrounding it. It does not account for any sources of later disruption, e.g. interactions with PBH clusters. For monochromatic mass functions and $a\to\infty$ we obtain $\bar{N}(y)=2$ consistent with the analysis of Ref. [45] which implied that $y$ should be of the order of the average distance between PBHs.

Using the simulations we can estimate if a binary with a nearest neighbour at distance $x_{\rm NN}$ gets disrupted. The results are summarised in Fig. 4. The numerical data is obtained from $2820$, $2423$, and $1382$ independent 70-body simulations with $f_{\rm PBH}=0.01$, $f_{\rm PBH}=0.1$ and $f_{\rm PBH}=1$, respectively. In all simulations, the central pair initially forms a binary with a coalescence time of roughly the age of the universe, as can be seen from Fig. 5. The blue region in Fig. 4 shows all initial configurations. As a consistency check we find that the initial distance of the nearest neighbour follows a Poisson distribution. Simulations in which at least one of the central BHs is in a binary system at $a=3a_{\rm eq}$ (upper panels) or $a=3a_{\rm eq}$ (lower panels), i.e. simulations where the initial binary may have swapped a BH, are shown by the yellow line, while the dashed yellow line shows simulations where the central BH remain bound to each other. More precisely, under the solid yellow line we included simulations where the energy of the system containing one of the central BHs and the BH closest to it was negative and we additionally imposed that the binding energy between the binary and its closest neighbour is less than 10% of the binary binding energy. For the dashed yellow line we only require that the total energy of the initial binary is negative. The green regions show simulations with central binaries with a coalescence time $\tau<10t_{0}$, which we use as the working definition for undisrupted binaries. This definition is justified by the fact that encounters with surrounding PBHs dominantly increase the coalescence time by several orders of magnitude, as can be clearly seen from the fist panel in Fig. 5.

Fig. 4 shows that the numerical result is consistent with the estimate (3.1) for the minimal value of $x_{\rm NN}$ and for its dependence on the scale factor. For $f_{\rm PBH}\ll 1$ the sharp transition at $y$ in Fig. 4 from almost all initial binaries being disrupted to almost all of them being non-perturbed confirms the underlying assumption in the computation of the rate (2.34) that $y$ is independent of $x_{0}$. For $f_{\rm PBH}=1$, however, half of the PBH with $x_{\rm NN}>y$ are disrupted already at $a=3a_{\rm eq}$. This implies an additional suppression factor on top of (2.37), because the latter accounts for disruption by the closest PBH only. Visual inspection of randomly selected simulations with $f_{\rm PBH}=1$ and $x_{NN}>y$ where the binary was disrupted indicates that the disruption is mostly, but not always, due to the interaction of the binary and a nearby $N$-body cluster as opposed to a direct collision with the nearest PBH. As the binaries will continue to be disrupted in the early clusters when $a>3a_{\rm eq}$, i.e. after the end of the simulation, it is expected that nearly all initial binaries will be disrupted within the age of the universe in case $f_{\rm PBH}\approx 1$. It is, however, not possible to draw definite conclusions from our numerical results beyond stating the need for a careful revision of early binary formation when $f_{\rm PBH}\approx 1$. We will return to the merger rate in that case in Sec. 3.5.

The expected coalescence time $\tau$ distribution shown in Fig. 5 for the $f_{\rm PBH}=0.1$ and $f_{\rm PBH}=0.01$ simulations is somewhat smaller than the age of the universe. This discrepancy is due to inaccuracies in the idealised analytic estimates (2.16) and (2.17) used to determine the initial separation $x_{0}$ for the central PBH pair. Notably, since $\tau\propto x_{0}^{37}$ small deviations in the initial separation estimate can lead to large deviations in $\tau$ from $t_{0}$. Based on test simulations we fixed the numerical parameters $c_{a}=0.129$ and $c_{j}=1.055$ to obtain merger times close to $t_{0}$ in the case $f_{\rm PBH}=1$. However, Fig. 5 shows that this choice overestimates $\tau$ for $f_{\rm PBH}=0.1$ and $f_{\rm PBH}=0.01$ by roughly a factor of 3 and 4, respectively. On the other hand, by Eq. (2.17) $\tau\propto c_{a}^{4}c_{j}^{7}$, so the obtained values $c_{a}=0.1$ and $c_{j}=0.95$ would have underestimated $\tau$ by a factor of 2. In all, in agreement with previous studies, we find that the analytic approach gives a reliable estimate of the coalescence time of the initial binary. Moreover, we stress that, because ${\rm d}R\propto c_{j}^{-21/37}c_{a}^{-12/37}$ by Eq. (2.35), the merger rate of initial PBH binaries is considerably less sensitive than $\tau$ to $\mathcal{O}(1)$ changes in $c_{j}$ and $c_{a}$.

Consider now the analytic estimate of the orbital parameters of binaries with a coalescence time $t_{0}$. The distribution (2.31) gives the conditional probability for a fixed initial separation $x_{0}$ and the size of the empty region $y$. Since the initial separation is distributed simply as $n\,{\rm d}V(x_{0})$, we obtain

[TABLE]

where we dropped overall factors independent of $j$ as they will be determined by properly normalising the distribution. In the last step we used ${\rm d}\tau/{\rm d}x_{0}\propto\tau/x_{0}$ and $j_{0}\propto x_{0}^{3}$. The initial separation $x_{0}$ is determined from the semimajor axis $r_{a}$ by Eq. (2.16) which in turn is fixed by the coalescence time (2.12) and the angular momentum $j$. This gives that $j_{0}\propto j^{-21/16}$. The distribution of $\ln j$, that is $j{\rm d}P/{\rm d}j$, depends on $j$ through the combination of $j/j_{0}$, thus we can define a similar characteristic angular momentum $j_{\tau}$ for fixed $\tau$ so that $\left.j{\rm d}P/{\rm d}j\right|_{\tau=t_{0}}$ will be a function of $j/j_{\tau}$ only, i.e. $j/j_{0}|_{\tau=t_{0}}=(j/j_{\tau})^{37/16}$. This gives777We used the analytic predictions $c_{a}=0.1$ and $c_{j}=0.95$. Since $j_{\tau}\propto c_{j}^{16/37}c_{a}^{-12/37}$, it is relatively insensitive to $\mathcal{O}(1)$ changes in these parameters.

[TABLE]

for binaries with total mass $M$ and mass asymmetry $\eta$. This quantity gives an order of magnitude estimate of the dimensionless angular momentum of initial binaries merging today.

In Fig. 6 we compare the distribution of angular momenta at $a=3a_{\rm eq}$ obtained numerically from the simulation and the corresponding analytic estimates based on Eqs. (2.31) and (3.2) with $\bar{N}(y)$ fixed using Eq. (3.1). The analytic prediction for unperturbed binaries is in good agreement with the numerical results (green). The $j$ distribution of perturbed binaries (yellow) is roughly uniform, i.e. $j{\rm d}P/{\rm d}j\propto j$.

Interestingly, the analytic prediction for the distribution of $j$ of the unperturbed binaries (solid line in Fig. 6) works well also when $f_{\rm PBH}=1$ although most of the initial binaries with $x_{NN}>y$ are disrupted. This observation can be used to shed some light on the cause of disruption of the binaries with $x_{NN}>y$. The choice (3.1) for the size of the empty region around the binary provides a rough estimate of the minimal $x_{\rm NN}$ for which the binary is almost certainly disrupted by the infall of the nearest PBH, but it does not predict what happens to the binaries for which $x_{\rm NN}>y$. As the analytic prediction for the distribution of $j$, which assumes that all initial binaries with $x_{NN}<y$ are disrupted while binaries with $x_{NN}>y$ remain undisrupted, matches the numerical one, the process disrupting the binaries with $x_{NN}>y$ must be statistically nearly independent of $j$ and thus also from $x_{\rm NN}$. This is because the nearest PBHs give the dominant contribution on the tidal torque which determines $j$. In conclusion, the disruption in simulations with $x_{\rm NN}>y$ is less likely to take place immediately after the surrounding PBHs have decoupled from expansion, but later, as they interact with clusters of PBHs. We have also observed this by visually studying individual simulations with $x_{\rm NN}>y$ that produce a disrupted central binary.

The distribution of semimajor axis is easily obtained by noting that $r_{a}\propto j^{-7/4}$ when the coalescence time is fixed, so $\left.j{\rm d}P/{\rm d}j\right|_{\tau=t_{0}}=(7/4)\left.r_{a}{\rm d}P/{\rm d}r_{a}\right|_{\tau=t_{0}}$. Analogously to the characteristic angular momentum $j_{\tau}$, we find the characteristic scale for the semimajor axis

[TABLE]

The semiminor axis $r_{b}=jr_{a}$ and the periapsis $r_{\rm per}\approx j^{2}r_{a}/2$ are therefore of the order of ${\rm au}$ and $10^{-2}{\rm au}$, respectively.

3.3 Properties of surrounding PBHs

Let us briefly examine the properties of the PBHs surrounding the central binary. By combing the data, excluding the central pair, from all simulations with a given $f_{\rm PBH}$, we estimate the two point function of the PBH spatial distribution of surrounding PBH. The result is shown in the left panel of Fig. 7. Initially, the two point function, shown by the grey dashed lines, is flat as expected for a Poisson distribution. The drop in $1+\xi(x)$ for large $x$ appears because of the finite size of the simulation. At $a=3a_{\rm eq}$ we observe a power law behaviour at small scales, $\xi(x)\approx\left(x/x_{c}\right)^{-\gamma}$, with $x_{c}=604(8){\rm pc}$ and $\gamma=2.38(1)$, for the data from $f_{\rm PBH}=1$ simulations. This shape does practically not depend on $f_{\rm PBH}$. Note that at $a=3a_{\rm eq}$ the flat region has almost disappeared for the $f_{\rm PBH}=1$ simulations and $1+\xi(x)$ has developed an extended tail at large $x$, as PBHs can be found outside the initial comoving volume. Thus, larger simulations are needed to probe the small scale structure beyond $a=3a_{\rm eq}$.

The distribution of velocities of the surrounding PBH is shown in the right panel of Fig. 7 for different values of $f_{\rm PBH}$. It has a peak which decreases as a power law on both sides and is followed by a roughly log-normal tail. The velocity dispersion in all simulations is approximately constant already after $0.1a_{\rm eq}$ and it takes the values 6 km/s, 1.5 km/s and 0.4 km/s for $f_{\rm PBH}=1$, $f_{\rm PBH}=0.1$ and $f_{\rm PBH}=0.01$, respectively.

3.4 Extended mass functions

To study the effect an extended mass function has on the survival of the binary, we simulated a simplified scenario where $f_{\rm PBH}=0.1$ and the PBH population consists of $3M_{\odot}$, $10M_{\odot}$ and $30M_{\odot}$ BH with mass density fractions $25\%$, $50\%$ and $25\%$, respectively. In detail, the simulation included 41, 25 and 4 surrounding PBHs of mass $3M_{\odot}$, $10M_{\odot}$ and $30M_{\odot}$ and a central pair with different combinations of these masses. The state of the central pair at $a=3a_{\rm eq}$ with a mass $(3+3)M_{\odot}$, $(10+10)M_{\odot}$, $(3+30)M_{\odot}$ and $(30+30)M_{\odot}$ is shown in Fig. 8. The histograms are based on 629, 636, 555 and 646 simulations, respectively. The definitions of the blue, yellow and green regions corresponding to all, bound and non-separated binaries are the same as for Fig. 4. Note that the dashed yellow line can exceed the solid one because the binding energy of the binary with its nearest neighbour is constrained only in the second case. In addition, the most massive, in this case $M=60M_{\odot}$, binaries can be bound to several light PBHs which provide only a small fraction of the binding energy of the system. The lighter PBHs will, in general, be eventually ejected from the system.

Fig. 8 shows that the estimate (3.1) for the disruption by the nearest neighbours works relatively well also for non-monochromatic mass functions with $f_{\rm PBH}\ll 1$. We checked that the angular momentum distribution follows the analytic prediction (2.31). The coalescence times of the binaries in these simulations were approximately $0.4t_{0}$, which is consistent with the monochromatic simulations (see the second panel of Fig. 5).

An important feature specific to extended mass functions is that, although heavy binaries are more easily disrupted, almost all of them will remain bound to each other. As a result, a population of less eccentric heavy binaries with a wide coalescence time distribution, peaked around $10^{10}t_{0}$, appears. On the other hand, the disruption rate is slightly increased when compared to the monochromatic case, as can also be seen from Fig. 9, because initial binaries containing light PBH are more easily disrupted.

3.5 Implications for the present merger rate

Let us now turn to the merger rate of initial PBH binaries in the late universe. To estimate which initial conditions will produce undisrupted binaries merging today we will rely on Eq. (3.1) for the choice of an appropriate value for $\bar{N}(y)$. The simulations indicate that this approach works when $f_{\rm PBH}\ll 1$ as for $f_{\rm PBH}\approx 1$ the disruption of initial binaries is not determined by the closest PBH, but their interaction with (or within) the compact $N$-body systems that form in the early universe. Extrapolating Eq. (3.1) to $a\gg 3a_{\rm eq}$ and considering only initial conditions that produce binaries expected to survive until the collapse of the first DM structures, i.e. until $a/a_{\rm eq}=1/\sigma_{M}$, we obtain

[TABLE]

By Eq.(2.37), the contribution of binaries much heavier than the average PBH mass will therefore be exponentially suppressed.

We stress that this is a rough approximation as it assumes that all initial PBH binaries with $x_{NN}\gtrsim y$ survive while all binaries with $x_{NN}\lesssim y$ are disrupted and will not merge within the age of the universe. This is certainly not the case when $f_{\rm PBH}\approx 1$. In the case of a monochromatic mass function, where $\bar{N}(y)\approx 2$ according to Eq. (3.5), we see from the left panel of Fig. 4 that about half of the binaries with $x_{NN}\gtrsim y$ are disrupted already after 0.4 Myr. The future of the remaining binaries will depend on their interaction with the surrounding PBH clusters and, given a half-life of the order of only 0.4 Myr, nearly all of these binaries are expected to be disrupted within the age of the universe. The fraction of unperturbed initial binaries will, of course, depend on the subsequent evolution of the PBH clusters, which is not accessible by our simulation. The disruption rate may decrease at later times when, for example, the dense clumps containing a few PBHs are dissolved within larger structures. By using Eq. (2.35) as a rough estimate, we note that when $f_{\rm PBH}\approx 1$ the initial binaries can produce a large enough merger rate to be consistent with LIGO even if only $0.1\%$ of the initial binaries remain unperturbed.

Most of the discussion focuses on a small fraction of all initial binaries – the ones that are expected to merge within the age of the universe. Even if nearly all binaries are disrupted, a large fraction of PBHs will still form binaries, as can be seen in Figs. 4 and 8. These binaries may contribute to the present merger rate.

The velocity dispersion estimates in Sec. 3.3 indicate that initial binaries merging within the age of the universe are hard – a hard binary is defined by having a binding energy, $E=m_{1}m_{2}/(2r_{a})$, that is larger than the average kinetic energy of the surrounding bodies. A well known result from studies of globular clusters is the Heggie-Hills law: hard binaries tend to get harder, while soft binaries get softer [66, 67]. The binding energy in binary-single body encounters is, on average, increased by an $\mathcal{O}(1)$ factor, which implies that the semimajor axis will decrease, but not by much. Binary-binary encounters, however, will generally lead to the ionisation of one of the binaries [68]. For wide mass spectra, the binary emerging from binary-single body collisions will generally comprise the heaviest PBH. Since the energy is roughly preserved, the final binary will expand so that $r_{a,\rm out}/r_{a,\rm in}\approx m_{\rm heaviest}/m_{\rm lightest}$.

The distribution of the relative change in the binding energy in different sets of simulations is shown Fig. 9. The data consists of simulations where the central binary was bound at $a=0.3a_{\rm eq}$. Hardened binaries are $\mathcal{O}(10)$ times more abundant than softened ones. In detail, the ratio of perturbed binaries (yellow) whose energy increased by at least $5\%$ to the number of perturbed binaries where the energy decreased by more than $5\%$ was 720/114, 493/40 and 622/92 in the panels of Fig. 9 counting from left to right, respectively. The frequency of hardened binaries in Fig. 9 drops exponentially fast as $E_{f}/E_{i}$ increases confirming that the binding energy will grow by a $\mathcal{O}(1)$ factor on average. In both sets of $f_{\rm PBH}=0.1$ simulations, the set of perturbed binaries (yellow) consist mainly of binaries that are expected to be perturbed at $a=3a_{\rm eq}$ by the estimate (3.1) as can be seen in Figs. 4 and 8.

Let us now attempt to derive a conservative estimate for the merger rate from perturbed initial binaries in the case $f_{\rm PBH}\approx 1$. For simplicity we will consider a monochromatic population of PBHs with mass $m_{\rm c}$. In that case, even if most of the initial binaries are disrupted, they will not be ionised. They will have roughly the same energy and semimajor axis as they had initially, but their eccentricity will be considerably decreased. Thus, in order to obtain a rough estimate for the merger rate of such binaries we may consider the following idealised scenario:

Most of the initial binaries remain bound. They may have exchanged a PBH. 2. 2.

The energy of the perturbed binaries is of the order of the initial binary. This means that we may estimate the semimajor axis from the initial conditions Eq. (2.14). To obtain the largest increase in coalescence time and thus the smallest rate we will assume that $r_{a}$ remains the same. 3. 3.

The eccentricity of the perturbed binaries is significantly decreased. Again, to obtain the smallest rate, we use $j=1$ as the final value. For final binaries with a coalescence time $\tau$, the initial comoving separation is then fixed by Eq. (2.12) and the previous assumption, $x_{0}\,\approx 10{\rm pc}\,(m_{\rm c}/M_{\odot})^{7/16}(\tau/t_{0})^{1/16}$. The horizon scale at PBH formation, $0.3{\rm pc}\,(m/M_{\odot})^{1/2}$, is about an order of magnitude smaller, thus such initial separations are viable. 4. 4.

The initial binary should be able to interact with the surrounding PBHs before it merges. Thus we assume that the coalescence time of the initial binary is larger than some time period $\tau_{c}$. When $r_{a}$ is fixed, the condition $\tau_{\rm init}>\tau_{c}$ implies that we need to consider initial binaries with $j_{\rm init}\geq\left(\tau_{c}/\tau\right)^{\frac{1}{7}}$. We do not constrain the distance of the closest PBH, thus we can use the initial $j$ distribution in the limit $y\to 0$, given by Eq. (2.32). Together with Eq. (2.21) we find that the probability of finding such angular momenta is

[TABLE]

for $\bar{N}(x_{0})\ll 1$. By the assumptions above $\bar{N}(x_{0})\approx 2\times 10^{-4}f_{\rm PBH}(m_{\rm c}/M_{\odot})^{5/16}$.

In all, with a final distribution of angular momenta peaked at $j=1$ and $\tau_{c}\approx 1\,\rm Myr$ we obtain from Eq. (2.19) that the present merger rate,888Formally, the evolution of the orbital parameters of initial binaries can expressed as

$P(j,r_{a})=\int{\rm d}j_{\rm init}{\rm d}r_{a,\rm init}\,P(j,r_{a}|j_{\rm init},r_{a,\rm init})P(j_{\rm init},r_{a,\rm init})\,,$

where according to the simplifying assumptions

$P(j,r_{a}|j_{\rm init},r_{a,\rm init})=\delta(j-1)\delta(r_{a}-r_{a,\rm init})\theta\left(j_{\rm init}-\left(\tau_{c}/\tau\right)^{\frac{1}{7}}\right).$

Note that the transition probability does not have to be normalised to unity since the number of binaries is not conserved. The rate can then be evaluated by replacing the distribution $P(j_{\rm init},r_{a,\rm init})$ with $P(j,r_{a})$ in Eq. (2.19).

[TABLE]

can be only slightly below the merger rate reported by the LIGO and Virgo collaborations, $9.7-101\,{\rm Gpc}^{-3}{\rm yr}^{-1}$ [6]. We stress that this is a very rough lower bound as more realistic scenarios will, on average, have final values of $j$ or $r_{a}$ which will increase this rate. Also the interactions with other surrounding matter will lead to an increase of $j$ and of the binding energy [49]. With $\tau_{c}=1\,\rm Myr$ we have that $j_{\rm initial}\geqslant 0.26$. The corresponding relatively high initial angular momentum will generally require the nearest neighbour to be initially much closer than the average separation between PBHs. Thus, it is expected that it couples to the binary early and, moreover, a more accurate description of the binary emerging from this configuration likely reduces to solving the full 3-body problem on an expanding background. In any case, the perturbative approach outlined in Sec. 2.1 has limited value in such situations. Another indicator that such binaries originate from early 3-body systems is the $f_{\rm PBH}^{3}$ dependence of the rate. To obtain an observable merger rate only a very small fraction of initial conditions have to produce binaries – the probability of the required initial 3-body configurations appear is roughly $\bar{N}(x_{0})^{2}\approx\mathcal{O}(10^{-8})$. The merger rate can decrease mainly by ionising these binaries. According to Fig. 9, this will take place with a small probability, and note that the binaries in our simulations are much softer than the ones contributing to the rate (3.7), i.e. they are easier to disrupt.

In conclusion, the perturbative estimate (2.34) becomes more accurate when $f_{\rm PBH}\ll 1$. For $f_{\rm PBH}\approx 1$ most of the initial binaries are disrupted999In contrast to Ref. [55], our results indicate that initial conditions for which the local PBH density is comparable or larger than the average density of DM likely lead to a decrease of the merger rate instead of a cascade of mergers. This is because such initial conditions tend to form compact $N$-body systems in which the initial binaries are rapidly disrupted., yet the disrupted binaries should still produce a merger rate consistent with LIGO. In both cases a more careful study of the formation and evolution of small scale structures and their interaction with the binaries is needed.

4 Phenomenology of PBH mergers

Now we turn to discuss the implications of the PBH mergers to the phenomenology of the PBH mergers. Using the results of the previous sections, we first study the temporal behaviour of the PBH merger rate, and by comparing it to the case of astrophysical BHs we confirm the earlier statements (see e.g. Ref. [24]) that the merger rates of primordial and astrophysical BHs are very different especially at large redshifts. Then, we perform a likelihood fit of the PBH mass function on the observed merger events, and finally derive potential constraints on the fraction of DM in PBHs from the rate of observed events and from the non-observation of the stochastic GW background.

4.1 Temporal behaviour of the merger rate

The $z$ dependence of the PBH merger rate (2.19) is shown by the red line in Fig. 10. For comparison, we show by the green line the $z$ dependence of the merger rate of astrophysical BHs. The properties of astrophysical BH binaries are assumed to inherit the properties of stellar binaries from which they are formed. In that respect, the $z$ dependence follows from the star formation rate and the delay from the formation of a stellar binary to merging of the BHs. The star formation rate is well approximated by the empirical formula [69, 8]

[TABLE]

and, as in Refs. [8, 70], we assume that the delay time distribution is $P_{d}(t)\propto t^{-1}$ for $t>50\,\rm{Myr}$ and zero otherwise. Then, the differential merger rate of BH binaries is given by

[TABLE]

where $\alpha$ and $\beta$ parametrise the mass dependence of the astrophysical binary formation. These are not relevant for the $z$ dependence that is given by the part in brackets.

The temporal behaviour of the merger rates for the primordial and astrophysical BHs are very different especially at large redshifts. Whereas the PBH merger rate monotonically increases as a function of $z$ due to the universal time dependence $R\propto\tau^{34/37}$ predicted by Eq. (2.35), the astrophysical BH merger rate starts to drop above $z\simeq 1$. LIGO is able to observe BH mergers only at low redshifts, z\mathrel{\raise 1.29167pt\hbox{<\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}0.4, but future GW observatories such as LISA or DECIGO can instead probe much higher redshifts, $z\simeq 10$. As pointed out e.g. in Ref. [24], the observations of high redshift BH mergers can reveal if a significant fraction of DM is in PBHs.

4.2 Likelihood fit to LIGO observations

Next we study the LIGO observations in the light of the derived merger rate. We start by finding a fit for the lognormal PBH mass function. We perform a maximum likelihood analysis considering the ten observed BH-BH merger events listed in Table III of Ref. [6]. The log-likelihood function reads

[TABLE]

where ${\rm d}P(m_{1},m_{2},z)\propto{\rm d}R(m_{1},m_{2},z){\rm d}V_{c}(z)$ is the differential probability of having a BH binary consisting of individual masses $m_{1}$ and $m_{2}$ merging at redshift $z$, and $V_{c}(z)$ is the comoving volume. The experimental uncertainties are accounted by $p_{j}(m_{j}|m)$ denoting the probability to observe a BH mass $m_{j}$ given that the BH has mass $m$, and $p_{j}(z_{j}|z)$, that is the probability to observe a BH-BH merger at redshift $z_{j}$ when it happens at redshift $z$. We take both $p_{j}$ to be Gaussian. The unit step function implements a detectability threshold based on the signal-to-noise ratio $\rho$ of the GW events.

Our estimate of the detectability of BH-BH mergers by LIGO is based on the signal-to-noise ratio of the GW events as in Refs. [71, 3]. We characterise the sensitivity of LIGO by the total strain noise, $S_{n}(f)$, of the detectors using the fit given in Ref. [15]. The actual strain sensitivity of the detectors fluctuates and has increased over time, whereas the total strain noise that we use is constant in time, and its amplitude is slightly above the best performance of the LIGO detectors. We therefore believe that our analysis may underestimate the detectability of the events, so a more careful treatment may produce slightly stronger constraints.

The signal-to-noise ratio is given by [71],

[TABLE]

where $\tilde{h}(\nu)$ is the Fourier transform of the signal [72, 73], that we compute using the results of Ref. [74]. For the detectability threshold we take $\rho>\rho_{c}=8$ [71, 3], implying roughly that mergers taking place at $z>0.4$ are not detectable.

Assuming that all LIGO BH-BH merger events are created by the early PBH binary formation mechanism, we can constrain the PBH mass function. The result of a maximum likelihood fit for a lognormal mass function (2.41) is shown in Fig. 11, where the red dashed and solid lines correspond to the $2\sigma$ and $3\sigma$ confidence levels, respectively, and the best fit, $m_{c}=20M_{\odot}$, $\sigma=0.6$, is indicated by the red dot. The best fit points to a relatively narrow mass functions, but the $3\sigma$ region extends to large widths. However, the constraints on the fraction of DM in PBHs disfavour very wide mass functions, as we will show in the next subsection.

4.3 Constraints on the PBH abundance

The prediction for the number of BH-BH merger events that LIGO should see within the period $\Delta t$ is

[TABLE]

During $\Delta t\simeq 165$ days, LIGO has observed $N_{\rm obs}=10$ BH-BH merger events [6]. We estimate the error of $N_{\rm obs}$ from Poisson statistics by taking the log-likelihood region $\ell-\ell_{\rm max}<-2$, which for Gaussian distribution would correspond to the $2\sigma$ confidence level. This gives $N_{\rm obs}=10\begin{subarray}{c}+7.7\\ -5.1\end{subarray}$, so by $2\sigma$ confidence level the number of observed events is smaller than $17.7$. This can be converted to an upper bound on the PBH abundance $f_{\rm PBH}$ for a given mass function using the merger rate estimate (2.34). In the case of lognormal mass function the resulting maximal value of $f_{\rm PBH}$ is shown by the colour coding in Fig. 11 in the $(m_{c},\sigma)$ -plane, and for $\sigma=0.6$ by the red solid line in Fig. 12. Moreover, the region between the red solid and red dashed lines in Fig. 12 corresponds $4.9<N<17.7$, that is, it corresponds to values of $f_{\rm PBH}$ for which the LIGO BH-BH merger event rate can be obtained at a $2\sigma$ confidence level. We remark that the PBH abundance inside the $2\sigma$ region in Fig. 11 is $f_{\rm PBH}<0.01$, confirming the validity of the perturbative merger rate estimate (2.34).

The non-observation of the stochastic GW background can further constrain the PBH abundance. The spectrum of the stochastic GW background from binary BH coalescences is [76, 77]

[TABLE]

where $\rho_{c}$ denotes the critical density, $\nu_{r}=(1+z)\nu$ is the redshifted source frequency, $H(z)$ is the Hubble parameter, and ${\rm d}E_{\rm GW}$ is the total GW energy in the frequency range $(\nu,\nu+{\rm d}\nu)$ emitted in the coalescence of BHs [72, 73, 74]. The factor $\theta(\rho_{c}-\rho(m_{1},m_{2},z))$ subtracts the contribution of events which can be observed individually. We find, however, that this factor can safely be neglected.

We calculate the stochastic GW background arising from coalescences of the early formed PBH binaries. By comparing its strength to the sensitivity of LIGO, we then constrain the fraction of DM in PBHs. The resulting upper bound arising from the latest LIGO observation run is shown by the black solid line in Fig. 12. The black dashed line shows instead the projected final sensitivity of LIGO. We see that in the range where LIGO is most sensitive the constraint from the observed merger rate is significantly stronger than the one from non-observation of the stochastic GW background. We remark that the constraint from the GW background is also weaker than the previous estimate given in Ref. [13]. This is due to the exclusion of the initial conditions where the initial binary is perturbed by its nearest neighbour. The constraint from the stochastic GW background is affected more strongly by the mass dependence of (3.5) and the resulting exponential suppression of the merger rate of heavy binaries. Note, however, that via the stochastic GW background it is possible to probe narrow mass functions with $m_{\rm c}<1M_{\odot}$ for which the individual mergers are not observable by LIGO. The stochastic GW background from early PBH binary mergers was also studied in [78].

The other constraints shown in Fig. 12 arising from lensing (EROS, Macho and SNe), accretion (Planck), observations of compact stellar systems (Eridanus II and Segue I) and wide binaries are evaluated for the lognormal mass function using the method introduced in Ref. [23]. The strongest of these (EROS and Segue I) exclude the hatched region above the black solid line in Fig. 11 for the fraction of DM in PBHs which would give $N=4.9$. Thus, the PBH origin of the LIGO merger events is excluded in that region as the required minimal PBH abundance is larger than allowed by EROS or Segue I.

We emphasise that whereas the early merger rate estimate for $f\simeq 1$ is subject to serious uncertainties as discussed in the previous sections, the merger rate (2.34) with the choice (3.5) works well for $f_{\rm PBH}\ll 1$. As the $f_{\rm PBH}\lesssim 0.1$ is already excluded by EROS or Segue I, we conclude that the observed merger rate significantly strengthens the constraint in the mass range $2-160M_{\odot}$. If the BHs observed by LIGO have a primordial origin, the fraction of DM in PBH and the peak mass of the distribution are contained within the region in Fig. 12 bound by the dashed red vertical lines and the red solid and dashed curves. These follow from the $2\sigma$ fit of the mass function and the merger rate, respectively.

5 Conclusions

We studied the formation and evolution of PBH binaries in the early universe, and revised the constraints from LIGO on the PBH abundance. Formed already before matter-radiation equality, the PBH binaries are the first gravitationally bound structures in the universe. We evaluated the characteristics of these binaries analytically allowing for extended PBH mass functions. The angular momentum was estimated perturbatively from the tidal force created by all surrounding matter, requiring that the PBH closest to the binary was sufficiently far. We then studied this scenario numerically using $N$-body simulations finding that the analytically obtained distributions of orbital characteristics of the binaries are in good agreement with the numerical result. We found that the angular momentum distribution is approximately Gaussian for the undisrupted PBH binaries.

The initially formed binary may be disrupted shortly after formation by a nearby PBH or a small cluster of PBHs. The latter may appear due to the Poissonian fluctuations, and we find that when $f_{\rm PBH}=1$ a large fraction of the initial binaries expected to merge within the age of the universe were absorbed into these clusters already before recombination. If $f_{\rm PBH}\ll 1$ the $N>2$ body systems form later, and in this case we found that single nearby PBHs are the dominant source of disruption. We estimated analytically which initial configurations avoid this disruption, and derived the suppression factor for the merger rate caused by this. We verified these estimates for the binary disruption by comparing them with the results of the $N$-body simulations.

We concluded that our analytical merger rate gives a good approximation for the merger rate of the initially formed PBH binaries in the case $f_{\rm PBH}\ll 1$, assuming that the disruption rate at later stages is not significant. For $f_{\rm PBH}=1$ the present merger rate estimate is greatly suppressed already at recombination, and our analytical results for the merger rate fail. Nevertheless, as most initially close PBH pairs will form binaries that are unlikely to be ionised, we argue that the population of perturbed binaries can still generate a merger rate consistent with LIGO. Moreover, since our results depend on the local PBH density, we conclude that, in case the PBH are initially clustered, the early binaries are more likely to be disrupted and the merger rate is therefore suppressed.

We used our merger rate results to re-evaluate the constraints on the PBH abundance from LIGO observations. We calculated the constraints from both the observed merger rate and the non-observation of the stochastic GW background. We found that the former gives the strongest constraint on the PBH abundance in the mass range $2-160M_{\odot}$. We also performed a likelihood fit for the lognormal PBH mass function on the 10 observed BH-BH merger events finding a best fit at $m_{c}=20M_{\odot}$ and $\sigma=0.6$. In that case the fraction of DM in PBHs needed to obtain the observed merger rate is $f_{\rm PBH}\simeq 0.002$.

Acknowledgments

The authors thank Yacine Ali-Haïmoud for useful discussions. This work was supported by the grants IUT23-6, PUT799, by EU through the ERDF CoE program grant TK133, and by the Estonian Research Council via the Mobilitas Plus grant MOBTT5. VV was supported by the United Kingdom STFC Grant ST/L000326/1.