Enhanced cosmological perturbations and the merger rate of PBH binaries

TL;DR

This paper explores how cosmological perturbations influence primordial black hole binary merger rates, showing that enhanced small-scale power spectra can significantly affect observational constraints and provide insights into early universe fluctuations.

Contribution

It introduces a detailed analysis of the impact of cosmological perturbations on PBH binary mergers, highlighting the role of small-scale power spectrum enhancements and their observational implications.

Findings

Cosmological perturbations can dominate the torque on PBH binaries, affecting their initial eccentricity.

Enhanced small-scale power spectra relax current bounds on PBH dark matter fraction.

Merger rates can reveal information about primordial perturbations at uncharted scales.

Abstract

The rate of merger events observed by LIGO/Virgo can be used in order to probe the fraction $f$ of dark mater in the form of Primordial Black Holes. Here, we consider the merger rate of PBH binaries, accounting for the effect of cosmological perturbations on their initial eccentricity $e$. The torque on the binaries receives significant contributions from a wide range of scales, that goes from the size of the horizon at the time when the binary forms, down to the co-moving size of the binary. In scenarios where PBH are formed from adiabatic perturbations, it is natural to expect an enhancement of the power spectrum $P_\Phi$ at small scales, where it is poorly constrained observationally. The effect can then be quite significant. For instance, a nearly flat spectrum with amplitude $P_\Phi \gtrsim 10^{-7}$ on scales smaller than $\sim 10 Mpc^{-1}$ gives a contribution $\langle j^2 \rangle \sim 10^3 P_\Phi$, where $j= (1-e^2)^{1/2}$ is the dimensionless angular momentum parameter of the binaries. This contribution can dominate over tidal torques from neighboring PBHs for any value of $f$. Current constraints allow for a power spectrum as large as $P_\Phi \sim 10^{-5}$ at the intermediate scales $10^3-10^5 Mpc^{-1}$, comparable to the co-moving size of the binaries at the time of formation. In particular, this can relax current bounds on the PBH abundance based on the observed LIGO/Virgo merger rate, allowing for a fraction $f\sim 10\%$ of dark matter in PBH of mass $\sim 30 M_\odot$. We investigate the differential merger rate $\Delta\Gamma(m_1,m_2)$, as a function of the masses of the binary components, and the corresponding ``universality" coefficient $\alpha = -(m_1+m_2)^2 \partial^2 \ln \Delta\Gamma/\partial m_1\partial m_2$, concluding that merger rates may provide valuable information on the spectrum of primordial cosmological perturbations at currently uncharted lengthscales.

Equations

$$dP=e^{-X}dX.$$

$$X=nV={\left(\frac{x}{\bar{x}}\right)}^3,$$

$$\bar{x}={\left(\frac{3\bar{m}}{2\pi f{\rho }_{eq}}\right)}^{1/3}.$$

$$\frac{1}{2}\mu H^2{s}^2{x}^2\sim \frac{G{m}_1{m}_2}{sx}.$$

$$\lambda \equiv \frac{X}{f}.$$

$$x=|\vec{x}|=\frac{\chi (s)}{s}{x}_0.$$

$$\chi (s;\lambda )=\lambda \chi (s/\lambda ;1).$$

$$s\approx s_b=\frac{\lambda }{3}\left(\frac{2\bar{m}}{M}\right),$$

$$a=3\beta s_{b}{x\over 2}={\beta\over 2}\left({2\bar{m}\over M}\right)\Big{(}\frac{3\bar{m}}{2\pi\rho_{eq}}\Big{)}^{1/3}\lambda^{4/3},$$

$$m\equiv \frac{\bar{m}}{M_{\odot }}.$$

$$a\approx 1.8\cdot 10^{-7}{\lambda }^{4/3}{m}^{1/3}\left(\frac{2\bar{m}}{M}\right)H_{eq}^{-1},$$

$$\Delta t\sim a^2(GMa{)}^{-1/2}\sim H_b^{-1}\lesssim t_{eq},$$

$$t=t[j,a]\equiv \frac{3}{85}\frac{a^4}{G^3{m}_1{m}_{2}M}{j}^7,$$

$$j\equiv \frac{\ell }{\sqrt{GMa}}.$$

$$j=j_0(\lambda )\approx 1.6\cdot 10^{-3}{\lambda }^{-16/21}{m}^{5/21}\left(\frac{M}{2\bar{m}}\right){\left(\frac{4m_1{m}_2}{M^2}\right)}^{1/7}.$$

$$d{n}_{bin}=d{n}_M{d}^2\mathcal{F}.$$

$$d^2\mathcal{F}\equiv d\mathcal{F}(m_1)d\mathcal{F}(m_2),$$

$$d{n}_M=\Theta \left(M-\lambda \bar{m}\right)f^2\frac{{\rho }_m(t_0)}{2\bar{m}}{e}^{-f\lambda }dP(j;\lambda )d\lambda ,$$

$$\lambda \lesssim \frac{M}{\bar{m}}.$$

$$d\Gamma (t_0,m_1,m_2)={\Gamma }_M(t_0)d^2\mathcal{F}$$

$${\Gamma }_M(t_0)=\int \delta \left(t_0-t[j,a]\right)d{n}_M.$$

$${\Gamma }_M(t_0)=f^2\frac{{\rho }_m}{14\bar{m}{t}_0}{\int }_{{\lambda }_{min}}^{M/\bar{m}}W(\lambda )d\lambda .$$

$$W(\lambda )=j{\frac{dP(j;\lambda )}{dj}|}_{j=j_0},$$

$${\lambda }_{min}<\lambda \lesssim M/\bar{m},$$

$$j=|{\vec{ȷ}}_{nb}+{\vec{ȷ}}_{cp}|$$

$$\vec{ȷ}\sim (\vec{x}\times \Delta \vec{g})\frac{\Delta t}{\sqrt{GMa}},$$

$$\Delta \vec{g}=\vec{\nabla }\Phi ({\vec{x}}_1)-\vec{\nabla }\Phi ({\vec{x}}_2),$$

$$j=|\vec{ȷ}|=\gamma \frac{j_X}{Y},$$

$$j_X=(\bar{m}/M)X.$$

$$dP(Y)=e^{-Y}dY.$$

Full text

1 Introduction

The detection of gravitational waves (GWs) from merging black hole binary systems has revived interest in the idea that primordial black holes (PBHs) may be a viable candidate for dark matter (DM). The abundance of PBHs is severely constrained for a wide range of masses (see e.g. [2] and references therein), but it could still be significant both for sublunar and stellar masses. In particular there is an active debate on wheter PBHs in the mass range recently detected by LIGO/Virgo collaboration could account for a sizable fraction of DM [3, 4].

The observed merger rate [5, 6, 7, 8, 9] $\Gamma\approx 10-100\ Gpc^{-3}yr^{-1}$ in the range $\sim 5-100M_{\odot}$ has recently been used in order to place limits on the PBH abundance [4, 10]. In such estimates, it has been assumed that PBH are spatially uncorrelated at the time of formation, and that the dominant contribution to the orbital angular momentum of the binaries originates from tidal forces exerted by other black holes in the neighborhood, around the time when the binary decouples from the Hubble flow [12, 13, 1]. With these assumptions, the observational upper bound on the merger rate limits the fraction of PBHs in DM to $f\lesssim 1\%$. Several refinements to this estimate have been considered, including initial spatial correlations of the PBHs [14, 15, 16, 17], tidal forces from non-relativistic matter perturbations [11, 10, 18], as well as the effect of a dark matter dress around the PBHs [19], with similar results for the bound on the PBH abundances. In Ref. [20], the effect of infalls of neighboring PBHs on the binary has been studied, with the conclusion that this may significantly reduce the observed merger rate. Also, N-body simulations for the formation and evolution of binaries [20] indicate that for high $f\sim 1$, the rate may be significantly reduced by disruption, through the interaction of binaries with compact N-body systems.

In this paper, we consider the effect of primordial cosmological perturbations on the angular momentum of PBH binaries during the radiation dominated era. We note that, even for a scale invariant spectrum of density perturbations, there is a wide range of scales contributing to the torque. Moreover, the amplitude of the power spectrum $P_{\Phi}$ is poorly constrained beyond the scale of $10Mpc^{-1}$, and it could be significantly larger than it is on cosmological scales. In the present context this possibility seems rather natural, since some scenarios for PBH formation111Not all scenarios for PBH formation require a bump in the power spectrum. For instance, PBH could be created by active seeds such as relic domain walls or false vacuum bubbles produced during inflation [29, 30, 31], rather than adiabatic perturbations. In such alternative scenarios, an enhancement in the spectrum of cosmological perturbations does not seem to be a necessary feature. The same is true for PBH formation at post-inflationary phase transitions (see e.g. [32, 33, 34, 35] and references therein). rely on a prominent enhancement or “bump” in the power spectrum at relatively short wavelengths, corresponding to the co-moving size of the horizon at the time when PBHs form. In the inflationary context, the height and location of the bump depend on specific features in the inflaton potential. For instance, in one field models, the field may undergo a short period of ultra-slow roll or constant roll as it encounters local extrema on its way down the potential [21, 22, 23, 24, 25, 26, 27, 28]. While the amplitude of perturbations on cosmological scales is of order $P_{\Phi}\sim 10^{-9}$, the r.m.s. amplitude at the bump should be much larger, $P_{\Phi}\sim 10^{-3}-10^{-2}$, so that PBHs can form in significant abundance. This is a strong departure from scale invariance, and it seems plausible that in generic models of this sort the power spectrum might be enhanced also at the scales interpolating from the cosmological plateau down to the PBH scale, including the intermediate scales comparable to the co-moving size of the binaries. Here, we shall be agnostic about the specific inflationary dynamics, and will simply explore the consequences of an enhanced spectrum which we shall model as a (piecewise) power law $P_{\Phi}(k)\propto k^{p}$. As we shall see, such an enhancement may have potentially observable consequences. In particular, it may affect the differential merger rate of binaries as a function of the component masses.

