{'fraction_non_alphanumeric': 0.039803570636176726, 'fraction_numerical': 0.014184921679707722, 'mean_word_length': 4.151085714285714, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An Improved Calculation of the Non-Gaussian Halo Mass Function

The abundance of collapsed objects in the universe, or halo mass function, is an important theoretical tool in 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 the smallness 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 saddle point techniques (which have been separately applied in previous approaches). Additionally, we carefully account for the various scale dependent combinations of small parameters which appear. Some of these combinations in fact become of order unity for large mass scales and at high redshifts, and must therefore be treated non-perturbatively. Our approach allows us to do this, and to also account for multi-scale density correlations which appear in the calculation. We thus derive an accurate expression for the mass function which is based on approximations that are valid over a larger range of mass scales and redshifts than those of other authors. By tracking the terms ignored in the analysis, we estimate theoretical errors for our result and also for the results of others. We also discuss the complications introduced by the choice of smoothing filter function, which we take to be a top-hat in real space, and which leads to the dominant errors in our expression. Finally, we present a detailed comparison between the various expressions for the mass functions, exploring the accuracy and range of validity of each.


Introduction
The primordial curvature inhomogeneities, generated by the inflationary mechanism, obey nearly Gaussian statistics. The deviations from Gaussianity, while expected to be small, provide a unique window into the physics 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 cosmic microwave background (CMB) temperature field, since inhomogeneities at the CMB epoch are small and their physics can be described by a perturbative treatment. In recent years, however, thanks to observations and developments in the theory, the large-scale structure (LSS) of the universe has emerged as a complementary probe to constrain primordial NG. While it is true that the n-point functions of the density field on small scales are dominated by the recent gravitational evolution, and do not reflect anymore the statistics of primordial perturbations, it turns out that the abundance of very massive objects, which form out of high peaks of the density perturbations, is a powerful probe of primordial NG. In this context, much attention has been given recently 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 sensitive probe of the magnitude of local NG f loc N L . From Ref. [3] one finds the following constraints: −29 < f loc N L < +69, already comparable with those obtained from CMB measurements in Ref. [4]: −10 < f loc N L < +74. The future is even more promising, with precisions of ∆f loc N L ∼ 10 [5] and ∆f loc N L ∼ 1 [6,7,8] being claimed for future surveys. The galaxy bispectrum is also a promising probe of NG as it could be more sensitive to other triangle configurations [9]. The mass function -which is the focus of this work, and which we discuss in detail below -has been used for example in Ref. [5] together with the scale dependent bias to produce forecasts for future surveys, and in Ref. [10] in an attempt to explain the presence of a very massive cluster at a large redshift as an indication of a large NG. For more references and information we refer the reader to reviews summarizing recent results on these topics [11,12].
The formation of bound dark matter halos from initially small density perturbations, as seen in numerical simulations, is a complicated and violent process. Some insight into the physics involved has been gained from the study of analytical models. The quantity of interest is the halo mass function, defined as the number density of dark matter halos with a mass between M and M + dM , whereρ is the average density of the universe, σ(M ) is the variance of the density contrast δ R filtered on some comoving scale R corresponding to the mass M , and the function f (σ) is to be computed. Throughout this work, we will refer to f (σ) itself as the mass function. A very useful tool in the analysis is the spherical collapse model [13], which predicts that the value of the linearly extrapolated density contrast of a spherical halo, at the time when the halo collapses, is δ c ≃ 1.686, with a weak cosmology dependence. This value serves as a collapse threshold for determining which inhomogeneous regions will end up as collapsed objects. Using this idea, Press & Schechter [14] (PS) first computed the mass function f (σ) in the case of Gaussian initial conditions. Their calculation however suffered from a problem of undercounting which affects the overall normalization -their approach does not count underdense regions embedded in larger overdense regions as eventually collapsed objects. To account for this discrepancy, PS introduced an ad-hoc factor of 2 by demanding that the mass function be correctly normalized, such that all the mass in the universe must be contained in collapsed objects. In the excursion set approach, Bond et al. [15] resolved this issue and derived a correctly normalized mass function, for Gaussian initial conditions. They argued that the filtered density contrast δ R follows a random walk as a function of the filtering scale, and the problem of computing f (σ) is translated 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, by modifying the probability density function for the (linearized) density contrast to describe non-Gaussian initial conditions. In their prescription, the relevant object is the ratio R ng of non-Gaussian to Gaussian mass functions. The full mass function is usually taken as the product of R ng and an appropriate Gaussian mass function as given by N -body simulations, e.g. the Sheth & Tormen mass function [18]. It is not clear however that 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] they pointed out that a PS-like prescription in fact misses some important non-Gaussian effects stemming from 3-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-point function of δ R , i.e. by linearizing in the non-Gaussian parameter f NL . Since the NG are assumed to be small, in the sense that the parameter ǫ = δ 3 /σ 3 satisfies ǫ ≪ 1, one might expect that such a perturbative treatment is valid. However, another crucial ingredient in the problem is that the length scales of interest are large, which leads to a second small parameter ν −1 where ν = δ c /σ. This is evident in the calculations of MR, who crucially use ν −2 ∝ σ 2 as a small parameter. Any perturbative treatment now depends not only on the smallness of ǫ and ν −1 individually, but also on the specific combinations of these parameters which appear in the calculation. It is known (and we will explicitly see below) that a natural combination that appears is ǫν 3 , which can become of order unity on scales of interest. The mass functions given by LMSV and MR therefore break down as valid series expansions when this occurs. Interestingly, MVJ's PS-like treatment on the other hand involved a saddle point approximation, allowing them to non-perturbatively account for the ǫν 3 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 calculation of the mass function. The goal of our paper is twofold. Firstly, we present a rigorous calculation of the mass function in the following way : (a) we use the techniques developed by MR in Refs. [19,20,21], which allow us to track the complex multi-scale correlations involved in the calculation, and (b) we demonstrate that MR's approach can be combined with saddle point techniques (used by MVJ), to non-perturbatively handle terms which can become of order unity. This leads to an expression for the mass function which is valid on much larger scales than those presented by MR and LMSV. Secondly, by keeping track of the terms ignored, we calculate theoretical error bars on the expressions for f (σ) resulting not only from our own calculations, but also for those of the other authors [16,17,21]. Since the terms ignored depend on ν in general, these error bars are clearly scale dependent. This allows us to estimate the validity of each of the expressions for the mass function at different scales, but importantly it also allows us to analytically compare between different expressions. 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 recent treatment of ellipsoidal collapse effects in the presence of non-Gaussianities on scales where ǫν 3 ≪ 1, see Lam & 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 and non-Gaussian N -body simulations. This probability distribution itself was approximated as a Gaussian, with mean and variance fit from simulations. While in practice this approach is easy to implement, it relies heavily on the output of N -body simulations. Our approach, on the other hand, allows us to compute a mass function almost entirely from first principles. This paper is organized as follows. In section 2 we fix some notation and briefly introduce the two most popular shapes of primordial NG, i.e. the local and equilateral ones. In section 3 we present our calculation of the mass function. In section 4 we discuss certain subtleties regarding the truncation of the perturbative series, and also compare with the other expressions for f (σ) mentioned above. In section 5 we discuss the effects induced by some additional complications introduced in the problem due to the specific choice of the filter function [19], which we take to be a top-hat in real space, and due to the inclusion of stochasticity in the value of the collapse threshold δ c (which is also expected to partially account for effects of ellipsoidal collapse) [20] . In section 6 we compare our final result Eqn. (43) with those of other authors, including theoretical 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.

