Initial clustering and the primordial black hole merger rate

If the primordial curvature perturbation followed a Gaussian distribution, primordial black holes (PBHs) will be Poisson distributed with no additional clustering. We consider local non-Gaussianity and its impact on the initial PBH clustering and mass function due to mode coupling between long and short wavelength modes. We show that even a small amount of non-Gaussianity results in a significant enhancement on the PBH initial clustering and subsequent merger rate and that the PBH mass function shifts to higher mass PBHs. However, as the clustering becomes strong, the local number density of PBHs becomes large, leading to a large theoretical uncertainty in the merger rate.


Introduction
Primordial black holes (PBHs) are black holes which may have formed very early in the history of the universe. There are several mechanisms by which they may have formed (for example, from cosmic strings [1] or bubble collisions [2]), but we will here focus on PBHs which form from the collapse of large density perturbations, as proposed in [3,4]. If a density perturbation has a large enough amplitude, it will collapse to form a PBH upon horizon entry. Such perturbations are sourced during cosmological inflation, where quantum fluctuations can become classical density perturbations as they exit the horizon, which then go on to re-enter the horizon following the end of inflation, and are also responsible for the growth of cosmological structure. The amplitude of such perturbations on small scales that form PBHs is required to be orders of magnitude larger than that observed on cosmological scales, although there are many models which do make this prediction (for example, [5][6][7][8][9][10][11], amongst many others).
The typical assumption is that the perturbations from which PBHs form have a Gaussian distribution, although recent papers have discussed the fact that, even if the curvature perturbation ζ is Gaussian, the density perturbations δ are not [12][13][14][15] -and it is the density which one should consider when investigating PBH formation [16]. There has been extensive research over the last decades to study the effects of primordial non-Gaussianity (PNG) on the abundance of PBHs [17][18][19][20][21][22][23][24][25][26] -finding that the presence of PNG can have contrasting effects on the abundance depending on the type of non-Gaussianity considered, but we note that local non-Gaussianity normally greatly enhances the abundance of PBHs, even with a negative skewness (f NL < 0) unless the power spectrum is narrowly peaked [27]. References [28][29][30][31] also studied the clustering of PBHs in non-Gaussian conditions, finding that even a small amount of PNG leads to significant clustering, which is otherwise absent.
In this paper, we will go beyond the work of previous papers, which have mostly focussed on the effect of PNG on the abundance of PBHs, and consider other key potential observables of PBHs as a function of the level of PNG: the PBH merger rate today and the PBH mass -1 -2 Primordial black hole abundance in a Gaussian universe Density perturbations above a certain threshold value will collapse upon horizon re-entry. The most robust criterion to use is the volume-averaged density contrast (also called the smoothed density contrast) [49]. Smoothing with a real-space top-hat function, the threshold value can take a value in the range 0.44 ≤ δ c ≤ 0.67, depending on the shape of the perturbation [50]. For a typical profile shape expected from the primordial perturbations generated, this takes a value δ c = 0.51 [49,51], which will be used throughout this paper. The effect of non-Gaussianity in the primordial distribution of matter and its effect on the profile shape, and its corresponding effect on the threshold value, was studied recently [52,53], finding that the critical value of the volume-averaged density contrast does not change significantly.
Perturbations which form PBHs in the early universe are necessarily rare. For example, in the case of solar mass PBHs, only one billionth of the universe needs to collapse into a PBH in order for PBHs to make up all of the dark matter today. Rare perturbations are expected to be approximately spherically symmetric [54]. Assuming spherical symmetry therefore, the smoothed density contrast δ R can be related to the curvature perturbation ζ in the comoving -2 -JCAP03(2020)004 synchronous gauge as Rζ (R) 2 + Rζ (R) , (2.1) where R is the smoothing scale (r m is typically used in the literature to represent the characteristic scale of a perturbation in the calculation of δ c , also representing the correct smoothing scale for a given perturbation), and the prime denotes a derivative with respect to the radial co-ordinate r. This equation allows us to relate the statistics of δ R , which determines PBH formation [16], to the statistics of ζ. The linear component of δ R , (which we will label as δ 1 ), is given by The mass of a PBH that forms from the perturbations depends upon the scale and amplitude of the perturbation as where K ≈ 4 for most profile shapes, and γ = 0.36 during radiation domination [12,[55][56][57][58], and M h is the horizon mass of a flat FRW universe at the time the Hubble horizon is the same scale as the perturbation. We will begin by assuming that ζ has a Gaussian distribution, before accounting for the effect of primordial non-Gaussiniaty in section 3. Following the method described in [12], the mass fraction of the universe collapsing to form PBHs at the time of formation is given by where δ c,1 is the critical value for the linear component, and finally σ j are moments of the power spectrum P δ 1 , given by where η is conformal time. During radiation domination, the power spectrum P δ 1 can be calculated from the curvature perturbation power spectrum P ζ as whereW (k, R) is the Fourier transform of the (real-space top-hat) smoothing function with a smoothing scale R = (aH) −1 , which is the Hubble scale, and T (k, η) is the linear transfer function at a time η,

JCAP03(2020)004
Importantly, on super-horizon scales, P δ 1 is proportional to the wave number to the fourth power, P δ 1 ∝ k 4 . We therefore conclude that (when perturbations follow a Gaussian distribution), super-horizon modes at the time of formation have a negligible effect on PBH formation. We will now discuss, briefly, the consequences of this in terms of the primordial clustering of PBHs arising from scale-independent bias. Scale-independent bias essentially happens because small-scale peaks can be situated inside large-scale peaks -meaning that many more small-scale peaks are likely to be above the threshold value in regions of the universe where there already exists a large-scale peaks, leading to the conclusion that PBHs preferentially form in large-scale overdense regions. However, because such large-scale overdensities are super-horizon at the time of PBH formation, they are negligibly small -and thus the effect of scale-independent bias is negligible on scales significantly larger than the scale of the PBH. This has been well documented in numerous papers [28-30, 45, 59-61], and also means that the spatial distribution of PBHs is expected to be Poissonian in the case of Gaussian statistics.
As described above, the perturbations which form PBHs are necessarily rare in order that their abundance does not quickly come to dominate the energy of the universe. For PBHs of around 1 solar mass, an initial abundance of β ∼ 10 −9 will cause them to dominate the energy fraction of the universe after the time of matter-radiation equality. Since the size of a PBH is approximately equal to the scale of the perturbation it formed from at horizon re-entry, the average separation between PBHs will be O(10 3 ) times greater than the horizon scale at the time of formation. If we take this scale as being the smallest scale relevant for the clustering of PBHs, then the amplitude of modes which may give rise to scale-independent bias is suppressed by a factor of at least O(10 −6 ) -and may safely be neglected.

Primordial black hole abundance in a non-Gaussian universe
We will now discuss the effects of non-Gaussianity on the abundance and initial clustering of PBHs. We will study the effect of local-type non-Gaussianity to second order, where the curvature perturbation ζ is related to the Gaussian-distributed ζ G as [62] where f NL is the non-linearity parameter (the δ 2 G term ensures that the mean of ζ remains zero while ζ G has a mean of zero). This model for non-Gaussianity is useful for several reasons: firstly, it allows us to make an analytic estimate of the effects of non-Gaussianity, 1 and secondly, local-type non-Gaussianity contains a strong coupling between the small-scale modes on which PBHs form, and the large-scale modes at which they cluster. It is also the typical type of non-Gaussianity generated during multiple-field inflation, e.g. [63] for a review.
It has been shown that local-type non-Gaussianity can lead to large dark-matter isocurvature modes, which are tightly constrained on CMB scales [28,29], leading to strong constraints on the non-Gaussianity in the scenario that a significant fraction of dark matter is composed of PBHs, f NL < O(10 −3 ). However, this constraint only applies to mode-coupling to scales large enough to be observable by the Planck satellite on the CMB and also does not apply in the case of single-field inflation. The constraints can be avoided if the bispectrum is assumed to be negligibly small on CMB scales, or uncorrelated to the PBH forming scales, JCAP03(2020)004 x ζ x δ Figure 1. A schematic plot of a universe containing exactly 2 modes. The x-axis represents spatial coordinates, whilst the y-axis represents the curvature perturbation ζ and the density contrast δ in the top and bottom plots respectively, at the time the small-scale wavelength enters the horizon. The red dashed line represents the threshold value in the density contrast in order for a PBH to form, whilst the black circles represent locations where a PBH will form. In a non-Gaussian, the two modes couple to each other, and the amplitude of the short-wavelength mode depends on the amplitude of the long-wavelength mode. This leads to larger perturbations in some regions of the universe -and consequently enhanced PBH formation in those regions.
but it may be larger on intermediate scales -which, as we will see, will result in significant clustering of PBHs on such scales.
The volume-averaged density contrast δ R is related to the curvature perturbation ζ during radiation domination by equation (2.1). Differentiating equation (3.1) to find ζ gives which gives an expression for δ R in terms of the Gaussian variable ζ G , This expression now depends not only upon derivatives of ζ, but on the absolute value of ζ itself. Therefore, the argument made in the previous section that large-scale modes do not affect PBH formation (and that therefore PBHs do not initially form in clusters) is no longer valid. This dependance of δ on the absolute value of ζ leads to modal coupling and initial clustering of PBHs, an example is shown schematically in figure 1. The consequences of this were discussed in detail in [27][28][29], which discussed how the abundance of PBHs and the amplitude of of isocurvature modes is affected. Over the course of this paper, we will extend the calculation to account for the non-linear relationship between ζ and δ and the critical scaling relationship, equation (2.3), as well as going on to calculate the effect on the mass function of PBHs, the primordial clustering, and the observed merger rate today.
In principle, ζ G and its derivative ζ G will be correlated. For example, a short distance away from the centre of a peak in ζ G , ζ G is likely to be positive, whilst ζ G is likely to be negative -and the opposite is true for troughs. However, the exact relation between ζ G and JCAP03(2020)004 ζ G depends on the profile shapes of the perturbation, which itself depends on the mechanism by which the perturbations were initially generated, and such a consideration goes beyond the scope of this work. 2 We will therefore make the assumption that ζ G is independent of ζ G , which can be treated as equivalent to assuming that the distribution of ζ is Gaussian in small regions of the universe, with the variance of perturbations varying from region to region. To validate this statement, we will split ζ G into long-and short-wavelength components, where only the short-wavelength modes will be relevant for PBH formation (the longerwavelength modes being super-horizon at the time of formation). Due to the Gaussianity of ζ G , ζ s will be uncorrelated with ζ l . Inserting this into equation (3.1) and taking the first derivative (as is necessary to calculate the formation criterion given by equation (2.1)) gives The derivatives of ζ l will be negligibly small compared to derivatives of ζ s , and so will be neglected. Finally, if we make the approximation that ζ s has a Gaussian distribution and therefore set f NL ζ s = 0, we obtain the relation which matches equation (3.2), with ζ G → ζ l and ζ G → ζ s , where the zeroth and first-order derivatives uncorrelated. For the sake of clarity in the rest of this paper, we will use the notation ζ l and ζ s . We note that the assumption of Gaussianity on the PBH forming scales will not significantly affect the main results of this paper: 3 firstly, whilst PBH abundance depends on the level of non-Gaussianity, it is degenerate with the amplitude of the small-scale power spectrum, which is treated here as a free parameter. Secondly, whilst the PBH merger rate does depend on the mass function of PBHs which itself will depend on the small-scale non-Gaussianity, this is a relatively small effect compared to the effect of the PBH abundance.
In a given region of the universe, with a constant ζ l , we can then simply treat ζ as a Gaussian variable, where the amplitude of the perturbation is simply modified by a factor (1 + 6 5 f NL ζ l ), which is essentially modifying the local power spectrum. Writing the variance of the small-scale perturbations σ 2 s as a function of f NL ζ l and the "background" (i.e. long wavelength) variance σ 2 b gives As described in section 2, the abundance of PBHs depends exponentially upon σ 2 , and so even small changes in σ 2 can mean a large change in the PBH abundance: it will be greatly amplified in regions of positive f NL ζ l , and greatly reduced for negative values -an effect often described as scale-dependant bias.

JCAP03(2020)004
In small regions of the universe, where ζ s can be treated as Gaussian, we follow the method of [12], and the local abundance of PBHs at the time of formation can be calculated with equation (2.4) as a function of f NL ζ l and σ 2 b by substituting σ 0 → σ s . The notation β local will be used to describe the local value of β in a small patch of the universe with constant ζ l . The fraction σ 1 /σ 0 appearing in equation (2.4) is not affected by the long wavelength mode because σ 1 and σ 0 are changed by the same factor.
For simplicity, we will here assume that the power spectrum has a Dirac-delta peak at some scale k * , and takes some smaller (but not necessarily constant) value at all other scales, (3.8) such that all PBH formation occurs at a single scale. The term P b (k) represents the (smaller) value of the power spectrum at all other scales, k = k * . We have made this assumption for simplicity, although our method can be extended to account for any power spectrum, which may have an extended PBH formation time.
The parameter δ β is introduced to describe the relative change in the abundance of PBHs, where β 0 (σ 2 b ) simply gives the "background" value for β in the absence of a large-scale ζ perturbation. 4 Figure 2 shows δ β as a function of f NL ζ l , with σ b = 0.015 and 0.010 for the solid blue and dotted red lines respectively. The right-hand plot shows a zoom of the central region. It can be immediately seen that f NL ζ l has an impact on the abundance of PBHs orders of magnitude larger than f NL ζ l , and a perturbative treatment will not provide accurate results except for very small values of f NL ζ l . 5 The total abundance of PBHs at formation may be obtained by integrating the local values of the abundance over the entire range of values of f NL ζ l , (3.10) It is found that the modal coupling due to non-Gaussianity always increases the number of PBHs form (see [27], where this was studied extensively).

Amplitude of non-Gaussianity
We will briefly discuss the amount of non-Gaussianity which may be expected in the primordial universe. In our model, the non-Gaussianity parameter f NL always appears multiplied by the long-wavelength component of the curvature perturbation, ζ l , and so the magnitude of the non-Gaussianity will be parameterised by the variance of this factor, f 2 NL ζ 2 l .
JCAP03(2020)004 Observational constraints by the Planck satellite require |f NL | 10 provided that it is scale independent (and we are focusing on the local model of non-Gaussianity). |f NL | ∼ 1 is a natural value in a range of multifield inflationary scenarios [64,65], although even in the multifield case |f NL | of order the slow-roll parameters is very common [66,67]. The longwavelength mode may be taken as the root mean squared value of the power spectrum (or larger by a factor of a few if the power spectrum takes a constant value over a large range of scales [27]), which is a scale dependent quantity. In the case of a very broad peak in the power spectrum (potentially giving rise to a broad PBH mass function), it may be realistic to take ζ l as large as 10 −1 , but given that the CMB µ-distortion constrains the primordial power spectrum to satisfy P ζ 10 −4 on scales not much larger than the scales required to generate LIGO mass PBHs, it is more reasonable to assume ζ l 10 −2 when studying the merger rate observed by LIGO and Virgo. As a lower bound, we should consider the observed CMB modes, which require ζ l 10 −4 . See figure 8 of [68] for an example plot showing the main power spectrum constraints. Combining all of this, values of f NL ζ l of order unity is possible but an upper bound, whilst for the motivated case that |f NL | ∼ 1 we should expect 10 −4 f 2 NL ζ 2 l 10 −2 . In the following sections, we will discuss the implications of the scale-dependent bias on the mass function of PBHs, the power spectrum of the PBH and dark matter density, and the observed black hole merger rate today.

Mass function
In addition to their abundance, a key (potential) observable of PBHs is their mass. We will now derive the mass function, and how it is affected by primordial non-Gaussianity. The mass function will be defined as where f PBH = Ω PBH /Ω DM is the ratio of PBHs to dark matter, and ψ(m) is defined such that ψ(m)dm = 1, where we have followed the definition in [46]. 6 The parameter f PBH is JCAP03(2020)004 the fraction of dark matter composed of PBHs, and can be calculated from β as where M H is the horizon mass at the time of horizon-entry, M eq = 2.8 × 10 17 M is the approximate horizon mass at matter-radiation equality [69]. The term accounts for the relative red-shift of the PBH density (which evolves like matter) relative to the dominant radiation density from the time of PBH formation until the time of matter-radiation equality. See [12,70] for further discussion. Assuming, for the time being, that ζ follows a Gaussian distribution, the mass function of PBHs is therefore defined completely by the primordial power spectrum P ζ .
In order to derive a mass function, we will, again, make the simplifying assumption that all PBH formation occurs at a single epoch, with a Dirac-delta form for the power spectrum, as in equation (3.8). It is noted that such a form for the power spectrum is not physical, and the fastest growth for the power spectrum as a function of the wavevector k has been shown to be ∼ k 4 [68], or slightly steeper in the case of a more contrived scenario [71] (see also [72]), which means that PBH formation at a single scale is not physical, but rather that it occurs over a range of scales. However, our numerical calculations assuming the narrowest possible peak in the power spectrum show that the mass function would only widen by a factor O(0.5%) compared to the Dirac-delta power spectrum, and so we will neglect it here in order to perform the calculations analytically. Again, the assumption of a Dirac-delta function for P ζ does not affect conclusions about the mass function, because PBH formation at each scale would be affected in the same manner. In this case, the mass function can be given as The mass, m, of a PBH is related to the amplitude of the perturbation from which it formed, by the well known critical scaling law, equation (2.3), which is inverted to give It is noteworthy here that there is a maximum value for the PBH mass which can form from perturbations of a given scale, due to the fact that equation (2.1) has a maximum at δ R = 2/3. This gives a maximum PBH mass of which gives M PBH,max = 2.05M H for the parameter choices considered here, see (2.3). Perturbations in ζ larger than that corresponding to the maximum δ R are possible, but these result in perturbations of a type for which the dynamics of PBH formation are not well understood [73]. In practice, this has a negligible effect on the calculation, since the abundance of such large ζ perturbations is exponentially suppressed.  Figure 3. The mass function ψ is shown for PBHs forming at a single time from a Dirac-delta power spectrum. The mass function generally peaks at approximately the horizon mass at the time of horizon crossing, shifting to higher masses as the variance of density perturbations, σ 2 0 , increases. Note that the mass function is normalised to integrate to unity, whilst the actual abundance of PBHs is larger by many orders of magnitude for larger values of σ 2 0 .

Substituting equation (4.4) into equation (4.3) gives the final expression for the mass function assuming a Gaussian distribution of
6) where σ 1 σ 0 = 1 for the Dirac-delta power spectrum considered here, and the subscript G denotes that this is valid for a Gaussian ζ. Figure 3 shows the mass function for PBHs generated at a single time as a function of the horizon mass. The mass function for several values of σ 2 0 is plotted, and it can be seen that larger values for σ 2 0 mean that the mass function peaks at higher masses, but the peak is broader and lower. Note that, whilst the mass function ψ is normalised to 1 when integrated over m, the actual abundance of PBHs is greatly suppressed for smaller values of σ 2 0 -such that more PBHs of all masses are actually produced for larger σ 2 0 . We will now turn our attention to the mass function in the presence of non-Gaussianity. As discussed in the previous section, this will be described by treating separate regions of the universe as having a Gaussian distribution, with the variance of perturbations modified by a long-wavelength perturbation, as in equation (3.7). In this case, the full mass function can be given in terms of ψ G as, (4.7) Figure 4 shows the mass function for various parameter choices. The left plot shows the mass function for σ 2 b = 0.005, for 3 choices of f 2 NL ζ 2 l . As f 2 NL ζ 2 l increases, the PBHs formed in the universe become more and more clustered -such that the PBH abundance becomes quickly dominated by PBHs formed in high-β regions. Therefore, the same behaviour for the mass spectrum is seen as for the Gaussian case: that the peak becomes lower, broader and shifted to the right as f 2 NL ζ 2 l increases, and again the total abundance of PBHs also increases dramatically. The right plot of figure 4 therefore shows a similar plot, except the total abundance of PBHs is held fixed at f PBH = 0.01 for different values of f 2 NL ζ 2 l . We now see that as f 2 NL ζ 2 l increases (which increases the PBH abundance), the background variance σ 2 b is decreased (which decreases the PBH abundance). The change in the mass function is therefore not as pronounced, due to these factors having opposite effects, although it is still seen that the peak broadens, lowers, and shifts to the right.
This predicts a unique mass function dependent on the abundance of PBHs and the level of non-Gaussianity, but does depend on the form we have specified for the power spectrum.
Given that the small-scale power spectrum is degenerate with the level of non-Gaussianity, it is not possible to use the mass function by itself (which is in principle measurable if PBHs are observed) to determine the level of non-Gaussianity or amplitude/shape of the power spectrum.