The paper is organized as follows. In Section 2 we briefly review the formation of PBH binaries and the distribution of orbital parameters, taking into consideration the effect of neighbouring black holes but ignoring cosmological perturbations. It is in this context that the universality coefficient $\alpha$, which charaterizes the dependence of the merger rate on the masses of the components, was first introduced [1]. Hence, this will be a useful reference case. We also comment on PBH infalls and their effect on $\alpha$.

In Section 3 we discuss the effect of cosmological perturbations on the dimensionless orbital angular momentum parameter $j$. In contrast with earlier analysis, here we include the perturbations in radiation, whose effect dominates over that of matter perturbations for binaries which decouple from the Hubble flow deep in the radiation era. For an enhanced $P_{\Phi}$ these tend to dominate the distribution of $j$.

In Section 4 we consider the merger rates in three different scenarios: the nearly scale invariant cosmological plateau (Case A), an enhanced spectrum at intermediate scales with a moderate spectral index $0<p\lesssim 2$ (Case B), and a rather steep power spectrum $p>2$, peaked at scales smaller than the binary size (Case C). Our conclusions are summarized in Section 5.

Throughout the paper $f$ will denote the fraction of dark matter in the form of PBHs, and $s$ will denote the cosmological scale factor, while $a$ will denote the semi-major axis of binaries. We adopt the convention that $s=1$ at the time of matter-radiation equality. The speed of light is set to $c=1$.

2 PBH binary formation and universality

In this Section, we briefly review the case where the angular momentum of binaries is due to the tidal torque from other PBHs in the vicinity, neglecting cosmological perturbations. We also introduce the universality coefficient $\alpha$ [1], and we discuss how this may be affected by the infall of neighboring PBH on binaries.

2.1 Initial orbital parameters and the life-time of binaries

Following [4, 12], let us assume a uniform distribution of PBHs, without any initial spatial correlations.222The effect of such an initial correlation has been discussed in Refs. [15, 17]. From a given PBH, the probability of finding the nearest neighbour at a certain distance is given by

[TABLE]

Here $X=nV$ is the product of the co-moving number density $n$ times the co-moving volume $V=(4/3)\pi x^{3}$, where $x$ is the co-moving distance. We adopt the convention that the cosmological scale factor is $s=1$ at the time of matter-radiation equality.

We shall also assume that the PBH mass function is not too broad 333This is expected when PBH are formed from very high peaks of a Gaussian random field of density perturbations [37, 36], even if the enhancement in the power spectrum has a sizable width. Unless the power spectrum involves different explicit scales, high peaks of the random field tend to have a well defined shape, which leads to a relatively narrow range of masses after gravitational collapse., allowing however for some spread in the masses within an order of magnitude or so. The co-moving number density takes the form $n=f\rho_{eq}/(2\bar{m})$, where $f$ is the fraction of DM in the form of PBHs, $\rho_{eq}$ is the density at the time of equality, and $\bar{m}$ is the average mass in the distribution. We may then write

[TABLE]

where

[TABLE]

In a spherical region of radius $\bar{x}$ we expect to find one PBH, on average, so the length scale $\bar{x}$ can also be thought of as a typical separation between PBHs.

A pair of black holes forms a binary when the relative kinetic energy due to the Hubble flow becomes comparable to the gravitational binding energy between the two objects [4],

[TABLE]

Here, $\mu=m_{1}m_{2}/M$ is the reduced mass, where $M=m_{1}+m_{2}$ is the total mass of the binary. The above relation has to be satisfied before the end of the radiation era, since both sides will scale as $s^{-1}$ during matter domination. For $x\lesssim\bar{x}$, and taking into account that $\rho\approx\rho_{eq}/(2s^{4})$ in the radiation era ($s\ll 1$), the relation (2.3) is satisfied when the cosmological scale factor $s$ is of order $s\sim\lambda(\bar{m}/M)\leq 1,$ where we have introduced

[TABLE]

More precisely, parametrizing the physical distance as $d=\chi(\eta;\lambda)x_{0}$, where $x_{0}$ is the initial co-moving separation, we may write

[TABLE]

The numerical analysis in Refs. [10, 20] shows that, that for binaries forming deep in the radiation era (i.e. $\lambda\ll 1$), the function $\chi$ is self-similar

[TABLE]

Initially, the two PBH are following the Hubble flow, so that $\chi\approx s$, $\vec{x}\approx\vec{x}_{0}$ is approximately constant, and the physical distance grows linearly in $s$. However, when the scale factor reaches the value

[TABLE]

the physical distance turns around and a bound system is formed with semi-major axis given by [10]

[TABLE]

where $\beta\approx 0.2$. Introducing the dimensionless average mass parameter

[TABLE]

we have

[TABLE]

where we have used $H_{eq}^{-1}\approx 0.9\cdot 10^{18}km$.

In an environment with no external forces and torques, two PBHs which are initially at rest would collide head-on due to gravitational attraction in a very short time-scale

[TABLE]

comparable to the Hubble radius $H_{b}^{-1}$ at the time when the binary forms. However, the binary system is immersed in a local tidal field, created by density perturbations and by other PBHs in the neighborhood. These forces will exert a torque on the binary, giving it an orbital angular momentum which avoids the head-on collision. The binary will then slowly radiate its energy by emitting gravitational waves in a much longer timescale, before the final merger occurs.

For a binary with initial orbital angular momentum $\ell$ per unit reduced mass, the life-time is given by Peters formula [38]

[TABLE]

where dimensionless parameter $j$, is defined as

[TABLE]

For an elliptic orbit with semi-minor axis $b$, we have $j=b/a=\sqrt{1-e^{2}}$, where $e$ is the eccentricity. Note that $j<1$. Using $t=t_{0}\approx 1.3\cdot{10^{23}}~{}km$, we find that the binaries which are merging today are characterized by

[TABLE]

where we have used Eq. (2.8) with $\beta\approx 0.2$. The distribution of $j$ in the ensemble of binaries depends on the specific mechanisms which give the binaries their angular momentum.

2.2 Merger rates

In general, the differential number density of binaries per unit volume is given by

[TABLE]

Here

[TABLE]

where $d{\cal F}$ is the PBH mass distribution function, and $dn_{M}(m_{1},m_{2},X,j)$ is the distribution of binaries with masses $m_{1}$ and $m_{2}$, initial separation of the partners characterized by $X$ and orbital angular momentum parameter $j$. The variable $X$ is distributed as (2.1), so using $X=f\lambda$, we have

[TABLE]

where $dP(j;\lambda)$ is the distribution of $j$ for given $\lambda$, $f\rho_{m}(t_{0})/\bar{m}$ is the number density of PBH at the present time, $\rho_{m}$ is the current matter density, and we have inserted a factor of $1/2$ to avoid double counting of binaries. The Heavyside function restrics the range of $\lambda$ since, according to our earlier discussion around Eq. (2.3), a given PBH will only be part of a binary if the distance to the nearest PBH satisfies

[TABLE]

Otherwise the Hubble flow velocity always remains larger than the binding energy. If the distribution of masses is not too wide, we have $M\sim 2\bar{m}$, and for $f\ll 1$ the exponential factor $e^{-f\lambda}$ can be approximated by $1$.

The intrinsic merger rate444Here, and for the rest of this paper, we consider intrinsic merger rates, ignoring effects due to time delay of events which occur at high redshift. These can be incorporated along the lines of Ref. [20]. of PBH binaries per unit time and volume can be written as

[TABLE]

where the rate at fixed total mass $M=m_{1}+m_{2}$ is given by

[TABLE]

Using (2.12) in the argument of the delta function, we can perform the $j$ integration to obtain

[TABLE]

Here ,we have dropped the factor $e^{-f\lambda}$, since as mentioned above this can be approximated by unity in the relevant range of parameters. Also, we have introduced

[TABLE]

where $j_{0}$ given by Eq. (2.14) is the value of $j$ for which the lifetime of the binary coincides with the present age of the universe $t_{0}$. The integral in (2.21) is in the range

[TABLE]

where the lower limit555In the expressions which are given only by order of magnitude, and in the interest of brevity, we will often omit the explicit dependence on the individual masses, assuming they are within one order of magnitude or so from each other. $\lambda_{min}\sim 2\cdot 10^{-4}m^{5/16}$ is determined from (2.14), taking into account that we must have $j_{0}\leq 1$.

The distribution $dP(j;\lambda)$ depends on the mechanism which gives angular momentum to the binaries. In general,

[TABLE]

is the added contribution from torques due to neighboring PBHs, and from torques due to cosmological perturbations. For the rest of this Section we concentrate on $j_{nb}$, while the effect of cosmological perturbations will be discussed in the following Sections.

2.3 Tidal torque due to neighboring black holes

Let us start by considering the effect of a single neighboring PBH, producing a tidal torque on the binary. By integrating the torque over time, we have,

[TABLE]

where $\Delta t$ is given in Eq. (2.11). Here, and for the rest of this section, we suppress the subscript $nb$ from $\vec{\jmath}$, since we are only dealing with the effect of neighboring PBHs. In Eq. (2.25), $\vec{x}=\vec{x}_{2}-\vec{x}_{1}$ is the relative co-moving separation between the members of the binary and

[TABLE]

is the tidal acceleration, expressed in terms of the Newtonian potential $\Phi$ created by the neighboring PBH. If the gradients vary on a length scale much larger than the separation $\vec{x}$, the tidal acceleration can be expanded in powers of $\vec{x}$ and the leading term is given by $\Delta g^{k}=-\Phi,_{kl}x^{l}$, where the spatial derivatives are with respect to the co-moving coordinates. If the 3rd black hole has a mass $m_{3}$, and is at a co-moving distance $y\gtrsim x$, then $\Phi\approx Gm_{3}/(s_{b}x)$, where $s_{b}$ is the scale factor around the time when the binary forms. The co-moving gradients can then be estimated as $|\Phi_{,ij}|\sim(Gm_{3}/s_{b}y^{3})(3y^{i}y^{j}-\delta_{ij})$, and Substituting in (2.25) with $a\sim s_{b}x$ we obtain

[TABLE]

Here we have introduced the variable $Y\equiv(y/\bar{x})^{3}$, characterizing the distance to the nearest third PBH, and

[TABLE]

