Long distance behavior of $O(N)$-model correlators in de Sitter space and the resummation of secular terms

We analyze the long distance behavior of the two-point functions for an interacting scalar $O(N)$ model in de Sitter spacetime. Following our previous work, this behavior is analyzed by analytic continuation of the Euclidean correlators, which are computed by treating the homogeneous zero mode exactly and using a partial resummation of the interactions between the zero and the non-zero modes. We focus on massless fields and present an alternative derivation of our method, which involves a double expansion in $1/N$ and the coupling constant of the theory. This derivation is simpler than the previous one and can be directly extended for fields with negative squared-mass. We extend our previous results by computing the long wavelength limit of the two-point functions at next-to-leading order in $1/N$ and at leading order in the coupling constant, which involves a further resummation of Feynman diagrams (needed when the two-point functions are analitically continued). We prove that, after this extra resummation, the two-point functions vanish in the long distance limit.


Introduction
Since its formulation, quantum field theory (QFT) in de Sitter (dS) spacetime has been a subject that has received a lot of interest and attention. One clear reason is that in a dS spacetime, being a maximally symmetric curved space, it is possible to obtain exact solutions and it is relatively simpler to quantize fields. Physical motivations increased when realizing its application to the early universe, with the emergence of the inflationary paradigm, and also because of the discovery of the late time accelerated expansion of the Universe. This subject is clearly worth studying further.
The exponential expansion of the metric in dS spacetime produces an effective growth in the couplings of the theory, that turns out to be crucial in the infrared (IR). For instance, when considering a self-interacting scalar field of mass m in a λφ 4 theory, the effective coupling becomes λH 2 /m 2 , and therefore the usual perturbative calculations break down for light fields [1]. For massless free fields, minimally coupled scalar fields do not admit a dS invariant vacuum state [2].
The understanding of IR effects in QFT in dS spacetime is still incomplete. For massless interacting scalar fields in dS spacetime, this requires the use of nonperturbative techniques. There has been a lot of progress in understanding the leading IR effects. Indeed, there is more than one framework in which nonperturbative results have been obtained. For instance, the so-called dynamical mass of a quantum scalar field φ with classical potential λφ 4 /8 has been computed using the well known stochastic approach [3,4] and also by formulating the theory on Euclidean dS space (a 4-sphere) [5], yielding which is clearly nonperturbative in λ. As was originally emphasized in [6], using the analogy with finite-size effects in condensed matter, the generation of a dynamical mass represents an important and characteristic IR effect of the theory. Indeed, by splitting the field into an IR part, consisting of the near-homogeneous sector of the field, φ IR , plus the remaining ultraviolet (UV) part, φ = φ IR + φ UV , it can be shown that the dynamical mass is related to the curvature of the effective potential in the symmetric phase, and to the leading IR order (leading order in √ λ) of the correlator of the IR part of the field, (1. 2) The above nonperturbative result is important to cure the divergence appearing in the free massless correlator. However, corrections beyond this constant IR contribution are necessary to understand the long distance limit of the correlators. This last claim is of course irrelevant for the sphere, which is compact, but it is crucial for the theory in dS spacetime. For a theory with only one scalar field, a nonperturbative method that allows a systematic calculations of such corrections is very difficult and, to our knowledge, it is still lacking. One possibility recently explored involves the use of nonperturbative renormalization group techniques [7]. There are also nonperturbative calculations based on stochastic methods from which it is possible to obtain a numerical solution that predicts how the correlator of the IR sector φ IR (x)φ IR (x ′ ) decays when the dS-invariant distance between the two points, r(x, x ′ ), approaches infinity. To achieve this, the method is based on an unequal treatment of the IR and UV parts of the field. It is not clear however how robust this result is against the way in which the field is split and it is unknown how to compute corrections to the leading order result in a systematic way. The relation between the different approaches (stochastic, Lorentzian, Euclidean), has been the subject of many works [8][9][10][11][12].
The O(N )-symmetric model has been considered before in the large N limit [13][14][15][16][17][18][19][20]. In particular, in [17,18] we have presented a systematic nonperturbative resummation scheme in the case the number of fields N is sufficiently large so that it is possible to assume a double expansion in 1/N and √ λ. It is based on the theory on the sphere. It consists on a reorganization of the perturbative expansion and on an analytic continuation to the dS. The main advantages are: it improves the IR behavior of correlation functions in dS; it can be systematically improved, and the renormalization process is well understood.
In this paper we consider the same O(N )-symmetric model as in [17,18]. We start by summarizing (in Sec. 2) the main relevant properties of QFT in dS spacetime, and its analytical continuation to the sphere. In Sec. 3 we present a simpler derivation of the reorganized perturbative (in 1/N and λ 1/2 ) expansion, and we provide diagramatic rules (which we call resummed Feynman rules) to systematically account for corrections. In Sec. 4 we analyze the long wavelength limit of the correlators in the large-N limit up to next-to-leading order in 1/N . We show that, when analytically continued to dS spacetime, an additional nonperturbative resummation is necessary to take into account all relevant diagrams at long distances. After this resummation, the correlators tend to zero as r(x, x ′ ) → ∞. The next to leading IR contribution, which determines the decay law as r(x, x ′ ) → ∞, is evaluated using the resummation of Sec. 3 for all values of N . We present our conclusions in Sec. 5. The Appendices provide additional details of the calculations and the extension of the resummation scheme to negative-squared-mass fields.

