The evolution of a primordial binary black hole due to interaction
with cold dark matter
and the formation rate of gravitational wave events

In what follows we neglect the influence of primordial black holes (PBHs) on the overall expansion of the Universe, assuming that it mainly consist of radiation and dust-like cold dark matter (CDM) particles. We are going to consider processes occuring at epochs before and after equipartition time, $t_{eq}\approx 23000yr$, defined by the usual condition that mass densities of matter and radiation are equal to each other at this time, and we denote these mass densities as $\rho_{eq}$. The scale factor, $s$, is assumed to be equal to one when $t=t_{eq}$. From the first Friedman equation we easily find a relation between $t_{eq}$ and $\rho_{eq}$

where $G$ is gravitational constant. Also, by integrating this equation we obtain

as well as useful relations between $s$ and its first and second derivatives with respect to the time $\tilde{t}$

where dot stands for the derivative with respect to $\tilde{t}$. In the asymptotic limits $\tilde{t}\rightarrow 0$ and $\tilde{t}\rightarrow\infty$ from (2) it follows that $s\approx{({\frac{4}{3}}(2-\sqrt{2})\tilde{t})}^{1/2}$ and $s\approx{((2-\sqrt{2})\tilde{t})}^{2/3}$, respectively.

In our analytical and numerical work we formally assume that a mass fraction of primordial black holes is $f<1$ and all black holes have the same masses $M=50M_{\odot}$. Our main results can be extended to other values of mass, therefore, we use dimensionless mass $M_{*}=\frac{M}{50M_{\odot}}$ in our expressions below as a parameter. A physical distance between two particular black holes is denoted by $R$, while the corresponding comoving distance, $R/s$, is $D$. It is convenient to introduce an average distance $\bar{R}=\sqrt[3]{\frac{3M}{4\pi\rho_{eq}}}$ and the associated dynamical time $t_{d}=\sqrt{\frac{{\bar{R}}^{3}}{GM}}$. From eq. (1) we have

Dimensionless physical and comoving distances are, respectively, $r=R/\bar{R}$ and $d=D/\bar{R}$. It turns out that the evolution of a separated binary black hole depends significantly on initial separation between binary components, $d_{0}=d(\tilde{t}\rightarrow 0)$. This quantity will be used to characterize different initial conditions for our numerical simulations.

When clusterization of CDM is not taken into account, one can use a simple ordinary differential equation to describe the evolution of $R$, its time derivative, and, accordingly, binary’s semimajor axis $a={\frac{GM}{|E|}}$, where $E$ is binary’s binding energy per unit of mass⁵⁵5Note that the energy is not conserved initially, when tidal forces determined by expansion of the Universe dominate over gravity of the binary. Nonetheless, we use the usual Keplerian definitions of binding energy and semimajor axis at all times., see e.g. [34]. In what follows we neglect orbital angular momentum of the binary, formally assuming that its eccentricity equals to one, see e.g. [26, 22, 27] for justification of this assumption. In this case the dynamical equation for the evolution of $R$ has the form

where the second term on l.h.s describes cosmological tidal forces, its explicit form follows from (3), the term on r.h.s is the usual Newtonian gravitational force per unit of mass, and we remind that dot stands for differentiation over $t/t_{eq}$⁶⁶6Note that one can also, in principle, include the dynamic friction term in this equation. However, it can be shown to be small during the radiation dominated stage due to the effect discussed in [54].. We introduce $r=R/\bar{R}$, use the definition of $\bar{R}$ and eq. (1) to bring (5) to the form

which should be solved subject to the condition that $r\rightarrow s(\tilde{t})d$ when $\tilde{t}\rightarrow 0$. Using the definitions of binding energy and semimajor axis as well as (1) we can easily express them in terms of the dimensionless distance $r$ and its time derivative