Since the position of the third black hole is random, the distribution for $Y$ is also given by

[TABLE]

The coefficient $\gamma$ is given by

[TABLE]

where $\theta_{3}$ is the angle between the relative coordinate $\vec{x}$ and the position of the third black hole $\vec{y}$, that is, $\cos\theta_{3}=\hat{x}\cdot\hat{y}$. The overall numerical factor in this expression is determined by taking into consideration the time dependence in Eq. (2.5) as the binary forms [10, 13, 1, 20]. Note that the average of $\gamma$ over angle and mass distribution function is

[TABLE]

Following [1], we keep $\gamma$ as an undetermined random variable of order one, over which we can integrate at the end of the computation, if needed.

Taking into consideration (2.27) and (2.29), we have

[TABLE]

Here, the superindex in the probability distribution indicates that, for the time being, we are considering the effect of the nearest neighboring PBH only. Then, we have

[TABLE]

where

[TABLE]

The characterisctic value

[TABLE]

is essentially the peak of the function $W^{(1)}(\lambda)$. The rate (2.21) is plotted in Fig. 1 for a range of values of $f$, and for different values of the mass $\bar{m}$, assuming that $m_{1}=m_{2}=\bar{m}$ and $\gamma\approx 1$. The curves have a knee which separates two different regimes with a power law behaviour in $f$. This can easily be understood analytically. The behaviour of the merger rate depends on whether $\lambda_{*}$ is large or small, and this in turn depends on the value of $f$.

The behaviour of the merger rate below the knee ($\lambda_{*}\gg 1$) corresponds to a low fraction of DM in PBH, $f\ll f_{*}$, where

[TABLE]

Since $\lambda\lesssim M/\bar{m}\sim 1\ll\lambda_{*}$, we have $Y_{0}\ll 1$, and we may then neglect the exponential dependence of the integrand in (2.21),

[TABLE]

Since the lower limit of integration $\lambda_{min}\sim 10^{-4}$ does not play a role, we have set it to zero for simplicity. This leads to the estimate

[TABLE]

where we have used

[TABLE]

For an approximately monochromatic PBH mass function, the observational bound $\Gamma\lesssim 10^{2}Gpc^{-3}yr^{-1}$ then leads to

[TABLE]

We conclude that, if the third black hole is the dominant source of orbital angular momentum, then solar mass black holes, with $m\sim 1$, can only account for a very small fraction of dark matter, with $f\lesssim 10^{-3}$. This is in agreement with the analysis of Refs. [4, 1]. Note that, even in the case $m\sim 1$ the upper limit of the observational bound (2.39) satisfies the condition $f\lesssim f_{*}$ only marginally.

Hence, let us now consider the complementary limit $f\gtrsim f_{*}$. For $\lambda\gg\lambda_{*}$ we have $Y_{0}(\lambda)=\gamma(\lambda/\lambda_{*})^{37/21}\gg 1$, and due to the factor $e^{-Y_{0}}$ the integral (2.21) is effectively cut-off at $\lambda=\lambda_{*}\ll 1$. Therefore it is a good approximation to remove the upper limit of integration, which doesn’t play a role, and then the integral scales as $\lambda_{*}$,

[TABLE]

The rate can then be approximated as

[TABLE]

The observational bound $\Gamma\lesssim 10^{2}Gpc^{-3}yr^{-1}$ then leads to the condition

[TABLE]

where, as in Eq. (2.39), in this inequality we assume a nearly monochromatic PBH mass function. In the mass range of LIGO/Virgo detections, $m\sim 30$, the bound on the fraction $f$ of DM in PBH is limited to $\sim 0.5\%$, again in good agreement with [4, 1]. The analytic estimate (2.41) is plotted in Fig. 1 as a dotted line for $m=1$.

The previous considerations can be extended to the case where we include the torque of all neighboring black holes, and not just the closest one. An expression for $dP^{(\infty)}(j)$ due to the cumulative effect of all PBHs in the neighborhood was derived in [10, 20, 1]. This has the form of a power law distribution with a break at $j_{X}\equiv(\bar{m}/M)X=(\bar{m}/M)f\lambda$:

[TABLE]

Note that at large $j$, the behaviour of (2.43) is similar to (2.31), where only the nearest PBH is considered. However, at small $j$ the distribution (2.31) vanishes exponentially in $1/j$, while (2.43) has the form $dP^{(\infty)}\propto jdj$. As pointed out in [1], the reason is that in the case of a single PBH, the only way to reduce the torque on the binary is to place the PBH sufficiently far. The probability for that decays exponentially in $1/j$ for large distance. On the other hand, when many neighboring PBHs are involved, their added torques may randomly produce a small effect, with a probability which is only phase space suppressed. The dimensionless angular momentum $\vec{\jmath}$ is in the plane orthogonal to the initial relative separation $\vec{x}$, so the corresponding measure is two dimensional $d^{2}\vec{\jmath}=2jdj$, and the behaviour $dP\propto jdj$ is expected. In conclusion, the distribution (2.31) does not provide a very good description at small $j$, even if it is true that the nearest PBH gives the dominant contribution to the torque.

Using (2.43) in (2.22) we have

[TABLE]

where

[TABLE]

is the same as $Y_{0}$ given in (2.33), with the coefficient $\gamma$ replaced by its averaged value over masses and directions, $\bar{\gamma}=1$.

At low $f$, where $\lambda_{*}\gg 1$, we have $\bar{Y}_{0}\ll 1$ throughout the range of integration in (2.21). Hence, $W\approx\bar{Y}_{0}$, and the rate will be given by Eq. (2.37). In other words, the inclusion of the effect of an infinite number of neighbours does not change the merger rate:

[TABLE]

Also, at high $f$, where $\lambda_{*}\ll 1$, the integral is dominated by $\lambda\sim\lambda_{*}\ll M/\bar{m}$, so we can approximate by extending the range of integration to infinity and evaluating in terms of Euler’s Gamma function. Then one finds

[TABLE]

and the difference between the two cases is only by a very small change in the overall numerical factor. In fact, the distributions (2.32) and (2.44) produce integrated merger rates which are almost indistiguishable from one another (also for intermediate values of $f\sim 10^{-3}$), in agreement with the notion that the nearest PBH gives the dominant contribution to the torque (See Fig. 1). More importantly, the dependence of the merger rate on binary masses is basically unaffected by the inclusion of an infinite number of neighbours. Let us now turn to the characterization of such mass dependence.

2.4 Universality in the mass dependence of the merger rates.

In principle, we cannot predict the mass dependence of the merger rates unless the initial mass distribution function $d{\cal F}(m_{i})$ is known. Unfortunately, the latter is model dependent. However, a very interesting observation was made in Ref. [1] which may bypass this difficulty. Noting that the rate in a given mass interval $\Delta m_{1}$, $\Delta m_{2}$ is given by $\Delta\Gamma(m_{1},m_{2})=\Gamma_{M}(m_{1},m_{2})\Delta{\cal F}(m_{1})\Delta{\cal F}(m_{2})$, the expression

[TABLE]

is independent of the unknown distribution funcion ${\cal F}$. It was argued in [1] that with a sufficiently large sample of PBH merger events, of order $10^{3}$, the coefficient $\alpha$ can be determined observationally with accuracy of order $15\%$. This makes it a very attractive observable, within reach of existing and upcoming gravitational wave detectors [39].

If neighboring PBHs are the only source of angular momentum for the binaries, the coefficient $\alpha$ can readily be found from the expressions (2.37) and (2.41). The powers of $m_{1}m_{2}$ in these expressions for the merger rate do not contribute to $\alpha$, since after taking the logarithm and the two derivatives with respect to $m_{1}$ and $m_{2}$ such terms drop out. This is the same reason why $\alpha$ does not depend on the initial mass distribution functions $\cal F$. The only contributions to $\alpha$ come from powers of the total mass $M$. Hence, from (2.37), we have

[TABLE]

and from (2.41) we have

[TABLE]

This leads to the prediction of a “hidden universality” in the merger rate [1], where the parameter $\alpha$ should be in the narrow range

[TABLE]

As we shall see, this universality coefficient can be altered by different effects [see Fig. 9]. This may convey useful information about the actual circumstances surrounding binary formation and evolution.

2.5 PBH infalls and the universality coefficient

Eq. (2.21) includes all binaries which have the appropriate initial conditions to merge at the present time, provided that they remain in isolation for the rest of cosmic history. It has been argued in [20] that this overestimates the rate, because some of the binaries may be affected by the infall of neighboring PBHs, which could disturb the eccentricity or even disrupt the binary. Note that binaries merging at the present time have a very low $j$, which is given by Eq. (2.14), $j\sim 10^{-3}$. PBH infalls would increase this value by a large factor, making the life-time of the binary much larger than the age of the universe. On the other hand, if the closest neighbour is at a distance such that $Y\gtrsim M/\bar{m}$, then the binary is in an underdense region, and the neighbouring PBHs will not decouple from the Hubble flow to fall onto it [13, 20].

To illustrate the potential impact of infalls as a function of $f$, let us start by using a crude approximation where only the closest PBH contributes to the torque. Within this approximation, we may compare the naive rate given by Eq. (2.21) with the merger rate of “pristine” binaries which are unaffected by infalls. This is achieved by restricting the integration to the range where, say,

[TABLE]

The result is illustrated in Fig. 2 (dashed lines). For that comparison, we assume a monochromatic mass spectrum, and the lower limit of the integral in Eq. (2.21) is taken to be the value of $\lambda$ that corresponds to $Y_{0}(\lambda)=M/\bar{m}=2$. For very low $f\sim 10^{-3}$, the lower limit of integration becomes of order one, and for lower $f$ the range of integration completely disappears. This causes the sudden drop of the dashed curves in Fig. 2. On the other hand the effect is not so dramatic for higher $f\gtrsim 10^{-3}$, amounting only to a moderate overall factor of order 10 or so. The reason is simple to understand. According to the analysis of Subsection 2.3, for $f$ above the knee in the curves of Fig. 1, the rates are dominated by binaries where the closest neighbor is at a distance $Y_{0}(\lambda_{*})\sim 1$. In this case the factor $e^{-Y_{0}(\lambda)}$ acts as the effective cut-off of the integral (2.21) at the value $\lambda_{*}\lesssim 1$. Hence, the effect of restricting the range of integration to $Y>Y_{min}=2$ is only a mild suppression by a factor of order $\exp[-(Y_{min}-Y(\lambda_{*}))]$, where the exponent is of order 1.