QFT in Euclidean de Sitter space
The line element of Euclidean de Sitter space can be obtained from the one of dS spacetime in the so-called global coordinates, by the analytical continuation to imaginary time t → i(τ − π/2H). The periodicity condition τ = τ + 2πH −1 must be imposed in order to avoid a multivalued metric. This leads to the metric of a d-sphere of radius H −1 with θ = Hτ . We consider a O(N ) scalar model with quartic self-interaction living in this metric with Euclidean action given by where φ a are the components of an element of the adjoint representation of the O(N ) group, with a = 1, .., N . The sum over repeated indices is implied. Following the properties of the metric in Eq. (2.2), the field can be expanded in a discrete set of modes It can be easily seen that the free propagator in the symmetric phase is then expressed in the following way where the superscript m stands for the mass, is the surface area of the d-sphere, and r is the dS invariant distance, with w and w ′ unit vectors on the d − 1-sphere. In the last equality of Eq. (2.5) we explicitly separated the zero-mode ( L = 0) part of the propagator, which has an IR divergence for m → 0, thereby definingĜ (m) (r) as the inhomogeneous (i.e. L > 0) part, which is instead finite in that limit.

Analytical continuation
After analytic continuation back to dS, the dS invariant distance can be immediately computed in the standard cosmological patch. Using conformal coordinates, we obtain There are no ambiguities with this procedure for invariant functions depending only on two points x and x ′ through r(x, x ′ ). For massive fields, the dS free propagator decays at large distances as Moreover, it was also shown [21][22][23] that loop corrections to the two-point function of a massive field also enjoy such decay at large distances. In this way, the standard perturbative expansion is well defined for massive fields in dS, by analytically continuation from the sphere.
In the case of massless fields the story is different. Firstly, we must work with the modified propagatorĜ (0) (r), as defined in Eq. (2.5) by substracting the (divergent) zero-mode contribution to the standard propagator. Then, although not divergent for m → 0, it still suffers from IR effects once taken to dS in the form of a secular behavior at long distances. Indeed, using Eq. (2.10), .
This secular behavior slips into each loop correction, worsening as the number of loops increases [1]. Then, unlike the massive case, the standard perturbative expansion breaks down at large distances/late times in dS. This demands a nonperturbative treatment that not only gives a mass to the zero modes, but also to the inhomogeneous modes.

Reorganizing the perturbative expansion
In the context of the QFT on the sphere, let us first consider the approach introduced in Ref. [5,6], and later extended in Ref. [24], which addresses the divergence for m → 0 in the full propagator G by means of a reformulation of the theory in terms ofĜ. The main point of the approach is to split the field as φ a (x) = φ 0a +φ a (x), and treat the constant zero modes φ 0a nonperturbatively, separately from the inhomogeneous partsφ a (x). The interaction part of the action (2.3) is separated as follows: where S (2) int contains at least two powers ofφ a (note that the term linear inφ a vanishes identically by orthogonality). Then, since the zero modes are constant, the path integral over them turns into an ordinary integral, which can be performed exactly (i.e. nonperturbatively in the coupling constant λ). The generating functional becomes where J 0a andĴ a are external sources and we introduced the shorthand notation The zero part Z 0 [J 0 ] is defined as the exact generating functional of the theory with the zero modes alone, The simplest example of application is to compute the variance of the zero modes. For a massless field m = 0, the well known finite result is obtained, which allows the identification of a dynamical mass m 2 dyn by analogy with the free massive field. This result is valid at LO in √ λ and for all N . Corrections coming from the inhomogeneous modes can then be computed by treating S (2) int in Eq. (3.2) perturbatively. Although now the correlator is finite for m → 0, the perturbation theory built with a masslesŝ G (0) is still ill defined when analytically continued to dS and at long distances/late times, due to the divergent behavior described in Eq. (2.11). Solving this issue requires further resummations of contributions that also involve the inhomogeneous modes. A subclass of such contributions come from the terms in S (2) int that are quadratic in bothφ a and φ 0a , which dressĜ with a nonperturbative mass. This treatment was introduced by us in Ref. [17] for the O(N )-model (while in Ref. [25] it was considered for the special case N = 1). Here we will present a simpler formulation, which has the added benefit of being more easily generalizable to other cases, such as for negative squared-mass.

Resummation of bi-quadratic terms
In the spirit of the separation of the interaction part of the action done in Eq. (3.1), we further isolate the bi-quadratic terms, where now S int contains terms with at least three powers ofφ a . The main idea is then to include the bi-quadratic terms in the definition of the propagatorĜ. The generating functional becomes where the φ 0 -dependent inverse propagator of theφ a fields is given bŷ with the following mass matrix where in the second line we have split the matrix into the parallel and transverse components with respect to the ǫ a ≡ φ 0a /| φ 0 | direction, by means of the projector P ab = δ ab − ǫ a ǫ b , and we defined This tells us there are (N − 1) inhomogeneous fields with mass m 2 1 and a single one with mass m 2 2 . Diagonalizing the mass matrix m 2 ab , we can factorize the freeφ a part of the path integral (last factor of Eq. (3.6)) x,yφ (2) ≡ ǫ aφa are the transverse and parallel components with respect to the φ 0 -direction respectively, and now the indexes i, j run from 1 to N − 1. We also use the shorthand G α ≡Ĝ (mα) , with α = 1, 2. Putting this back into Eq. (3.6) we obtain i ,Ĵ (2) ], (3.12) where we are introducing a new "free" generating functional, whereẐ α ≡Ẑ f Ĵ (α) , m 2 α is the free generating functional of a single inhomogeneous field of mass m α , normalized toẐ α [0] = 1, and we have defined the notation choosing the normalization N such that Z[0] = 1. The superindex J 0 indicates that the0expectation value is taken over the zero modes in the presence of an external source J 0 . Notice that the functional derivatives with respect to the fullĴ a can be split in terms of derivatives with respect toĴ where we are taking advantage of the already manifest O(N − 1)-symmetry of the fields of mass m 1 , in the plane orthogonal to φ 0 . Afterwards, when applying these derivatives either toẐ 1 orẐ 2 , we will drop the superindices (1) and (2) as there will be no ambiguity on whichĴ is which. Later on, when setting J 0 = 0, the expectation value taken over the zero modes has the effect of averaging over the direction of φ 0 , giving rise to manifestly O(N )-symmetric expressions for the correlators.

