Rare Events Are Nonperturbative: Primordial Black Holes From Heavy-Tailed Distributions

In recent years it has been noted that the perturbative treatment of the statistics of fluctuations may fail to make correct predictions for the abundance of primordial black holes (PBHs). Moreover, it has been shown in some explicit single-field examples that the nonperturbative effects may lead to an exponential tail for the probability distribution function (PDF) of fluctuations responsible for PBH formation -- in contrast to the PDF being Gaussian, as suggested by perturbation theory. In this paper, we advocate that the so-called $\delta N$ formalism can be considered as a simple, yet effective, tool for the nonperturbative estimate of the tail of the PDF. We discuss the criteria a model needs to satisfy so that the results of the classical $\delta N$ formalism can be trusted and most possible complications due to the quantum nature of fluctuations can be avoided. As a proof of concept, we then apply this method to a simple example and show that the tail of the PDF can be even {\it heavier} than exponential, leading to a significant enhancement of the PBH formation probability, compared with the predictions of the perturbation theory. Our results, along with other related findings, motivate the invention of new, nonperturbative methods for the problem and open up new ideas on generating PBHs with notable abundance.


I. INTRODUCTION
The observation of gravitational waves from binary systems by LIGO/Virgo [1], LIGO/VIRGO/KAGRA [2] and the possibility that the PBHs are the source of these events has revived interest in the study of PBHs [3][4][5][6][7][8][9][10][11][12]. The origin of the PBHs is large curvature perturbations ζ > ∼ 1 that enter the horizon during the radiation era in the early universe and collapse an entire horizon-sized region into a black hole [13][14][15][16]. Such fluctuations are known to be generated during inflation with a PDF ρ ζ , to be determined by the physics governing inflation. The abundance of produced PBHs is then characterized by where ζ c is some critical threshold that has to be calculated from collapse models [12]. Here we ignore the subtleties associated with the determination of ζ c and simply assume ζ c = 1. A major obstacle in computing β is the difficulty in calculating ρ ζ for large values of ζ, i.e., at its tail. At a naïve level, one may assume ρ ζ is Gaussian (as predicted for small perturbations and as observed on the CMB) and extrapolate the observed power spectrum of curvature perturbations on CMB scales to PBH scales. This leads to minuscule PBH abundance for most inflationary models, unless exotic features are added to the potential to enhance the power at shorter scales. However, the Gaussian approximation is clearly unjustified away from * sina.hooshangi@ipm.ir † mh.namjoo@ipm.ir ‡ mnoorbala@ut.ac.ir the peak of the PDF and a reliable calculation of ρ ζ at its tail is necessarily of nonperturbative nature. There have been a number of nonperturbative investigations taking into account various effects that give rise to an enhancement in the PBH abundance, sometimes by several orders of magnitude [17][18][19][20][21][22][23][24][25][26]. This enhancement is due to lifting the Gaussian fall-off of ρ ζ to a slower type of decay; the slowest single-field case being an exponential fall-off e −Λζ . This raises a natural question whether a heavier tail for ρ ζ , slower than any exponential, is possible. We show that a heavy tail is indeed possible. To achieve this, we employ the δN formalism [27][28][29][30][31], which is much simpler than other approaches used elsewhere. Unlike the conventional usage of the δN formalism in other applications, which is based on Taylor expanding N in small perturbations, we work with the full nonlinear N and study how this affects the PDF of ζ and modifies its tail. This nonperturbative view of the δN formalism is permitted since it holds nonlinearly for large-scale fluctuations. We shall argue that there will be conditions under which our approach is valid and other effects-not captured by the classical δN formalism-can be neglected. We will discuss these conditions in general and for a toy model example that we use to demonstrate our idea. Throughout this paper, we set M P = 1/ √ 8πG N = 1.

