$O(N)$ model in Euclidean de Sitter space: beyond the leading infrared approximation

We consider an $O(N)$ scalar field model with quartic interaction in $d$-dimensional Euclidean de Sitter space. In order to avoid the problems of the standard perturbative calculations for light and massless fields, we generalize to the $O(N)$ theory a systematic method introduced previously for a single field, which treats the zero modes exactly and the nonzero modes perturbatively. We compute the two-point functions taking into account not only the leading infrared contribution, coming from the self-interaction of the zero modes, but also corrections due to the interaction of the ultraviolet modes. For the model defined in the corresponding Lorentzian de Sitter spacetime, we obtain the two-point functions by analytical continuation. We point out that a partial resummation of the leading secular terms (which necessarily involves nonzero modes) is required to obtain a decay at large distances for massless fields. We implement this resummation along with a systematic double expansion in an effective coupling constant $\sqrt\lambda$ and in 1/N. We explicitly perform the calculation up to the next-to-next-to-leading order in $\sqrt\lambda$ and up to next-to-leading order in 1/N. The results reduce to those known in the leading infrared approximation. We also show that they coincide with the ones obtained directly in Lorentzian de Sitter spacetime in the large N limit, provided the same renormalization scheme is used.


Introduction
The study of interacting quantum fields in de Sitter geometry is of interest for a variety of reasons. In inflationary models, interactions could lead to non-Gaussianities in the cosmic microwave background. Quantum effects could also contribute to the dark energy, and explain, at least partially, the present acceleration of the universe. From a conceptual point of view, being maximally symmetric, de Sitter geometry allows for a number of explicit analytic calculations that illustrate the role of the curved background on the dynamics of the quantum field, and also the backreaction of the quantum fields on the geometry.
Although at first sight quantum field theory (QFT) in de Sitter geometry should be simpler than on other cosmological backgrounds, the exponential expansion of the metric produces an effective growth in the couplings. Indeed, when considering a λφ 4 theory, a diagram with L loops will be proportional to (λH 2 /m 2 ) L , where m is the mass of the quantum field and H the Hubble constant [1]. Therefore, the usual perturbative calculations break down for light (m 2 H 2 ) quantum fields. In the massless limit, minimally coupled scalar fields do not admit a de Sitter invariant vacuum state due to infrared (IR) divergences, and the two-point functions do not respect the symmetries of the classical theory [2].
There have been several attempts to cure the IR problem, all of them introducing some sort of nonperturbative approach. The well known stochastic approach [3] assumes that the long distance fluctuations of the quantum field behave as a random classical field, that satisfies a Langevin equation. From the corresponding Fokker-Planck equation one can obtain finite and de Sitter invariant correlation functions even for minimally coupled massless fields and fields with negative mass-squared, provided there are interaction terms that stabilize the potential. The formalism has proven very useful in understanding the IR effects and, in particular, the spontaneous symmetry breaking phenomena (see for instance [3,4]). However, it does not provide a systematic way of accounting for the interactions of the short distance fluctuations, which are generically neglected, nor a justification for the classical treatment of the field beyond the leading IR approximation.
In the context of QFT in Lorentzian de Sitter spacetime, many of the nonperturbative approaches are based in mean field (Hartree) [5,6] or large N approximations [7], which in turn can be formulated through the 2PI effective action (2PIEA) [8]. Other approaches consider the analysis of the Schwinger-Dyson equations (or their nonequilibrium counterpart, the Kadanoff-Baym equations), combined with a separation of long and short wavelengths using a physical momentum decomposition [9]. There are also calculations based on the dynamical renormalization group, and on the exact renormalization group equation for the effective potential [10]. In relation to the development of a method for computing correlation functions in a systematically improvable way, the difficulties one finds with these approaches are mainly of technical nature.
An alternative and simpler approach emerges for quantum fields in Euclidean de Sitter space. As this is a compact space (a sphere of radius H −1 ) the modes of a quantum field are discrete, and the origin of the IR problems can be traced back to the zero mode [11]. This important observation suggests itself the solution: the zero mode should be treated exactly while the nonzero modes (UV modes in what follows) could be treated perturbatively. Using this idea, it has been shown that in the massless λφ 4 theory the interaction turns on a dynamical squared mass for the field, that in the leading IR limit is proportional to √ λH 2 . Indeed, in the Euclidean space, this mass cures the IR problems. Moreover, a systematic perturbative procedure for calculating the n-point functions has been delineated in Ref. [12], where it was shown that, for massless fields, the effective coupling is √ λ instead of λ. It has been pointed out that this procedure together with an analytical continuation could be used to cure some IR problems also in the Lorentzian de Sitter spacetime, and in particular to obtain n-point functions that respect the de Sitter symmetries. However, so far explicit calculations have been restricted to obtaining corrections to the variance of the zero mode, which has no analog quantity in the Lorentzian de Sitter spacetime. In particular, an explicit calculation of the inhomogeneous two-point functions of the scalar field that can lead to two-point functions respecting the de Sitter symmetries after analytical continuation is still lacking. To perform that explicit calculation is one of the main goals of this paper.
The relation between the different approaches (Lorentzian, Euclidean, stochastic), has been the subject of several works. For example, recently, it has been shown that the stochastic and the in-in Lorentzian calculations are equivalent at the level of Feynman rules for massive fields, when computing equal-time correlation functions in the leading IR approximation [13], in agreement with previous arguments [14]. The evaluation of the field variance in the large N limit gives the same result in the three approaches, up to next to leading order (NLO) in 1/N , in the IR limit [9]. The connection between the Feynman diagrams for computing correlation functions in the in-in formalism and the analytically continued ones obtained in the Euclidean space has been described in detail in [15,16] for massive fields. Similar arguments have been used in [17] to generalize the connection to the massless case when the zero mode is treated nonperturbatively as proposed in [11]. In the latter case, the inhomogeneous IR behavior of the the correlation functions is still unclear (they could grow at large distances [17]), and to understand this it is necessary to go beyond the leading IR approximation. It is worth to highlight that beyond the leading IR approximation a comparison of the results obtained by different approaches necessary requires the use of equivalent regularization methods and renormalization schemes.
In this paper we will pursue the approach of Ref. [12], providing a generalization to the case of O(N ) scalar field theory, and including a discussion of the renormalization process. We will present a detailed calculation of the corrections to the two-point functions up to second order in the parameter √ λ. As we will see, the zero mode part of the twopoint functions, that also receives UV corrections, determines the quadratic part of the effective potential. Being positive for all values of d and N , it implies that spontaneous symmetry breaking does not occur in the O(N ) models for λ 1. We will check that, in the leading large N limit, the Euclidean results are fully consistent with the Lorentzian ones, even beyond the leading IR approximation. We will also show that the corrections of order 1/N of our result coincide in the leading IR approximation with the ones of Ref. [9], which are the most precise results known so far for this model and were obtained working in the leading IR approximation and directly in the framework of the QFT in Lorentzian spacetime. Our results improve on those by including systematically the corrections coming from the interactions of both IR and UV sectors.
The paper is organized as follows. In Section II we present our model and describe the systematic perturbative calculation of Ref. [12]. We compute the homogeneous part of the two-point functions, related to the quadratic part of the effective potential, in the leading IR approximation (i.e. neglecting the contributions coming from the UV modes). In Section III we evaluate the full two-point functions of the theory, up to corrections of second order from the interaction with the UV modes, including the counterterms needed for renormalization. We find out that, although this framework deals with the IR divergences by giving a dynamical mass to the zero modes, the behavior of the two-point functions at large distances still has IR problems for massless fields. In Section IV we analyze in detail the massless case. We perform a large N expansion of the results of Section III and compare the Euclidean and Lorentzian results. Section V is concerned with an extension of the nonperturbative treatment to perform a resummation of the leading IR secular terms, in order to recover the proper decay of the two-point functions at large distances. This resummation is combined with a systematic expansion both in √ λ and 1/N . Finally, in Section VI we describe the main conclusions of our work. Several Appendices contain some details of the calculations.

