PBH in single field inflation: the effect of shape dispersion and non-Gaussianities

Primordial black holes (PBHs) may result from high peaks in a random field of cosmological perturbations. In single field inflationary models, such perturbations can be seeded as the inflaton overshoots a small barrier on its way down the potential. PBHs are then produced through two distinct mechanisms, during the radiation era. The first one is the familiar collapse of large adiabatic overdensities. The second one is the collapse induced by relic bubbles where the inflaton field is trapped in a false vacuum, due to large backward fluctuations which prevented horizon sized regions from overshooting the barrier. We consider (numerically and analytically) the effect of non-Gaussianities on the threshold for overdensities to collapse into a PBH. Since typical high peaks have some dispersion in their shape or profile, we also consider the effect of such dispersion on the corresponding threshold for collapse. With these results we estimate the most likely channel for PBH production as a function of the non-Gaussianity parameter $f_{\rm NL}$. We also compare the threshold for collapse coming from the perturbative versus the non perturbative template for the non-Gaussianity arising in this model. We show that i) for $f_{\rm NL}\gtrsim 3.5$, the population of PBH coming from false vacuum regions dominates over that which comes from the collapse of large adiabatic overdensities, ii) the non-perturbative template of the non-Gaussianities is important to get accurate results. iii) the effect of the dispersion is small in determining the threshold for the compaction function, although it is appreciable in determining the threshold amplitude for the curvature perturbation at low $f_{\rm NL}$. We also confirm that the volume averaged compaction function provides a very accurate universal estimator for the threshold.


