An Improved Calculation of the Non-Gaussian Halo Mass Function
Guido D'Amico, Marcello Musso, Jorge Noreña, Aseem Paranjape
aa r X i v : . [ a s t r o - ph . C O ] F e b An Improved Calculation of the Non-Gaussian Halo Mass Function
Guido D’Amico a , b , Marcello Musso c , Jorge Nore˜na a , b , Aseem Paranjape c a SISSA, via Bonomea 265, 34136 Trieste, Italy b INFN - Sezione di Trieste, via Bonomea 265, 34136 Trieste, Italy c Abdus Salam International Centre for Theoretical PhysicsStrada Costiera 11, 34151, Trieste, Italy
Abstract
The abundance of collapsed objects in the universe, or halo mass function, is an important theoretical toolin studying the effects of primordially generated non-Gaussianities on the large scale structure. The non-Gaussian mass function has been calculated by several authors in different ways, typically by exploiting thesmallness of certain parameters which naturally appear in the calculation, to set up a perturbative expansion.We improve upon the existing results for the mass function by combining path integral methods and saddlepoint techniques (which have been separately applied in previous approaches). Additionally, we carefullyaccount for the various scale dependent combinations of small parameters which appear. Some of these com-binations in fact become of order unity for large mass scales and at high redshifts, and must therefore betreated non-perturbatively. Our approach allows us to do this, and also to account for multi-scale densitycorrelations which appear in the calculation. We thus derive an accurate expression for the mass functionwhich is based on approximations that are valid over a larger range of mass scales and redshifts than thoseof other authors. By tracking the terms ignored in the analysis, we estimate theoretical errors for our resultand also for the results of others. We also discuss the complications introduced by the choice of smoothingfilter function, which we take to be a top-hat in real space, and which leads to the dominant errors in ourexpression. Finally, we present a detailed comparison between the various expressions for the mass functions,exploring the accuracy and range of validity of each.
The primordial curvature inhomogeneities, generated by the inflationary mechanism, obey nearly Gaussianstatistics. The deviations from Gaussianity, while expected to be small, provide a unique window into thephysics of inflation. For example, single-field slow-roll models of inflation lead to a small level of non-Gaussianity (NG), so that an observation of a large NG would indicate a deviation from this paradigm.Until a few years ago, the main tool to constrain NG was considered to be the statistics of the cosmicmicrowave background (CMB) temperature field, since inhomogeneities at the CMB epoch are small and theirphysics can be described by a perturbative treatment. In recent years, however, thanks to observations anddevelopments in the theory, the large-scale structure (LSS) of the universe has emerged as a complementaryprobe to constrain primordial NG. While it is true that the n -point functions of the density field on small scalesare dominated by the recent gravitational evolution, and do not reflect anymore the statistics of primordialperturbations, it turns out that the abundance of very massive objects, which form out of high peaks of thedensity perturbations, is a powerful probe of primordial NG. In this context, much attention has been givenrecently to three possible methods of constraining the magnitude and shape of the primordial NG with the LSS:the galaxy power spectrum, the galaxy bispectrum and the mass function. It was pointed out in Refs. [1, 2]that a NG of a local type induces a scale dependence on the galaxy power spectrum, thus making it a sensitiveprobe of the magnitude of local NG f locNL . From Ref. [3] one finds the following constraints: − < f locNL < +69, E-mail: [email protected], [email protected], [email protected], [email protected] lready comparable with those obtained from CMB measurements in Ref. [4]: − < f locNL < +74. The futureis even more promising, with precisions of ∆ f locNL ∼
10 [5] and ∆ f locNL ∼ M and M + dM , dndM = ¯ ρM f ( σ ) (cid:12)(cid:12)(cid:12)(cid:12) d ln σd ln M (cid:12)(cid:12)(cid:12)(cid:12) , (1)where ¯ ρ is the average density of the universe, σ ( M ) is the variance of the density contrast δ R filtered onsome comoving scale R corresponding to the mass M , and the function f ( σ ) is to be computed. Throughoutthis work, we will refer to f ( σ ) itself as the mass function. A very useful tool in the analysis is the sphericalcollapse model [13], which predicts that the value of the linearly extrapolated density contrast of a sphericalhalo, at the time when the halo collapses, is δ c ≃ . f ( σ ) in the case of Gaussian initialconditions. Their calculation however suffered from a problem of undercounting which affects the overallnormalization – their approach does not count underdense regions embedded in larger overdense regionsas eventually collapsed objects. To account for this discrepancy, PS introduced an ad-hoc factor of 2 bydemanding that the mass function be correctly normalized, such that all the mass in the universe must becontained in collapsed objects. In the excursion set approach, Bond et al . [15] resolved this issue and deriveda correctly normalized mass function, for Gaussian initial conditions. They argued that the filtered densitycontrast δ R follows a random walk as a function of the filtering scale, and the problem of computing f ( σ ) istranslated into the problem of finding the rate of “first crossing” of the barrier δ c , whose solution is well-known.We will study this formalism in detail in section 3 for the more general non-Gaussian case.Turning to non-Gaussianities, the most popular non-Gaussian mass functions are those due to Matarrese,Verde and Jimenez [16] (MVJ) and LoVerde et al . [17] (LMSV). Both groups used the PS approach, bymodifying the probability density function for the (linearized) density contrast to describe non-Gaussianinitial conditions. In their prescription, the relevant object is the ratio R ng of non-Gaussian to Gaussian massfunctions. The full mass function is usually taken as the product of R ng and an appropriate Gaussian massfunction as given by N -body simulations, e.g. the Sheth & Tormen mass function [18]. It is not clear howeverthat this is the correct way to proceed. Indeed, in a series of papers [19, 20, 21], Maggiore & Riotto (MR)presented a rigorous approach to the first-passage problem in terms of path integrals, and in Ref. [21] theypointed out that a PS-like prescription in fact misses some important non-Gaussian effects stemming from3-point correlations between different scales (so-called “unequal time” correlators).On the other hand, MR treated non-Gaussian contributions to f ( σ ) by simply linearizing in the 3-pointfunction of δ R , i.e. by linearizing in the non-Gaussian parameter f NL . Since the NG are assumed to besmall, in the sense that the parameter ǫ = h δ i /σ satisfies ǫ ≪
1, one might expect that such a perturbativetreatment is valid. However, another crucial ingredient in the problem is that the length scales of interest arelarge, which leads to a second small parameter ν − where ν = δ c /σ . This is evident in the calculations of MR,who crucially use ν − ∝ σ as a small parameter. Any perturbative treatment now depends not only on thesmallness of ǫ and ν − individually, but also on the specific combinations of these parameters which appearin the calculation. It is known (and we will explicitly see below) that a natural combination that appearsis ǫν , which can become of order unity on scales of interest. The mass functions given by LMSV and MRtherefore break down as valid series expansions when this occurs. Interestingly, MVJ’s PS-like treatment onthe other hand involved a saddle point approximation, allowing them to non-perturbatively account for the ν term (which appears in an exponential in their approach). For a discussion, see Ref. [22].It appears to us therefore, that there is considerable room for improvement in the theoretical calculationof the mass function. The goal of our paper is twofold. Firstly, we present a rigorous calculation of the massfunction in the following way : (a) we use the techniques developed by MR in Refs. [19, 20, 21], which allow usto track the complex multi-scale correlations involved in the calculation, and (b) we demonstrate that MR’sapproach can be combined with saddle point techniques (used by MVJ), to non-perturbatively handle termswhich can become of order unity. This leads to an expression for the mass function which is valid on muchlarger scales than those presented by MR and LMSV. Secondly, by keeping track of the terms ignored, wecalculate theoretical error bars on the expressions for f ( σ ) resulting not only from our own calculations, butalso for those of the other authors [16, 17, 21]. Since the terms ignored depend on ν in general, these errorbars are clearly scale dependent. This allows us to estimate the validity of each of the expressions for themass function at different scales, but importantly it also allows us to analytically compare between differentexpressions. In this paper we will not explicitly account for effects of the ellipsoidal collapse model [23, 24],since these are expected to be negligible on the very large scales which are of interest to us. For a recenttreatment of ellipsoidal collapse effects in the presence of non-Gaussianities on scales where ǫν ≪
1, seeLam & Sheth [25]. For a different approach to computing the non-Gaussian halo mass function, see Ref. [1],where the authors proposed that this mass function can be approximated as a convolution of the
Gaussian mass function with a probability distribution function that maps between halos identified in Gaussian andnon-Gaussian N -body simulations. This probability distribution itself was approximated as a Gaussian, withmean and variance fit from simulations. While in practice this approach is easy to implement, it relies heavilyon the output of N -body simulations. Our approach, on the other hand, allows us to compute a mass functionalmost entirely from first principles.This paper is organized as follows. In section 2 we fix some notation and briefly introduce the two mostpopular shapes of primordial NG, i.e. the local and equilateral ones. In section 3 we present our calculationof the mass function. In section 4 we discuss certain subtleties regarding the truncation of the perturbativeseries, and also compare with the other expressions for f ( σ ) mentioned above. In section 5 we discuss theeffects induced by some additional complications introduced in the problem due to the specific choice of thefilter function [19], which we take to be a top-hat in real space, and due to the inclusion of stochasticity inthe value of the collapse threshold δ c (which is also expected to partially account for effects of ellipsoidalcollapse) [20] . In section 6 we compare our final result Eqn. (43) with those of other authors, includingtheoretical errors for each, and conclude with a brief discussion of the results and directions for future work.Some technical asides have been relegated to the Appendices. We need to relate the linearly evolved density field to the primordial curvature perturbation, which carries theinformation of the non-linearities produced during and after inflation. We start from the Bardeen potentialΦ on subhorizon scales, given by Φ( k , z ) = − T ( k ) D ( z ) a R ( k ) , (2)where R ( k ) is the (comoving) curvature perturbation, which stays constant on superhorizon scales; T ( k ) isthe transfer function of perturbations, normalized to unity as k →
0, which describes the suppression of powerfor modes that entered the horizon before the matter-radiation equality; and D ( z ) is the linear growth factorof density fluctuations, normalized such that D ( z ) = (1+ z ) − in the matter dominated era. Then, the densitycontrast field is related to the potential by the Poisson equation, which in Fourier space reads δ ( k , z ) = − ak m H Φ( k , z ) = 2 k m H T ( k ) D ( z ) R ( k ) ≡ M ( k, z ) R ( k ) , (3)where we substituted Eqn. (2). Here, Ω m is the present time fractional density of matter (cold dark matterand baryons), and H = 100 h km s − Mpc − is the present time Hubble constant. The redshift dependence s trivially accounted for by the linear growth factor D ( z ) and in the following, for notational simplicity, wewill often suppress it. All our calculations will use a reference ΛCDM cosmology compatible with WMAP7data [4], using parameters h = 0 . m = 0 . b = 0 . n s = 0 .
961 and σ = 0 . σ is the variance of the density field smoothed on a length scaleof 8 h − Mpc. For simplicity, for the transfer function T ( k ) we use the BBKS form, proposed in Bardeen etal . [23]: T BBKS ( x ) ≡ . x ln (1 + 2 . x ) (cid:0) . x + (16 . x ) + (5 . x ) + (6 . x ) (cid:1) − / , (4)where x ≡ k ( h Mpc − ) / Γ with a shape parameter Γ = Ω m h exp h − Ω b (1 + √ h/ Ω m ) i that accounts for bary-onic effects as described in Ref. [26]. For more accurate results, one could use a numerical transfer function,as obtained by codes like CMBFAST [27] or CAMB [28]; the results are not expected to be qualitativelydifferent.In order to study halos, which form where an extended region of space has an average overdensity whichis above threshold, it is useful to introduce a filter function W R ( | x | ), and consider the smoothed density field(around one point, which we take as the origin), δ R = Z d k (2 π ) f W ( kR ) δ ( k ) , (5)where f W ( kR ) is the Fourier transform of the filter function. For all numerical calculations we will use thespherical top-hat filter in real space, whose Fourier transform f W ( kR ) is given by f W ( y ) = 3 y (sin y − y cos y ) . (6)This choice allows us to have a well-defined relation between length scales and masses, namely M = (4 π/ m ρ c R with ρ c = 3 H / (8 πG ) = 2 . · h − M sol ( h − Mpc) − . However it introduces some complexities in the anal-ysis, which we will comment on later. By using Eqns. (5) and (3) we have, for the 3-point function, h δ R δ R δ R i c = Z d k (2 π ) d k (2 π ) d k (2 π ) f W ( k R ) f W ( k R ) f W ( k R ) M ( k ) M ( k ) M ( k ) h R ( k ) R ( k ) R ( k ) i c , (7)where the subscript c denotes the connected part, and analogous formulae are valid for the higher ordercorrelations. The function h R ( k ) R ( k ) R ( k ) i c encodes information about the physics of the inflationary epoch. Bytranslational invariance, it is proportional to a momentum-conserving delta function: h R ( k ) R ( k ) R ( k ) i c = (2 π ) δ D ( k + k + k ) B R ( k , k , k ) , (8)where the (reduced) bispectrum B R ( k , k , k ) depends only on the magnitude of the k ’s by rotational invari-ance. According to the particular model of inflation, the bispectrum will be peaked about a particular shapeof the triangle. The two most common cases are the squeezed (or local) NG, peaked on squeezed triangles k ≪ k ≃ k , and the equilateral NG, peaked on equilateral triangles k ≃ k ≃ k . Indeed, one can definea scalar product of bispectra, which describes how sensitive one is to a NG of a given type if the analysis isperformed using some template form for the bispectrum. As expected, the local and equilateral shapes areapproximately orthogonal with respect to this scalar product [29]. We will now describe these two models inmore detail. The local model: he local bispectrum is produced when the NG is generated outside the horizon, for instance in the curvatonmodel [30, 31] or in the inhomogeneous reheating scenario [32]. In these models, the curvature perturbationcan be written in the following form, R ( x ) = R g ( x ) + 35 f locNL (cid:0) R g ( x ) − h R g i (cid:1) + 925 g NL R g ( x ) , (9)where R g is the linear, Gaussian field. We have included also a cubic term, which will generate the trispectrumat leading order. The bispectrum is given by B R ( k , k , k ) = 65 f locNL [ P R ( k ) P R ( k ) + cycl . ] , (10)where “cycl.” denotes the 2 cyclic permutations of the wavenumbers, and P R ( k ) is the power spectrum givenby P R ( k ) = Ak n s − . The trispectrum is given by h R ( k ) R ( k ) R ( k ) R ( k ) i c = (2 π ) δ D ( k + k + k + k ) × f X b Models with derivative interactions of the inflaton field [33, 34, 35] give a bispectrum which is peaked aroundequilateral configurations, whose specific functional form is model dependent. Moreover, the form of thebispectrum is usually not convenient to use in numerical analyses. This is why, when dealing with equilateralNG, it is convenient to use the following parametrization, given in Ref. [36], B R ( k , k , k ) = 185 f equilNL A h k − n s k − n s + 13( k k k ) − n s ) / − k k k ) (4 − n s ) / + 5 perms. i . (12)This is peaked on equilateral configurations, and its scalar product with the bispectra produced by the realisticmodels cited above is very close to one. Therefore, being a sum of factorizable functions, it is the standardtemplate used in data analyses. We now turn to the main calculation of the paper. The non-Gaussian halo mass function can be obtained bycalculating the barrier first crossing rate F of a random walk with non-Gaussian noise, in the presence of anabsorbing barrier. This can be done perturbatively, starting from a path integral approach as prescribed byMR [19, 21] and the mass function can be shown to be f ( σ ) = 2 σ F ( σ ). As discussed by MR, the calculationof f involves certain assumptions regarding the type of filter used and also the location of the barrier. Inparticular, the formalism is simplest for a sharp filter in k -space, and using the spherical top-hat of Eqn.(6) introduces complications in the form of non-Markovian effects. Further, in order to make the sphericalcollapse ansatz more realistic and obtain better agreement with N -body simulations, MR show that it isuseful to treat the location of the barrier δ c as a stochastic variable itself, and allow it to diffuse. For the timebeing, we will ignore these complications, and will return to their effects in section 5.To make the paper self-contained, we begin with a brief review of the path integral approach to thecalculation of the mass function. The reader is referred to Ref. [19] for a more pedagogical introduction. Inthe path integral approach, one treats the variance σ R ≡ h ˆ δ R i as a “time” parameter, t ≡ σ R , and considersthe random walk followed by the smoothed density field ˆ δ R as this “time” is increased in discrete steps startingfrom small values (equivalently, as R is decreased from very large values). Here ˆ δ R ( ~x ) is a stochastic quantityin real space due to the stochasticity inherent in the initial conditions. We use the notation ˆ δ R to distinguish he stochastic variable from the values it takes, which will be noted by δ i below. We probe this stochasticityby changing the smoothing scale at a fixed location ~x = 0, thus making the variable perform a random walk,which obeys a Langevin equation ∂ ˆ δ∂t = ˆ η , (13)with a stochastic noise ˆ η whose statistical properties depend on the choice of filter used. In particular, for atop hat filter in k -space, the noise is white, i.e. its 2-point function is a Dirac delta [15], h ˆ η ( t )ˆ η ( t ) i = δ D ( t − t ) (14)The random walk can be described as a trajectory { δ , δ , . . . , δ n } which starts with ˆ δ taking the value δ = 0at t = 0 (or R → ∞ which is the homogeneous limit), then taking values δ i at times t i , finally arriving at δ n at time t n , with a discrete timestep ∆ t = t k +1 − t k = t n /n . The probability P ( t ) that the trajectory crossesthe barrier at δ c at a time larger than some t (i.e. at scales smaller than the corresponding R or M ), is thesame as the probability that the trajectory did not cross the barrier at any time smaller than t , so that P ( t ) = Z δ c −∞ d δ . . . d δ n W ( { δ j } ; t ) , (15)where the probability density over the space of trajectories, W ( { δ j } ; t ) is defined as W ( { δ j } ; t ) ≡ h δ D (ˆ δ ( t ) − δ ) . . . δ D (ˆ δ ( t n ) − δ n ) i , (16)where δ D is the Dirac delta distribution. The first crossing rate is given by the negative time derivative of P , F = − ∂ t P , and the mass function is then f = 2 t F ( t ). In Eqn. (16) one can write the Dirac deltas using theintegral representation δ D ( x ) = R ∞−∞ dλe − iλx / π , to obtain W ( { δ j } ; t ) = Z ∞−∞ dλ π . . . dλ n π h e − i P j λ j ˆ δ ( t j ) i e i P j λ j δ j . (17)The object h e − i P j λ j ˆ δ j i is the exponential of the generating functional of the connected Green’s functions,and can be shown to reduce to [37] h e − i P j λ j ˆ δ j i = exp ∞ X p =2 ( − i ) p p ! n X j ,..,j p =1 λ j . . . λ j p h ˆ δ j . . . ˆ δ j p i c , (18)where h ˆ δ j . . . ˆ δ j p i c is the connected p -point function of ˆ δ , with the short-hand notation ˆ δ j = ˆ δ ( t j ). k filter In the Gaussian case, all connected n -point correlators vanish except for n = 2, and in the Markovian (sharp- k filter) case which we are considering, the 2-point function becomes h ˆ δ j ˆ δ k i = min( t j , t k ), where min( t j , t k ) isthe minimum of t j and t k . The resulting n -dimensional Gaussian integral can be handled in a straightforwardway to obtain W gm = n − Y k =0 Ψ ∆ t ( δ k +1 − δ k ) ; Ψ ∆ t ( x ) = (2 π ∆ t ) − / e − x / (2∆ t ) , (19)where we follow MR’s notation and use the superscript “gm” to denote “Gaussian Markovian”. As MR haveshown [19], the resulting expression for P gauss ( t ) in the continuum limit ∆ t → P gauss = Z δ c −∞ dδ . . . dδ n W gm = erf (cid:18) ν √ (cid:19) , (20)where we use the notation ν ≡ δ c /σ . (This in principle also includes the redshift dependence of the collapsethreshold δ c , see below.) This expression for the continuum limit probability P gauss is of course a well-known esult going back to Chandrasekhar [38]. This leads to the standard excursion set result for the Gaussianmass function f PS = − t∂ t | δ c P gauss , f PS ( ν ) = r π ν e − ν / , (21)where we use the subscript PS (for Press-Schechter), to conform with the conventional notation for this object. k filter In the non-Gaussian case (but still retaining the sharp- k filter), the probability density W ( { δ j } ; t ) also getscontributions from connected n -point correlators with n ≥ 3, since these in general do not vanish. These canbe handled by using the relation λ k e i P j λ j δ j = − i∂ k e i P j λ j δ j , with ∂ j ≡ ∂/∂δ j . A straightforward calculationthen shows the mass function to be f = − t ∂∂t (cid:12)(cid:12)(cid:12)(cid:12) δ c Z δ c −∞ d δ . . . d δ n exp (cid:20) − n X j,k,l =1 h ˆ δ j ˆ δ k ˆ δ l i c ∂ j ∂ k ∂ l + 14! n X j,k,l,m =1 h ˆ δ j ˆ δ k ˆ δ l ˆ δ m i c ∂ j ∂ k ∂ l ∂ m + . . . (cid:21) W gm , (22)where it is understood that one takes the continuum limit ∆ t → t . We will find it useful to change variables from ( δ c , t ) to ( ν , t ), in which case the partialderivative becomes − t ( ∂/∂t ) | δ c = ν ( ∂/∂ν ) | t − t ( ∂/∂t ) | ν ≡ ν∂ ν − t∂ t . (23)It is also useful at this stage to take a small detour and introduce some notation which we will use throughoutthe rest of the paper. We define the scale dependent “equal time” functions ε n − ≡ h ˆ δ nR i c σ nR ; n ≥ , (24)which as we will see, remain approximately constant over the scales of interest. We assume the ordering ε n − ∼ O ( ǫ n − ) with ǫ ≪ 1, which can be motivated from their origin in inflationary physics, where one finds ε ∼ f NL A / , ε ∼ g NL A , etc . Typically we expect ǫ . − for f NL . ε and ε in the local and equilateral models, as a function of t = σ R . We see e.g. that ε in the local modelis comparable to ε . In the literature one usually encounters the reduced cumulants S n , which are related tothe ε n − by ε = σ S , ε = σ S and so on. The motivation for using the S n comes from the study of NGinduced by nonlinear gravitational effects. However, as we see from Fig. 1, when studying primordial NG itis more meaningful to consider the ε n which are approximately scale-independent and perturbatively ordered.We will soon see that a natural expansion parameter that arises in the calculation has the form ∼ ǫν , andwe therefore require that the mass scales under scrutiny are not large enough to spoil the relation ǫν ≪ 1. Itturns out that observationally interesting mass scales can nevertheless be large enough to satisfy ǫν ∼ O (1).Fig. 2 shows the behaviour of ε ν and ε ν at different redshifts, as a function of mass, in our reference ΛCDMmodel for local type NG, with f locNL = 100. The behaviour for the equilateral NG is qualitatively similar. Theredshift dependence of these quantities comes from the definition of ν , ν ( M, z ) ≡ √ a δ c σ ( M ) D (0) D ( z ) ≡ δ c ( z ) σ ( M ) , (25)where we denote the usual spherical collapse threshold as δ c = (3 / π/ / ≃ . δ c for thefull, redshift dependent quantity, and a is a parameter accounting for deviations from the simplest collapsemodel. In the standard spherical collapse picture we have a = 1. A value of a different from unity (specifically Notationally we distinguish the order parameter ǫ from the specific NG functions ε and ε . .00010.0010.010.10.1 1 ε σ LocalEquilateral (a) σ ε (f NL2 =10 ) ε (g NL =10 ) ε (b) Figure 1: Scale dependence of the ε n . Panel (a) : Behaviour of ε vs. σ in the local and equilateral models, for f NL = 100 in each case. Panel (b) : Behaviour of 2nd order ( ∼ ǫ ) terms. We show ε for the local model with f NL = 100 and g NL = 10 . The terms proportional to f and g NL are shown separately. Also shown is ε for thesame model. √ a ≃ . 89) can be motivated by allowing the collapse threshold to vary stochastically [20], as we will discussin section 5. We will soon see that the object ǫν appears naturally in the calculation, and to be definite wewill assume ǫν ∼ O (1) for now. In section 4 we will discuss the effects of relaxing this condition and probingsmaller length scales.We now turn to the “unequal time” correlators appearing in Eqn. (22). Since we are concerned withlarge scales, we are in the small t limit, and following MR we expand the n -point correlators around the “finaltime” t . We can define the Taylor coefficients G ( p,q,r )3 ( t ) ≡ " d p dt pj d q dt qk d r dt rl h ˆ δ ( t j )ˆ δ ( t k )ˆ δ ( t l ) i c t j = t k = t l = t , (26)and then expand h ˆ δ j ˆ δ k ˆ δ l i c = ∞ X p,q,r =0 ( − p + q + r p ! q ! r ! G ( p,q,r )3 ( t )( t − t j ) p ( t − t k ) q ( t − t l ) r . (27)For the 4-point function we will have an analogous expression involving coefficients G ( p,q,r,s )4 .Since calculations involving a general set of coefficients G , G , etc. are algebraically rather involved, wefind it useful to first consider an example in which these coefficients take simple forms. In this toy model weassume that the ε n are exactly constant, and moreover that the n -point correlators take the form h ˆ δ j ˆ δ k ˆ δ l i c = ε ( t j t k t l ) / ; h ˆ δ j ˆ δ k ˆ δ l ˆ δ m i c = ε ( t j t k t l t m ) / . (28)For clarity, we will display details of the calculation only for this model. We have relegated most of thetechnical details of our calculation to Appendix A. In Appendix A.1 we show that the mass function for thismodel can be brought to the form f = (cid:18) π (cid:19) / ν e − ( ε / ∂ ν +( ε / ∂ ν + ... (cid:20) e − ν / − ε ν e − ν / + 516 ε (cid:16) π (cid:17) / erfc (cid:18) ν √ (cid:19) + 18 (cid:18) ε − ε (cid:19) e − ν / (cid:0) ν − (cid:1) + O ( ǫ ν ) (cid:21) , (29) Throughout the paper we will consider at most 4-point correlators. This truncation is justified given our assump-tions, as we will see later. .010.11 10 ε ν / M (h -1 M sol ) z = 1z = 0.5z = 0 (a) ε ν M (h -1 M sol )z = 1z = 0.5z = 0 (b) Figure 2: Panel (a) : Behaviour of ε ν / f NL = 100. The threecurves correspond to different redshifts. The horizontal line corresponds to ε ν / ε ν with the same setup as in panel (a). where we ignore terms like ε ν O ( ν − ). The remaining exponentiated derivative can be handled using a saddlepoint approximation. We show how to do this in Appendix A.2. The expression for the mass function f ( ν )works out to f ( ν ) = (cid:18) π (cid:19) / ν exp (cid:20) − ν (cid:18) − ε ν + 14 (cid:16) ε − ε (cid:17) ν + O ( ǫ ν ) (cid:19)(cid:21) × (cid:18) − ε ν (cid:18) − ν (cid:19) + (cid:16) ε − ε (cid:17) ν + O ( ǫ ν ) (cid:19) , (30)which superficially at least, is comprised of two series expansions, one in the exponential and one as apolynomial, both based on the small parameter ǫν (see however the next section).This derivation assumed that ε and ε are constant, and that the unequal time correlations are givenby Eqn. (28). It is straightforward to relax these assumptions and perform the calculation with the exactstructure of the correlations. Appendix A.3 shows how to do this, and the result is f ( ν, t ) = (cid:18) π (cid:19) / ν exp (cid:20) − ν (cid:18) − ε ν + 14 (cid:16) ε − ε (cid:17) ν + O ( ǫ ν ) (cid:19)(cid:21) × (cid:26) − ε ν (cid:18) (4 − c ) + 1 ν (cid:18) c − c − (cid:19)(cid:19) + 18 ν (cid:18) ε (11 − c ) − ε (cid:18) c − d ln ε d ln t (cid:19)(cid:19) + O ( ǫ ν ) (cid:27) , (31)where the functions c n ( t ) are defined in Eqn. (A.17) and characterize the behaviour of the unequal timecorrelations. In our toy model above, the c n reduce to unity. Indeed, as a check we see that the expression in(31) reduces to Eqn. (30) if we take ε , ε to be constant and set the c n to unity.One issue which we have ignored so far, is that the definition of ν involves the variance t = σ of the non-Gaussian field. Computationally it is more convenient to work with the variance σ of the Gaussian field in terms of which cosmological NG are typically defined. We should then ask whether this differencewill require changes in our expressions for f . We start by noting that this difference in variances is of order ∼ ǫ . For example, in the local model one has σ ( R ) = Ad ( R ) + A ( Af ) d ( R ) where A ∼ − is an overallnormalization constant, d and d are scale dependent functions of comparable magnitude on all relevantscales, and ǫ is estimated as ǫ ∼ f NL A / . We therefore have ν = δ c /σ = ( δ c /σ g )(1 + O ( ǫ )). However, withour assumption that ǫν ∼ O (1), we see that this correction is actually of order ∼ ( ǫ ν ) ν − ∼ ǫ ν , whichwe have been consistently ignoring. We will see that even when we relax the assumption ǫν ∼ O (1) and robe smaller scales where ǫν ≪ 1, this correction can still be consistently ignored. Hence we can safely set ν = δ c /σ g in all of our expressions. Now that all the derivative operators which we consider important have been accounted for, we can checkwhether our final result is consistently truncated, i.e. whether we have retained all terms at any given orderin the expansion. Symbolically, our current result for the mass function can be written as f ∼ e − ν ( ǫν + ǫ ν + O ( ǫ ,ǫ ν ) ) h ǫν + ǫν + ǫ ν + O ( ǫν − , ǫ , ǫ ν ) i , (32)with the understanding that coefficients are computed (but not displayed) for all terms except those indicatedby the O () symbols. Also, ǫ refers to both ε and ε .Since the expansions involve two parameters, ǫν and ν − , they make sense only if we additionally prescribea relation between these parameters. So far we assumed that ǫ is fixed and ν is such that ǫν ≃ 1, which wasbased on the observation that the term ǫν naturally appears in the exponent and is not restricted in principleto small values. In Appendix B.1 we discuss this condition in more detail, and also analyse the consequencesof relaxing this condition and probing smaller mass scales. We find that for observationally accessible massscales larger than the scale where ǫν ≃ ν − , the single expression f ∼ e − ν ( ǫν + ǫ ν ) h ǫν + ǫν + O ( ǫ ν , ǫ ν , ǫν − ) i , (33)is parametrically consistent as it stands – the terms ignored are smaller than the smallest terms retained – andin fact it remains a very good approximation even when ǫν ≃ 1, since the only “inconsistent” term then is ǫν − , whose effect reduces as ν increases. On scales where ǫν ≃ ν − and lower, the theoretical error becomescomparable to or larger than the quadratic term in the exponential. Plugging back all the coefficients, wehave the following result for the mass function (excluding filter effects, see section 5), f ( ν, t ) = f PS ( ν ) exp (cid:18) ε ν − (cid:16) ε − ε (cid:17) ν (cid:19) × (cid:26) − ε ν (cid:18) (4 − c ) + 1 ν (cid:18) c − c − (cid:19)(cid:19) + O ( ǫ ν , ǫ ν , ǫν − ) (cid:27) . (34) In this subsection we compare our results with previous work on the non-Gaussian mass function. As men-tioned in the introduction, this quantity has been computed by several authors in different ways [16, 17, 21].If one considers the range of validity of the perturbative expansion, the strongest result so far has been dueto MVJ [16], who explicitly retain the exponential dependence on ε . Their expression for f can be writtenas f MVJ = f PS ( ν ) e ε ν / (1 − ε ν/ / (cid:18) − ε ν (cid:18) − d ln ε d ln t (cid:19)(cid:19) . (35)The major shortcoming of their result is that it is based on a Press-Schechter like prescription, and musttherefore be normalized by an appropriate Gaussian mass function, typically taken to be the one due to The analysis presented by MVJ in fact allows one to retain terms like ∼ ǫ ν in the exponential as well, and wehave seen that when ǫν ≃ 1, these terms are as important as the polynomial ǫν term retained by MVJ. However,since the MVJ expression misses unequal time effects of order ∼ ǫν anyway, it is reasonable to compare our resultswith the expression (35), which is also the one used by most other authors (see e.g. Refs. [39, 40]). heth & Tormen [18]. Additionally, it always misses the contributions due to the unequal time correlators,which contribute to the terms ∼ ǫν , ǫν − , etc. in Eqn. (34). When one considers formal correctness onthe other hand, MR have presented a result based on explicit path integrals, which accounts for the unequaltime contributions, and which also does not need any ad hoc normalizations (in this context see also Lam &Sheth [25]). Indeed, our calculations in section 3 were based on techniques discussed by MR in Refs. [19, 21].As we discuss below however, the fact that MR do not explicitly retain the exponential dependence of ε ν ,means that their result is subject to significant constraints on the range of its validity. Their expression for f , ignoring filter effects, is f MR = f PS ( ν ) (cid:18) ε ν (cid:26) − ν (4 − c ) − ν (cid:18) c − c − (cid:19)(cid:27)(cid:19) . (36)This expression is precisely what one obtains by linearizing our expression (34) in ε , which serves as a checkon our calculations. LMSV [17] present a result based on an Edgeworth expansion of the type encounteredwhen studying NG generated by nonlinear gravitational effects [41]. The result most often quoted in theliterature is their expression linear in ε (and hence in ε ν ), which is f LMSV , lin = f PS ( ν ) (cid:18) ε ν (cid:26) − ν (cid:18) − d ln ε d ln t (cid:19) − ν d ln ε d ln t (cid:27)(cid:19) . (37)In Appendix B.3 of Ref. [17], LMSV also give an expression involving ε and ε , which can be written as f LMSV , quad = f PS ( ν ) (cid:20) ε (cid:18) H ( ν ) + 2 ν d ln ε d ln t H ( ν ) (cid:19) + 172 ε (cid:18) H ( ν ) + 4 ν d ln ε d ln t H ( ν ) (cid:19) + 124 ε (cid:18) H ( ν ) + 2 ν d ln ε d ln t H ( ν ) (cid:19) (cid:21) , (38)where the H n ( ν ) are the Hermite polynomials of order n . This expression was used by LMSV only as a checkon the validity of their linear expression. By comparing with our expression which is non-perturbative in ε ν , we will see below that these quadratic terms in fact significantly improve LMSV’s prediction.Sticking to the linearized results, we see that the expressions of both MR and LMSV have the symbolicform f ∼ e − ν / h ǫν + ǫν + ǫν + . . . i , (39)where the ellipsis denotes all terms of the type ǫν − , ǫν − , etc., as well as all terms containing ǫ . As we haveseen, deciding where to truncate the expression for f is not trivial, and using our more detailed expression wecan ask whether the expression (39) is consistent at all the relevant length scales. Immediately, we see thatthis expression cannot be correct once ǫν becomes close to unity. In Appendix B.2 we discuss the situationfor MR and LMSV at smaller mass scales. In Ref. [20], MR showed that the agreement between a Gaussian mass function calculated using the statisticsof random walks, and mass functions observed in numerical simulations with Gaussian initial conditions, canbe dramatically improved by allowing the barrier itself to perform a random walk. This approach is motivatedby the fact that the ignorance of the details of the collapse introduces a scatter in the value of the collapsethreshold for different virialized objects. The width of this scatter was found by Robertson et al . [42] to bea growing function of σ ( M ), which is consistent with the physical expectation that deviations from sphericalcollapse become relevant at small scales. The barrier can thus be treated (at least on a first approximation)as a stochastic variable whose probability density function obeys a Fokker-Planck equation with a diffusion We have corrected a typo in MR’s result [21] : the object they define as V should appear with an overall positivecoefficient in their Eqns. (85), (87) and (92). oefficient D B , which can be estimated numerically in a given N -body simulation. In particular, MR found D B ≃ . 25 using the simulations of Ref. [42].Conceptually, the variation of the value of the barrier is due to two types of effects, one intrinsicallyphysical and one more inherent to the way in which one interprets the results of simulations. From a physicalpoint of view, the dispersion accounts for deviations from the simple model of spherical collapse, for instancethe effects of ellipsoidal collapse, baryonic physics, etc. On the other hand, the details of the distribution ofthe barrier (and therefore the precise value of D B ) will depend on the halo finder algorithm used to identifyhalos in a particular simulation, since different halo finders identify collapsed objects with different properties.MR concluded that the final effect of this barrier diffusion on large scales can be accounted for in a simpleway, by changing δ c → √ aδ c where a = (1 + D B ) − . In practice this change is identical to the one proposedby Sheth et al . [24] . As MR argue in section 3.4 of Ref. [21], this barrier diffusion effect can also be accountedfor in the non -Gaussian case, again by the simple replacement of δ c → √ aδ c . It is easy to see that theirarguments go through for all our calculations as well, and we have implemented this change in our definitionof ν in Eqn. (25), setting √ a = 0 . k filter for which the results of section (3) apply. This is done by writing the 2-point func-tion h ˆ δ ( R ) ˆ δ ( R ) i calculated using the real space top-hat filter, as the Markovian value plus a correction, h ˆ δ ( R j ) ˆ δ ( R k ) i = min( t j , t k ) + ∆ jk , and noting that the correction ∆ jk remains small over the interestingrange of length scales. In fact, MR show that a very good analytical approximation for the symmetric object∆ jk , is ∆ jk ≃ κ min( t j , t k ) (cid:18) − min( t j , t k )max( t j , t k ) (cid:19) , (40)where in our case we find κ ( R ) ≃ . 464 + 0 . R , with R measured in h − Mpc. The mass function is thenobtained by perturbatively expanding in ∆ ij , with the leading effect being due to the integral Z δ c −∞ dδ . . . dδ n n X j,k =1 ∆ jk ∂ j ∂ k W gm , which on evaluation leads to f g , sharp − x ( ν, t ) = (cid:18) π (cid:19) / ν (cid:20) (1 − κ ) e − ν / + κ (cid:18) , ν (cid:19) + O ( κ ) (cid:21) , (41)where the subscript stands for Gaussian noise with the top-hat filter in real space, and κ introduces a weakexplicit t (= σ ) dependence. In Ref. [21] MR proposed an extension of this result to the non-Gaussian case,by assuming that all the non-Gaussian terms that they computed with the sharp- k filter, would simply getrescaled by the factor (1 − κ ) at the lowest order, but otherwise retain their coefficients. Symbolically, theirresult (Eqn. 88 of Ref. [21]) is f ng , sharp − x ( ν, t ) ∼ ν (cid:20) (1 − κ ) e − ν / (cid:0) ǫν + ǫν + ǫν − (cid:1) + κ (cid:18) , ν (cid:19)(cid:21) , (42)with the specific coefficients of the ǫν , ǫν and ǫν − terms being identical to those in Eqn. (36). However,the coefficient of e.g. the κǫν term arises from the action of an operator ∼ P j,k ∆ jk ∂ j ∂ k combining with thefirst unequal time operator ∼ ε t / P j ( t − t j ) ∂ j P k,l ∂ k ∂ l , and there is no simple way of predicting its exact A potential issue in this argument lies in the assumption of a linear Langevin equation for the stochastic barrier B , resulting in a simple Fokker-Planck equation with a constant D B like the one in MR, while the distribution of B was found to be approximately log-normal (and therefore non-Gaussian) in Ref. [42]. One can see that a Langevinequation of the type ∂ t B = Bξ (which would produce a log-normal distribution) can be approximated by ∂ t B = h B i ξ ,whenever the fluctuations around h B i are small, and gives a constant diffusion coefficient as long as h B i is constant.Although both approximations are reasonable on the scales of interest, non-Gaussian and scale dependent correctionsto the barrier diffusion should be studied, since in principle they could be of the same order as the other correctionsretained here. This investigation is left for future work. alue beforehand. Since MR explicitly neglect such “mixed” terms, their formula is not strictly inconsistent,as long as one keeps in mind that the theoretical error in their expression is of the same order as the terms ∼ κǫν that they include. However, if one wants to consistently retain such terms, a detailed calculation isneeded . Our calculations (not displayed) indicate that the coefficient of the κǫν term depends on certaindetails of the continuum limit of the path integral near the barrier, which require a more careful study. Weare currently investigating methods of computing these effects. At present however, we conclude that themixed terms involving both filter effects and NG, must be truncated at order ∼ κǫν .Finally, the filter-corrected mass function is also subject to effects of barrier diffusion. Here we make thesame assumptions as MR do in Ref. [20], namely that the barrier location satisfies a Langevin equation withwhite noise and diffusion constant D B , which can be accounted for by replacing κ → ˜ κ = κ/ (1 + D B ) = aκ .However, it is difficult to theoretically predict the unequal time behaviour of the barrier correlations, andthese simple assumptions must also be tested, perhaps by suitably comparing with the detailed results ofRobertson et al . [42]. We leave this for future work. Our final expression for the mass function, corrected foreffects of the diffusing barrier and the top-hat real space filter, is f ( ν, t ) = f PS ( ν ) (cid:18) − ˜ κ + O (˜ κ ) (cid:19) exp (cid:20) ε ν − (cid:16) ε − ε (cid:17) ν (cid:21) × (cid:26) − d ln ˜ κ/d ln t )1 − ˜ κ ˜ κν − (cid:0) − ν − (cid:1) − ε ν (cid:18) c − ˜ κ + 4 − c (cid:19) − ε ν − (cid:18) c − c − (cid:19) + O (˜ κ ν − , ˜ κǫν, ˜ κν − ) + O ( ǫ ν , ǫ ν , ǫν − ) (cid:27) , (43)where we have chosen to account for the scale independent O (˜ κ ) error arising from filter effects, as an overallnormalization uncertainty, and have explicitly displayed the orders of the various terms we ignore. Here f PS ( ν ) is given by Eqn. (21) with ν ( M, z ) defined in Eqn. (25).To summarize, Eqn. (43) gives an analytical expression for the non-Gaussian mass function. This expres-sion is based on approximations that are valid over a larger range of length scales than the ones presented byMR and LMSV, and incorporates effects which are ignored in the expression presented by MVJ and LMSV.Like all these other mass functions, it suffers from the errors introduced by filter effects. However, the largestof these can be accounted for as an overall normalization constant, which can be fixed using, for instance,results of a Gaussian simulation. Among the NG functions ε , ε and the c n which appear in the mass func-tion, the most important ones over the mass range of interest are ε , ε and c which are nearly constant. Inprinciple though, all these functions must be computed numerically for every mass scale of interest, and indeedall the plots in this paper use the results of such numerical calculations. However, since this is somewhattedious to do in practice, in Table 1 we provide analytical approximations for ε , ε , c , c and c , for thelocal and equilateral case as a function of σ . As mentioned earlier, all these quantities are independent ofredshift, although they depend on the choice of cosmological parameters in a complicated way in general dueto the presence of the transfer function in their definitions. However, the dependence on σ is simple, andone can check that we have ε ∝ σ , ε ∝ σ and that the c n are independent of σ . Recall that the c n are also independent of f NL and g NL . For completeness, in Table 1 we also give approximations for the filterparameters ˜ κ and d ln ˜ κ/d ln t which appear in the mass function, as functions of σ . In this section we conclude with our final results for the non-Gaussian halo mass function, comparing ourapproach with previous work. In principle, we should compare the full expressions for the mass functionsof various authors with ours. However, recall that for MVJ and LMSV one has to multiply an analytically Notice that this issue is completely decoupled from the subtleties in truncation discussed in section 4 – this problemremains even at scales where MR’s expression is formally consistent. arameter Fitting form b + c t n Local NG b c nε . . 015 0 . c . 98 0 . 073 0 . c . 15 0 . 79 0 . c . 15 0 . 45 0 . ε ( f ) − . . . ε ( g NL ) 7 . · − . . 25 Parameter Fitting form b + c t n Equilateral NG b c nε . − · − . c . − . 052 0 . c . 32 0 . 93 0 . c . 72 0 . 36 0 . κ . 36 0 . − . d ln ˜ κ/d ln t . − . − . Table 1: Analytical approximations for the various NG parameters, in the local and equilateral cases, as a functionof t = σ , in the range 2 · < M/ ( h − M sol ) < · , for f NL = 100 and g NL = 10 . We have ε ∝ f NL in bothcases, and for ε in the local case we give separate approximations for the terms proportional to f and g NL . Wedo not consider ε in the equilateral case, since the trispectrum in this case is highly model dependent. We also giveapproximations for the filter parameters ˜ κ and d ln ˜ κ/d ln t as functions of t , in the same mass range. The errors on allthe approximations are less than 1%, except for ε ( f ) where the error is ∼ predicted ratio R ng = f ( ν, M, f NL ) /f ( ν, M, f NL = 0) with a suitable Gaussian mass function based on fitsto simulations. It is not clear how to compute theoretical error bars on the latter. On the other hand, theobject R ng itself is an unambiguous theoretical prediction of every approach, that is MVJ, LMSV, MR andour work, and we can compute theoretical errors on it. In this work, we will restrict ourselves to comparingthe different expressions for R ng . In future work, we hope to compare both R ng and the full mass functionwith the results of N -body simulations.In Fig. 3 and Fig. 4 we plot the ratio R ng , respectively without and with the filter effects, at redshift z = 1. In this way we can explicitly disentangle the errors due to an approximate treatment of non-Gaussianeffects, from those due to the filter effects. We compare our expression (43) with the expressions of MR (36),LMSV (37) and (38), and of MVJ (35). Notice that, when considering the filter effects, the Gaussian functionthat enters in the ratio R ng is defined to be the function with f NL = 0 (i.e. without NG but with filtereffects when present). We use the local model, setting f NL = 100 and g NL = 0, and use the reference ΛCDMcosmology described in section 2. We do not explicitly show the final results for the equilateral model, butthey are qualitatively similar. As is commonly done in the literature, we modify the LMSV and MVJ curvesby applying the Sheth et al . correction of δ c → √ aδ c . An identical correction is already present in theexpressions (43) and (36) due to the barrier diffusion. We set √ a = 0 . 89, which is the value inferred by MR inRef. [20] using the simulations of Robertson et al . [42]. We wish to emphasize a feature that our calculationshares with that of MR, which is that the constant a is the only parameter whose value depends on the outputof N -body simulations. The rest of the calculation for the mass function is completely analytical and fromfirst principles.To make the comparison meaningful, we introduce theoretical error bars on the curves. These error barshave no intrinsic statistical meaning – they simply keep track of the absolute magnitude of the terms thatare ignored in any given prescription for computing the mass function. As we have discussed at length insection 4, these theoretical errors are scale dependent. The estimated error magnitude for each point is themaximum among the terms ignored in the expression. More explicitly, the errors for the linearized LMSVexpression (37) are estimated as the maximum of ( ǫν ) which comes from the expansion of the exponential, ǫν which is the order of the largest unequal time terms missing, and ˜ κν − which comes from the filter effects.The errors for the LMSV expression (38) are similarly estimated as the maximum of ( ǫν ) , ǫν and ˜ κν − .The largest error for the MVJ expression (35) is the maximum of ǫν (unequal time terms) and ˜ κν − (filtereffects). Finally, the error for the MR expression (36) is the maximum of ( ǫν ) from the expansion of theexponential, ǫν − from the largest unequal time terms ignored, and ˜ κ ν − and ˜ κǫν from the filter effects. Weinclude the filter effects and the associated errors only in Fig. 4. .511.522.533.5 10 R ng M (h -1 M sol )This workMRMVJLMSV quadLMSV lin Figure 3: Theoretical comparison of the different mass functions at z = 1, without the filter effects, i.e. setting˜ κ = 0. We plot the ratio R ng of the non-Gaussian and Gaussian mass functions, in the local model with f NL = 100and g NL = 0. See main text for a discussion of the error bars. The arrow indicates the mass scale where ε ν / From these figures, we can draw some interesting conclusions. First of all, we see that it is importantto retain terms which are quadratic in the NG, either with a saddle point method like in MVJ and in ourformula, or by expanding the exponential up to second order, like in LMSV. Actually, we argue that it iscorrect to keep the exponential, otherwise the expansion breaks down when ǫν is of order unity. We noticein passing that the term proportional to ε which comes from the trispectrum, partially cancels with the ε term. Secondly, comparing our expression with MVJ’s, we can observe that keeping the unequal time termsallows us to sensibly reduce the theoretical errors due to the approximate treatment of NG. In fact, if theseterms are missing, they provide the largest theoretical error on large scales. Instead, the largest theoreticalerror on small scales comes from the approximations involved in dealing with a real space top-hat filter, as isapparent from Fig. 4.To conclude, in this work we have calculated the non-Gaussian halo mass function in the excursion setframework, improving over previous calculations. We started from a path integral formulation of the randomwalk of the smoothed density field, following Ref. [19]. This allows us to take into account effects due tomulti-scale correlations of the smoothed density field (“unequal time” correlations), and due to the real spacetop-hat filter, which generates non-Markovianities in the random walk. We recognize two small parametersin which we perturb: ǫ , defined below Eqn. (24), which measures the magnitude of the primordial NG; and ν − = σ R /δ c , which is small on very large scales. In order to do a consistent expansion and to estimate thetheoretical errors, one must study the (scale dependent) relation between these two parameters, which wehave discussed in Sec. 4. We then used saddle point techniques which allowed us to non-perturbatively retainthe dependence on ǫν , which naturally appears in the calculation and whose magnitude becomes of orderunity at high masses and high redshift. Finally, we included effects due to the choice of filter function anddue to deviations from spherical collapse, as explained in Sec. 5. Our final result is presented in Eqn. (43), .511.522.533.5 10 R ng M (h -1 M sol )This workMRMVJLMSV quadLMSV lin Figure 4: Same as Fig. 3, but including filter effects. These affect only the error bars for MVJ and LMSV, and theyaffect both the curve and the error bars for MR and our result. For MR and our result, the Gaussian mass functionused to construct the ratio R ng , is taken as the non-Gaussian result at f NL = 0, and hence includes filter effects. which we reproduce here: f ( ν, t ) = f PS ( ν ) (cid:18) − ˜ κ + O (˜ κ ) (cid:19) exp (cid:20) ε ν − (cid:16) ε − ε (cid:17) ν (cid:21) × (cid:26) − d ln ˜ κ/d ln t )1 − ˜ κ ˜ κν − (cid:0) − ν − (cid:1) − ε ν (cid:18) c − ˜ κ + 4 − c (cid:19) − ε ν − (cid:18) c − c − (cid:19) + O (˜ κ ν − , ˜ κǫν, ˜ κν − ) + O ( ǫ ν , ǫ ν , ǫν − ) (cid:27) . (44)In Table 1 we provide analytical approximations for the various parameters that appear in this expression,which in general must be computed numerically. We also considered other expressions for the mass functionfound in the literature, which use different expansion methods but do not estimate the theoretical errors.We estimated the theoretical errors for each formula, and we show comparative plots in Fig. 3 and Fig.4. In our work we have improved over the calculations of MVJ [16] and LMSV [17] (who ignore unequaltime correlations) and of MR [21] (who do not retain the exponential dependence on ǫν ). We have alsodemonstrated that the (linearized) result of LMSV can be significantly improved by retaining the quadraticterms of their calculation which are usually ignored in the literature. We find that at large scales and highredshifts, the biggest theoretical errors are introduced by ignoring the exponential dependence on ǫν , followedby the neglect of unequal time correlations. The errors on our expression (44) are therefore significantly smallerthan those of the others. The strength of our approach lies in the combination of path integral methods aslaid out by MR [21], and the saddle point approximation as used by MVJ [16].Our work can be continued in several directions. First, a thorough calculation of the effects due to thechoice of the filter should be performed, since they lead to significant uncertainties in our final expression.This would include a study of the details of the continuum limit of the path integral near the barrier, and alsoa study of the statistics of the barrier diffusion process in the presence of filter effects. Second, a comparisonwith N -body simulations should be performed, in order to quantitatively assess the possibility of constrainingNG using our work. Third, it would be interesting to study how to account for the effects of ellipsoidal ollapse, in a framewrok such as the one employed in this paper. Finally, an application to the void statisticsalong the same lines should be feasible. The problem here is made more interesting by the presence of two barriers, as discussed by Sheth & van de Weygaert [43], and since voids probe larger length/mass scales thanhalos, they constitute a promising future probe of primordial NG [44]. Acknowledgements It is a pleasure to thank Stefano Borgani, Paolo Creminelli, Francesco Pace, Emiliano Sefusatti, Ravi Sheth,Licia Verde and Filippo Vernizzi for useful discussions. AppendixA Mass function calculation A.1 Equal time vs. unequal time terms In this appendix we show how the exponentiated derivative operator in the path integral can be handled byseparating the contributions of the equal time and unequal time correlations. While this calculation assumesthe toy model introduced in Eqn. (28), it easily generalizes to the more realistic case as we discuss later.Using the first few terms of the unequal time expansions, in our toy model one can write n X j,k,l =1 h ˆ δ j ˆ δ k ˆ δ l i c ∂ j ∂ k ∂ l = ε t / (cid:18) n X j,k,l =1 ∂ j ∂ k ∂ l − n X j =1 (1 − t j t ) ∂ j n X k,l =1 ∂ k ∂ l − n X j =1 (1 − t j t ) ∂ j n X k,l =1 ∂ k ∂ l + 34 n X j,k =1 (1 − t j t )(1 − t k t ) ∂ j ∂ k n X l =1 ∂ l + . . . (cid:19) , (A.1) n X j,k,l,m =1 h ˆ δ j ˆ δ k ˆ δ l ˆ δ m i c ∂ j ∂ k ∂ l ∂ m = ε t (cid:18) n X j,k,l,m =1 ∂ j ∂ k ∂ l ∂ m − n X j =1 (1 − t j t ) ∂ j n X k,l,m =1 ∂ k ∂ l ∂ m + . . . (cid:19) . (A.2)These derivative operators are exponentiated in the path integral, and act on W gm . One simplification thatoccurs in our toy model, is that the path integral in Eqn. (22) becomes a function only of ν (although thisis not obvious at this stage), and hence eventually only the ν∂ ν part of the overall derivative contributes.However, the structure of the exponentiated derivatives is still rather formidable. Moreover, the truncationof the series at this stage is based more on the intuition that higher order terms should somehow be smaller,rather than on a strict identification of the small parameters. In fact, we will see in detail in section 4 thatthe issue of truncation involves several subtleties.To make progress, it helps to analyze the effect on W gm of each of the terms in the above series, before exponentiation. The leading term in Eqn. (22) involves the multiple integral of W gm , which is just thequantity P gauss encountered in Eqn. (20). The operator ν∂ ν acts on the error function to give the Gaussianrate of Eqn. (21). Next, notice that the action of the operator P nj =1 ∂ j on any function g ( δ , . . . δ n ) underthe multiple integral, is simply Z δ c −∞ dδ . . . dδ n n X j =1 ∂ j g = ∂∂δ c Z δ c −∞ dδ . . . dδ n g , (A.3)Using this, and the fact that t / ( ∂/∂δ c ) | t = ∂ ν | t , we see that the leading term in Eqn. (A.1) (i.e. the termwith no powers of (1 − t j /t )), leads to a term involving ε ν∂ ν ( ∂ ν ) erf (cid:16) ν/ √ (cid:17) ∼ f PS ε ν (1 + O ( ν − )) , The problem with this term is that the quantity ε ν can be of order unity, and hence cannot be treatedperturbatively. To be consistent, we should keep all terms involving powers of ε ν . Luckily, this can be one in a straightforward way due to the result in Eqn. (A.3). We see that the entire exponential operatorexp[ − ( ε t / / P nj,k,l =1 ∂ j ∂ k ∂ l ] in Eqn. (22) can be pulled across the multiple integral and converted toexp [ − ( ε / ∂ ν ] acting on the remaining integral. Similarly, the operator exp [( ε t / P nj,k,l,m =1 ∂ j ∂ k ∂ l ∂ m ]can be pulled out and converted to exp [( ε / ∂ ν ], and the same applies for all such equal time operators. Wewill see later that the action of these operators can be easily accounted for, using a saddle-point approximation.To summarize, the function f at this stage is given by f = ν e − ( ε / ∂ ν +( ε / ∂ ν + ... ∂ ν Z δ c −∞ dδ . . . dδ n exp (cid:20) ε t / (cid:18) n X j =1 (1 − t j t ) ∂ j n X k,l =1 ∂ k ∂ l + 38 n X j =1 (1 − t j t ) ∂ j n X k,l =1 ∂ k ∂ l − n X j,k =1 (1 − t j t )(1 − t k t ) ∂ j ∂ k n X l =1 ∂ l + . . . (cid:19) − ε t (cid:18) n X j =1 (1 − t j t ) ∂ j n X k,l,m =1 ∂ k ∂ l ∂ m + . . . (cid:19)(cid:21) W gm . (A.4)Now consider the action of the individual terms in the remaining exponential under the integrals, but withoutexponentiation. From MR [21], we have the following results , n X j =1 (1 − t j t ) n X k,l =1 Z δ c −∞ dδ . . . dδ n ∂ j ∂ k ∂ l W gm = (cid:18) π (cid:19) / t / e − ν / , (A.5a) n X j =1 (1 − t j t ) n X k,l =1 Z δ c −∞ dδ . . . dδ n ∂ j ∂ k ∂ l W gm = (cid:18) π (cid:19) / t / h ( ν ) , (A.5b) n X j,k =1 (1 − t j t )(1 − t k t ) n X l =1 Z δ c −∞ dδ . . . dδ n ∂ j ∂ k ∂ l W gm = (cid:18) π (cid:19) / t / h ( ν ) , (A.5c)where we have defined h ( ν ) ≡ e − ν / − (cid:16) π (cid:17) / ν erfc (cid:18) ν √ (cid:19) = ν / Γ (cid:18) − , ν (cid:19) , (A.6)where Γ( − / , ν / 2) is an incomplete gamma function. Let us focus on the term in Eqn. (A.5a). If welinearize in ε in Eqn. (A.4), then this term appears with ε t / ∂ ν acting on it, leading to ∼ f PS ε ν ≪ f PS .This term can therefore be treated perturbatively. Similarly, one can check that the terms given by Eqns.(A.5b) and (A.5c) also lead to perturbatively small quantities, which are in fact further suppressed comparedto ε ν by powers of ν − . Specifically, one obtains terms involving ε erfc (cid:0) ν/ √ (cid:1) which, for large ν , reducesto ∼ f PS · ε ν · ν − (1 + O ( ν − )).A few comments are in order at this stage. First, this ordering in powers of ν − is a generic feature ofintegrals involving an increasing number of powers of (1 − t j /t ) being summed. This can be understood ina simple way from the asymptotic properties of the incomplete gamma function, as we show in Appendix C.We are therefore justified in truncating the Taylor expansion of the unequal time correlators, even thoughsuperficially (on dimensional grounds) each term in the series appears to be equally important. Secondly, wehave not yet accounted for the effect of the exponential derivatives. In fact we will see in the next sectionthat when ǫν ∼ O (1), it is these terms that impose stricter conditions on the series truncations. For now,however, we have no guidance other than the fact that if we account for one term of order ∼ ǫ n ν n , then weshould account for all terms at this order. Given this, note that for ǫν ∼ O (1) we have ν − ∼ ǫν , and hencethe terms arising from Eqns. (A.5b) and (A.5c) are of order ∼ ǫ ν . To consistently retain them, we musttherefore also retain the term linear in ε and the one quadratic in ε , when expanding the exponential. These The terms in Eqns. (A.5a), (A.5b) and (A.5c) are, upto prefactors, the integrals of what MR denote as Π (3 , NL) ,Π (3 , NNLa) and Π (3 , NNLb) respectively in Ref. [21]. nvolve the following quantities: n X j =1 (1 − t j t ) n X k,l,m =1 Z δ c −∞ dδ . . . dδ n ∂ j ∂ k ∂ l ∂ m W gm = − (cid:18) π (cid:19) / t ν e − ν / , (A.7a) n X j,k =1 (1 − t j t )(1 − t k t ) n X l,l ,l ,l =1 Z δ c −∞ dδ . . . dδ n ∂ j ∂ k ∂ l ∂ l ∂ l ∂ l W gm = − (cid:18) π (cid:19) / t ν e − ν / , (A.7b)where we have used the result (A.3), and in Eqn. (A.7b) also the identity ∂ ν h ( ν ) = − ν e − ν / . (A.8)We now see that the result of the path integral depends only on ν . Putting things together and computingthe overall ν derivative, we find the result in Eqn. (29). A.2 Saddle point calculation To compute the action of the exponentiated derivative operators, we start by writing the expression in squarebrackets in Eqn. (29) in terms of its Fourier transform, using the relations e − ν / = Z ∞−∞ dλ √ π e iλν e − λ / , − ν e − ν / = Z ∞−∞ dλ √ π ( iλ ) e iλν e − λ / ,ν e − ν / = − Z ∞−∞ dλ √ π ( λ − e iλν e − λ / , (cid:16) π (cid:17) / erfc (cid:18) ν √ (cid:19) = Z ∞−∞ dλ √ π iλ e iλν e − λ / . (A.9)Together with the identity e A ( − d/dν ) n e iλν = e A ( − iλ ) n e iλν , for constant A and B , this gives f ( ν ) = (cid:18) π (cid:19) / ν Z ∞−∞ dλ √ π e iλν e − λ / − iλ ) ε / − iλ ) ε / ... P ( λ ) (A.10)where P ( λ ) is the truncated series given by P ( λ ) = 1 + 14 iε λ + 516 iε λ − λ (cid:18) ε − ε (cid:19) + . . . (A.11)The integral in eq. (A.10) can be performed using the saddle point approximation. We write it as f ( ν ) = (cid:18) π (cid:19) / ν Z ∞−∞ dλ √ π e φ ( λ ) , (A.12)where φ ( λ ) ≡ iλν − λ + iε λ + ε λ + ln P ( λ ) + . . . (A.13)The location of the saddle point, λ = λ ∗ , is the solution of φ ′ ( λ ∗ ) = 0, and the saddle point approximationthen tells us that Z ∞−∞ dλ √ π e φ ( λ ) = e φ ( λ ∗ ) ( | φ ′′ ( λ ∗ ) | ) − / , (A.14) We are using a regulator which shifts the pole at λ = 0 in the last expression in Eqn. (A.9), to λ = − iα where α is real, positive and small. see Appendix D for a discussion of the errors introduced by this approximation). It turns out that in orderto obtain f ( ν ) correctly up to order ∼ ǫ ν , we only need λ ∗ correct up to order ∼ ǫν . The expression for φ ′ at the relevant order is, φ ′ ( λ ) = iν − λ + iε λ + . . . , (A.15)and solving for λ ∗ perturbatively up to order ǫν , we find λ ∗ = iν (cid:20) − ε ν + O ( ǫ ν ) (cid:21) . (A.16)This leads to the expression for f ( ν ) in Eqn. (30). A.3 Result with full unequal time terms In the more realistic case of slowly-varying ε n , we choose to parametrize the coefficients G and G (see Eqn.(26)) in a convenient way as follows : G (1 , , = 12 ε ( t ) c ( t ) t / ; G (2 , , = − ε ( t ) c ( t ) t − / , G (1 , , = 14 ε ( t ) c ( t ) t − / ; G (1 , , , = 12 ε ( t ) c ( t ) t , (A.17)where the coefficients c n ( t ) are smoothly varying functions and depend on the NG model. They are definedin such a way that they all reduce to unity in the toy model defined by Eqn. (28). Fig. 5 shows the behaviourof c , c and c with σ , for the local and equilateral models. The ε n and c n are independent of redshift byconstruction, since the linear growth rate D ( z ) always drops out in their definitions. Further, the c n do notdepend on the values of f NL and g NL . One can then use the definitions of ε and the c n to prove the followinguseful relations d ln ε d ln t = 32 ( c − 1) ; d ln c d ln t = 1 − c + 1 c (cid:18) c − c (cid:19) . (A.18)The calculation of the mass function for this general case proceeds completely analogously to that for the toymodel, apart from a few subtleties which we will discuss later. In this case Eqn. (A.4) is replaced with f = ( ν∂ ν − t∂ t ) e − ( ε ( t ) / ∂ ν +( ε ( t ) / ∂ ν + ... g ( ν, t )= (cid:20) ν + 13 d ln ε d ln t ε ∂ ν − d ln ε d ln t ε ∂ ν (cid:21) e − ( ε ( t ) / ∂ ν +( ε ( t ) / ∂ ν + ... ∂ ν g ( ν, t ) − te − ( ε ( t ) / ∂ ν +( ε ( t ) / ∂ ν + ... ∂ t g ( ν, t ) , (A.19)wher the function g ( ν, t ) can be shown to be g ( ν, t ) = (cid:18) π (cid:19) / (cid:20) (cid:16) π (cid:17) / erf (cid:18) ν √ (cid:19) + 14 ε c e − ν / + ε (cid:18) c − c (cid:19) h ( ν ) − ε c νe − ν / + 112 ε c νe − ν / + . . . (cid:21) , (A.20)The expression in Eqn. (A.19) can be evaluated analogously to Eqn. (29), since the additional derivatives poseno conceptual difficulty. The result of the saddle point calculation, correct up to quadratic order assuming ǫν ∼ O (1), and after using the relations (A.18), is given in Eqn. (31). .40.5480.7521.031.411.942.663.6550.1 1 c σ LocalEquilateral (a) c σ LocalEquilateral (b) c σ LocalEquilateral (c) Figure 5: The derivative coefficients c (panel (a)), c (panel (b)) and c (panel (c)), as a function of σ , for localand equilateral NG models. These quantities are independent of redshift and the NG amplitudes f NL and g NL . Theaxes are logscale. B Truncation of the perturbative series B.1 Analysis of “transition points” In this appendix we analyse the consistency of our series truncation at various mass scales. When ǫν ≃ 1, inthe polynomial in (32) we retain the terms ǫν ≃ ν − , ( ǫν − , ǫ ν ) ≃ ν − , and we discard ( ǫν − , ǫ , ǫ ν ) ≃ ν − .It would seem that our expression is then correct upto order ∼ ν − . However, the terms discarded in theexponential have the form exp( O ( ǫ ν )) ∼ exp( O ( ν − )) ∼ O ( ν − ). The error we are making is thus of thesame order as the smallest terms we are retaining, and it therefore makes sense to also ignore all the terms oforder ∼ ν − which we computed in the polynomial. The consistent expression when ǫν ≃ f ∼ e − ν ( ǫν + ǫ ν ) (cid:2) ǫν + O (cid:0) ν − (cid:1)(cid:3) . (B.1)Clearly, similar arguments can be applied at smaller scales where, e.g. one might have ǫν ≃ ν − , ν − , etc. Itis then important to ask which mass scales correspond to these “transition points”. In Fig. 6 we plot ν ( M, z )given by Eqn. (25) in an observationally interesting mass range, for three different redshifts. The horizontallines mark the transition points where ǫν becomes equal to (from top to bottom) 1, ν − , ν − , ν − , ν − and ν − . We fix ǫ = 1 / 300 which follows from the fact that in the local model with f NL = 100 we have ε ≃ . 02 (see Fig. 1), and the expression for f ( ν, M ) contains the quantity ε / ´ ´ ´ ´ ´ ´ H h - M sol L Ν Figure 6: ν ≡ δ c ( z ) /σ ( M ) in the range 5 · < ( M/h − M sol ) < · for three different redshifts, with ǫ = 1 / z = 1, 0 . ǫν becomes equal to (from top to bottom) 1, ν − , ν − , ν − , ν − and ν − . at different redshifts, and their locations also obviously depend on the value of ǫ . For example, we find thatthe transition point where ǫν ≃ ν − , remains accessible even when ǫ is an order of magnitude smaller (with ǫ ≃ / ν ≃ . ǫν ≃ , ν − on the other hand, are notaccessible for this level of NG. The transition at ǫν ≃ ν − is therefore observationally very interesting.We will now discuss in some detail the truncation of our expression for f , at various transition points.The goal is to try and settle on a single expression which is valid over a wide range of scales (i.e. acrossseveral transition points). This can then be applied without worrying about truncation inconsistencies. Ofcourse, the order of the discarded terms will then depend on the particular transition point being considered,leading to a scale dependent theoretical error. B.1.1 ǫν ≃ ν − At this transition point, the terms we retain in the exponential are ǫν ≃ ν − ; ǫ ν ≃ ν − , while discarding O ( ǫ ν ) = O ( ν − ). In the polynomial meanwhile, we retain ǫν ≃ ν − ; ǫν − ≃ ν − ; ǫ ν ≃ ν − , while discarding O ( ǫν − ) = O ( ν − ) ; O ( ǫ ) = O ( ν − ) ; O ( ǫ ν ) = O ( ν − ) . Our expression (32) therefore retains all terms correctly up to order ∼ ν − , and is consistent. With someforesight however, it turns out to be more convenient to degrade this expression somewhat by also discardingthe polynomial quadratic term ǫ ν ≃ ν − . The remaining expression, f ∼ e − ν ( ǫν + ǫ ν ) h ǫν + ǫν + O (cid:0) ν − (cid:1)i , (B.2)is also consistent at this transition point, and has a form which is identical to the ones we will see next. B.1.2 ǫν ≃ ν − As we mentioned earlier, this transition point is observationally quite interesting. The terms we retain in theexponential are ǫν ≃ ν − ; ǫ ν ≃ ν − , hile discarding O ( ǫ ν ) = O ( ν − ), and in the polynomial we retain ǫν ≃ ν − ; ǫν − ≃ ν − ; ǫ ν ≃ ν − , while discarding O ( ǫν − ) = O ( ν − ) ; O ( ǫ ) = O ( ν − ) ; O ( ǫ ν ) = O ( ν − ) . This time we see that the term ǫν − has become as important as the quadratic term ǫ ν in the polynomial,and to be consistent we should discard the quadratic term. The expansion should read f ∼ e − ν ( ǫν + ǫ ν ) h ǫν + ǫν + O (cid:0) ν − (cid:1)i . (B.3) B.1.3 ǫν ≃ ν − A similar analysis as above shows that at this stage ǫν − ≃ ν − > ǫ ν , and a consistent expression againrequires dropping the quadratic term in the polynomial, leaving f ∼ e − ν ( ǫν + ǫ ν ) h ǫν + ǫν + O (cid:0) ν − (cid:1)i . (B.4) B.1.4 ǫν ≃ ν − and smaller Beyond this point, the term ǫν − which we discard in the polynomial, becomes comparable or larger than thequadratic term of the exponential as well, and a consistent expression becomes f ∼ e − ν (1+ ǫν ) h ǫν + ǫν + . . . i (B.5)The parametric order of the terms now discarded, depends on the exact relation between ǫν and ν − .Finally, note that the error introduced by setting ν → ν g where ν g is defined using the variance of a Gaussian field, was estimated in section 3 as O ( ǫ ). When ǫν ≃ 1, this error is of order O ( ǫ ν ) and cantherefore be consistently ignored. It is not hard to see that at all lower transition points, this error continuesto be comparable to or smaller than the largest terms being discarded, and can hence be consistently ignored.This finally leads to the conclusion stated in the main text. B.2 Truncation in MR and LMSV results Since the mass scale where ǫν ≃ ǫν ≃ ν − which, as we saw, isaccessible over a wide range of redshifts for ǫ ∼ − , and at high redshifts also for ǫ ∼ − . In this case theterms MR and LMSV retain have magnitudes ǫν ≃ ν − ; ǫν ≃ ν − ; ǫν − ≃ ν − , and terms like ǫν − ≃ ν − are discarded. We know from our expression however, that ǫν appears inthe exponential, and therefore leads to terms like ( ǫν ) ≃ ν − and ( ǫν ) ≃ ν − when the exponential isexpanded, which are of the same order as the terms retained in (39). The exponential also contributes a term ǫ ν ≃ ν − , which in fact involves the tri spectrum of NG, again at the order retained by MR and LMSV.The error in the expression (39) when ǫν ≃ ν − , is therefore O ( ǫν ). (A similar analysis shows that the errorat transition point where ǫν ≃ ν − , is O ( ν − ) > O ( ǫν ).)From a purely parametric point of view, the situation for MR and LMSV improves as ν is decreasedfurther, and the expression (39) as it stands, becomes exactly consistent (in the sense discussed in the previoussubsection, see below Eqn. (33)) when ǫν ≃ ν − , because at this stage ǫν − ≃ ν − while ( ǫν ) ≃ ν − and ǫ ν ≃ ν − , and hence the exponential only contributes a single linear term ǫν . More importantly, LMSV’sexpression also has errors due to the absence of the unequal time terms discussed earlier, which are of order ǫν and can be dominant over the others. For the intermediate transitions, the analysis shows that when ǫν ≃ ν − , the error in (39) is O ( ν − ) > O ( ǫν − ), and when ǫν ≃ ν − , the error is O ( ǫν − ). This should becompared with our result (34), in which the error (at least on large scales) is always parametrically smaller than the smallest terms we retain. C Hierarchy of terms in Eqn. (A.4) Here we argue why the hierarchy of terms ordered by powers of ν − emerges on expanding the exponentiatedderivative operators in Eqn. (A.4). Focusing on terms involving the 3-point correlator, one sees that a genericterm in the expansion contains some powers of ( ε t / ), multiplying an n -dimensional integral containingsome summations ∼ P nj ,j ,... =1 (1 − t j /t ) p (1 − t j /t ) p . . . ∂ j ∂ j . . . , and also some summations over “free”derivatives ∼ P nk ,k ... =1 ∂ k ∂ k . . . , all of this acting on W gm . More precisely, the structure of the terms is ∼ ( ε t / ) m X j ,..,j m Z δ c −∞ dδ . . . dδ n [(1 − t j /t ) . . . (1 − t j m /t )] p (cid:2)(cid:0) − t j m +1 /t (cid:1) . . . (1 − t j m /t ) (cid:3) q × (cid:2)(cid:0) − t j m +1 /t (cid:1) . . . (1 − t j m /t ) (cid:3) r ∂ j . . . ∂ j m W gm , (C.1)for m ≥ p, q, r such that not all three are zero. The terms we have considered in the textare ( m, p, q, r ) = (1 , , , , , , , , , 0) and (2 , , , ∂ ν . For the “non-free” derivatives, we see thatwhat is important is the total number of (1 − t j /t ) factors accompanying these derivatives. For example,the (1 , , , 0) term in Eqn. (A.5c) has the same structure as the (1 , , , 0) term in Eqn. (A.5b) – the effectof P j,k (1 − t j /t )(1 − t k /t ) ∂ j ∂ k , up to numerical factors, is identical to that of P j,k (1 − t j /t ) ∂ j ∂ k . This isexpected to be true also with higher numbers of non-free derivatives.It is then possible to understand the hierarchy of terms by only considering terms containing P j (1 − t j /t ) p ∂ j , and no other non-free derivatives. The basic object to study now becomes X j (1 − t j /t ) Z dδ . . . dδ n ∂ j W gm , which in the continuum limit can be shown to reduce to the integral g (0) (cid:18) ν (cid:19) ≡ Z dyy / (1 − y ) / e − ν / y = √ π (cid:18) − , ν (cid:19) . (C.2)Notice the similarity with the function h ( ν ) in Eqn. (A.8), which of course is not accidental given thedefinitions of these objects. It is now easy to check that increasing the powers of (1 − t j /t ) in the summationamounts to increasing the powers of (1 − y ) in g (0) . We are then comparing (with A = ν / g (0) ( A ) with g ( p ) ( A ) where g ( p ) ( A ) ≡ Z dyy / (1 − y ) / p e − A/y . (C.3)Starting with p = 1 and manipulating the integrals, it is straightforward to establish the recurrence g ( p +1) ( A ) = g ( p ) ( A ) − Z ∞ A d ˜ A g ( p ) ( ˜ A ) . (C.4)The argument is now almost complete. We know that for large A = ν / 2, we have Γ( n, A ) = e − A A n − (1 + O ( A − )). Hence g (0) ( A ) = ( √ π/ A − / e − A (1 + O ( A − )), and its integral from A to ∞ gives a leading termproportional to Γ( − / , A ) = e − A A − / (1 + O ( A − )). The pattern is now clear: g ( p ) ( A ) ∼ A − / − p e − A (1 + O ( A − )), and since A = ν / 2, this explains the hierarchy of terms in powers of ν − , in Eqn. (A.4). The saddle point approximation In this appendix we discuss the saddle point approximation of the integrals of the type appearing in section 3.1,and estimate the error it induces. We will argue that the errors introduced by the saddle point approximationare much smaller than those due to truncating the perturbative series in the small parameters ǫ and ν − . Foran introduction to the saddle point approximation see Ref. [45]. Since we only wish to discuss the saddle pointmethod in this appendix, we will ignore here the complications introduced by the unequal time correlators,i.e. in Eqn. (A.10) we set P ( λ ) = 1. We will also work here to first order in ǫν . The extension to a moregeneral case is straightforward and the result is given by (31) as described in section 3.1. We begin withexpression (A.10): f ( ν ) = (cid:18) π (cid:19) / ν Z ∞−∞ d λ √ π e g ( λ ) , (D.1)where g ( λ ) ≡ iνλ − λ / − iλ ) ε / O ( ǫ λ ).We first find the location of a saddle point λ ∗ of the function g ( λ ), by perturbatively solving g ′ ( λ ∗ ) = 0using ǫν as the small parameter and demanding g ′′ ( λ ∗ ) < 0. The first-order solution is λ ∗ = iν (cid:0) − ε ν/ O ( ǫ ν ) (cid:1) , (D.2) g ( λ ∗ ) = − ν − ε ν + O ( ǫ ν )) , (D.3) g ′′ ( λ ∗ ) = − − ε ν + O ( ǫ ν ) . (D.4)The saddle point approximation consists roughly of performing a Taylor expansion of g ( λ ) to second orderaround λ ∗ in the integrand of (D.1) and performing the resulting Gaussian integral. We will carry this outexplicitly below. The saddle point prescription will give a good approximation to the integral as long as g ( λ )attains a global maximum at λ ∗ (along the contour of integration); this is indeed our case since the integrandin Eqn. (D.1) will be nearly a Gaussian centered at λ ∗ in the complex plane.Notice that Im λ ∗ = 0, requiring a deformation of the contour of integration such that it passes through λ ∗ . The deformation of the path of integration can be performed by taking a closed contour formed by fourpieces: The real axis C , the line Im λ = Im λ ∗ which we call here − C , and the closures of this contour atpossitive and negative infinity. The integral in this closed contour must be zero, and since the integral on theclosures of the contour at infinity can be assumed to vanish, we have R C = R C . Therefore C is the desireddeformation of the contour which passes through λ ∗ . We can then make a series of approximations in theintegral (D.1), which we discuss below, Z ∞−∞ dλ √ π e g ( λ ) ≈ Z ∞−∞ dλ √ π e g ( λ ∗ )+ g ′′ ( λ ∗ )( λ − λ ∗ ) / = e g ( λ ∗ ) ( − g ′′ ( λ ∗ )) − / = e − ν ( − ε ν/ O ( ǫ ν ) ) (cid:0) ε ν + O ( ǫ ν ) (cid:1) − / . (D.5)Here the integrations are performed along the deformed contour.In order to estimate the errors induced by the approximation done in equation (D.5), one can keep higher Technically, one should also require that Im g ( λ ) be nearly constant along the deformed contour for the saddlepoint approximation to work. In our case one can show that Im g will be suppressed by ǫ . This and all errrors inducedby the saddle point are accounted for in equation (D.6). fr ac ti on a l e rr o r M (h -1 M sol ) f NL =30 f NL =100 (a) fr ac ti on a l e rr o r M (h -1 M sol ) f NL =30 f NL =100 (b) Figure 7: Panel (a) : Fractional difference between the saddle point approximation on the r.h.s. of Eqn. (D.7) andthe numerical integration of the l.h.s. of the same equation. Panel (b) : Total error induced on the result of the toymodel (D.7) by both the saddle point approximation and the perturbative expansion to leading order in ǫν . We plotthe fractional difference between the numerical integration of the l.h.s. of Eqn. (D.7) and the approximation (D.5).Both panels show the results for a local NG with two values of f NL . orders in the Taylor expansion of the function in the exponential: Z ∞−∞ dλ √ π e g ( λ ) ≈ Z ∞−∞ dλ √ π e g ( λ ∗ )+ g ′′ ( λ ∗ )( λ − λ ∗ ) / g (3) ( λ ∗ )( λ − λ ∗ ) / g (4) ( λ ∗ )( λ − λ ∗ ) / ... ≈ e g ( λ ∗ ) ( − g ′′ ( λ ∗ )) − / + Z ∞−∞ d z √ π (cid:26) g (3) ( λ ∗ ) z + 172 (cid:2) g (3) ( λ ∗ ) (cid:3) z + 124 g (4) ( λ ∗ ) z + . . . (cid:27) e g ( λ ∗ )+ g ′′ ( λ ∗ ) z / = e g ( λ ∗ ) ( − g ′′ ( λ ∗ )) − / (1 + O ( ǫ )) . (D.6)Here we used the fact that g (3) ( λ ∗ ) = O ( ǫ ) and g (4) ( λ ∗ ) = O ( ǫ ). The integrals in the second equality of thisderivation can be computed analytically, which allows one to go to arbitrary accuracy with the saddle pointtechnique. Notice that the results of these integrations are of higher order than the terms we retain. In themain text, where the integral contains also a polynomial P ( λ ) one can again compute the errors via similarTaylor expansions. These errors can be shown to be of order O ( ǫ ), comparable to other terms which weignore.One can also estimate the errors introduced by our opproximations by using the following toy model inwhich everything is computable: Take the 3-point cumulant ε to be different from zero and all higher ordercumulants ε n for n ≥ . For such a model the integral is Z ∞−∞ d λ e iνλ − λ / − iλ ) ε / ≈ π √ ε ν ! / exp − √ ε ν + ε ν (cid:0) − √ ε ν (cid:1) ε ! . (D.7)In the r.h.s of this equation we have used the saddle point approximation but have made no expansion in ǫν .By comparing the numerical integration of the l.h.s. with the expression on the r.h.s. (panel (a) of Fig. 7),one can see that the errors introduced by the saddle point approximation are indeed of order ǫ as indicatedby (D.6). On the other hand, one can use the numerical integration of the left hand side of this equation and This toy model is inconsistent because if the third cumulant is different from zero, then all higher cumulants mustalso be different from zero. We use it here only to estimate how good the saddle point prescription is in approximatingan integral, and compare it with errors induced by a perturbative expansion in ǫν . ompare it with the approximation (D.5) (panel (b) of Fig. 7), to see that the biggest error is of order ǫ ν induced by the fact that we perform a perturbative expansion in ǫν . Notice that here we considered only theleading order in ǫν and ignored unequal time correlators, while in the main text we present a result which ismore precise (to next to leading order in ǫν ) and complete (using the excursion set formalism rigorously). References [1] N. Dalal, O. Dore, D. Huterer and A. Shirokov, “The imprints of primordial non-Gaussianities on large-scale structure: scale dependent bias and abundance of virialized objects,” Phys. Rev. D (2008)123514 [arXiv:0710.4560 [astro-ph]].[2] S. Matarrese and L. Verde, “The effect of primordial non-Gaussianity on halo bias,” Astrophys. J. (2008) L77 [arXiv:0801.4826 [astro-ph]].[3] A. Slosar, C. Hirata, U. Seljak, S. Ho and N. Padmanabhan, “Constraints on local primordial non-Gaussianity from large scale structure,” JCAP (2008) 031 [arXiv:0805.3580 [astro-ph]].[4] E. Komatsu et al. , “Seven-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Cos-mological Interpretation,” arXiv:1001.4538 [astro-ph.CO].[5] B. Sartoris, S. Borgani, C. Fedeli, S. Matarrese, L. Moscardini, P. Rosati and J. Weller, “The potentialof X-ray cluster surveys to constrain primordial non-Gaussianity,” arXiv:1003.0841 [astro-ph.CO].[6] C. Carbone, L. Verde and S. Matarrese, “Non-Gaussian halo bias and future galaxy surveys,” Astrophys.J. , L1 (2008) [arXiv:0806.1950 [astro-ph]].[7] C. Carbone, O. Mena and L. Verde, “Cosmological Parameters Degeneracies and Non-Gaussian HaloBias,” JCAP , 020 (2010) [arXiv:1003.0456 [astro-ph.CO]].[8] C. Cunha, D. Huterer and O. Dore, “Primordial non-Gaussianity from the covariance of galaxy clustercounts,” arXiv:1003.2416 [astro-ph.CO].[9] E. Sefusatti, “1-loop Perturbative Corrections to the Matter and Galaxy Bispectrum with non-GaussianInitial Conditions,” Phys. Rev. D (2009) 123002 [arXiv:0905.0717 [astro-ph.CO]].[10] R. Jimenez and L. Verde, “Implications for Primordial Non-Gaussianity ( f NL ) from weak lensing massesof high-z galaxy clusters,” Phys. Rev. D , 127302 (2009) [arXiv:0909.0403 [astro-ph.CO]].[11] L. Verde, “Non-Gaussianity from Large-Scale Structure Surveys,” arXiv:1001.5217 [astro-ph.CO].[12] V. Desjacques and U. Seljak, “Primordial non-Gaussianity from the large scale structure,”arXiv:1003.5020 [astro-ph.CO].[13] J. E. Gunn and J. R. I. Gott, “On the infall of matter into cluster of galaxies and some effects on theirevolution,” Astrophys. J. , 1 (1972).[14] W. H. Press and P. Schechter, “Formation of galaxies and clusters of galaxies by selfsimilar gravitationalcondensation,” Astrophys. J. , 425 (1974).[15] J. R. Bond, S. Cole, G. Efstathiou and N. Kaiser, “Excursion set mass functions for hierarchical Gaussianfluctuations,” Astrophys. J. , 440 (1991).[16] S. Matarrese, L. Verde and R. Jimenez, “The abundance of high-redshift objects as a probe of non-Gaussian initial conditions,” Astrophys. J. (2000) 10 [arXiv:astro-ph/0001366].[17] M. LoVerde, A. Miller, S. Shandera and L. Verde, “Effects of Scale-Dependent Non-Gaussianity onCosmological Structures,” JCAP (2008) 014 [arXiv:0711.4126 [astro-ph]].[18] R. K. Sheth and G. Tormen, “An Excursion Set Model Of Hierarchical Clustering : Ellipsoidal CollapseAnd The Moving Barrier,” Mon. Not. Roy. Astron. Soc. , 61 (2002) [arXiv:astro-ph/0105113]. 19] M. Maggiore and A. Riotto, “The Halo Mass Function from the Excursion Set Method. I. First principlederivation for the non-Markovian case of Gaussian fluctuations and generic filter,” arXiv:0903.1249 [astro-ph.CO].[20] M. Maggiore and A. Riotto, “The halo mass function from the excursion set method. II. The diffusingbarrier,” arXiv:0903.1250 [astro-ph.CO].[21] M. Maggiore and A. Riotto, “The halo mass function from the excursion set method. III. First principlederivation for non-Gaussian theories,” arXiv:0903.1251 [astro-ph.CO].[22] P. Valageas, “Mass function and bias of dark matter halos for non-Gaussian initial conditions,”arXiv:0906.1042 [astro-ph.CO].[23] J. M. Bardeen, J. R. Bond, N. Kaiser and A. S. Szalay, “The Statistics Of Peaks Of Gaussian RandomFields,” Astrophys. J. , 15 (1986).[24] R. K. Sheth, H. J. Mo and G. Tormen, “Ellipsoidal collapse and an improved model for the numberand spatial distribution of dark matter haloes,” Mon. Not. Roy. Astron. Soc. , 1 (2001) [arXiv:astro-ph/9907024].[25] T. Y. Lam and R. K. Sheth, “Halo abundances in the f NL model,” arXiv:0905.1702 [astro-ph.CO].[26] N. Sugiyama, “Cosmic background anistropies in CDM cosmology,” Astrophys. J. Suppl. , 281 (1995)[arXiv:astro-ph/9412025].[27] U. Seljak and M. Zaldarriaga, “A Line of Sight Approach to Cosmic Microwave Background Anisotropies,”Astrophys. J. (1996) 437 [arXiv:astro-ph/9603033].[28] A. Lewis, A. Challinor and A. Lasenby, “Efficient Computation of CMB anisotropies in closed FRWmodels,” Astrophys. J. (2000) 473 [arXiv:astro-ph/9911177].[29] D. Babich, P. Creminelli and M. Zaldarriaga, “The shape of non-Gaussianities,” JCAP (2004) 009[arXiv:astro-ph/0405356].[30] D. H. Lyth, C. Ungarelli and D. Wands, “The primordial density perturbation in the curvaton scenario,”Phys. Rev. D (2003) 023503 [arXiv:astro-ph/0208055].[31] N. Bartolo, S. Matarrese and A. Riotto, “On non-Gaussianity in the curvaton scenario,” Phys. Rev. D (2004) 043503 [arXiv:hep-ph/0309033].[32] G. Dvali, A. Gruzinov and M. Zaldarriaga, “A new mechanism for generating density perturbations frominflation,” Phys. Rev. D (2004) 023505 [arXiv:astro-ph/0303591].[33] M. Alishahiha, E. Silverstein and D. Tong, “DBI in the sky,” Phys. Rev. D , 123505 (2004) [arXiv:hep-th/0404084].[34] N. Arkani-Hamed, P. Creminelli, S. Mukohyama and M. Zaldarriaga, “Ghost Inflation,” JCAP ,001 (2004) [arXiv:hep-th/0312100].[35] P. Creminelli, “On non-Gaussianities in single-field inflation,” JCAP , 003 (2003) [arXiv:astro-ph/0306122].[36] P. Creminelli, A. Nicolis, L. Senatore, M. Tegmark and M. Zaldarriaga, “Limits on non-Gaussianitiesfrom WMAP data,” JCAP (2006) 004 [arXiv:astro-ph/0509029].[37] J. Zinn-Justin, “Quantum field theory and critical phenomena,” Int. Ser. Monogr. Phys. (2002) 1.[38] S. Chandrasekhar, “Stochastic problems in physics and astronomy,” Rev. Mod. Phys. , 1 (1943).[39] T. Giannantonio and C. Porciani, “Structure formation from non-Gaussian initial conditions: multivari-ate biasing, statistics, and comparison with N-body simulations,” arXiv:0911.0017 [astro-ph.CO].[40] C. Carbone et al. , “The properties of the dark matter halo distribution in non-Gaussian scenarios,” Nucl.Phys. Proc. Suppl. , 22 (2009). 41] F. Bernardeau, S. Colombi, E. Gaztanaga and R. Scoccimarro, “Large-scale structure of the universeand cosmological perturbation theory,” Phys. Rept. , 1 (2002) [arXiv:astro-ph/0112551].[42] B. Robertson, A. Kravtsov, J. Tinker and A. Zentner, “Collapse Barriers and Halo Abundance: Testingthe Excursion Set Ansatz,” Astrophys. J. (2009) 636 [arXiv:0812.3148 [astro-ph]].[43] R. K. Sheth and R. van de Weygaert, “A hierarchy of voids: Much ado about nothing,” Mon. Not. Roy.Astron. Soc. , 517 (2004) [arXiv:astro-ph/0311260].[44] M. Kamionkowski, L. Verde and R. Jimenez, “The Void Abundance with Non-Gaussian PrimordialPerturbations,” JCAP , 010 (2009) [arXiv:0809.0506 [astro-ph]].[45] A. Erd´elyi, “Asymptotic Expansions,”