In general, equation (5) should be solved numerically, and its solutions could be parametrized by values of the initial comoving distance, $d_{0}$. These solutions have the following qualitative behaviour. Initially, at sufficiently small values of $\tilde{t}$, the pair is unbound and the distance follows the Hubble expansion, $r\approx s(\tilde{t})d$. However, the relative importance of gravitational attraction of the companion black hole grows with time, and, at a certain time $t_{bind}$ defined by the condition that $E(t_{bound})=0$, the pair gets bound, with its binding energy larger than zero. During the following stage the semimajor axis monotonically decreases up to some approximately constant value $a_{fin}$. At later times, in the absence of interaction with accreting dark matter, $a_{fin}$ is assumed to be changed only by emission of gravitational waves and we are going to compare our numerical results on the evolution of $a$ due to interaction with CDM with $a_{fin}$. The stage of a significant decrease of $a$ terminates when $r$ becomes formally equal to zero at first time, and we denote this moment of time as $t_{fin}$.

When $\tilde{t}$ is small it is easy to see that ${\frac{\ddot{s}}{s}}\approx{\frac{1}{4\tilde{t}}}$ on r.h.s. of (5). In this case (5) is invariant with respect to a linear change of variables $r\rightarrow\beta r_{n}$ and $\tilde{t}\rightarrow\beta r_{n}$ provided that $\gamma^{2}=\beta^{3}$. The remaining freedom in choosing $\beta$ and $\gamma$ can be used to make the asymptotic behaviour of $r_{n}(\tilde{t}_{n}\rightarrow 0)\propto{\tilde{t}}^{1/2}d_{0}$ independent of $d_{0}$. It turns out that in this case one have to use $\beta\propto d_{0}^{4}$ and $\gamma\propto d_{0}^{6}$. After such a choice of $\beta$ and $\gamma$ is made, being considered in terms of the new variables $r_{n}$ and $\tilde{t}_{n}$ both equation (5) and the initial condition to it do not depend on any parameters, and, therefore, the solution to this equation is uniquely fixed. The scaling of all quantities characterising the solution expressed in terms of the old variables, such as $s_{fin}$, $s_{bind}$ and $a_{fin}$ with $d_{0}$, should follow from the scaling of $\beta$ and $\gamma$. The scaling of $a_{fin}$ expressed in units of $\bar{R}$ should be the same as the distance $r$, and, accordingly, $a_{fin}/\bar{R}\propto d_{0}^{4}$ while the scalings of $s_{fin}$ and $s_{bind}$ should be the same as the scaling of the scale factor $\propto\sqrt{\tilde{t}}\propto d_{0}^{3}$. The coefficients of proportionality in front of these scaling relations can be chosen by comparing them to the results of a numerical solution of (5). In this way we obtain that

when $d_{0}$ is sufficiently small (see, also e.g. [27]).

We plot $a_{fin}$ as a function of $d_{0}$ in Fig. 1 and $s_{bound}=s(t=t_{bound})$ and $s_{fin}=s(t=t_{fin})$ in Fig 2 together with the corresponding approximate relations (8). One can see from this Fig. that the relations (8) are indeed quite close to the numerical solutions at $d_{0}<1$. Note that at large values of $d_{0}$ the effect of dark matter clusterization is very important and must be necessarily taken into account.

In order to simulate an interaction of PBHs and surrounding dark matter we use a direct summation N-body 4-th order Hermite code ph4 which is a part of Amuse suite [36, 37]. This code is MPI-parallel and we made it also OpenMP parallel. Since we simulate systems at times close to the matter-radiation equality, we implemented radiation as a background. Our system also follows Hubble expansion at the beginning. In cosmological N-body codes this is usually implemented by solving particle motion equations in the comoving coordinates and using periodic boundary conditions. Since we are more interested in the behaviour of a gravitationally bound part of the system, we use integration in physical coordinates. Periodic boundary conditions will introduce a distortion to the simulation of a spherically symmetric system, so we use vacuum boundary conditions and our system is a sphere with finite radius, large enough to avoid effects of density discontinuity at the sphere boundary on the central part of the system. The sphere is centered at the origin of coordinates. Thus, in order to implement the cosmological background in the simulations, we introduce a radial acceleration, acting on particle at a position $\mathbf{r}$:

We use a cubic grid for seeding the sphere with DM particles. Particles are assigned velocities according to the Hubble law. In case of a single black hole, it is placed close to the center with a tiny offset $<10^{-3}$ of the interparticle distance. If the BH is placed at (0, 0, 0) exactly, the trajectories of the closest particles pass directly through the BH and the timesteps of the simulation become very tiny, making it impossible to run. We do not take into account the absorption of DM particles by BHs in our simulations.

