Influence of early dark matter haloes on the primordial black holes merger rate

ABSTRACT

Primordial black hole (PBH) binaries forming in the early Universe may contribute to the merger events observed by the LIGO–Virgo–KAGRA collaborations. Moreover, the inferred merger rate constraints the fraction of PBH with masses |$m \sim 10 \, {\rm M}_{\odot }$| in the dark matter (DM) to f_PBH ≲ 10⁻³. This constraint assumes that after the formation of PBH binaries, they do not get destroyed or their parameters are not perturbed until the merger. However, PBHs themselves contribute to the formation of early DM structures in which the interactions between PBHs take place actively. This leads to the fact that the binaries can be perturbed in such a way that their lifetime becomes longer than the Hubble time t_H. In this work, we consider the effect of the initial spatial Poisson distribution of PBHs on the structure formation at the high redshifts z ≳ 10. Next, we explore the evolution of such haloes due to the interaction of PBHs with each other and with DM particles. We show that the early haloes evolve on time-scales much shorter than the age of the Universe. Furthermore, for fractions of PBHs f_PBH < 1, the internal dynamics of a halo is significantly accelerated due to the dynamical friction of PBHs against DM particles. As a result, a significant fraction of binaries will be perturbed in such structures, and the gravitational waves constraints on PBHs with masses |$m \sim 10 \, {\rm M}_{\odot }$| can be weakened to f_PBH ∼ 0.1.

black hole mergers, gravitational waves, dark matter

1 INTRODUCTION

Merging primordial black holes (PBHs) binaries are actively considered as sources of gravitational wave events observed by the LIGO–Virgo–KAGRA collaborations (Bird et al. 2016; Blinnikov et al. 2016; Sasaki et al. 2016, 2018; Ali-Haïmoud, Kovetz & Kamionkowski 2017; Clesse & García-Bellido 2017; Raidal, Vaskonen & Veermäe 2017; Raidal et al. 2019; Dolgov et al. 2020; Hütsi et al. 2021). However, in addition to gravitational waves, PBHs can explain a number of other modern problems of cosmology and astrophysics (Clesse & García-Bellido 2018; Carr & Kuhnel 2022; Carr et al. 2023). Among them are the observation of supermassive black holes (quasars) at high redshifts z ≳ 6 (Banados et al. 2018; Inayoshi, Visbal & Haiman 2020; Yang et al. 2020; Wang et al. 2021), possibility of explaining some part or even all of the dark matter (DM) depending on the PBH mass (Green & Kavanagh 2021; Carr & Kuhnel 2022; Carr et al. 2023). Galaxies discovered by the JWST telescope in the very young Universe (Castellano et al. 2022; Finkelstein et al. 2022; Atek et al. 2023; Ferrara, Pallottini & Dayal 2023; Labbé et al. 2023) can be also explained in terms of PBH (Liu & Bromm 2022; Dolgov 2023; Hütsi et al. 2023).

Despite the fact that PBHs can potentially be part of gravitational wave events, the observed merger rate imposes the strongest constraints on the fraction of PBHs in DM f_PBH ≲ 10⁻³ at tens of solar masses (Sasaki et al. 2016; Ali-Haïmoud et al. 2017; Kavanagh, Gaggero & Bertone 2018; Pilipenko, Tkachev & Ivanov 2022; Jangra, Kavanagh & Diego 2023; Postnov & Mitichkin 2023). This constraint suggests that in the early Universe, a pair of PBHs decouple from the Hubble flow and form a binary system (Nakamura et al. 1997; Ioka et al. 1998). Further, these binaries gradually undergo the stage of inspiral due to the emission of gravitational waves and eventually merge. However, from the moment of formation to the merger, the parameters of the binary system can be significantly perturbed so that the lifetime of the binary exceeds the Hubble time t_H, as a result of which the merger rate can be reduced (De Luca et al. 2020; Jedamzik 2020; Vaskonen & Veermäe 2020). In addition, if PBHs are initially strongly clustered, the constraint on gravitational waves can be relaxed to f_PBH ∼ 1 (Eroshenko & Stasenko 2023).

The discrete nature of PBHs leads to the early formation of structures due to their initial Poisson distribution in space (Meszaros 1975; Afshordi, McDonald & Spergel 2003; Inman & Ali-Haïmoud 2019). The forming haloes have a much larger mass than predicted in purely adiabatic inflationary fluctuations with nearly scale invariant power spectrum. The evolution of such DM haloes containing some amount of PBHs is, in many respects, similar to that of globular star cluster; in particular, for very early haloes, the phenomenon of core collapse takes place (see for review Spitzer 1987), time of which is much less than the age of the Universe t_H. If a PBH binary finds itself in a halo that experiences core collapse, then it is very likely that it will not contribute to the PBH merger rate (Vaskonen & Veermäe 2020). This is due to the fact that in the process of halo evolution, the probability of perturbing a pair of PBHs increases significantly. However, the analysis of Vaskonen & Veermäe (2020) is somewhat simplified, because it does not take into account the halo density profile and an influence of DM particles on the dynamics of PBHs in the halo for the case f_PBH < 1. PBHs will experience dynamical friction against DM particles, which leads to a decrease in the core collapse time t_cc. The main aim of this article is to find the fraction of binaries that will be in haloes that experience core collapse in a time less than the age of the Universe. In this paper, to study the listed effects, the kinetic Fokker–Planck equation is solved. This equation describes the evolution of a self-gravitating system and is often used to study the dynamics of globular star clusters and galactic nuclei (Vasiliev 2017).

