Non-perturbative non-Gaussianity and primordial black holes

We present a non-perturbative method for calculating the abundance of primordial black holes given an arbitrary one-point probability distribution function for the primordial curvature perturbation, $P(\zeta)$. A non-perturbative method is essential when considering non-Gaussianities that cannot be treated using a conventional perturbative expansion. To determine the full statistics of the density field, we relate $\zeta$ to a Gaussian field by equating the cumulative distribution functions. We consider two examples: a specific local-type non-Gaussian distribution arising from ultra slow roll models, and a general piecewise model for $P(\zeta)$ with an exponential tail. We demonstrate that the enhancement of primordial black hole formation is due to the intermediate regime, rather than the far tail. We also show that non-Gaussianity can have a significant impact on the shape of the primordial black hole mass distribution.


I. INTRODUCTION
One of the unanswered questions in ΛCDM cosmology is the nature of dark matter. One suggestion that has attracted renewed interest recently is that of primordial black holes (PBHs), which could have formed from the collapse of very overdense regions in the early universe [1][2][3]. If they survive to the present day, they would constitute some or all of the dark matter content of the universe. This is characterised through the parameter f PBH , the fraction of dark matter in the form of PBHs, i.e. f PBH = 1 indicates that PBHs make up the entirety of the dark matter. For a recent review of constraints on f PBH , see Ref. [4]. Even if they are not the main constituent of dark matter, the merger of PBH binaries will also generate gravitational waves, and may contribute to the LIGO-Virgo-KAGRA detections [5][6][7][8][9].
The simplest formation mechanism for PBHs is from overdensities seeded by inflation, a period of accelerated expansion in the very early universe that also explains the seeding of large scale structure. The perturbations from inflation are typically characterised by the primordial power spectrum of the curvature perturbation, P ζ (k), which is observed on cosmic microwave background (CMB) scales to be ∼ 2.1 × 10 −9 [10]. It is commonly assumed that to form PBHs requires increasing P ζ (k) by several orders of magnitude on small scales [11]. However, this assumes that the probability distribution function (PDF) of the curvature perturbation is Gaussian. PBHs are formed in the tail of the PDF, and so non-Gaussianity that enhances the occurrence of rare events could have a significant impact on PBH formation, potentially reducing the power required to form them.
Previous studies of PBH abundance in the presence of primordial non-Gaussianity have usually assumed a local, perturbative expansion in terms of a first-order Gaussian curvature perturbation, ζ G with higher-order parameters f NL , etc. [12][13][14][15][16][17]. However, a non-perturbative calculation can be made using the δN -method [18]. It has been argued that by fully incorporating quantum diffusion, in the so-called stochastic δN formalism [19,20], we expect the far tail of the probability distribution for primordial curvature perturbations to generically display a non-Gaussian, exponential tail [21][22][23][24][25][26][27][28][29][30] (such heavy tails were also found in Refs. [31][32][33] using different methods). Therefore it is important to develop and apply techniques which enable us to explore more general forms for the PDF, and for the tail of the distribution in particular, which encompass a wider range of non-Gaussianity. In this paper, we present for the first time a fully non-perturbative treatment of primordial non-Gaussianity that makes no assumptions about the model other than that ζ is related to a Gaussian field ζ G by what we will call a "generalised local" transformation.