The primordial black hole number density power spectrum
In this section, we will derive a simple analytic estimate for the power spectrum of the PBH density perturbations, P PBH . To begin, we will greatly simplify the calculation of β. The dominant term in the integral in equation (2.4) to calculate β is the exponential term. Neglecting the other terms gives the standard result from the Press-Schechter calculation, where ν c = δ c,1 /σ 0 . Extending this to include the modal coupling to large-scale modes, the expression for the local PBH abundance at formation gives We will use the same definition for the PBH density perturbation δ β as previously, equation (3.9), assuming that f 2 NL ζ 2 l is small and that the total PBH abundance is not changed significantly from the background value. Expanding the resulting expression to first order in ν in the ν → ∞ limit (retaining only the leading order terms), and to first order in f NL ζ l around zero gives a linear expression for δ β in terms of ζ l , or to second-order, Note that, while we have ignored the subdominant terms in equation (2.4) it is possible to include them, and the full numerical result relating δ β to ζ l is shown in figure 2. From this linear expression, it is simple to relate the power spectrum of ζ to the power spectrum of δ β , which is consistent with the expression more rigorously obtained in [30], if we have the standard picture that 6 5 f NL 2 = τ NL from single-source inflation (a scenario in which any one field is responsible for generating the curvature perturbation, such as the standard curvaton or modulated reheating models) [63,74]. Including the second-order term gives Figure 5 shows a comparison of the value of δ β obtained via the full numerically integrated expression given by equation (3.9), and the first-and second-order expansions obtained above. For the numerically obtained result, we have assumed a Dirac-delta form for the small-scale power spectrum as before, as in (3.8), with an amplitude of A s ≈ 0.0938 (such -12 -JCAP03(2020)004 that f PBH = 1 in the absence of large scale modes), giving ν c = 5.77 used in the first-and second-order expansions.
Here, contributions to the power spectrum from shot noise and adiabatic perturbations have been neglected, although they are well understood, as neither will contribute to the calculated merger rate. The shot noise has no effect because it is a consequence of the random positions of PBHs, which is already accounted for in the calculation of the expected merger rate. However, one should consider the signal arising in the power spectrum from shot noise if searching for evidence of non-Gaussianity [31]. The adiabatic term is due to the different densities in different regions of the universe. However, at early times when PBHs form, the much larger modes under consideration here are in the super-horizon regime, and simply correspond to time-shifts in a smaller region of the universe -as described in the separate universe approach (e.g. [75,76]).