In this work, we generally consider the two-component DM, which consists of a PBH with mass |$m = 10 \, {\rm M}_{\odot }$| (we assume a monochromatic mass spectrum) and some massive unknown particles that gravitationally interact with the PBHs. First, we show how PBHs influence the formation of early DM haloes at high redshifts (z ≳ 10) due to Poisson noise. Then, the evolution of such haloes is considered by solving the Fokker–Planck equation. We show that the time-scale of the core collapse depends significantly on the PBH fraction in the DM composition f_PBH. Next, we qualitatively consider the dynamics of PBH binaries that form in the early Universe and show that their parameters can be significantly perturbed during halo evolution. Eventually, we show how the modern PBH merger rate is changing and conclude that the constraints can be relaxed to f_PBH ∼ 0.1.

2 EARLY DM HALOES

The discrete nature of PBHs induces Poisson fluctuations with amplitude |$\delta \sim f_{PBH}/\sqrt{N}$|⁠, where N is the average number of PBHs in the considered volume (Carr & Silk 2018). Note that these fluctuations are isocurvature perturbations and can have a much larger amplitude than the adiabatic inflationary fluctuations δ_ad ∼ 10⁻⁵ ÷ 10⁻⁴. As a result, on small scales, the matter power spectrum is modified as follows:

$$\begin{eqnarray} P(k) = P_{PBH} + T^2_{ad}(k) P_{ad}(k). \end{eqnarray}$$

(1)

Here, P_ad = Akⁿ with n = 0.96¹ (Aghanim et al. 2020) and T_ad are the matter power spectrum and the transfer function for adiabatic inflationary fluctuations (Mo & White 2002) and P_PBH is the PBHs contribution (Hütsi, Raidal & Veermäe 2019; Mena et al. 2019):

$$\begin{eqnarray} P_{PBH} = \frac{f^2_{PBH}}{n_{PBH}} T^2_{iso}, \end{eqnarray}$$

(2)

where n_PBH is the comoving PBH number density and T_iso is the transfer function for isocurvature perturbations, which describes the linear growth of fluctuations (Peacock 1999; Gorbunov & Rubakov 2011)

$$\begin{eqnarray} T_{iso} = \frac{3}{2} (1 + z_{eq}) g(0), \end{eqnarray}$$

(3)

where g(0) ≈ 0.74 reflects the suppression of density fluctuation growth due to the Λ term (Mo & White 2002) and z_eq ≈ 3402 is the redshift of matter–radiation equality. Thus, the power spectrum induced by the Poisson noise of PBHs has the form

$$\begin{eqnarray} P_{PBH} = \frac{9}{4} g^2(0) (1 + z_{eq})^2 \frac{m f_{PBH}}{\Omega _{DM} \rho _{crit}}, \end{eqnarray}$$

(4)

where |$\rho _{crit} \approx 126 \, {\rm M}_{\odot }$| kpc⁻³ is the critical density.

The variance of fluctuations on the mass scale M = 4πR³ρ_M/3 as the function of the redshift z is given by standard expression

$$\begin{eqnarray} \sigma ^2_M(R,z) = D^2(z)\int \frac{dk}{2 \pi ^2} \, k^2 P(k) W^2 (kR), \end{eqnarray}$$

(5)

where function D(z) = g(z)/[g(0)(1 + z)] describes the growth of fluctuations and W(kR) is the window function in the Fourier space (Top-hat spherical filter)

$$\begin{eqnarray} W(kR) = \frac{3 \, (\sin kR-kR \cos kR)}{(kR)^3}. \end{eqnarray}$$

(6)

On the scales where the main contribution to the power spectrum is due to PBHs (particularly at high redshifts z), the variance of fluctuations will be

$$\begin{eqnarray} \sigma ^2_M(M,z) = \frac{9 \, m f_{PBH}}{4 \, M} \left(\frac{1 + z_{eq}}{1 + z} \right)^{2}. \end{eqnarray}$$

(7)

Here, for simplicity, the contribution of baryons is neglected, i.e. it is assumed that Ω_M = Ω_DM and the value of the integral from the square window function

$$\begin{eqnarray} \int _{0}^{\infty } dx \, x^2 W^2(x) = \frac{3}{2} \pi \end{eqnarray}$$

(8)

is used.