Since PBH infalls have a sizeable impact on the merger rates, it important to consider the overall effect of all neighboring PBH under the assumption that all of them are outside the basin of attraction of the binary. This issue was considered in generality in Ref. [20], where the distribution $dP(j;\lambda,Y_{min})$ was calculated, by taking into account all neighboring PBHs, and assuming that these are at a distance larger than $Y_{min}$ from the binary. The result was found in closed form in terms of integrals of Hypergeometric functions, and it is somewhat cumbersome in general. However, the expression greatly simplifies in the limits $Y_{min}\to 0$ and $Y_{min}\gg 1$, which are of our primary interest. For $Y_{min}\to 0$, Eq. (2.43) is recovered, as expected, since in this case we do not exclude any of the binaries from the count. On the other hand, for $Y_{min}\gg 1$, it was found that the distribution can be approximated as a Gaussian

[TABLE]

with variance

[TABLE]

Here, $K=6\langle m^{2}\rangle/(5\bar{m}^{2})\sim 1$, where the brackets indicate average over the mass distribution. The approximation is already quite accurate for $Y_{min}\gtrsim 2$, which is the range of our interest.

The differential merger rate is therefore given by

[TABLE]

where the factor $e^{-Y_{min}}$ accounts for the probability that the closest neighboring PBH is further than $Y_{min}$. Here, and in what follows, we restrict attention to the regime where $f\gtrsim f_{*}$, which seems most relevant for observations. Then, we can approximate (2.21) by removing the upper limit of integration and we have

[TABLE]

where we have used (2.45). The corresponding merger rate is plotted in Fig. 2 (solid lines), for different values of the mass.

Let us now consider the universality coefficient $\alpha$. If $Y_{min}$ were independent of the masses, then (2.57) would scale like $\lambda_{*}$, just like in Eq. (2.40). In that case, we would recover the value $\alpha=36/37$. However, since heavier binaries have a larger basin of attraction, $Y_{min}$ scales as (2.52), and we have666An analogous computation considering only the effect of the nearest PBH produces a somewhat different answer,

$\int Wd\lambda\approx\int Y_{0}e^{-Y_{0}}d\lambda\approx{21\over 37}\gamma^{-{21\over 37}}\lambda_{*}\int_{Y_{min}}^{\infty}Y_{0}^{21\over 37}e^{-Y_{0}}dY_{0}={21\over 37}\gamma^{-{21\over 37}}\lambda_{*}\Gamma(58/37,Y_{min}).$

(2.57)

Approximating $\Gamma(58/37,Y_{min})\approx 1.25\ Y_{min}^{21/37}e^{-Y_{min}}$ for $Y_{min}\gtrsim 2$, and using (2.52) we find $\alpha=(36+21)/37\approx 1.54$. However, it should be noted that here we are considering large $Y$, which corresponds to small $j$, and in this regime it is not a good approximation to neglect the contribution from all other PBHs, as explained in the paragraph following Eq. (2.43).

[TABLE]

The second term comes from $M^{21/74}$ in (2.56). It is interesting to note that the factor $e^{-Y_{min}}$, which is exponential in $M$, does not contribute to the parameter $\alpha$, since it drops out after taking two derivatives of $\log\Gamma$ with respect to the masses. We conclude that the effect of infalls of neighboring PBHs onto binaries produces a significant shift of $\alpha$ towards a higher value. 777Aside from infalls, the simulations in Ref [20] also indicate that, for $f\gtrsim 0.1$, binaries can be disrupted during the matter dominated era by interaction with compact N-body systems. This effect can be particularly important for $f\sim 1$, where a sizable fraction of the binaries undergo interactions even before the time of recombination. This effect is also likely to suppress the rates at high $f$, and further work is needed to assess what fraction of the binaries may ultimately remain unaffected. In what follows, we shall simply ignore this possibility, assuming that $f$ is low enough for this effect to be unimportant. Let us now turn our attention to the effect of cosmological perturbations. As we shall see, the parameter $\alpha$ can also be sensitive to these.

3 Cosmological perturbations

In this Section we consider the effect of cosmological perturbations on the eccentricity of PBH binaries. In contrast with previous work, we include the effect of density perturbations in radiation, which can be dominant for binaries with $\lambda\ll 1$, decoupling from the Hubble flow deep in the radiation era.

Let us assume that primordial density perturbations are adiabatic and Gaussian. These are completely characterized by the primordial power spectrum of a single scalar variable, such as the temporal component of the metric perturbation in the longitudinal gauge [40], often denoted by $\Phi$. This variable plays the role of the Newtonian potential in the non-relativistic limit. In Fourier space the gravitational potential is expressed as

[TABLE]

For perturbations with $k\eta_{eq}\gg 1$, entering the horizon well before equality($\eta\ll\eta_{eq}$), radiation dominates over dark matter and baryons. Neglecting the decaying mode on supercurvature scales, the time dependence of such modes for $\eta\ll\eta_{eq}$ is then given by [40]

[TABLE]

where $\kappa=\Omega_{DM}/\Omega_{M}\approx 0.84$ is the fraction of non-relativistic matter in the form of dark matter, and the time dependence is given in terms of

[TABLE]

Here, ${\bf C}\approx 0.577$. The initial amplitudes $\Phi_{\vec{k}}^{0}$ of the gravitational potential on superhorizon scales are Gaussian distributed, with variance given by

[TABLE]

The functions $G$ and $H$ represent the contribution of radiation and matter density perturbations, respectively. Which one dominates the torque will depend on the time $\eta_{b}$ when binaries decouple from the Hubble flow, which in turn is related to the co-moving binary size [see Eq. (2.7)]. At early times, matter is subdominant, and the contribution of matter perturbations to $\Phi$ is suppressed by the scale factor $s=(\eta/\eta_{eq})$ in front of $H$ in Eq. (3.3). The factor $\Theta(x-1)$ in $H$ should not be taken too literally, it is just meant to indicate that the expression is only valid after the modes cross the horizon, $k\eta_{eq}\gg 1$, and the logarithmic growth begins.

The angular momentum per unit reduced mass of the binary, $\vec{\ell}$, can be written as a time integral of the tidal torque exerted by the gravitational potential. Denoting by $\vec{x}_{1}(\eta)$ amd $\vec{x}_{2}(\eta)$ the co-moving positions of the two members of the binary, we have:

[TABLE]

where $\vec{x}(\eta)=\vec{x}_{2}-\vec{x}_{1}$ is the relative co-moving coordinate, which is time dependent from the time $\eta_{b}$ when the binary decouples from the Hubble flow [see Eq. (2.5)].

In order to calculate the variance of the angular momentum, we will work at lowest order in the gravitational potential, so that inside the integrand in Eq. (3.5) we can use the unperturbed head-on trajectory, which we shall take along the $z$ axis, with relative coordinate:

[TABLE]

Assuming that the center of mass is at the origin of coordinates, the positions of the two PBHs are given by $\vec{x}_{2}=(m_{1}/M)\vec{x}$ and $\vec{x}_{1}=-(m_{2}/M)\vec{x}$. Using (3.4) and (3.5) we have

[TABLE]

where the indices $i,j,m,n$ can only take values $x$ or $y$ and

[TABLE]

Introducing spherical coordinates in momentum space, $k_{x}=k\sin\theta\cos\phi$, $k_{y}=k\sin\theta\sin\phi$, $k_{z}=k\cos\theta$, and integrating over $\phi$, we have

[TABLE]

where we have used $\langle\ell^{2}\rangle=\langle\ell^{i}\ell^{j}\rangle\delta^{ij}=2\langle\ell^{x}\ell^{x}\rangle$ and we have introduced the change of variable $w=\cos\theta$.

To estimate the integral $F$ we first note that the radiation $G$, and matter $H$ terms make their contribution at very different times. Consider a binary with initial separation $x_{0}$ that decouples from the Hubble flow at the conformal time $\eta_{b}$. From (2.8) and (2.10), these two scales are widely separated, and parametrically related by

[TABLE]

In general, all perturbations make their contribution to the torque at times $\eta\lesssim\eta_{b}$. After that, the binary starts oscillating, its co-moving size shrinks, and tidal gradients decay in inverse proportion to the scale factor. Perturbations in radiation start oscillating once they enter the horizon, and make most of their contribution to the torque at $\eta\sim k^{-1}\lesssim\eta_{b}$, while matter perturbations make their contribution near the time $\eta\sim\eta_{b}$, regardless of their wavelength.

3.1 Perturbations in radiation

Let us start by considering small binaries, for which $\lambda\ln(k_{0}\eta_{b})\ll 1$. These form deep in the radiation dominated era, at the time when $s=s_{b}(\lambda)\sim\lambda\ll 1/\ln(k_{0}\eta)$. In this case, perturbations in radiation dominate over matter perturbations. From Eq. (3.10), and assuming masses in the stellar range, the logarithm is of order 10, and so this condition requires $\lambda\ll 0.1$. In this regime, we can neglect $s(\eta)H(k\eta/3)$ relative to $G(k\eta/3)$ for $\eta\lesssim\eta_{b}$, and we have

[TABLE]

Due to the oscillating nature of $G$ for $\eta\gg k^{-1}$, the integral is dominated by early times $\eta\ll k^{-1}\lesssim\eta_{b}$, where $x=x_{0}\chi(\eta)/s$ is approximately constant, $x\approx x_{0}$, and can be taken out of the integral (3.8). Noting that

[TABLE]

it is clear that the contribution of modes outside the horizon at the time of binary formation, $k\eta_{b}\ll 1$, is suppressed. For modes with $k\eta_{b}\gg 1$ we can approximate

[TABLE]

Substituting in (3.9), and using

[TABLE]

where the function $G$ is defined in (3.3), we have888In determining the numerical coefficient in front of (3.16), we have used, from (2.2) and (2.8),