Models of non-Gaussianity
We need to relate the linearly evolved density field to the primordial curvature perturbation, which carries the information of the non-linearities produced during and after inflation. We start from the Bardeen potential Φ on subhorizon scales, given by where R(k) is the (comoving) curvature perturbation, which stays constant on superhorizon scales; T (k) is the transfer function of perturbations, normalized to unity as k → 0, which describes the suppression of power for modes that entered the horizon before the matter-radiation equality; and D(z) is the linear growth factor of density fluctuations, normalized such that D(z) = (1+z) −1 in the matter dominated era. Then, the density contrast field is related to the potential by the Poisson equation, which in Fourier space reads where we substituted Eqn. (2). Here, Ω m is the present time fractional density of matter (cold dark matter and baryons), and H 0 = 100h km s −1 Mpc −1 is the present time Hubble constant. The redshift dependence is trivially accounted for by the linear growth factor D(z) and in the following, for notational simplicity, we will often suppress it. All our calculations will use a reference ΛCDM cosmology compatible with WMAP7 data [4], using parameters h = 0.702, Ω m = 0.272, present baryon density Ω b = 0.0455, scalar spectral index n s = 0.961 and σ 8 = 0.809, where σ 2 8 is the variance of the density field smoothed on a length scale of 8h −1 Mpc. For simplicity, for the transfer function T (k) we use the BBKS form, proposed in Bardeen et al. [23]: where x ≡ k(hMpc −1 )/Γ with a shape parameter Γ = Ω m h exp −Ω b (1 + √ 2h/Ω m ) that accounts for baryonic 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 qualitatively different.
In order to study halos, which form where an extended region of space has an average overdensity which is 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), where W (kR) is the Fourier transform of the filter function. For all numerical calculations we will use the spherical top-hat filter in real space, whose Fourier transform W (kR) is given by This choice allows us to have a well-defined relation between length scales and masses, namely However it introduces some complexities in the analysis, which we will comment on later. By using Eqns. (5) and (3) we have, for the 3-point function, where the subscript c denotes the connected part, and analogous formulae are valid for the higher order correlations.

Shapes of non-Gaussianity
The function R(k 1 )R(k 2 )R(k 3 ) c encodes information about the physics of the inflationary epoch. By translational invariance, it is proportional to a momentum-conserving delta function: where the (reduced) bispectrum B R (k 1 , k 2 , k 3 ) depends only on the magnitude of the k's by rotational invariance. According to the particular model of inflation, the bispectrum will be peaked about a particular shape of the triangle. The two most common cases are the squeezed (or local) NG, peaked on squeezed triangles k 1 ≪ k 2 ≃ k 3 , and the equilateral NG, peaked on equilateral triangles k 1 ≃ k 2 ≃ k 3 . Indeed, one can define a scalar product of bispectra, which describes how sensitive one is to a NG of a given type if the analysis is performed using some template form for the bispectrum. As expected, the local and equilateral shapes are approximately orthogonal with respect to this scalar product [29]. We will now describe these two models in more detail.

The local model:
The local bispectrum is produced when the NG is generated outside the horizon, for instance in the curvaton model [30,31] or in the inhomogeneous reheating scenario [32]. In these models, the curvature perturbation can be written in the following form, where R g is the linear, Gaussian field. We have included also a cubic term, which will generate the trispectrum at leading order. The bispectrum is given by where "cycl." denotes the 2 cyclic permutations of the wavenumbers, and P R (k) is the power spectrum given by P R (k) = Ak ns−4 . The trispectrum is given by The equilateral model: Models with derivative interactions of the inflaton field [33,34,35] give a bispectrum which is peaked around equilateral configurations, whose specific functional form is model dependent. Moreover, the form of the bispectrum is usually not convenient to use in numerical analyses. This is why, when dealing with equilateral NG, it is convenient to use the following parametrization, given in Ref. [36], This is peaked on equilateral configurations, and its scalar product with the bispectra produced by the realistic models cited above is very close to one. Therefore, being a sum of factorizable functions, it is the standard template used in data analyses.

