Correlation Functions in Stochastic Inflation

Combining the stochastic and $\delta N$ formalisms, we derive non perturbative analytical expressions for all correlation functions of scalar perturbations in single-field, slow-roll inflation. The standard, classical formulas are recovered as saddle-point limits of the full results. This yields a classicality criterion that shows that stochastic effects are small only if the potential is sub-Planckian and not too flat. The saddle-point approximation also provides an expansion scheme for calculating stochastic corrections to observable quantities perturbatively in this regime. In the opposite regime, we show that a strong suppression in the power spectrum is generically obtained, and comment on the physical implications of this effect.


Introduction
Inflation is one of the leading paradigms describing the physical conditions that prevailed in the very early Universe [1][2][3][4][5][6]. It is a phase of accelerated expansion that solves the puzzles of the standard hot Big Bang model, and it provides a causal mechanism for generating scalar [7][8][9][10][11] and tensor [12] inhomogeneous perturbations on cosmological scales. These inhomogeneities result from the parametric amplification of the vacuum quantum fluctuations of the gravitational and matter fields during the accelerated expansion.
The transition from these quantum fluctuations to classical but stochastic density perturbations [13][14][15][16] gives rise to the stochastic inflation formalism [17][18][19][20][21][22][23][24][25]. 1 It consists of an effective theory for the long-wavelength parts of the quantum fields, which are "coarse grained" at a fixed physical scale (i.e. non-expanding), somewhat larger than the Hubble radius during the whole inflationary period. 2 The non-commutative parts of this coarse grained field ϕ are small, and at this scale, short-wavelength quantum fluctuations have negligible non-commutative parts too. In this framework, they act as a classical noise on the dynamics 1 This formalism was, in fact, first used in Ref. [9] at the level of the Langevin equation, from which results lying beyond the one-loop approximation for the inflaton field were obtained. 2 More precisely, the coarse grained part of the field consists of the modes k for which k σaH. Here, σ is a cutoff parameter satisfying [25] e −1/(3 1 ) σ 1, where 1 is the first slow-roll parameter. Under this condition, the physical results are independent of σ.
of the super-Hubble scales, and ϕ can thus be described by a stochastic classical theory, following the Langevin equation (1.1) our choice of combining the stochastic and δN formalisms. We then settle our computational strategy and proceed with the calculation itself. Results are presented in section 4; see in particular Eqs. (4.1) and (4.6). We show that the standard formulas are recovered in a "classical" limit that we carefully define, and discuss the regimes where they are not valid. Finally, in section 5, we summarize our main results and conclude.

Time Variable Issue
Because of the Friedmann equation, the Hubble parameter appearing in Eq. (1.1) is sourced by the inflaton field itself, through the slow-roll function H(ϕ). At leading order in the noise, one simply has H(ϕ) = H(φ c ), which is a classical quantity. Beyond the leading order, however, H is dependent on the full coarse grained field and is therefore a stochastic quantity. 4 This has two consequences. The first one is that starting from a classical time label, any other time variable defined through a or H is a stochastic quantity, and cannot be used to label the Langevin equation, otherwise one would describe a physically different process. The time label must therefore be carefully specified. The second one is that, since H is related to the curvature of space-time, its stochasticity has to do with the one of space-time itself. We are thus a priori describing effective quantum gravitational effects, corresponding to the gravitational-and self-interactions of the inflaton field. The corresponding corrections should therefore remain small as long as the energy density of the inflaton field is small compared to the Planck scale. For this reason, it is convenient to define the dimensionless potential which we will make use of extensively in the following. Before turning to the calculation of the correlation functions, in this section, we show on a simple example why different time labels in the Langevin equation typically yield results that differ by ∝ v corrections.
First, let us recast the stochastic process (1.1) through a Fokker-Planck equation, which governs the time evolution of the probability density P (φ, N ) that ϕ = φ at time N . In the Itô interpretation 5 [17,35,36], it reads 6 3) 4 Hereafter, by "stochastic quantities", we simply refer to realization dependent quantities, as opposed to quantities that are fixed for all realizations. 5 More generally, the last term in Eq. (2.3) can be written in the form with 0 ≤ α ≤ 1, where α = 0 corresponds to the Itô interpretation and α = 1/2 to the Stratonovich one [34]. However, analysis shows that keeping terms explicitly depending on α exceeds the accuracy of the stochastic approach in its leading approximation (1.1). In particular, corrections to the noise term due to self-interactions of small-scale fluctuations (if they exist) are at least of the same order or even larger. 6 Note also that we never use the "volume weighted" variant of Eq. (2.3) proposed as an alternative in Ref. [24] since then the resulting distribution is not normalizable: its integral over dφ is time-or N -dependent. Thus, it leads to probability non-conservation. Neither is it justified from the physical point of view, since it is based on the assumption that all Hubble physical volumes ("observers") emerging from the expansion of a previous inflationary patch are clones of each other, while they are strongly correlated. Now let us compare this equation with the one that would have been obtained if the Langevin equation had been written in terms of cosmic time t. Performing the simple change of time variable dN = Hdt in Eq. (1.1), this is given by Here we use the notationφ to stress the fact that, a priori,φ does not describe the same stochastic process as ϕ. The Fokker-Planck equation corresponding to Eq. (2.4) is given by If H is taken to be a function of time only, independent of ϕ, the H factors can be taken out of the derivatives with respect to φ in Eqs. (2.3) and (2.5). In this case, it is straightforward to see that these two are perfectly equivalent through the change of time variable dN = Hdt, and that they describe the same stochastic process. On the contrary, if H explicitly depends on ϕ, this is obviously no longer the case and P =P . This can be better illustrated by calculating the stationary distributions associated with these processes. Let P stat (φ) denote a stationary probability distribution for the stochastic process (1.1), or equivalently, (2.3). By definition, which defines the probability current J. This current thus needs to be independent of φ for a stationary distribution. In most interesting situations, it is actually 0. This is notably the case when the allowed values for φ are unbounded. For example, if V (φ) is defined up to φ = ∞, the normalization condition P stat (φ)dφ = 1 requires that P stat (φ) decreases at infinity strictly faster than |φ| −1 . In this case, both P stat (φ) and ∂P stat (φ)/∂φ vanish at infinity. From Eq. (2.6), J vanishes at infinity also, hence everywhere. This yields a simple differential equation to solve for P stat (φ), and one obtains .
Here, an overall integration constant, which makes the distribution normalized, P stat (φ)dφ = 1, is omitted. Similarly, Eq. (2.5) can be written as ∂P /∂t = ∂J/∂φ, and requiring that the currentJ vanishes gives rise to a differential equation for the stationary distributionP stat (φ), which can easily be solved. One obtains .
The two distributions are close, and the effects coming from the H(ϕ) dependence are small, only in the regions of the potential where v 1. At this point, we are left with the issue of identifying the right time variable to work with. Actually, one can explicitly show [26,27,37] that N is the correct answer, and that it is the only time variable that allows the stochastic formalism to reproduce a number of results from QFT on curved space-times. We leave this discussion to appendix A, where we elaborate on existing results and show why, since we deal with metric perturbations, we must work with N .