$GMa=G\beta\lambda\bar{m}x=G\beta\bar{m}{x_{0}^{4}\over f\bar{x}^{3}}=\beta x_{0}^{4}{2\pi G\over 3}\rho_{eq}={\beta\over 4}x_{0}^{4}H^{2}_{eq}\approx{x_{0}^{4}H^{2}_{eq}\over 20}.$

(3.15)

Also, we have used $H_{eq}^{2}\eta_{eq}^{2}\approx 1$. Note that we are using the convention where the scale factor is equal to unity at the time of equality.

[TABLE]

The subindex in $j_{cp}$ indicates that this is due to cosmological perturbations, as opposed to the neighboring black holes which we considered in the previous Section. Here we have introduced the standard expression for the primordial power spectrum

[TABLE]

In standard slow-roll inflationary scenarios, $P_{\Phi}(k)$ is nearly independent of $k$. In that case, the two factors in square brackets in Eq. (3.16) play the role of the infrared and ultraviolet cut-off which regulate the logarithmic behaviour of the integral. Note that

[TABLE]

while

[TABLE]

Therefore, for a nearly flat power spectrum, the integral will be dominated by the range $\eta_{b}^{-1}\lesssim k\lesssim k_{0}=x_{0}^{-1}$,

[TABLE]

where

[TABLE]

In the last approximate equality, we have used (3.10).

It may seem counterintuitive that the dispersion in $j$ receives contributions from a wide range of scales, since at the time $\eta_{b}$ when the binary forms, the amplitude of the gravitational potential for modes within the horizon falls off as $k^{-2}$. However, the torque at wavelenths larger than $x_{0}$ depends on second derivatives of the potential, which brings in a factor of $k^{2}$. As a result, the contribution is independent of scale in the range we are considering. Physically, the effect takes place well before the binary decouples from the Hubble flow, at the time $\eta\lesssim k^{-1}$. Hence, it seems appropriate to refer to this as the contribution of the peculiar velocities of the PBHs to the orbital angular momentum at the time when the binary forms.

It is also worth noting that the prefactor in front of the logarithm is independent of the parameters characterizing the binary. Each decade in wavelength gives the same contribution to $\langle j^{2}\rangle$, and the dependence on parameters such as massess and semi-major axis, is only through the range of scales contributing to the logarithm. This is in contrast with the contribution from matter perturbations, which we now review.

3.2 Adding matter perturbations

Matter perturbations can be included along similar lines. One difference is that their effect on the binary occurs near the time $\sim\eta_{b}$ when the binary starts oscillating, and we cannot ignore the time dependence of the separation $x$ in the integral (3.8) which gives $F(k,x_{0},w)$. For $kx_{0}\gg 1$, this has a dependence on $m_{1}$ and $m_{2}$ which, unlike the radiation case, is hard to disentangle in general.

For $kx_{0}\ll 1$, we may use the approximation

[TABLE]

in the integral (3.8). With this approximation, the dependence on masses disappears and we obtain the total contribution of radiation and matter perturbations as 999In order to obtain the numerical coefficient in front of the matter contribution, we have used [10, 20] $\int_{0}^{1}ds(\chi^{2}/s^{2})\approx 0.3\lambda(2\bar{m}/M)$ to do the intergral of the second term in (3.8). Also, we have used $\kappa\approx 0.84$ for the ratio of dark matter density to the total non-relativistic matter density, and we have ignored the slow logarithmic dependence in $\eta$. Since the integral is dominated by $\eta\sim\eta_{b}$, we have used the value $\eta=\eta_{b}$ inside the logarithm in the mode function $H$.

[TABLE]

where

[TABLE]

Here, we have used (3.21), neglecting the small logarithmic dependence in $\lambda$ and $m$.

Substituting (3.23) in (3.9) and performing the $w$ integration, we immediately find the total variance of the orbital parameter due to long wavelength cosmological perturbations:

[TABLE]

The first term in the square brackets in (3.23) corresponds to radiation, while the second one, accompanied by the factor of $\lambda$, corresponds to matter perturbations. The latter become subdominant for sufficiently small $\lambda\lesssim 0.27$. In view of our earlier discussion in Subsection 2.3, for $f\gg 3\cdot 10^{-3}$ the rates are dominated by small binaries, with $\lambda\lesssim\lambda_{*}\ll 1$. Hence, it appears that perturbations in radiation may be as relevant for observations as the matter perturbations which have been considered in earlier analysis.

Let us now turn to a discussion of the effect of such perturbations on the merger rates.

4 Effect of perturbations on the merger rates

The size of binaries at the time of formation (shaded in gray in Fig. 3) is at intermediate co-moving scales which are much smaller than those probed by CMB temperature anisotropies or large scale structure. Indeed, from Eqs. (2.8) and (2.10), we have

[TABLE]

Taking into account that $\lambda$ is in the range (2.23), for stellar mass black holes $m\sim 1-100$ we have $k_{0}\sim(10^{5}-10^{7})k_{eq}$, which corresponds to the present co-moving scale in the range101010Throughout this paper, we adopt the convention that the scale factor $s$ is equal to $1$ at the time of equality. Thus, to avoid confusion, we will refer to the present day wave number by $\bar{k}=z_{eq}k$. Note that relations such as Eq.(4.1) between $k_{0}$ and $k_{eq}$ are valid in both conventions, since the factor of $z_{eq}$ applies to both sides of the equation. $\bar{k}\sim(10^{3}-10^{5})Mpc^{-1}$.

The angular momentum of PBH binaries may be affected by the power spectrum $P_{\Phi}$ of cosmological perturbations over a very wide range of scales. As illustrated in Fig. 3, such power spectrum is poorly constrained on scales smaller than $3Mpc^{-1}$, and here we would like to explore the consequences this uncharted territory might have on the merger rate of binaries.111111In Fig. 3, we are ignoring bounds which are related to the abundance of ultra-compact mini-halos. Such constraints depend on the nature of dark matter, and could be absent in certain models (e.g. if dark matter belongs to a hidden sector). For a recent discussion, see [41] and references therein. For this purpose, let us consider three distinct behaviours which may capture the generic effect of an enhanced power spectrum at small scales. These are labeled case A, B and C in Fig. 3. Let us consider them in turn.

4.1 Case A: Nearly scale invariant cosmological perturbations

Consider a nearly scale invariant power spectrum of the form

[TABLE]

This is consistent with observations of the CMB and large scale structure on cosmological scales, down to $\bar{k}\sim 3Mpc^{-1}$, with $n_{s}\approx.97$, $A_{\Phi}\approx.97\cdot 10^{-9}$ and $\bar{k}_{*}\approx 0.05Mpc^{-1}$ [43]. A minimal assumption we can make, consistent with standard slow roll inflationary models, is that the nearly flat spectrum can be extrapolated down to the co-moving size of binaries. Since the tilt is rather small, we may approximate $P_{\Phi}\approx const.$ over the range of interest.

Using (3.25) we then find

[TABLE]

Assuming the CMB normalization $P_{\Phi}\approx 10^{-9}$, we have $\sigma^{2}_{cp}\sim(2-9)10^{-6}$. Hence, for moderate values of $\lambda$, cosmological perturbations may easily provide an angular momentum comparable to $j_{0}$, given in (2.14), necessary for binaries to have a life-time comparable to the age of the universe. Still, neighboring PBH contribute to $j$ with variance given by (2.54) [20],

[TABLE]

where we are excluding binaries which may be disrupted by a PBH at distances smaller than $Y_{min}$. Thus, at low $f$ the effect of neighboring PBHs on the angular momentum is negligible compared to that of cosmological perturbations, while the latter effect can become important only at higher $f$. This is illustrated in Fig. 4.

To better understand the relative importance of the two effects, it is illustrative to look at the differential merger rate in the integrand in Eq. (2.21):

[TABLE]

where $\sigma^{2}=\sigma^{2}_{nb}+\sigma^{2}_{cp}.$

The function $W(\lambda)$ for Case A is plotted in thick lines in the left pannel of Fig. 5, for different values of $f$, and $m=30$. At low $f\lesssim 10^{-2}$, the effect of neighboring PBHs is negligible. This corresponds to the values below the knee in Fig. 4. In this case Fig. 5 shows that for the low amplitude cosmological plateau with $P_{\Phi}\approx 10^{-9}$ the dominant contribution to $W$ is at $\lambda\sim 1$. This value corresponds to a time of binary formation $\eta_{b}\sim\eta_{eq}$. Even at this relatively late time the effect of relativistic matter perturbations has a noticeable effect on the position of the peaks. For comparison, we plot in dotted lines the case where perturbations in radiation are ignored. At higher $f\sim 10^{-1}$ we see from Fig. 5 that the rate is dominated by lower $\lambda\sim 0.1$. For such values of $\lambda$, binaries decouple earlier in time, when perturbations in radiation would be more important than those in non-relativistic matter. But in fact, both are negligible, since for such low values of $f$, the variance in $j$ is dominated instead by the torque from neighboring PBHs. The situation is different when we consider an enhanced power spectrum (Case B), which is plotted in the right pannel of Fig. 5.

Before moving on to a discussion of enhanced power spectra, let us consider the parameter $\alpha$ in the presence of the standard scale invariant spectrum of cosmological perturbations (4.2). Note that the variance (4.4) dominates at high $f$, whereas the last term in (4.3) dominates at low $f$. Both terms have exactly the same dependence on $\lambda$ and on the masses. They only differ in the explicit dependence of the variance on $Y_{min}$, which is absent at low $f$. Hence, by the same argument which lead us to Eq. (2.58), it is straightforward to conclude that

[TABLE]

where the lower and upper bounds correspond to the limiting behaviours for low and high $f$ respectively.

4.2 Case B: Enhancement of the power spectrum at intermediate scales.

Bounds on the amplitude of the power spectrum above $\bar{k}\sim 10Mpc^{-1}$ are rather loose. The absence of CMB spectral distortions gives an upper bound at intermediate scales in the range $\bar{k}\sim 10-10^{5}Mpc^{-1}$ [46, 41]

[TABLE]

several orders of magnitude higher than the nearly scale invariant plateau observed on cosmological scales.