Random walks and the halo mass function
We now turn to the main calculation of the paper. The non-Gaussian halo mass function can be obtained by calculating the barrier first crossing rate F of a random walk with non-Gaussian noise, in the presence of an absorbing barrier. This can be done perturbatively, starting from a path integral approach as prescribed by MR [19,21] and the mass function can be shown to be f (σ) = 2σ 2 F (σ). As discussed by MR, the calculation of f involves certain assumptions regarding the type of filter used and also the location of the barrier. In particular, 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 spherical collapse ansatz more realistic and obtain better agreement with N -body simulations, MR show that it is useful to treat the location of the barrier δ c as a stochastic variable itself, and allow it to diffuse. For the time being, 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 the calculation of the mass function. The reader is referred to Ref. [19] for a more pedagogical introduction. In the path integral approach, one treats the variance σ 2 R ≡ δ2 R as a "time" parameter, t ≡ σ 2 R , and considers the random walk followed by the smoothed density fieldδ R as this "time" is increased in discrete steps starting from small values (equivalently, as R is decreased from very large values). Hereδ R ( x) is a stochastic quantity in real space due to the stochasticity inherent in the initial conditions. We use the notationδ R to distinguish the stochastic variable from the values it takes, which will be noted by δ i below. We probe this stochasticity by changing the smoothing scale at a fixed location x = 0, thus making the variable perform a random walk, which obeys a Langevin equation with a stochastic noiseη whose statistical properties depend on the choice of filter used. In particular, for a top hat filter in k-space, the noise is white, i.e. its 2-point function is a Dirac delta [15], The random walk can be described as a trajectory {δ 0 , δ 1 , . . . , δ n } which starts withδ taking the value δ 0 = 0 at 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 crosses the barrier at δ c at a time larger than some t (i.e. at scales smaller than the corresponding R or M ), is the same as the probability that the trajectory did not cross the barrier at any time smaller than t, so that where the probability density over the space of trajectories, W ({δ j }; t) is defined as 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 = 2tF (t). In Eqn. (16) one can write the Dirac deltas using the The object e −i j λjδj is the exponential of the generating functional of the connected Green's functions, and can be shown to reduce to [37] e −i j λjδj = exp where δ j1 . . .δ jp c is the connected p-point function ofδ, with the short-hand notationδ j =δ(t j ).

Halo mass function: Gaussian case, sharp-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 δ jδk = min(t j , t k ), where min(t j , t k ) is the minimum of t j and t k . The resulting n-dimensional Gaussian integral can be handled in a straightforward way to obtain where we follow MR's notation and use the superscript "gm" to denote "Gaussian Markovian". As MR have shown [19], the resulting expression for P gauss (t) in the continuum limit ∆t → 0 is simply where we use the notation ν ≡ δ c /σ. (This in principle also includes the redshift dependence of the collapse threshold δ c , see below.) This expression for the continuum limit probability P gauss is of course a well-known result going back to Chandrasekhar [38]. This leads to the standard excursion set result for the Gaussian mass function f PS = −2t∂ t | δc P gauss , where we use the subscript PS (for Press-Schechter), to conform with the conventional notation for this object.

Halo mass function: non-Gaussian case, sharp-k filter
In the non-Gaussian case (but still retaining the sharp-k filter), the probability density W ({δ j }; t) also gets contributions from connected n-point correlators with n ≥ 3, since these in general do not vanish. These can be handled by using the relation λ k e i j λj δj = −i∂ k e i j λj δj , with ∂ j ≡ ∂/∂δ j . A straightforward calculation then shows the mass function to be where it is understood that one takes the continuum limit ∆t → 0 before computing the overall derivative with respect to t. We will find it useful to change variables from (δ c , t) to (ν, t), in which case the partial derivative becomes It is also useful at this stage to take a small detour and introduce some notation which we will use throughout the rest of the paper. We define the scale dependent "equal time" functions which as we will see, remain approximately constant over the scales of interest. We assume the ordering ε n−2 ∼ O(ǫ n−2 ) with ǫ ≪ 1, which can be motivated from their origin in inflationary physics, where one finds Typically we expect ǫ 10 −2 for f NL 100. Fig. 1 shows the behaviour of ε 1 and ε 2 in the local and equilateral models, as a function of t = σ 2 R . We see e.g. that ε 2 in the local model is comparable to ε 2 1 . In the literature one usually encounters the reduced cumulants S n , which are related to the ε n−2 by ε 1 = σS 3 , ε 2 = σ 2 S 4 and so on. The motivation for using the S n comes from the study of NG induced by nonlinear gravitational effects. However, as we see from Fig. 1, when studying primordial NG it is 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 ∼ ǫν, and we therefore require that the mass scales under scrutiny are not large enough to spoil the relation ǫν ≪ 1. It turns out that observationally interesting mass scales can nevertheless be large enough to satisfy ǫν 3 ∼ O(1). Fig. 2 shows the behaviour of ε 1 ν 3 and ε 1 ν at different redshifts, as a function of mass, in our reference ΛCDM model for local type NG, with f loc NL = 100. The behaviour for the equilateral NG is qualitatively similar. The redshift dependence of these quantities comes from the definition of ν, where we denote the usual spherical collapse threshold as δ c0 = (3/5)(3π/2) 2/3 ≃ 1.686, reserving δ c for the full, redshift dependent quantity, and a is a parameter accounting for deviations from the simplest collapse model. In the standard spherical collapse picture we have a = 1. A value of a different from unity (specifically   a ≃ 0.89) can be motivated by allowing the collapse threshold to vary stochastically [20], as we will discuss in section 5. We will soon see that the object ǫν 3 appears naturally in the calculation, and to be definite we will assume ǫν 3 ∼ O(1) for now. In section 4 we will discuss the effects of relaxing this condition and probing smaller length scales.
We now turn to the "unequal time" correlators appearing in Eqn. (22). Since we are concerned with large scales, we are in the small t limit, and following MR we expand the n-point correlators around the "final time" t. We can define the Taylor coefficients and then expand 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 3 , G 4 , etc. are algebraically rather involved, we find it useful to first consider an example in which these coefficients take simple forms. In this toy model we assume that the ε n are exactly constant, and moreover that the n-point correlators take the form 2 For clarity, we will display details of the calculation only for this model. We have relegated most of the technical details of our calculation to Appendix A. In Appendix A.1 we show that the mass function for this model can be brought to the form  where we ignore terms like ε 1 νO(ν −4 ). The remaining exponentiated derivative can be handled using a saddle point approximation. We show how to do this in Appendix A.2. The expression for the mass function f (ν) works out to which superficially at least, is comprised of two series expansions, one in the exponential and one as a polynomial, both based on the small parameter ǫν (see however the next section). This derivation assumed that ε 1 and ε 2 are constant, and that the unequal time correlations are given by Eqn. (28). It is straightforward to relax these assumptions and perform the calculation with the exact structure of the correlations. Appendix A.3 shows how to do this, and the result is where the functions c n (t) are defined in Eqn. (A.17) and characterize the behaviour of the unequal time correlations. 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 ε 1 , ε 2 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 = σ 2 of the non-Gaussian field. Computationally it is more convenient to work with the variance σ 2 g of the Gaussian field in terms of which cosmological NG are typically defined. We should then ask whether this difference will require changes in our expressions for f . We start by noting that this difference in variances is of order ∼ ǫ 2 . For example, in the local model one has σ 2 (R) = Ad 1 (R) + A(Af 2 NL )d 2 (R) where A ∼ 10 −9 is an overall normalization constant, d 1 and d 2 are scale dependent functions of comparable magnitude on all relevant scales, and ǫ is estimated as ǫ ∼ f NL A 1/2 . We therefore have ν = δ c /σ = (δ c /σ g )(1 + O(ǫ 2 )). However, with our assumption that ǫν 3 ∼ O(1), we see that this correction is actually of order ∼ (ǫ 2 ν 2 )ν −2 ∼ ǫ 3 ν 3 , which we have been consistently ignoring. We will see that even when we relax the assumption ǫν 3 ∼ O(1) and probe smaller scales where ǫν 3 ≪ 1, this correction can still be consistently ignored. Hence we can safely set ν = δ c /σ g in all of our expressions.