Resummed Feynman rules
After having performed the above resummation of the bi-quadratic interactions, we now have a new "free" generating functional Eq. (3.13) where the propagatorsĜ α are massive, with masses that depend on | φ 0 |. If the remaining interaction terms in S (3) int are treated perturbatively, the new type of perturbative corrections can be computed using new Feynman rules, which we call Resummed Feynman rules, or R-Feynman rules for short. The R-Feynman diagrams are built with the dressed propagatorsĜ 1 andĜ 2 as internal lines. The last step is to evaluate the weighted average over the zero modes. The derivation of the new rules is straightforward, but with an important new ingredient that we will discuss below. First, let us write the interactions explicitly (according to the decomposition in the fieldsφ 1 andφ 2 ), (3. 16) In what concerns the Feynman rules, the | φ 0 | factor in front of the first two terms shall not be regarded as a "leg" of those vertices, but rather just as a coefficient 1 . The rules are summarized in Figs. 1 and 2.
x, i x ′ , ĵ i j Figure 2. Vertices. The factor A ijkl /3 = (δij δ kl + δ ik δ jl + δ il δ jk ) /3 in the first vertex accounts for the possible permutations of the indices.
There is an important difference with traditional Feynman rules to consider here. As a consequence of the definition of the new "free" generating functional in Eq. (3.13) as a weighted average over the zero modes, there is no direct cancellation of disconnected graphs when computing perturbative corrections. Indeed, consider a correction ∆Ẑ to the generating functional of the inhomogeneous modes, which at LO is justẐ 1Ẑ2 , prior to the zero mode average. The corrected complete 2 generating functional Z ′ now reads which corrects Eq. (3.13). Now consider a corrected n-point function of inhomogeneous fields 1 In the terminology of the standard perturbative expansion, the integration over the zero modes in . . . 0 takes care of summing over all possible ways of connecting the legs associated to φ 0a , both the ones that here appear explicitly in the vertices in Fig. 2 as well as those that are implicit in the masses of the resummed propagators m 2 1 ( φ 0 ) and m 2 2 ( φ 0 ). 2 It is complete in the sense that both zero and inhomogeneous modes are taken into account.
computed from the previous expression, where we used the normalizationẐ α [0] = 1, with α = 1, 2, and we also treated the correction ∆Ẑ perturbatively. The first term in the right-hand side is the leading contribution obtained from Eq. (3.13), while the second and third terms are the corrections. In the usual case, the second term contains both connected and disconnected contributions, the latter of which are cancelled by the third term. However, in the current situation this does not occur due to the weighting over the zero modes. Indeed, By adding and subtracting the left-hand side of the above equation to Eq. (3.18), we can identify two contributions to the correction of the n-point function as follows: a connected part which is built in the standard way with the connected R-Feynman diagrams using the above rules; and a 0-connected part which accounts for graphs that are disconnected according to the R-Feynman rules (that is, graphs that are not connected by lines associated to the propagatorsĜ α ), but when written in terms of the original perturbation theory they are actually connected by the lines that would correspond to the zero-modes. This means that in general, one has to include both contributions (3.20) and (3.21) when computing corrections to the correlators using this formalism. However, as we show in Appendix A, the 0-connected parts (3.21) are suppressed by extra powers of λ with respect to the connected parts (3.20). This will become clearer next, after considering an example in which we compute the 0-connected parts and verify they are of higher order in λ. For this reason, we will not need to compute the 0-connected parts for the computations we will consider here.

Resummed two-point functions
Let us apply the formalism just described to compute the two-point functions. We have to consider two distinct contributions, int , we only need to take derivatives of Z[ J 0 ,Ĵ (1) i ,Ĵ (2) ] in Eq. (3.13). For the constant part, i.e. the first term of Eq. (3.22), we take two derivatives with respect to J 0a and then set all external sources to zero. Exploiting also the O(N )-symmetry, it can be expressed as where now it is straightforward to perform the integration over the directions of φ 0 . Indeed, using which is now explicitly O(N )-symmetric. The remaining step is to compute the zero-mode expectation value. Prior to elaborating on how to do this, let us briefly discuss some perturbative corrections to these expressions.

Corrections
As a usage example of the R-Feynman rules, we compute the corrections coming from some of the interactions in S int . The simplest ones are local diagrams with at most one of the quartic vertices in Eq. (3.16).
We start by noticing that according to the R-Feynman rules we cannot build connected diagrams that correct φ 0a φ 0b , ∆ φ 0a φ 0b C ≡ 0. (3.27) There are, however, corrections in the form of 0-connected diagrams from the fact that there are bubble diagrams (i.e. corrections toẐ 1Ẑ2 ) that do not cancel in general. Indeed, where ∆Ẑ is given by the bubble diagrams shown in Fig. 3, Let us now consider corrections to the inhomogeneous part of the two-point function, Eq. (3.24). The needed connected R-Feynman diagrams are shown in Fig. 4. Notice that from Eq. (3.15), the diagrams with two external lines associated toφ 1 contribute proportional to P ab , while those where the external lines correspond toφ 2 go with a factor of ǫ a ǫ b . Therefore, using the R-Feynman rules the connected part reads where we are using the property Further simplification using properties (3.25) and dropping terms of order O(1/N 2 ) leads to Finally, for the 0-connected part, according to the definition (3.21) for n = 2, we have In order to evaluate the corrected two-point functions, we need to compute explicitly the integrals over the zero-modes that define the expectation values on the right hand sides of Eqs. (3.23), (3.26), (3.28), (3.32), and (3.33). This can be systematically performed order by order assuming a double expansion in λ and 1/N , as we show next.