Hence, let us consider a generic power spectrum of the form

[TABLE]

where for convenience we choose the pivot scale to be $k_{B}=5\cdot 10^{5}k_{eq}$, corresponding to the present co-moving scale of $5\cdot 10^{3}Mpc^{-1}$. From Eq. (4.1), this is by order of magnitude the co-moving size of the binaries, and we see from Fig. 3 that the power spectrum can be rather high at those scales, allowing for an enhanced amplitude up to $B_{\Phi}\lesssim 10^{-5}$. Let us now discuss the impact of this enhancement as a function of $p$. Here we focus on $p<2$, leaving the discussion of steeper power spectra for the next Subsection.

The case $p=0$ is special, since the variance $\sigma^{2}_{cp}$ depends logarithmically on the width of the plateau in the power spectrum. For illustration, let us consider a scale invariant spectrum of amplitude $B_{\Phi}$ between the scale $k_{min}=10^{3}k_{eq}$ and $k_{max}\gtrsim 10^{7}k_{eq}$. Then, using (3.25), we find121212In this estimate, we neglect a small subleading dependence in $\lambda$ and $m$, assuming $|\ln\left[(\lambda/0.3)(m/30)\right]|\ll 16$.

[TABLE]

Comparing with (4.4), it is clear that for $f\lesssim 0.1$ cosmological perturbations will dominate over the effect of neighboring PBH for amplitudes $B_{\phi}\gtrsim 10^{-7}$, which are well within the range allowed by observational constraints. This is illustrated in the right pannel in Fig. 5, where we see that the rate is dominated by $\lambda\lesssim 0.1$ Hence, we may approximate

[TABLE]

Using (2.14) we then have

[TABLE]

where

[TABLE]

For sufficiently high $B_{\phi}\gtrsim 10^{-7}$ we have $\lambda_{B}^{(p=0)}\ll 1$, and we can ignore the upper limit of integration in (2.21). Then we have:

[TABLE]

and using (2.38) and (4.10) we have

[TABLE]

We note that the same parametric dependence is obtained if we simply assume that all binaries are born equal, with the same initial dimensionless angular momentum 131313Indeed, if all binaries have the same $j=j_{*}$, then the rate is determined by those binaries whose semi-major axis is such that $j_{0}(\lambda)=j_{*}$. Assuming, for simplicity, a monochromatic PBH mass spectrum, this determines $\lambda=\lambda_{0}=1.3\cdot 10^{-4}j_{*}^{-21/16}m^{5/16}.$ Hence, we have

$\Gamma(t_{0};j_{*})={f\rho_{m}\over 4\bar{m}}\int\delta(t_{0}-t[a;j_{*}])e^{-X}dX\approx 1.8\cdot 10^{4}f^{2}{m^{-11/16}\over j_{*}^{21/16}}Gpc^{-3}yr^{-1}e^{-Y_{min}},$

(4.15)

where we used $t\propto a^{4}\propto\lambda^{16/3}$ and $X=\lambda f$ in order to do the integration over $X$ by using the Dirac delta function. This matches Eq. (4.14) for $j_{*}\approx 0.33\ \sigma$. $j_{*}\sim\sigma$. From (4.14) we have $\Gamma_{M}\propto M^{15/16}e^{-Y_{min}}$, and therefore

[TABLE]

which is significantly lower than the values in the range (4.6), corresponding to the cosmological plateau. The difference comes from the fact that for the enhanced spectrum, perturbations in radiation dominate over matter perturbations, and this changes the mass dependence of merger rates. For $P_{\Phi}\approx 10^{-5}$, and for a monochromatic spectrum, we have

[TABLE]

This turns out to be an excellent approximation to the numerical result which is plotted in Fig. 6. In the mass range $30-100$ solar masses, the LIGO bound on the merger rate requires $f\sim 5-9\%$, where we have used $Y_{min}\approx 2$. Hence, the bound on $f$ is significantly relaxed compared to the case without cosmological perturbations.

For $0<p<2$, the integral in (3.25) is dominated by $k\approx k_{0}$. As we have just seen for the case of $p=0$, when the amplitude of cosmological perturbations is enhanced, it is the perturbations in radiation that dominate, and in this case we know that the contribution of modes with $k\gtrsim k_{0}$ can also be important. This is given by Eq. (3.16). Defining

[TABLE]

and using (4.8) in (3.16) we have

[TABLE]

For the purpose of keeping track of the relative importance of radiation versus matter perturbations, here we have included the effect of matter perturbations up to the time $\eta=\eta_{b}$, which corresponds to the second term in square brackets in Eq. (4.19). Matter perturbations can contribute to the angular momentum for some time after $\eta_{b}$, when the two PBHs forming the binary start falling towards each other. Here we have not included this contribution, which in principle can have a more complicated mass dependence. However, by order of magnitude, this effect is comparable to the contribution up to the time $\eta_{b}$. This characterization suffices in order to check whether perturbations in radiation are dominant or not. The numerical coefficient $\psi(p)$ is represented in Fig. 7, where we can see that $\psi>0.44$ for all values of $p$, and that it diverges near both ends of the interval $0<p<2$. Such divergences correspond to the infrared and ultraviolet logarithmic divergences in the integral (3.16), respectively. In practice, these are cut-off by the finite range where the power spectrum has the specified behaviour, as we discussed for the case $p=0$. Here, we shall simply consider a generic $p\sim 1$, which is not too close to [math] or $2$, so that $\psi\sim 1$. In this case, we obtain

[TABLE]

In order to estimate the merger rates, let us first start by neglecting matter perturbations in (4.20). Using (2.14), we have

[TABLE]

with

[TABLE]

where we have set $k_{B}=5\cdot 10^{5}k_{eq}$. For $B_{\Phi}\gtrsim 10^{6}$ and $m\lesssim 100$, we find that $\lambda_{B}^{(p)}\ll 1$. This is consistent with the assumption that perturbations in radiation dominate over perturbations in matter, since the peak of the differential rate $W(\lambda)$, given by (4.5), occurs at $\lambda\sim\lambda_{B}^{(p)}\ll 1$.

For $0\lesssim p\lesssim 1$, the integral in (2.21) quickly converges for $\lambda\gg\lambda_{B}^{(p)}$ and the upper limit of integration becomes irrelevant. In this case we have

[TABLE]

As $p$ approaches the value $11/7$, the ratio given in (4.21) approaches the behaviour $(j_{0}^{2}/\sigma^{2})\propto\lambda^{-1}$, and the integral (2.21) becomes logarithmically divergent at large $\lambda$. This is regulated by the finite range of $\lambda$. Indeed, for $1\lesssim p\lesssim 2$, the integral (2.21) is dominated by the interval $\lambda_{B}^{(p)}\lesssim\lambda<\lambda_{m}$. Here, $\lambda_{m}\sim 0.3$ is the value for which matter perturbations start becoming important. At that point, the behaviour of the integrand switches to $(j_{0}^{2}/\sigma^{2})\propto\lambda^{-3+{7p-11\over 21}}$, which rapidly converges. Hence, in that case we can approximate

[TABLE]

Let us now consider the ”universality” coefficient $\alpha$, defined in Eq. (2.48). From (4.23) or (4.24), we immediately obtain

[TABLE]

where in the case (4.24) we have ignored the subleading logarithmic dependence on $M$ (we will comment on this subleading correction below).

Factoring out the $f^{2}$ and $Y_{min}$ dependence in $\Gamma_{M}$, as

[TABLE]

we have represented the approximations (4.23) and (4.24) in Fig. 8 as dashed lines (for the case $m=1$). For comparison, we numerically calculate the rate which is obtained when we include, in addition to the effect of perturbations in radiation, the torque exerted by non-relativistic matter perturbations up to the time $\eta\approx\eta_{b}$ when the binary decouples from the Hubble flow. The latter torque corresponds to the second term in square brackets in Eq. (4.19). The numerical result for $m=1$ and $m=30$ is depicted in thick lines in Fig. 8. Note that the approximations (4.23) and (4.24) can be imprecise, by up to a factor of 2 or so in the region $p\sim 1$, due to the fact that we have neglected the effect of matter perturbations, which can affect the merger rates by a sizable fraction even if they are subdominant.

A much better analytic approximation to the numerical result can be obtained by keeping the effect of matter perturbations in $\sigma^{2}_{cp}$. Using $\lambda_{B}^{(p)}\ll 1$, we may approximate

[TABLE]

Here, we have neglected the integrand for $\lambda\lesssim\lambda_{B}^{(p)}$, since this is a small interval where $W$ is suppressed relative to its peak value at $\lambda\sim\lambda_{B}^{(p)}$, and we have also neglected the exponential factor in $W$ for $\lambda\gtrsim\lambda_{B}^{(p)}$, since the exponent is small in this range. The integral in the right hand side can be calculated in terms of the incomplete Euler’s $\beta$ function. Expanding this function for small $\lambda_{B}^{(p)}$ we find

[TABLE]

The analytic approximation (4.26) with (4.28) is plotted in Fig 8 in dotted lines for $m=1$ and $m=30$. We find that it reproduces the full numerical result plotted in thick lines, with very good accuracy in the full range of $p$. The term in round brackets contains additional dependence on $m_{1}$ and $m_{2}$ which is not present in (4.23), and this will contribute a correction $\Delta\alpha$ to the expression (4.25), given by

[TABLE]

This correction turns out to be rather small, and typically $\Delta\alpha\lesssim 0.05-0.1$, unless the mass ratio is hierarchical. The reason can be understood as follows. The mass dependence comes from the term in square brackets in (4.29),

[TABLE]

The exponents $\beta$ and $\gamma$ happen to be very small $\beta,\gamma\lesssim 0.1$ in the range of interest. Then, unless one of the two masses is hierarchically smaller than the total mass $M$, we have $\Delta\alpha\propto\beta+O(\beta^{2},\beta\gamma,\gamma^{2})\sim\beta$. The correction $\Delta\alpha$ is plotted in Fig. 10 for $p=1.8$, $\bar{m}=10M_{\odot}$, and a range of values of $m_{1}$ and $m_{2}$. Note that $|\Delta\alpha|\lesssim 0.08$, unless we consider the region where $m_{1}\ll m_{2}$. We have checked that for the case of similar masses $m\lesssim 30$, we have $|\Delta\alpha|\lesssim 0.08$ in the full range $0.2\lesssim p\lesssim 1.8$

