Impact of dark matter spikes on the merger rates of Primordial Black Holes

Mergers of Primordial Black Holes (PBHs) may contribute to the gravitational wave mergers detected by the LIGO-Virgo-KAGRA (LVK) Collaboration. We study the dynamics of PBH binaries dressed with dark matter (DM) spikes, for PBHs with extended mass functions. We analyze the impact of DM spikes on the orbital parameters of the PBH binaries formed in the early Universe and calculate their merger rates at the age of the Universe today. We consider two possible scenarios for the dynamics of the dressed binaries: assuming that either the DM spikes are completely evaporated from the binaries before merger or they remain static until the merger. Contrary to previous studies, we find that the presence of spikes may increase or decrease the present-day PBH merger rates, in some cases dramatically. Comparing with merger rates reported by the LVK Collaboration in the third Gravitational Wave Transient Catalog (GWTC-3), we derive approximate constraints on the fraction of Solar-mass PBHs in cold dark matter as f pbh ≤ 𝒪(10-5–10-3), depending on the mass function. Our calculations are valid only for the idealized scenarios in which the DM spikes are either evaporated or static. However, they suggest that the impact of DM spikes on PBH merger rates may be more complicated than previously thought and motivate the development of a more general description of the merger dynamics, including feedback of the DM spikes in highly eccentric PBH binaries.