Euclidean de Sitter Space
In this section we describe the methods of QFT in the Euclidean de Sitter space developed in [11] and [12] and generalize them for the case of a O(N )-symmetric model with Euclidean action Euclidean de Sitter space is obtained from Lorentzian de Sitter space in global coordinates by performing an analytical continuation t → −i(τ − π/2H) and a compactification in imaginary time τ = τ + 2πH −1 . The resulting metric is that of a d-sphere of radius H −1 where θ = Hτ . Due to the symmetries and compactness of this space, the field can be expanded in d-dimensional spherical harmonics and then the free scalar propagator of mass m is, in the symmetric phase, where the superscript indicates the mass. The L = 0 contribution G (m) 0 = |Y 0 | 2 H d /m 2 is clearly responsible for the infrared divergence in the correlation functions of the scalar field for m 2 → 0. We split φ a (x) = φ 0a +φ a (x) in order to treat the constant zero modes φ 0a separately from the inhomogeneous partsφ a (x). This prompts to separate the propagator as well, where nowĜ (m) has the property of being finite in the infrared (m 2 → 0). The interaction part of the action takes the following form: (2.6) Here |φ 0 | 2 = φ 0a φ 0a and V d is the total volume of Euclidean de Sitter space in d-dimensions, which thanks to the compactification is finite and equal to the hypersurface area of a dsphere where Γ is Euler's Gamma function. The explicit form ofS int will be written below.
In order to compute the quantum correlation functions of the theory we define the generating functional in the presence of sources J 0a andĴ a , where we introduced the shorthand notation is defined as the generating functional of the theory with the zero modes alone. This part gives the leading infrared contribution and can be exactly computed in several interesting cases following Ref. [11]. Note that, as the zero modes are constant on the sphere, their kinetic terms vanish, and Z 0 [J 0 ] involves only ordinary integrals.
The effective potential gives valuable information about how the quantum fluctuations around a background fieldφ influence its behavior. We are interested in particular in the generation of a dynamical mass due to quantum effects. Up to quadratic order it can be shown that (see Appendix A), This is an exact property of the Euclidean theory valid for all N and λ, which shows that the dynamical mass is related to the inverse of the variance of the zero modes as (2.10) At the leading infrared order the interaction between the infrared and ultraviolet modes in Eq. (2.8) can be neglected, and For a vanishing tree level mass m = 0, the integrals on the right-hand side can be computed exactly leading to a dynamical mass, (2.12) For N = 1, we recover the result of [11], which is also the one from the stochastic approach [3], where we evaluated in d = 4, V 4 = 8π 2 /3H 4 . In the following section, we will compute the two-point functions of the full scalar field including up to the second perturbative correction coming from the UV modes.

