<figure><img src="image_1.png" /></figure>

HAL
open science

# Non-perturbative non-Gaussianity and primordial black holes

Andrew D. Gow, Hooshyar Assadullahi, Joseph H.P. Jackson, Kazuya

Koyama, Vincent Vennin, David Wands

► To cite this version:

Andrew D. Gow, Hooshyar Assadullahi, Joseph H.P. Jackson, Kazuya Koyama, Vincent Vennin, et al.. Non-perturbative non-Gaussianity and primordial black holes. EPL - Europhysics Letters, 2023, 142 (4), pp.49001. 10.1209/0295-5075/acd417. hal-03877906

HAL Id: hal-03877906

https://hal.science/hal-03877906v1

Submitted on 2 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
LETTER • OPEN ACCESS

Non-perturbative non-Gaussianity and primordial
black holes

To cite this article: A. D. Gow *et al* 2023 *EPL* **142*** 49001

View the [article online](#) for updates and enhancements.

You may also like

- *Primordial black holes from a curvaton scenario with strongly non-Gaussian perturbations*
  Andrew Gow, Tays Miranda and Sami Nurmi

- *Revisiting small-scale fluctuations in -attractor models of inflation*
  Laura Iacconi, Hooshyar Assadullahi, Matteo Fasiello et al.

- *Search for Gravitational Waves from the Coalescence of Subsolar Mass and Eccentric Compact Binaries*
  Alexander H. Nitz, Yi-Fan Wang and
# Non-perturbative non-Gaussianity and primordial black holes

A. D. GOW<fnref n="1" />(a) ID, H. ASSADULLAHI<fnref n="1" />,<fnref n="2" />(b), J. H. P. JACKSON<fnref n="1" />(c), K. KOYAMA<fnref n="1" />(d), V. VENNIN<fnref n="3" />,<fnref n="4" />,<fnref n="1" />(e)
and D. WANDS<fnref n="1" />(f)

received 25 January 2023; accepted in final form 10 May 2023
published online 30 May 2023

**Abstract** – 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.

open access

Copyright © 2023 The author(s)

Published by the EPLA under the terms of the [Creative Commons Attribution 4.0 International License](https://creativecommons.org/licenses/by/4.0/)
(CC BY). Further distribution of this work must maintain attribution to the author(s) and the published
article's title, journal citation, and DOI.

**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–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_{\text{PBH}}$, the fraction of dark matter in the form of PBHs, i.e., $f_{\text{PBH}} = 1$ indicates that PBHs make up the entirety of the dark matter. For a recent review of constraints on $f_{\text{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–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_\zeta(k)$, which is observed on cosmic microwave background (CMB) scales to be $\sim 2.1 \times 10^{-9}$ [10]. It is commonly assumed that to form PBHs requires increasing $P_\zeta(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, $\zeta_G$ with higher-order parameters $f_{\text{NL}}$, etc. [12–17]. However, a non-perturbative calculation

<footnote n="1">¹ Institute of Cosmology & Gravitation, University of Portsmouth - Dennis Sciama Building, Burnaby Road, Portsmouth, PO1 3FX, UK</footnote>

<footnote n="2">² School of Mathematics and Physics, University of Portsmouth - Lion Gate Building, Lion Terrace, Portsmouth, PO1 3HF, UK</footnote>

<footnote n="3">³ Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, CNRS, Université PSL, Sorbonne Université, Université Paris Cité - 75005 Paris, France</footnote>

<footnote n="4">⁴ Laboratoire Astroparticule et Cosmologie, CNRS, Université Paris Cité - 75013 Paris, France</footnote>

<footnote>(a)E-mail: andrew.gow@port.ac.uk (corresponding author)</footnote>

<footnote>(b)E-mail: hooshyar.assadullahi@port.ac.uk</footnote>

<footnote>(c)E-mail: joseph.jackson@port.ac.uk</footnote>

<footnote>(d)E-mail: kazuya.koyama@port.ac.uk</footnote>

<footnote>(e)E-mail: vincent.vennin@ens.fr</footnote>

<footnote>(f)E-mail: david.wands@port.ac.uk</footnote>
can be made using the $\delta N$-method [18]. It has been argued that by fully incorporating quantum diffusion, in the so-called stochastic $\delta 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–30] (such heavy tails were also found in refs. [31–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 $\zeta$ is related to a Gaussian field $\zeta_G$ by what we will call a “generalised local” transformation.

**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[\zeta(\mathbf{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 top-hat window function [34,35]. The relationship between $\zeta$ and $C$ is non-linear, but the compaction during radiation domination can be written in terms of its linear part $C_\ell$ as

$$C = C_{\ell} - \frac{3}{8} C_{\ell}^{2}, \qquad (1)$$

with $C_\ell$ related to $\zeta$ by

$$C_{\ell} = -\frac{4}{3} r \zeta', \qquad (2)$$

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_\ell)$ instead of $P(C)$. Equation (1) has a maximum of $C = 2/3$ for $C_\ell = 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_\ell \le 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 $\zeta$ in the same region must also be spherical. In general, we can expand the field in spherical harmonics,

$$\zeta(\mathbf{x}) = \zeta_{00}(r) + \sum_{\ell>0, m \in [-\ell, \ell]} Y_{\ell m}(\theta, \phi) \zeta_{\ell m}(r). \qquad (3)$$

From now on we will consider only the spherically symmetric monopole term in this expansion and use the shorthand $\zeta \equiv \zeta(r) \equiv \zeta_{00}(r)$.

The reliance on $\zeta'$ rather than $\zeta$ itself means that $P(\zeta)$ is not sufficient to calculate the PBH mass distribution and abundance, and instead the joint probability $P[\zeta(r_1), \dots, \zeta(r_n)]$ is required [37]. Therefore, in our analysis, we will assume that $\zeta(r)$ can be related to a Gaussian field $\zeta_G(r)$ whose statistical properties are well understood. We introduce a “generalised local” transformation,

$$\zeta(r) = \zeta[\zeta_G(r), r], \qquad (4)$$

where we have allowed for an explicit $r$-dependence independent of the Gaussian field $\zeta_G$. This extra dependence is required to consider arbitrary $P(\zeta)$ since the Gaussian PDF will always carry an $r$ dependence through the variance $\Sigma_{YY}$ as well as $\zeta_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 $\zeta$ is a local function of a Gaussian field, $\zeta = \zeta[\zeta_G(r)]$. For a monotonic transformation, high peaks in $\zeta$ necessarily correspond to high peaks in $\zeta_G$, and hence this relation is consistent with our spherical assumption.

We want to connect an arbitrary non-Gaussian $\zeta$ with PDF $P(\zeta)$ to a Gaussian $\zeta_G$ with PDF

$$P_G(\zeta_G, \Sigma_{YY}) = \frac{1}{\sqrt{2\pi\Sigma_{YY}}} \exp\left(-\frac{\zeta_G^2}{2\Sigma_{YY}}\right), \qquad (5)$$

where we write the variance as $\Sigma_{YY}$ 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 $\zeta$ and $\zeta_G$,

$$F[\zeta(r), r] = F_G[\zeta_G(r), \Sigma_{YY}(r)], \qquad (6)$$

which can be applied for any $P(\zeta)$, whether analytical or numerical.

The simplest generalised local transformation has all the $r$-dependence coming from $\zeta_G(r)$ and the variance $\Sigma_{YY}(r)$, for which the non-Gaussian field can be written as

$$\zeta(\zeta_G, \Sigma_{YY}) = F^{-1}[F_G(\zeta_G, \Sigma_{YY})], \qquad (7)$$

where $F^{-1}$ is the inverse CDF for $\zeta$. In this case, $\zeta'$ in eq. (2) includes two terms, one from the dependence of $\zeta_G$ on $r$, and one from the dependence of $\Sigma_{YY}$ on $r$. We can therefore write the linear compaction (2) as

$$C_{\ell} = -\frac{4}{3} r [\mathcal{J}_1(\zeta_G)\zeta'_{G} + \mathcal{J}_2(\zeta_G)\Sigma'_{YY}], \qquad (8)$$

with Jacobian factors given by

$$\mathcal{J}_1(\zeta_G) = \frac{\partial \zeta}{\partial \zeta_G}, \qquad \mathcal{J}_2(\zeta_G) = \frac{\partial \zeta}{\partial \Sigma_{YY}}. \qquad (9)$$
Note that by construction, from eq. (7), we have $J_1 = P_G(\zeta_G)/P(\zeta)$. In the case of a local transformation where $\zeta$ depends only on $\zeta_G$, $J_2$ will be zero.

We can see from eq. (8) that the linear compaction depends on both $\zeta_G$ and $\zeta'_G$, and so $P(C_\ell)$ will be written in general as an integration over a 2D PDF,

$$P(C_{\ell}) = \int dX \int dY P(X, Y) \delta_{\text{D}}[C_{\ell} - C_{\ell}(X, Y)], \quad (10)$$

where $C_\ell(X, Y)$ is given by eq. (8), $\delta_D$ is a Dirac delta function, and X and Y are given by

$$X = r\zeta'_{G}, \quad Y = \zeta_{G}. \qquad (11)$$

The 2D Gaussian distribution, $P(X, Y)$, is given by

$$P(X, Y) = \frac{1}{2\pi\sqrt{\det(\Sigma)}} \exp\left(-\frac{V^T \Sigma^{-1} V}{2}\right), \quad (12)$$

with

$$V^T = (X, Y), \quad \Sigma = \begin{pmatrix} \Sigma_{XX} & \Sigma_{XY} \\ \Sigma_{XY} & \Sigma_{YY} \end{pmatrix}. \quad (13)$$

The correlators can then be written in terms of the power spectrum of the Gaussian field, $\mathcal{P}_{\zeta_G}(k)$,

$$\Sigma_{XX} = \int d(\ln k)(kr)^2 \left[ \frac{dj_0}{dz}(kr) \right]^2 \mathcal{P}_{\zeta_G}(k), \quad (14)$$

$$\Sigma_{XY} = \int d(\ln k)(kr) j_0(kr) \frac{dj_0}{dz}(kr) \mathcal{P}_{\zeta_G}(k), \quad (15)$$

$$\Sigma_{YY} = \int d(\ln k) j_0^2(kr) \mathcal{P}_{\zeta_G}(k), \quad (16)$$

where $j_0 = (\sin z)/z$ is the spherical Bessel function. Here we differ from ref. [38] by retaining the $\Sigma_{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 $\Sigma_{XY}^2 = \Sigma_{XX}\Sigma_{YY}$, 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,

$$\mathcal{P}_{\zeta_G}(k) = \frac{A}{\sqrt{2\pi\Delta}} \exp\left[-\frac{\ln^2(k/k_*)}{2\Delta^2}\right], \quad (17)$$

with amplitude A, width $\Delta$, and peak scale $k_*$. This is normalised such that it approaches the monochromatic case of a Dirac delta in $\ln k$ in the limit $\Delta \to 0$.

We can then rewrite the linear compaction in terms of X and Y as

$$C_{\ell} = -\frac{4}{3} [\mathcal{J}_1(Y)X + \mathcal{J}_2(Y)r\Sigma'_{YY}] \quad (18)$$

$$= -\frac{4}{3} [\mathcal{J}_1(Y)X + 2\mathcal{J}_2(Y)\Sigma_{XY}]. \quad (19)$$

Finally, carrying out the marginalisation in eq. (10) with $C_\ell(X, Y)$ given by eq. (19),

$$P(C_{\ell}) = \int d\zeta_{G} \frac{3}{4|\mathcal{J}_{1}(\zeta_{G})|} \times P\left\{-\frac{1}{\mathcal{J}_{1}(\zeta_{G})} \left[\frac{3}{4}C_{\ell} + 2\Sigma_{XY}\mathcal{J}_{2}(\zeta_{G})\right], \zeta_{G}\right\}. \quad (20)$$

The PBH mass m is related to the linear compaction $C_\ell$ by the critical collapse relation [40–42]

$$m = KM_H \left(C_{\ell} - \frac{3}{8}C_{\ell}^2 - C_c\right)^{\gamma}, \quad (21)$$

where K, $C_c$, and $\gamma$ are parameters that control the collapse, which in general are dependent on the shape of the perturbations [43–47]. The mass distribution is then given by [39]

$$f(m) = Z \frac{K (C_{\ell} - \frac{3}{8}C_{\ell}^2 - C_c)^{\gamma+1}}{\gamma (1 - \frac{3}{4}C_{\ell})} P(C_{\ell}), \quad (22)$$

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_\ell$ 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_\ell = 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

$$f_{\text{PBH}} = \int_{0}^{m_{\text{max}}} d(\ln m) f(m), \quad (23)$$

where $m_{\text{max}}$ corresponds to $C_\ell = 4/3$.

**Results.** – We now apply the method described in the previous section to two examples of non-Gaussian distributions, $P(\zeta)$. We choose the power spectrum width $\Delta = 0.3$ in eq. (17), and take the same peak scale as [39], $k_* = 1.56 \times 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$, $\gamma = 0.36$ (see footnote 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 \times 10^{16}$.

The first non-Gaussian distribution $P(\zeta)$ is an example of a local transformation arising from the classical $\delta N$ approach for an ultra slow roll (USR) phase in the early universe. This results in the transformation [38,39,49–51]

$$\zeta = -\frac{1}{3} \ln(1 - 3\zeta_G), \quad (24)$$

<footnote n="1">¹We note that the maximum scale and critical collapse parameters are valid for a monochromatic power spectrum, and will differ for our case of a lognormal with finite width [43–47]. However, the results of the non-Gaussian calculation presented here will not be strongly sensitive to these modifications.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1: Comparison of the piecewise exponential and the classical USR δN transformation. The slope, α = 3, is taken to match the classical USR δN tail. The two values of ζt are chosen to approximately match the classical USR δN P(ζ) tail (ζt = 1/3) and fPBH curve (ζt = 0.15). In the first three panels, A ~ 4 × 10⁻³, for which the classical USR δN distribution gives fPBH = 1. The blue shaded region in the top right panel represents the region relevant for PBH formation. The grey shaded region in the bottom right panel indicates overproduction of PBHs (fPBH > 1).</figcaption></figure>

and in the non-Gaussian PDF

$$P(\zeta) = \frac{1}{\sqrt{2\pi\Sigma_{YY}}} \exp \left[ -\frac{(1 - e^{-3\zeta})^2}{18\Sigma_{YY}} - 3\zeta \right], \quad (25)$$

where ΣYY 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,

$$P(\zeta) = N \begin{cases} P_G(\zeta), & \zeta \le \zeta_t, \\ P_G(\zeta_t) \exp[-\alpha(\zeta - \zeta_t)], & \zeta > \zeta_t. \end{cases} \quad (26)$$

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 fPBH 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 fPBH 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 fPBH. 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 > Mkpc. 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/4Cℓ) 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
<figure><img src="image_4.png" /><figcaption>Fig. 2: Comparison of $f_{\text{PBH}}(A)$ (top) and $f(m)$ (bottom) for different values of the piecewise distribution parameters $\zeta_t$ and $\alpha$, with $\alpha = 3$ (left) and $\zeta_t = 0.15$ (right). The grey shaded region in the top row indicates the overproduction of PBHs ($f_{\text{PBH}} > 1$). The mass distributions $f(m)$ in the bottom row are normalised to $f_{\text{PBH}} = 1$.</figcaption></figure>

decaying very slowly, so does not beat the diverging term. The mass distribution $f(m)$ arising from the classical USR $\delta N$ model does decay for $m > M_{k_*}$, despite having a far tail slope of $\alpha = 3$ in $P(\zeta)$ like the piecewise model. This is because, as stated previously, the relevant part of $P(\zeta)$ is not the far tail, but the transition, during which the effective value of $\alpha$ is significantly larger, leading to a decaying $P(C_\ell)$ 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_\ell = 4/3$, but is expected to break down before this point [52].

Figure 2 shows the effect of changing the piecewise transition point $\zeta_t$ and exponential slope $\alpha$ in eq. (26) on the total PBH abundance $f_{\text{PBH}}$ and the mass distribution $f(m)$. As in fig. 1, the largest values of $\zeta_t$ provide no enhancement over the Gaussian PDF except at the highest amplitudes, for which $f_{\text{PBH}} \gg 1$. For smaller values of $\zeta_t$, the exponential tail does provide a significant enhancement over the Gaussian case, greatly reducing the amplitude required for $f_{\text{PBH}} = 1$. In contrast, changing $\alpha$ only moves the $f_{\text{PBH}}$ curve by a (relatively) small amount, e.g., by factors of 10 rather than the $10^{10}$ coming from changing $\zeta_t$. The steeper values of $\alpha$ are suppressed compared to the shallower slopes, as expected. For large amplitudes, $f_{\text{PBH}}$ is suppressed (even compared to the Gaussian). This suppression comes from the normalisation of $P(\zeta)$ when including a large tail. The details of this suppression are unimportant, since they correspond to $f_{\text{PBH}} \gg 1$ and so

<table><caption>Table 1: Approximate values of $\nu_t = \zeta_t/\sqrt{\Sigma_{YY}}$ and $P(\zeta_t)$ required for a specific $f_{\text{PBH}}$, with a given $\alpha$. There is a linear relationship between $P(\zeta_t)$ and $f_{\text{PBH}}$ for constant $\alpha$, and an inverse linear relationship between $f_{\text{PBH}}$ and $\alpha$ for fixed $\nu_t$, which comes about from the integration of the far tail exponential behaviour.</caption><thead><tr><th>α</th><th>ν<sub>t</sub></th><th>P(ζ<sub>t</sub>)</th><th>f<sub>PBH</sub></th></tr></thead><tbody><tr><td rowspan="2">3</td><td>8.2</td><td>10<sup>-10</sup></td><td>1</td></tr><tr><td>10.7</td><td>10<sup>-20</sup></td><td>10<sup>-10</sup></td></tr><tr><td>0.3</td><td>8.2</td><td>10<sup>-10</sup></td><td>10</td></tr><tr><td>30</td><td>8.2</td><td>10<sup>-10</sup></td><td>10<sup>-1</sup></td></tr></tbody></table>

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 $\zeta_t$ and $\alpha$ on the mass distribution $f(m)$. For $\zeta_t$, there are three cases that are visible. The $\zeta_t = 0.3$ curve lies on top of the Gaussian, as for the $\zeta_t = 1/3$ case in fig. 1. For $\zeta_t = 0.25$, the function $f(m)$ starts to decrease for $m > M_{k_*}$ as $P(C_\ell)$ decreases, before growing again due to the diverging term discussed previously. Finally, all values of $\zeta_t$ smaller than 0.25 show no drop-off for $m > M_{k_*}$, instead being dominated by the diverging factor. Comparing the $f_{\text{PBH}}$ and $f(m)$ plots for varying $\zeta_t$, it seems clear that the amount of non-Gaussianity required to significantly reduce the power spectrum amplitude corresponding to
a fixed $f_{\text{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 $\mu$-distortions [53]. The right panel shows that, as with $f_{\text{PBH}}$, the mass distribution is less affected by $\alpha$ than $\zeta_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.

**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 $\zeta$ to the compaction function $C$. A non-perturbative method is essential when considering non-Gaussianities that are not described by the usual $f_{\text{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, $\zeta(\zeta_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 $\delta N$ relation for a local-type non-Gaussian $\zeta(\zeta_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 $\delta N$ model is due to the transition regime, rather than the far tail. Additionally, we have examined the transition value, $\zeta_t$, required for $f_{\text{PBH}} = 1$ using our piecewise model. We find that $\nu_t = \zeta_t/\sqrt{\sum_{YY}}$ is approximately independent of the peak amplitude of the Gaussian power spectrum, $A$, for fixed $f_{\text{PBH}}$ and $\alpha$. We also see a linear relationship between $P(\zeta_t)$ and $f_{\text{PBH}}$, and when changing the exponential slope $\alpha$ we find an inverse linear relation between $f_{\text{PBH}}$ and $\alpha$, see table 1. 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_{\text{PBH}}$, e.g., to evade the $\mu$-distortion constraints for supermassive PBHs.

***

The authors would like to thank ILIA MUSCO and SHI PI for helpful discussions. This work is supported by the Science and Technology Facilities Council (grant Nos. ST/S000550/1, ST/W001225/1, and ST/T506345/1). For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising. Supporting research data are available on reasonable request from the corresponding author, ADG.

*Additional remark:* During the writing up of this paper, ref. [54] appeared on the arXiv. The authors utilise the

same mechanism of writing $\zeta$ in terms of $\zeta_G$ and applying the product rule to the derivative in the compaction function, but restrict themselves to local transformations $\zeta(\zeta_G)$, rather than allowing $P(\zeta)$ to be arbitrary.

*Data availability statement:* The data that support the findings of this study are available upon reasonable request from the authors.

## REFERENCES

[1] ZEL'DOVICH Ya. B. and NOVIKOV I. D., Sov. Astron., **10** (1967) 602.

[2] HAWKING S., Mon. Not. R. Astron. Soc., **152** (1971) 75.

[3] CARR B. J. and HAWKING S. W., Mon. Not. R. Astron. Soc., **168** (1974) 399.

[4] CARR B., KOHRI K., SENDOUDA Y. and YOKOYAMA J., Rep. Prog. Phys., **84** (2021) 116902.

[5] BIRD S., CHOLIS I., MUÑOZ J. B., ALI-HAÏMOUD Y., KAMIONKOWSKI M., KOVETZ E. D., RACCANELLI A. and RIESS A. G., Phys. Rev. Lett., **116** (2016) 201301.

[6] SASAKI M., SUYAMA T., TANAKA T. and YOKOYAMA S., Phys. Rev. Lett., **117** (2016) 061101.

[7] CLESSE S. and GARCÍA-BELLIDO J., Phys. Dark Univ., **15** (2017) 142.

[8] GOW A. D., BYRNES C. T., HALL A. and PEACOCK J. A., JCAP, **2020** (2020) 031.

[9] FRANCIOLINI G., BAIBHAV V., DE LUCA V., NG K. Y., WONG K. W. K., BERTI E., PANI P., RIOTTO A. and VITALE S., Phys. Rev. D, **105** (2022) 083526.

[10] AGHANIM N. et al., Astron. Astrophys., **641** (2020) A6.

[11] GOW A. D., BYRNES C. T., COLE P. S. and YOUNG S., JCAP, **2021** (2021) 002.

[12] BULLOCK J. S. and PRIMACK J. R., Phys. Rev. D, **55** (1997) 7423.

[13] YOUNG S., REGAN D. and BYRNES C. T., JCAP, **2016** (2016) 029.

[14] YOO C.-M., GONG J.-O. and YOKOYAMA S., JCAP, **2019** (2019) 033.

[15] PALMA G. A., SCHEIHING HITSCHFELD B. and SYPSAS S., JCAP, **2020** (2020) 027.

[16] TAOSO M. and URBANO A., JCAP, **2021** (2021) 016.

[17] YOUNG S., JCAP, **2022** (2022) 037.

[18] LYTH D. H. and RODRÍGUEZ Y., Phys. Rev. Lett., **95** (2005) 121302.

[19] FUJITA T., KAWASAKI M., TADA Y. and TAKESAKO T., JCAP, **2013** (2013) 036.

[20] VENNIN V. and STAROBINSKY A. A., Eur. Phys. J. C, **75** (2015) 413.

[21] PATTISON C., VENNIN V., ASSADULLAHI H. and WANDS D., JCAP, **2017** (2017) 046.

[22] EZQUIGA J. M., GARCÍA-BELLIDO J. and VENNIN V., JCAP, **2020** (2020) 029.

[23] VENNIN V., Stochastic inflation and primordial black holes, arXiv:2009.08715 [astro-ph.CO] (2020).

[24] FIGUEROA D. G., RAATIKAINEN S., RÄSÄNEN S. and TOMBERG E., Phys. Rev. Lett., **127** (2021) 101302.

[25] ANDO K. and VENNIN V., JCAP, **2021** (2021) 057.

[26] PATTISON C., VENNIN V., WANDS D. and ASSADULLAHI H., JCAP, **2021** (2021) 080.
[27] RIGOPOULOS G. and WILKINS A., JCAP, **2021** (2021) 027.

[28] TADA Y. and VENNIN V., JCAP, **2022** (2022) 021.

[29] JACKSON J. H. P., ASSADULLAHI H., KOYAMA K., VENNIN V. and WANDS D., JCAP, **2022** (2022) 067.

[30] ANIMALI C. and VENNIN V., *Primordial black holes from stochastic tunneling*, arXiv:2210.03812 (2022).

[31] PANAGOPOULOS G. and SILVERSTEIN E., *Primordial Black Holes from non-Gaussian tails*, arXiv:1906.02827 (2019).

[32] ACHÚCARRO A., CÉSPedes S., DAVIS A.-C. and PALMA G. A., *J. High Energy Phys.*, **2022** (2022) 052.

[33] CAI Y.-F., MA X.-H., SASAKI M., WANG D.-G. and ZHOU Z., JCAP, **2022** (2022) 034.

[34] SHIBATA M. and SASAKI M., Phys. Rev. D, **60** (1999) 084002.

[35] HARADA T., YOO C.-M., NAKAMA T. and KOGA Y., Phys. Rev. D, **91** (2015) 084057.

[36] KOPP M., HOFMANN S. and WELLER J., Phys. Rev. D, **83** (2011) 124025.

[37] DE LUCA V. and RIOTTO A., Phys. Lett. B, **828** (2022) 137035.

[38] BIAGETTI M., DE LUCA V., FRANCIOLINI G., KEHAGIAS A. and RIOTTO A., Phys. Lett. B, **820** (2021) 136602.

[39] KITAJIMA N., TADA Y., YOKOYAMA S. and YOO C.-M., JCAP, **2021** (2021) 053.

[40] CHOPTUIK M. W., Phys. Rev. Lett., **70** (1993) 9.

[41] EVANS C. R. and COLEMAN J. S., Phys. Rev. Lett., **72** (1994) 1782.

[42] NIEMEYER J. C. and JEDAMZIK K., Phys. Rev. Lett., **80** (1998) 5481.

[43] MUSCO I., Phys. Rev. D, **100** (2019) 123524.

[44] YOUNG S., Int. J. Mod. Phys. D, **29** (2019) 2030002.

[45] GERMANI C. and MUSCO I., Phys. Rev. Lett., **122** (2019) 141302.

[46] GERMANI C. and SHETH R. K., Phys. Rev. D, **101** (2020) 063520.

[47] ESCRIVÀ A., GERMANI C. and SHETH R. K., JCAP, **2021** (2021) 030.

[48] MUSCO I., DE LUCA V., FRANCIOLINI G. and RIOTTO A., Phys. Rev. D, **103** (2021) 063538.

[49] CAI Y.-F., CHEN X., NAMJOO M. H., SASAKI M., WANG D.-G. and WANG Z., JCAP, **2018** (2018) 012.

[50] ATAL V., GARRIGA J. and MARCOS-CABALLERO A., JCAP, **2019** (2019) 073.

[51] PI S. and SASAKI M., *Primordial Black Hole Formation in Non-Minimal Curvaton Scenario*, arXiv:2112.12680 (2021).

[52] ESCRIVÀ A., Phys. Dark Univ., **27** (2020) 100466.

[53] CHLUBA J., ERICKCEK A. L. and BEN-DAYAN I., *Astro-ph. J.*, **758** (2012) 76.

[54] FERRANTE G., FRANCIOLINI G., IOVINO A. J. and URBANO A., Phys. Rev. D, **107** (2023) 043520.