Let us define the characteristic halo mass (such haloes are also called 1σ haloes) as σ_M(M*) = δ_c, where δ_c = 1.69 is the critical fluctuation value for spherical collapse calculated in linear theory. Fig. 1 shows the redshift dependencies of the characteristic halo mass calculated for different PBH contributions to DM f_PBH, where the dotted lines are obtained with the formula

$The characteristic halo mass as a function of the redshift for different fractions of PBHs. The solid lines correspond to the solution of the equation σM(M, z) = δc, where σM is defined by equation (5). Dashed lines are given by the approximate equation (9).$

Figure 1.

The characteristic halo mass as a function of the redshift for different fractions of PBHs. The solid lines correspond to the solution of the equation σ_M(M, z) = δ_c, where σ_M is defined by equation (5). Dashed lines are given by the approximate equation (9).

Open in new tab Download slide

$$\begin{eqnarray} M^{*} = \frac{9 \, m f_{PBH} }{4 \, \delta _c^2} \left(\frac{1 + z_{eq}}{1 + z} \right)^{2}, \end{eqnarray}$$

(9)

which approximates well at the redshifts z ≳ 10. It can be seen that the presence of PBHs leads to the production of more massive dark haloes at high redshifts as compared with the case of purely adiabatic inflationary perturbations.

In the spherical top hat collapse model, after the formation of a structure with mass M, the average density of a halo will be |$\rho _H = \Delta _c \, \rho _{M}(z)$|⁠, where Δ_c = 18π²

$$\begin{eqnarray} \rho _H = \frac{3 M}{4 \pi R_{vir}^3} = 18 \pi ^2 \rho _{crit} \Omega _M (1+z)^{3}, \end{eqnarray}$$

(10)

where R_vir is the virial radius. At high redshifts z ≳ 10 sufficiently dense haloes are formed with a small velocity dispersion σ ∼ 1 km s⁻¹, which leads to active interaction of PBHs with each other. The dynamical evolution of such haloes is similar to that of globular star clusters, in particular, the core collapse occurs for them (Lynden-Bell & Wood 1968; Cohn 1980) that is shrinking of the central region with increasing in density, it further increases the rate of PBH interactions. However, the structures are formed hierarchically from the bottom up, which can lead to the destruction of small haloes in the process of their absorption by large structures. It is a difficult task to quantitatively take into account the dynamic internal evolution of haloes and their interaction with each other during the structure formations. In this paper, we accept the common assumption that the dark haloes are not destroyed with time and focus only on the internal dynamics of the halo and its influence on the PBH merger rate.

3 DYNAMICS OF PBH IN THE EARLY STRUCTURES

Two processes affect the dynamics of PBHs in a halo: the interaction of PBHs with each other and with DM particles. Both these processes can be characterized by diffusion coefficients using the Fokker–Planck equation (Binney & Tremaine 2008). They show the average rate of change in the body’s velocity (PBH in our case) as a result of many weak gravitational encounters with other objects. The interaction between PBHs is characterized by the following coefficient

$$\begin{eqnarray} \langle (\Delta v)^2 \rangle = \frac{4 \sqrt{2} \pi G^2 \rho _{BH} m \ln {\Lambda }}{\sigma }, \end{eqnarray}$$

(11)

where σ is the one-dimensional velocity dispersion, ρ_BH = ρ_Hf_PBH is the density of PBHs in the halo center and ln Λ ∼ 10 is the Coulomb logarithm. The characteristic time-scale of PBH interactions is called the relaxation time t_r, which is defined as t_r = σ²/〈(Δv)²〉 (see details in Binney & Tremaine 2008)

$$\begin{eqnarray} t_r = \frac{0.34 \, \sigma ^3}{G^2 m \rho _{BH} \ln {\Lambda }}. \end{eqnarray}$$

(12)

If the halo consisted only of PBHs, then its evolution is not different from the dynamics of the simplest globular star cluster and proceeds according to the scenario of a gravitational catastrophe (Lynden-Bell & Wood 1968): under the influence of pairwise interactions of PBHs, the central region of the halo (hereinafter core) is compressed to a sufficiently large density. The core collapse time is a certain number of relaxation times t_cc = βt_r, where the proportionality constant β is the number of those times, which depends on the choice of the initial density profile (Quinlan 1996). Evolution after collapse is driven by three-body interactions with the formation of binary systems, which stop the collapse of the core and lead to the self-similar expansion of the cluster (Lee 1987; Takahashi 1996; Binney & Tremaine 2008). The merger rate of such binaries was estimated in Franciolini et al. (2022), but this channel is less efficient than mergers of binaries formed in the early Universe.

On the other hand, if a halo contains DM particles in addition to PBHs, then the interaction of PBHs with DM particles can be characterized by the dynamical friction coefficient (Merritt 2013)

$$\begin{eqnarray} \langle \Delta v \rangle = - \frac{4 \pi G^2 \rho _{DM} m \ln {\Lambda }}{\sigma ^2}, \end{eqnarray}$$

(13)