Consistency of the truncation 4.1 Comparative sizes of terms in the mass function
Now that all the derivative operators which we consider important have been accounted for, we can check whether our final result is consistently truncated, i.e. whether we have retained all terms at any given order in the expansion. Symbolically, our current result for the mass function can be written as with the understanding that coefficients are computed (but not displayed) for all terms except those indicated by the O() symbols. Also, ǫ 2 refers to both ε 2 1 and ε 2 . Since the expansions involve two parameters, ǫν and ν −2 , they make sense only if we additionally prescribe a relation between these parameters. So far we assumed that ǫ is fixed and ν is such that ǫν 3 ≃ 1, which was based on the observation that the term ǫν 3 naturally appears in the exponent and is not restricted in principle to small values. In Appendix B.1 we discuss this condition in more detail, and also analyse the consequences of relaxing this condition and probing smaller mass scales. We find that for observationally accessible mass scales larger than the scale where ǫν 3 ≃ ν −3 , the single expression is parametrically consistent as it stands -the terms ignored are smaller than the smallest terms retained -and in fact it remains a very good approximation even when ǫν 3 ≃ 1, since the only "inconsistent" term then is ǫν −1 , whose effect reduces as ν increases. On scales where ǫν 3 ≃ ν −4 and lower, the theoretical error becomes comparable to or larger than the quadratic term in the exponential. Plugging back all the coefficients, we have the following result for the mass function (excluding filter effects, see section 5),

Comparing with previous work
In this subsection we compare our results with previous work on the non-Gaussian mass function. As mentioned 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 due to MVJ [16], who explicitly retain the exponential dependence on ε 1 . Their expression for f can be written as 3 The major shortcoming of their result is that it is based on a Press-Schechter like prescription, and must therefore be normalized by an appropriate Gaussian mass function, typically taken to be the one due to Sheth & Tormen [18]. Additionally, it always misses the contributions due to the unequal time correlators, which contribute to the terms ∼ ǫν, ǫν −1 , etc. in Eqn. (34). When one considers formal correctness on the other hand, MR have presented a result based on explicit path integrals, which accounts for the unequal time 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 ε 1 ν 3 , means that their result is subject to significant constraints on the range of its validity. Their expression for f , ignoring filter effects, is 4 This expression is precisely what one obtains by linearizing our expression (34) in ε 1 , which serves as a check on our calculations. LMSV [17] present a result based on an Edgeworth expansion of the type encountered when studying NG generated by nonlinear gravitational effects [41]. The result most often quoted in the literature is their expression linear in ε 1 (and hence in ε 1 ν 3 ), which is In Appendix B.3 of Ref. [17], LMSV also give an expression involving ε 2 1 and ε 2 , which can be written as where the H n (ν) are the Hermite polynomials of order n. This expression was used by LMSV only as a check on the validity of their linear expression. By comparing with our expression which is non-perturbative in ε 1 ν 3 , 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 symbolic form f ∼ e −ν 2 /2 1 + ǫν 3 + ǫν + ǫ ν + . . . , where the ellipsis denotes all terms of the type ǫν −3 , ǫν −5 , etc., as well as all terms containing ǫ 2 . As we have seen, deciding where to truncate the expression for f is not trivial, and using our more detailed expression we can ask whether the expression (39) is consistent at all the relevant length scales. Immediately, we see that this expression cannot be correct once ǫν 3 becomes close to unity. In Appendix B.2 we discuss the situation for MR and LMSV at smaller mass scales.