Introduction
Primordial Black Holes (PBH) may have formed during the radiation dominated era due to unusually high peaks in the distribution of cosmological density perturbations [1,2]. There are strong observational constraints on the abundance of PBH over a wide range of mass scales [3]. Nonetheless, these still allow for several phenomenologically interesting possibilities. For instance, PBH of sublunar [4] or stellar mass [5,6] may constitute a sizable fraction of all dark matter in the universe 1 . Also, the origin of supermassive black holes at the center of galaxies is not very well understood at present, and one possibility is that they may have formed by accretion from a smaller intermediate mass PBH seed (for a recent review, see [9]).
In order to make accurate predictions on the statistical properties of PBHs, it is necessary to be specific about their formation process. One of the simplest mechanisms is the collapse of large adiabatic perturbations seeded during a period of single-field inflation [10][11][12][13][14][15][16]. While fluctuations must be predominantly Gaussian at the cosmic microwave background scales, those leading to PBH formation at smaller scales are typically non-Gaussian [17][18][19] 2 . Sufficiently large amplification of the perturbations are induced while the inflation overshoots a small barrier on its way down the potential, undergoing a periond of "constant roll" (see Fig. 1). In this context, PBH can be formed not only from the collapse of a large adiabatic overdensity, but also from false vacuum bubbles which continue inflating in the ambient radiation dominated universe, and eventually pinch off from it. This results in a black hole which separates the ambient universe from an inflating baby universe [18,25,26]. 3 A question of practical interest is to determine the abundance of PBHs. Several works have already treated the influence of non-Gaussianities in the abundance of PBHs [18,[27][28][29][30][31][32][33][34][35][36]. Since this turns to be large, it is important to i) predict the amplitude and shape of the non-Gaussianities for a given model of PBH formation, and ii) consider their influence beyond perturbation theory.
When PBHs are formed from rare overdensities, their abundance will depend on the threshold for the amplitude of the overdensity to collapse once it reenters the horizon. This threshold notoriously depends on the shape (or profile) for the overdensity [37][38][39][40][41]. For a Gaussian random field, the typical shape of high peaks is determined from the power spectrum, but if the distribution is non-Gaussian, the shape will also depend on the nature of the non-Gaussianity [18,34,35]. Furthermore, since fluctuations are drawn from a statistical distribution, the shapes of perturbations susceptible of collapsing will inherit a dispersion. While the mean profile is usually taken to be representative of the typical shape, it seems important to consider how the threshold may vary due to the dispersion of shapes. This point is particularly relevant when a mean profile for the perturbations cannot be defined, as it is the case for large overdensities coming from the model of single-field inflation with a barrier 4 [18].
In this work we study the dependence of the threshold on the dispersion of the profiles, including the non-Gaussianity resulting from the physics of single-field inflation. The non-Gaussianity is entirely due to the non-linear relation between the Gaussian variable and the non-Gaussian gauge-invariant curvature perturbation ζ. Here δφ is the inflaton field perturbation in the flat slicing, evaluated at the onset of the slow roll attractor behaviour past the top of the barrier, and H is the expansion rate during inflation 5 . For the nonlinear relation between ζ and ζ g , we will compare the non-perturbative expression which follows from the single field model where the inflaton overshoots a small barrier [18], with the more widely used perturbative Taylor expansion of ζ to second order in ζ g (parametrized by the standard coefficient f NL ). These non-perturbative and perturbative versions of local non-Gaussianity are given, respectively, in Eqs. (2.13) and (2.15) below. We will find the thresholds for collapse into a PBH under the assumption of spherical symmetry, by using a recently developed numerical code [45]. This solves the Misner-Sharp (MS) partial differential equations by using spectral methods. We will also compare the results obtained by numerical evolution with the results which can be obtained from a recently 3 These are sometimes refered to as black holes with a baby universe inside. Note, however, that the baby universe is not in the trapped region, or "interior" of the black hole. Rather, the trapped region separates two normal regions, one in the parent ambient universe and the other in the baby universe, which were once causally connected but are not anymore, after the trapped region forms. 4 In a nutshell, the problems is that ζ diverges when the amplitude of ζg reaches a critical value µ * , and it is not even defined for larger amplitude of ζg, for which there is a finite probability. 5 Refs. [42][43][44] consider the non-Gaussianity in the density perturbation δ due to the non-linear relation between δ and ζ. Note that such discussion would be redundant in our approach, where the initial conditions for numerical evolution, as well as the threshold estimators for gravitational collapse, are expressed directly in terms of ζ.
proposed universal estimator for the strength of a perturbation [46]. This is given by a suitable spatial averageC of the so-called initial compaction function C[ζ(r)] [37], out to a certain optimal radius r m . The threshold value forC which triggers gravitational collapse turns out to be extremely robust, in the sense that it is nearly independent of the radial profile of the perturbation ζ(r).
The plan of the paper is the following: In section 2, we consider the typical shapes of a high peak in the curvature perturbation profile, within one standard deviation of the median profile, and we introduce the non-perturbative relation between the curvature perturbation ζ and the Gaussian variable ζ g . In Section 3, we present the Misner-Sharp equations and we review the criteria for the formation of PBHs. The results of the numerical simulation and the analytical estimates, together with their interpretation are presented in section 4.
2 Large and rare peaks from single field inflation At cosmological scales, the power spectrum of primordial perturbations must be of the order of 10 −9 , in accordance with observations of the cosmic microwave background. However, in order for PBH formation to be significant, the power must be of the order of 10 −3 − 10 −2 at the PBH scale. This jump in the amplitude can be achieved if the inflaton field passes trough a period of constant-roll. 6 Parametrically, the fraction of dark matter in PBH is Ω P BH ∼ 10 9 (M /M P BH ) 1/2 β 0 , where the probability of PBH formation at the time when a large perturbation crosses the horizon can be roughly estimated as β 0 ∼ exp[−ζ 2 th /2σ 2 0 ], for some threshold value ζ th ∼ 1. The remaining factors in the estimate of Ω P BH account for the dilution of radiation relative to PBH, from the time of their formation until the time t eq . For M P BH in the broad range 10 −13 − 10 2 M , the threshold for the perturbations to undergo gravitational collapse must be in the range ζ th ∼ (6 − 8)σ 0 , sizably larger than the standard deviation, in order to obtain a significant Ω P BH ∼ 1. Because these perturbations are very rare, we can use the theory of high peaks to describe them.