Computation of zero-mode expectation values in a double λ and 1/N -expansion
In order to perform a systematic evaluation of the expectation values over the zero modes . . . 0, order by order in 1/N , we can use the saddle-point approximation (Laplace method). Indeed, by a simple change of variables u = λ| φ 0 | 2 /2N , the entire N dependence can be collected in an overall factor in the exponential, so that for a generic function of u, g(u), the zero-mode expectation values can be written as where we have introduced the following functions We explicitly emphasize the dependence of the propagatorsĜ 1 andĜ 2 with u through their masses m 2 1 = m 2 + u and m 2 2 = m 2 + 3u. Assuming an expansion for large values of N , as described in the Appendix A, we can approximate where a prime means a derivative with respect to u. The coefficients B (1) i (u), as well as the higher order ones, are built from the functions h(u) and D(u) and their derivatives, whileū is the solution where we have approximated [Ĝ 1 ] for smallū at linear order, assuming thatū ≪ H 2 . The coincidence limit of the UV propagator is divergent and shall be renormalized in the standard way, however we defer this discussion to the Appendix B. After the renormalization procedure, the expressions of these quantities turn out to be the same, but with the coincidence limit of the propagators and their derivatives replaced by the corresponding finite counterpart, and the constants m and λ by the renormalized quantities. From now on, we assume the replacements have been done, but for the sake of simplicity, we use the same notation for the finite quantities. The relevant (positive) solution is then Another property, shown in Appendix A, that will be useful for the computation of 0-connected parts is where g(u) and k(u) are arbitrary (although sufficiently smooth) functions of u , and all quantities on the right-hand side are evaluated atū. The coefficients C 2 (u) and C (2) i (u) depend of h(u), D(u) and their derivatives and are given in the Appendix A. Notice that in this case it is important to keep the N −2 terms due to the typical extra overall factor of N in the bubble diagrams (see for example Eq. (3.29)).