Corrections from the UV modes to the two-point functions
Corrections to the leading order result come from expanding the exponential withS int in Eq. (2.8). The explicit expression for the interaction part of the action that involves the UV modes is where A abcd is the totally symmetric 4-rank tensor The terms linear inφ that would appear when splitting both the mass and the interaction terms vanish, since d d x √ gY L (x) = 0 for L > 0.
As a guiding principle for computing the perturbative corrections in this section, we recall that for a massless minimally coupled field, φ 2p 0 0 ∼ λ −p/2 (see Eqs. (2.10) and (2.12), and Eq. (4.1) below), and therefore we will have a perturbative expansion in powers of √ λ (in contrast, the correlation functions of the zero modes start at λ 0 when m = 0, and therefore the order at which each perturbative term contributes changes with respect to the massless case). With this in mind, here we will keep terms that will be at most of order λ when m = 0.
The first correction to the generating functional comes from expanding the exponential linearly and keeping the first term of Eq. (3.1). Following the standard procedure, the generating functional at NLO reads .
The next to NLO (NNLO) correction has two contributions. The first one is given by the square of the interaction term considered before, as it comes from expanding the exponential up to quadratic order. The second one is the last term in Eq. (3.1) at linear order. These NNLO contributions to the generating functional are given by , (3.4) and , (3.5) respectively. Here, and in what follows, we use the notation as a shorthand in the functional derivatives. We split the two-point functions of the total fields φ a into UV and IR parts where the cross-terms vanish by orthogonality. In what follows we compute each part separately.

UV part of the two-point functions
Starting with the UV part, we calculate the two point functions ofφ a by taking two func- where the factor Z[0, 0] −1 takes care of the normalization of the interacting theory: Here [Ĝ (m) ] denotes the coincidence limit of the free UV propagator with mass m, which is independent of x by de Sitter invariance, and we used A abcd δ cd = (N + 2)δ ab . Furthermore, to arrive at this expression we have used Eqs. (B.1) and (B.2) to write the derivatives of Z f [Ĵ] in terms of free propagators, relying on the fact that it is a free generating functional. On the other hand, for the derivatives of Z 0 [J 0 ], we have with φ 0a φ 0b 0 the exact two-point functions of the zero modes in the absence of the UV modes. By exact we mean that we include the self-interaction nonperturbatively. For convenience, in the last equality of Eq. (3.9) we expressed the result in terms of its trace with respect to the internal O(N ) indices. In general, any of the traced exact n-point functions of zero modes at leading order can be expressed in terms of ordinary integrals. For even n = 2p we have while they vanish for odd n.
To evaluate the second derivative of Z[J 0 ,Ĵ], we can set J 0 = 0 in Eq. (3.7) and compute .
it can be seen that Eq. (3.13) can be made finite by a mass counterterm of the form The expression for the UV part of the propagator can be simplified considerably using that the integrals of free UV propagators in Euclidean space in Eq. (3.13) can be expressed in terms of derivatives of a single propagator with respect to its mass. This is shown in the Appendix D: and thus Eq. (3.13) reads, after renormalization, Naïvely, this result looks like a Taylor expansion of a free UV propagator with respect to the mass squared, around the mass parameter m 2 . However, as we will see later, while this is indeed the case in the large-N limit, it does not hold anymore when the 1/N corrections are included.

