Multiple Fields in Stochastic Inflation

Stochastic effects in multi-field inflationary scenarios are investigated. A hierarchy of diffusion equations is derived, the solutions of which yield moments of the numbers of inflationary $e$-folds. Solving the resulting partial differential equations in multi-dimensional field space is more challenging than the single-field case. A few tractable examples are discussed, which show that the number of fields is, in general, a critical parameter. When more than two fields are present for instance, the probability to explore arbitrarily large-field regions of the potential, otherwise inaccessible to single-field dynamics, becomes non-zero. In some configurations, this gives rise to an infinite mean number of $e$-folds, regardless of the initial conditions. Another difference with respect to single-field scenarios is that multi-field stochastic effects can be large even at sub-Planckian energy. This opens interesting new possibilities for probing quantum effects in inflationary dynamics, since the moments of the numbers of $e$-folds can be used to calculate the distribution of primordial density perturbations in the stochastic-$\delta N$ formalism.


Introduction
Inflation is the leading paradigm to describe 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 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, that later seed the large scale structure of the Universe.

JCAP06(2016)043
grained" at a fixed physical scale (i.e. non-expanding), larger than the Hubble radius during the whole inflationary period. The non-commutative parts of these coarse grained fields are small (compared to their anti-commutative parts), and at this scale, short-wavelength quantum fluctuations have negligible non-commutative parts too. In this framework, they behave as a classical noise acting on the dynamics of the super-Hubble scales, and the coarse grained fields can thus be described by a stochastic classical theory, following Langevin equations.
The stochastic formalism accounts for the quantum modification of the super-Hubble scales dynamics, and allows us to study how quantum effects modify inflationary observable predictions. Stochastic inflation is indeed a powerful tool for calculating correlation functions of quantum fields during inflation [27][28][29][30][31][32][33][34][35][36][37][38][39][40] (see also ref. [41,42] for the case of scalar electrodynamics during inflation and ref. [43] for the case of derivative interactions and constrained fields). In practice, it can be connected to cosmological observations through the δN formalism [9,[44][45][46][47][48][49][50], which relates the statistical properties of scalar curvature perturbations ζ to the distribution of the number of e-folds N among a family of homogeneous universes. In the stochastic formalism, this number of e-folds (realised between an initial flat slice of space-time and a final slice of uniform energy density) is a stochastic quantity that we denote N , and its statistical moments directly give rise to correlation functions of cosmological perturbations. For example, the power spectrum of curvature perturbations can be expressed as where N is the number of e-folds realised between the time when the scale at which P ζ is calculated exits the Hubble radius during inflation and the end of inflation, and similar expressions can be derived for higher correlation functions. This is the so-called "stochastic-δN formalism" [51][52][53][54][55]. In fact, this approach may also be called "stochastic-N formalism" since it does not rely on an expansion in δN and in the metric perturbation ζ (for instance, these two quantities are not small in the so-called regime of "eternal inflation"). The problem therefore boils down to calculating statistical moments of the realised number of e-folds in stochastic inflation. In ref. [54], this was done for single-field setups where it was shown that in most potentials, stochastic corrections to inflationary predictions remain small within the observational window as long as the inflationary energy scale is sub-Planckian. In this paper, we extend the analysis to multiple field scenarios. Indeed, stochastic dynamics in multiple field potentials can be highly non-trivial [53,[55][56][57][58], and has a priori the potential to substantially extend the single-field phenomenology.
This work designs a generic stochastic-δN formalism for multi-field inflation and is organised as follows. In section 2, techniques from "first passage time analysis" are employed to calculate the statistical moments of the inflationary numbers of e-folds in multiple field setups. The results obtained in single-field potentials are reviewed, and first indications are given why including more than one scalar field can lead to fundamental differences. In section 3, "harmonic potentials" are introduced, for which the problem can be addressed analytically. In section 4, the calculation is carried out for a subclass of harmonic potentials called "v(r) potentials", where it is notably shown that, in some cases, including more than one scalar field leads to an infinite mean e-folds number. Implications of this "infinite inflation" phenomenon are discussed, in particular studying the probability of probing arbitrarily large-field regions of the potential that are classically impossible to reach otherwise. In section 5, the case where extra fields modulate the surface of end of inflation is investigated.

JCAP06(2016)043
It is shown that depending on the details of the end-surface, stochastic corrections may be large even at sub-Planckian energy densities, and tend to smear out the effect of the modulating fields. Finally, in section 6, we present a few concluding remarks, and we end the paper with several appendices containing various technical aspects.

Moments of inflationary e-folds
We investigate the situation where inflation is driven by D canonical scalar fields φ i , where 1 ≤ i ≤ D, slowly rolling down the potential V (φ i ). Their coarse-grained parts evolve according to the Langevin equations [44] where ξ i are D independent normalised white Gaussian noises, so that ξ i (N 1 )ξ j (N 2 ) = δ i,j δ(N 1 − N 2 ), and V φ i denotes the derivative of the potential V with respect to φ i . Time is labeled by the number of e-folds N ≡ ln a, where a is the scale factor. At leading order in slow roll, the Hubble parameter H is related to the potential energy through the Friedmann where M Pl is the reduced Planck mass. In this section, we explain how the statistical moments of the inflationary e-folds generated by eq. (2.1) can be calculated. We briefly review the results [54] obtained in single-field setups, and provide hints why including more than one scalar field may a priori lead to fundamental differences.

Fokker-Planck equation
Let us first recast these stochastic processes through a Fokker-Planck equation, which governs the time evolution of the probability density P (φ i , N ) that the fields take value φ i at time N . Introducing the dimensionless potential in the Itô interpretation 1 [19,60,61], it reads This equation can be written as ∂P/∂N ≡ L FP · P ≡ −∇ · J , where L FP is the Fokker-Planck differential operator defined in eq. (2.4), ∇ denotes the vector differential operator ∇ i = ∂/∂φ i in field space, which, for simplicity, we assume to be flat, and J is the probability current, In a stationary distribution P stat , by definition, the 1 More generally, the last term in eq. (2.4) can be written in the form with 0 ≤ α ≤ 1, where α = 0 corresponds to the Itô interpretation and α = 1/2 to the Stratonovich one [59]. One can show that keeping terms explicitly dependent on α exceeds the accuracy of the stochastic approach in its leading approximation (2.1). In particular, corrections to the noise term due to self-interactions of smallscale fluctuations (if they exist) are at least of the same order or even larger. This is why explicit dependence on α is dropped here. divergence of J vanishes, corresponding to incompressible flows. If only one field is present, this means that J is uniform in field space, and that, in most interesting situations, the probability current itself vanishes. For example, if field space is unbounded, the normalisation 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, hence everywhere, and this gives If more than one field is present, one can already see that the analysis is much less trivial. The distribution (2.5) is still a solution of the stationarity problem, but non-uniform probability currents are also allowed (since only their divergence must vanish), yielding other solutions. As will be seen, going from one to several fields always renders the problems more difficult, but makes the phenomenology of their solutions richer.

First passage time analysis
For the reasons explained in section 1, we are interested in the duration of the inflationary dynamics described by eq. (2.1), or equivalently, by eq. (2.4). More precisely, the problem we want to solve is the one depicted in figure 1. The inflationary domain in field space is denoted Ω. In the standard situation, it is given by the set of points such that the first slow- However, in order to also consider situations where inflation does not end by slow-roll violation [62][63][64][65][66][67], in the following, Ω is left unspecified. Starting from φ in i at time N in , 2 let N denote the number of e-folds realised before ∂Ω − , the boundary of Ω where inflation stops, JCAP06(2016)043 is crossed. Obviously, N is different for each realisation of the stochastic process under consideration.
In appendix A (see also appendix B), it is shown how the first passage time formalism [68,69] allows one to express the moments of N as the solutions of deterministic differential equations. Introducing the set of functions it is found that, at leading order in slow roll, they satisfy the hierarchy of partial differential see eq. (A.21). This equation is valid for n ≥ 1, where we have defined f 0 = 1. It is the main result of this section and its properties will be discussed in details in what follows. At this stage, let us simply notice that if only one field is present, eq. (2.9) is an ordinary differential equation, the solution of which is given by [54] f whereφ n andφ n are integration constants (see section 2.3). All moments can therefore be calculated analytically, for any potential v(φ). When more than one field is present however, eq. (2.9) is a partial differential equation, for which no generic analytical solution exists. Partial differential equations are notoriously more difficult to study (even at the numerical level) than ordinary differential equations. This is why multiple field scenarios are more challenging from a technical point of view.

Boundary conditions
Single-and multi-field scenarios also differ when considering the boundary conditions according to which eq. (2.9) must be solved. Such conditions need to be specified, for instance to fix the integration constantsφ n andφ n appearing in eq. (2.10). A first requirement is that f n = 0 on ∂Ω − , which simply encodes the fact that, if one starts off the evolution on ∂Ω − , one instantaneously crosses ∂Ω − by definition, so that N = 0 and all its moments vanish. In case Ω is compact in field space (as in hilltop models) so that the inflationary domain is bounded by ∂Ω − in all directions, this defines a Cauchy problem (i.e. there exists a single solution once the boundary condition f n | ∂Ω − = 0 is imposed). In all other situations however, where inflation can proceed at arbitrarily large field value (see figure 1) as in large-field or plateau models, this single absorbing condition on ∂Ω − is not 3 For n ≥ 2, it is sometimes be more convenient to work with the centred moments They satisfy a similar hierarchy of partial differential equations, namely valid for n ≥ 2 if one defines g0 = 1 and g1 = 0.

JCAP06(2016)043
sufficient and one needs to add a second boundary condition. This is why in what follows, a second reflecting (or absorbing) boundary condition is placed at large field value on ∂Ω + . The impact of this extra boundary condition is in fact very much dependent on whether one or several fields are present. Indeed, in appendix C, it is shown that the probability p + (φ i ) to bounce again (or being absorbed into) the extra boundary ∂Ω + before exiting inflation through ∂Ω − , starting from φ in i = φ i , is given by the solution of the ordinary differential equation with boundary conditions p + = 1 on ∂Ω + and p + = 0 on ∂Ω − , see eq. (C.8). If only one scalar field is present, the solution of this equation is given by [54] where one has noted ∂Ω ± = {φ ± }. From this expression, one can see that as soon as v is not vanishing at large-field value [or if it does, as soon as it does not go to 0 faster than 1/ ln(φ)], p + → 0 when φ + → ∞. This means that the system never (i.e. with probability 0) explores large-field regions of the potential. As will be shown in section 4, this implies that, in most single-field setups, the exact location of the second boundary cancels out in all physical results when removed to infinity. In section 4.3 however, it will be shown that this property does not generalise to multi-field scenarios. This is one of the reasons why different results can be obtained in these setups.

Harmonic potentials
In the previous section, the partial differential equations (2.9) governing the moments of the number of inflationary e-folds were derived. For a fully generic multi-field potential, these equations have no analytical solutions and one needs to resort to numerical analysis. Alternatively, in this section we identify a subclass of inflationary potentials for which eq. (2.9) can be solved exactly and the effects associated with the inclusion of multiple fields can be studied analytically.

Polar coordinates
Let us first note that eqs. (2.9) are diffusion equations, akin to the Laplace equation. This suggests that some insight may be gained by reparameterising field space with polar-type coordinates. Since the slow-roll trajectory follows the gradient of the potential at the classical level [i.e. without including the diffusion term in eq. (2.1)], a natural choice 4 is to take the potential v itself for the radial coordinate, completed by D − 1 angular coordinates θ j (with In this expression, recall that the vectorial notation (and the differential operators ∇ and ∆ = |∇ 2 |) refer to field space. For example, where {e φ i } stands for the field space basis. One can always choose the angular variables θ j to form a system of orthogonal variables 6 and one obtains

Harmonic potentials
We now restrict the analysis to potentials for which separable solutions (in the basis {v, θ j }) of eqs. (2.9) exist. An important remark is that purely radial (i.e. independent of θ j ) solutions of eqs. (2.9) can be found if the coefficient in front of ∂/∂v is a function of v only. For this reason, we define "harmonic potentials" as being such that is a function of v only. In order to understand to which extent harmonic potentials allow one to proceed analytically, let us discuss the case where D = 2 fields are present. For two-field potentials, one has a single angular variable θ, and the orthogonality condition ∇(v) ⊥ ∇(θ) mentioned between eqs. (3.1) and (3.2) , where h is an overall factor that is left unspecified at this stage. Let us simply note that, in order for θ to be globally defined [72], the curl of ∇(θ) must vanish, This is the only condition h needs to satisfy, and for harmonic potentials where g depends on v only, it is interesting to notice that it can be fulfilled if h is taken as depending on v only as well, according 5 Strictly speaking, this procedure is well-defined only if the level lines of v(φi) form simply connected hyper surfaces in field space. This is implicitly assumed in what follows, even if more complicated situations can also be studied, either making use of symmetries in the potential function v(φi) as in hybrid inflation, or paving field space with several maps. 6 For example [71], one can start from ∇(v) and use Gram-Schmidt orthogonalisation procedure to iteratively derive ∇(θ1), ∇(θ2), etc. 7 Indeed, in this case, one has

JCAP06(2016)043
One can also check that in this case, ∆θ = h φ 2 v φ 1 − h φ 1 v φ 2 = 0. The differential operator of eq. (2.9) then takes the simple form In this case, up to the overall |∇(v)| 2 factor, all coefficients of the differential operator of eq. (2.9) are explicit functions of the radial coordinate v only, and the problem boils down to solving ordinary differential equations after Fourier transforming the angular coordinate, as will be exemplified in sections 4 and 5. Let us give a few concrete examples of harmonic potentials.

v(r) potentials
A subclass of harmonic potentials is provided by potentials v(r) that depend on The angular coordinates can be chosen to match the ones of the usual spherical coordinates system in D dimension, given in appendix D where their gradients and Laplacians are also derived. This gives rise to where, for simplicity, r has been used as the radial coordinate instead of v. It is interesting to notice that, compared to the single-field case where D = 1, angular terms involving ∂/∂θ j and ∂ 2 /∂θ 2 j are obviously present, but the radial term proportional to ∂/∂r also receives a new contribution. One can also check that when D = 2, eq. (3.4) is recovered. These v(r) potentials are further studied in section 4.

Linear potentials
Another subclass of harmonic potentials is provided by potentials v(u) that depend on a linear combination of the fields only. Here, one can choose the α i constants to be normalised so that α 2 i = 1. In this case, one has |∇(v)| 2 = v 2 (u) and ∆v = v (u), so that g = v (u)/v 2 (u). The "angular" coordinates can then be defined with constant gradients so that {α, ∇(θ j )} form an orthonormal basis of field space (here, the θ j variables are unbounded and should not be viewed as geometrical angles, and "angular" must be understood in a generic way). In this case, one obtains where, for simplicity, u has been used as the radial coordinate instead of v. From this expression, it is clear that the situation is very close to a single-field setup, since inflation is only driven by the "scalar field" u. The only difference with a purely single-field setup arises if the boundary conditions discussed in section 2.3 depend on the other fields, and introduce some "angular" dependence in the solutions of eqs. (2.9). This situation is further investigated in section 5.

Straight potentials
In sections 3.3 and 3.4, we showed that v(r) and linear potentials are "harmonic" in the sense defined in section 3.2. In fact, as we are now going to see, such potentials share the property that the slow-roll classical trajectories are straight lines in field space. We call such potentials "straight potentials". On large scales, entropy perturbations can source adiabatic perturbations only if the background solution follows a curved trajectory in field space [70], which is why these potentials are exactly the ones for which adiabatic perturbations are conserved on large scales, at least at the classical level. We therefore expect them to play a special role in the present analysis. For this reason, in this section, we try to better characterise them.
Since the slow-roll classical trajectories follow the local gradients of the potential, starting from some point φ, the next point on the classical trajectory has coordinates φ + ∇(v), being an infinitesimal number. The gradients evaluated at these two points must be parallel for straight potentials. In other words, the variation in the gradients between φ and φ + ∇(v) must be aligned with the gradient at φ, that is to say By expanding this relation into its components e φ i , it is easy to show that it leads to stands for the i th Poisson bracket between a and b. Straight potentials are therefore such that all Poisson brackets between |∇(v)| and v vanish, meaning that |∇(v)| depends on v only. 8 Since this quantity appears in various places in eq. (3.2), we understand why straight potentials play a special role in the present context. In particular, in section 3.3, it was shown that for v(r) potentials, |∇(v)| 2 = v 2 (r) is a function of v only, which confirms that v(r) potentials are straight potentials. Similarly in section 3.4, it was shown that for linear potentials, is a function of v only, and linear potentials also are straight potentials, as announced above.
Reciprocally, one can show that straight potentials can only be of one of these two types: v(r) potentials or linear potentials. Indeed, let us consider a straight potential v and its (straight) gradient lines in dimension D = 2. We first assume that its gradient lines never intersect in field space. This means that they all are parallel, and one can write footnote 8) and is therefore linear. Let us now assume that there is exactly one intersection point in the gradient lines of v, which, after performing a constant field shift, we set at the origin of field space. It is easy to see that any gradient line not passing through the origin would produce a second intersection point at least, hence all gradient lines go through the origin and one can write ∇(v) = a(φ)e r , where e r is the unit vector pointing to the radial direction r = hence v depends on r only and is of the v(r) type. Finally, let us assume that the gradient lines of v intersect at two or more points. Then, one can convince oneself that an infinite number of other intersection points can be obtained, that fill the entire (or a dense subset of the) field space. Since the gradient of v must vanish when two non-parallel lines intersect (otherwise its direction would be ill-defined), this means that v is constant, and this case is in fact trivial. This result can be generalised to D > 2 where one finds that the potential is of the v(r) type within the field subspace that is orthogonal to the one containing the fields of which v is independent.
The situation is schematically summarised in figure 2. Straight potentials are a specific class of harmonic potentials. They are either linear or v(r) potentials, and singlefield potentials lie at the intersection between these two. Let us finally notice that not all harmonic potentials are straight. For example, let us consider a "loop corrected" potential of the form The function g defined in eq. (3.3) is constant, g = −1/(v 0 α), and such potentials are therefore harmonic. However, one has which is not a vanishing function, hence loop corrected potentials are not straight. More generally, this is the case for all potentials v(w) that are functions of w = i φ i only, for which g = v (w)/v 2 (w).

v(r) potentials and infinite inflation
In section 3.3, it was shown that v(r) potentials provide a subclass of harmonic potentials, for which analytical solutions of the diffusion equations (2.9) can be found. In this section, {v, θ1, · · · θD−1}, the Poisson bracket {a, b} v,θ j is given by Therefore, if all Poisson brackets between a and b vanish, then {a, b} v,θ j vanishes as well. If one takes b = v, this means that ∂a/∂θj = 0, for all 1 ≤ j ≤ D − 1, hence a depends on v only, as is the case for |∇(v)| in eq. (3.10). we derive such solutions and use these potentials to illustrate the physical implications of including more than one scalar field in the analysis. If one sets boundary conditions to be angular independent, that is to say if one assumes that inflation ends at r = r − and that a reflecting wall is placed at r = r + , angular independent solutions of eq. (2.9) can be obtained. More precisely, combining eqs. (2.9) and (3.6), one obtains where we recall that f 0 = 1. One can check that f n (r − ) = 0 (absorbing boundary condition) and that f (r + ) = 0 (reflecting boundary condition). If one wanted to place an absorbing boundary condition at r + instead, one would have to change the upper bound of the second integral to a smaller value [54], but this would not affect the following considerations. In the same manner, if, instead of the situation depicted in figure 1 where the fields classically decrease during inflation, one considered a hilltop potential symmetric about r = 0, one would have to replace r + by 0 in the above formula, and this case is also discussed in what follows. A preliminary remark is that both the number of fields D and the location of the reflecting boundary condition r + explicitly appear in eq. (4.1), and are therefore expected to play a role. For illustration, in the left panel of figure 3, the mean number of e-folds (4.1), N = f 1 (r in ), is displayed as a function of the initial condition v in for a quartic potential v ∝ r 4 , for different values of r + , and for D = 1 (solid lines) and D = 7 (dashed lines). One can check that, in some regimes at least, the result strongly depends on r + and D indeed (see also the right panel of figure 3), in a way that we now analyse in more details.   1. In the opposite regime, N deviates from N cl in a way that depends on v + and D, and that is further discussed in the main text. One should note that, in principle, v in > 1 lies outside the validity range of the present calculation since it corresponds to initial super-Planckian energy density, but it is displayed to make clear the asymptotic behaviour of the mean number of e-folds at large initial field value. In the right panel, v in = 1 is fixed, but v + varies, and different values of D are displayed. One can check that when D < p, a finite asymptotic value is reached when v + → ∞, while when D ≥ p, N diverges in this limit.

Classical limit
A first important consistency check consists of verifying that the correct classical limit is recovered. In the classical picture (i.e. without the stochastic terms), the slow-roll equation of motion (2.1) for the fields is given by dφ i /dN = −M 2 Pl v φ i /(vr). This gives rise to dr/dN = −M 2 Pl v /v, and the classical number of e-folds reads Let us see how this formula can be recovered from eq. (4.1). In the classical limit, energy densities are sub-Planckian v 1 and the integral over r in eq. (4.1) is dominated by its contribution close to the lower bound r − , near which one can Taylor expand 1/v at first order, 1/v(r ) 1/v(r ) − v (r )/v 2 (r )(r − r ). This gives rise to  Table 1. Mean number of e-folds realised in v(r) potentials and probability of exploring arbitrarily large-field regions of the potential, when D fields are present.
By plugging this formula into eq. (4.1), one obtains in the classical limit. In the left panel of figure 3, the classical formula (4.2) is displayed and one can check that, when v in 1, it provides a good approximation to the full stochastic results indeed.
The fact that the classical trajectory arises as a saddle-point limit of the full quantum dynamics is not so surprising, since it is a common feature of path integral calculations in quantum field theory. Let us also notice that the classical limit (4.2) depends neither on the number of fields D nor on the location of the upper boundary condition r + , and matches the single-field result. However, as we are now going to see, stochastic corrections break this classical D and r + invariance and introduce dependence on both the number of fields and the location of the upper boundary condition.

Infinite inflation
The validity of the classical limit relies on the assumption that the second integral in eq. (4.1) is dominated by its contribution close to the lower bound r . If this is correct, this means that the upper bound, r + , can be removed to infinity without affecting the leading order result, providing a well-defined regularisation procedure. In the left panel of figure 3, one can see that for D = 1, the curves saturate to an asymptotic behaviour when v + increases, and such a procedure seems therefore to exist. However, for D = 7, the result does not seem to converge when v + increases. This is why in the right panel of figure 3, the mean number of e-folds is displayed as a function of v + for quartic v(r)-potentials, for a fixed v in = 1 and a few values of D. One can see that when D < 4, N goes to a constant value when v + → ∞, while when D ≥ 4, it diverges. This confirms that the number of fields plays an important role in determining whether the limit r + → ∞ is finite or not.
More precisely, the mean number of e-folds is finite if the integrand of eq. (4.1) is integrable when r → ∞, that is to say if r D−1 /v(r) is an integrable function. This criterion depends on the number of fields D, as already noticed, but also on the large-field behaviour of the potential. Let us distinguish the three following cases: • If, at large-field value, the potential is of the "Plateau" type and v goes to a constant value v ∞ > 0, then r D−1 /v(r) is never integrable and an infinite mean number of e-folds is always realised, regardless of the number of fields.

JCAP06(2016)043
• If, at large-field value, the potential is of the monomial type v ∝ r p , then r D−1 /v(r) is integrable only when D < p, and an infinite mean number of e-folds is realised as soon as more than p fields are present. This is consistent with the previous discussion about the right panel of figure 3.
• If the potential is of the "hilltop" type and symmetric around 0, as explained above, r + has to be replaced by 0 in eq. (4.1). In this case, the integrability of r D−1 /v(r) needs to be checked around 0 instead of infinity. If v is finite at r = 0 this is always the case, hence the mean number of e-folds is never infinite in such potentials.
The situation is summarised in table 1. In a large class of potentials (plateau potentials and some monomial potentials), the mean number of e-folds is infinite, and we call this phenomenon "infinite inflation". Let us notice that this is different from "eternal inflation" [73][74][75][76] where volume weighting is included and the diverging quantity is the physical volume of the inflating part of the Universe, not N . Infinite inflation implies eternal inflation but is a stronger statement. For example, eternal inflation can be realised in hilltop models [74,77] while, as we have just shown, infinite inflation never occurs in such potentials. Another important remark is that for monomial potentials, whether infinite inflation occurs crucially depends on the number of fields, which therefore plays the role of an "order parameter" (as illustrated below in figure 4). The number of dimensions is a critical parameter for many stochastic processes (for instance in recurrence problems [78]) and this may therefore not be so surprising. The key feature is that the more fields, the larger the volume in field space to realise inflation and the more common infinite inflation.
Beyond the physical implications related to the possibility of realising arbitrarily large number of e-folds [79][80][81], infinite inflation raises the issue of practical calculability of observables. Indeed, since the correlation functions of scalar adiabatic fluctuations are related to the moments of the number of e-folds, see section 1, it is not clear what the predictions for these observables are when those moments are infinite. How these infinities regularise or not is in fact a non-trivial question that we investigate separately in ref. [82]. At this stage however, let us notice that infinite inflation may suggest that the system explores regions of the potential that are far away from what its classical trajectory would allow it to reach, and that observables may be sensitive to the physics at play in these remote regions. For this reason, we now study how likely it is to explore large-field regimes in stochastic multiple field inflation.

Large field exploration
In this section, we derive the probability p + (r) that, starting from r, the system bounces at least once against the reflecting wall located at r + before exiting inflation at r − (or alternatively, if an absorbing wall is located at r + , the probability that the system exits inflation at r + rather than r − ). In section 2.3, it was explained that p + is given by the solution of eq. (2.11) with boundary conditions p + (r − ) = 0 and p + (r + ) = 1. Making use of eq. (3.6), one obtains When the upper boundary condition r + is removed to infinity, one obtains a non-vanishing probability p + if the function r 1−D is integrable (assuming that v has a positive limit at   infinity). Contrary to the case of infinite inflation in section 4.2, this condition is independent of the shape of the potential at large-field value, and p + > 0 as soon as strictly more than 2 fields are present. This information is added in table 1, and in figure 4, eq. (4.6) is displayed for monomial potentials v ∝ r p (left panel, for different values of p and v in ) and a plateau potential, the Starobinsky model [1], v = v ∞ [1 − exp(− 2/3r/M Pl )] 2 (right panel, for different values of v ∞ ), when r + is removed to infinity. One can check that when D ≤ 2, p + = 0. When D > 2, strictly speaking, p + > 0, but one can see that p + is non-negligible only when D is larger than some value that depends on the parameters of the potential and on the initial field vev. Schematically, this value is realised when the integrand of the integrals in eq. (4.6) is maximal at v in . In monomial potentials, this leads to the conclusion that p + is non-negligible when while in plateau potentials, one obtains the condition One can numerically check that, indeed, these expressions provide good estimates of the point where p + starts to be non-negligible. This shows that including more fields increases the -15 -

JCAP06(2016)043
probability to explore large-field regions of the potential, but for sub-Planckian energy scales, one needs a very large number of fields to obtain a substantial probability. For example, if one normalises the overall mass scale of the potentials to fit the measured amplitude of the scalar power spectrum [83] and starts the evolution 50 (classical) e-folds before the end of inflation, one finds v in ∼ 10 −11 p for monomial models and v ∞ ∼ 10 −12 for the Starobinsky potential, so that 10 11 fields would be required to obtain appreciable values of p + , a very large number indeed. Of course, from a model building perspective, the shape of the potential may be very different at very large field outside the observational window than what cosmological observations constrain at smaller field values (for example [84,85], the potential may be of the Plateau type where the scales probed in the CMB cross the Hubble radius, but of the monomial type at larger field), and if inflation starts high enough in the potential, large-field exploration, enhanced by the presence of multiple fields, may become likely. But the above results suggest that, in the simplest setups, cosmological observations at small (i.e. sub-Planckian) energies carry limited information about the physics taking place at much higher energy (at least through stochastic effects). We further investigate this question in ref. [82].

Inhomogeneous end of inflation
In section 4, we have considered the case of v(r) potentials where the dynamics is governed by the "radial" field v only, and angular independent solutions can be found. In this section, we study situations where both v and θ j play a role, and study the simple setup where the inflationary dynamics is effectively driven by a single field φ while extra fields χ j (to be identified with the "angles" θ j ) only appear at the surface defining the end of inflation. This notably allows one to describe "inhomogeneous end of inflation" [86] where inhomogeneities induced from the additional light fields on the end-surface make additional contributions to curvature perturbations on super-Hubble scales.
In the terminology introduced in section 3, this case is called "linear potential" (i.e. v depends on a linear combination of the fields only) and is described in section 3.4. In eq. (3.8), one can see that the "angular" sector is only affected by a pure diffusion term, which suggests that some insight may be gained by Fourier transforming the functions Plugging this ansatz into eq. (3.8), eq. (2.9) gives rise to the set of recursive ordinary differential equations where k 2 ≡ |k| 2 . If the end-surface ∂Ω − is parametrised by the function φ − (χ), and if the inflationary domain is limited from above at the reflecting surface ∂Ω + defined by φ = φ + , these equations need to be solved with the boundary conditions 3) The procedure one needs to carry out is therefore the following: solve eq. to set all integration constants simultaneously. In practice, such a calculation may need to partly rely on numerical analysis, but it is still more straightforward (and numerically less expensive) than having to solve the full partial differential equations (2.9).

The example of exponential potentials
In order to illustrate how the above procedure works in practice, let us consider the case of "power-law inflation" [87] where the potential is of the exponential type v ∝ e αφ/M Pl . To be explicit, we study the situation where one extra field χ modulates the end-surface through where µ is a modulation parameter (when µ → 0, one recovers the standard single-field setup), and χ 0 is the scale over which the modulation takes place. In this case, solutions of eq. (5.2) that are 2πχ 0 -periodic in χ can be found, and one can replace the continuous Fourier transform of eq. (5.1) by a discrete Fourier sum over integer numbers k, f n (φ, χ) = k e −ikχ/χ 0 f k n (φ). Let us also note that since the exponential potential is conformally invariant, shift symmetry in the inflaton field value allows us to take φ end = 0 without loss of generality, and write v = v end e αφ/M Pl . For the mean number of e-folds f 1 , recalling that f 0 = 1, eq. (5.2) then gives rise to .

(5.5)
When k = 0, the solution can be expressed in terms of the exponential integral function Ei, while when k = 0, the solution is given in terms of confluent hypergeometric functions 1 F 1 .
Requiring that f 1 is real, one obtains where A 0 , A k , B 0 , B k are integration constants, that must be fixed making use of eqs. (5.3). At this stage, one has to proceed numerically. In practice, if the summation over k in eq. (5.6) is truncated at order k max , one has 2(k max + 1) integration constants to fix. One can choose k max + 1 values of χ uniformly distributed in [0, πχ 0 ], and evaluate both parts of eqs. (5.3) at these values. This gives rise to 2(k max + 1) linear equations for the 2(k max + 1) integration constants, that one can solve with standard matrix inversion methods. One then increases k max until a sufficient level of convergence is reached.
The result of such a procedure is displayed in the left panel of figure 5, where the mean number of e-folds is plotted as a function of the initial values of φ and χ, with α = 0.1, χ 0 /M Pl = 1, µ/M Pl = 0.5 and v end = 0.05 (this last value does not lead to the right scalar power spectrum amplitude [83], but it is used to make clearer the effects we want to comment on). The value of φ + has been taken to be sufficiently large so that it does not play any role, which is possible since the model is effectively single-field and non-plateau during inflation, as follows from the discussion in section 4. One has taken k max = 100, but the result is already very well converged when k max 8. The black dashed lines are various level lines of N and have been superimposed to guide the eye. The white region at the bottom left corresponds to φ < φ − (χ), which lies outside of the inflationary domain. One may also note that only the region 0 ≤ χ ≤ πχ 0 is displayed, since the result for other values of χ can easily be inferred using the symmetry and periodicity properties of f 1 .

Smearing out the modulating field
In the left panel of figure 5, one can notice that, going from the left to the right, the level lines of N at first follow the modulation of the end-surface, and then become more and more straight. This result can be understood in terms of the two limits v 1 and v 1. In the classical limit v 1, the diffusion term acting on the extra field χ in eq. (3.8) is negligible, and χ freezes to its initial value. In this limit, the point where the system exits inflation in field space becomes deterministic and is simply given by φ − (χ in ). A similar expression to eq. (4.2) can therefore be obtained, except that the lower bound now explicitly depends on χ, For the power-law model under consideration, this gives rise to N cl = [φ−µ cos(χ/χ 0 )]/(M Pl α). This quantity is displayed in the right panel of figure 5 where one can check that, at small -18 -

JCAP06(2016)043
field where v is not too large, it provides a good approximation to the full stochastic result given in the left panel.
In the stochastic dominated limit where v 1, the diffusion term acting on χ becomes large, and quickly randomises the vev of this extra field. In this regime, memory of the initial conditions on χ is quickly erased, which explains why the result becomes dependent on φ only and the level lines of N in the left panel of figure 5 tend to being merely vertical. Technically, one can check that the first line in eq. (5.6), which is the 0 th (i.e. χ-independent) mode, provides the dominant contribution in the limit v 1. This term exactly matches the one obtained in eq. (2.10) in a purely single-field setup.
Therefore, stochastic effects tend to erase the presence of modulating fields by blurring their initial values and averaging the exit point over the end-surface. In practice, the size of the effect depends on how the scale over which modulation takes place (denoted χ 0 in the present model) compares to the dispersion acquired by the modulating field at the end of inflation, but not so much on the size of the modulation itself (here denoted µ). In the simple toy model discussed in this section, if one sets v 10 −10 , corresponding to the value that would fit the measured scalar power spectrum amplitude, one finds 9 that the effect is large if χ 0 /M Pl 10 −4 . The situation is therefore different from the purely single-field case [54] where, in the simplest setups, the stochastic effects must be subdominant since v is observationally constrained to be small. Here, large stochastic corrections can be obtained even if v 1, depending on the scales of the features in the end-surface.

Conclusion
Let us now summarise our main results. In this paper, we have investigated the inclusion of multiple fields in the stochastic inflationary framework. In particular, we have derived a hierarchy of diffusion equations that can be used to calculate the moments of the numbers of e-folds. This sets up the formalism that will be required to compute correlation functions of primordial curvature perturbations in the stochastic δN approach in multi-field inflation, which we will investigate in ref. [82]. Compared to the single-field case, two main differences have been noted. Firstly, the volume available in field space is strongly dependent on the number of fields D. For instance, for a given mass scale µ, the volume associated to the region |φ| < µ is given by π D/2 µ D /Γ(D/2 + 1). This is why the dimension often plays the role of a critical parameter, or an order parameter, in stochastic processes. This is for example the case in the well-known recurrence problems where it can be shown that, for a random walk on a D-dimensional lattice, the probability to return to the origin is 1 for D = 1 or D = 2, and strictly smaller than 1 for D > 2 (and is then given by the so-called Pólya's random walk constants [78]). Thus D = 2 can be seen as a critical dimension for the recurrence problems. A similar effect takes place in stochastic inflation, where we have shown that the probability to explore arbitrarily large field regions of the potential vanishes for D ≤ 2 and is strictly non-zero when D > 2, regardless of the potential. In practice, depending on the large-field profile of the potential, this probability is often much smaller than one (although not strictly vanishing) for initial sub-Planckian field values, even when a large number of 9 Here, the dispersion acquired by the freely diffusing χ field is given by √ 50H/2π √ 100vMPl, where 50 is taken to be the number of e-folds between Hubble exit time of the scales probed in the CMB and the end of inflation, and H/2π is the amplitude of the noise term in eq. (2.1) which we assume to be roughly constant over these last 50 e-folds of inflation.

JCAP06(2016)043
fields is present. However, this illustrates why the number of fields plays a critical role in stochastic inflation and why including more fields tends to make the system explore regions of the potential that would otherwise be inaccessible with a single field. This notably translates into "infinite inflation", i.e. the fact that the mean number of realised e-folds is infinite in all plateau models, and in monomial potentials v ∝ φ p if p or more fields are present. Infinite inflation has important consequences for problems related to the total duration of inflation (which can be a crucial parameter in deriving the distribution of curvaton-type fields [81,88] or other light spectator fields [80] at the end of inflation, but also to determine whether pre-inflationary imprints on large scales can arise in "just enough inflation" types of scenarios [79]). Since the correlation functions of cosmological perturbations are related to the moments of e-folds number, this also means that infinite quantities appear in cosmological observables. Whether these infinities can be regularised and whether one needs to modify the large-field sector of the theory to make the results finite is therefore an important, number of fields dependent, question, that we address separately in a companion paper [82].
Secondly, including more than one scalar field offers more directions (than the classical one, aligned with the gradient of the potential in field space) along which the system can fluctuate. If quantum diffusion spreads the fields distribution along these extra directions on a distance (in field space) over which the potential does not vary much, the associated effect remains small. However, in the opposite case, large stochastic corrections can be produced. This situation was exemplified in section 5 in the case where inflation ends inhomogeneously, i.e. for a value of the inflaton field φ that depends on a modulating field χ. If the dispersion in the χ field direction at the end of inflation is larger than the χ vev scale over which the surface of end of inflation is modulated, large stochastic effects are obtained. Let us stress that this can happen even at sub-Planckian energies, contrary to what is found in purely single-field setups. In ref. [54] indeed, it was shown that, even if stochastic effects can shift the location of the observational window along the inflationary potential, once the observational window is fixed, single-field stochastic corrections scale with v, which is aways a very small quantity. In multiple field scenarios, we have found that this property does not hold, and that stochastic corrections to cosmological observables can a priori be large within the observational window [82]. This opens interesting possibilities for probing quantum effects on inflationary dynamics.

JCAP06(2016)043
equation (2.4), valid for any initial condition, gives, in this specific case, We now want to write a similar equation but where the time derivative acts on the first time argument, N in . Let us start from the Chapman-Kolmogorov relation which simply states that any process starting at φ in i at N in and ending up at φ i at N goes through some pointφ i at some intermediate timeN , and we integrate over all intermediary pointsφ i . When one differentiates this relation with respect toN , the left hand side vanishes, and one obtains where, in the second line, we have used the Fokker-Planck equation (A.1). The second term in the integral of eq. (A.4) can be integrated by parts making use of the adjoint Fokker-Planck operator L † FP defined as and one obtains the adjoint Fokker-Planck equation The next step is to introduce the survival probability S(N ) of having not yet crossed ∂Ω at time N , This corresponds to the probability of having N > N . If P (N ) denotes the probability distribution associated with N , this means that By differentiating this expression with respect to N , one obtains

JCAP06(2016)043
The n th moment of N can therefore be expressed as where integration by parts has been performed, with the requirement that ρ vanishes on ∂Ω (absorbing boundary conditions). By applying the adjoint Fokker-Planck operator L † FP (φ in i ) on this relation and making use of the adjoint Fokker-Planck equation (A.6), one obtains At this stage, let us recall that one is dealing with a Markovian process, for which the transition probability depends on N − N in only. One then has When n = 1, this gives rise to Here, we have used the fact that, in the infinite future, all realisations have crossed ∂Ω, 10 so that ρ(φ i , ∞|φ in i , N in ) = 0. When n ≥ 2, from eq. (A.10), one notices that N n−1 appears in the right hand side of eq. (A.15), giving rise to 20) and recalling that the adjoint Fokker-Planck operator is given by eq. (A.5), the previous analysis shows that the f n functions satisfy the hierarchy of partial differential equations for n ≥ 1, where f 0 = 1.

JCAP06(2016)043 B First passage time from Langevin equation
The same results as the ones derived in appendix A can also be obtained starting from the Langevin equation (2.1). In this section, we quickly sketch such a derivation for illustrative purpose, and to make comparison with ref. [54] easier. Let f (φ) be a generic function of the fields value φ = (φ 1 , φ 2 , · · · , φ D ). If φ is a realisation of the stochastic process (2.1), its variation is given by Integrating this relation between N = 0 where φ = φ in and N = N where φ = φ end ∈ ∂Ω − , one obtains the Itô's lemma [89] f Let us now apply this lemma to f 1 , defined in section 2 as being the solution of eq. (2.9) which vanishes along ∂Ω − . By definition, the first term in the left hand side of eq. (B.3) vanishes, and the integrand of the second integral of the right hand side is −1. This gives rise to By taking the stochastic average of this equation, one is led to Note that the fact that the stochastic average of the integral term in eq. (B.4) vanishes is not trivial a priori since not only the integrand but the upper bound of the integral itself is stochastic, but because the noises ξ i are uncorrelated at different times, this can be shown rigorously [90]. This demonstrates eq. (A.20) for n = 1. Larger values of n can be dealt with in a similar manner. Indeed, by squaring eq. (B.4) and taking the stochastic average of it, one obtains Let us now apply Itô's lemma (B.3) to g 2 ≡ f 2 − f 2 1 , where f 2 is defined in section 2 as being the solution of eq. (2.9) which vanishes along ∂Ω − . By definition, the first term in the left hand side of eq. (B.3) vanishes, and the integrand of the second integral of the right hand side is given by −2 i f 2 1,φ i v. One obtains where eq. (B.6) has been used in the last equality. Since g 2 = f 2 − f 2 1 , this gives rise to f 2 (φ in ) = N 2 , i.e. eq. (A.20) for n = 2. Applying the same method, one can iteratively proceed and extend the result to any n. In this section, we consider the case where ∂Ω is made of two (or more) disconnected pieces (say ∂Ω − and ∂Ω + ) and one wants to determine with which probability p + the system exits Ω crossing ∂Ω + (or respectively, with which probability p − = 1 − p + the system exits Ω crossing ∂Ω − ).
The same techniques as the ones employed in appendix A can be employed to determine this quantity as a function of the initial conditions φ in where the system starts off its evolution. By definition, p + corresponds to the probability that the system is somewhere along ∂Ω + at time N , where N is integrated over all possible values between N in and ∞, Let us now apply the adjoint Fokker-Planck operator defined in eq. (A.5) to this relation. Making use of eq. (A.6), one obtains In eq. (C.5), similarly to what was performed in eq. (A.13), one has used the fact that the stochastic process under consideration is Markovian, hence the transition probability depends on N − N in only. To obtain the final result (C.7), one has also used the fact that, as mentioned below eq. (A.18), all realisations have crossed ∂Ω in the infinite future hence ρ(φ i , ∞|φ in i , N in = 0), together with the initial condition ρ φ i , N in |φ in i , N in = δ D (φ i − φ in i ) (and the assumption that φ in i ∈ ∂Ω + , otherwise we already know that p + = 1 by definition). The probability p + (φ i ) that the system first reaches ∂Ω + starting from φ in i = φ i is therefore given by the solution of the ordinary differential equation with boundary conditions p + = 1 on ∂Ω + and p + = 0 on ∂Ω − . In these expressions, we have defined From here, the derivatives of the angular coordinates can be calculated, and one obtains The norm of the gradient of the angular coordinates appearing in eq. (3.2) is therefore given by (D. 13) In particular, |∇(θ 1 )| depends on r only. In the same manner, the second derivatives of the angular coordinates can be calculated, from which their Laplacian, appearing in eq. (3.2), are found to be sin −2 (θ ) . (D.14) In particular, this implies that ∆θ D−1 = 0 and one recovers the fact that, when D = 2, the only angular coordinate has vanishing Laplacian.