Solving peak theory in the presence of local non-gaussianities
SSolving peak theory in the presence of local non-gaussianities
Flavio Riccardi a,b , ∗ Marco Taoso c , † and Alfredo Urbano d,e ‡ a SISSA, via Bonomea 265, I-34132 Trieste, Italy b I.N.F.N. sezione di Trieste, SISSA, via Bonomea 265, I-34132 Trieste, Italy c I.N.F.N. sezione di Torino, via P. Giuria 1, I-10125 Torino, Italy d Dipartimento di Fisica, “Sapienza” Università di Roma, Piazzale Aldo Moro 5, 00185, Roma, Italy and e I.F.P.U., Institute for Fundamental Physics of the Universe, via Beirut 2, I-34014 Trieste, Italy. (Dated: February 9, 2021)We compute the probability density distribution of maxima for a scalar random field in the presence of localnon-gaussianities. The physics outcome of this analysis is the following. If we focus on maxima whosecurvature is larger than a certain threshold for gravitational collapse, our calculations illustrate how thefraction of the Universe’s mass in the form of primordial black holes (PBHs) changes in the presence of localnon-gaussianities. We find that previous literature on the subject exponentially overestimate, by manyorders of magnitude, the impact of local non-gaussianities on the PBH abundance. We explain the origin ofthis discrepancy, and conclude that, in realistic single-field inflationary models with ultra slow-roll, one canobtain the same abundance found with the gaussian approximation simply changing the peak amplitudeof the curvature power spectrum by no more than a factor of two. We comment about the relevance ofnon-gaussianities for second-order gravitational waves.
I. INTRODUCTION
The possibility that the totality of dark matter in the Universe consists of primordial black holes (PBHs) still holdsthe stage even though almost half-a-century has passed after the pioneering proposal of Hawking and Carr [1, 2]. Thisis especially true in the mass range (cid:46) M PBH [g] (cid:46) in which black holes are neither too light (otherwise theywould have evaporated in the past through Hawking radiation [3]) or too heavy (otherwise they would distort space-time in a way that contradicts present bounds from lensing experiments [4, 5]).PBHs could have formed in the very early Universe during the radiation dominated era. The key ingredient thattriggers the formation of a PBH is the presence of an over-fluctuation in the density of the Universe which, if largeenough, gravitationally collapses dragging down any matter within its horizon, that is the parcel of space around anypoint reachable at the speed of light.The theory of inflation provides an elegant mechanism that explains the origin of density perturbations in the Uni-verse. In the inflationary picture, space-time fluctuates quantum mechanically around a background that is expandingexponentially fast. After the end of inflation, these curvature fluctuations are transferred to the radiation field, creatingslightly overdense and under-dense regions. It is, therefore, fascinating to ask whether the formation of PBHs fits in theinflationary picture of structure formation. To answer this question, two (related) aspects need to be addressed. i) The inflaton dynamics should give rise to a peak in the power spectrum of curvature perturbations. This trans-lates into a large variance for density perturbations that, in turn, enhances the chance to create overdense regionsabove the threshold for gravitational collapse. ii)
The abundance of such collapsing regions should be large enough to explain the totality of dark matter.In this paper we focus on simple single-field inflationary models. It is known that in order to fulfil point i) slow-rollconditions must be violated. The simplest option is to introduce, few e -folds before the end of inflation, an approximatestationary inflection point in the inflaton potential (see refs. [11–13] for the earliest proposal in this direction). Whenthe inflaton, during its classical dynamics, crosses this region (during the so-called “ultra slow-roll phase”), curvatureperturbations, due to the presence of negative friction, get exponentially enhanced [14, 15]. Point ii) is more subtle. ∗ Electronic address: fl[email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] It is also possible to have PBH formation during matter domination [6]. This is not the only option. The formation of PBHs may have been independent of inflationary physics; PBHs may have been originated fromtopological defects formed during symmetry breaking phase transition, for instance from the collapse of string loops [7–9]. An exception, where no amplification of the power spectrum is needed, are models where PBHs form from the collapse of domain walls createdduring inflation [10]. Alternatively, a parametric amplification of curvature perturbations could be caused by resonance with oscillations in the sound speed of theirpropagation [16]. Another possibility is that, after the inflationary phase, the inflaton begins to oscillate near the minimum of the potential andfragments into oscillons which, in turn, lead to copious production of PBHs [17]. a r X i v : . [ a s t r o - ph . C O ] F e b Due of their intrinsic quantum-mechanical origin, the way in which quantum fluctuations lead to a classical patternof perturbations can be described only in a probabilistic sense. Consequently, the computation of the abundance ofcollapsing regions requires informations about the statistical distribution of density perturbations. Most of the time, forsimplicity, the gaussian approximation is assumed. However, the very same fact that slow-roll conditions are violatedas a consequence of point i) suggests that non-gaussianities may play a relevant role. Refs. [18, 19] indeed find thatduring an ultra slow-roll phase sizable non-gaussianities of local type are generated. In the rest of this paper we willdub these non-gaussianities “primordial” to distinguish them from non-gaussianities that arise from the non-linearrelation between curvature and density perturbations.What is the impact of primordial non-gaussianities on the gaussian approximation when computing the PBH abun-dance? Ref. [20] addressed this question in the context of threshold statistics. The main result of ref. [20] is that theabundance of PBHs is exponentially sensitive to primordial non-gaussianities. Based on this result, ref. [18] claimsthat, in the context of single-field inflationary models which feature an approximate stationary inflection point, thegaussian approximation is hardly trustable when computing the PBH abundance.The goal of this work is to address the same question using a different computational strategy inspired by peak the-ory [21]. More precisely, we associate regions where the overdensity field takes values above the threshold for gravita-tional collapse with spiky local maxima of the curvature perturbation field, and compute the number density of thelatter using peak theory that we extend to include local non-gaussianities.Our main conclusion is that the impact of local non-gaussianities on the PBH abundance is far less important com-pared to what previously thought. We confirm that in models for PBH production (at least the class of models that weare going to consider), local non-gaussianities are sizeable enough to invalid the use of the the gaussian approximationto estimate their abundance. However we find that their impact is modest when translated in terms of the amplitudeof curvature power spectrum, namely it is enough to change it by a factor (cid:39) or smaller to obtain the same PBH abun-dance predicted by the gaussian calculation. This shift can be obtained by a small change of the parameters of theinflationary model.As a phenomenological application, we consider the impact of local non-gaussianities on the computation of theamplitude of the induced second-order gravitational-wave signal. Our main conclusion is that a careful treatment ofnon-gaussianities is needed in order to provide a reliable comparison with the expected experimental sensitivities offuture gravitational-wave interferometers. En route, we discuss the difference between threshold statistics and peaktheory, and we explain under which conditions (and why) peak theory gives a PBH abundance which is larger than theone computed by means of threshold statistics.The structure of this paper is as follows. ∗ In section II we introduce the problem and present our solution strategy. This section is paired with appendix Awhere we explain in more detail the cosmological interpretation of all quantities involved. ∗ In section III we discuss our main results and we explain the discrepancy with the previous literature. This sectionis paired with appendix B-G where we collect all relevant technical details. ∗ We conclude in section IV.
II. PROBLEM SETUP AND SOLUTION STRATEGY
Consider in position space h ( (cid:126)x ) = R ( (cid:126)x ) + α R ( (cid:126)x ) , (1)where α is a constant, R ( (cid:126)x ) is a gaussian scalar random field while h ( (cid:126)x ) is non-gaussian because of the presence of thenon-linear term on the right-hand side. In this case, non-gaussianities are called of local type because for a given (cid:126)x atwhich we evaluate h the amount of non-gaussianity is localized at the same position. We briefly discuss in appendix Athe physical interpretation of eq. (1).The first observation is that ∂ i h ( (cid:126)x ) = [ ∂ i R ( (cid:126)x )][1 + 2 α R ( (cid:126)x )] , (2)meaning that stationary points of R are also stationary points of h ( ∂ i R = 0 implies ∂ i h = 0 ).What is crucial, however, is that the nature (saddle points, maxima or minima) of these “shared” stationary pointsdepends on the sign of the factor α R . Consider the matrix of second derivatives evaluated at a stationary point (cid:126)x st (we will further indicate a minimum with (cid:126)x m and a maximum with (cid:126)x M ). One finds the Hessian matrix h ij ( (cid:126)x st ) = [ R ij ( (cid:126)x st )](1 + 2 α R st ) . To fix ideas, consider the simple case of two spatial dimensions (cid:126)x = { x, y } and the case in which the stationary point (cid:126)x m is a minimum of R .Minima of R are identified by two conditions. The first one, R xx ( (cid:126)x m ) R yy ( (cid:126)x m ) − R xy ( (cid:126)x m ) > separates extremafrom saddle points. The second one, R xx ( (cid:126)x m ) > and R yy ( (cid:126)x m ) > , separates minima from maxima. Since wehave h xx ( (cid:126)x m ) h yy ( (cid:126)x m ) − h xy ( (cid:126)x m ) = (1 + 2 α R m ) [ R xx ( (cid:126)x m ) R yy ( (cid:126)x m ) − R xy ( (cid:126)x m ) ] it is obvious that the condition R xx ( (cid:126)x m ) R yy ( (cid:126)x m ) − R xy ( (cid:126)x m ) > is also satisfied by h .On the contrary, since h xx ( (cid:126)x m ) = (1 + 2 α R m ) R xx ( (cid:126)x m ) and h yy ( (cid:126)x m ) = (1 + 2 α R m ) R yy ( (cid:126)x m ) , it is possible that aminimum of R becomes a maximum of h if α R m < . Viceversa, a maximum of R can become a minimum of h .The argument trivially generalizes to the more realistic case of three spatial dimensions.In peak theory, one computes the number density of maxima [21]. We are interested in the number density of maximaof the non-gaussian variable h . As argued before, identifying this quantity with the number density of maxima of R (based on the observation that h and R have the same stationary points) is not correct. Let us give a quantitativeargument to support this claim. From the previous discussion, it is clear that counting the maxima of R might be notenough. On the contrary, a simple modification could be the following. One should i) Count the maxima of R ; ii) Add the minima of R that, depending on the value of (1 + 2 α R m ) , become maxima of h ; iii) Subtract the maxima of R that, depending on the value of (1 + 2 α R M ) , become minima of h .Ref. [22] assumes i) . However, the two operations ii) and iii) do not balance between each others, and a sizablecorrection to i) will be introduced if α is large enough. In fig. 1 we show how the number density of maxima of thegaussian variable R changes (in percentage) as a function of α when ii) and iii) are implemented. Schematically, wecompute ∆ n max = ( R → maxima of h ) − ( R → minima of h ) R . (3)We obtain fig. 1 using gaussian peak theory (implementing the results of ref. [21], see appendix B). If we take α (cid:28) , ii) and iii) do not alter the estimate of i) . However, for sizable α (cid:38) . the change in the number density of maxima of R becomes evident.In situations of cosmological interest, the issue is further complicated by the fact that we are not really interestedin all maxima of h but only in those which are “spiky enough.” The reason is that the quantity which is relevant is thedensity contrast δ ( (cid:126)x, t ) (also dubbed overdensity field in the following) whose relation with h ( (cid:126)x ) (assuming radiationdominated epoch) reads [23] δ ( (cid:126)x, t ) = − (cid:18) aH (cid:19) e − h ( (cid:126)x ) (cid:20) (cid:52) h ( (cid:126)x ) + 12 h i ( (cid:126)x ) h i ( (cid:126)x ) (cid:21) , (4)where the time dependence comes from the scale factor a = a ( t ) and the Hubble rate H = H ( t ) while h does not de-pend on time because eq. (4) assumes perturbations to be on super-horizon scales. Eq. (4) can be thought as a Poissonequation in which h plays the role of gravitational potential while the density contrast can be written more precisely as δ ( (cid:126)x, t ) ≡ δρ ( (cid:126)x, t ) /ρ b ( t ) where ρ b ( t ) is the average background radiation energy density and δρ ( (cid:126)x, t ) = ρ ( (cid:126)x, t ) − ρ b ( t ) its perturbation. The physics-case that is relevant for the present study is the one in which the density contrast hasa peak localized in some region of space that is high enough to trigger the gravitational collapse into a black hole. Ifthe number of these peaks above threshold is large enough, these black holes can be part of dark matter. In the range (cid:46) M PBH [g] (cid:46) , a population of PBHs may account for the totality of dark matter observed in the Universetoday.Consider a peak of the overdensity field, located at some spatial point (cid:126)y pk . δ ( (cid:126)y pk , t ) = − (cid:18) aH (cid:19) e − h ( (cid:126)y pk ) (cid:20) (cid:52) h ( (cid:126)y pk ) + 12 h i ( (cid:126)y pk ) h i ( (cid:126)y pk ) (cid:21) (cid:39) − (cid:18) aH (cid:19) (cid:52) h ( (cid:126)y pk ) , (5)where in the second step we linearized in h . We follow here the approach of refs. [24, 25] in which the linear approxima-tion was adopted. Since the peak amplitude of the overdensity must be larger than some critical value δ c , we deducethe condition −(cid:52) h ( (cid:126)y pk ) (cid:38)
94 ( aH ) δ c , (6) We use the short-hand notation f k ≡ f ( (cid:126)x k ) , f i ( (cid:126)x k ) ≡ ∂ i f ( (cid:126)x ) evaluated at (cid:126)x k and f ij ( (cid:126)x k ) ≡ ∂ ij f ( (cid:126)x ) evaluated at (cid:126)x k for some generic function f . The flat spatial Laplacian is (cid:52) f ( (cid:126)x ) ≡ (cid:80) i f ii ( (cid:126)x ) . ��� ��� ��� ��� ��� ������������ α Δ [ % ] FIG. 1:
Percentage increase in the number density of maxima of R when we ii) add the minima of R that become maxima of h andiii) subtract the maxima of R that become minima of h (see eq. (3)). To make this (illustrative) plot we set σ = 1 and γ = 3 / (seeappendix B for definitions). on the curvature of h at the peak of δ .If we assume that local maxima of h coincide with peaks of δ (that is (cid:126)y pk (cid:39) (cid:126)x M ), then the condition in eq. (6) tellsthat only maxima of h which are “spiky enough” contribute to the formation of black holes. Of course, the assumptionthat local maxima of h coincide with peaks of δ requires some care. Ref. [26] argues, both analytically and numerically,that this assumption is well justified. However, ref. [26] only considers the case in which h is gaussian, that is, h = R with α = 0 in our case (but they include the presence of the non-linearities in eq. (5)). Since stationary points of R arealso stationary points of h , we tend to believe that the same conclusion holds true in the case with α (cid:54) = 0 but of coursethis is an important point that has to be checked explicitly.All in all, the strategy we shall follow in the course of this work is the following.First, we will compute the number density of maxima of the non-gaussian random field h that are “spiky enough”according to the condition in eq. (6). This requires a generalization of the work in ref. [21] such as to implement localnon-gaussianities. Second, we will check that these maxima are also peaks of the overdensity field. If this last point willturn out to be true, our computation of the number density of “spiky enough” maxima of h will provide the abundanceof peaks of the overdensity field that are large enough to form black holes. III. RESULTS AND DISCUSSION
We present in this section the main results of our analysis. In section III A we discuss how primordial non-gaussianities of local type alter the abundance of PBHs. In section III B we compare with the existing literature. Insection III C we (partially) include the effect of non-linearities in the relation between curvature and density perturba-tions.All technical details are collected in appendix A (where we discuss the origin of eq. (1) from a cosmological view-point), appendix B (where we discuss the gaussian limit), appendixes C and D (where we discuss how to construct thenon-gaussian part and the approximations that are involved), appendix E (where we give formulas for computing cu-mulants of generic order), appendix F (where we discuss how to compute the threshold value for collapse into blackholes) and appendix G (where we discuss non-linearities).
A. The abundance of PBHs in the presence of primordial non-gaussianities of local type
The quantity of central interest is the fraction of the Universe’s mass in the form of PBHs at the time of their formation.As customary in the literature, we indicate this quantity with β . The present-day fractional abundance of dark matterin the form of PBHs is given by (for a review, see ref. [27]) Ω PBH Ω DM = O (1) × (cid:18) β − (cid:19) (cid:20) g ∗ ( t f )106 . (cid:21) − / (cid:18) M PBH g (cid:19) − / , (7)where g ∗ ( t f ) is the number of relativistic degrees of freedom at the time of black hole formation (that we normalizeand set to its standard model value). Eq. (7) is defined modulo an overall O (1) factor whose precise value dependson the detail of the gravitational collapse that leads to black hole formation. In this paper we consider M PBH (cid:39) g; consequently, as an order-of-magnitude estimate, β (cid:38) − is excluded since it would imply overclosure of thepresent-day Universe, Ω PBH > Ω DM .We find the following formulafraction of the Universe’s mass in PBH in the presence of primordial local non-gaussianities β (cid:39) (cid:112) π (1 − γ ) (cid:34)(cid:90) ∞− ασ d ¯ ν (cid:90) ∞ x δ (¯ ν ) dx e − ¯ ν / f ( x ) e − ( x − x ∗ )22(1 − γ + (cid:90) − ασ −∞ d ¯ ν (cid:90) x δ (¯ ν ) −∞ dx e − ¯ ν / f ( x ) e − ( x − x ∗ )22(1 − γ (cid:35) (8)that we derive in detail in appendix B (as far as the gaussian limit is concerned), appendix C (where we discuss howto construct the non-gaussian part and the approximations that are involved) and appendix F (where we discuss howto compute the threshold value for collapse into black holes). In short: ∗ The parameter α , already defined in eq. (1), indicates the presence of local non-gaussianities (of quadratic type).The limit α → reproduces the gaussian result. A more physical interpretation of this parameter is given inappendix A. In concrete models of inflation which generate a sizable abundance of dark matter in the form ofPBHs (see, for instance, ref. [28]), we expect α (cid:39) [0 . ÷ . [18, 19].We derive our result based on peak theory. More precisely, we associate regions where the overdensity field takeslarge values with spiky local maxima of the comoving curvature perturbation, and compute the number densityof the latter using peak theory that we extend to include local non-gaussianities. Within this approach, eq. (8)represents an original result. ∗ The spectral moments σ j are defined by (see eq. (B8) and discussion in appendix B) σ j ≡ (cid:90) dkk P R ( k ) k j , (9)where P R ( k ) is the dimensionless power spectrum of the gaussian random field R . The a-dimensional parameter γ is defined as γ = σ /σ σ and takes values < γ < . In this paper we analyze two possible cases. In order toelucidate some intermediate results of our computational strategy, we use in appendix B a simple toy-model forthe power spectrum given by the log-normal function (see eq. (B24) and related discussion) P R ( k ) = A g √ πv exp (cid:20) − log ( k/k (cid:63) )2 v (cid:21) , (10)since in this case the spectral moments can be computed analytically and they are given by σ j = A g k j(cid:63) e j v .The three parameters { k (cid:63) , A g , v } in eq. (10) control, respectively, the position of the peak of the power spectrum,the peak amplitude of the power spectrum and its width. However, we remark that in single-field inflationarymodels the value of α that defines the amount of local non-gaussianities and the shape of the power spectrumare intimately related, and in general one can not take α as a free parameter and fix the power spectrum to aspecific functional form like the one introduced in eq. (10). A more realistic example is the following. Considerthe power spectrum defined by the piecewise function realistic power spectrum : P R ( k ) = P R ( k (cid:63) ) × (cid:16) kk (cid:17) n exp (cid:104) − log ( k /k (cid:63) )2 v (cid:105) for k < k exp (cid:104) − log ( k/k (cid:63) )2 v (cid:105) for k (cid:54) k (cid:54) k (cid:16) kk (cid:17) n exp (cid:104) − log ( k /k (cid:63) )2 v (cid:105) for k > k (11)with two power-law behaviors for k (cid:28) k (cid:63) and k (cid:29) k (cid:63) that, for k ≈ k (cid:63) , are connected by the log-normal functionin eq. (10). In this case, it is possible to show that the spectral index of the fall-off of the power spectrum afterthe peak at k = k (cid:63) is related to α by the relation n ≈ − α that is twice the value of the Hubble parameter η after the end of the ultra slow-roll phase. In our numerical analysis, therefore, we use the realistic power spec-trum in eq. (11) with the condition n = − α for fixed α . This provides the above-mentioned relation between This can be understood as follows. The modes that constitute the fall-off of the power spectrum after the peak are those for which the horizon- �� � ����� - � �� - � ����� - � �� - � �� - � �� - � / ★ ℛ ( ) / ℛ ( ★ ) - ���� - � � �� - ���� - � - � � � � � � � - ������ � - � η η � = - ���� η � = - ���� η � = - ���� FIG. 2:
Power spectra for the model in ref. [28] (red, with label k − . ), ref. [29] (blue, with label k − . ) and ref. [30] (green, with label k − . ) computed numerically by solving the Mukhanov-Sasaki equation (left panel). The slope of the power-law falloff after the peakis k η where η is the value of the Hubble parameter η (whose evolution, as function of the e -fold time N with N in the beginning ofthe ultra slow-roll phase, is shown in the right panel) after the end of the ultra slow-roll phase; η , in turn, is related to the parameter α that controls the size of local non-gaussianities via α = − η / (see appendix A). the amount of local non-gaussianities and the shape of the power spectrum. We remark that this point is oftenoverlooked in the literature. As a numerical check, we show in fig. 2 (left panel) the power spectra of comovingcurvature perturbations computed numerically for three inflationary models in which a ultra slow-roll phasetakes place. All power spectra are well described by the analytical ansatz in eq. (11). The slope of the power spec-tra after the peak is related to the value of the Hubble parameter η (shown in the right panel of the same figure)after the end of the ultra slow-roll phase which, in turn, controls the size of local non-gaussianities (see captionand appendix A for more details).As far as the values of the other parameters in eq. (11) are concerned, we use k (cid:63) = 1 . × Mpc − , n = 3 . , k = k (cid:63) / , k = 3 k (cid:63) / and v = 0 . ; the values of k , , v and n are motivated by a fit of eq. (11) done with respectto the numerical power spectrum obtained in the context of the explicit models studied in ref. [28]. In particular,notice that the spectral index n describes the growth of the power spectrum that leads to the formation of thepeak; its value is related to the value of the Hubble parameter η during the ultra slow-roll phase and the durationof the latter. Semi-analytical arguments (see ref. [31]) suggest that n < . The value k (cid:63) = 1 . × Mpc − ischosen because it implies M PBH (cid:39) (for which Ω PBH (cid:39) Ω DM is possible, see ref. [19]). We consider the peakamplitude of the power spectrum P R ( k (cid:63) ) as a free parameter. ∗ The function f ( x ) is given in eq. (C17) and the quantity x ∗ in eq. (C12). The function x δ (¯ ν ) is defined by therelation (see eq. (C14)) (1 + 2 ασ ¯ ν ) x δ (¯ ν ) = 9( a m H m ) σ δ c , (12)where δ c = O (1) is a threshold value above which a peak of the overdensity field collapses to form a black hole.The left hand side of eq. (12) corresponds to the critical curvature of h in eq. (6), see section C. Formally, eq. (12)depends on time via the comoving Hubble radius /aH . We evaluate eq. (12) at the time t m when curvatureperturbations re-enter the horizon and become causally connected (see refs. [24, 25] and appendix F). In eq. (12)we use the short-hand notation a m H m ≡ a ( t m ) H ( t m ) . crossing condition k = aH happens after the end of the ultra slow-roll phase [28]; during this part of the dynamics the power spectrum can beapproximated by means of the conventional slow-roll relation P R ( k ) = H / π (cid:15) where the Hubble parameter (cid:15) evolves in time according to (cid:15) ( N ) ∝ e − η N where η (which is a negative number, η < , see appendix A) is the value of the Hubble parameter η after the end of the ultraslow-roll phase. We neglect the contribution coming from the time-evolution of H , which is sub-leading. This means that we have P R ( k ) ∝ e η N . From k = aH we have dk/k = dN , and we can convert the e -fold time-dependence into a k -dependence, P R ( k ) ∝ k η = k − α where we used that α = − η / (see appendix A). However, see ref. [32] for a special case (derived in the context of non-attractor inflation) in which the growth of the power spectrum can be steeper,although for a limited range of k . ∗ An important comment concerns the so-called smoothing procedure. Consider the case in which one takes avery narrow power spectrum, like the toy-model introduced in eq. (10) with a small value of v , say v = 0 . . In thiscase it is not strictly necessary to introduce a smoothing procedure because the power spectrum is characterizedby a well-defined scale in momentum space, k = k (cid:63) . The realistic case introduced in eq. (11), on the contrary,requires more care. Although the power spectrum peaks at k = k (cid:63) , the peak is broadened by the relatively largevalue of v , and it also possesses a pronounced power-law tail at large k (cid:29) k (cid:63) . In this situation we can not blindlyapply eq. (8) to compute the PBH abundance because the spectral moment σ is formally ultraviolet-divergentunless the power spectrum decays fast enough, which is however not the case in the realistic model. The solutionto this issue (discussed in appendix F) is to smooth-out small scales by introducing an appropriate cut-off. At theoperative level, we use, instead of eq. (11), the power spectrum P cut R ( k ) ≡ P R ( k ) exp( − k /k ) and we choose k cut such as to minimize the threshold value in the right-hand side of eq. (12).Physically, the fact that the power spectrum in eq. (11) does not possess a well-defined scale means that the PBHsit generates will be characterized by a relatively broad mass distribution (rather than sharply peaked at the valueassociated to k (cid:63) ). The cut-off procedure described before selects the scale (and, therefore, the value of the mass M PBH ) at which the abundance of PBH will be the largest. ∗ We use the linear approximation in eq. (4). We discuss the role of non-linearities in section III C.We now discuss the implications of eq. (8). In the left panel of fig. 3 we show the fraction of Universe’s mass in the �� - � �� - � �� - �� �� - �� �� - �� �� - �� �� - �� �� - � �� - � � ℛ ( ★ ) β ℯℓ ℛ ( ) α = � ( ℊ ) α = ���� α = ���� α = ���� �� �� - � �� - � �� - �� �� - �� �� - �� �� - �� �� - �� �� - � �� - � � ℛ ( ★ ) β ℯℓ ℛ ( ) α = � ( ℊ ) α = ���� �� ℯ ℯ ℴ ℴ ℓ � FIG. 3:
Left panel: Fraction of the Universe’s mass in PBHs at the time of their formation computed with (solid lines with colors)and without (dashed black line) non-gaussianities as a function of the peak amplitude of the power spectrum P R ( k (cid:63) ) . We adopt therealistic model for the power spectrum introduced in eq. (11) and the abundance is computed using eq. (8). We show the impact ofnon-gaussianities for different benchmark values of α . In the hatched region we have β > − and the Universe is overclosed (seeeq. (7) and related discussion). Right panel: we consider two spatial dimensions, and compare the exact computation of β with theapproximation obtained including only the third-order cumulants. We also show the value of β obtained by means of the exponentialapproximation in eq. (D101); the latter gives an estimate of the abundance off by many orders of magnitude compared with the actualresult. form of PBH computed according to eq. (8) for increasing values of the parameter α starting from the gaussian casewith α = 0 . What values of α are expected in concrete models? In popular single-field models of inflation that generatea sizable abundance of dark matter in the form of PBHs, we find α (cid:39) . (ref. [28]), α (cid:39) . (ref. [33]), α (cid:39) . (ref. [34]), α (cid:39) . (ref. [30]), α (cid:39) . (ref. [29]). The plot shows that including local non-gaussianities of primordialorigin makes the formation of PBHs easier, and the value of P R ( k (cid:63) ) required to reproduce the benchmark abundance β = 10 − turns out to be smaller than the gaussian one. We find that the rescaling of the peak amplitude implied by Of course, any power spectrum generated by the inflationary dynamics has an intrinsic cut-off set by the smallest scale (largest k ) that exits theHubble horizon before inflation ends. More precisely, therefore, with the words “ultraviolet-divergent integral” we mean that σ is dominated bysmall scales. the presence of local non-gaussianities is modest, a factor a few in realistic models. Moreover, let us mention that itcan be obtained at a price of an even smaller retuning of the parameters of the inflationary models. Similar results havebeen obtained in [19] for the calculation of the PBH abundance with threshold statistics.In addition to the exact result presented in eq. (8), we also consider a “perturbative” approach based on a power-series α -expansion around the gaussian distribution. To this end, we work in two (instead of three) spatial dimensions.This simplifying assumption allows to perform most of the computations in appendix B and appendix D analytically,and makes possible to visualize and check numerically a number of intermediate results by means of simple two-dimensional plots.Let us first clarify the exact meaning of the word “perturbative.” From a statistical viewpoint, the α -expansion cor-responds to an expansion in cumulants of the joint non-gaussian probability distribution according to the schematicsummarized in table I (see appendix B, appendix D and appendix E for details). In the gaussian approximation, all cu- Gaussian O ( α ) O ( α ) O ( α ) O ( α ) O ( α ) . . .α = 0 second-order cumulants C (cid:51) (cid:55) (cid:51) (cid:55) (cid:55) (cid:55) . . . third-order cumulants C (cid:55) (cid:51) (cid:55) (cid:51) (cid:55) (cid:55) . . . fourth-order cumulants C (cid:55) (cid:55) (cid:51) (cid:55) (cid:51) (cid:55) . . . fifth-order cumulants C (cid:55) (cid:55) (cid:55) (cid:51) (cid:55) (cid:51) . . .. . . . . . . . . . . . . . . . . . . . . . . . TABLE I:
We organize the cumulants C n of the joint six-dimensional probability distribution P ( h, h x , h y , h xx , h xy , h yy ) as a seriesexpansion in α . In the case α = 0 , only second-order cumulants are non-vanishing (eqs. (D74-D77) with α = 0 ), and we reconstructthe gaussian limit. At order O ( α ) , the leading correction is given by third-order cumulants (eqs. (D82-D89)). At order O ( α ) , weinclude corrections to the second-order cumulants (eqs. (D74-D77)) and the leading pieces in the expression of fourth-order cumulants(eqs. (D102-D116)). mulants C n (cid:62) vanish and the second-order ones correspond to the entries of the covariance matrix of the distribution.This result is valid in the limit α → . If α (cid:54) = 0 , all cumulants C n (cid:62) are generated. However, the non-zero cumulants canbe organized in terms of an α -expansion as shown in table I (see caption, and appendix D and appendix E for details).Crudely speaking, we have C n (cid:62) ∼ O ( α n − ) + O ( α n ) . (13)At order O ( α ) , only the leading part of the third-order cumulants C appears. We can, therefore, consider an expansionaround the gaussian distribution including only the leading part of the third-order cumulants C . The rationale for thisapproximation is twofold. i) If we compare the perturbative approach with the exact computation (downgraded in two spatial dimensions)we can estimate the validity of the approximation in which only the leading part of the third-order cumulant C is included.This exercise is useful for the following reason. In the approach based on threshold statistics, one usually com-putes, using the tools of cosmological perturbation theory, the cumulants in the form of (the connected part of)correlators of the density perturbation field. In the presence of ultra slow-roll, however, this computation is notsimple (see refs. [18, 19]), and one typically includes only the leading term in the so-called bispectrum (that isthe three-point correlator). This approximation precisely corresponds to the one in which only the leading partof the third-order cumulant C is included (that is the one proportional to α in eq. (13)). In the right panel offig. 3 we show the comparison between the exact computation of β and the approximation in which only the More precisely, for each cumulant the expansion is controlled—considering for simplicity the log-normal power spectrum in eq. (10)—by thedimensionless parameter αA g (cid:28) . leading part of C is included; as mentioned before, we work in two spatial dimensions and we fix α = 0 . . Thecomparison shows that truncating the expansion at the first order in the non-gaussian corrections does not fullycapture the impact of non-gaussianities on β .Based on this result, we pose the attention on the fact that a similar conclusion is likely to be valid also whencomputing correlators in cosmological perturbation theory. The local non-gaussianity which is present at thelevel of the three-point correlator (computed, for instance, in refs. [18, 19] and used in ref. [18] to estimate theimpact of non-gaussianities on β ) induces corrections also at higher orders. In terms of (cosmological) Feynmandiagrams, the situation can be sketched as follows third-order cumulant fourth-order cumulant from (cubic interactions) pure fourth-order cumulantend of inflation (14)where the square of the third-order local interaction enters at the fourth-order (as well as at higher ones). Thecomparison shown in the right panel of fig. 3 suggests that, without including the contributions that third-orderlocal interactions induce in higher-order cumulants, a precise computation of β cannot be claimed.In addition to this observation there is a second, and by far more problematic, issue that we shall discuss next. ii) When computing non-gaussian corrections in the form of a series expansion around the gaussian distribution,some care must be taken. Seemingly harmless approximations, indeed, may lead to erroneous conclusions. Con-sider the simplified case discussed before in which only the leading part of the third-order cumulants C is in-cluded. A wrong approximation in the evaluation of the non-gaussian probability distribution leads to the ex-pression of the abundance that we report in eq. (D101). In this compact expression, the non-gaussian correctionexponentiates and alters the argument of the exponential function in the gaussian distribution. If this approx-imation were true, the effect of non-gaussianities would be exponentially large. For illustration, we add in theright panel of fig. 3 the value of β that one gets by means of eq. (D101). Clearly, eq. (D101) overestimates the ac-tual value of β by many orders of magnitude. In appendix D we explain why this approximation is wrong andwhat is the correct procedure to follow.The reader may wonder why we are wasting time discussing a wrong result. The reason is that it rings a bell.Previous studies on the impact of primordial non-gaussianities found, although in a slightly different statisti-cal context, that non-gaussian corrections alter exponentially the gaussian value of β , a conclusion that we justbranded inaccurate. A closer look at this literature, therefore, is mandatory. This will be the subject of the nextsection. B. Comparison with the literature
The impact of primordial non-gaussianities on the computation of PBH abundance was discussed in ref. [20] in thecontext of the so-called threshold statistics. Ref. [20] finds that primordial non-gaussianities may play a very relevantrole because they alter the argument of the exponential function that sets the value of the PBH abundance in the gaus-sian case. Ref. [18] applies the results of ref. [20] to the case of single-field inflationary models which feature the pres-ence of an approximate stationary inflection point, and concludes that the gaussian approximation does not give thecorrect estimate of the PBHs abundance precisely because the non-gaussian correction exponentiates and drasticallychanges the gaussian result.However, this conclusion clashes with our result. To explain the discrepancy, let us first re-derive the main result ofref. [20] in a simplified way.In a nutshell, in the context of threshold statistics one computes, assuming a gaussian distribution, the probabilityto find regions where the overdensity field δ takes values above a given threshold. It is known that threshold statisticsgives a smaller PBH abundance if compared with peak theory. Let us briefly explain the origin of this difference. Thispoint is not crucial to understand the discrepancy between our conclusions and ref. [20] but it will play an important Of course, higher-order correlators generated by pure higher-order interactions—for instance, a pure quartic interaction in the violet vertex of theexample above—could also be present but, as discussed in appendix A, we do not consider them explicitly in this paper. δ . Furthermore, we start considering the gaussian limit. The gaussian probability density distribution of δ with variance σ δ and spectral moments ¯ σ j (where ¯ σ = σ δ ) is given in the two cases by the expressions P threshold ( δ ) = 1 √ πσ δ e − δ / σ δ , (15) P peak ( δ ) = 1(2 π ) σ δ R ∗ (cid:26) (cid:90) ∞ dx f ( x ) (cid:112) π (1 − γ ) exp (cid:20) − ( x − δγ/σ δ ) − γ ) (cid:21) (cid:27)(cid:124) (cid:123)(cid:122) (cid:125) ≡ G ( γ,δ ) e − δ / σ δ ≡ π ) σ δ R ∗ (cid:90) ∞ dxP peak ( x, δ ) , (16)where the function f ( x ) is given explicitly in eq. (C17). The two quantities (cid:54) γ ≡ ¯ σ / ¯ σ σ δ < and R ∗ ≡ √ σ / ¯ σ are factors depending on the power spectrum of δ (notice in particular that γ is the same defined below eq. (9) butnow written in terms of the spectral moments of δ instead of R ). The variable x is defined by x ≡ −(cid:52) δ/ ¯ σ . The twoexpressions are similar but there are two differences. First, P peak ( δ ) has dimension of inverse spatial volume (becauseof the factor /R ∗ ). This is because P peak ( δ ) is defined in peak theory as a number density of maxima. This impliesthat the true comparison is between the a-dimensional quantities P threshold ( δ ) and R ∗ P peak ( δ ) . Second, P peak ( δ ) con-tains an extra factor (the one in curly brackets) with respect to P threshold ( δ ) . This factor arises because in peak the-ory one starts from the ten-dimensional gaussian joint probability density distribution P ( δ, δ i , δ ij ) of δ and its first( δ i ) and second ( δ ij ) spatial derivatives and imposes a number of conditions that select among the stationary pointsthose that are maxima [21]. In threshold statistics, on the contrary, one simply integrates out all spatial informationsby reducing P ( δ, δ i , δ ij ) to the one-dimensional gaussian distribution in eq. (15). The crucial aspect is that the func-tion G ( γ, δ ) in eq. (16) depends on δ in a way which is proportional to the parameter γ . The latter controls the de-gree of correlation between δ and x . When γ = 0 , the two variables are completely uncorrelated, and we find that G (0 , δ ) = (29 − √ / √ π . In this case, G (0 , δ ) does not depend on δ : peak theory and threshold statistics givethe same qualitative answer (in the sense that the only δ -dependence is encoded in the gaussian exponential whichis in common between the two). When γ → , δ and x are strongly correlated. This fact deforms the shape of the - � - � � � ������� ≈ �� - � γ = � - � - � � � ������� ≈ �� - � γ = ���� FIG. 4:
Isocontours of constant probability density distribution P peak ( x, δ ) defined in eq. (16). We use ν ≡ δ/σ δ and we set (forillustrative purposes only) the threshold at ν c = 3 . The probability above the threshold is obtained integrating P peak ( x, δ ) for x ∈ [0 , ∞ ) and ν ∈ [ ν c , ∞ ) . We show the case with γ = 0 (no correlation between x and δ , left panel) and γ = 0 . (strong correlationbetween x and δ , right panel). probability density distribution of x and δ and gives more weight to the region in which δ crosses the threshold forcollapse. This is illustrated schematically in fig. 4 (using ν ≡ δ/σ δ ). In the case of strong correlation between δ and x ,therefore, it is well expected that after integrating above the threshold ν > ν c peak theory gives a result which is largerthan threshold statistics. In threshold statistics, this information about the correlation between δ and x is completely Notice that, for the sake of simplicity, we are implicitly assuming in this example the same threshold value ν c in the case of peak theory andthreshold statistics in order to highlight the main difference between the two statistical approaches. δ is treated independently from the spatial configuration. We remark that the difference between thresholdstatistics and peak theory is more and more relevant as we approach the limit γ → . As already pointed out, the valueof γ is dictated by the properties of the power spectrum (in this case, strictly speaking, the power spectrum of δ whichis however related to the power spectrum of comoving curvature perturbations) with a very peaked power spectrumthat corresponds to the limit γ → .After this digression, we can back to the main point of this section. In the context of threshold statistics, non-gaussianities are included in the form of a Gram–Charlier A series around the gaussian ansatz [35] P NG ( δ ) = P G ( δ ) (cid:34) ∞ (cid:88) n =3 n !2 n/ B n (0 , , C , . . . , C n ) H n (cid:18) ν √ (cid:19)(cid:35) (17) = P G ( δ ) (cid:20) C / H (cid:18) ν √ (cid:19) + C / H (cid:18) ν √ (cid:19) + C / H (cid:18) ν √ (cid:19) + (10 C + C )6!2 / H (cid:18) ν √ (cid:19) + . . . (cid:21) , where P G ( δ ) = (1 / √ πσ δ ) exp( − δ / σ δ ) is the gaussian probability density distribution of δ with variance σ δ (thatis P threshold ( δ ) introduced in eq. (15)) and H n are the physicists’ Hermite polynomials. We define ν ≡ δ/σ δ and weintroduce the n -th complete exponential Bell polynomial B n ( x , . . . , x n ) = (cid:80) nk =1 B n,k ( x , . . . , x n − k +1 ) with B n,k the partial exponential Bell polynomials. We indicate with C n the n th normalized cumulant defined as the connectedpart of the n -point correlator (evaluated at the same point) of the overdensity field normalized by the n th power of thestandard deviation C n ≡ (cid:104) n times (cid:122) (cid:125)(cid:124) (cid:123) δ ( (cid:126)x ) . . . δ ( (cid:126)x ) (cid:105) conn σ nδ . (19)The explicit computation of eq. (19) in the presence of an ultra-slow roll phase is discussed, in the context of cosmo-logical perturbation theory, in ref. [19] (see also ref. [18]). Notice that the quantity (cid:104) δ ( (cid:126)x ) . . . δ ( (cid:126)x ) (cid:105) is dimensionless, andso is C n . We now integrate P NG ( δ ) over some threshold δ c = ν c σ δ . We define the abundance β ( ν c ) = (cid:82) ∞ δ c dδ P NG ( δ ) .We find β ( ν c ) = 12 Erfc (cid:18) ν c √ (cid:19) + 1 √ πν c e − ν c / ∞ (cid:88) n =3 √ ν c n !2 n/ B n (0 , , C , . . . , C n ) H n − (cid:18) ν c √ (cid:19) (21) (cid:39) √ πν c e − ν c / (cid:124) (cid:123)(cid:122) (cid:125) gaussian approx β G ( ν c ) (cid:34) ∞ (cid:88) n =3 √ ν c n !2 n/ B n (0 , , C , . . . , C n ) H n − (cid:18) ν c √ (cid:19)(cid:35)(cid:124) (cid:123)(cid:122) (cid:125) non gaussian correction , (22)where in the last step we use (1 / x/ √ (cid:39) (1 / √ πx ) e − x / for x (cid:29) .We now adopt the approach suggested in ref. [20] and approximate H n (cid:18) ν c √ (cid:19) = 2 n/ ν nc + O ( ν n − c ) , (23) We remind that d n dx n e − x / σ = ( − n n/ σ n e − x / σ H n (cid:18) x √ σ (cid:19) , with H n ( x ) = ( − n e x d n dx n e − x , (18)so that the reader can immediately recognize in eq. (17) the structure of a derivative expansion. We use the property √ πσ (cid:90) ∞ δ c dδe − δ / σ H n (cid:18) δ √ σ (cid:19) = 1 √ π e − ν c / H n − (cid:18) ν c √ (cid:19) . (20) ν c (cid:29) . In this case eq. (22) becomes β ( ν c ) (cid:39) √ πν c exp (cid:32) − ν c ∞ (cid:88) n =3 C n ν nc n ! (cid:33) . (27)We find, therefore, the main result of ref. [20]: The non-gaussian correction exponentiates, and changes the argumentof the exponential in the gaussian distribution.Consider the simplified case in which only the third-order normalized cumulant (a.k.a. skewness) is non-vanishing, C (cid:54) = 0 and C n> = 0 . This is the approximation studied in ref. [18] in which the third-order cumulant is computed inthe case of local non-gaussianities. We find β ( ν c ) (cid:39) √ πν c exp (cid:18) − ν c C ν c (cid:19) , with C (cid:54) = 0 and C n> = 0 . (29)Ref. [20] and ref. [18] used the above equations to conclude that the gaussian estimate of the PBH abundance is hardlytrustable. This result is based on the approximation in eq. (23). However, as we shall now discuss, the applicability ofthis approximation is not as straightforward as one may think. Let us critically inspect the issue.We define, for ease of reading, the quantity b n ( ν c ) ≡ ( √ ν c /n !2 n/ ) B n (0 , , C , . . . , C n ) H n − ( ν c / √ which entersin the non-gaussian correction in eq. (22). As a consequence of eq. (23), we have (slashed terms are neglected if we The complete exponential Bell polynomial is defined by the exponential generating function [36] exp ∞ (cid:88) j =1 x j t j j ! = ∞ (cid:88) n =0 B n ( x , . . . , x n ) t n n ! , (24)and the first few complete Bell polynomials are B = 1 , B ( x ) = x , B ( x , x ) = x + x , B ( x , x , x ) = x + 3 x x + x . In our casewe have x = x = 0 and the previous definition takes the form exp ∞ (cid:88) j =3 x j t j j ! = 1 + ∞ (cid:88) n =3 B n (0 , , x , . . . , x n ) t n n ! , (25)which applies to our case with t = ν c and x i (cid:62) = C i (cid:62) since the application of eq. (23) to eq. (22) gives β ( ν c ) (cid:39) √ πν c e − ν c / (cid:34) ∞ (cid:88) n =3 B n (0 , , C , . . . , C n ) ν nc n ! (cid:35) . (26)Notice that, compared to our result, ref. [20] finds an additional factor ( − n in eq. (27) which, however, does not appear in our computation. Notice, however, that assuming C (cid:54) = 0 and C n> = 0 is not fully consistent with the local non-gaussianities computed in ref. [18]. Indeed, localnon-gaussianities automatically generate non-vanishing cumulants of any order. This is evident in our approach, in which α (cid:54) = 0 generates awhole tower of non-zero cumulants. From a more mathematical viewpoint, consider the following theorem. The sequence { µ n , n = 0 , , , . . . } corresponds to moments of a non-negative probability density function if and only if the determinants D n +1 ≡ det µ µ µ · · · µ n µ µ µ · · · µ n +1 µ µ µ · · · µ n +2 ... ... ... ... ... µ n µ n +1 µ n +2 · · · µ n , n = 0 , , , . . . , (28) are all non-negative [37] . Consider the illustrative case with C n = 0 for n (cid:62) . We have µ = 1 , µ = 0 , µ = σ δ , µ = (cid:104) δ (cid:105) and µ n = σ nδ ( n − for n (cid:62) even ( µ n = 0 for n (cid:62) odd). The last condition simply means that cumulants of order higher than three vanish (and the correspondingmoments purely gaussian). The first determinant D > already sets a non-trivial condition on the skewness, that is −√ < C < √ .However, there is more than this. It is indeed trivial to see that higher-order conditions D n > for n > impose increasingly strong bound on C so that only C → is allowed if all the infinite number of constraints are implemented. In other words, the theorem above implies that skewnessalone can not consistently parameterize a non-gaussian probability density function . This result resonates with what we discussed at point ii) insection III A. Computing only the three-point correlator of the overdensity field does not fully describe non-gaussianities but one should at leastinclude the contributions that the three-point function generates at higher orders. b = C ν c (cid:18) − (cid:1)(cid:1)(cid:1) ν c (cid:19) , (30) b = C ν c (cid:18) − (cid:26)(cid:26)(cid:26)(cid:26) ν c + 15 ν c (cid:19) , (31) b = C ν c (cid:18) − (cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24) ν c + 210 ν c − ν c + 105 ν c (cid:19) , (32) b = . . . (33)and so on. We see, for instance, that we neglected the term − C ν c / in eq. (32) but we kept the term C ν c / ineq. (30). This is hardly justifiable given that ν c (cid:29) and we expect C ∼ O (1) (see ref. [19] for a careful computation of C ). This is enough to question the validity of the result based on eq. (23). Furthermore, from this simple comparison,we also see that the approximation in eq. (23)—contrary to what naïvely expected—gets worse for larger ν c . � � � � � � � � �� ���� - �� �� - �� �� - �� �� - � �� - � � ν β ( ν ) ℊℯ ℴ ℯ ℯ ℴ � = ���� �� - � �� - � �� - �� �� - �� �� - �� �� - �� �� - �� �� - � �� - � � ℛ ( ★ ) β ℊℯ ℴ ℯ ℯ ℴ FIG. 5:
Left panel. Numerical comparison between i) the gaussian result β G ( ν c ) defined in eq. (22), ii) the exponential approximationin eq. (29) and iii) the exact expansion of β ( ν c ) = β G ( ν c )[1 + (cid:80) n max n =3 b n ( ν c )] for increasing values of n max . We take C = 0 . andwe show our results as function of ν c . Right panel. Same as in the left panel but using the scaling in eq. (37) to draw β as function of P R ( k (cid:63) ) . We take δ c = 1 . . In particular, in this specific case, we can rewrite the sum (21) by using the recurrence property of the Bell polyno-mials: B n +1 ( x , . . . , x n +1 ) = n (cid:88) k =0 (cid:18) nk (cid:19) B n − k ( x , . . . , x n − k ) x k +1 . (34)In our particular case x k +1 = C δ k +1 , , therefore B k +1) (0 , , C , , . . . ,
0) = (3 k +1)(3 k +2) C B k (0 , , C , , . . . , .This can be easily recasted in the form: B n (0 , , C , , . . . ,
0) = (cid:40) (cid:0) C (cid:1) k (3 k − k − , if n = 3 k, , if n (cid:54) = 3 k. (35)Using this identity inside the eq. (21) we get: β ( ν c ) = 12 Erfc (cid:18) ν c √ (cid:19) + 1 √ π e − ν c / ∞ (cid:88) k =1 k ! (cid:18) C √ (cid:19) k H k − (cid:18) ν c √ (cid:19) . (36)4To show more explicitly the error that one makes by taking the exponential approximation in eq. (29) we plot inthe left panel of fig. 5 the comparison between i) the gaussian result β G ( ν c ) defined in eq. (22), ii) the exponentialapproximation in eq. (29) and iii) the exact expansion of β ( ν c ) = β G ( ν c )[1 + (cid:80) n max n =3 b n ( ν c )] for increasing values of n max . We take C = 0 . and we set C n> = 0 . The approximation gives a wrong estimate of the actual magnitude ofthe non-gaussianities. As expected, the exponential approximation diverges from the actual result for larger values of ν c .There is another point that is worth emphasizing. Taking C fixed and changing only ν c , as done in the left panel offig. 5 is not fully consistent since both C and ν c are functions of the power spectrum P R . More explicitly, we have thescaling ν c ∝ P R ( k (cid:63) ) − / and C ∝ P R ( k (cid:63) ) / (see appendix A). This means that reducing P R ( k (cid:63) ) to decrease C doesnot improve on the applicability of the exponential approximation since ν c , in turn, increases. To better visualize thispoint let us consider the following benchmark scalings σ δ = 0 . (cid:20) P R ( k (cid:63) )0 . (cid:21) / , C = 0 . (cid:20) P R ( k (cid:63) )0 . (cid:21) / , (37)and take δ c = 1 . . We can now redo the comparison we did before but now as function of P R ( k (cid:63) ) . The result is shownin the right panel of fig. 5. We conclude again that the exponential approximation overestimates the impact of localnon-gaussianities on the PBH abundance by many orders of magnitude compared with the actual result. Using theexponential approximation, one would wrongly conclude that, in order to fit the reference value β (cid:39) − , an order-of-magnitude decrease in the peak amplitude of the power spectrum is needed compared to the gaussian result.A precise determination of the value of P R ( k (cid:63) ) that is needed to get Ω DM (cid:39) Ω PBH has important phenomenologicalimplications. Let us discuss a specific example. PBHs with mass M PBH (cid:39) g have approximatively the size of theatomic nucleus (remember that M (cid:12) (cid:39) .
48 km in the Planck unit system). Testing the nature of such small objectsin the form of dark matter is extremely challenging. An interesting prospect is the following. The peak in the powerspectrum of scalar perturbations that is responsible for the formation of PBHs also generates (as a second-order effect)a gravitational wave signal. The position of the peak amplitude of the power spectrum of curvature perturbations, thepeak height in the mass distribution of PBHs and the frequency of the peak of the induced gravitational wave signal ( f (cid:63) in the following) are related by the approximate relation (cid:18) M PBH g (cid:19) − / (cid:39) k (cid:63) × Mpc − (cid:39) f (cid:63) . , (38)meaning that a peak amplitude around k (cid:63) = 10 Mpc − corresponds to a frequency range detectable by futuregravitational-wave interferometers like LISA, DECIGO and MAGIS-100. This is, in principle, a powerful probe. Bymeans of the condition Ω DM (cid:39) Ω PBH , one can fix the parameters of a given inflationary model that produces a sizableabundance of dark matter in the form of PBHs; in turn, this “predicts” the induced signal of gravity waves since the lat-ter is completely determined once the former step is taken. This is because the current energy density of gravitationalwaves induced by scalar perturbations is given by Ω GW = c g Ω r (cid:90) √ dt (cid:90) ∞ √ ds (cid:20) ( t − / s − / t − s (cid:21) (cid:2) I c ( t, s ) + I s ( t, s ) (cid:3) P R (cid:34) k √
32 ( s + t ) (cid:35) P R (cid:34) k √
32 ( s − t ) (cid:35) , (39)where c g ≈ . , Ω r is the current energy density of radiation and I c and I s are two functions that can be computedanalytically (see, for instance, ref. [38] and references therein); Ω GW is completely fixed once P R ( k ) is given. In mostof the analysis the power spectrum P R ( k ) is fixed to a specific functional form, and the abundance of PBHs computedassuming gaussian statistics. Typical functional forms for P R ( k ) used in the literature are the delta-function powerspectrum (see, e.g., ref. [39]), the log-normal power spectrum (see, e.g., refs. [40, 41]) and the more realistic brokenpower-law in eq. (11) (see, e.g., ref. [42]); all these studies assume gaussian statistics for PBH formation.The first remark that we make is that, as already anticipated, taking a specific functional form for P R ( k ) and assuminggaussian statistics can be conceptually wrong; for instance, for single-field inflationary models of phenomenologicalrelevance, eq. (11) with a fixed value of n = − α implies the presence of local non-gaussianities, and the value of β should be computed accordingly. The second remark is that if one takes the exponential approximation in eq. (29) One can argue that the power series (cid:80) ∞ k =1 z k k ! H k − ( x ) contained in eq. (36) is not convergent. Nevertheless it should be regarded as an asymp-totic series that gives a trustable value after an optimal truncation. Moreover it is Borel summable, and its value can be computed numericallyevalutating the Borel sum (cid:82) ∞ e − t (cid:20)(cid:80) ∞ k =1 ( tz ) k ( k !) H k − ( x ) (cid:21) dt . This check has been done and the value obtained agrees with the one got by theoptimal truncation. P R ( k (cid:63) ) . In reality (see fig. 5), primordial non-gaussianities reduce the value of P R ( k (cid:63) ) by no morethan a factor of 2 compared to the gaussian result. We will come back in more detail to this point in the next section. C. Non gaussianities from the non-linear relation between density and curvature perturbations.
In the previous sections we have assumed a linear relation between the density and curvature perturbations. This isbecause our main focus was the study of the non-gaussianities of primordial origin. Now we are going to include alsothe effect of the non-linearities in eq. (4). In our approach, they enter in two ways. First, they modify the relation be-tween the overdensity field and the curvature perturbation; this modification translates into a sort of “renormalization”of the value of δ c whose net effect is to make the production of black holes harder.Importantly, non-linearities also change the shape of the profile of the overdensity peaks which eventually collapseinto black holes. Since the threshold value depends on the shape of the overdensity, non-linearities also affect theway in which δ c is computed. Ref. [43] concludes that δ c is robust against non-linearities since its value changes at thepercent level. However, ref. [43] only considers the idealized case of a δ -function power spectrum. Deviations fromthe δ -function limit are considered in refs. [44–47] using different procedures whose consistency with each other is notyet fully understood. We do not aim to clarify this issue in the present work; however, we point out that refs. [44–47]assume that the curvature perturbation field is gaussian while, at least for the class of models which are relevant for thepresent analysis, primordial non-gaussianities should be included as well.A comprehensive study of the interplay between non-gaussianities of primordial origin and non-linearities in eq. (4)is left for future work. Here we shall focus on the first effect above, and assume that the threshold δ c is not signifi-cantly affected by the presence of non-linearities. Under this assumption, in appendix G we show how to compute thePBH abundance accounting for primordial non-gaussianities of local type and non-linearities between density andcurvature perturbations. We find that one can still use eq. (8), but replacing the threshold x δ (¯ ν ) with x NL δ (¯ ν ) , definedas x NL δ (¯ ν ) ≡ a m H m ) σ δ c ασ ¯ ν exp (cid:8) σ [¯ ν + ασ (¯ ν − (cid:9) . (40) ∗ Consider first the gaussian case α → . We have σ x NL δ (¯ ν ) | α =0 = 9( a m H m ) δ c e σ ¯ ν to be compared withits linear limit σ x δ (¯ ν ) | α =0 = 9( a m H m ) δ c . The main difference, as already emphasized in ref. [26], is thatnon-linearities “renormalize” the value of δ c by means of the exponential factor e σ ¯ ν . This means that whenconsidering maxima of the curvature field their height, in addition to their curvature, becomes important at thenon-linear level. We remark that these two quantities are statistically correlated, and the amount of correlationis controlled by the power spectrum via the a-dimensional parameter γ (see appendix B). If the power spectrumis very narrow, maxima with large curvature are likely to have also large ¯ ν ; consequently, the threshold value δ c effectively increases because of the factor e σ ¯ ν . This effect decreases the abundance with compared to thegaussian case. However, if the power spectrum is not very narrow, the correlation mentioned before is less ac-centuated, and maxima with large curvature are not necessarily characterized by large ¯ ν ; in this case, δ c increasesby a smaller factor and the reduction of the abundance turns out to be less important. Our realistic power spec-trum in eq. (11) belongs to the last case. In the left panel of fig. 6 we compare the gaussian computation of β inthe case without (dashed black line) and with (dot-dashed black line, label “NL”) non-linearities. As expected, β decreases but the change in the peak amplitude of the power spectrum that is needed to reproduce the referencevalue β (cid:39) − is less than a factor of 2 (while very narrow power spectra may increase the peak amplitude ofthe power spectrum up to one order of magnitude [26]). ∗ We take α (cid:54) = 0 . As discussed in section III A, primordial non-gaussianities modify the statistics of the overden-sity peaks. More peaks form, and, at constants threshold for collapse, the abundance of PBHs increases. Non-linearities, however, tend to enhance the value of δ c in the sense discussed in the previous paragraph, and thisgives the opposite effect of reducing the abundance. Moreover, when both non-linearities and primordial non-gaussianities are present the renormalization of δ c also depends on the value of α , as shown in eq. (40). We have,therefore, two competing effects: More overdensity peaks form but their collapse becomes less probable. Wefind that the result is a net decrease of the PBH abundance. This is shown in the left panel of fig. 6 (solid colorlines, label “PNG+NL”). Quantitatively, in order to fit the reference value β (cid:39) − the combination of primor-dial non-gaussianities and non-linearities require a peak amplitude of the power spectrum larger compared toits gaussian value. However, we find that P R ( k (cid:63) ) only increases by a factor which is at most 2 in all the realisticcases analyzed in this paper.Before concluding, let us explain the reason why it is important to have control on the computation of the PBHabundance at the level of the effects discussed in the present work.6 �� - � �� - � �� - �� �� - �� �� - �� �� - �� �� - �� �� - � �� - � � ℛ ( ★ ) β ℯℓ ℛ ( ) α = � ( ℊ ) α = ���� α = ���� �� �� ��� + �� �� - � �� - � �� - � �� - � � �� �� � �� - �� �� - �� �� - �� �� - �� �� - �� �� - �� �� - �� �� - � �� - � �� - � �� - � �� �� �� �� �� �� ��������� [ �� ] � Ω � � [ ��� - � ] ���� � � � � � - ��� �� ������� + �� ������ ������ ( ��� ) � Ω � � � + � � � � � � � � � + � � α = ���� ℯ ℴ FIG. 6:
Left panel. Same as in the left panel of fig. 3; we compare the gaussian result (dashed black line) with i) the computation thatincludes non-linearities between curvature and density perturbations (dot-dashed line, label “NL”) and ii) the computation that in-cludes both primordial non-gaussianities and non-linearities (solid clored lines, label “PNG+NL”). Right panel. Fraction of the energydensity in gravitational waves relative to the critical energy density of the Universe as a function of the frequency (see eq. (39)). Wecompare the gaussian result (dashed black line) with the computation that includes non-linearities between curvature and densityperturbations and the computation that includes both primordial non-gaussianities and non-linearities. We also consider for illus-tration the comparison with the exponential approximation (dot-dashed line, see the right panel of fig. 3). We take α = 0 . . Wesuperimpose the signal on the expected sensitivity curves of the future gravitational wave detectors LISA (assuming the C1 configura-tion [48]), DECIGO [49], MAGIS-100 [50] and the Einstein Telescope [51]. We also show present bound [52] and future prospect [53] forthe aLIGO experiment. The diagonal gray band is the stochastic gravitational wave background from binary black holes (BBH) andbinary neutron stars (BNS), which has been computed following ref. [54]. The only condition Ω DM (cid:39) Ω PBH (which is β (cid:39) − under the assumptions stated at the beginning of this section)is not enough to make the presence of non-gaussianities “phenomenologically observable”. With this term in quotationmarks we mean that the presence (or absence) of non-gaussianities can be simply reabsorbed in a re-tuning of P R ( k (cid:63) ) (which is equivalent to a very modest re-tuning of model parameters, see ref. [19]) if the only goal is to reproduce β (cid:39) − .What changes the rules of the game is the fact that, once P R ( k (cid:63) ) is fixed, the amplitude of the induced gravity wavesignal in eq. (39) is also predicted. Importantly, this implies that partial or inaccurate treatments of non-gaussianitieswill result in different values of P R ( k (cid:63) ) (required by imposing the same condition β (cid:39) − ) and, consequently, differ-ent gravitational wave signals. The (inaccurate) exponential approximation in eq. (27) provides a paradigmatic exam-ple. As discussed, if taken at face value the exponential approximation in eq. (27) may imply a one order-of-magnitudereduction in the value of P R ( k (cid:63) ) which gives β (cid:39) − (see, e.g., fig. 3 or fig. 5). Consequently, the amplitude of theinduced gravity wave signal in eq. (39) would be suppressed by two orders-of-magnitude.In the right panel of fig. 6 we compare the value of h Ω GW that one gets using peak theory in the gaussian limit(dashed black line) and the exponential approximation (dot-dashed black line; we take α = 0 . , and consider non-gaussianities at the level of the leading part of third-order cumulants). This example is tailored to show that in this casethe gravitational wave signal would drop below the expected stochastic background due to mergers of astrophysicalblack holes and neutron stars (the diagonal gray band). In the same plot we also show the value of h Ω GW that we get inthe presence of primordial non-gaussianities (label “PNG” this time computed according to eq. (8)) and in the presenceof both primordial non-gaussianities and non-linearities (label “PNG+NL” computed according to our discussion insection III C). As expected, primordial non-gaussianities (the combination of primordial non-gaussianities and non-linearities) decrease (increases) the amplitude of the gravitational wave signal compared to the gaussian result but thenet effect in both cases does not exceed the order of magnitude. As an order-of-magnitude approximation, therefore,the gaussian limit can be considered trustable (while using the exponential approximation would hide the signal belowthe stochastic background).There is one last point that is worth emphasizing. In the gaussian approximation, peak theory gives a PBH abundancewhich is systematically larger than the one provided by threshold statistics (this is in agreement with the expectation of7ref. [55] and discussed in section III B). The amplitude of the gravitational wave signal in fig. 6 is indeed smaller than theone computed with threshold statistics (see ref. [28]) and dangerously closer to the expected stochastic background.For this reason, we argue that keeping under control the impact of non-gaussianities could play an important part inthe interpretation of possible future detection. IV. CONCLUSIONS AND OUTLOOK
We conclude summarizing the main results and novelties of our work. ◦ We have studied the impact of primordial non-gaussianities of local type on the abundance of PBHs. For thispurpose we have focused on the maxima of the the comoving curvature perturbations. We have shown thatthose peaks “spiky enough”, i.e. with a curvature larger than a certain threshold, are good proxies for the maximaof the density field which undergo gravitational collapse, and produce PBHs. Exploiting this observation, we haveobtained that the cosmological abundance of PBHs can be computed using eq. (8). Our calculations extend thegaussian peak theory formalism [21] to include the effect of local non-gaussianities.We have examined the consequences of our results for models of single-field inflation with an ultra slow-rollphase. These scenarios have been investigated for the production of PBHs. In this context, the impact of primor-dial non-gaussianities is not dramatic. In fact, the desired PBH abundance can be obtained with an amplitudeof the power spectrum of curvature perturbations which is only a factor (cid:46) smaller from the value inferred withthe gaussian approximation. ◦ In parallel to the exact result presented in eq. (8), we have developed a computational strategy that approximatesthe effect of primordial non-gaussianities in the form of an expansion in cumulants around the gaussian result.In order to have full analytical control, we worked out this part in two (instead of three) spatial dimensions. Thepresence of local non-gaussianities affect cumulants at any order, and we have shown what is their structure.Furthermore, we have shown that the only inclusion of third-order cumulants is not sufficient to fully capturethe effect of local non-gaussianities. ◦ Previous works have studied primordial non-gaussianities in the context of threshold statistics, finding that non-gaussianities exponentially affect the PBH abundance as given by eq. (27). We have found that that result is notcorrect and largely overestimate the PBH abundance. ◦ In addition to non-gaussianities of primordial origin, we have also considered non-gaussianities that originatefrom the non-linear relation between curvature and density perturbations. In this case we have two competingeffects. On the one hand, primordial non-gaussianities alter the statistics of peaks enhancing the abundanceof PBHs as described by eq. (8). On the other one, primordial non-gaussianities and non-linearities change thecondition for collapse. At the linear level and without accounting for primordial non-gaussianities, the conditionfor collapse only involves the laplacian of the curvature perturbation at the position of the peak of the overdensityfield while in the presence of non-gaussianities it takes the form of eq. (40). This alteration makes the formation ofPBHs harder. When combining the two effects, we find that non-gaussianities suppress PBH production but thedesired PBH abundance can be obtained with an amplitude of the power spectrum of curvature perturbationswhich is only a factor (cid:46) larger from the value inferred with the gaussian approximation. ◦ En route, we discuss the difference between threshold statistics and peak theory, and we explain under whichconditions (and why) peak theory gives a PBH abundance which is larger than the one computed by means ofthreshold statistics.Finally, let us mention that our treatment of non-gaussianities can be further improved. We have computed thethreshold for collapse without including the effect of non-gaussianities (both primordial and non-linear) on the shapeof the collapsing peak (technically speaking, in eq. (F2) we used the averaged density profile—while in principle peaktheory also gives informations about the shape of the peak—and the linear approximation in eq. (F1)). Including theseeffects could give an additional modification of δ c . In light of the motivations discussed before, it would be importantto keep under control these additional effects. We will address this issue in a forthcoming work focused on the interplaybetween the computation of the black hole abundance and the induced gravitational wave signal. 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.8
Appendix A: Scalar perturbations during inflation in the presence of ultra slow-roll
In this paper we focus on single-field models of inflation. We indicate with φ the canonically normalized inflatonfield and with U = U ( φ ) its potential. The dynamics of φ can be obtained solving the equation of motion d φdN + 3 dφdN − (cid:18) dφdN (cid:19) + (cid:34) − (cid:18) dφdN (cid:19) (cid:35) d log Udφ = 0 , (A1)with slow-roll initial conditions; N indicates the number of e -folds defined by dN = Hdt , where H ≡ ˙ a/a is the Hubblerate, a the scale factor of the Friedmann-Lemaître-Robertson-Walker metric, and t the cosmic time (with ˙ ≡ d/dt ). Wewill also use in the following the conformal time τ defined by means of dt/dτ = a or, equivalently, dN/dτ = aH (noticethat, in the limit in which H is constant, we have the relation τ = − /aH ; the conformal time, therefore, is negative,and late times towards the end of inflation can be formally identify with the limit τ → − ).Scalar fluctuations can be efficiently described in terms of the Mukhanov-Sasaki field variable u ( τ, (cid:126)x ) which isa gauge-invariant combination of both fluctuations of the inflaton field and scalar fluctuations of the backgroundFriedmann-Lemaître-Robertson-Walker geometry. The Mukhanov-Sasaki field variable solves the differential equa-tion (for a review, see ref. [56]) d udτ = (cid:18) (cid:52) + 1 z d zdτ (cid:19) u , (A2) z d zdτ = a H (cid:20) (1 + (cid:15) − η )(2 − η ) + 1 aH (cid:18) d(cid:15)dτ − dηdτ (cid:19)(cid:21) , (A3)with z ≡ (1 /H )( dφ/dτ ) and (cid:52) the laplacian acting on spatial coordinates. The Hubble parameters are defined by (cid:15) ≡ − ˙ H/H and η ≡ − ¨ H/ H ˙ H .The field u can be quantized by defining the operator ˆ u ( τ, (cid:126)x ) = (cid:90) d (cid:126)k (2 π ) (cid:104) u k ( τ ) a (cid:126)k e + i(cid:126)k · (cid:126)x + u ∗ k ( τ ) a † (cid:126)k e − i(cid:126)k · (cid:126)x (cid:105) , (A4)with the annihilation and creation operators that satisfy the commutation relations of bosonic fields [ a (cid:126)k , a (cid:126)k (cid:48) ] = [ a † (cid:126)k , a † (cid:126)k (cid:48) ] = 0 , [ a (cid:126)k , a † (cid:126)k (cid:48) ] = (2 π ) δ (3) ( (cid:126)k − (cid:126)k (cid:48) ) , a (cid:126)k | (cid:105) = 0 , (A5)where the last condition defines the vacuum. The equation of motion for each mode u k ( τ ) takes the form of aSchrödinger equation d u k dτ + (cid:18) k − z d zdτ (cid:19) u k = 0 , (A6)which can be solved imposing Bunch-Davies initial conditions at some initial time when k (cid:29) aH . This choice definesthe initial vacuum state as the minimum energy eigenstate for an harmonic oscillator with time-independent frequencyin a space-time which is locally flat (because we are at length scales much smaller than the de Sitter curvature radius).Quantization of scalar perturbations is easier in terms of the Mukhanov-Sasaki field variable u ( τ, (cid:126)x ) since the quadraticaction corresponding to the equation of motion in eq. (A2) is the action describing a canonically normalized free fieldwith an effective time-dependent mass. However, the dynamics of the modes u k in Fourier space is more transparentif we introduce the so-called comoving curvature perturbation R = u/z ; in analogy with eq. (A4), R can be promotedto a quantum operator ˆ R ( τ, (cid:126)x ) = (cid:90) d (cid:126)k (2 π ) (cid:104) R k ( τ ) a (cid:126)k e + i(cid:126)k · (cid:126)x + R ∗ k ( τ ) a † (cid:126)k e − i(cid:126)k · (cid:126)x (cid:105) , (A7)with R k = u k /z that satisfies the time-evolution equation d R k dN + (3 + (cid:15) − η ) d R k dN + k a H R k = 0 , (A8)which is the differential equation of a damped harmonic oscillator.In the canonical picture of slow-roll inflation (that is for small Hubble parameters (cid:15), η (cid:28) ), eq. (A8) can be solvedin two complementary regimes divided by the so-called horizon-crossing condition k = aH . At early times, when k (cid:29) aH (in such case the modes are called sub-horizon), the last term dominates over the friction one, and the so-lution oscillates. As time passes by during inflation, the comoving Hubble radius /aH shrinks, the horizon-crossing9condition is met, and one eventually enters in the regime characterized by k (cid:28) aH (in such case the modes are calledsuper-horizon); the last term in eq. (A8) can be neglected and the latter admits the constant solution d R k /dN = 0 . Inother words, after horizon crossing the mode R k freezes to a constant value that it maintains in time until the end ofinflation (more precisely, until the mode re-enters the conformal Hubble horizon after the end of inflation).In the canonical picture of slow-roll inflation, after horizon crossing quantum fluctuations can be regarded as clas-sical, which motivates the description of cosmological perturbations in terms of classical random fields [57]. We willgive a more precise definition of random fields in appendix B. At the conceptual level, the previous statement impliesthat one has the schematic relation lim k/aH (cid:28) ˆ R ( τ, (cid:126)x ) = R ( (cid:126)x ) meaning that for super-horizon modes the quantumoperator ˆ R can be interpreted as a classical random field R ( (cid:126)x ) ; notice that the latter is time-independent because themodes R k are frozen in time. The previous relation can be formulated in a more precise form by saying that vacuumexpectation values on the quantum side are interpreted as statistical averages on the classical side. Formally, in theclassical picture we still have the Fourier decomposition R ( (cid:126)x ) = (cid:90) d (cid:126)k (2 π ) (cid:16) R k a (cid:126)k e + i(cid:126)k · (cid:126)x + R ∗ k a † (cid:126)k e − i(cid:126)k · (cid:126)x (cid:17) , (A9)but now creation and annihilation operators are no longer q -numbers but stochastic c -numbers which are defined bythe statistical averages (cid:104) a (cid:126)k a † (cid:126)k (cid:48) (cid:105) = 12 (2 π ) δ (3) ( (cid:126)k − (cid:126)k (cid:48) ) = (cid:104) a † (cid:126)k (cid:48) a (cid:126)k (cid:105) , (cid:104) a (cid:126)k a (cid:126)k (cid:48) (cid:105) = (cid:104) a † (cid:126)k a † (cid:126)k (cid:48) (cid:105) = 0 . (A10)Notice that (cid:104) a (cid:126)k a † (cid:126)k (cid:48) (cid:105) = (cid:104) a † (cid:126)k (cid:48) a (cid:126)k (cid:105) is only valid for classical fluctuations, since it implies that the stochastic parameterscommute. Furthermore, the stochastic parameters a (cid:126)k , a † (cid:126)k have zero mean value (consequently, (cid:104)R ( (cid:126)x ) (cid:105) = 0 ) in orderto match the fact that the quantum operator ˆ R ( τ, (cid:126)x ) in eq. (A7) has zero vacuum expectation value.The random field R ( (cid:126)x ) is fully specified by the entire hierarchy of its correlation functions. The simplest one isthe two-point correlation function. The latter can be defined by introducing the idea of power spectrum. The powerspectrum ∆ R ( k ) is defined by the Fourier transform of the two-point correlation function (cid:104)R ( (cid:126)x ) R ( (cid:126)x + (cid:126)r ) (cid:105) = (cid:90) d k (2 π ) e i(cid:126)k · (cid:126)r ∆ R ( k ) = (cid:90) ∞ dk sin( kr ) kr k π ∆ R ( k ) . (A11)The fact that ∆ R ( k ) depends only on k ≡ | (cid:126)k | and the explicit angular integrations that we performed in eq. (A11) areconsequences of the assumptions of spatial homogeneity and isotropy (equivalently, spatial homogeneity and isotropyimply that (cid:104)R ( (cid:126)x ) R ( (cid:126)x + (cid:126)r ) (cid:105) only depends on the relative distance r ≡ | (cid:126)r | ). The limit r → in eq. (A11) defines thevariance σ of the comoving curvature perturbation σ ≡ (cid:104)R ( (cid:126)x ) R ( (cid:126)x ) (cid:105) = lim r → (cid:104)R ( (cid:126)x ) R ( (cid:126)x + (cid:126)r ) (cid:105) = (cid:90) ∞ dkk k π ∆ R ( k ) ≡ (cid:90) ∞ dkk P R ( k ) , (A12)where we defined the dimensionless power spectrum P R ( k ) ≡ ( k / π )∆ R ( k ) . This equation gives to the powerspectrum an intuitive statistical meaning; P R ( k ) represents the contribution to the variance of the field per unit loga-rithmic bin around the comoving wavenumber k . As stated above, we can also compute eq. (A12) by taking the vacuumexpectation value of the quantum operator ˆ R ( τ, (cid:126)x ) . We have (using eq. (A5)) lim k/aH (cid:28) (cid:104) ˆ R ( τ, (cid:126)x ) ˆ R ( τ, (cid:126)x ) (cid:105) = lim k/aH (cid:28) (cid:90) ∞ dkk k π |R k ( τ ) | = (cid:90) ∞ dkk k π |R k | , (A13)where the last step means that we are considering a sufficiently late time (the limit lim k/aH (cid:28) , for fixed k , representsa time-limit since aH depends on time) such that the mode with comoving wavenumber k is frozen to its constantvalue after horizon crossing (and the time-dependence drops in the last equality). Eq. (A13) gives an operative defi-nition of the power spectrum. For a given model of inflation, we can solve eq. (A8) (equivalently, eq. (A6)) for each k and take a “late-time limit” in the sense specified above (that is we evaluate R k at some late time after it freezes to aconstant value). The power spectrum is given by P R ( k ) ≡ ( k / π ) |R k | and fully specifies the two-point correlatorof the random field R ( (cid:126)x ) . Equivalently, eq. (A13) can be derived from the computation of (cid:104)R ( (cid:126)x ) R ( (cid:126)x ) (cid:105) by means of thedecomposition given in eq. (A9).To proceed further, we need to specify the structure of higher-order correlators. Under the assumption that therandom field R ( (cid:126)x ) is gaussian, however, the power spectrum is enough to fully reconstruct higher-order correlators.This is because the N -point correlation function either vanishes (for odd N ) or can be expressed in terms of the powerspectrum as a consequence of the Isserlis’ theorem (for even N ). The assumption that the statistics of the randomfield R ( (cid:126)x ) is gaussian is well-motivated in the context of the canonical picture of slow-roll inflation. This is because if0one tries to compute, on the quantum-side, the three-point correlator lim k/aH (cid:28) (cid:104) ˆ R ( τ, (cid:126)x ) ˆ R ( τ, (cid:126)x ) ˆ R ( τ, (cid:126)x ) (cid:105) the resultingexpression turns out to be suppressed—compared to the two-point correlator—by additional powers of the Hubbleparameters (cid:15) and η [58]. This means that non-gaussianities are not relevant during conventional slow-roll dynamicsduring which the Hubble parameters take O ( (cid:28) values (said differently, this means that any detection of sizablenon-gaussianities at CMB scales will rule out all single field slow-roll models of inflation).However, we are interested in a situation which deviates from standard slow-roll dynamics. We refer to ref. [28] for adetailed (both analytical and numerical) study, and we summarize here the main points.Standard slow-roll dynamics takes place at large field values. However, few e -folds before the end of inflation (whichends at the absolute minimum of the potential located at the origin) the inflaton field crosses an approximate stationaryinflection point. The inflaton field almost stops but it possesses just enough inertia to overcome the approximatestationary inflection point. During this part of the dynamics the Hubble parameter η transits from η (cid:39) (that is typicalof slow-roll) to a large positive value that is maintained for few e -folds until the field crosses the approximate stationaryinflection point. If η (cid:38) / (typically one has η (cid:38) ), the friction term in eq. (A8) becomes negative. This part of thedynamics characterized by the presence of negative friction is dubbed ultra slow-roll. During the negative frictionphase, the modes R k —more precisely, their modulus |R k | —change exponentially fast, and can be either enhanced orsuppressed depending on the specific value of k .After the end of ultra slow-roll, η transits to a phase during which it takes negative O (1) values. The friction termin eq. (A8) turns positive, and the modes R k —after being enhanced or suppressed by the negative friction phase—arenow free to freeze to their final constant value. Since P R ( k ) ≡ ( k / π ) |R k | , the negative friction phase that modifiesexponentially |R k | produces a distinctive peak in the power spectrum of curvature perturbations.This peculiar dynamics has important consequences as far as non-gaussianities are concerned. Because of the pres-ence of the negative friction phase, classicalization of the modes do not happens after their horizon crossing but isdelayed after the end of ultra slow-roll [59]. This means that the three-point correlator (cid:104) ˆ R ( τ, (cid:126)x ) ˆ R ( τ, (cid:126)x ) ˆ R ( τ, (cid:126)x ) (cid:105) hasto be evaluated after the end of ultra slow-roll. Crucially, after the end of ultra slow-roll η takes sizable negative O (1) values (while we expect (cid:15) (cid:28) ). Let us indicate this value with η (which is a negative number). This implies thatnon-gaussianities are no longer negligible since the expected slow-roll suppression is not valid anymore. The explicitcomputation of the three-point correlator in Fourier space gives, for a triad of comoving wavenumbers k , k , k , theso-called local bispectrum [18] B R ( k , k , k ) (cid:39) − η [∆ R ( k )∆ R ( k ) + ∆ R ( k )∆ R ( k ) + ∆ R ( k )∆ R ( k )] , (A14)where the three-point correlator is given by the Fourier transform lim end of USR (cid:104) ˆ R ( τ, (cid:126)x ) ˆ R ( τ, (cid:126)x ) ˆ R ( τ, (cid:126)x ) (cid:105) = (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 ) B R ( k , k , k ) , (A15)and where lim end of USR indicates explicitly that the three-point correlator has to be evaluated at some late time afterthe end of ultra slow-roll, when the modes R k finally set to their final constant value (so that the right-hand side ofeq. (A15) does not depend on time).The analysis carried out in refs. [18] shows that non-gaussianities in the presence of ultra slow-roll are expected to benon-negligible. Eqs. (A14, A15) represent the quantum side of the story, and eq. (A14) can be obtained by means of theso-called “ in - in ” formalism [58]. It is important to understand the implications from the point of view of the randomfield R ( (cid:126)x ) . Consider the non-gaussian random field R NG ( (cid:126)x ) defined by R NG ( (cid:126)x ) ≡ R G ( (cid:126)x ) + ( − η )2 (cid:2) R G ( (cid:126)x ) − (cid:104)R G ( (cid:126)x ) (cid:105) (cid:3) , (A16)where R G ( (cid:126)x ) is a gaussian random field—with variance σ (see eq. (A12))—which admits the decomposition givenin eq. (A9); it is a simple exercise to show that the three-point correlator (cid:104)R NG ( (cid:126)x ) R NG ( (cid:126)x ) R NG ( (cid:126)x ) (cid:105) has precisely theform given in eq. (A15). Notice that in eq. (A16) the presence of the constant piece (cid:104)R G ( (cid:126)x ) (cid:105) guarantees that the non-gaussian random field R NG ( (cid:126)x ) has zero mean. This property is physically motivated by the fact that the backgroundsolution is stable.The physics-case sketched in this appendix motivates the non-gaussianities studied in this paper which are of theform given in eq. (A16). To avoid cluttering the notation, in the main part of this work we indicate simply with R thegaussian random field R G in eq. (A16) and with h the non-gaussian one R NG . Furthermore, we set α ≡ ( − η ) / .Before proceeding, an important comment is in order. As we have discussed, the structure of the non-gaussian ran-dom field given in eq. (A16) is motivated by the explicit computation of the three-point correlator in eqs. (A14, A15) inthe presence of ultra slow-roll. This computation is based on the “ in - in ” formalism in which one expands the interac-tion Hamiltonian up to the cubic order.However, computing only the bispectrum is not the end of the story. In principle, one should compute also thetrispectrum (that is the connected part of the four-point correlator) and check that its contribution does not alter the1form of non-gaussianities derived including only cubic interactions. Needless to say, the computation of the trispec-trum is anything but simple, and we are not aware of explicit results in the context of ultra slow-roll. From a morepragmatic phenomenological perspective, one possible way to proceed—widely used for the analysis of CMB non-Gaussianity, see ref. [60]—is the following. Instead of eq. (A16), one takes the more general ansatz R NG ( (cid:126)x ) ≡ R G ( (cid:126)x ) + f (cid:2) R G ( (cid:126)x ) − (cid:104)R G ( (cid:126)x ) (cid:105) (cid:3) + f R G ( (cid:126)x ) , (A17)where f , are free coefficients that parametrize quadratic and cubic deviations from the gaussian limit. Notice thateq. (A17) generalizes eq. (A16) in the sense that it introduces cubic corrections but it preserves locality (in the sense thatdeviations from exact gaussianity at (cid:126)x are located at the same spatial position). As discussed before, in the absence ofan explicit computation there is no guarantee that ultra slow-roll generates deviations from eq. (A16) that have theform given in eq. (A17). Nevertheless, eq. (A17) can be considered as a phenomenological parametrization to studydeviations from from eq. (A16), as done for instance in ref. [61]. For this reason, in our analysis we will try to set theformalism considering a generic deviation from gaussianity (but always assuming locality) even though we will presentour results for the motivated case of quadratic non-gaussianities given in eq. (A16). Appendix B: CliffsNotes on gaussian peak theory
As a warm-up, consider the case of a n -dimensional scalar gaussian random field R ( (cid:126)x ) that we identify with therandom field associated to curvature perturbations. Peak theory in the case of a scalar gaussian random field is well-known [21]. However, in this appendix we will give a detailed discussion. The reason is that in our approach there is animportant conceptual difference compared to the standard results of ref. [21]. Ref. [21] computes the number densityof peaks of the overdensity field working directly with δ ( (cid:126)x, t ) , and without any reference to curvature perturbations. Inthis paper, on the contrary, we aim to compute the same quantity but starting from the distribution of local maxima ofthe curvature perturbation. Following this alternative route, we will be able to find (see appendix D) a generalizationthat accounts for the case in which local non-gaussianities in the definition of R ( (cid:126)x ) are present.Let us start from basics. A n -dimensional scalar random field R ( (cid:126)x ) is a set of random variables, one for each point (cid:126)x in the n -dimensional real space, equipped with a probability distribution p [ R ( (cid:126)x ) , . . . , R ( (cid:126)x m )] d R ( (cid:126)x ) . . . d R ( (cid:126)x m ) which measures the probability that the function R has values in the range R ( (cid:126)x j ) to R ( (cid:126)x j ) + d R ( (cid:126)x j ) for each of the j = 1 , . . . , m , with m an arbitrary integer and (cid:126)x , . . . , (cid:126)x m arbitrary points. We are interested in the behavior of the random field for a point in space that is stationary, and we consider m = 1 with (cid:126)x = (cid:126)x st . We can expand around this point according to R ( (cid:126)x ) = R ( (cid:126)x st ) + 12 n (cid:88) i,j =1 R ij ( (cid:126)x st )( (cid:126)x − (cid:126)x st ) i ( (cid:126)x − (cid:126)x st ) j , (B1)from which we have R i ( (cid:126)x ) = n (cid:88) j =1 R ij ( (cid:126)x st )( (cid:126)x − (cid:126)x st ) j . (B2)The goal is to obtain the number density of these stationary points in the n -dimensional space. Not to violate thecosmological principle, we only want to consider random fields which are statistically homogeneous and isotropic.Consequently, the specific value of (cid:126)x st is irrelevant, and we can always shift to (cid:126)x st = (cid:126) . Let us, therefore, drop theexplicit dependence on (cid:126)x st .The quantity of central interest for the computation of the number density of stationary points of R is the joint prob-ability density distribution of the field R , its first and second derivatives. This is intuitively obvious, since identifyingmaxima (or minima) implies a set of conditions on field derivatives, and it is thus mandatory to know what is theirprobability distribution (derivatives of a random field are also random variables themselves).Consider the realistic case with n = 3 . We indicate with P ( R , R i , R ij ) d R d R i d R ij the joint probability distribu-tion for the field being in the range R to R + d R , the field gradient being in the range R i to R i + d R i and the secondderivative matrix elements being in the range R ij to R ij + d R ij , all at the same point in space. P ( R , R i , R ij ) is thejoint ten-dimensional probability density distribution.If the point is stationary, we can write the joint probability distribution as P ( R , R i = 0 , R ij ) | det( R ij ) | d R d (cid:126)xd R ij .We set R i = 0 since the point is stationary, and we used eq. (B1) to change variables in the gradient volume element. Notice that, as in ref. [21], all spatial separations and length scales are described in comoving coordinates in the cosmological background. Thismeans that the number density that we shall compute at the end of this section in eq. (B33) must be understood as a comoving number density. d (cid:126)x with the field being inthe range R to R + d R as n st ( R ) d R d (cid:126)x ≡ d R d (cid:126)x (cid:90) P ( R , R i = 0 , R ij ) | det( R ij ) | d R ij (cid:124) (cid:123)(cid:122) (cid:125) ≡ n st ( R ) , (B3)where the integral is extended to the whole range of variability of the second derivatives since we are consideringgeneric stationary points. The probability distribution to have a maximum in a volume d (cid:126)x with the field being inthe range R to R + d R is n max ( R ) d R d (cid:126)x ≡ d R d (cid:126)x (cid:90) max P ( R , R i = 0 , R ij ) | det( R ij ) | d R ij (cid:124) (cid:123)(cid:122) (cid:125) ≡ n max ( R ) , (B4)where now the integral is subject to the conditions on the Hessian matrix that define a maximum. An equivalent defi-nition holds in the case of a minimum. The probability density distribution n max ( R ) represents the number density oflocal maxima (where “number density” is defined in a probabilistic sense) with field value in the range R to R + d R .In order to extract quantitative informations, we need to compute P ( R , R i , R ij ) , set R i = 0 , and integrate. Thisstrategy does not depend on the specific statistics of the random field.The computation of P ( R , R i , R ij ) drastically simplifies in the gaussian case. This is because the joint probabilitydistribution of a gaussian field and its derivatives is a multivariate normal distribution. Consider the simplified casewith n = 2 in which we have more control on analytical formulas; we have six random variables that we collect in thecolumn vector R ≡ ( R , R x , R y , R xx , R xy , R yy ) T . We have P ( R , R i , R ij ) = 1(2 π ) k/ √ det C exp (cid:20) −
12 ( R − (cid:104) R (cid:105) ) T ( C − )( R − (cid:104) R (cid:105) ) (cid:21) , (B5)where k = 6 is the dimension of R , (cid:104) R (cid:105) is the column vector of the expectation values of R and C is the covariancematrix defined by C ≡ (cid:104) ( R − (cid:104) R (cid:105) )( R − (cid:104) R (cid:105) ) T (cid:105) with elements C ij = (cid:104) ( R − (cid:104) R (cid:105) ) i ( R − (cid:104) R (cid:105) ) j (cid:105) = (cid:104) R i R j (cid:105) − (cid:104) R i (cid:105)(cid:104) R j (cid:105) . (B6)From the computation of the covariance matrix one can fully reconstruct the joint probability distribution.We restrict the analysis to zero-mean random fields since this assumption is physically motivated (see discussion inappendix A). From the explicit computation of the two-point correlators (cid:104) R i R j (cid:105) , one finds that the covariance matrixtakes the form C = R R x R y R xx R xy R yy R σ − σ / − σ / R x σ / R y σ / R xx − σ / σ / σ / R xy σ / R yy − σ / σ / σ / , (B7)and we find that R x , R y and R xy are completely uncorrelated while R , R xx and R yy are correlated. We introduce thespectral moments σ j ≡ (cid:90) dkk P R ( k ) k j , (B8)where P R ( k ) is the dimensionless power spectrum of R . Notice that in two spatial dimensions the dimensionlesspower spectrum of R is related to the power spectrum by means of P R ( k ) = ( k / π )∆ R ( k ) with ∆ R ( k ) = |R k | . The simplest possibility would be n = 1 . However, the case n = 2 is the simplest setup in which a non-trivial discussion about spatial isotropy ispossible. (cid:104)R xx R yy (cid:105) . We use the explicit form of R given in eq. (A9).We find (cid:104)R xx R yy (cid:105) = (cid:90) d (cid:126)k (2 π ) d (cid:126)k (2 π ) (cid:104) (cid:16) − k ,x R k a (cid:126)k e i(cid:126)k · (cid:126)x − k ,x R ∗ k a † (cid:126)k e − i(cid:126)k · (cid:126)x (cid:17) (cid:16) − k ,y R k a (cid:126)k e i(cid:126)k · (cid:126)x − k ,y R ∗ k a † (cid:126)k e − i(cid:126)k · (cid:126)x (cid:17) (cid:105) = (cid:90) d (cid:126)k (2 π ) d (cid:126)k (2 π ) k ,x k ,y (cid:104) e i ( (cid:126)k − (cid:126)k ) · (cid:126)x R k R ∗ k (cid:104) a (cid:126)k a † (cid:126)k (cid:105) + e − i ( (cid:126)k − (cid:126)k ) · (cid:126)x R k R ∗ k (cid:104) a † (cid:126)k a (cid:126)k (cid:105) (cid:105) = (cid:90) d (cid:126)k (2 π ) k x k y |R k | = 14 π (cid:90) dkdϕk cos ϕ sin ϕ ∆ R ( k ) = 18 (cid:90) dkk k P R ( k ) , (B9)where in the last line we just introduced polar coordinates. All the entries in eq. (B7) can be computed in a similar way.Interestingly, the pattern of zeros in eq. (B7) and the relations among different non-zero entries—obtained before bymeans of a direct computation—can be understood as a consequence of homogeneity and isotropy.Homogeneity, that is translational invariance, implies that correlators do not depend on the specific spatial positionat which they are computed. For instance, this means that the spatial derivative of (cid:104)RR(cid:105) should vanish (this mustbe true for a generic correlator evaluated at a given spatial point); from this condition, one finds ∂ i ( (cid:104)RR(cid:105) ) = 0 →(cid:104)RR i (cid:105) = 0 so that R , R x and R y are uncorrelated. Similarly, from ∂ x ( (cid:104)R x R x (cid:105) ) = 0 one gets (cid:104)R x R xx (cid:105) = 0 (withsimilar relations along other directions) so that first and second derivatives are uncorrelated. On the contrary, from ∂ x ( (cid:104)RR x (cid:105) ) = 0 it follows that (cid:104)RR xx (cid:105) = −(cid:104)R x R x (cid:105) as indeed obtained in eq. (B7). Similarly, (cid:104)RR yy (cid:105) = −(cid:104)R y R y (cid:105) .Isotropy, that is rotational invariance, implies that correlators do not depend on a particular direction in space. Thesimplest consequence of isotropy is that (cid:104)R x R x (cid:105) = (cid:104)R y R y (cid:105) and (cid:104)R x R y (cid:105) = 0 . In order to derive these two conditions,a very useful trick (that we shall use also in the non-gaussian computation) is to introduce—instead of the two compo-nents x and y of the two-dimensional vector (cid:126)x —the complex conjugated variables z ≡ x + iy and z ∗ = x − iy from whichwe have ∂ z = ( ∂ x − i∂ y ) / and ∂ z ∗ = ( ∂ x + i∂ y ) / . If we now rotate the two-dimensional vector (cid:126)x (without changingits length) the complex number z changes by a phase factor e iτ , that is z → e iτ z (and z ∗ → e − iτ z ∗ ). Consequently,the derivatives with respect to z and z ∗ change according to ∂ z → e − iτ ∂ z and ∂ z ∗ → e iτ ∂ z ∗ . If we now consider thecorrelators (cid:104)R z R z (cid:105) and (cid:104)R z ∗ R z ∗ (cid:105) they rotate according to (cid:104)R z R z (cid:105) → e − iτ (cid:104)R z R z (cid:105) and (cid:104)R z ∗ R z ∗ (cid:105) → e iτ (cid:104)R z ∗ R z ∗ (cid:105) .Because of isotropy of the two-dimensional space, (cid:104)R z R z (cid:105) and (cid:104)R z ∗ R z ∗ (cid:105) can not depend on τ and, therefore, theymust vanish. Consequently, the system of equations (cid:104)R z R z (cid:105) = 14 ( (cid:104)R x R x (cid:105) − i (cid:104)R x R y (cid:105) − (cid:104)R y R y (cid:105) ) = 0 , (B10) (cid:104)R z ∗ R z ∗ (cid:105) = 14 ( (cid:104)R x R x (cid:105) + 2 i (cid:104)R x R y (cid:105) − (cid:104)R y R y (cid:105) ) = 0 , (B11)admits the solution (cid:104)R x R x (cid:105) = (cid:104)R y R y (cid:105) and (cid:104)R x R y (cid:105) = 0 which are precisely the relations we were looking for. Morein general, if we consider the rotation ∂ z → e − iτ ∂ z and ∂ z ∗ → e iτ ∂ z ∗ a generic correlator will take a phase factor e iκτ where κ ≡ ( z ∗ derivatives) − ( z derivatives) . If κ (cid:54) = 0 , then the correlator must be equal to zero as a consequenceof isotropy. For instance, we have (cid:104)R zz R z ∗ (cid:105) = 0 (because κ = − ) but (cid:104)R z R z ∗ (cid:105) (cid:54) = 0 (because κ = 0 ).Combining homogeneity and isotropy validates the remaining entries in eq. (B7). For instance, from ∂ y ( (cid:104)RR x (cid:105) ) = (cid:104)R x R y (cid:105) + (cid:104)RR xy (cid:105) = 0 (homogeneity) we find (cid:104)RR xy (cid:105) = 0 since isotropy implies (cid:104)R x R y (cid:105) = 0 . As a final check, weconsider the block of the second derivatives in eq. (B7). From the previous argument, isotropy implies that (cid:104)R zz R zz (cid:105) =0 and (cid:104)R zz ∗ R zz (cid:105) = 0 (together with their complex conjugated (cid:104)R z ∗ z ∗ R z ∗ z ∗ (cid:105) = 0 and (cid:104)R zz ∗ R z ∗ z ∗ (cid:105) = 0 ). These tworelations imply (cid:104)R xx R xy (cid:105) = (cid:104)R yy R xy (cid:105) = 0 and (cid:104)R xx R xx (cid:105) = (cid:104)R xx R yy (cid:105) + 2 (cid:104)R xy R xy (cid:105) . Both these relations areverified by the entries in eq. (B7). We can actually do more since it is possible to show that (cid:104)R xx R yy (cid:105) = (cid:104)R xy R xy (cid:105) . Letus write (cid:104)R xx R yy (cid:105) = (cid:104) ∂ xx R ( (cid:126)x ) ∂ yy R ( (cid:126)x ) (cid:105) = (cid:104) ∂ x x R ( (cid:126)x ) ∂ y y R ( (cid:126)x ) (cid:105)| (cid:126)x = (cid:126)x = (cid:126)x = ∂ x x ∂ y y (cid:104)R ( (cid:126)x ) R ( (cid:126)x ) (cid:105)| (cid:126)x = (cid:126)x = (cid:126)x = ∂ x ∂ y (cid:104)R x ( (cid:126)x ) R y ( (cid:126)x ) (cid:105)| (cid:126)x = (cid:126)x = (cid:126)x . (B12)where the first step is just a more explicit definition of (cid:104)R xx R yy (cid:105) while in the following ones we consider two distinctpoint (cid:126)x = ( x , y ) and (cid:126)x = ( x , y ) that we later set equal again. Now the point is that because of homogeneity ofspace the correlator (cid:104)R x ( (cid:126)x ) R y ( (cid:126)x ) (cid:105) depends only on the distance | (cid:126)x − (cid:126)x | ; we can, therefore, exchange ↔ ob-taining (cid:104)R x ( (cid:126)x ) R y ( (cid:126)x ) (cid:105) without altering the result. From eq. (B12) this means that we have (cid:104)R xx R yy (cid:105) = (cid:104)R xy R xy (cid:105) as indeed verified in eq. (B7). If we combine this result with the previous relation (cid:104)R xx R xx (cid:105) = (cid:104)R xx R yy (cid:105) + 2 (cid:104)R xy R xy (cid:105) we find (cid:104)R xx R xx (cid:105) = 3 (cid:104)R xy R xy (cid:105) which is again verified in eq. (B7). Finally, we also note that the condition (cid:104)R xx R yy (cid:105) = (cid:104)R xy R xy (cid:105) implies that (cid:104)R zz ∗ R zz ∗ (cid:105) = (cid:104)R zz R z ∗ z ∗ (cid:105) . These kind of relations based on homogeneity and isotropy will beuseful later in the non-gaussian case.Using the properties of the exponential function, eq. (B5) takes the form P ( R , R i , R ij ) = P ( R x ) P ( R y ) P ( R xy ) P ( R , R xx , R yy ) , (B13)4where P ( R x ) = 1 (cid:112) πσ exp (cid:18) − R x σ (cid:19) , P ( R y ) = 1 (cid:112) πσ exp (cid:32) − R y σ (cid:33) , P ( R xy ) = 2 (cid:112) πσ exp (cid:32) − R xy σ (cid:33) , (B14)and P ( R , R xx , R yy ) = 1(2 π ) / (cid:112) det ˜ C exp (cid:18) −
12 ˜ R T ˜ C − ˜ R (cid:19) , ˜ C ≡ R R xx R yy R σ − σ / − σ / R xx − σ / σ / σ / R yy − σ / σ / σ / , (B15)with ˜ R ≡ ( R , R xx , R yy ) T . As customary in the gaussian case, from the knowledge of the power spectrum it is possibleto fully reconstruct the statistics of the random field. From eq. (B4) we get the number density of maxima n max ( R ) = 1 πσ (cid:90) max d R xx d R xy d R yy (cid:12)(cid:12) R xx R yy − R xy (cid:12)(cid:12) P ( R xy ) P ( R , R xx , R yy ) , (B16)where we used P ( R x = 0) = P ( R y = 0) = 1 / (cid:112) πσ . The integration region is defined by the conditions max = {R xx R yy − R xy > ∧ R xx < ∧ R yy < } .The change of variables {R xx , R yy , R xy } → { r, s, θ } defined by r cos θ ≡
12 ( R xx − R yy ) , r sin θ ≡ R xy , s ≡ −
12 ( R xx + R yy ) , (B17)turns out the be useful. The Jacobian of the transformation is J = 2 r , and we have R xx R yy − R xy = s − r . Thecondition s − r > becomes − s < r < s with s > since R xx < and R yy < . Furthermore, we restrict to r > if < θ < π . All in all, we have max = { < θ < π ∧ s > ∧ < r < s } .The parametrization in terms of { r, s, θ } is useful because we have s = − ( R xx + R yy ) = −(cid:52)R . (B18)This implies that any condition that restricts the value of the curvature −(cid:52)R can be implemented through s by im-posing s > s min instead of s > . This is the case of eq. (6) with α = 0 (thus h = R ), which reads s min σ = 98 ( a m H m ) σ δ c , (B19)Eq. (B16) becomes n max ( R , s min ) = 8 π σ σ (cid:112) σ σ − σ (cid:90) ∞ s min ds (cid:90) s dr r ( s − r ) exp (cid:20) − ( R σ + 4 s σ − s R σ )2( σ σ − σ ) − r σ (cid:21) , (B20)where, according to the previous argument, we set to s min the lower limit of integration over s and define n max ( R , s min ) such that n max ( R , s min = 0) = n max ( R ) . The integration over r gives n max ( R , s min ) = (B21) σ π σ (cid:112) σ σ − σ (cid:90) ∞ s min ds (cid:18) s σ + e − s σ − (cid:19) exp (cid:26) − σ σ σ σ − σ ) (cid:20) s σ − sσ (cid:18) σ σ σ (cid:19) R σ + R σ (cid:21)(cid:27)(cid:124) (cid:123)(cid:122) (cid:125) ≡ (cid:82) ∞ s min ds ¯ n max ( R ,s ) , where ¯ n max ( R , s ) ≡ σ π σ (cid:112) σ σ − σ (cid:18) s σ + e − s σ − (cid:19) exp (cid:26) − σ σ σ σ − σ ) (cid:20) s σ − sσ (cid:18) σ σ σ (cid:19) R σ + R σ (cid:21)(cid:27) , (B22)can be interpreted as the number density of maxima with field value in the range R to R + d R and curvature −(cid:52)R inthe range s to s + ds ) . We note that the argument of the exponential function in ¯ n max ( R , s ) is invariant under theexchange R /σ ↔ s/σ .5From the structure of the argument of the exponential function in ¯ n max ( R , s ) , it is natural to introduce the dimen-sionless parameter γ ≡ σ /σ σ which is completely determined, as we shall explain in a moment, by the properties ofthe power spectrum. Furthermore, it is easy to see that we have < γ < . In turn, this condition implies σ σ − σ > so that the square root in eq. (B22) is always real valued. Let us verify the non-trivial condition σ σ − σ > ; from thedefinition in eq. (B8), we have σ σ − σ = (cid:20)(cid:90) dkk P R ( k ) (cid:21) (cid:20)(cid:90) dk (cid:48) k (cid:48) P R ( k (cid:48) ) k (cid:48) (cid:21) − (cid:20)(cid:90) dkk P R ( k ) k (cid:21) = (cid:90) dkk dk (cid:48) k (cid:48) P R ( k ) P R ( k (cid:48) ) (cid:0) k (cid:48) − k (cid:48) k (cid:1) . (B23)To conclude that σ σ − σ > , we only need to show that ( k (cid:48) − k (cid:48) k ) > since the power spectrum is positivedefinite and the integrals over k cover the positive real axis. In the first double-integral, we can change k ↔ k (cid:48) withoutaltering the result, and this means that we can also substitute k (cid:48) → ( k + k (cid:48) ) / ; if we do this transformation, thefactor ( k (cid:48) − k (cid:48) k ) becomes ( k (cid:48) − k ) / which is always positive. This concludes the proof.From eq. (B22) we see that the value of γ controls the amount of correlation between R and s . If we take γ = 0 , thetwo variables are completely uncorrelated (because σ = 0 ). On the contrary, γ = 1 corresponds to the case in whichthey are maximally correlated.To proceed further, we assume the following analytical expression for the power spectrum P R ( k ) = A g √ πv exp (cid:20) − log ( k/k (cid:63) )2 v (cid:21) . (B24)As far as the spectral moments in eq. (B8) are concerned, we find the analytical result σ j = A g k j(cid:63) e j v , (B25)which implies σ = A g and γ = e − v , where we see that < γ < as expected. The three parameters { A g , v, k (cid:63) } fully specify our problem. The physical picture is the following. ◦ The scale k (cid:63) represents the comoving wavenumber at which the power spectrum peaks. This quantity is relatedto the mass of the black holes produced after the collapse of the regions where the overdensity field is abovethreshold. We take k (cid:63) = O (10 ) Mpc − corresponding to M PBH = O (10 ) g. ◦ The parameter v controls the broadness of the power spectrum and, in turn, the broadness of the mass distribu-tion of the black holes. The case v = 0 . , for instance, corresponds to a very narrow power spectrum. This, inturn, will generate a very narrow mass distribution of black holes. ◦ The amplitude of the power spectrum A g is related to the abundance of dark matter in the present-day Universein the form of black holes, and typical values are of order A g / √ πv = O (10 − ) . To fix ideas, to get P R ( k (cid:63) ) = 10 − one needs A g = 2 . × − for v = 0 . .Let us now pause for a moment to clarify the rationale of the computation that we are doing. We are interested in regionsof space where the field R has large curvature −(cid:52)R since in this case the overdensity field (which is proportional to −(cid:52)R ) takes large values. In eq. (B22) we can select regions with large curvature by implementing (as done in eq. (B21))a lower limit of integration over s . However, eq. (B22) always associates, by construction, regions with large curvature(consequently, peaks of the overdensity field) with local maxima of R . This association can be analytically justified asfollows. Consider the Taylor expansion in eq. (B1) which we rewrite in two spatial dimensions taking a local maximumas stationary point of R R ( (cid:126)x ) = R M + 12 (cid:88) i,j =1 R ij ( (cid:126)x M )( (cid:126)x − (cid:126)x M ) i ( (cid:126)x − (cid:126)x M ) j . (B26)This quadratic equation can be written in a canonical form (that is without cross terms) if we do a coordinate transfor-mation (cid:126)x → (cid:126)x (cid:48) ≡ ( x (cid:48) , y (cid:48) ) such that the new axes are aligned along the eigenvectors of the matrix R ij ( (cid:126)x M ) . We are freeto do this rotation because of isotropy. In such case, eq. (B26) takes the form R ( (cid:126)x (cid:48) ) = R M −
12 ( λ x (cid:48) + λ y (cid:48) ) , (B27)where λ i =1 , > are the (minus) eigenvalues of R ij ( (cid:126)x M ) . Notice that in eq. (B27) we used homogeneity to shift theposition of the maximum (cid:126)x M to the origin of the new coordinate system. From eq. (B27) we see that an iso-densitysurface with constant R (cid:126)x (cid:48) ≡ R ( (cid:126)x (cid:48) ) is an ellipse with canonical equation λ R M − R (cid:126)x (cid:48) ) x (cid:48) + λ R M − R (cid:126)x (cid:48) ) y (cid:48) = 1 , a i ≡ (cid:20) R M − R (cid:126)x (cid:48) ) λ i (cid:21) / , (B28)6and semi-axes a i =1 , . The key point is that the actual magnitude of the eigenvalues λ i depends on the steepness of thefield R around the position of its maximum. This follows from eq. (B27) if we apply the Laplacian with respect to thecoordinates (cid:126)x (cid:48) since we find λ + λ = −(cid:52)R ( (cid:126)x (cid:48) ) . (B29)Suppose now that the point (cid:126)x (cid:48) coincides with a peak of the overdensity field, (cid:126)y pk . Using eq. (6), we find λ + λ = −(cid:52)R ( (cid:126)y pk ) (cid:38)
94 ( aH ) δ c = 2 σ
98 ( aH ) σ δ c (cid:124) (cid:123)(cid:122) (cid:125) (cid:29) = ⇒ λ i σ (cid:29) . (B30)In eq. (B30) we used the fact that we typically expect ( aH ) /σ (cid:29) ; we will comment in more detail about this estimateat the end of this section. Notice that in eq. (B30) we assumed that the two eigenvalues λ i have the same magnitude.This is because we have separately λ = −R xx ( (cid:126)y pk ) and λ = −R yy ( (cid:126)y pk ) , and the second derivatives R xx ad R yy have the same covariance (see eq. (B7)). We can now do the same expansion in eq. (B26) but with respect to the firstderivatives of R at (cid:126)x M . The stationary condition R i ( (cid:126)x M ) = 0 reads R i ( (cid:126)x ) − (cid:88) j =1 R ij ( (cid:126)x M )( (cid:126)x − (cid:126)x M ) j = 0 = ⇒ R x ( (cid:126)x (cid:48) ) + λ x (cid:48) = 0 R y ( (cid:126)x (cid:48) ) + λ y (cid:48) = 0 (B31)where in the last line we introduced, as done before, the eigenvalues λ i =1 , . We identify again the point (cid:126)x (cid:48) with a peakof the overdensity field so that we can use the estimate in eq. (B30). Furthermore, (cid:126)x (cid:48) represents now, by construction,the distance between the local maximum of R and the peak of the overdensity field. We can estimate this distance bymeans of eq. (B30) and eq. (B31). In eq. (B31), R x ( (cid:126)x (cid:48) ) and R y ( (cid:126)x (cid:48) ) are not equal to zero (since we moved away fromthe local maximum of R ) and their magnitude can be estimated (in a probabilistic sense) by means of the covariancematrix. We have R x ( (cid:126)x (cid:48) ) ≈ R y ( (cid:126)x (cid:48) ) ≈ (cid:104)R x R x (cid:105) / ≈ σ . We find | (cid:126)x (cid:48) | ≈ σ σ (cid:18) λ i /σ (cid:19) (cid:28) σ σ = e − v k (cid:63) . (B32)As we will see in a moment (see also appendix F), we have that /k (cid:63) is typically of the order of the comoving horizonlength /aH at the time when the perturbations re-enter the horizon and become causally connected. Therefore, wefind | (cid:126)x (cid:48) | (cid:28) /aH . This means that high peaks of the overdensity field lie “close” (that is within an Hubble radius) tolocal maxima of the curvature perturbation. We can validate the analytical approximation by means of a numericalcheck. To this end, we use the full joint probability density distribution in eq. (B5) to generate a sample of randomvalues that we distribute on a two-dimensional grid (see caption of fig. 7). In other words, at each point on the spatialgrid corresponds a value of R = ( R , R x , R y , R xx , R xy , R yy ) T randomly generated from eq. (B5). In the left panel offig. 7 we show the spatial distribution of the random variable s/σ = −(cid:52)R /σ . We focus on a region in which −(cid:52)R takes a large value (in units of σ ). In the left panel of fig. 7, we indicate this region with a red contour. We zoom inthis part of the plot in the right panel of fig. 7. The analytical argument explained before suggests that we should find amaximum close to the point where the curvature field peaks. We look for this maximum numerically by looking at thebehavior of the gradient field {R x , R y } that we plot using blue arrows. We indeed find a local maximum that we markwith a yellow star. We checked numerically that at the position of the yellow star where the gradient field vanishes theconditions on the second derivatives that define a maximum are verified. We note that at the position of the peak of theoverdensity field the gradient field {R x , R y } does not vanish but we find that its magnitude is of order O (1) (in unitsof σ ) as we argued in eq. (B32). Furthermore, we checked that the local maximum lies closer to the peak of −(cid:52)R /σ for increasing higher values of the latter.This numerical result corroborates the validity of the analytical argument in eq. (B32), and we conclude, therefore,that eq. (B22) is the right distribution to consider: We count the peaks of δ by looking at the maxima of R with largecurvature.Next, we ask if there exists some relation between the curvature of a local maximum of R and the value of R at themaximum. We can use the probability density distribution ¯ n max ( R , s ) defined in eq. (B22) to generate numerically,in position space, a sample of maxima by extracting randomly the value of R and curvature s = −(cid:52)R . We can thenuse eq. (5) to extract from the distribution of s the distribution of the overdensity field (of course, with R instead of h The result of this exercise will be useful later to discuss the impact of non-linearities in eq. (5), see appendix G. ★★ ★★ ��� ��� ������������ � � FIG. 7:
Numerical simulation of the random field R together with its first and second derivatives. At each point in space (discretizedin steps ∆ x = 5 and ∆ y = 5 ) we associate a vector of values {R , R x , R y , R xx , R xy , R yy } randomly generated from eq. (B5). We usethe power spectrum in eq. (B24) to compute the correlation matrix. We set v = 0 . , k (cid:63) = 1 . × Mpc − and A g = 2 . × − .Left panel. We show the density plot of the random variable s/σ = −(cid:52)R /σ which is related to the overdensity field. Right panel.Same as in the left panel but zoomed in the region where the random variable −(cid:52)R /σ has a pronounced peak. We superimpose(blue arrows) a vector plot that keeps track of the gradient field {R x , R y } . The yellow star marks the position of the local maximumof R that lies close to the peak of −(cid:52)R /σ . since we are considering here the gaussian case). The outcome of this exercise is shown in fig. 8, fig. 9 and fig. 10 (seecaptions for details).Let us consider the case in which we take v = 0 . . We remind that this choice corresponds to a very narrow powerspectrum. We have γ (cid:39) . . As noticed before, in this case we expect a strong correlation between R and s . Thismeans that regions with large R /σ are likely to be also regions with large s/σ (and, consequently, regions where theoverdensity field peaks). This is evident in fig. 8, fig. 9 and fig. 10 (left panel). In particular, in fig. 9 we see that maximawith large values of R (left panel) maps precisely regions where the overdensity field δ peaks (right panel). As statedbefore, this is a consequence of the fact that R and s (hence δ ) are highly correlated. Furthermore, in the left panel offig. 10 the symmetry of the randomly generated points under the exchange R /σ ↔ s/σ is evident.It is instructive to consider what happens if we take a different value of v . Consider, for instance, the case with v = 0 . .We have γ (cid:39) . , and R and s are now much less correlated compared to the case with v = 0 . . This is shown in theright panel of fig. 10 (see caption for details). For v = 0 . , spiky maxima (that is maxima with large curvature) areseldom characterized also by a large value of R . This means that if we aim at deriving a generic formula for the numberdensity of spiky maxima we can not restrict eq. (B21) to special values of R .Let us summarize our findings so far. We have that the number density of local maxima of R with large curvaturegives the number density of regions where the overdensity field peaks. Physically, we are only interested in local maximawith large curvature irrespectively on their value of R . Since there is no special relation—for generic values of v (hence γ )—between s and R , eq. (B21) must be integrated over the entire range of variability of R .If we consider the case in which s min does not depend on R , we can integrate eq. (B21) analytically. We find thenumber density N max ( s min ) ≡ (cid:90) ∞−∞ d R (cid:90) ∞ s min ds ¯ n max ( R , s ) = σ √ π / σ (cid:34) s min σ exp (cid:18) − s σ (cid:19) + 12 (cid:114) π (cid:32) √ s min σ (cid:33)(cid:35) . (B33)We can define the dimensionful quantity R ∗ ≡ √ d σ /σ , where d is the number of spatial dimensions. From eq. (B25)we find R ∗ = √ σ σ = √ k (cid:63) e − v . (B34)8 FIG. 8:
Numerical simulation of maxima of R in two spatial dimensions. At each point ( x, y ) in space (discretized in steps ∆ x = 10 and ∆ y = 10 ), we extract randomly from the density distribution ¯ n max ( R , s ) (defined in eq. (B22)) the value of R and s = −(cid:52)R / .The former are shown in the left panel. We remark that, by construction, all point generated by means of ¯ n max ( R , s ) are maxima of R . As far as the values of s are concerned, we plot on the right panel the corresponding values of δ = (8 / /aH ) s . We use thepower spectrum in eq. (B24), and we set v = 0 . , k (cid:63) = 1 . × Mpc − and A g = 2 . × − . Furthermore, we take ( k (cid:63) /aH ) (cid:39) . as suggested by numerical simulations (see appendix F). For illustrative purposes, we set the threshold δ c = 0 . (which issignificantly smaller compared to the value expected from numerical simulation of gravitational collapse into black holes that is δ c (cid:39) . , see appendix F; we use a smaller value of δ c otherwise events over the threshold would be too rare to be simulated in our simplifiednumerical analysis). Points in the simulation with δ > δ c are marked with a black dot. FIG. 9:
Same as in fig. 8 but zoomed in the region delimited by dashed lines and displayed in the form of a tri-dimensional densityplot. We see that spiky maxima with large values of R (left panel) coincide with peaks of δ (right panel). In this simplified two-dimensional set-up, we can estimate the mass fraction of black holes by means of the dimen-sionless quantity (for more details, see appendix F) R ∗ N max ( s min ) = 1 √ π / (cid:34) s min σ exp (cid:18) − s σ (cid:19) + 12 (cid:114) π (cid:32) √ s min σ (cid:33)(cid:35) (cid:39) √ π / (cid:18) s min σ (cid:19) exp (cid:18) − s σ (cid:19) , (B35)where the last approximation is valid if s min /σ (cid:38) . The mass fraction is controlled by an exponential decaying9 � � � � ������������������������������������� ℛ / σ � - Δ ℛ / � σ � - � - � � � � � � ������������������������������������� ℛ / σ � - Δ ℛ / � σ � �� % ���� = ��� = ��� FIG. 10:
Left panel. Simulated points in figs. 8, 9 shown in the plane {R /σ , −(cid:52)R / σ } . The horizontal dashed line correspondsto the threshold value δ c (see caption of fig. 8). Points above threshold collapse into black holes (marked with black dots here andin figs. 8, 9). This simulation corresponds to v = 0 . , and we see that R /σ and −(cid:52)R / σ are strongly correlated. Right panel.Comparison between two numerical simulations with v = 0 . (green) and v = 0 . (red). The two ellipses represent the 95% confidenceregions. function with argument s min σ = 98 ( aH ) σ δ c = 9 e − v A / g (cid:18) aHk (cid:63) (cid:19) δ c . (B36)The value of s min /σ , which is crucial for the determination of the correct order-of-magnitude of the mass fraction, de-pends on the properties of the power spectrum via the factor e − v /A / g k (cid:63) and the details of the gravitational collapsethat leads to black hole formation via the factor ( aH ) δ c . The comoving horizon length /aH depends on time. Thecomputation of the threshold for black hole production introduces the time t m that is defined by the time when thecurvature perturbations cross the horizon and become causally connected. This is the time at which eq. (B36) has to becomputed. The time t m defines implicitly, by means of the condition a ( t m ) H ( t m ) r m = 1 , the length scale r m . Thislength scale enters in the numerical evaluation of eq. (B36), and one typically gets ( a m H m /k (cid:63) ) δ c = O (1) . A precisenumerical evaluation of this factor is needed in order to set the size of the exponential suppression in eq. (B35). Wepostpone to appendix F a more detailed discussion about this point.It is important to check the validity of eq. (B35) against the standard result. In ref. [21], the number density of peaksof the overdensity field is controlled by the exponential function exp( − δ c / σ δ ) where σ δ refers to the variance of theoverdensity field. By means of eq. (5) (linearized, and with R instead of h ), one finds σ δ = (16 / /aH ) σ . Conse-quently, exp( − δ c / σ δ ) matches precisely the exponential function in eq. (B35). This is another indication that comput-ing the number density of peaks of the overdensity field via the number density of maxima of R that are spiky enoughleads to sensible results; importantly, this is in spite of the fact that we derived eq. (B35) in two spatial dimensions.After this long and detailed discussion about the gaussian case, we are ready to move to the more interesting situationin which local non-gaussianities are present. Actually, we are already in the position to make an interesting comment.Suppose that we compute the analogue of eq. (B22) with local non-gaussianities with of course now h and −(cid:52) h insteadof R and −(cid:52)R . It is crucial to answer the same question that we asked in the gaussian case: Is the number density ofmaxima of h with large curvature −(cid:52) h a good proxy for the number density of peaks of the overdensity field? The sameanalytical argument discussed in eqs. (B26-B32) can be repeated with h and −(cid:52) h instead of R and −(cid:52)R . The estimate More precisely, the process of black hole formation involves three different times. First of all, the number density of peaks that eventually formblack holes has to be calculated at some initial time when perturbations are still super-horizon. All the analysis done so far is based on this time(even though we do not specify it explicitly) since we are considering comoving curvature perturbations which are constant (see appendix A).Second, we have the time t m when the perturbations re-enter the horizon and become causally connected. Finally, we have the time t f at whichblack holes form. In general t f (cid:54) = t m , and from simulations of gravitational collapse in numerical relativity we have ( t f /t m ) / (cid:39) [25]; eq. (6)is usually evaluated at time t m while the factor ( t f /t m ) / (cid:39) can be included in the computation of the black hole abundance (see, e.g.,discussion in ref. [62]). λ i =1 , in eq. (B30) remains the same. In eq. (B32), the only difference is that local non-gaussianitiesalter the entries of the covariance matrix that we used to estimate the magnitude of the gradient field. However, weanticipate that we find (see eq. (D73)) (cid:104) h x h x (cid:105) = (cid:104) h y h y (cid:105) = σ (1 + 4 α σ ) / and (cid:104) h x h y (cid:105) = 0 ; we conclude that thenon-gaussian correction in this case is negligible since we have α σ = α A g (cid:28) . Remember indeed that the size of σ is controlled by the amplitude of the power spectrum which is σ = A g (cid:28) . This is particularly clear for the simplechoice of P R in eq. (B24) but remains true in general.We conclude that the same argument discussed in eqs. (B26-B32) is valid also in the presence of local non-gaussianities. In appendix C and appendix D we will, therefore, move to compute the number density of maxima of h with curvature above the threshold for black hole production that is needed to generalize eq. (B35) to the case α (cid:54) = 0 .We consider two different approaches. ∗ In appendix C we come back to the formulation of the problem in three spatial dimensions. We will derive anexpression analogue to eq. (B35) but valid in the case α (cid:54) = 0 and three spatial dimensions. We dub the resultof appendix C “exact” because no approximations will be used throughout the computation (apart from the lin-earization in eq. (5); the case in which also non-linearities are included will be discussed in appendix G). ∗ In appendix D we follow a different route. We will consider an expansion in cumulants around the gaussianprobability density distribution. To make this approach more transparent from the analytic point of view, we willwork again in two spatial dimensions.
Appendix C: Peak statistics with local non-gaussianities
In appendix B we have shown that the number density of peaks of the overdensity field can be approximated withgood accuracy with the number density of maxima of the comoving density perturbation which are spiky enough, i.e.with a Laplacian smaller than a threshold. The analysis has been performed assuming two spatial dimensions, to ahave a better control of the analytic expressions, and for ease of visualization of our simulations. Nevertheless, thesame conclusion is also valid in three dimensions, which is the case of physical interest that we are going to considerin this section. We shall now use these results to compute the number density of peaks of the overdensity field inpresence of local non-gaussianities. Our starting point is the analogous of eq. (B3) for the non-gaussian random field h ( (cid:126)x ) . Furthermore, as explained above, we are interested into maxima of h which are spiky enough. Explicitly, we have n pk ( h ) dh d x = dh d x (cid:90) spiky max d h i d h ij P NG ( h, h i , h ij ) δ ( h i ) | det( h ij ) | , (C1)where P NG ( h, h i , h ij ) is the joint probability distribution function of the non-gaussian field h, the field gradient h i andthe second derivatives h ij . The Dirac delta δ ( h i ) enforces the condition that the point is stationary. We will show ina moment how to restrict the integration volume to the field configurations which represent spiky maxima of h . Thecomoving number density of peaks is obtained integrating over all the heights of the peaks N pk = (cid:90) ∞ h min dh n pk ( h ) = (cid:90) spiky max dh d h i d h ij P NG ( h, h i , h ij ) δ ( h i ) | det( h ij ) | . (C2)Notice that for α > , which is the case relevant for the PBH production (see sec. A), and from the relation h ( (cid:126)x ) = R ( (cid:126)x ) + α (cid:0) R ( (cid:126)x ) − σ (cid:1) , one realizes that h ( x ) attains a minimum h min = − (1 + 4 α σ ) / α for R = − / α. Toproceed, it turns out to be useful to consider the gaussian variables R , R i and R ij instead of the non-gaussian fields h , h i and h ij . The reason is that the joint probability distribution function of the former variables, that we denotewith P ( R , R i , R ij ) , is known and it has a simple analytic expression, see sec. B for the explicit formula in the caseof two dimensions. Using the conservation of the probability in a differential volume P NG ( h, h i , h ij ) dh d h i d h ij = P ( R , R i , R ij ) d R d R i d R ij , we can perform a change a variable and write N pk = (cid:90) spiky max d R d R i d R ij P ( R , R i , R ij ) δ [ h i ( R , R i , α )] | det [ h ij ( R , R i , R ij , α )] | . (C3)where the field gradient and the second derivatives of the non gaussian variables are written in terms of the gaussianfields as h i ( R , R i , α ) = R i (1 + 2 α R ) , h ij ( R , R i , R ij , α ) = R ij (1 + 2 αR ) + 2 α R i R j . (C4)As noticed in sec. II, stationary points of R are also stationary points of h. Moreover, from the expressions above, onecan see that also the configuration (1 + 2 αR ) = 0 is a stationary point of h. However, as mentioned above, it is a1minimum of h (for α > ), therefore we need to focus only on the configurations with R i = 0 . This means that onecan write eq. (C3) as N pk = (cid:90) spiky max d R d R ij P ( R , R i = 0 , R ij ) | det ( h ij ( R , R i = 0 , R ij , α )) | | α R| = (cid:90) spiky max d R d R ij P ( R , R i = 0 , R ij ) | det ( R ij ) | . (C5)In the equation above all the information about the non-gaussianities is confined in spiky max , i.e. in the condition toimpose on the integration volume in order to restrict on maxima of h spiky enough. Let us see how these constraintscan be implemented explicitly. The Hessian matrix, for R i = 0 , is h ij = R ij (1 + 2 α R ) . (C6)We should, therefore, select maxima (minima) of R for positive (negative) values of (1 + 2 α R ) . The calculation can besimplified aligning the coordinate axes along the eigenvectors of the matrix −R ij . The corresponding eigenvalues aredenoted as λ i with i = 1 , , (the same strategy has been adopted in sec. B in the case of two spatial dimensions). Thethree eigenvectors, and the three Euler angles defining their orientation, can be used to parametrize the six indepen-dent random fields of R ij . Let us also define the following variables ¯ ν = R /σ , σ x = −(cid:52)R = λ + λ + λ σ y = ( λ − λ ) / σ z = ( λ − λ + λ ) / (C7)Changing variables and integrating over the Euler angles, from eq. (C5) one obtains [21] N pk = (cid:90) spiky max ¯ n pk (¯ ν, x, y, z ) d ¯ ν dx dy dz with ¯ n pk (¯ ν, x, y, z ) = A e − Q | F ( x, y, z ) | , (C8)and A ≡ (cid:114) σ σ π (cid:112) − γ , (C9) Q ≡ ¯ ν x − x ∗ ) − γ ) + 52 (cid:0) y + z (cid:1) , (C10) F ( x, y, z ) ≡ y (cid:0) y − z (cid:1) (cid:104) ( x + z ) − y (cid:105) ( x − z ) , (C11) x ∗ ≡ γ ¯ ν . (C12)We can divide the integration on ¯ ν in two parts. ◦ α R > .This implies that maxima of R are also maxima of h , and that we are considering ¯ ν > − / σ α. One can choosean ordering for the eigenvalues of the matrix −R ij : λ ≥ λ ≥ λ . Therefore, we can select maxima of R requiringthat λ > . Under these conditions, and using eqs. (C7), the domain of integration for the variables y and z reads: (cid:90) x/ dy (cid:90) y − y dz + (cid:90) x/ x/ dy (cid:90) y y − x dz . (C13)Then, working within the linear approximation in eq. (5), eq. (6) implies that only spiky maxima should be se-lected x > a m H m ) σ δ c ασ ¯ ν ≡ x δ (¯ ν ) . (C14)In the equation above, the horizon scale /a H has been evaluated at the time t m when the curvature perturba-tions cross the horizon. Having specified the appropriate domain of integration for all the variables, we can nowintegrate eq. (C8) and multiply by a factor 6 to take into account all the possible orderings of the eigenvalues λ i .We find N (I)pk = (cid:90) ∞− ασ d ¯ ν (cid:90) ∞ x δ (¯ ν ) dx ¯ n pk (¯ ν, x ) , (C15)2where ¯ n pk (¯ ν, x ) = e − ¯ ν / (2 π ) R ∗ f ( x ) e − ( x − x ∗ )22(1 − γ (cid:112) π (1 − γ ) , (C16) f ( x ) = x − x (cid:34) Erf (cid:32)(cid:114) x (cid:33) + Erf (cid:32)(cid:114) x (cid:33)(cid:35) + (cid:114) π (cid:20)(cid:18) x (cid:19) e − x + (cid:18) x − (cid:19) e − x (cid:21) , (C17)and we remind that R ∗ ≡ √ σ /σ . ◦ α R < .In this case minima of R are maxima of h. Proceeding analogously as before, we have: N (II)pk = (cid:90) − ασ −∞ d ¯ ν (cid:90) x δ (¯ ν ) −∞ dx ¯ n pk (¯ ν, x ) . (C18)The comoving number density of peaks of the overdensity field is approximated as the sum of the two terms above: N pk = N (I)pk + N (II)pk . One can interpret these two contributions along the lines exposed in section II. The term N (II)pk counts then minima of R which are maxima of h. Instead N (I)pk corresponds to maxima of R , and the restriction in theintegration range of ¯ ν (the lower limit) is designed to subtract those maxima of R which are minima of h . Obviously,the gaussian result is recovered for α → . In this limit we can compare with the standard expression in ref. [21].Numerically, and for the power spectra under consideration, we found that the two calculations agree within a factor (cid:39) . As already mentioned in section B for the case of two spatial dimensions, this confirms that peaks of the overdensityfield are well approximated by peaks of the comoving density perturbation which are spiky enough. In appendix F wewill discuss how to translate N pk into the primordial black hole abundance in eq. (8). Appendix D: Peak theory with local non-gaussianities: a perturbative approach in 2D
The idea is to derive the joint probability density distribution for the variables h , h x , h y , h xx , h xy , h yy (consideringagain, for simplicity, the two-dimensional case). Once we get this probability density distribution, we set h x = h y =0 and we integrate over h xx , h xy and h yy in the domain defining maxima. This strategy was simple to implementin the case of the gaussian variable R because the joint probability density distribution was a multivariate normaldistribution. The non-gaussian case is more complicated.In order to tackle the problem, we shall use the approach based on the characteristic function. In full generality, fora set of N correlated random variables ξ i the characteristic function is the Fourier transform of their joint probabilitydensity distribution χ ( λ , . . . , λ N ) ≡ (cid:90) dξ . . . dξ N P ( ξ , . . . , ξ N ) exp [ i ( ξ λ + · · · + ξ N λ N )]= (cid:90) dξ . . . dξ N P ( ξ , . . . , ξ N ) (cid:20) i ( ξ λ + · · · + ξ N λ N ) + i
2! ( ξ λ + · · · + ξ N λ N ) + . . . (cid:21) = 1 + i (cid:88) j (cid:104) ξ j (cid:105) λ j + i (cid:88) j ,j (cid:104) ξ j ξ j (cid:105) λ j λ j + i (cid:88) j ,j ,j (cid:104) ξ j ξ j ξ j (cid:105) λ j λ j λ j + . . . (D1)where, after a Taylor expansion, we introduced the moments of the joint distribution by means of the integrals (cid:104) ξ j . . . ξ j k (cid:105) = (cid:90) dξ . . . dξ N P ( ξ , . . . , ξ N ) ξ j . . . ξ j k . (D2)If we take the natural log of the characteristic function and Taylor expand, we define the cumulants log χ ( λ , . . . , λ N ) ≡ i (cid:88) j C ( ξ j ) λ j + i (cid:88) j ,j C ( ξ j , ξ j ) λ j λ j + . . . . (D3)3The relation between moments and cumulants follows from the comparison of the two Taylor expansions. We find thewell-known relations C ( ξ j ) = (cid:104) ξ j (cid:105) , (D4) C ( ξ j , ξ j ) = (cid:104) ξ j ξ j (cid:105) − (cid:104) ξ j (cid:105)(cid:104) ξ j (cid:105) , (D5) C ( ξ j , ξ j , ξ j ) = (cid:104) ξ j ξ j ξ j (cid:105) − (cid:104) ξ j (cid:105)(cid:104) ξ j ξ j (cid:105) − (cid:104) ξ j (cid:105)(cid:104) ξ j ξ j (cid:105) − (cid:104) ξ j (cid:105)(cid:104) ξ j ξ j (cid:105) + 2 (cid:104) ξ j (cid:105)(cid:104) ξ j (cid:105)(cid:104) ξ j (cid:105) , (D6) C ( ξ j , ξ j , ξ j , ξ j ) = . . . (D7)and so on. Working with cumulants instead of moments is more efficient. This is particularly true in the gaussian case.For a set of correlated gaussian variables, all cumulants of order higher than two vanish. In the gaussian case, eq. (D3)gives χ ( λ , . . . , λ N ) = exp i (cid:88) j (cid:104) ξ j (cid:105) λ j − (cid:88) j ,j C ( ξ j , ξ j ) λ j λ j , (D8)and the inverse Fourier transform that gives the probability density distribution reads P ( ξ , . . . , ξ N ) = (cid:90) dλ (2 π ) . . . dλ N (2 π ) exp i (cid:88) j ( (cid:104) ξ j (cid:105) − ξ j ) λ j − (cid:88) j ,j C ( ξ j , ξ j ) λ j λ j . (D9)If we complete the square inside the integrand and compute the resulting multivariate gaussian integral, we find pre-cisely eq. (B5) where the second-order cumulants reconstruct the covariance matrix elements in eq. (B6). This wasprecisely the strategy that we followed in the previous section: we computed the elements of the covariance matrix(that are the second-order cumulants) and we (implicitly) performed an inverse Fourier transform to get back theprobability density distribution.In the non-gaussian case, we can try to apply the same logic. First, we compute the cumulants; second, we recon-struct the probability density distribution by means of an inverse Fourier transform.The computation of the cumulants require some mathematical tricks that we shall explain in the following. In fullgenerality, we need to compute the n th -order cumulant C n ( ∂ j h, . . . , ∂ l h ) where each ∂ j h represents a certain numberof spatial derivatives (zero, one or two) acting on h . Remember also that the random field h (and its derivatives) iscomputed at a specific spatial position (say, (cid:126)x ) that will be later identified with a stationary point. We can write C n ( ∂ j h, . . . , ∂ l h ) = C n [ ∂ j h ( (cid:126)x ) , . . . , ∂ l h ( (cid:126)x )]= C n [ ∂ j h ( (cid:126)x ) , . . . , ∂ l n h ( (cid:126)x n )] | (cid:126)x = ··· = (cid:126)x n = (cid:126)x = ∂ j . . . ∂ l n C n [ h ( (cid:126)x ) , . . . , h ( (cid:126)x n )] | (cid:126)x = ··· = (cid:126)x n = (cid:126)x . (D10)In the first step of eq. (D10) we consider each ∂ j k h to act at a different point (cid:126)x k , and later we set all points equal again (wealready used a similar trick in eq. (B12)). The advantage of this step is that since now each derivative acts at a differentspatial point, we can bring them outside the cumulant (as done in the last step of eq. (D10)). It is a very simple exerciseto check explicitly (for instance, by computing second-order cumulants of a gaussian variable with its derivatives) thatthis procedure is completely legitimate. Using our notation h ( (cid:126)x n ) = h n , the problem reduces to the computation of C n ( h , . . . , h n ) where now the random fields are evaluated at different spatial points. After computing C n ( h , . . . , h n ) ,we will take the spatial derivatives and finally set all points equal according to the prescription in eq. (D10).To compute C n ( h , . . . , h n ) , we consider the relation h = R + α R . We write in general h = R + f NL ( R ) where f NL ( R ) can be a generic non-linear function controlled by the parameter α . Let us consider the term α R as a pertur-bation and expand in α . We have C n ( h , . . . , h n ) = C n ( R , . . . , R n ) + C n [ f NL ( R ) , . . . , R n ] + · · · + C n [ R , . . . , f NL ( R n )] (D11) + { C n [ f NL ( R ) , f NL ( R ) , . . . , R n ] + . . . } + O ( α ) . (D12)This deconstruction makes clear the fact that we can compute the cumulants C n ( h , . . . , h n ) based on the cumulantscomputed for the gaussian random field R .The problem is now the computation of C n [ f NL ( R ) , . . . , R n ] where we remind that R n = R ( (cid:126)x n ) .In order to compute C n [ f NL ( R ) , . . . , R n ] , we need to work out an intermediate result. First, remember the defi-nition of random field that we gave at the beginning of appendix B: The scalar random field R ( (cid:126)x ) is a set of randomvariables, one for each point (cid:126)x in space, equipped with a joint probability density distribution p ( R , . . . , R n ) . Whenthe random field is gaussian, p ( R , . . . , R n ) is a multivariate Gaussian distribution and we can write p ( R , . . . , R n ) = 1(2 π ) n/ √ det σ exp − (cid:88) i,j ( σ − ) ij R i R j , (D13)4where σ ij = (cid:104)R i R j (cid:105) (as before, we consider zero-mean random field).Consider now a new set of random variables {Y , . . . , Y n } related to the previous one via the transformations Y = g ( R , . . . , R n ) , . . . , Y n = g n ( R , . . . , R n ) with inverse R = g − ( Y , . . . , Y n ) , . . . , R n = g − n ( Y , . . . , Y n ) . The jointprobability density distribution for the transformed variables is given by p Y ( Y , . . . , Y n ) = p [ g − ( Y , . . . , Y n ) , . . . , g − n ( Y , . . . , Y n )] | det J ( Y , . . . , Y n ) | , (D14)where the Jacobian matrix is J ( Y , . . . , Y n ) ≡ ∂ R ∂ Y . . . ∂ R ∂ Y n ... . . . ... ∂ R n ∂ Y . . . ∂ R n ∂ Y n , with ∂ R j ∂ Y k = ∂∂ Y k (cid:2) g − j ( Y , . . . , Y n ) (cid:3) . (D15)Eq. (D14) is known as Jabobi’s multivariate theorem.We note that in eq. (D14) we can use the relation | det J ( Y , . . . , Y n ) | = 1 / | det J ( R , . . . , R n ) | where J ( R , . . . , R n ) is the Jacobian matrix of the inverse transformation with elements ∂ Y j /∂ R k = ∂g j ( R , . . . , R n ) /∂ R k . This impliesthat we can rewrite p Y ( Y , . . . , Y n ) = (cid:90) d R . . . d R n p ( R , . . . , R n ) δ [ Y − g ( R , . . . , R n )] . . . δ [ Y n − g n ( R , . . . , R n )] , (D16)and eq. (D14) follows from eq. (D16) if one applies the transformation property of the multi-dimensional delta function δ [ Y − g ( R , . . . , R n )] . . . δ [ Y n − g n ( R , . . . , R n )] = δ [ R − g − ( Y , . . . , Y n )] . . . δ [ R n − g − n ( Y , . . . , Y n )] | det J ( R , . . . , R n ) | , (D17)and integrates over R i . Eq. (D16) is a very useful formula. Consider the characteristic function of p Y ( Y , . . . , Y n ) de-fined by the Fourier transform χ Y ( λ , . . . , λ n ) = (cid:90) d Y . . . d Y n p Y ( Y , . . . , Y n ) exp [ i ( Y λ + · · · + Y n λ n )] . (D18)Using eq. (D16) and integrating over Y i thanks to the delta functions, we find χ Y ( λ , . . . , λ n ) = (cid:90) d R . . . d R n e i [ g ( R ,..., R n ) λ + ··· + g n ( R , ... , R n ) λ n ] p ( R , . . . , R n ) . (D19)From the definition of cumulants with respect to the characteristic function χ Y ( λ , . . . , λ n ) , we finally find (seeeq. (D3)) C n ( Y , . . . , Y n ) = ( − i ) n ∂∂λ . . . ∂∂λ n log χ Y ( λ , . . . , λ n ) (cid:12)(cid:12)(cid:12)(cid:12) λ = ... = λ n =0 = ( − i ) n ∂∂λ . . . ∂∂λ n log (cid:90) d R . . . d R n e i [ g ( R ,..., R n ) λ + ··· + g n ( R ,..., R n ) λ n ] p ( R , . . . , R n ) (cid:12)(cid:12)(cid:12)(cid:12) λ = ... = λ n =0 . (D20)This is a generic formula that can be applied to our case.Consider the linear order in α . From eq. (D20), the cumulant C n [ f NL ( R ) , . . . , R n ] is given by C n [ f NL ( R ) , . . . , R n ] = ( − i ) n ∂∂λ . . . ∂∂λ n log (cid:90) d R . . . d R n e i [ f NL ( R ) λ + R λ + ··· + R n λ n ] p ( R , . . . , R n ) (cid:12)(cid:12)(cid:12)(cid:12) λ = ... = λ n =0 . (D21)A similar formula is applicable to all remaining first-order terms in eq. (D11). Similarly, at order O ( α ) one needs tocompute, for instance, cumulants like C n [ f NL ( R ) , f NL ( R ) , . . . , R n ] =( − i ) n ∂∂λ . . . ∂∂λ n log (cid:90) d R . . . d R n e i [ f NL ( R ) λ + f NL ( R ) λ + ··· + R n λ n ] p ( R , . . . , R n ) (cid:12)(cid:12)(cid:12)(cid:12) λ = ... = λ n =0 . (D22)5All we need to do is computing the integral and taking derivatives with respect to λ i . The computation of the integralsis simplified by the fact that the variables R i are gaussian with joint distribution given by eq. (D13).Before proceeding, an important remark is in order. For the sake of simplicity, we derived eq. (D14) assuming thatthe functions Y i =1 ,...,n = g i =1 ,...,n ( R , . . . , R n ) define one-to-one mappings. In this case, there exist unique inversefunctions so that R i =1 ,...,n = g − i =1 ,...,n ( Y , . . . , Y n ) . The case we are interested in, however, is not exactly of this form.If we solve h = R + α R for R , we indeed find R = ( − ± √ αh ) / α which has two roots.However, this is not an insurmountable problem since Jabobi’s multivariate theorem in eq. (D14) can be generalizedto the case in which the system Y i =1 ,...,n = g i =1 ,...,n ( R , . . . , R n ) admits at most a countable number of roots. Let usindicate these q = 1 , . . . , Q roots in the form R q,i =1 ,...,n = g − q,i =1 ,...,n ( Y , . . . , Y n ) . Eq. (D14) becomes p Y ( Y , . . . , Y n ) = Q (cid:88) q =1 p [ g − q, ( Y , . . . , Y n ) , . . . , g − q,n ( Y , . . . , Y n )] | det J q ( Y , . . . , Y n ) | , (D23)where now J q ( Y , . . . , Y n ) is the Jacobian corresponding to the q th root J q ( Y , . . . , Y n ) ≡ ∂ R q, ∂ Y . . . ∂ R q, ∂ Y n ... . . . ... ∂ R q,n ∂ Y . . . ∂ R q,n ∂ Y n , with ∂ R q,j ∂ Y k = ∂∂ Y k (cid:2) g − q,j ( Y , . . . , Y n ) (cid:3) . (D24)The multi-dimensional delta-function identity in eq. (D17) changes accordingly, and the final result in eq. (D20) re-mains unaltered.After this digression, we are ready to compute the integral in eq. (D20). Instead of looking for a generic expression,let us consider specific cases organized for increasing level of difficulty. ◦ The simplest possibility is to truncate the analysis at the linear order in α , eq. (D21). Consider the integral ineq. (D21) which we rewrite as I ( λ , . . . , λ n ) ≡ (cid:90) d R e if NL ( R ) λ (cid:20) (cid:90) d R . . . d R n e i [ R λ + ... + R n λ n ] p ( R , . . . , R n ) (cid:124) (cid:123)(cid:122) (cid:125) ≡ I (cid:48) ( R ,λ ,...,λ n ) (cid:21) . (D25)The key observation is that the integral inside the square brackets would precisely match the definition of thecharacteristic function χ ( λ , . . . , λ n ) of p ( R , . . . , R n ) if it were completed by an additional integration over R (see the definition in the first line of eq. (D1)). In such case, we could simply use (see eq. (D8)) (cid:90) d R e iλ R I (cid:48) ( R , λ , . . . , λ n ) = χ ( λ , . . . , λ n ) = exp (cid:20) − (cid:88) ij σ ij λ i λ j (cid:21) , (D26)since {R , . . . , R n } are gaussian. However, the fact that integration over R is missing implies that the inte-gral I (cid:48) ( R , λ , . . . , λ n ) is actually equal to the inverse Fourier transform of χ ( λ , . . . , λ n ) with respect to λ . Informulas, eq. (D25) becomes I ( λ , . . . , λ n ) = (cid:90) d R e if NL ( R ) λ (cid:90) dλ (2 π ) e − iλ R exp (cid:20) − (cid:88) ij σ ij λ i λ j (cid:21) = (cid:90) d R e if NL ( R ) λ exp (cid:20) − (cid:88) ij (cid:62) σ ij λ i λ j (cid:21) (cid:90) dλ (2 π ) exp (cid:20) − σ λ − (cid:18) i R + (cid:88) j (cid:62) σ j λ j (cid:19) λ (cid:21) = (cid:90) d R e if NL ( R ) λ √ πσ exp (cid:20) − (cid:88) ij (cid:62) σ ij λ i λ j + 12 σ (cid:18) i R + (cid:88) j (cid:62) σ j λ j (cid:19) (cid:21) , (D27)where we just reorganized terms before integrating over λ in the last step. If we take the natural log, we find log[ I ( λ , . . . , λ n )] = − (cid:88) ij (cid:62) σ ij λ i λ j + log (cid:90) d R e if NL ( R ) λ √ πσ exp (cid:20) σ (cid:18) i R + (cid:88) j (cid:62) σ j λ j (cid:19) (cid:21) . (D28)6According to eq. (D21), we now need to compute derivatives with respect to λ i =1 , ... ,n and finally set λ i =1 , ... ,n =0 . We find C n [ f NL ( R ) , . . . , R n ] = 1 √ πσ (cid:18) n (cid:89) k =2 σ k (cid:19) (cid:90) d R f NL ( R ) g n ( R n − ) e −R / σ , (D29)where g n ( R n − ) is a polynomial of order R n − whose explicit expression is given by g n ( R n − ) = ( − n − e R / σ d n − d R n − e −R / σ . (D31)This result implies that we can just integrate by parts n − times in eq. (D29). We find C n [ f NL ( R ) , . . . , R n ] = (cid:18) n (cid:89) k =2 σ k (cid:19) √ πσ (cid:90) d R f ( n − ( R ) e −R / σ (cid:124) (cid:123)(cid:122) (cid:125) = (cid:104) f ( n − ( R ) (cid:105) . (D32)The leftover integral is nothing but the expectation value of f ( n − ( R ) . All in all, we find (remember that σ ij = (cid:104)R i R j (cid:105) ) C n [ f NL ( R ) , . . . , R n ] = (cid:104) f ( n − ( R ) (cid:105)(cid:104)R R (cid:105) . . . (cid:104)R R n (cid:105) (D33) ◦ At order O ( α ) , consider eq. (D22) and the integral J ( λ , . . . , λ n ) ≡ (cid:90) d R d R e if NL ( R ) λ + if NL ( R ) λ (cid:20) (cid:90) d R . . . d R n e i [ R λ + ... + R n λ n ] p ( R , . . . , R n ) (cid:21) . (D34)The computation goes as before with the only difference that now we need to consider the inverse Fourier trans-form of the characteristic function with respect to both λ and λ . For the natural log of J ( λ , . . . , λ n ) we find log[ J ( λ , . . . , λ n )] = − (cid:88) ij (cid:62) σ ij λ i λ j (D35) + log (cid:90) d R d R e if NL ( R ) λ + if NL ( R ) λ π (cid:112) σ σ − σ exp (cid:20) σ σ − σ ) (cid:18) A σ + A σ − A A σ (cid:19)(cid:21) , where A ≡ i R + (cid:80) j (cid:62) σ j λ j and A ≡ i R + (cid:80) j (cid:62) σ j λ j . This equation is the analogue of eq. (D28) at order O ( α ) . As before, we need to to compute derivatives with respect to λ i =1 , ... ,n and finally set λ i =1 , ... ,n = 0 . Letus consider first the derivatives with respect to λ and λ . We find ( − i ) ∂∂λ ∂∂λ log[ J ( λ , . . . , λ n )] (cid:12)(cid:12)(cid:12)(cid:12) λ = λ =0 = (D36) + 12 π (cid:112) σ σ − σ (cid:124) (cid:123)(cid:122) (cid:125) π √ det ˜ σ (cid:90) d R d R f NL ( R ) f NL ( R ) exp (cid:20) σ σ − σ ) (cid:18) A σ + A σ − A A σ (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) = (1 / (cid:80) i,j =1 (˜ σ − ) ij A i A j (cid:21) − (cid:20) √ πσ (cid:90) d R f NL ( R ) e A / σ (cid:21)(cid:20) √ πσ (cid:90) d R f NL ( R ) e A / σ (cid:21) , where ˜ σ is the two-by-two sub-matrix of σ formed by the first two rows and columns. Before proceeding, let usconsider a simple example. If we are interested to the computation of the second-order cumulant, these two We can notice that the g n polynomials are nothing but the Hermite polynomials, except for a pre-factor. In particular: g n ( R n − ) = 1(2 σ ) n/ H n − (cid:18) R √ σ (cid:19) (D30) A j =1 , = i R j =1 , . We findthe simple result C [ f NL ( R ) , f NL ( R )] = (cid:104) f NL ( R ) f NL ( R ) (cid:105) − (cid:104) f NL ( R ) (cid:105)(cid:104) f NL ( R ) (cid:105) (D37)where the expectation values (cid:104) f NL ( R ) (cid:105) and (cid:104) f NL ( R ) (cid:105) are defined as in eq. (D32) while (cid:104) f NL ( R ) f NL ( R ) (cid:105) iscomputed by means of the joint probability distribution of R and R ; in formulas, we have (cid:104) f NL ( R ) f NL ( R ) (cid:105) = 12 π √ det ˜ σ (cid:90) d R d R f NL ( R ) f NL ( R ) exp (cid:20) − (cid:88) i,j =1 (˜ σ − ) ij R i R j (cid:21) . (D38)For n (cid:62) , we need to compute in eq. (D36) derivatives with respect to λ i =3 , ... ,n and finally set λ i =3 , ... ,n = 0 .These derivatives act on the exponential functions in the integrands of the two terms on the right-hand side ofeq. (D36). Let us start from the second term; if we compute the λ -derivatives on the exponential function, wefind the following result ( − i ) n − ∂∂λ n . . . ∂∂λ (cid:18) e A / σ e A / σ (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) λ = ... λ n =0 = (D39) ( − n − (cid:18) n (cid:89) k =3 σ k (cid:19)(cid:18) d n − d R n − e −R / σ (cid:19) e −R / σ + ( − n − (cid:18) n (cid:89) k =3 σ k (cid:19)(cid:18) d n − d R n − e −R / σ (cid:19) e −R / σ + ( − n − (cid:20) a + b = n − (cid:88) a,b> (cid:18) n (cid:89) k =3+ b σ k (cid:19)(cid:18) n − a (cid:89) j =3 σ j (cid:19)(cid:18) d a d R a e −R / σ (cid:19)(cid:18) d b d R b e −R / σ (cid:19) + perms of (3 , . . . , n ) (cid:21) . To fully understand the content of this equation, few comments are in order. Inside the square brackets, thesum (cid:80) a,b runs over all possible combinations of a > and b > such that a + b = n − . For in-stance, if we take n = 4 we have only { a = 1 , b = 1 } while for n = 5 we have two combinations, namely { a = 1 , b = 2 } and { a = 2 , b = 1 } (notice that if n = 3 , there are no combinations that are allowed, and theterm in the square brackets does not contribute to the third-order cumulant). Furthermore, after computingthe sum over a and b , a further sum over all permutations of the indices (3 , . . . , n ) that give distinct results isrequired. Let us clarify this point with an explicit example. Consider n = 4 . The sum over { a = 1 , b = 1 } gives ( e −R / σ e −R / σ R R /σ σ ) σ σ ; we now need include the term with (3 , exchanged so that the finalresult is ( e −R / σ e −R / σ R R /σ σ )( σ σ + σ σ ) . Only permutations that give distinct results needto be included. For instance, take for n = 5 the term proportional to σ σ σ ; in this case, the sum over permu-tations of (3 , , gives only three distinct terms, σ σ σ → σ σ σ + σ σ σ + σ σ σ (instead ofsix—that is the number of permutations of three objects—since the remaining three do not give distinct results).As far as the first term on the right-hand side of eq. (D36) is concerned, we find ( − i ) n − ∂∂λ n . . . ∂∂λ (cid:20) e A σ A σ − A A σ σ σ − σ (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) λ = ... λ n =0 = (D40) ( − n − (cid:18) n (cid:89) k =3 σ k (cid:19)(cid:20) d n − d R n − e − R σ R σ − R R σ σ σ − σ (cid:21) + ( − n − (cid:18) n (cid:89) k =3 σ k (cid:19)(cid:20) d n − d R n − e − R σ R σ − R R σ σ σ − σ (cid:21) + ( − n − (cid:26) a + b = n − (cid:88) a,b> (cid:18) n (cid:89) k =3+ b σ k (cid:19)(cid:18) n − a (cid:89) j =3 σ j (cid:19)(cid:20) ∂ n − ∂ R a ∂ R b e − R σ R σ − R R σ σ σ − σ (cid:21) + perms of (3 , . . . , n ) (cid:27) , where the sum (cid:80) a,b and the subsequent one over different permutations have the same meaning explainedbefore.Eqs. (D39, D40) allow to compute the integrals in eq. (D36) by means of n − integration by parts. We find thefinal result (for n (cid:62) )8 C n [ f NL ( R ) , f NL ( R ) , . . . , R n ] = (D41) (cid:104) (cid:104) f ( n − ( R ) f NL ( R ) (cid:105) − (cid:104) f ( n − ( R ) (cid:105)(cid:104) f NL ( R ) (cid:105) (cid:105) (cid:104)R R (cid:105) . . . (cid:104)R R n (cid:105) + (cid:104) (cid:104) f NL ( R ) f ( n − ( R ) (cid:105) − (cid:104) f NL ( R ) (cid:105)(cid:104) f ( n − ( R ) (cid:105) (cid:105) (cid:104)R R (cid:105) . . . (cid:104)R R n (cid:105) + (cid:26) a + b = n − (cid:88) a,b> (cid:104) (cid:104) f ( a )NL ( R ) f ( b )NL ( R ) (cid:105) − (cid:104) f ( a )NL ( R ) (cid:105)(cid:104) f ( b )NL ( R ) (cid:105) (cid:105) n (cid:89) k =3+ b (cid:104)R R k (cid:105) n − a (cid:89) j =3 (cid:104)R R j (cid:105) + perms of (3 , . . . , n ) (cid:27) which generalizes eq. (D33) at order O ( α ) . We remark that eqs. (D33, D41) are completely generic and do notdepend on the functional form of f NL ( R ) . ◦ One can in principle continue the computation and include higher-order α -corrections to the cumulants. Need-less to say, the corresponding expressions quickly become quite unwieldy (even though the computation remainsconceptually simple). One possible strategy is to include only O ( α ) corrections, and check that O ( α ) termsremain sub-leading. In order for this criterium to be satisfactory, one should correctly identify the expansionparameter. As we shall see, in the case with no derivatives involved in eq. (D10) the latter turns out to be α σ . ◦ We can conclude this part with few checks of eqs. (D33, D41). For simplicity, let us consider the case in whichthere are no derivatives in eq. (D10). n = 2 We compute C ( h , h ) . At order O ( α ) , we have C ( h , h ) = C ( R , R ) + C [ f NL ( R ) , R ] + C [ R , f NL ( R )] (D42) = (cid:104)R R (cid:105) + (cid:104) f (1)NL ( R ) (cid:105)(cid:104)R R (cid:105) + (cid:104) f (1)NL ( R ) (cid:105)(cid:104)R R (cid:105) , (D43)where we use eq. (D5) to rewrite the first term (which is the gaussian one) and eq. (D33) to rewrite the last twoones. Since we have f NL ( R ) = α R , we have (cid:104) f (1)NL ( R i =1 , ) (cid:105) = 0 since we are considering zero-mean randomfields. This means that there are no corrections at order O ( α ) . We now move to consider corrections at order O ( α ) . We use eq. (D37). At order O ( α ) , we have C ( h , h ) = (cid:104)R R (cid:105) + (cid:104) f NL ( R ) f NL ( R ) (cid:105) − (cid:104) f NL ( R ) (cid:105)(cid:104) f NL ( R ) (cid:105) = (cid:104)R R (cid:105) + α (cid:104)R R (cid:105) − α (cid:104)R (cid:105)(cid:104)R (cid:105) . (D44)We now set the two points equal (since we are not considering derivatives of cumulants, there is no need ofconsidering different spatial point). We have (cid:104)R (cid:105) = σ and (cid:104)R (cid:105) = 3 σ . In conclusion, we find C ( h, h ) = σ (cid:0) α σ (cid:1) . (D45)This is an exact result since corrections of order O ( α ) enter only at higher-orders. As a rule of thumb, for a cumu-lant C n of order n , only corrections up to order O ( α n ) are possible. Eq. (D45) shows that the correct expansionparameter is α σ . In realistic models of inflation we expect α = O (1) while we have σ = A g for the specificpower spectrum in eq. (B24). Since one typically has A g = O (10 − ) , we expect α σ (cid:28) . Let us also notice thatin applications of cosmological interest it is customary to use f NL ( R ) = α ( R − (cid:104)R (cid:105) ) since in this way one hasa zero-mean non-gaussian field h , (cid:104) h (cid:105) = 0 . It is simple to see that the constant shift in f NL ( R ) that is implied bythis particular choice does not affect the computation of cumulants of order equal or higher than two. However,it does change the first order cumulant since C ( h ) = (cid:104) h (cid:105) (see eq. (D4)). The choice f NL ( R ) = α ( R − (cid:104)R (cid:105) ) ,therefore, leads to C ( h ) = 0 . n = 3 Consider now C ( h , h , h ) . This cumulant vanishes in the gaussian limit. At order O ( α ) , we use eq. (D33) andcompute C ( h , h , h ) = C [ f NL ( R ) , R , R ] + C [ R , f NL ( R ) , R ] + C [ R , R , f NL ( R )] (D46) = 2 α ( (cid:104)R R (cid:105)(cid:104)R R (cid:105) + (cid:104)R R (cid:105)(cid:104)R R (cid:105) + (cid:104)R R (cid:105)(cid:104)R R (cid:105) ) . (D47)If we set all points equal, we find C ( h, h, h ) = 6 ασ . (D48)9We now move to consider the possible presence of a correction at order O ( α ) that can be computed by meansof eq. (D41). As already noticed, for n = 3 the sum (cid:80) a,b does not contribute to the final result. We are left withthe two terms (cid:104) (cid:104) f (1)NL ( R ) f NL ( R ) (cid:105) − (cid:104) f (1)NL ( R ) (cid:105)(cid:104) f NL ( R ) (cid:105) (cid:105) (cid:104)R R (cid:105) + (cid:104) (cid:104) f NL ( R ) f (1)NL ( R ) (cid:105) − (cid:104) f NL ( R ) (cid:105)(cid:104) f (1)NL ( R ) (cid:105) (cid:105) (cid:104)R R (cid:105) . (D49)If we set the three points equal and consider explicitly f NL ( R ) = α R , we find that eq. (D49) reduces to α (cid:104)R (cid:105) σ which is zero because (cid:104)R (cid:105) = 0 . Eq. (D48), however, is not an exact result since for n = 3 correctionsof order O ( α ) are possible. Mimicking eq. (D45), we expect C ( h, h, h ) = 6 ασ [1 + O ( α σ )] . (D50)The correction of order O ( α ) can not be computed by means of eq. (D41) but it is expected to be sub-leadingsince we work under the assumption that α σ (cid:28) . n = 4 Consider the fourth-order cumulant C ( h , h , h , h ) . If we limit the analysis to quadratic non-gaussianities asin f NL ( R ) = α R , it is clear that there are no corrections of order O ( α ) since f (3)NL ( R ) = 0 . The first non-trivialcorrection arises at order O ( α ) C ( h , h , h , h ) = C [ f NL ( R ) , f NL ( R ) , R , R ] + 5 combinations . (D51)If we take eq. (D41) it is simple to see that the first two terms on the right-hand side give vanishing contribution.In the last line, we have one contribution corresponding to a = b = 1 that gives two terms since one has to sum ↔ exchange. We have C [ f NL ( R ) , f NL ( R ) , R , R ] = (cid:104) (cid:104) f (1)NL ( R ) f (1)NL ( R ) (cid:105) − (cid:104) f (1)NL ( R ) (cid:105)(cid:104) f (1)NL ( R ) (cid:105) (cid:105) ( (cid:104)R R (cid:105)(cid:104)R R (cid:105) + 3 ↔ (cid:105) )= 4 α (cid:104)R R (cid:105) ( (cid:104)R R (cid:105)(cid:104)R R (cid:105) + (cid:104)R R (cid:105)(cid:104)R R (cid:105)(cid:105) ) . (D52)If we set all points equal and multiply by six because of eq. (D51), we find C ( h, h, h, h ) = 48 α σ [1 + O ( α σ )] , (D53)where we introduced again a sub-leading correction of order O ( α ) .We can check the validity of eqs. (D45, D50, D53) explicitly. The reason is that if we limit the analysis—as donehere—to random variables without derivatives , the computation of the cumulants of h admits an exact analyticalsolution. We can indeed extract the probability density distribution of the random variable h —at a given spatialpoint, that is what we need in order to compare with eqs. (D45, D50, D53)—by means of the Jacobi’s multivariatetheorem in eq. (D23) using the fact that R is a gaussian variable with zero mean and variance σ . We find p ( h ) = e −R / σ + e −R − / σ √ πσ √ αh , where R ± = − ± √ αh α , (D54)with h ∈ [ − / α, ∞ ) . By means of the Faà di Bruno’s formula, we find the generic n th -order cumulant C n ( h, . . . , h ) = 2 n − σ ( ασ ) n − ( n + 4 α σ )( n − , (D55)One can check that eqs. (D45, D50, D53) are correctly reproduced in the appropriate limits. ◦ Since in the case with no derivatives the probability density distribution can be obtained analytically (seeeq. (D54)), it is instructive to do the following exercise. First of all, let us rewrite eq. (D54) for the case with f NL ( R ) = α ( R − (cid:104)R (cid:105) ) with (cid:104)R (cid:105) = σ . We find p ( h ) = e −R / σ + e −R − / σ √ πσ (cid:112) α ( h + ασ ) , where R ± = − ± (cid:112) α ( h + ασ )2 α . (D56) From the probability density distribution we compute the moments µ n = (cid:82) ∞− / α dh h n p ( h ) . The explicit expression for the n th cumulant interms of the first n moments can be obtained by using Faà di Bruno’s formula for higher derivatives of composite functions; explicitly, for n (cid:62) ,we have C n = (cid:80) nk =1 ( − k − ( k − B n,k (0 , µ , . . . , µ n − k +1 ) , where B n,k are the incomplete Bell polynomials. C ( h ) = 0 while C n ( h, . . . , h ) ≡ C n ( h ) are still given by eq. (D55). We can compute analyt-ically the characteristic function log χ ( λ ) = (cid:80) ∞ n =2 ( i n /n !) C n ( h ) λ n . We find χ ( λ ) = exp (cid:20) − λ σ − iαλσ ) − iαλσ (cid:21) (cid:112) − iαλσ , (D57)and one can check for consistency that p ( h ) in eq. (D56) is given by the inverse Fourier transform p ( h ) = 12 π (cid:90) ∞−∞ dλ χ ( λ ) e − ihλ . (D58)Once we know p ( h ) , we can compute the integral ¯ p ( h c ) ≡ (cid:82) ∞ h c dh p ( h ) , which is the probability to find h abovethe threshold h c . We find ¯ p NG ( h c ) = 12 (cid:20) − Erf (cid:18) (cid:112) h c α + 4 α σ √ ασ (cid:19) − Erf (cid:18) − (cid:112) h c α + 4 α σ √ ασ (cid:19)(cid:21) , (D59)with Erf( x ) the error function. This simple case, therefore, can be solved exactly. This means that we can use itas a playground to test the following approximation. Let us organize the sum over cumulants in a power-seriesexpansion in α . We have ∞ (cid:88) n =2 i n n ! C n ( h ) λ n = − σ λ − iαλ σ + α σ ( − λ + 2 λ σ ) + O ( α ) , (D60)where the first term on the right-hand side corresponds to the quadratic cumulant that reproduces the gaus-sian limit. We can truncate eq. (D60) at some order in α , compute the corresponding ¯ p ( h c ) and compare witheq. (D59). At order α , we find the gaussian result ¯ p G ( h c ) = 12 Erf (cid:18) h c √ σ (cid:19) (cid:39) √ πv c exp (cid:18) − v c (cid:19) , (D61)where we define v c ≡ h c /σ , and the last approximation corresponds to v c (cid:29) which is the limit that is relevantfor our analysis. Consider now the order O ( α ) . Let us write eq. (D59) at order O ( α ) as ¯ p α ( h c ) = (cid:82) ∞ h c dh p α ( h ) with p α ( h ) = 12 π (cid:90) ∞−∞ dλ exp (cid:20) − σ λ − i ( ασ ) λ σ (cid:21) e − ihλ = 12 π (cid:90) ∞−∞ dλ e − σ λ / (cid:20) ∞ (cid:88) m =0 ( − iασ ) m m ! (cid:21) λ m e − ihλ (cid:124) (cid:123)(cid:122) (cid:125) i m ∂ me − ihλ∂h m = ∞ (cid:88) m =0 ( − iασ ) m m ! i m ∂ m ∂h m π (cid:90) ∞−∞ dλ e − σ λ / e − ihλ (cid:124) (cid:123)(cid:122) (cid:125) gaussian integral √ πσ e − h / σ ≡ ∞ (cid:88) m =0 p ( m ) α ( h ) . (D62)In practice, we can approximate ¯ p α ( h c ) as a sum of terms, ¯ p α ( h c ) = (cid:80) ∞ m =0 ¯ p ( m ) α ( h c ) each one of them defined bymeans of eq. (D62), that is ¯ p ( m ) α ( h c ) = (cid:82) ∞ h c dh p ( m ) α ( h ) . A similar decomposition can be defined at higher ordersin α . At order O ( α ) , we find ¯ p (1) α ( h c ) = 1 √ πv c exp (cid:18) − v c (cid:19) ( ασ ) v c (cid:18) − v c (cid:19) , (D63) ¯ p (2) α ( h c ) = 1 √ πv c exp (cid:18) − v c (cid:19) ( ασ ) v c (cid:18) − v c + 15 v c (cid:19) , (D64) ¯ p (3) α ( h c ) = 1 √ πv c exp (cid:18) − v c (cid:19) ( ασ ) v c (cid:18) − v c + 210 v c − v c + 105 v c (cid:19) (D65) ¯ p (4) α ( h c ) = . . . Since v c (cid:29) , we are tempted to approximate (by neglecting all sub-leading terms in each of the round brackets) ¯ p ( m ) α ( h c ) = 1 √ πv c exp (cid:18) − v c (cid:19) ( ασ v c ) m m ! = ⇒ ¯ p α ( h c ) = 1 √ πv c exp (cid:18) − v c ασ v c (cid:19) . (D66)1We are now in the position to compare i) the gaussian result in eq. (D61), ii) the exact non-gaussian result ineq. (D59), iii) the order O ( α ) approximation ¯ p α ( h c ) = (cid:80) ∞ m =0 ¯ p ( m ) α ( h c ) obtained by truncating the series forincreasing values of m and iv) the order O ( α ) approximation resummed as in eq. (D66). This comparison isshown in the left panel of fig. 11.We can now consider the order O ( α ) . We find ¯ p (1) α ( h c ) = 1 √ πv c exp (cid:18) − v c (cid:19) ( ασ ) (cid:20) v c (cid:18) − v c (cid:19) + 2 v c ασ (cid:18) − v c (cid:19)(cid:21) , (D67) ¯ p (2) α ( h c ) = 1 √ πv c exp (cid:18) − v c (cid:19) ( ασ ) × (cid:20) v c (cid:18) − v c + 15 v c (cid:19) + 2 ασ v c (cid:18) − v c + 42 v c − v c (cid:19) + 2( ασ ) v c (cid:18) − v c + 3814 v c − v c (cid:19)(cid:21) , (D68) ¯ p (3) α ( h c ) = . . . We are again tempted to consider the limit v c (cid:29) in which we keep only the leading term inside the roundbrackets, meaning that we write ¯ p (1) α ( h c ) = 1 √ πv c exp (cid:18) − v c (cid:19) ( ασ ) (cid:0) v c + 2 v c ασ (cid:1) , (D69) ¯ p (2) α ( h c ) = 1 √ πv c exp (cid:18) − v c (cid:19) ( ασ ) (cid:20) v c ασ v c + 2( ασ ) v c (cid:21) , (D70) ¯ p (3) α ( h c ) = . . . In this case, we find that the series ¯ p α ( h c ) = (cid:80) ∞ m =0 ¯ p ( m ) α ( h c ) can be resummed, and we get the neat expression ¯ p α ( h c ) (cid:39) √ πv c exp (cid:20) − v c ασ v c (1 + 2 ασ v c ) (cid:21) . (D71)We can do the same comparison we did before, and compare i) the gaussian result in eq. (D61), ii) the exactnon-gaussian result in eq. (D59), iii) the order O ( α ) approximation ¯ p α ( h c ) = (cid:80) ∞ m =0 ¯ p ( m ) α ( h c ) obtained bytruncating the series at increasing values of m and iv) the order O ( α ) approximation resummed as in eq. (D71).This comparison is shown in the right panel of fig. 11. � � � � � ��� - �� �� - �� �� - � �� - � �� - � � � ( ) ℊ ℴ ℊ ℯ ℴ � ( α ) � � � � � ��� - �� �� - �� �� - � �� - � �� - � � � ( ) ℊ ℴ ℊ ℯ ℴ � ( α � ) FIG. 11:
Left panel. We compare at order O ( α ) the exact non-gaussian expression for ¯ p ( h c ) (defined below eq. (D58)) given in eq. (D59)with: i) its gaussian limit (eq. (D61)), ii) its power-series expansion ¯ p α ( h c ) = (cid:80) ∞ m =0 ¯ p ( m ) α ( h c ) for increasing values of m (green region)and iii) the resummed expression (labelled “exp approx”) given in eq. (D66). For illustration, we take σ = 0 . and α = 0 . . Rightpanel. Same as in the left panel but at order O ( α ) . h without any information about itsderivatives) but, as we shall discuss later, we will use an approximation scheme very similar to the one adoptedabove. It is, therefore, instructive to draw some conclusions (this is particularly important in light of the fact thatin the case without derivatives we know the exact form of the probability density distribution). Consider firstthe O ( α ) approximation illustrated in the left panel of fig. 11. The sum ¯ p α ( h c ) = (cid:80) ∞ m =0 ¯ p ( m ) α ( h c ) (truncatedat some finite m ; in the plot we show the cases with m (cid:54) ) already gives a good—even though not optimal—approximation of ¯ p NG ( h c ) . On the contrary, the resummed expression in eq. (D66), more simple to handle, tendsto overestimate the true result. Consider now the O ( α ) approximation illustrated in the right panel of fig. 11. Thesum ¯ p α ( h c ) = (cid:80) ∞ m =0 ¯ p ( m ) α ( h c ) (truncated at some finite m ; in the plot we show the cases with m (cid:54) ) gives anexcellent approximation of ¯ p NG ( h c ) . On the contrary, it is evident that the resummed expression in eq. (D71)quickly diverges from the true result and can not be trusted.The reason why the exponential approximation deviates from the exact result can be identified already at order O ( α ) . Consider (in units of ¯ p G ( h c ) ) the first term in eq. (D63) that we kept in our expansion, that is ( ασ ) v c , andcompare it with one of the terms that we neglected in eq. (D65), for instance the fourth in the round brackets − ασ ) v c . These two terms have the same power of v c but the latter is parametrically suppressed by ( ασ ) .It is true that we expect ( ασ ) ∼ − (cid:28) but we also have an extra numerical factor − that partially com-pensate the suppression. On similar ground, neglecting the term − (14 / ασ ) v c in eq. (D65) compared to theterms that we kept in eq. (D63) and eq. (D64) seems not justifiable.We conclude that the exponential approximations in eq. (D66) and eq. (D71) can not be considered as acceptablygood proxies for the exact result since they may overestimate it by many orders of magnitude. However, in orderto have a good approximation of the exact result it is more proper to decompose ¯ p α ( h c ) = (cid:80) ∞ m =0 ¯ p ( m ) α ( h c ) (oreven better at order O ( α ) ) and sum over m up to some finite value according to the desired accuracy.After this digression, we are finally ready to compute the non-gaussian cumulants for the random variables h , h i , h ij . For simplicity, we start again from the two-dimensional case, and we have six random variables { h, h x , h y , h xx , h xy , h yy } . We consider the non-gaussian function f NL ( R ) = α ( R − (cid:104)R (cid:105) ) since in this case we have C ( h ) = (cid:104) h (cid:105) = 0 . Furthermore, in order to fully exploit the properties of homogeneity and isotropy (see discussionbelow eq. (B8)), it is useful to change basis to the complex conjugated random variables { h, h z , h z ∗ , h zz , h zz ∗ , h z ∗ z ∗ } . All first-order cumulants vanish. This is true for C ( h ) = (cid:104) h (cid:105) = 0 as discussed before. As far as the other first-ordercumulants are concerned, we have for instance C ( h z ) = (cid:104) h z (cid:105) = ∂ z (cid:104) h (cid:105) = 0 which vanishes (even without imposing C ( h ) = 0 ) because (cid:104) h (cid:105) does not depend, as a consequence of homogeneity, on the specific spatial point at which it iscomputed. Similarly, we have C also for the remaining random variables.As far as the second-order cumulants are concerned, we find the following non-zero entries h h z h z ∗ h zz h zz ∗ h z ∗ z ∗ h C ( h, h ) 0 0 0 C ( h, h zz ∗ ) 0 h z C ( h z , h z ∗ ) 0 0 0 h z ∗ C ( h z , h z ∗ ) 0 0 0 0 h zz C ( h zz , h z ∗ z ∗ ) h zz ∗ C ( h, h zz ∗ ) 0 0 0 C ( h zz ∗ , h zz ∗ ) 0 h z ∗ z ∗ C ( h zz , h z ∗ z ∗ ) 0 0 , (D73)with C ( h, h ) = σ (1 + 2 α σ ) , (D74) C ( h z , h z ∗ ) = σ α σ ) , (D75) C ( h, h zz ∗ ) = − σ α σ ) , (D76) C ( h zz , h z ∗ z ∗ ) = C ( h zz ∗ , h zz ∗ ) = σ (cid:20) α (cid:18) σ σ + σ (cid:19)(cid:21) . (D77) The following relations turn out to be useful. For a generic function f ( x, y ) : R → R we have f z = ( f x − if y ) / , f z ∗ = ( f x + if y ) / and f zz ∗ = 14 ( f xx + f yy ) , f zz = 14 ( f xx − f yy − if xy ) , f z ∗ z ∗ = 14 ( f xx − f yy + 2 if xy ) . (D72)The stationary point condition f x = f y = 0 becomes f z = 0 (hence f z ∗ = 0 ). The condition f xx f yy − f xy > that separates extrema fromsaddle points becomes f zz ∗ > f zz f z ∗ z ∗ = | f zz | that is f zz ∗ < −| f zz | ∨ f zz ∗ > | f zz | (equivalently, | f zz ∗ | > | f zz | ). Notice that f zz ∗ is real.Maxima correspond to the condition f zz ∗ < while minima are identified by f zz ∗ > . O ( α ) , all non-gaussian corrections vanish. We already computed C ( h, h ) in eq. (D45). It is instructive toconsider explicitly few more cases. The relevant equations are eq. (D10), eq. (D33), eq. (D37) and eq. (D41).Consider for instance the computation of C ( h z , h z ∗ ) . We have C ( h z , h z ∗ ) eq.(D10) = ∂ z ∂ z ∗ C ( h , h ) (cid:12)(cid:12) (cid:126)x = (cid:126)x (D78) = ∂ z ∂ z ∗ { C ( R , R ) + C [ f NL ( R ) , R ] + C [ R , f NL ( R )] + C [ f NL ( R ) , f NL ( R )] } (cid:12)(cid:12) (cid:126)x = (cid:126)x = ∂ z ∂ z ∗ (cid:2) (cid:104)R R (cid:105) + (cid:104) f (1)NL ( R ) (cid:105)(cid:104)R R (cid:105) + (cid:104) f (1)NL ( R ) (cid:105)(cid:104)R R (cid:105) (cid:124) (cid:123)(cid:122) (cid:125) O ( α ) , eq . (D33) + (cid:104) f NL ( R ) f NL ( R ) (cid:105) − (cid:104) f NL ( R ) (cid:105)(cid:104) f NL ( R ) (cid:105) (cid:124) (cid:123)(cid:122) (cid:125) O ( α ) , eq . (D37) (cid:3)(cid:12)(cid:12) (cid:126)x = (cid:126)x = (cid:104)R z R z ∗ (cid:105) + 2 (cid:104) f (1)NL ( R ) (cid:105)(cid:104)R z R z ∗ (cid:105) + α ∂ z ∂ z ∗ (cid:104)R R (cid:105) (cid:12)(cid:12) (cid:126)x = (cid:126)x (D79) = (cid:104)R z R z ∗ (cid:105) + 4 α (cid:104)R R z R z ∗ (cid:105) . (D80)Eq. (D79) follows from the fact that the expectation values (cid:104) f NL ( R ) (cid:105) and (cid:104) f (1)NL ( R ) (cid:105) do not depend, because of ho-mogeneity, on the specific spatial point at which they are evaluated and, therefore, their spatial derivatives van-ish. Eq. (D80) follows from the fact that (cid:104) f (1)NL ( R ) (cid:105) = 2 α (cid:104)R(cid:105) = 0 ; at order O ( α ) , we moved the derivatives insidethe statistical average and set (cid:126)x = (cid:126)x . In eq. (D80) we have (cid:104)R z R z ∗ (cid:105) = ( (cid:104)R x R x (cid:105) + (cid:104)R y R y (cid:105) ) / σ / and (cid:104)R R z R z ∗ (cid:105) = (cid:104)R ( R x + R y ) (cid:105) / . This statistical average can be computed by means of the multivariate normaldistribution in eq. (B13) (cid:104)R ( R x + R y ) (cid:105) = (cid:90) d R d R x d R y d R xx d R xy d R yy R ( R x + R y ) P ( R x ) P ( R y ) P ( R xy ) P ( R , R xx , R yy ) (D81) = (cid:90) d R x d R y ( R x + R y ) P ( R x ) P ( R y ) (cid:124) (cid:123)(cid:122) (cid:125) = σ (cid:90) d R xy P ( R xy ) (cid:124) (cid:123)(cid:122) (cid:125) = 1 (cid:90) d RR (cid:90) d R xx d R yy P ( R , R xx , R yy ) (cid:124) (cid:123)(cid:122) (cid:125) = √ πσ exp( −R / σ ) , and, from the explicit computation of the integrals, we get (cid:104)R ( R x + R y ) (cid:105) = σ σ . All in all, we find C ( h z , h z ∗ ) = σ (1 + 4 α σ ) / . All the remaining non-zero entries in eq. (D73) can be derived in the same way.In eq. (D73) the zero entries vanish as a consequence of isotropy. We can introduce again the parameter κ ≡ ( z ∗ derivatives in C ) − ( z derivatives in C ) and set to zero all cumulants with κ (cid:54) = 0 since not invariant underspatial rotations. The reasoning that we followed in the gaussian case for the computations of the two-point correla-tors (see discussion below eq. (B8)) applies also in the case of generic n th -order cumulants. This is because cumulantsare functions of moments (see eqs. (D4-D7)). Furthermore, the relation between C ( h z , h z ∗ ) and C ( h, h zz ∗ ) can beunderstood as a consequence of homogeneity.We now move to consider the third-order cumulants. Based on isotropy, we expect only a handful of non-zero cu-mulants (that are those with κ = 0 ) that we list in table II. No derivatives 2 derivatives 4 derivatives 6 derivatives C ( h, h, h ) C ( h, h, h zz ∗ ) C ( h, h zz ∗ , h zz ∗ ) C ( h zz ∗ , h zz ∗ , h zz ∗ ) × × × × C ( h, h z , h z ∗ ) C ( h z , h z ∗ , h zz ∗ ) C ( h zz , h z ∗ z ∗ , h zz ∗ ) × × × C ( h, h zz , h z ∗ z ∗ ) × C ( h z , h z , h z ∗ z ∗ ) × C ( h z ∗ , h z ∗ , h zz ) × TABLE II:
Third-order cumulants for the random variables { h, h z , h z ∗ , h zz , h zz ∗ , h z ∗ z ∗ } that are non-zero based on isotropy. For thesecumulants, we have κ ≡ ( z ∗ derivatives in C ) − ( z derivatives in C ) = 0 . We gather together in each column cumulants withthe same number of spatial derivatives (cumulants without derivatives in the first column, with two derivatives in the second, four inthe third and so on). For each entry, the tiny numbers in the second row indicate the multiplicity of the corresponding cumulant dueto distinct permutations of its arguments. C ( h, h, h ) = 6 ασ in eq. (D50). As far as the remaining cumulants are concerned, we find C ( h, h z , h z ∗ ) = ασ σ , (D82) C ( h, h, h zz ∗ ) = − ασ σ , (D83) C ( h, h zz ∗ , h zz ∗ ) = α (cid:0) σ + 2 σ σ (cid:1) , (D84) C ( h, h zz , h z ∗ z ∗ ) = α σ σ , (D85) C ( h z , h z ∗ , h zz ∗ ) = − α σ , (D86) C ( h z , h z , h z ∗ z ∗ ) = C ( h z ∗ , h z ∗ , h zz ) = α σ , (D87) C ( h zz ∗ , h zz ∗ , h zz ∗ ) = − α σ σ , (D88) C ( h zz , h z ∗ z ∗ , h zz ∗ ) = − α σ σ . (D89)We find that corrections at order O ( α ) vanish for all third-order cumulants listed before . However, as already dis-cussed in the case of eq. (D50), eqs. (D82-D89) are not exact because we expect the presence of O ( α ) correctionswe did not include. It is simple to check that properties of homogeneity are respected by the explicit expressions ineqs. (D82-D89). For instance, from ∂ z ∗ C ( h z , h z , h z ∗ ) = 0 we find C ( h z , h z ∗ , h zz ∗ ) = − C ( h z , h z , h z ∗ z ∗ ) that is in-deed verified by eq. (D86) and eq. (D87). Similarly, from ∂ z ∗ C ( h, h, h z ) = 0 we have C ( h, h, h zz ∗ ) = − C ( h, h z , h z ∗ ) that is indeed verified by eq. (D82) and eq. (D83).After computing the cumulants at the desired order in α , we are finally ready to take the last step in our computation.From the cumulants, we will reconstruct the characteristic function and, via an inverse Fourier transform, the prob-ability density distribution. Few technical remarks are in order because it is crucial to properly identify the variablesparticipating to the Fourier transforms. The inverse Fourier transform of the first line in eq. (D1) is given by P ( ξ , . . . , ξ N ) = (cid:90) dλ (2 π ) . . . dλ N (2 π ) χ ( λ , . . . , λ N ) exp [ − i ( ξ λ + · · · + ξ N λ N )] . (D90)Both eq. (D1) and eq. (D90) are valid for real variables. The simplest identification, therefore, would be { ξ i =1 ,..., } = { h, h x , h y , h xx , h xy , h yy } . However, we found more efficient to work with { h, h z , h z ∗ , h zz , h zz ∗ , h z ∗ z ∗ } since propertieslike isotropy become more transparent. Among these variables, h and h zz ∗ are real while h z and h zz are complex (withconjugated variables given by h z ∗ and h z ∗ z ∗ , respectively). In this case a suitable choice of real variables, therefore, is { ξ i =1 ,..., } = { h, Re h z , Im h z , Re h zz , h zz ∗ , Im h zz } with the obvious relations Re h z = ( h z + h z ∗ ) / , Im h z = − i ( h z − h z ∗ ) / , Re h zz = ( h zz + h z ∗ z ∗ ) / , Im h zz = − i ( h zz − h z ∗ z ∗ ) / . (D91)Consequently, we introduce their Fourier counterparts { λ i =1 ,..., } ≡ { λ, λ Re h z , λ Im h z , λ Re h zz , λ zz ∗ , λ Im h zz } such thatthe characteristic function in eq. (D1) is given by the integral (cid:90) dh d Re h z d Im h z d Re h zz dh zz ∗ d Im h zz P ( h, h z , h z ∗ , h zz , h zz ∗ , h z ∗ z ∗ ) exp[ i ( hλ + Re h z λ Re h z + Im h z λ Im h z + . . . )] . (D92)Consider, for instance, the two terms Re h z λ Re h z + Im h z λ Im h z in the argument of the exponential. By means ofeqs. (D91) we can simply rewrite Re h z λ Re h z + Im h z λ Im h z = h z λ ∗ z + h z ∗ λ z if we define λ z ≡ ( λ Re h z + iλ Im h z ) / .Notice that in this case we have λ Re h z = 2Re λ z and λ Im h z = 2Im λ z . We can, therefore, use λ z as the complex Fouriervariable associated to h z ∗ (and, correspondingly, its conjugated λ ∗ z associated to h z ). Similarly, we introduce λ zz ≡ ( λ Re h zz + iλ Im h zz ) / , and the exponential function in eq. (D92) reads exp { i [ hλ + h zz ∗ λ zz ∗ + ( h z λ ∗ z + h zz λ ∗ zz + c.c. )] } .We are now in the position to consider explicitly the inverse Fourier transform of eq. (D92). As integration variables, This happens because the O ( α ) contribute to C ( h , h , h ) is: ∆ C ( h , h , h ) = 2 α (cid:2) (cid:104)R R (cid:105)(cid:104)R ( R + R ) (cid:105) + (cid:104)R R (cid:105)(cid:104)R ( R + R ) (cid:105) + (cid:104)R R (cid:105)(cid:104)R ( R + R ) (cid:105) (cid:3) Thus, when we derive and then we set (cid:126)x = (cid:126)x = (cid:126)x , following the prescription (D10) and we return to the { x, y } basis, we see that all the termshave the form (cid:104) ξ j ξ j ξ j (cid:105) , where ξ j k ∈ {R , R x , R y , R xx , R xy , R yy } . These are third order moments of a gaussian multivariate distribution,which are vanishing. λ Re h z , λ Im h z , λ Re h zz , λ Im h zz we use λ z and λ zz introduced before by means of the relations λ Re h z = 2Re λ z , λ Im h z = 2Im λ z , λ Re h zz = 2Re λ zz and λ Im h zz = 2Im λ zz . We have P ( h, h z , h z ∗ , h zz , h zz ∗ , h z ∗ z ∗ ) = (D93) (cid:90) dλ (2 π ) d Re λ z π d Im λ z π dλ zz ∗ (2 π ) d Re λ zz π d Im λ zz π χ ( λ, λ z , λ zz , λ zz ∗ ) exp {− i [ hλ + h zz ∗ λ zz ∗ + ( h z λ ∗ z + h zz λ ∗ zz + c.c. )] } , with of course λ z = Re λ z + i Im λ z and λ zz = Re λ zz + i Im λ zz . The natural log of the characteristic function log χ ( λ, λ z , λ zz , λ zz ∗ ) is written, according to eq. (D3), as a series expansion in terms of the cumulants with respectto the Fourier variables { λ, λ z , λ zz , λ zz ∗ } (remember that λ z and λ zz are complex variables while λ and λ zz ∗ are real).The first non-gaussian correction arises at order O ( α ) and reads log χ ( λ,λ z , λ zz , λ zz ∗ ) = (D94) − C ( h, h ) λ − C ( h z , h z ∗ ) | λ z | − C ( h, h zz ∗ ) λλ zz ∗ − C ( h zz , h z ∗ z ∗ ) | λ zz | − C ( h zz ∗ , h zz ∗ ) λ zz ∗ − i C ( h, h, h ) λ − i C ( h, h, h zz ∗ ) λ λ zz ∗ − iC ( h, h z , h z ∗ ) λ | λ z | − i C ( h, h zz ∗ , h zz ∗ ) λλ zz ∗ − iC ( h, h zz , h z ∗ z ∗ ) λ | λ zz | − iC ( h z , h z ∗ , h zz ∗ ) | λ z | λ zz ∗ − i C ( h z , h z , h z ∗ z ∗ )( λ ∗ z ) λ zz − i C ( h z ∗ , h z ∗ , h zz ) λ z λ ∗ zz − i C ( h zz ∗ , h zz ∗ , h zz ∗ ) λ zz ∗ − iC ( h zz , h z ∗ z ∗ , h zz ∗ ) | λ zz | λ zz ∗ , where the second-order cumulants given in eq. (D73) are taken at order O ( α ) and, therefore, coincide with their gaus-sian limit. Notice that each term in the series expansion enters with a coefficient that counts its multiplicity (explicitlywritten in table II). This follows from the fact that in the last line in eq. (D1) there are different terms that give the samecontribution. The integration in eq. (D93) gives the probability density distribution P ( h, h z , h z ∗ , h zz , h zz ∗ , h z ∗ z ∗ ) . Af-ter computing P ( h, h z , h z ∗ , h zz , h zz ∗ , h z ∗ z ∗ ) , it is also possible to reconstruct P ( h, h x , h y , h xx , h xy , h yy ) by means ofthe Jacobi’s multivariate theorem in eq. (D14) using the transformation in eq. (D72).It is possible to check this procedure in the gaussian limit α → . In the gaussian limit only the second-order cumu-lants survive (with α → taken in the corresponding expressions given in eq. (D73)). The inverse Fourier transform ineq. (D93) can be computed analytically in the gaussian limit, and one gets P ( R , R z , R z ∗ , R zz , R zz ∗ , R z ∗ z ∗ ) where weused R instead of h since α → . The determinant of the Jacobian matrix is | det J ( R , R x , R y , R xx , R xy , R yy ) | = 1 / ,and Jacobi’s multivariate theorem gives precisely eq. (B13). This is a non-trivial check that our procedure is technicallycorrect.In the presence of local non-gaussianities, the computation of eq. (D93) is much more complicated because χ ( λ, λ z , λ zz , λ zz ∗ ) in eq. (D94) is an exponential function whose argument is a polynomial of cubic order. A possi-ble solution strategy is the following. Let us write eq. (D93) in the schematic form (we use for illustration a generic setof variables { λ i =1 ,...,N } that in the actual computation must be replaced with { λ, λ z , λ zz , λ zz ∗ } ) χ ( λ i ) = exp[ p ( λ i ) + p ( λ i )] = exp[ p ( λ i )] exp[ p ( λ i )] = exp[ p ( λ i )] ∞ (cid:88) m =0 m ! p ( λ i ) m , (D95)where p ( λ i ) is the quadratic polynomial in { λ i } that corresponds to the gaussian limit α → while p ( λ i ) is the cubicpolynomial that contains O ( α ) deviations. At a given order m in the series expansion defined by eq. (D95), the inverseFourier integral in eq. (D93) can be computed analytically. This is because we can make use of the following identity(the sum over k is understood) (cid:90) dλ . . . λ N ( λ ni . . . λ mj ) exp[ p ( λ i )] exp( − iλ k h k ) = ( i ) n . . . ( i ) m ∂ n ∂h ni . . . ∂ m ∂h mj (cid:90) dλ . . . λ N exp[ p ( λ i )] exp( − iλ k h k ) , (D96)where the last integral is nothing but the gaussian integral that can be easily computed. We go along with this strategy,that we summarize for ease of reading in the following. i) For fixed m , we compute the probability density distribution in eq. (D93) performing the inverse Fourier trans-form by means of the trick in eq. (D95) and eq. (D96). ii) We transform P ( h, h z , h z ∗ , h zz , h zz ∗ , h z ∗ z ∗ ) into P ( h, h x , h y , h xx , h xy , h yy ) by means of Jacobi’s multivariate the-orem (as illustrated before for the gaussian case). iii) We set h x = h y = 0 , and compute the number density of local maxima n max ( h ) according to the general defini-tion in eq. (B4) (with of course h instead of R , and working in two spatial dimensions).6 iv) We use the same change of variables proposed in eq. (B17) (with again h instead of R ), and we integrate over θ , r , h and s in the domain defining maxima, as discussed in the gaussian case; similarly, we implement the samelower bound of integration over s given by the threshold condition for black hole formation in eq. (B36). Noticethat in the following we will set the lower limit of integration over h to be h min = −∞ despite the fact that we have h min = − (1 + 4 α σ ) / α for h = R + α ( R − σ ) and R ∈ ( −∞ , + ∞ ) . Approximating h min with its gaussianvalue allows to obtain close analytical formulas. We shall justify the validity of this approximation at the end ofthis section. v) Finally, we define the mass fraction of black holes as in eq. (B35).We now illustrate the result of this procedure. Needless to say, the term with m = 0 reproduces the gaussian result ineq. (B35). The first correction is given by the term with m = 1 . Let us introduce for the sake of simplicity the notation β = β G + (cid:80) ∞ m =1 β ( α,m )NG with β G given by eq. (B35). We find β ( α, = (cid:18) ασ σ (cid:19) √ π / e − s /σ (cid:20) s σ + e s /σ (cid:18) − s σ + 4 s σ (cid:19)(cid:21) . (D97)For m = 2 , we find β ( α, = (cid:18) ασ σ (cid:19) e − s /σ π / (cid:26) √ πe s /σ Erfc (cid:18) √ s min σ (cid:19) + 3 √ s min σ (cid:20) − s σ + 48 s σ + e s /σ (cid:18) − s σ − s σ + 64 s σ (cid:19)(cid:21)(cid:27) . (D98)For increasing values of m the corresponding expressions of β ( m )NG become more lengthy and less transparent, and we donot report them explicitly. On the contrary, we make use of the same expansion introduced in eq. (B35) for s min /σ (cid:38) .Interestingly, in this limit we find that the generic term β ( α,m )NG can be written as β ( α,m )NG (cid:39) (cid:20)(cid:18) s σ (cid:19) ασ σ (cid:21) m m − / π / m ! (cid:18) s min σ (cid:19) e − s /σ , (D99)and the sum (cid:80) ∞ m =1 β ( α,m )NG admits the following analytical form ∞ (cid:88) m =1 β ( α,m )NG (cid:39) √ π / (cid:18) s min σ (cid:19) e − s /σ (cid:26) − (cid:20) (cid:18) s σ (cid:19) ασ σ (cid:21)(cid:27) . (D100)The sum β = β G + (cid:80) ∞ m =1 β ( α,m )NG gives β (cid:39) √ π / (cid:18) s min σ (cid:19) exp (cid:20) − s σ + 16 (cid:18) s σ (cid:19) ασ σ (cid:21) . (D101)Notice that the structure of this equation is analogue to the one that we found in eq. (D66).Before proceeding, let us come back to the issue of the lower limit of integration over h . As anticipated, from the verysame definition h = R + α ( R − σ ) it follows that h > h min = − (1 + α σ ) / α . To be precise, therefore, this lowerlimit of integration should be implemented when computing the non-gaussian contributions β ( α,m )NG . However, imple-menting the gaussian lower limit of integration h min = −∞ does not change qualitatively our results, and the reasonis the following. When computing β , we only focus on maxima with large curvature −(cid:52) h/σ since they correspond topeak of the overdensity field. The value of −(cid:52) h/σ is positively correlated with the value of h/σ , and the amount ofcorrelation is controlled by the parameter γ . This fact was illustrated in appendix B in the gaussian case, and it remainsvalid also in the presence of non-gaussianities (for the values of α that are relevant for our analysis). When γ → ,the correlation is maximal and regions with large −(cid:52) h/σ are also regions with positive h/σ (cid:38) . This means thatregions where h takes negative values are completely irrelevant for our purposes. In the cases that are relevant for ouranalysis, we typically have γ ∼ . ; we checked numerically (using for simplicity the log-normal power spectrum, andintegrating numerically over s ) that the value of β ( α,m )NG computed with the proper lower bound of integration over h re-main equal to those obtained with the gaussian limit h min = −∞ . In order to see some deviation, one should considerthe opposite limit γ → . In this limit −(cid:52) h/σ and h/σ are only weakly correlated and it is therefore possible to findregions with large −(cid:52) h/σ in which h takes negative values. In such situation, of course, implementing the correctlower limit of integration over h becomes important. However, cases with γ → fall outside the class of models we areconsidering in this paper.7 � ��� � ��� � ����� - �� �� - �� �� - �� �� - �� �� - �� �� - �� �� - � / σ � β ( ) ℊ ℴ ℊ ℯ ℴ � ( α ) � ��� � ��� � ����� - �� �� - �� �� - �� �� - �� �� - �� �� - �� �� - � / σ � β ( ) ℊ ℴ ℊ ℯ ℴ � ( α � ) FIG. 12:
Left panel. We compare the gaussian approximation for β given in eq. (B35) with i) its power-series expansion β = β G + (cid:80) ∞ m =1 β ( α,m )NG computed at order O ( α ) for increasing values of m (green region, with the thin solid green lines that label the case m = 1 , . . . , ) and ii) the resummed expression (labelled “exp approx”) given in eq. (D101). Right panel. Same as in the left panel butat order O ( α ) and up to m = 4 . In both panels we consider the log-normal model for the power spectrum, and we use the benchmarkvalues v = 0 . and P R ( k (cid:63) ) = 5 × − . We take α = 0 . . We compare in the left panel of fig. 12 the series β = β G + (cid:80) m max m =1 β ( α,m )NG , truncated at some finite m max , withthe resummed expression in eq. (D101). We show the cases m max = 1 , . . . , and we use the log-normal powerspectrum in eq. (B24) to compute the spectral parameters. The situation is very similar to what already discussed infig. 11, and we can exploit the knowledge that we gained from that case to argue the following conclusions. The series β = β G + (cid:80) m max m =1 β ( α,m )NG shows a convergence towards the order O ( α ) approximation (that is, we recall, the approxi-mation in which we only include the leading part of the third-order cumulants) of the exact non-gaussian result whilethe exponential approximation overestimates the abundance.We move to consider corrections at order O ( α ) . We already computed in eqs. (D74-D77) the second-order cumu-lants at order O ( α ) . In addition, we need to compute the fourth-order cumulants at order O ( α ) . We computed ineq. (D53) the simplest one without derivatives, C ( h, h, h, h ) = 48 α σ . The remaining non-zero cumulants are sum- No derivatives 2 derivatives 4 derivatives 6 derivatives 8 derivatives C ( h, h, h, h ) C ( h, h, h, h zz ∗ ) C ( h, h, h zz , h z ∗ z ∗ ) C ( h z , h z ∗ , h zz ∗ , h zz ∗ ) C ( h zz , h zz , h z ∗ z ∗ , h z ∗ z ∗ ) × × × × × C ( h, h, h z , h z ∗ ) C ( h, h, h zz ∗ , h zz ∗ ) C ( h z , h z , h z ∗ z ∗ , h zz ∗ ) C ( h zz ∗ , h zz ∗ , h zz ∗ , h zz ∗ ) × × × × C ( h, h z , h z ∗ , h zz ∗ ) C ( h z ∗ , h z ∗ , h zz , h zz ∗ ) C ( h zz , h z ∗ z ∗ , h zz ∗ , h zz ∗ ) × × × C ( h, h z , h z , h z ∗ z ∗ ) C ( h, h zz ∗ , h zz ∗ , h zz ∗ ) × × C ( h, h z ∗ , h z ∗ , h zz ) C ( h, h zz , h z ∗ z ∗ , h zz ∗ ) × × C ( h z , h z , h z ∗ , h z ∗ ) C ( h z , h z ∗ , h zz , h z ∗ z ∗ ) × × TABLE III:
Fourth-order cumulants for the random variables { h, h z , h z ∗ , h zz , h zz ∗ , h z ∗ z ∗ } that are non-zero based on isotropy. Forthese cumulants, we have κ ≡ ( z ∗ derivatives in C ) − ( z derivatives in C ) = 0 . As in table II, we gather together in eachcolumn cumulants with the same number of spatial derivatives (cumulants without derivatives in the first column, with two deriva-tives in the second, four in the third and so on). For each entry, the tiny numbers in the second row indicate the multiplicity of thecorresponding cumulant due to distinct permutations of its arguments. C ( h, h, h, h zz ∗ ) = − α σ σ , (D102) C ( h, h, h z , h z ∗ ) = 6 α σ σ , (D103) C ( h, h, h zz , h z ∗ z ∗ ) = 3 α σ σ , (D104) C ( h, h, h zz ∗ , h zz ∗ ) = α σ (10 σ + 3 σ σ ) , (D105) C ( h, h z , h z ∗ , h zz ∗ ) = − α σ σ , (D106) C ( h, h z , h z , h z ∗ z ∗ ) = C ( h, h z ∗ , h z ∗ , h zz ) = α σ σ , (D107) C ( h z , h z , h z ∗ , h z ∗ ) = 2 α σ σ , (D108) C ( h z , h z ∗ , h zz ∗ , h zz ∗ ) = α σ ( σ + σ σ ) , (D109) C ( h z , h z , h z ∗ z ∗ , h zz ∗ ) = C ( h z ∗ , h z ∗ , h zz , h zz ∗ ) = 0 , (D110) C ( h, h zz ∗ , h zz ∗ , h zz ∗ ) = − α σ ( σ + 2 σ σ ) , (D111) C ( h, h zz , h z ∗ z ∗ , h zz ∗ ) = − α σ σ σ , (D112) C ( h z , h z ∗ , h zz , h z ∗ z ∗ ) = α σ ( σ + σ σ ) , (D113) C ( h zz , h zz , h z ∗ z ∗ , h z ∗ z ∗ ) = α σ σ , (D114) C ( h zz ∗ , h zz ∗ , h zz ∗ , h zz ∗ ) = 3 α σ (3 σ + σ σ ) , (D115) C ( h zz , h z ∗ z ∗ , h zz ∗ , h zz ∗ ) = α σ (3 σ + 2 σ σ ) , (D116)that we derive as before using eq. (D10) and eq. (D41). Conceptually, the rest of the computation follows again points i) - v) discussed below eq. (D96). At the technical level, instead of eq. (D95) we now have χ ( λ i ) = exp[ p ( λ i ) + p ( λ i )] = exp[ p ( λ i )] exp[ p ( λ i )] = exp[ p ( λ i )] ∞ (cid:88) m =0 m ! p ( λ i ) m , (D117)where p ( λ i ) is the quartic polynomial that contains O ( α ) and O ( α ) deviations. The computation of the abundancefollows the same prescription discussed before, and we introduce the series expansion β = β G + (cid:80) m max m =1 β ( α ,m )NG with β G given by eq. (B35) and β ( α ,m )NG which corresponds to the series in eq. (D117) truncated at some m > . For instance,we find (with ¯ s min ≡ s min /σ and γ = σ /σ σ ) β ( α , = (D118) √ ασ ¯ s min π / e − s (cid:2) α (2 + γ ) σ − γ ¯ s min − α (5 + 8 γ ) σ ¯ s + 8 γ ¯ s + 16 α (1 + 3 γ ) σ ¯ s (cid:3) + O ( e − s ) . The exact analytic expressions for β ( α ,m )NG are quite lengthly and we do not report them here explicitly. The sum (cid:80) ∞ m =1 β ( α ,m )NG admits an analytical expression only if we retain, in each one of the β ( α ,m )NG , the highest power of ¯ s min in the polynomial that multiplies the leading exponential suppression. We find β (cid:39) √ π / (cid:18) s min σ (cid:19) exp (cid:40) − s σ + 16 (cid:18) s σ (cid:19) ασ σ + 32 (cid:18) s σ (cid:19) α (cid:34) (cid:18) σ σ (cid:19) + σ (cid:35)(cid:41) . (D119)The structure of this expression is analogue to eq. (D71). As shown in fig. 12, the exponential approximation containedin eqs. (D101) and (D119) gets worse for increasing s min /σ and diverges from the actual result, overstimating theabundance. The failure of the exponential approximation is analogue to the one discussed in section III B in the contextof threshold statistics.9 Appendix E: Towards an exact computation
The generic n -th order cumulant can be exactly computed if we specify the explicit form of the non-gaussianity. Inour case f NL ( R ) = α ( R − (cid:104)R (cid:105) ) .The translation only affects the first-order cumulant: C [ R + α ( R − (cid:104)R (cid:105) )] = C ( R + α R ) − α (cid:104)R (cid:105) , (E1) C n [ R + α ( R − (cid:104)R (cid:105) ) , . . . , R n + α ( R n − (cid:104)R (cid:105) )] = C n ( R + α R , . . . , R n + α R n ) , ∀ n > . (E2)Therefore we can set, without loss of generality, f NL ( R ) = α R in the computation of the cumulants of order higherthan one. The quadratic form of the non-linearity allows us to use the Gaussian integral extended to the complex planethanks to analytic continuation.We know that, in general C n ( h , . . . , h n ) is sum of a series of contributes of different orders of which the last onemust be O ( α n ) and it is represented by C n ( α R , . . . , α R n ) .
1. Cancellation of the contributes up to O ( α n − ) As a first step, we want to prove that C n ( α R , . . . , α R k , R k +1 , . . . , R n ) = 0 for k < n − . This implies that C n ( h , . . . , h n ) starts at order O ( α n − ) . Using eq. (D20): C n ( α R , . . . , α R k , R k +1 , . . . , R n ) = (E3) = ( − i ) n ∂∂λ . . . ∂∂λ n log (cid:90) d R . . . d R n e i (cid:80) ki =1 α R i λ i + i (cid:80) ni = k +1 R i λ i p ( R , . . . , R n ) (cid:12)(cid:12)(cid:12)(cid:12) λ = ... = λ n =0 = (E4) = ( − i ) n ∂∂λ . . . ∂∂λ n log (cid:90) d R . . . d R k e i (cid:80) ki =1 α R i λ i (cid:20)(cid:90) d R k +1 . . . d R n e i (cid:80) ni = k +1 R i λ i p ( R , . . . , R n ) (cid:21)(cid:124) (cid:123)(cid:122) (cid:125) ≡I (cid:48) ( R ,..., R k ,λ k +1 ,...,λ n ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) λ =0 . (E5)We know that (cid:90) d R . . . d R k e i (cid:80) ki =1 R i λ i I (cid:48) ( R , . . . , R k , λ k +1 , . . . , λ n ) = χ ( λ , . . . , λ n ) = e − (cid:80) ni,j =1 σ ij λ i λ j , (E6)and so I (cid:48) ( R , . . . , R k , λ k +1 , . . . , λ n ) = (cid:90) dλ π . . . dλ k π exp (cid:20) − i k (cid:88) i =1 λ i R i − (cid:18) k (cid:88) i,j =1 σ ij λ i λ j + +2 k (cid:88) i =1 n (cid:88) j = k +1 σ ij λ i λ j + n (cid:88) i,j = k +1 σ ij λ i λ j (cid:19)(cid:21) . (E7)Let us use the matrix notation, defining the following objects: λ ≡ ( λ , . . . , λ k ) T = ( λ i ) ki =1 , R ≡ ( R , . . . , R k ) T = ( R i ) ki =1 , (E8) ˆ σ ≡ ( σ ij ) ki,j =1 , σ j ≡ ( σ j , . . . , σ kj ) T = ( σ ij ) ki =1 with j ∈ { k + 1 , . . . , n } . (E9)In this way I (cid:48) ( R , λ k +1 , . . . , λ n ) = e − (cid:80) ni,j = k +1 σ ij λ i λ j (cid:90) d k λ (2 π ) k e − λ T ˆ σ λ + λ · ( − i R − (cid:80) nj = k +1 λ j σ j ) (E10) = 1 (cid:112) (2 π ) k det ˆ σ exp i R + n (cid:88) j = k +1 λ j σ j T ˆ σ − (cid:32) i R + n (cid:88) l = k +1 λ l σ l (cid:33) − n (cid:88) i,j = k +1 σ ij λ i λ j . Now, let us consider the integral I ( λ , . . . , λ n ) = (cid:90) d R . . . d R k e i (cid:80) ki =1 α R i λ i I (cid:48) ( R , λ k +1 , . . . , λ n ) ; (E11)0we can write k (cid:88) i =1 α R i λ i = k (cid:88) i,j =1 α R i λ i δ ij R j = α R T Λ R , (E12)where we defined the diagonal matrix Λ = diag( λ , . . . , λ k ) . In this way, we have I ( λ , . . . , λ n )= 1 (cid:112) (2 π ) k det ˆ σ (cid:90) d k R exp iα R T Λ R + 12 i R + n (cid:88) j = k +1 λ j σ j T ˆ σ − (cid:32) i R + n (cid:88) l = k +1 λ l σ l (cid:33) − n (cid:88) i,j = k +1 σ ij λ i λ j (E13) = 1 (cid:112) (2 π ) k det ˆ σ e ( (cid:80) nj = k +1 λ j σ j ˆ σ − (cid:80) nl = k +1 λ l σ l − (cid:80) ni,j = k +1 σ ij λ i λ j ) (cid:90) d k R e − R T ( ˆ σ − − iα Λ ) R + i R T ˆ σ − (cid:80) nj = k +1 λ j σ j (E14) = 1 (cid:112) det ( − iα ˆ σ Λ) exp − n (cid:88) j = k +1 λ j σ T j ˆ σ − (cid:104) ( − iα ˆ σ Λ) − − (cid:105) n (cid:88) l = k +1 λ l σ l − n (cid:88) i,j = k +1 σ ij λ i λ j , (E15)where we used that det(ˆ σ − − iα Λ) = det ˆ σ − det( − iασ Λ) and ˆ σ T = ˆ σ . Using the following property of thedeterminant of a matrix: det( − iα ˆ σ Λ) = exp[tr log( − iα ˆ σ Λ)] = exp (cid:34) − tr ∞ (cid:88) m =1 m (2 iα ˆ σ Λ) m (cid:35) , (E16)and also the geometric series [( − iα ˆ σ Λ) − − ] = ∞ (cid:88) m =1 (2 iα ˆ σ Λ) m , (E17)we find I ( λ , . . . , λ n ) = exp (cid:20)
12 tr ∞ (cid:88) m =1 m (2 iα ˆ σ Λ) m − n (cid:88) j = k +1 λ j σ T j ˆ σ − ∞ (cid:88) m =1 (2 iα ˆ σ Λ) m n (cid:88) l = k +1 λ l σ l − n (cid:88) i,j = k +1 σ ij λ i λ j (cid:21) . Thus, the cumulant is C n ( α R , . . . , α R k , R k +1 . . . , R n ) = ( − i ) n ∂∂λ . . . ∂∂λ n log I ( λ , . . . , λ n ) (cid:12)(cid:12)(cid:12)(cid:12) λ =0 (E18) = ( − i ) n ∂∂λ . . . ∂∂λ n (cid:20)
12 tr ∞ (cid:88) m =1 m (2 iα ˆ σ Λ) m − n (cid:88) j = k +1 λ j σ T j ˆ σ − ∞ (cid:88) m =1 (2 iα ˆ σ Λ) m n (cid:88) l = k +1 λ l σ l − n (cid:88) i,j = k +1 σ ij λ i λ j (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) λ =0 . A term is non zero only if it contains all the λ i ’s. For k < n the first term vanishes because it does not contain the λ i with i > k . Instead the last term is always vanishing for n > , because the number of derivatives is higher than thenumber of λ i ’s present. So, until now C n ( α R , . . . , α R k , R k +1 . . . , R n ) = ( − i ) n ∂∂λ . . . ∂∂λ n − n (cid:88) j = k +1 λ j σ T j ˆ σ − ∞ (cid:88) m =1 (2 iα ˆ σ Λ) m n (cid:88) l = k +1 λ l σ l (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) λ =0 . (E19)Now, let us make the derivative with respect to λ n C n ( α R , . . . , α R k , R k +1 . . . , R n ) = (E20) ( − i ) n ∂∂λ . . . ∂∂λ n − (cid:20) − n (cid:88) j = k +1 δ jn σ T j ˆ σ − ∞ (cid:88) m =1 (2 iα ˆ σ Λ) m n (cid:88) l = k +1 λ l σ l − n (cid:88) j = k +1 λ j σ T j ˆ σ − ∞ (cid:88) m =1 (2 iα ˆ σ Λ) m n (cid:88) l = k +1 δ ln σ l (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) λ =0 . ˆ σ − (ˆ σ Λ) m +1 = Λ(ˆ σ Λ) m and this matrix is symmetric: [Λ(ˆ σ Λ) m ] T = [(ˆ σ Λ) m ] T Λ T = [(ˆ σ Λ) m ] T Λ = (cid:2) (ˆ σ Λ) T (cid:3) m Λ = (cid:0) Λ T ˆ σ T (cid:1) m Λ = (Λˆ σ ) m Λ = Λ(ˆ σ Λ) m ; (E21)so, deriving with respect to λ n − C n ( α R , . . . , α R k , R k +1 . . . , R n ) = ( − i ) n ∂∂λ . . . ∂∂λ n − (cid:20) − σ T n Λ ∞ (cid:88) m =1 (2 iα ˆ σ Λ) m − n (cid:88) l = k +1 λ l σ l (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) λ =0 = ( − i ) n ∂∂λ . . . ∂∂λ n − (cid:20) − σ T n Λ ∞ (cid:88) m =1 (2 iα ˆ σ Λ) m − σ n − (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) λ =0 . (E22)It is clear that, if k < n − there are still the derivatives ∂∂λ k +1 . . . ∂∂λ n − to apply, but no λ k +1 , . . . , λ n − to meet.Therefore: C n ( α R , . . . , α R k , R k +1 , . . . , R n ) = 0 , ∀ k < n − (E23)Now let us move to compute the only three orders left.
2. Computation of the O ( α n − ) For the leading part of the cumulant, it is enough to set k = n − in the eq. (E22) (let us consider the case n > , sowe can get rid of the last term that we canceled in eq. E18): C n ( α R , . . . , α R n − , R n − , R n ) = ( − i ) n ∂∂λ . . . ∂∂λ n − (cid:20) − σ T n Λ ∞ (cid:88) m =1 (2 iα ˆ σ Λ) m − σ n − (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) λ =0 , (E24)where now Λ = diag ( λ , . . . λ n − ) .The first observation we can do is that only the O ( α n − ) matters in the sum over m . In fact, for m < n − thederivatives are too much and they annihilate the term. Instead, for m > n − the derivatives are not enough and whenwe put λ = 0 the term vanishes. C n ( α R , . . . , α R n − , R n − , R n ) = (2 α ) n − ∂∂λ . . . ∂∂λ n − (cid:20) σ T n − Λ(ˆ σ Λ) n − σ n (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) λ =0 . (E25)Now we can see that: σ T n − (cid:2) Λ(ˆ σ Λ) n − (cid:3) σ n = n − (cid:88) i,j =1 σ n − ,i n − (cid:88) i ,...,i n − =1 λ i σ ii λ i σ i i λ i . . . λ i n − σ i n − j λ j σ jn . (E26)Of course, because of the presence of the derivatives, in the sum, only the terms with all the { λ i } n − i =1 different matters,i.e. we only keep the permutations of { , , . . . , n − } . So, putting all together C n [ α R , . . . , α R n − , R n − , R n ] = (2 α ) n − ∂∂λ . . . ∂∂λ n − n − (cid:88) i ,...,i n − =1 λ i . . . λ i n − σ n − ,i σ i i . . . σ i n − i n − σ i n − n (cid:12)(cid:12)(cid:12)(cid:12) λ =0 = (2 α ) n − (cid:88) { i ,...,i n − } ==perms { ,...,n − } σ n − ,i σ i i . . . σ i n − i n − σ i n − n . (E27)So, the formula for the O ( α n − ) contribute to the n -th order cumulant is: C n ( α R , . . . , α R n − , R n − , R n ) = (2 α ) n − (cid:88) { i ,...,i n − } ==perms { ,...,n − } (cid:104)R n − R i (cid:105)(cid:104)R i R i (cid:105) . . . (cid:104)R i n − R i n − (cid:105)(cid:104)R i n − R n (cid:105) (E28)2For example, we can check this formula for n = 4 C ( α R , α R , R , R ) = (2 α ) (cid:88) { i,j } =perms { , } (cid:104)R R i (cid:105)(cid:104)R i R j (cid:105)(cid:104)R j R (cid:105) (E29) = 4 α ( (cid:104)R R (cid:105)(cid:104)R R (cid:105)(cid:104)R R (cid:105) + (cid:104)R R (cid:105)(cid:104)R R (cid:105)(cid:104)R R (cid:105) ) = 4 α (cid:104)R R (cid:105) ( (cid:104)R R (cid:105)(cid:104)R R (cid:105) + (cid:104)R R (cid:105)(cid:104)R R (cid:105) ) , that precisely matches the result obtained in eq. (D52).
3. Computation of the O ( α n − ) For k = n − we have: I ( λ , . . . , λ n ) = 1 (cid:112) det( − iα ˆ σ Λ) exp (cid:26) − λ n σ T n ˆ σ − (cid:2) ( − iα ˆ σ Λ) − − (cid:3) σ n − σ nn λ n (cid:27) , (E30)and so: C n ( α R , . . . , α R n − , R n ) = ( − i ) n ∂∂λ . . . ∂∂λ n log I ( λ , . . . , λ n ) (cid:12)(cid:12)(cid:12)(cid:12) λ =0 = (E31) = ( − i ) n ∂∂λ . . . ∂∂λ n (cid:20)
12 tr ∞ (cid:88) m =1 m (2 iα ˆ σ Λ) m − λ n σ T n ˆ σ − ∞ (cid:88) m =1 (2 iα ˆ σ Λ) m σ n − σ nn λ n (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) λ =0 = (E32) = ( − i ) n ∂∂λ . . . ∂∂λ n − (cid:34) − λ n (cid:32) σ T n ˆ σ − ∞ (cid:88) m =1 (2 iα ˆ σ Λ) m σ n + σ nn (cid:33)(cid:35)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) λ =0 . (E33)In the last step, everything vanishes since it remains a λ n with no more derivatives with respect to it and in the end wemust set λ = 0 , so there is no O ( α n − ) contribution: C n ( α R , . . . , α R n − , R n ) = 0 (E34)
4. Computation of the O ( α n ) Now we have
Λ = diag ( λ , . . . , λ n ) and the trace part is the only term left, so: I ( λ , . . . , λ n ) = 1 (cid:112) (2 π ) n det σ (cid:90) d n R e − R T ( σ − − iα Λ ) R = exp (cid:20) −
12 tr log ( − iασ Λ) (cid:21) , (E35)where now σ = ( σ ij ) ni,j =1 and Λ = diag ( λ , . . . , λ n ) . With the usual steps, knowing that only the O ( α n ) matters: C n ( α R , . . . , α R n ) = ( − i ) n ∂∂λ . . . ∂∂λ n (cid:34)
12 tr ∞ (cid:88) m =0 m (2 iασ Λ) m (cid:35)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) λ =0 = 2 n − α n n ∂∂λ . . . ∂∂λ n tr( σ Λ) n (cid:12)(cid:12)(cid:12)(cid:12) λ =0 . (E36)The trace is: tr( σ Λ) n = n (cid:88) i ,...,i n =1 λ i λ i . . . λ i n σ i n i σ i i . . . σ i n − i n , (E37)and after the derivation only the permutations of { , , . . . , n } survive, so: C n ( α R , . . . , α R n ) = 2 n − α n n (cid:88) { i ,...,i n } ==perms { ,...,n } (cid:104)R i n R i (cid:105) . . . (cid:104)R i n − R i n (cid:105) (E38)3
5. Computation of the generic cumulant
Now that we have all the ingredients, we can compute the generic n -th order cumulant up to all orders. First of all,let us focus on the O ( α n − ) contribute: C O ( α n − ) n ( h , . . . , h n ) = (2 α ) n − (cid:88) { i,j }∈ A (cid:88) { i ,...,i n − } ==perms { ,...,n }\{ i,j } (cid:104)R i R i (cid:105)(cid:104)R i R i (cid:105) . . . (cid:104)R i n − R j (cid:105) (E39)where A is the set in which all the non-ordered couples of numbers between { , . . . , n } are contained. E.g. for n = 4 , A = (cid:8) { , } , { , } , { , } , { , } , { , } , { , } (cid:9) . From combinatorics we know that the permutations of n − differentobjects are ( n − , while the cardinality of A is the number of ways in which we can take 2 objects among n differentobjects without taking care of their order, which is (cid:0) n (cid:1) = n ( n − . So, putting these 2 together, the total sum is composedof n !2 elements. If we compute also the exchanged couples, we double the sum, having n ! elements, that is the cardinalityof perms { , . . . , n } . So we can put the two sums together dividing by 2: C O ( α n − ) n ( h , . . . , h n ) = 2 n − α n − (cid:88) { i ,...,i n } ==perms { ,...,n } (cid:104)R i R i (cid:105)(cid:104)R i R i (cid:105) . . . (cid:104)R i n − R i n (cid:105) (E40)Now, let us add also the O ( α n ) contribution, that we can write as C n ( α R , . . . , α R n ) = 2 n − α n n (cid:88) { i ,...,i n } ==perms { ,...,n } (cid:104)R i R i (cid:105) . . . (cid:104)R i n R i (cid:105) (E41)And so, summing the two C n ( h , . . . , h n ) = 2 n − α n − (cid:88) { i ,...,i n } ==perms { ,...,n } (cid:104)R i R i (cid:105) . . . (cid:104)R i n − R i n (cid:105) (cid:18) α n (cid:104)R i R i n (cid:105) (cid:19) (E42)And we can easily see that, putting (cid:126)x = · · · = (cid:126)x n : C n ( h, . . . , h ) = 2 n − α n − n ! (cid:104)R (cid:105) n − (cid:18) α n (cid:104)R (cid:105) (cid:19) = 2 n − ( n − σ (cid:0) ασ (cid:1) n − (cid:0) n + 4 α σ (cid:1) , (E43)and this result exactly reproduces the formula given in eq. (D55). Appendix F: Threshold for gravitational collapse into black holes
We follow the approach of refs. [24, 25] that we generalize to include local non-gaussianities in the curvature pertur-bation field. The main steps of the computation are the following i) Consider eq. (4) at the linear order in h δ ( (cid:126)x, t ) = − (cid:18) aH (cid:19) (cid:52) h ( (cid:126)x ) . (F1)As usual, δ ( (cid:126)x, t ) is a statistical variable and only statistical averages can be used for practical purposes. Insteadof using δ ( (cid:126)x, t ) , one then defines the averaged density radial profile ¯ δ ( r, t ) by means of ¯ δ ( r, t ) ≡ F ( aH ) ψ ( r ) , with ψ ( r ) = (cid:104) δ ( (cid:126)x, t ) δ ( (cid:126)x + (cid:126)r, t ) (cid:105)(cid:104) δ ( (cid:126)x, t ) δ ( (cid:126)x, t ) (cid:105) , (F2) Recently, ref. [63] proposed a simplified analytical prescription for the computation of the threshold for PBH formation. We note that in ref. [63]the quantity δ c corresponds to δ th in our notation. (cid:104) δ ( (cid:126)x, t ) δ ( (cid:126)x + (cid:126)r, t ) (cid:105) = (F3) aH ) (cid:90) dk k sin( kr ) kr (cid:20) P R ( k ) + α k (cid:90) dqq d cos θ ( k + q − kq cos θ ) / P R ( q ) P R ( (cid:112) k + q − kq cos θ ) (cid:124) (cid:123)(cid:122) (cid:125) ≡P NG ( k ) (cid:21) , while (cid:104) δ ( (cid:126)x, t ) δ ( (cid:126)x, t ) (cid:105) is given by the limit lim r → (cid:104) δ ( (cid:126)x, t ) δ ( (cid:126)x + (cid:126)r, t ) (cid:105) . Eq. (F3) follows from the explicit evaluation(going first in Fourier space) of the two-point correlator; P NG ( k ) represents the correction to the averaged densityprofile due to local non-gaussianities, and the integral over q arises because of the two convolutions that areneeded to compute the term α (cid:104)R ( (cid:126)x ) R ( (cid:126)x + (cid:126)r ) (cid:105) . Remember indeed that the Fourier transform of R is not thesquare of the Fourier transform of R but it involves a convolution.Notice that ψ ( r ) in eq. (F2) is time-independent (assuming of course the perturbations to be fully super-horizonthus constant in time), and the only time-dependence of ¯ δ ( r, t ) comes from the overall factor / ( aH ) . Further-more, spherical symmetry is assumed. This means that we are neglecting interactions between adjacent over-densities [21]. Intuitively, this is a very good assumption since the formation of a PBH is already a rare event,and having two or more black holes forming on the same site should be extremely unlikely. In eq. (F2) F is adimensionful parameter which is related to the amplitude of the over-density at the center of the radial profile.Finally, notice that, by construction, we have ψ (0) = 1 and, consequently, ¯ δ (0 , t ) = F / ( aH ) . ii) The next step consists in the evaluation of the so-called compaction function C ( r ) . This is the most importantquantity because, as explained in refs. [24, 25], PBHs form when local maxima of the compaction function exceeda certain threshold value [64, 65]. The main point of the analysis is the identification of such threshold.Physically, the compaction function is defined as twice the local excess-mass over the comoving areal radius, thatis (we use explicitly, for the sake of clarity, the reduced Planck mass ¯ M = 1 / πG N ) C ( r, t ) ≡ G N δM ( r, t ) R ( r, t ) , (F4)where the areal radius is R ( r, t ) ≡ a ( t ) r and the local excess-mass due to the overdensity field δ ( r, t ) within aspherical region of radius R is given by δM ( r, t ) M b ( r, t ) ≡ V b ( r, t ) (cid:90) R πδ ( r (cid:48) , t ) R (cid:48) dR (cid:48) , (F5)where V b ( r, t ) ≡ (4 π/ R ( r, t ) , the background radiation energy density is ρ b ( t ) = 3 ¯ M H ( t ) and M b ( r, t ) = V b ( r, t ) ρ b ( t ) . Using these relation, eq. (F4) can be recast in the form C ( r, t ) = 3( aH ) r (cid:90) r dr (cid:48) r (cid:48) δ ( r (cid:48) , t ) . (F6)This is a generic expression. If we now use for the overdensity field the averaged density radial profile ¯ δ ( r, t ) defined in eq. (F2), we find C ( r ) (cid:39) F r (cid:90) r dr (cid:48) r (cid:48) ψ ( r (cid:48) ) , (F7)which is valid in the so-called long-wavelength approximation of ref. [25]. Notice that in this approximation C ( r ) does not depend on time. From C ( r ) , one computes r m which is the position of its maximum. The precisecomputation of r m for a given power spectrum of curvature perturbation will be discussed below. iii) After computing r m , we use the fact that, as shown in ref. [25], we have δ th = 3¯ δ ( r m , t m ) where δ th is definedas the threshold value for the local excess-mass (defined in eq. (F5)) associated to the the averaged curvatureperturbations ¯ δ ( r, t ) and evaluated at position r m and time t m . At super-horizon scales, the same computationthat led to eq. (F7) gives δM ( r, t ) M b ( r, t ) = 1 V b ( r, t ) (cid:90) R π ¯ δ ( r (cid:48) , t ) R (cid:48) dR (cid:48) (cid:39) (cid:18) aHr (cid:19) C ( r ) , (F8)and, known r m , the time t m is defined implicitly by a ( t m ) H ( t m ) r m = 1 . From the previous definition, it followsthat δ th = δM ( r m , t m ) /M b ( r m , t m ) (cid:39) C c ( r m ) .5If we knew δ th it would be possible to use δ th = 3¯ δ ( r m , t m ) and eq. (F2) to extract the critical value of F , named F c in the following. Correspondingly, we indicate with C c ( r ) the critical compaction function, that is eq. (F7)with F = F c .The quantity δ th is usually extracted from simulations in numerical relativity, and has been found to range in theinterval . (cid:46) δ th (cid:46) / . The left-side of this interval corresponds to the so-called Harada-Yoo-Kohri limit thatis the threshold for which a very sharply peaked over-density profile collapses into a zero-mass black hole [66]. iv) To proceed further, the simplest way to go is to fix the value of δ th within the range . (cid:46) δ th (cid:46) / and extract F c from δ th = 3¯ δ ( r m , t m ) . We can actually do better and exploit the result of ref. [67]. By using the parameter q defined as q ≡ −C (cid:48)(cid:48) ( r m ) r m / C ( r m ) , the value of δ th can be analytically computed as δ th = 415 e − /q (cid:20) q − / q Γ(5 / q ) − Γ(5 / q, /q ) (cid:21) . (F9)We refer to ref. [67] for a detailed discussion. Eq. (F9) is an extremely neat result and shows that the threshold forthe compaction function is only sensitive to its curvature at the maximum.The strategy, therefore, is the following. We compute q from C ( r ) and r m , then δ th from eq. (F9). From δ th , using δ th = 3¯ δ ( r m , t m ) , we can solve for F c . From eq. (F2), we indeed have F c = δ th / r m ψ ( r m ) . v) Finally, we are in the position to extract the value of δ c that enters in eq.(6), leading to eq. (12) and eq. (B19).We simply have δ c = δ th ψ ( r m ) = F c r m , (F10)and, from eq. (F2), we see that δ c corresponds to ¯ δ c (0 , t m ) that is the value of the critical averaged density pro-file at the center. Furthermore, in eq. (B19) we evaluated the comoving horizon length /aH at time t m so that ( a m H m ) = 1 /r m . We remark that the threshold value for the overdensity field δ is not set by δ th (that is thethreshold on the compaction function) but rather by F c which is the threshold at the center of the averagedover-density. This point will become more clear in our numerical analysis.We apply this procedure to the case relevant for our analysis.First of all, let us comment about the non-gaussian correction in eq. (F3). For illustration, we consider the powerspectrum in eq. (B24). The advantage is that we can do the angular integration that appears in P NG ( k ) semi-analytically. We find P NG ( k ) = α k (cid:90) dqq A g √ πvkk (cid:63) exp (cid:32) v − log qk (cid:63) v (cid:33) (cid:34) Erf (cid:32) v + log k + qk (cid:63) √ v (cid:33) − Erf (cid:32) v + log | k − q | k (cid:63) √ v (cid:33)(cid:35) , (F11)where − (cid:54) Erf( x ) (cid:54) is the error function. We can now compare P NG ( k ) with the leading gaussian term P R ( k ) ineq. (B24). We show this comparison in fig. 13 for a narrow ( v = 0 . ) and a broad ( v = 0 . ) spectrum. The outcomeof this comparison is that the presence of the non-gaussian correction P NG ( k ) can, in principle, modify the shape ofthe leading order power spectrum. For instance, we note that in the expression for P NG ( k ) the peak in k is shiftedtowards slightly higher k compared to k (cid:63) as a consequence of the convolution. However, the correction P NG ( k ) isnaturally suppressed, in amplitude, since proportional to A g with A g (cid:28) . This means that, for P NG ( k ) to be relevant,one needs to overcome this suppression by taking α (cid:29) , and this is not what typically happens in the class of modelsthat we have in mind for our analysis (where α (cid:46) ). For this reason, in the following we will neglect the non-gaussiancorrection P NG ( k ) in eq. (F3). Furthermore, since the computation of ψ ( r ) in eq. (F2) will be in any case numerical, weput aside the toy-model in eq. (B24) and move to consider the more realistic case of the power spectrum introduced ineq. (11).As already anticipated in section III, in realistic situations in which the power spectrum is not sharply peaked thecomputation of δ c requires some care. We start by introducing the function P cut R ( k ) ≡ P R ( k ) e − k /k with some cut-off wavenumber k cut . The reason is that we are interested in the local maximum of the compaction function for which This can be easily understood if we look at eq. (F11). The dominant contribution to the q -integral comes from the region where the argument ofthe overall exponential vanishes (for q (cid:29) k (cid:63) and q (cid:28) k (cid:63) we have an exponential suppression), that is for q (cid:39) e v k (cid:63) . If we select this contributionfrom the q -integral, it is easy to see that the resulting expression P NG ( k ) is maximized when k (cid:39) k (cid:63) e v . This is because when q = e v k (cid:63) thesecond error function in eq. (F11) takes the value − Erf[ v + log( k − k (cid:63) e v k (cid:63) )] which is maximized to when k = k (cid:63) e v , that is − Erf[ −∞ ] = 1 . ����� - � ��� - � �� - � �� - � �� - � �� - � �� � / ★ � ( ) / ℛ ( ★ ) � = ℛ � = �� ( α = � ) = ��� ����� - � ��� - � �� - � �� - � �� - � �� - � �� � / ★ � ( ) / ℛ ( ★ ) � = ℛ � = �� ( α = � ) = ��� FIG. 13: P R ( k ) and P NG ( k ) for the power spectrum in eq. (B24) as a function of the comoving wavenumber k in units of k (cid:63) . Wenormalize the y -axis with respect to P R ( k (cid:63) ) . We set α = 1 while A g is chosen in such a way that P R ( k (cid:63) ) = 10 − . � � � � ��������������������� ��� / ★ � � � / σ � α = � � � α = � � � α = � � � � FIG. 14:
Ratio s min /σ defined in eq. (B19) as a function of the cutoff wavenumber k cut in units of k (cid:63) . More precisely, s min /σ is computed following the procedure sketched in i)-v) in appendix F and, for a given k cut , we used the power spectrum P cut R ( k ) ≡P R ( k ) e − k /k with P R ( k ) given in eq. (11) for three different benchmark values of α . Notice that σ ∝ P R ( k (cid:63) ) ≡ p (cid:63) and, conse-quently, for the ratio on the y -axis we have s min /σ ∝ / √ p (cid:63) . We show the value of s min /σ in units of / √ p (cid:63) meaning that for agiven peak amplitude P R ( k (cid:63) ) and a given k cut the corresponding threshold for collapse s min /σ is given by the value on the y -axistimes / √ p (cid:63) . Since P R ( k (cid:63) ) (cid:28) , we have s min /σ (cid:29) . the threshold for collapse in eq. (12) and eq. (B19) is minimized, thus giving the largest probability for collapse. We canlook for this maximum by means of a numerical scan over the value k cut . We show our numerical results in fig. 14 (seecaption for details). We find that the value of s min /σ is minimized at around k cut /k (cid:63) (cid:39) . . The precise value dependson the value of α . This is because, as explained below eq. (11), the spectral index of the power-law falloff of the powerspectrum after the peak is controlled by α , and different values of the latter give different shapes.In the left panel of fig. 15, we show the critical compaction function (which is defined by eq. (F7) with F = F c ) forthe case with k cut = 1 . k (cid:63) and α = 0 . . In the right panel of the same figure, we show the radial profile of the criticalaveraged over-density which is given by eq. (F2) with F = F c at time t = t m . This figure makes clear the differencebetween δ th and δ c = F c r m in eq. (F10). The former is the threshold that refers to the value of the compaction functionat its maximum while the latter refers to the over-density at the center.After clarifying the computation of the threshold for gravitational collapse, let us now turn to the computation of the7 � � � � � ����������������� � / � � � ( � ) δ �� ≃ ���� � � � � � ����������������������� � / � � δ ( � � ) δ �� / � δ � = ℱ � � � FIG. 15:
Left panel. Critical compaction function (that is eq. (F7) evaluated for F = F c ) obtained for P cut R ( k ) with k cut = 1 . k (cid:63) and α = 0 . ; δ th (cid:39) . is the threshold computed with eq. (F9), and corresponds to the maximum of the critical compaction function.Right panel. For the same value of k cut and α = 0 . , we show the radial profile of the critical averaged over-density. PBH abundance. At the time t f of their formation the latter is defined by β = 1 ρ b ( t f ) (cid:90) ∞ ν c dνρ PBH ( ν ) , (F12)and it simply corresponds to the ratio between the mass density in PBHs and the background energy density at for-mation time, ρ b ( t f ) = 3 H f / π . In eq. (F12) we use ν ≡ δ/σ δ = 2 s/σ . The parameter ν controls the hight of theoverdensity peak that collapses into a PBH, and the integration in eq. (F12) tells us that only peaks above the thresholdcontribute to the PBH abundance (while perturbations with ν < ν c disperse into the expanding Universe). The massdensity of PBHs takes the form ρ PBH ( ν ) = M PBH ( ν ) n physmax ( ν ) , (F13)where M PBH ( ν ) is the PBH mass (usually calculated at horizon crossing time t m ) and n physmax ( ν ) is the number density ofpeaks of height ν in the physical space at formation time. The latter is related to the comoving number density of peaks n max ( ν ) by the scaling (in three spatial dimensions) n physmax ( ν ) = n max ( ν ) /a f . We consider for the moment the realisticcase of three spatial dimensions. The PBH mass M PBH ( ν ) is proportional to the horizon mass at horizon-crossing time, M H ( t m ) , according to the relation [24] M PBH ( ν ) = K M H ( t m ) (cid:18) ¯ σ a H a m H m (cid:19) ˜ γ ( ν − ν c ) ˜ γ . (F14)In this equation the horizon mass at horizon-crossing time is given by the equation M H ( t m ) = 1 / H m while the factor ( ν − ν c ) γ is the scaling law for critical collapse (extracted from numerical simulation, and with ˜ γ (cid:39) . in the case ofradiation [68]) while K = O (1) . In our notation, ν = 2 s/σ and ν c = 2 s min /σ . Notice that in eq. (F14) we introduce,following ref. [24], the spectral moments referred to the density perturbation ¯ σ j = 1681 1( aH ) (cid:90) dkk P R ( k ) k j +4 . (F15)All in all, eq. (F12) translates into β = 4 π K (cid:18) ¯ σ a H a m H m (cid:19) ˜ γ (cid:18) a m H m (cid:19) a f a m (cid:90) ∞ ν c dν ( ν − ν c ) ˜ γ n max ( ν ) . (F16)8We need to specify the peak number density and integrate. Let us consider first the gaussian case, in which wehave [21] n max ( ν ) = 14 π (cid:18) γR ∗ (cid:19) ν exp( − ν / , with γ ≡ ¯ σ ¯ σ ¯ σ , R ∗ ≡ √ σ ¯ σ . (F17)We, therefore, have β = K π (cid:18) ¯ σ a H a m H m (cid:19) ˜ γ (cid:18) a f a m (cid:19) (cid:18) γr m R ∗ (cid:19) ν γc (cid:90) ∞ dx ( x − ˜ γ x e − x ν c / . (F18)Ref. [24] computed the above integral by means of a saddle point approximation. We, instead, note that it admits thefollowing exact expression in terms of generalized hypergeometric functions (cid:90) ∞ dx ( x − ˜ γ x e − x ν c / = 2 ˜ γ/ ν ˜ γ +4 c (cid:26) γ (cid:18) γ (cid:19) F (cid:20)(cid:26)
12 (1 − ˜ γ ) , − ˜ γ (cid:27) , (cid:26) , − − ˜ γ (cid:27) , − ν c (cid:21) (F19) − √ ν c ˜ γ (cid:18) γ (cid:19) F (cid:20)(cid:26)
12 (1 − ˜ γ ) , − ˜ γ (cid:27) , (cid:26) , − (cid:18) − ˜ γ (cid:19)(cid:27) , − ν c (cid:21) (cid:27) . This result, in turn, allows us to extract the simple approximation (cid:90) ∞ dx ( x − ˜ γ x e − x ν c / (cid:39) e − ν c / ν − γ ) c Γ(1 + ˜ γ ) . (F20)which is valid for ν c / ˜ γ (cid:29) . We checked numerically that this approximation works exquisitely well (at the % level for ν c = 10 with increasing precision for larger thresholds). We find β (cid:39) (cid:34) K π (cid:18) ¯ σ a H a m H m ν c (cid:19) ˜ γ (cid:18) a f a m (cid:19) (cid:18) γr m R ∗ (cid:19) Γ(1 + ˜ γ ) (cid:35)(cid:124) (cid:123)(cid:122) (cid:125) ∼ O (1) ν c e − ν c / , �� - � �� - �� �� - �� �� - �� �� - � ℛ ( ★ ) β (F21)The numerical value of β is controlled by the exponential function e − ν c / while the pre-factor inside the square bracketsin eq. (F21) is a dimensionless O (1) number whose exact value only plays a sub-leading role in the determination of β . To corroborate this statement, in the figure attached to eq. (F21) we compare (as function of the peak amplitude ofthe power spectrum) the abundance β computed with (black solid line) and without (red dashed line) the pre-factorin the square brackets, and we show that the two computations almost coincide.The computation in eq. (F21) shows that in order to estimate the value of β in a pragmatic way we can just integratethe comoving number density of maxima of the overdensity field and make it dimensionless β (cid:39) π R ∗ (cid:90) ∞ ν c dν n max ( ν ) . (F22)Eq. (F22) gives a perfect approximation of eq. (F21).In presence of local non-gaussianities the abundance can be obtained following the same logic than leads toeq. (F16), but computing the number density of peaks as explained in sec. C, namely identifying this quantity withthe number density of spiky maxima of the comoving density perturbations, the sum of eq. (C15) and eq. (C18). Weobtain β = 4 π K (cid:18) ¯ σ a H a m H m (cid:19) ˜ γ (cid:18) a m H m (cid:19) a f a m × (cid:34)(cid:90) ∞− ασ d ¯ ν (cid:90) ∞ x δ (¯ ν ) dx [˜ ν (¯ ν, x ) − ν c ] ˜ γ ¯ n max (¯ ν, x ) + (cid:90) ασ −∞ d ¯ ν (cid:90) x δ (¯ ν ) −∞ dx [˜ ν (¯ ν, x ) − ν c ] ˜ γ ¯ n max (¯ ν, x ) (cid:35) , (F23) Notice that R ∗ is now defined in terms of the spectral moments ¯ σ j , differently from what done in section B where we have an analogous expressionin terms of σ j . ¯ n max (¯ ν, x ) is defined in eqs. (C16, C17). Notice that ¯ ν = R /σ , while in eq. (F14) we have ν = δ/σ δ . In theexpression above we have written the latter quantity as ν = ˜ ν (¯ ν, x ) = x (1 + 2 ασ ¯ ν ) . Once again, we find that theabundance is well approximated by a simpler expression, obtained by taking the number density of peaks, N (I)pk + N (II)pk (eqs. (C15, C18) ), and making it dimensionless β (cid:39) π R ∗ (cid:104) N (I)pk + N (II)pk (cid:105) . (F24)This equation reproduces eq. (8).Having looked at the exact definition of β , let us now consider our simplified two-dimensional model. We followthe rationale that led to eq. (F22). The number density of maxima above threshold of the overdensity field is given by N max ( s min ) (that is, for instance, eq. (B33) in the gaussian case) and, in two spatial dimensions, to make it dimension-less we need to multiply it times R ∗ β (cid:39) R ∗ N max ( s min ) . (F25) Appendix G: Non-gaussianities from non-linearities
Consider the relation between the overdensity field and the primordial curvature perturbation δ ( (cid:126)x, t ) = − (cid:18) aH (cid:19) e − h ( (cid:126)x ) (cid:20) (cid:52) h ( (cid:126)x ) + 12 h i ( (cid:126)x ) h i ( (cid:126)x ) (cid:21) ≡ − (cid:18) aH (cid:19) δ r ( (cid:126)x ) . (G1)We have two kinds of non-gaussianities. On the one hand, non-gaussianities of primordial origin—which are themain subject of this paper—that arise from the non-gaussian nature of the random field h . On the other one, non-gaussianities also arise from the presence of non-linear terms in the relation between h and δ in eq. (G1). The lattertype of non-gaussianities is ineludible in the sense that the overdensity field is non-gaussian even if the curvatureperturbation is gaussian [26]. However, the presence of non-gaussianities of primordial origin is also expected in thecontext of realistic models of single-field inflation which lead to black hole formation; this is because in these modelsslow-roll conditions (which usually suppress the amount of non-gaussianities) are badly violated. A comprehensivestudy in which both effects are included is needed (see ref. [22] for an attempt in this direction). In this appendix, welimit our discussion to a number of simple qualitative considerations. Furthermore, as already stated in section III, wedo not consider the effect of non-linearities in eq. (G1) on the shape of the peak of the overdensity field (in principle,this may affect the way in which δ c is computed; as explained in appendix F, we compute δ c by taking the linear ap-proximation in eq. (G1), and we use the average density profile to describe the shape of the peaks). In the following, wewill only discuss how non-linearities in eq. (G1) are expected to change the probability distribution of the peaks of theoverdensity field. a. The case without primordial non-gaussianities Let us start from the gaussian case in which h = R (that is α = 0 ). Compared to the analysis carried out in ap-pendix B, the first thing to check is whether peaks of the overdensity field are still identifiable with local maxima of R even in the presence of non-linearities. To answer this question, we repeat the numerical analysis that led to fig. 7, thistime using the full relation in eq. (G1). A positive answer is found, in agreement with the result of ref. [26]. On this basis,we can identify (cid:126)y pk (the position of the peak of the overdensity field) with (cid:126)x M (the position of the associated local max-imum of R ); consequently, the condition that the peak amplitude of the overdensity must be larger than some criticalvalue δ c becomes −(cid:52)R ( (cid:126)x M ) (cid:38)
94 ( aH ) e R M δ c , (G2)with δ c that is now “renormalized” by the factor e R M . The bottom line is the following. If the value of R M is so largethat e R M (cid:29) then it may sizably change the threshold condition on −(cid:52)R ( (cid:126)x M ) . This point was already emphasizedin ref. [26].In addition, the important remark that we want to make is the following. At the quantitative level, the importanceof this effect strongly depends on the power spectrum. If the power spectrum is very peaked R M and −(cid:52)R ( (cid:126)x M ) are In eqs. (F24, F25) R ∗ is defined as in section B, i.e. R ∗ ≡ √ d σ /σ . −(cid:52)R ( x M ) are likely to be also regions with large R M . In turn, this enhances (because of thefactor e R M ) the threshold condition on the curvature −(cid:52)R ( (cid:126)x M ) compared to the linear case, and we expect a smallerblack hole abundance since higher threshold for black hole formation means rarer objects. The net result is that, inorder to keep the value of β constant, one needs to take larger values of P R ( k (cid:63) ) to compensate the above effect. Thisexpectation is confirmed in ref. [26] for the case of a very narrow log-normal power spectrum. If we take eq. (B24) with v = 0 . , in particular, ref. [26] finds that one needs to increase P R ( k (cid:63) ) by more than one order of magnitude.However, in single-field models of inflation that are relevant in the present work the power spectrum is broader thana very narrow log-normal function. This has an important impact. As argued in appendix B (see fig. 10 and related dis-cussion), for a broad power spectrum the correlation between R M and −(cid:52)R ( x M ) is much less pronounced (becauseit is controlled by the parameter < γ ≡ σ /σ σ < with γ → for a narrow power spectrum and γ → fora broad one, see eq. (B22)). We expect, therefore, that in the realistic case the effect of non-linearities will be far lessrelevant if compared to what expected based on the case of a narrow power spectrum: in the realistic case large valuesof −(cid:52)R ( x M ) are not necessarily combined with large values of R M , and the enhancement effect in eq. (G2) will be lessimportant. This intuition is confirmed in the left panel of fig. 6. The dot-dashed black line (with label “NL”) representsthe computation of β in three spatial dimensions in the gaussian limit but including the condition in eq. (G2). In theanalyzed case, we find γ (cid:39) . which is significantly far from the narrow limit γ → . Non-linearities suppress, asexpected, the PBH abundance, but the peak amplitude of the power spectrum which is required to get the referencevalue β (cid:39) − increases only by a factor smaller than 2. b. The case with primordial non-gaussianities We now consider eq. (G1) with h ( (cid:126)x ) = R ( (cid:126)x ) + α [ R ( (cid:126)x ) − σ ] . Eq. (G1) takes the form −(cid:52)R ( (cid:126)x M ) (cid:38)
94 ( aH ) (1 + 2 α R M ) exp (cid:8) R M + α ( R − σ )] (cid:9) δ c . (G3)The number density of peaks of the overdensity field can be computed following the same procedure outlined in sec-tion C. The only difference is that we should now impose eq. (G3) instead of eq. (6). This means that eq. (C14) is replacedby the following expression: x > a m H m ) σ δ c ασ ¯ ν exp (cid:8) σ [¯ ν + ασ (¯ ν − (cid:9) ≡ x NL δ (¯ ν ) , (G4)and the function x NL δ (¯ ν ) replaces x δ (¯ ν ) in the integrals in eqs. (C15, C18). We show our results for the PBH abundancein the left panel of fig. 6. We find that the additional presence of primordial non-gaussianities makes the formation ofPBHs even harder than the expectation based solely on the presence of non-linearities. However, the enhancementin the value of P R ( k (cid:63) ) that is required (compared to the gaussian case) to fit the reference value β = 10 − remainsextremely small (less than a factor of even for α = 0 . ). [1] S. Hawking, “ Gravitationally collapsed objects of very low mass, ” Mon. Not. Roy. Astron. Soc. (1971), 75[2] B. J. Carr and S. Hawking, “
Black holes in the early Universe, ” Mon. Not. Roy. Astron. Soc. (1974), 399-415[3] 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]].[4] 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]].[5] A. Katz, J. Kopp, S. Sibiryakov and W. Xue, “ Femtolensing by Dark Matter Revisited, ” JCAP (2018), 005 [arXiv:1807.11495[astro-ph.CO]].[6] T. Harada, C. M. Yoo, K. Kohri, K. i. Nakao and S. Jhingan, “ Primordial black hole formation in the matter-dominated phase ofthe Universe, ” Astrophys. J. (2016) no.1, 61 [arXiv:1609.01588 [astro-ph.CO]].[7] A. Vilenkin, “Cosmological Density Fluctuations Produced by Vacuum Strings,” Phys. Rev. Lett. (1981), 1169-1172 [erratum:Phys. Rev. Lett. (1981), 1496] This limit corresponds to α → in appendix C. However, notice that, because of our definition of the realistic power spectrum in eq. (11), some(very mild) α -dependence is present (via the shape of the power spectrum) even if we take the gaussian limit in the computation of the abundance.The dot-dashed line in the left panel of fig. 6 corresponds to the case α = 0 . . [8] S. W. Hawking, “ Gravitational radiation from collapsing cosmic string loops, ” Phys. Lett. B (1990), 36-38[9] J. Fort and T. Vachaspati, “
Do global string loops collapse to form black holes?, ” Phys. Lett. B (1993), 41-46 [arXiv:hep-th/9305081 [hep-th]].[10] H. Deng, J. Garriga and A. Vilenkin, “
Primordial black hole and wormhole formation by domain walls, ” JCAP (2017), 050[arXiv:1612.03753 [gr-qc]].[11] A. A. Starobinsky, “ Spectrum of adiabatic perturbations in the universe when there are singularities in the inflation potential, ”JETP Lett. (1992) 489-494.[12] P. Ivanov, P. Naselsky and I. Novikov, “ Inflation and primordial black holes as dark matter, ” Phys. Rev. D (1994) 7173.[13] 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]].[14] S. M. Leach and A. R. Liddle, “ Inflationary perturbations near horizon crossing, ” Phys. Rev. D (2001), 043508 [arXiv:astro-ph/0010082].[15] 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].[16] Y. F. Cai, X. Tong, D. G. Wang and S. F. Yan, “ Primordial Black Holes from Sound Speed Resonance during Inflation, ” Phys. Rev.Lett. (2018) no.8, 081306 [arXiv:1805.03639 [astro-ph.CO]].[17] E. Cotner, A. Kusenko and V. Takhistov, “
Primordial Black Holes from Inflaton Fragmentation into Oscillons, ” Phys. Rev. D (2018) no.8, 083513 [arXiv:1801.03321 [astro-ph.CO]].[18] V. Atal and C. Germani, “ The role of non-gaussianities in Primordial Black Hole formation, ” Phys. Dark Univ. (2019), 100275[arXiv:1811.07857 [astro-ph.CO]].[19] M. Taoso, A. Urbano, “ Non-gaussianities for primordial black hole formation, ” in preparation.[20] 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]].[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. 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]].[23] 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]].[24] 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]].[25] 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]].[26] 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]].[27] M. Sasaki, T. Suyama, T. Tanaka and S. Yokoyama, “ Primordial black holes—perspectives in gravitational wave astronomy, ” Class.Quant. Grav. (2018) no.6, 063001 [arXiv:1801.05235 [astro-ph.CO]].[28] G. Ballesteros, J. Rey, M. Taoso and A. Urbano, “ Primordial black holes as dark matter and gravitational waves from single-fieldpolynomial inflation, ” JCAP (2020), 025 [arXiv:2001.08220 [astro-ph.CO]].[29] 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]].[30] I. Dalianis, A. Kehagias and G. Tringas, “ Primordial black holes from α -attractors, ” JCAP (2019), 037 [arXiv:1805.09483 [astro-ph.CO]].[31] C. T. Byrnes, P. S. Cole and S. P. Patil, “ Steepest growth of the power spectrum and primordial black holes, ” JCAP (2019), 028[arXiv:1811.11158 [astro-ph.CO]].[32] O. Özsoy and G. Tasinato, “ On the slope of the curvature power spectrum in non-attractor inflation, ” JCAP (2020), 048[arXiv:1912.01061 [astro-ph.CO]].[33] O. Özsoy, S. Parameswaran, G. Tasinato and I. Zavala, “ Mechanisms for Primordial Black Hole Production in String Theory, ”JCAP (2018), 005 [arXiv:1803.07626 [hep-th]].[34] M. Cicoli, V. A. Diaz and F. G. Pedro, “ Primordial Black Holes from String Inflation, ” JCAP (2018), 034 [arXiv:1803.02837 [hep-th]].[35] H. Cramer, “ Mathematical methods of statistics, ” Princeton University Press, 1999.[36] L. Comtet, “
Advanced Combinatorics, ” D. Reidel publishing company, 1974; pg. 134.[37] J. A. Shohat and J. D. Tamarkin, “
The problem of moments, ” Mathematical Surveys and Monographs, (1943).[38] J. R. Espinosa, D. Racco and A. Riotto, “ A Cosmological Signature of the SM Higgs Instability: Gravitational Waves, ” JCAP (2018), 012 [arXiv:1804.07732 [hep-ph]].[39] C. Yuan, Z. C. Chen and Q. G. Huang, “ Probing primordial–black-hole dark matter with scalar induced gravitational waves, ”Phys. Rev. D (2019) no.8, 081301 [arXiv:1906.11549 [astro-ph.CO]].[40] S. J. Kapadia, K. L. Pandey, T. Suyama and P. Ajith, “
Prospects for probing ultralight primordial black holes using the stochas-tic gravitational-wave background induced by primordial curvature perturbations, ” Phys. Rev. D (2020) no.12, 123535[arXiv:2005.05693 [astro-ph.CO]].[41] K. Inomata and T. Nakama, “
Gravitational waves induced by scalar perturbations as probes of the small-scale primordial spec-trum, ” Phys. Rev. D (2019) no.4, 043511 [arXiv:1812.00674 [astro-ph.CO]].[42] N. Bhaumik and R. K. Jain, “ Stochastic induced gravitational waves and lowest mass limit of primordial black holes with theeffects of reheating, ” [arXiv:2009.10424 [astro-ph.CO]].[43] 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]].[44] 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]].[45] 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]].[46] 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]].[47] 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]].[48] C. Caprini, M. Hindmarsh, S. Huber, T. Konstandin, J. Kozaczuk, G. Nardini, J. M. No, A. Petiteau, P. Schwaller and G. Servant, et al. “ Science with the space-based interferometer eLISA. II: Gravitational waves from cosmological phase transitions, ” JCAP (2016), 001 [arXiv:1512.06239 [astro-ph.CO]].[49] K. Yagi and N. Seto, “ Detector configuration of DECIGO/BBO and identification of cosmological neutron-star binaries, ” Phys. Rev.D (2011), 044011 [erratum: Phys. Rev. D (2017) no.10, 109901] [arXiv:1101.3940 [astro-ph.CO]].[50] J. Coleman [MAGIS-100], “ Matter-wave Atomic Gradiometer InterferometricSensor (MAGIS-100) at Fermilab, ” PoS
ICHEP2018 (2019), 021 [arXiv:1812.00482 [physics.ins-det]].[51] M. Maggiore, C. Van Den Broeck, N. Bartolo, E. Belgacem, D. Bertacca, M. A. Bizouard, M. Branchesi, S. Clesse, S. Foffa andJ. García-Bellido, et al. “ Science Case for the Einstein Telescope, ” JCAP (2020), 050 [arXiv:1912.02622 [astro-ph.CO]].[52] B. P. Abbott et al. [LIGO Scientific and Virgo], “ Search for the isotropic stochastic background using data from Advanced LIGO’ssecond observing run, ” Phys. Rev. D (2019) no.6, 061101 [arXiv:1903.02886 [gr-qc]].[53] L. Barsotti, P. Fritschel, M. Evans and S. Gras, Updated Advanced LIGO sensitivity design curve.[54] Z. C. Chen, F. Huang and Q. G. Huang, “
Stochastic Gravitational-wave Background from Binary Black Holes and Binary NeutronStars and Implications for LISA, ” Astrophys. J. (2019) no.1, 97 [arXiv:1809.10360 [gr-qc]].[55] S. Young, C. T. Byrnes and M. Sasaki, “
Calculating the mass fraction of primordial black holes, ” JCAP (2014), 045[arXiv:1405.7023 [gr-qc]].[56] A. Riotto, “ Inflation and the theory of cosmological perturbations, ” ICTP Lect. Notes Ser. (2003), 317-413 [arXiv:hep-ph/0210162 [hep-ph]].[57] C. Kiefer and D. Polarski, Annalen Phys. (1998), 137-158 [arXiv:gr-qc/9805014 [gr-qc]].[58] J. M. Maldacena, “ Non-Gaussian features of primordial fluctuations in single field inflationary models, ” JHEP (2003), 013[arXiv:astro-ph/0210603 [astro-ph]].[59] G. Ballesteros, J. Rey, M. Taoso and A. Urbano, “ Stochastic inflationary dynamics beyond slow-roll and consequences for primor-dial black hole formation, ” [arXiv:2006.14597 [astro-ph.CO]].[60] N. Kogo and E. Komatsu, “
Angular trispectrum of cmb temperature anisotropy from primordial non-gaussianity with the fullradiation transfer function, ” Phys. Rev. D (2006), 083007 [arXiv:astro-ph/0602099 [astro-ph]].[61] S. Young and C. T. Byrnes, “ Primordial black holes in non-Gaussian regimes, ” JCAP (2013), 052 [arXiv:1307.4995 [astro-ph.CO]].[62] 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]].[63] I. Musco, V. De Luca, G. Franciolini and A. Riotto, “ The Threshold for Primordial Black Hole Formation: a Simple Analytic Pre-scription, ” [arXiv:2011.03014 [astro-ph.CO]].[64] M. Shibata and M. Sasaki, “
Black hole formation in the Friedmann universe: Formulation and computation in numerical rela-tivity, ” Phys. Rev. D (1999), 084002 [arXiv:gr-qc/9905064 [gr-qc]].[65] A. Helou, I. Musco and J. C. Miller, “ Causal Nature and Dynamics of Trapping Horizons in Black Hole Collapse, ” Class. Quant.Grav. (2017) no.13, 135012 [arXiv:1601.05109 [gr-qc]].[66] T. Harada, C. M. Yoo and K. Kohri, “ Threshold of primordial black hole formation, ” Phys. Rev. D (2013) no.8, 084051[arXiv:1309.4201 [astro-ph.CO]].[67] A. Escrivà, C. Germani and R. K. Sheth, “ Universal threshold for primordial black hole formation, ” Phys. Rev. D (2020) no.4,044022 [arXiv:1907.13311 [gr-qc]].[68] D. W. Neilsen and M. W. Choptuik, “
Critical phenomena in perfect fluids, ” Class. Quant. Grav.17