In case of a pair of BHs, we place them at comoving separation between BHs $d_{0}$, and assign tangential velocities assuming they are at apocenters of pure Keplerian orbits with eccentricity $e=0.8$.

We have made several test runs to check the simulation parameters and a series of production runs, which are summarized in a Table 1.

We took caution of choosing a proper value of gravitational softening length $\varepsilon$. On the one hand, the DM is a collisionless fluid, and it is desirable to have $\varepsilon$ larger than the mean interparticle distance. On the other hand, we are interested in the behaviour of a pair of BHs, for which scattering and dynamical friction is important, so $\epsilon$ should not be too high. Having this in mind, we have conducted tests to check how our results depend on softening. We use the same value of softening length for both DM and BH particles. For a single BH, we have run 5 simulations with 50k particles and the values of $\varepsilon=4.3\times 10^{-6}$, $8.5\times 10^{-5}$, $8.5\times 10^{-3}$, $2.5\times 10^{-2}$ and $5\times 10^{-2}$, respectively. The halo mass growth rate and density profile of these simulations is similar, while the anisotropy and distribution of angular momenta differ significantly. Both anisotropy and angular momentum are larger for simulations with higher $\varepsilon$ when we consider $\varepsilon=8.5\times 10^{-5}$, $8.5\times 10^{-3}$, $2.5\times 10^{-2}$ and $5\times 10^{-2}$ respectively. The simulation with $\varepsilon=4.3\times 10^{-6}$ demonstrates the same distributions of particle angular momentum and anisotropy as the simulation with $\varepsilon=8.5\times 10^{-5}$. As we discuss later on, angular momentum of infalling particles is a crucial parameter determining the evolution of a pair of BHs. From this analysis we conclude that values of $\varepsilon\leq 8.5\times 10^{-5}$ give stable results which are not affected by softening. For a pair of BHs we have conducted a similar investigation and found that the evolution of the pair semimajor axis is stable for $\varepsilon\leq 8.5\times 10^{-5}$. One should note that the mean interparticle distance at the halo center, estimated as $n^{-1/3}$, where $n$ is particle number density, always exceeds $10^{-3}$ for all our simulations.

Having a set of simulations with identical initial conditions and gravitational softening, but different number of DM particles, we investigate how well our results converge with the change of $N$, or particle mass $m_{p}$. As seen from Table 1, our simulations of a pair with $d_{0}=0.5$ cover almost 2 orders of magnitude in the particle mass, $m_{p}$.

It is well known that in a homogeneous expanding Universe a halo is formed around a single point mass (see, e.g., Gott 1975 [38] & Gunn 1977 [39]). The halo boundary can be characterized by a turn around radius: at this radius a shell of matter turns around its radial velocity from Hubble expansion to infall. As has been shown in [40], the radius and mass of the halo grow with time as:

We check that our simulations reproduce equations (10) and (11). At the same time we measure halo density profile from our numerical results. To do this, we construct a set of spherical shells centered at the black hole (or at the mass center of the pair, in the case of 2 BHs), each shell containing the same amount of particles $N_{\mathrm{shell}}=25$. The turn around shell is the furthest shell from the center which has mean radial velocity less than zero. The turn around radius is a radius at which the linearly interpolated mean shell velocity between the turn around shell and the next shell with positive radial velocity equals to zero. The turn around mass is the mass of all particles, except the black hole(s) inside the turn around radius.

In Fig. 3 we show the evolution of $r_{\mathrm{ta}}$ and $M_{\mathrm{ta}}$ in several simulations. It is clearly seen that in simulations with a single BH the evolution of these characteristics very well reproduces the equations (10) and (11). Simulations with a pair also show a perfect match until some specific time which is close to the first pericenter passage of the pair. This passage happens at $s=1.15$ for simulations with $d_{0}=0.95$ and at $s=0.075$ for simulations with $d_{0}=0.5$. After that time, the halo radius becomes somewhat smaller, while the mass drops significantly. The turn around mass evolution of a halo around a pair of black holes is close to $M_{\mathrm{ta}}\propto s^{3/4}$ at late times. This demonstrates that the pair throws away a significant amount of matter.