IR part of the two-point functions
We now calculate the IR part of the two-point functions by taking two derivatives of the It is useful to first setĴ = 0 and take the derivatives afterwards: where we made use of Eqs. (3.9) and (3.10), and that In the last term of Eq. (3.20) we can make the replacement . This shows that this expression has, besides [Ĝ (m) ], another divergent quantity ∂[Ĝ (m) ]/∂m 2 . We split the latter as before and following the details of the Appendix C we obtain the renormalized expression. The renormalization involves a new counterterm to compensate for this divergence, namely Finally, we can write down the renormalized two-point functions for the zero modes . (3.24) The last equality follows after interpreting the corrections as a modification to the mass m 2 dyn (IR) of the zero modes which, as mentioned before, determines the curvature of the effective potential. Eqs.

Massless fields
A case of great interest is when the fields are massless at tree level, m = 0, as it is in this case in which the perturbative treatment becomes ill-defined. The nonperturbative treatment of the zero modes ensures that these modes acquire a dynamical mass, avoiding the IR divergence associated to the free two-point functions in the massless limit. This can be verified by checking that Eq. (3.24) remains finite in this limit. Indeed, the n-point functions of the zero modes can be exactly computed to be which exhibit no IR divergences. This equation shows a scaling of the form φ 0 ∼ λ −1/4 , making the perturbative expansion of the UV modes to be in powers of √ λ. It is worth to note that, when considering the two-point functions of the UV modes, Eq. (3.18), the free UV propagators that build up this expression will now become massless. After performing the analytical continuation to the Lorentzian spacetime, this leads to an IR enhanced behavior at large distances, i.e. secular growth (see Appendix E). The reason for this is that so far we have only given mass to the zero modes, while the UV modes remain massless, as sketched diagramatically in Appendix F. This situation can be dealt with by improving the result Eq. (3.18) by further resumming a given subset of diagrams that give mass to the UV propagators appearing in that expression. We will focus on this point in Section 5.

The 1/N expansion
In Euclidean de Sitter space, after dealing with the zero modes, the calculations of the npoint functions can be done by means of a perturbative expansion in powers of √ λ. When λ is sufficiently small, being a compact space, this perturbative expansion is valid for any set of points and for any value of N . So an expansion in 1/N is certainly unnecessary in that case. However, we are ultimately interested in computing quantities in the Lorentzian de Sitter spacetime. For this, as will become clear in the next section, a double expansion in √ λ and 1/N turns out to be crucially convenient in order to obtain a tractable perturbative expansion that remains valid at large distances. With this in mind, we perform an expansion in 1/N of the results obtained in the previous section. In order to compare with known results appearing in the literature, it is sufficient to remain at NLO in 1/N . At the same time, the known Lorentzian results that are nonperturbative in the coupling constant will have to be expanded in powers of √ λ to bring them to the same precision. The Euclidean results for the two-point functions at order λ and order 1/N are obtained by inserting Eq. (4.1) into Eqs. (3.18) and (3.24), and then expanding in powers of 1/N . We arrive at the following expressions: and for the UV and IR contributions, respectively. It is worth to stress that the leading contribution in the limit N → ∞ of Eq. (4.2) is compatible with a Taylor expansion of a massive UV propagator, with an UV dynamical mass given by up to order λ. Furthermore, the IR dynamical mass m 2 dyn (IR) read from Eq. (4.3) also coincides up to that order with m 2 dyn (U V ), and therefore the whole propagator can be interpreted as a free de Sitter propagator with a dynamical mass m 2 dyn . Beyond the LO contribution in 1/N , this is not true, since in Eq. (4.2) the coefficient of the second derivative is no longer half the square of that of the first derivative, and the two-point functions has a more complicated structure.