Effects of the diffusing barrier and the filter
In Ref. [20], MR showed that the agreement between a Gaussian mass function calculated using the statistics of random walks, and mass functions observed in numerical simulations with Gaussian initial conditions, can be dramatically improved by allowing the barrier itself to perform a random walk. This approach is motivated by the fact that the ignorance of the details of the collapse introduces a scatter in the value of the collapse threshold for different virialized objects. The width of this scatter was found by Robertson et al. [42] to be a growing function of σ(M ), which is consistent with the physical expectation that deviations from spherical collapse 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 coefficient D B , which can be estimated numerically in a given N -body simulation. In particular, MR found D B ≃ 0.25 using the simulations of Ref. [42]. Conceptually, the variation of the value of the barrier is due to two types of effects, one intrinsically physical and one more inherent to the way in which one interprets the results of simulations. From a physical point of view, the dispersion accounts for deviations from the simple model of spherical collapse, for instance the effects of ellipsoidal collapse, baryonic physics, etc. On the other hand, the details of the distribution of the barrier (and therefore the precise value of D B ) will depend on the halo finder algorithm used to identify halos 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 simple way, by changing δ c0 → √ aδ c0 where a = (1 + D B ) −1 . In practice this change is identical to the one proposed by Sheth et al. [24] 5 . As MR argue in section 3.4 of Ref. [21], this barrier diffusion effect can also be accounted for in the non-Gaussian case, again by the simple replacement of δ c0 → √ aδ c0 . It is easy to see that their arguments go through for all our calculations as well, and we have implemented this change in our definition of ν in Eqn. (25), setting √ a = 0.89.
In Ref. [19], MR also accounted for the non-Markovian effects of the real space top-hat filter, as opposed to the sharp-k filter for which the results of section (3) apply. This is done by writing the 2-point function δ (R 1 )δ(R 2 ) calculated using the real space top-hat filter, as the Markovian value plus a correction, δ (R j )δ(R k ) = min(t j , t k ) + ∆ jk , and noting that the correction ∆ jk remains small over the interesting range of length scales. In fact, MR show that a very good analytical approximation for the symmetric object ∆ jk , is where in our case we find κ(R) ≃ 0.464 + 0.002R, with R measured in h −1 Mpc. The mass function is then obtained by perturbatively expanding in ∆ ij , with the leading effect being due to the integral which on evaluation leads to where the subscript stands for Gaussian noise with the top-hat filter in real space, and κ introduces a weak explicit t(= σ 2 ) 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 get rescaled by the factor (1 − κ) at the lowest order, but otherwise retain their coefficients. Symbolically, their result (Eqn. 88 of Ref. [21]) is with the specific coefficients of the ǫν 3 , ǫν and ǫν −1 terms being identical to those in Eqn. (36). However, the coefficient of e.g. the κǫν term arises from the action of an operator ∼ j,k ∆ jk ∂ j ∂ k combining with the first unequal time operator ∼ ε 1 t 1/2 j (t − t j )∂ j k,l ∂ k ∂ l , and there is no simple way of predicting its exact 5 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 DB 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 Langevin equation of the type ∂tB = Bξ (which would produce a log-normal distribution) can be approximated by ∂tB = B ξ, whenever the fluctuations around B are small, and gives a constant diffusion coefficient as long as B is constant. Although both approximations are reasonable on the scales of interest, non-Gaussian and scale dependent corrections to the barrier diffusion should be studied, since in principle they could be of the same order as the other corrections retained here. This investigation is left for future work. value 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 is needed 6 . Our calculations (not displayed) indicate that the coefficient of the κǫν term depends on certain details of the continuum limit of the path integral near the barrier, which require a more careful study. We are currently investigating methods of computing these effects. At present however, we conclude that the mixed 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 the same assumptions as MR do in Ref. [20], namely that the barrier location satisfies a Langevin equation with white 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, and these simple assumptions must also be tested, perhaps by suitably comparing with the detailed results of Robertson et al. [42]. We leave this for future work. Our final expression for the mass function, corrected for effects of the diffusing barrier and the top-hat real space filter, is where we have chosen to account for the scale independent O(κ 2 ) error arising from filter effects, as an overall normalization 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 expression is based on approximations that are valid over a larger range of length scales than the ones presented by MR 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 largest of 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 ε 1 , ε 2 and the c n which appear in the mass function, the most important ones over the mass range of interest are ε 1 , ε 2 and c 1 which are nearly constant. In principle though, all these functions must be computed numerically for every mass scale of interest, and indeed all the plots in this paper use the results of such numerical calculations. However, since this is somewhat tedious to do in practice, in Table 1 we provide analytical approximations for ε 1 , ε 2 , c 1 , c 2 and c 3 , for the local and equilateral case as a function of σ 2 . As mentioned earlier, all these quantities are independent of redshift, although they depend on the choice of cosmological parameters in a complicated way in general due to the presence of the transfer function in their definitions. However, the dependence on σ 8 is simple, and one can check that we have ε 1 ∝ σ 8 , ε 2 ∝ σ 2 8 and that the c n are independent of σ 8 . 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 filter parametersκ and d lnκ/d ln t which appear in the mass function, as functions of σ 2 .