The typical high peak profiles
Since the non-Gaussian curvature perturbation ζ is a local function of the Gaussian field ζ g , let us start by reviewing the latter [47]. This will be the basis to describe the non-Gaussian realisations. Fluctuations of ζ g are characterized by the power spectrum P ζ (k), representing the variance of the random field per logarithmic interval in k, Introducing the normalized two point correlation function of ζ g ( x) as peaks of the Gaussian random field of given amplitude µ = νσ 0 at the origin, have a mean profile given by where the last term can be neglected in the limit of high peaks ν 1. Note that ψ(0) = 1. The above expectation is calculated by using the number density distribution of peaks. This distribution is almost Gaussian, except for a Jacobian prefactor which relates the condition of being an extremum with the condition for the peak to be at a certain location. If we simply condition the field value to be at a certain height, the distribution is Gaussian, an leads to the simpler expression [47] ζ g (r)|ν = σ 0 νψ(r), (2.4) which coincides with the large ν limit of (2.3). For a Gaussian distribution, the mean and the median coincide, and therefore for the rest of this paper we shall refer to (2.4) as the median Gaussian profile.
Still, there will be some deviations around the median, so that the typical profile will be of the form ζ g (r) = µψ(r) ± ∆ζ, (2.5) where the variance of the shape is given by [47] (∆ζ(r)) 2 where the gradient moments of the power spectrum are given by In what follows, we are going to consider two different forms for the enhancement of the power spectrum at the PBH scale. Monochromatic power spectrum: This is simply an idealized a delta function enhancement, such that the power spectrum is given by In this case the median shape is given by ζ g (r) = µ sinc(k 0 r), (2.9) while the dispersion takes the following form (∆ζ(r)) 2 Note that in this case, we have γ = 1, and the general expression (2.6) contains indeterminate ratios. In order to obtain (2.10) we have regularized the delta function by using a normalized distribution which is constant in an interval of radius ε around k = k 0 , and vanishes outside of this interval, taking the limit ε → 0 at the end. Sharply peaked power spectrum: We are also going to consider a more realistic case, in which the enhancement follows a power law growth k n . Models of the type considered here tend to have n in the range 3 − 4, and for definiteness we shall consider n = 4, which seems to be the maximum possible value in canonical single field scenarios [48]. This upper bound is actually saturated in the example considered in [18]. We also consider a rapid fall of the power spectrum after the peak. In the single field model, such fall-off behaves as k − 12 5 f NL [18]. Note that at low f NL < 5/3, the fall-off is not sharp enough to make σ 2 [given in (2.7)] indepent on the ultraviolet details of the spectrum. In other words, peak theory cannot be blindly used in this case to find the number density of PBH at the scale of the peak k p , because the distribution of peaks is dominated by smaller scales. In what follows, we will simply introduce a sharp cut-off at the peak value k p . This amounts to a top hat window function in momentum space, which filters out the smaller scales. 7 The spectrum is then given by (2.11) In this case, the correlation function determining the shape of the peak is given by where we have further assumed that k 0 k p . For its dispersion we can take directly Eq. (2.6), since in this case γ = 1. We now discuss the effect of non-Gaussianities.

Non-Gaussianity
In single-field inflation, when the inflaton passes through a period of constant-roll as it overshoots a barrier, the non-Gaussian curvature perturbation ζ is related to the Gaussian field ζ g defined in (1.1) as [18] (2.13) The parameter µ * can be written as a function of the potential as with η ≡ V /V , evaluated at the local maximum of the barrier. The relation (2.13) is only defined for perturbations with ζ g < µ * . Perturbations with ζ g > µ * are so large that they prevent the inflaton field from overshooting the local maximum [18]. The regions where the inflaton is trapped in the false vacuum are localized false vacuum bubbles which, from the point of view outside observers, end up forming a black hole, while from the point of view of internal observers they continue inflating. That is the reason why such PBH are said to carry a baby universe inside [25,26]. In this context, black holes can be formed in two different ways. If ζ g is larger than a certain threshold µ th , whith µ th < ζ g < µ * , then standard black holes will be created by the gravitational collapse of the adiabatic overdensity. On the other hand, regions where ζ g > µ * , will lead to false vacuum bubbles. 8 considered in reference [30], by using a slightly di↵erent (although essentially equivalent) method for generating a dispersion of shapes [18,31]. In [30], however, a fiducial value C th ⇡ 0.268 was used throughout, independently of the value of f NL . Here, we show that C th is actually a decreasing function of f NL , which contributes to enhancing the abundance of PBH with growing f NL . The total e↵ect can be inferred directly from the left panel of Fig. 2, where we see that µ th changes by roughly a factor of 2 as f NL varies from 0 to 6. Since the abundance of PBHs is exponential in µ 2 / 2 0 , it follows that for the larger values of non-Gaussianity the power spectrum can be about a factor of 4 or so smaller than it is for the Gaussian case.

Case B: Non-perturbative template
In the single-field case, we note that the threshold for collapse approaches the critical value for jumping to the false minimum for f NL ⇠ 3. Around this value of f NL , this means that black holes would rather be formed though to the creation of false vacuum regions. We -11 -non-Gaussianity parameter f NL , from µ th ⇡ 0.6 for f NL = 0 to µ th ⇡ 0.33 for f NL = 6.
Recently, the e↵ect of non-Gaussianity with the quadratic template (2.14) was also considered in reference [30], by using a slightly di↵erent (although essentially equivalent) method for generating a dispersion of shapes [18,31]. In [30], however, a fiducial value C th ⇡ 0.268 was used throughout, independently of the value of f NL . Here, we show that C th is actually a decreasing function of f NL , which contributes to enhancing the abundance of PBH with growing f NL . The total e↵ect can be inferred directly from the left panel of Fig. 2, where we see that µ th changes by roughly a factor of 2 as f NL varies from 0 to 6. Since the abundance of PBHs is exponential in µ 2 / 2 0 , it follows that for the larger values of non-Gaussianity the power spectrum can be about a factor of 4 or so smaller than it is for the Gaussian case.