Another test of simulations is the density profile of dark matter halos. Theory predicts the formation of a power law profile $\rho\propto r^{-9/4}$ [40]. In Fig. 4 we demonstrate how our simulated halos follow this predicted law⁷⁷7Note that this result is different from that obtained in [29], who claimed that $\rho\propto r^{-3}$. It is seen that there is a good agreement, but the profile for a single black hole is flattening to $\rho\propto r^{-1.5}$ at $r<0.1$. For the pair of black holes this flattening is more pronounced, and also the whole profile has smaller values of density than the one corresponding to a single black hole at $r\lesssim 2$, which is of the order of the radius of influence of the binary. The effect of flattening of the density profile is clearly determined by domination of gravitational field of the black hole or the binary at small distances. The evolution of halo mass and the dark matter profiles show a remarkable stability with the change of the numerical resolution, i.e. the particle mass $m_{p}$ over almost two orders of magnitude.

For the pairs of BHs with a large initial distances, $d_{0}=1.5$ and $d_{0}=2$, there is some notable differences in halo evolution. The two halos around each BH grow in isolation for a long period of time, gaining mass much larger than the mass of BHs they host. After these two halos merge, the resulting single halo shows growth of mass and radius quite similar with the case of a single BH and it is well approximated as $M_{ta}\propto s^{0.9}$, see Fig. 5. However, the halo itself becomes severely anisotropic. We demonstrate it by measuring the major to minor axis ratio of a reduced gyration tensor, defined as:

where $x^{k}_{i}$ is $i$-th coordinate of $k$-th particle. We compute $G_{ij}$ for all particles with $r<r_{ta}$, excluding the BH particles, then we find eigenvalues of $G_{ij}$. The axis ratio is defined as the square root of the ratio of maximal and minimal eigenvalues of $G_{ij}$. In Fig. 5 we subtract 1 from this axis ratio, i.e. value of 0 corresponds to a spherically symmetric distribution. The halos around BH pairs that have merged earlier are much more spherical, as seen from the evolution of the $G_{ij}$ axis ratio in 2BH_50k_d0.95 simulation in Fig. 5. But the halo around a pair with $d_{0}=0.95$ is less spherical than a halo around a single BH at times $s<5$.

We characterize the pair of BHs at every simulation snapshot by two main parameters: the semimajor axis of the orbit $\tilde{a}=a/\bar{R}$, where $a$ is defined in (7) and the dimensionless angular momentum $l$, which is computed as

where $v_{t}$ is the relative tangential velocity of BHs (measured in units of $\bar{R}/t_{eq}$).

The results are shown in Fig. 6. The top row demonstrates the numerical convergence over almost two orders of magnitude in test particle mass in case of simulations with $d_{0}=0.5$. All simulations with the total number of particles from $5\times 10^{3}$ to $3\times 10^{5}$ show a steady decrease of $\tilde{a}$ after the first passage of pericenter, which happens at $s\approx 0.08$. Note that this value is very close to the corresponding value $s_{fin}$ obtained from the solution of eq. (6), see Fig. 2. One should note that due to the decrease of $\tilde{a}$ with time, at some moment there is less than 1 DM particle between the BHs which affects the pair evolution. For the lowest resolution simulation with $m_{p}=3\times 10^{-3}$ this happens at $s=0.2$, while for simulations with $m_{p}=3\times 10^{-4}$ and $m_{p}=9\times 10^{-5}$ this happens at $s=1$ and $s=2$, respectively. Our highest resolution simulation with $m_{p}=3.9\times 10^{-5}$ is stopped before the regime with less than 1 particle between BHs is reached.

In the bottom left panel of Fig. 6 it is seen that simulations with $d_{0}<1.5$ show a steady decrease of $\tilde{a}$ after the first pericenter passage. Simulations with $d_{0}\geq 1.5$ show a rapid drop of $\tilde{a}$ at the pericenter passage, determined by the merging of two massive halos.

The evolution of the dimensionless angular momentum, $l$, in Fig 6 shows two sets of lines with different initial $l$. One set corresponds to simulations which start at $s_{\mathrm{start}}=0.001$ and another one corresponds to $s_{\mathrm{start}}=0.1$. This is caused by the preparation of the initial conditions with fixed eccentricity described in Section III.1. Our results show an approximately tenfold drop of angular momentum during the evolution of the system.