Method
Let us now review how correlation functions of curvature fluctuations can be calculated in stochastic inflation, and see which approach is best suited to the issue we are interested in.
The problem can first be treated at linear order [38][39][40] by expanding the coarse grained field ϕ about its classical counterpart φ cl , ϕ = φ cl + δφ (1) . Here, recall that φ cl is the solution of the Langevin equation (1.1) without the noise term. The quadratic moment of δφ (1) can be calculated as in appendix A.2, see Eq. (A.23). It corresponds to the integrated power spectrum of the field fluctuations on super-Hubble scales, and can therefore be related [40] to the power spectrum P ζ of curvature perturbations thanks to the relation In this expression, the right hand side needs to be evaluated when the scale associated with the wavenumber k (at which the power spectrum is calculated) exits the Hubble radius. If one plugs the expression (A.23) obtained in appendix A for δφ (1) 2 using N as the time variable into Eq. (3.1), one obtains .

(3.2)
As before, φ cl needs to be evaluated when the scale associated with the wavenumber at which the power spectrum is calculated exits the Hubble radius. The quantity 1 ≡ −dH/(H 2 dt) is the first slow-roll parameter. At leading order in slow roll, it verifies 1 = M 2 Pl /2(v /v) 2 . The above expression exactly matches the standard result [41,42]. In order to get the first corrections to this standard result, one thus needs to go to higher orders in δφ. Actually, one can show that no contributions arise at next-to-leading order, and that one needs to go at least to next-to-next-to-leading order. This renders the calculation technically difficult. This is why we will prefer to make use of non-perturbative techniques. In passing, let us stress that in Ref. [40], the Langevin equation is written and solved with t, whereas, as already said, the number of e-folds N must be used instead. This has important consequences. Indeed, if one makes use of cosmic time t and plugs the associated expression (A.29) for the quadratic moment of δφ (1) into Eq. (3.1), one obtains Here, we have adopted the same notation as in section 2 where a tilde recalls that not the same quantity is worked out andζ is not ζ. This result matches Eq. (2.11) of Ref. [40]. However, in this work, it is concluded that, because of the second term in the braces of Eq. (3.3), which is always negative, the amplitude of the power spectrum in the stochastic approach is in general reduced with respect to the standard result. One can see that such a statement is incorrect, since the extra term in Eq. (3.3) is simply due to not working with the correct time variable. This is why, if such an approach were to be followed and extended to higher orders, it would again be crucial to work with N as the time variable.
Another strategy is followed in Refs. [43][44][45], where methods of statistical physics, such as replica field theory, are employed in a stochastic inflationary context. However, only the case of a free test field evolving in a de Sitter or power-law background is investigated, while we need to go beyond the fixed background assumption in order to study the effects of the explicit H(ϕ) dependence. This is why we cannot directly make use of this computational scheme in the present work.
Finally, in Refs. [33,46,47], the δN formalism is used to relate the curvature perturbations to the number of e-folds statistics. This is this last route that we chose to follow here, since it does not rely on any perturbative expansion scheme, and since it does not prevent us from implementing the explicit H(ϕ) dependence. In Ref. [47], numerical solutions are obtained for quadratic and hybrid potentials. In the present work, we derive fully analytical and non-perturbative results that apply to any single-field potential, and which do not require a numerical solution of the Langevin equation. As a by-product, this allows us to prove, for the first time, that the standard results are always recovered in the classical limit, for any potential.