Case B: Non-perturbative template
An inflaton potential with a small barrier on its slope. As the background field goes over the barrier, it undergoes a period of constant-roll, which strongly amplifies the power spectrum of adiabatic perturbations to the amplitude required for significant PBH production. Large backward fluctuations may also prevent some horizon sized regions from overshooting the barrier, generating false vacuum bubble relics. From the point of view of internal observers, these continue inflating at a high rate, while from the external point of view, these bubbles will form PBH once they enter the horizon during the radiation dominated era.
By Taylor expanding the non-perturbatuve relation (2.13) to quadratic order in ζ g , we obtain the widely used perturbative template of local type non-Gaussianity, The parameter µ * is related to f NL through It is clear, however, that this truncated expansion is far from accurate, since PBH formation occurs in the regime where ζ g is not small. Furthermore, the perturbative template does not capture the existence of a second channel for PBH production from regions trapped in a false vacuum, since ζ in Eq.(2.15) is well defined for any amplitude of ζ g . Nonetheless, because of its prevalence in the literature, and in order to compare with other approaches, it seems of some interest to also consider this quadratic template. Hence, in the following we will consider two different cases. Case A corresponds to an idealized Dirac delta function power spectrum for the Gaussian variable ζ g , as in Eq. (2.8), where we will consider the "vanilla" perturbative local template (2.15) in order to obtain the non-Gaussian curvature perturbation. Case B is a more realistic scenario based on the single field model of [18], where the logarithmic template for non-Gaussianity will be combined with the power spectrum (2.11) in order to determine the range of typical shapes for ζ.

The formation of PBH
In this section we describe the relevant equations for the evolution of spherically symmetric perturbations, which can be solved with the help of the numerical code recently developed in [45]. We also describe the criteria for creation of a BH.

The Misner-Sharp equations
The Misner-Sharp equations (MS) are the Einstein's equations for a spherically symmetric spacetime, in a frame comoving with a perfect fluid [52]. In this case, the metric can be written in the general diagonal form where dΩ = dθ 2 + sin 2 θdφ 2 is the metric on the unit 2-sphere. The fluid is at rest relative to the radial coordinate r. The metric components A(r, t), B(r, t) and R(r, t) are all positive, the latter one corresponding to the areal radius of the 2-spheres. Following [52], we define the partial derivatives with respect to proper time and proper distance as Applying the last two operators to R we can also define We may now introduce the Misner-Sharp mass M (r, t) through the equation We will also use the form of the stress energy tensor of a perfect fluid where u µ = (A −1 , 0, 0, 0) is the four-velocity field of the fluid, ρ is the energy density, and p is the pressure. During the radiation era, the equation of state is p = 1 3 ρ. In terms of these variables, the MS equations take the form Let us now discuss the initial conditions for evolution in terms of the random field ζ(r) of primordial curvature perturbations.

The long wavelength approximation and initial conditions
Initially, at early times, perturbations have a physical wavelength L much larger than the Hubble radius H −1 [37]. Hence, we are going to consider the long wavelength approximation to determine the form of our initial metric and hydrodynamical variables. This is based in expanding the exact solutions in a power series of a parameter , (3.13) to the lowest non-vanishing order in (t) 1. In the limit → 0, the metric of a perturbed FRW model can be written in the form (3.14) This is in a coordinate system where the energy density of the fluid is used as a clock, so that t = const. surfaces coincide with ρ = const. surfaces. He have also restricted to spherical symmetry, which excludes the presence of tensor modes (gravitational waves). In terms of ζ(r), the long wavelength solution of the MS equations reads [41] U = H(t)R(1 + 2Ũ ), where the functionsρ,Ũ represent the energy density and velocity perturbation, given bỹ Here r k is the comoving lengthscale of the perturbation associated to the wavenumber k, i.e. r k e ζ(r k ) = [H(t)a(t) ] −1 . As mentioned below Eq. (3.11), all remaining variables can be obtained from the set (U, ρ, M, R), and so it is not necessary to specify any additional initial conditions.