We conclude that, in the range $0\lesssim p\lesssim 2$ the spectral index $p$ determines the parameter $\alpha$, which may range from $15/16\approx.93$ up to $5/3\approx 1.66$. This is represented in Fig. 9. A measurement of $\alpha$ may therefore provide valuable information on the primordial perturbation power spectrum at intermediate scales. Note that the value of $\alpha$ is independent of the amplitude of perturbations, as long as this amplitude is sufficiently large 141414At low amplitudes, $B\lesssim 10^{-7}$ the torque from matter perturbations becomes as important as that from radiation perturbations. In general, in this case, the dependence of $\Gamma_{M}$ on $M$ does not necessarily factor out as a power of $M$. Still, if we make the assumption that the torque from matter perturbations is dominated by the contribution from times up to $\eta_{b}$, then using (4.19) and neglecting the first term in square brackets (which is due to radiation perturbations), we find $\alpha=72/(74-7p)$, which is larger than $36/37\approx 0.97$ and smaller than $6/5=1.2$. This is a narrower range than allowed by perturbations in radiation. $B_{\Phi}\gtrsim 10^{-6}$.

It should also be pointed out that the bound on $f$ from the observed merger rates can be considerably relaxed by the enhanced power spectrum at intermediate scales. For instance, if we take $m=30$ and $B_{\Phi}=10^{-5}$, then from Fig. 8 we have $\tilde{\Gamma}(p)\lesssim 10^{5}$. Consequently, for a monochromatic PBH mass function, the bound

[TABLE]

with $Y_{min}\approx 2$ is satisfied for $f\lesssim 0.1$. 151515For smaller values of the power spectrum, in the range $10^{-7}<B_{\Phi}<10^{-5}$, the merger rate scales approximately in proportion to $\lambda_{B}^{(p)}\propto B_{\Phi}^{-21/(32-7p)}$. This translates into the bound $f\lesssim 0.1(B_{\phi}/10^{-5})^{21/(64-14p)}$.

4.3 Case C: Steep spectrum at very small scales.

Consider a power spectrum with spectral index $p\gtrsim 2$ up to some high scale $k_{C}$. In this case, the contribution of perturbations to the angular momentum of binaries will be dominated by $k\sim k_{C}$, rather than $k\sim k_{0}$, and it is convenient to parametrize by using $k_{C}$ as the pivot scale:

[TABLE]

From the constraints in Fig. (3), and assuming $k_{C}\geq 10^{7}k_{eq}$, we require

[TABLE]

The case with $p<2$ was discussed in the previous subsection, with the correspondence

[TABLE]

between the prefactors in (4.8) and (4.32). Here we concentrate in the case $p>2$. The specific form of $P_{\Phi}(k)$ for $k\gtrsim k_{C}$ will not be important, as long as it grows slower than $k^{2}$, or that it decays for $k>k_{C}$.

From (3.16), the contribution of radiation perturbations to the variance of $j$ is dominated by wavelengths shorter than the binary size $x_{0}$, and is given by

[TABLE]

In the last step we have neglected the contribution from the lower limit of integration, assuming that $p$ is not too close to $2$. As we discussed in the previous subsection, we can easily incorporate the torque of matter perturbations up to the time $\eta=\eta_{b}$, which gives a combined total of

[TABLE]

Here we have used $k_{0}/k_{B}=(m\lambda)^{-1/3}$, with $k_{B}=5\cdot 10^{5}k_{eq}$.

Given that $C_{\phi}$ is bounded above by (4.33), the factor $(k_{B}/k_{C})^{2}$ will considerably suppress the effect of cosmological perturbations unless $k_{C}$ is not too far from the intermediate scale $k_{B}$. Eq. (4.36) can be compared with de contribution of neighboring PBH to the torque, given in (4.4)

[TABLE]

where in the last step we use the fiducial values $K\sim 1$, $M\sim 2\bar{m}$ and $Y_{min}\approx 2$. For $p$ not too close to 2, and $m\sim 30$, the effect of cosmological perturbations can only be dominant for a narrow range of $k_{C}$:161616In the case of cosmological perturbations, the merger rate for $p>2$ is dominated by the values of $\lambda$ where the behaviour of the variance (4.36) changes from $\sigma^{2}_{cp}\propto\lambda^{-2/3}$ to $\sigma^{2}_{cp}\propto\lambda^{4/3}$. This happens at $\lambda\sim 1/3$, and we use this in the relation $\sigma_{cp}\gtrsim\sigma_{nb}$ in order to estimate the upper limit of the range (4.38).

[TABLE]

Hence, the existence of a window for $k_{C}$ where $\sigma_{nb}\ll\sigma_{cp}$ requires $f\ll 1$.

Following similar steps as in Subsection 4.2, we find that the expected merger rate of binaries will be given by (4.26), with

[TABLE]

Here,

[TABLE]

This leads to a value of the universality coefficient which is given by (4.25) with $p=2$. That is

[TABLE]

Eq. (4.41) corresponds to the thick horizontal line in the range $p>2$ in Fig. 9. Note, however, that the approximation (4.39) requires $\lambda_{C}\ll 1$ (since it is obtained by Taylor expanding an incomplete Euler beta function in small $\lambda_{C}$). The merger rates for $k_{C}=10^{7}k_{eq}$ are plotted in Fig. 11, where it can be seen that the approximation (4.39) plotted as dashed lines, agrees very well with the numerical result plotted in thick lines, for $m\lesssim 30$. For higher values of $k_{C}$ or higher values of $m$, we are outside the regime $\lambda_{C}\ll 1$, and the approximation is not valid. Nonetheless, the behaviour of the merger rates from the numerical result is similar for all masses up to $m\lesssim 100$. In particular, we find that $f$ should be less than a few percent in the whole range $m=1-100$. This is in contrast with the results for Case B ($p<2$), where values of $f$ as large as $10\%$ can be consistent with the observed merger rates. In spite of the high amplitude of the power spectrum, the effect in Case C is not as significant as in case B, since perturbations are now on scales smaller than the binary size.

4.4 Effect of a peak at the PBH scale

The mechanism of PBH formation from adiabatic perturbations relies on a prominent enhancement in $P_{\phi}$, with a maximum at some co-moving scale $k_{PBH}$ with amplitude

[TABLE]

This is enough to produce PBHs in significant abundance, in regions where high peaks in the Gaussian random field of cosmological perturbations reach non-linear values above a certain threshold. As discussed in Refs. [36, 37], PBH formation does not occur when $k_{PBH}$ crosses the horizon, but slightly later, when the scale

[TABLE]

crosses the horizon171717The scale $r_{m}$ corresponds to the maximum of the so-called compaction function ${\cal C}(r)$, which characterizes the averaged overdensity as a function of the distance to its center. The approximate factor of 3 in (4.43) depends somewhat on the shape of the overdensity profile, which for amplitudes well above the standard deviation is in turn determined by the form of $P_{\Phi}(k)$. For instance, for a monochromatic power spectrum, $P_{\Phi}(k)\propto\delta(k-k_{PBH})$, we have $r_{m}\approx 2.7\ k_{PBH}^{-1}$. On the other hand, the strong overdensity causes a non-linear distortion of the spatial metric which affects the relation between co-moving and physical distance. Generically, the combination of these non-linear effects gives a mass of the black hole which is $\sim 10$ times bigger than the mass within the unperturbed horizon at the time when $k_{PBH}$ crosses it. Hence the factor of 3 in (4.44).. The mass of the black hole is typically equal to the mass within the horizon at that time, so that $2GM=H_{r_{m}}^{-1}$, and we have

[TABLE]

Comparing with Eq. (4.1) this corresponds to lengthscales which are shorter than the co-moving size of the binary, by three orders of magnitude or so, corresponding to the present co-moving wave-number of order

[TABLE]

On a logarithmic scale, this is not too far from the intermediate scales discussed in Subsection 4.2, and hence it is natural to ask about the consequences of this peak on the merger rates, in the light of our earlier results.

For illustration, in Fig. 3 we include an example extracted from [42], where the enhancement of the power spectrum is obtained from a period where the inflaton goes from slow roll to ultra-slow roll and back to slow roll, through a discrete set of values of the second slow roll parameter $\eta$. This model may be somewhat artificial, but it will be useful for illustrative purposes since it includes several features which may be present in more general cases. The curve includes a steep growth with $P_{\Phi}\propto k^{p}$ from the cosmological plateau up to a very high value $P_{\Phi}\sim 10^{-3}$, at $k\sim k_{C}$, as in Case C. This is followed by a less steep part up to some maximum value $k_{PBH}$, corresponding to the scale of PBH formation. Finally, the spectrum decays for $k>k_{PBH}$.

It has been argued that $p=4$ corresponds to the theoretical upper bound for the spectral index for a wide class of inflationary one field models [42]. Values in the range $3\lesssim p\lesssim 4$ are easily obtained in the case where the enhancement is due to a short period of ultra-slow roll and constant roll inflation. The upper bound can be saturated in relatively simple models where the constant roll phase has a sufficiently negative value of the slow roll parameter $\eta$ (see e.g. [45]). For $p\gtrsim 2$, the regime where cosmological perturbations dominate over the effect of neighboring PBHs is given by Eq. (4.38), which requires

[TABLE]

Here we have used (4.44). For a sizable value of the fraction of dark matter in the form of PBHs, say $f\gtrsim 10^{-2}$, cosmological perturbations will only play a role if the peak of the power spectrum is rather broad, with $k_{PBH}\gtrsim 10\ k_{C}$. Conversely, for a narrow peak with $k_{PBH}\sim k_{C}$ and sizable $f\gtrsim 10^{-2}$, the eccentricity of binaries is determined by the tidal forces from neighboring PBHs. Broad peaks in the power spectrum are not necessarily generic, but as shown by the example in Fig. 3, they can in principle be obtained in phenomenological models.