Comparison with Lorentzian QFT: dynamical mass
There are several nonperturbative approaches in the Lorentzian QFT. Here we consider the 2PIEA formulation, in which both the mean value of the fieldφ(x) and the exact twopoint functions G (m dyn ) (x, x ) are treated as independent degrees of freedom. A detailed description of this method is given elsewhere. In this framework it is possible to obtain the exact result for the two-point functions in the large-N limit. For this, we follow [6,18] and write down the exact equations of motion in the large-N limit, (4.6) In this limit, the equations become local. Eq. (4.6) corresponds to that of a free propagator in de Sitter spacetime with mass m 2 dyn satisfying a self-consistent gap equation, where the counterterms have to be suitably chosen to cancel the divergences of the coincidence limit of the propagator in the right hand side. We expand [G (m dyn ) ] in powers of where the dots stand for terms of order O(m 4 dyn /H 4 ) and can be neglected since we are interested in the small mass case m 2 dyn H 2 . We take advantage of the fact that the coincident propagator is exactly the same for both the Lorentzian and Euclidean theories. Therefore we use the same notation as for our Euclidean calculations for the different parts in the expansion. The necessary counterterms are (4.10) When expanded at the leading order in λ, these counterterms coincide with those used in the Euclidean calculation, Eqs. (3.15) and (3.23), when the latter are expanded for small masses and the large-N limit is taken. The renormalized gap equation is then, in the symmetric phaseφ = 0: This is a quadratic algebraic equation for m 2 dyn , whose positive and physically relevant solution is , (4.12) where we have defined (4.14) The dynamical mass is finite for m = 0, and exact in the large-N limit.
In order to draw a comparison with the Euclidean results, we expand the Lorentzian expression up to order √ λ, The corresponding Euclidean calculation is given in Eq. (4.3). We see that, at leading order in 1/N , the dynamical masses computed in both approaches coincide. When going beyond the leading order in 1/N , the previous approach is no longer valid, since the propagator cannot be described as a free propagator with a dynamical mass given by the gap equation, Eq. (4.6), and the full Schwinger-Dyson equations must be solved instead. The complete Lorentzian results that take into account the UV modes and the renormalization process are technically involved and still unknown. However, in Ref. [9] the authors were able to obtain results up to the NLO in the 1/N expansion which are valid at the leading order in the IR. We refer the reader to their paper for the details. Basically, exploiting the de Sitter symmetries, the scalar fields are split into two sectors: an IR one, formed by modes with physical momenta smaller than a critical scale µ IR H, and the one containing the remaining modes. The interactions of the IR modes are taken into account, but the ones of the remaining modes are neglected. In this approximation, the authors have obtained a self-consistent solution of the Schwinger-Dyson equations for the IR two-point functions that is valid at NLO in 1/N and at leading order in the IR. In particular, the result for the dynamical mass can be written as [9]: where the dots stand for corrections that are higher order in either √ λ or 1/N (which cannot be unambiguously computed in the approximation considered, since they depend on the arbitrary IR scale µ IR ). Clearly, this result agrees with our Eq. (4.3) at leading order in √ λ.

Resumming the leading IR secular terms to the two-point functions
In section 4 we discussed the case of massless fields and took notice that the calculated two-point functions of the UV modes, given in Eq. (3.18), would be expressed in terms of massless UV propagators. The problem is then that the analytically continued correlation function does not decay at large distances, as happens for massive fields. The reason for this being that the UV modes remained massless in this framework. In this section we will extend the nonperturbative treatment, with the aim of resumming the leading IR secular terms. In order to achieve this, we need to perform a resummation of diagrams to give mass to the UV propagators present in Eq. (3.18). As we will see, it will be enough to resum only a subclass of diagrams: those coming from the interaction term that is quadratic in both φ 0 andφ, which will be treated nonperturbatively, while still perturbing on the remaining terms. We start by rewritting the generating functional Eq. (2.8) by grouping this term with the other terms quadratic inφ, as part of the "free" generating functional of the UV modes, where now the "free" UV propagatorĜ ab (φ 0 ) has a φ 0 -dependent mass, with m 2 ab (φ 0 ) = (λ/2N )A abcd φ 0c φ 0d and whereS int =S int − S (2) int has the remaining interaction terms that should be treated perturbatively. The normalization factor N ensures that Z[0, 0] = 1.
In order to compare with the results of the previous sections, it will be enough to keep terms up to order λ. Therefore, it is necessary to include perturbatively only the first correction coming from the term, and so the generating functional reduces to Now we proceed to calculate the connected two-point functions of the UV modes, which we split in two contributions, according to the interaction term that we are treating perturbatively. Approximating Z (1) in Eq. (5.6) by the first term of Eq. (5.5) we obtain for the first contribution, Here it is important to note that both square roots of the determinant of the propagator will not cancel each other out due to the integral over φ 0 . This is crucial to obtain the correct result. In order to evaluate this formal expression, we must expand both the numerator and the denominator in powers of m 2 ab (δ/δJ 0 ). We start with the numerator: The key point is that our resummation needs only include the (infinite) subset of contributions that modify the UV propagator at separate points, while the determinant ofĜ ab has no IR problems and can be safely evaluated at m = 0. Therefore, in the above series it is enough to consider only terms with as many derivatives acting on the determinant as needed for a given precision in λ, since each time we increase the number of those derivatives we pick up a factor of λ φ 2 0 ∼ √ λ. In our case, we shall keep terms with zero and one derivatives of the determinant.
If first we consider the subset of terms with no derivatives of the determinant, we have a series for the connected UV propagator, We see that each term of the series corresponds to a connected diagram with p insertions of (λ/2N )A ijkl φ 0j φ 0k and two external legs, as shown in Fig. 1. The contractions of the O(N ) indices amount to an overall δ ab and a factor depending both on N and p. To calculate this factor, we recall that each A ijkl has 3 terms made of pairs of Kronecker deltas, so for each p there will be 3 p terms. Of these, one term will be proportional to while the other 3 p − 1 terms will leave two powers of φ 0 untraced, Therefore, the series in Eq. (5.9) becomes This series can be resummed order by order in 1/N , as we will see in the next section. Before doing this, let us first deal with the other subset of terms that we need to consider up to order λ 1 , namely, those in which there is one derivative acting on the determinant, ∂ detĜ In this case the UV propagators in the diagrams are not all connected among themselves, but rather through their interaction with the zero modes. Indeed, these diagrams are composed of a single bubble with one insertion of (λ/2N )A ijkl φ 0j φ 0k (factorized in the previous expression) times a connected part with p − 1 insertions, as depicted in Fig. 2. The prefactor is simply proportional to δ kl , while the series now contains only connected diagrams, and the same argument as before applies. Therefore, relabeling the summation index p = l + 1, we obtain ∂ detĜ Next we consider the denominator of Eq. (5.7). As we discussed, there is no need to resum the determinant since it has no external legs. Therefore we can treat it perturbatively in λ. To be consistent with what we did with the numerator, we should keep the first two terms, that is, with zero and one derivatives acting on the determinant, The combination in big parentheses can be easily computed (see Appendix G) to give where, although there is a massless propagator, there are no infrared problems since it is evaluated at its coincidence limit. To deal with the UV divergence, the renormalization can be carried out by introducing (at this order) the mass counterterm δm 2 given in Eq. (3.15) at m = 0, both in the action of the zero modes as well as for the UV modes. Then, when performing the above mass expansions, the diagrams now can be constructed with both lines with insertions of (λ/2N )A ijkl φ 0j φ 0k and of δm 2 δ ij . As we did with the derivatives of the determinant, we need only keep a finite number of δm 2 δ ij insertions, given that δm 2 ∼ λ. In our case, we need only one, but in general we will need as many as derivatives of the determinant we have. At higher orders in λ it is also necessary to include a δλ counterterm.
Putting everything together, distributing, and keeping the terms up to the corresponding order in λ, we arrive at . Now we look at the remaining term coming form the perturbative correction, Once again, we must expand this in powers of the mass in order to evaluate the J 0derivatives acting on Z 0 [J 0 ] an resum order by order in 1/N . As before, the only part we should be concerned with resumming is the first derivative of the UV propagator at separate points, while the other factors are well behaved for m = 0. The renormalization of the coincidence limit comes from the inclusion of the δm 2 counterterm. In contrast with the previous calculation, we have no need to keep any terms with derivatives of the determinant, since they would only contribute at higher orders in λ. Then, keeping only the contributions with derivatives over ∂Ĝ (m) (x,x ) where once again the counting of the N -dependent factor relies on the fact that each term in the series corresponds to a connected diagram.