II. METHOD
Primordial black holes form from rare, large fluctuations in the primordial density field. The power spectrum is no longer sufficient to determine the abundance of large fluctuations if we wish to go beyond the simplest Gaussian distribution. Instead we will start from an arbitrary PDF for the primordial curvature perturbation P [ζ(x)], e.g., obtained from an inflationary model.
The criteria for formation of PBHs can be given in terms of the compaction function, C(r). This is a volume-averaged version of the density contrast, equivalent to smoothing with a real-space tophat window function [34,35]. The relationship between ζ and C is non-linear, but the compaction during radiation domination can be written in terms of its linear part C ℓ as with C ℓ related to ζ by where a prime ′ denotes a derivative with respect to r. This simple relation makes it useful to calculate the PBH properties in terms of P (C ℓ ) instead of P (C). Eq. (1) has a maximum of C = 2/3 for C ℓ = 4/3, which corresponds to the boundary between type I and type II perturbations. It is unclear what the properties of a PBH resulting from a type II perturbation would be [36], so we restrict to C ℓ ⩽ 4/3 in this work.
For a Gaussian field, the high peaks in the density field that form PBHs can be approximated as spherical. We assume that this also holds in our non-Gaussian case, and we will show that this is a self-consistent assumption in our subsequent calculations. For the linear compaction in eq. (2) to be spherical, the curvature field ζ in the same region must also be spherical. In general, we can expand the field in spherical harmonics, From now on we will consider only the spherically symmetric monopole term in this expansion and use the shorthand ζ ≡ ζ(r) ≡ ζ 00 (r). The reliance on ζ ′ rather than ζ itself means that P (ζ) is not sufficient to calculate the PBH mass distribution and abundance, and instead the joint probability P [ζ(r 1 ), · · · , ζ(r n )] is required [37]. Therefore, in our analysis, we will assume that ζ(r) can be related to a Gaussian field ζ G (r) whose statistical properties are well-understood. We introduce a "generalised local" transformation, where we have allowed for an explicit r-dependence independent of the Gaussian field ζ G . This extra dependence is required to consider arbitrary P (ζ) since the Gaussian PDF will always carry an r dependence through the variance Σ Y Y as well as ζ G , as can be seen in eqs. (6) and (7). This transformation is a generalisation of a local transformation, in the sense that it includes the local case, where ζ is a local function of a Gaussian field, ζ = ζ[ζ G (r)]. For a monotonic transformation, high peaks in ζ necessarily correspond to high peaks in ζ G , and hence this relation is consistent with our spherical assumption. We want to connect an arbitrary non-Gaussian ζ with PDF P (ζ) to a Gaussian ζ G with PDF where we write the variance as Σ Y Y for reasons that shall become apparent. To do this, we utilise the fact that a variable defined by F X (X) is always uniformly distributed, where F X is the cumulative distribution function (CDF) for the variable X. This means we can equate the CDFs for ζ and ζ G , which can be applied for any P (ζ), whether analytical or numerical. The simplest generalised local transformation has all the r-dependence coming from ζ G (r) and the variance Σ Y Y (r), for which the non-Gaussian field can be written as where F −1 is the inverse CDF for ζ. In this case, ζ ′ in eq. (2) includes two terms, one from the dependence of ζ G on r, and one from the dependence of Σ Y Y on r. We can therefore write the linear compaction (2) as with Jacobian factors given by Note that by construction, from eq. (7), we have J 1 = P G (ζ G )/P (ζ). In the case of a local transformation where ζ depends only on ζ G , J 2 will be zero. We can see from eq. (8) that the linear compaction depends on both ζ G and ζ ′ G , and so P (C ℓ ) will be written in general as an integration over a 2D PDF, where C ℓ (X, Y ) is given by eq. (8), δ D is a Dirac delta function, and X and Y are given by The 2D Gaussian distribution, P (X, Y ), is given by with The correlators can then be written in terms of the power spectrum of the Gaussian field, P ζ G (k), where j 0 = (sin z)/z is the spherical Bessel function. Here we differ from Ref. [38] by retaining the Σ XY term, which is not zero for isotropy as claimed. Indeed, for a monochromatic Gaussian power spectrum, studied in Ref. [39], the fields are fully correlated, with Σ 2 XY = Σ XX Σ Y Y , and the joint PDF in eq. (12) collapses to a 1D case. In this work, to retain generality, we keep the cross-correlation term but assume a lognormal shape for the Gaussian power spectrum, with amplitude A, width ∆, and peak scale k * . This is normalised such that it approaches the monochromatic case of a Dirac delta in ln k in the limit ∆ → 0. We can then rewrite the linear compaction in terms of X and Y as Finally, carrying out the marginalisation in eq. (10) with C ℓ (X, Y ) given by eq. (19), The PBH mass m is related to the linear compaction C ℓ by the critical collapse relation [40][41][42] where K, C c , and γ are parameters that control the collapse, which in general are dependent on the shape of the perturbations [43][44][45][46][47]. The mass distribution is then given by [39] f where Z depends on the peak scale in the power spectrum and the background cosmology. The diverging term in the denominator arises from the Jacobian factor when changing variable from C ℓ to m, and is ultimately inherited from the non-linear form of C in eq. (1). The boundary between type I and type II perturbations at C ℓ = 4/3 corresponds to the maximum value of C = 2/3, and hence a zero derivative. The PBH abundance is determined from the mass distribution f (m) through where m max corresponds to C ℓ = 4/3.