Massless case
We now focus on the case of a massless, minimally coupled field we are interested in. It can be readily seen that setting m = 0 it leads to an expansion in powers of √ λ. Indeed, expanding Eq. (3.40) up to order λ we havē Then, evaluating the coefficients of the expansions Eqs. (3.38) and (3.41) at thisū, allows them to be expanded in a similar fashion With all these ingredients we can now compute the two-point functions up to next-to-next-to LO (NNLO) in √ λ and next-to LO (NLO) in 1/N . At such order, it is enough to evaluate the sum of the contributions in Eqs. (3.23) and (3.26), and the connected correction to the inhomogeneous part Eq. (3.32) (recall there is no such correction for the constant part). For completion, in what follows, we also show separately the 0-connected corrections, both the constant Eq. (3.28) and the inhomogeneous Eq. (3.33) corrections, and we verify that they are of higher order in √ λ with respect to the parts we keep.
The results are the following. For the constant part while for the inhomogeneous part while for the inhomogeneous part In both cases, with these explicit expressions we can now confirm that the 0-connected contributions are suppressed by extra factors of λ with respect to Eqs. (3.44) and (3.45).
It is interesting to analyze the large distance behavior of resummed two-point functions obtained after our procedure. First, notice that in the final expression for the inhomogeneous part, Eq. (3.45), the free inhomogeneous propagatorsĜ 1 (r) andĜ 2 (r) are now evaluated for positive masses M 2 1 ≡ m 2 1 (ū) and M 2 2 ≡ m 2 2 (ū), Therefore, upon analytic continuation to dS, both approach a constant exponentially at large distances, as discussed in subsection 2.1 for massive fields (Eq. (2.10) minus a constant part), rather than the logarithmic divergence of massless fields (as in Eq. (2.11)), that plagued the perturbative expansion prior to our resummation. However, notice that the result Eq. (3.45) also contains derivatives ofĜ 1 (r) andĜ 2 (r) with respect to their masses, and moreover, the number of derivatives increases with the orders in the double √ λ and 1/N -expansion. This is due to the treatment given to the zero-mode weighted integrals (see Eq. (3.38)). The behavior of these derivatives in dS and at large distances is easily 4 Notice that in Ref. [17] a further approximation was assumed to be valid, namely λ ∂Ĝ ( . This is correct on the sphere, but it actually fails in dS at sufficiently long distances. obtained from Eq. (2.10) and found to be (3.49) meaning that, although there is an overall exponential decay, for large enough distances higher order terms become as relevant as the lower order ones, thus breaking the expansion. The origin of this problem is in the assumption used when computing (3.45) in the double expansion in √ λ and 1/N , which requires √ λ log(r) ≪ N . A proper analysis of the long distance behavior must then consider the opposite hierarchy √ λ log(r) ≫ N , and will be discussed in the next section. Nevertheless, this does not invalidate the result (3.45), as long as one is only interested in points not too far apart (or in the Euclidean correlators).
Let us close this section by pointing out a related aspect to be considered when using the results on the sphere to study the large distance (IR) behavior by analytical continuation to dS, which is the use of power counting rules to select the relevant diagrams. Notice that the power counting in the coupling constant λ we use here 5 assumes there is no additional nonperturbative IR enhancement. This is certainly the case for the sphere, but is not when the correlators are analytically continued to dS. Indeed, due to the non compactness of the dS spacetime, it turns out that the counting used here applies only up to a certain maximum distance, but not beyond. This demands going beyond perturbative R-Feynman diagrams to understand the far IR limit by collecting all contributions at a given order in λ that are relevant in this limit. In order to make this feasible, an organizing principle such as powers of 1/N is necessary. We will discuss more on this in the following section.

Leading IR contribution: a consistent calculation beyond N → +∞
Let us first discuss the limiting value of the two-point functions, once analytically continued back to dS, for large distances/late times, r → +∞. Contributions to this limit may come either from the zero-mode part, or from the limiting value of the inhomogeneous part. Indeed, according to Eqs. (2.5) and (2.10), in this limit the massive UV propagatorsĜ α go tô with m 2 2 = 3m 2 1 = 3u. However, once all the constant contributions are put together, one expects the full two-point function to approach a vanishing value.
After the resummation of the bi-quadratic terms, the full two-point function without any corrections from S Combining the constant contributions, we have where we used Laplace's method at NLO in 1/N to compute the integrals 6 .
First observe that this result vanishes at LO in the large N limit, N → +∞. This is because the family of diagrams that contribute at large-N is that of the daisy and superdaisy type, which add a local part to the self-energy, and whose leading contribution in √ λ is already completely taken into account by the exact treatment of the zero modes.
Starting at NLO in 1/N we encounter a nonvanishing value. The reason is that at NLO in 1/N there are diagrams that contribute at LO in √ λ that are not of the type included in the resummation above. Indeed, the relevant diagrams that must be added to the self-energy are non local and have the form of a bubble-chain. These are only partially accounted for in our treatment so far, even at the LO in √ λ. This stems from the fact that, although on the sphere there is a hierarchy between interactions with φ 0 and φ in terms of powers of √ λ, upon analytical continuation to dS, when r → ∞ the correlators of inhomogeneous modes receive an enhancement that makes them as relevant as the zero mode correlators, Therefore, a consistent calculation of the two-point function at large distances at in dS, at LO in √ λ and NLO in 1/N , requires the addition of diagrams correcting Eq. (3.24). To show this explicitly, let us now evaluate the limiting value of the diagrams when the distance between the two points goes to infinity. Doing this, we will get a non-vanishing result that should be added to Eq. (4.3).
In order to identify the relevant diagrams that potentially contribute to the limiting value at LO in √ λ we use the scaling in Eq. (4.4) as a new way of counting the corresponding order in √ λ. An strategy to perform the calculation of the diagrams consists in the following two steps: I) use the R-Feynman rules described in Sec. 3.3 to compute the corrections to the propagatorŝ G 1 andĜ 2 (with | φ 0 |-dependent masses) and replace the corrected propagators into Eq. (3.24); II) evaluate the two point functions using (3.25), as for Eq. (3.26), and taking the expectation value over | φ 0 | assuming the double expansion in √ λ and 1/N . Recall that the 0-connected diagrams are suppressed by extra factors of √ λ. From Eq. (3.26) it is immediate to conclude that forĜ 1 it is necessary to compute the NLO correction in 1/N , while forĜ 2 the LO is enough. We illustrate the calculation with the diagram in Fig. 5a, that corrects the propagatorĜ 2 . To begin with, it is easy to see that the diagram does not vanish at N → ∞.
Our method to obtain the contribution of a given diagram in the large-distance limit is based on general theorems proved in Refs. [21,22]. The theorems imply that a calculation of a two-point Feynman diagram using full massive propagators G α (i.e. the standard ones) instead ofĜ α (both inside the loops and in the external legs) would give a vanishing result at large distances. The R-Feynman rules involve the modified propagatorŝ that differ from the standard massive propagators by the constant contributions of the Euclidean zero modes. Therefore we can make the substitution above into the R-Feynman diagrams and keep, in the large distance limit, only the contributions that are not connected by G α . This is so because after making the substitution (4.5), any two-point diagram connected by full massive propagators G α , up to constant factors (possibly involving factors of G (0) α ), can be thought as a particular two-point Feynman diagram which are known to vanish at large distances [21,22].
There is an infinite number of diagrams that contribute to this order, that can be obtained by dressing the bubble in Fig. 5a. The dressed diagram in Fig. 5a ′ contains the chains of bubbles given in Fig. 6. Fortunately, this set of diagrams can be computed and resummed, following the same method applied for the one with a single bubble, that is, writing the propagatorĜ α in terms  At higher orders, there are even more additional diagrams contributing to the constant limiting value. This can be seen from the fact that the exact treatment of the zero modes contains contributions at all orders in 1/N , which however are not compensated by the limiting value of the diagrams included so far for the inhomogeneous part, even at the LO in √ λ. Indeed, while the diagrams in Fig. 6 are enough at NLO in 1/N , using Eq.(4.4) one can see that a diagram of the type shown in Fig. 8, which is NNLO in 1/N , should also be taken into account at LO in √ λ. Despite the added complexity to perform a consistent calculation at large distances beyond N → +∞, the situation now is nevertheless much better than before resumming the bi-quadratic terms, where the correlators of the inhomogeneous modes were massless and outright divergent in the IR (see Eq. (2.11)).