Results and Discussion
In this section we conclude with our final results for the non-Gaussian halo mass function, comparing our approach with previous work. In principle, we should compare the full expressions for the mass functions of various authors with ours. However, recall that for MVJ and LMSV one has to multiply an analytically  Table 1: Analytical approximations for the various NG parameters, in the local and equilateral cases, as a function of t = σ 2 , in the range 2 · 10 13 < M/(h −1 M sol ) < 5 · 10 15 , for fNL = 100 and gNL = 10 4 . We have ε1 ∝ fNL in both cases, and for ε2 in the local case we give separate approximations for the terms proportional to f 2 NL and gNL. We do not consider ε2 in the equilateral case, since the trispectrum in this case is highly model dependent. We also give approximations for the filter parametersκ and d lnκ/d ln t as functions of t, in the same mass range. The errors on all the approximations are less than 1%, except for ε2(f 2 NL ) where the error is ∼ 6%. This was due to numerical difficulties in calculating this object. These approximations of course depend on our choice of cosmological parameters.
predicted ratio R ng = f (ν, M, f NL )/f (ν, M, f NL = 0) with a suitable Gaussian mass function based on fits to simulations. It is not clear how to compute theoretical error bars on the latter. On the other hand, the object R ng itself is an unambiguous theoretical prediction of every approach, that is MVJ, LMSV, MR and our work, and we can compute theoretical errors on it. In this work, we will restrict ourselves to comparing the different expressions for R ng . In future work, we hope to compare both R ng and the full mass function with 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-Gaussian effects, 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 function that enters in the ratio R ng is defined to be the function with f NL = 0 (i.e. without NG but with filter effects when present). We use the local model, setting f NL = 100 and g NL = 0, and use the reference ΛCDM cosmology described in section 2. We do not explicitly show the final results for the equilateral model, but they are qualitatively similar. As is commonly done in the literature, we modify the LMSV and MVJ curves by applying the Sheth et al. correction of δ c0 → √ aδ c0 . An identical correction is already present in the expressions (43) and (36) due to the barrier diffusion. We set √ a = 0.89, which is the value inferred by MR in Ref. [20] using the simulations of Robertson et al. [42]. We wish to emphasize a feature that our calculation shares with that of MR, which is that the constant a is the only parameter whose value depends on the output of N -body simulations. The rest of the calculation for the mass function is completely analytical and from first principles.
To make the comparison meaningful, we introduce theoretical error bars on the curves. These error bars have no intrinsic statistical meaning -they simply keep track of the absolute magnitude of the terms that are ignored in any given prescription for computing the mass function. As we have discussed at length in section 4, these theoretical errors are scale dependent. The estimated error magnitude for each point is the maximum among the terms ignored in the expression. More explicitly, the errors for the linearized LMSV expression (37) are estimated as the maximum of (ǫν 3 ) 2 which comes from the expansion of the exponential, ǫν which is the order of the largest unequal time terms missing, andκν −2 which comes from the filter effects. The errors for the LMSV expression (38) are similarly estimated as the maximum of (ǫν 3 ) 3 , ǫν andκν −2 . The largest error for the MVJ expression (35) is the maximum of ǫν (unequal time terms) andκν −2 (filter effects). Finally, the error for the MR expression (36) is the maximum of (ǫν 3 ) 2 from the expansion of the exponential, ǫν −3 from the largest unequal time terms ignored, andκ 2 ν −2 andκǫν from the filter effects. We include the filter effects and the associated errors only in Fig. 4. From these figures, we can draw some interesting conclusions. First of all, we see that it is important to retain terms which are quadratic in the NG, either with a saddle point method like in MVJ and in our formula, or by expanding the exponential up to second order, like in LMSV. Actually, we argue that it is correct to keep the exponential, otherwise the expansion breaks down when ǫν 3 is of order unity. We notice in passing that the term proportional to ε 2 which comes from the trispectrum, partially cancels with the ε 2 1 term. Secondly, comparing our expression with MVJ's, we can observe that keeping the unequal time terms allows us to sensibly reduce the theoretical errors due to the approximate treatment of NG. In fact, if these terms are missing, they provide the largest theoretical error on large scales. Instead, the largest theoretical error on small scales comes from the approximations involved in dealing with a real space top-hat filter, as is apparent from Fig. 4.
To conclude, in this work we have calculated the non-Gaussian halo mass function in the excursion set framework, improving over previous calculations. We started from a path integral formulation of the random walk of the smoothed density field, following Ref. [19]. This allows us to take into account effects due to multi-scale correlations of the smoothed density field ("unequal time" correlations), and due to the real space top-hat filter, which generates non-Markovianities in the random walk. We recognize two small parameters in which we perturb: ǫ, defined below Eqn. (24), which measures the magnitude of the primordial NG; and ν −1 = σ R /δ c , which is small on very large scales. In order to do a consistent expansion and to estimate the theoretical errors, one must study the (scale dependent) relation between these two parameters, which we have discussed in Sec. 4. We then used saddle point techniques which allowed us to non-perturbatively retain the dependence on ǫν 3 , which naturally appears in the calculation and whose magnitude becomes of order unity at high masses and high redshift. Finally, we included effects due to the choice of filter function and due to deviations from spherical collapse, as explained in Sec. 5. Our final result is presented in Eqn. which we reproduce here: 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 function found 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 unequal time correlations) and of MR [21] (who do not retain the exponential dependence on ǫν 3 ). We have also demonstrated that the (linearized) result of LMSV can be significantly improved by retaining the quadratic terms of their calculation which are usually ignored in the literature. We find that at large scales and high redshifts, the biggest theoretical errors are introduced by ignoring the exponential dependence on ǫν 3 , followed by the neglect of unequal time correlations. The errors on our expression (44) are therefore significantly smaller than those of the others. The strength of our approach lies in the combination of path integral methods as laid 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 the choice 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 also a study of the statistics of the barrier diffusion process in the presence of filter effects. Second, a comparison with N -body simulations should be performed, in order to quantitatively assess the possibility of constraining NG using our work. Third, it would be interesting to study how to account for the effects of ellipsoidal collapse, in a framewrok such as the one employed in this paper. Finally, an application to the void statistics along 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 than halos, they constitute a promising future probe of primordial NG [44].