The criterion for BH production
The formation of a black hole for a given initial condition can be inferred from the behaviour of perturbations which do not dissipate after entering the horizon but continue growing until a trapped surface [53] is formed. This signals the onset of gravitational collapse. To identify the trapped surfaces, we consider the expansion Θ ± ≡ h µν ∇ µ k ± ν of null geodesic congruences k ± µ orthogonal to a spherical surface Σ. Here h µν is the metric induced on Σ. There are two such congruences, which we may call inward and outward directed, with components k ± µ = (A, ±B, 0, 0), such that k + · k − = −2. In flat space, Θ − < 0, while Θ + > 0. Surfaces Σ with this property are called "normal". If both expansions are negative, the surface is called "trapped", while if both are positive, the surface is "anti-trapped". In terms of the MS variables [54], we have In a spherically symmetric spacetime, any point in the (r, t) plane can be thought of as a closed surface Σ of proper radius R(r, t). We can classify such points into normal, trapped and anti-trapped. In the transition from a normal region to a trapped region, we must go through a boundary where Θ − < 0 and Θ + = 0. This is a marginally trapped surface which is usually called the apparent horizon. Since Θ + Θ − ∝ U 2 − Γ 2 = 0 and using Eq. (3.5), the condition for the formation of an apparent horizon is simply This could be marginally trapped, as it occurs for black holes, or marginally anti-trapped, as is the case for a cosmological horizon. If the condition R < 2M is satisfied in the vicinity of the apparent horizon, this means that we have trapped surfaces, and a PBH will be formed in the subsequent evolution. A useful estimator for the strength of a spherically symmetric perturbation is the socalled compaction function, which is the mass excess δM (R) = M −M b enclosed in the aereal radius R(r, t) relative to the FLRW background M b , divided by the areal radius 9 [37] C(r, t) ≡ δM R .
(3.17) From (3.11), the MS mass out to the aereal radius R is given by and a similar expression for M b (R), where ρ is replaced by the background density ρ b . Here we have introduced the volume element 10 where dR = R dr is evaluated on t = const. hypersurfaces. In the last step whe consider the long wavelength limit, which is valid at sufficiently early times, when ζ(r) is time independent 9 Here we use the definition of the compaction function given originally in [37], which has also been used in most of the subsequent literature. Note, however, that some recent papers use a convention which differs by a factor of 2.
10 Note that this differs from the proper volume element on t = const. hypersurfaces dVMS = ΓdVproper. and R ≈ a(t)re ζ(r) . It is straigthforward to check that the compaction function can also be expressed as where we have introduced the "volume" averaged density perturbation Let us note that he condition R < 2M for the formation of a trapped surface is related the criterion C max ≈ 1/2 suggested in [45], where C max (r, t) is the maximum value of the compaction function at some given moment of time. Note that, indeed, for C max ≥ 1/2 we have 2M/R = 2C max +2M b /R > 1, guaranteeing that the surface is trapped. From a practical point of view, both criteria perform with similar efficiency in the simulations we have run.

A universal threshold for collapse
In the long wavelength limit, the compaction function is time independent, and can be expressed in terms of the curvature perturbation as [38] C(r) = 1 3 (1 − (1 + rζ (r)) 2 ), (3.22) where the prime denotes the partial derivative with respect to the radial coordinate. It has long been recognized that the initial compaction function C(r) is a useful tool for predicting whether a pertrubation will end up collapsing into a PBH. If C(r) has a maximum at r = r m which satisfies C(r m ) > C th , then a PBH will be formed after R(r m , t) enters the horizon [37]. An interesting feature of the threshold value C th , is that its possible range is rather limited. Indeed, from Eq. (3.22) it is easy to see that C th cannot be larger than 1/3. Also, it was argued in [38] that on physical grounds C th 0.21. The precise value of this lower bound is not very tightly determined by the argument, but recently it has been shown by numerical studies that slightly lower values are possible [46], and that the threshold lies in the range for a broad class of shapes. Since this window is relatively narrow, spanning less than a factor of 2, the use of a threshold C th has been popular in phenomenological studies of PBH production. Nonetheless, the precise value of C th still depends on the profile of the perturbation, and so in this approach we cannot completely dispense with numerical simulations of collapse in order to obtain accurate results.
Remarkably, a universal estimator has recently been proposed, whose threshold value for PBH formation seems to be independent on the shape of the high peak overdensity [46]. This consists of a spatial average of the compaction function out to the optimal radius r m corresponding to the maximum of C, where again, the MS volume element (3.19) is used for averaging 11C where R m = R(r m , t). The shape-independent threshold value for gravitational collapse is then given byC th ≈ 1/5. (3.25) This universal behaviour has been tested in [46] for a very broad class of shapes. Here we shall further confirm its validity by checking that it holds to very good accuracy in the class of profiles that we will study.

