We review the theory behind the formation of primordial black hole binaries and their merger rates. We consider the binary formation in the early and late Universe emphasising the former as it gives the dominant contribution of the present primordial black hole merger rate. The binaries formed in the early Universe are highly eccentric so their interactions with other primordial black holes can significantly increase their coalescence times and thereby suppress the merger rate. We discuss in detail how the suppression of the merger rate arising from such interactions can be estimated and how such interactions lead to the formation of another much harder binary population that contributes to the present merger rate if more than 10 of dark matter consists of primordial black holes with a relatively narrow mass distribution. When the primordial abundance is below 1 encounters between primordial black holes are unlikely and their effect on the merger rate becomes negligible.
['Astrophysics - Cosmology and Nongalactic Astrophysics', 'Astrophysics - High Energy Astrophysical Phenomena', 'High Energy Physics - Phenomenology']
Formation of primordial black hole binaries and their merger rates
{'Formation of primordial black hole binaries and their merger rates': "Martti Raidal, Ville Vaskonen and Hardi Veermae  \nAbstract We review the theory behind the formation of primordial black hole binaries and their merger rates. We consider the binary formation in the early and late Universe, emphasising the former as it gives the dominant contribution of the present primordial black hole merger rate. The binaries formed in the early Universe are highly eccentric so their interactions with other primordial black holes can significantly increase their coalescence times and thereby suppress the merger rate. We discuss in detail how the suppression of the merger rate arising from such interactions can be estimated and how such interactions lead to the formation of another, much harder, binary population that contributes to the present merger rate if more than 10% of dark matter consists of primordial black holes with a relatively narrow mass distribution. When the primordial abundance is below 1%, encounters between primordial black holes are unlikely and their effect on the merger rate becomes negligible.  \nMartti Raidal  \nKeemilise ja Bioloogilise Fuusika Instituut, Ravala pst. 10, 10143 Tallinn, Estonia e-mail: martti.raidal@cern.ch  \nVille Vaskonen  \nKeemilise ja Bioloogilise Fuusika Instituut, Ravala pst. 10, 10143 Tallinn, Estonia, Dipartimento di Fisica e Astronomia, Universit'a degli Studi di Padova, Via Marzolo 8, 35131 Padova, Italy  \nIstituto Nazionale di Fisica Nucleare, Sezione di Padova, Via Marzolo 8, 35131 Padova, Italy, e-mail: ville.vaskonen@pd.infn.it  \nHardi Veerm¨ae  \nKeemilise ja Bioloogilise F¨u¨usika Instituut, R¨avala pst. 10, 10143 Tallinn, Estonia  \ne-mail: hardi.veermae@cern.ch", '1 Introduction': 'Gravitational wave (GW) signals created by primordial black hole (PBH) binaries provide a direct probe of hypothetical PBH populations. The PBH populations are expected to differ from astrophysical black hole (ABH) populations in aspects such as mass distribution, redshift dependence of the merger rate, spins or even spatial distribution. A good theoretical understanding of both populations is crucial for employing the full potential of GW observations. In the following, we will focus on describing the PBH binary populations and the resulting PBH merger rates.  \nThe characteristics of the present-day PBH populations are moderately modeldependent. They can depend on the numerous ways how the PBHs could have been produced in the early Universe and potentially also on deviations from the standard cosmology, which can affect the evolution of the PBH populations after their formation. To our luck, it is not necessary to specify the details of PBH formation to estimate the properties of the present-day PBH binary populations. Only the PBH abundance, the distribution of their masses and spins as well as their spatial distribution before matter-radiation equality should be specified. This sets the initial stage of PBH structure formation as nearby PBH pairs decouple from the expansion and form the first bound structures in the Universe.  \nIn certain scenarios, a population of merging PBH binaries can be probed by GWexperiments. With enough data, these mergers can help us reconstruct the initial characteristics of PBHs in the early Universe and begin probing the early Universe by reconstructing scenarios that would produce such initial conditions. The goal of this chapter is to outline the first part of this exercise and show how an initial population of PBHs coupled to the Hubble flow in the very early Universe will give rise to a population of PBH binaries merging in the present day.', '1.1 Binary formation channels': 'There are various channels for forming PBH binaries. These can be classified by 1) whether the binaries were formed in the early or late Universe, and 2) whether the formation included dynamics of two or more bodies. The PBH binary characteristics of each combination can be outlined as follows:  \n- 1. In the early two-body channel , binaries are formed from nearby PBH pairs that decouple the Hubble flow. Such binaries generally dominate the PBH binary merger rate when 𝑓 PBH - the fraction of DM in PBHs - is small. However, they are very eccentric since most of their angular momentum comes from tiny tidal forces induced by fluctuations in the density of surrounding matter. This makes the merger rate susceptible to perturbations by encounters with other PBHs, so their subsequent evolution should be accounted for when estimating their present merger rate. In particular, when 𝑓 PBH ≳ 0 . 1, small-scale structure formation is enhanced. This leads to greater disruption of these binaries and suppresses the corresponding merger rate.\n- 2. In the early three-body channel , bound three-body systems are formed after initial compact 3-PBH configurations decouple from the Hubble flow. In such a system, a binary remains if one of the PBHs is ejected or if two of the PBHs merge. The resulting binaries are much less eccentric and, therefore, much harder /one.sup than the binaries from the early two-body channel. As a result, their merger rate is not significantly affected by encounters with other PBHs. However, since compact 3-PBH systems are less likely to occur in the early Universe, this channel is generally subdominant to the early two-body channel for small 𝑓 PBH. However, it can dominate if 𝑓 PBH ≈ 1 or if the PBHs are formed in clusters.\n- 3. In the late two-body channel , PBH binaries are formed from hyperbolic encounters of PBHs in DM haloes if they lose a sufficient amount of energy via GWemission. The resulting binaries are typically eccentric and merge relatively fast. Although the resulting merger rate depends on the PBH structure formation and may be boosted by spiky DM profiles or PBH clustering, such binaries tend to have a smaller merger rate than PBH binaries formed in the early Universe. Although this channel could also produce ABH-PBH binaries, their contribution is expected to be marginal [1].\n- 4. The late three-body channel can be active in the environment with large PBH densities. The binaries are formed in 2+1 PBH encounters, where the 3rd PBH carries away enough energy for the remaining 2 to become bound. Analogously to the early three-body channel, these binaries are hard and not too eccentric. Although their formation implies an environment with relatively frequent PBHPBHbinary encounters, such encounters will not significantly affect their merger rate. This channel is generally weaker than the ones listed above. It can be  \nenhanced if 𝑓 PBH ≈ 1 or if the formation of PBH structures is enhanced, i.e., by initial clustering.  \nIn the following, we will consider each channel in detail.', '1.2 Preliminaries': "We will use geometric units 𝑐 = 𝐺 = 1 throughout this chapter. We consider a generic PBH mass function, that we define as  \n𝜓 ( 𝑚 ) ≡ 𝑚 𝜌 PBH d 𝑛 PBH d ln 𝑚 , (1)  \nwhere 𝜌 PBH and 𝑛 PBH denote the total energy and comoving number densities of PBHs. With this definition we get ∫ 𝜓 ( 𝑚 ) d ln 𝑚 = 1 and the average of a quantity 𝑋 is given by  \n⟨ 𝑋 ⟩ ≡ ∫ d 𝑛 PBH 𝑛 PBH 𝑋 = ∫ d ln 𝑚 ⟨ 𝑚 ⟩ 𝑚 𝜓 ( 𝑚 ) 𝑋 , (2)  \nwhere ⟨ 𝑚 ⟩ = 𝜌 PBH / 𝑛 PBH = 1 / ∫ d ln 𝑚𝑚 -1 𝜓 ( 𝑚 ) .  \nWedenote the masses of individual PBHs by a lowercase 𝑚 and the PBHs making up the binary by 𝑚 1 and 𝑚 2. The total mass and the reduced mass of the binary are 𝑀 ≡ 𝑚 1 + 𝑚 2 and 𝜇 ≡ 𝑚 1 𝑚 2 / 𝑀 and we often use the symmetric mass ratio 𝜂 ≡ 𝜇 / 𝑀 to quantify unequal mass binaries.  \nWe work in the Newtonian limit. Thus, the binding energy of the binary is  \n𝐸 bin = 𝜇 / r 2 2 -𝑀𝜇 𝑟 = 𝑀𝜇 2 𝑟 𝑎 , (3)  \nwhere 𝑟 𝑎 denotes the semimajor axis of the binary, r = r 1 -r 2 is the separation of the PBHs positioned at r 1 and r 2, and the dot denotes time derivation. The binary's angular momentum L = 𝜇 r × / r will be quantified via the dimensionless quantity  \nj ≡ L / 𝜇 √ 𝑟 𝑎 𝑀 , (4)  \nwhich satisfies 𝑗 ≡ | j | ≤ 1 and is related to the eccentricity as  \n𝑒 = √︃ 1 -𝑗 2 . (5)  \nThus, for our purposes, it is sufficient to characterize PBH binaries by four quantities, e.g., 𝑚 1, 𝑚 2, 𝑟 𝑎 and 𝑗 .  \nTo estimate the merger rate, we need to consider the formation of PBH binaries and their disruption between formation and merger. Given the PBH binary formation rate density d 𝑅 b ( 𝑡 ' )/( d 𝑚 1d 𝑚 2d 𝑗 d 𝑟 𝑎 ) at time 𝑡 ' , the merger rate density at time 𝑡 is  \nFormation of primordial black hole binaries and their merger rates  \nd 𝑅 ( 𝑡 ) d 𝑚 1d 𝑚 2 = ∫ d 𝑡 ' d 𝑗 d 𝑟 𝑎 d 𝑅 b ( 𝑡 ' ) d 𝑚 1d 𝑚 2d 𝑗 d 𝑟 𝑎 𝛿 ( 𝑡 -𝑡 ' -𝜏 ( 𝑚 1 , 𝑚 2 , 𝑟 𝑎 , 𝑗 )) , (6)  \nwhere 𝜏 denotes the coalescence time. The coalescence time can be estimated assuming the binary will evolve only via GW emission after its formation. In the Newtonian regime, the coalescence time of eccentric binaries is [2]  \n𝜏 = 3 85 𝑟 4 𝑎 𝑗 7 𝜂𝑀 3 . (7)  \nThis expression would overestimate the coalescence time mostly by a factor of 1.81 when applied to almost circular binaries. /two.sup If the binary is disrupted /three.sup , or evolves between formation and merger due to effects other than GW emission, e.g., accretion, then (6) must be modified accordingly.  \n𝜏 = 3 𝑟 4 𝑎 𝑗 7 85 𝜂𝑀 3 𝐹 ( 𝑗 ) -𝐹 ( 1 ) 𝑗 GLYPH<16> 1 -121 425 𝑗 2 GLYPH<17> 3480 / 2299 GLYPH<0> 1 -𝑗 2 GLYPH<1> 24 / 19 ≈ 3 𝑟 4 𝑎 𝑗 7 85 𝜂𝑀 3 1 -0 . 70 𝑗 + 0 . 62 𝑗 2 1 + 0 . 67 𝑗 , (8)  \nwhere 𝐹 ( 𝑗 ) ≡ 𝐹 1 GLYPH<0> -1 / 2; -1181 / 2299 , -5 / 19; 1 / 2; 𝑗 2 121 / 425 , 𝑗 2 GLYPH<1> with 𝐹 1 the first Appell hypergeometric function. The second approximation deviated from the exact one by less than 0.5%. The limiting case of circular binaries, 𝑗 → 1, gives 𝜏 = ( 5 / 256 ) 𝑟 4 𝑎 / 𝜂𝑀 3 .", '2 PBH binary formation in the early Universe': 'In the early Universe, PBH binaries can form when nearby PBHs decouple from the Hubble flow [3]. The formation process and the resulting merger rate depend on whether a third PBH falls into the pair. In this section, we will first lay out the derivation for the two-body formation channel, and look at phenomena that may perturb these binaries and suppress this rate. Finally, we will consider the merger rate stemming from the three-body channel.', '2.1 Dynamics of early binary formation': 'To set up the initial configuration for the early Universe two-body channel, let us consider a pair of PBHs with masses 𝑚 1 and 𝑚 2 at a comoving separation 𝑥 0 which initially follow the Hubble flow. We assume that there are no PBHs within a comoving radius 𝑦 of the pair, to guarantee that the pair evolves as a two-body system after decoupling. This initial setup is shown in Fig. 1. The size of the empty region will be determined later by requiring that the PBHs closest to the pair will not become gravitationally bound to the pair after decoupling from the Hubble flow and, therefore, form a bound system consisting of 3 or more bodies. Such systems will be addressed later in Sec. 2.3. We further assume that the positions of PBHs are uncorrelated, that is, they follow a Poisson distribution. In this case, the comoving number density of PBH configurations described above is  \nd 𝑛 pairs = 1 2 𝑒 -¯ 𝑁 ( 𝑦 ) d 𝑛 ( 𝑚 1 ) d 𝑛 ( 𝑚 2 ) d 𝑉 ( 𝑥 0 ) , (9)  \nwhere d 𝑛 ( 𝑚 ) is the comoving number density of PBHs in the mass range ( 𝑚, 𝑚 + d 𝑚 ) , and  \n¯ 𝑁 ( 𝑦 ) ≡ 𝑛𝑉 ( 𝑦 ) (10)  \nis the expected number of PBHs in a spherical volume 𝑉 ( 𝑦 ) ≡ ( 4 𝜋 / 3 ) 𝑦 3 with a comoving radius 𝑦 . The factor 1 / 2 avoids overcounting.  \nThe population of PBH binaries created from such initial conditions is described by the formation rate  \nd 𝑅 b ( 𝑡 ) d 𝑚 1d 𝑚 2d 𝑗 d 𝑟 𝑎 = ∫ d 𝑥 0 d 𝑛 p d 𝑚 1d 𝑚 2d 𝑥 0 d 𝑃 ( 𝑗 , 𝑟 𝑎 | 𝑥 0 , 𝑦 ) d 𝑗 d 𝑟 𝑎 𝛿 ( 𝑡 -𝑡 dc ( 𝑥 0 )) , (11)  \nwhere 𝑡 dc ( 𝑥 0 ) denotes the time of decoupling of the pair from the Hubble flow and d 𝑃 ( 𝑗 , 𝑟 𝑎 | 𝑥 0 , 𝑦 )/ d 𝑗 d 𝑟 𝑎 is the distribution orbital parameters given 𝑥 0 and 𝑦 . /four.sup A few simplifications can be made. Since the early binaries are formed before the matterradiation equality, we can omit the dependence on formation time and assume that  \nFig. 1 A initial PBH configuration and its parametrization.  \n<!-- image -->  \n𝑡 dc ( 𝑥 0 ) ≈ 𝑡 0, where 𝑡 0 is the age of the Universe. /five.sup Additionally, 𝑃 ( 𝑗 , 𝑟 𝑎 | 𝑥 0 , 𝑦 ) ∝ 𝛿 ( 𝑟 𝑎 -𝑟 𝑎 ( 𝑥 0 )) d 𝑃 ( 𝑗 | 𝑥 0 , 𝑦 )/ d 𝑗 since 𝑟 𝑎 can be predicted from 𝑥 0 alone. Given the coalescence time (7), we find that the merger rate density (6) can be recast as  \nd 𝑅 ( 𝑦 ) E2 ( 𝑡 ) = ∫ d 𝑛 pairsd 𝑗 d 𝑃 ( 𝑗 | 𝑥 0 , 𝑦 ) d 𝑗 𝛿 ( 𝑡 -𝜏 ( 𝑟 𝑎 ( 𝑥 0 ) , 𝑗 )) = 1 14 𝜏 d 𝑛 ( 𝑚 1 ) d 𝑛 ( 𝑚 2 ) ∫ d 𝑉 ( 𝑥 0 ) 𝑒 -¯ 𝑁 ( 𝑦 ) 𝑗 d 𝑃 ( 𝑗 | 𝑥 0 , 𝑦 ) d 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> 𝑗 = 𝑗 ( 𝜏,𝑟 𝑎 ( 𝑥 0 ) ) , (12)  \nwhere, in the second line, 𝑗 ( 𝜏, 𝑟 𝑎 ( 𝑥 0 )) is obtained by inverting Eq. (7).  \nThis estimate must, however, be modified to account for the disruption of the early binaries via collisions of the binary with other PBHs in DM clusters between the binary formation and merger. We will describe this effect by the survival probability 𝑆 L, which we take to be independent of the initial conditions. Furthermore, the exclusion of initial conditions in which the binary gets disrupted by nearby PBHs shortly after formation can be characterised by another survival probability 𝑆 E. So, the merger rate density Eq. (12) takes the final form  \nd 𝑅 E2 ( 𝑡 ) = 𝑆 L 𝑆 E × d 𝑅 ( 0 ) E2 ( 𝑡 ) , (13)  \nwhere 𝑆 E ≡ d 𝑅 ( 𝑦 ) E2 / d 𝑅 ( 0 ) E2 ≤ 1. /six.sup Note that 𝑦 is generally not an independent variable and can depend on various model parameters.  \nAll in all, to proceed, we need to compute ( i ) how the semimajor axis 𝑟 𝑎 ( 𝑥 0 ) is formed during the decoupling from the Hubble flow (see section 2.1.1), ( ii ) how spatial inhomogeneities set the distribution of initial angular momenta of the binary (see section 2.1.2), ( iii ) what is the minimal distance 𝑦 between the pair and other PBHs that would avoid disrupting the binary due to the infall of other close PBHs into the binary and how does it affect 𝑆 E (see section 2.1.3), ( iv ) estimate disruption  \nin late DM haloes by computing the probability 𝑆 L that the binaries do not collide with other PBHs during before the merger (see section 2.2).', '2.1.1 Decoupling': "The line element for spacetimes with small inhomogeneities in the comoving Newtonian gauge and the absence of anisotropic stress is  \nd 𝑠 2 = -( 1 + 2 𝜙 ( x )) d 𝑡 2 + ( 1 -2 𝜙 ( x )) 𝑎 ( 𝑡 ) 2 d 𝑥 2 , (14)  \nwhere the potential 𝜙 is determined by the total matter energy density 𝜌 M ( x ) as Δ 𝜙 = 4 𝜋𝑎 2 𝜌 𝑀 ( x ) . Treating the PBHs as point masses, we obtain  \n𝜙 ( x ) = -∑︁ 𝑖 𝑚 𝑖 𝑎 | x 𝑖 -x | -4 𝜋 ¯ 𝜌 M ∫ d 3 𝑘 ( 2 𝜋 ) 3 𝑎 2 𝑘 2 𝑒 -𝑖 k · x 𝛿 M ( k ) , (15)  \nwhere 𝑚 𝑖 denote the masses and x 𝑖 the comoving positions of the PBHs, ¯ 𝜌 M the average comoving matter energy density and 𝛿 M ( k ) the matter density perturbation. A system of 𝑁 non-relativistic PBHs, 𝑎 d 𝑥 / d 𝑡 ≪ 1, is described by the action  \n𝑆 ( 𝑁 ) = ∫ d 𝑡 GLYPH<20> ∑︁ 𝑖 𝑚 𝑖 GLYPH<18> 1 2 / r 2 𝑖 + 1 2 ¥ 𝑎 𝑎 r 2 𝑖 -𝜙 ex ( x ) GLYPH<19> + ∑︁ 𝑖>𝑗 𝑚 𝑖 𝑚 𝑗 | r 𝑖 -r 𝑗 | GLYPH<21> , (16)  \nwhere r 𝑖 ≡ 𝑎 x 𝑖 denotes the proper distance, the last term accounts for the pairwise interaction of the PBHs and the potential 𝜙 ex ( x ) describes the effect of the surrounding matter on the 𝑁 -body system. For a PBH pair, the action can be approximated as  \n𝑆 ( 2 ) ≈ ∫ d 𝑡 GLYPH<20> 𝑀 2 GLYPH<18> / r 2 𝑐 + ¥ 𝑎 𝑎 r 2 𝑐 GLYPH<19> -𝑀𝜙 ex ( r 𝑐 ) + 𝜇 2 GLYPH<18> / r 2 + ¥ 𝑎 𝑎 r 2 + 2 𝑀 𝑟 -r · T · r GLYPH<19> GLYPH<21> , (17)  \nwhere r denotes the proper separation of the PBH, r 𝑐 the centre of mass of the PBH pair, and T 𝑖 𝑗 ≡ 𝜕 𝑖 𝜕 𝑗 𝜙 ex ( r 𝑐 ) is the tidal tensor due to the external forces. In summary, the PBH pair is subject to the following forces  \nF / 𝜇 = r ¥ 𝑎 / 𝑎 |{z} Hubble flow -𝑀 ˆ r / 𝑟 2 | {z } self-gravity + ( ˆ r · T · r ) ˆ r | {z } radial tidal forces + ( r × ( T · r )) | {z } tidal torque ×( ˆ r / 𝑟 ) , (18)  \nwhere ˆ r ≡ r / 𝑟 . The first three forces are radial, whereas the last term provides the torque that prevents the head-on collision of the PBH pair.  \nThe angular momentum L of the pair is generated by the tidal torque,  \nL = ∫ d 𝑡 r × F = 𝜇 ∫ d 𝑡 r × ( T · r ) . (19)  \nDuring radiation domination, the PBHs surrounding the pair are assumed to follow the Hubble flow and the matter density perturbations are constant. This implies that the potential (15) scales roughly as 𝜙 ∝ 𝑎 -1 . Consequently, T ∝ 𝑎 -3 and the tidal forces get damped fast after the pair decouples from the Hubble flow. The pair has then formed a binary with angular momentum L and semi-major axis 𝑟 𝑎 .  \nThebinary decouples when the Hubble flow becomes subdominant and the second term in (18) takes over. To estimate when this happens, consider the overdensity in matter due to the PBH pair of comoving size 𝑥 0,  \n𝛿 pair ≡ 𝑀 / 2 𝜌 M 𝑉 ( 𝑥 0 ) , (20)  \nwhich should exceed the average matter density 𝛿 pair ≫ 1 if the PBH pair is to form a binary. Such a region begins to collapse when its density becomes comparable to the density of radiation, 𝜌 R 𝑎 -4 ≈ 𝛿 pair 𝜌 M 𝑎 -3 , that is when  \n𝑎 ≈ 𝑎 dc ≡ 𝑎 eq / 𝛿 pair , (21)  \nwhere 𝑎 eq denotes the scale factor during matter-radiation equality. Since 𝛿 pair ≫ 1, we expect the binary to be formed during radiation domination. To make this more precise, consider first the idealized case of a head-on collision for which, the radial equation of motion of the pair, ¥ 𝑟 -𝑟 ¥ 𝑎 / 𝑎 + 𝑀𝑟 -2 = 0, can be recast in terms of the comoving separation 𝑥 ≡ 𝑟 / 𝑎 as  \n𝑥 '' + 𝑎 𝑎 dc 𝑥 3 0 𝑥 2 = 0 , (22)  \nwhere the prime denotes derivation with respect to ln ( 𝑎 ) and the initial conditions are given by 𝑥 ( 𝑎 0 ) = 𝑥 0, 𝑥 ' ( 𝑎 0 ) = 0 with 𝑎 0 ≪ 𝑎 dc. The solution of the equation of motion can be expressed as 𝑥 ( 𝑎 ) ≈ 𝑥 0 𝜒 ( 𝑎 / 𝑎 dc ) , where 𝜒 ( 𝑦 ) solves ( 𝑦𝜕 𝑦 ) 2 𝜒 + 𝑦𝜒 -2 = 0 with 𝜒 ( 𝑦 → 0 ) = 1, 𝜒 ' ( 𝑦 → 0 ) = 0. The function 𝜒 ( 𝑦 ) oscillates with an amplitude that decreases asymptotically as 0 . 2 / 𝑦 , so the semi-major axis of the binary is  \n𝑟 𝑎 = 𝑎𝑥 𝑎 / 2 ≈ 0 . 1 𝑎 dc 𝑥 0 ≈ 0 . 84 𝜌 R 𝑥 4 0 / 𝑀. (23)  \nThe first root of 𝜒 ( 𝑦 ) lies at 𝑦 ≈ 0 . 54 corresponding to the first close encounter of the pair at 𝑎 ≈ 0 . 54 𝑎 dc and indicating that the binary decouples from the Hubble flow before that. In Fig. 2 we indicate by the dashed vertical line the moment when the self-gravity of the pair starts to dominate over the Hubble flow.  \nThe idealised case discussed above is in excellent agreement with the radial evolution in the presence of a weak tidal field illustrated in Fig. 2 as well as with 𝑁 -body simulations in which the tidal effects are created by dynamical PBHs surrounding the pair [5] suggesting that the angular momentum can be studied perturbatively as long as the orbit remains eccentric ( 𝑗 ≪ 1) [6]. Plugging the solution of (22) into (19), we find that the dimensionless angular momentum (4) is  \n)  \na  \nr  \n2  \n(  \n/  \nr  \nFig. 2 Decoupling of the PBH pair. The left and right panels show, respectively, the evolution of the pair's separation and angular momentum in a weak tidal field. The vertical dashed lines show the moment when the self-gravity of the pair starts to dominate over the Hubble flow (see Eq. (18)). The horizontal dashed line in the right panel shows the approximation in Eq. (24).  \n<!-- image -->  \na  \n/  \na  \n0.1  \n1  \ndc  \nj = 1 √ 𝑟 𝑎 𝑀 ∫ d 𝑎 𝑎 2 𝐻 x × ( T 𝑎 3 · x ) = 0 . 95 𝑥 3 0 𝑀 ˆ r × ( T 𝑎 3 · ˆ r ) , (24)  \nwhere we used ∫ ∞ 0 d 𝑦 𝜒 ( 𝑦 ) 2 ≈ 0 . 3. Note that under the above assumptions T 𝑎 3 is constant.", '2.1.2 Distribution of angular momenta': 'The tidal torque is determined by the distribution of surrounding PBHs as well as the inhomogeneities of the rest of matter. The resulting dimensionless angular momentum can thus be expressed as  \nj = j PBH + j M , j PBH ≡ ∑︁ 𝑖 j 1 ( x 𝑖 , 𝑚 𝑖 ) , (25)  \nwhere j 1 is the angular momentum generated by single PBH with mass 𝑚 𝑖 at a comoving distance 𝑥 𝑖 from the pair, and j M is the angular momentum generated by matter fluctuations. Using Eq. (24) with T 𝑖 = 𝑚 𝑖 ( 1 -3ˆ x 𝑖 ⊗ ˆ x 𝑖 )/ 𝑥 3 𝑖 𝑎 -3 and defining  \n𝑗 0 ≡ 0 . 95 ¯ 𝑁 ( 𝑥 0 ) ⟨ 𝑚 ⟩ 𝑀 ≈ 0 . 4 𝑓 PBH 𝛿 pair , (26)  \nwhere ⟨ 𝑚 ⟩ = 𝜌 PBH / 𝑛 is the average PBH mass and we have used Ω M / Ω DM = 1 . 2, we get  \nj 1 ( x 𝑖 , 𝑚 𝑖 ) = 𝑗 0 𝑚 𝑖 ⟨ 𝑚 ⟩ 3 ¯ 𝑁 ( 𝑥 𝑖 ) ˆ x 𝑖 × ˆ r ( ˆ x 𝑖 · ˆ r ) . (27)  \nThe quantity 𝑗 0 is a reference scale for the dimensionless angular momentum. It captures its order of magnitude because torque is dominantly generated by the closest PBHs, for which ¯ 𝑁 ( 𝑥 𝑖 ) ≈ 1, and the remaining terms in Eq. (27) will contribute an  \n1.0  \n0.8  \n0.6  \n0.4  \n0.2  \n0.0  \n-  \n2  \n10  \nO( 1 ) on average. Binaries form from PBH pairs which satisfy 𝛿 pair ≫ 1, thus we naturally expect that 𝑗 0 ≪ 1.  \nTo compute the distribution of angular momenta, we will consider the most common scenario in which the masses 𝑚 𝑖 and positions x 𝑖 as well as the matter density perturbations are statistically independent. In this case, it is convenient to use the additive property of the cumulant-generating function of j  \n𝐾 ( k ) ≡ ln GLYPH<10> 𝑒 𝑖 k · j GLYPH<11> = 𝐾 M ( k ) + 𝐾 PBH ( k ) , (28)  \nwhere 𝐾 M ( k ) ≡ ln GLYPH<10> 𝑒 𝑖 k · j M GLYPH<11> and 𝐾 PBH ( k ) ≡ ln GLYPH<10> 𝑒 𝑖 k · j PBH GLYPH<11> are cumulant generating functions of j M, j PBH, respectively. The average is taken over all random quantities contained in the expression. The probability distribution of j can then be computed as  \nd 𝑃 d 3 𝑗 ≡ ⟨ 𝛿 ( j -j M -j PBH )⟩ = ∫ d 3 𝑘 ( 2 𝜋 ) 3 𝑒 -𝑖 k · j + 𝐾 ( k ) . (29)  \nNote that the above distribution is two-dimensional since the angular momentum is perpendicular to ˆ r : the cumulant generating function is a function of k ⊥ ≡ k × ˆ r . Statistical isotropy of the surrounding matter further constrains it to be a function of 𝑘 ⊥ only.  \nLet us first consider the matter density fluctuations. The generating function 𝐾 M ( k ) is determined by the first two cumulants if the matter density fluctuations are Gaussian, that is  \n𝐾 M ( k ) = -1 2 GLYPH<10> ( k · j ) 2 GLYPH<11> = -1 4 𝜎 2 𝑗, M 𝑘 2 ⊥ , (30)  \nsince the mean ⟨ j ⟩ = 0 must vanish due to isotropy arguments. The variance 𝜎 𝑗, M can be estimated by noting that matter density fluctuations induce a tidal field which in Fourier space reads T M = 𝑎 -3 ˆ q ⊗ ˆ q 4 𝜋 ¯ 𝜌 M 𝛿 M ( 𝑞 ) by Eq. (15). Using Eq. (24) and averaging over orientations then gives [5, 6]  \n𝜎 2 𝑗, M = GLYPH<10> j 2 M GLYPH<11> = 0 . 95 𝑥 3 0 𝑀 ! 2 𝑎 6 5 GLYPH<28> tr ( T M · T M ) -1 3 tr ( T M ) 2 GLYPH<29> = 6 5 𝑗 2 0 𝜎 2 M 𝑓 2 PBH , (31)  \nwhere 𝜎 2 M ≡ ( Ω M / Ω DM ) 2 GLYPH<10> 𝛿 2 M GLYPH<11> is the rescaled variance of matter density perturbations at the time the binary is formed. /seven.sup  \nThe cumulant generating function 𝐾 PBH accounts for identical and statistically independent contributions from all surrounding PBHs. To send the number of surrounding PBHs to infinity, we will first consider a finite volume and then take the limit 𝑉 →∞ while keeping 𝑁 / 𝑉 = 𝑛 fixed: /eight.sup  \n𝐾 PBH ( k ) = ∫ d 𝑛 ( 𝑚 ) ∫ | x | >𝑦 d Ω 4 𝜋 d 𝑉 ( 𝑥 ) GLYPH<20> exp GLYPH<18> 𝑖 3 𝑚𝑗 0 𝜌 PBH 𝑉 ( 𝑥 ) ( k ⊥ · ˆ x ) ( ˆ x · ˆ r ) GLYPH<19> -1 GLYPH<21> . (32)  \n𝐾 PBH ( k ) = lim 𝑉 →∞ 𝑁 / 𝑉 = 𝑛 𝑁 ln GLYPH<18> ∫ | x | >𝑦 d 3 𝑥 𝑉 d 𝑛 ( 𝑚 ) 𝑛 𝑒 𝑖 k · j 1 ( x ,𝑚 ) GLYPH<19> = ∫ d 3 𝑥 d 𝑛 ( 𝑚 ) GLYPH<16> 𝑒 𝑖 k · j 1 ( x ,𝑚 ) -1 GLYPH<17> = -¯ 𝑁 ( 𝑦 ) ∫ d 𝑛 ( 𝑚 ) 𝑛 𝐹 GLYPH<18> 𝑚 ⟨ 𝑚 ⟩ 1 ¯ 𝑁 ( 𝑦 ) 𝑗 0 𝑘 ⊥ GLYPH<19> , (33)  \nwhere 𝐹 ( 𝑧 ) = 1 𝐹 2 GLYPH<0> -1 / 2; 3 / 4 , 5 / 4; -9 𝑧 2 / 16 GLYPH<1> -1 and 1 𝐹 2 is the generalised hypergeometric function. It asymptotes to 𝐹 ( 𝑧 ) ∼ 𝑧 -1 when 𝑧 →∞ and 𝐹 ( 𝑧 ) ∼ 3 𝑧 2 / 10 when 𝑧 → 0. These limits correspond to ¯ 𝑁 ( 𝑦 ) ≪ 1 and ¯ 𝑁 ( 𝑦 ) ≫ 1, respectively. In particular, the ¯ 𝑁 ( 𝑦 ) ≫ 1 case recovers the Gaussian limit 𝐾 PBH ( k ) ∼ -1 4 𝜎 2 𝑗, PBH 𝑘 2 ⊥ , with the variance given by  \n𝜎 2 𝑗, PBH = 6 5 𝑗 2 0 ⟨ 𝑚 2 ⟩ ⟨ 𝑚 ⟩ 2 1 ¯ 𝑁 ( 𝑦 ) . (34)  \nThe variance diverges when ¯ 𝑁 ( 𝑦 ) → 0. For this reason, it is not possible to use the central limit theorem to compute the angular momentum distribution even though it arises from an infinite sum of independent identical random variables. We also remark that a simple parallel can be drawn between the variances (31) and (34) by noting that the variance of PBH density fluctuations in a volume 𝑉 ( 𝑦 ) is ⟨ 𝑚 2 ⟩/(⟨ 𝑚 ⟩ 2 ¯ 𝑁 ( 𝑦 )) and has thus an analogous interpretation to the factor 𝜎 2 M / 𝑓 PBH in Eq. (31).  \nDue to isotropy, it is sufficient to consider the distribution of 𝑗 instead of j . The cumulant generating functions imply the angular momentum distribution (29) is given by  \n𝑗 d 𝑃 d 𝑗 = ∫ ∞ 0 d 𝑢 𝑢𝐽 0 ( 𝑢 ) exp GLYPH<20> -¯ 𝑁 ( 𝑦 ) ∫ d 𝑚 𝑚 𝜓 ( 𝑚 ) 𝐹 GLYPH<18> 𝑢 𝑚 ⟨ 𝑚 ⟩ 1 ¯ 𝑁 ( 𝑦 ) 𝑗 0 𝑗 GLYPH<19> -𝑢 2 3 10 𝜎 2 M 𝑓 2 PBH 𝑗 2 0 𝑗 2 GLYPH<21> . (35)  \nDepending on the size of the empty region surrounding the PBH pair, this distribution interpolates between two asymptotic cases  \nThe spatial integral can be evaluated by making the change of variables 𝑢 ≡ 𝑧𝑉 ( 𝑦 )/ 𝑉 ( 𝑥 ) , with 𝑧 ≡ 𝑘 ⊥ 𝑗 0 𝑚 /( 𝜌 PBH 𝑉 ( 𝑦 ) ) , and taking the following steps:  \n-𝑧 ∫ d 2 Ω 4 𝜋 ∫ 𝑧 0 d 𝑢 𝑢 2 h exp GLYPH<16> 𝑖 3 𝑢 ( ˆ k ⊥ · ˆ x ) ( ˆ x · ˆ r ) GLYPH<17> -1 i = -𝑧 ∫ 𝑧 0 d 𝑢 𝑢 2 ∫ 𝜋 0 ∫ 2 𝜋 0 d cos ( 𝜃 ) d 𝜙 4 𝜋 GLYPH<20> exp GLYPH<18> 𝑖 3 𝑢 2 sin 2 ( 𝜃 ) sin ( 2 𝜙 ) GLYPH<19> -1 GLYPH<21> = -𝑧 ∫ 𝑧 0 d 𝑢 𝑢 2 GLYPH<20> 𝜋 2 √ 2 𝐽 -1 4 GLYPH<18> 3 𝑢 4 GLYPH<19> 𝐽 1 4 GLYPH<18> 3 𝑢 4 GLYPH<19> -1 GLYPH<21> = 1 𝐹 2 GLYPH<18> -1 2 ; 3 4 , 5 4 ; -9 𝑧 2 16 GLYPH<19> -1 .  \nFig. 3 Angular momentum distributions of PBH binaries after formation (35) for different ¯ 𝑁 ( 𝑦 ) are shown by solid lines. The dashed lines show the corresponding Gaussian approximation (36). A monochromatic mass function is assumed in the figure.  \n<!-- image -->  \n𝑗 d 𝑃 d 𝑗 =              𝑗 2 / 𝑗 2 0 ( 1 + 𝑗 2 / 𝑗 2 0 ) 3 / 2 , ¯ 𝑁 ( 𝑦 ) ≪ 1 (other PBHs arbitrarily close) , 2 𝑗 2 𝜎 2 𝑗 𝑒 -𝑗 2 / 𝜎 2 𝑗 , ¯ 𝑁 ( 𝑦 ) ≫ 1 (Gaussian limit) , (36)  \nwhere  \n𝜎 2 𝑗 ≡ 𝜎 2 𝑗, M + 𝜎 2 𝑗, PBH = 6 5 𝑗 2 0 1 ¯ 𝑁 ( 𝑦 ) ⟨ 𝑚 2 ⟩ ⟨ 𝑚 ⟩ 2 + 𝜎 2 M 𝑓 2 PBH ! (37)  \nis the total variance of the distribution. The angular momentum distribution (35) together with the Gaussian approximation is shown in Fig. 3 for a monochromatic mass function. It can be seen that the Gaussian approximation works well already for ¯ 𝑁 ( 𝑦 ) = 2, while the ¯ 𝑁 ( 𝑦 ) → 0 limit becomes a good approximation when ¯ 𝑁 ( 𝑦 ) ≲ 0 . 2.  \nThe appearance of 𝜓 ( 𝑚 ) in (35) implies a non-trivial dependence on the mass function. However, the dependence on 𝜓 ( 𝑚 ) vanishes completely (36) in the ¯ 𝑁 ( 𝑦 ) ≪ 1 limit and can be reduced to the relative width of the mass distribution in the Gaussian ¯ 𝑁 ( 𝑦 ) ≫ 1 case.  \nEq. (35) is in excellent agreement with numerical simulations as illustrated by Fig. 4. It shows the angular momentum distributions obtained from 𝑁 -body simulations of PBH binaries expected to merge within the age of the Universe 𝑡 0 (green region) and compares them to the distribution predicted by Eq. (35) (black dashed curve). Note that the distribution (35) gives a conditional probability assuming a fixed initial separation 𝑥 0. The corresponding distribution conditioned on a fixed coalescence time 𝑡 is  \nd 𝑃 ( 𝑗 | 𝜏 = 𝑡 ) d 𝑗 ∝ 𝑗 0 ( 𝑥 0 ) d 𝑃 ( 𝑗 | 𝑥 0 ) d 𝑗 GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> 𝜏 ( 𝑟 𝑎 ( 𝑥 0 ) , 𝑗 ) = 𝑡 , (38)  \n1  \nFig. 4 Distribution of angular momenta at 𝑎 = 3 𝑎 eq of PBH binaries for different 𝑓 PBH and a monochromatic mass function with 𝑚 𝑐 = 30 𝑀 ⊙ . The initial configurations were chosen in a way that the PBH pair would form a binary that merges today if unperturbed. The perturbed and unperturbed binaries are coloured in yellow and green respectively. The thick dashed line shows the distribution (38) with ¯ 𝑁 ( 𝑦 ) estimated from Eq. (40). The thin dashed curves show limiting cases. The thick grey dashed curve compares the distribution of disrupted binaries to d 𝑃 / d 𝑗 ∝ 1. Based on simulation data from Ref. [5].  \n<!-- image -->  \nwhere we dropped prefactors independent of 𝑗 . Analogously to the distribution 𝑗 d 𝑃 ( 𝑗 | 𝑥 0 )/ d 𝑗 that depends on 𝑗 only through the combination 𝑗 / 𝑗 0, it is possible to define 𝑗 𝜏 so that d 𝑃 ( 𝑗 | 𝜏 )/ d 𝑗 would depend on 𝑗 through the combination 𝑗 / 𝑗 𝜏 only, where 𝑗 𝜏 gives an order of magnitude estimate of the dimensionless angular momentum of initial binaries merging today. Since 𝑗 / 𝑗 0 | 𝜏 = 𝑡 = ( 𝑗 / 𝑗 𝜏 ) 37 / 16 , we have that  \n𝑗 𝜏 = 1 . 7 × 10 -2 GLYPH<20> 𝑡 𝑡 0 GLYPH<21> 3 37 GLYPH<20> 𝑀 𝑀 ⊙ GLYPH<21> 5 37 ( 4 𝜂 ) 3 37 𝑓 16 37 PBH . (39)  \nThe dependence on both the mass and the coalescence time is very mild, so it is safe to assume that PBHs around the stellar mass range had initial eccentricities around O( 10 -2 ) or lower when 𝑓 PBH ≪ 1.', '2.1.3 Merger rate and binary disruption by nearest PBHs': 'The binary may be disrupted shortly after formation by becoming bound to a nearby PBHor later by interacting with PBHs in compact DM haloes, which are more likely to form in PBH cosmologies. The first possibility is accounted by an appropriate choice of ¯ 𝑁 ( 𝑦 ) and it is typically the dominant effect when 𝑓 PBH ≪ 1. When 𝑓 PBH ≈ 1, however, the binary is likely to be perturbed by later encounters due to the enhanced small-scale structure of PBH DM. This effect will be addressed separately in Sec. 2.2.  \nSince 𝑗 0 = O( 10 -2 ) for binaries merging today, the early PBH binaries tend to be extremely eccentric at formation. Since the coalescence time scales as 𝜏 ∝ 𝑗 7 , even a O( 1 ) increase in 𝑗 will increase the coalescence time by several orders of magnitude. Therefore, even mild interactions can effectively remove the binary from the population of binaries merging today. It is certainly possible that some of the dis-  \nFig. 5 Binary disruption at 𝑎 = 3 𝑎 eq depending on the initial comoving distance 𝑥 NN of the PBH nearest to the binary for 30 𝑀 ⊙ PBHs and different 𝑓 PBH. The blue region indicates simulations in which the initial pair did not form a binary. The yellow region shows pairs that evolved into a binary that was later perturbed but remained bound. Binaries above the yellow dashed line swapped at least one of their components. The green region shows undisrupted binaries. The vertical dashed line indicates the estimate Eq. (40). The numerical simulations evolve a PBH pair estimated to evolve into a binary that would merge at the present time if not disrupted. Based on simulation data from Ref. [5].  \n<!-- image -->  \ns may contribute to the present merger rate, but they must evolve from initially more compact PBH systems to merge within the age of the Universe after being perturbed. This scenario constitutes an entirely different formation channel the early three-body channel - discussed in Sec. 2.3.  \nTo determine ¯ 𝑁 ( 𝑦 ) , consider collisions between the binary and the PBH initially closest to it. Let 𝑥 NN denote the initial comoving distance of the nearest neighbour. Since 𝑦 sets a lower bound on viable values of 𝑥 NN for binary formation, we must find the minimal 𝑥 NN which does not lead to the disruption of the binary after the nearest neighbour decouples from the Hubble flow. This estimate can be made following the reasoning of Sec. 2.1.1. The region surrounding the binary corresponds to an effective matter density fluctuation 𝛿 NN ≈ ( 𝑀 -𝜌 PBH 𝑉 ( 𝑥 NN ))/( 𝜌 M 𝑉 ( 𝑥 NN )) . Given that such configurations are expected to collapse at 𝑎 ≈ 𝑎 eq / 𝛿 NN we estimate that, at a given redshift, the binary has collided with its neighbour if 𝑥 NN < 𝑦 , when  \n¯ 𝑁 ( 𝑦 ) ≈ 𝑀 ⟨ 𝑚 ⟩ 𝑓 PBH 𝑓 PBH + 𝑎 eq / 𝑎 . (40)  \nThe condition (40) agrees well with numerical 𝑁 -body simulations of PBH binary formation in the early Universe [5]. This is illustrated in Fig. 5, which shows that (40) predicts the interface between perturbed (yellow region) and unperturbed (green region) binaries very well when 𝑓 PBH ≪ 1. When 𝑓 PBH = 1, about half of the binaries with 𝑥 NN > 𝑦 are seen disrupted. As we will show in Sec. 2.2, the early formation of PBH clusters in such PBH cosmologies [7] can account for this effect. The estimate Eq. (40) was also used when comparing numerically obtained angular momentum distributions to our analytic predictions in Fig. 4. In that case, an excellent agreement can be observed also when 𝑓 PBH = 1.  \nWhen 𝑎 ≫ 𝑎 eq and for narrow mass functions, we have that ¯ 𝑁 ( 𝑦 ) ≈ 2 implying that the PBH pair should be positioned in an underdense region. This condition is, however, too strong when the PBH abundance is small 𝑓 PBH ≪ 1. In this case, it is sufficient to demand that the binary survives until the non-linear structure formation of CDM haloes sets in. Since this happens when 𝑎 / 𝑎 eq = 1 / 𝜎 𝑀 , we have that  \n¯ 𝑁 ( 𝑦 ) ≈ 𝑀 ⟨ 𝑚 ⟩ 𝑓 PBH 𝑓 PBH + 𝜎 M . (41)  \nThe binary is expected to evolve without much interference after that point since the disruption of PBH binaries in the standard late CDM haloes is very unlikely [6]. This point is illustrated further in section 2.2, where we show that, even when considering the enhanced small-scale structures in PBH cosmologies, PBH encounters are expected to take place in relatively small haloes containing less than O( 10 4 ) PBHs which are formed at high redshifts.  \nA potential shortcoming of the condition (41) is that it accounts for the masses of the surrounding PBHs on average. Thus it can fail when the mass distribution of PBHs spans several orders of magnitude. For example, it can underestimate the suppression when the nearest PBH is much heavier than the total mass of the binary since heavier PBHs are more likely to disrupt light binaries even if their distance is large. On the other hand, it can overestimate the disruption of heavy binaries by light PBHs surrounding it. In extreme cases, the surrounding binaries can be too light to significantly alter the orbit of the binary. In that case, excluding initial configurations where a light PBH is close to binary is unjustified.  \nHaving determined the characteristics of the first binaries and constraints to their environment, we can proceed to evaluate the merger rate density (13) d 𝑅 𝐸 2 = 𝑆 E 𝑆 Ld 𝑅 ( 0 ) 𝐸 2 . The unsuppressed contribution d 𝑅 ( 0 ) 𝐸 2 is evaluated by integrating Eq. (12) in the 𝜎 𝑀 → 0 limit and using 𝑦 → 0 asymptotic (35) of the angular momentum distribution,  \nd 𝑅 ( 0 ) 𝐸 2 ( 𝑡 ) d 𝑚 1d 𝑚 2 = 1 . 4 𝑡 𝑡𝜂𝑀 14 𝑓 7 PBH 𝜌 11 𝑀 ! 3 37 d 𝑛 ( 𝑚 1 ) d 𝑛 ( 𝑚 2 ) ≈ 1 . 6 × 10 6 Gpc 3 yr 𝑓 53 37 PBH 𝜂 -34 37 GLYPH<18> 𝑀 𝑀 ⊙ GLYPH<19> -32 37 GLYPH<18> 𝑡 𝑡 0 GLYPH<19> -34 37 𝜓 ( 𝑚 1 ) 𝜓 ( 𝑚 2 ) 𝑚 1 𝑚 2 . (42)  \nThe suppression factor 𝑆 E gives the remaining contribution /nine.sup  \n𝑆 E = 𝑒 -¯ 𝑁 ( 𝑦 ) Γ ( 21 / 37 ) ∫ d 𝑣 𝑣 -16 37 exp GLYPH<20> -¯ 𝑁 ( 𝑦 )⟨ 𝑚 ⟩ ∫ d 𝑚 𝑚 𝜓 ( 𝑚 ) 𝐹 GLYPH<18> 𝑚 ⟨ 𝑚 ⟩ 𝑣 ¯ 𝑁 ( 𝑦 ) GLYPH<19> -3 𝜎 2 M 𝑣 2 10 𝑓 2 PBH GLYPH<21> , (43)  \nwith ¯ 𝑁 ( 𝑦 ) estimated from (41).  \nThe suppression factor includes the effect of matter perturbations. In particular, when 𝑓 PBH ≪ 𝜎 𝑀 , the angular momentum of the early binaries is generated mainly by matter perturbations (compare Eqs. (31) and (34)) instead of the surrounding PBHs. As a result, the merger rate scales 𝑓 2 PBH instead of 𝑓 53 / 27 PBH . This slightly stronger dependence causes the merger rate to drop faster when 𝑓 PBH is decreased. Since in this regime d 𝑅 ( 0 ) 𝐸 2 ( 𝑡 ) ∝ 𝜎 -21 / 74 M a larger 𝜎 M, that is, more pronounced matter perturbations, will lead to a smaller merger rate.  \nIn practice (43) must be evaluated numerically, which can be quite computationally expensive. A good analytic approximation can be constructed by considering the ¯ 𝑁 ( 𝑦 ) → 0 and ¯ 𝑁 ( 𝑦 ) → ∞ asymptotic of the angular momentum and noting that the following inequalities hold [5]  \n𝑆 E , min ≤ 𝑆 E ≤ 𝑆 E , max < 1 , (44)  \nwhere the ¯ 𝑁 ( 𝑦 ) → 0 and ¯ 𝑁 ( 𝑦 ) → ∞ limits are  \n𝑆 E , max = 5 𝑓 2 PBH 6 𝜎 2 M ! 21 74 𝑈 21 74 , 1 2 , 5 𝑓 2 PBH 6 𝜎 2 M ! , 𝑆 E , min = √ 𝜋 ( 5 / 6 ) 21 / 74 Γ ( 29 / 37 ) " ⟨ 𝑚 2 ⟩/⟨ 𝑚 ⟩ 2 ¯ 𝑁 ( 𝑦 ) + 𝜎 2 M 𝑓 2 PBH # -21 74 𝑒 -¯ 𝑁 ( 𝑦 ) , (45)  \nwhere 𝑈 denotes the confluent hypergeometric function. The factor 𝑆 E can be approximated as  \n𝑆 E ≈ √ 𝜋 ( 5 / 6 ) 21 / 74 Γ ( 29 / 37 ) " ⟨ 𝑚 2 ⟩/⟨ 𝑚 ⟩ 2 ¯ 𝑁 ( 𝑦 ) + 𝐶 + 𝜎 2 M 𝑓 2 PBH # -21 74 𝑒 -¯ 𝑁 ( 𝑦 ) , (46)  \nwhere the factor 𝐶 is fixed by demanding that 𝑆 E → 𝑆 E , max in the limit ¯ 𝑁 ( 𝑦 ) → 0,  \n𝐶 = 𝑓 2 PBH ⟨ 𝑚 2 ⟩/⟨ 𝑚 ⟩ 2 𝜎 2 M       " Γ ( 29 / 37 ) √ 𝜋 𝑈 21 74 , 1 2 , 5 𝑓 2 PBH 6 𝜎 2 M !# -74 21 -1       -1 . (47)  \nThis approximation is accurate within 7% for log-normal mass function with widths 𝜎 ≤ 2. Another approximation, that is accurate within 2% but retains one numerical integral, has been constructed in Ref. [8].', '2.2 Disruption of early binaries in DM haloes': 'After its formation, the binary may get perturbed by encounters with a third PBH. Such encounters are likely in very dense haloes. We can estimate the probability of this by computing the rate of encounters that significantly alter the angular momentum of the binary. However, the small haloes are also unstable and experience gravothermal instabilities within the age of the Universe. This can lead to core collapse and the eventual evaporation of the PBH cluster. In the former case, the density can increase significantly and all binaries in such haloes are likely to get perturbed. In the following, we first estimate the rate of encounters that perturb the binaries and then the probability that the binary is within a halo that goes through gravothermal collapse before the binary merges.  \nConsider an early binary whose coalescence time initially equals the age of the Universe, 𝜏 = 𝑡 0 and a third PBH of mass 𝑚 3 approaching the binary. In a hyperbolic encounter, /one.sup/zero.sup the velocity of the third body of mass 𝑚 3 at the distance of closest approach 𝑟 𝑝 is 𝑣 𝑝 = 𝑏𝑣 rel / 𝑟 𝑝 , where 𝑣 rel denotes the relative velocity of the third body and 𝑏 the impact parameter that is related to 𝑟 𝑝 by 𝑏 2 = 𝑟 2 𝑐 + 2 ( 𝑀 + 𝑚 3 ) 𝑟 𝑐 / 𝑣 2 rel . We can estimate the timescale of the interaction between the binary and the third body as 𝑡 𝑝 ∼ 𝑟 𝑝 / 𝑣 𝑝 . The tidal torque caused by the third body is given by T 𝑝 = 𝑚 3 ( 1 -3ˆ r 𝑝 ⊗ ˆ r 𝑝 )/ 𝑟 3 𝑝 and the resulting change in the binary angular momentum can be estimated as Δ L ≃ 𝑡 𝑝 r 𝑎 × ( T 𝑝 · r 𝑎 ) . This gives  \nΔ 𝑗 ≃ 𝑚 3 𝑟 3 / 2 𝑎 √ 𝑀𝑏𝑣 rel 𝑟 𝑝 . (48)  \nThe binary gets perturbed if Δ 𝑗 is comparable to the initial angular momentum. For the early binaries, we can estimate the initial angular momentum by 𝑗 𝜏 given in Eq. (39). So, encounters with the impact parameter  \n𝑏 ≲ [ 𝑚 3 ( 𝑀 + 𝑚 3 )] 1 / 3 𝑟 1 / 2 𝑎 𝑗 1 / 3 𝜏 𝑀 1 / 6 𝑣 rel ≡ 𝑏 Δ 𝑗>𝑗 𝜏 (49)  \nwould perturb the binary. The corresponding interaction cross section is 𝜎 Δ 𝑗>𝑗 𝜏 = 𝜋𝑏 2 Δ 𝑗>𝑗 𝜏 .  \nTo estimate how large the rate of such processes is, let us consider early PBH haloes of 𝑁 PBHs. Such haloes form when 𝑎 𝑐 ≡ ( 1 + 𝑧 𝑐 ) -1 ≈ 𝑎 eq √ 𝑁 / 𝑓 PBH, where 𝑎 eq denotes the matter-radiation equality. If such haloes are virialized, the velocity dispersion is 𝜎 2 𝑣 ≈ 𝑀 H / 𝑅 where 𝑅 denotes the virial radius and 𝑀 H = ⟨ 𝑚 ⟩ 𝑁 / 𝑓 PBH the mass of the halo. /one.sup/one.sup The average matter density in these haloes is, following the Press-Schechter theory, given by 𝜌 = 3 𝑀 H /( 4 𝜋𝑅 3 ) ≈ 18 𝜋 2 𝜌 𝑐 𝑎 -3 𝑐 , where 𝜌 𝑐 is the critical comoving density. So, the timescale of processes that perturb the binary in  \nsuch haloes can be estimated as  \n𝑡 𝑝 ≃ GLYPH<20> 𝜌 ⟨ 𝑚 ⟩ ⟨ 𝑣 rel 𝜎 Δ 𝑗>𝑗 𝜏 ⟩ GLYPH<21> -1 ≃ 1 Gyr GLYPH<20> 𝑁 1000 GLYPH<21> 19 12 𝑓 -397 222 PBH GLYPH<20> ⟨ 𝑚 ⟩ 𝑀 ⊙ GLYPH<21> -10 111 , (50)  \nwhere used Eq. (7) to translate 𝑟 𝑎 into 𝜏 , took 𝜏 = 𝑡 0, replaced all masses by ⟨ 𝑚 ⟩ and the relative velocity by the velocity dispersion 𝜎 𝑣 . With this rough estimate, we find that for ⟨ 𝑚 ⟩ = O( 10 𝑀 ⊙ ) and 𝑓 PBH = 1 binaries in haloes with 𝑁 ≲ 5800 get perturbed before today. Moreover, by definition, the haloes need at least 𝑁 = 3, implying that binaries are likely not perturbed when 𝑓 PBH ≲ 10 -3 .  \nNext, let us estimate the effect of gravothermal instabilities of small haloes. We estimate this using the characteristic time scale of core collapse, 𝑡 cc ≥ 18 𝑡 r, where 𝑡 r denotes the relaxation time [9]:  \n𝑡 r = 0 . 065 𝜎 3 𝑣 𝑚𝜌 ln Λ ≈ 2kyr 𝑁 7 / 4 𝑓 5 / 2 PBH ln Λ . (51)  \nBy approximating the Coulomb logarithm as ln Λ ≈ ln ( 𝑁 / 𝑓 PBH ) and requiring 18 𝑡 r < 𝑡 0 we get that haloes with  \n𝑁 ≤ 1500 𝑓 10 / 7 PBH ln Λ 4 / 7 ≡ 𝑁 𝑐 (52)  \nexperience gravothermal collapse before today (see also [10]). This indicates that, if 𝑓 PBH = 1, then binaries are perturbed in haloes with 𝑁 < 5300. All binaries survive if 𝑓 PBH ≲ 0 . 005. These numbers are very similar to those found above by estimating the rate of the encounters that would perturb the binaries assuming that the haloes survive.  \nLet us briefly consider the hardness of the early binaries. In the largest haloes in which binaries are likely disrupted according to Eq. (52), the velocity dispersion is  \n𝜎 𝑣 ≈ 0 . 8 km s 𝑓 1 6 PBH GLYPH<20> ⟨ 𝑚 ⟩ 𝑀 ⊙ GLYPH<21> 1 3 𝑁 1 12 ≈ 1 . 8 km s 𝑓 2 7 PBH GLYPH<20> ⟨ 𝑚 ⟩ 𝑀 ⊙ GLYPH<21> 1 3 . (53)  \nOn the other hand, the characteristic binding energy of binaries merging within the age of the Universe is  \n𝐸 bin 𝜇 ≈ 3 . 4 km s 𝜂 -1 8 GLYPH<20> 𝑗 0 . 01 GLYPH<21> 7 8 GLYPH<20> 𝑀 𝑀 ⊙ GLYPH<21> 1 8 ! 2 . (54)  \nTherefore, the binaries can be considered hard as encounters with other PBHs are unlikely to ionise them within environments where encounters are frequent enough.  \nAssuming that all of the binaries in haloes of size 𝑁 ≤ 𝑁 𝑐 ( 𝑧 ) get perturbed, we can estimate the suppression of the merger rate by estimating the probability that the binary is in such halo. Including both isolated haloes and subhaloes, the suppression factor, shown in Fig. 6, is given by [11]  \nFig. 6 The suppression factor of merger rate for the early two-body channel as a function of PBH abundance at different redshifts. The dashed curves show the approximation (60) and the dotted curves show the suppression factor neglecting the subhalos.  \n<!-- image -->  \n𝑆 L ( 𝑧 ) = 1 -𝑁 𝑐 ( 𝑧 ) ∑︁ 𝑁 = 3      𝑝 𝑁 ( 𝑧 𝑐 ) + ∑︁ 𝑁 \' >𝑁 𝑐 ( 𝑧 ) 𝑝 𝑁,𝑁 \' ( 𝑧 𝑐 )      , (55)  \nwhere 𝑧 𝑐 = 𝑓 PBH 𝑧 eq / √ 𝑁 𝑐 denotes the redshift at which the haloes of 𝑁 𝑐 PBHs form and 𝑧 eq denotes the redshift of matter-radiation equality. The probability of the binary being part of a halo containing 𝑁 PBHs at redshift 𝑧 is given by  \n𝑝 𝑁 ( 𝑧 ) = 𝑁 -1 / 2 𝑒 -𝑁 / 𝑁 ∗ ( 𝑧 ) " 𝑁> 2 𝑁 -1 / 2 𝑒 -𝑁 / 𝑁 ∗ ( 𝑧 ) , (56)  \nand the probability of the binary being part of a subhalo containing 𝑁 PBHs inside a halo of 𝑁 \' by  \n𝑝 𝑁,𝑁 \' ( 𝑧 ) = 𝑝 𝑁 \' ( 𝑧 ) 𝑁 -1 / 2 𝑒 -𝑁 / 𝑁 ∗ ( 𝑧 ) " 𝑁 \' 𝑁 = 2 𝑁 -1 / 2 𝑒 -𝑁 / 𝑁 ∗ ( 𝑧 ) . (57)  \nHere 𝑁 ∗ ( 𝑧 ) is the characteristic number of PBHs in a halo at redshift 𝑧 , and can be estimated by [7]  \n𝑁 ∗ ( 𝑧 ) = [ ln ( 1 + 𝛿 ∗ ) -𝛿 ∗ /( 1 + 𝛿 ∗ )] -1 (58)  \nwith  \n𝛿 ∗ ≈ 1 . 69 h 𝑓 PBH 1 𝐹 2 GLYPH<16> ( 1 -√ 21 )/ 4 , ( 1 + √ 21 )/ 4; 1; -( 𝑧 eq + 1 )/( 𝑧 + 1 ) GLYPH<17> i -1 . (59)  \nThe numerically evaluated suppression factor 𝑆 L ( 𝑧 ) can be approximated by [11]  \n𝑆 L ( 𝑡 ) ≈ min       1 , 0 . 01 " GLYPH<18> 𝑡 𝑡 0 GLYPH<19> 0 . 44 𝑓 PBH # -0 . 65 𝑒 0 . 03 ln 2 GLYPH<20> GLYPH<16> 𝑡 𝑡 0 GLYPH<17> 0 . 44 𝑓 PBH GLYPH<21>       , (60)  \nwhich, as seen in Fig. 6, fits the numerical estimate well. The plot also shows that for large 𝑓 PBH (e.g. 𝑓 PBH > 0 . 03 at 𝑧 = 0) the dominant contribution to the suppression arises from binary disruption in subhalos. We remark that this suppression factor might be affected by phenomena, such as tidal stripping, which can modify subhalo evolution. So far, DM subhaloes in PBH cosmologies have not been studied in detail and thus this term might be subject to unknown systematic uncertainties.', '2.3 The early three-body channel: merger rate of perturbed binaries': "Asdiscussed above, a sizeable fraction of PBHs will form binaries after they decouple from the Hubble flow. However, many of them can be perturbed by interactions with other PBHs. The population of such binaries will not share the characteristics of the unperturbed binaries evolving in isolation. This was sharply illustrated by the simulated angular momentum distribution shown in Fig. 4, which revealed a peaked distribution for unperturbed binaries, but a nearly flat one (d 𝑃 / d 𝑗 ≈ 1) for perturbed ones. The perturbed binaries tend to be hard so the interactions with other PBHs are likely not break them apart. Since they are expected to be still around, many of them could contribute to the present merger rate.  \nMost perturbed binaries in the early Universe that would merge today would have formed from PBH pairs that initially evolved into an eccentric binary with a very short coalescence time. Let us first make a few order-of-magnitude estimates to understand the formation mechanism. Consider a compact PBH pair that decouples from the Hubble flow. Its expected initial angular momentum due to tidal forces is 𝑗 = O( 10 -2 ) . After being perturbed we expect that 𝑗 = O( 1 ) , that is, the angular momentum grows by about two orders of magnitude. Assuming that this process formed a binary that merges today and given that 𝜏 ∝ 𝑗 7 , the initial coalescence time would have been 14 orders of magnitude shorter than the age of the Universe, that is, about an hour. Such binaries would not have had time to interact with other PBHs before the merger.  \nIt follows that perturbed binaries contributing to the present PBH merger rate are formed from PBH pairs residing in a tidal field that induces angular momenta 𝑗 ≫ O( 10 -2 ) , which tend to be larger than for unperturbed binaries discussed in the last section. Such a tidal field implies the presence of a third PBH at a distance closer than the average separation. The presence of a nearby 3rd PBH means that we are looking at a 3-body system that eventually evolves into a binary after the third body is ejected. The binding energy of the binary is comparable to the binding energy of the 3-body system and its angular momentum can be estimated to follow the distribution  \nd 𝑃 ( 𝑗 ) d 𝑗 = 𝛾 𝑗 𝛾 -1 , (61)  \nwith 𝛾 ∈ [ 1 , 2 ] . The limiting case 𝛾 = 1 is seen in numerical simulations of early Universe binary formation [5] and is shown in Fig. 4 by the yellow histogram. The other limit, 𝛾 = 2, corresponds to the thermal distribution [12]. However, we note  \nthat gravitating many-body systems are unstable and will not thermalize completely, so it is not clear if a thermal 𝛾 = 2 distribution can be realized. Also, assuming that the angular momentum distribution 𝑃 ( 𝑗 ) is independent of the initial angular momentum agrees with numerical studies of binary-single body collisions [13, 14].  \nWith these considerations in mind, we can estimate the merger rate of perturbed binaries. Consider three PBHs with masses 𝑚 𝑖 , 𝑖 ∈ { 1 , 2 , 3 } that are of the similar. We assume that the Hubble decoupling of the 3-body system is hierarchical /one.sup/two.sup so that the PBHs '1' and '2' at a comoving separation 𝑥 will decouple first at  \n𝑎 ( 𝑡 12 ) ≈ 𝑎 eq 𝑁 ( 𝑥 ) 𝑓 PBH ⟨ 𝑚 ⟩ 𝑚 1 + 𝑚 2 , (62)  \nand the third PBH at a comoving distance 𝑦 from the centre of mass of the pair will decouple at  \n𝑎 ( 𝑡 𝑝 ) ≈ 𝑎 eq 𝑁 ( 𝑦 ) 𝑓 PBH ⟨ 𝑚 ⟩ 𝑚 1 + 𝑚 2 + 𝑚 3 . (63)  \nDuring radiation domination  \n𝑡 p = 65 . 3 kyr × GLYPH<18> ⟨ 𝑚 ⟩ 𝑚 1 + 𝑚 2 + 𝑚 3 GLYPH<19> 2 𝑁 ( 𝑦 ) 2 𝑓 2 PBH . (64)  \nAs discussed above, the coalescence time 𝜏 i of the binary '1+2' can be short. In the following, we will consider systems with  \n𝜏 i > 20 𝑡 p . (65)  \nThe binding energy of eccentric binaries scales as 𝐸 ( 𝑡 ) = 𝐸 0 ( 1 -𝑡 / 𝜏 ) -2 [2], so this condition is equivalent to demanding that the binding energy is affected by less than 11% by GW emission. We also note that, during matter-domination, Eq. (64) overestimates 𝑡 p, so (65) would lead to a lower merger rate estimate.  \nWith hierarchical collapse, the coalescence time of the initial '1+2' binary depends on its angular momentum 𝑗 12 which is set mostly by the tidal torque from the third PBH. By Eqs. (27) and (26), it is given by  \n𝑗 12 = 1 . 4 𝑚 3 𝑚 1 + 𝑚 2 ¯ 𝑁 ( 𝑥 ) ¯ 𝑁 ( 𝑦 ) | sin ( 2 𝜃 )| , (66)  \nwhere 𝜃 denotes the angle between the line joining the '1+2' pair and the line joining the center of mass of the pair and the third PBH. The order of magnitude of the angular momentum generated by the remaining PBHs is of the order of the width of the distribution (37), that is, 𝜎 𝑗 ≈ 0 . 5 ¯ 𝑁 ( 𝑥 )/ √︁ ¯ 𝑁 ( 𝑦 ) . If ¯ 𝑁 ( 𝑦 ) ≪ 1, then 𝑗 12 ≫ 𝜎 𝑗 and we can neglect the contribution from the remaining PBHs.  \nIn a hierarchical setup, the binding energy of the 3-body system is determined mostly by the binding energy of the initial binary. The comoving density of initial '1+2' binaries with per binding energy interval ( 𝐸, 𝐸 + d 𝐸 ) is  \nd 𝑛 𝑖 ( 𝐸, 𝑚 𝑖 ) d 𝐸 = ∫ d 𝑛 3 𝑃 ( 𝐸, 𝜏 𝑖 ≥ 20 𝑡 𝑝 | 𝑥, 𝑦, 𝜃, 𝑚 𝑖 ) , (67)  \nwhere  \nd 𝑛 3 = 𝑛 2 d 𝑛 ( 𝑚 1 ) 𝑛 d 𝑛 ( 𝑚 2 ) 𝑛 d 𝑛 ( 𝑚 3 ) 𝑛 d ¯ 𝑁 ( 𝑥 ) d ¯ 𝑁 ( 𝑦 ) 𝑒 -¯ 𝑁 ( 𝑦 ) d cos 𝜃 2 (68)  \ngives the density of 3-body configurations yielding the initial binary. It generalizes Eq. (9). The probability that the '1+2' pair produces a binary with binding energy 𝐸 that satisfies the bound (65) is  \n𝑃 ( 𝐸, 𝜏 𝑖 ≥ 20 𝑡 𝑝 | 𝑥, 𝑦, 𝜃, 𝑚 𝑖 ) = 𝛿 ( 𝐸 -𝐸 ( 𝑥 )) Θ ( 𝜏 i ( 𝐸, 𝑗 12 ) -20 𝑡 p ) , (69)  \nwhere Θ denotes the step function. The binding energy 𝐸 = 𝑚 1 𝑚 2 / 2 𝑟 𝑎 (3) is determined from the initial separation 𝑥 by Eq. (23). It gives 𝑥 4 ≈ 0 . 6 ( 𝑚 1 + 𝑚 2 ) 𝑚 1 𝑚 2 /( 𝜌 𝑅 𝐸 ) .  \nConsider now the integrals over 𝑥 , 𝑦 and 𝜃 in (67). The 𝑥 integral can be eliminated using the 𝛿 -function. The 𝑦 integral can be replaced with the integral over ¯ 𝑁 ( 𝑦 ) and the condition (65) can now be recast as an upper bound on ¯ 𝑁 ( 𝑦 ) ,  \n𝑁  \n¯  \n(  \n𝑦  \n) ≤  \n¯  \n𝑁  \n(  \n𝑦  \n)  \nmax  \n≡ 3800 𝑓 PBH 𝐸 -37 36 𝜌 1 4 𝑅 𝑚 11 12 1 𝑚 11 12 2 𝑚 7 9 3 ( 𝑚 1 + 𝑚 2 + 𝑚 3 ) 2 9 ( 𝑚 1 + 𝑚 2 ) 11 36 ⟨ 𝑚 ⟩ | sin ( 2 𝜃 )| 7 9 . (70)  \nFor consistency, we will further require that ¯ 𝑁 ( 𝑦 ) max < 1 at every orientation, which yields a lower bound on the binding energy  \n𝐸 > 𝐸 min ≡ ⟨ 𝑚 ⟩ GLYPH<20> 11 km s GLYPH<21> 2 GLYPH<20> ⟨ 𝑚 ⟩ 𝑀 ⊙ GLYPH<21> 18 37 𝑓 36 37 PBH G( 𝑚 1 , 𝑚 2 , 𝑚 3 ) (71)  \nwhere we defined the dimensionless factor  \nG( 𝑚 1 , 𝑚 2 , 𝑚 3 ) = 2 𝑚 3 1 𝑚 3 2 𝑚 1 + 𝑚 2 ! 11 37 𝑚 28 37 3 GLYPH<16> 𝑚 1 + 𝑚 2 + 𝑚 3 3 GLYPH<17> 8 37 ⟨ 𝑚 ⟩ -91 37 (72)  \nthat becomes unity for monochromatic mass functions. With this cut-off, we can simplify the computation and omit the exponential suppression in (68), so the integral over ¯ 𝑁 ( 𝑦 ) can be approximated by ¯ 𝑁 ( 𝑦 ) max.  \nAveraging over 𝜃 then gives the initial binding energy distribution  \nwhere  \nd 𝑛 𝑖 ( 𝐸, 𝑚 𝑖 ) d 𝐸 d 𝑚 1d 𝑚 2d 𝑚 3 = d 𝑛 𝑖, mono ( 𝐸 ) d 𝐸 F( 𝑚 1 , 𝑚 2 , 𝑚 3 ) 3 GLYPH<214> 𝑖 = 1 d 𝑛 𝑛 d 𝑚 𝑖 , (73)  \nd 𝑛 𝑖, mono ( 𝐸 ) d 𝐸 = 4 . 0 × 10 10 𝐸 -25 9 𝑓 3 PBH ⟨ 𝑚 ⟩ 16 9 𝜌 3 2 𝑅 (74)  \nand we defined the dimensionless factor  \nF( 𝑚 1 , 𝑚 2 , 𝑚 3 ) = 𝑚 5 3 1 𝑚 5 3 2 𝑚 7 9 3 GLYPH<16> 𝑚 1 + 𝑚 2 2 GLYPH<17> 4 9 GLYPH<16> 𝑚 1 + 𝑚 2 + 𝑚 3 3 GLYPH<17> 2 9 ⟨ 𝑚 ⟩ -43 9 (75)  \nthat reduces to unity for monochromatic mass functions.  \nTo convert Eq. (73) into the energy distribution of binary parameters after the 3rd body is ejected, we assume that a 3-body system with binding energy 𝐸 ' yields a binary with binding energy 𝐸  \n𝜕𝐾 ( 𝐸 | 𝐸 ' ) 𝜕𝐸 = 𝛼 𝐸 ' 𝑒 -𝛼 ( 𝐸 / 𝐸 ' -1 ) Θ ( 𝐸 -𝐸 ' ) . (76)  \nThis approximates the behaviour observed in numerical simulations of PBH formation in the early universe for which 𝛼 = 1 gives a decent fit [5, 11]. The limiting case in which the initial binding energy does not change corresponds to the limit 𝛼 →∞ . The Heggie-Hills law states that hard binaries tend to get harder, while soft binaries get softer [15, 16]. Thus, as the early binaries tend to be hard, decreasing 𝛼 will support more sizeable energy changes and can thus model the effect of subsequent collisions with other PBHs in PBH clusters. Since we work with the assumption of narrow mass functions, we do not consider the dependence on the masses of the three interacting PBHs in Eq. (76).  \nTypically the lightest PBH will be ejected from the 3-body system, as it is more likely to be accelerated beyond the escape velocity. So, the binary will be typically composedofthetwoheavier PBHs. To model this, we will assume that the probability of ejecting a PBH from the system is exponentially suppressed by the relative mass 𝑝 𝑖 ( 𝑚 𝑖 ) ∝ exp (-𝜆𝑚 𝑖 /( 𝑚 1 + 𝑚 2 + 𝑚 3 )) , where 𝜆 ≥ 0 is a free parameter. 𝜆 = 0 corresponds to equal probability of ejection, while for 𝜆 → ∞ the lightest PBH is always expelled. The exact shape of this function must be determined numerically. Given 𝑚 1, 𝑚 2, 𝑚 3, the binary constituent mass distribution is  \n𝑃 𝑚 ( 𝑚 ' 1 , 𝑚 ' 2 | 𝑚 1 , 𝑚 2 , 𝑚 3 ) = 𝑝 1 𝛿 ( 𝑚 ' 1 -𝑚 2 ) 𝛿 ( 𝑚 ' 2 -𝑚 3 ) + 𝑝 2 𝛿 ( 𝑚 ' 1 -𝑚 1 ) 𝛿 ( 𝑚 ' 2 -𝑚 3 ) + 𝑝 3 𝛿 ( 𝑚 ' 1 -𝑚 1 ) 𝛿 ( 𝑚 ' 2 -𝑚 2 ) . (77)  \nThe distribution of orbital parameters of the early binaries from the 3-body channel is  \nd 𝑛 E3 ( 𝑗 , 𝐸, 𝑚 ' 1 , 𝑚 ' 2 ) d 𝑗 d 𝐸 d 𝑚 ' 1 d 𝑚 ' 2 = d 𝑃 ( 𝑗 ) d 𝑗 ∫ d 𝐸 ' d 𝑚 1d 𝑚 2d 𝑚 3 𝜕𝐾 ( 𝐸 | 𝐸 ' ) 𝜕𝐸 × d 𝑛 𝑖 ( 𝐸 ' ) 𝐸 ' 𝑚 𝑚 𝑚 𝑃 𝑚 ( 𝑚 ' 1 , 𝑚 ' 2 | 𝑚 1 , 𝑚 2 , 𝑚 3 ) , (78)  \nd d 1d 2d 3  \nwith the angular momentum distribution given by Eq. (61). Since 𝜕𝑛 𝑖 ( 𝐸 )/ 𝜕𝐸 ∝ 𝐸 -25 / 9 , it is possible to evaluate the 𝐸 ' integral  \nd 𝑛 E3 ( 𝑗 , 𝐸, 𝑚 1 , 𝑚 2 ) d 𝑗 d 𝐸 = d 𝑃 ( 𝑗 ) d 𝑗 d 𝑛 𝑖, mono ( 𝐸 ) d 𝐸 K( 𝛼 ) ¯ F( 𝑚 ' 𝑖 ) d 𝑛 ( 𝑚 1 ) 𝑛 d 𝑛 ( 𝑚 2 ) 𝑛 , (79)  \nwhere  \n¯ F( 𝑚 1 , 𝑚 2 ) ≡ ∫ d 𝑛 ( 𝑚 ) 𝑛 GLYPH<20> ( 𝑝 1 ( 𝑚, 𝑚 1 , 𝑚 2 ) + 𝑝 1 ( 𝑚 1 , 𝑚, 𝑚 2 ))F ( 𝑚, 𝑚 1 , 𝑚 2 ) + 𝑝 3 ( 𝑚 1 , 𝑚 2 , 𝑚 )F ( 𝑚 1 , 𝑚 2 , 𝑚 ) GLYPH<21> (80)  \naccounts for the variation in masses and the overall factor  \nK( 𝛼 ) ≡ 𝛼 -25 / 9 𝑒 𝛼 Γ ( 25 / 9 , 𝛼 ) (81)  \ncontains the effect of hardening in binary-single PBH collisions. Γ ( 𝑎, 𝑏 ) denotes the incomplete gamma function. We will use 𝛼 = 1, which corresponds to K( 1 ) = 4 . 0. As was discussed above, accounting for hardening via later encounter would reduce 𝛼 and increase the merger rate. However, as the effect of later encounters was estimated to be small [17], we will not consider this contribution. The enhancing effect of hardening might, however, become more relevant in scenarios in which PBHs are initially clustered because the probability of binary-PBH collisions is expected to increase in such cases.  \nWhen integrating over masses in Eq. (80) we neglected the lower bound 𝐸 min (71) on the binding energy 𝐸 . To estimate 𝐸 min for the binary after the 3rd PBH is ejected, we make the substitution in Eq. (71),  \nG( 𝑚 1 , 𝑚 2 , 𝑚 3 ) → ¯ G( 𝑀 ) = GLYPH<18> 𝑀 2 ⟨ 𝑚 ⟩ GLYPH<19> 91 / 37 , (82)  \nwhere 𝑀 is the total mass of the binary. As the binary is more likely to consist of the heaviest PBHs, ¯ G will give a stricter lower bound than what would follow from G , i.e., ¯ G ≳ G . This will in turn soften the resulting merger rate. However, we expect the effect to be small. As with G , we have that ¯ G = 1 for monochromatic mass functions.  \nThe merger rate implied by the distribution of binaries (79) is given by Eq. (6). As with 𝑅 E2 in Eq. (12), we can assume that the binaries formed at the Big Bang, so that  \nd 𝑅 E3 d ln 𝑚 1d ln 𝑚 2 = ∫ d 𝑗 d 𝐸 d 𝑛 E3 ( 𝑗 , 𝐸 ) d 𝑗 d 𝐸 d ln 𝑚 1d ln 𝑚 2 𝛿 ( 𝑡 -𝜏 ( 𝑗 , 𝐸 )) = 𝑅 E3 , mono GLYPH<18> 𝑀 2 ⟨ 𝑚 ⟩ GLYPH<19> 179 𝛾 259 -2122 333 ( 4 𝜂 ) -3 𝛾 7 -1 ¯ F( 𝑚 1 , 𝑚 2 ) 𝜓 ( 𝑚 1 ) 𝜓 ( 𝑚 2 ) , (83)  \nwhere the prefactor corresponds to the merger rate for a monochromatic mass function  \n𝑅 E3 , mono = 7 . 9 × 10 4 Gpc 3 yr 𝑓 144 𝛾 259 + 47 37 PBH GLYPH<20> 𝑡 𝑡 0 GLYPH<21> 𝛾 7 -1 GLYPH<20> ⟨ 𝑚 ⟩ 𝑀 ⊙ GLYPH<21> 5 𝛾 -32 37 𝑒 3 . 2 ( 1 -𝛾 ) 𝛾 28 / 9 -𝛾 K . (84)  \nThis merger rate depends weakly on the exact form of the condition (65). For instance, imposing 𝜏 𝑖 > 𝑞𝑡 𝑝 ( 𝑦 ) , we get 𝑅 E3 ∝ 𝑞 16 ( 28 -9 𝛾 ) 2331 .  \nWestress that the three-body merger rate estimate (83) holds for relatively narrow mass functions . A hierarchical decoupling, i.e., 𝑎 ( 𝑡 12 ) ≪ 𝑎 ( 𝑡 p ) implies a relatively narrow mass function which guarantees that 𝑚 1 + 𝑚 2 is of the same order as 𝑚 3. If 𝑚 3 ≫ 𝑚 1 + 𝑚 2, then the gravitational potential of the third body could dominate the '1+2' system before its decoupling from the Hubble flow. On the other hand, if 𝑚 3 ≪ 𝑚 1 + 𝑚 2, the third body may not be heavy enough to disrupt the initial binary. So, extending these estimates to wider mass functions requires considering the cases involving large mass ratios separately.", '3 PBH binary formation in the late Universe': 'In the late Universe, PBH binaries can form dynamically through two- or three-body interactions in dense structures that have decoupled from the Hubble flow. Binary formation in two-body encounters happens through the emission of GWs while in a three-body system, one of the bodies carries away kinetic energy from the system leaving the other two in a bound state.  \nWecharacterize the PBH number overdensities with an effective parameter 𝛿 eff > 1. Consider first a structure that includes 𝑁 PBHs. In such a structure, the average PBH number overdensity is 𝛿 𝑁 = 𝑁 /( 𝑉 𝑁 𝑛 PBH ) > 1 where 𝑛 PBH = " 𝑁 𝑛 𝑁 𝑁 is the comoving average PBH number density and 𝑉 𝑁 is the volume containing the 𝑁 PBHs. The PBH merger rate is 𝑅 𝑁 = 𝛿 𝑞 𝑁 ( 𝜎 𝑣,𝑁 / 𝜎 𝑣 ) 𝑏 𝑅 1 where 𝑞 = 2 for the two-body and 𝑞 = 3 for the three-body channel, 𝜎 𝑣,𝑁 is the velocity dispersion of PBHs in a halo with 𝑁 PBHs, 𝜎 𝑣 is a reference value for the velocity dispersion, 𝑏 = -11 / 7 for the two-body and 𝑏 ∈ (-55 , -47 )/ 7 for the three-body channel (as shown in the following) and 𝑅 1 is the merger rate assuming the background value for the PBH number density and the velocity dispersion 𝜎 𝑣 . The total PBH merger rate is then obtained by summing over different size halos:  \n𝑅 = ∑︁ 𝑁 𝑛 𝑁 𝑉 𝑁 𝑅 𝑁 = " 𝑁 𝑛 𝑁 𝑁𝛿 𝑞 -1 𝑁 ( 𝜎 𝑣,𝑁 / 𝜎 𝑣 ) 𝑏 𝑛 PBH 𝑅 1 ≡ 𝛿 𝑞 -1 eff 𝑅 1 , (85)  \nwhere in the last step we defined the effective PBH overdensity 𝛿 eff that includes also the effect of the differences in the velocity dispersions in different halos. We assume that the mass distribution of PBHs in DM haloes matches their overall mass distribution.', '3.1 Two-body encounters': 'Consider a hyperbolic encounter of two bodies of masses 𝑚 1 and 𝑚 2 with initial relative velocity 𝑣 rel and impact parameter 𝑏 , as shown in Fig. 7. In the centre of mass frame, the energy and angular momentum of the system can be written in terms of the polar coordinates of the separation r 1 -r 2 ≡ ( 𝑟 cos 𝜙, 𝑟 sin 𝜙, 0 ) as  \n𝐸 = 1 2 𝜇 GLYPH<16> / 𝑟 2 + 𝑟 2 / 𝜙 2 GLYPH<17> -𝜇𝑀 𝑟 , L = 𝜇𝑟 2 / 𝜙 ˆ 𝑧 , (86)  \nwhere 𝑀 = 𝑚 1 + 𝑚 2 and 𝜇 = 𝑚 1 𝑚 2 / 𝑀 . Due to conservation of energy and angular momentum, the distance 𝑟 between the bodies evolves as  \n𝑟 = 𝑏 sin 𝜙 0 cos ( 𝜙 -𝜙 0 ) -cos 𝜙 0 , / 𝜙 = -𝑏𝑣 rel 𝑟 2 , (87)  \nwhere 𝜙 0 = 𝜋 -tan -1 ( 𝑏𝑣 rel / 𝐺𝑀 ) corresponds to the distance of closest approach.  \nFig. 7 A hyperbolic encounter of two PBHs.  \n<!-- image -->  \nThe system loses energy by emitting GWs. The energy emitted by unit time is given by the third time derivative of the quadrupole moment tensor 𝑄 𝑖 𝑗 = 𝑚 1 𝑟 1 ,𝑖 𝑟 1 , 𝑗 + 𝑚 2 𝑟 2 ,𝑖 𝑟 2 , 𝑗 as / 𝐸 GW = ¤ 𝑄 𝑖 𝑗 ¤ 𝑄 𝑖 𝑗 / 45. By integrating over time and taking the limit of small 𝑣 rel we get  \n𝐸 GW = 19 𝜋𝜂 2 𝑀 8 6 𝑏 7 𝑣 7 rel , (88)  \nwhere 𝜂 = 𝜇 / 𝑀 .  \nThe bodies form a binary if the energy emitted in GWs is larger than the initial kinetic energy of the pair, 𝐸 0 = 𝜇𝑣 2 rel / 2. This corresponds to the maximal impact parameter  \n𝑏 max = GLYPH<18> 19 𝜋𝜂 3 GLYPH<19> 1 / 7 𝑀𝑣 -9 / 7 rel . (89)  \nThe semimajor axis and eccentricity of the resulting binary are  \n𝑟 𝑎, 0 = -𝜂𝑀 2 2 𝐸 𝑓 , 𝑒 0 = √︄ 1 + 2 𝐸 𝑓 𝑣 2 rel 𝑏 2 𝜂𝑀 3 , (90)  \nwhere 𝐸 𝑓 ≡ 𝐸 0 -𝐸 GW, and by Eq. (8) the binary coalescence time is  \n𝜏 = 19 𝜋𝑀 85 𝑣 3 rel √︄ ( 𝑏 / 𝑏 max ) 21 1 - ( 𝑏 / 𝑏 max ) 7 ≈ 1 day 𝑀 𝑀 ⊙ GLYPH<18> 𝑣 rel 100 km / s GLYPH<19> -3 √︄ ( 𝑏 / 𝑏 max ) 21 1 - ( 𝑏 / 𝑏 max ) 7 . (91)  \nThis expression diverges for 𝑏 = 𝑏 max since in that case the orbit is a parabola, but already for 𝑏 < 0 . 99 𝑏 max the square-root factor is smaller than 4 and the binaries merge very quickly.  \nThe capture cross section is then given by 𝜎 = 𝜋𝑏 2 max . Assuming a MaxwellBoltzmann velocity distribution 𝑃 ( 𝑣 ) ∝ 𝑣 2 exp GLYPH<0> -𝑣 2 / 𝜎 2 𝑣 GLYPH<1> , that matches reasonably  \nwell the results of N-body simulations [18], we get [19, 20, 21]  \n⟨ 𝑣 rel 𝜎 2 ⟩ = ∫ d 𝑣 rel 𝑣 rel 𝑃 ( 𝑣 rel ) 𝜋𝑏 2 max ≈ 2 √ 𝜋 GLYPH<18> 19 𝜋 𝜂 3 GLYPH<19> 2 / 7 𝑀 2 𝜎 11 / 7 𝑣 . (92)  \nThe velocity dispersion 𝜎 𝑣 typically takes a value in the 1 -100 km / s range. Since the coalescence times of the binaries formed through the two-body capture process are very short, we can directly translate the capture rate into the PBH merger rate:  \nd 𝑅 L2 d ln 𝑚 1d ln 𝑚 2 = ⟨ 𝑣 rel 𝜎 mer ⟩ 𝛿 eff d 𝑛 PBH d ln 𝑚 1 d 𝑛 PBH d ln 𝑚 2 = 2 √ 𝜋 GLYPH<18> 19 𝜋 𝜂 3 GLYPH<19> 2 / 7 𝑀 2 𝜎 11 / 7 𝑣 𝛿 𝑓 2 PBH 𝜌 2 DM 𝑚 1 𝑚 2 𝜓 ( 𝑚 1 ) 𝜓 ( 𝑚 2 ) ≈ 3 . 4 × 10 -6 Gpc 3 yr 𝑓 2 PBH 𝛿 eff GLYPH<18> 𝜎 𝑣 km / s GLYPH<19> -11 7 𝜂 -5 7 𝜓 ( 𝑚 1 ) 𝜓 ( 𝑚 2 ) . (93)  \nClustered PBHs tend to have also a smaller velocity dispersion which further enhances the rate (93). Moreover, clustering tends to expel light PBHs from the regions of large 𝛿 eff which effectively causes spatial dependence of the PBH mass function 𝜓 ( 𝑚 ) .', '3.2 Three-body encounters': 'The merger rate of binaries formed from in three-body encounters is estimated by multiplying the rate of two-body encounters by the probability that a third body happens to be close to the pair and estimating the coalescence times of these binaries [22, 23, 24, 25, 26, 1, 27, 17].  \nUsing (87) at 𝜙 = 𝜙 0, we find that in a hyperbolic encounter, the distance between the two bodies gets smaller than 𝑟 𝑎 if the impact parameter is 𝑏 2 < 𝑟 2 𝑎 ( 1 + 2 𝑀 /( 𝑟 𝑎 𝑣 2 rel )) ≡ 𝑏 2 𝑎 . Assuming a Maxwell-Boltzmann velocity distribution, the velocity averaged cross section for the PBHs to reach within distance 𝑟 𝑎 is  \n⟨ 𝑣 rel 𝜎 2 ⟩ = ∫ d 𝑣 rel 𝑣 rel 𝑃 ( 𝑣 rel ) 𝜋𝑏 2 𝑎 ≈ 2 √ 𝜋𝑟 𝑎 𝜎 𝑣 GLYPH<18> 1 + 2 𝑀 𝑟 𝑎 𝜎 2 𝑣 GLYPH<19> . (94)  \nThe time the PBHs spend at a distance smaller than 𝑟 𝑎 can be integrated from the solution (87). In the limit of small 𝑣 rel this gives 𝜏 𝑎 ≈ √︁ 8 𝑟 3 𝑎 /( 9 𝑀 ) . Then, the probability that a third PBH passes the PBH pair within distance 𝑟 𝑎 from its center of mass is  \n𝑃 3 = ∫ d 𝑣 3 𝑃 ( 𝑣 3 ) 𝜋𝑟 2 𝑎 ( 𝑣 3 𝜏 𝑎 + 𝑟 𝑎 ) 𝛿𝑛 PBH =      4 3 √︄ 2 𝜋𝑟 7 𝑎 𝑀 𝜎 𝑣 + 𝜋𝑟 3 𝑎      𝛿 eff 𝑛 PBH . (95)  \nThe rate of three-body encounters can be approximated as ∝ 𝑃 3 ⟨ 𝑣 rel 𝜎 2 ⟩ . The encounter leads to the formation of a binary if the semimajor axis 𝑟 𝑎 is sufficiently small, 𝑟 𝑎 < 𝑟 𝑎, max. The upper limit on 𝑟 𝑎 is estimated by setting a lower limit on the binary hardness factor 𝜅 defined by comparing the binary binding energy to the average kinetic energy of the surrounding PBHs:  \n𝜅 ≡ 𝐸 bin ⟨ 𝑚 ⟩ 𝜎 2 𝑣 / 2 = 𝜇𝑀 𝑟 𝑎 ⟨ 𝑚 ⟩ 𝜎 2 𝑣 , (96)  \nwhere 𝐸 bin is the binding energy given by Eq. (3). The threshold value of 𝜅 > 𝜅 min = 5 is often used in the literature (see e.g. [25, 27, 17]).  \nUnlike the two-body encounters, the three-body encounters produce wide binaries whose coalescence time can be even longer than the age of the Universe. To estimate the merger rate, we need to first estimate the angular momenta of the binaries. Analogously to the early universe 3-body scenarios (61), we take  \nd 𝑃 d 𝑗 = 𝛾 𝑗 𝛾 -1 , (97)  \nwhere 𝛾 ∈ [ 1 , 2 ] . The value 𝛾 = 2 corresponds to the thermal distribution [12] but Refs. [28, 5] find that the distribution could be superthermal with 𝛾 < 2.  \nWith the above ingredients, the binary formation rate in the three-body encounters is  \nd 𝑅 𝑏 d 𝑚 1d 𝑚 2d 𝑗 d 𝑟 𝑎 = d 𝑃 d 𝑗 d 𝑃 3 ⟨ 𝑣 rel 𝜎 2 ⟩ d 𝑟 𝑎 Θ ( 𝑟 𝑎 -𝑟 𝑎, max ) 𝛿 eff d 𝑛 PBH d 𝑚 1 d 𝑛 PBH d 𝑚 2 . (98)  \nThen, using Eq. (6), we get the merger rate of the binaries formed in the three-body encounters:  \nd 𝑅 L3 d ln 𝑚 1d ln 𝑚 2 = 7 𝜋 24 𝑓 3 PBH 𝛿 2 eff GLYPH<18> 1360 𝜂𝜎 8 𝑣 𝑡 3 𝑀 GLYPH<19> 𝛾 7 𝑀 3 𝜌 3 DM 𝜂 ⟨ 𝑚 ⟩ 𝜎 9 𝑣 F GLYPH<18> ⟨ 𝑚 ⟩ 2 𝜂𝑀 𝜅 min GLYPH<19> 𝜓 ( 𝑚 1 ) 𝜓 ( 𝑚 2 ) ≈ 1 . 3 × 10 -16 𝑒 -6 . 0 ( 𝛾 -1 ) Gpc 3 yr 𝑓 3 PBH 𝛿 2 eff GLYPH<18> 𝜎 𝑣 km / s GLYPH<19> -9 + 8 𝛾 7 𝜂 -1 + 𝛾 7 GLYPH<18> 𝑀 𝑀 ⊙ GLYPH<19> 3 -𝛾 7 GLYPH<18> 𝑡 𝑡 0 GLYPH<19> 𝛾 7 × F GLYPH<18> ⟨ 𝑚 ⟩ 2 𝜂𝑀 𝜅 min GLYPH<19> 𝜓 ( 𝑚 1 ) 𝜓 ( 𝑚 2 ) , (99)  \nwhere  \nFormation of primordial black hole binaries and their merger rates  \nF( 𝜅 ) ≡ 𝜅 -4 + 4 𝛾 7 GLYPH<20> 6 √ 𝜋 7 -𝛾 + 72 63 -8 𝛾 𝜅 -1 2 + 15 √ 𝜋 70 -8 𝛾 𝜅 -1 + 22 77 -8 𝛾 𝜅 -3 2 GLYPH<21> . (100)  \nSimilarly, as in the case of two-body encounters, the rate is larger in the higher density environments with smaller velocity dispersion. However, the three-body channel is even more sensitive to the clustering. Moreover, the coalescence time of the binaries formed through the three-body encounters is typically very large, so the binaries that would merge in the present Universe should have formed at a high redshift in compact PBH clusters induced by their Poisson distribution.', '4 Summary': 'In conclusion, the merger rate of PBH binaries consists of four components which we report here for standard PBH cosmologies:', '1. Early two-body channel': 'd 𝑅 𝐸 2 d ln 𝑚 1d ln 𝑚 2 ≈ 1 . 6 × 10 6 Gpc 3 yr 𝑓 53 37 PBH 𝜂 -34 37 GLYPH<18> 𝑀 𝑀 ⊙ GLYPH<19> -32 37 GLYPH<18> 𝑡 𝑡 0 GLYPH<19> -34 37 𝑆 L 𝑆 E 𝜓 ( 𝑚 1 ) 𝜓 ( 𝑚 2 ) , (101)  \nwhere 𝑆 L and 𝑆 E are suppression factors given in Eqs. (46) and (60).', '2. Early three-body channel': 'd 𝑅 E3 d ln 𝑚 1d ln 𝑚 2 = 7 . 9 × 10 4 Gpc 3 yr 𝑓 144 𝛾 259 + 47 37 PBH GLYPH<20> 𝑡 𝑡 0 GLYPH<21> 𝛾 7 -1 GLYPH<20> ⟨ 𝑚 ⟩ 𝑀 ⊙ GLYPH<21> 5 𝛾 -32 37 𝑒 3 . 2 ( 1 -𝛾 ) 𝛾 28 / 9 -𝛾 K × GLYPH<18> 𝑀 2 ⟨ 𝑚 ⟩ GLYPH<19> 179 𝛾 259 -2122 333 ( 4 𝜂 ) -3 𝛾 7 -1 ¯ F( 𝑚 1 , 𝑚 2 ) 𝜓 ( 𝑚 1 ) 𝜓 ( 𝑚 2 ) , (102)  \nwhere K and ¯ F are given by Eqs. (81,80) and account for 3-body dynamics. Numerical simulations suggest that K ≈ 4 . 0. The second line contains the effect of extended mass functions.', '3. Late two-body channel': 'd 𝑅 L2 d ln 𝑚 1d ln 𝑚 2 ≈ 3 . 4 × 10 -6 Gpc 3 yr 𝑓 2 PBH 𝛿 eff GLYPH<18> 𝜎 𝑣 km / s GLYPH<19> -11 7 𝜂 -5 7 𝜓 ( 𝑚 1 ) 𝜓 ( 𝑚 2 ) , (103)  \nwhere 𝛿 eff and 𝜎 𝑣 characterize the PBH density contrast and velocity dispersion in the structures.', '4. Late three-body channel': "d 𝑅 L3 d ln 𝑚 1d ln 𝑚 2 ≈ 1 . 3 × 10 -16 𝑒 -6 . 0 ( 𝛾 -1 ) Gpc 3 yr 𝑓 3 PBH 𝛿 2 eff GLYPH<18> 𝜎 𝑣 km / s GLYPH<19> -9 + 8 𝛾 7 × 𝜂 -1 + 𝛾 7 GLYPH<18> 𝑀 𝑀 ⊙ GLYPH<19> 3 -𝛾 7 GLYPH<18> 𝑡 𝑡 0 GLYPH<19> 𝛾 7 F GLYPH<18> ⟨ 𝑚 ⟩ 2 𝜂𝑀 𝜅 min GLYPH<19> 𝜓 ( 𝑚 1 ) 𝜓 ( 𝑚 2 ) , (104)  \nwhere 𝛿 eff and 𝜎 𝑣 characterize the PBH density contrast and velocity dispersion in the structures, F( 𝜅 ) is given in Eq. (100) and 𝛾 determines the initial angular momentum distribution (97) of the binaries.  \nTo demonstrate the relative importance of different PBH merger rate channels for different values of 𝑓 PBH, we present numerical examples of the merger rates (101), (102), (103), (104) in Fig. 8. For concreteness, we have chosen a monochromatic mass function with 𝑚 = 30 𝑀 ⊙ . It is evident that for 𝑓 PBH ≲ 0 . 1  \nFig. 8 Numerical estimates of the four channels contributing to the PBH merger rate as functions of the fraction of DM in PBHs. We assume a monochromatic mass function with 𝑚 = 30 𝑀 ⊙ (upper panel) and 𝑚 = 0 . 1 𝑀 ⊙ (lower panel). The bands show variations due to uncertainties in the angular momentum distribution. We show optimistic late universe merger rates with 𝛿 eff = 10 5 for the late two-body and 𝛿 eff = 10 7 for the late three-body channels. The gray dashed lines show the LVK upper bounds from [29] and [30].  \n<!-- image -->  \nthe early two-body channel (101) dominates. This channel is also best understood and has the smallest uncertainties. The situation may change when 𝑓 PBH ≳ 0 . 1, that is, when a large fraction of DM of the Universe should be in the form of PBHs. In this case, the PHB density is large and the early three-body channel (102) may be the dominant one. However, the uncertainties associated with this solution are larger than in the first case. The upper (solid) and lower (dashed) boundaries of the three-body channels correspond to the merger rate with 𝛾 = 1 and K = 4 and 𝛾 = 2 and K = 1, respectively. The late Universe channels (103) and (104) are always expected to be subdominant. Optimistically, in Fig. 8 we have fixed 𝜎 𝑣 = 1 km / s and  \nused 𝛿 eff = 10 5 for the late two-body and 𝛿 eff = 10 7 for the late three-body channel where we also show 𝛾 = 1 (upper solid) and 𝛾 = 2 (lower dashed) cases. Yet the late-time channels cannot be enhanced enough to be dominant.  \nThe current upper bound from GW observations by LIGO-Virgo-KAGRA (LVK) observational run O3 is indicated by the gray dashed lines in Fig. 8. Comparing it with the predicted merger rate will implies a constraint on 𝑓 PBH. However, the observational bound on the merger rate is model dependent - it is affected by the characteristics of the binary population, e.g., the mass and redshift dependencies in the merger rates. Thus a dedicated analysis is needed for accurate GW constraints on 𝑓 PBH and Fig. 8 can at best indicate the order of magnitude. Given the present sensitivity to 𝑓 PBH, both early channels can be relevant to observations. The early two-body channel constraints 𝑓 PBH < O( 10 -3 ) for heavier 𝑚 ∈ 1 -100 𝑀 ⊙ PBHs, while the early three-body channel implies that the abundance sub-solar mass PBHs should not exceed O( 10% ) [31, 32, 33].  \nInitial clustering of PBHs affects /one.sup/three.sup , in particular, the early Universe binary formation channels. A rough estimate of the effect of initial clustering can be obtained by considering a constant 2-point function at the scales relevant for binary formation, 1 + 𝜉 ( x ) = 𝛿 . As this is equivalent to a local change in the PBH number density, its effect on the merger rate is a simple rescaling: 𝑅 ( 𝑓 PBH ) → 𝛿 -1 𝑅 ( 𝛿 𝑓 PBH ) . However, this scaling does not hold for later processes, such as disruption in haloes, that do not depend on the local PBH density at their formation. The merger rate of the early three-body channel is proportional to a higher power of 𝑓 PBH than that of the early two-body channel so it is also more sensitive to the initial clustering. Moreover, the initial clustering would enhance the disruption probability of the early binaries, suppressing 𝑅 𝐸 2. On the other hand, 𝑅 𝐸 3 gets enhanced even more as the three-body encounters are more likely.  \nAcknowledgments. This work was supported by the European Regional Development Fund through the CoE program grant TK202 and by the Estonian Research Council grants PRG803, PSG869, RVTT3 and RVTT7. The work of V.V. was partially supported by the European Union's Horizon Europe research and innovation program under the Marie Sklodowska-Curie grant agreement No. 101065736.", 'References': "- 1. K. Kritos, V. De Luca, G. Franciolini, A. Kehagias, and A. Riotto. The Astro-Primordial Black Hole Merger Rates: a Reappraisal. JCAP , 05:039, 2021.\n- 2. P. C. Peters. Gravitational Radiation and the Motion of Two Point Masses. Phys. Rev. , 136:B1224-B1232, 1964.\n- 3. Takashi Nakamura, Misao Sasaki, Takahiro Tanaka, and Kip S. Thorne. Gravitational waves from coalescing black hole MACHO binaries. Astrophys. J. Lett. , 487:L139-L142, 1997.\n- 4. N. Aghanim et al. Planck 2018 results. VI. Cosmological parameters. Astron. Astrophys. , 641:A6, 2020. [Erratum: Astron.Astrophys. 652, C4 (2021)].\n- 5. Martti Raidal, Christian Spethmann, Ville Vaskonen, and Hardi Veermae. Formation and Evolution of Primordial Black Hole Binaries in the Early Universe. JCAP , 02:018, 2019.\n- 6. Yacine Ali-Ha¨ımoud, Ely D. Kovetz, and Marc Kamionkowski. Merger rate of primordial black-hole binaries. Phys. Rev. D , 96(12):123523, 2017.\n- 7. Derek Inman and Yacine Ali-Ha¨ımoud. Early structure formation in primordial black hole cosmologies. Phys. Rev. D , 100(8):083528, 2019.\n- 8. Alex Hall, Andrew D. Gow, and Christian T. Byrnes. Bayesian analysis of LIGO-Virgo mergers: Primordial vs. astrophysical black hole populations. Phys. Rev. D , 102:123524, 2020.\n- 9. Gerald D. Quinlan. The Time scale for core collapse in spherical star clusters. New Astron. , 1:255, 1996.\n- 10. N. Afshordi, P. McDonald, and D. N. Spergel. Primordial black holes as dark matter: The Power spectrum and evaporation of early structures. Astrophys. J. Lett. , 594:L71-L74, 2003.\n- 11. Ville Vaskonen and Hardi Veermae. Lower bound on the primordial black hole merger rate. Phys. Rev. D , 101(4):043015, 2020.\n- 12. J. H. Jeans. The origin of binary systems. MNRAS, 79:408, April 1919.\n- 13. Steinn Sigurdsson and E. S. Phinney. Binary-Single Star Interactions in Globular Clusters. ApJ, 415:631, October 1993.\n- 14. John M. Fregeau, P. Cheung, S. F. Portegies Zwart, and F. A. Rasio. Stellar collisions during binary - binary and binary - single star interactions. Mon. Not. Roy. Astron. Soc. , 352:1, 2004.\n- 15. D. C. Heggie. Binary evolution in stellar dynamics. MNRAS, 173:729-787, December 1975.\n- 16. J. G. Hills. Encounters between binary and single stars and their effect on the dynamical evolution of stellar systems. AJ, 80:809-825, October 1975.\n- 17. Gabriele Franciolini, Konstantinos Kritos, Emanuele Berti, and Joseph Silk. Primordial black hole mergers from three-body interactions. Phys. Rev. D , 106(8):083529, 2022.\n- 18. Yao-Yuan Mao, Louis E. Strigari, Risa H. Wechsler, Hao-Yi Wu, and Oliver Hahn. Halo-toHalo Similarity and Scatter in the Velocity Distribution of Dark Matter. Astrophys. J. , 764:35, 2013.\n- 19. Gerald D. Quinlan and Stuart L. Shapiro. Dynamical Evolution of Dense Clusters of Compact Stars. ApJ, 343:725, August 1989.\n- 20. Hideaki Mouri and Yoshiaki Taniguchi. Runaway merging of black holes: analytical constraint on the timescale. Astrophys. J. Lett. , 566:L17-L20, 2002.\n- 21. Simeon Bird, Ilias Cholis, Julian B. Mu˜noz, Yacine Ali-Ha¨ımoud, Marc Kamionkowski, Ely D. Kovetz, Alvise Raccanelli, and Adam G. Riess. Did LIGO detect dark matter? Phys. Rev. Lett. , 116(20):201301, 2016.\n- 22. S. J. Aarseth and D. C. Heggie. The probability of binary formation by three-body encounters. A&A, 53:259-265, December 1976.\n- 23. Natalia Ivanova, K. Belczynski, J. M. Fregeau, and F. A. Rasio. The Evolution of binary fractions in globular clusters. Mon. Not. Roy. Astron. Soc. , 358:572-584, 2005.\n- 24. N. Ivanova, S. Chaichenets, J. Fregeau, C. O. Heinke, J. C. Lombardi, Jr., and T. Woods. Formation of black-hole X-ray binaries in globular clusters. Astrophys. J. , 717:948-957, 2010.\n- 25. Meagan Morscher, Bharath Pattabiraman, Carl Rodriguez, Frederic A. Rasio, and Stefan Umbreit. The Dynamical Evolution of Stellar Black Holes in Globular Clusters. ApJ, 800(1):9, February 2015.  \n- 26. Valeriya Korol, Ilya Mandel, M. Coleman Miller, Ross P. Church, and Melvyn B. Davies. Merger rates in primordial black hole clusters without initial binaries. Mon. Not. Roy. Astron. Soc. , 496(1):994-1000, 2020.\n- 27. Carl L. Rodriguez et al. Modeling Dense Star Clusters in the Milky Way and beyond with the Cluster Monte Carlo Code. Astrophys. J. Supp. , 258(2):22, 2022.\n- 28. Nicholas C. Stone and Nathan W. C. Leigh. A statistical solution to the chaotic, non-hierarchical three-body problem. Nature, 576(7787):406-410, December 2019.\n- 29. R. Abbott et al. Population of Merging Compact Binaries Inferred Using Gravitational Waves through GWTC-3. Phys. Rev. X , 13(1):011048, 2023.\n- 30. R. Abbott et al. Search for Subsolar-Mass Binaries in the First Half of Advanced LIGO's and Advanced Virgo's Third Observing Run. Phys. Rev. Lett. , 129(6):061104, 2022.\n- 31. Gert Hutsi, Martti Raidal, Ville Vaskonen, and Hardi Veermae. Two populations of LIGO-Virgo black holes. JCAP , 03:068, 2021.\n- 32. Gabriele Franciolini, Ilia Musco, Paolo Pani, and Alfredo Urbano. From inflation to black hole mergers and back again: Gravitational-wave data-driven constraints on inflationary scenarios with a first-principle model of primordial black holes across the QCD epoch. Phys. Rev. D , 106(12):123526, 2022.\n- 33. M. Andr'es-Carcasona, A. J. Iovino, V. Vaskonen, H. Veermae, M. Mart'ınez, O. Pujol'as, and Ll. M. Mir. Constraints on primordial black holes from LIGO-Virgo-KAGRA O3 events. 5 2024.\n- 34. Valerio De Luca, Gabriele Franciolini, Antonio Riotto, and Hardi Veermae. Ruling Out Initially Clustered Primordial Black Holes as Dark Matter. Phys. Rev. Lett. , 129(19):191302, 2022."}