The merger rate of binaries is plotted in the left pannel of Fig. 12 for different values of the mass, for the case of a broad peak which raises steeply with $p=4$ up to the scale $k_{C}\approx 10^{8}k_{eq}$, and then proceeds with moderate slope to the maximum at $k_{PBH}=10k_{C}\approx 10^{9}k_{eq}$. In this case the analytic approximation (4.39) does not apply, since $\lambda_{C}\gtrsim 1$ in the relevant range of $f$. The effect of cosmological perturbations is only significant for low masses $m\lesssim 1$.

For the case of a narrow peak with $p=4$, which raises from low values with a steep spectral index to all the way to the PBH scale $k_{PBH}=k_{C}\approx 10^{9}k_{eq}$, the merger rate is plotted in the right pannel of Fig. 12, where we see that the effect of the narrow peak is completely insignificant.

For less steep spectral index $0<p\lesssim 2$, the growing part of the spectrum corresponds to Case B. The condition that cosmological perturbations provide the dominant source of torque on binaries can be derived along similar lines, and is given by

[TABLE]

This is much less restrictive on $f$ than (4.46). Indeed, for $p\lesssim 1$, the condition does not significantly restrict the fraction $f$, even for the case where $k_{C}\approx k_{PBH}$. 181818Note, however, that in order to interpolate between the cosmological value $P_{\Phi}\sim 10^{-9}$ and the peak value $P_{\Phi}\sim 10^{-3}$ in the span of wavelengths which goes from $10Mpc^{-1}$ to $10^{7}Mpc^{-1}$, we need $p\gtrsim 1$. Incidentally, for $p=1$ the value of the “universality” parameter $\alpha$ which is obtained when cosmological perturbations provide the dominant torque (see Fig. 9) is quite similar to the value $1.26$ which corresponds to a dominant torque from neighboring PBHs (at a distance larger than the infall radius $Y_{min}$), so for this particular value of $p$ it seems harder to distinguish one mechanism from the other. The effect of such a feature in the power spectrum is described in Subsection 4.2.

5 Summary and conclusions

We have studied the effect of cosmological perturbations on the merger rate of PBHs in the stellar mass range. The effect can be quite significant depending on the amplitude of the power spectrum at different scales.

For a scale invariant spectrum with amplitude $P_{\Phi}\approx 10^{-9}$ (Case A), matter perturbations have a dominant effect on the eccentricity of binaries for $f\lesssim 10^{-2}$, while for larger $f$ the distribution of $j$ is actually dominated by tidal torques from PBHs in the vicinity. In this case, the merger rates would be greater than the current LIGO/Virgo bounds unless the fraction of dark matter in the form of PBH is rather small, of the order $f\lesssim 10^{-2}$ for PBHs with $m\sim 30$.

On the other hand, $P_{\Phi}$ could be much larger at scales $10-10^{5}Mpc^{-1}$ (Case B). For instance, we find that a nearly flat spectrum with amplitude $P_{\Phi}\gtrsim 10^{-7}$ within such scales leads to a dimensionless angular momentum with mean squared value

[TABLE]

which is mostly due to perturbations in the radiation fluid. For a nearly flat $P_{\Phi}$, the variance of $j$ is almost independent of the mass of the PBH and the size of the binaries. There is only a mild subleading logarithmic dependence on such parameters, which accounts for the range of scales contributing to the effect, from the co-moving size of the binary to the co-moving size of the horizon at the time when the binary forms. The situation is different when we have a tilted spectrum, of the form $P_{\Phi}\propto k^{p}$, ($p\lesssim 2$), in which case the variance of $j$ depends on binary size and masses. In this case the effect is determined by the amplitude of $P_{\Phi}$ near the co-moving scale of the binaries $k\sim k_{0}$, which is in the intermediate range $10^{3}-10^{5}Mpc^{-1}$. The bound from spectral distorsions in the CMB caused by dissipation of acoustic modes requires $P_{\Phi}\lesssim 10^{-5}$ at such scales [46]. In the situation where this bound is saturated, the torque can be large enough to significantly suppress the merger rate to a level consistent with LIGO/Virgo observations even for $f\sim 10\%$.

An even stronger enhancement of the power spectrum may occur at scales beyond $10^{5}Mpc^{-1}$. These are smaller than the co-moving size of the binaries, but they can contribute to the peculiar velocities of the PBHs and hence to the initial orbital angular momentum. For a steep power spectrum (Case C), with spectral index $p\gtrsim 2$ up to some short wavelength $k_{C}\gtrsim 10^{7}k_{eq}$, the contribution to $\langle j^{2}\rangle$ is dominated by the shortest scales $k\sim k_{C}$, rather than $k\sim k_{0}$. The observational upper bound $P_{\Phi}(k_{C})\lesssim 10^{-3}-10^{-2}$ is comparable to the amplitude which is necessary to produce PBH in significant abundances by the collapse of adiabatic perturbations. For generality, we may assume that beyond $k_{C}$, the power spectrum can still grow slightly (say, with a much lower spectral index $p^{\prime}\lesssim 2$), by a factor of a few, up to a maximum value at the scale which we may call $k_{PBH}$, corresponding to the scale which dominates PBH formation (see Fig. 3). In this case, and assuming $P_{\Phi}(k_{C})\sim 10^{-3}$, the effect on the orbital parameter of binaries is of order $\langle j^{2}\rangle\sim(k_{0}/k_{C})^{2}$, and becomes irrelevant unless $k_{C}\lesssim 10^{3}k_{0}$. Since $k_{0}<10^{5}k_{eq}$ while $k_{PBH}\sim 10^{9}k_{eq}$ for PBH in the stellar mass range, the effect of the steeply rising power spectrum on the eccentricity of binaries will only be important if this is crowned by a broad peak with a maximum at $k_{PBH}\gtrsim 10k_{C}$. The corresponding effect on the merger rates is somewhat intricate in general, but it is only significant at low masses. For $m\lesssim 100$, the bound on the fraction of dark matter can relaxed to $f\lesssim 2\cdot 10^{-2}$ due to this effect, for sufficiently low $k_{C}\sim 10^{7}k_{eq}$. Interestingly, in this regime the merger rate becomes almost independent of $m$ (see Fig. 11). This is in contrast with the standard situation where the angular momentum is supplied by a neighboring PBH, where we have $\Gamma\propto m^{-32/37}$.

We have also investigated the dependence of merger rates on the constituent masses, with particular attention to the universality coefficient $\alpha$ [1]. This is rather insensitive to the unknown initial PBH mass distribution function, and can be determined observationally with some precision (of order $15\%$ given a sufficiently large number of PBH merger events $\sim 10^{3}$ [1, 39]). It seems therefore of great empirical relevance for PBH scenarios. Our results for $\alpha$ are summarized in Fig. 9. In the case where cosmological perturbations are subdominant relative to the torque from neighboring PBH, it was argued in [1] that $\alpha\approx 1\pm 0.05$ (this is displayed as a horizontal shading in Fig. 9). We point out, however, that this narrow range shifts to $\alpha\approx 1.26$ once the effect of PBH infall onto binaries is taken into consideration, by excluding disturbed binaries from the merger count (after the infall these binaries are likely to have a much larger lifetime than the present age of the universe).191919The shift of $\alpha$ in the case where cosmological perturbations are subdominant is due to the factor $Y_{min}^{21/74}$ in Eq. (2.56). Note that the merger rates have an additional exponential dependence on the infall radius $Y_{min}$. Such dependence drops out from the universality coefficient $\alpha$ due to the linearity of $Y_{min}$ in the total mass $M$ [see Eq. (2.52)]. The assumption of linear behaviour seems very reasonable [13, 20], but may require further validation from numerical simulations, since any departures from it may have a significant effect on $\alpha$. This study is outside the scope of the present work and is left for further research. On the other hand, in the case where cosmological perturbations dominate the variance of $j$, the coefficient $\alpha$ ranges from $15/16$ to $5/3$ depending on the value of the spectral index $p$ [see Fig. 9]. In this case, the value of $\alpha$ is unaffected by infalls. We conclude that, as a matter of principle, an accurate measurement of merger rates of PBH might carry some information on the circumstances surrounding PBH binary formation, including the amplitude and spectral index of primordial perturbations on very short wavelengths or the effect of PBH infalls.

In turn, any information on the primordial power spectrum which may be obtained through a measument of $\alpha$ would constrain the underlying inflationary dynamics, from the scale of binaries down to the scale of PBH. This might complement other possible probes on the amplitude of $P_{\Phi}$ at small scales, such as upcoming searches for spectral distortions in the CMB, or pulsar timing array constraints [44, 41, 42].

There are several directions in which our study could be extended. It was recently pointed out in [19] that for $f\ll 1$ a dark matter dress may develop around PBH before they form binaries, and that this would affect the semi-major axis and the eccentricity of binaries after the first few oscilations which shake off most of the dark matter cloud. It was also shown that the effect on the merger rates is nonetheless small. In the present context, the effect may be even smaller, since in the presence of enhanced cosmological perturbations the merger rate is dominated by much smaller binaries forming deeper in the radiation era, when the halo around each PBH has had less time to accrete. Nonetheless, it might be interesting to be more quantitative about this effect, taking also into consideration the case of constituent masses which differ by a sizable factor, in order to assess the possible impact on $\alpha$. Finally, it was pointed out in [20] that for $f\gtrsim 0.1$ binaries can be disrupted by collision with compact N-body systems, which may substantially deplete the population of pristine binaries. Investigation of these issues seems to require further simulations, and remains an interesting direction of research.

Acknowledgements

We thank V. Atal, N. Bellomo, J. Bernal, C. Byrnes, J. Miralda, M. Sasaki and T. Tanaka for useful discussions. This work has been partically supported by FPA2016-76005-C2-2-P, MDM-2014-0369 of ICCUB (Unidad de Excelencia Maria de Maeztu), and AGAUR2017-SGR-754. N.T. is supported by an INPhINIT grant from la Caixa Foundation (ID 100010434) code LCF/BQ/IN17/11620034. This project has also received funding from the European Union s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 713673. We have used the publicly available WebPlotDigitizer [47] to digitise plots.