Introduction
Primordial Black Holes (PBHs) may form shortly after the Big Bang due to the collapse of very large density fluctuations on horizon scales [1,2].Other possible channels of PBH formation include the collapse of domain walls, the collapse of cosmic string loops, and bubble collisions [3,4].The mass spectrum of PBHs is dependent on the scale of perturbations leading to their origin [2].PBHs lighter than ∼ 10 −19 M ⊙ would have already evaporated due to Hawking radiation but PBHs heavier than 10 −19 M ⊙ could in principle be very interesting as they can be a possible non-particle candidate for cold dark matter [4][5][6][7].PBHs of mass ≳ 1 M ⊙ may contribute to the binary black hole mergers detected by LIGO-Virgo-KAGRA (LVK) Collaboration [8][9][10].Solar-mass PBHs might also explain the origin of supermassive black holes, which without lighter seed black holes, may struggle to grow rapidly enough to reach their present observed masses [11][12][13][14].The existence of bright galaxies observed by the James Web Space Telescope (JWST) at high redshift of z ≳ 10 are also difficult to explain as per the standard ΛCDM model [15][16][17][18].But if PBHs with mass around ≳ 100 M ⊙ are formed in the early Universe then they can seed the formation of massive structures at such high redshift [19].
As per the Gravitational Wave Transient Catalog-3 (GWTC-3), the LVK Collaboration has reported the detection of ∼ 90 compact binary black hole (BBH) mergers in the component mass range of 5 − 105 M ⊙ [20,21].It has been proposed that some of the mergers detected by LVK Collaboration could originate from the collisions of black holes that might be primordial in nature [8][9][10].It is predicted that PBHs can decouple from the Hubble expansion very early and form bound systems such as binaries which can ultimately merge to produce ripples in the fabric of space-time.The detection of gravitational waves (GWs) from such mergers can lead to constrain the abundance of PBHs which in turn can impose constraints on various models of inflation leading to their origin (e.g.[22][23][24][25]).
In the mass range probed by the LVK Collaboration, we know that the fractional abundance f pbh of PBHs in cold dark matter (CDM) should not be greater than ∼ 10 −3 [10,[26][27][28][29][30][31][32][33][34][35][36][37][38][39].This implies the existence of other candidates contributing to the dark matter, which may have a particle nature.The presence of such dark particles can lead to the formation of dense spikes around all PBHs [40,41].PBHs with DM spikes around them are often referred to as 'dressed' PBHs.In principle, spikes around the PBHs can be formed by any particle-like DM.However, gamma ray observations exclude this possibility if DM particles have a large annihilation cross section, as in the case of classical WIMPs [42][43][44][45][46][47][48].We therefore consider the scenario in which the dominant particle DM is not WIMP-like or is WIMP-like with a very small annihilation cross section [45,47].The merger rate of PBH binaries without DM spikes has been extensively studied, for monochromatic as well as extended PBH mass functions (MF) (see e.g.[26,[49][50][51][52][53][54][55]).The effects of DM spikes on the orbital parameters of PBH binaries with monochromatic PBH mass functions was taken in to account by Ref. [56] and it was shown that the merger rates of PBH binaries dressed with DM spikes may be around twice that of binaries without spikes.The formation of a common DM spike around the whole of the PBH binary was neglected in Ref. [56].Taking this into account, Ref. [57] showed that the merger rate of PBH binaries with DM spikes may be as much as 6 − 8 times larger than PBH binaries without DM spikes.
Since the BBH mergers observed by LVK show a range of component masses, so it becomes important to extend these studies to dressed PBHs with extended mass functions.This is crucial in order to quantify how frequently the mergers of PBHs with realistic mass functions could take place in the Universe.In addition, dressed PBHs with unequal masses can lead to intermediate mass ratio inspirals (IMRIs) or extreme mass ratio inspirals (EMRIs) in which the lighter PBH in the binary is hundreds to millions of times less massive than its companion.In such systems, it may be possible to see the effects of DM spikes directly on the GW signal [58][59][60][61][62][63][64].The possibility to distinguish the gravitational waves originating from the merger of PBH binaries with DM spikes and without DM spikes was recently explored in Ref. [65].In order to accurately study the GW signals from these PBH binaries, however, it is important to understand their formation mechanism and time evolution up to the merger.
In this paper, then, we consider PBH binaries with DM spikes having a wide range of PBH masses and estimate their merger rates by extending the approach of Ref. [56].The simultaneous evolution of PBHs and DM spikes in the binaries can be a complex problem in full generality [61,63,65].We therefore calculate the merger rates of these binaries assuming that either the DM spikes are completely evaporated before the merger (as in Ref. [56]) or that they remain static until the merger.These assumptions provide a starting point for understanding how DM spikes affect the orbital properties and evolution of unequal-mass PBH binaries.We compare our results with the merger rates reported by LVK in GWTC-3, in order to derive approximate constraints on the abundance of PBHs in dark matter f pbh .
The paper is organized in the following manner.In Sec. 2, we review the accretion process of DM spikes around isolated PBHs by closely following the formalism of Ref. [46].In this section, which can be skipped by the expert reader, we describe the density profile of dark matter spikes expected around PBHs.Then, in Sec. 3, we describe the formation mechanism of PBH binaries with and without DM spikes taking into account the impact of DM spikes.This set up is an improved and generalized description of binary formation given in Refs.[52,56], for extended mass functions with the addition of DM spikes into the picture.In Sec. 4, we analyze the variation of the final merger time of PBH binaries with and without DM spikes as a function of the initial semi-major axis and eccentricity, which also includes an updated review Ref. [56].After that in Sec. 5, we calculate the present-day merger rates of PBH binaries surrounded by DM spikes by convolving the final merger time with the probability distributions of initial semi-major axis and eccentricity of the systems.We present these merger rates for three different extended PBH mass functions as illustrative examples.Then, in Sec.5.4, we discuss how the gap between the merger rates of PBH binaries with and without DM spikes varies with the abundance of PBHs in CDM f pbh , a novel aspect in comparison to previous studies.We also highlight the mass ratio distribution of PBHs merging today, as described in Appendix C. Finally in Sec.6, we conclude and point out the possible future directions.

Formation of DM spikes around PBHs
In this section, we review the formation mechanism and density profile of DM spikes around isolated PBHs proposed originally in Ref. [40,41], following closely the treatment in Ref. [46].We start with the assumption that in the early Universe, non-annihilating CDM particles (such as axions or WIMPs with almost negligible annihilation cross-sections) can encounter isolated PBHs in their vicinity, due to which their motion becomes influenced by the combined effect of the gravitational attraction of the PBHs and Hubble expansion.Then, the equation of motion of a DM shell of radius r around a PBH of mass m pbh in the FLRW metric is given by: where the first term corresponds to the gravitational attraction of the PBH and the second term denotes the decelerating Hubble expansion of the Universe.The evolution of these DM shells starts deep in the radiation era, corresponding to the initial conditions r = r i and ṙ = At some point, the gravitational attraction term in Eq. (2.1) starts to dominate over the Hubble expansion, i.e. ṙ becomes zero, due to which the DM shell starts to move towards the PBH and becomes gravitationally bound to it.This is known as the turnaround of the DM shell and occurs at a time t ta .The size of the DM shell at turnaround is known as the turnaround radius r ta of the shell.Shells within the turnaround radius are considered bound to the PBH.The turnaround of subsequent DM shells leads to the formation of a DM spike around the PBH which is shown pictorially in the left panel of Fig. 1.During radiation domination, the turnaround radius of the DM shell can be estimated analytically by substituting ṙ = 0 and H = 1/(2 t) in Eq. (2.1).We have verified that in radiation domination, the numerical estimate of the turnaround radius, holds better than its analytic estimate, in agreement with Ref. [46].So, we shall use the numerical estimate given by Eq. (2.2) as the turnaround radius of DM shells in the rest of the paper.As per Eq.(2.1), the variation of the distance r of different DM shells from the center of a PBH of mass m pbh = 100 M ⊙ are depicted in the right panel of Fig. 1.This figure shows that at an initial distance r i (defined at some arbitrary early time) smaller than 10 −5 pc , the gravitational attraction of the PBH strongly dominates over the Hubble expansion and the shells of the dark matter become quickly gravitationally bound to it.If r i lies in the range 10 −5 − 10 −4 pc, DM shells initially evolve under the simultaneous effect of gravitational attraction of the PBH and Hubble expansion.After the turnaround point (where gravitational attraction of the PBH becomes stronger than the Hubble expansion), they start moving towards the PBH and become bound to it.For r i greater than 10 −4 pc, the Hubble expansion dominates over the gravitational attraction of the PBH and the DM shells do not get gravitationally bound to the PBH before matter-radiation equality.

Density profile of DM spikes
Based on the PBH mass and time scale at which the kinetic decoupling of the DM particles occurs [66], the density profile of the DM spikes around the PBHs can be a power law with different radial indices, as discussed in detail in Ref. [67].Here, we assume that the formation of the PBHs takes place after the kinetic decoupling of the DM particles such that there is no previous accretion of DM particles around the PBHs [68].Since, in radiation domination, the density of matter varies as ρ ∝ a −3 and a ∝ t1/2 then the density of dark matter can be written as [46]: (2.3) In the above equation, the factor of one half comes from the fact that at matter-radiation equality, density of matter is half of the total energy density.We know that the size of the DM spike at any time is given by the turnaround radius r ta of the DM shells at that time.
Then, substituting Eq. (2.2) in Eq. (2.3), the density of the DM spike in radiation domination at a distance r from the center of a PBH of mass m pbh can be written as: with ρ eq being the total energy density of the Universe at matter-radiation equality (t eq ) and Ω cdm /Ω m ≈ 0.85 as the fraction of CDM in the matter density of the Universe [69].Here, the subscript 'sp' denotes the DM spike accreted around an isolated PBH.The density profile of ρ(r) ∝ r −9/4 given by Eq. (2.4) was originally introduced for the gravitational collapse of collisionless non-relativistic matter in an Einstein-de Sitter universe [70,71].We note that this profile is steeper than the density profile of ρ(r) ∝ r −3/2 for the accretion of DM particles around PBHs studied in Ref. [56].Moreover, Eq. (2.4) slightly differs from the density profiles of DM spikes calculated in Ref. [46] by a factor of Ω cdm /Ω m , which we have included here to account for the fact that DM does not comprise all of the matter content of the Universe.Assuming spherical symmetry, the mass of the DM spike around the PBH can be estimated as: r 3/4 for r < r ta r 3/4 ta for r ≥ r ta . ( The total mass of the DM spike at any moment is given by the mass enclosed within all the DM shells within the turnaround radius.Now, let us define s = a/a eq as the scale factor normalized to unity at matter-radiation equality.Then using H(s) = 4 π G ρ eq /3 h(s) with h(s) = √ s −3 + s −4 , the turnaround time of DM shells in terms of turnaround scale factor s ta is 1 : Substituting the values of the turnaround radius given by Eq. (2.2) and turnaround time given by Eq. (2.6) in Eq. (2.5), one can easily verify that the mass of the DM spike being accreted around the PBH up to matter-radiation equality is: This value is around half of the mass of the corresponding PBH and is in rough agreement with previous calculations of [40,56] which suggest that at matter-radiation equality, the mass of the DM spike becomes comparable to the mass of the PBH itself.Small differences compared to previous calculations arise from minor numerical factors, including the extra factor of Ω cdm /Ω m in the density profile of the DM spike given by Eq.(2.4).

Formation of PBH Binaries
In this section, we set up a model for the formation of PBH binaries in the presence of DM spikes around each PBH.In doing so, we first provide an improvement over Ref. [56] by self-consistently taking into account the growth of DM spikes during the process of binary formation.Then we extend the formalisms of Refs.[49,52] for the initial angular momentum of binaries by including the presence of DM spikes.The relevant technical details are given in the following sub-sections, but we here briefly summarize them.In Sec.3.1, we show that the impact of DM spikes on the initial semi-major axis of the PBH binaries depends on the scale factor at which the formation of the binaries takes place.As the scale factor of binary formation increases, the initial separation of dressed PBH binaries becomes increasingly smaller than the same PBH binaries without DM spikes.In Sec.3.2 we show that the change in initial angular momentum due the presence of DM spikes is negligible for PBH binaries decoupling early in radiation domination and its value is similar to Ref. [49].Similar to Ref. [52], we write f as the total abundance of PBHs in non-relativistic matter such that the fraction of PBHs in CDM is f pbh = Ω pbh /Ω cdm ≈ f /0.85.Depending on the formation mechanism, PBHs will possess a mass function P (m), which can be described by a normalized probability distribution function (PDF).Common parametrizations for the mass function include Power Law, Broken Power Law and Lognormal [4].In our calculations, we adopt a few specific extended PBH mass functions discussed later in Sec. 5. We assume that the PBHs of masses m i and m j have the discrete resolution of ∆ i and ∆ j respectively [52] 2 such that the equation of normalization for any PDF in the mass range from m min to m max can be written as: Since P (m i ) describes the probabilistic distribution of PBHs having mass m i then the fraction of matter in the forms of PBHs with mass m i becomes [52]: We assume that the PBHs are uniformly distributed in the Universe such that the average comoving separation3 between the PBHs of mass m i can be written as: Similarly, PBHs of mass m j have an average comoving separation xj in between them.Then, using the values of comoving separations xi and xj , we can calculate the average comoving separation between the PBHs of masses m i and m j as [52]: with and In the rest of the paper, the subscript 'ij' is omitted in the definition of µ ij and the comoving separation ⟨x ij ⟩ or xij .Note that, unlike in Ref. [52], we have allowed for different discrete spacing i.e. ∆ i ̸ = ∆ j for PBHs of mass m i and m j .In Sec.3.1 and Sec.3.2, we mention all relevant technical details as an extension of Ref. [52], which can be skipped by the knowledgeable reader.

Initial semi-major axis of PBH binaries
Now, we estimate how the presence of DM spikes with density profile given by Eq. (2.4) influences the formation of PBH binaries.Since all the neighbouring PBHs and matter-density fluctuations of the Universe also affect the dynamics of the binaries formation, so we include their impact too, which is mentioned in Section 3.2.The formation of a binary takes place between the two nearest neighbouring PBHs of masses m i and m j with x being the comoving separation between them.We assume that no other PBH lies inside a spherical volume of radius equal to the comoving separation x.The two PBHs decouple from the Hubble flow due to their mutual gravitational attraction and form a bound system.However, every PBH in the Universe may have multiple PBHs in its neighbourhood which also might decouple at a similar time, leading to hierarchical interactions of systems comprising more than two PBHs.These N-body interactions of PBHs have been studied in Ref. [26] where they conclude that if PBHs are a very tiny fraction of the total DM i.e. f pbh ≪ 1 then the possibility of hierarchical PBHs interactions becomes very low.In this work, we consider small values of f pbh ≲ 10 −2 and we therefore neglect the disruption of individual binaries via multiple PBH interactions in their vicinity.Now, similar to Eq. (2.1), in a binary, the equation of motion of the proper separation r b between the PBHs with masses m i and m j having DM spikes of masses m sp,i (t) and m sp,j (t) respectively is given as: where the dot represents the derivative with respect to proper time t.Here, we assume that each PBH and its DM spike can be treated as a point mass which depends on time.As the dressed PBHs evolve under the simultaneous effect of their mutual gravitational attraction and Hubble expansion, mass of each DM spike increases continuously with time t.Eventually at some point, the gravitational attraction term in Eq. (3.7) dominates over the expansion term and the pair of PBHs reach a maximum separation before turning around to approach each other again.This point is referred to as the point of binary decoupling and we consider this to be the time of formation of the binary.After binary decoupling, we no longer consider the growth of DM spikes around the isolated PBHs or around the binary as a whole 4 .So, we assign the initial parameters of the binary at the point of decoupling.Then, similar to [49] with the addition of DM spikes around PBHs in picture, we define a dimensionless separation, χ = r b /x and express Eq.(3.7) as: where primes represent the derivative with respect to scale factor s and h ≡ h(s) = √ s −3 + s −4 .Also, with a dimensionless parameter which quantitatively describes how the mass of the PBHs present in the volume occupied by the binary differs from the mass of the background matter lying in the same volume. 5Large values of λ correspond to PBH pairs having a large initial separation, due to which the PBHs decouple late as the background density dominates over the mutual density of the two PBHs for a longer time.Here, the only term signifying the effect of DM spikes in the PBH binary is: This explicit expression shows that the impact of DM spikes on the dynamics of binary formation is independent of the PBHs masses in the binary.This is applicable to any PBH binary dressed with DM spikes having the density profile of ρ(r) ∝ r −9/4 .Following Eq. (3.8), the formation of PBH binaries with and without DM spikes for λ = 1 is illustrated in the left panel of Fig. 2.This panel shows that PBHs with DM spikes decouple earlier than the PBHs without spikes due to the increased gravitational attraction with the presence of DM spikes6 .We estimate the initial semi-major axis a i of the binaries as half of the proper separation of the PBHs at the point of binary decoupling because this is the maximum separation shared by the PBHs of the binaries i.e. r b,max ≈ (1 + e i ) a i with e i ∼ 1.Hence, we take the initial dimensionless separation of the PBHs at binary decoupling i.e. when the curve χ/λ first turns over such that |χ| becomes equal to 2a i /x.Now, for different values of λ, the impact of DM spikes on the initial semi-major axis of PBH binaries can be seen in the the numerical solution of Eq. (3.8) shown in the right panel of Fig. 2. From this panel, it is clear that for λ ≪ 1 values, the behaviour of a PBH binary having DM spikes is not significantly different from the same binary having no DM spikes.We find that the scale factor s dec of binary decoupling is directly proportional to the value of λ.For λ ≪ 1, binaries decouple very early and hence do not have enough time to accrete substantial DM spikes.But for λ ≳ 1, the initial separation of the binary with DM spikes is comparatively smaller because the binary has enough time to accrete DM spikes.With the accretion of DM spikes, the gravitational attraction between the PBHs in the binary enhances and prevents their further evolution due to Hubble expansion.This ultimately makes the initial separation of the PBH binaries smaller than the PBH binaries without DM spikes.This behaviour of the initial semi-major axis of dressed PBH binaries qualitatively agrees  (3.8).The right panel shows the variation of the initial semi-major axis a i of the PBH binaries with and without DM spikes as a function of λ as per the numerical solution of Eq. (3.8).Here, λ is the dimensionless parameter given by (3.10) which is a proxy for the PBH separation and which quantifies how much larger the background density is compared to the mutual density of the two PBHs.In addition, x is the comoving separation of the PBHs in the binaries and s dec is the scale factor of the binary decoupling.
with the initial semi-major axis of PBH binaries having DM spikes with the density profile of ρ(r) ∝ r −3/2 shown in Fig. 3 of Ref. [56].From this panel, the initial semi-major axis of the PBH binaries without DM spikes for a given value of λ can be approximated as: with the initial semi-major axis for PBH binaries having DM spikes as: The subscript '0' signifies PBH binaries without DM spikes.

Initial angular momentum of PBH binaries
During the formation of binaries, apart from their mutual gravitational attraction, PBHs also experience tidal torque due the gravitational force of all other PBHs and large scale density perturbations in the Universe [72,73].This tidal force provides an initial angular momentum to the PBHs which avoids their head-on collision, and leads to the formation of very eccentric binaries.Now in Sec.3.2.1 and Sec.3.2.2,we calculate the probability distribution of this initial angular momentum of PBH binaries by including the presence of DM spikes around PBHs in the formalism of Ref. [49].

Tidal Torque due to all PBHs outside the binaries
The gravitational attraction of all the PBHs outside the binary leads to a net tidal force on the PBHs in the binary.The local tidal field produced by any PBH outside the binary is defined as T ij = −∂ i ∂ j ϕ with ϕ being the gravitational potential [49].Due to this tidal field, the PBHs in the binary experience a force given by [49,52]: Here, F is the tidal force per unit reduced mass of the binary.This perturbative force provides an angular momentum per unit reduced mass of the PBH binary as: Then, using T ∝ m total,d (s)/s 3 and the above equation, the tidal field experienced by the PBH binary at any scale factor s can be written as: Now substituting the value of the tidal torque given by Eq. (3.17) in Eq. (3.15), the angular momentum of the PBH binary becomes: with numerical solution Here, s i ≪ 1 is some initial scale factor corresponding to radiation domination and s dec is the scale factor at which the binary decouples from the Hubble expansion.For λ ≪ 1, the value of the angular momentum given by Eq. (3.19) for the PBH binary with DM spikes agrees with the value for PBH binaries without DM spikes calculated in Ref. [49].And from Eq. (3.18), we see that the presence of DM spikes around the distant PBHs impacts the tidal torque substantially only for s dec ≥ s eq corresponding to λ > 1.This is because the masses of DM spikes are lighter for PBH binaries which decouple much earlier than matter-radiation equality, meaning that the enhanced tidal gravitational force due to the presence of spikes is not enough to substantially change the initial angular momentum.Now, the dimensionless angular momentum of the PBH binary due to the tidal torques is [49,52]: with m b = [m i + m j + m sp,i (s) + m sp,j (s)] being the total mass of the binary at scale factor s.
For the tidal field generated by the distant PBH situated at comoving separation y ≫ x [49] and semi-major axis a given by Eq. (3.13), we can rewrite the dimensionless angular momentum of the binary as: The dimensionless constant characterizes the configuration of the PBH binary based on the PBH masses and their comoving separation x.Note that in the limit λ ≪ 1, the value of C(λ) becomes unity and matches with the case without DM spikes.The total reduced angular momentum of the PBH binary with DM spikes resulting from all distant PBHs in the Universe at comoving separation y ≫ x is: Now, applying the formalism described in Appendix 1 of Ref. [49], the probability distribution of the initial angular momentum of PBH binaries is given as: with γ λ = j/j λ and Here, j λ is the characteristic value of the initial angular momentum j for a given value of the dimensionless variable λ which specifies the comoving separation of the PBHs in the binaries.This distribution of initial angular momentum given by Eq. (3.24) for PBH binaries with DM spikes is valid for 0 < j < 1 and we set P (j) = 0 for j > 1.

Tidal torque due to large scale density perturbations
In the approximation that PBHs do not make up all of the dark matter, large scale matter density perturbations in the Universe also exert tidal forces on the PBH binaries.Since the dynamics can become very complex if the scale of density perturbations is either smaller or of the order of the separation of PBH binaries, so we are only considering the tidal torque generated by density fluctuations having comparatively larger length scales.Then, the tidal torque exerted on PBH binaries by such large scale matter density perturbations δ m is [49]: with ϕ p as the gravitational potential due to the large scale density perturbations.Since the tidal field generated by the perturbations δ m is Gaussian so the tidal tensor T eq ij is Gaussian, resulting in reduced angular momentum of the PBH binaries with DM spikes as Now as per the reduced angular momentum given by the above equation and Appendix 2 of Ref. [49], the variance of the distribution of j in the plane perpendicular to x is: with σ eq = δ 2 eq 1/2 ≈ 0.005.
In general, the total initial angular momentum of the PBH binary is the sum of the angular momenta due to all distant PBHs and large scale density perturbations with characteristic value: The characteristic angular momentum of PBH binaries without DM spikes can be calculated by substituting m b = (m i + m j ) and C(λ) → 1 (in the limit λ ≪ 1) in the above equation.Then, similar to [49], we assume that the probability distribution of the total angular momentum has the same expression as given in Eq. (3.24), but with the characteristic angular momentum given by Eq. (3.29).We find that the impact of DM spikes is larger for higher values of λ, since binaries which decouple late have enough time to accrete massive DM spikes.However, comparing with the case without DM spikes, we find that even for large λ the impact of DM spikes on the angular momentum is not significant.

Merger time of PBH binaries
Since a lot of detailed analysis has already been done on PBH binaries without DM spikes [49,51,52,74,75], our main focus will be on PBH binaries which contain DM spikes.Now, in Sec.4.1 and Sec.4.2 we calculate the final merger time of PBH binaries by taking into account the impact of DM spikes on their merger.We consider two possible scenarios in which the DM spikes can alter the merger dynamics of PBH binaries, listed as follows: 1. Evaporated DM spikes -We assume that for very eccentric orbits, every time the PBHs in the binaries come close to each other there is an exchange of energy between the PBHs and the DM spikes which lead to the ejection of DM particles from the binaries.This continuous loss of energy from the PBHs leads to a decrease in the semi-major axis and circularizes the binaries until the DM spikes are completely evaporated before the merger, as studied in Ref. [56].In the rest of this work, PBH binaries merging after complete evaporation of DM spikes are labelled as PBH binaries with evaporated DM spikes.We define the final orbital parameters of these PBH binaries at the stage where the masses of the DM spikes become zero after which we assume that the binaries lose energy only via emission of gravitational waves.
2. Static DM spikes -We consider PBH binaries with DM spikes in analogy to Intermediate Mass Ratio Inspiral (IMRI) systems in which a lighter object inspirals into another object of intermediate mass.The impact of DM spikes in the orbital evolution of IMRIs is an on-going field of study (see e.g.[60][61][62]).We consider highly eccentric PBH binaries in which the primary PBH of mass m i is at the center of its DM spike and the secondary PBH with mass m j ≪ m i inspirals into the primary PBH.Due to the existence of a strong gravitational potential between the massive primary PBH and its DM spike, a lot of energy is required to remove the DM particles from the spike.For a very light secondary PBH spiralling down into the primary PBH, it is therefore difficult to eject DM particles.Hence, under this assumption, the density of the DM spikes can be treated as static up to the merger of the PBH binary.In rest of work, PBH binaries with such type of spikes as denoted as PBH binaries with static DM spikes.
The actual dynamics of a realistic PBH binary will lie in between these two cases of evaporated DM spikes and static DM spikes.The assumption of evaporated spikes should be an accurate description when the PBH masses are close to equal (as shown in Ref. [56]).The assumption of static spike should be valid in the limit of very large mass ratios.In general, the DM spike of the primary PBH will respond to the inspiral of the secondary PBH but we expect the spike to survive in this process.We will apply these two scenarios regardless of the PBH mass ratio, though we discuss in Sec.4.3 and Sec.6 more precisely when we expect these assumptions to fail and what steps are needed for a more complete treatment.

Evaporated DM spikes
Now, we calculate the final merger time of PBH binaries with evaporated DM spikes by extending the analytical approach of Ref. [56] for extended mass functions and the DM density profile of ρ sp (r) ∝ r −9/4 mentioned in Sec.2.1.In highly eccentric binaries, the orbit is almost radial, due to which the exchange of angular momenta between the PBHs and the DM spikes is negligible [56].Hence, we make the assumption that the angular momentum of the PBHs remains conserved irrespective of the angular momentum of the DM spikes.So, the binaries circularize only due to the loss of the binding energy of DM spikes in every close encounter of PBHs.In that case, the angular momentum of the PBHs in the binary is given as: with µ 0 = m i m j / (m i + m j ) being the reduced mass of the binary in the absence of DM spikes.Then, applying the conservation of angular momentum before and after the evaporation of the DM spikes, we get: with e being the eccentricity of the PBH binary orbit.Now, the binding energy of each DM spike is [56]: where r s is the Schwarzschild radius of the PBH.In the approximation that PBHs with DM spikes can be treated as point objects, the initial orbital energy of the binary with PBH masses m i and m j having DM spikes of masses m sp,i (r) and m sp,j (r) respectively, can be given as: with m total,k = (m k + m sp,k (r)).The final orbital energy of the PBH binary after the complete evaporation of DM spikes is: Then, using the conservation of energy, the final semi-major axis a f of the binary after the DM spikes are totally evaporated [56] can be written as: Here, β is the ratio of the final to initial semi-major axis of the PBH binary, taking into account the complete evaporation of the DM spikes.The variation of the final semi-major axis a f with the initial semi-major axis a i for binaries of dressed PBHs with masses m i and m j is shown in Fig. 3.In all cases, we see that a f remains close to a i for binaries with small initial sizes.Binaries with small values of a i decouple early, meaning that the DM spike has not yet had time to grow.The effect of DM spikes comes into play more and more as the semi-major axis increases, because in that case the binary decouples later.The heavier the DM spikes, the more the energy ejected by the PBH binaries during the evaporation of the spike, which ultimately makes the final separation of the binary comparatively smaller.
In the left panel of this figure (for equal mass PBHs), we see that for a density profile of ρ sp (r) ∝ r −9/4 , the final semi major axis a f of the binaries becomes constant at large values of the initial semi major axis a i .This is because the binding energy of DM spikes tends to a constant value with the growth of the DM spikes.In contrast, the final semi-major axis a f decreases with increase in the initial semi major axis a i for a DM density profile of ρ sp (r) ∝ r −3/2 [56], because the binding energies of the DM spikes keep increasing for this shallower density profile as the initial separation of the binaries increases.In the right panel (for un-equal mass PBHs), we see that as the primary PBH mass m i increases, the mass of its DM spike increases which leads to an increase its binding energy.This enhances the loss of energy from the binary during the evaporation of the DM spike, shrinking the binary more in comparison to other binaries with smaller value of m i .Also, this decrease in the final separation of binaries is much larger for the DM density profile of ρ sp (r) ∝ r −9/4 compared to r −3/2 .Now, the merger time of highly eccentric (e ∼ 1) isolated binaries emitting energy only via emission of gravitational waves is [76]: Then, substituting the value of the final angular momentum given by Eq. (4.2) in the above equation, the final merger time of the PBH binaries with evaporated DM spikes can be written as: with t i being the initial merger time predicted by Eq. (4.7) right after the formation of PBH binary.From the left panel of Fig. 3, it is evident that the final merger time t f of PBH binaries with evaporated DM spikes is smaller than their initially predicted merger time t i .The reason for this behaviour is that while the evaporation of the DM spikes tends to circularize the binary orbit (tending to slow down the merger), it also shrinks the orbit, which tends to speed up the merger, and which dominates over the circularization effect.Hence, PBH binaries with evaporated DM spikes merge faster than the PBH binaries without DM spikes.

Static DM spikes
To find out the merger time of the PBH binaries with static DM spikes, we consider PBH binaries merging analogous to IMRI systems using the same formalism as that of Ref. [62].These PBH binaries lose energy via the emission of gravitational waves (GW) and dynamical friction (DF) over timescales which are much larger than the period of the binary orbit.We assume that the path of the secondary PBH around the central PBH can be considered as a classical Keplerian orbit.During the evolution of this orbit the secondary PBH loses energy and angular momentum leading to merger of the binary.The total orbital energy E is given by: with angular momentum L of the PBHs given by Eq. (4.1).Here, we have assumed that the mass of the DM spikes is much less than the reduced mass of the binary.Then, similar to Ref. [62], using these values of the orbital energy E and angular momentum L, the equation of evolution for the semi-major axis a and eccentricity e of the binary orbit can be respectively written as:  Here, m i is the central/primary PBH and m j is the secondary PBH which spirals down in to m i .The solid part of the curves shows the final merger time t f of these binaries obtained numerically from code imripy [62] and the dashed part shows the power law extrapolation of t f for j i < 10 −3 with power index equal to 0.46.Then, the final merger time of PBH binaries with static DM spikes calculated by solving Eqs.(4.10) and (4.11) using the imripy code7 [62] is shown in Fig. 4. In this figure, the value of the final merger time t f is power-law extrapolated for j i < 10 −3 (i.e.t f ∝ j i ν ) with power-law index ν = 0.46.For a small range of j i , the merger time is approximately constant shown by the flattened part of the curves in Fig. 4.
This behaviour of PBH binaries with static DM spikes shown in Fig. 4 can be explained on the basis of Fig. 5. Here, the terms (1/E) dE/dt and (1/L) dL/dt, which are opposite in sign, contribute to the evolution of the eccentricity given by Eq. (4.11) for PBH binaries with static DM spikes.The rate of energy loss is largest for very eccentric binaries (j ≪ 1).These therefore merge the fastest, as shown in the right panel of Fig. 5.However, in the left panel of Fig. 5, we see that for very small values of initial angular momentum j i , the orbital energy term dominates over the angular momentum term.This dominance as per Eq.(4.11) leads to decrease the eccentricity of the binary with time.This circularization of initially highly eccentric orbits tends to slow their merger compared to the case with a fixed eccentricity.Thus, the merger time grows only slowly with j.In the left panel of this figure, we see that for j i ≳ 0.02, the angular momentum term dominates over the contribution of orbital energy in Eq. (4.11) due to which the eccentricity of the PBH binaries increases with time.In general, the merger time of PBH binaries should increase with increase in j i .However, due to this eccentrification, the increase in the merger time of PBH binaries is smaller than expected.This can be seen in the right panel of Fig. 5 and also explains the flattening of the curves in Fig. 4 close to j i ∼ 0.02.Finally, the merger time of PBH binaries increases quickly as j i → 1, as the process of eccentrification becomes less and less efficient (for j i ∼ 1, this eccentrification occurs only very close to the merger, as seen in the top most curve in the right panel of Fig. 5).The values of the total energy E and angular momentum L of the binary orbit is given by Eq. (4.9) and Eq.(4.1) respectively.We checked that for the PBH binary mentioned in the left panel, the loss of orbital energy due to gravitational waves dE GW /dt is ≈ O(28) smaller than the loss of energy due to dynamical friction dE DF /dt.The right panel shows the variation of angular momentum j of the same PBH binary as a function of time t for different values of the initial angular momentum j i .
By studying the variation of the final merger time t f with initial angular momentum j i shown in Fig. 4, and using the parameters of the reference binary, we find the following approximate expression for final merger time t f of PBH binaries with static DM spikes: where The function f (j i ) is an interpolation function which gives the value of j i corresponding to a given value of t f , following the shape of the curves in Fig. 4. Here, C ref (a i , m i , m j ) is a constant factor which normalizes the final merger time of PBH binaries with static DM spikes with respect to the reference binary defined earlier.In Appendix A, we provide analytic arguments for the scaling of the merger time with the component masses and initial orbital elements of the PBH binaries with static DM spikes.These scalings derive from the fact that in most of the PBH binaries we are interested in, energy losses due to dynamical friction dominate over gravitational wave emission.We note that for PBH binaries with static DM spikes having very small values of initial semi-major axis and initial angular angular momentum, the loss of orbital energy is dominated by the emission of gravitational waves.So, in such cases, we make the assumption that the merger time of PBH binaries with static DM spikes is the same as that of PBH binaries without DM spikes.

Impact of DM spikes on merger time
Now, the variation of final merger time t f of PBH binaries without DM spikes (Eq.(4.7)), with evaporated DM spikes (Eq.(4.8)) and with static DM spikes (Eq.(4.12)) is shown in Fig. 6.In both panels of the figure, we see that the PBH binaries without DM spikes take longest time to merge.The merger time of PBH binaries without DM spikes decreases with increase in m j as t f ≈ 1/m j for m j ≪ m i .Comparing the left and right panels, we see that the merger time of PBH binaries without DM spikes increases as per Eq.(4.7) with increase in the initial angular momentum j i .
In the case of PBH binaries with evaporated DM spikes, the final semi-major axis decreases with decrease in m j i.e. β ∝ m j in Eq. (4.6).This implies that for m j ≪ m i , the loss of energy during the evaporation of DM spikes is much larger than the initial orbital energy of the PBH binaries.Hence, most of the energy of PBH binaries with evaporated DM spikes lies in the form of the binding energy of DM spikes.Then, using β ∝ m j in Eq. (4.8) we get t f ∝ 1/ √ m j due to which the final merger time of PBH binaries with evaporated DM spikes decreases with increase in m j (but more slowly than the case without spikes).
In the case of PBH binaries with static DM spikes, the final merger time decreases with increase in m j as t f ∝ m −0.89 j (see Appendix A).As the value of the initial eccentricity increases then PBH binaries with static DM spikes merge faster which is already explained as per Fig. 4 in the previous section.But, for m i = m j , the merger time of PBH binaries with static DM spikes rises rapidly and coincides with the merger time of PBH binaries without DM spikes.This is because for m i = m j , the value of Coulomb logarithm log Λ = m i /m j in dynamical friction becomes zero.This implies that the loss of orbital energy happens via emission of gravitational waves only due to which its dynamics becomes similar to PBH binaries without DM spikes.Hence, the final merger time of equal mass PBH binaries with static DM spikes coincides with the merger time of equal mass PBH binaries without DM spikes.
In Fig. 6 we also see that there are scenarios in which the system with a static DM spike merges faster than the same system assuming an evaporated spike.The loss of energy which tends to accelerate the merger comes from the binding energy of the DM spike.Since all the DM particles get ejected before the merger, so the evaporation of DM spikes can be interpreted as the maximum possible amount of energy which can be extracted from the PBH binaries.Hence, we would expect that PBH binaries with evaporated DM spikes should have the minimum possible merger time.When systems with static spikes merge faster than this, it suggests that more energy is being extracted by dynamical friction than is available in the binding energy of the spike, signalling the breakdown of the static spike formalism described above.As we have noted already, the actual physics of PBH binaries and spikes lies somewhere in between these two scenarios and it is the ratio of PBHs masses q = m j /m i which plays the key role to decide which description of the dynamics is closer to the truth.In particular, from Fig. 6, we see that the merger time for binaries with static spikes becomes increasingly smaller than the merger time with evaporated spikes, as q → 1, matching our intuition that the static spike assumption should be poor for equal-mass binaries.Of course, in the general case, we require a more complete model to describe the dynamics of dressed PBH binaries which includes the feedback of DM spikes in the picture.We address this issue in more detail in Sec. 6.

Merger rates of PBH binaries
We assume that the PBHs are uniformly distributed in the Universe and there is no other PBH lying inside the comoving volume V = 4π 3 x 3 occupied by the binary with PBHs of masses m i and m j .We note that if f pbh ≈ 1, then the binaries can get trapped in a PBH cluster.The dense environment of the PBH cluster can disrupt the binaries leading to the suppression of their current merger rates [26,77,78]. 8We neglect this suppression factor in our calculations by avoiding f pbh ≈ 1 because we are focusing on subdominant PBHs.Moreover, an abundance of f pbh ≈ 1 for Solar-mass PBHs is disfavoured by LVK estimates of BBH merger rates.The dimensionless variable λ given by Eq. (3.10) quantifies the initial separation of the PBHs in the binary.Then, the probability distribution of λ for a PBH binary takes the form [52]: with n T = f ρ eq ∞ 0 m dm as the total number density of PBHs.Here, P (λ) is the probability of finding a pair of nearest neighbouring PBHs of mass m i and m j with a dimensionless separation λ ∝ x 3 , with no other PBHs within the comoving distance x.Now, we assume that the combined probability distribution function P (λ, j i ) is separable i.e.P (λ, j i ) = P (λ) • P (j i ).Then, the combined probability distribution of dimensionless variable λ and final merger time t f of PBH binaries is: Substituting the value of P (j i ) given by Eq. (3.24) in the above equation, the probability distribution of the final merger time t f of PBH binaries at a given value of λ becomes: (5.4) Then, using Eqs.(5.1) and 5.3, the combined probability distribution of (λ, t f ) can be written as: (5.5) Using Eq. ( 5.2), the distribution of the final merger time can be calculated as: Then, similar to Ref. [52], for PBH binaries having PBHs of masses m i and m j , the merger rate at time t f is equal to the number of mergers N merger per unit merger time t f per unit comoving volume V m , written as: where ρ 0 m ≃ 4 × 10 19 M ⊙ Gpc −3 is the matter density of the Universe today.The final merger time used in our further calculations is equal to the age of the Universe today i.e. t f = 13.78Gyr.Here, R ij is the merger rate of PBH binaries having extended PBH mass functions but for monochromatic PBH mass functions, R ij can be calculated using m i = m j = m i.e.P (m i ) = P (m j ) = P (m)/2.

Evaporated DM spikes
For PBH binaries with evaporated DM spikes, the dependence of angular momentum on merger time is given by Eq. (4.7) which along with Eq. (4.2) gives: Substituting the above equation in Eq. (5.5), the probability distribution of (λ, t f ) for binaries with evaporated DM spikes becomes: We note that the functional form of P (λ, t f ) is the same for PBH binaries with evaporated DM spikes and without DM spikes (though in that case the relationship between j i and t f is different).Now using Eq.(5.6), the probability distribution of the final merger time for PBH binaries with evaporated DM spikes is: Substituting the value of µ given by Eq. (3.5) and probability distribution of final merger time given by Eq. (5.10) in Eq. (5.7), the merger rate of PBH binaries with evaporated DM spikes at time t f is: with Here, is a dimensionless constant and c is the speed of light.The merger rate R ij, 0 of PBH binaries without DM spikes [49,52] at time t f can be calculated by setting the mass of DM spikes equal to zero and β = 1 in the merger rate of PBH binaries with evaporated DM spikes, given by Eq. (5.11).

Static DM spikes
The probability distribution of (λ, t f ) for PBH binaries having static DM spikes can be calculated using Eq.(4.12) in Eq. (5.5) as: From Eq. (5.6), the probability distribution of final time of merger for PBH binaries with static DM spikes can be written as: (5.15) Again, substituting the value of µ given by Eq. (3.5) and the probability distribution of the final merger time given by Eq. (5.15) in Eq. (5.7), the merger rate of PBH binaries with static DM spikes is given as:

Component mass dependence of merger rates
To calculate the merger rates of PBH binaries with evaporated and static DM spikes given by Eq. (5.11) and Eq.(5.16) respectively, we consider three specific examples of PBH mass functions (MFs): Lognormal [79], Power-Law [52,80] and Multipeak [65] illustrated in Fig. 7.The details of these three MFs are as follows: • Multipeak in mass range of 10 −4 − 100 M ⊙ .This MF is extracted from a piece-wise primordial power spectrum (PPS) which is in agreement with all the current constraints on the abundance of PBHs in the Universe [65].The variation of the threshold for PBH formation with the equation of state in the early Universe leads to the multipeak structure of this MF [81].
In the rest of the work, we shall denote the lower mass limit of these MFs as m min and the upper mass limit as m max .
For the Lognormal MF, we fix f pbh = 1.61 × 10 −3 which corresponds to a merger rate of R = 44 Gpc −3 yr −1 for PBH binaries without DM spikes and roughly matches with the merger rate observed by the LVK Collaboration.Then, we use this value of f pbh to calculate the merger rates of PBH binaries with and without DM spikes having Lognormal MF.For the Power Law and Multipeak MFs, we use f pbh = 2.53 × 10 −2 and f pbh = 0.066 respectively which are in agreement with all the constraints on PBHs having extended mass functions.To simplify the structure of the paper, we shall focus more on the Power Law and Multipeak MFs in the main text and mention some of the further details about the Lognormal MF in Appendix B. In the following, we will present the merger rates as a function of the component masses m i and m j .In addition, in Appendix C, we present the same information as a function of the mass ratio q = m j /m i .This is particularly instructive for understanding whether mergers with large or small values of q dominate the present day merger rate and therefore which formalism (evaporated or static spikes) might be more applicable, following the discussion in Sec.4.3.Now, the present-day merger rates of PBH binaries without DM spikes are shown in Fig. 8.We see that for both the MFs, the merger rates are largest for binaries composed of the most abundance PBHs (compare with Fig. 7).If the mass fraction of PBHs P (m)/m were constant with m, the merger rate would be expected to decrease roughly proportional to 1/m [51].However, the bin widths ∆ we have used in Fig. 8 grow with m.Thus, the binary population and therefore the binned merger rates simply trace the underlying abundance of PBHs with masses m i and m j , as expected from Eq. (5.7).
In Fig. 9, we show the ratio of the present-day merger rates of PBH binaries with evaporated DM spikes given by Eq. (5.11) and binaries without DM spikes.In this figure, we see that for our Power Law and Multipeak MF examples, the presence of evaporated DM spikes may enhance or reduce the present-day merger rate by an O(1) factor, depending on the component masses.
This behaviour can be explained on the basis of how the initial properties of PBH binaries merging today is affected by the presence of the DM spike.In Fig. 10, we show the probability distribution of the final merger time t f as a function of the dimensionless PBH separation λ, i.e.P (t f |λ) given by Eq. (5.3), for two different choices of component PBH masses.The black dashed line shows the case without DM spikes.The peak of this distribution represents the value of λ for binaries which are most likely to merge today, which we denote λ peak .Both panels of Fig. 10 show that in comparison to the PBH binaries without DM spikes, the value of λ peak is shifted towards higher values in the presence of evaporated DM spikes (shown by the orange colored curves).This can be explained as per Eq.(4.8) which signifies that for fixed values of a i ∼ λ 4/3 and initial angular momentum j i , the presence of evaporated DM spikes reduces the final merger time of the PBH binaries.Hence, one to λ ∼ λ peak .The shift to larger values of λ leads to the availability of more binaries merging today (without much penalty from larger values of λ being rarer).In this case, then, the merger rate is increased by the presence of DM spikes.
In order to distinguish these two possibilities, we can rewrite P (λ), given by Eq. (5.1), as: with the typical value of λ: Here, V avg = 1/n T is the average comoving volume occupied by any PBH in the Universe (regardless of its mass) and V ij = 4π 3 x3 ij is the average comoving volume occupied by a pair of PBHs of mass m i and m j .As the mass of the component PBHs increases, their number density decreases and V ij increases, leading to a decrease in λ typ .Then, Eq. (5.17) signifies that the formation of PBH binaries is pushed to smaller values of λ < λ typ (as in the left panel of Fig. 10).Physically, as V ij becomes larger, the probability that another PBH (with mass m k ̸ = m i , m j ) will be in the intervening volume increases (preventing the formation of the m i − m j binary).This in turn pushes the typical separation of binaries to smaller values.We should then compare λ typ to λ peak to understand whether the merger rate increases or decreases with the addition of DM spikes.The value of λ peak can be estimated by fixing the initial angular momentum of the binary to its most likely value j i ∼ j λ ∝ λ, given by Eq. (3.25), and choosing then λ = λ peak to give the desired merger time.If λ peak > λ typ , then P (λ) decreases rapidly close to λ peak , as in the left panel of Fig. 10, leading to a decrease in the present-day merger rate in the presence of DM spikes.Instead, If λ peak < λ typ , then P (λ) is roughly constant close to λ peak , as in the right panel of Fig. 10, and the merger rate increases.
In Fig. 11, we show the ratios of the current merger rates of PBH binaries with static DM spikes, given by Eq. (5.16), and binaries without DM spikes.From this figure, we see that the merger rates for most of the PBH binaries with static DM spikes are massively suppressed for Power Law MF.This behaviour can be explained as in the case of evaporated spikes, though the addition of static DM spikes typically reduced the merger time substantially more than static DM spikes (see e.g.Fig. 6).This pushes λ peak to much higher values, meaning that binaries wide enough to be observed merging today are very rare (λ peak ≫ λ typ ).This in turn leads to a huge suppression of the merger rate.For the Multipeak MF, we note that the presence of DM spikes enhances the merger rate for PBH binaries composed of the lightest PBHs (bottom left corner).For a fixed mass fraction, very light PBHs have a much larger number density, decreasing the typical volume V ij in between PBHs.In addition, light PBHs are the most abundant in the Multipeak MF, meaning that the probability of finding heavier PBHs in the volume V ij is small.This means that the typical dimensionless separation λ typ is large (λ peak < λ typ ).The shift in λ peak to larger values increases the availability of binaries merging today, which in turn can enhance the merger rates for these light, abundance PBH binaries.
We note that if the DM density profile is shallower (for example, ρ sp (r) ∝ r −3/2 ) the behaviour of the PBH binaries is expected to be similar to the case with ρ sp (r) ∝ r −9/4 .That is, DM spikes of density profile of ρ sp (r) ∝ r −3/2 should also lead to an increase in the merger  The typical value of λ given by Eq. (5.18) decreases as we increase f pbh .9However, as we increase f pbh , the most probable value of λ peak for binaries merging today also decreases.This can be seen by noting that the distribution P (t f |λ), given by Eq. (5.1), peaks at γ λ = √ 2. The typical initial angular momentum of PBH binaries scales as j λ ∼ λ(f 2 + σ 2 eq ) 1/2 , so setting j i ∼ √ 2j λ and fixing λ peak choosing the final merger time, we obtain the estimate: where we have neglected some O(1) factors inside the brackets.Here, the value of the exponent κ for different PBH binaries can be calculated as: for no DM spikes, 14 21 for evaporated DM spikes, 2 for static DM spikes. (5.20) The different values of κ arise because of the different scaling of the merger time with j i and λ in each scenario.The presence of evaporated spikes causes large binaries to shrink as they eject DM, meaning that the merger time scales more slowly with the initial semi-major axis (and therefore λ) than in the absence of DM spikes.This translates to a smaller value of κ.
In contrast, in the static spike case, the inspiral is driven by dynamical friction.The merger time then scales much more slowly with j i than in the absence of DM spikes.This translates to a larger value of κ.Thus, λ peak decreases as we increase f pbh , but the precise scaling with f pbh varies depending on our assumptions about the presence and behaviour of DM spikes.We see from the above discussion that as we decrease f pbh , the value of λ peak grows more slowly for binaries without DM spikes than for binaries with evaporated DM spikes.This means that the relative shift in λ peak induced by evaporated spikes grows as we decrease f pbh .In the mass range of interest in Fig. 12, we are typically in the scenario where λ typ > λ peak , leading to an increase in the merger rate as we decrease f pbh .However, the impact of evaporated DM spikes on the merger rate grows only slowly (the two values of κ are similar in Eq. (5.20)) and eventually saturates as f pbh ≪ σ eq .In contrast, as we decrease f pbh , the value of λ peak grows more quickly for binaries without DM spikes than for binaries with static DM spikes.At large values of f pbh , we find that λ peak,static ≪ λ typ , strongly suppressing the merger rate.However, as we decrease f pbh , λ peak in the static case becomes closer to the case without DM spikes.Eventually, for f pbh ≪ σ eq , we reach the point where λ peak,static > λ typ and the merger rate can be substantially enhanced.
In Fig. 12, we also highlight with a grey band the merger rates of BH binaries detected by the LVK Collaboration in mass range of 5 − 100 M ⊙ : R = (17.9− 44) Gpc −3 yr −1 , as reported in GWTC-3 [82].We note that these rates can be well explained by the mergers of PBH binaries with and without DM spikes having either Lognormal or Power Law MF (though of course the contribution of astrophysical binaries should also be taken into account).To emphasise the effect that DM spikes would have on an inference of the PBH abundance, we list the values of f pbh which would be compatible with the reported LVK merger rates in Table 1.

Conclusions
In this work, we have reviewed the accretion of dark matter (DM) spikes around isolated primordial black holes (PBHs).Then we studied the effect of the growth of DM spikes on the initial orbital parameters of PBH binaries formed in the early Universe.To probe the merger dynamics of PBH binaries, we made assumptions about the ways in which these spikes can impact the evolution of PBH binaries.We have calculated the current merger rates of PBH binaries with DM spikes for three different extended PBH mass functions under the assumption that either the DM spikes are completely evaporated before merger ("evaporated" spikes) or they remain static up to the merger ("static" spikes).
Our calculations verify that for the density profile of ρ(r) ∝ r −9/4 , the mass of the DM spike being accreted around an isolated PBH is directly proportional to the mass of the PBH itself.We have demonstrated for the first time that the impact of DM spikes on the formation of PBH binaries is independent of the masses of the component PBHs.We find that the effect of DM spikes on the orbital parameters increases with the scale factor at which the formation of the binary takes place, as this allows more time for the DM spike to grow.The impact of DM spikes on the initial semi-major axis can be substantial for binaries forming around matter-radiation equality or later, but the impact on the initial binary eccentricity is typically negligible.
Our results for the merger rates of PBH binaries suggest a subtle behaviour for extended mass functions.We found that in the scenario of evaporated DM spikes, the merger rates can be increased or decreased by around a factor of two.This is due to the fact that the presence of DM spikes speeds up the merger of PBH binaries.This in turn means that in order to find a PBH binary which merges today, the initial PBH separation required is shifted to larger values compared to the case without DM spikes.If this required separation is larger than the typical PBH separation, this corresponds to a shift to increasingly rare PBH binaries, leading to a reduction in the merger rate.Instead, if this required separation is smaller than the typical separation, this shift leads to an increase in the availability of binaries merging today, leading to an increase in the merger rate.The same logic applies when the DM spikes are considered static until merger, though in that case the merger rates of PBH binaries can be reduced by many orders of magnitude.This is especially true for very asymmetric binaries, which potentially has implications for searches for dressed PBH IMRIs and EMRIs.These results strongly disagree with [56,57] in the sense that the presence of DM spikes around PBHs might only leads to increase the merger rates of PBH binaries.We would also like to point out that the focus of our numerical results was on the present-day merger rates of PBH binaries.However, our formalism still holds for PBH binaries merging at arbitrary redshift, and can be used to study the merger rates over cosmic time.
Our calculations suggest that the current data of the LVK Collaborations could include a contribution from the mergers of PBH binaries having DM spikes.As shown in Table 1, this constrains the fraction of PBHs in CDM f pbh ≤ O(10 −5 − 10 −3 ), depending on the PBH mass function, consistent with previous results.We have also highlighted how the gap between the merger rates of PBH binaries with and without DM spikes increases with decrease in f pbh .This is also an important aspect which has not been considered so far and must be taken into account in order to derive accurate limits from future data of the LVK Collaboration, as existing constraints are pushed to smaller and smaller values of f pbh .In addition, these results may also be relevant for studying the rapid and enhanced mergers of PBH binaries in the early Universe due to the presence of DM spikes.Such PBHs may provide the seeds for the growth of supermassive black holes and early-forming galaxies [11,12].
In our description of DM spikes around PBHs, we have ignored the annihilation, mutual gravitational attraction and growth of DM spikes after decoupling of PBH binaries.We have also ignored processes such as the capture of binaries in PBH clusters or the interaction of PBHs lying very close to the binaries.We have argued that for sufficiently small f pbh some of these should be negligible (such as the formation of PBH clusters).However, a more complete picture of these factors should be taken into account for careful analysis of the accretion of DM spikes around PBHs.
In this work, we have considered two extreme scenarios for the behaviour of DM spikes around PBHs binaries.We expect our description of evaporated spikes to be accurate for binaries with close-to-equal mass ratios, in which the injection of energy leads to the rapid unbinding of the spike.Instead, the static spike scenario may be more accurate when the mass ratio is very large and the lighter PBH does not substantially perturb the spike of the larger companion.Of course, the true dynamics will be somewhere in between these extremes.This true dynamics can be very complex, including the co-evolution of the binaries and the DM spike.To accurately probe the merger rates of PBH binaries dressed with DM spikes, the feedback of the PBHs on DM spikes also need to be considered in a way similar to Ref. [61].
Even so, our results highlight that the presence of DM spikes around PBH binaries may enhance or suppress their merger rates, depending strongly on their abundance and mass functions.In some cases, our results suggest that this suppression can be extreme.The development of a prescription for describing the feedback of DM spikes in highly eccentric systems should therefore be considered a priority.With this, it will be possible to develop a more general and complete approach which can probe the real merger dynamics of the PBH binaries dressed with DM spikes.section, numerical calculations using imripy [62] suggested that the final merger time t f should scale with the initial semi-major axis, a i and the component masses, m i and m j as: with α = 0.75, γ = 0.65 and δ = −0.89.We begin by noting that for large enough initial separations, gravitational wave energy losses are typically subdominant to dynamical friction energy losses (see e.g.Ref. [63]).This roughly matches the empirical scaling observed in Eq. (A.1), with small differences coming from the dependence of log Λ on m i and m j and on the effects of eccentricity in the orbit.Lognormal mass function

Figure 1 :
Figure 1: Dynamics of different DM shells around a PBH of mass m pbh = 100 M ⊙ from deep radiation domination up to matter-radiation equality, t eq .The distance r of all the DM shells are measured from the center of the PBH.In the left panel, the blue colored region represents the DM spike around the PBH consisting of various DM shells shown by colored concentric circles.In the right panel, the colored lines show the evolution of different DM shells around the PBH, in accordance with Eq. (2.1).The black dotted line shows the turnaround radii, r ta of DM shells given by Eq. (2.2) which grows with time.

Figure 2 :
Figure2: The left panel shows the dynamics of PBH binaries with and without DM spikes up to the first cross-over of the PBHs using λ = 1 in Eq.(3.8).The right panel shows the variation of the initial semi-major axis a i of the PBH binaries with and without DM spikes as a function of λ as per the numerical solution of Eq. (3.8).Here, λ is the dimensionless parameter given by (3.10) which is a proxy for the PBH separation and which quantifies how much larger the background density is compared to the mutual density of the two PBHs.In addition, x is the comoving separation of the PBHs in the binaries and s dec is the scale factor of the binary decoupling.

. 15 )
Now, consider a dressed PBH lying outside the binary, designated as the distant PBH.The total mass of the distant PBH is m total,d (s) = m d + m sp,d (s) with m d being the mass of the PBH itself and m sp,d as the mass of its DM spike at scale factor s. The tidal torque will grow as the total mass of the distant PBH grows with s.So, if the comoving separation y of the distant PBH from the center of the PBH binary is approximately constant and much larger than the separation of the binary then T ∝ m total,d (s)/s 3 .The tidal field exerted on the PBH binary at matter-radiation equality by the distant PBH at comoving separation y ≫ x is[49]:

10 − 8 Figure 3 :
Figure3: Variation of the final semi-major axis a f with initial semi major axis a i for PBH binaries with evaporated DM spikes.In both panels, the black dotted line shows the case in which the final semi-major axis is equal to the initial semi-major axis i.e. a f = a i .The colored solid lines represent results for the density profile of ρ sp (r) ∝ r −9/4 given by Eq. (4.6) and the colored dashed lines correspond to ρ sp (r) ∝ r −3/2[56,71].
/dt signifies the sum of the loss of energy due to the emission of gravitational waves dE GW /dt and due to dynamical friction dE DF /dt which are further discussed in Appendix A. Now, we define a reference PBH binary with static DM spikes having m i,ref = 1 M ⊙ , m j,ref = 10 −3 M ⊙ and a i,ref = 1.67 × 10 2 pc.This value of a i,ref is calculated using Eq.(3.13) assuming that the comoving separation x of the PBHs is equal to its average value i.e. x = x.10 −6 10 −5 10 −4 10 −3 10 −2 10 −1 10 0

Figure 4 :
Figure4: Variation of the final merger time t f with initial angular momentum j i for different PBH binaries having static DM spikes.Here, m i is the central/primary PBH and m j is the secondary PBH which spirals down in to m i .The solid part of the curves shows the final merger time t f of these binaries obtained numerically from code imripy[62] and the dashed part shows the power law extrapolation of t f for j i < 10 −3 with power index equal to 0.46.

Figure 5 :
Figure5: The left panel shows the variation of terms (1/E) dE/dt (dark magenta line) and (1/L) dL/dt (olive colored line) with angular momentum j, mentioned in the equation of evolution of the eccentricity e given by Eq. (4.11), for a PBH binary with static DM spikes.The values of the total energy E and angular momentum L of the binary orbit is given by Eq. (4.9) and Eq.(4.1) respectively.We checked that for the PBH binary mentioned in the left panel, the loss of orbital energy due to gravitational waves dE GW /dt is ≈ O(28) smaller than the loss of energy due to dynamical friction dE DF /dt.The right panel shows the variation of angular momentum j of the same PBH binary as a function of time t for different values of the initial angular momentum j i .

10 −4 10 − 3 10 −2 10 −Figure 6 :
Figure 6: Variation of final merger time t f as a function of the PBH mass m j for PBH binaries with evaporated DM spikes and static DM spikes.The left and right panels show different values of the initial angular momentum j i of the binary.Here, we have chosen the values of initial semi-major axis as a i = 0.08 pc.

Figure 7 :
Figure 7: Variation of Lognormal, Power Law and Multipeak mass functions with PBH mass m pbh in mass range from m min to m max .

Figure 8 :
Figure 8: Current merger rates of PBH binaries without DM spikes R ij, 0 for Power Law and Multipeak MFs with PBH mass range from m min to m max .

Figure 12 :
Figure 12: Variation of current merger rates of PBH binaries with and without DM spikes as a function of f pbh .Here, f pbh = f /0.85 is the fractional abundance of PBHs in cold dark matter.The shaded region corresponds to the latest merger range of (17.9 − 44) Gpc −3 yr −1 seen by LVK Collaboration for BBHs mergers in mass range of 5 − 100 M ⊙ such that m i ≥ m j .

) 7/ 2 . 4 i a − 9 / 4 6 )
(A.4)    with µ 0 = m i m j / (m i + m j ) as the reduced mass of the PBH binary in the absence of DM spikes.The rate of loss of energy of the PBH binary due to the dynamical friction is[62]:dE DF dt (r, v) = 4π (Gm j ) 2 ρ sp (r)v −1 log Λ ,(A.5)such that Coulomb logarithm log Λ = log m i /m j .Counting powers of the masses and semi-major axis, we find that ρ sp ∼ m 3/and v ∼ (m i /a) 1/2 .We thus find that:We can now estimate the scaling of the typical final merger time as:

Figure 13 shows
Figure13shows the current merger rates of PBH binaries with and without DM spikes in the mass range of 5 − 100 M ⊙ for our Lognormal benchmark mass function.The top left panel signifies that in the absence of DM spikes the merger rates are for PBH binaries which consist of the most abundant PBHs as shown in Fig.7.The top right panel shows that the presence of evaporated DM spikes increases the merger rate of PBH binaries by roughly a factor of 2

Figure 13 :
Figure 13: Current merger rates of PBH with and without DM spikes having Lognormal mass function with PBH mass range from m min to m max .Here, R ij, 0 represents the current merger rate of PBH binaries having no DM spikes around the PBHs.R ij, evaporated denotes the current merger rate of PBH binaries consisting of DM spikes which get completely evaporated before merger.And R ij, static signifies the current merger rate of PBH binaries containing DM spikes which remain static until the merger.

Table 1 :
Values of f pbh constrained by latest LVK merger data as per the current merger rates of PBH binaries with and without DM spikes shown in Fig.12.
−4 − 1.10 × 10 −3 [83]e the orbital energy is E = −G m i m j /2 a.The loss of total energy E of the binary orbit occurs via the emission of gravitational waves (GW) and due to dynamical friction (DF), Similar to Ref.[83], the rate of loss of orbital energy of the binary due to the emission of gravitational waves is given as: (m i + m j ) 3 1 +7325 e 2 + 37 96 e 4 a 5 (1 − e 2