Non-gaussianities for primordial black hole formation
NNon-gaussianities for primordial black hole formation
Marco Taoso a and Alfredo Urbano b,c a I.N.F.N. sezione di Torino, via P. Giuria 1, I-10125 Torino, Italy b Dipartimento di Fisica, “Sapienza” Università di Roma, Piazzale Aldo Moro 5, 00185, Roma, Italy and c I.F.P.U., Institute for Fundamental Physics of the Universe, via Beirut 2, I-34014 Trieste, Italy. (Dated: February 9, 2021)We analyze primordial non-gaussianities in presence of an ultra-slow phase during the inflationary dynam-ics, focusing on scenarios relevant for the production of primordial black holes. We compute the three-point correlation function of comoving curvature perturbations finding that non-gaussianities are sizable,and predominantly local. In the context of threshold statistics, we analyze their impact for the abundanceof primordial black holes, and their interplay with the non-gaussianities arising from the non-linear rela-tion between density and curvature perturbations. We find that non-gaussianities significantly modify theestimate of the primordial black holes abundance obtained with the gaussian approximation. However,we show that this effect can be compensated by a small change, of a factor ÷ at most, of the amplitudeof the primordial power spectrum of curvature perturbations. This is obtained with a small tuning of theparameters of the inflationary model. I. MOTIVATIONS
The theory of inflation provides an elegant mechanism that explains the origin of structures in the Universe. In theinflationary picture, space-time fluctuates quantum mechanically around a background that is expanding exponen-tially fast. After the end of inflation, these fluctuations are transferred to the radiation field, creating slightly overdenseand under-dense regions. These tiny ripples are the initial conditions that start the process of gravitational collapsewhich eventually forms the intricate architecture of galaxies that populate our observable Universe today.If the amplitude of some of the overdense regions exceeds a specific threshold value, the pull of gravity becomesunsustainable and the gravitational collapse directly forms primordial black holes (PBHs) [1]. Interestingly, for massesin the range (cid:46) M PBH [g] (cid:46) , a population of such PBHs may account for the totality of dark matter observedin the Universe today, as envisaged by Hawking and Carr in their pioneering research [2, 3].Because of their intrinsic quantum-mechanical origin, the way in which quantum fluctuations lead to a classicalpattern of perturbations can be described only in a probabilistic sense, and it is often assumed that primordial per-turbations follow a gaussian statistics. It is important to stress that there is no fundamental reason to believe that thisassumption is true, and it should be rather considered to be an approximation. This approximation seems to workexquisitely well at the typical scales which are relevant for cosmic microwave background (CMB) anisotropy measure-ments, . (cid:46) k [Mpc − ] (cid:46) . , and this is because deviations from the gaussian pattern are expected to be sup-pressed by powers of the slow-roll parameters [4] which during conventional slow-roll dynamics take O ( (cid:28) values.As far as the formation of PBHs is concerned, however, the story might be drastically different. The process of PBHformation involves completely different scales compared to those probed by CMB observations. Typically, in orderto be compatible with observational bounds on PBH abundance in the present-day Universe, one needs to considerscales (cid:46) k [Mpc − ] (cid:46) . At such small scales, the physics involved can be completely different compared towhat expected at CMB scales. Indeed, the formation of fluctuations large enough to trigger gravitational collapse intoPBHs requires a departure from conventional slow-roll dynamics, and the slow-roll parameters can easily reach O (1) values. Consequently, the parametric suppression that limits the presence of non-gaussianities at CMB scales is notguaranteed anymore.For illustration, we summarize the time-evolution of physical scales probed by CMB observations and scales relevantfor PBH formation in fig. 1.The introductory discussion above makes clear that the possible presence of non-gaussianities plays a fundamentalrole in the precise computation of the PBH abundance. The formation of PBHs requires perturbations above somethreshold, located in the tail of their probability distribution. Small deviations from the gaussian shape of the tail,therefore, may produce exponentially different results.In this work we investigate the role of non-gaussianities in the context of inflationary models featuring an ultraslow roll (USR) phase [5, 6]. This part of the inflationary dynamics has been invoked for the production of PBHs (seerefs. [7, 8] for the earliest proposal in this direction). This is because, in these scenarios, an USR phase allows to boostthe primordial fluctuations at small scales, leading to a sizable population of PBHs. Crucially, conventional slow-rolldynamics is violated during such a period, potentially enabling the generation of large non-gaussianities. In our anal-ysis, we focus on single-field inflation models in which the USR phase arises when the inflaton, rolling down its poten-tial, crosses an approximate stationary inflection point before the end of inflation. Specifically, we consider both a toymodel, which allows an analytical exposition, and a more realistic parametrization of the inflaton potential.There exist a number of preceding studies devoted to the understanding of non-gaussianities in the context of PBHformation, and we will properly acknowledge them in the course of our work. Before diving into the details of our a r X i v : . [ a s t r o - ph . C O ] F e b analysis, let us sketch briefly the structure of this paper.In section II, we will frame in a more formal way the computation of the abundance of PBHs and the connectionwith inflationary dynamics. In section III, we will compute the so-called bispectrum of comoving curvature perturba-tions which controls the first non-trivial (i.e. connected) deviation from exact gaussianity. The bispectrum vanishesidentically for a gaussian distribution, and a large non-zero result will introduce sizable deviations from the gaussianstatistics. In section IV, we will move to consider the connected part of the trispectrum of comoving curvature pertur-bations which controls deviation from exact gaussianity at the level of fourth-order statistical correlators. We will notcompute the full connected trispectrum but we will argue that the non-zero contribution which is generated by cubicinteractions plays a crucial role in the computation of the PBH abundance. In section V we discuss the computation ofPBH abundance in the presence of non-gaussianities. We will show that their impact is not dramatic. We find that thesize of non-gaussianities is large enough to significantly affect the PBH abundance. On the other hand, it is enough tochange the amplitude of the power spectrum of comoving curvature perturbations of a factor ÷ to obtain the sameabundance obtained with the gaussian approximation.We present our conclusions in section VI. Appendix A complements the computations illustrated in the main textwith further technical details. � * (cid:1) ℯ (cid:1) (cid:1) � (cid:1) (cid:2) � �� � ��� �� - �� �� - �� �� - �� �� - �� �� - �� �� - �� � ��� ( (cid:1) ) �� � � � � � �� � � � � � � � � � � � [ � � - � ] ≈ �� ℯ � (cid:4)ℴℓ(cid:7)(cid:8) ℴ(cid:4) (cid:9)(cid:10)(cid:4)ℓ(cid:1)(cid:11)(cid:9)ℴ(cid:10) (cid:12)(cid:1)(cid:7)(cid:9)(cid:1)(cid:11)(cid:9)ℴ(cid:10) ℯ(cid:13)ℴ(cid:14)(cid:15) Λ ℴ (cid:2) (cid:3) ℯ (cid:5) (cid:6) (cid:7) (cid:8) (cid:9) ℴ (cid:10) (cid:7) ℓ (cid:12) (cid:9) (cid:10) (cid:13) ℴ (cid:12) (cid:1) ℬ ℋ (cid:4) (cid:5) (cid:6) ℓ ℯ (cid:4) (cid:1) ℯ (cid:3) ℴ (cid:5) (cid:6) (cid:7) (cid:8) (cid:9) (cid:10) (cid:7) ℴ (cid:8) � � � (cid:11) ℓ (cid:10) (cid:1) (cid:9) (cid:13) ℓ ℴ (cid:14) � (cid:1) ℴ ℓℓ (cid:15) (cid:1) ℯ (cid:13) ℯ (cid:8) (cid:10) � (cid:16) (cid:9) (cid:17) (cid:18) (cid:8) (cid:7) (cid:19) ℯ (cid:1) (cid:13) ℯ (cid:1) ℯ (cid:20) ℯ (cid:9) (cid:10) (cid:7) (cid:8) ℊ (cid:20) ℴ (cid:1) (cid:7) (cid:22) ℴ (cid:8) (cid:3) (cid:1) ℴ (cid:13)(cid:13) (cid:7) (cid:8) ℊ (cid:23) ℴ (cid:1) (cid:24) * (cid:16) � ℯ � (cid:1) (cid:2) (cid:3) (cid:4) (cid:5) (cid:6) (cid:7) ℓ ℋ (cid:10) (cid:11)(cid:11) ℓ ℯ ℓ ℯ (cid:13) ℊ (cid:15) (cid:2) (cid:3)ℬℋ(cid:6) (cid:7)ℯ(cid:9)(cid:6)(cid:10)(cid:11)(cid:2)(cid:12) (cid:13)(cid:6) (cid:14)(cid:13)(cid:12)(cid:12)ℯ(cid:7) (cid:1) ℬ ℋ (cid:4) (cid:5) ℴ (cid:7) (cid:8) NG suppressedNG enhanced
FIG. 1:
Time-evolution of the physical Hubble length /H (red line), written in units of its present-day value /H (cid:39) . × Mpc (with H (cid:39) − GeV in natural units), as function of the logarithm of the scale factor a of the Friedmann-Lemaître-Robertson-Walker line element. We follow the evolution of /H throughout the history of our observable Universe. The inflationary phase iscomputed according to the model in ref. [27]. After (instantaneous) reheating, the evolution of /H during the radiation epoch, mat-ter epoch (m.e.) and the present-day epoch ( Λ ) is computed according to the standard Λ CDM model. We superimpose the time-evolution of physical length scales λ = a/k , with k the comoving wavenumber, that are relevant for CMB observations (the greenband labelled “observational window” with . (cid:46) k [Mpc − ] (cid:46) . and k ∗ = 0 . Mpc − ). We also show the time-evolution ofphysical length scales associated to perturbations that are responsible for PBH formation (the blue band labelled “PBH scales” with (cid:46) k [Mpc − ] (cid:46) ). Scales that cross the ultra slow-roll phase when they are deep in the super-Hubble regime are not sizablyaffected by the presence of negative friction [27]. For these scales, non-gaussianities are slow-roll suppressed. This is the case of scalesrelevant for CMB observations. On the contrary, scales that cross the Hubble radius close to the ultra slow-roll regime can be exponen-tially enhanced [27]. The interval of e -folds during which the ultra slow-roll phase takes place is indicated with a vertical strip. At thesescales, non-gaussianities are enhanced. PBH formation takes place at time t f after the PBH scales re-enter the Hubble radius duringthe radiation epoch. Once formed, we consider the simple scenario in which the PBH energy density redshifts as matter. II. ON THE ABUNDANCE OF PRIMORDIAL BLACK HOLES: FROM THEORY TO OBSERVATIONS AND BACK
In cosmology and in the theory of structure formation, we assume that fluctuations in the energy density or in thegravitational potential form a three-dimensional random field. A three-dimensional random field δ ( (cid:126)x ) is a set of ran-dom variables, one for each point in the three-dimensional real space, defined by a probability functional P [ δ ( (cid:126)x )] whichspecifies the probability for the occurrence of a particular realization of the field over the ensemble. In cosmology, therandom field δ ( (cid:126)x ) is the three-dimensional mass overdensity field which gives the fractional overdensity in a given re-gion of space (and at a given cosmic time) with respect to the unperturbed background (a.k.a. density contrast). Arandom field is fully specified by the entire hierarchy of its correlation functions, (cid:104) δ ( (cid:126)x ) δ ( (cid:126)x ) . . . δ ( (cid:126)x n ) (cid:105) . The simplesttype of random field is the gaussian random field for which the probability functional P is a gaussian distribution. Thestatistical properties of a gaussian (subscript G hereafter) random field are completely specified by the two-point cor-relation function or, working in Fourier space, by the power spectrum. The n -point correlation functions either vanish(for odd n ) or can be expressed in terms of the power spectrum as a consequence of the Isserlis’ theorem (for even n ).In Fourier space, we have (assuming a vanishing mean value as well as homogeneity and isotropy of space) (cid:104) δ k (cid:105) = 0 , (1) (cid:104) δ k δ k (cid:48) (cid:105) = (2 π ) δ (3) ( (cid:126)k + (cid:126)k (cid:48) )∆ δ ( k ) , (2) (cid:104) δ k δ k δ k (cid:105) G = 0 , (3) (cid:104) δ k δ k δ k δ k (cid:105) G = (cid:104) δ k δ k (cid:105)(cid:104) δ k δ k (cid:105) + (cid:104) δ k δ k (cid:105)(cid:104) δ k δ k (cid:105) + (cid:104) δ k δ k (cid:105)(cid:104) δ k δ k (cid:105) = (2 π ) δ (3) ( (cid:126)k + (cid:126)k ) δ (3) ( (cid:126)k + (cid:126)k )∆ δ ( k )∆ δ ( k ) + two permutations , (4)and so on. The variance of the gaussian distribution is σ = lim (cid:126)y → (cid:126)x (cid:104) δ ( (cid:126)x ) δ ( (cid:126)y ) (cid:105) = (cid:104) δ ( (cid:126)x ) δ ( (cid:126)x ) (cid:105) = (cid:90) dkk P δ ( k ) , ∆ δ ( k ) = 2 π k P δ ( k ) , (5)from which we see that the power spectrum of density fluctuations P δ ( k ) is the contribution to the variance of the fieldper unit logarithmic interval of k , and saying that P δ ( k ) ∼ means that the Fourier modes in a unit logarithmic binaround wavenumber k generate fluctuations δρ/ρ of order unity. Notice that σ does not depend on the position, byvirtue of the delta function eq. (2) which is indeed a consequence of homogeneity. We are interested in the situation inwhich there is a sizable probability to form in a given region of typical size R density fluctuations δρ/ρ large enough tocollapse and form black holes under the pull of their own gravitational weight. In order to study PBH formation andmake contact with the power spectrum of comoving curvature perturbations, eq. (5) has to be modified in two ways.First, we convolve the random field δ ( (cid:126)x ) with some window function W ( | (cid:126)x − (cid:126)y | , R ) to obtain δ R ( (cid:126)x ) = (cid:82) d (cid:126)y W ( | (cid:126)x − (cid:126)y | , R ) δ ( (cid:126)y ) so that, at every point (cid:126)x , the smoothed field represents the weighted average of δ over a spherical region ofcharacteristic dimension R centred in (cid:126)x . Second, we relate the density contrast to the comoving curvature perturbation R by means of the linear-order relation (in Fourier space) [9] δ k = 2(1 + ω )(5 + 3 ω ) (cid:18) kaH (cid:19) R k , (6)where ω = p/ρ defines the equation of state, and we shall take ω = 1 / for radiation domination. It is important toremark that eq. (6) is not an exact relation but it is only valid at the linear order, and non-linear corrections to eq. (6)unavoidably source non-gaussianities, see [10–18] for related studies. We will come back to this point in section V. Forthe moment, we shall work assuming the linear relation in eq. (6).Eq. (5) becomes (cid:104) δ (cid:105) R = σ R ≡ ω ) (5 + 3 ω ) (cid:90) dkk W ( kR ) (cid:18) kaH (cid:19) P R ( k ) , (7)where P R ( k ) is the power spectrum of the the comoving curvature perturbation. W ( kR ) is the Fourier trans-form of the window function introduced before, and we shall use the volume-normalized gaussian window func-tion W ( kR ) = exp( − k R / . Assuming that PHBs are formed when the density perturbation exceeds the thresh-old value δ th , the fraction of the Universe ending up in PHBs is given by the tail of the density contrast distribution β ( M PBH ) = γ (cid:82) ∞ δ th P ( δ ) dδ , and in the gaussian approximation we have β G ( M PBH ) = γ (cid:90) ∞ δ th dδ √ πσ R ( M PBH ) exp (cid:20) − δ σ R ( M PBH ) (cid:21) , (8)where we convert the smoothed variance σ R to a function of the PBH mass taking into account that the size R is relatedto the mass M PBH as R ≈ G N M PBH /γa form . The parameter γ depends on the details of the gravitational collapseand a form is the scale factor at the time at which the PBH is formed.The distribution of PHBs as function of their mass is strongly peaked at the value that maximizes σ ( M PBH ) , in viewof the exponential factor in eq. (8). In terms of the power spectrum P R this means that the PBH abundance roughlyscales as e − / P R . Finally, the present-day fractional abundance of dark matter in the form of PBH is given by f PBH = (cid:16) γ . (cid:17) / (cid:20) β ( M PBH )1 . × − (cid:21) (cid:20) g ∗ ( T form )106 . (cid:21) − / (cid:18) M PBH g (cid:19) − / . (9)The computation of β by means of eq. (8) is known as threshold statistics. The value of δ th is subject to some uncer-tainty (varying in the range . (cid:46) δ th (cid:46) / ), and also depends on the specific shape of the power spectrum. In ouranalysis, for simplicity, we adopt the fixed values δ th = 0 . and γ = 0 . [19].An alternative way of computing the PBH abundance relies on the use of peak theory [20], in which one computesthe number density of peaks of the overdensity field δ (which in general contains many peaks, i.e. local maxima, andvalleys, i.e. local minima) and integrate above some threshold δ c [21]. The latter in general differs from δ th and dependson the shape of the power spectrum [22] (see also ref. [10] in which the full non-linear relation between density andcurvature perturbation is explored). Using again the gaussian approximation, one finds [11, 20, 21] β pk G ( M PBH ) (cid:39) γ π (cid:18) (cid:104) k (cid:105) (cid:19) / r m (cid:18) δ c σ (cid:19) e − δ c / σ , with (cid:104) k (cid:105) = 1 σ (cid:90) d (cid:126)k (2 π ) k ∆ δ ( k ) , (10)where δ c = δ th / ψ ( r m ) with ψ the average density profile and r m is the comoving distance where the so-called com-paction function is maximized (we refer to ref. [22, 23] for a comprehensive discussion and more detailed definitions). Notice also that, following refs. [11, 22, 23], we do not use the smoothing in eq. (10) when dealing with peak theory.The reason is that when applying peak theory, and in particular when computing the threshold for collapse, the scale r m in position space associated to the peak of the power spectrum arises naturally.On general ground, the abundance obtained by means of peak theory is larger than the one provided by thresholdstatistics [20, 22]. In this paper, however, our interest is not a comparison between the two approaches. Instead, theissue that we want to explore is whether the gaussian approximation in eqs. (8, 10) is justified. It is clear that the for-mation of PBH is a tail effect and, as such, strongly depends on possible deviations from the gaussian approximation.Investigating the possible presence of non-gaussianities in the context of realistic inflationary models in which one gen-erates a sizable abundance of PBH is precisely the main goal of this work. After clarifying whether non-gaussianitiesare present and what are their main properties, we will quantify their impact on the PBH abundance focussing onthreshold statistics. We will consider the case of peak theory and present a comparison with threshold statistics in aforthcoming work [26].Let us, therefore, illustrate in more details what are the conditions that the inflationary dynamics should possess inorder to generate the large density fluctuations needed for PBH formation.In the gaussian approximation, a sizable PBH abundance is obtained for P R (cid:39) − – − at the PBH scales, re-quiring therefore a huge enhancement of P R with respect to its value at the CMB scales, P R (cid:39) − . This can happenin presence of an USR phase during inflation. In fact, during such a transitory part of the inflationary dynamics, theamplitude of the comoving curvature perturbations can be exponentially enhanced. Let us introduce the Hubble pa-rameters (cid:15) and η : (cid:15) ≡ − ˙ HH η ≡ − ¨ H H ˙ H , (12)where H is the Hubble function and dot is the derivative with respect to time. During standard slow roll inflation | (cid:15) | (cid:28) , | η | (cid:28) and, in the super-horizon limit, the evolution of R k is a linear combination of a constant mode and anexponentially decaying one. The latter becomes quickly negligible and the amplitude of R k freezes to a constant valueuntil the end of inflation. Instead, if (cid:15) − η (cid:39) − η < , a condition which we take as a definition of the USRphase, the decaying mode becomes a growing one. In practice R k evolves as dumped harmonic oscillator, with thesign of the friction term controlled by the combination (cid:15) − η. During the USR phase this sign flips and the frictionterm becomes a driving force. This is responsible for the exponential growth of R k and the enhancement of P R at smallscales. Notice however that for certain comoving wavelengths the amplitude of the growing mode (once it freezes to aconstant value after the USR phase) is opposite to the one of the constant mode. This cancellation produces a dip inthe power spectrum. The outcome of this intricate dynamics can produce a power spectrum of comoving curvatureperturbations as the one in fig. 3. Eq. (10) is obtained as follows. Introducing ν ≡ δ/σ, the comoving number density of peaks is [21] n pk G , com = 1(2 π ) (cid:18) (cid:104) k (cid:105) (cid:19) / ν e − ν , (11)and in physical space one has n pk G , phys = n pk G , com /a f , with a f is the scale factor at the formation time t f . The PBH mass reads M PBH ( ν ) = K M PBH ( t m ) σ ˜ γ ( ν − ν c ) ˜ γ where ˜ γ is a critical exponent that depends on the equation of state in the formation era (with ˜ γ (cid:39) . for radia-tion [24]) and K ∼ O (1) . M PBH ( t m ) is the horizon mass at horizon-crossing time t m defined implicitly by the relation a ( t m ) H ( t m ) r m = 1 . One can write β pk G = (cid:82) ∞ ν c dν M PBH ( ν ) n pk G , phys ( ν ) /ρ f where ρ f is the background density at the formation time. Using the result from numer-ical simulation a f /a m ∼ [23, 25] and in the limit ν c ≡ δ c /σ (cid:29) one obtains eq. (10) with γ = 3 K Γ(1 + ˜ γ ) ν − ˜ γc σ ˜ γ . To make a quantitative analysis we shall consider an explicit inflationary model where the USR phase is realized. Wefocus on the single field model discussed in [27], described by the following action in the Jordan frame S = (cid:90) d x √− g (cid:20) − (cid:0) M P + ξ φ (cid:1) R + 12 g µν ∂ µ φ ∂ ν φ − V ( φ ) (cid:21) . (13)where φ is the inflaton field, M P is the reduced Planck mass and ξ the non-minimal coupling to gravity. The potentialis a polyonomial V ( φ ) = a φ + a φ + a φ + N (cid:88) n =5 a n φ n . (14)and in full generality contains also higher order operators ( n > ). These are not mandatory for the production ofPBH but, as noted in [27], they improve the compatibility of the model with the CMB observables. After a judicioustuning of the coefficients of the quadratic and cubic terms, the potential presents an approximate stationary inflectionpoint. This key features induces an USR phase in the inflationary dynamics. The dashed blue line in fig. 2 shows theevolution of η as a function of the e -fold time for a specific choice of the parameters of the model. As it is evident, thecondition η > / , and therefore the USR phase, is realized for a while. This solution, depicted with a cyan star in [27],is in agreement with CMB observables and present 100% of DM in the form of PBHs (the abundance is computed witheq. (8)). For concreteness, in the rest of the paper, we focus on this benchmark point for our numerical analysis. Werefer to it as numerical model in the following. However, let us stress that our results will not be limited to this specificsolution. On the contrary, as we will see, our conclusion will be of general validity for all the single-field inflationarymodels in which PBH production is induced by the presence of an approximate stationary inflection point before theend of inflation.As shown by the black line in fig. 2, the evolution of η can be approximated, at least at the qualitative level, by a piece-wise function. This observation has prompted a simple analytical model to describe the inflationary dynamics [27, 28].The inflationary phase is splitted in the three periods, each of them with constant values of η (and (cid:15) (cid:28) ). In this wayone can obtain an analytical solution for the evolution of R k and consequently its power spectrum. The advantage isthat one can gain more insight into the features of the USR dynamics, as the enhancement or decrease of R k for certaincomoving wavelenghts (see [27] for an in-depth discussion). Moreover, the analytical model does not rely on a specificinflationary model, rather, changing its parameters, it can describe different USR phases. The drawback is that it is notrealistic enough for a detailed comparison with observables. In the rest of the paper we will compare the results basedon the analytical model with the outcome of a numerical analysis performed in the context of the inflationary modeldescribed above.Let us briefly summarize the main features of the analytical model. The three regions discussed above are as follows. ◦ Region I . It corresponds to a standard slow roll phase. We assume η I = 0 and therefore we have a constant (cid:15) I (cid:28) . ◦ Region II . It corresponds to the USR phase, which extends from the the e -fold time N in to N end . We define ∆ N ≡ N end − N in and take η II (cid:29) / . Considering (cid:15) (cid:28) η, its evolution is (cid:15) II ( N ) = (cid:15) I e − η II ( N − N in ) . ◦ Region III . It corresponds to the last part of the inflationary dynamics. We assume a negative value for η III . (cid:15) evolves according to (cid:15)
III ( N ) = (cid:15) I e − η II ∆ N e − η III ( N − N end ) . In each region the evolution of the R k can be obtained analytically by solving the Mukhanov-Sasaki equation. Switchingto conformal time τ (we remind that dN e /dτ = a H ), for each region (i) we have R (i) k ( τ ) = H √ π (cid:112) (cid:15) (i) ( τ ) ( − τ ) / (cid:104) α (i) k H (1) ν i ( − kτ ) e i ( ν i +1 / π/ + β (i) k H (2) ν i ( − kτ ) e − i ( ν i +1 / π/ (cid:105) , (15)where H is taken constant since we always have (cid:15) (cid:28) . H (1) ν i and H (2) ν i are the Hankel function of first and second kindrespectively and ν i = (cid:112) / − η i (3 − η i ) . Imposing Bunch-Davies initial conditions, one has α I k = 1 and β I k = 0 in thefirst region. The complex coefficients α (i) k and β (i) k in the regions II and III are obtained requiring continuity of R k andits first derivative across the three regions. Finally, the power spectrum P R ( k ) is computed at super-Hubble scales as P R ( k ) ≡ lim kτ → − k π |R k ( τ ) | . We can now try to relate the parameters of the model to observable quantities (see table I). For this purpose in the leftpanel of fig. 3 we show the power spectrum of comoving curvature perturbations as a function of the comoving scale k for a choice of the parameters of the analytical model. The large enhancement at small scales allows for a sizablePBH production. As anticipated, the analytical model reproduces, at least qualitatively, the main features obtainedwith a more accurate numerical calculation, namely the presence of a peak and a dip in the power spectrum. Thecorrespondence between what we compute and what we observe can be explained with the following points. � �� � ��� � η ������ � �� ��� Δ � η = � η = � / � η = - � η = � η �� η ��� FIG. 2:
The blue dashed line shows the evolution of η as function of the number of e -folds in the model of ref. [27]. Specifically, withinthis solution (depicted with a cyan star in [27]) all the dark matter is in form of PBHs. The black line shows the behaviour of η in ourtoy-model. Model parameter H / π (cid:15) I N in ∆ N η II η III
Description P R in region I beginning of USR duration of USR value of η during USR value of η after USR Features in the normalization of P R k peak k dip , P R ( k peak ) k dip , P R ( k peak ) P R ( k ) power spectrum at large scales Observable A s M PBH n s f PBH total number of e -foldsTABLE I: Correspondence between the parameters of the model introduced in section II and the cosmological observables relevant forPBH formation. i) We fix the position of the peak of P R ( k ) at k peak = 10 Mpc − . This choice, in turn, gives a PBH distribu-tion peaked around M PBH = 10 g. These values are motivated by the fact that such population of PBHs cansuccessfully account for the totality of dark matter in the present-day Universe. Furthermore, from the value M PBH = 10 g we know that the comoving wavenumber k peak crosses the horizon approximately e -foldsafter the perturbation associated to the pivot wavenumber associated with the CMB scale k ∗ = 0 . Mpc − doesit. From that, we infer the value of N peak (and N in ). ii) For a given value of ∆ N and η II , we can compute the position of the dip and, consequently, the number of e -folds between horizon crossing of the pivot scale k ∗ and k dip . We can, therefore, compute the normalization ofthe power spectrum H / π (cid:15) I by rescaling the observed value at the pivot scale, namely A s (cid:39) × − , up to k dip by means of the slow-roll evolution P R ( k ) = A s ( k/aH ) n s − where we take n s = 0 . . iii) Finally, the values of ∆ N , η II and η III are readjusted–recomputing for each choice the normalization describedat point ii) –by requiring f PBH = 1 . The value of η III mostly controls the dynamics after the USR phase (althoughit also enters in the computation of the power spectrum), and its value can be gauged in order to guarantee thatinflation lasts approximately 60 e -folds.Typical values of ∆ N , η II and η III which generate the correct PBH abundance of dark matter are shown in fig. 3. Despitethe simplicity of the model, it is nice to see that these numbers are in full agreement with those obtained, by means ofa careful numerical analysis, in the context of more complete models [27].However, this discussion makes also clear that the picture risks being grossly incomplete without including the effectof non-gaussianities. This is mostly because the integral in eq. (8) is dominated by the tail of the distribution where �� �� �� �� �� �� �� �� ��� - � �� - � �� - � �� - �� [ ��� - � ] ℛ ( ) η �� = ���� η ��� = - ���� Δ � = ���� �� �� �� �� �� �� �� �� �� - � �� - � �� - � �� - � �� - � �� - � � � ��� [ � ] � � � � ≡ Ω � � � / Ω � � � � � � � � � � � � � � � � � � � � �� �� ����� ������ ������������ �������������� ���������� FIG. 3:
Power spectrum of comoving curvature perturbations and PBH abundance in the gaussian approximation. Both quantitiesrefer to the toy analytical model. The abundance is computed by means of threshold statistics, eq. (9). We compare the abundance ofPBH with existing bound from Hawking evaporation using both the extragalactic background radiation (EGBR, ref. [29]) and the 511keV gamma-ray line (keV, ref. [30]) and micro-lensing constraints from observation of the Andromeda galaxy M31 (HSC M31, ref. [31]).We also show future detection prospects using femto-lensing of gamma-ray bursts (FL, ref. [32]). For completeness, we also show boundsfrom neutron star disruption (NS, ref. [33]; see also refs. [34–36] for a more detailed discussion), and white dwarf explosions (WD,ref. [37]). deviation from the gaussian approximation may give a sizable effect.On general ground, deviations from the gaussian approximation are expected in realistic models of inflation. At CMBscales, non-gaussianities are typically negligible since suppressed by the small values of the Hubble parameters (cid:15) and η (see ref. [38] for the latest experimental constraints on primordial non-gaussianities at the scales probed by the Planckexperiment). However, we are interested in the opposite situation in which slow-roll is violated, and in particular η takes sizable values. It is, therefore, plausible to expect sizable non-gaussianities at scales much smaller that thoseprobed by CMB measurements where the USR dynamics might take place.Deviation from normality can be explored by considering skewness and kurtosis. Skewness ( S ) and kurtosis ( K )measure deviations with respect to a gaussian normal distribution. Gaussian distributions are symmetric around theirmean while, for a generic distribution, skewness measures the degree of departure from the symmetric shape. Kurtosis,on the contrary, is a measure of how differently shaped are the tails of a distribution with respect to the tails of thegaussian distribution. In the presence of non-gaussianities, eq. (8) is modified as [39, 40] β N G ( M PBH ) = β G ( M PBH ) + γ e − δ / σ R √ π (cid:20) S R / H (cid:18) δ th √ σ R (cid:19) + K R H (cid:18) δ th √ σ R (cid:19) + . . . (cid:21) , (16)which can be considered as an application of the Gram-Charlier series, that is the expansion of the probability den-sity function (PDF) in terms of its cumulant (sometimes also called Edgeworth expansion, the latter being the sameseries but with a different truncation and ordering of terms). In eq. (16), S R and K R are the skewness and the kurtosissmoothed on the scale R (in analogy with eq. (7)) and their explicit definitions will be discussed in the next sections. Inshort, skewness is related to the three-point correlation function which vanishes in the gaussian limit (see, e.g., eq. (3))while kurtosis is related to the connected part of the four-point correlation function (that is the part of the four-pointcorrelator that remains after subtracting the gaussian contribution, see, e.g., eq. (4)). The dots in eq. (8) indicate contri-butions from higher order correlators ( n (cid:62) ) and H n are the Hermite polynomials. In the following, we shall proceedwith the discussion of non-gaussianities assuming eq. (16). In the next section we will compute the skewness and in sec-tion IV we will focus on the kurtosis. With these results, in section V we will investigate the impact of non-gaussianitieson the abundance of PBHs. Notice that we use the so-called “physicists’ Hermite polynomials” which are defined by H n ( x ) = ( − n e x d n e − x /dx n instead of the so-called“probabilists’ Hermite polynomials” He n ( x ) often preferred in statistics theory. The relation among the two is He n ( x ) = 2 − n/ H n ( x ) . III. THE BISPECTRUM AND THE SKEWNESS
The three-point correlation function in momentum space is (cid:104)R k R k R k (cid:105) = (2 π ) δ (3) ( (cid:126)k + (cid:126)k + (cid:126)k ) B R ( k , k , k ) , (cid:126)k (cid:126)k (cid:126)k (17)where B R ( k , k , k ) defines the bispectrum. In the gaussian approximation, the bispectrum vanishes. While thepower spectrum is a function of k only, the bispectrum is a three-dimensional function defined on a tetrahedral re-gion for which (cid:126)k , (cid:126)k and (cid:126)k satisfy the triangle condition enforced–as illustrated in the schematic picture above–by thedelta function in eq. (17). It is convenient to define the reduced bispectrum B R ( k , k , k ) by means of B R ( k , k , k ) = B R ( k , k , k ) [∆ R ( k )∆ R ( k ) + ∆ R ( k )∆ R ( k ) + ∆ R ( k )∆ R ( k )] . (18)This is a useful definition since the reduced bispectrum B R is dimensionless (roughly speaking, the scaling B R ∼ /k is absorbed by ∆ R ∼ /k ). In position space, the skewness is S = 1( σ ) / (cid:90) d (cid:126)k (2 π ) d (cid:126)k (2 π ) d (cid:126)k (2 π ) (cid:104)R k R k R k (cid:105) = 1( σ ) / (cid:90) d (cid:126)k (2 π ) d (cid:126)k (2 π ) B R ( k , k , k ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:126)k + (cid:126)k + (cid:126)k =0 . (19)We take this opportunity to give a more general definition. We define the n -th cumulant C ( n ) ≡ (cid:104) n times (cid:122) (cid:125)(cid:124) (cid:123) R ( (cid:126)x ) . . . R ( (cid:126)x ) (cid:105) ( σ ) n/ = (cid:104)R n (cid:105) σ n , (20)as the connected part of the n -point correlator (evaluated at the same point) normalized by the n -th power of thestandard deviation. The connected part of the n -point correlator is obtained by subtracting from the full expression ofthe n -point correlator the gaussian contribution (which is zero for the 3-point correlator but, as we shall see in the nextsection, does not vanish for n (cid:62) , see e.g. eq. (4)). Given the definition in eq. (20), we have C (3) = S .The integration in R is limited by the tetrahedral region illustrated below (with semi-perimeter K ≡ ( k + k + k ) / ) V B ≡ k = K (1 − α + β ) / k = K (1 + α + β ) / k = K (1 − β ) (cid:126)k (cid:126)k (cid:126)k K k = k , k = k = k , k = αβ with slices β = 1 β = α = 0 α = − α = 1 k = k = K , k = 0 k = k = Kk = 0 k = k = Kk = 0 constant K (21)and ranges K ∈ [0 , ∞ ) , α ∈ [ − (1 − β ) , − β ] and β ∈ [0 , . Qualitatively, if we write (assuming P R and B R arescale-independent) ( σ ) / ∼ P / R and B R ∼ B R P R , we see that one expects S ∼ B R P / R . As a rule of thumb, weexpect an almost gaussian perturbation if B R (cid:28) P − / R . This means that B R is an interesting quantity to compute tohave at least a first idea about the impact of the three-point correlator on non-gaussianities. To perform the integral ineq. (19) some care is needed for the angular part since the triangular condition enforced by the delta function involvesalso angles. We use the exponential integral form of the delta function, by means of which we write (cid:90) d (cid:126)k (2 π ) d (cid:126)k (2 π ) d (cid:126)k (2 π ) (2 π ) δ (3) ( (cid:126)k + (cid:126)k + (cid:126)k ) = 12 π (cid:90) d log x (cid:90) dk dk dk ( k k k ) sin( k x ) sin( k x ) sin( k x )= 18 π (cid:90) V B dk dk dk ( k k k ) , (22) The delta function corresponds to invariance under translations, and the fact that the bispectrum depends only on the lengths of the three sides ofthe triangle formed by the momenta corresponds to invariance under rotations. Homogeneity and isotropy, therefore, reduce the nine parametersof the triad (cid:126)k , , down to three. where, crucially, the second equality in which we performed the integral over x is valid only inside the tetrahedrondefined in eq. (21) (as one can check by using the parametrization of k , , given in eq. (21) and performing the integralover x after imposing the correct constraints for α and β ). Of course, this way of rewriting the angular integral worksonly because we are assuming isotropy for the bispectrum (i.e. the omitted integrand in eq. (22) does not depend onangles).In analogy with eq. (7), we define the smoothed skewness as S R = 1( σ R ) / (cid:90) d (cid:126)k (2 π ) d (cid:126)k (2 π ) d (cid:126)k (2 π ) W ( k R ) W ( k R ) W ( k R ) 8(1 + ω ) (5 + 3 ω ) (cid:18) k aH (cid:19) (cid:18) k aH (cid:19) (cid:18) k aH (cid:19) (cid:104)R k R k R k (cid:105) = 1( σ R ) / π (cid:90) V B dk dk dk ( k k k ) W ( k R ) W ( k R ) W ( k R ) (cid:18) k aH (cid:19) (cid:18) k aH (cid:19) (cid:18) k aH (cid:19) B R ( k , k , k ) (cid:124) (cid:123)(cid:122) (cid:125) ≡(cid:104) δ (cid:105) R , (23)computed in the region defined by eq. (21). We define the n -th cumulant (in analogy with eq. (20) but now in terms of δ , and including the smoothing) C ( n ) R ≡ (cid:104) n times (cid:122) (cid:125)(cid:124) (cid:123) δ ( (cid:126)x ) . . . δ ( (cid:126)x ) (cid:105) ( σ R ) n/ = (cid:104) δ n (cid:105) R σ nR . (24)By using eq. (23), the goal is to compute the skewness as function of the PBH mass. To this end, we need to compute thethree-point function (cid:104)R k R k R k (cid:105) in Fourier space. Notice that, in order to compute the integral in eq. (23) withoutany approximation, we would like to keep the full shape in momentum space of B R ( k , k , k ) into account.We compute (cid:104)R k R k R k (cid:105) by means of the standard background+perturbations splitting and using the standard in-in formalism. We report directly the results and refer to appendix A 1 for details. One can isolate two dominant con-tributions in the calculation of three-point function. The first one (the sum of expressions in and end below) originatesfrom the following operator in the cubic action of R : S ⊃ − (cid:82) dτ d x a (cid:15) η (cid:48) R R (cid:48) . The second contribution (dubbed red ) originates by a convenient redefinition of R , see eq. (A7).We consider first the analytical model discussed in section II. Our result for B R ( k , k , k ) is B R ( k , k , k ) = B (in) R ( k , k , k ) + B (end) R ( k , k , k ) + B (red) R ( k , k , k ) (cid:124) (cid:123)(cid:122) (cid:125) local , (25)where k B (in) R = − Im H (cid:15) Γ( ν III ) π / e η II − η III )∆ N / − ν III ( x x x ) η III − η II R F in ( x , x , x ) , (26) R ≡ (cid:89) i =1 (cid:104) α III k i e iπ ( ν III + ) − β III k i e − iπ ( ν III + ) (cid:105) , (27) k B (end) R = − Im H (cid:15) Γ( ν III ) π / e η II − η III )∆ N / − ν III ( x x x ) η III − ( − η II + η III ) e N (1 − η II ) R F end ( x , x , x ) , (28) k B (red) R = − π η III ( x x x ) (cid:2) x P R ( k ) P R ( k ) + x P R ( k ) P R ( k ) + x P R ( k ) P R ( k ) (cid:3) . (29)In the expressions above we defined x i ≡ k i /k in , with k in the comoving wavenumber associated to the start of the USRphase at conformal time τ in by the relation k in τ in = − . The two functions F in ( x , x , x ) and F end ( x , x , x ) aredefined by the general expression F i ( x , x , x ) ≡ i (cid:15) / k / H (cid:18) R ∗ k R ∗ k d R ∗ k dτ + R ∗ k R ∗ k d R ∗ k dτ + R ∗ k R ∗ k d R ∗ k dτ (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) τ = τ i , (30)and we have F in ( x , x , x ) = e − i ( x + x + x ) ( x x x ) / (cid:2) x ( x x − ix − ix −
1) + x ( x x − ix − ix −
1) + x ( x x − ix − ix − (cid:3) , (31)0while the corresponding expression for F end ( x , x , x ) is more complicated and must be evaluated numerically fornon-integers η II , III . Armed with these expressions, we are ready to compute the reduced bispectrum in eq. (18) andtackle the integral in eq. (23).In fig. 4 we show the bispectrum as function of the PBH mass. The results obtained in the context of the analyticalmodel are shown in the left panel of fig. 4. Notice that the ratio S R = (cid:104) δ (cid:105) R /σ R is proportional to the factor A / s = �� �� �� �� �� �� �� �� �� �� �� - � �� - � � �� � �� � �� � � ��� [ � ] � � ≡ � [ � � / � ] η �� = ���� η ��� = - ���� Δ � = ���� � ℛ � � � � ℛ � �� ���������� ����� �� �� �� �� �� �� �� �� �� �� �� - � �� - � � � ��� [ � ] � [ ℴ - η � ] ��������� ����� � ℛ ��� FIG. 4:
Smoothed skewness in eq. (23) as function of the PBH mass computed for the model in which f PBH = 1 in the gaussian approx-imation (see fig. 3). The yellow line represents the sum of the three contributions in eq. (25). Separately, we also show the contributionfrom B (end) R and B (red) R (that is the local contribution); B (in) R gives a negligible contribution not shown in the plot. The right-side ofthe y -axis shows the value of S R in units of ( A s ) / with A s = H / π (cid:15) I . On the left-side y -axis, A s is fixed by the normalizationchosen in fig. 3. Right panel. Smoothed skewness in eq. (23) as function of the PBH mass computed for the numerical model in ref. [27].In the numerical model, only the contribution B (red) R survives (solid yellow line). H/ √ π (cid:15) I . On the left-side of the y -axis, we fix this value by means of the procedure outlined in the introduction ofsection II (see table I and related discussions). On the right-side y -axis, on the contrary, we indicate the size of the ratio S R /A / s such that this quantity is not affected by the details of the fit of the power spectrum at large scales (where ouranalytical model is not accurate since (cid:15) I is taken to be constant).We find that B (in) R is negligible while B (end) R gives the dominant contribution compared to the local term B (red) R . Inthe analytical model, the two contributions B (in) R and B (end) R come, respectively, from the two delta function transitionin the time evolution of η dηdτ = η II δ ( τ − τ in ) + ( − η II + η III ) δ ( τ − τ end ) . (32)Modeling the transition at τ = τ end as a sharp step function can be misleading. This was already noticed in ref. [41](even though in the context of non-attractor inflation models) and further analyzed in ref. [42]. The outcome of theanalysis carried out in ref. [41] is that if one takes—in place of the sharp step function—a smooth transition at τ end , thecontribution B (end) R gets drastically reduced.The evolution of η was already shown in fig. 2 where it is indeed clear that the step-function approximation at τ end is the crude approximation of a smoother transition. When the smooth transition is implemented, we indeed find that B (end) R is completely negligible in the computation of the bispectrum. In short, this is because in the presence of asmooth transition B (end) R receives—instead of the instantaneous contribution of the delta function—two integratedcontributions (respectively, before and after the transition at τ end ) with opposite signs that almost cancel betweeneach other (at the ‰ level). Let us give here just the idea of the proof (which is different compared to the one given inref. [41]). Instead of the step-function approximation, we introduce the hyperbolic tangent parametrization η ( N ) = 12 (cid:20) − η II + η II tanh (cid:18) N − N in δN (cid:19)(cid:21) + 12 (cid:20) η II + η III + ( η III − η II ) tanh (cid:18) N − N end δN (cid:19)(cid:21) , (33)1 � �� � ��� η �� � η ��� � η ℐ ���� ℐ ����� ��� ��� ��� ����� - � ����� � δ � ( ℐ � � �� + ℐ � � � � � ) / ℐ � � �� [ % ] ℯ - ℴ ℓ ℯ ℊδ � → � FIG. 5:
Left panel. Hyperbolic tangent parametrization of η in eq. (33) for different δN (from darker to lighter red, δN = 0 . ÷ . in steps of . ). The black line is the step-function limit δN → . The two integrals I left and I right refer to the intervals [ N end − . , N end ] and [ N end , N end + 1 . , respectively. Right panel. Sum of the two integrals I left + I right (normalized to I left , and expressedin percentage) as function of δN . When the ratio ( I left + I right ) / I left is above 100%, the two integrals I left and I right have the samesign and they sum constructively. Below this value, a cancellation occurs with increasing precision. where the parameter δN controls the width of the two transitions at N in and N end . This parametrization can be con-sidered as a proxy for the actual numerical result (see fig. 2). On the other hand, the limit δN → reproduces the step-function approximation. Notice that we also introduced an hyperbolic tangent to model the transition at N = N in ,although not strictly needed for the validity of our considerations at N = N end . The hyperbolic tangent parametriza-tion in eq. (33) is shown in the left panel of fig. 5 for different values of δN . Consider the transition at N = N end . Inthe cubic action − (cid:82) dτ d x a (cid:15) η (cid:48) R R (cid:48) we now have to integrate over time instead of just picking up a delta functionterm at N = N end . One ends up, as anticipated before, with the sum of two integrals (respectively, I left and I right ) cov-ering the two regions indicated in the left panel of fig. 5. For small δN , when the transition is very sharp, the two terms I left and I right have the same sign, and their sum reconstructs the large bispectrum found in eq. (28). For larger δN , I right flips sign and the cancellation occurs with increasing precision. In the right panel of fig. 5 we show the quantity ( I left + I right ) / I left as function of the width of the transition δN . For δN (cid:38) . , the sum I left + I right quickly gives acancellation with the ‰ accuracy.In conclusion, we find that the analytical toy model in which the evolution of η in the presence of USR is modeled bymeans of a step-function approximation is not best-suited for the computation of non-gaussianities since it misses theimportant cancellation that is at work when the transition between adjacent regions is implemented in a more physicalway. For this reason, we switch from now on to the numerical model introduced in section II. By direct computation,we confirm that B (end) R is negligible in the numerical model, where we verified the same cancellation discussed abovefor the hyperbolic tangent parametrization. This is because the evolution of η in the numerical model can be indeedapproximated very well by eq. (33) with width parameter δN (cid:39) . (for which the cancellation is optimal, see fig. 5).As a final comment before proceeding our analysis, notice that the cancellation of B (end) R is not an exact mathematicalidentity but, as discussed before, its precision depends on the specific value of δN . In the numerical model introducedin section II, the latter happens to be large enough that one can safely consider B (end) R (cid:39) when compared to B (red) R .However, if one takes a model in which, say, δN (cid:39) . then our result shows that the cancellation works only at the ∼ % level; this means that B (end) R gets reduced by, say, a factor of 2 but remains very sizable. Since the momentumstructure of B (end) R is not of local type, in such a situation it would become interesting to include it. The relevant ques-tion becomes the following: Is a model with δN (cid:39) . realistic? To answer this question, we can reverse-engineeringthe inflaton equation of motion and get V ( N ) = V ( N ref ) exp (cid:40) − (cid:90) NN ref dN (cid:48) (cid:20) (cid:15) (3 − η )3 − (cid:15) (cid:21)(cid:41) , (34) φ ( N ) = φ ( N ref ) ± (cid:90) NN ref dN (cid:48) √ (cid:15) , (35)as far as inflaton potential and field profile are concerned. These two equations are exact. We can now use the2parametrization of η given in eq. (33) and derive the time-evolution of (cid:15) . Using η (cid:39) − (1 / d log (cid:15)/dN , we find thefollowing (cumbersome but analytical) expression (cid:15) ( N ) = (cid:15) I e − η III ( N − N ref ) (cid:20) cosh (cid:18) N − N end δN (cid:19) cosh (cid:18) N − N in δN (cid:19)(cid:21) − δNη III2 (cid:20) cosh (cid:18) N ref − N end δN (cid:19) cosh (cid:18) N ref − N in δN (cid:19)(cid:21) δNη III2 (cid:20) cosh (cid:18) N − N end δN (cid:19) sech (cid:18) N − N in δN (cid:19)(cid:21) δN ( η II − η III2 ) (cid:20) cosh (cid:18) N ref − N end δN (cid:19) sech (cid:18) N ref − N in δN (cid:19)(cid:21) δN ( − η II + η III ) , (36)where (cid:15) I (cid:28) is the value of (cid:15) at some reference N ref in region I. We now integrate eqs. (34, 35) and get the inflatonpotential as function of φ for different δN . We show our result in fig. 6. We find that the width of the transition region ϕ - ϕ ��� � ( ϕ ) / � ( ϕ � � � ) δ � = � δ � = ���� δ � = ��� η �� = ���� η ��� = - ���� Δ � = � FIG. 6:
Solution of eqs. (34, 35) with η -evolution given in eq. (33) and (cid:15) -evolution given in eq. (36) for different δN with φ ref ≡ φ ( N ref ) . of η is related to the deformation of the stationary inflection point of the potential. Small values δN → deepenthe local minimum and lead to a situation in which the inflaton risks being trapped into it (thus preventing a gracefulexit). This is clear as far as the limit δN → is concerned. One is tempted to conclude that the same result holds truewhen comparing the cases with δN = 0 . and δN = 0 . . However, any intuition must be taken with a grain of saltsince we know that small changes in the way one deforms the inflection point from its stationary configuration havea huge impact on the amplitude of the power spectrum (see discussion in ref. [27]). For this reason, we postpone aquantitative estimate of the relation between the power spectrum and δN to the end of this section, and we move onwith our discussion of non-gaussianities in the context of the numerical model in its benchmark configuration where B (end) R can be neglected.The smoothed skewness in the context of the numerical model is shown in the right panel of fig. 4. The contribution B (red) R survives also in the numerical model, and it is shown in yellow in the right panel of fig. 4. In the numericalmodel, the contribution B (red) R is evaluated at some later time after the end of the USR phase when the relevant modescease evolving (that is region III in the analytical model). We indicate with η (which is de facto equivalent to η III ) the(negative) value of η at this transition time. Since η enters as an overall multiplicative factor in the computation ofthe local contribution of the bispectrum (the dominant one, as explained above), in the right panel of fig. 4 we plot thesmoothed skewness in units of − η , that is S R / ( − η ) . As far as this contribution is concerned, the numerical modelgives a result that is in good agreement with the analytical model (despite the poor description of the latter at largescales).In conclusion, we find that the presence of USR leaves a sizable imprint on the bispectrum of comoving curvatureperturbations (and, consequently, on the density contrast) at the typical scales which are relevant for PBH formation.Non-gaussianities are predominantly of local type, meaning that the bispectrum of comoving curvature perturbationsin Fourier space has the form (see eq. (29)) B loc R ( k , k , k ) = (6 f NL / R ( k )∆ R ( k ) + ∆ R ( k )∆ R ( k ) + ∆ R ( k )∆ R ( k )] . (37) Locality here is just the statement that if the curvature is evaluated at some point (cid:126)x then the non-gaussianity is localized at same position. To beclear, a non-local type of non-gaussianity takes the form R ( (cid:126)x ) = R G ( (cid:126)x ) + (cid:82) d (cid:126)yd (cid:126)zK ( (cid:126)x ; (cid:126)y, (cid:126)z )[ R G ( (cid:126)y ) R G ( (cid:126)z ) − (cid:104)R G ( (cid:126)y ) R G ( (cid:126)z ) (cid:105) ] with the localcase corresponding to the kernel K ( (cid:126)x ; (cid:126)y, (cid:126)z ) = 3 f NL / δ (3) ( (cid:126)y − (cid:126)x ) δ (3) ( (cid:126)z − (cid:126)x ) . f NL with its historical normalization that became conventional, which for theclass of models that we are studying reads f NL = − / η .The presence of USR dynamics is crucial to generate potentially sizable non-gaussianities. Conceptually, the reasonis the following. The modes that contribute to the formation of the peak in the power spectrum of comoving curvatureperturbations are those that are (in modulus) enhanced (or suppressed, as far as the modes that contribute to thedip of the power spectrum are concerned) by the negative friction phase that takes place during USR. These modescease evolving only after the end of USR, a handful of e -folds before the end of inflation. During this last part of theinflationary dynamics, the Hubble parameter η takes O (1) negative values which are needed to connect the end ofthe USR phase with the end of inflation. In fact the inflaton almost stops on the top of the local maximum that formswhen an approximate stationary inflection point is present in the potential, and accelerates afterwards to reach thesubsequent absolute minimum where inflation ends. Since we are forced by the presence of the USR dynamics toevaluate correlators for comoving scales k (cid:38) Mpc − in a region where | η | ∼ O (1) , the usual slow-roll suppressionwhich affects non-gaussianities at CMB scales is not at work anymore, thus opening the possibility of having potentiallylarge effects.Let us now close this section with a final comment related to the cancellation of B (end) R . Consider the scalar potentialin eq. (14) that we write in the Einstein frame in the form V ( φ ) = λφ ξφ ) (cid:20) ξ φ − c ) φ φ + 2(1 + c )(3 + ξφ ) φ φ (cid:21) = ��� ��� ��� ��� ��������������������� ϕ � ! � ( ϕ ) / λ where φ = 1 indicates the position of the inflection point and c , describe deformations from its stationary con-figurations (that corresponds to c , = 0 ). Changing the value of the non-minimal coupling (gray region in the plotabove) keeping c , (cid:54) = 0 fixed does not alter the shape of the potential around the inflection point. This means that allthese solutions will be described by a similar value of δN. However, we find that a ‰ variation of ξ changes the peakof the power spectrum by many orders of magnitude. The smaller ξ the lower P R ( k peak ) . In practice, in the context ofthe numerical model, the peak of the power spectrum is controlled by the size of the non-minimal coupling ξ. A smalltuning of ξ leaves δN almost unchanged but allows to select the desired value of P R ( k peak ) . We will see this explicitlyin section V.On the other hand, it is clear that changing the values of c , directly alters, by construction, the shape of the approx-imate stationary inflection point (red curves in the plot above, see ref. [27]). For instance, we found good inflationarysolutions (that is solutions of the inflaton equation of motion in which the inflaton successfully overcomes the inflec-tion point and reaches the absolute minimum at the origin) with δN (cid:39) . . According to the right panel of fig. 5, thisis enough to get a cancellation that works only at the 10% level. In such a case, neglecting B (end) R is not justifiable. Ofcourse, more work is needed to validate these solutions against CMB data and bound from Hawking evaporation (inparticular, all solutions found with δN (cid:39) . tend to generate PBHs that are too light to be dark matter). The pointthat we want to make here is that the condition B (end) R (cid:39) is not guaranteed but one should always confront the η -evolution of a particular solution with our result in the right panel of our fig. 5 before drawing any conclusion aboutthe size of non-gaussianities generated by the transition at N = N end . The recipe to follow is simple. Compute thetime-evolution of η by solving the inflaton equation of motion and compare it with eq. (33) to extract the value of δN that best suits the shape of the transition at N = N end . From the right panel of fig. 5, one then reads the accuracy ofthe cancellation.In the rest of this work we will focus on cases with B (end) R (cid:39) . Notice that here φ is not canonically normalized but this fact does not alter the validity of the qualitative argument below. See ref. [27] for details.We also neglect here for simplicity the presence of higher-dimensional operators. IV. THE TRISPECTRUM AND THE KURTOSIS
The four-point correlation function in momentum space has the structure (cid:104)R k R k R k R k (cid:105) = (cid:104) (2 π ) δ (3) ( (cid:126)k + (cid:126)k ) δ (3) ( (cid:126)k + (cid:126)k )∆ R ( k )∆ R ( k ) + two permutations (cid:105)(cid:124) (cid:123)(cid:122) (cid:125) disconnected contribution +(2 π ) δ (3) ( (cid:126)k + (cid:126)k + (cid:126)k + (cid:126)k ) T R ( k , k , k , k , k , k ) (cid:124) (cid:123)(cid:122) (cid:125) connected contribution (cid:126)k (cid:126)k (cid:126)k (cid:126)k (38)The trispectrum T R depends on the six scalars specifying the quadrilateral formed by the wave-vectors, as illustrated inthe schematic picture above where the two dashed lines correspond to k ≡ | (cid:126)k − (cid:126)k | and k ≡ | (cid:126)k − (cid:126)k | . Notice thatin the computation of (cid:104)R k R k R k R k (cid:105) we need to isolate the connected part since the disconnected part is alreadypresent at the gaussian level, see eq. (4). In analogy with the reduced bispectrum introduced in eq. (18), it is useful todefine the reduced trispectrum T R which is given by T R ( k , k , k , k , k , k ) = T R (cid:110) ∆ R ( k )∆ R ( k ) (cid:104) ∆ R ( | (cid:126)k + (cid:126)k | ) + ∆ R ( | (cid:126)k + (cid:126)k | ) (cid:105) + 11 permutations (cid:111) . (39)As done in the previous section, we define the smoothed kurtosis as K R = 1( σ R ) (cid:90) (cid:89) i =1 d (cid:126)k i (2 π ) W ( k i R ) 16(1 + ω ) (5 + 3 ω ) (cid:18) k i aH (cid:19) (cid:104)R k R k R k R k (cid:105) (cid:124) (cid:123)(cid:122) (cid:125) ≡(cid:104) δ (cid:105) R . (40)The exact computation of K R is beyond the scope of this paper. However, an important remark is in order. When non-gaussianities are of local type, they generate—in addition to the bispectrum of the form defined by eq. (37) and obtainedin eq. (29) as a consequence of our explicit computation—also a contribution to the reduced trispectrum in eq. (39). Forthe sake of simplicity, we ignore in the following the smoothing procedure and the relation between density contrastand curvature perturbations, and we just consider the curvature perturbation R . It is a simple exercise to show that if wetake R = R G + (3 f NL / (cid:0) R G − (cid:104)R G (cid:105) (cid:1) we find a trispectrum of the form given in eq. (39) with T R = 18 f / [43, 44].This can be understood in the context of the in-in formalism as shown in appendix A 2 (see derivation of eq. (A37)).We do not attempt here to compute K R carefully but we just limit the discussion to the simple (and well-known) ob-servation that a non-zero local-type bispectrum of order O ( f NL ) also implies a non-zero contribution to the reducedtrispectrum of order O ( f ) . Statistically speaking, this means that we expect a non-zero fourth-order cumulant re-lated to the third-order one. This comment actually applies also to higher-order cumulants as well.Of course, we stress again that in addition to this O ( f ) contribution other terms are expected which are generated,for instance, by a pure quartic interaction. In the following, we shall work under the assumption that all these additionalterms are negligible, and that the dominant contribution to the four-point correlator is the one generated by local cubicinteractions. V. ON THE PBH ABUNDANCE IN THE PRESENCE OF NON-GAUSSIANITIES: THE LOCAL CASE.
We are now in the position, using eq. (16), to analyze the impact of non-gaussianities for the PBH abundance. Tomake our discussion more clear, let us start by considering the PDF for the fractional overdensity δ whose integralabove the threshold gives the mass fraction β (see eq. (8)). We indicate with P N G ( δ ) the generic non-gaussian PDF, andwith P G ( δ ) the PDF in the gaussian approximation. As anticipated, eq. (16) follows from the Gram-Charlier series ofthe PDF P N G ( δ ) = P G ( δ ) ∞ (cid:88) n =0 n !2 n/ B n [0 , , C (3) R , . . . , C ( n ) R ] H n (cid:18) ν √ (cid:19) (41) = P G ( δ ) (cid:26) C (3) R / H (cid:18) ν √ (cid:19) + C (4) R / H (cid:18) ν √ (cid:19) + C (5) R / H (cid:18) ν √ (cid:19) + [10 C (3) 2 R + C (6) R ]6!2 / H (cid:18) ν √ (cid:19) + . . . (cid:27) , In general, T R depends on four three-momenta, namely twelve parameters; however, as discussed before in the case of the bispectrum, assumingisotropy and homogeneity the number of parameters reduces to six. ν ≡ δ/σ and where we introduced the n -th complete exponential Bell polynomial B n ( x , . . . , x n ) = n (cid:88) k =1 B n,k ( x , . . . , x n − k +1 ) , (42)with B n,k the partial exponential Bell polynomials. C ( n ) R is the n -th cumulant defined in eq. (24). In eq. (41), we carriedout explicitly the expansion up to the sixth order to show that, from this order on, the coefficient of the n -th term isnot just given by C ( n ) R . This is a well-known result in the theory of statistics [45, 46]. Without entering to deep into themathematical detail, eq. (41) can be thought as a Taylor expansion of P N G ( δ ) around a reference PDF, taken to be thegaussian distribution P G ( δ ) . For completeness, if we integrate eq. (41) to compute β N G = γ (cid:82) ∞ δ th dδP N G ( δ ) and usethe property √ πσ (cid:90) ∞ δ th dδe − δ / σ H n (cid:18) δ √ σ (cid:19) = 1 √ π e − δ / σ H n − (cid:18) δ th √ σ (cid:19) , (44)we find precisely eq. (16). From this simple discussion, we therefore learn that eq. (16) follows from a series expansionaround the gaussian PDF.In principle, all the cumulants are needed to compute the PBH abundance with eq. (16). However, to examine therelevance of non-gaussianities, we first start focusing only on the third order cumulant C (3) R , which we computed insection III (see fig. 4), and set the remaining ones to zero. A complete calculation is performed later in this section. Weshow our results in the left panel of fig. 7. The analysis is performed in the context of the numerical model, and wefocus on the PBH mass corresponding to the peak of the mass distribution (i.e. the most abundant PBHs). The dottedline (labelled GC ) shows the PBH abundance computed integrating eq. (41) above the threshold as a function of thepower spectrum of curvature perturbations at the peak position. Let us stress again that in eq. (41) C (3) R does not enteronly at third order in the expansion, but also affects higher ones. For definiteness, we include terms up to n = 15 ineq. (41). This means that C (3) R enters at n = 3 , , , , . We do not attempt here to include more terms and test theconvergence of eq. (41) since, as we shall discuss in a moment, we will eventually adopt a different approach for thecomputation of the abundance (see [26] for a critical discussion on the truncation of the Gram-Charlier serie). On thetop x -axis we also show the relative change of the non-minimal coupling (with respect to the value ξ f PBH =1 that givesthe totality of dark matter in PBHs for the gaussian model) that is needed to modify the peak of the power spectrumaccording to the corresponding point on the bottom x -axis. The dashed black line corresponds instead to the gaussianapproximation in eq. 8. Two main messages can be extracted from this plot. Non-gaussianities are relevant, and modifythe abundance computed with gaussian statistics by several orders of magnitude. On the other hand, the abundance isvery sensitive to the amplitude of the power spectrum, and it is enough to change the latter quantity by a small amountto obtain the same abundance predicted by the gaussian approximation.As mentioned above, this calculation is incomplete, because it misses the contributions from higher orders in theexpansion and from the cumulants C ( n ) R with n > (which are non-zero because, as discussed in section IV for thespecific case of the fourth-order cumulant, they at least receive a non-zero “induced” contribution from the local three-point function). Baring cancellations among the different contributions, one would expect that, at qualitative level,the results presented above would still hold, once all the other terms are included. In other words, the purpose of theestimate performed above is to highlight the relevance of non-gaussianities for the PBHs abundance. On the otherhand, for a proper calculation one could in principle use eq. (16), and compute all the cumulants needed to reachthe convergence of the serie at the desired accuracy. However, exploiting the results of section III, we follow a simplerstrategy, which will also allow us to take into account the non-gaussianities arising from the non-linear relation betweencurvature and density perturbations. The idea is the following. In section III we have shown that the three-pointcorrelator in the presence of USR generates non-gaussianities that are predominantly of local type, since they give abispectrum of the form given in eq. (37). This means that we can write in position space R = R G + 3 f NL (cid:0) R G − (cid:104)R G (cid:105) (cid:1) , (45) We remind that d n dx n e − x / σ = ( − n n/ σ n e − x / σ H n (cid:18) x √ σ (cid:19) , (43)so that the reader can immediately recognize in eq. (41) a derivative expansion. A similar strategy has been presented in [11] for the study of the non-gaussianities produced by the non-linear relation among density and curva-ture perturbations, see below. R G is gaussian and the term (cid:104)R G (cid:105) is added to ensure that (cid:104)R(cid:105) = 0 . The statistical properties of the non-gaussianvariable R are therefore completely specified in terms of those of the gaussian variable R G and the parameter f NL , which controls the amount of primordial non-gaussianities. As explained in sec. II, in threshold statistics the abun-dance of PBHs is obtained computing the probability that the density perturbation exceeds the threshold value δ th . Here, the strategy is to relate the density perturbation to the comoving curvature perturbation R and its spatial deriva-tives ∂ i R and (cid:52)R , where ∂ i denotes the derivative with respect to the spatial direction i. Using eq. (45), these non-gaussian variables can be expressed in terms of the analogous ones for the gaussian variable R G , whose statisticalproperties are known. In practice, exploiting the relation in eq.(45), we bypass the need of computing all the cumu-lants in eq. (16).To see explicitly how this procedure can be implemented, let us first define the following quantity: σ j = (cid:90) dkk W ( kR ) P R ( k ) k j , γ = σ σ σ . (46)For ease of notation, we also define the following gaussian variables: ν = R G σ , η i = ∂ i R G , x = − (cid:52)R G σ , (47)The joint PDF of the gaussian variables ν , (cid:126)η and x is [21]: P ( ν, (cid:126)η, x ) d ν d(cid:126)η dx = A e − Q dν d (cid:126)η dx , with Q = ν x − x ∗ ) − γ ) + 3 (cid:126)η · (cid:126)η σ , (48)and A = 6 √ π / (cid:112) − γ ) σ , x ∗ = γ ν. (49)To proceed, we should relate the comoving curvature perturbation with the density perturbation. As anticipated,this map is non-linear, and therefore it hides some unavoidable non-gaussianities. The full relation can be written inthe form [47] δ ( (cid:126)x, t ) = − ω )5 + 3 ω (cid:18) aH (cid:19) e − R ( (cid:126)x ) / (cid:52) e R ( (cid:126)x ) / = − ω )5 + 3 ω (cid:18) aH (cid:19) e − R ( (cid:126)x ) (cid:124) (cid:123)(cid:122) (cid:125) i ) (cid:20) (cid:52)R ( (cid:126)x ) + 12 ∂ i R ( (cid:126)x ) ∂ i R ( (cid:126)x ) (cid:124) (cid:123)(cid:122) (cid:125) ii ) (cid:21) , (50)We have two sources of non-linearities (and thus non-gaussianities) that come from i) terms of order higher than twoin the expansion of the exponential and ii) the quadratic term in the square brackets. Notice that eq. (50) is written incomoving gauge (defined by requiring both comoving slicing and comoving threading) whose choice seems appropri-ate since often used to perform simulations that describe the formation of PBHs in the context of numerical relativity.Furthermore, eq. (50) is not exact but includes terms up to the second order in the so-called spatial gradient expan-sion [47]. Notice that in eq. (50) R denotes the non-gaussian comoving curvature perturbation. Using eq. (45) one canrecast the expression above in terms of the gaussian variables ν , (cid:126)η, x in eq. (47). Then, solving for the variable x onefinds: x δ ( δ, ν, (cid:126)η ) = 9( aH ) e σ [ ν + ( ν − ) σ α ] δ + (cid:0) α + 8 νασ + 8 ν α σ (cid:1) (cid:126)η · (cid:126)η ν α σ ) σ . (51)Above we have defined α = 3 f NL / and taken w = 1 / , corresponding to the radiation domination era.Finally, the fraction of the Universe collapsing into PBHs is obtained integrating the joint PDF in eq. (48) rewrittenin terms of ν , (cid:126)η and δ, and imposing that density perturbation is above the threshold δ th : β N G = γ (cid:90) ∞ δ th dδ (cid:90) d (cid:126)η (cid:90) ∞−∞ dν J P ( ν, (cid:126)η, x δ ) , (52)where J is the Jacobian of the transformation in eq. (51): J = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) aH ) e σ [ ν + ( ν − ) σ α ]4 (1 + 2 ν α σ ) σ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (53)7Eq. (52) generalizes the gaussian expression in eq. (8) including the effect of primordial non-gaussianities of localtype, and of non-gaussianities induced by the non-linear relation in eq. (50). For the latter source of non-gaussianities,an analogous expression has been obtained in [11]. One can check that in the limit of f NL → and using eq. (50) atlinear order in R , eq. (52) reduces to eq. (8).We are now in position to quantify the impact of the non-gaussianities on the PBH abundance. We show in theleft panel of fig. 7 our results for β ( M PBH ) as a function of the power spectrum at peak position for both the gaussianand non-gaussian cases. As before, we consider the numerical model, for which η (cid:39) − . and therefore f NL = − / η (cid:39) . . To isolate the effect of primordial local non-gaussianites, we first perform the calculation employing �� - � �� - � �� - �� �� - �� �� - �� �� - �� �� - �� �� - �� �� - �� �� - �� �� - � �� - � � ξ � ��� = � - ���� ℛ ( ���� ) β ( � � � � ) Δξ [ � ] ℯℯ ℴℯℓℴℯ β β β + ℴℓ ��� ��� ��� ��� ������������������ - η � Δ ℛ ℊℴ ℊℯℴ ℊℯ + ℴ ℓℯℯ FIG. 7:
Left panel: PBH abundance as a function of the the power spectrum of curvature perturbations at peak position. The dashedblack line is for the gaussian approximation, eq. (8). The dotted black line refers to the abundance computed integrating the Gram-Charlier series in eq. (41), see text for details. The red line takes into account the effect of primordial non-gaussianities, and the greenline further includes the non-linearities between density and curvature perturbations. Right panel: Rescaling factor of the powerspectrum needed, in presence of non-gaussianities, in order to get the same PBH abundance predicted by the gaussian approximation.This quantity is shown as a function of − η , the parameter controlling the amount of local non-gaussianities. The color code of thelines is the same as in the left panel. The calculations in both panels are for the numerical model used in section III. the linear relation among density and the comoving curvature perturbations, eq. (6). The result is shown with a red linein fig. 7. Local non-gaussianities (with positive f NL ) tend to enhance the abundance of PBHs, and their effect can bereabsorbed by decreasing the amplitude of the power spectrum by a factor ÷ . This can be achieved at the prize ofa small re-tuning of the parameter that controls the height of the peak of the power spectrum. In the case discussed inthis paper, such parameter is the non-minimal coupling to gravity, and the re-tuning that one needs to gauge away theeffect of non-gaussianities is of ∼ O (10 − ) .We notice that using the Gram-Charlier series with the only inclusion of the third order cumulant C (3) R and few termsin the expansion, leads to an increase of the PBH abundance and therefore, at the qualitative level, it captures the maineffect of primordial non-gaussianities (in the limit in which non-linearities are neglected). However, we reiterate thatsome care is needed in the truncation of the serie, and a sufficient number of terms should be included in order toobtain a good convergence. We do not discuss further the Gram-Charlier series, and we focus on the approach that ledto eq. (52). A critical discussion about the validity of the approach based on the Gram-Charlier series can be found inref. [26].Then, we introduce in the calculation also the non-linearities among density and comoving curvature perturbations.Using eq. (52) we obtain the green line in fig. 7. We find that their effect is opposite to the one of primordial local non-gaussianities, making the formation of PBHs harder, in agreement with what found in [11] (but differently from thisreference we also include primordial non-gaussianities). The reason for that can be read from eq. (51): non-linearitiestends to increase the threshold x δ relevant for PBH formation. From fig. 7 one can see that, once again, the effect ofnon-gaussianities on the abundance can be reabsorbed by rescaling (in this case decreasing) the height of the power8spectrum by a modest amount.One could wonder how these results change for different values of η . To address this question, we show in the rightpanel of fig. 7, as a function of − η , the rescaling factor of the amplitude of power spectrum needed to get the same PBHabundance predicted by the gaussian approximation. For comparison, let us mention that the models of single fieldinflation with an USR phase in refs. [27, 48, 49], which are designed to obtained a sizable PBH abundance, predict − η (cid:39) . − . (see also [50] for a survey of different models in the literature). Our results show that non-gaussianities donot play a dramatic role for the computation of the abundance in this class of models, meaning that their effect can becompensated by a small tuning of the model parameters.Before concluding, we shall mention that in our analysis we focused on how non-gaussianties, primordial or fromthe non-linear relation, change the statistics of the density field. However, they also modify the threshold δ th . This isbecause δ th depends on the shape of the density perturbation, which is affected by the presence of non-gaussianties.These modifications has been studied in the literature following different approaches, focusing on the effect of non-linearities [10, 12, 13, 17, 18], or including also primordial non-gaussianities [15] (for the latter case the power spectrumwas assumed to be a Dirac-delta distribution; our power spectrum is significantly different). A complete investigationof the role of non-gaussianities should include also these effects. We leave it for future work. VI. SUMMARY AND CONCLUSIONS
Let us conclude with a brief summary of the main novelties and results of our work. We have focused on the studyof non-gaussianities which are generated by the presence of USR dynamics during inflation. The latter occur when theHubble parameter η takes values > / . An USR phase is crucial for the production of a sizable number of PBHs, andprevious studies often limited the computation of their abundance to the gaussian approximation.We have computed the bispectrum of the comoving curvature perturbations R in the context of the in-in formalism.From the operators appearing in the cubic action of R we have isolated two dominant terms. The contribution fromone of those drastically depends on the transition between the USR phase and the conventional slow-roll phase. Wehave shown that unless this transition is very sharp, a cancellation occurs in the calculation, dramatically reducing theimprint of this operator on the bispectrum. We have then offered a procedure, based on a simple analysis of the timeevolution of η , to estimate how sharp the transition should be in order for this cancellation to apply. Our results extendthe findings of ref. [41], which focused on a scenario of non-attractor inflation where the USR phase, taken with η = 3 , is connected to a phase with negligible η. On the contrary, in the context of models for PBHs production, the USR phase(generically with η (cid:54) = 3 ) is followed by a phase where η is negative and O (1) . This is a generic property in the presenceof USR, required to connect the end of the latter with the end of inflation.We have found that the remaining contributions to the bispectrum generate non-gaussianities predominantly of thelocal type, with f NL = 5 / − η ) . (54)where η is the value of the parameter η after the end of the USR phase. We worked out the details of our analysis bothsemi-analytically and numerically, in the context of the model introduced in eqs. (13, 14), in which the inflaton field isa real scalar singlet equipped with a generic polynomial potential and non-minimally coupled to gravity in the Jordanframe. For this class of models, we extend the results presented in [50]. Recently it has been claimed that at the end ofan USR phase non-gaussianities are slow-roll suppressed and thus negligible [51]. However, as we have shown with anexplicit example, this is not necessarily true. On the contrary, in the class of models that we have studied, and whichare motivated by the PBHs production, η is generically negative and O (1) . Then, using the formalism of threshold statistics, we have investigated the impact of non-gaussianities found here onthe PBH abundance. To evaluate their relevance, or in other words the validity of the gaussian approximation, we haveconsidered eq. (16), which allows to compute the PBH abundance once the cumulants of non-gaussian distributionare known. Focusing only on the first non-gaussian correction, i.e. the third-order cumulant C (3) R , we have found thatprimordial non-gaussianities sizeably modify the estimation of the PBH abundance based on gaussian statistics. Quali-tatively, the same conclusion has been obtained in [40, 50], where an approximation of eq. (16) was adopted to quantifythe impact of non-gaussianities. A quantitative comparison with these analysis and the approximations adopted therewill be presented in [26]. Then, we have taken a step forward and we have shown how to compute the PBH abundancein presence of local non-gaussianities, extending the gaussian threshold statistics in eq. (8). This leads to eq. (52). Inour calculations we have kept into account the non-linear relation between the density and curvature perturbations,which introduces another source of non-gaussianities in the addition to the primordial ones. A similar formalism hasbeen presented in [11], but focused only to the study of non-linearities. Using eq. (52), we have found that the impactof non-gaussianities can be reabsorbed by a change of the amplitude of the power spectrum of curvature perturbationsof a factor ÷ at most. Namely, after such shift one recovers the same PBH abundance obtained with the gaussianapproximation. For the inflationary model that we have analyzed, this modification is achieved by a minor (below the‰ level) tuning of one parameter (the non-minimal coupling to gravity).9Finally, let us stress that these results do not necessarily imply that non-gaussianities are irrelevant for the phe-nomenology of PBHs. A notable example of the contrary are the primordial gravitational wave induced by scalar cur-vature perturbations. This signal is an important route to test PBHs with future gravitational-wave interferometers. Theenergy density of the gravitational waves scales quadratically with the amplitude of the power spectrum of curvatureperturbations, and therefore there exists a correlation between the PBH abundance and the gravitational wave signal.One shall be interested in scenarios where PBHs account for all the dark matter, or such that their abundance is inagreement with observational constraints in the relevant mass range. The gravitational signal would then be a predic-tion, whose theoretical uncertainty is tied to the ones involved in the calculations of the PBHs abundance. The effectsdiscussed here are therefore relevant in this context. The interplay between gravitational signal and PBH abundancewill be studied in a future work. Acknowledgments
The research of A.U. is supported in part by the MIUR under contract 2017 FMJFMW (“New Avenues in Strong Dy-namics,” PRIN 2017) and by the INFN grant “SESAMO – SinergiE di SApore e Materia Oscura.” M.T. acknowledges sup-port from the INFN grant “LINDARK,” the research grant “The Dark Universe: A Synergic Multimessenger ApproachNo. 2017X7X85” funded by MIUR, and the project “Theoretical Astroparticle Physics (TAsP)” funded by the INFN.
Appendix A1. Computation of the three-point correlator (cid:104)R k R k R k (cid:105) In this appendix we sketch the derivation of eqs. (26, 28, 29). Even though the formalism used is quite standard,we shall explain in detail our analysis in order to better highlight the difference between our result and the existingliterature.The obtain the power spectrum of the comoving curvature perturbation, computed in section II, one has to solve theequation of motion obtained from the second-order action for the comoving curvature perturbation.This is not enough for the computation of non-gaussianities (equivalently, for the computation of higher-order cor-relators) since the latter can be thought as an imprint left by interactions. More in detail, we can split the full Hamilto-nian in three parts H full = H cl + H + H int , where H cl is the background Hamiltonian constructed with the classicalbackground field, H is the Hamiltonian constructed with scalar field perturbations up to quadratic order and H int isconstructed with perturbations at order equal or higher than three. The formalism to compute correlation functionsin cosmological settings is called the in-in formalism (introduced in ref. [52, 53] for condensed matter systems and theapplied to cosmology in refs. [4, 54–56]). The three-point correlator can be obtained by applying the master formula (cid:104) ˆ R ( τ, (cid:126)k ) ˆ R ( τ, (cid:126)k ) ˆ R ( τ, (cid:126)k ) (cid:105) = 2 Im (cid:34)(cid:90) τ −∞ (1 − i(cid:15) ) dτ (cid:48) a ( τ (cid:48) ) (cid:104) ˆ R I ( τ, (cid:126)k ) ˆ R I ( τ, (cid:126)k ) ˆ R I ( τ, (cid:126)k ) ˆ H int ( τ (cid:48) ) (cid:105) (cid:35) , (A3)where the subscript I denotes interaction picture operators, and the interaction Hamiltonian ˆ H int is also constructedwith the interaction picture operators. The i(cid:15) prescription (where (cid:15) , not to be confused with the Hubble parameter, is asmall real positive number) is introduced in order to regularize the time integral in the early time limit. This prescriptionassures that the integrand vanishes at early times, and hence the integral remains finite. At this stage, it is important tokeep in mind that the vacuum expectation value on the left-hand side is evaluated on the interaction-picture vacuumwhile the the right-hand side is evaluated on the vacuum of the free theory. The interaction Hamiltonian ˆ H int evolvesthe vacuum of the free theory to the interaction vacuum at the time τ at which we evaluate the three-point function.Eq. (A3) makes clear that non-gaussianities, encoded in three- and higher-point correlators, are a manifestation ofinteractions. Notice also that eq. (A3) is valid at the first order in ˆ H int . We remark that our goal is to compute, using For the sake of clarity, and to connect with section III, let us be clear with our notation. In the following we use ˆ R I ( τ, (cid:126)x ) = (cid:90) d (cid:126)k (2 π ) ˆ R I ( τ,(cid:126)k ) e i(cid:126)x · (cid:126)k , ˆ R I ( τ,(cid:126)k ) = R k ( τ ) a (cid:126)k ( τ (cid:63) ) + R ∗ k ( τ ) a †− (cid:126)k ( τ (cid:63) ) . (A1)The subscript I is used here to distinguish between the left- and right-hand side of eq. (A3). Furthermore, we have the contraction (cid:104) | ˆ R I ( τ ,(cid:126)k ) ˆ R I ( τ ,(cid:126)k ) | (cid:105) = (2 π ) δ (3) ( (cid:126)k + (cid:126)k ) R k ( τ ) R ∗ k ( τ ) . (A2) τ close to the end of inflation when we can safely take the super-horizon curvature perturbations as constants.The cubic action, in terms of the cosmic time, is given by S = (cid:82) dtd (cid:126)x L , with [4] L = a (cid:15) R ˙ R + a(cid:15) R ( ∂ R ) − a(cid:15) ˙ R ( ∂ R )( ∂χ ) + a (cid:15) ˙ (cid:15) R ˙ R + (cid:15) a ( ∂ R )( ∂χ ) ∂ χ + (cid:15) a ( ∂ R )( ∂χ ) + 2 f ( R ) δ S δ R , (A4) f ( R ) = (cid:15) R + 1 H R ˙ R + terms quadratic in R with spatial derivatives acting on R , (A5)where (cid:15) = ˙ (cid:15)/H(cid:15) = 2 (cid:15) − η , ∂ χ = a (cid:15) ˙ R and δ S /δ R the variation of the quadratic action with respect to theperturbation. Notice that the cubic action S is exact for arbitrary (cid:15) and η . As shown in ref. [4], the last term in eq. (A4)can be absorbed by the field redefinition R → R n + f ( R n ) so that, after this redefinition, the interaction Hamiltonianin conformal time τ is H int ( τ ) = − (cid:90) d (cid:126)x (cid:20) a(cid:15) R n ( R (cid:48) n ) + a(cid:15) R n ( ∂ R n ) − (cid:15) R (cid:48) n ( ∂ R n )( ∂χ n ) + a(cid:15)(cid:15) (cid:48) R n R (cid:48) n + (cid:15) a ( ∂ R n )( ∂χ n )( ∂ χ n ) + (cid:15) a ( ∂ R n )( ∂χ n ) (cid:21) , (A6)with ∂ χ n = a (cid:15) ˙ R n and (cid:48) denotes the derivative with respect to τ . This is not the only effect that is introduced by thethe field redefinition since it also gives some extra terms in the three-point correlation function of R . In Fourier space,we have (cid:104) ˆ R ( τ, (cid:126)k ) ˆ R ( τ, (cid:126)k ) ˆ R ( τ, (cid:126)k ) (cid:105) (A7) = (cid:104) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) (cid:105) + (cid:15) (cid:104) (cid:104) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) (cid:105) + two perms (cid:105) + higher order terms , where one should keep in mind that in eq. (A7) ˆ R n ( τ, (cid:126)k ) is not just the square of ˆ R n ( τ, (cid:126)k ) but it originates from theFourier transform of the square of ˆ R n ( τ, (cid:126)x ) (since the field redefinition is defined at the level of the action in positionspace) meaning that, by means of the convolution theorem, we actually have ˆ R n ( τ, (cid:126)k ) = (cid:90) d (cid:126)q (2 π ) ˆ R n ( τ, (cid:126)q ) ˆ R n ( τ, (cid:126)k − (cid:126)q ) , (A8)and similarely for the other two permutations. In eq. (A7) we neglected terms with spatial derivatives in eq. (A5) sincewe are interested in super-horizon modes and the contribution from the second term in eq. (A5) since we are interestedin the computation of the three-point function at late time. All in all, the strategy of the computation is the following. i) We compute the three-point function (cid:104) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) (cid:105) , that is the first term on the right-hand sideof eq. (A7), by means of eq. (A3) using the interaction Hamiltonian given in eq. (A6); ii) we shift back from the bispectrum of R n to that of R by including the leading correction given by the secondterm on the right-hand side of eq. (A7).For notation simplicity, we drop from now on the subscript n . We start from point i) . In the interaction Hamiltonianmost of the terms are (cid:15) -suppressed. The only exception which is relevant for us is the term that contains (cid:15) (cid:48) . In ourmodel, we have dηdτ = η II δ ( τ − τ in ) + ( − η II + η III ) δ ( τ − τ end ) , (A9)meaning that in the integral in eq. (A3) we have two contributions corresponding to the transition at τ = τ in and τ = τ end . Let us compute explicitly only the first contribution. The integral on the right-hand side of eq. (A3) gives(after promoting classical fields to operators) Im η II a ( τ in ) (cid:15) I (cid:90) d (cid:126)x (cid:104) ˆ R I ( τ, (cid:126)k ) ˆ R I ( τ, (cid:126)k ) ˆ R I ( τ, (cid:126)k ) ˆ R I ( τ in , (cid:126)x ) ˆ R I ( τ in , (cid:126)x ) (cid:105) = 2 Im η II a ( τ in ) (cid:15) I × (A10) (cid:90) d (cid:126)x (cid:90) d (cid:126)q (2 π ) d (cid:126)q (2 π ) d (cid:126)q (2 π ) (cid:104) ˆ R I ( τ, (cid:126)k ) ˆ R I ( τ, (cid:126)k ) ˆ R I ( τ, (cid:126)k ) ˆ R I ( τ in , (cid:126)q ) ˆ R I ( τ in , (cid:126)q − (cid:126)q ) ˆ R (cid:48) I ( τ in , (cid:126)q ) (cid:105) e i(cid:126)x · ( (cid:126)q + (cid:126)q ) , | (cid:105) of the free theory can be computed by means of the Wick theorem. We only have three relevantcontractions, each of them that counts twice since, for instance, the two contractions (cid:90) d (cid:126)x (cid:90) d (cid:126)q (2 π ) d (cid:126)q (2 π ) d (cid:126)q (2 π ) (cid:104) ˆ R I ( τ, (cid:126)k ) ˆ R I ( τ, (cid:126)k ) ˆ R I ( τ, (cid:126)k ) ˆ R I ( τ in , (cid:126)q ) ˆ R I ( τ in , (cid:126)q − (cid:126)q ) ˆ R (cid:48) I ( τ in , (cid:126)q ) (cid:105) e i(cid:126)x · ( (cid:126)q + (cid:126)q ) = (A11) (cid:90) d (cid:126)xe − i(cid:126)x · ( (cid:126)k + (cid:126)k + (cid:126)k ) (cid:124) (cid:123)(cid:122) (cid:125) =(2 π ) δ (3) ( (cid:126)k + (cid:126)k + (cid:126)k ) R k ( τ ) R k ( τ ) R k ( τ ) R ∗ k ( τ in ) R ∗ k ( τ in ) R (cid:48) ∗ k ( τ in ) , (A12)and (cid:90) d (cid:126)x (cid:90) d (cid:126)q (2 π ) d (cid:126)q (2 π ) d (cid:126)q (2 π ) ˆ R I ( τ, (cid:126)k ) ˆ R I ( τ, (cid:126)k ) ˆ R I ( τ, (cid:126)k ) ˆ R I ( τ in , (cid:126)q ) ˆ R I ( τ in , (cid:126)q − (cid:126)q ) ˆ R (cid:48) I ( τ in , (cid:126)q ) (cid:105) e i(cid:126)x · ( (cid:126)q + (cid:126)q ) , (A13)give the same result. If we factor out (2 π ) δ (3) ( (cid:126)k + (cid:126)k + (cid:126)k ) , the contribution of eq. (A10) to the bispectrum becomes B (in) R ( k , k , k ) = 4 Im (cid:26) η II a ( τ in ) (cid:15) I R k ( τ ) R k ( τ ) R k ( τ ) × (A14) (cid:2) R ∗ k ( τ in ) R ∗ k ( τ in ) R (cid:48) ∗ k ( τ in ) + R ∗ k ( τ in ) R ∗ k ( τ in ) R (cid:48) ∗ k ( τ in ) + R ∗ k ( τ in ) R ∗ k ( τ in ) R (cid:48) ∗ k ( τ in ) (cid:3) (cid:27) , taken at late time τ → − (so that the term R k ( τ ) R k ( τ ) R k ( τ ) has to be computed in region III). Using eq. (15) andthe expansion of the Hankel functions for small argument, we find eq. (26) used in the text.The computation of the contribution at τ = τ end follows the same steps but now we pick up the second delta functionin eq. (A9). We find the result B (end) R ( k , k , k ) = 4 Im (cid:26) ( − η II + η III ) a ( τ end ) (cid:15) II ( τ end ) R k ( τ ) R k ( τ ) R k ( τ ) × (A15) (cid:2) R ∗ k ( τ end ) R ∗ k ( τ end ) R (cid:48) ∗ k ( τ end ) + R ∗ k ( τ end ) R ∗ k ( τ end ) R (cid:48) ∗ k ( τ end ) + R ∗ k ( τ end ) R ∗ k ( τ end ) R (cid:48) ∗ k ( τ end ) (cid:3) (cid:27) . Finally, we move to consider the contribution ii) that comes from the field redefinition in eq. (A7). If we neglect (cid:15) , wehave (cid:15) = − η III since this term must be evaluated towards the end of inflation in region III. The correction in eq. (A7)reads − η III (cid:90) d (cid:126)q (2 π ) (cid:104) (cid:104) ˆ R I ( τ, (cid:126)q ) ˆ R I ( τ, (cid:126)k − (cid:126)q ) ˆ R I ( τ, (cid:126)k ) ˆ R I ( τ, (cid:126)k ) (cid:105) + two permutations (cid:105) (A16)Consider the first term in the square brackets. If we apply the Wick theorem, we have again that the contraction − η III (cid:90) d (cid:126)q (2 π ) (cid:104) ˆ R I ( τ, (cid:126)q ) ˆ R I ( τ, (cid:126)k − (cid:126)q ) ˆ R I ( τ, (cid:126)k ) ˆ R I ( τ, (cid:126)k ) (cid:105) = − η III π ) δ (3) ( (cid:126)k + (cid:126)k + (cid:126)k ) |R k ( τ ) | |R k ( τ ) | (A17)counts twice since it gives the same result of − η III (cid:90) d (cid:126)q (2 π ) (cid:104) ˆ R I ( τ, (cid:126)q ) ˆ R I ( τ, (cid:126)k − (cid:126)q ) ˆ R I ( τ, (cid:126)k ) ˆ R I ( τ, (cid:126)k ) (cid:105) . (A18) Notice that in principle there are 15 ways of contracting pairs of operators (divided in three classes: pairings of external fields with fields at theinteraction vertex, pairings of external fields among themselves and vertex fields among themselves). However, we are only interested in theconnected contractions, i.e. contractions where one external operator in the string ˆ R I ( τ,(cid:126)k ) ˆ R I ( τ,(cid:126)k ) ˆ R I ( τ,(cid:126)k ) is contracted with one operatorof the interaction vertex H int ; this reduces the number of contractions down to 6. Physically, the disconnected contractions which are not includedare proportional to tadpole diagrams but, since we are assuming that the background solution is stable (or, equivalently, that (cid:104) ˆ R(cid:105) = 0 ) all suchamplitudes vanish. (2 π ) δ (3) ( (cid:126)k + (cid:126)k + (cid:126)k ) fromthe definition of the bispectrum, we find B (red) R ( k , k , k ) = − η III (cid:0) |R k ( τ ) | |R k ( τ ) | + |R k ( τ ) | |R k ( τ ) | + |R k ( τ ) | |R k ( τ ) | (cid:1) = − π η III ( k k k ) (cid:2) k P R ( k ) P R ( k ) + k P R ( k ) P R ( k ) + k P R ( k ) P R ( k ) (cid:3) , (A19)which reproduces eq. (29) (and agrees, if one takes P R ( k ) = H / π (cid:15) , with the standard result obtained in single-fieldslow-roll inflation, see e.g. ref. [57]).We are now in the position to compare our result with refs. [41, 58]. Ref. [58] computes the amount of non-gaussianities in the case of non-attractor inflation models. In the specific realization considered, a canonically nor-malized scalar field rolls under a potential V ( φ ) = (cid:40) V for φ < φ c V ( φ ) for φ > φ c (A20)which is constant below some critical field value φ c and where V ( φ ) supports a second phase of slow-roll inflationas in conventional models. In this model, which assumes the Minkowski vacuum deep inside the horizon, η transitsfrom η = 3 (dubbed non-attractor phase) to η = 0 during the second slow-roll phase. Since η = 0 during the secondinflationary phase, the only contribution to the three-point function comes from the delta function at transition timefor η . Let us indicate this transition time with τ e . Ref. [58] also assumes that, immediately after the end of the non-attractor phase, the super-horizon curvature perturbations cease evolving (meaning that the comoving time τ at whicheq. (A3) has to be evaluated is taken to be the transition time τ e ). As a result, ref. [58] finds a local bispectrum (seeeq. (45)) with amplitude which, in the squeezed limit, is given by f NL = 5 / . Notice that, based on this result, ref. [40]evaluated the amount of non-gaussianities expected in the case of PBH formation in single-field inflationary modelsfeaturing a USR phase. However, the model of ref. [58], designed for the case of non-attractor inflation, cannot beapplied to the case relevant for PBH formation.Ref. [41] improves the analysis of ref. [58] by assuming that the non-attractor inflation model consists of at leastthree different phases: the non-attractor phase (with η = 3 ), the transition phase (which is modeled in two waysdubbed smooth and sharp), and the final slow-roll phase (where η eventually relaxes to η = 0 ). We reproduce, for � � ��� - � - �� η � � � � � � ��� - � - �� η FIG. 8:
Time evolution of η for the two models discussed in ref. [41]. The left panel shows a typical smooth transition from the non-attractor to the slow-roll phase. We indicate with N e the e -fold time at which the non-attractor phase ends. The right panel showsthree typical sharp transitions from the non-attractor to the slow-roll phase (the three lines have different inflaton velocity at the endof the non-attractor phase, see ref. [41] for details; furthermore we chose, for the sake of clarity, three different values of N e ). ease of reading, in fig. 8 the qualitative behavior of η as function of time (with N e the e -fold time at which the non-attractor phase ends) for the two cases discussed in ref. [41] (smooth and sharp transition). In the case of a smoothtransition, ref. [41] finds that the non-gaussianity computed in ref. [58] by means of the instantaneous approximationis almost cancelled. In the case of a sharp transition, ref. [41] finds that the amplitude of local non-Gaussianity is mainlydetermined by the value of the inflaton velocity at the end of the non-attractor phase.3Without going into detail, we remark here that none of these models, designed for the case of non-attractor infla-tion, can be considered as an optimal proxy for the case of PBH production. This is evident if we compare the timeevolution of η shown in fig. 2 with those in fig. 8. In realistic single-field inflationary models of PBH production whichare compatible with observations, the field starts from slow-roll conditions ( η ≈ ), goes through a USR phase ( η > )and approaches the end of inflation via a third phase during which η < . This dynamics leaves a specific imprints onthe three-point correlator which departs, both qualitatively and quantitatively, from the previous works, as discussedin this appendix.
2. Computation of the three-point correlator (cid:104)R k R k R k (cid:105) : further remarks Let us add few additional comments to the computation performed in appendix A 1. a) We use the cubic action in eq. (A4) as opposed to the expression [4, 57] S = (cid:90) dtd (cid:126)x (cid:26) − a(cid:15) R ( ∂ R ) − a (cid:15)H ˙ R + 3 a (cid:15) R ˙ R + 12 a (cid:18) R − ˙ R H (cid:19) (cid:2) ( ∂ i ∂ j ψ )( ∂ i ∂ j ψ ) − ( ∂ ψ )( ∂ ψ ) (cid:3) − a ( ∂ i ψ )( ∂ i R )( ∂ ψ ) (cid:27) , (A21)with ψ = −R /H + a (cid:15)∂ − R because all terms in eq. (A4) are explicitly proportional to the Hubble parame-ters and, therefore, more easy to handle. Furthermore, in eq. (A4) the cubic vertex is simpler since the last termcan be eliminated via a field redefinition [4]. The action obtained by means of the Lagrangian density in eq. (A4)can be derived from the one in eq. (A21) after integrating by parts, discarding total derivatives, and using thefield equation that follows from the quadratic action [4]. Notice that eq. (A21) was used in ref. [8] to computenon-gaussianities related to PBH production in the context of the chaotic new inflation model. Ref. [8] finds �� �� �� �� ����� � ϕ ϕ � ���� �� �� �� �� �� ���� - � �� - � �� - � ����� � �� � � η η = � FIG. 9:
Inflaton dynamics for the Coleman-Weinberg model studied in ref. [8]. In the left panel we show the time evolution of theinflaton field and in the right panel the time evolution of the Hubble parameter η . that non-gaussianities give small corrections to the PBH abundance. However, the model presented in ref. [8]differs—both qualitatively and quantitatively—from the class of models studied in this paper. The inflaton fol-lows a Coleman-Weinberg potential (see inset plot in the left panel of fig. 9) starting from large field values androlling down towards the origin. Initial conditions and model parameters are tuned in such a way that the infla-ton crosses the positive potential minimum and the barrier at the origin, falls in the negative potential minimum,bounces back, crosses a second time the origin, bounces back again (this time on the positive side of the poten-tial), climbs again the central barrier and it has just enough inertia to spend N (cid:38) e -folds moving slowly in theneighborhood of the origin after the oscillation. For ease of reading, we reproduce the classical dynamics of theinflaton in the left panel of fig. 9 where the oscillation and the subsequent second phase of inflation is evident.The latter ends when the inflaton falls again in the negative potential minimum, where it finally settles down. Inthis model, PBHs are produced during the second phase of inflation when the inflaton moves slowly close to theorigin. In the right panel of fig. 9, we computed the time evolution of the η for the same model parameters andinitial conditions used in the left panel. The solid (dashed) line indicates positive (negative) values. The time4profile of η is considerably more complicated than the one expected in our class of models mostly because ofthe presence of the oscillating phase during which η changes twice from large positive to large negative valuesbefore taking for ∆ N (cid:39) e -folds the constant value η (cid:39) close to the origin (where the potential is flat). Thisis a clear example of a model that differs from the class of models that we analyze in this paper. However, we donot proceed further in the explicit analysis of this setup since one gets n s (cid:39) . and r (cid:39) . for, respectively,the spectral index and the tensor-to-scalar ratio. These values are now grossly excluded by the most recent CMBmeasurements. b) Consider the field redefinition
R → R n + f ( R n ) and the quadratic action for RS [ R ] = (cid:90) dτ d (cid:126)x a (cid:15) (cid:2) ( R (cid:48) ) − ( ∂ R ) (cid:3) . (A22)Since f ( R n ) is a quadratic function in R n , the shift in the quadratic action generates cubic and quartic termsfor R n but not additional quadratic ones; schematically, we have S [ R ] → S [ R n ] + ˜ S [ R n ] + O ( R n ) , and thequadratic action for R n has, therefore, exactly the same form as for R . This is the reason why even if we formallywork with R n in eqs. (A6, A7) we can still use the free-field solutions. For the very same reason, the shift does notchange L in eq. (A4) but it just relabels R → R n ; however, the utility of the field redefinition introduced in [4] isthat what cancels the last term in eq. (A4) is the cubic correction ˜ S [ R n ] introduced before and originating fromthe quadratic action, as one can explicitly check. c) Consider now some of the interaction terms in eq. (A6) that we neglected. Let us focus on the interaction term H ( c )int ( τ ) ≡ − (cid:90) d (cid:126)xa(cid:15) R ( R (cid:48) ) , (A23)which generates a contribution to the bispectrum of the form B ( c ) R = − Im R k ( τ ) R k ( τ ) R k ( τ ) (cid:90) τ −∞ (1 − i(cid:15) ) d ¯ τ a (¯ τ ) (cid:15) (¯ τ ) (cid:2) R ∗ k (¯ τ ) R (cid:48) ∗ k (¯ τ ) R (cid:48) ∗ k (¯ τ ) + two perms (cid:3) . (A24)At late time the integrand eventually goes to zero because the perturbations settle to constant values. However,the integral is computed along the whole inflationary history, and, in particular, it crosses the USR phase duringwhich the perturbations grow exponentially. It is simple to see that the contribution to the bispectrum fromeq. (A24) scales according to B ( c ) R = (cid:15) I (cid:18) H (cid:15) I k (cid:19) F c ( x , x , x ) , (A25)where F c is a shape function which depends on the integral in eq. (A24). As expected, eq. (A25) pays a sup-pression factor (cid:15) I (cid:28) if compared to the contributions we included in our computation. For this reason, weneglected this contribution. Of course, in more realistic models in which it is well possible to reach (cid:15) I = O (1) during the pre-USR dynamics, the same argument can be invalidated. To proceed further, one has to evaluatethe time integral and compute the shape function explicitly. In our model, we can split the integral in three parts.In region I, the integral is analytical. The contribution at τ = −∞ (1 − i(cid:15) ) vanishes thanks to the (cid:15) prescription(physically, this means that interaction are switched off in the Bunch-Davis vacuum), and only the contribution at τ = τ in survives. In region II and region III, the integral can be computed numerically. In region III perturbationswith k (cid:46) k peak settle exponentially fast to their final constant value, and the integral gives a small contributionsince the integrand quickly vanishes. On the contrary, as discussed in ref. [27], perturbations with k (cid:38) k peak evolves in region III until their horizon crossing time N k > N end (after which they become constant) and theycontribute to the integral in eq. (A24). We checked numerically that the contribution in eq. (A25) to the bispec-trum is largely subdominant.Consider now the interaction term H ( c )int ( τ ) ≡ (cid:90) d (cid:126)x(cid:15) R (cid:48) ( ∂ i R )( ∂ i χ ) , χ = (cid:15)a ∂ − ˙ R , (A26)with ∂ − the inverse laplacian. The corresponding contribution to the bispectrum reads B ( c ) R = 4 Im R k ( τ ) R k ( τ ) R k ( τ ) (cid:90) τ −∞ (1 − i(cid:15) ) d ¯ τ a (¯ τ ) (cid:15) (¯ τ ) (cid:34) (cid:126)k · (cid:126)k k R (cid:48) ∗ k (¯ τ ) R ∗ k (¯ τ ) R (cid:48) ∗ k (¯ τ ) + six perms (cid:35) , (A27)5where the scalar products can be expressed in terms of the magnitudes k i using the conservation enforced by thetriangle identity, for instance (cid:126)k + (cid:126)k = − (cid:126)k = ⇒ (cid:126)k · (cid:126)k = ( k − k − k ) / . As before, eq. (A27) features thescaling B ( c ) R = (cid:15) I (cid:18) H (cid:15) I k (cid:19) F c ( x , x , x ) , (A28)with the extra suppression given by (cid:15) I (cid:28) . The explicit computation of F c ( x , x , x ) is not conceptuallydifferent from that of F c ( x , x , x ) and leads to comparable results.The few examples considered here confirm that the bispectrum computed from the time variation of η capturesthe leading contribution coming from the interaction Hamiltonian. We also showed that organizing cubic inter-actions according to eq. (A6) is indeed useful since it allows to exploit a power counting in terms of the Hubbleparameter (cid:15) (assuming the latter remains small during the inflationary dynamics). d) A comment on eq. (A7) is needed. As anticipated, the three-point function is meant to be evaluated in the vacuumstate of the full interacting theory. More formally, one should write (cid:104) ˆ R ( τ, (cid:126)k ) ˆ R ( τ, (cid:126)k ) ˆ R ( τ, (cid:126)k ) (cid:105) = (cid:104) Ω( τ ) | ˆ R I ( τ, (cid:126)k ) ˆ R I ( τ, (cid:126)k ) ˆ R I ( τ, (cid:126)k ) | Ω( τ ) (cid:105) , (A29)with the vacuum of the interacting theory given by | Ω( τ ) (cid:105) = T e − i (cid:82) τ −∞ (1 − i(cid:15) ) d ¯ τa (¯ τ ) H int (¯ τ ) | (cid:105) where T is the time-ordering operator. At initial time, when interactions are turned off (thanks to the i(cid:15) prescription) the initial vac-uum reduces to the Bunch-Davis vacuum | (cid:105) . Notice that the three-point function (cid:104) ˆ R ( τ, (cid:126)k ) ˆ R ( τ, (cid:126)k ) ˆ R ( τ, (cid:126)k ) (cid:105) would vanish if evaluated in | (cid:105) (since it contains an odd numbers of operators) so that its non-zero value is agenuine imprint left by interactions.The same comment applies to all terms on the right-hand side of eq. (A7) which should be evaluated in the vac-uum state of the full interacting theory. This is indeed the way in which we computed the first contribution onthe right-hand side, (cid:104) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) (cid:105) . Consider now the other term (cid:15) (cid:104) (cid:104) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) (cid:105) + two permutations (cid:105) . (A30)We computed this term in the vacuum of the free theory. The difference is that in this case the operator is quarticin ˆ R n , and we have a non-zero leading contribution already in the vacuum | (cid:105) (the first term in the Dyson seriesof | Ω( τ ) (cid:105) ). This is the contribution that has been computed in eq. (A19). Further corrections to this contributioncan be discussed (dubbed “higher order terms” in eq. (A7)). One possibility is given by the term (cid:15) (cid:104) (cid:104) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) (cid:105) + permutations (cid:105) , (A31)which is also non-zero evaluated on the free-field vacuum. This term generates a contribution to the bispectrumwhich has the form B ( d ) R = − η (cid:90) d Ω q (4 π ) (cid:90) dqq (2 π ) ( k + q ) ( k − q ) P R ( k + q ) P R ( q ) P R ( k − q ) + permutations . (A32)If we approximate P R ( q ) as a spiky power spectrum so that the integration over q picks up only the contributionat k peak , a quick estimate of the previous term compared with the one in eq. (29) would be B ( d ) R B (red) R = O ( η P R ( k peak )) , with P R ( k peak ) ∼ − . (A33)The point to make here is that terms like the one in eq. (A32) are usually tossed away in slow-roll inflation modelssince highly suppressed by P R ∼ − . In the presence of USR, one can not take lightly the same conclusionsince the power spectrum can be much larger than the value suggested by CMB measurements. However, thesuppression in eq. (A33) might still allow to neglect these contributions. We can compute numerically the ratioin eq. (A33) and check our estimate. It is worth doing this check explicitly since this is the largest correction to ourapproximation. We show our results, computed in the context of the numerical model, in fig. 10. The suppressionestimated in eq. (A33) is confirmed.6 �� �� �� �� �� �� �� �� �� - � �� - � �� - � [ ��� - � ] � ℛ ( � ) / � ℛ ( ℯ ) � = � = � ≡ ��������� ����� - ��� - ��� - ��� - ��� � ������ ��� ����� - � �� - � �� - � α ∈ [-( � - β ) � ( � - β )] � ℛ ( � ) / � ℛ ( ℯ ) β = ����������� � = �� �� ��� - � FIG. 10:
Ratio in eq. (A33) computed numerically. We consider the numerical model, and we fix η III = η = − . . In the left panel,we show the ratio in the equilateral configuration. In the right panel, we use the parametrization in eq. (21) for fixed K and threedifferent choices for β . e) Consider the effect of the field redefinition on the four-point correlator (cid:104) ˆ R ( τ, (cid:126)k ) ˆ R ( τ, (cid:126)k ) ˆ R ( τ, (cid:126)k ) ˆ R ( τ, (cid:126)k ) (cid:105) . InFourier space, the field redefinition generates, at order O ( (cid:15) ) , the following 6 terms (cid:15) (cid:20) (cid:104) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) (cid:105) + (cid:104) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) (cid:105) + (cid:104) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) (cid:105) + (cid:104) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) (cid:105) + (cid:104) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) (cid:105) + (cid:104) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) ˆ R n ( τ, (cid:126)k ) (cid:105) (cid:21) . (A34)Each one of them can be computed in the vacuum of the free theory. Consider for simplicity only the first one,which we write in the form (cid:15) (cid:90) d (cid:126)q (2 π ) d (cid:126)q (2 π ) (cid:104) ˆ R I ( τ, (cid:126)k ) ˆ R I ( τ, (cid:126)k ) ˆ R I ( τ, (cid:126)q ) ˆ R I ( τ, (cid:126)k − (cid:126)q ) ˆ R I ( τ, (cid:126)q ) ˆ R I ( τ, (cid:126)k − (cid:126)q ) (cid:105) . (A35)It is easy to see that there are 8 possible contractions, and one obtains (cid:15) (cid:110) ∆( k )∆( k ) (cid:104) ∆( | (cid:126)k + (cid:126)k | ) + ∆( | (cid:126)k + (cid:126)k | ) + ∆( | (cid:126)k + (cid:126)k | ) + ∆( | (cid:126)k + (cid:126)k | ) (cid:105)(cid:111) , (A36)where the four terms in the curly brackets counted twice. Putting all together, eq. (A34) reads (cid:15) (cid:104) ∆( k )∆( k )∆( | (cid:126)k + (cid:126)k | ) + 23 permutations (cid:105) , (A37)where the
4! = 24 terms in the square brackets counts the permutations of the indices . This is the samestructure that we introduced for the reduced trispectrum in eq. (39). If we identify (cid:15) / f NL / , the coefficientin front of eq. (A37) becomes f / .7 [1] Ya. B. Zel’dovich and I. D. Novikov, “ The Hypothesis of Cores Retarded during Expansion and the Hot Cosmological Model, ” As-tronomicheskii Zhurnal, (1966), 758.[2] S. Hawking, “ Gravitationally collapsed objects of very low mass, ” Mon. Not. Roy. Astron. Soc. (1971), 75[3] B. J. Carr and S. Hawking, “
Black holes in the early Universe, ” Mon. Not. Roy. Astron. Soc. (1974), 399-415[4] J. M. Maldacena, “
Non-Gaussian features of primordial fluctuations in single field inflationary models, ” JHEP (2003), 013[arXiv:astro-ph/0210603 [astro-ph]].[5] S. M. Leach and A. R. Liddle, “ Inflationary perturbations near horizon crossing, ” Phys. Rev. D (2001), 043508 [arXiv:astro-ph/0010082].[6] S. M. Leach, M. Sasaki, D. Wands and A. R. Liddle, “ Enhancement of superhorizon scale inflationary curvature perturbations, ”Phys. Rev. D (2001), 023512 [arXiv:astro-ph/0101406].[7] P. Ivanov, P. Naselsky and I. Novikov, “ Inflation and primordial black holes as dark matter, ” Phys. Rev. D (1994) 7173.[8] R. Saito, J. Yokoyama and R. Nagata, “ Single-field inflation, anomalous enhancement of superhorizon fluctuations, and non-Gaussianity in primordial black hole formation, ” JCAP (2008), 024 [arXiv:0804.3470 [astro-ph]].[9] A. M. Green, A. R. Liddle, K. A. Malik and M. Sasaki, “ A New calculation of the mass fraction of primordial black holes, ” Phys.Rev. D (2004), 041502 [arXiv:astro-ph/0403181 [astro-ph]].[10] C. Germani and R. K. Sheth, “ Nonlinear statistics of primordial black holes from Gaussian curvature perturbations, ” Phys. Rev.D (2020) no.6, 063520 [arXiv:1912.07072 [astro-ph.CO]].[11] V. De Luca, G. Franciolini, A. Kehagias, M. Peloso, A. Riotto and C. Ünal, “
The Ineludible non-Gaussianity of the Primordial BlackHole Abundance, ” JCAP (2019), 048 [arXiv:1904.00970 [astro-ph.CO]].[12] C. M. Yoo, T. Harada, J. Garriga and K. Kohri, “ Primordial black hole abundance from random Gaussian curvature perturbationsand a local density threshold, ” PTEP (2018) no.12, 123E01 [arXiv:1805.03946 [astro-ph.CO]].[13] M. Kawasaki and H. Nakatsuka, “
Effect of nonlinearity between density and curvature perturbations on the primordial black holeformation, ” Phys. Rev. D (2019) no.12, 123501 [arXiv:1903.02994 [astro-ph.CO]].[14] S. Young, I. Musco and C. T. Byrnes, “ Primordial black hole formation and abundance: contribution from the non-linear relationbetween the density and curvature perturbation, ” JCAP (2019), 012 [arXiv:1904.00984 [astro-ph.CO]].[15] A. Kehagias, I. Musco and A. Riotto, “ Non-Gaussian Formation of Primordial Black Holes: Effects on the Threshold, ” JCAP (2019), 029 [arXiv:1906.07135 [astro-ph.CO]].[16] C. M. Yoo, J. O. Gong and S. Yokoyama, “ Abundance of primordial black holes with local non-Gaussianity in peak theory, ” JCAP (2019), 033 [arXiv:1906.06790 [astro-ph.CO]].[17] C. M. Yoo, T. Harada, S. Hirano and K. Kohri, “ Abundance of Primordial Black Holes in Peak Theory for an Arbitrary PowerSpectrum, ” [arXiv:2008.02425 [astro-ph.CO]].[18] I. Musco, V. De Luca, G. Franciolini and A. Riotto, [arXiv:2011.03014 [astro-ph.CO]].[19] I. Musco, J. C. Miller and L. Rezzolla, “
Computations of primordial black hole formation, ” Class. Quant. Grav. (2005), 1405-1424 [arXiv:gr-qc/0412063 [gr-qc]].[20] S. Young, C. T. Byrnes and M. Sasaki, “ Calculating the mass fraction of primordial black holes, ” JCAP (2014), 045[arXiv:1405.7023 [gr-qc]].[21] J. M. Bardeen, J. Bond, N. Kaiser and A. Szalay, “ The Statistics of Peaks of Gaussian Random Fields, ” Astrophys. J. (1986),15-61[22] C. Germani and I. Musco, “
Abundance of Primordial Black Holes Depends on the Shape of the Inflationary Power Spectrum, ”Phys. Rev. Lett. (2019) no.14, 141302 [arXiv:1805.04087 [astro-ph.CO]].[23] I. Musco, “
Threshold for primordial black holes: Dependence on the shape of the cosmological perturbations, ” Phys. Rev. D (2019) no.12, 123524 [arXiv:1809.02127 [gr-qc]].[24] D. W. Neilsen and M. W. Choptuik, “
Critical phenomena in perfect fluids, ” Class. Quant. Grav. (2000), 761-782 [arXiv:gr-qc/9812053 [gr-qc]].[25] A. Kalaja, N. Bellomo, N. Bartolo, D. Bertacca, S. Matarrese, I. Musco, A. Raccanelli and L. Verde, “ From Primordial Black HolesAbundance to Primordial Curvature Power Spectrum (and back), ” JCAP (2019), 031 [arXiv:1908.03596 [astro-ph.CO]].[26] F. Riccardi, M. Taoso, A. Urbano, “ Solving peak theory in the presence of local non-gaussianities, ” to appear.[27] G. Ballesteros, J. Rey, M. Taoso and A. Urbano, “
Primordial black holes as dark matter and gravitational waves from single-fieldpolynomial inflation, ” JCAP , 025 (2020) arXiv:2001.08220 [astro-ph.CO].[28] G. Ballesteros, J. Rey, M. Taoso and A. Urbano, ‘ ‘Stochastic inflationary dynamics beyond slow-roll and consequences for primor-dial black hole formation, ” JCAP , 043 (2020) [arXiv:2006.14597 [astro-ph.CO]].[29] B. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, “ New cosmological constraints on primordial black holes, ” Phys. Rev. D (2010),104019 [arXiv:0912.5297 [astro-ph.CO]].[30] R. Laha, “ Primordial Black Holes as a Dark Matter Candidate Are Severely Constrained by the Galactic Center 511 keV γ -RayLine, ” Phys. Rev. Lett. (2019) no.25, 251101 [arXiv:1906.09994 [astro-ph.HE]].[31] H. Niikura, M. Takada, N. Yasuda, R. H. Lupton, T. Sumi, S. More, T. Kurita, S. Sugiyama, A. More, M. Oguri and M. Chiba,“ Microlensing constraints on primordial black holes with Subaru/HSC Andromeda observations, ” Nature Astron. (2019) no.6,524-534 [arXiv:1701.02151 [astro-ph.CO]].[32] A. Katz, J. Kopp, S. Sibiryakov and W. Xue, “ Femtolensing by Dark Matter Revisited, ” JCAP (2018), 005 [arXiv:1807.11495[astro-ph.CO]].[33] F. Capela, M. Pshirkov and P. Tinyakov, “ Constraints on primordial black holes as dark matter candidates from capture by neutronstars, ” Phys. Rev. D (2013) no.12, 123524 [arXiv:1301.4984 [astro-ph.CO]].[34] P. Pani and A. Loeb, “ Tidal capture of a primordial black hole by a neutron star: implications for constraints on dark matter, ” JCAP (2014), 026 [arXiv:1401.3025 [astro-ph.CO]].[35] F. Capela, M. Pshirkov and P. Tinyakov, “ A comment on “Exclusion of the remaining mass window for primordial black holes ...”,arXiv:1401.3025, ” [arXiv:1402.4671 [astro-ph.CO]].[36] G. Defillon, E. Granet, P. Tinyakov and M. H. Tytgat, “
Tidal capture of primordial black holes by neutron stars, ” Phys. Rev. D (2014) no.10, 103522 [arXiv:1409.0469 [gr-qc]].[37] P. W. Graham, S. Rajendran and J. Varela, “ Dark Matter Triggers of Supernovae, ” Phys. Rev. D (2015) no.6, 063007[arXiv:1505.04444 [hep-ph]].[38] Y. Akrami et al. [Planck], “ Planck 2018 results. IX. Constraints on primordial non-Gaussianity, ” [arXiv:1905.05697 [astro-ph.CO]].[39] S. Matarrese, F. Lucchin and S. A. Bonometto, “
A path integral approach to large scale matter distribution originated by non-gaussian fluctuations, ” Astrophys. J. (1986), L21-L26[40] G. Franciolini, A. Kehagias, S. Matarrese and A. Riotto, “
Primordial Black Holes from Inflation and non-Gaussianity, ” JCAP (2018), 016 [arXiv:1801.09415 [astro-ph.CO]].[41] Y. Cai, X. Chen, M. H. Namjoo, M. Sasaki, D. Wang and Z. Wang, “ Revisiting non-Gaussianity from non-attractor inflation mod-els, ” JCAP (2018), 012 [arXiv:1712.09998 [astro-ph.CO]].[42] S. Passaglia, W. Hu and H. Motohashi, “ Primordial black holes and local non-Gaussianity in canonical inflation, ” Phys. Rev. D (2019) no.4, 043536 [arXiv:1812.08243 [astro-ph.CO]].[43] E. Komatsu, “ Hunting for Primordial Non-Gaussianity in the Cosmic Microwave Background, ” Class. Quant. Grav. (2010),124010 [arXiv:1003.6097 [astro-ph.CO]].[44] L. Boubekeur and D. Lyth, “ Detecting a small perturbation through its non-Gaussianity, ” Phys. Rev. D (2006), 021301[arXiv:astro-ph/0504046 [astro-ph]].[45] A. Stuart and K. Ord, “ Kendall’s Advanced Theory of Statistics. Volume I: Distribution Theory, ” Oxford University Press, 1987.[46] J. E. Kolassa, “
Series Approximation Methods in Statistics, ” Lecture Notes in Statistics, Springer, 2006.[47] T. Harada, C. M. Yoo, T. Nakama and Y. Koga, “
Cosmological long-wavelength solutions and primordial black hole formation, ”Phys. Rev. D (2015) no.8, 084057 [arXiv:1503.03934 [gr-qc]].[48] G. Ballesteros and M. Taoso, “ Primordial black hole dark matter from single field inflation, ” Phys. Rev. D (2018) no.2, 023501[arXiv:1709.05565 [hep-ph]].[49] I. Dalianis, A. Kehagias and G. Tringas, “ Primordial black holes from α -attractors, ” JCAP (2019), 037 [arXiv:1805.09483 [astro-ph.CO]].[50] V. Atal and C. Germani, Phys. Dark Univ. , 100275 (2019) doi:10.1016/j.dark.2019.100275 [arXiv:1811.07857 [astro-ph.CO]].[51] R. Bravo and G. A. Palma, [arXiv:2009.03369 [hep-th]].[52] J. S. Schwinger, “ Brownian motion of a quantum oscillator, ” J. Math. Phys. (1961), 407-432.[53] L. Keldysh, “ Diagram technique for non-equilibrium processes, ” Zh. Eksp. Teor. Fiz. (1964), 1515-1527.[54] R. Jordan, “ Effective Field Equations for Expectation Values, ” Phys. Rev. D (1986), 444-454.[55] E. Calzetta and B. Hu, “ Closed Time Path Functional Formalism in Curved Space-Time: Application to Cosmological Back Reac-tion Problems, ” Phys. Rev. D (1987), 495[56] S. Weinberg, “ Quantum contributions to cosmological correlations, ” Phys. Rev. D (2005), 043514 [arXiv:hep-th/0506236 [hep-th]].[57] X. Chen, M. Huang, S. Kachru and G. Shiu, “ Observational signatures and non-Gaussianities of general single field inflation, ”JCAP (2007), 002 [arXiv:hep-th/0605045 [hep-th]].[58] M. H. Namjoo, H. Firouzjahi and M. Sasaki, “ Violation of non-Gaussianity consistency relation in a single field inflationarymodel, ” EPL101