II. OUTLINE OF THE METHOD
According to the δN formalism, the curvature perturbation ζ on a final uniform density surface is related to the perturbation δφ on an initial surface by where N (φ) is the total number of e-folds (in an unperturbed universe) between the initial hypersurface, labeled arXiv:2112.04520v2 [astro-ph.CO] 2 Sep 2022 by the field value φ, and the final uniform density hypersurface;φ is the unperturbed initial value of the field, and the mode under consideration is super-horizon throughout this evolution. The conventional approach is to treat this equation perturbatively and make Taylor expansions to first order or second order The PDF of the field perturbations δφ on the initial surface is often assumed to be Gaussian: where σ δφ is equal to H/2π at the initial time. It then follows that, to first order, the PDF of ζ is also Gaussian with σ ζ = |N (φ)|σ δφ . Including the second order corrections, one finds a non-Gaussianity ∝ N /N 2 . We don't appeal to such Taylor expansions in this paper; instead use the fully nonlinear relation between ζ and δφ.
In particular, the PDF of ζ reads where, for simplicity, we have assumed that N (φ) is a one-to-one function. An important consequence of this nonlinearity is that large values of ζ need not be mapped to large values of δφ. Therefore, the tail of ρ ζ does not need to inherit the Gaussian nature of the tail of ρ δφ . This is a key point that plays a central role in our approach. Given a potential V (φ), our strategy to calculate the PDF of ζ at a givenφ, corresponding to a given PBH scale, is as follows: We assume that ρ δφ is Gaussian on the initial surface atφ. Next we calculate N without resort to any slow-roll approximation. 1 To do so, we use the equation of motion where prime denotes derivative with respect to the number of e-folds defined by dN = Hdt. 2 We define the final hypersurface (after which ζ remains intact until the horizon reentry during the radiation era) by φ = φ e . We can 1 Since the nontrivial evolution of ζ in single-field models takes place in a non-attractor phase where the slow-roll conditions are violated, it is essential to avoid slow-roll approximations. 2 Throughout this paper, a prime on a function denotes differentiation with respect to its natural argument. So N means dN /dφ whereas φ means dφ/dN . then find N (φ + δφ) by solving Eq. (7) with initial conditionsφ + δφ andπ =φ for a range of different values of δφ and then requiring that φ(N = N (φ+δφ)) = φ e . This yields N as a function of the initial valuesφ + δφ, which can then be inserted in Eq. (6), in conjunction with (5), to provide the PDF of ζ. Notice that the initial field velocity will not be perturbed (while, in principle, N does depend on both initial conditions). This will be justified later after implementing the method in a simple model. In the next section we work out the details for a particular potential leading to a heavy-tailed PDF. Here, we clarify what we mean by a heavy tail (which is consistent with the standard nomenclature in mathematics) [32]: dζ beats any exponentially decaying function, i.e., if for every λ > 0, we have lim sup ζ→+∞ e λζF (ζ) = +∞. This is not very easy to check in a numerical calculation, so we employ a related definition that is more practical: Note that both definitions exclude exponential tails (where D approaches a positive constant) and of course Gaussian tails, but admit power-law tails, as well as those falling like exp[−kζ p ] with 0 < p < 1 and k > 0.