1/N expansion at NLO
As we mentioned, the series above can be resummed order by order in 1/N , for which we need to expand the summands of the various series. Up to order N −1 we have Using the first of these we can resum the first Taylor series of Eq. (5.17): where m 2 dyn,0 = λ 2V d , and the second series: Inserting these expressions into Eq. (5.17) we get, All the UV propagators at separated points now have a mass squared of order √ λ, giving them the correct behavior at large separations. A noteworthy observation is the presence of some propagators whose squared mass is three times that of the others, something which could not have been anticipated from the perturbative result. Also, the only instance of a massless UV propagator in the previous expression has its coincidence limit taken, and it is therefore just a finite constant factor with no IR issues. Now we turn to the remaining part, Eq. (5.19), where the series can be resummed as done for Eq. (5.17). However, since the contribution at NLO in 1/N is also higher order in λ with respect to the LO one, we can just keep the latter, The full result for the connected two-point functions of the UV modes up to order λ and N −1 , with a partial resummation of the infinite subset of diagrams, is obtained by combining Eqs. (5.24) and (5.25): It was verified as a cross-check that this result reduces to the perturbative one Eq. (3.18), upon expanding the latter at NLO in 1/N (see Appendix H). The large distance behaviour of the two-point functions ultimately depends on the masses of the free propagators that build up the expression, m 2 dyn,0 and 3m 2 dyn,0 , which determine how fast it decays. However, the true exponent could be other, the difference lying beyond the precision of this result.

Comparison with Lorentzian QFT: two-point functions
The Lorentzian calculations of Ref. [9] for the two-point functions are expressed differently, as a sum of two free (full) propagators, with masses which are different to ours. The coefficients are given by c − = 5/16N and c + + c − = 1. Therefore, in order to draw a comparison between both results, they are better expressed in terms of UV propagators and its derivatives around a common mass parameter, which we choose to be m 2 dyn,0 . Also, the full propagators of Eq. (5.27) must be first separated into their IR (constant) and UV parts, for which it is better to think of that expression as being analytically continued to Euclidean space.
Expanding the Euclidean result Eq. (5.26) in this way and dropping higher order terms we obtain while doing the same with the inhomogeneous part of the Lorentzian result Eq. (5.27) yields all but those terms proportional to [Ĝ (m) ] ren . However, it is expected for these contributions to be missing in the Lorentzian result of Ref. [9], given the nonsystematic way the interactions among IR and UV sectors are treated there. Indeed, we can only trust that result up to the leading IR order, √ λ, and up to NLO in 1/N , in which case both results are compatible.
The convenience of the Euclidean approach becomes evident as it provides a systematic, order by order expansion in both √ λ and 1/N , while also allowing for the inclusion of partial resummations that cure the IR effects at large separations.