At a first glance our numerical results for the evolution of $\tilde{a}$ and $l$ are at odds with the results of [31]. To investigate the cause of the difference, we first run with our code a simulation with the initial conditions from [31], which has the PBH mass 30 M_⊙, the initial semimajor axis $a_{i}=0.01$ pc and eccentricity $e_{i}=0.995$. These initial conditions have only 6.2 M_⊙ mass in the dark matter, which is about 1/10 of the binary mass. We obtain results in a fine agreement with that shown in [31], so the difference with our simulations is not due to the fact that we use a different N-body code. We also run a simulation with the same initial conditions, but with the additional force caused by the radiation background. This gives almost no effect on the evolution of the binary.

Next we check the influence of accretion of DM particles from the Hubble flow. In simulations of [31] there is no such accretion, but according to a simple estimate, the halo around a binary should grow by ten times during the time period shown in Fig. 6 of [31]. We run a test simulation with a Hubble sphere filled by DM particles, as described in Section III.1, but the sphere radius is artificially chosen small enough to end the accretion in the middle of the simulation (this simulation is labelled as 2BH_d0.5_small in Table 1). We show in Fig. 7 that when the accretion stops, the angular momentum stops decreasing. The same happens for the semimajor axis. So, we conclude that the accretion of DM particles is the main factor responsible for the late evolution of the PBH pair, after the formation of the common halo. The infalling DM particles have near radial trajectories and thus persistently refill loss cone of the binary, see also the following Section.

We can see from Fig. 6 when $d<0.95$ both semimajor axis and angular momentum experience a steady decrease. As discussed above, the most likely explanation for this behaviour is due to ’hardening’ of the binary due to its interaction with CDM particles having their periastron distances, $r_{p}$, order of or smaller than binary’s semimajor axis $a$. As discussed in e.g. [42], the contribution of a given particle to the evolution of both quantities is determined by the ratio $r_{p}/a$ (or, alternatively, the ratio of particle specific angular momentum to that of the binary), and only the particles with $r_{p}/a\sim 1$ can influence it significantly. In principle, the origin of particle’s angular momentum could be due to the tidal torque exerted by the binary on the particles. However, we checked that a value of angular momentum expected from this effect is too small to explain what is observed in simulations. A relatively large value of the angular momentum may be explained by torques acting on the particles from the side of non-spherical part of density distribution of the halos related to the observed halo’s anisotropy. The anisotropy itself may be originated from the action of the radial orbit instability (see e.g. [43]) with initial perturbations provided by the action of binary’s quadrupole gravitation field. In this picture the smaller is $d_{0}$ the smaller would be the anisotropy and this is indeed observed in the simulations. An analytical calculation of the corresponding distribution of angular momentum is a non-trivial problem, which is left for a future work. In this paper we use an average value of $r_{p}$, $\bar{r}_{p}$, evaluated numerically at the ’influence radius’ determined by the condition that halo’s mass inside this radius is equal to that of the binary. We show the dependency of $\bar{r}_{p}$ on $s$ in Fig. 8, for the case $d_{0}=0.95$.

Sesana et al. [42] showed numerically that, a contribution of a particle with a given value of $r_{p}$ to the hardening rate of semimajor axis is determined by a function $C(x)$ of $x=\sqrt{r_{p}/a}$, while the same contribution to the evolution of angular momentum is given by a similar function $B(x)$. We used their numerical results presented in their Figs 1 and 2, for an equal mass binary with eccentricity $e=0.9$ and fitted $C(x)$ by a simple function, $C(x)=2.647(1+x^{2})^{-2}$, which is in a quite good agreement with the numerical curve. $B(x)$ is a non-monotonic function of $x$, and, a simple fit of it does not appear to be possible. However, for our purposes we only need to know the integral $I(x)=\int_{0}^{x}dyyB(y)$, which can be fitted as $I(x)=2.86(1-\frac{1}{{(1+x^{3})}^{3/5}})$ remarkably well.

When $C(x)$ and $B(x)$ are known, it follows from the discussion in Sesana et al. [42] the time evolution of $a$ and $L$ are given by the simple expressions