III. A TOY MODEL
As a proof of concept we carry out the above procedure for a simple toy model. Let us suggest the potential and choose the unperturbed initial value to beφ = 0 (equivalently, one may replace φ by φ −φ throughout). We could work with a number of other inflationary potentials, but we can offer a heuristic derivation for this particular V based on requiring a power-law tail for the PDF of ζ, that turns out to be of the form ρ ζ ∼ ζ −2 (hence D = 2/ζ), starting from a Gaussian PDF for δφ. See Appendix A for the details of this derivation. We now follow the method outlined in the previous section to obtain the PDF associated with the potential (8) by choosing numerical values forφ > φ e and for the range δφ min < δφ < δφ max . Importantly, notice that we assume, for simplicity, an abrupt transition to a slowroll phase of inflation immediately after φ e . As argued in [33] and further investigated in [34] this sharp transition leads to the conservation of ζ just after the transition so that following its evolution up to the end of inflation will not be needed. It would be interesting to study the effect of a milder transition to our results but it is beyond the scope of this paper. We require Eq. (8) to be valid only in the range φ e < φ <φ + δφ max ; beyond that, it suffices that a slow-roll phase be realized. Furthermore, notice that we choose an initial velocity that leads to a non-attractor phase of inflation from φ to near φ e . This phase leads to the nontrivial evolution of ζ on super-horizon scales, leading to a non-standard PDF for ζ.
In Fig. 1 we show ρ ζ as a result of the procedure outlined above. In addition to plotting the fully nonlinear result ρ ζ , we have computed N (φ) and N (φ) numerically, inserted them in Eqs. (3) and (4), and followed the rest of the numerical calculation again using Eq. (6) but now with these linear and quadratic relations between ζ and δφ. This gives ρ (1) ζ and ρ (2) ζ , respectively, which depart from ρ ζ around ζ ∼ 1 and are markedly suppressed afterwards. The values of the parameters used in this calculation are quoted in the caption of Fig. 1. Fig. 2 is intended to show the heaviness of the tail by plotting D, defined in Def. 2. Again, we have demonstrated the result for D (1) and D (2) , too. The nonstop fall-off trend in D is a plausible witness for the heaviness of its tail. We are further assured by the perfect agreement between the predicted analytic behavior D = 2/ζ The parameter β in Eq. (1), which is related to the PBH abundance, is plotted in Fig. 3 as a function of σ δφ = H/2π. We have also plotted β (1) and β (2) , which are derived from ρ (1) and ρ (2) , respectively. We see that taking into account the full nonlinear nature of the δN formula (2) amounts to considerable enhancement in the PBH abundance compared to the perturbative calculation β (1) based on the Gaussian PDF ρ (1) ζ . One may argue that maybe it was unnecessary to incorporate the full nonlinear form of Eq. (2) and that if one included the second order term as in Eq. (4), the same level of enhancement could be reached by β (2) . However this is not the case in the present model, as Fig. 3 clearly shows that β (2) β by several orders of magnitude. We conclude that the result, including the enhancement, is intrinsically nonperturbative. We further note that the large ζ behavior of ρ ζ while consistent with analytical predictions, does not suffice to estimate β since the integral in Eq. (1) is mainly supported around ζ ∼ O(1). As a result-just like the ζ 1 limit-the methods that only predict the behavior of the PDF in the ζ 1 limit may not reliably estimate β. It is indeed the transient regime between the two limits that contributes to β, which might be considered as a serious technical difficulty. This is not an issue here as we were able to construct the full PDF using the δN formalism. However, note that our approach relies on some hidden assumptions that need to be satisfied. We will briefly discuss these assumptions in the next section and confirm that they are indeed met in our presented example.

