The $1/N$ expansion for stochastic fields in de Sitter spacetime

We propose a $1/N$ expansion of Starobinsky and Yokoyama's effective stochastic approach for light quantum fields on superhorizon scales in de Sitter spacetime. We explicitly compute the spectrum and the eigenfunctions of the Fokker-Planck operator for a O($N$)-symmetric theory with quartic selfinteraction at leading and next-to-leading orders in this expansion. We obtain simple analytical expressions valid in various nonperturbative regimes in terms of the interaction coupling constant.

Focusing on the case of spectator fields with a standard kinetic term 1 in pure de Sitter spacetime, the relevant stochastic dynamics possesses a late time equilibrium state described by a stationary probability distribution whose form is known for an arbitrary potential. This allows one to compute a variety of one point correlators often with analytic control. Higher order correlators exhibit nontrivial spacetime dependences which can be conveniently expressed in terms of a spectral decomposition involving the eigenvalues and eigenfunctions of the associated Fokker-Planck operator. If those can be computed numerically [1,2,26,27], it is often of interest to also have some analytic control, for instance, as checks of numerical results, or for comparison with direct QFT calculations [4,6,17,19,24].
To date, only few explicit analytical results are known concerning the spectrum of the Fokker-Planck operator, even in the case of a simple quartic potential. Of course, when the relevant coupling constant (see below) is small, a systematic perturbative treatment of the eigenvalue problem is feasible. This has been implemented at the first nontrivial orders both in the case of a positive [2] and of a negative [26] square mass for a single scalar field theory. Such perturbative results, however, are not valid in the (phenomenologically relevant) cases of essentially massless fields. Nonperturbative expressions of the three lowest eigenvalues have been obtained from the calculation of various correlators in a 1/N expansion for a O(N )-symmetric theory [3].
In the present work, we setup a proper 1/N expansion directly at the level of the Fokker-Planck eigenvalue equation for systems with O(N ) symmetry. In the case of a quartic potential, we obtain simple analytical expressions of all eigenvalues and eigenfunctions both at leading and next-to-leading orders, which reproduce and encompass the results mentioned above. These provide benchmark results, valid for arbitrary value of the coupling (within the validity of the stochastic approach, i.e., for light fields), for various quantities of physical interest, such as correlation lengths and times, relaxation and decoherence timescales, or various spectral indices, relevant for phenomenological applications [2,[27][28][29].
In Sec. II, we briefly review the stochastic approach for the O(N ) theory and its formulation in terms of an eigenvalue problem for the associated Fokker-Planck operator. We setup the 1/N expansion of the problem and present the solution to the eigenvalue problem at leading and next-to-leading orders in Sec. III. We discuss our findings together with their physical interpretation and the comparison with previously existing results in Sec. IV. Sec. V summarizes our conclusions. The details of the next-toleading-order calculation are given in the Appendix A and we present some comparison with numerical results in Appendix B.

II. STOCHASTIC FORMALISM AND FOKKER PLANCK EQUATION
We consider a N -component scalar fieldφ a embedded in the expanding Poincaré patch of de Sitter spacetime in D = d + 1 dimensions, with metric ds 2 = g µν dx µ dx ν = −dt 2 + e 2Ht d x 2 , in terms of the cosmological time t and the comoving spatial coordinates x. We set the Hubble scale H = 1 in the following. With standard notations, the microscopic action reads For light fields in units of H, 2 the quantum fluctuations on superhorizon scales can be described as those of an effective stochastic variable driven by the subhorizon degrees of freedom. On such scales, spatial gradient are negligible and one can treat the problem as a collection of independent Hubble patches described by an appropriate Langevin equation [1]. Absorbing unimportant numerical factors through the redefinitions ϕ a =φ a dΩ D+1 /2 and V =V Ω D+1 /2, with Ω n = 2π n/2 /Γ(n/2), the latter reads where ∂ a = ∂/∂ϕ a . Here, the field ϕ a denotes a spatially averaged field over a Hubble patch and the noise ξ a reflects the effect of the subhorizon (quantum) fluctuations, which constantly feed this coarse-grained degree of freedom as a result of the gravitational redshift. We refer the reader to the literature [1,5] for details on the derivation of Eq. (2). Treating the subhorizon sector in the linear approximation and assuming the Bunch-Davies vacuum yields a white Gaussian noise, entirely characterized by the two-point correlator Following standard procedures [30], this stochastic dynamics is equivalently formulated in terms of the following Fokker-Planck equation for the field probability distribution function (PDF) P ≡ P (ϕ a , t) which can, itself, be reduced to an eigenvalue problem, as we now recall for the case of a potential with O(N ) symmetry. First, introduce the reduced PDF p(ϕ a , t), defined as P (ϕ a , t) = e −V (ϕa) p(ϕ a , t), in terms of which Eq. (4) takes the form of the Schrödinger-like equation where ∆ ϕ = ∂ a ∂ a and For a O(N )-symmetric potential, it is convenient to use spherical coordinates in field space and to decompose the angular dependence onto generalized spherical harmonics Y i (θ i ), where θ i=1,...,N −1 denote the N − 1 angular variables in field space and the i are integers such that | 1 | ≤ 2 ≤ . . . ≤ N −1 . These harmonics diagonalise the angular part of the operator ∆ ϕ = where we have noted = N −1 and ρ = ϕ 2 . For the purpose of the 1/N expansion below, it is convenient to introduce the scaled radial variable and potentials x = ρ/ √ N , v = V /N , and w = W/N . We have where the prime denotes a derivative with respect to x. Seeking solutions to Eq.
To set up an appropriate 1/N expansion we first need to properly control the limit N → ∞ of the theory. This requires one to understand how the various quantities in Eq. (9) scale with N at large N . To this aim, it is instructive to consider the exactly solvable case of a purely quadratic potential v(x) = m 2 x 2 /2, or, equivalently, In that case, Eq. (9) is nothing but the radial Schrödinger equation for a symmetric N -dimensional harmonic oscillator with unit mass and pulsation ω = m 2 and whose energy levels are shifted by −m 2 /2. The spectrum is degenerate in the "quantum numbers" i and labeled by a nonnegative integer n, and the eigenfunctions are easily obtained in Cartesian coordinates (in field space) as products of Hermite polynomials. In radial coordinates, they can be written as finite polynomials, where n − = 2k is bound to be a nonnegative even integer and is a finite polynomial in x whose coefficients a q are determined by the recursion relation The latter and, hence, the polynomial in Eq. (13) has a well-defined limit when N → ∞ at fixed n and . In this limit, Eq. (14) becomes which is solved as a q = a 0 C q k (−2m 2 ) q , with C q k the binomial coefficient, yielding the leading-order radial eigenfunctions, up to a normalization constant, A few comments are in order here. First, the eigenvalues (11) do not scale with N . Second, the appropriate radial variable to work with in order to obtain a nontrivial large-N limit is the scaled variable x. Finally, taking the limit N → ∞ directly at the level of the eigenvalue equation (9) yields the result R = 0 which, although consistent with the naive large-N limit of Eq. (12) at fixed x, is clearly too harsh. To avoid this caveat, it thus appears important to factor out the exponential factor in (12). We now apply these lessons to the case of a more general potential.

B. Interacting case
Following the previous discussion, we introduce the reduced radial function R = e −N v r, with v the relevant potential. The eigenvalue equation (9) becomes which possesses a well-defined large-N limit. Setting N → ∞, we get the following first order equation This can be easily integrated for polynomial potentials in terms of the roots of (1 − 2xv ). We will focus on the case of a quartic potential which provides simple analytical formulas. Using the identity where the right-hand side of Eq. (18) can be decomposed in simple fractions as with Integrating Eq. (22) is now elementary and yields the leading-order radial function The obtained eigenfunctions lead to normalizable PDFs thanks to the exponential factors we have extracted, P ∝ e −V R ∝ e −2V r. Requiring the solutions to be regular for all x selects a discrete subset, as expected from the analogous quantum mechanical problem. Using the fact that m 2 ± ≥ 0, we see that regularity imposes α + = k ∈ N. In turns, this implies that the eigenvalues are indexed by to nonnegative integers n and such that n − = 2k (so that ≤ n) and are given by The corresponding eigenfunctions thus read Notice that, as expected, the lowest eigenstate of the system has a vanishing eigenvalue, Λ 0,0 = 0. This is, in fact, guaranteed by the symmetries of the Fokker-Planck operator and simply corresponds to existence of a latetime, equilibrium state of the system, with PDF P ∝ e −2V . Equations (25) and (26) completely solve the eigenvalue problem in the large-N limit and provide the leading order of a systematic 1/N expansion. As an illustration, we explicitly compute the next-to-leading order corrections in Appendix A. We report here the results for the eigenvalues where a n, = n(3n − 2) − ( − 2) and b n, = a −n, .

IV. DISCUSSION
As already mentioned, all physical information in the stochastic approach can be expressed in terms of the eigenvalues and eigenfunctions. The results of the previous section provide analytical expressions, nonperturbative in the coupling constant, which allow one to discuss The leading-order eigenvalues Λ n, as functions of m 2 at fixed coupling λ = 1. We show three groups corresponding to (from bottom to top on the axis m 2 = 0) n = , n = + 2, and n = + 4 with = 0, 1, 2, 3 for each group.
various physically relevant regimes. Before doing so, let us quickly mention that the general results of the previous section reproduce the findings of Refs. [3,19], where the eigenvalues Λ 1,1 , Λ 2,0 , and Λ 3,1 had been obtained by other means at leading and next-to-leading-for Λ 1,1orders. Also worth mentioning is the fact that Eqs. (25) and (26) trivially reduce to the Gaussian results (11) and (16) when λ = 0, which corresponds to m 2 − = 0 and m 2 + = m 2 . We show, in Fig. 1, the evolution of the spectrum of the theory as a function of the parameters of the potential. The Gaussian limit is controlled by the dimensionless coupling λ/m 4 . The lifting of the Gaussian degeneracy, where the eigenvalues Λ n, are independent of , is given by, in the large-N limit, For m 2 → 0, the perturbative treatment is invalid. However, the nonperturbative expressions of the previous section remain valid and give The square mass m 2 dyn = λ/2 is dynamically generated by the self-interactions and is of gravitational origin. It is the so-called dynamical mass, which quantifies the local field fluctuation ϕ 2 = 1/(2m 2 dyn ). We see that the spectrum (29) consists in multiples of this square mass and is thus analog to that of a Gaussian potential with pulsation m 2 dyn (although the degeneracies of the eigenvalues are different from the Gaussian case).
Finally, another interesting regime is that of a doublewell potential with m 2 < 0, which is also strongly nonperturbative due to the flat (Goldstone) directions in the potential. We find that the deep-well limit λ/m 4 1 mirrors the near Gaussian case in that m 2 + → 0 and  FIG. 2. The leading-order (lines) and next-to-leading-order (dashed) eigenfunctions R n, (x) for some of the lowest levels in the case m 2 = 0 and λ = 1. We take N = 2 in this figure to amplify the difference between the two curves in each case. Here, the normalization are chosen such that either the function or its first nonzero derivative at x = 0 is fixed to 1. In practice, this means that a0 = 1 and the function C(x) in Appendix A is chosen such that C(0) = 0. m 2 − → |m 2 |. In this regime, The eigenvalues now only depend on the even integer n − and are thus multiples of 2|m 2 |. This is a simple consequence of the fact that in the deep-well regime, the lowest excitations are those of the approximately Gaussian well of pulsation 2|m 2 | near the nontrivial minimum. The increased (infinite) degeneracy of each level as compared to the free-field case reflects the presence of flat directions in the potential. In particular, there are infinitely many states with almost zero eigenvalue Λ n,n ≈ nλ/(2|m 2 |) |m 2 | 1, which results in large correlation times and lengths for operators in arbitrary nontrivial (vector, tensor, i.e., = 0) representation of the symmetry group. The scalar ( = 0) sector is particular in that the only contribution it receives from this light multiplet is the ground state level Λ 0,0 = 0, which corresponds to the equilibrium PDF and describes the disconnected piece of correlators. The nontrivial time dependence of correlators of scalar operators is thus entirely dictated by the higher levels Λ 2n,0 ≈ 2n|m 2 |, with n > 0, corresponding to the heavy radial directions in the potential, with square mass 2|m 2 |, and thus by small correlation times and lengths relative to the = 0 sectors.
We also show, in Fig. 2, the eigenfunctions R n, corresponding to some of the lowest eigenstates in the case m 2 = 0. We compare the leading-and next-to-leadingorder results (see Appendix B) with N = 2 in order to maximise the difference. The explicit expression of the eigenfunctions at next-to-leading order is given in Appendix A.
We end this Section with some comments concerning the applicability of the present results to arbitrary values of N . It is often the case that the large-N expansion provides a good-qualitative if not quantitative-guide down to small values of N , in particular, in the case m 2 ≥ 0. The case N = 1 has been studied in great detail in the literature [1,2,26] but only very few results exist for N ≥ 2 [27], to which we compare our findings in the Appendix. In Fig. 3 , we show the result of numerically computing the lowest eigenvalue Λ 1,1 for various N in the case m 2 = 0 and we compare with the results of the 1/N expansion. The leading-order result gives a good qualitative description down to rather low values of N and the first, 1/N correction improves the agreement to a quantitative level. We refer the reader to Appendix B for more details and more comparisons.
Double-well potentials, with m 2 < 0, needs to be discussed separately. First, the case N = 1 is qualitatively different because the symmetry at work is discrete. In particular, there are no flat directions in the potential and the relevant physics is governed by tunnelling effects, not Goldstone modes. Another aspect which plays an important role for small values of N is the fact that the relevant radial potential for the present eigenvalue problem is not directly V , but rather W ; see Eq. (5). Indeed, to reformulate the radial eigenvalue equation (9) in terms of a standard one-dimensional problem with an effective potential W eff , one eliminates the single derivative term ∝ R by means of the redefinition ψ(x) = x N −1 2 R(x). This yields, for arbitrary N , where Equation (31) is the standard form of the onedimensional Schrödinger equation with potential W eff , except for the factor N in the first term, which can be absorbed in a rescaling of x. As Markkanen and Rajantie [26] pointed out in the case N = 1 (where = 0, 1 and thus W eff = W ), in the deep double-well limit, W in fact exhibits a three-well structure as a function of ϕ with, in addition to the symmetric wells at ϕ = 0, a third well around the origin ϕ = 0. The resulting spectrum is thus, up to exponentially small splittings due to tunnelling effects, a superposition of the Gaussian spectra from the wells at x = 0 and x = 0, with pulsation |m 2 | and 2|m 2 |, respectively, with the bottom of the central well being upshifted by 3|m 2 |/2 relative to that of the external wells.
The central well remains for arbitrary N > 1 and the potential (32) receives additional centrifugal and geometrical contributions ∝ 1/x 2 . Because of the latter and the centrifugal barrier, the minimum of the central well is slightly shifted away from x = 0 (for N ≥ 3). For increasing N , the potential rapidly approaches the asymptotic form For the quartic potential (19), the position of the central and external wells in the radial direction are given by, respectively, x 2 − = 1/(2|m 2 |) and x 2 + = |m 2 |/λ and the corresponding values of the potential are W eff (x − ) = 3N |m 2 |/4 and W eff (x + ) = 0. We conclude that the excitations of the central well are rapidly lifted relative to that of the external well for increasing N and, hence, decouple in the large-N limit. We show, in Fig. 4, the effective potential (32) together with some of the large-N wavefunctions R n, = e −N v r n, , with r n, given by Eq. (26).

V. CONCLUSIONS
We have proposed a systematic 1/N -expansion of the stochastic approach for quantum fields in de Sitter spacetime which we have applied to O(N )-symmetric models. In its Fokker-Planck formulation, the stochastic approach amounts to solving an equivalent quantum mechanical eigenvalue problem for a single degree of freedom in N dimensions. Various large-N limits of this problem have been considered before in the literature [31][32][33] but, to our knowledge, not the one we have proposed here.
We have performed explicit calculations in the case of a quartic potential, for which we have obtained simple analytical expressions of the eigenvalues and eigenfunctions of the Fokker-Planck operator at leading and nextto-leading orders. These reproduce and generalize our previous partial results in Ref. [3], where a small subset of eigenvalues could be extracted from the calculation of various correlators in the large-N expansion. The eigenvalues and eigenfunctions obtained here can be used to compute a variety of correlators as well as various time and length scales relevant for both phenomenological and fundamental questions, such as spectral indices, relaxation and decoherence times, etc. [28,29,34,35]. The expressions obtained here are nonperturbative in the relevant coupling constant and are thus useful to analyze the various regimes where a perturbative expansion is unavailable, namely, the cases of light interacting fields, m 4 /H 4 λ, and of potentials, with m 2 < 0. Of course the relevant eigenvalue problem can be exactly solved for a wide variety of potentials by numerical means. First results for multifield systems with continuous symmetries have been presented in Ref. [27] and in the present work for massless fields with a purely quartic potential. A detailed investigation of more general potentials for various values of N -in the spirit of Refs. [2,26] for the case N = 1-would be of great interest.
in Eq. (9). The leading-order equation is given by Eq. (18), and was solved in Sec. III. The leading-order eigenfunctions and eigenvalues depend on two quantum numbers n and , and are given in Eqs. (25) and (26).
To keep the formulas simple, we will no write explicitly the dependence in the quantum numbers in the following. We define The next-to-leading-order equation reads The right-hand side can be written in terms of g using the following relations together with the factorization (20), and we end up with where Using the method of variation of constants, we take the following ansatz, r 1 (x) = C(x)r 0 (x), into Eq. (A7), which yields the equation For the quartic potential (19), the function h is a polynomial fraction which can, again, be decomposed into partial fractions. Introducing the notations with the definition (21), the functions r 0 and g can be expressed as and with α + = n− 2 and α − = − n 2 . Note also that We obtain, after some calculations, with the coefficients and where we defined M 2 = m 2 + + m 2 − and a n, = n(3n Note that b n, = a −n, . We verify explicitly the + ↔ − symmetry, obvious from Eqs. (A12) and (A13). In particular, we check that the coefficients α k ↔ β k under the exchange m 2 + ↔ −m 2 − and n − ↔ −n. With the decomposition (A14), Eq. (A9) is readily integrated as with a 1 a free integration constant to be fixed, e.g., by a normalization condition at next-to-leading order. As before, at leading order, possible singularities are related to the zero of the polynomial p + . Remembering that the solution we seek is r 1 = Cr 0 , we see that the last two terms in the first line of Eq. (A23) contribute as α 2 p α+−1 + and α 3 p α+−2 + and are thus potentially singular for n− = 0 and n − = 0, 2, respectively. This singularities are, in fact, absent thanks to the fact that the coefficients α 2 and α 3 vanish for these values of n − . The only possible singular behavior comes from the term ln p + and regularity thus imposes α 1 = 0. This fixes Λ 1 as The expression (27) is obtained using the identities M 2 = √ m 4 + 2λ and m 2 + m 2 − = λ/2. Finally, the corresponding eigenfunction reads The above expressions are valid for all values of the parameters m 2 and λ. To end this Section, we present the explicit formulas for the case m 2 = 0, where m 2 + = m 2 − = λ/2. We have and the various coefficients in the function C(x) read As an illustration, the corresponding eigenfunctions are plotted against the leading-order ones in Fig. 2 for N = 2. In practice, we observe that the next-to-leading-order eigenfunctions provide a pretty good approximation of the numerical results down to N = 2 for the eigenstates we have computed numerically here, namely, R 1,1 and R 2,0 .

Appendix B: Comparison with numerical results
In order to test the validity and convergence of the 1/N expansion in a simple-but nontrivial-case, we solve numerically the eigenvalue equation (17) for a purely quartic potential (19), with m 2 = 0. We first compute the lowest nonzero eigenvalue Λ 1,1 as a function of N and compare with the leading and next-to-leading-order predictions (A26). This is presented in Fig. 3. The first observation is that the leading-order result gives a reasonable estimate of the exact result down to rather low values of N . Furthermore, the next-to-leading-order approximation neatly improves the matters and gives a fairly accurate description of the exact results down to N = 1, where the relative error is 8%. From Eq. (A26), we also observe that the vector ( = 1) sector is the one with the smallest 1/N correction.
To test further the present expansion scheme, we do the same analysis for the lowest nonzero eigenvalue in the scalar ( = 0) sector, namely Λ 2,0 . This is presented in Fig. 5. The leading-order result gives, again, a good estimate of the exact result down to low values of N . As before, the next-to-leading-order approximation quantitatively improves the description, however, for not too small values of N , for which it becomes worse. The relative error reaches 8% for N = 4 and increases up to 25% for N = 1. We have computed the 1/N 2 correction in that case, which reads  Fig. 3 for the lowest nonzero eigenvalue in the scalar sector, Λ2,0. Here, we also show the next-to-nextto-leading-order result (NNLO).
and is also shown in Fig. 5. We see that it greatly improves the description at small N , with a relative error of 10% for N = 1. We end this Section by comparing, when possible, our numerical results to existing ones in the literature. The case N = 1 has been studied in great detail in Refs. [1,2,26] and, recently, Adshead et al have presented first results for continuous symmetries, with N = 2, 3 in the case of a purely quartic potential [27]. The quar-tic couplingλ in that Reference is related to ours as λ = d 2 Ω D+1 λ/(2N ) and the authors use the quantum numbers k = (n − )/2 and to label the eigenstates. Their definition of the eigenvaluesΛ k, is related to ours asΛ with d = 3 and Ω 5 = 8π 2 /3. This holds for N ≥ 2. In the case N = 1, there is no angular momentum in field space and only = 0, 1 are permitted. One has For N = 1, we find Λ 1,1 / √ λ = 0.9693 and Λ 2,0 / √ λ = 3.1532. This translates intoΛ 1 = 0.0891 andΛ 2 = 0.2897, which agrees with the known results [1,2,26,27]. For cases with continuous symmetries, we find Λ 2,0 / √ λ = 2.8133 for N = 2 and Λ 2,0 / √ λ = 2.7296 for N = 3, which giveΛ N =2 1,0 = 0.3656 andΛ N =3 1,0 = 0.4344, respectively, in agreement with the results of Ref. [27].