where ${\it P}_{r_{p}}(r_{p})$ is a probability of finding the periaston radius with a value smaller than $r_{p}$. We model the corresponding probability density as a step function, $d{\it P}_{r_{p}}(r_{p})/dr_{p}=\frac{1}{2\bar{r}_{p}}\Theta(r_{2}-2\bar{r}_{p})$, where the factor $\frac{1}{2\bar{r}_{p}}$ in the front of the step function $\Theta(x)$ is chosen in such a way that the average value of $r_{p}$ coincides with $\bar{r}_{p}$. $\dot{M}$ in (14) is the mass flow through the sphere of influence of the pair (i.e. a sphere centered at the center of mass of the binary with the radius equal to the influence radius, $r_{inf}\approx 2$, see the discussion above). We assume that it is approximately the same as the mass flow through the turn around radius, $\dot{M}\approx M\frac{ds}{dt}$.

Combining all ingredients discussed above together, we can write down two evolution equations for $a$ and $L$:

As an initial condition we use values $a_{fin}$ and $s_{fin}$ obtained from the solution of eq. (6), $s_{fin}\approx 0.6$ and $a_{fin}(s_{fin})\approx 0.1$ Results of the solution of eqns (15) are shown in Figs 9 and 10. One can show from these Figs. that, indeed, at sufficiently large values of $s$ the fully numerical results are quite close to the solutions of (15).

The formation rate $R$ is determined, to a large extent, by a characteristic time scale of the evolution of $a$, $t_{GW}$, for a bound PBH pair having a large value of initial eccentricity $e\sim 1$

where $j=\sqrt{1-e^{2}}$, see e.g [27].

The pairs have a time to merge only when $t_{GW}<t_{H}$, where $t_{H}\approx 1.38\cdot 10^{10}yr$ is the present age of the Universe. We introduce $\tilde{t}_{GW}=t_{GW}/t_{H}$ and express it in terms of the quantities appropriate for the problem at hand. We obtain

where $\tilde{a}=a/\bar{R}$ and $M_{*}=M/50M_{\odot}$. From the condition $\tilde{t}_{GW}<1$ we obtain

In what follows we assume that both $\tilde{a}$ and $j$ (or, equivalently, the dimensionless angular momentum $l=\sqrt{\tilde{a}}j$) are certain functions of the initial comoving distance $d_{0}$, and a possible form of these functions follow from a particular scenario of the evolution of the pair before the process of its evolution due to emission of gravitational waves sets in. Note that it is different from the approach undertaken by AKK, who considered $j$ as a stochastic variable, but we are going to show that both approaches give quite similar results in the scenario when the effects of accreting CDM are neglected. When these functions are specified the condition $t=t_{GW}$ provide a certain value of $d_{0}$ (for given $f$ and $M_{*}$), $d_{*}(t)$, and only the pairs with $d<d_{*}(t)$ can merge before time $t$. When $t=t_{H}$ we set $d_{*}\equiv d_{*}(t_{H})$.

Assuming that spatial distribution of PBHs obey the Poisson statistics and introducing the variable $X=fd_{0}^{3}$ instead of $d_{0}$, one can easily show that probability density of distribution over a distance to the closest neighbor having some value of $X$, $\frac{d{\it P}}{dX}$ is $\frac{d{\it P}}{dX}=e^{-X}$. Accordingly, the probability to have a pair with $d_{0}<d_{*}$, ${\it P}_{*}$, is ${\frac{1}{2}}X_{*}e^{-X_{*}}$, where the factor ${\frac{1}{2}}$ is inserted to avoid double counting of the components, and $X_{*}=fd_{*}^{3}$, see also AKK. Taking into account that we expect $X_{*}\ll 1$ we can set ${\it P}_{*}=X_{*}$.

In order to estimate the merger rate we need to know the present day number density, $n_{pbh}$. This follows from our definition of the PBH mass fraction $f$ and the average radius $\bar{R}$

where $s_{PT}\approx 5.88\cdot 10^{3}$ is the present day value of the scale factor and we remind that we normalize it by the condition that it is equal to one at the equipartition time. A number density of PBHs, which have been merged before the present time, $n_{merge}$, follows from (19) and the definition of our probability ${\it P}_{*}=X_{*}$