where ρ_DM = (1 − f_PBH)ρ_H is the density of DM particles. Dynamical friction leads to the energy loss of PBHs, as a result of which they settle in the centre of the halo. Similarly to the case considered above, the characteristic dynamical friction time can be defined as follows:

$$\begin{eqnarray} t_{df} = \frac{3}{8} \sqrt{\frac{2}{\pi }} \frac{\sigma ^3}{G^2 m \rho _{DM} \ln {\Lambda }}. \end{eqnarray}$$

(14)

It is physically obvious that, for a small fraction of PBHs in the DM, their dynamics is determined only by the interaction with DM particles, and the scattering of PBHs on each other is negligible. In the intermediate case, the picture is as follows: in the early stages of evolution, dynamical friction is the dominant process, as a result of which PBHs will sink into the centre of the halo. When a sufficiently high concentration of PBHs is reached, their further dynamics will already be determined by pair interactions with each other. From the condition that the diffusion coefficients are equal σ〈Δv〉 ∼ 〈(Δv)²〉, one can find the density of PBHs, starting from which the dynamics will be determined only by the pairwise interaction of PBHs. It can be seen that this occurs at ρ_BH ∼ ρ_DM, which was obvious from physical considerations. PBH settling occurs on the scale of dynamical friction time equation (14). Further, as it was outlined, the evolution will occur due to interactions of PBHs with each other on the relaxation time-scale equation (12), where ρ_BH ≈ ρ_DM. Therefore it is convenient to choose the following quantity

$$\begin{eqnarray} t_{ch} = \frac{\sigma ^3}{G^2 m \rho _{c} \ln {\Lambda }} \end{eqnarray}$$

(15)

as a characteristic time for the evolution of a halo consisting of both PBHs and DM. Here, ρ_c is the central density of halo and σ is calculated using the Jeans formula (Binney & Tremaine 2008)

$$\begin{eqnarray} \sigma ^2 = \frac{1}{\rho _c} \int _0^{\infty } \rho _H(r^{\prime }) \frac{G M(r^{\prime })}{r^{\prime 2}} \, dr^{\prime }, \end{eqnarray}$$

(16)

where ρ_H = ρ_BH + ρ_DM is the DM halo density profile.

For further analysis, it is necessary to determine the halo core collapse time for different fractions of PBHs f_PBH in the DM composition. For this, the Fokker–Planck kinetic equation is used, the procedure for numerical solution of which is described in Cohn (1980), Vasiliev (2017), and Stasenko, Kirillov & Belotsky (2022). In this work for definiteness, the Burkert profile is used as the initial density distribution (Burkert 1995)

$$\begin{eqnarray} \rho _H = \frac{\rho _c}{\left(1 + r/r_0 \right) \left(1 + r^2/r_0^2 \right)}, \end{eqnarray}$$

(17)

where r₀ is the radius that determines the size of the halo core. We also assume that both PBHs and DM particles are equally distributed in space. That is, it is assumed that the fraction of PBHs matches f_PBH. However, during the process of halo formation in the central region, the PBH number density can be higher, but this requires careful numerical simulation. In this case, the gravitational waves constraints will be relaxed slightly stronger.

As an example, Fig. 2 shows the evolution of PBHs distribution in the halo formed at z_f = 20 (it corresponds t_f ≈ 180 Myr). The PBH fraction was chosen to be f_PBH = 0.1, and the parameter r₀ was chosen so that R_vir/r₀ = 5 and the halo mass |$M = 3 \cdot 10^4 \, {\rm M}_{\odot }$|⁠, which corresponds to a typical halo forming during this epoch, see Fig. 1. It can be seen that during the time t_cc ∼ 1.2 Gyr the core of the halo shrinks to r_c ∼ 1 pc while the density increases to |$\rho \sim 100 \, M_{ \odot }$| pc⁻³. We stopped the calculations when the number of PBHs in the central region is N ≈ 20. Subsequent evolution within the framework of the Fokker–Planck equation would lead to further compression of the core to an infinite density, that is unphysical. In a realistic scenario, as noted above, the core collapse is terminated due to the formation of binaries, which act as a heat source and lead to the expansion of the cluster. However, post-collapse evolution may be more complex due to the dominant role of DM particles composing the dark halo. Never the less, it can be noted that the natural formation of PBH clusters is possible, namely, small dark structures with a significant concentration of PBHs in the central regions.

$PBH density profile evolution (solid lines) in the halo with mass $M = 3 \cdot 10^4\, {\rm M}_{\odot }$ formed at zf = 20 and the PBH fraction fPBH = 0.1. The dashed line is the density profile of DM particles.$

Figure 2.

PBH density profile evolution (solid lines) in the halo with mass |$M = 3 \cdot 10^4\, {\rm M}_{\odot }$| formed at z_f = 20 and the PBH fraction f_PBH = 0.1. The dashed line is the density profile of DM particles.

Open in new tab Download slide