Conclusions
In this paper we considered an interacting O(N ) scalar field model in d-dimensional Euclidean de Sitter spacetime, paying particular attention to the IR problems that appear for massless and light fields. We extended the approach of Refs. [11,12] to the O(N ) model. The zero modes are treated exactly while the corrections due to the interactions with the UV modes are computed perturbatively. The calculation of the two-point functions of the field shows that the exact treatment of the zero modes cures the IR divergences of the usual massless propagator: the two-point functions becomes de Sitter invariant. However, the NLO contains derivatives of the free propagator of the UV modes. Although the massless UV propagator is de Sitter invariant, its Lorentzian counterpart exhibits a growing behavior at large distances, invalidating the perturbative expansion in this limit. This problem can be fixed in the leading order large N limit by resumming the higher order corrections: one can show that the final result corresponds to a two-point functions of a free field with a self-consistent mass.
In order to alleviate the behavior of the correlation functions in the IR limit, we performed a resummation of a class of diagrams that give mass to the UV propagator. After this resummation, higher order corrections can be systematically computed in a perturbative expansion in powers of both √ λ and 1/N . We presented explicit results up to second order in √ λ and NLO in 1/N . In the leading large N limit, the results can be computed exactly, both in the Euclidean approach and also working directly in the Lorentzian spacetime (which corresponds to the Hartree approximation [5,6]). We derived the corresponding results in the Lorentzian spacetime paying particular attention to the renormalization process. We showed that the Euclidean approach reproduces the Lorentzian results in this large N limit, provided the same renormalization scheme is used. As we mentioned before, to our knowledge, the most precise results previously obtained for this model are those presented in Ref. [9], which are also valid up to the NLO in the large N expansion, but only at the leading IR order. We have also shown that our results coincide with the ones of Ref. [9], when expanded up to the corresponding order. Beyond the leading IR order, as we have emphasized along the main text, a consistent treatment of the UV sector becomes necessary. The use of the Euclidean path integral (which is simpler than its in-in counterpart) together with the double perturbative expansion (in √ λ and 1/N ) performed in our calculations, allowed us to further include the contribution of the UV modes and to consistently take into account the renormalization process. Moreover, in this framework, the precision of the calculation can be systematically improved by computing higher order corrections.
There are many directions in which the present work can be extended. An interesting generalization of our calculations that we leave for future work is the study of the corrections in the case of a tree level double-well potential. As customarily done in Minkowski spacetime, it would be definitively interesting to use the Euclidean approach, together with the appropriate analytical continuation, to perform calculations in nonstationary situations in a fixed de Sitter background spacetime. For instance, to develop a systematic way of computing corrections to the effective action or to the expectation value of the energy-momentum tensor of the quantum scalar fields, which is important for studying the backreaction of the quantum fields on the dynamics of the spacetime geometry. Furthermore, it would be valuable to extend the Euclidean techniques to other models, such as theories with derivative interactions, with fermionic and/or gauge fields, and to study metric perturbations around a de Sitter (or quasi de Sitter) background.

Acknowledgments
The work of FDM and LGT has been supported by CONICET and ANPCyT.
LGT was also partially supported by the ICTP/IAEA Sandwich Training Educational Programme and UBA. DNL was supported by ICTP during the initial stage of this work. FDM would like to thank ICTP for hospitality during part of this work. Some algebraic computations were assisted by the software Cadabra [19].

A Effective Potential
In this Appendix we examine the quadratic part of the effective potential and relate it to the variance of the zero modes. We start by defining the effective action for this theory, ) the generating functional of connected diagrams, and wherē define the "classical" fields. The Effective Potential is obtained by evaluating the effective action at a constant field, that isφ = 0, which in turn demands thatĴ = 0, and then dividing by the space volume V d : With the purpose of calculating the quadratic term of V ef f (φ 0 ) as a function ofφ 0 , we perform the following expansion, where the linear term vanishes forφ 0 = 0, as it can be seen by differientiating Eq. (A.4) with respect toφ 0 , δΓ[φ 0 , 0] δφ 0a = −J 0a , (A. 6) and taking into account that in the symmetric phase, the mean fieldφ 0 vanishes if and only if J 0 = 0. Taking another derivative of the previous expression but now with respect to J 0 , we obtain where we used the chain rule for the second equality. Now, differientiating Eq. (A.2) with respect to J 0 gives, which, inserted in the previous expression leads to the conclusion that We now have to evaluate forφ 0 = 0 (J 0 = 0), allowing us to identify the exact two-point functions of the zero modes φ 0a φ 0b . In the symmetric phase we expect any rank-2 tensor with respect to the internal O(N ) indeces to be proportional to the identity δ ab . Therefore inverting the previous quantity is straightforward, where we have expressed the result in terms of the variance of |φ 0 |, i.e. φ 2 0 = δ ab φ 0a φ 0b . Finally, we replace this last expression in Eq. (A.5) and we divide by V d , in order to obtain the effective potential up to quadratic order, Eq. (2.9).