Asymptotic IR behaviour of the two-point functions
The next step is to study the decay of the full two-point function at large distances. To do this properly at a given order in 1/N , again it would be necessary to include the most relevant rdependent parts of all the diagrams that contribute at that order. We intend to pursue this in a systematic way in future work. Our aim here is to just analyze what kind of decay is expected after resumming the bi-quadratic terms and the exact treatment of the zero modes.
Here, being interested in the asymptotic behavior at large distances and/or late times, we work out a different approach to compute the r-dependent part of the two-point function in the r → +∞ limit for all N . For this, we will assume that this limit can be taken before integrating over the zero modes.
For r → +∞, we can approximate the analytically-continued (to dS) propagators Here we have defined y = log(r)/(d − 1). Looking only at the r-dependent part, having already discussed the constant part, the two-point functions at long distances are given by the following expression In the second line we have expanded the propagators for small u up to quadratic order. Here we can see that the integrands of Eq. (4.7) are exponentially suppressed for u ≫ λ/2V d , and therefore, for small enough λ the approximation of Eqs. (4.6) is justified. Notice that for the cases N ≤ 2 there is a IR divergence in the integral of Eq. (4.7). This divergence is spurious and is actually absent when one considers both the r-dependent and constant parts together. The needed integrals are of the form where U (a, b, z) is the Tricomi confluent hypergeometric function, and A and B are the linear and quadratic in u coefficients ofh(u) respectively, defined in Eq.(4.9). With these expressions, we can evaluate Eq. (4.7), where we defined U t (x) = U N +t as a shorthand. Although we have assumed a large value of y when using Eqs. (4.6), nothing has been said about N . For N → +∞, the limit must be taken carefully and before y → +∞, since the parameter N enters in several arguments of the hypergeometric functions. Doing this, as crosscheck, we recover the known result for the asymptotic behavior at large distances of the two-point function, with m 2 dyn = λ/(2V d ). We are now in position to study the y → +∞ asymptotic behavior for fixed N . For this we can use that U 0 (x) ≃ B N/4 x −N/2 at large x, and therefore the resummed UV two-point function takes the following form, where we kept the most relevant term for y → +∞. The coefficient is defined as with their LO dependence in λ explicitly shown. The behavior of the resummed two-point functions at large distances changes drastically with N , as seen from Eq. (4.13). On the one hand, for N > 2 they decay as a power law for y → +∞. On the other hand, for the cases N ≤ 2, as already mentioned, both the r-dependent and constant parts must be kept during the computation. This allows to show that for N = 1, 2 the two-point function diverges as √ y and log(y) respectively 7 . These distinct asymptotic behaviors for different values of N are summarized in Table 1, where also the standard perturbative results (prior to the resummation of bi-quadratic interactions) are shown for comparison.
We now have a more complete picture of the asymptotic behavior at large distances/late times of the resummed two-point functions. Indeed, for large but not strictly infinite N , the behavior can initialy be well described by Eq. (3.45) as quasi-exponential. But as y approaches and exceeds ∼ N/ √ λ, the behavior changes to a power-law of log(r), as shown in this section, with an exponent strongly dependent on the value of N .
In summary, in this Section we have analyzed the effects of the resummation on the IR behavior of the two-point functions. Given that the bi-quadratic interactions between the zero and UV modes are treated nonperturbatively, the UV modes acquire a | φ 0 |-dependent mass, and therefore the IR behavior is affected. It is convergent or less divergent, depending on the value of N . To understand these results, we go back to Eq. (3.26), that shows that the interacting two-point function is a weighted average of free and massive two-point functions. As the free correlators decay exponentially, one could naively expect the interacting correlator to have the same behavior. However, for massless fields m = 0 the bi-quadratic interaction induces a mass for the UV modes that is proportional to u = λ| φ 0 | 2 /2N and vanishes at u = 0. Therefore, the IR behavior is dominated by the contributions coming from the region u ≪ H 2 in the integral of Eq. (4.7). Due to the volume factor u N/2−1 in that equation, these contributions are suppressed for large values of N , and the correlators decay more strongly as N increases.
It is unclear whether additional corrections coming from the remaining (infinitely many) R-Feynman diagrams that contribute at each order in 1/N , when treated nonperturbatively, might change this behavior. This can be particularly important for small N , for which the decay given by the current resummation is milder, or absent altogether. For instance, if the UV modes were to have an additional dynamical mass coming from the interactions in S (3) int , this mechanism could dominate the IR behavior.