A 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 by separating the contributions of the equal time and unequal time correlations. While this calculation assumes the 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 j,k,l=1 These derivative operators are exponentiated in the path integral, and act on W gm . One simplification that occurs in our toy model, is that the path integral in Eqn. (22) becomes a function only of ν (although this is 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 truncation of 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 that the 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 the quantity P gauss encountered in Eqn. (20). The operator ν∂ ν acts on the error function to give the Gaussian rate of Eqn. (21). Next, notice that the action of the operator n j=1 ∂ j on any function g(δ 1 , . . . δ n ) under the multiple integral, is simply Using this, and the fact that t 1/2 (∂/∂δ c )| t = ∂ ν | t , we see that the leading term in Eqn. (A.1) (i.e. the term with no powers of (1 − t j /t)), leads to a term involving The problem with this term is that the quantity ε 1 ν 3 can be of order unity, and hence cannot be treated perturbatively.
To be consistent, we should keep all terms involving powers of ε 1 ν 3 . Luckily, this can be done in a straightforward way due to the result in Eqn. (A.3). We see that the entire exponential operator exp[−(ε 1 t 3/2 /3!) n j,k,l=1 ∂ j ∂ k ∂ l ] in Eqn. (22) can be pulled across the multiple integral and converted to exp [−(ε 1 /3!)∂ 3 ν ] acting on the remaining integral. Similarly, the operator exp [(ε 2 t 2 /4!) n j,k,l,m=1 ∂ j ∂ k ∂ l ∂ m ] can be pulled out and converted to exp [(ε 2 /4!)∂ 4 ν ], and the same applies for all such equal time operators. We will 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 Now consider the action of the individual terms in the remaining exponential under the integrals, but without exponentiation. From MR [21], we have the following results 7 , where Γ(−1/2, ν 2 /2) is an incomplete gamma function. Let us focus on the term in Eqn. (A.5a). If we linearize in ε 1 in Eqn. (A.4), then this term appears with ε 1 t 3/2 ∂ ν acting on it, leading to ∼ f PS ε 1 ν ≪ f PS . This term can therefore be treated perturbatively. Similarly, one can check that the terms given by Eqns.
A few comments are in order at this stage. First, this ordering in powers of ν −2 is a generic feature of integrals involving an increasing number of powers of (1 − t j /t) being summed. This can be understood in a 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 though superficially (on dimensional grounds) each term in the series appears to be equally important. Secondly, we have not yet accounted for the effect of the exponential derivatives. In fact we will see in the next section that when ǫν 3 ∼ 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 we should account for all terms at this order. Given this, note that for ǫν 3 ∼ O(1) we have ν −2 ∼ ǫν, and hence the terms arising from Eqns. (A.5b) and (A.5c) are of order ∼ ǫ 2 ν 2 . To consistently retain them, we must therefore also retain the term linear in ε 2 and the one quadratic in ε 1 , when expanding the exponential. These involve the following quantities: where we have used the result (A.3), and in Eqn. (A.7b) also the identity We now see that the result of the path integral depends only on ν. Putting things together and computing the 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 square brackets in Eqn. (29) in terms of its Fourier transform, using the relations 8 Together with the identity e A(−d/dν) n e iλν = e A(−iλ) n e iλν , for constant A and B, this gives where P(λ) is the truncated series given by P(λ) = 1 + 1 4 iε 1 λ + 5 16 The integral in eq. (A.10) can be performed using the saddle point approximation. We write it as The location of the saddle point, λ = λ * , is the solution of φ ′ (λ * ) = 0, and the saddle point approximation then tells us that (see Appendix D for a discussion of the errors introduced by this approximation). It turns out that in order to obtain f (ν) correctly up to order ∼ ǫ 2 ν 2 , we only need λ * correct up to order ∼ ǫν. The expression for φ ′ at the relevant order is, and solving for λ * perturbatively up to order ǫν, we find 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 3 and G 4 (see Eqn. (26)) in a convenient way as follows : where the coefficients c n (t) are smoothly varying functions and depend on the NG model. They are defined in such a way that they all reduce to unity in the toy model defined by Eqn. (28). Fig. 5 shows the behaviour of c 1 , c 2 and c 3 with σ 2 , for the local and equilateral models. The ε n and c n are independent of redshift by construction, since the linear growth rate D(z) always drops out in their definitions. Further, the c n do not depend on the values of f NL and g NL . One can then use the definitions of ε 1 and the c n to prove the following useful relations The calculation of the mass function for this general case proceeds completely analogously to that for the toy model, apart from a few subtleties which we will discuss later. In this case Eqn. (A.4) is replaced with wher the function g(ν, t) can be shown to be In this appendix we analyse the consistency of our series truncation at various mass scales. When ǫν 3 ≃ 1, in the polynomial in (32) we retain the terms ǫν ≃ ν −2 , (ǫν −1 , ǫ 2 ν 2 ) ≃ ν −4 , and we discard (ǫν −3 , ǫ 2 , ǫ 3 ν 3 ) ≃ ν −6 . It would seem that our expression is then correct upto order ∼ ν −4 . However, the terms discarded in the exponential have the form exp . The error we are making is thus of the same order as the smallest terms we are retaining, and it therefore makes sense to also ignore all the terms of order ∼ ν −4 which we computed in the polynomial. The consistent expression when ǫν 3 ≃ 1 is then given by Clearly, similar arguments can be applied at smaller scales where, e.g. one might have ǫν 3 ≃ ν −1 , ν −2 , etc. It is 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 horizontal lines mark the transition points where ǫν 3 becomes equal to (from top to bottom) 1, ν −1 , ν −2 , ν −3 , ν −4 and ν −5 . We fix ǫ = 1/300 which follows from the fact that in the local model with f NL = 100 we have ε 1 ≃ 0.02 (see Fig. 1 at different redshifts, and their locations also obviously depend on the value of ǫ. For example, we find that the transition point where ǫν 3 ≃ ν −2 , remains accessible even when ǫ is an order of magnitude smaller (with ǫ ≃ 1/3000, this transition occurs at ν ≃ 4.96). The transitions at ǫν 3 ≃ 1, ν −1 on the other hand, are not accessible for this level of NG. The transition at ǫν 3 ≃ ν −2 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. across several transition points). This can then be applied without worrying about truncation inconsistencies. Of course, the order of the discarded terms will then depend on the particular transition point being considered, leading to a scale dependent theoretical error.
At this transition point, the terms we retain in the exponential are . In the polynomial meanwhile, we retain Our expression (32) therefore retains all terms correctly up to order ∼ ν −6 , and is consistent. With some foresight however, it turns out to be more convenient to degrade this expression somewhat by also discarding the polynomial quadratic term ǫ 2 ν 2 ≃ ν −6 . The remaining expression, is also consistent at this transition point, and has a form which is identical to the ones we will see next.
B.1.2 ǫν 3 ≃ ν −2 As we mentioned earlier, this transition point is observationally quite interesting. The terms we retain in the exponential are while discarding O(ǫ 3 ν 5 ) = O(ν −10 ), and in the polynomial we retain This time we see that the term ǫν −3 has become as important as the quadratic term ǫ 2 ν 2 in the polynomial, and to be consistent we should discard the quadratic term. The expansion should read A similar analysis as above shows that at this stage ǫν −3 ≃ ν −9 > ǫ 2 ν 2 , and a consistent expression again requires dropping the quadratic term in the polynomial, leaving

B.1.4 ǫν 3 ≃ ν −4 and smaller
Beyond this point, the term ǫν −3 which we discard in the polynomial, becomes comparable or larger than the quadratic term of the exponential as well, and a consistent expression becomes The parametric order of the terms now discarded, depends on the exact relation between ǫν 3 and ν −1 .
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(ǫ 2 ). When ǫν 3 ≃ 1, this error is of order O(ǫ 3 ν 3 ) and can therefore be consistently ignored. It is not hard to see that at all lower transition points, this error continues to 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 ǫν 3 ≃ 1 is on the border of the observed mass window (for galaxy cluster observations), even at high redshifts, let us therefore directly look at the case ǫν 3 ≃ ν −2 which, as we saw, is accessible over a wide range of redshifts for ǫ ∼ 10 −2 , and at high redshifts also for ǫ ∼ 10 −3 . In this case the terms MR and LMSV retain have magnitudes and terms like ǫν −3 ≃ ν −8 are discarded. We know from our expression however, that ǫν 3 appears in the exponential, and therefore leads to terms like (ǫν 3 ) 2 ≃ ν −4 and (ǫν 3 ) 3 ≃ ν −6 when the exponential is expanded, which are of the same order as the terms retained in (39). The exponential also contributes a term ǫ 2 ν 4 ≃ ν −6 , which in fact involves the trispectrum of NG, again at the order retained by MR and LMSV. The error in the expression (39) when ǫν 3 ≃ ν −2 , is therefore O(ǫν). (A similar analysis shows that the error at transition point where ǫν 3 From a purely parametric point of view, the situation for MR and LMSV improves as ν is decreased further, and the expression (39) as it stands, becomes exactly consistent (in the sense discussed in the previous subsection, see below Eqn. (33)) when ǫν 3 ≃ ν −5 , because at this stage ǫν −1 ≃ ν −9 while (ǫν 3 ) 2 ≃ ν −10 and ǫ 2 ν 4 ≃ ν −12 , and hence the exponential only contributes a single linear term ǫν 3 . More importantly, LMSV's expression 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 ǫν 3 ≃ ν −3 , the error in (39) is O(ν −6 ) > O(ǫν −1 ), and when ǫν 3 ≃ ν −4 , the error is O(ǫν −1 ). This should be compared 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 ν −2 emerges on expanding the exponentiated derivative operators in Eqn. (A.4). Focusing on terms involving the 3-point correlator, one sees that a generic term in the expansion contains some powers of (ε 1 t 3/2 ), multiplying an n-dimensional integral containing some summations ∼ n j1,j2,.
. ., and also some summations over "free" derivatives ∼ n k1,k2...=1 ∂ k1 ∂ k2 . . ., all of this acting on W gm . More precisely, the structure of the terms is for m ≥ 1 and non-negative p, q, r such that not all three are zero. The terms we have considered in the text are (m, p, q, r) = (1, 1, 0, 0), (1, 1, 1, 0), (1, 2, 0, 0) and (2, 1, 0, 0). We have already discussed how the "free" derivatives can be pulled out of the integral and converted to ∂ ν . For the "non-free" derivatives, we see that what is important is the total number of (1 − t j /t) factors accompanying these derivatives. For example, the (1, 1, 1, 0) term in Eqn. (A.5c) has the same structure as the (1, 2, 0, 0) term in Eqn. (A.5b) -the effect This is expected 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 j (1 − t j /t) p ∂ j , and no other non-free derivatives. The basic object to study now becomes j (1 − t j /t) dδ 1 . . . dδ n ∂ j W gm , which in the continuum limit can be shown to reduce to the integral Notice the similarity with the function h(ν) in Eqn. (A.8), which of course is not accidental given the definitions of these objects. It is now easy to check that increasing the powers of (1 − t j /t) in the summation amounts to increasing the powers of (1 − y) in g (0) . We are then comparing (with A = ν 2 /2) g (0) (A) with g (p) (A) where

D 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 approximation are much smaller than those due to truncating the perturbative series in the small parameters ǫ and ν −1 . For an introduction to the saddle point approximation see Ref. [45]. Since we only wish to discuss the saddle point method 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 more general case is straightforward and the result is given by (31) as described in section 3.1. We begin with expression (A.10): where g(λ) ≡ iνλ − λ 2 /2 + (−iλ) 3 ε 1 /6 + O(ǫ 2 λ 4 ). We first find the location of a saddle point λ * of the function g(λ), by perturbatively solving g ′ (λ * ) = 0 using ǫν as the small parameter and demanding g ′′ (λ * ) < 0. The first-order solution is The saddle point approximation consists roughly of performing a Taylor expansion of g(λ) to second order around λ * in the integrand of (D.1) and performing the resulting Gaussian integral. We will carry this out explicitly 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 integrand in 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 four pieces: The real axis C 1 , the line Im λ = Im λ * which we call here −C 2 , and the closures of this contour at possitive and negative infinity. The integral in this closed contour must be zero, and since the integral on the closures of the contour at infinity can be assumed to vanish, we have C1 = C2 . Therefore C 2 is the desired deformation of the contour which passes through λ * 9 . We can then make a series of approximations in the integral (D.1), which we discuss below,  Here we used the fact that g (3) (λ * ) = O(ǫ) and g (4) (λ * ) = O(ǫ 2 ). The integrals in the second equality of this derivation can be computed analytically, which allows one to go to arbitrary accuracy with the saddle point technique. Notice that the results of these integrations are of higher order than the terms we retain. In the main text, where the integral contains also a polynomial P(λ) one can again compute the errors via similar Taylor expansions. These errors can be shown to be of order O(ǫ 2 ), comparable to other terms which we ignore. One can also estimate the errors introduced by our opproximations by using the following toy model in which everything is computable: Take the 3-point cumulant ε 1 to be different from zero and all higher order cumulants ε n for n ≥ 2 to be zero 10 . For such a model the integral is  Fig. 7), to see that the biggest error is of order ǫ 2 ν 2 induced by the fact that we perform a perturbative expansion in ǫν. Notice that here we considered only the leading order in ǫν and ignored unequal time correlators, while in the main text we present a result which is more precise (to next to leading order in ǫν) and complete (using the excursion set formalism rigorously).