III. RESULTS
We now apply the method described in sec. II to two examples of non-Gaussian distributions, P (ζ). We choose the power spectrum width ∆ = 0.3 in eq. (17), and take the same peak scale as [39], k * = 1.56 × 10 13 Mpc −1 . The scale r is chosen such that k * r = 2.74, which maximises the compaction function for a monochromatic power spectrum [48]. Additionally we choose the same critical collapse parameters as Ref. [39], K = 1, C c = 0.587, γ = 0.36 1 , and write the PBH mass in terms of the horizon mass M H evaluated at the peak scale in the power spectrum, M k * . For these choices, the factor in the mass distribution (22) is Z = 6.88 × 10 16 .
The first non-Gaussian distribution P (ζ) is an example of a local transformation arising from the classical δN approach for an ultra slow roll (USR) phase in the early universe. This results in the transformation [38,39,[49][50][51] and in the non-Gaussian PDF where Σ Y Y is the variance of the Gaussian field ζ G . From now on we shall refer to this as the classical USR δN approach. In this case, and assuming a monochromatic power spectrum, our general method reduces to the ratio distribution approach in Ref. [39].
As an example of a minimal generalised local transformation governed by eq. (7), we also introduce a test model for P (ζ) consisting of a piecewise matching between a Gaussian distribution and a simple exponential tail, The width of the Gaussian part is chosen to be identical to the distribution for ζ G , and the factor N is included to give the correct overall normalisation.  Figure 1 shows a comparison of the piecewise exponential model in eq. (26) with the Gaussian and classical USR δN cases. We choose α = 3 in order to match the slope in the exponential tail of the classical USR δN distribution (25). We select two values for the piecewise transition: ζ t = 1/3, to match the far tail of the classical USR δN distribution as can be seen in the top-left panel, and ζ t = 0.15, to approximately match the f PBH curve in the bottom-right panel. The top-right panel shows P (C ℓ ), calculated using eq. (20). It is clear that for ζ t = 1/3, there is no enhancement in the piecewise model over the Gaussian case, indicating that the exponential tail in the classical USR δN model does not actually enhance the production of PBHs, with the enhancement instead coming from the transition region between ζ ≈ 0.2 and 1.6 in the top left panel. This is corroborated by the f PBH curves, where the ζ t = 1/3 case is identical to the Gaussian. It is also clear from the top right panel that the classical USR δN model, and the piecewise model with ζ t = 0.15 chosen to produce a similar PBH enhancement, have significant support for C ℓ values larger than the cutoff of 4/3. This suggests that if type II perturbations do collapse to form PBHs, they will contribute significantly to the total abundance in these non-Gaussian models.
The bottom left panel shows the mass distribution f (m), normalised by the total dark matter fraction f PBH . As expected, the Gaussian and ζ t = 1/3 distributions are identical. However, it is clear that the ζ t = 0.15 curve is significantly different, without the decaying behaviour for m > M k * . This is because the mass distribution defined in eq. (22) has a divergence at the mass corresponding to C ℓ = 4/3, coming from the (1 − 3 4 C ℓ ) term in the denominator. All of the mass distributions in fig. 1 have this divergence at the right hand limit of the mass range, but it is not visible for the other curves, because P (C ℓ ) is decaying faster than the term in the denominator. However, for the ζ t = 0.15 case, the top-right panel shows that P (C ℓ ) is decaying very slowly, so does not beat the diverging term. The mass distribution f (m) arising from the classical USR δN model does decay for m > M k * , despite having a far tail slope of α = 3 in P (ζ) like the piecewise model. This is because, as stated previously, the relevant part of P (ζ) is not the far tail, but the transition, during which the effective value of α is significantly larger, leading to a decaying P (C ℓ ) and a tail in the mass distribution. It should also be noted that the critical collapse relation in eq. (21) has been used right up to C ℓ = 4/3, but is expected to break down before this point [52]. Figure 2 shows the effect of changing the piecewise transition point ζ t and exponential slope α in eq. (26) on the total PBH abundance f PBH and the mass distribution f (m). As in fig. 1, the largest values of ζ t provide no enhancement over the Gaussian PDF except at the highest amplitudes, for which f PBH ≫ 1. For smaller values of ζ t , the exponential tail does provide a significant enhancement over the Gaussian case, greatly reducing the amplitude required for f PBH = 1. In contrast, changing α only moves the f PBH curve by a (relatively) small amount, e.g. by factors of 10 rather than the 10 10 coming from changing ζ t . The steeper values of α are suppressed compared to the shallower slopes, as expected. For large amplitudes, f PBH is suppressed (even compared to the Gaussian). This suppression comes from the normalisation of P (ζ) when including a large tail. The details of this suppression are unimportant, since they correspond to f PBH ≫ 1 and so these power spectrum amplitudes are disallowed anyway for PBH masses that survive to the present day.
In the bottom two panels we can see the effect of changing ζ t and α on the mass distribution f (m). For ζ t , there are three cases that are visible. The ζ t = 0.3 curve lies on top of the Gaussian, as for the ζ t = 1/3 case in fig. 1. For ζ t = 0.25, the function f (m) starts to decrease for m > M k * as P (C ℓ ) decreases, before growing again due to the diverging term discussed previously. Finally, all values of ζ t smaller than 0.25 show no drop-off for m > M k * , instead being dominated by the diverging factor. Comparing the f PBH and f (m) plots for varying ζ t , it seems clear that the amount of non-Gaussianity required to significantly reduce the power spectrum amplitude corresponding to a fixed f PBH will also have a large impact on the mass distribution. This is important in the case of evading tight constraints such as those from µ-distortions [53]. The right panel shows that, as with f PBH , the mass distribution is less affected by α than ζ t , with all the curves appearing basically identical. This is consistent with the conclusion from fig. 1 that the transition regime provides the enhancement, rather than the far tail.
Σ Y Y and P (ζ t ) required for a specific f PBH , with a given α. There is a linear relationship between P (ζ t ) and f PBH for constant α, and an inverse linear relationship between f PBH and α for fixed ν t , which comes about from the integration of the far tail exponential behaviour. α ν t P (ζ t ) f PBH

IV. CONCLUSIONS
In this work, we have developed a non-perturbative treatment of primordial non-Gaussianity, when it is of a generalised local type, and applied it to the formation of primordial black holes by relating the curvature perturbation ζ to the compaction function C. A non-perturbative method is essential when considering non-Gaussianities that are not described by the usual f NL expansion, such as those appearing in the stochastic treatment of inflation. Our method assumes there exists a generalised local transformation from a Gaussian field, ζ(ζ G , r), and it remains to be seen whether this is a good description of the non-Gaussianity generated by stochastic diffusion. Previous works have used the classical USR δN relation for a local-type non-Gaussian ζ(ζ G ) with a PDF that smoothly transitions from a Gaussian to an exponential tail. By comparing to a simple piecewise model with a sudden transition between a Gaussian and an exponential with the same slope, we have shown that the enhanced PBH formation in the classical USR δN model is due to the transition regime, rather than the far tail. Additionally, we have examined the transition value, ζ t , required for f PBH = 1 using our piecewise model. We find that ν t = ζ t / √ Σ Y Y is approximately independent of the peak amplitude of the Gaussian power spectrum, A, for fixed f PBH and α. We also see a linear relationship between P (ζ t ) and f PBH , and when changing the exponential slope α we find an inverse linear relation between f PBH and α, see table I. We also show that, for sufficiently shallow exponential tail slopes, the PBH mass distribution is dominated by a diverging term for high masses. This effect will be important when significant non-Gaussianity is required to reduce the power spectrum amplitude by orders of magnitude at fixed f PBH , e.g. to evade the µ-distortion constraints for supermassive PBHs.

NOTE ADDED
During the writing up of this paper, Ref. [54] appeared on the arXiv. The authors utilise the same mechanism of writing ζ in terms of ζ G and applying the product rule to the derivative in the compaction function, but restrict themselves to local transformations ζ(ζ G ), rather than allowing P (ζ) to be arbitrary.