IV. UNDERLYING ASSUMPTIONS
In this section we justify the various underlying assumptions and approximations made in our work.
Non-Gaussianities in δφ: We have assumed that δφ has a Gaussian PDF. Although it is straightforward to incorporate a non-Gaussian PDF into our numerical calculations, the precise form of such a PDF requires an in-in calculation and depends on the exact form of the interactions dictated by higher order terms L n (n > 2) in the Lagrangian of δφ. If these higher order (interaction) terms are negligible compared to the second order (free) term L 2 , i.e., if L n L 2 for n > 2, we can safely ignore the non-Gaussianity of ρ δφ for the range of δφ under consideration. To this end, we estimateδ φ and 1 a ∇δφ by Hδφ, so that a typical term in the Lagrangian is of order L n ∼ c n (t)δφ n max where δφ max corresponds to ζ max , the largest fluctuation that is relevant to our calculation of the PBH abundance. The ratio of interest is therefore L n /L 2 ∼ (c n /c 2 )δφ n−2 max which we need to ensure remains small. In practice, it often suffices to check the size of L 3 /L 2 and in most situations it is the matter sectorrather than the gravity sector-that contributes most to the interactions, reducing this ratio to |V /V |δφ max .
Horizon Crossing Effects: There is a common lore that the horizon crossing effects cause inaccuracies in the δN formalism of the order of the slow-roll parameters. Accordingly, a naïve application of the δN formalism does not yield the non-Gaussianity consistency relation f NL ∝ (n s − 1) of single-field inflation [35]. It may therefore be a source of concern that such effects may as well spoil our nonperturbative calculation. However, it is not true that these horizon crossing effects render the δN formalism invalid. They do cause corrections, but those corrections appear in the correlations of δφ at the initial surface, whose deviation from Gaussianity is controlled by interaction terms L n as we mentioned in the previous paragraph. 3 Therefore, the smallness of L n /L 2 can simultaneously resolve such concerns.
Stochastic Effects: We have neglected stochastic effects throughout this paper. Indeed there are two sources of randomness that give rise to a PDF for ζ. First, there is the randomness of the initial conditions δφ of the field at the initial surface, which is incorporated as ρ δφ in our calculation. The second source is the stochastic noise that the field encounters along the way. Let us consider how a rare event ζ > ∼ 1 can be generated under the effect of these sources of randomness. One possibility is to start with a rare initial condition (from the tail of δφ) but then proceed on a typical trajectory (close to the deterministic evolution of the field) toward the end of inflation. Another possibility is to start with a typical initial field value (around the peak of ρ δφ ) but then ζmax = 1 ζmax = 10 5 L3/L2 0.14 0.71 S 0.0060 0.75 follow a rare path (highly fluctuating and far from the average trajectory) until the end of inflation. It is highly unlikely that both extreme events occur simultaneously in a single realization, so it is fair to say that the probability of a rare event ζ > ∼ 1 is the sum of the probabilities of these two alternative possibilities.
We have only considered the first possibility in this work, so strictly speaking, we have a lower bound on the PDF of ζ and hence on the PBH abundance. This is already a considerably enhanced probability compared to the Gaussian or exponential suppression. To justify the approximation that the inflaton proceeds on a typical path, let us estimate the size of the stochastic noise when diffusion is small. It is well known that the stochastic noise corresponds to random jumps of magnitude H/2π per e-fold that are superimposed on top of the classical equation of motion of φ. Let ∆φ be the total field excursion in a given path from the initial valueφ + δφ to φ e . The expected displacement due to the accumulated jumps is then, just as in the conventional random walk, equal to (H/2π) √ N , where N is the total number of efolds during this excursion. Therefore, a fair measure of relative importance of the stochastic effects is given by so we need to check that S is small. We have not considered the second possibility (following a rare path) here. But other studies [17,19,21,37] claim no more than an exponential enhancement due to the stochastic effects. It is therefore reasonable to assume that the main contribution to the PBH abundance comes from our approach.
Finally, note thatπ must be chosen such that the inflaton classically reaches φ e for all δφ ∈ [δφ min , δφ max ].
Contribution of δπ: We have ignored the fluctuations to the conjugate momentum π = φ , so we work with a fixed value of φ 0 set equal toφ . This is justified whenever the mode function δφ freezes and hence δπ decays exponentially. Since δφ is the perturbation of an almost massless field during inflation, this assumption is well justified.
Before closing this section, let us verify the validity of these assumptions in the toy model of the previous section. Let us first consider the set of parameters used in Figs. 1 and 2. In the single-field model at hand, there are 11 terms L 3 in the cubic Lagrangian in the flat gauge each with its own c 3 [35]. We quote the maximum of the ratios L 3 /L 2 for two values of ζ max in Table I. For each individual path with initial condition δφ, we calculate ∆φ and N . The maximum S among all paths with ζ < ζ max in our numerical calculation is given in Table I. It is remarkable that both L 3 /L 2 and S are small, even for ζ as large as 10 5 . Thus ρ ζ is reliable over the entire visible range in Figs. 1 and 2-including the heavy tail. Fig. 3 uses a different set of parameters for which our ρ ζ at the tail is not as reliable. However, those values are not relevant for the purpose of calculating β; instead we can confirm that L 3 /L 2 ≈ 0.22, and S varies between 0.02 to 0.31 for ζ ∼ 1 which is the region relevant to β. So we see that in all cases the approximations are justified.