Conclusions
In this paper we have considered the same O(N )-symmetric theory and the same reorganization of the perturbation theory on the sphere (in λ and 1/N ), as in [17]. We have presented a procedure that allows the systematic perturbative expansion we developed in [17] for massless fields to be performed in a much simpler way. Another advantage of the present procedure is that it can be immediately extended to fields with negative mass squared, as we show in Appendix C.
The results presented here also extend the previous ones, providing an explicit evaluation of the long wavelength limit of the two-point functions for arbitrary values of the number of fields N . The resummed two-point function is given in Eq. (3.26), that shows that it can be written as a weighted average of free, massive two-point functions. This result can be interpreted as a spectral representation of the two-point correlation function of an interacting field, analogous to the Källen-Lehmann representation in Minkowski spacetime [25]. Being an average of functions that decay exponentially at large distances, the IR behaviour of the resummed two-point function is improved: it tends to a constant as r → ∞ for all N > 2. The aymptotic r-dependence summarized in Table  1 shows how this constant value is approached. These results should be contrasted with the usual divergent results obtained in perturbation theory for interacting massless fields. In the IR, the weighted average is dominated by the contributions of free two-point functions of small masses, that become more relevant for low values of N . Because of this reason, for N = 1, 2 the resummed two-point functions do not tend to a constant as r → ∞, but diverge with a milder divergence than in perturbation theory.
The reorganization of the perturbation theory, implemented by treating exactly the bi-quadratic interaction terms between the constant zero modes and the inhomogeneous fields, is not enough to describe the far IR behaviour of the correlation functions. The corrections computed within a combined double expansion in 1/N and √ λ show that higher order terms become as relevant as lower order ones when the correlation functions are analytically continued to dS spacetime, and evaluated at very long distances √ λ log(r) N (note that this is not a problem when computing the correlator functions on the sphere, which is compact). We have shown that in order to analyze the behaviour of the two-point functions for such large distances, an additional resummation is mandatory, since an infinite number of diagrams that correct the inhomogeneous propagators contribute at order O √ λ, 1/N in the double expansion. We have identified these diagrams, broadly speaking, as those containing chains of an arbitrary number of bubbles. The fact that the R-Feynman rules needed in the additional resummation involve the free massive propagators that come from the first resummation, allowed us to use general results about the IR behaviour of interacting massive propagators in dS spacetime, and to perform explicitly the resummation of the bubble chains in the large-r limit. Within the double expansion in 1/N and √ λ we have shown that the contribution of the bubble chains changes the behaviour of the two-point functions, that now tend to zero (instead of a constant) as r → ∞ up to O √ λ, 1/N . We expect the resummation of bubble chains to modify the decays in Table 1 for large values of N as well, when the subleading IR parts are included. When the number of fields is small, the situation is less clear, since even more diagrams would be needed to understand the large-r limit. We consider this is out of the scope of the present paper, but we hope to address this issue in the future.
Finally, let us comment about the applicability of our method to physical situations that involve n-point functions with n > 2. This is in principle non trivial, due to the fact that the method directly applies to the sphere rather than to dS. Although it is immediate to obtain the different two-point functions in dS spacetime given the unique two-point function on the sphere (which is well-known [26] and is completely analogous to the same problem in Minkowski space), it is unclear how to proceed in general for n-point functions. Notice however that this problem arises due to causality, there is no such a problem as far as one is interested in computing n-point functions for spatiallyseparated points. Therefore, the calculation of n-point functions using our resummation could have several cosmological applications. We plan to consider this problem, starting with the calculation analogous to that in Ref. [27] for comparison of the methods.

Acknowledgments
This work has been supported by CONICET , ANPCyT and UNCuyo. We would like to thank Bei-Lok Hu and Stefan Hollands for useful discussions on related matters. The diagrams were done with the tikz-feynman package [28]. which is approximated up to NLO in 1/N by the following formula where it is understood that all the functions are evaluated in the saddle-pointū, i.e. h ′ (ū) = 0, h ′′ (ū) < 0. If there were more than one point that satisfies these conditions, the solution would involve a sum over them of the previous expression. We can now apply this result to expressions of the kind by choosing f (u) = D(u)g(u) and expanding the denominator in 1/N as well. This last step will ensure the cancelation of any NLO terms without derivatives of g(u). The resulting expression is Another useful combination that appears when computing perturbative corrections is the following, For this last expression it was necessary to extend Eq. (A.2) up to NNLO in 1/N . In the massless case,ū ∼ √ λ, and the coefficients up to order λ are expanded as As it turns out, all of the C coefficients start at least at order λ. When the combination (A.6) comes from a perturbative calculation, there will be another λ in front, meaning that up to higher order terms, it is valid to factorize a mean value over the zero modes of a product, in terms of products of mean values, i.e. g k 0 ≃ g 0 k 0, even beyond the large-N limit. This greatly simplifies calculations of perturbative corrections. Indeed, the application of this property to the 0-connected parts defined in Eq. (3.21), implies that they are suppressed by extra powers of √ λ with respect to the connected parts (3.20).

B Renormalization
The functions h(u) andh(u), defined in Eqs. (3.36) and (4.8) respectively, contain divergent contributions from the coincidence limit of the UV propagator [Ĝ (m) ] and its first derivative (second and higher order derivatives are finite in this limit). These can be absorbed in the usual way by the introduction of counterterms associated with each of the parameters m 2 and λ. Consider for instance the functionh(u) for an arbitrary mass m 2 , where the labels B have been added to denote the bare quantities. Splitting these into a finite part and a conterterm, λ B = λ + δλ and m 2 B = m 2 + δm 2 , it is straightforward to see that the following choice of counterterms allows for perturbative renormalization, Although the first term is still divergent, it actually drops from any physical quantity due to the normalization of the generating functional. For the sake of clarity, we avoid the ren and f in labels in the main text.