The primordial black hole merger rate
Over recent years, multiple observations of gravitational wave signals from the merging of binary black hole systems has been observed by LIGO [77,78], Since the paper "Did LIGO detect dark matter?" [37], there has been an extensive amount of literature dedicated to the discussion of the formation of binary PBHs and the gravitational wave signal from merging PBHs [38,41,43,44,46,[79][80][81][82][83][84][85][86]. The general consensus is that, were PBHs to make up the entirety of dark matter, the rate of observed BH merger events would be significantly higher than actually detected -thereby ruling out PBHs in the approximate mass range 1-100M from composing the entirety of dark matter.
In this section, we will consider the effect on the merger rate due to primordial clustering caused by PNG. We consider a scenario where PBH binaries form shortly after the formation of the PBHs themselves. At the time of formation, PBHs are coupled to the expansion of the universe, and will typically become gravitationally bound to each other when their local density becomes greater than the surrounding radiation density, normally after the time of matter-radiation equality. However, due to the randomness of the Poisson fluctuations, some PBHs will form much closer to each other than average, and will therefore decouple significantly earlier. PBH pairs may therefore start falling towards each shortly after formation, with a direct collision avoided by torque provided by gravitational forces from nearby matter perturbations and other PBHs. The pair of PBHs then forms a binary system, which will eventually merge. The merger may be observable by gravitational wave detectors.
By calculating the distribution of orbital parameters from the distribution of initial conditions leading to PBH formation, reference [46] gives an analytic expression for the merger rate of PBHs today, which reads where dR 0 is approximated as with τ the coalescence time, t 0 = 13.8 Gyr the current age of the universe, m 1 and m 2 the masses of the PBHs in a binary pair, and η is the symmetric mass fraction,

JCAP03(2020)004
We will take the maximum value for the suppression factor S (valid in the regime where the expected number of PBHs in a spherical volume with a radius equal to the initial separation between the PBH binary is small), given by where U is the confluent hypergeometric function, and σ 2 M = (Ω M /Ω DM ) 2 σ 2 f is the rescaled variance of matter density perturbations, evaluated here at the time of matter-radiation equality, which is close to the time when binaries that merge today originally formed. We will here follow previous works [42,46,87] and take σ f = 0.005. In a Gaussian universe, all of the remaining factors affecting the merger rate can therefore be derived purely from the variance of density perturbations, R = R σ 2 0 , again assuming a Dirac-delta form for the enhanced part of the power spectrum.
The expression derived by reference [46] assumed that PBHs were distributed across the universe with a uniform Poisson distribution and a constant mass function. Extending this to the non-Gaussian distribution considered here is straightforward, requiring the merger rate to be integrated over different regions of the universe in the same manner as the PBH abundance or the mass function. Calculating the local value for the merger rate also requires the calculation of the local value of the PBH DM fraction, f PBH , and the local PBH mass function, ψ G . Both of these factors depend on the local variance, σ 2 s , given in equation (3.7) as a function of f NL ζ l , which follows a Gaussian distribution.
The expression for the total merger rate R T is then given as, (6.5) However, reference [46] also tested their analytic prediction for the merger rate with N -body simulations, and found that, when the PBH density becomes large, f PBH 0.1, the initial binaries are likely to be disrupted by nearby PBHs which is likely to result in a reduced merger rate. As a result, equation (6.1) used for the merger rate cannot be trusted to be accurate for such regions. To account for this uncertainty as much as is currently possible, we have posited several different models. We have introduced a "hard" and "soft" cut-off in the integral above, where the merger rate is set to zero or to the value corresponding to f PBH = 0.1 in regions where f PBH > 0.1, respectively.
The merger rate today will then depend on 2 factors: the abundance of PBHs, and the level of non-Gaussianity, parameterised by f PBH and f 2 NL ζ 2 l respectively. Figure 6 shows the predicted merger rate today as a function of f 2 NL ζ 2 l for f PBH = 0.01 and 0.001. Where f 2 NL ζ 2 l = 0, our expression is reduced to that given by [46]. It can be seen that as the level of non-Gaussianity increases, the clustering also increases, which leads to a rise in the predicted merger rate in all cases. However, as the level of non-Gaussianity continues to rise, the PBH abundance in the universe starts being dominated by regions where the local PBH density is very high -which means equation (6.1) can no longer be trusted. In this case, the merger rate may start dropping -as shown by the dashed-and dotted-lines, representing hard and soft cut-offs respectively. Above the point at which the lines diverge, there is still a great deal of uncertainty in what the merger rate may be. The dashed and solid lines can be considered as lower and upper bounds respectively, with a "more realistic" estimate given by the dotted lines.   Figure 6. The effect of local-type non-Gaussianity is plotted against the predicted merger rate for PBHs. The merger rate is plotted for 2 different total abundances of PBHs, f PBH = 0.01 and f PBH = 0.001 in blue and red respectively. For small amounts of clustering, the merger rate is expected to increase, but there is significant uncertainty when the clustering becomes very large -due to the uncertainty in the evolution of binary systems in PBH-dense regions.

Smallest value of f PBH which might produce the observed merger rate
In this section, we will provide a simple order-of-magnitude estimate for the lower limit of PBH abundance which may give rise to the observed BH-BH mergers observed by LIGO. Whilst there are significant uncertainties arising in the calculation from the current uncertainty in the merger rate when the local PBH abundance becomes large, in this section we will neglect this in order to give a rough estimate for the lower bound. We will assume a simple picture where PBHs form in extremely high abundance in certain regions (hereafter referred to high PBH density (HPD) regions), and in negligible amounts elsewhere. 7 Such a scenario may be realised within the framework of our model if the variance of f NL ζ l is extremely large, leading to a small number regions with a very large small-scale power spectrum with high PBH abundance, and a large number of regions with small power spectrum and low PBH abundance (although the full calculation described above breaks down when PBHs can no longer be considered rare events).
In order to provide the required number of PBH mergers, we will take that the minimum merger rate to be 10 Gpc −3 yr −1 , that PBHs form when the horizon mass is equal to 20 M , and that the maximum formation rate in any HPD region is β = 0.1. In such HPD regions, the PBH energy density would be O(10 7 ) times greater than the dark matter density in the universe, f PBH ∼ 10 7 . The method described in section 4 breaks down when PBHs are not rare events, and so we will assume a lognormal mass function given by,

JCAP03(2020)004
where we take values m c = 20 M and σ m = 0.1 (the result is not especially sensitive to the exact values). The merger rate predicted by equation (6.1) is then R = O(10 15 ). In order to produce the minimum merger rate of 10 Gpc −3 yr −1 , such regions must therefore occupy approximately 10/10 15 = 10 −14 of the entire volume of the universe. Therefore, the minimum fraction of dark matter composed of PBHs in order to be responsible for the BH-BH merger events observed by LIGO is We stress that this number is intended to be a crude order-of-magnitude estimate only, and assumes that the analytic prediction for the merger rate holds for (locally) very large PBH abundances. In reality, the lower bound is likely to be significantly higher.

Summary
We have studied the way in which primordial non-Gaussianity may affect key observables related to the production of PBHs in the early universe. Unlike for a purely Gaussian distribution, primordial non-Gaussianity arising from inflation, especially local-type, is expected to result in significant modal coupling. We have studied the effect that such modal coupling would have on the abundance of PBHs, the mass function, the power spectrum of PBH number density perturbations on larger scales at the time of formation, and the merger rate of PBHs observable today from binary PBH systems. The effect of modal coupling is found to always increase the abundance of PBHs, regardless of the sign of f NL . This is due to the exponential dependance of the PBH abundance on the amplitude of the power spectrum. In the case of positive (negative) f NL , PBH formation is enhanced in regions of positive (negative) ζ l , and vice-versa. This leads to a significant variation in the formation rate of PBHs in different regions of the universe, δ β . In section 5, we calculated the power spectrum P PBH to second order in terms of the power spectrum of the curvature perturbation P ζ . A possible future observation of such perturbations could therefore provide insights into the early universe.
We have discussed the mass function of PBHs which form, assuming a very narrow peak in the power spectrum -in our case, a Dirac-delta peak. The mass function peaks at approximately the horizon mass at the time perturbations enter the horizon. The location of the peak shifts to higher masses when the amplitude of the power spectrum is larger (or in regions where the local amplitude is higher) -corresponding to an increase in the formation rate of PBHs. Therefore, in addition to the variance in the abundance of PBHs in different regions of the universe, variations in the mass function of PBHs may, in the future, provide information about primordial physics.
The effect of modal coupling from non-Gaussianity on the possible observed GW signal today from merging binary PBH systems has also been considered. In the case that the modal coupling is relatively small, f 2 NL ζ 2 l 4 × 10 −4 , the effect of modal coupling can be safely stated to increase the predicted merger rate. This means that much smaller abundances of PBHs than previously thought could still produce a high enough present day merger rate. We have estimated that the minimum amount of dark matter composed of PBHs which could still be responsible the merger events observed by LIGO is f PBH ∼ 10 −7 . A lower PBH abundance is therefore excluded from producing the observed GW signals.
However, the formation rate of binary systems in regions of high PBH abundance (i.e. where the PBH density is greater than 10% of the background dark matter density) JCAP03(2020)004 is currently not well understood -and this means that there is significant uncertainty in our calculation of the merger rate when the initial clustering becomes too large. It may even have the effect of reducing the merger rate sufficiently such that the merger rate constraints on f PBH may be weakened if there is significant PBH clustering. We note that a paper which appeared while we were completing this work [47] have estimate the weakest possible constraint from the merger rate by considering the merger rate of PBHs from binary pairs which have been disrupted. PBH abundances previously thought ruled out by the lack of observed merger events may not be ruled out.
Throughout this paper, we have ignored the effect of non-Gaussianity on PBH scales -although the results presented would be qualitatively unchanged by its inclusion. Full consideration of this will require consideration of the the shape of the small-scale power spectrum, typical profile shapes and the shape of the bispectrum, amongst other factors, and goes beyond the scope of this paper. However, much work has been completed on this topic in the past [15,[19][20][21][22][23][24][25]27], and a more complete analysis using the methods presented in this paper would yield qualitatively similar results. We have also only considered an expansion to second order in local-type non-Gaussianity, although it has been shown that higher-order terms can also be important [21].