The δN Formalism
The δN formalism [9,[48][49][50][51][52] is very well suited to addressing the calculation of correlation functions in stochastic inflation, since it relates the statistical properties of curvature perturbations to the distribution of the number of e-folds among a family of homogeneous universes. Let us first recall where this correspondence comes from and, as an example, how the scalar power spectrum is usually calculated in the associated formalism.
Starting from the unperturbed flat Friedmann-Lemaître-Robertson-Walker (FLRW) line element, ds 2 = −dt 2 + a 2 (t)δ ij dx i dx j , deviations from homogeneity and isotropy can be included in a more general, perturbed metric, which contains some gauge redundancy. A specific gauge choice consists in requiring that fixed t slices of space-time have uniform energy density, and that fixed x worldlines be comoving. When doing so, and including scalar perturbations only, the perturbed metric in this gauge (which coincides in the super-Hubble regime with the synchronous gauge supplemented by some additional conditions fixing it uniquely) becomes [9,53,54] ds 2 = −dt 2 + a 2 (t)e 2ζ(x) δ ij dx i dx j , up to small terms proportional to gradients of ζ. Here, ζ is the adiabatic (curvature) perturbation, which is time-independent in single-field inflation once the decaying mode can be neglected. The omission of tensor perturbations is justified by the fact that their amplitude is suppressed compared to the scalar ones by the small slow-roll parameter 1 . This allows us to define a local scale factorã(t, x) = a(t)e ζ(x) . Starting from an initial flat slice of space-time at time t in , the amount of expansion N (t, x) ≡ ln [ã(t, x)/a(t in )] to a final slice of uniform energy density is then related to the curvature perturbation through where N 0 (t) ≡ ln [a(t)/a(t in )] is the unperturbed amount of expansion. From this, an important simplification arises on large scales where anisotropy and spatial gradients can be neglected, and the local density, expansion rate, etc., obey the same evolution equations as a homogeneous FLRW universe. Thus we can use the homogeneous FLRW solutions to describe the local evolution, which is known as the "quasi-isotropic" [55][56][57][58] or "separate universe" approach [51,59,60]. It implies that N (t, x) is the amount of expansion in unperturbed, homogeneous universes, so that ζ can be calculated from the knowledge of the evolution of a family of such universes. Written in terms of the inflaton field φ(x) = φ + δφ(x), consisting of an unperturbed, homogeneous piece φ and of a perturbation δφ originating from quantum fluctuations, Eq. (3.4) gives rise to Here, N is to be evaluated in unperturbed universes from an initial epoch when the inflaton field has an assigned value φ to a final epoch when the energy density has an assigned value ρ. Since the observed curvature perturbations are almost Gaussian, at leading order in perturbation theory, one has Here, N (φ) is usually evaluated with the slow-roll, classical formula . It can be expressed in terms of the power spectrum of δφ (defined by similar relations) thanks to Eq. (3.6). For quasi-de Sitter inflation, and when the curvature of the inflaton potential is much smaller than H, on super-Hubble scales, the latter is given by [61] P δφ (k) H 2 (k)/4π 2 , where H(k) means H evaluated at the time when the k mode crosses the Hubble radius, i.e. when aH = k. Together with Eq. (3.7), one therefore obtains which is the same as Eq. (3.2) and which matches the standard result [41,42].
A fundamental remark is that in the above usual calculation, the quasi-isotropic (separate universe) approximation is assorted with the assumption that on super-Hubble scales, the evolution of the inflaton field is governed by its classical equation of motion (3.7). The stochastic dispersion in the number of e-folds thus only comes from the field dispersion at Hubble crossing δφ * . In most cases, this is a good approximation for the following reason. From the Langevin equation (1.1), one can see that during the typical time scale of one e-fold, the classical drift of the inflaton field is of the order ∆φ cl = V /(3H 2 ) = √ 2 1 M Pl , while the quantum kick is of the order ∆φ qu = H/(2π). This allows us to define a rough "classicality" criterion ∆φ qu /∆φ cl that assesses the amplitude of the stochastic corrections to the classical trajectory. Making use of Eqs. (3.8), this ratio can be expressed as which is valid for single-field slow-roll models of inflation with canonical kinetic terms. Since P ∼ 2×10 −9 for the modes observed in the Cosmic Microwave Background (CMB), stochastic effects are already small when these modes cross the Hubble radius. If one further assumes that 1 monotonously grows toward 1 during the last stages of inflation, P ζ ∝ H 2 / 1 decreases (since H can only decrease) and one is therefore ensured that the stochastic corrections to the inflaton trajectory remain small. However, they are two caveats to this line of reasoning. The first one is that, as we will show below, ∆φ qu /∆φ cl is not the correct way to assess the importance of stochastic effects and one should use instead another classicality criterion that we will derive. The second one is that, in some situations, 1 becomes tiny or even vanishes in some transient phase between the Hubble exit time of the observed modes and the end of inflation. This is the case, for example, when the potential has a flat inflection point, such as in MSSM inflation [62][63][64] or as in punctuated inflation [65,66]. Another situation of interest is when inflation does not have a graceful exit but ends due to tachyonic instability involving an auxiliary field, like in hybrid inflation [67,68], or by brane annihilation in string-theoretical setups [69,70]. In such cases, 1 can decrease and the last e-folds of inflation may be dominated by the quantum noise. It is therefore important to study the dispersion δN arising not only from δφ * but from the complete subsequent stochastic history of the coarse grained field.
Note also that in these expressions, ζ need not be small as was shown in Refs. [9,51,71] [note, however, that ζ is defined up to a constant due to an arbitrary possible rescaling of a(t)], thus, δN need not be small, too. As follows from the quasi-isotropic (separate universe) approach, the condition for inflation to proceed is only that H M Pl . On the other hand, if P ζ (k) ∼ H/(M Pl √ 1 ) exceeds unity (the so called regime of "eternal inflation"), then the Universe loses its local homogeneity and isotropy after the end of inflation, but not immediately. This occurs much later than the comoving scale a(t)/k at which this inhomogeneity occurs crosses the Hubble radius H −1 second time. Thus, in the scope of the inflationary scenario P ζ may well exceed unity at scales much exceeding the present Hubble radius. The stochastic inflation approach provides us with a possibility to obtain quantitatively correct results in this non-linear regime, too.

Computational Programme
This is why we now generalize this approach to a fully stochastic framework. For a given wavenumber k, let φ * (k) be the mean value of the coarse grained field when k crosses the Hubble radius. If inflation terminates at φ end , let N (k) denote the number of e-folds realized between φ * (k) and φ end . Obviously, N is a stochastic quantity, and we can define its variance (3.10) It is related with the curvature perturbation δN of Eq. (3.6) in the following manner. Since δN is computed between two fixed points φ * (k) and φ end , it receives an integrated contribution from all the modes crossing the Hubble radius between these two points, and one has Here we have used the relation N (k) = ln(a end /a * ) = ln(k end /k)(1 + 1 * + · · · ), where 1 * + · · · stand for slow-roll corrections that we do not need to take into account at leading order in slow roll. One then has In the same manner, the third moment of the number of e-folds distribution,  the power spectrum squared, is then given by (3.14) where 5/72 is a conventional historical factor. Analogously, the trispectrum is related to the third derivative of δN 4 with respect to N , and so on and so forth. The computational programme we must follow is now clear. For a given mode k, we first calculate φ * (k) (this sets the location of the observational window). We then consider stochastic realizations of Eq. (1.1) that satisfy ϕ = φ * (k) at some initial time, 7 and denote by N the number of e-folds that is realized before reaching φ end . Among these realizations, we calculate the first moments of this stochastic quantity, N , N 2 , N 3 , etc. We finally apply relations such as Eqs. (3.12) and (3.14) to obtain the power spectrum, the non-Gaussianity local parameter, or any higher order correlation function.

First Passage Time Analysis
In what follows, this calculation is performed using the techniques developed in "first passage time analysis" [72,73], which was applied to stochastic inflation in Ref. [17]. We consider the situation sketched in Fig. 1, where the inflaton is initially located at φ * and evolves in some potential V (φ) under Eq. (1.1). Because any part of the potential can a priori be explored, here we consider two possible ending points, φ 1 and φ 2 , located on each side of φ * . If the potential is, say, of the hilltop type (left panel), φ 1 and φ 2 can be taken at the two values where inflation has a graceful exit, on each side of the maximum of the potential. If, on the other hand, a flat potential extends up to φ = ∞ (right panel), one of these points, say φ 2 , can be taken where V becomes super-Planckian and inhomogeneities prevent inflation from occurring. In such cases, the precise value of φ 2 plays a negligible role, as we will show in section 3.3.1. Let N be the number of e-folds realized during this process.
Before proceeding with the calculation of the N moments, a first useful result to establish is the Itô lemma, which is a relation verified by any smooth function f of ϕ. The Taylor expansion of such a function at second order is given by ϕ is a realization of the stochastic process under study, dϕ is given by Eq. (1.1) and at first order in dN , one obtains which we will repeatedly make use of in the following.

Ending Point Probability
As a first warm-up, let us calculate the probability p 1 that the inflaton field first reaches the ending point located at φ 1 [i.e. φ (N ) = φ 1 ], or, equivalently the probability p 2 = 1 − p 1 that the inflaton field first reaches the ending point located at φ 2 [i.e. φ (N ) = φ 2 ]. This will also allow us to determine when the ending point located at φ 2 plays a negligible role. First of all, let ψ (ϕ) be a function of the coarse grained field that can be expressed as where h (ϕ) will be specified later. By construction, one has ψ (φ 1 ) = 1 and ψ (φ 2 ) = 0. This implies that the mean value of ψ evaluated at ϕ (N ) is given by The idea is then to find an h (ϕ) function that makes easy the evaluation of the left hand side of the previous relation, so that we can deduce p 1 . In order to do so, let us apply the Itô lemma (3.16) to h (ϕ). If one requires that the integral of the second line of Eq. (3.16) vanishes, that is, Because h and ψ are linearly related, see Eq. (3.17), the same equation is satisfied by ψ. When averaged over all realizations, 8 its right hand side vanishes. One then obtains ψ [ϕ (N )] = ψ (φ * ), which is the probability p 1 one is seeking for. All one needs to do is therefore to solve Eq. (3.18) to obtain h (ϕ), to plug the obtained expression into Eq. (3.17) to derive ψ(ϕ), and finally to evaluate this function at φ * . A formal solution to Eq. (3.18) is given by where A and B are two integration constants that play no role, since they cancel out when calculating ψ thanks to Eq. (3.17). Indeed, the latter gives rise to and a symmetric expression for p 2 . 9 A few remarks are in order about this result. First, one can check that, since φ * lies between φ 1 and φ 2 , the probability (3.20) is ensured to be between 0 and 1. Second, one can also verify that when φ * = φ 1 , p 1 = 1, and when φ * = φ 2 , p 1 = 0, as one would expect. Third, in the case depicted in the right panel of Fig. 1, in the limit where φ 2 → ∞, one is sure to first reach the ending point located at φ 1 , that is, Indeed, the numerator of the expression for p 2 is finite, since a bounded function is integrated over a bounded interval. If the potential is maximal at φ 2 , and if it is monotonous over an interval of the type [φ 0 , φ 2 [, its denominator is on the contrary larger than the integral of a function bounded from below by a strictly positive number, over an unbounded interval [φ 0 , φ ∞ [. This is why it diverges, and why p 2 vanishes. This means that if φ 2 is sufficiently large, its precise value plays no role, since inflation always terminates at φ 1 = φ end .

Mean Number of e-folds
Let us now turn to the calculation of the mean number of e-folds N . As above, we want to make use of the Itô lemma (3.16). To do so, let us define f (ϕ) as the solution of the differential equation with boundary conditions f (φ 1 ) = f (φ 2 ) = 0. Such a solution will be explicitly calculated in due time. For now, it is interesting to notice that when this is plugged into the Itô lemma (3.16), the first term of the left hand side, f (φ 1 or φ 2 ), vanishes, and the second integrand of the right hand side is −1. Thus, the Itô equation can be rewritten as Integration domain of Eq. (3.24) when evaluated at ϕ = φ 2 , in the case φ 1 < φ 2 (the opposite case proceeds the same way). The discrete parameter x is integrated between φ 1 and φ 2 , while y varies between x andφ. The resulting integration domain is displayed in green. When x <φ, one has dxdy > 0 and one integrates a positive contribution to the mean number of e-folds. Conversely, when x >φ, one has dxdy < 0 and one integrates a negative contribution. This is necessary in order for the overall integral to vanish. This is whyφ must lie between φ 1 and φ 2 .
By averaging over realizations, one obtains 10 What one needs to do is therefore to solve the deterministic differential equation (3.21) with the associated boundary conditions, and to evaluate the solution at φ * . One obtains whereφ is an integration constant set to satisfy the condition f (φ 2 ) = 0. There is no generic expression for it, 11 but one can be more specific. First of all, as can be seen in Fig. 2,φ must be such that, when f is evaluated at φ 2 , the integration domain of Eq. (3.24) possesses a positive part and a negative part, which are able to compensate for each other. This implies thatφ must lie between φ 1 and φ 2 . A second generic condition comes from splitting the The first integral vanishes because 10 Here again, since both the integrand and the upper bound are stochastic quantities, it is non-trivial that the integral in the right hand side of Eq. (3.22) vanishes when averaged, but it can be shown rigorously. 11 Alternatively, one can write Eq. (3.24) in the explicit form [17] f where p1 is given by Eq. (3.20) and, in the configuration of Fig. 1, θ(x − x * ) = 1 when x > x * and 0 otherwise.
f (φ 2 ) = 0, which means that in order for f to be symmetrical in φ 1 ↔ φ 2 ,φ(φ 1 , φ 2 ) must satisfy this symmetry too, that is to say,φ (φ 1 , φ 2 ) =φ (φ 2 , φ 1 ). Third, in the case where the potential is symmetric about a local maximum φ max close to which inflation proceeds, the integrand in Eq. (3.24) is symmetric with respect to the first bisector in Fig. 2. The two green triangles must therefore have the same surface, which readily leads toφ = φ max . Fourth, finally, in the case displayed in the right panel of Fig. 1, if φ 2 is sufficiently large, we have established in section 3.3.1 that p 2 = 0 and the quantity we compute is the mean number of e-folds between φ * and φ 1 = φ end . For explicitness, let us assume that v > 0 (the same line of arguments applies in the case v < 0). Inflation proceeds at φ < φ 2 . In the domain of negative contribution in Fig. 2, the argument of the exponential in Eq. (3.24) is positive. As a consequence, ifφ is finite and φ 2 → ∞, the negative contribution to the integral is infinite while the positive one remains finite, which is impossible. In order to avoid this, one must then haveφ = φ 2 . In practice, almost all cases boil down to one of the two previous ones andφ is specified accordingly. Combining Eqs. (3.23) and (3.24), one finally has 12 This quantity is plotted for large and small field potentials in Fig. 3, where it is compared with the results of a numerical integration of the Langevin equation (1.1) for a large number of realizations over which the mean value of N is computed. One can check that the agreement is excellent.

Classical Limit
Let us now verify that the above formula boils down to the classical result (3.7) in some "classical limit". This can be done by performing a saddle-point approximation of the integrals appearing in Eq. (3.25). Let us first work out the y-integral, that is to say, φ x dy/v(y) exp[1/v(y)]. Since the integrand varies exponentially with the potential, the strategy is to evaluate it close to its maximum, i.e. where the potential is minimum. The potential being maximal atφ in most cases (see the discussion above), the integrand is clearly maximal 13 which exactly matches the classical result (3.7). The classical trajectory thus appears as a saddle-point limit of the mean stochastic trajectory, analogously to what happens e.g. in the context of path integral calculations. This calculation also allows us to identify under which conditions the classical limit is recovered. A priori, the Taylor expansion of 1/v can be trusted as long as the difference between 1/v(x) and 1/v(y) is not too large, say |1/v(y) − 1/v(x)| < R, where R is some  small number. If one uses the Taylor expansion of 1/v at first order, this means that |y − x| < Rv 2 /v . Requiring that the second order term of the Taylor expansion is small at the boundary of this domain yields the condition |2v − v v 2 /v 2 | 1. For this reason, we define the classicality criterion This quantity is displayed in the top axes in Fig. 3 and one can check that indeed, the classical trajectory is a good approximation to the mean stochastic one if and only if η cl 1. In the following, we will see that η cl is the relevant quantity to discuss the strength of the stochastic effects in general and in section 4.3, we will further discuss the physical implications of Eq. (3.27).
For now, and for future use, let us give the first correction to the classical trajectory. This can be obtained going one order higher in the saddle-point approximation, that is to say, using a Taylor expansion of 1/v at second order. One obtains where the dots stand for higher order terms. In the brackets of Eq. (3.28), the two last terms stand for the first stochastic correction and one should not be surprised that, in general, when η cl is small, it is small. It is also interesting to notice that it is directly proportional to d 1 /dN . When 1 increases as inflation proceeds, the stochastic leading correction is therefore positive and the stochastic effects tend to increase the realized number of e-folds, while when 1 decreases as inflation proceeds, the correction is negative and the stochastic effects tend to decrease the number of e-folds, at least at linear order.

Number of e-folds Variance
Let us now move on with the calculation of the dispersion in the number of e-folds, defined in Eq. (3.10). If one squares Eq. (3.22), and takes the stochastic average of it, one obtains 15 In order to make use of the Itô lemma, let then g(φ) be the function defined by where f is the function defined in Eq. (3.24). When the Itô lemma (3.16) is applied to g (φ), if one further sets g (φ 1 ) = g (φ 2 ) = 0, one obtains where the second equality is just a consequence of Eq. (3.29) and where the third equality is just a consequence of Eq. (3.23). Therefore, one just needs to solve Eq. (3.30) with boundary conditions g (φ 1 ) = g (φ 2 ) = 0 and to evaluate the resulting function at φ * in order to obtain δN 2 = g(φ * ). The differential equation (3.30) can formally be integrated, and one obtains is an integration constant that must be chosen in order to have g (φ 2 ) = 0. One can show that it satisfies the four properties listed in section 3.3.2 forφ and can therefore be specified in the same manner. With φ 1 = φ end , one then has . (3.33)

Classical Limit
As was done for the mean number of e-folds in section 3.3.2, let us derive the classical limit of Eq. (3.33). Obviously, in the classical setup the trajectories are not stochastic and δN 2 = 0, and what we are interested in here is the non-vanishing leading order contribution to δN 2 in the limit η cl 1. As before, the y-integral can be worked out with a saddle-point approximation, and one obtains 16 Plugging back this expression into Eq. (3.33), one obtains Finally, and for future use again, let us give the first correction to this classical limit. Going one order higher in the saddle-point approximation, one obtains (3.35)

Number of e-folds Skewness and Higher Moments
In the same manner, if one denotes the third moment of the distribution of number of efolds by m (φ * ) = δN 3 defined in Eq. (3.13), one can show that m(φ) is the solution of the differential equation m − m v /v 2 = −6f g that obeys m(φ 1 ) = m(φ 2 ) = 0. As before, taking φ 1 = φ end and φ 2 = φ ∞ , one obtains whereφ 3 can be set asφ. Similarly to above, making use of Eqs. (3.28) and (3.35), a saddle-point approximation of this integral leads to the classical limit (3.37) Let us finally explain how the same procedure can be iterated and higher order moments can be calculated. Let us denote the p th momentum of the number of e-folds distribution by (3.38) where, by convention, we set σ 0 = 1 and σ 1 = 0. As above, one can recursively show that σ p is the solution of the differential equation .
(3.40) When p = 2, this yields the variance (3.33); when p = 3, the skewness (3.36) is obtained; when p = 4, the kurtosis could be derived as well, and so on and so forth for any moment. 16 In the η cl 1 limit, f is close to the classical trajectory (3.26) as shown in section 3.3.2, and one can take f (y) v (y) /v (y) M −2 Pl .

Results
We are now in a position where we can combine the intermediary results of the previous sections to give explicit, non-perturbative and fully generic expressions for the first correlation functions of curvature perturbations in stochastic inflation. We first derive the relevant formulas and their classical limits, before commenting on their physical implications in section 4.3.

Power Spectrum
Following the programme we settled in section 3.2, if one plugs Eqs. (3.25) and (3.33) into Eq. (3.12), one obtains P ζ (φ * ) = g (φ * )/f (φ * ), that is, . (4.1) In this expression, P ζ (φ * ) stands for the power spectrum calculated at a scale k such that when it crosses the Hubble radius, the mean inflaton field value is φ * . This formula provides, for the first time, a complete expression of the curvature perturbations power spectrum calculated in stochastic inflation. It is plotted for large and small field potentials in Fig. 4.
From this, a generic expression for the spectral index can also be given. Since, at leading order in slow roll, ∂/∂ ln(k) −∂φ/∂ N × ∂/∂φ, one has Here, for conciseness, we do not expand this expression in terms of integrals of the potential, but it is straightforward to do so with Eqs. (3.24) and (3.32).

Classical Limit
Before commenting further on the physical implications of the above result, let us make sure that in the classical limit, η cl 1, the standard formula is recovered. Combining Eqs. (3.12), (3.26) and (3.34), one has which exactly matches the usual result (3.8) at leading order in slow roll. This fully generalizes the work of Ref. [46], where this agreement is shown but only in the case where the Hubble parameter varies linearly with φ and if the noise has constant amplitude. Here we have extended this result to any potential is single-field slow-roll inflation, and included dependence of the noise amplitude on the coarse grained field. In the same manner, making use of Eqs. (3.26), (3.34) and (4.2) together, one obtains n S | cl = 1 − M 2 Pl 3 (v /v) 2 − 2v /v , which again matches the standard slow-roll result n S = 1 − 2 1 − 2 where 2 ≡ d ln 1 /dN is the second slow-roll parameter, since at leading order in slow roll, one has 2 = 2M 2 Pl (v 2 /v 2 − v /v), and 1 has been given above.
Let us now derive the leading order corrections to these "classical" results. This can be done making use of Eqs. (3.28) and (3.35) in Eqs. (4.1) and (4.2). For the power spectrum, one obtains while for the spectral index one gets (4.5)

Non-Gaussianity and Higher Moments
The local non-Gaussianity parameter can be calculated in the same manner, and Eq. (3.14) gives rise to For conciseness, this expression is not expanded in terms of integrals of the potential, but it is straightforward to do so with Eqs. (3.24), (3.32) and (3.36).
Here also, we need to make sure that in the classical limit, the standard result is recovered. Combining Eqs. (4.6), (3.28), (3.35) and (3.37), one obtains The first two terms in the brackets match the usual result [75]. In contrast, it is important to stress that within the usual δN formalism, the standard result cannot be obtained because of the intrinsic non-Gaussianity of the fields at Hubble exit [75,76]. Such effects are automatically taken into account in our formalism, which readily gives rise to the correct formula.
Obviously, one can go on and calculate any higher order correlation function with Eq. (3.40). However, with the power spectrum and non-Gaussianity local parameter at hand, we already are in a position where we can draw important physical conclusions.

Discussion
A first important consequence of Eqs. (4.1) and (4.6) is the correctness of their classical limits. They show the validity of our computational programme for calculating correlation functions in general. This may be particularly useful for investigating other cases than singlefield slow-roll inflation, especially when the standard procedure is difficult to follow. Indeed, our method can easily be numerically implemented, and it could then be applied to more complicated scenarios such as multi-field inflation where it has been shown [77] that the δN formalism retains reliable, modified kinetic terms where the stochastic inflation formalism has been generalized [78][79][80][81], etc. In particular, it is well suited to situations where stochastic effects dominate the inflationary dynamics in some parts of the potential [47,68,82] and where one must take the stochastic effects into account.
However, even if the stochastic effects within the CMB observable window need to be small, let us stress that the location of the observable window along the inflationary potential can be largely affected. This notably happens when the potential has a flat region between the location where the observed modes exit the Hubble radius and the end of inflation, as is the case e.g. in hybrid inflation or in potentials with flat inflection points.
Another point to note is that, contrary to what one may have expected, the corrections we obtained are not controlled by the ratio ∆φ qu /∆φ cl extensively used in the literature, but by the classicality criterion η cl derived in Eq. (3.27). This has two main consequences.
First, η cl has dimension v, which means that it is Planck suppressed. 17 This makes sense, since, as noted in section 2, some of the corrections we obtained physically correspond to the self-and gravitational interactions of the inflaton field. 18 This is why it can be Table 1. Classicality criterion η cl defined in Eq. (3.27) for a few types of inflationary potentials. Except for "large field", the expression given for η cl is valid close to the flat point of the potential. useful to compare our results with loop calculations performed in the literature by means of other techniques. In particular, the self-loop correction to the power spectrum is derived in Ref. [86], and graviton loop corrections are obtained in Ref. [87] (for a nice review, see also Ref. [88]). A diagrammatic approach based on the δN formalism is also presented in Ref. [89] where the power spectrum and the bispectrum are calculated up to two loops. In all these cases, the obtained corrections are of the form P 1−loop α is a numerical factor of order one that depends on the kind of loops one considers, and 2 stands for second order combinations of slow-roll parameters. When the number of e-folds N is of the order 1/ , this is exactly the kind of leading corrections we obtained. This feature is therefore somewhat generic. Obviously, it remains to understand which loops exactly our approach allows one to calculate, and how our results relate to the above mentioned ones. We leave it for future work. Second, η cl contains 1/v 2 terms. This means that, even if v needs to be very small, 19 if the potential is sufficiently flat, η cl may be large. In table 1, we have summarized the shape of η cl for different prototypical inflationary potentials. For large field potentials, η cl is directly proportional to v. This is why, in the left panels of Figs. 3 and 4, departure of the stochastic results from the standard formulas occur only when v 1, in a regime where our calculation cannot be trusted anyway. However, for potentials with flat points, different results are obtained. If the flat point is of the hilltop type, η cl diverges at the maximum of the potential. This is why, in the right panels of Figs. 3 and 4, even if v saturates to a small maximal value, the stochastic result differs from the classical one close to the maximum of the potential. However, this happens many e-folds before the scales probed in the CMB cross the Hubble radius, that is to say, at extremely large, non-observable scales. The same conclusion holds for plateau potentials (either of the polynomial or exponential type) where stochastic effects lead to non-trivial modifications in far, non observable regions of the plateau. On the other hand, if the potential has a flat inflection point, η cl can be large at intermediate wavelengths, too small to lie in the CMB observable window but still of astrophysical interest. This could have important consequences in possible non-linear effects at those small scales, such as the formation of primordial black holes (PBHs) [90]. In such models, the production of PBHs is calculated making use of the standard classical formulas for the amount of scalar perturbations. However, we have shown that in such regimes, stochastic effects largely modify its value. An important question is therefore how this changes the production of PBHs in these models. In particular, it is interesting to notice that if the potential is concave (v < 0), which is the case favored by observations [91][92][93], the leading correction in Eq. (4.4) is an enhancement of the power spectrum amplitude. However, as can be seen in the right panel of Fig. 4, as soon as one leaves the perturbative regime, this is replaced by the opposite trend: at the flat point, the classical result accounts for a diverging power spectrum while the stochastic effects make it finite. If generic, this effect may be important for the calculation of PBHs formation, and we plan to address this issue in a future publication.

Conclusion
Let us now summarize our main results. Making use of the δN formalism, we have shown how curvature perturbations can be related to fluctuations in the realized amount of inflationary e-folds in stochastic inflation trajectories. We have then applied "first passage time analysis" techniques to derive all the statistical moments of the number of e-folds, hence all scalar correlation functions in stochastic inflation.
We have shown that the standard results can be recovered as saddle-point limits of the full expressions. The situation is therefore analogous to, e.g., path integral calculations. A new simple classicality criterion has been derived, which should replace the common estimate based on the ratio between the mean quantum kick and the classical drift during one efold. It shows that quantum corrections to inflationary observables are Planck suppressed in general (that is to say, they are proportional to V /M 4 Pl ), but can be large if the potential is flat enough, even at sub-Planckian scales. For simple inflationary models where |dV /dφ|/V increases monotonously as inflation proceeds, the corresponding effects play a non-trivial role only at extremely large, non-observable scales. However, models containing a flat point in the potential between the Hubble exit location of the modes currently observed in the CMB and the end of inflation behave differently. First, the stochastic effects change the mean total number of inflationary e-folds and can therefore largely modify the location of the observational window along the inflationary potential. Second, the amount of scalar perturbations produced around the flat point is strongly modified by stochastic effects. This may be crucially important for a number of non-linear effects computed at these small scales, such as the formation of PBHs, or non-Gaussianity.
Together with the case of tensor perturbations, which we have not addressed in this paper, we plan to study these issues in future publications.
A Why must we use the number of e-folds in the Langevin equation?
In section 2, we have made explicit that different choices of time variable in the Langevin equation account for different stochastic processes. In this appendix, we explain why the correct time variable to work with is the number of e-folds N , elaborating on already existing results. We first present a generic argument, based on a perturbation of the background equations, before turning to an explicit comparison of stochastic inflation and QFT predictions, in order to identify the correct time variable.

A.1 Perturbations Equation derived from the Background Equation
A heuristic derivation [17] of the Langevin equation relies on splitting the full quantum inflaton field into a coarse grained, large scales part ϕ, and a short-wavelength component φ > , and on performing an expansion of the equation of motion in φ > . As suggested in Ref. [37], the correct time variable should therefore be the one such that the equations for the perturbations, which must be established at the action level, can correctly be obtained from varying the equation of motion for the background itself, when written in terms of this time variable. In this section, we establish that this condition selects out N as the time variable.
In the case where inflation is driven by a single scalar field φ, the action we start from is given by From this action (and this action only), we first want to derive equations of motion for the scalar perturbations that can be compared with what will be obtained below from varying the background equation of motion itself. To make our point even more convincing, we go up to second order in the perturbations. This is why we expand the background fields {φ, g µν } at second order in the scalar perturbations. 20 When the time variable in the metric is the conformal time η, one has The degrees of freedom introduced above are partially redundant and in absence of anisotropic stress, the scalar sector can be described in terms of a single gauge invariant variable. One possible choice is the Mukhanov-Sasaki variable [7,42,94] v, which can be defined, order by order, as the scalar field fluctuation φ (n) on uniform curvature hypersurfaces [95]. To first and second orders, after a lengthy but straightforward calculation, one obtains [96] v (1)  The equation of motion for the scalar perturbations φ (1) and φ (2) is therefore given by the one for v (1) and v (2) in this gauge. At leading order in the slow-roll approximation, and in the long-wavelength limit, they read 21 (A.7) We now need to compare these equations with the ones that arise when varying the equation of motion for the background, and to find for which time variable they match.
If t is used When cosmic time t is used, the leading order of the slow-roll approximation for the Klein-Gordon equation for the background is given by where we take H 2 V /(3M 2 Pl ) at leading order in slow roll. When plugging φ = φ (0) + φ (1) + φ (2) in this equation, one obtains at first and second order in the perturbations One should stress that these equations do not apply to φ (1) and φ (2) since they are different from Eqs. (A.6) and (A.7), which is why we use the notationφ (1,2) instead of φ (1,2) . The differences with Eqs. (A.6) and (A.7) are displayed in red. One can see that several factors do not match. This is because in general, the equations for the perturbations must be derived from the action itself and cannot be obtained by simply varying the equation of motion for the background.
If ds = H p a q dt is used For this reason, let us look for a time variable s which is such that the equations for the perturbations arise from varying the equation of motion of the background when written in terms of s. Let us assume that s is related to t thanks to a relation of the form ds = H p (φ) a q (φ) dt , (A.11) where p and q are power indices that we try to determine. For example, when p = 0 and q = 0, s is the cosmic time t, when p = 1 and q = 0, s is the number of e-folds N , while when p = 0 and q = −1, s is the conformal time η. In terms of s, the equation of motion for the background is given by where again we take H 2 V /(3M 2 Pl ) at leading order in slow roll. When plugging in φ = φ (0) + φ (1) + φ (2) , one obtains at first and second order in the perturbations (A.14) Again, these equations do not apply to φ (1) and φ (2) in general, since the correct ones are given by Eqs. (A.6) and (A.7) which is why we use the notationφ (1,2) . The differences between these two sets of equations are displayed in red. In order for the above to match Eqs. (A.6) and (A.7), one must have q = 0 and (p + 1)/6 = 1/3, which gives p = 1, (p + 1)/2 = 1, which also gives p = 1, and (p + 1)(p + 3)/36 = 2/9, which gives p = 1 or p = −5. As a conclusion, with p = 1 and q = 0 only, the equations for the perturbations (from what is shown here, up to second order in perturbation theory) can be seen as if they were derived from varying the equation of motion for the background. This choice corresponds to the number of e-folds N .

A.2 Stochastic Inflation and Quantum Field Theory on Curved Space-Times
To go beyond this generic argument, one can explicitly show [26,27] that N is the time variable which allows one to consistently connect stochastic inflation with results from QFT on curved space-times. For example, let us consider the leading order of the fluctuations δφ = ϕ − φ cl in the coarse grained inflaton field about its classical background value φ cl . By "classical", recall that we mean that φ cl is the solution of the equation of motion without the noise term. We want to compute the mean square value of δφ and compare what we obtain