B Functional derivatives ofẐ f [Ĵ]
Using thatẐ f [Ĵ] is a free generating functional, its functional derivatives evaluated atĴ = 0 can be easily expressed in terms of the free UV propagatorĜ(x, x ). The only somewhat tricky part is to keep track of the O(N )-indeces. The useful expressions are In the case of the sixth derivative, it is not necessary to write down the most general expression for six different points, since we only need particular cases with some of them evaluated in coincidence. The two particular cases we need are and

C Renormalization
The renormalization process is performed by the addition of two counterterms in the action, . It is safe to assume that, as in the usual perturbative case, δm 2 ∼ λ and δλ ∼ λ 2 , therefore at NNLO we need to consider terms with δm 2 , (δm 2 ) 2 and δλ. This leads to the following new contributions to the generating functional, Tracking these terms into the calculation of both the UV and IR two-point functions, produce the following contributions: to be added to Eq. (3.13), and to be added to a previous step of Eq. (3.24) (not shown), which is just like Eq. (3.24) with the ren and f in labels removed. With these additions coming from the counterterms, it is straightforward to see that the choices Eqs.

D Integrals of the Euclidean UV propagator
The UV propagator can be written aŝ (D.1) where we have used the orthogonality relation of the spherical harmonics showing that when r → +∞Ĝ Analogously, it can be shown that the derivatives with respect to the mass pick up further powers of log(r): (E.8)

F Double expansion and diagrammatics (the massless case)
In this Appendix we focus on the corrections to the UV propagator and the associated diagrammatics for the case of fields with vanishing tree-level masses. The zero modes are treated nonperturbatively using the exact Euclidean path integrals in the absence of the nonzero modes. On the one hand, from these path integrals we have learnt that each power of φ 0,a scales as λ −1/4 , while eachφ a does not add any factor of λ (i.e., it counts as λ 0 ). On the other hand, we know that in the massless limit (m 2 → 0) the free propagators in the Lorentzian spacetime increase at large distances, and that a derivative with respect to m 2 of the propagators adds a logarithmic-growth factor (see Appendix E). For assessing the importance of each diagram after the analytical continuation to the Lorentzian spacetime we need to understand how the behavior of each of them at large distances would be. For this, it is very useful to note that equations like (3.16) and (3.17) hold in general, namely ... Therefore, by representing each logarithmic-growth factor by y (i.e., using the notation of Appendix E, y ∼ log r) we see the right hand side scales as y k−1 .
In order to draw Feynman diagrams, we use a dash line for the variance of the zero modes (which is always computed nonperturbatively using the exact path integral for the zero modes in isolation and scales as 1/ √ λ) and a solid line for the free UV propagator (which depending on the diagram may contribute with a factor of y). Hence, using that for each vertex we have a factor of λ/N , one can conclude that the diagram in Fig. 3 scales as y( √ λy) and becomes of the same order as the free UV propagator (which scales as y) when √ λy ∼ 1, indicating a break down of the perturbation theory when the distance is sufficiently large. Actually, there are more diagrams that also become of the same order in that case, as, for instance, those shown in Fig. 4. Indeed, it is not difficult to see that the one with n-vertices scales as y( √ λy) n , representing a relative correction that goes as ( √ λy) n . We were able to perform the resummation of these diagrams in section 5. Now, in order to see that the rest of the contributions can be treated perturbatively, let us analyze each of them separately, according to the type of vertex. Let us start with other corrections associated to the bi-quadratic interaction term of Eq. (5.1). The leading order diagram involving this vertex which is not included in the resummation is the one on the left in Fig. 5, and gives a relative correction that scale as √ λ( √ λy). Then, for λ sufficiently small it can be considered as a small correction to the ones of Fig. 3 at all times. The other type of correction that also contributes up to the NLO in 1/N involves the last interaction term in Eq. (3.1) and the corresponding diagram is the one in Fig. 5 to the right and gives a correction relative to the free propagator that also goes as √ λ( √ λy).  . Feynman diagrams contributing to the UV propagator, which become of the same order as the one in Fig. 3 when √ λy ∼ 1. The first one scales as y( √ λy) 2 , the second one as y( √ λy) 3 and so on.