V. CONCLUDING REMARKS
Our primary point has been to give a proof of concept. There are several issues that need to be addressed to calculate the PBH abundance reliably. To begin, the Eq. (1) for β is only a crude estimate. A more accurate treatment requires more rigorous determination of the threshold [38][39][40]. Actually, the physically relevant quantity is the density contrast δ rather than ζ. These subtleties modify our estimate of β, but not the fact that it differs significantly from the predictions of perturbation theory. The second issue is the contributions of the stochastic effects, which as we discussed, can further lift the PDF, albeit not as heavily as we have found. As we mentioned before, we have not considered these effects by working in the drift-dominated regime, but they deserve further study. In particular, it would be interesting to explore scenarios where both classical and stochastic effects (corresponding, respectively, to the initial and onthe-way fluctuations) are equally important.
There are various extensions that can be made to this work. In a future publication, we will apply the nonperturbative δN to other single-field models to see how generic the heavy tail might be and also to investigate how β may be modified compared to the perturbative treatments. As another extension, one can apply the same methods to multi-field models of inflation [41]. Furthermore, the possibility of obtaining heavy tails in multi-field scenarios seems to be fairly unexplored (see, however, Ref. [18]). It is also interesting to investigate if the contribution of the tail region to β can be made more pronounced.

ACKNOWLEDGMENTS
We thank Eva Silverstein for useful comments. M.N. acknowledges financial support from the research council of University of Tehran.
Appendix A: Speculating a potential leading to a heavy tail Here we briefly outline the analysis that suggests the potential (8), leading to a power-law tail for ρ ζ . For large values of ζ, we like to avoid the Gaussian fall-off due to ρ δφ in Eq. (6). To achieve this, it would be plausible to have a situation where the ρ δφ factor approaches a constant in the same limit. A simple relation like the following fulfills our requirement: where α, ζ 0 and ζ c are constants. From the last relation it is evident that δφ and hence ρ δφ approaches a constant, as desired. In fact, we have achieved more: The remaining factor, i.e., dδφ/dζ gives a power-law decay of the PDF scaling as 1/ζ 2 for large ζ, which is indeed heavy-tailed with D = 2/ζ → 0. We need to show that there is indeed a potential that supports this solution.
Let us denote the number of e-folds required to reach φ 2 from the initial configuration φ 1 by N (φ 1 , φ 2 ). 4 Inverting Eq. (A1) and invoking Eq. (2), this implies From this, we may identify Therefore, which suggests that Notice that these expressions were-and still areapproximate (for large ζ) and may not be taken literally. In particular, Eq. (A5) is not well defined when φ 1 = φ 2 . Nonetheless, we can use them to guess a potential that may work. Assuming that the inflaton field evolves monotonically, we can write Eq. (7) in the following equivalent form Using Eq. (A5) for φ 1 =φ and φ 2 = φ(N ), a simple integration results in We have made two approximations in deriving this potential. First, we have ignored the first slow-roll parameter in Eq. (A6) ( 1 2 N −2 ,φ = 1 2 φ 2 = 1 which is valid in an inflationary background) without assuming anything about higher slow-roll parameters. Secondly, we have assumed α|φ −φ| 1 and α|φ −φ| 3 1 to obtain a simpler form for the potential, although this was unnecessary. Notice that the potential depends onφ (and, in principle, onπ) which are the initial conditions at the initial hypersurface. This is not problematic; it just means that for a given set of initial conditions a potential of the form (A7) can lead to a heavy tail. This concludes our motivation for considering the potential (A7). We emphasize again that this is based on some guesswork and a dubious reader may choose to accept it only as a suggestion and wait to see if it indeed leads to a heavy tail.