Results
Here we consider the set of typical profiles corresponding to Case A and Case B described at the end of Subsection 2.2 (with their corresponding dispersion in shapes). We determine the corresponding thresholds for collapse µ th as a function of the non-Gaussianity parameter, in the range 0 < f NL < 6, and for the family of profiles in the range ζ g = µψ±∆ζ, corresponding to one standard deviation from the median profile. We have done so by using two different methods. Namely, by numerical evolution with the code developed in [45], and by using the universal criterion based onC = 1/5. In Fig. 3 we show the results for the thresholds µ th and C th evaluated from these two methods, in the case where we do not include the dispersion ∆ζ = 0, for the perturbative and non perturbative template. We see a good agreement between both, within a deviation of ∼ 2%, as was reported in [46]. Let us now comment on the more qualitative features of the results and their physical implications.

Case A: Perturbative template
This case corresponds to the Dirac delta function power spectrum (2.8), together with the perturbative local template (2.15) for the relation between ζ g and the curvature perturbation ζ.
In Fig. 2 we display the time evolution of the mean profile (2.9) for the Gaussian case (f NL = 0). The "sinc" profile (2.9) is somewhat peculiar, in that the initial compaction function (represented as a blue line in the figure) has a dominant peak at r = r m ≈ 2.7k −1 0 , and then an infinite number of nearly equally spaced secondary peaks of nearly equal height at r r m . The threshold for gravitational collapse of the dominant peak once it enters the horizon is determined numerically to be C th ≈ 0.29. This raises the somewhat naive question of what happens to the secondary peaks if the compaction function exceeds C th also at the secondary peaks. Will these also trigger the gravitational collapse of bigger PBHs once they enter the horizon? It is clear from the figure that this will not be the case. As soon as the dominant peak enters the horizon, at the time t H , the width of the secondary peaks will also be within the horizon, and we see that these secondary structures disipate due to pressure gradients. 12 By contrast, the dominant peak continues to grow and in a time-scale t ∼ 10t H , it reaches C > 1/2, signaling the existence of a trapped region with 2M > R.
In fact, once we consider the dispersion of profiles, ζ ± g = µψ ± ∆ζ, with ∆ζ g given in (2.10), we find that the initial compaction function for ζ + g can be lower at the first peak than it is at the subsequent "secondary" ones. Still, the first peak is the one that grows under Figure 2. Time evolution of the compaction function C(r, t) for the Gaussian profile (2.9), with µ = 0.64, slightly larger than the threshold value µ th ≈ 0.61. For reference, the threshold value C th is indicated as a dashed line. The radial coordinate is in units of the initial time t i , which we take to be much smaller than the time t H at which r m crosses the horizon, t H = 100 t i . The size of the grid is actually somewhat larger than displayed, with r max = 200 t i , much larger than the initial Hubble radius H −1 i = 2t i . After the time t H the secondary peaks in the compaction function dissipate due to pressure gradients. The dominant peak, on the other hand, continues to grow. By the time t = 16t H , the compaction function has reached values significantly larger than 1/2, indicating that a trapped region has already formed. time evolution, until a trapped surface forms, whereas the secondary ones dissipate, This is important, because it highlights the fact that the relevant optimal radius r m at which we evaluate the compaction function in order to determine the threshold, and which also enters the universal estimator Eq. (3.24), is not the absolute maximum of the compaction function, but the local maximum which is closest to the origin.
We have determined the threshold amplitude µ th and the threshold compaction function C th for different values of the non-Gaussianity parameter in the range 0 ≤ f NL ≤ 6. The numerical results are shown in Fig. 3 and Fig. 4.
In particular, we find that the threshold for the compaction function reaches a constant as we increase the non-linear parameter f NL . To gain some insight into the origin of this behaviour, let us note that at sufficiently large f NL the overdensity is dominated by the non-linear term. Indeed, for µf NL 1/ψ(r m ) ≈ 1.85, (4.1) the median shape can be approximated as ζ(r) ≈ f NL µ 2 ψ 2 (r) out to the radius r m where in the last step in (4.1) we use ψ(r) = sinc(r), and r m ≈ 1.8 is the maximum of the compaction function for the profile ψ 2 (r). In this regime, the shape of the perturbation is independent of f NL , and hence, we expect C th to be independent of f NL : Here the numerical value is calculated by evolving the profile ζ ∝ sinc 2 (k 0 r). From the right panel in Fig. 3 we see that C th is indeed nearly constant for f NL in the range from 2 to 6. This is, however, somewhat coincidental, since the condition µ th f NL 2 is only satisfied for f NL 10. In the same regime, from r m ζ (r m ; µ th )) = √ 1 − 3C th − 1, we expect Note that this overestimates the actual values of µ th in the interval 1 < f NL < 6, by 30% or so (See Fig. 3). The reason is that for f NL 10, the value of r m and, more importantly ψ(r m )ψ (r m ), changes appreciably with f NL .
By contrast with C th , we find that the threshold amplitude µ th decreases quite significantly with f NL in the 0 < f NL < 6 interval. We also note that the dispersion of the shapes accounts for a very small dispersion of C th . On the other hand, the threshold for the amplitude µ th has a larger variability, in particular at low f NL .
Qualitatively, this behaviour is consistent with the recent work by [35], for the same monochromatic power spectrum. It should be noted, however, that there are some differences in the approach, and that we are considering a different set of profiles. Ref. [35] develops a perturbative method in order to calculate the average profile of δρ, while here we consider instead a full family of profiles for ζ within one standard deviation from the median at fixed ν. For f NL = 0, our results for C th corresponding to the median profile for ζ coincide. However, even in this case Ref. [35] finds that the average profile for the density perturbation δρ has a somewhat smaller value of C th , differing by ∼ 10%. Also, the asymptotic value of C th for large f NL which we find seems to differ from the threshold for the profiles considered in [35] by a similar amount. Although the average profile for δρ can be somewhat different from the median, it seems unlikely that the former is not within a standard deviation from the latter, corresponding to 68 % of all realisations. As shown in the left panel of Fig. 4, the dispersion of C th around the value 0.28 is quite small for large f NL . A more plausible origin for the discrepancy may be a certain inaccuracy of the lowest order perturbative approach which was used in [35]. In particular, the lowest order nonlinear correction to the mean density profile given in the left panel of Fig. 1 in [35] seems to be comparable to the zero-th order linear term, indicating the absence of any obvious small parameter in the perturbative expansion.
Recently, the effect of non-Gaussianity with the quadratic template (2.15) was also considered in reference [34], by using a slightly different method for generating the dispersion of shapes [39], which is essentially equivalent to the one we use here [18]. In [34], a fiducial value C th ≈ 0.267 was used independently of the value of f NL . Here, we find that C th is slightly decreasing with f NL , in the range 0 < f NL 2, which contributes to enhancing the abundance of PBH with growing f NL relative to the Gaussian case.
The total effect of non-Gaussianity in the perturbative template can be inferred directly from the left panel of Fig. 3, where we see that µ th changes by roughly a factor of 2 as f NL varies from 0 to 6. Since the abundance of PBHs is exponential in µ 2 /σ 2 0 , it follows that for the larger values of non-Gaussianity the power spectrum can be about a factor of 4 or so smaller than it is for the Gaussian case.