We obtain the merger rate differentiating $X$ in (20) and remembering that $X$ is a function of time due to the condition $t=t_{GM}$ and our assumption that $\tilde{a}$ and $j$ entering $t_{GW}$ are functions of $d$. We have

When the influence of accreting CDM on the evolution of semimajor axis and angular momentum of the binary is neglected the problem can be analyzed with the help of solutions to equation (5) and calculation of some integrals related to those solutions, see e.g. AKK, [27]. We closely follow AKK in this Paper for an estimate of the event rate in this case and use this estimate as a ’reference value’ for the problem, which takes into account the effect of clusterization.

When the accretion is neglected one can use $a_{fin}$ represented in Fig. 1 as the corresponding value of semimajor axis entering (17). In the same approach $j$ does change significantly at time $t>t_{eq}$ and its value is assumed to be determined by two factors - 1) gravitational torques on the pair from the side of other black holes, and 2) gravitational torques determined by the cosmological density perturbations. In the first case a distribution over possible values $j$ is given by equation (19) of AKK, which can be used to calculate a mean value of $j$. However, the corresponding integral is logarithmically divergent, which is related to the fact that nonphysical values of $j>1$ are allowed. We set, accordingly, the upper integration limit $j=1$ there, thus obtaining the average value

where we remind that $f$ is the ratio of PBH mass density to CDM mass density and it is assumed that $f\ll 1$. When $f\geq 10^{-3}$ $\Lambda\sim 10$. The average value of $j$ due to the influence of cosmological perturbations, $j_{cp}$ is given by eq. (21) of AKK. Being rewritten with the help of our notations it has the form

where we assume that $\sigma_{eq}=5\cdot 10^{-3}$ in eq. (21) of AKK. In order to find a typical value of $t_{GW}$ in eq. (16) we use $j=j_{bh}$ when $j_{bh}>j_{cp}$ and $j=j_{cp}$ otherwise.

When the result of solution of eq. (5) and $j$ given by the expressions above are substituted in (18) the condition $\tilde{a}_{crit}=\tilde{a}_{fin}$ can be used to express $d_{*}$ as a function of $f$, for a fixed value of $M_{*}$.

We show this function in Fig. 11 for the case $M_{*}=1$. Note that the saw-like behaviour of the curve is simply due to relatively small numbers of our grid points used to represent our quantity $d_{0}$ when solving (5). This is clearly nonphysical, but, since it does not influence any our results it is left unchanged. Also note that $d_{*}$ is constant at small values of $f$. This is because when $f$ is small $j=j_{cp}$ which does not depend on $f$, see eq. (23)

From eq.(8) it follows that $\tilde{a}_{fin}\propto d_{*}^{4}$ when $d_{*}$ is not too large, while from the discussion above we see that we may approximate the dependence of $j$ on $d_{8}$ as $j\propto d^{3}$. Taking into account that $X_{*}=fd_{*}^{3}$ we see from eq. (16) that $t_{GW}\propto X_{*}^{37/3}$ and, accordingly, $\frac{t_{GW}}{X_{*}}\frac{\partial X_{*}}{\partial t_{GW}}=\frac{3}{37}$ in (21) when accretion is neglected. We use $\frac{t_{GW}}{X_{*}}\frac{\partial X_{*}}{\partial t_{GW}}=\frac{3}{37}$ in eq. (21), thus obtaining the event rate in the model case of the absence of accretion, which is shown in Fig. 14 below together with an estimate of the event rate, which takes into account the effects of CDM accretion. From this Fig. it is seen that, in the case when accretion is neglected, in order to have $R$ in the interval 10-100 $Gpc^{-3}yr^{-1}$ the mass fraction $f$ should be in the interval $2.5\cdot 10^{-3}-1.2\cdot 10^{-2}$. It agrees perfectly with the corresponding result of AKK ⁸⁸8Note that the value of $M_{*}$ in their Fig. 5, which is the closest to our value $M_{*}=1$, is $M_{*}=3/5$. This shouldn’t, however, influence the agreement in any significant way, see their Fig. 5. From Fig. 11 it also follows that when $f$ is approximately in the range $10^{-3}-10^{-2}$, $d_{*}$ should approximately be in the range $0.5-1$. This justifies the range of $d_{0}$ chosen for our numerical simulations.