Fig. 3 shows the halo core collapse time in units of characteristic times equation (15) formed at z_f = 10 for different PBH fractions f_PBH. The red dots in the graph are obtained as the result of the numerical solution of the Fokker–Planck equation and are well-approximated by the following formula represented by a solid line on the graph

Figure 3.

Core collapse time of the characteristic halo mass formed at z_f = 10. The red dots are given by the direct numerical solution of the Fokker–Planck equation. The solid curve is given by equation (18).

Open in new tab Download slide

$$\begin{eqnarray} \frac{t_{cc}}{t_{ch}} = 15.9 \, \Big (1.3 \, e^{2.1 f_{PBH}}-1\Big). \end{eqnarray}$$

(18)

Due to computational difficulties, we considered only haloes containing N_PBH ≥ 30, so the Fokker–Planck equation was not solved in the region f_PBH ≲ 10⁻². In fact, for such fractions f_PBH, the number of PBHs in dark haloes will be less than 30 in this epoch, see equation (9). It was found that the core collapse time, expressed in characteristic times equation (15), is universal and does not depend on the profile parameter r₀ and the moment of halo formation, which will be also reflected in Section 5. However, it may depend on the specific choice of the initial density profile (Quinlan 1996).

The analysis of Vaskonen & Veermäe (2020) suggests that the core collapse time t_cc does not depend on the PBH fraction f_PBH and corresponds to a halo consisting only of PBHs, which in our case is |$t_{cc} \approx 150 \, t_{ch}$|⁠. However, as can be seen from Fig. 3 for the PBH fraction f_PBH ≲ 0.1, such estimation is strongly overestimated. Haloes dominated by DM evolve faster: as f_PBH decreases, the core collapse time becomes constant |$t_{cc} \approx 5 \, t_{ch}$|⁠, which corresponds to the fact that the PBH dynamics in the halo is mainly determined by dynamical friction.

Let us estimate the characteristic time equation (15) for haloes formed at different z. To do this, we set |$\sigma \sim \sqrt{G M/R_{vir}}$| and assume that the halo central density is α times bigger than the average matter density in the Universe ρ_c = αρ_M(z)

$$\begin{eqnarray} t_{ch} \sim 0.7 \, \frac{10^4}{\alpha } \left(\frac{N}{10^3} \right) \left(\frac{1+z}{10}\right)^{-3/2} \, \text{Gyr}, \end{eqnarray}$$

(19)

where we set ln Λ = 10 and N – the number of PBHs in the halo. The choice of α ≈ 10⁴ corresponds to R_vir/r₀ = 5 in the density profile of equation (17). It can be seen that at the PBH fraction in DM f_PBH ≈ 0.1, haloes formed at z ≈ 10 experience core collapse within the time t_cc ≲ t_H.

4 PBH BINARIES

In the early Universe, two PBHs at some moment decouple from the Hubble expansion and form a bound binary system, which will subsequently experience the inspiral motion due to the emission of gravitational waves and eventually merge. The lifetime of such a binary is given by Peters (1964)

$$\begin{eqnarray} t_{gw}=\frac{3 \, c^5 a^4 j^7}{170 \, G^3m^3}, \end{eqnarray}$$

(20)

where m – mass of each PBH, a and |$j = \sqrt{1-e^2}$| are the semimajor axis and dimensionless angular momentum of the binary system and e is the eccentricity. The distribution of parameters a and e (which is easily converted to j) of binaries can be found in Nakamura et al. (1997), Sasaki et al. (2016), and Eroshenko (2018).

The merger rate of the unperturbed PBH binaries is given by Raidal et al. (2019) and Vaskonen & Veermäe (2020)

$$\begin{eqnarray} \mathcal {R}_0 = \frac{3.1 \times 10^6}{\text{Gpc}^{3} \, \text{yr}} S \, f_{PBH}^{53/37} \left(\frac{m}{{\rm M}_{\odot }} \right)^{-32/37}, \end{eqnarray}$$

(21)

where S is the suppression factor arising due to PBH interactions during the formation of a binary

$$\begin{eqnarray} S \approx 0.24 \left(1 + \frac{2.3 \, \sigma _M^2}{f^2_{PBH}} \right)^{-21/74}, \end{eqnarray}$$

(22)

where σ_M ≈ 0.005 is the variance of matter fluctuations at the moment of PBH pair formation.

As shown in the previous section, early DM haloes evolve on time-scales smaller than the age of the Universe. Therefore, it can be expected that the parameters of the binaries will be perturbed as a result of scattering with other PBHs. Moreover, these PBH binaries formed in the early Universe are highly eccentric j ≪ 1. The scattering of such binaries with a single PBH will lead to a decrease in eccentricity (Jedamzik 2020), which ultimately leads to a significant increase of the lifetime. Taking this into account, equation of the merger rate (21) will be modified

$$\begin{eqnarray} \mathcal {R} = \mathcal {R}_0 S_f, \end{eqnarray}$$

(23)