Case B: Non-perturbative template
Let us now consider the single field model where the inflaton overshoots a barrier, as in Fig.  1. In this case the, the curvature perturbation is related to the Gaussian variable through Figure 3. Results with ∆ζ = 0. The orange and blue points represents the values got using the perturbative ζ A and the non perturbative template ζ B with the corresponding error bars. The red points are those computed using the universal law of (3.24). The inner plot represents the deviation d =| µ N th − µ A th | /µ N th between the numerical µ N th and the analytical values µ A th (the same is applied for C th ). We also show in dashed line the critical amplitude ζ * ≡ µ * , such that a perturbation jumps into the false local minimum of the potential. For values of f NL ∼ 3 − 4, the thresholds for collapse approaches this limit. left) Variation of the threshold for the amplitude µ th with respect to the non-Gaussian parameter f NL . right) Variation of the threshold for the maximum of the compaction function C th with respect to the non-Gaussian parameter f NL . Figure 4. Results with ∆ζ = 0 including the dispersion term of (2.6). left) Variation of the threshold for the amplitude µ th with respect to the non-Gaussian parameter f NL , for both the perturbative template ζ A (orange) and the non perturbative template ζ B (blue). The shaded region indicates the dispersion in the numerical results from the dispersion of shapes. right) Variation of the threshold for the maximum of the compaction function C th with respect to the non-Gaussian parameter f NL . While for the perturbative template, the threshold for the compaction function is constant for large f NL , for the non perturbative template the threshold keeps evolving with increasing f NL . the non-perturbative relation (2.13).
For f NL 1 we expect the results of Case B to be very similar to Case A, and indeed this can be seen in Figs. 3 and 4. Even though the power spectra are slightly different in both cases, the Dirac delta spectrum seems a good approximation to the sharp spike (2.11) which follows from the one-field model. In Figs. 3 and 4 we also plot the curve f NL ≡ 5/(6µ * ) as a dashed line. Note that for ζ g ∼ µ * non-linearities are very important, and in fact for ζ g > µ * the backward fluctuation in the inflaton potential causes a horizon sized region to remain stuck in the false vacumm [18]. We find that the threshold µ th for adiabatic perturbations to collapse into PBHs approaches the critical value µ * for f NL ∼ 3.5. Around this value of f NL , black holes will actually be more likely to be formed though to the creation of false vacuum regions than by adiabatic perturbations. Indeed, we can calculate the abundance of black holes produced by the latter mechanism relative to the abundance of black holes with a baby universe in their interior, given by where we have used the peak theory prescription for computing number density of high peaks [47], for ν = µ th /σ 0 1. In the same limit, their ratio is then simply given by where we have also used the fact that µ th < µ * . Note that PBHs created from large overdensities follow the critical collapse scaling, and therefore their mass can range from zero up to the mass contained within the horizon at the time of their formation. On the other hand, PBHs formed from false vacuum bubbles will have a mass which is a fixed (order one) fraction of the mass of radiation contained within a horizon sized region [25]. The ratio (4.4) is then an upper bound on the dark matter fraction in the form of standard PBHs relative to that in the form of PBHs containing a baby universe. In Fig. 5 we show this ratio as a function of the non-Gaussian parameter f NL , for different values of ν. We see that for f NL < 3 standard black holes dominate, for 3 < f NL < 4 both types of black holes are produced with a comparable abundance and for f NL > 4, black holes with baby universe in their interior dominate. In principle, as mentioned above, both populations could be distinguished if we could measure the mass distribution of PBH accurately enough to tell whether it follows the critical collapse distribution or it is instead very monochromatic. Whether this can be done realistically is an interesting open question.