C Negative squared-mass case
In this Appendix we consider an extension of the resummation to cases with negative-squared-mass fields, and we show that our method can be applied without any substantial modification. Indeed, we can directly evaluate Eq. (3.40) for a negative squared-mass m 2 = −µ 2 < 0, and then expand in λū Then, the constant part of the two-point function (3.23) now reads, at the leading order in λ, which is a well known result. The positive value of the dynamical mass shows that the symmetry is restored in the effective potential due to the IR effects. For the two-point function at separate points, we find that the effective masses of the propagators are It is worth remarking that both masses are strictly positive, showing that this resummation is enough to overcome the divergences appearing for the tree-level propagators. When analytically continuing the results to dS, again the same discussion as for massless fields applies, that is, the results can only be trusted for points not too far apart. Let us briefly comment on the perturbative corrections. The main aspect to take into account is that now the effective coupling is different than in the massless case. From Eq. (C.2) we have the scaling | φ 0 | 2 0 ∼ λ −1/2 , and therefore the IR enhancement is stronger than for the massless case ( | φ 0 | 2 0 ∼ λ −1/4 ). Hence, all the interaction terms in S int contribute about the same. Indeed, taking into account that the terms of the form λ( φ 0 · φ )| φ | 2 must always be considered in pairs (by symmetry considerations, correlators with an odd number of factors of φ 0a will vanish), their contribution will go as λ 2 | φ 0 | 2 0 ∼ λ. Given that the effective coupling in this case is also λ, the perturbative corrections contribute already at the NLO level. Of course, a calculation of all contributions at that order will be much more involved than for the massless case, but, in principle, it can be performed using the R-Feynman rules.
In this Appendix we provide the details of the calculation of the long distance limit, r(x, x ′ ) → ∞, of the two-point functions at leading order in √ λ and next-to-leading (NLO) order in 1/N . As mentioned in Sec. 4.1, there is an infinite set of Feynman diagrams contributing to the two-point functions in this limit.
We need to compute all Feynman diagrams which contribute at NLO in 1/N and are in turn enhanced at long distances, in the sense that they give a correction at LO in √ λ in the limit r(x, x ′ ) → ∞. The first diagram given in Fig. 4 is an example that contribute at NLO in 1/N (and also at LO), but that is clearly NLO in √ λ, since the propagator in the loop is evaluated at coincident points and therefore cannot give an enhancement for r(x, x ′ ) → ∞. Hence, note the diagrams we need contain at least two vertexes. To proceed, we use the R-Feynman rules described in Sec. 3.3 to compute corrections to both propagatorsĜ 1 andĜ 2 whose masses depend on | φ 0 |. Then, the two point functions can be evaluated using Eq. (3.25). Hence, forĜ 1 it is necessary to compute the NLO correction in 1/N , while forĜ 2 the LO is enough. To count the power of 1/N at leading order in the limit N → ∞, it is useful to note that | φ 0 | must be considered as counting as √ N . In other words | φ 0 | = 2N u/λ and u, as it turns into an integration variable, is considered to be independent of N . As usual, each trace with solid lines gives a factor of N .
Let us first consider the diagrams correctingĜ 2 . Taking the above considerations into account, it is simple to see the diagram in Fig. 5a does not vanish at N → ∞. There are two ways of adding corrections to this diagrams that can be thought as dressing the bubble in the middle in two steps: First, the bubble can be corrected by adding another bubble with a vertex with four solid lines (this gives a trace and therefore a factor N that compensate the one suppressing the vertex). A sum of the chain of bubbles corrected in this way gives a partially-dressed bubble, which we draw as a solid gray circle, = + + + . . . ≡B , (D.1) and we will denote asB. Second, each partially-dressed bubble can be further dressed by adding two vertexes with two solid lines and a discontinuous one, as follows: The dressed bubble (which is plotted as a black solid circle and we will denote asB) is obtained by summing this chain of diagrams. Therefore, the diagram we need to compute can be represented as in Fig. 5a ′ .
We now consider the diagrams correctingĜ 1 . We start by writing all the different diagrams that contribute at order 1/N and have the minimum amount of vertexes. These are given in the first line of Fig. 9. Then, as forĜ 2 , there is actually a set of diagrams contributing at the same order, and the sum of them can be expressed in terms of the same resummed bubbles as in the second line of Fig. 9. Of course, a brute-force evaluation of the diagrams we have drawn is too difficult, however, our aim here is to compute the leading order in √ λ in the long distance limit. To achieve this, we will explote the knowledge of perturbation theory for massive fields [21,22]. The main observation is that the R-Feynman rules are built from modified "massive" propagatorsĜ 1 andĜ 2 , which differ from standard massive propagators only by an homogeneous (constant) contribution (see Fig. 10 for the corresponding diagrammatic representation): where G α is the standard propagator with mass m a and (D.4) With the use of Eq. (D.3), the two-point diagrams we need to compute can be split into diagrams where the two points are connected by at least one massive propagator G α and disconnected diagrams. It has been shown that any connected two-point diagram built with massive propagators and standard cubic or quadratic vertexes vanishes at long distances [21,22]. Therefore, we need to keep and compute only the two-point diagrams that once split according to Eq. (D.3) have the two points disconnected. All of these facts will become clearer in a moment, after considering a few examples. Let us start with the diagram in Fig. 5a ′ . Notice that the bubble itself,B(z, z ′ ), can be split into a connected bubble, B(z, z ′ ) (represented as a black square), plus an homogeneous disconnected part, B d (drawn as a cut bubble): Then, the diagram in Fig. 5a ′ can be split as where in the last line we discarded the connected contributions, which vanish at long distances. In formulas: In order to obtain the final result we need to perform the integration in u, which in the 1/N expansion can be done systematically using the Laplace method, as described in Appendix A. The result cancels the constant in Eq. (4.3) in the main text, and makes the two-point functions to vanish in the long distance limit.
We finally mention here that the inhomogeneous 0-connected contribution (3.47), which we dropped in the calculation on the sphere, will also be enhanced in dS. This enhancement is stronger the more derivatives with respect to m 2 act on the propagatorsĜ 1 andĜ 2 , and therefore it will be more important for the 0-connected part than for the part we kept (3.45). This is why we kept those terms with the most derivatives acting on the r-dependent functions in (3.47). However, it is straightforward to check that even with the stronger enhancement, the 0-connected part is still at least suppressed by one factor of √ λ with respect to the other parts. Moreover, this is true beyond the specific diagrams computed in previous section, and in particular it remains true for the infinite set of diagrams that we need to resum here in order to deal with the limiting value for r → ∞. This follows directly from counting powers of λ in the expansion (3.41) with coefficients (3.43), and taking into account that only one of the two functions g(u) and k(u) depends on r and therefore receives an enhancement in dS of a factor λ −1/2 for each derivative acting on it. Therefore, we can still ignore the 0-connected parts in dS in the current discussion, as it does not affect the limiting value of the two-point functions for r → ∞, at least up to NNLO in √ λ.