where S_f ≲ 1 is the suppression factor that shows the fraction of binaries remaining unperturbed and is calculated in the next section. Scatterings of a binary system with other PBHs can be divided into the following: encounters with large impact parameters, which are of a tidal nature; and close hard scatterings.

Let us estimate the first. The change in energy in a binary system as a result of such tidal interaction will be (Binney & Tremaine 2008)

$$\begin{eqnarray} \Delta E_{tid} \approx \frac{8 \, G^2 m^3 a^2}{3 \, v^2 b^4}, \end{eqnarray}$$

(24)

where b is the impact parameter and v is the relative velocity (further, for estimates, it is assumed that v = σ). The rate of energy change is obtained by the standard procedure of multiplying equation (24) by a factor |$n_{BH} \sigma \, 2 \pi b \, db$| and integration over the impact parameter from b_min = a to ∞

$$\begin{eqnarray} \frac{d E_{tid}}{dt} = \frac{8 \, \pi G^2 \rho _{BH} m^2}{3 \, \sigma }. \end{eqnarray}$$

(25)

The lower limit of integration was chosen as a based on the consideration that for b < a, it is assumed that the scattering is already strong. It can be shown that an order of magnitude estimate of the change of j is (Ali-Haïmoud et al. 2017)

$$\begin{eqnarray} \Delta j \sim \frac{\Delta E_{tid}}{E}, \end{eqnarray}$$

(26)

then, the rate of change in j will be

$$\begin{eqnarray} \frac{d j}{d t} \sim \frac{8 \, \pi G \rho _{BH} a}{3 \, \sigma }. \end{eqnarray}$$

(27)

Let us estimate the change time of the dimensionless angular momentum by a value of the order of j itself (then the lifetime of binary will increase by the factor ∼2⁷, see equation 20) as follows:

$$\begin{eqnarray} t_{j} \sim \frac{j}{ dj/dt} \sim \frac{j \sigma }{ a G \rho _{BH} }. \end{eqnarray}$$

(28)

At the beginning of the halo evolution, this time is very small (due to the smallness of ρ_BH), but then the PBH density begins to dominate and the quantity |$1/\sqrt{G \rho _{BH}}$| becomes about the dynamical time in the centre of the halo t_dyn ∼ r_c/σ, which for ρ_BH = 100 M_⊙ pc⁻³ (see Fig. 2) is t_dyn ∼ 1 Myr. Thus, on the time scales of halo evolution, the characteristic time of the change of j can be estimated as follows:

$$\begin{eqnarray} t_j \sim j t_{dyn} \frac{r_c}{a}, \end{eqnarray}$$

(29)

that is several dynamical times, since j ∼ 10⁻² and a ∼ 100 a.u. (Ali-Haïmoud et al. 2017) and r_c ∼ 1 pc. Thus, the efficiency of tidal perturbations of binaries regarding changing j becomes important on the time-scale of halo evolution.

Another process that may be responsible for the perturbation of the binary parameters is strong scattering with a single PBH. We will assume that the pericentre (the distance of closest approach) for such scattering is r_p ∼ a. Then, the cross-section of such a process will be

$$\begin{eqnarray} \Sigma \sim \pi a^2 \left(1 + \frac{2 \, Gm }{a \sigma ^2} \right). \end{eqnarray}$$

(30)

Since the binding energy of the binaries is greater than the characteristic kinetic energy of the PBHs in the halo, the gravitational focusing approximation is valid (i.e. the second term dominates in equation 30), then the rate of close scattering can be estimated as follows:

$$\begin{eqnarray} \Gamma &\sim & n_{BH} \Sigma \sigma \sim \frac{2 \pi a \rho _{BH} G}{ \sigma } \\ &\sim & 0.1 \left(\frac{a}{10 \text{au}} \right) \left(\frac{\rho _{BH}}{100 {\rm M}_{\odot } \text{pc}^{-3}} \right) \left(\frac{1 \text{km}\, \text{s}^{-1}}{\sigma } \right) \, \text{Gyr}^{-1}. \end{eqnarray}$$

(31)

So, it can be seen that strong scatterings are less efficient than long-range tidal interactions, but, never the less, they can also make some contribution to the perturbation of binary parameters. Thus, on the time-scales of halo core collapse t_cc, the parameters of the binaries can be significantly perturbed, leading to an increase in the lifetime of the PBH pairs and a reduction in the merger rate.

An important remark should be made. In addition to the merger of binaries formed in the early Universe, direct mergers of PBHs in the dark halo are also possible (Bird et al. 2016; Clesse & García-Bellido 2017). Since halo evolution leads to an increase in the PBHs number density in the central region, direct mergers of PBHs can make a significant contribution to gravitational wave events. However, accurate accounting for these mergers in the modern era requires knowledge of the abundance of such clusters. As an example, it can be seen from Fig. 2 that the maximum PBH number density in the halo is reached by redshift z ∼ 4. However, in the further process of structures formation, such clusters can be destroyed, which will reduce the rate of direct PBH mergers. We leave this analysis for the next work.