Summary and conclusions
In this paper we have investigated the effect of non-Gaussianities, and of the statistical dispersion in the shape of high peaks, on the threshold for PBH formation.
We assume that the fluctuations δφ of the inflaton field at the time of horizon crossing are Gaussian distributed, so that ζ g given in Eq. (1.1) is a Gaussian random field. This variable is non-linearly related to the standard curvature perturbation through a local relation ζ = ζ(ζ g ). In cosmological perturbation theory, where ζ is typically very small, it is customary to expand the local relation to second order in ζ g . The parameter f NL is then defined as the coefficient of the quadratic term. However, in the context of PBH formation, the curvature perturbation ζ is sizable, and it is important to consider the full non-perturbative relation between ζ and ζ g . In particular, this reveals a new regime for PBH formation through the collapse of false vacuum bubbles. These formed at places where a large fluctuation prevented the inflaton from overshooting a small barrier on the slope of the potential [18].
For the evolution of large adiabatic perturbations, we have used the numerical code developed in [45]. We have investigated the threshold amplitude µ th for the curvature perturbation to trigger gravitational collapse, and the corresponding threshold C th in terms of the compaction function, in two different scenarios. In Case A, we used the standard template for perturbative non-Gaussianity, parametrized by f NL , and a monochromatic power spectrum. Case B is based on a more realistic scenario where the inflaton overshoots a barrier, and we use the non-perturbative template for non-Gaussianity. The results of numerical evolution have also been compared successfully with the universal threshold (3.24), for a broad range of f NL . Both methods agree within a deviation of ∼ 2% as was reported in [46]. For the median profiles, the results of this comparison are plotted in Fig. 3.
The results which include the dispersion of shapes are summarized in Fig. 4. We find that the effect of the dispersion of shapes on C th is small, while its larger on the threshold for the amplitude of fluctuations µ th , particularly at low f NL . We find that the impact of non-Gaussianity on the thresholds is more substantial in the non-perturbative treatment. For instance, while C th saturates to a constant for large f NL in the perturbative template, we find that in the non-perturbative template it decays approximately linearly as C th ≈ 0.29 − 0.03 × f NL for f NL 3.75. Numerically, it is hard to probe larger values of f NL because µ th approaches µ * , and the profiles become extremely peaked near the origin. Nonetheless, we expect the linear behaviour to saturate to its lowest possible value C th ≈ 1/5 for f NL 4.
Finally, we have computed the relative abundance of PBHs coming from the normal collapse of an overdensity with respect to those coming from false vacuum regions. We conclude that false vacuum regions dominate the production of PBHs for f NL 3.5. PBHs created from large overdensities have a distribution of masses which follows from the critical collapse scaling and the dispersion in shapes, whereas those created from false vacuum bubbles have a fairly monochromatic spectrum. Prospects for observational discrimination of these two possibilities remain an interesting direction for further reseach. Another possi-ble phenomenological application of our results may be in the study of gravitational waves induced by non-Gaussian scalar perturbations [55].