5 SUPPRESSION FACTOR AND THE PBH MERGER RATE

To estimate the suppression factor S_f in equation (23), we use the formalism of Vaskonen & Veermäe (2020), the idea of which is to calculate the fraction of PBH binaries that are in haloes experiencing core collapse. As shown in the previous section, it is assumed that the lifetime of binaries becomes longer than the Hubble t_H and, therefore, they will not contribute to the merger rate. If at some redshift z_f a halo is formed containing N_c PBHs and experiencing core collapse to the redshift z_c, then a halo containing a smaller number of PBHs N < N_c (but forming at the same redshift z_f) will collapse to redshifts z > z_c. The fraction of PBH binaries that will not be perturbed as a result of the evolution of haloes formed at z_f is estimated as follows:

$$\begin{eqnarray} S_f(z) &= &1-\sum _{N = 3}^{N_c} p_N(z_f) \\ && - \sum _{N^{\prime } \gt N_c} \left(\sum _{N = 3}^{N_c} p^{\prime }_N(z_f) \right) p_N^{\prime }(z_f). \end{eqnarray}$$

(32)

Here, p_N is the distribution function of the number of PBHs N in the halo

$$\begin{eqnarray} p_N(z) \propto \frac{1}{\sqrt{N}} e^{-N/N^{*}(z)}, \end{eqnarray}$$

(33)

which follows from the Press–Schechter mass function and it is valid at high redshifts (Hütsi et al. 2019) and N*(z) = M*(z)/m is the characteristic number of PBHs in the halo. The physical meaning of the terms in equation (32) is as follows: the second term is the probability that the binary is in a halo containing N PBHs, the third term is the probability that the binary is in substructures of the halo with the number of PBHs N′ > N_c (but inside this large halo there will be substructures with N < N_c). As noted earlier, binaries which are in haloes collapsing to the redshift z will not contribute to the merger rate. Therefore, the fraction of these binaries should be subtracted, which is reflected in equation (32).

The further idea is to find the halo formation redshift z_f at which the suppression factor will be minimal to the redshift z. This is done as follows: some redshift is taken at which haloes are formed. Then, the critical value of PBHs N_c in the halo, which will experience the core collapse to the redshift z, is found. Next, we calculate the suppression factor using equation (32), then we shift along axis z_f in the direction of decreasing S_f until its minimum value is reached. To implement this algorithm, it is necessary to calculate the halo core collapse time. This can be done either directly by solving the Fokker–Planck equation or using equation (18). Since we are interested in the modern merger rate, it is further assumed that z = 0.

Fig. 4 shows the modern suppression factor calculated for the Burkert density profile (17) for different magnitudes of the parameter R_vir/r₀. The thick dots correspond to the numerical solution of the Fokker–Planck equation, while the solid curves were obtained using the expression for the core collapse time in equation (18). As expected, with the increase of the halo density (a decrease of the parameter r₀), the suppression factor S_f decreases, which is due to the fact that the evolution of the halo proceeds more rapidly, see equation (19). It should be repeated that we do not take into account the influence of hierarchical structure formation, which can lead to tidal disruption of the halo. However, in such process, the central part of the halo will probably continue to evolve, and, in addition, acceleration of the core collapse is possible (Quinlan 1996; Nishikawa, Boddy & Kaplinghat 2020).

The suppression factor at z = 0 calculated for different values Rvir/r0: curves from bottom to top correspond to a decrease of the r0. The solid lines and dots are obtained, respectively, by estimating the lifetime from equation (18) and by calculating directly the core collapse time using the Fokker–Planck equation.

Figure 4.

The suppression factor at z = 0 calculated for different values R_vir/r₀: curves from bottom to top correspond to a decrease of the r₀. The solid lines and dots are obtained, respectively, by estimating the lifetime from equation (18) and by calculating directly the core collapse time using the Fokker–Planck equation.

Open in new tab Download slide

It is also important to note that the suppression factor S_f does not depend on the PBH mass m. This is due to the fact that the characteristic time t_ch depends only on the number of PBHs, see equation (19). On the other hand, the characteristic number of PBHs in the halo N* also does not depend on the mass of PBH, as can be seen from equation (9), and also from the fact that Poisson fluctuations are determined only by the number of PBHs |$\delta \sim 1/\sqrt{N}$|⁠.

Fig. 5 shows the PBH merger rate with and without the suppression factor (given by equation 21), where the parameter r₀ was chosen as R_vir/r₀ = 5 for density profile in equation (17). The shaded area corresponds to the measurements of the LIGO–Virgo–KAGRA collaborations. It can be seen that the fraction of PBHs in DM allowed from the point of view of observations of the BH merger rate is relaxed weakened to f_PBH ∼ 0.1, which corresponds to the constraints on PBHs at masses |$\sim 10 \, {\rm M}_{\odot }$| in the modern Universe on dwarf galaxies (Brandt 2016; Koushiappas & Loeb 2017) and lensing (Alcock et al. 2000; Oguri et al. 2018).

$The PBH merger rate with the suppression factor (solid line) and without it (dashed line). The shaded area is the PBH merger rate inferred by the LIGO–Virgo–KAGRA collaboration $\mathcal {R} = 17.9 \div 44$ Gpc−3 yr−1 (Abbott et al. 2023).$

Figure 5.

The PBH merger rate with the suppression factor (solid line) and without it (dashed line). The shaded area is the PBH merger rate inferred by the LIGO–Virgo–KAGRA collaboration |$\mathcal {R} = 17.9 \div 44$| Gpc⁻³ yr⁻¹ (Abbott et al. 2023).

Open in new tab Download slide

In the case of PBH clustering stronger than the Poisson noise predicted in works Rubin, Sakharov & Khlopov (2001) and Khlopov, Rubin & Sakharov (2005) (see also the review of Belotsky et al. 2019), the constraints are likely to be weakened more. Since in this case, the PBH structures will be formed in a younger Universe. These clusters will experience core collapse on smaller time-scales (Stasenko et al. 2022), as a result of which most of the binaries will be perturbed compared to pure Poisson clustering. However, in this case, PBH binaries will already be formed in these clusters through dynamical channels: due to the emission of gravitational waves during close approaches (Bird et al. 2016; Clesse & García-Bellido 2017; Stasenko & Kirillov 2021; García-Bellido, Jaraba & Kuroyanagi 2022) and as a result of three-body interactions (Franciolini et al. 2022). The merger rate of these binaries may dominate over early ones, which requires a separate analysis.

6 CONCLUSION

PBHs with the masses |$\sim 10 \, {\rm M}_{\odot }$| are the subject of active discussion regarding strongest constraints on their contribution to the composition of DM (f_PBH ≲ 10⁻³) due to the observation of gravitational wave signals. This constraint assumes that binaries forming in the early Universe are not perturbed and/or destroyed with time. However, the Poisson initial space distribution of PBHs leads to the early formation of a dark haloes. In such haloes, the probability of perturbing and/or destroying a pair of PBHs is significant. In this work, the dynamics of PBHs in early DM structures and its influence on the merger rate were considered.

To study the dynamics of PBHs in early dark haloes, the Fokker–Planck kinetic equation was solved numerically. We have shown that the halo core collapse time essentially depends on the fraction of PBHs in the DM composition f_PBH. For f_PBH < 1, the halo evolves much faster than in the case when the entire DM consists of PBHs, which is due to dynamical friction against DM particles. Halo evolution leads to an increase in the density of PBHs in the central region of the halo, which leads to the perturbation of PBH binaries, as a result of which their lifetime can increase significantly.

The result was obtained under the assumption that early haloes are not destroyed during the structure formation and their internal dynamics leads to core collapse due to the interaction of PBHs both with each other and with DM particles. By calculating the core collapse time, we estimated the suppression factor for the modern merger rate. Namely, it is assumed that if a binary finds itself in a halo that experiences core collapse to the redshift z = 0, then it will no longer contribute to the merger rate. This is due to the fact that in the processes of interactions of a binary with other PBHs, its parameters will be significantly perturbed, leading to a binary lifetime will exceed the age of the Universe.

Ultimately, we showed that the constraints on the fraction of PBHs in DM can be relaxed to f_PBH ≲ 0.1, which is compatible with the constraints obtained in the modern Universe from dwarf galaxies and lensing. Moreover, in the case of initial PBH clustering, the constraints will probably be weakened more strongly due to the fact that an even larger fraction of binaries will be perturbed.

ACKNOWLEDGEMENTS

The authors are grateful for useful discussions by Yu.N. Eroshenko and S.G. Rubin. The work was supported by Russian Science Foundation (RSF) grant number 23-42-00066, https://rscf.ru/project/23-42-00066/.

DATA AVAILABILITY

No new data were generated or analysed in support of this research.

Footnotes

Throughout this article, we use the standard cosmological ΛCDM model with parameters Ω_CDM = 0.27, Ω_B = 0.05, Ω_Λ = 0.68.

References

Abbott

et al. ,

2023

Phys. Rev. X

011048

Month:	Total Views:
October 2023	64
November 2023	55
December 2023	23
January 2024	44
February 2024	19
March 2024	26
April 2024	27
May 2024	16
June 2024	14
July 2024	14
August 2024	7
September 2024	12
October 2024	14
November 2024	15
December 2024	12
January 2025	8
February 2025	4
March 2025	7
April 2025	6
May 2025	3

Article Contents

Influence of early dark matter haloes on the primordial black holes merger rate

ABSTRACT

1 INTRODUCTION

2 EARLY DM HALOES

3 DYNAMICS OF PBH IN THE EARLY STRUCTURES

4 PBH BINARIES

5 SUPPRESSION FACTOR AND THE PBH MERGER RATE

6 CONCLUSION

ACKNOWLEDGEMENTS

DATA AVAILABILITY

Footnotes

References

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only