e1e-mail<EMAIL_ADDRESS>e2e-mail<EMAIL_ADDRESS>e3e-mail:
<EMAIL_ADDRESS>11institutetext: Dipartimento di Fisica e
Astronomia, Università di Bologna and INFN, via Irnerio 46, 40126 Bologna,
Italy 22institutetext: L. D. Landau Institute for Theoretical Physics, Moscow,
119334 Russia 33institutetext: Bogolyubov Laboratory of Theoretical Physics,
Joint Institute for Nuclear Research, 141980 Dubna, Russia 44institutetext:
Dipartimento di Fisica e Chimica, Università di L’Aquila, 67100 Coppito,
L’Aquila, Italy
and INFN, Laboratori Nazionali del Gran Sasso, 67010 Assergi, L’Aquila, Italy

# Massive scalar field in de Sitter spacetime:
a two-loop calculation and a comparison with the stochastic approach

Alexander Yu. Kamenshchike1,addr1,addr2  Alexei A. Starobinskye2,addr2,addr3
Tereza Vardanyane3,addr4

(Received: date / Accepted: date)

###### Abstract

We examine long-wavelength correlation functions of massive scalar fields in
de Sitter spacetime. For the theory with a quartic self-interaction, the two-
point function is calculated up to two loops. Comparing our results with the
Hartree-Fock approximation and with the stochastic approach shows that the
former resums only the cactus type diagrams, whereas the latter contains the
sunset diagram as well and produces the correct result. We also demonstrate
that the long-wavelength expectation value of the commutator of two fields is
equal to zero both for spacelike and timelike separated points.

††journal: Eur. Phys. J. C

##  1 Introduction

Quantum field theory in expanding spacetimes is crucial for exploring the
origin and the evolution of the universe. Scalar fields in de Sitter spacetime
are of particular importance for understanding the period of inflation and the
growth of quantum fluctuations during this period.

Free scalar fields on the de Sitter background—both massive and massless—have
been studied extensively (see, e.g., ChT ; Tagirov ; BD ; Linde:1982uu ;
Starobinsky:1982ee ; Vilenkin:1982wt ; Allen:1985ux ; Allen:1987tz ). It has
been shown that for massive fields there exists a one-parameter family of de
Sitter-invariant states; this includes the Bunch-Davies vacuum, whose modes in
the Poincare patch match with the adiabatic in-vacuum modes in the remote
past. From this family, only the Bunch-Davies vacuum has finite renormalized
values of $\langle\phi^{2}\rangle$ and the energy-momentum tensor of the
scalar field.

The massless minimally coupled scalar field also has the Bunch-Davies state,
but it is no longer de Sitter invariant. It is, however, invariant under the
$\rm E(3)$ subgroup, which is the group of translations and rotations in the
flat spatial sections of de Sitter spacetime. One of the manifestations of the
breakdown of de Sitter invariance is the presence of an infrared divergence in
the theory. Hence, to evaluate correlation functions, an infrared cutoff
should be introduced in momentum space. As a result, expectation values of
products of fields are explicitly time-dependent; for example,

$\langle\phi^{2}(\vec{x},t)\rangle_{\rm ren}={H^{3}t\over 4\pi^{2}}+{\rm
Constant}\,,$

where $H$ is the constant Hubble parameter in the Poincare patch of de Sitter
space (see (3) below) and $t=0$ denotes the beginning of de Sitter stage
Linde:1982uu ; Starobinsky:1982ee ; Vilenkin:1982wt  (see also Allen:1985ux ),
e.g., arising from a previous generic anisotropic and inhomogeneous stage with
curvature much exceeding $H^{2}$ that does not require fine tuning of pre-
inflationary initial conditions. At late times, the first term on the right
side dominates. This time-dependent term arises solely from long-wavelength
modes, that is, modes with physical wavelengths much greater than the Hubble
scale.

A similar result is obtained for small gauge-invariant metric tensor
perturbations— gravitational waves (GW)—of the de Sitter background
Starobinsky:1979ty :

$\langle h_{ij}h^{ij}\rangle_{\rm ren}=\frac{16GH^{3}t}{\pi}+{\rm
Constant}\,,$

where $G$ is Newton’s gravitational constant and $i,j=1,2,3$ (in
Starobinsky:1979ty  this is calculated in the synchronous gauge, but it does
not matter, since the tensor perturbations $h_{ij}$ are gauge invariant). In
the context of the inflationary scenario where the duration of de Sitter
(inflationary) stage is finite, i.e. our past light cone is finite, this
formula just corresponds to the standard expression for the primordial power
spectrum of tensor perturbations generated during inflation, summed over
polarizations: $P_{g}=\frac{16GH^{2}}{\pi}$, where slowly varying $H$ is
estimated at the moment of the first Hubble radius crossing, $k/a(t)=H(t)$,
for each mode with a comoving momentum k. In Starobinsky:1979ty , these
results were explicitly presented in another equivalent form—as the energy
density of primordial GW after the second Hubble radius crossing at the
subsequent radiation dominated stage.

In the massive case there is no infrared divergence, but if $0<m^{2}\ll
H^{2}$, the leading contribution to $\langle\phi^{2}\rangle_{\rm ren}$, which
is given by Linde:1982uu ; Starobinsky:1982ee ; Vilenkin:1982wt ; Allen:1987tz

$\langle\phi^{2}(\vec{x},t)\rangle_{\rm ren}={3H^{4}\over
8\pi^{2}m^{2}}+\mathcal{O}\left(\left(m^{2}/H^{2}\right)^{0}\right),$

still derives entirely from the long-wavelength modes.

When there is a self-interaction, each successive term in the weak coupling
perturbative expansion contains higher and higher powers of $Ht$ or
$H^{2}/m^{2}$, depending on whether the theory is massless or massive. This
means that perturbation theory breaks down when the value of $Ht$ or
$H^{2}/m^{2}$ overwhelms the smallness of the coupling constant Woodard ;
Prokopec:2007ak ; we . A non-perturbative method for calculating the
expectation values of the coarse-grained theory, containing only the long-
wavelength fluctuations of a scalar field, was developed by Starobinsky in
Star  and further developed in many papers, in particular in Star-Yoko ;
Vennin:2015hra . An important result of Starobinsky’s approach is that
expectation values can be determined by using a probability distribution
function that is a solution to a simple Fokker-Planck equation. The emergence
of the Fokker-Planck equation from the full quantum evolution of the theory
was presented in more recent works Tereza ; Burgess:2014eoa ; Burgess:2015ajz
; Tokuda:2017fdh ; Pinol:2020cdp ; Baum-Sun ; Gorb-Sen .

Due to the non-perturbative nature of Starobinsky’s stochastic approach, it is
important to check it by comparing its results with those obtained from the
standard field-theoretic perturbative approach whenever the latter is
possible. This was done in the one-loop (Hartree-Fock, or Gaussian)
approximation in Finelli:2008zg ; Finelli:2010sh . Here we consider the two-
loop case in the regime when the scalar field mass and its quartic interaction
are sufficiently small. A number of two-loop calculations for scalar fields
with quartic self-interaction in de Sitter spacetime were made using the
stochastic approach and different quantum field theoretical methods
Marolf:2010nz ; Gat-Ser ; Hollands:2010pr ; Akhmedov:2013xka ; Garb-Rig-Zhu ;
Gautier:2015pca ; Onemli ; Onemli1 ; Mark-et ; Kitamoto ; Mor-Ser ; Mark-Raj .

In our previous paper we , we considered a massless minimally coupled scalar
field with a quartic self-interaction. We calculated the long-wavelength part
of $\langle\phi^{2}(\vec{x},t)\rangle$ and $\langle\phi^{4}(\vec{x},t)\rangle$
to second order in the coupling constant and developed a method for taking the
late-time limit of these perturbative series. This method was inspired by the
renormalization group and is based on the assumption that expectation values
of products of fields at coinciding spacetime points satisfy autonomous first-
order differential equations. We then compared our results with the results of
the stochastic approach and found them to be in good agreement.

In this paper we turn to the study of the massive case. In Sec. 2 we review
the quantization of a free massive scalar field in de Sitter space. Section 3
presents a perturbative calculation of the long-wavelength part of the two-
point function. Section 3.1 contains the zeroth-order computation and a
comparison with the untruncated two-point function of the free theory. In Sec.
3.2 we outline the “in-in” formalism and use it to calculate the one- and two-
loop corrections in Sections 3.3 and 3.4. For coinciding spacetime points our
perturbative expression agrees with those obtained in Gat-Ser ; Garb-Rig-Zhu ;
Mor-Ser . We compare our results with the Hartree-Fock approximation in Sec. 4
and with the stochastic approach in Sec. 5. In Sec. 6 we argue that the
perturbative expression for the two-point function can be reorganized into a
sum of exponential functions that depend on the two given points in a de
Sitter invariant way.

##  2 Scalar field in de Sitter spacetime

We will study a massive scalar field with a quartic self-interaction,

$S=\int
d^{4}x\sqrt{-g}\left(\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi-{1\over
2}{m^{2}}\phi^{2}-{\lambda\over 4}\phi^{4}\right)\;,$ (1)

propagating in the de Sitter background represented as an expanding spatially
flat Friedmann universe (the Poincare patch) with the metric

$ds^{2}=dt^{2}-a^{2}(t)\delta_{ij}dx^{i}dx^{j}\;,$ (2)

where the scale factor $a(t)$ is

$a(t)=e^{Ht}\;.$ (3)

Here $t$ is a cosmic time coordinate and $H$ is the Hubble constant or the
inverse of the de Sitter radius. The cosmic time in an expanding de Sitter
universe runs in the interval $-\infty<t<\infty$. It will also be convenient
to use a conformal time coordinate $\eta$, which is related to the cosmic time
$t$ by the condition $dt=a(\eta)d\eta$. Expressed in terms of the conformal
time, the metric is

$ds^{2}=a^{2}(\eta)(d\eta^{2}-\delta_{ij}dx^{i}dx^{j})\;,$ (4)

where

$a(\eta)=-\frac{1}{H\eta}\;,$ (5)

and $\eta$ runs from $-\infty$ to $0$.

We will assume that the coupling constant and the mass are small: $\lambda\ll
1$, $m^{2}\ll H^{2}$. However, as we will see from the perturbative loop
expansion, it has to be assumed additionally that the mass is not too small:
$\sqrt{\lambda}H^{2}\ll m^{2}\ll H^{2}$. Note that there is no need in the
last condition for the stochastic approach, but we have to impose it since we
want to compare the stochastic approach with the results of perturbation
theory.

When $\lambda=0$, the equation of motion for the rescaled field $\chi\equiv
a(\eta)\phi$ is

$\chi^{\prime\prime}-\nabla^{2}\chi-{1\over\eta^{2}}\left(2-{m^{2}\over
H^{2}}\right)\chi=0\;,$

where primes denote derivatives with respect to the conformal time and
$\nabla^{2}$ is the three-dimensional Laplacian.

To quantize the field $\phi(\vec{x},t)$, we expand it in terms of creation and
annihilation operators

$\phi(\vec{x},t)=\int{d^{3}\vec{k}\over(2\pi)^{3}}\left\\{\phi_{k}(\eta)e^{i\vec{k}\cdot\vec{x}}a_{\vec{k}}+\phi_{k}^{*}(\eta)e^{-i\vec{k}\cdot\vec{x}}a_{\vec{k}}^{\dagger}\right\\}\;,$

where $a_{\vec{k}}^{\dagger}$ and $a_{\vec{k}}$ satisfy the usual commutation
relations. The mode functions $\chi_{k}\equiv a(\eta)\phi_{k}$ then obey the
differential equation

$\displaystyle\chi_{k}^{\prime\prime}+k^{2}\left[1-{1\over
k^{2}\eta^{2}}\left(2-{m^{2}\over H^{2}}\right)\right]\chi_{k}=0\;,$

where $k=|\vec{k}\rvert$. The general solution of this equation can be
expressed as a linear combination of Hankel functions:

$\chi_{k}(\eta)=\sqrt{-k\eta}\left[A_{k}\,{\mathcal{H}}_{\nu}^{(1)}(-k\eta)+B_{k}\,{\mathcal{H}}_{\nu}^{(2)}(-k\eta)\right]\;,$

where the index $\nu$ is

$\nu=\sqrt{{9\over 4}-\frac{m^{2}}{H^{2}}}\;.$ (6)

The choice of the coefficients $A_{k}$ and $B_{k}$ defines a vacuum state
$|0\rangle$ annihilated by $a_{\vec{k}}$:

$\displaystyle a_{\vec{k}}|0\rangle=0\;\;\text{for any $\vec{k}$}\;.$

If one wants to have a vacuum that in the remote past $\eta\to-\infty$ (or,
equivalently, for modes with very short physical wavelengths, $-kH\eta\gg H$)
behaves like the vacuum in Minkowski spacetime,

$\displaystyle\chi_{k}(\eta)\to{e^{-ik\eta}\over\sqrt{2k}}\;,$

one should choose

$\chi_{k}(\eta)=-{\sqrt{\pi}\over
2}\sqrt{-\eta}\,{\mathcal{H}}_{\nu}^{(1)}\left(-k\eta\right)\;.$

Such a choice is called the Bunch-Davies vacuum BD . The mode functions of the
original field, $\phi_{k}=a^{-1}\chi_{k}$, are then given by

$\phi_{k}(\eta)=-{\sqrt{\pi}\over
2}H(-\eta)^{3/2}{\mathcal{H}}_{\nu}^{(1)}(-k\eta)\;.$ (7)

As long as $m>0$, this state is de Sitter invariant Allen:1985ux .

Throughout this paper we will be assuming that $m^{2}\ll H^{2}$, in which case
the order (6) of the Hankel function becomes

$\nu\approx{3\over 2}-u,\text{ with $u\equiv{m^{2}\over 3H^{2}}\ll 1$}\;.$ (8)

##  3 Perturbative calculation of the two-point correlation function

In this section we present a perturbative calculation of the long-wavelength
part of the two-point function. There are two reasons why it is meaningful to
consider exclusively the long-wavelength modes.

The first reason is physical. The fluctuations relevant for the formation of
the observed large-scale structure of the universe are those whose wavelength,
by the end of inflation, has been stretched to a size much larger than the
Hubble horizon.

The second reason is mathematical. Calculations are much simpler if instead of
the modes (7), one uses their long-wavelength limit. At the same time, in many
cases the results can reflect the behavior of the untruncated theory. In the
small mass limit, the long-wavelength two-point function matches with the
untruncated one for large separations and for coinciding spacetime points. For
the free theory, this is explicitly shown at the end of Sec. 3.1. Perturbative
corrections (which we calculate in Sections 3.3 and 3.4) consist of products
of Wightman functions and retarded Green’s functions, integrated over momenta
and intermediate time coordinates. Here too the leading contribution in small
mass derives entirely from the long-wavelength limit of these functions.

###  3.1 The free theory

At the level of the free theory, the two-point function in a vacuum state
$|0\rangle$ defined by a set of modes $\phi_{k}(\eta)$ is given by

$\displaystyle\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{0}}$
$\displaystyle\equiv$ $\displaystyle\langle
0|\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})|0\rangle$ $\displaystyle=$
$\displaystyle\int\frac{d^{3}{\vec{k}}}{(2\pi)^{3}}\phi_{k}(\eta_{1})\phi_{k}^{*}(\eta_{2})e^{i\vec{k}\cdot(\vec{x}-\vec{y})}$
$\displaystyle=$ $\displaystyle{1\over
2\pi^{2}}\int_{0}^{\infty}dk\;k^{2}{\sin{kr}\over
kr}\phi_{k}(\eta_{1})\phi_{k}^{*}(\eta_{2})\,,$

where on the last line the integration over the angles was performed, and
$r\equiv|\vec{x}-\vec{y}\rvert$.

The long-wavelength part of (LABEL:2point) consists of modes with physical
momenta much less than $H$:

$\displaystyle{k\over a(\eta_{1})H}=-k\eta_{1}<\epsilon\;,\;\;{k\over
a(\eta_{2})H}=-k\eta_{2}<\epsilon\;,$

where $\epsilon$ is a small constant parameter, $\epsilon\ll 1$. Expanding the
mode functions (7) in this limit yields

$\displaystyle\phi_{k}(\eta)$ $\displaystyle\approx$
$\displaystyle{iH\over\sqrt{2}}(-\eta)^{3/2}(-k\eta)^{-\nu}={iH\over\sqrt{2k^{3}}}(-k\eta)^{u}.$
(10)

The long-wavelength part of the two-point function is then given by

$\displaystyle\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{0},L}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}={H^{2}(\eta_{1}\eta_{2})^{u}\over
4\pi^{2}}\int_{0}^{-\epsilon/\eta_{\rm m}}{dk\over k}{k^{2u}\sin(kr)\over
kr}\;,$ (11)

where $\eta_{\rm m}$ is the earliest time that accompanies the momentum
$\vec{k}$: $\eta_{\rm m}\equiv{\rm min}(\eta_{1},\eta_{2})$, and the subscript
$L$ stands for “long-wavelength”.

In the case of coinciding spacetime points, Eq. (11) gives

$\displaystyle\langle\phi^{2}(\vec{x},t)\rangle_{\lambda^{0},L}={H^{2}(-\eta)^{2u}\over
4\pi^{2}}\int_{0}^{-\epsilon/\eta}{dk\over k^{1-2u}}={H^{2}\over
8\pi^{2}}{\epsilon^{2u}\over u}.$ (12)

By expanding $\epsilon^{2u}$ for small values of $u$

$\epsilon^{2u}=1+2u\ln\epsilon+\mathcal{O}(u^{2}\ln^{2}\epsilon),$

we see that if we choose $\epsilon$ in such a way that it satisfies the
following conditions

$\exp\left(-u^{-1}\right)\ll\epsilon\ll 1,$ (13)

then $\epsilon^{2u}$ in Eq. (12) may be replaced by 1:

$\displaystyle\langle\phi^{2}(\vec{x},t)\rangle_{\lambda^{0},L}={H^{2}\over
8\pi^{2}u}={3H^{4}\over 8\pi^{2}m^{2}}\;.$ (14)

It is very important that $\langle\phi^{2}\rangle_{\lambda^{0},L}$ appears to
be anomalously large—much larger than $H^{2}$. This makes
$\langle\phi^{2}\rangle_{L}$ a very good approximation for the total
renormalized value of $\langle\phi^{2}\rangle_{\rm ren}$, both here and at
higher loops, since small-scale and renormalization contributions to this
quantity are of the order of $H^{2}$ (possibly up to a logarithmic
multiplier).

To consider a more general case, it is convenient to introduce a new
integration variable $w\equiv kr$ and rewrite Eq. (11) as

$\displaystyle\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{0},L}={H^{2}\over
4\pi^{2}}\left({\eta_{1}\eta_{2}\over r^{2}}\right)^{u}\int_{0}^{-\epsilon
r/\eta_{\rm m}}{\sin w\over w^{2-2u}}{dw}\;.$ (15)

In order to analyze this expression, let us look at the following quantity:
given two points in de Sitter spacetime, there is a de Sitter invariant
function associated with them Allen:1985ux ; Allen:1987tz :

$Z(X,Y)=-H^{2}\eta_{\mu\nu}X^{\mu}Y^{\nu},$ (16)

where $X$ and $Y$ represent coordinates in five-dimensional Minkowski
embedding spacetime with the metric $\eta_{\mu\nu}={\rm diag}(1,-1,-1,-1,-1)$.
We describe this function in more details in Appendix A. The explicit
expression for (16) depends on the coordinate system. In spatially flat
coordinates (4) that we use throughout this paper it has the form:

$\displaystyle
Z(\vec{x}_{1},\eta_{1};\vec{y},\eta_{2})={\eta_{1}^{2}+\eta_{2}^{2}-|\vec{x}-\vec{y}\rvert^{2}\over
2\eta_{1}\eta_{2}}\;.$ (17)

The important property of this quantity is that it allows us to distinguish
between timelike and spacelike related points. If points are timelike
separated, then $Z>1$; if points are lightlike separated, then $Z=1$, and if
points are spacelike separated, then $Z<1$.

If $(\vec{x},t_{1})$ and $(\vec{y},t_{2})$ are related in such a way that

$Z>1-{1\over 2\epsilon^{2}}\;,$ (18)

then it can be deduced from

$\displaystyle{\eta_{1}^{2}+\eta_{2}^{2}-|\vec{x}-\vec{y}\rvert^{2}\over
2\eta_{1}\eta_{2}}>1-{1\over 2\epsilon^{2}}\;$

that

${-r/\eta_{\rm m}}<{1/\epsilon}\;,$

where $\eta_{\rm m}\equiv{\rm min}(\eta_{1},\eta_{2})$ and
$r\equiv|\vec{x}-\vec{y}\rvert$. This means that the upper limit of the
integral in (15) is smaller than $1$, so we can use the approximation $\sin
w\approx w$ and obtain

$\displaystyle\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{0},L}$
$\displaystyle=$ $\displaystyle{H^{2}\over
4\pi^{2}}\left({\eta_{1}\eta_{2}\over r^{2}}\right)^{u}\int_{0}^{-\epsilon
r/\eta_{\rm m}}{dw\over w^{1-2u}}$ (19) $\displaystyle=$
$\displaystyle{H^{2}\over 8\pi^{2}}{\epsilon^{2u}\over
u}\bigg{(}{\eta_{1}\eta_{2}\over\eta_{\rm m}^{2}}\bigg{)}^{u}\;.$

As explained earlier, $\epsilon^{2u}$ can be replaced by $1$ if it satisfies
the condition (13). Hence in terms of the cosmic time we have

$\displaystyle\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{0},L}$
$\displaystyle=$ $\displaystyle{H^{2}\over
8\pi^{2}u}e^{-uH|t_{1}-t_{2}\rvert}\;.$ (20)

There are two important subcases of (18) for which the two-point function is
given by the above expression. One is when points $(\vec{x},t_{1})$ and
$(\vec{y},t_{2})$ are timelike or lightlike separated; the other is when these
points have coinciding time coordinates and the physical spatial distance
between them satisfies $a(t)r<(\epsilon H)^{-1}$. In the latter instance, from
Eq. (20) we obtain the same result as in Eq. (14). This means that as far as
the long-wavelength correlation function is concerned, there is no difference
between coinciding spacetime points and points on a constant time hypersurface
that are separated by a proper distance less than $(\epsilon H)^{-1}$.

Turning to the case where the upper limit of the integral in (15) is bigger
than $1$, that is,

$(-r/\eta_{\rm m})>1/\epsilon\;,$ (21)

it can be shown by using (17) that

$Z<1-{1\over 2\epsilon^{2}}\ll-1\;;$

therefore, this corresponds to the regime of large spacelike separation
between points. Let us split the integral in (15) into two parts

$\displaystyle\int_{0}^{-\epsilon r/\eta_{\rm m}}{\sin w\over
w^{2-2u}}{dw}=\\!\\!\int_{0}^{\infty}{\sin w\over
w^{2-2u}}{dw}-\int_{-\epsilon r/\eta_{\rm m}}^{\infty}{\sin w\over
w^{2-2u}}{dw}.$ (22)

The first integral can be calculated exactly

$\displaystyle\int_{0}^{\infty}{\sin w\over w^{2-2u}}{dw}=-\cos(\pi
u)\Gamma(2u-1)={1\over 2u}\bigg{(}1+\mathcal{O}(u)\bigg{)}\;,$

while the second integral can be shown to be smaller than $(1-2u)^{-1}$:

$\displaystyle\left|\int_{-\epsilon r/\eta_{\rm m}}^{\infty}{\sin w\over
w^{2-2u}}{dw}\right\rvert<\int_{-\epsilon r/\eta_{\rm
m}}^{\infty}{w^{2u-2}}{dw}<{1\over 1-2u}\;,$

where the second inequality follows from (21). Hence in the small mass limit
($u\ll 1$), the second term on the right-hand side of (22) can be neglected,
and we conclude that in the regime of large spacelike separation the two-point
function equals

$\displaystyle\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{0},L}$
$\displaystyle=$ $\displaystyle{H^{2}\over
8\pi^{2}u}\left({\eta_{1}\eta_{2}\over r^{2}}\right)^{u}$ $\displaystyle=$
$\displaystyle{H^{2}\over 8\pi^{2}u}e^{-uH(t_{1}+t_{2})}(rH)^{-2u}\;.$

As we can see, in this regime the equal-time two-point function,

$\displaystyle\langle\phi(\vec{x},t)\phi(\vec{y},t)\rangle_{\lambda^{0},L}={3H^{4}\over
8\pi^{2}m^{2}}(RH)^{-{2m^{2}\over 3H^{2}}}\;,$ (24)

depends only on the physical spatial distance $R\equiv re^{Ht}$.

Let us compare the long-wavelength results we have obtained so far with the
untruncated two-point function. This two-point function is calculated by using
Eq. (LABEL:2point) with mode functions given by (7). The result is well known
BD ; Allen:1987tz : for spacelike separated points, it is equal to

$\displaystyle\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{0}}$
$\displaystyle~{}~{}~{}~{}={H^{2}(1-u)(2-u)\over 16\pi\sin\pi
u}F\left(3-u,u;2;{1+Z\over 2}\right),$ (25)

where $F$ is a hypergeometric function, which is real for $Z<1$ and has a
branch cut that runs from $Z=1$ to $+\infty$. To study (25) in the limit of
coinciding spacetime points, we expand it around $Z=1$:

$\displaystyle\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{0}}\to{H^{2}\over
8\pi^{2}}\bigg{(}{1\over u}$ $\displaystyle-$ $\displaystyle 2+{1\over 1-Z}$
$\displaystyle-$ $\displaystyle\ln{(1-Z)\over 2}+\mathcal{O}(u)\bigg{)}$

as $(\vec{x},t_{1})\to(\vec{y},t_{2})$; $\mathcal{O}(u)$ contains terms that
are finite as well as terms that diverge as $Z\to 1$. It is evident from the
above expression that the term leading in $u$ is free of ultraviolet
divergences and coincides with the result (14), obtained by the long-
wavelength truncation.

When points are timelike separated, the untruncated correlation function
equals

$\displaystyle\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{0}}$
$\displaystyle~{}~{}~{}~{}={H^{2}(1-u)(2-u)\over 16\pi\sin\pi
u}F\left(3-u,u;2;{1+Z\pm i\epsilon\over 2}\right),$ (27)

where $\epsilon$ is a positive real infinitesimal. The plus sign in front of
$\epsilon$ refers to the case when $(\vec{x},t_{1})$ is in the past light cone
of $(\vec{y},t_{2})$, and the minus sign to the case when $(\vec{x},t_{1})$ is
in the forward light cone of $(\vec{y},t_{2})$. The $+i\epsilon$ and
$-i\epsilon$ prescriptions push the branch cut in the complex $(1-Z)$-plane
slightly above and slightly below the negative real axis, respectively.
Consequently, the sign of the imaginary part of (27) depends on the sign of
$t_{1}-t_{2}$: when $(t_{1}-t_{2})>0$, the imaginary part of (27) is negative,
and when $(t_{1}-t_{2})<0$, the imaginary part of (27) is positive.

It follows from Eqs. (25) and (27) that for $|Z|\gg 1$,

$\displaystyle\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{0}}\to{H^{2}\over
8\pi^{2}u}\bigg{(}1+\mathcal{O}(u^{0})\bigg{)}|Z|^{-u}\;.$ (28)

For large timelike separations, $Z$ is much greater than unity and can be
approximated as $Z\approx 0.5{e^{H|t_{1}-t_{2}\rvert}}$. Therefore, (28)
becomes

$\displaystyle\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{0}}\to{H^{2}\over
8\pi^{2}u}e^{-uH|t_{1}-t_{2}\rvert}\;,$

which matches with its long-wavelength counterpart (20). For large spacelike
separations, $Z\ll-1$ and, with the additional assumption that $-r/\eta_{\rm
m}\gg 1$, it can be approximated as
$Z\approx-0.5H^{2}r^{2}e^{H(t_{1}+t_{2})}$. Hence (28) reduces to

$\displaystyle\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{0}}\to{H^{2}\over
8\pi^{2}u}e^{-uH(t_{1}+t_{2})}(rH)^{-2u}\;,$

which also coincides with its long-wavelength analog (LABEL:tps).

Just as in flat spacetime, the expectation value of the commutator of two
fields vanishes for spacelike separated points and is nonzero for timelike
separated points. For example, if $(\vec{x},t_{1})$ is in the future of
$(\vec{y},t_{2})$, then

$\displaystyle\\!\\!\\!\\!\langle[\phi(\vec{x},t_{1}),\phi(\vec{y},t_{2})]\rangle_{\lambda^{0}}$
$\displaystyle={iH^{2}(1-u)(2-u)\over 8\pi\sin\pi u}{\rm
Im}\left\\{F\left(3-u,u;2;{1+Z-i\epsilon\over 2}\right)\right\\}.$

On the other hand, it follows from Eqs. (20) and (LABEL:tps) that

$\displaystyle\langle[\phi(\vec{x},t_{1}),\phi(\vec{y},t_{2})]\rangle_{\lambda^{0},L}=0$
(29)

both for timelike and spacelike related points. The vanishing of the
commutator in (29) indicates that the long-wavelength part of the field in a
sense behaves like a classical quantity. As we will see, the perturbative
corrections don’t change the value of the long-wavelength commutator: it
remains equal to zero.

###  3.2 The “in-in” formalism

To calculate perturbative corrections to the two-point function, we will work
in the interaction picture and use the “in-in” formalism Schwinger ; Keldysh ;
Bakshi:1962dv ; Jordan ; Weinberg:2005vy . In this formalism the two-point
function can be written as

$\displaystyle\\!\\!\\!\\!\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle$
$\displaystyle\quad\quad\quad\quad=\bigl{\langle}0\big{|}U^{\dagger}(t_{1},-\infty)\phi_{I}(\vec{x},t_{1})U(t_{1},-\infty)$
$\displaystyle\quad\quad\quad\quad\quad~{}~{}~{}\times
U^{\dagger}(t_{2},-\infty)\phi_{I}(\vec{y},t_{2})U(t_{2},-\infty)\big{|}0\bigr{\rangle}\;.$
(30)

Here, $|0\rangle$ is the vacuum state of the free theory and $U(t,-\infty)$ is
the unitary time-evolution operator:

$\displaystyle U(t,-\infty)$ $\displaystyle=$ $\displaystyle
Te^{-i\int_{-\infty}^{t}dt^{\prime}\,H_{I}(t^{\prime})}\;,$ $\displaystyle
U^{\dagger}(t,-\infty)$ $\displaystyle=$
$\displaystyle\overline{T}e^{i\int_{-\infty}^{t}dt^{\prime}\,H_{I}(t^{\prime})}\;,$

where $T$ stands for time-ordering, $\overline{T}$ — for anti-time-ordering.
Just as in “in-out” formalism, $\phi_{I}$-s are interaction picture fields
with time evolution governed by the free theory Hamiltonian and $H_{I}$ is the
interaction Hamiltonian in the interaction picture. In what follows we
suppress the subscript $I$ of the interaction picture fields.

Unlike in “in-out” formalism, conditions on fields and states are not imposed
at both very early and very late times, but only at very early times: both the
vacuum state $|0\rangle$ and its Hermitian conjugate are defined at $-\infty$.

For practical calculations it will be convenient to rewrite Eq. (30) in a
slightly different form. Let us represent the time-evolution operator and its
Hermitian conjugate as

$\displaystyle U(t_{1},-\infty)$ $\displaystyle=$ $\displaystyle
U^{\dagger}(\infty,t_{1})U(\infty,-\infty)\;,$ $\displaystyle
U^{\dagger}(t_{2},-\infty)$ $\displaystyle=$ $\displaystyle
U^{\dagger}(\infty,-\infty)U(\infty,t_{2})\;.$

Using the right-hand side of these equalities to replace the product
$U(t_{1},-\infty)U^{\dagger}(t_{2},-\infty)$ in (30) gives

$\displaystyle\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle$
$\displaystyle=$
$\displaystyle\bigl{\langle}0\big{|}U^{\dagger}(t_{1},-\infty)\phi(\vec{x},t_{1})U^{\dagger}(\infty,t_{1})$
$\displaystyle~{}~{}~{}\times
U(\infty,t_{2})\phi(\vec{y},t_{2})U(t_{2},-\infty)\big{|}0\bigr{\rangle}$
$\displaystyle=$
$\displaystyle\bigl{\langle}0\big{|}\overline{T}\bigl{(}\phi(\vec{x},t_{1})e^{i\int_{-\infty}^{\infty}dt^{\prime}\,H_{I}(t^{\prime})}\bigr{)}$
$\displaystyle~{}~{}~{}\times
T\bigl{(}\phi(\vec{y},t_{2})e^{-i\int_{-\infty}^{\infty}dt^{\prime}\,H_{I}(t^{\prime})}\bigr{)}\big{|}0\bigr{\rangle}\,.$

The second equality is obtained by noticing that on the right-hand side of the
first equality the first three operators are anti-time-ordered and the last
three operators are time-ordered.

The right-hand side of (LABEL:ctp1) can be interpreted in the following way:
we start with the initial state at $-\infty$, evolve forward to $t_{2}$ where
the operator $\phi(\vec{y},t_{2})$ is inserted, then continue to evolve
forward to $\infty$; afterwards, we evolve backward to $t_{1}$ where the
operator $\phi(\vec{x},t_{1})$ occurs, then continue to evolve back to
$-\infty$. This is why the “in-in” formalism is also called “closed-time-path”
formalism. This interpretation makes it possible to write (LABEL:ctp1) in
terms of a single time-ordered expression Jordan ; Collins:2011mz : label the
fields on the forward-flowing part of the path with a “$+$” superscript, and
the fields on the backward-flowing part of the path with a “$-$” superscript;
thereby, (LABEL:ctp1) can be written as

$\displaystyle\\!\\!\\!\\!\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle$
$\displaystyle=\bigl{\langle}0\big{|}T\bigl{(}\phi^{-}(\vec{x},t_{1})\phi^{+}(\vec{y},t_{2})e^{-i\int_{-\infty}^{\infty}dt^{\prime}\;[H_{I}^{+}(t^{\prime})-H_{I}^{-}(t^{\prime})]}\bigr{)}\big{|}0\bigr{\rangle}\,,$
(32)

where $H_{I}^{\pm}(t)\equiv H_{I}[\phi^{\pm}(t,\vec{x})]$, and the time-
ordering operation is extended in the following way: two “$+$” fields are
ordered as usual,

$\displaystyle T(\phi^{+}(\vec{x},t)\phi^{+}(\vec{y},t^{\prime}))$
$\displaystyle=$
$\displaystyle\theta(t-t^{\prime})\phi^{+}(\vec{x},t)\phi^{+}(\vec{y},t^{\prime})$
$\displaystyle+\theta(t^{\prime}-t)\phi^{+}(\vec{y},t^{\prime})\phi^{+}(\vec{x},t)\;,$

“$-$” fields always occur after “$+$” fields,

$\displaystyle T(\phi^{+}(\vec{x},t)\phi^{-}(\vec{y},t^{\prime}))$
$\displaystyle=$
$\displaystyle\phi^{-}(\vec{y},t^{\prime})\phi^{+}(\vec{x},t)\;,$
$\displaystyle T(\phi^{-}(\vec{x},t)\phi^{+}(\vec{y},t^{\prime}))$
$\displaystyle=$
$\displaystyle\phi^{-}(\vec{x},t)\phi^{+}(\vec{y},t^{\prime})\;,$

and two “$-$” fields are ordered in the opposite of the usual sense,

$\displaystyle T(\phi^{-}(\vec{x},t)\phi^{-}(\vec{y},t^{\prime}))$
$\displaystyle=$
$\displaystyle\theta(t^{\prime}-t)\phi^{-}(\vec{x},t)\phi^{-}(\vec{y},t^{\prime})$
$\displaystyle+\theta(t-t^{\prime})\phi^{-}(\vec{y},t^{\prime})\phi^{-}(\vec{x},t)\;.$

We can use Wick’s theorem to express the time-ordered products in (32) in
terms of field contractions, but unlike in “in-out” formalism, there are four
types of Wick contractions (and hence, four propagators)

$\displaystyle\langle
0|T(\phi^{+}(\vec{x},t)\phi^{+}(\vec{y},t^{\prime}))|0\rangle$
$\displaystyle~{}~{}=\theta(t-t^{\prime})G^{>}(\vec{x},t;\vec{y},t^{\prime})+\theta(t^{\prime}-t)G^{<}(\vec{x},t;\vec{y},t^{\prime})\;,$
$\displaystyle\langle
0|T(\phi^{+}(\vec{x},t^{\prime})\phi^{-}(\vec{y},t^{\prime}))|0\rangle=G^{<}(\vec{x},t;\vec{y},t^{\prime})\;,$
$\displaystyle\langle
0|T(\phi^{-}(\vec{x},t)\phi^{+}(\vec{y},t^{\prime}))|0\rangle=G^{>}(\vec{x},t;\vec{y},t^{\prime})\;,$
$\displaystyle\langle
0|T(\phi^{-}(\vec{x},t)\phi^{-}(\vec{y},t^{\prime}))|0\rangle$
$\displaystyle~{}~{}=\theta(t^{\prime}-t)G^{>}(\vec{x},t;\vec{y},t^{\prime})+\theta(t-t^{\prime})G^{<}(\vec{x},t;\vec{y},t^{\prime})\;,$

where $G^{>}(\vec{x},t;\vec{y},t^{\prime})$ and
$G^{<}(\vec{x},t;\vec{y},t^{\prime})$ are the Wightman functions

$\displaystyle G^{>}(\vec{x},t;\vec{y},t^{\prime})$ $\displaystyle\equiv$
$\displaystyle\langle 0|\phi(\vec{x},t)\phi(\vec{y},t^{\prime})|0\rangle$
$\displaystyle=$
$\displaystyle\int\frac{d^{3}\vec{k}}{(2\pi)^{3}}\;e^{i{\vec{k}(\vec{x}-\vec{y})}}G_{k}^{>}(t,t^{\prime})\;,$
$\displaystyle G^{<}(\vec{x},t;\vec{y},t^{\prime})$ $\displaystyle\equiv$
$\displaystyle\langle 0|\phi(\vec{y},t^{\prime})\phi(\vec{x},t)|0\rangle$
$\displaystyle=$
$\displaystyle\int\frac{d^{3}\vec{k}}{(2\pi)^{3}}\;e^{i{\vec{k}(\vec{x}-\vec{y})}}G_{k}^{<}(t,t^{\prime})\;.$

In terms of the conformal time, the momentum representations of the Wightman
functions are given by

$\displaystyle\\!\\!\\!\\!G_{k}^{>}(\eta,\eta^{\prime})=\phi_{k}(\eta)\phi_{k}^{*}(\eta^{\prime})$
$\displaystyle={\pi H^{2}\over
4}(-\eta)^{3/2}(-\eta^{\prime})^{3/2}{\mathcal{H}}_{\nu}^{(1)}(-k\eta){\mathcal{H}}_{\nu}^{(2)}(-k\eta^{\prime})\;,$
(33) $\displaystyle
G_{k}^{<}(\eta,\eta^{\prime})=\phi_{k}^{*}(\eta)\phi_{k}(\eta^{\prime})=G_{k}^{>}(\eta^{\prime},\eta)\;,$
(34)

where we have used the mode functions (7) associated with the Bunch-Davies
state.

###  3.3 The one-loop correction

The interaction Hamiltonian of the scalar field with quartic self-interaction
is given by

$\displaystyle H_{I}(t)=\int d^{3}\vec{x}\;\left(a^{3}(t){\lambda\over
4}\phi^{4}\right)\;.$

The one-loop correction to the two-point function can be obtained by expanding
Eq. (32) to first order in the coupling constant $\lambda$:

$\displaystyle\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda}$
$\displaystyle~{}~{}=-i\int_{-\infty}^{\infty}dt^{\prime}\;\langle
0|T\bigg{(}\phi^{-}(\vec{x},t_{1})\phi^{+}(\vec{y},t_{2})$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad~{}\times\Big{[}H_{I}^{+}(t^{\prime})-H_{I}^{-}(t^{\prime})\Big{]}\bigg{)}|0\rangle\;.$
$(\vec{x},t_{1})$$t^{\prime}$$(\vec{y},t_{2})$Fig.1

Taking the contractions that lead to the connected diagram (Fig.1) and
switching from the cosmic time $t$ to the conformal time $\eta$ yields

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda}$
$\displaystyle=$
$\displaystyle-{3i\lambda}\int_{-\infty}^{0}\\!\\!{d\eta^{\prime}}a^{4}(\eta^{\prime})\int\\!\\!\frac{d^{3}{\vec{k}}}{(2\pi)^{3}}{e^{i\vec{k}\cdot(\vec{x}-\vec{y})}}$
$\displaystyle\times\bigg{[}\theta(\eta_{1}-\eta^{\prime})G_{k}^{C}(\eta_{1},\eta^{\prime})G_{k}^{<}(\eta_{2},\eta^{\prime})$
$\displaystyle+\theta(\eta_{2}-\eta^{\prime})G_{k}^{C}(\eta_{2},\eta^{\prime})G_{k}^{>}(\eta_{1},\eta^{\prime})\bigg{]}\\!\\!\int\\!\\!\frac{d^{3}{\vec{\ell}}}{(2\pi)^{3}}G_{\ell}^{>}(\eta^{\prime},\eta^{\prime}),$

where $G_{k}^{C}(\eta,\eta^{\prime})\equiv
G_{k}^{>}(\eta,\eta^{\prime})-G_{k}^{<}(\eta,\eta^{\prime})$. Note that the
product $\theta(\eta-\eta^{\prime})G_{k}^{C}(\eta,\eta^{\prime})$ is nothing
but the momentum representation of the retarded Green’s function.

To evaluate the long-wavelength part of (LABEL:1l), let us look at the small
(physical) momenta expansion of the Wightman functions. It follows from Eq.
(10) that

$G_{k}^{>}(\eta,\eta^{\prime})\approx{H^{2}\over
2k^{3}}(-k\eta)^{u}(-k\eta^{\prime})^{u}\approx
G_{k}^{<}(\eta,\eta^{\prime})\;.$ (36)

The long-wavelength part of the integral over $\vec{\ell}$ in Eq. (LABEL:1l)
is equal to $\langle\phi^{2}\rangle_{\lambda^{0},L}$ (cf. Eq (12)):

$\displaystyle\int^{-\epsilon/\eta^{\prime}}\\!\\!\frac{d^{3}{\vec{\ell}}}{(2\pi)^{3}}G_{\ell}^{>}(\eta^{\prime},\eta^{\prime})$
$\displaystyle\approx$ $\displaystyle{H^{2}\over
4\pi^{2}}\int_{0}^{-\epsilon/\eta^{\prime}}\\!\\!{d\ell\over\ell}(-\ell\eta^{\prime})^{2u}$
(37) $\displaystyle=$ $\displaystyle{H^{2}\over 8\pi^{2}}{\epsilon^{2u}\over
u}\;.$

In order to extract the leading terms from the difference of the Wightman
functions, it is convenient to rewrite it in terms of the Bessel functions of
the first and second kind; from Eqs. (33) and (34) we have

$\displaystyle G_{k}^{C}(\eta,\eta^{\prime})$ $\displaystyle=$ $\displaystyle
i{\pi H^{2}\over 2}(-\eta)^{3/2}(-\eta^{\prime})^{3/2}$
$\displaystyle\times\biggl{[}Y_{\nu}(-k\eta)J_{\nu}(-k\eta^{\prime})-J_{\nu}(-k\eta)Y_{\nu}(-k\eta^{\prime})\biggr{]},$

which in the long-wavelength limit gives

$\displaystyle G_{k}^{C}(\eta,\eta^{\prime})$ $\displaystyle\approx$
$\displaystyle-{i\over 2\nu}H^{2}(-\eta)^{3/2}(-\eta^{\prime})^{3/2}$
$\displaystyle\times\biggl{[}(-k\eta)^{-\nu}(-k\eta^{\prime})^{\nu}-(-k\eta)^{\nu}(-k\eta^{\prime})^{-\nu}\biggr{]}$
$\displaystyle\approx$ $\displaystyle-{i\over
3}H^{2}(-\eta)^{u}(-\eta^{\prime})^{3-u}+{i\over
3}H^{2}(-\eta)^{3-u}(-\eta^{\prime})^{u}.$

In what follows we will only keep the first term of this approximation,

$\displaystyle G_{k}^{C}(\eta,\eta^{\prime})$ $\displaystyle\approx$
$\displaystyle-{i\over 3}H^{2}(-\eta)^{u}(-\eta^{\prime})^{3-u}\;,$ (38)

since (as it will become clear from the calculations below) the part of the
two-point function that would arise due to the second term would be suppressed
by an extra factor of $u$.

Substituting Eqs. (36), (37) and (38) into Eq. (LABEL:1l) gives

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda,L}$
$\displaystyle=$ $\displaystyle{\lambda H^{2}\over
16\pi^{2}}{\epsilon^{2u}\over
u}(\eta_{1}\eta_{2})^{u}\bigg{\\{}\int_{-\infty}^{\eta_{1}}{d\eta^{\prime}\over\eta^{\prime}}\int^{-\epsilon/\eta_{\rm
m}}\frac{d^{3}{\vec{k}}}{(2\pi)^{3}}{e^{i\vec{k}\cdot(\vec{x}-\vec{y})}\over
k^{3-2u}}$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad+\int_{-\infty}^{\eta_{2}}{d\eta^{\prime}\over\eta^{\prime}}\int^{-\epsilon/\eta_{\rm
m}}\frac{d^{3}{\vec{k}}}{(2\pi)^{3}}{e^{i\vec{k}\cdot(\vec{x}-\vec{y})}\over
k^{3-2u}}\bigg{\\}}$ $\displaystyle=$ $\displaystyle{\lambda H^{2}\over
32\pi^{4}}{\epsilon^{2u}\over u}\left({\eta_{1}\eta_{2}\over
r^{2}}\right)^{u}\\!\bigg{\\{}\int_{-\infty}^{\eta_{1}}{d\eta^{\prime}\over\eta^{\prime}}\int_{0}^{-\epsilon
r/\eta_{\rm m}}{\sin w\over w^{2-2u}}{dw}$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad+\int_{-\infty}^{\eta_{2}}{d\eta^{\prime}\over\eta^{\prime}}\int_{0}^{-\epsilon
r/\eta_{\rm m}}{\sin w\over w^{2-2u}}{dw}\bigg{\\}}\;,$

where $\eta_{\rm m}$ is the earliest time that accompanies the momentum
$\vec{k}$: $\eta_{\rm m}\equiv{\rm min}(\eta_{1},\eta_{2},\eta^{\prime})$.
This momentum cutoff is chosen in such a way as to ensure that all the
physical scales associated with the comoving scale $k$ are smaller than
$\epsilon H$ (that is, $-k\eta_{1}<\epsilon$, $-k\eta_{2}<\epsilon$ and
$-k\eta^{\prime}<\epsilon$). To make the calculations more tractable, let us
assume that $t_{1}\geq t_{2}$ (that is, $\eta_{1}\geq\eta_{2}$), in which case
$\eta_{\rm m}={\rm min}(\eta_{2},\eta^{\prime})$. In the second term in curly
brackets, $\eta_{\rm m}$ is then equal to $\eta^{\prime}$. In the first term,
$\eta_{\rm m}$ is determined by splitting the time integral into two parts

$\displaystyle\int_{-\infty}^{\eta_{1}}{d\eta^{\prime}\over\eta^{\prime}}\int_{0}^{-\epsilon
r/\eta_{\rm m}}{\sin w\over w^{2-2u}}{dw}$ (40)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!=\int_{-\infty}^{\eta_{2}}{d\eta^{\prime}\over\eta^{\prime}}\\!\int_{0}^{-\epsilon
r/\eta^{\prime}}\\!{\sin w\over w^{2-2u}}{dw}$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+\int_{\eta_{2}}^{\eta_{1}}{d\eta^{\prime}\over\eta^{\prime}}\\!\int_{0}^{-\epsilon
r/\eta_{2}}\\!{\sin w\over w^{2-2u}}{dw}\;;$

hence, Eq. (LABEL:1ll) becomes

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda,L}$
$\displaystyle=$ $\displaystyle{\lambda H^{2}\over
32\pi^{4}}{\epsilon^{2u}\over u}\left({\eta_{1}\eta_{2}\over
r^{2}}\right)^{u}\bigg{\\{}2\int_{-\infty}^{\eta_{2}}{d\eta^{\prime}\over\eta^{\prime}}\\!\int_{0}^{-\epsilon
r/\eta^{\prime}}\\!{\sin w\over w^{2-2u}}{dw}$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad~{}+\int_{\eta_{2}}^{\eta_{1}}{d\eta^{\prime}\over\eta^{\prime}}\\!\int_{0}^{-\epsilon
r/\eta_{2}}\\!{\sin w\over w^{2-2u}}{dw}\bigg{\\}}\;.$

If the points $(\vec{x},t_{1})$ and $(\vec{y},t_{2})$ are related in such a
way that $Z>1-(2\epsilon^{2})^{-1}$, the upper limit of the integral over $w$
in the second term in curly brackets is less than unity. Since $\eta^{\prime}$
in the first term is less than $\eta_{2}$, the upper limit of the integral
over $w$ in here is also less than unity. Hence in both terms $\sin{w}$ can be
replaced by $w$,

$\int_{0}^{-\epsilon r/\eta}{\sin w\over
w^{2-2u}}{dw}\approx\int_{0}^{-\epsilon r/\eta}{dw\over
w^{1-2u}}={\epsilon^{2u}\over 2u}\left({r\over-\eta}\right)^{2u}\;,$

after which we obtain

$\displaystyle\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda,L}$
$\displaystyle~{}~{}~{}~{}=-{\lambda H^{2}\over 64\pi^{4}}{\epsilon^{4u}\over
u^{3}}\left({\eta_{1}\over\eta_{2}}\right)^{u}\left(1-u\ln{\eta_{1}\over\eta_{2}}\right)$
$\displaystyle~{}~{}~{}~{}=-{\lambda H^{2}\over 64\pi^{4}}{\epsilon^{4u}\over
u^{3}}\bigg{(}1+uH(t_{1}-t_{2})\bigg{)}e^{-uH(t_{1}-t_{2})}\;.$

Had we assumed that $t_{2}>t_{1}$, this result would still be true but with
$t_{1}$ and $t_{2}$ exchanged. Therefore it can be generalized as

$\displaystyle\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda,L}$
$\displaystyle~{}~{}~{}~{}=-{\lambda H^{2}\over
64\pi^{4}u^{3}}\bigg{(}1+uH|t_{1}-t_{2}\rvert\bigg{)}e^{-uH|t_{1}-t_{2}\rvert}\;.$
(42)

Specifically, in the case of coinciding spacetime points one finds

$\langle\phi^{2}(\vec{x},t)\rangle_{\lambda,L}=-{27\lambda H^{8}\over
64\pi^{4}m^{6}}\;.$ (43)

Note that even though the second term in the parenthesis of (42) has an
additional factor of $u$, this term cannot be neglected because it is greater
than unity when $|t_{1}-t_{2}\rvert>1/u$.

As explained in Sec. 3.1, when $-\epsilon r/\eta_{2}>1$ (and consequently
$Z<1-{1\over 2\epsilon^{2}}$), the integral over $w$ in the second term in the
curly brackets of Eq (LABEL:1loop) can be approximated as

$\displaystyle\int_{0}^{-\epsilon r/\eta_{2}}{\sin w\over
w^{2-2u}}{dw}\approx\int_{0}^{\infty}{\sin w\over w^{2-2u}}{dw}={1\over
2u}\bigg{(}1+\mathcal{O}(u)\bigg{)}\;.$ (44)

Performing the $\eta^{\prime}$ integration then produces

$\displaystyle\int_{\eta_{2}}^{\eta_{1}}{d\eta^{\prime}\over\eta^{\prime}}\\!\int_{0}^{-\epsilon
r/\eta_{2}}\\!{\sin w\over w^{2-2u}}{dw}\approx{1\over
2u}\ln\left({\eta_{1}\over\eta_{2}}\right)\;.$ (45)

To evaluate the first term, we split the time integral as follows:

$\displaystyle
2\int_{-\infty}^{\eta_{2}}{d\eta^{\prime}\over\eta^{\prime}}\\!\int_{0}^{-\epsilon
r/\eta^{\prime}}\\!{\sin w\over w^{2-2u}}{dw}$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!=2\int_{-\infty}^{-\epsilon
r}{d\eta^{\prime}\over\eta^{\prime}}\\!\int_{0}^{-\epsilon
r/\eta^{\prime}}\\!{\sin w\over w^{2-2u}}{dw}$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+2\int_{-\epsilon
r}^{\eta_{2}}{d\eta^{\prime}\over\eta^{\prime}}\\!\int_{0}^{-\epsilon
r/\eta^{\prime}}\\!{\sin w\over w^{2-2u}}{dw}\;.$

The second term on the right side of the above equality should be treated the
same way as integrals (44) and (45),

$\displaystyle 2\int_{-\epsilon
r}^{\eta_{2}}{d\eta^{\prime}\over\eta^{\prime}}\\!\int_{0}^{-\epsilon
r/\eta^{\prime}}\\!{\sin w\over w^{2-2u}}{dw}\approx{1\over
u}\ln{\left(-\eta_{2}\over\epsilon r\right)}\;,$ (46)

since the ratio $-\epsilon r/\eta^{\prime}$ in here stays greater than unity.
In the first term this ration stays smaller than unity, and hence this term
can be calculated as

$\displaystyle 2\int_{-\infty}^{-\epsilon
r}{d\eta^{\prime}\over\eta^{\prime}}\\!\int_{0}^{-\epsilon
r/\eta^{\prime}}\\!{\sin w\over w^{2-2u}}{dw}$
$\displaystyle\quad\quad\quad\approx 2\int_{-\infty}^{-\epsilon
r}{d\eta^{\prime}\over\eta^{\prime}}\\!\int_{0}^{-\epsilon
r/\eta^{\prime}}\\!{dw\over w^{1-2u}}=-{1\over 2u^{2}}\;.$ (47)

Upon substituting Eqs. (45), (46) and (47) into Eq. (LABEL:1loop), we obtain
the long-wavelength correlation function at large spacelike separations:

$\displaystyle\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda,L}$
$\displaystyle\quad\quad=-{\lambda H^{2}\over 64\pi^{4}}{\epsilon^{2u}\over
u^{3}}\left({\eta_{1}\eta_{2}\over
r^{2}}\right)^{u}\left[1+u\ln\left({\epsilon^{2}r^{2}\over\eta_{1}\eta_{2}}\right)\right]\;,$

which in terms of cosmic time becomes

$\displaystyle\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda,L}$
$\displaystyle=$ $\displaystyle-{\lambda H^{2}\over
64\pi^{4}u^{3}}\left(r^{2}H^{2}e^{H(t_{1}+t_{2})}\right)^{-u}$
$\displaystyle\times\bigg{\\{}1+u\ln\left({r^{2}H^{2}e^{H(t_{1}+t_{2})}}\right)\bigg{\\}}\;.$

If we considered the case $t_{2}>t_{1}$, this result would still be true.

The comparison of Eqs. (42) and (LABEL:1loops) with (20) and (LABEL:tps)
reveals that the effective parameter of the perturbative expansion is not
$\lambda$ but $\lambda/u^{2}$, so the perturbation theory is valid as long as
$\lambda\ll m^{4}/H^{4}$.

###  3.4 The two-loop correction

Upon expanding (32) to second order in $\lambda$ and taking all of the
possible field contractions that produce connected graphs, we obtain several
diagrams with different topologies that contribute to the two-point function
at $\lambda^{2}$-order.

####  3.4.1 Two independent loops

$(\vec{x},t_{1})$$t^{\prime}$$t^{\prime\prime}$$(\vec{y},t_{2})$Fig.2

Taking the field contractions that correspond to the diagram with two
independent loops (Fig.2), one finds

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{2}}^{(1)}$
$\displaystyle=$
$\displaystyle-9\lambda^{2}\bigg{\\{}\int_{-\infty}^{0}\\!\\!d\eta^{\prime}a^{4}(\eta^{\prime})\int_{-\infty}^{\eta^{\prime}}\\!\\!d\eta^{\prime\prime}a^{4}(\eta^{\prime\prime})$
$\displaystyle\times\int\frac{d^{3}{\vec{k}}}{(2\pi)^{3}}{e^{i\vec{k}\cdot(\vec{x}-\vec{y})}}\bigg{[}\theta(\eta_{1}-\eta^{\prime})G_{k}^{C}(\eta_{1},\eta^{\prime})G_{k}^{<}(\eta_{2},\eta^{\prime\prime})$
$\displaystyle~{}~{}~{}~{}+\theta(\eta_{2}-\eta^{\prime})G_{k}^{C}(\eta_{2},\eta^{\prime})G_{k}^{>}(\eta_{1},\eta^{\prime\prime})\bigg{]}G_{k}^{C}(\eta^{\prime},\eta^{\prime\prime})$
$\displaystyle+\int_{-\infty}^{\eta_{1}}\\!\\!{d\eta^{\prime}}a^{4}(\eta^{\prime})\int_{-\infty}^{\eta_{2}}\\!\\!{d\eta^{\prime\prime}}a^{4}(\eta^{\prime\prime})\int\frac{d^{3}{\vec{k}}}{(2\pi)^{3}}{e^{i\vec{k}\cdot(\vec{x}-\vec{y})}}$
$\displaystyle\times
G_{k}^{C}(\eta_{1},\eta^{\prime})G_{k}^{C}(\eta_{2},\eta^{\prime\prime})G_{k}^{>}(\eta^{\prime},\eta^{\prime\prime})\bigg{\\}}$
$\displaystyle\times\\!\int\\!\frac{d^{3}{\vec{\ell}}}{(2\pi)^{3}}G_{\ell}^{>}(\eta^{\prime},\eta^{\prime})\\!\int\\!\frac{d^{3}{\vec{p}}}{(2\pi)^{3}}G_{p}^{>}(\eta^{\prime\prime},\eta^{\prime\prime}).$

Using the small momenta approximations (36), (37) and (38) gives

$\displaystyle\\!\\!\\!\\!\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{2},L}^{(1)}={\lambda^{2}H^{2}\over
256\pi^{6}}{\epsilon^{4u}\over u^{2}}\left({\eta_{1}\eta_{2}\over
r^{2}}\right)^{u}$
$\displaystyle~{}~{}\times\bigg{\\{}\int_{-\infty}^{\eta_{1}}{d\eta^{\prime}\over\eta^{\prime}}\int_{-\infty}^{\eta^{\prime}}{d\eta^{\prime\prime}\over\eta^{\prime\prime}}+\int_{-\infty}^{\eta_{2}}{d\eta^{\prime}\over\eta^{\prime}}\int_{-\infty}^{\eta^{\prime}}{d\eta^{\prime\prime}\over\eta^{\prime\prime}}$
$\displaystyle~{}~{}\quad\quad+\int_{-\infty}^{\eta_{1}}{d\eta^{\prime}\over\eta^{\prime}}\int_{-\infty}^{\eta_{2}}{d\eta^{\prime\prime}\over\eta^{\prime\prime}}\bigg{\\}}\int_{0}^{-\epsilon
r/\eta_{\rm m}}{\sin w\over w^{2-2u}}{dw}\,,$ (49)

where $\eta_{\rm m}$ is the earliest time that accompanies the momentum
$\vec{k}$: $\eta_{\rm m}\equiv{\rm
min}(\eta_{1},\eta_{2},\eta^{\prime},\eta^{\prime\prime})$. We can evaluate
(49) by going through steps similar to the ones in the previous subsection; it
is slightly more tedious here since splitting the time integrals produces more
terms.

When $Z>1-(2\epsilon^{2})^{-1}$, we obtain

$\displaystyle\\!\\!\\!\\!\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{2},L}^{(1)}={\lambda^{2}H^{2}\over
512\pi^{6}u^{5}}e^{-uH|t_{1}-t_{2}\rvert}$
$\displaystyle\quad\quad\times\bigg{(}1+uH|t_{1}-t_{2}\rvert+{1\over
2}u^{2}H^{2}|t_{1}-t_{2}\rvert^{2}\bigg{)}\,,$ (50)

which for coinciding spacetime points becomes

$\displaystyle\langle\phi^{2}(\vec{x},t)\rangle_{\lambda^{2},L}^{(1)}={243\lambda^{2}H^{12}\over
512\pi^{6}m^{10}}\;.$ (51)

In $Z<1-(2\epsilon^{2})^{-1}$ regime, with the additional assumption that
$-\epsilon r/{\rm min}(\eta_{1},\eta_{2})>1$, evaluating (49) yields

$\displaystyle\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{2},L}^{(1)}={\lambda^{2}H^{2}\over
512\pi^{6}u^{5}}\left(r^{2}H^{2}e^{H(t_{1}+t_{2})}\right)^{-u}$
$\displaystyle\quad\quad\quad\quad\times\bigg{\\{}1+u\ln\left({r^{2}H^{2}e^{H(t_{1}+t_{2})}}\right)$
$\displaystyle\quad\quad\quad\quad\quad~{}+{1\over
2}u^{2}\ln^{2}\left({r^{2}H^{2}e^{H(t_{1}+t_{2})}}\right)\bigg{\\}}.$ (52)

####  3.4.2 Snowman diagram

$(\vec{x},t_{1})$$t^{\prime}$$t^{\prime\prime}$$(\vec{y},t_{2})$Fig. 3

Next, we consider the snowman diagram (Fig.3). Contractions corresponding to
this diagram produce the following result:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{2}}^{(2)}$
$\displaystyle=$
$\displaystyle-9\lambda^{2}\int_{-\infty}^{0}\\!\\!{d\eta^{\prime}}a^{4}(\eta^{\prime})\int\frac{d^{3}{\vec{k}}}{(2\pi)^{3}}{e^{i\vec{k}\cdot(\vec{x}-\vec{y})}}$
$\displaystyle\times\bigg{[}\theta(\eta_{1}-\eta^{\prime})G_{k}^{C}(\eta_{1},\eta^{\prime})G_{k}^{<}(\eta_{2},\eta^{\prime})$
$\displaystyle~{}~{}+\theta(\eta_{2}-\eta^{\prime})G_{k}^{C}(\eta_{2},\eta^{\prime})G_{k}^{>}(\eta_{1},\eta^{\prime})\bigg{]}\int_{-\infty}^{\eta^{\prime}}{d\eta^{\prime\prime}}a^{4}(\eta^{\prime\prime})$
$\displaystyle\times\int\frac{d^{3}{\vec{p}}}{(2\pi)^{3}}\bigg{(}G^{>}_{p}(\eta^{\prime},\eta^{\prime\prime})+G^{<}_{p}(\eta^{\prime},\eta^{\prime\prime})\bigg{)}G_{p}^{C}(\eta^{\prime},\eta^{\prime\prime})$
$\displaystyle\times\int\frac{d^{3}{\vec{\ell}}}{(2\pi)^{3}}G_{\ell}^{>}(\eta^{\prime\prime},\eta^{\prime\prime})\;.$

The long-wavelength approximation of this expression has the form:

$\displaystyle\\!\\!\\!\\!\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{2},L}^{(2)}={\lambda^{2}H^{2}\over
64\pi^{6}}{\epsilon^{2u}\over u}\left(\eta_{1}\eta_{2}\over r^{2}\right)^{u}$
$\displaystyle~{}~{}\times\bigg{[}\int_{-\infty}^{\eta_{1}}{d\eta^{\prime}\over\eta^{\prime}}(-\eta^{\prime})^{2u}+\int_{-\infty}^{\eta_{2}}{d\eta^{\prime}\over\eta^{\prime}}(-\eta^{\prime})^{2u}\bigg{]}$
$\displaystyle~{}~{}\times\int_{0}^{-\epsilon r/\eta_{\rm m}}{\sin w\over
w^{2-2u}}{dw}\int_{-\infty}^{\eta^{\prime}}{d\eta^{\prime\prime}\over\eta^{\prime\prime}}\int_{0}^{-\epsilon/\eta^{\prime\prime}}{dp\over
p^{1-2u}}\;,$ (53)

where, as before, $\eta_{\rm m}$ is the earliest time that accompanies the
momentum $\vec{k}$: $\eta_{\rm m}\equiv{\rm
min}(\eta_{1},\eta_{2},\eta^{\prime})$. Evaluating this in the case
$Z>1-(2\epsilon^{2})^{-1}$ gives

$\displaystyle\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{2},L}^{(2)}$
$\displaystyle~{}~{}~{}~{}={\lambda^{2}H^{2}\over
512\pi^{6}u^{5}}\bigg{(}1+uH|t_{1}-t_{2}\rvert\bigg{)}e^{-uH|t_{1}-t_{2}\rvert}\;,$
(54)

which for coinciding spacetime points reduces to

$\displaystyle\langle\phi^{2}(\vec{x},t)\rangle_{\lambda^{2},L}^{(2)}={243\lambda^{2}H^{12}\over
512\pi^{6}m^{10}}\;.$ (55)

On the other hand, when $Z<1-(2\epsilon^{2})^{-1}$, Eq. (53) gives

$\displaystyle\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{2},L}^{(2)}$
$\displaystyle=$ $\displaystyle{\lambda^{2}H^{2}\over
512\pi^{6}u^{5}}\left(r^{2}H^{2}e^{H(t_{1}+t_{2})}\right)^{-u}$
$\displaystyle\times\left\\{1+u\ln\left({r^{2}H^{2}e^{H(t_{1}+t_{2})}}\right)\right\\}\;.$

####  3.4.3 Sunset diagram

$(\vec{x},t_{1})$$t^{\prime}$$t^{\prime\prime}$$(\vec{y},t_{2})$Fig. 4

For the field contractions that correspond to the sunset diagram (Fig.4) we
obtain the following expression:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{2}}^{(3)}$
(57) $\displaystyle=$
$\displaystyle-6\lambda^{2}\bigg{\\{}\int_{-\infty}^{0}{d\eta^{\prime}}a^{4}(\eta^{\prime})\int_{-\infty}^{\eta^{\prime}}{d\eta^{\prime\prime}}a^{4}(\eta^{\prime\prime})$
$\displaystyle\times\int\frac{d^{3}{\vec{k}}}{(2\pi)^{3}}\frac{d^{3}{\vec{q}}}{(2\pi)^{3}}\frac{d^{3}{\vec{p}}}{(2\pi)^{3}}\frac{d^{3}{\vec{\ell}}}{(2\pi)^{3}}$
$\displaystyle\times(2\pi)^{3}\delta^{3}(\vec{k}+\vec{p}+\vec{l}+\vec{q}){e^{i\vec{k}\cdot(\vec{x}-\vec{y})}}$
$\displaystyle\times\bigg{[}\theta(\eta_{1}-\eta^{\prime})G_{k}^{C}(\eta_{1},\eta^{\prime})G_{k}^{<}(\eta_{2},\eta^{\prime\prime})$
$\displaystyle~{}~{}~{}+\theta(\eta_{2}-\eta^{\prime})G_{k}^{C}(\eta_{2},\eta^{\prime})G_{k}^{>}(\eta_{1},\eta^{\prime\prime})\bigg{]}$
$\displaystyle\times\bigg{[}G_{p}^{>}(\eta^{\prime},\eta^{\prime\prime})G_{q}^{>}(\eta^{\prime},\eta^{\prime\prime})+G_{p}^{>}(\eta^{\prime},\eta^{\prime\prime})G_{q}^{<}(\eta^{\prime},\eta^{\prime\prime})$
$\displaystyle~{}~{}~{}+G_{p}^{<}(\eta^{\prime},\eta^{\prime\prime})G_{q}^{<}(\eta^{\prime},\eta^{\prime\prime})\bigg{]}G_{\ell}^{C}(\eta^{\prime},\eta^{\prime\prime})\\!\\!\\!\\!\\!\\!$
$\displaystyle+\int_{-\infty}^{\eta_{1}}{d\eta^{\prime}}a^{4}(\eta^{\prime})\int_{-\infty}^{\eta_{2}}{d\eta^{\prime\prime}}a^{4}(\eta^{\prime\prime})$
$\displaystyle\times\int\frac{d^{3}{\vec{k}}}{(2\pi)^{3}}\frac{d^{3}{\vec{q}}}{(2\pi)^{3}}\frac{d^{3}{\vec{p}}}{(2\pi)^{3}}\frac{d^{3}{\vec{\ell}}}{(2\pi)^{3}}$
$\displaystyle\times(2\pi)^{3}\delta^{3}(\vec{k}+\vec{p}+\vec{l}+\vec{q}){e^{i\vec{k}\cdot(\vec{x}-\vec{y})}}G_{k}^{C}(\eta_{1},\eta^{\prime})$
$\displaystyle\times
G_{k}^{C}(\eta_{2},\eta^{\prime\prime})G_{\ell}^{>}(\eta^{\prime},\eta^{\prime\prime})G_{p}^{>}(\eta^{\prime},\eta^{\prime\prime})G_{q}^{>}(\eta^{\prime},\eta^{\prime\prime})\bigg{\\}}\;.$

Let us use the delta function on the fourth line to perform the integral over
$\vec{\ell}$ and the delta function on the second to last line to perform the
integral over $\vec{k}$:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{2}}^{(3)}$
$\displaystyle=$
$\displaystyle-6\lambda^{2}\bigg{\\{}\int_{-\infty}^{0}{d\eta^{\prime}}a^{4}(\eta^{\prime})\int_{-\infty}^{\eta^{\prime}}{d\eta^{\prime\prime}}a^{4}(\eta^{\prime\prime})$
$\displaystyle\times\int\frac{d^{3}{\vec{k}}}{(2\pi)^{3}}\frac{d^{3}{\vec{q}}}{(2\pi)^{3}}\frac{d^{3}{\vec{p}}}{(2\pi)^{3}}{e^{i\vec{k}\cdot(\vec{x}-\vec{y})}}$
$\displaystyle\times\bigg{[}\theta(\eta_{1}-\eta^{\prime})G_{k}^{C}(\eta_{1},\eta^{\prime})G_{k}^{<}(\eta_{2},\eta^{\prime\prime})$
$\displaystyle~{}~{}~{}+\theta(\eta_{2}-\eta^{\prime})G_{k}^{C}(\eta_{2},\eta^{\prime})G_{k}^{>}(\eta_{1},\eta^{\prime\prime})\bigg{]}$
$\displaystyle\times\bigg{[}G_{p}^{>}(\eta^{\prime},\eta^{\prime\prime})G_{q}^{>}(\eta^{\prime},\eta^{\prime\prime})+G_{p}^{>}(\eta^{\prime},\eta^{\prime\prime})G_{q}^{<}(\eta^{\prime},\eta^{\prime\prime})$
$\displaystyle~{}~{}~{}+G_{p}^{<}(\eta^{\prime},\eta^{\prime\prime})G_{q}^{<}(\eta^{\prime},\eta^{\prime\prime})\bigg{]}G_{\ell}^{C}(\eta^{\prime},\eta^{\prime\prime})\\!\\!\\!\\!\\!\\!$
$\displaystyle+\int_{-\infty}^{\eta_{1}}{d\eta^{\prime}}a^{4}(\eta^{\prime})\int_{-\infty}^{\eta_{2}}{d\eta^{\prime\prime}}a^{4}(\eta^{\prime\prime})$
$\displaystyle\times\int\frac{d^{3}{\vec{q}}}{(2\pi)^{3}}\frac{d^{3}{\vec{p}}}{(2\pi)^{3}}\frac{d^{3}{\vec{\ell}}}{(2\pi)^{3}}{e^{i(\vec{\ell}+\vec{p}+\vec{q})\cdot(\vec{x}-\vec{y})}}G_{k}^{C}(\eta_{1},\eta^{\prime})$
$\displaystyle\times
G_{k}^{C}(\eta_{2},\eta^{\prime\prime})G_{\ell}^{>}(\eta^{\prime},\eta^{\prime\prime})G_{p}^{>}(\eta^{\prime},\eta^{\prime\prime})G_{q}^{>}(\eta^{\prime},\eta^{\prime\prime})\bigg{\\}}\;,$

where on the seventh line $\ell=|\vec{k}+\vec{p}+\vec{q}\rvert$, and on the
last two lines $k=|\vec{\ell}+\vec{p}+\vec{q}\rvert$. Because in the long-
wavelength limit $G^{C}$-s are momentum-independent, the momentum integrals
decouple in this limit and one obtains

$\displaystyle\\!\\!\\!\\!\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{2},L}^{(3)}={\lambda^{2}H^{2}\over
12}(\eta_{1}\eta_{2})^{u}$ $\displaystyle~{}~{}\times\Bigg{\\{}{3\over
r^{2u}}\bigg{[}\int_{-\infty}^{\eta_{1}}{d\eta^{\prime}\over\eta^{\prime}}(-\eta^{\prime})^{2u}\int_{-\infty}^{\eta^{\prime}}{d\eta^{\prime\prime}\over\eta^{\prime\prime}}(-\eta^{\prime\prime})^{2u}$
$\displaystyle\quad\quad\quad\quad\quad+\int_{-\infty}^{\eta_{2}}{d\eta^{\prime}\over\eta^{\prime}}(-\eta^{\prime})^{2u}\int_{-\infty}^{\eta^{\prime}}{d\eta^{\prime\prime}\over\eta^{\prime\prime}}(-\eta^{\prime\prime})^{2u}\bigg{]}$
$\displaystyle\quad\quad\quad\times\int_{0}^{-\epsilon r/\eta_{\rm
m}}\\!\\!{\sin w\over
w^{2-2u}}{dw}\bigg{(}\int_{0}^{-\epsilon/\eta^{\prime\prime}}{dp\over
p^{1-2u}}\bigg{)}^{2}$ $\displaystyle\quad\quad\quad+{1\over
r^{6u}}\int_{-\infty}^{\eta_{1}}{d\eta^{\prime}\over\eta^{\prime}}(-\eta^{\prime})^{2u}\int_{-\infty}^{\eta_{2}}{d\eta^{\prime\prime}\over\eta^{\prime\prime}}(-\eta^{\prime\prime})^{2u}$
$\displaystyle\quad\quad\quad\quad\times\bigg{(}\int_{0}^{-\epsilon
r/\eta^{\prime}_{\rm m}}{\sin w\over w^{2-2u}}{dw}\bigg{)}^{3}\Bigg{\\}}\;,$
(58)

where $\eta_{\rm m}$ is the earliest time that accompanies the momentum
$\vec{k}$: $\eta_{\rm m}\equiv{\rm
min}(\eta_{1},\eta_{2},\eta^{\prime},\eta^{\prime\prime})$, and
$\eta^{\prime}_{\rm m}$ is the earliest time that accompanies the momenta
$\vec{\ell}$, $\vec{p}$ and $\vec{q}$: $\eta^{\prime}_{\rm m}\equiv{\rm
min}(\eta^{\prime},\eta^{\prime\prime})$. For $Z>1-(2\epsilon^{2})^{-1}$,
evaluating (58) gives

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{2},L}^{(3)}$
(59) $\displaystyle=$ $\displaystyle{\lambda^{2}H^{2}\over
1024\pi^{6}u^{5}}\bigg{(}1+2uH|t_{1}-t_{2}\rvert\bigg{)}e^{-uH|t_{1}-t_{2}\rvert}$
$\displaystyle+{\lambda^{2}H^{2}\over
3072\pi^{6}u^{5}}e^{-3uH|t_{1}-t_{2}\rvert}\;.$

For coinciding spacetime points, (59) reduces to

$\displaystyle\langle\phi^{2}(\vec{x},t)\rangle_{\lambda^{2},L}^{(3)}={81\lambda^{2}H^{12}\over
256\pi^{6}m^{10}}\;.$ (60)

In the opposite regime, $Z<1-(2\epsilon^{2})^{-1}$, Eq. (58) yields

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{\lambda^{2},L}^{(3)}$
(61) $\displaystyle=$ $\displaystyle{\lambda^{2}H^{2}\over
1024\pi^{6}u^{5}}\left(r^{2}H^{2}e^{H(t_{1}+t_{2})}\right)^{-u}$
$\displaystyle\quad\quad\quad~{}\times\left\\{1+2u\ln\left({r^{2}H^{2}e^{H(t_{1}+t_{2})}}\right)\right\\}$
$\displaystyle+{\lambda^{2}H^{2}\over
3072\pi^{6}u^{5}}\left(r^{2}H^{2}e^{H(t_{1}+t_{2})}\right)^{-3u}\;.$

Could the ultraviolet divergences and renormalization change the results
obtained in this section? To answer this question, let us look, for example,
at one of the terms in (57) rewritten in coordinate space:

$\displaystyle\\!\\!\\!\\!\int_{-\infty}^{\eta_{2}}\\!{d\eta^{\prime}}a^{4}(\eta^{\prime})\\!\int_{-\infty}^{\eta^{\prime}}\\!{d\eta^{\prime\prime}}a^{4}(\eta^{\prime\prime})\\!\int\\!d^{3}\vec{v}\\!\\!\int\\!d^{3}\vec{w}\,G^{C}(\vec{y},t_{2};\vec{v},t^{\prime})$
$\displaystyle\times
G^{>}(\vec{x},t_{1};\vec{w},t^{\prime\prime})G^{C}(\vec{v},t^{\prime};\vec{w},t^{\prime\prime})\Big{[}G^{>}(\vec{v},t^{\prime};\vec{w},t^{\prime\prime})\Big{]}^{2}\;.$

Ultraviolet divergences in this expression occur when
$(\vec{v},t^{\prime})=(\vec{w},t^{\prime\prime})$. Recall that
$G^{>}\equiv\langle\phi(\vec{v},t^{\prime})\phi(\vec{w},t^{\prime\prime})\rangle_{\lambda^{0}}$;
hence, around $Z=1$ it is given by (LABEL:UV) for spacelike related points and
by

$\displaystyle
G^{>}(\vec{v},t^{\prime};\vec{w},t^{\prime\prime})\to{H^{2}\over
8\pi^{2}}\bigg{(}{1\over u}$ $\displaystyle-$ $\displaystyle 2-{1\over Z-1\pm
i\epsilon}$ $\displaystyle-$ $\displaystyle\ln{(Z-1)\over 2}\pm
i\pi+\mathcal{O}(u)\bigg{)}$

for timelike related points. As shown in Sec. 3.1, the first term in this
expression arises solely from modes with small physical momenta.

To write the function $G^{C}\equiv G^{>}-G^{<}=2\,{\rm Im}(G^{>})$, we can use
the identity

$\displaystyle{1\over x\pm i\epsilon}=P{1\over x}\mp i\pi\delta(x)\;,$

where $P$ means the principal value. For timelike related points,

$\displaystyle
G^{C}(\vec{v},t^{\prime};\vec{w},t^{\prime\prime})\to{H^{2}\over
4\pi^{2}}\Big{(}\pm i\pi\delta(Z-1)\pm i\pi\Big{)}\;,$

and for spacelike related points it is equal to zero.

In the product
$G^{C}(\vec{v},t^{\prime};\vec{w},t^{\prime\prime})\Big{[}G^{>}(\vec{v},t^{\prime};\vec{w},t^{\prime\prime})\Big{]}^{2}$,
the leading term in $u$ that does not produce UV divergences is of the order
$H^{6}/u^{2}$ and derives entirely from the long-wavelength modes. On the
other hand, the leading term containing UV divergences is of the order
$H^{6}/u$. Therefore, the renormalization does not affect the value of the
two-point function at leading order in small mass.

Now, taking this into account and combining (14), (43), (51), (55) and (60),
we get the final expression for the average value of $\langle\phi^{2}\rangle$
at coinciding spacetime points in the equilibrium Bunch-Davies state:

$\langle\phi^{2}\rangle={3H^{4}\over 8\pi^{2}m^{2}}-{27\lambda H^{8}\over
64\pi^{4}m^{6}}+{81\lambda^{2}H^{12}\over
64\pi^{6}m^{10}}+\mathcal{O}(\lambda^{3})\;.$ (62)

We have not obtained any instability at the two-loop level conjectured in
Polyakov:2012uc . This result is in qualitative agreement with general
considerations in Marolf:2010nz ; Hollands:2010pr ; Gautier:2015pca . Note,
however, that the constancy of $\langle\phi^{2}\rangle$ does not mean that
there is no particle creation in this state. Although the notion of particles
is not well defined in the long-wavelength (super-Hubble, anti-WKB) regime and
$\langle\phi^{2}\rangle$ includes contributions from both created particles
and vacuum polarization, from the fact that the effective number of particles
in each Fourier mode ${\vec{k}}$ grows with time in the super-Hubble regime,
one may say that strong particle creation occurs in the Bunch-Davies state,
but this growth is exactly compensated by the universe expansion. This also
shows that it is misleading to call the Bunch-Davies state a vacuum, at least
if $m^{2}\ll H^{2}$. However, it is certainly the equilibrium state.

Our expression (62) agrees with the results in Gat-Ser ; Garb-Rig-Zhu , which
were also obtained by quantum field theory calculation.

##  4 Comparison with the Hartree-Fock approximation

Let us follow paper Star-Yoko  and write the equation of motion for the scalar
field with the action (1),

$\displaystyle\Box\phi=-m^{2}\phi-\lambda\phi^{3}\;,$

where $\Box\equiv g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}$ is the covariant
d’Alembertian. Multiplying both sides of this equation by $\phi$, using the
identity

$\displaystyle\phi\Box\phi={1\over
2}\Box{\phi^{2}}-\partial_{\mu}\phi\partial^{\mu}\phi\;,$

and taking expectation values of the field operators results in the following
equation

$\frac{1}{2}\Box\langle\phi^{2}\rangle-\langle\partial_{\mu}\phi\,\partial^{\mu}\phi\rangle=-m^{2}\langle\phi^{2}\rangle-\lambda\langle\phi^{4}\rangle\;.$

If we use the Hartree-Fock (Gaussian) approximation,
$\langle\phi^{4}\rangle=3\langle\phi^{2}\rangle^{2}$, this equation can be
rewritten as

$\frac{1}{2}\Box\langle\phi^{2}\rangle-\langle\partial_{\mu}\phi\,\partial^{\mu}\phi\rangle=-m^{2}\langle\phi^{2}\rangle-3\lambda\langle\phi^{2}\rangle^{2}\;.$
(63)

In the case of a massless noninteracting field, the long-wavelength part of
$\langle\phi^{2}\rangle$ is given by

$\displaystyle\langle\phi^{2}(\vec{x},t)\rangle_{L}$ $\displaystyle=$
$\displaystyle\int^{-\epsilon/\eta}\frac{d^{3}{\vec{k}}}{(2\pi)^{3}}\phi_{k}(\eta)\phi_{k}^{*}(\eta)={H^{2}\over
4\pi^{2}}\int_{\kappa}^{-\epsilon/\eta}{dk\over k}$ (64) $\displaystyle=$
$\displaystyle-{H^{2}\over
4\pi^{2}}\ln\left({-\kappa\eta\over\epsilon}\right)\;,$

where Eq. (10) with $u=0$ was used to express the mode functions. We
introduced an infrared cutoff $\kappa$ for the comoving momentum $k$ since the
integral is divergent at $k=0$. In terms of cosmic time (64) becomes

$\displaystyle\langle\phi^{2}\rangle_{L}={H^{3}\over 4\pi^{2}}(t-t_{0})\;,$
(65)

where $t_{0}\equiv(1/H)\ln(\kappa/H\epsilon)$; thus, it grows linearly with
time Allen:1985ux ; Allen:1987tz ; Starobinsky:1982ee ; Linde:1982uu ;
Vilenkin:1982wt . This means that the only term that survives in the
d’Alembertian operator is the first time derivative:

$\displaystyle\Box\langle\phi^{2}\rangle_{L}=3H{\partial\over\partial
t}\langle\phi^{2}\rangle_{L}\;.$ (66)

The second term on the left side of Eq. (63) contains two derivatives and,
therefore, is a constant proportional to $\epsilon^{2}$,

$\langle\partial_{\mu}\phi\,\partial^{\mu}\phi\rangle_{L}\propto\epsilon^{2}\;;$
(67)

hence, it can be neglected.

When $\lambda$ and $m$ are small, expressions (66) and (67) are still
approximately correct, so Eq. (63) reduces to

${\partial\over\partial
t}\langle\phi^{2}\rangle_{L}=\frac{H^{3}}{4\pi^{2}}-{2m^{2}\over
3H}\langle\phi^{2}\rangle_{L}-\frac{2\lambda}{H}\langle\phi^{2}\rangle_{L}^{2}\;,$
(68)

where the first term on the right side is needed to recover (65) in the limit
of zero mass and coupling. As $t\to\infty$, all of the solutions to this
equation approach an equilibrium value that satisfies

$\frac{H^{3}}{4\pi^{2}}-{2m^{2}\over
3H}\langle\phi^{2}\rangle_{L}-\frac{2\lambda}{H}\langle\phi^{2}\rangle_{L}^{2}=0\;.$
(69)

For $\lambda=0$, we have

$\langle\phi^{2}\rangle_{L}={3H^{4}\over 8\pi^{2}m^{2}}\;.$ (70)

When $\lambda\neq 0$, Eq (69) gives

$\langle\phi^{2}\rangle_{L}={m^{2}\over 6\lambda}\left(\sqrt{1+{9\lambda
H^{4}\over 2\pi^{2}m^{4}}}-1\right);$ (71)

we chose the root that coincides with (70) in the limit $\lambda\to 0$.
Assuming that $\lambda H^{4}/m^{4}\ll 1$, and expanding (71) yields

$\langle\phi^{2}\rangle_{L}={3H^{4}\over 8\pi^{2}m^{2}}-{27\lambda H^{8}\over
64\pi^{4}m^{6}}+{243\lambda^{2}H^{12}\over
256\pi^{6}m^{10}}+\mathcal{O}(\lambda^{3})\;.$ (72)

Comparing this expansion with the two-loop results of the previous section
(with (62) in particular), we see that they match at zeroth- and first-order
in $\lambda$, but there is a mismatch at second order: the $\lambda^{2}$-term
in (72) omits the contribution of the sunset diagram and is equal to the sum
of (51) and (55). As a result, in the Hartree-Fock approximation this term is
$4/3$ times smaller. Also, it has to be concluded that Eq. (71) resumes all
cactus type diagrams of the perturbation theory.

##  5 Comparison with the stochastic approach

The stochastic approach argues Star ; Star-Yoko  that the behavior of the
long-wavelength part of the quantum field $\phi(\vec{x},t)$ in de Sitter
spacetime can be modeled by an auxiliary classical stochastic variable
$\varphi$ with a probability distribution $\rho(\varphi,t)$ that satisfies the
Fokker-Planck equation

${\partial\rho\over\partial t}={H^{3}\over
8\pi^{2}}{\partial^{2}\rho\over\partial\varphi^{2}}+{1\over
3H}{\partial\over\partial\varphi}\biggl{(}{\partial
V\over\partial\varphi}\rho(t,\varphi)\biggr{)}\;;$ (73)

that is, the expectation value of any quantity constructed from the long-
wavelength part of $\phi(\vec{x},t)$ is equal to the expectation value of the
same quantity constructed from the variable $\varphi$. $V(\varphi)$ is a
function of the stochastic variable and has the same functional form as the
corresponding potential of the quantum field theory: for the theory with the
action (1) this potential has the form

$V(\varphi)={1\over 2}{m^{2}}\varphi^{2}+{\lambda\over 4}\varphi^{4}\;.$ (74)

At late times any solution of the equation (73) approaches the static
equilibrium solution

$\rho_{\rm
eq}(\varphi)=N^{-1}\exp\left(-\frac{8\pi^{2}}{3H^{4}}V(\varphi)\right),$ (75)

where $N$ is the normalization fixed by the condition

$\displaystyle\int_{-\infty}^{\infty}\rho_{\rm eq}(\varphi)\,d\varphi=1\;.$

In particular, with the potential (74) we can calculate this normalization
explicitly

$\displaystyle N$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\exp\left[-{8\pi^{2}\over
3H^{4}}\left({\lambda\varphi^{4}\over 4}+{m^{2}\varphi^{2}\over
2}\right)\right]d\varphi$ (76) $\displaystyle=$
$\displaystyle{m\over{\sqrt{2\lambda}}}\exp(z)\mathcal{K}_{1\over 4}(z)\;,$

where $\mathcal{K}_{1\over 4}(z)$ is a modified Bessel function of the second
kind, and $z\equiv{\pi^{2}m^{4}\over 3\lambda H^{4}}$. The distribution (75)
can be used to compute the expectation value of $\varphi^{2}$:

$\displaystyle\langle\varphi^{2}\rangle$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\varphi^{2}\rho_{\rm
eq}(\varphi)d\varphi$ $\displaystyle=$ $\displaystyle
N^{-1}\int_{-\infty}^{\infty}\varphi^{2}\exp\left[-{8\pi^{2}\over
3H^{4}}\left({\lambda\varphi^{4}\over 4}+{m^{2}\varphi^{2}\over
2}\right)\right]d\varphi\;.$

The integral above can be expressed in the following way:

$\displaystyle\\!\\!\\!\\!\int_{-\infty}^{\infty}\\!\\!\varphi^{2}\exp\left[-{8\pi^{2}\over
3H^{4}}\left({\lambda\varphi^{4}\over 4}+{m^{2}\varphi^{2}\over
2}\right)\right]$ $\displaystyle={\pi m^{3}\over
4\lambda^{3/2}}\exp(z)\left[{\cal I}_{1\over 4}(z)-{\cal I}_{-{1\over
4}}(z)+{\cal I}_{5\over 4}(z)-{\cal I}_{3\over 4}(z)\right]$
$\displaystyle~{}~{}+{3H^{4}\over 8\pi m}{\exp(z)\over\lambda^{1/2}}{\cal
I}_{1\over 4}(z)\,,$ (78)

where $\mathcal{I}_{\nu}(z)$ are modified Bessel functions of the first kind.
Combining the right sides of (76) and (78), and using the connection formula

$\displaystyle\mathcal{K}_{\nu}(z)={\pi\over
2}{\mathcal{I}_{-\nu}(z)-\mathcal{I}_{\nu}(z)\over\sin(\nu\pi)}$

and the recurrence relation

$\displaystyle\mathcal{I}_{\nu-1}(z)-\mathcal{I}_{\nu+1}(z)={2\nu\over
z}\mathcal{I}_{\nu}(z),$

we obtain a fairly concise expression for (LABEL:two)

$\displaystyle\langle\varphi^{2}\rangle={m^{2}\over
2\lambda}{\mathcal{K}_{3\over 4}(z)\over\mathcal{K}_{1\over 4}(z)}-{m^{2}\over
2\lambda}.$ (79)

Interestingly, this result is in agreement with the expression obtained in
Beneke:2012kn  for the amplitude of the zero mode in Euclidean de Sitter
space.

Expanding (79) in the limit $\lambda H^{4}/m^{4}\ll 1$ (which corresponds to
$z\gg 1$) gives

$\langle\varphi^{2}\rangle={3H^{4}\over 8\pi^{2}m^{2}}-{27\lambda H^{8}\over
64\pi^{4}m^{6}}+{81\lambda^{2}H^{12}\over
64\pi^{6}m^{10}}-\frac{24057\lambda^{3}H^{16}}{4096\pi^{8}m^{14}}+\mathcal{O}(\lambda^{4})\;.$
(80)

As we can see, this result coincides exactly with the result (62) of the
quantum field theory calculations of Sec. 3, and unlike the Hartree-Fock
approximation (72), it includes the contribution of the sunset diagram; that
is, the $\lambda^{2}$-term is equal to the sum of (51), (55) and (60). In Eq.
(80) we have also given the expression for the term proportional to
$\lambda^{3}$. The corresponding expression obtained by the field-theoretical
methods will be presented in the future work future .

Note also that in the opposite limit $z\to 0$, the known expression for the
massless self-interacting field arises Star-Yoko  that would require summation
of all loops in the standard field-theoretic approach:

$\langle\varphi^{2}\rangle=\sqrt{\frac{3}{2\pi^{2}}}\,\frac{\Gamma(3/4)}{\Gamma(1/4)}\,\frac{H^{2}}{\sqrt{\lambda}}\approx
0.132\,\frac{H^{2}}{\sqrt{\lambda}}\,.$ (81)

In this case, the Hartree-Fock approximation (71) produces the result that is
only 1.17 times smaller than (81) (see Star-Yoko ).

Interestingly, by differentiating the formula (76) with respect to $m^{2}$ and
using the recurrence relations for the modified Bessel functions of the second
kind (see e.g. Grad ),

$z\frac{d}{dz}\mathcal{K}_{\nu}(z)+\nu\mathcal{K}_{\nu}(z)=-z\mathcal{K}_{\nu-1}(z),$

and the fact that $\mathcal{K}_{-\nu}(z)=\mathcal{K}_{\nu}(z)$, one can
express an arbitrary average value $\langle\varphi^{2n}\rangle$ in terms of
the ratio $\mathcal{K}_{3\over 4}(z)/\mathcal{K}_{1\over 4}(z)$ . For example,
for $\langle\varphi^{4}\rangle$, we obtain

$\langle\varphi^{4}\rangle=\frac{3H^{4}}{8\pi^{2}\lambda}+\frac{m^{4}}{2\lambda^{2}}\left(1-\frac{\mathcal{K}_{3\over
4}(z)}{\mathcal{K}_{1\over 4}(z)}\right).$

As discussed in Star-Yoko , the long-wavelength two-point function of
$\phi(\vec{x},t)$ too can be calculated by using the classical stochastic
variable $\varphi$: if the points $(\vec{x},t_{1})$ and $(\vec{y},t_{2})$ are
timelike or lightlike related, this correlation function is given by

$\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{L}=\langle\varphi(t_{1})\varphi(t_{2})\rangle\;.$
(82)

If the correlation function $\langle\varphi(t_{1})\varphi(t_{2})\rangle$
depends only on the absolute value of the time difference
$T\equiv|t_{1}-t_{2}\rvert$, it can be expressed as Star-Yoko ; Risken

$\langle\varphi(t_{1})\varphi(t_{2})\rangle=\int_{-\infty}^{\infty}\varphi\,\Xi(\varphi,T)d\varphi\;,$
(83)

where the function $\Xi(\varphi,T)$ satisfies the Fokker-Planck equation (73),

${\partial\Xi\over\partial T}={H^{3}\over
8\pi^{2}}{\partial^{2}\Xi\over\partial\varphi^{2}}+{1\over
3H}{\partial\over\partial\varphi}\biggl{(}{\partial
V\over\partial\varphi}\Xi(\varphi,T)\biggr{)}\;,$ (84)

with the initial condition

$\Xi(\varphi,0)=\varphi\rho_{\rm eq}(\varphi)\;,$ (85)

and $\rho_{\rm eq}(\varphi)$ is the equilibrium solution (75). Derivatives of
$\langle\varphi(t_{1})\varphi(t_{2})\rangle$ at $T=0$ can be computed by using
Eqs. (83)–(85):

${\partial\over\partial
T}\langle\varphi(t_{1})\varphi(t_{2})\rangle\bigg{|}_{T=0}=-{H^{3}\over
8\pi^{2}}\;,$ (86) ${\partial^{2}\over\partial
T^{2}}\langle\varphi(t_{1})\varphi(t_{2})\rangle\bigg{|}_{T=0}={H^{2}\over
24\pi^{2}}\Big{(}3\lambda\langle\varphi^{2}\rangle+m^{2}\Big{)}\;,$ (87)

and so on. It is easy to confirm that the $T$-derivatives of the two-point
correlation function computed in Sec. 3 (for $Z>1-{1\over 2\epsilon^{2}}$
case) satisfy these equalities as well.

##  6 Exponentiation of the perturbative series

Let us put together the calculations of Sec. 3: up to two loops, for
$Z>1-(2\epsilon^{2})^{-1}$ the long-wavelentgth two-point correlation function
is given by

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{L}$
$\displaystyle=$ $\displaystyle{H^{2}\over 8\pi^{2}u}e^{-uHT}-{\lambda
H^{2}\over 64\pi^{4}u^{3}}\bigg{(}1+uHT\bigg{)}e^{-uHT}$
$\displaystyle+{\lambda^{2}H^{2}\over 512\pi^{6}u^{5}}\bigg{(}1+uHT+{1\over
2}u^{2}H^{2}{T}^{2}\bigg{)}e^{-uHT}$ $\displaystyle+{\lambda^{2}H^{2}\over
512\pi^{6}u^{5}}\bigg{(}1+uHT\bigg{)}e^{-uHT}$
$\displaystyle+{\lambda^{2}H^{2}\over
1024\pi^{6}u^{5}}\bigg{(}1+2uHT\bigg{)}e^{-uHT}$
$\displaystyle+{\lambda^{2}H^{2}\over
3072\pi^{6}u^{5}}e^{-3uHT}+\mathcal{O}(\lambda^{3})\;,$

where $T\equiv|t_{1}-t_{2}\rvert$. It will be convenient to reorganize this
expression in the following way:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{L}$
(88) $\displaystyle=$ $\displaystyle{H^{2}\over
8\pi^{2}u}\left(1-{\lambda\over 8\pi^{2}u^{2}}+{5\lambda^{2}\over
128\pi^{4}u^{4}}\right)e^{-uHT}$ $\displaystyle-{\lambda H^{3}T\over
64\pi^{4}u^{2}}\left(1-{3\lambda\over 8\pi^{2}u^{2}}\right)e^{-uHT}$
$\displaystyle+{\lambda^{2}H^{4}{T}^{2}\over
1024\pi^{6}u^{3}}e^{-uHT}+{\lambda^{2}H^{2}\over
3072\pi^{6}u^{5}}e^{-3uHT}+\mathcal{O}(\lambda^{3})$ $\displaystyle=$
$\displaystyle{H^{2}\over 8\pi^{2}u}\left(1-{\lambda\over
8\pi^{2}u^{2}}+{5\lambda^{2}\over
128\pi^{4}u^{4}}+\mathcal{O}(\lambda^{3})\right)$
$\displaystyle\times\left[1-C_{1}{\lambda HT\over 8\pi^{2}u}+{C_{2}\over
2}\left({\lambda HT\over
8\pi^{2}u}\right)^{2}+\mathcal{O}(\lambda^{3})\right]e^{-uHT}$
$\displaystyle+{\lambda^{2}H^{2}\over 3072\pi^{6}u^{5}}e^{-3uHT}+\dots\;,$

where $C_{1}\equiv(1-{3\lambda\over 8\pi^{2}u^{2}})(1-{\lambda\over
8\pi^{2}u^{2}}+{5\lambda^{2}\over 128\pi^{4}u^{4}})^{-1}$ and
$C_{2}\equiv(1-{\lambda\over 8\pi^{2}u^{2}}+{5\lambda^{2}\over
128\pi^{4}u^{4}})^{-1}$. Upon expanding $C_{1}$ and $C_{2}$ in powers of
$\lambda$, the second term in the square brackets of (88) becomes

$\displaystyle-C_{1}{\lambda HT\over 8\pi^{2}u}=-{\lambda HT\over
8\pi^{2}u}+{\lambda^{2}HT\over 32\pi^{4}u^{3}}+\mathcal{O}(\lambda^{3})\;,$

and the third term

$\displaystyle{C_{2}\over 2}\left({\lambda HT\over
8\pi^{2}u}\right)^{2}={1\over 2}\left({\lambda HT\over
8\pi^{2}u}\right)^{2}+\mathcal{O}(\lambda^{3}).$

Hence the expression in square brackets reduces to

$1-{\lambda HT\over 8\pi^{2}u}+{\lambda^{2}HT\over 32\pi^{4}u^{3}}+{1\over
2}\left({\lambda HT\over 8\pi^{2}u}\right)^{2}+\mathcal{O}(\lambda^{3})\;.$

To second order in $\lambda$, this matches with the first three terms in the
Taylor series of the exponential function

$\displaystyle\exp\left[{-{\lambda HT\over 8\pi^{2}u}\left(1-{\lambda\over
4\pi^{2}u^{2}}\right)}\right],$ (89)

so it is plausible that an infinite series of diagrams may be resummed into
this exponent. With this assumption, we arrive at

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{L}$
(90) $\displaystyle=$ $\displaystyle{H^{2}\over
8\pi^{2}u}\left(1-{\lambda\over 8\pi^{2}u^{2}}+{5\lambda^{2}\over
128\pi^{4}u^{4}}+\mathcal{O}(\lambda^{3})\right)$
$\displaystyle\times\exp\bigg{[}{-uHT\left(1+{\lambda\over
8\pi^{2}u^{2}}-{\lambda^{2}\over
32\pi^{4}u^{4}}+\mathcal{O}(\lambda^{3})\right)}\bigg{]}$
$\displaystyle+{\lambda^{2}H^{2}\over
3072\pi^{6}u^{5}}\exp\Big{[}{-3uHT}\Big{]}+\dots\;.$

Similarly, for $Z<1-(2\epsilon^{2})^{-1}$ we obtain

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{L}$
(91) $\displaystyle=$ $\displaystyle{H^{2}\over
8\pi^{2}u}\left(1-{\lambda\over 8\pi^{2}u^{2}}+{5\lambda^{2}\over
128\pi^{4}u^{4}}+\mathcal{O}(\lambda^{3})\right)$
$\displaystyle\times\left(r^{2}H^{2}e^{H(t_{1}+t_{2})}\right)^{-u\left(1+{\lambda\over
8\pi^{2}u^{2}}-{\lambda^{2}\over
32\pi^{4}u^{4}}+\mathcal{O}(\lambda^{3})\right)}$
$\displaystyle+{\lambda^{2}H^{2}\over
3072\pi^{6}u^{5}}\left(r^{2}H^{2}e^{H(t_{1}+t_{2})}\right)^{-3u}+\dots\;,$

and it is evident that in this regime the equal-time correlation function
depends only on the physical spatial distance $R\equiv re^{Ht}$.

From the expressions (90) and (91), we see that the perturbative corrections
don’t change the long-wavelength part of the commutator: just as in the free
theory case, it is equal to zero both for timelike and spacelike related
points.

When $|Z|\gg 1$, both (90) and (91) can be expressed in a de Sitter invariant
form. As was noted in Sec. 3.1, for large timelike separations $2Z\approx
e^{HT}$, and for large spacelike separations
$2Z\approx-H^{2}r^{2}e^{H(t_{1}+t_{2})}$. Therefore,

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{L}$
(92) $\displaystyle=$ $\displaystyle{H^{2}\over
8\pi^{2}u}\left(1-{\lambda\over 8\pi^{2}u^{2}}+{5\lambda^{2}\over
128\pi^{4}u^{4}}+\mathcal{O}(\lambda^{3})\right)$
$\displaystyle\times|Z|^{-u\left(1+{\lambda\over
8\pi^{2}u^{2}}-{\lambda^{2}\over
32\pi^{4}u^{4}}+\mathcal{O}(\lambda^{3})\right)}$
$\displaystyle+{\lambda^{2}H^{2}\over 3072\pi^{6}u^{5}}|Z|^{-3u}+\dots\;$

both for large timelike and large spacelike separations. This expression is
also correct for coinciding spacetime points, since in this case $Z=1$ and
(92) gives the same result as (90) with $T$ set equal to zero. In Gat-Ser ,
momentum representation of the unequal-time correlation function is calculated
by using the Schwinger-Dyson equations to resum the infinite series of self-
energy insertions. A comparison of its results with (90)–(92) is presented in
Appendix B.

As was shown in Star-Yoko , correlation functions deduced by the stochastic
approach have a form of an infinite sum of exponentials; e.g., if points are
timelike related, the two-point function is given by

$\langle\phi(\vec{x},t_{1})\phi(\vec{y},t_{2})\rangle_{L}=\langle\varphi(t_{1})\varphi(t_{2})\rangle=N\sum_{n=1}^{\infty}A_{n}^{2}e^{-\Lambda_{n}T},$

where $N$ is the normalization constant (76). By comparing this with Eq. (90),
we can identify the first three coefficients $A_{n}$ and $\Lambda_{n}$:

$\displaystyle NA_{1}^{2}$ $\displaystyle=$ $\displaystyle{H^{2}\over
8\pi^{2}u}\left(1-{\lambda\over 8\pi^{2}u^{2}}+{5\lambda^{2}\over
128\pi^{4}u^{4}}+\mathcal{O}(\lambda^{3})\right)\;,$ $\displaystyle A_{2}$
$\displaystyle=$ $\displaystyle 0\,,\,NA_{3}^{2}={\lambda^{2}H^{2}\over
3072\pi^{6}u^{5}}+\mathcal{O}(\lambda^{3})\;,$ $\displaystyle\Lambda_{1}$
$\displaystyle=$ $\displaystyle\left(1+{\lambda\over
8\pi^{2}u^{2}}-{\lambda^{2}\over
32\pi^{4}u^{4}}+\mathcal{O}(\lambda^{3})\right)uH\;,$
$\displaystyle\Lambda_{3}$ $\displaystyle=$ $\displaystyle
3uH+\mathcal{O}(\lambda)\;.$

These coefficients coincide with the results of Mark-et ; Mor-Ser  where they
are derived using some versions of the stochastic approach (e.g., Mor-Ser
uses the formulation of the stochastic theory in terms of a supersymmetric
one-dimensional field theory).

As $T\to\infty$, the two-point function (90) decays with the characteristic
correlation time

$T_{c}\sim{1\over uH}={3H\over m^{2}}\gg{1\over H}\;.$

Similarly, it follows from (91) that as $R\to\infty$, the equal-time
correlation function decays with the characteristic correlation length

$R_{c}\sim{1\over H}\exp\left({3H^{2}\over 2m^{2}}\right)\;.$

The behavior of (91) differs from a much faster exponential decay of the
equal-time correlation function in flat spacetime:
$\langle\phi(\vec{x},t)\phi(\vec{y},t)\rangle_{\rm
flat}\sim\sqrt{m/r^{3}}e^{-mr}$ as $r\to\infty$.

##  7 Conclusions

In Minkowski spacetime, quantum fluctuations decay quickly with scale and are
significant only at short distances. The situation is different in de Sitter
spacetime. Here the exponential expansion stretches short-wavelength modes to
super-Hubble scales without decreasing its amplitude too much. On super-Hubble
scales, the amplitude of fluctuations is scale-independent for the massless
minimally coupled scalar field, and only weakly decays with scale when
$0<m^{2}\ll H^{2}$.

In this article we have calculated—up to two loops—the long-wavelength two-
point function for a scalar theory with a small mass and a quartic
interaction. It has been shown that it is de Sitter invariant for coinciding
points as well as at large spacelike and large timelike separations. We have
also demonstrated that the commutator of the long-wavelength part of the field
is equal to zero both at the free theory level and with the perturbative
corrections.

Our results are in agreement with Starobinsky’s stochastic approach in which
the coarse-grained quantum field is equivalent to a classical stochastic
quantity.

###### Acknowledgements.

The work of A. K. and A. S. was partially supported by the RFBR Grant No.
20-02-00411. The work of T. V. was partially supported by the research grant
“The Dark Universe: A Synergic Multimessenger Approach” No. 2017X7X85K under
the program PRIN 2017 funded by the Ministero dell’Istruzione, delle
Università e della Ricerca (MIUR).

##  Appendix A Geodesics in de Sitter spacetime

De Sitter spacetime can be represented as the four-dimensional hyperboloid

$(X^{1})^{2}+(X^{2})^{2}+(X^{3})^{2}+(X^{4})^{2}-(X^{0})^{2}=H^{-2}$ (93)

embedded in five-dimensional Minkowski spacetime. In this appendix we set
$H=1$.

For two points in de Sitter spacetime, the only invariant quantity associated
with these points is the function

$Z(X,Y)\equiv X^{1}Y^{1}+X^{2}Y^{2}+X^{3}Y^{3}+X^{4}Y^{4}-X^{0}Y^{0}\;.$ (94)

Let us consider the spacetime interval between these points:

$\displaystyle s^{2}(X,Y)$ $\displaystyle\equiv$
$\displaystyle(X^{0}-Y^{0})^{2}-(X^{1}-Y^{1})^{2}-(X^{2}-Y^{2})^{2}$ (95)
$\displaystyle-(X^{3}-Y^{3})^{2}-(X^{4}-Y^{4})^{2}.$

Using the hyperboloid equation (93) (with $H=1$) and the definition of $Z$
from Eq. (94), we obtain

$Z(X,Y)=1+{1\over 2}s^{2}(X,Y).$ (96)

Thus, if the interval between $X$ and $Y$ is timelike, then $Z(X,Y)>1$; if the
interval is spacelike, then $Z(X,Y)<1$, and if the interval is lightlike, then
$Z(X,Y)=1$. Note the if $\overline{Y}$ is the antipodal point of $Y$ (that is,
$\overline{Y}=-Y$), then $Z(X,\overline{Y})=-Z(X,Y)$; hence, if $Z(X,Y)<-1$,
the interval between $X$ and $\overline{Y}$ is timelike.

Let us show that the geodesics in de Sitter spacetime are intersections of the
hyperboloid with hyperplanes passing through the origin of the ambient
spacetime just like the geodesics on the sphere are the great circles. To
simplify notations, we will consider two-dimensional de Sitter spacetime
embedded into 3-dimensional flat spacetime. In terms of the spherical
coordinates

$\displaystyle X^{0}=\sinh t,$ $\displaystyle X^{1}=\cosh t\cos\phi,$
$\displaystyle X^{2}=\cosh t\sin\phi\;,$ (97)

which cover all of the de Sitter hyperboloid, the metric has the form

$ds^{2}=dt^{2}-\cosh^{2}td\phi^{2}\;.$ (98)

The Christoffel symbols for this metric are

$\Gamma_{\phi t}^{\phi}=\tanh t,\ \ \Gamma_{\phi\phi}^{t}=\sinh t\cosh t\;.$
(99)

The equations for the geodesics are

$\frac{d^{2}\phi}{d\tau^{2}}+2\Gamma_{\phi
t}^{\phi}\frac{d\phi}{d\tau}\frac{dt}{d\tau}=0\;,$ (100)
$\frac{d^{2}t}{d\tau^{2}}+\Gamma_{\phi\phi}^{t}\left(\frac{d\phi}{d\tau}\right)^{2}=0\;.$
(101)

Upon substituting the expressions for the Christoffel symbols (99) into Eqs.
(100)-(101) , we obtain

$\frac{d^{2}\phi}{d\tau^{2}}+2\tanh t\frac{d\phi}{d\tau}\frac{dt}{d\tau}=0\;,$
(102) $\frac{d^{2}t}{d\tau^{2}}+\sinh t\cosh
t\left(\frac{d\phi}{d\tau}\right)^{2}=0\;.$ (103)

Equation (102) gives

$\frac{d\phi}{d\tau}=\frac{A}{\cosh^{2}t}\;,$ (104)

where $A$ is an integration constant. Substituting the expression (104) into
Eq. (103) yields

$\frac{d^{2}t}{d\tau^{2}}=-\frac{A^{2}\sinh t}{\cosh^{3}t}\;.$ (105)

Multiplying Eq. (105) by $2(dt/d\tau)$ produces

$\frac{d}{d\tau}\left(\frac{dt}{d\tau}\right)^{2}=A^{2}\frac{d}{d\tau}\left(\frac{1}{\cosh^{2}t}\right).$
(106)

After integrating Eq. (106), we obtain

$\frac{dt}{d\tau}=\sqrt{\frac{A^{2}}{\cosh^{2}t}+B},$ (107)

where $B$ is an integration constant. Comparing Eqs. (104) and (107) gives the
expression for the geodesics in terms of the variables $\phi$ and $t$:

$\frac{d\phi}{dt}=\frac{A}{\cosh t\sqrt{A^{2}+B\cosh^{2}t}}\;.$ (108)

The equation for a plane passing through the origin of the ambient spacetime
can be written as

$\alpha X_{0}-\beta X_{1}-\gamma X_{2}\;.$ (109)

To obtain the expression for the intersection of (109) with the de Sitter
hyperboloid, one should substitute the expressions (97) into Eq. (109):

$\alpha\sinh t-\beta\cosh t\cos\phi-\gamma\cosh t\sin\phi=0.$ (110)

We can rewrite Eq. (110) as follows:

$\alpha\tanh t=\sqrt{\alpha^{2}+\beta^{2}}\sin(\phi+\phi_{0})\;.$ (111)

Eq. (111) gives the differential equation

$\frac{\alpha}{\cosh^{2}t}dt=\sqrt{\alpha^{2}+\beta^{2}}\cos(\phi+\phi_{0})d\phi\;.$
(112)

From Eq. (111) we find that

$\cos(\phi+\phi_{0})=\sqrt{1-\frac{\alpha^{2}}{\beta^{2}+\gamma^{2}}\tanh^{2}t}\;,$
(113)

and using Eqs. (112) and (113), we obtain

$\frac{d\phi}{dt}=\frac{\alpha}{\cosh
t\sqrt{(\beta^{2}+\gamma^{2}-\alpha^{2})\cosh^{2}t+\alpha^{2}}}\;.$ (114)

The equation (114) coincides with Eq. (108), and hence we have shown that
intersections of planes through the origin with the de Sitter hyperboloid are
geodesics.

This result should not depend on the choice of the coordinate system. However,
for completeness we will reproduce it in the flat coordinates

$\displaystyle X^{0}=\sinh t+\frac{1}{2}x^{2}e^{t},$ $\displaystyle
X^{1}=xe^{t},$ $\displaystyle X^{2}=\cosh t-\frac{1}{2}x^{2}e^{t}\;.$ (115)

In these coordinates the metric is

$ds^{2}=dt^{2}-e^{2t}dx^{2}\;,$ (116)

and the equations for the geodesics are

$\frac{d^{2}x}{d\tau^{2}}+2\frac{dx}{d\tau}\frac{d\phi}{d\tau}=0\;,$ (117)
$\frac{d^{2}t}{d\tau^{2}}+e^{2t}\left(\frac{dx}{d\tau}\right)^{2}\;.$ (118)

Integrating Eq. (117) gives

$\frac{dx}{d\tau}=Ce^{-2t}\;,$ (119)

where $C$ is an integration constant. Substituting the expression (119) into
Eq. (118) and proceeding in a manner similar to that used for the spherical
coordinates yields

$\frac{dt}{d\tau}=\sqrt{C^{2}e^{-2t}+D}\;.$ (120)

Then,

$\frac{dx}{dt}=\frac{Ce^{-2t}}{\sqrt{C^{2}e^{-2t}+D}}\;.$ (121)

This equation can be integrated explicitly:

$x(t)=F-\frac{1}{C}\sqrt{D+C^{2}e^{-2t}}\;.$ (122)

In terms of the coordinates (115), intersections of planes through the origin
of the ambient spacetime with the de Sitter hyperboloid are given by

$\displaystyle\delta\left(\sinh t+\frac{1}{2}x^{2}e^{t}\right)+\varepsilon
xe^{t}+\eta\left(\cosh t-\frac{1}{2}x^{2}e^{t}\right)=0\;.$ (123)

Solving this quadratic equation with respect to $x$, we find

$x(t)=-\frac{\varepsilon}{\delta-\eta}+\sqrt{\frac{\varepsilon^{2}-\delta^{2}+\eta^{2}}{(\delta-\eta)^{2}}+e^{-2t}}.$
(124)

The form of Eq. (124) coincides with Eq. (122).

For any two points in Minkowski spacetime, there always exist a geodesic that
connects these points. This geodesics can be timelike, lightlike or spacelike,
but it always exists. In de Sitter spacetime the situation is different. Let
us take two arbitrary points $X$ and $Y$ on the de Sitter hyperboloid. When is
it possible to connect them by a geodesic?

To answer this question, consider the plane that passes through these two
points and the origin of the coordinates. Note that if some point $U$ belongs
to the intersection of this plane with the de Sitter hyperboloid, then its
antipodal point $\overline{U}=-U$ also belongs to this intersection. If the
angle between this plane and the plane $X^{0}=0$ is less than $\pi/4$, its
intersection with the hyperboloid is an ellipse. In this case, all vectors
tangent to this curve are spacelike, and we can call this geodesic spacelike.
Therefore, $Z(X,Y)<1$ for any two points on the ellipse, and since their
antipodal points are also on this ellipse, $Z(X,{\overline{Y}})<1$ as well.
Hence, any two points that can be connected by an elliptical geodesic should
satisfy the condition that $-1<Z(X,Y)<1$.

If the angle between a plane through a pair of points and the origin and the
plane $X^{0}=0$ is $\geq\pi/4$, its intersection with the hyperboloid gives a
hyperbola with two branches. Any two points on the same branch are timelike
related (or lightlike related if the angle is equal to $\pi/4$), which means
that $Z(X,Y)\geq 1$. The antipodal point of any point on a given branch is
located on the other branch, and hence for any two points lying on different
branches $Z(X,Y)\leq-1$. The latter points, obviously, cannot be not connected
by any geodesic.

##  Appendix B A comparison with the results in $p$-representation

For completeness of presentation and for convenience of readers, in this
Appendix we will present a detailed comparison of our results (90)–(92) with
the two-point function calculated in Gat-Ser  in the framework of momentum
representation Par-Ser . We will use formulas from Gat-Ser  and considerations
of Appendix B of Gautier:2015pca  concerning the relation between Green’s
functions in $p$-representation and in coordinate space.

In $d+1$ dimensional flat de Sitter space, where the Hubble parameter is set
to be equal to 1, the metric is

$ds^{2}=\eta^{-2}(-d\eta^{2}+d{\bf X}\cdot d{\bf X}).$ (125)

A Green’s function can be represented as follows:

$G(x,x^{\prime})=\int\frac{d^{d}K}{(2\pi)^{d}}e^{i{\bf K}\cdot({\bf X}-{\bf
X^{\prime}})}\tilde{G}(\eta,\eta^{\prime},K).$ (126)

De Sitter symmetries ensure that the correlation function has the following
$p$-representation Par-Ser :

$\tilde{G}(\eta,\eta^{\prime},K)=\frac{(\eta\eta^{\prime})^{\frac{d-1}{2}}}{K}\hat{G}(p,p^{\prime}),$
(127)

where

$p=-K\eta,\ \ p^{\prime}=-K\eta^{\prime}.$

Now, let us suppose that for the Green’s function which we would like to
calculate the expression $\hat{G}$ is proportional to

$\frac{\sqrt{pp^{\prime}}}{(pp^{\prime})^{\nu}}.$ (128)

Then, one can show Gautier:2015pca  that the substitution of the expression
(128) into the formula (126) gives the following expression for the case of
the large spacelike interval between the points $x$ and $x^{\prime}$:

$\frac{1}{(4\pi)^{d/2}}\frac{1}{\Gamma(d/2)}\frac{(\eta\eta^{\prime})^{\varepsilon}}{\varepsilon|{\bf
X}-{\bf X^{\prime}}|^{2\varepsilon}},$ (129)

where the parameter $\varepsilon$ is connected with the $d,\nu$ and the scalar
field mass $m$ as follows:

$\nu=\sqrt{\frac{d^{2}}{4}-m^{2}}=\frac{d}{2}-\varepsilon,\ \
\varepsilon\approx\frac{m^{2}}{d}.$ (130)

Using the de Sitter invariance and the formula (130), we can rewrite the
formula (129) as follows:

$\frac{1}{(4\pi)^{d/2}}\frac{d(2|Z|)^{-\varepsilon}}{m^{2}\Gamma(d/2)},$ (131)

where $Z$ is, as usual, the de Sitter invariant function (17).

The symmetrized two-point function calculated in Gat-Ser  is given by (see Eq.
(46) of Gat-Ser )

$\hat{F}(p,p^{\prime})=\sqrt{pp^{\prime}}\tilde{F}_{\nu}\left(\frac{c_{+}}{(pp^{\prime})^{\bar{\nu}_{+}-\varepsilon}}+\frac{c_{-}}{(pp^{\prime})^{\bar{\nu}_{-}-\varepsilon}}\right).$
(132)

Correspondingly, its coordinate representation $F(x,x^{\prime})$ is related to
(132) in the same way that (131) is related to (128):

$F(x,x^{\prime})=\frac{d}{(4\pi)^{d/2}\Gamma(d/2)}\tilde{F}_{\nu}\left(\frac{c_{+}}{m_{+}^{2}}|Z|^{-\varepsilon_{+}}+\frac{c_{m}}{m_{-}^{2}}|Z|^{-\varepsilon_{-}}\right),$
(133)

where

$\tilde{F}_{\nu}=\frac{[2^{\nu}\Gamma(\nu)]^{2}}{4\pi}.$ (134)

Here the exponents $\varepsilon_{+}$ and $\varepsilon_{-}$ are defined as

$\varepsilon_{\pm}=\frac{d}{2}-\bar{\nu}_{\pm}+\varepsilon,$ (135)

where

$\bar{\nu}_{\pm}=\sqrt{\nu^{2}\pm 2\varepsilon\tilde{\nu}+\varepsilon^{2}},$
(136) $\tilde{\nu}=\sqrt{\nu^{2}+\frac{\sigma_{\rho}}{4\varepsilon^{2}}},$
(137)
$\sigma_{\rho}=\frac{\lambda^{2}}{\varepsilon^{2}}\frac{N+2}{24N^{2}}\frac{\Gamma^{2}(d/2)}{4\pi^{d+2}},$
(138)

where $N$ is the number of scalar fields. The quantities $\nu$ and
$\varepsilon$ in Eqs. (133)-(138) are defined as

$\nu=\sqrt{\frac{d^{2}}{4}-M^{2}}=\frac{d}{2}-\varepsilon,\
\varepsilon\approx\frac{M^{2}}{d},$ (139)

where the mass $M$ is obtained by taking into account the local contribution
to the self-energy Gat-Ser . The relation between the mass $M$ and the mass
parameter $m$ of the Lagrangian is given by the equation

$M^{2}\approx m^{2}+\frac{c_{d}\lambda}{M^{2}},$ (140)

where

$c_{d}=\frac{(N+2)d\Gamma\left(\frac{d}{2}\right)}{12N\cdot
2\pi^{\frac{d}{2}+1}}.$ (141)

Solving Eq. (140) perturbatively up to quadratic terms in $\lambda$, one finds

$M^{2}=\frac{m^{2}}{2}+\sqrt{\frac{m^{4}}{4}+c_{d}\lambda}=m^{2}\left(1+\frac{c_{d}\lambda}{m^{4}}-\frac{c_{d}^{2}\lambda^{2}}{m^{8}}\right).$
(142)

Now we are in a position to calculate the exponents $\varepsilon_{\pm}$. First
of all, from Eq. (137) it follows that up to quadratic terms in $\lambda$,

$\tilde{\nu}=\nu+\frac{\sigma_{\rho}}{8\varepsilon^{2}\nu}.$ (143)

Then, from Eq. (136) it follows that

$\bar{\nu}_{\pm}=\nu\pm\varepsilon\pm\frac{\sigma_{\rho}}{8\varepsilon\nu(\nu\pm\varepsilon)},$
(144)

and hence,

$\varepsilon_{+}=\varepsilon-\frac{\sigma_{\rho}}{8\varepsilon\nu^{2}},$ (145)
$\varepsilon_{-}=3\varepsilon+\frac{\sigma_{\rho}}{8\varepsilon\nu^{2}}.$
(146)

Using the formulas written above we obtain the following expression for
$\varepsilon_{+}$:

$\displaystyle\varepsilon_{+}=\frac{m^{2}}{d}\bigg{(}1$ $\displaystyle+$
$\displaystyle\frac{\lambda d(N+2)\Gamma(d/2)}{24N\pi^{d/2+1}m^{4}}$ (147)
$\displaystyle-$
$\displaystyle\frac{\lambda^{2}d^{2}(N+2)(N+5)\Gamma^{2}(d/2)}{576N^{2}\pi^{d+2}m^{8}}\bigg{)}.$

Substituting the values $d=3$ and $N=1$ into Eq. (147) and taking into account
that the coupling constant $\lambda$ in Gat-Ser  is 6 times greater than the
coupling constant in the present paper, we can see that $\varepsilon_{+}$
coincides with the exponent of the first term in our Eqs. (90)–(92).

The expression for $\varepsilon_{-}$ also coincides with the corresponding
exponent in our Eqs. (90)–(92), where its is calculated only at $\lambda^{0}$
order, because the coefficient in front of the exponential is already
proportional to $\lambda^{2}$.

Let us now calculate the coefficients in front of the exponentials in Eq.
(133). The coefficients $c_{\pm}$ are given by Gat-Ser

$c_{\pm}=\frac{\tilde{\nu}\pm\nu}{2\tilde{\nu}}.$ (148)

Using Eqs. (137) and (138), we obtain

$c_{+}=1-\frac{\lambda^{2}d^{2}(N+2)\Gamma^{2}(d/2)}{384N^{2}\pi^{d+2}m^{8}}.$
(149)
$c_{-}=\frac{\lambda^{2}d^{2}(N+2)\Gamma^{2}(d/2)}{384N^{2}\pi^{d+2}m^{8}}.$
(150)

We should also find the “masses” $m_{\pm}$ arising in Eq. (133). They are
given by

$\bar{\nu}_{\pm}-\varepsilon=\sqrt{\frac{d^{2}}{4}-m^{2}_{\pm}}.$ (151)

Using the definitions given above, we find that

$m_{+}^{2}=M^{2}-\frac{\lambda^{2}d^{2}(N+2)\Gamma^{2}(d/2)}{192N^{2}\pi^{d+2}}$
(152)

and

$m_{-}^{2}=3M^{2}+\frac{\lambda^{2}d^{2}(N+2)\Gamma^{2}(d/2)}{192N^{2}\pi^{d+2}}.$
(153)

Recall that in Eq. (133) the expressions $\frac{c_{+}}{m_{+}^{2}}$ and
$\frac{c_{-}}{m_{-}^{2}}$ are multiplied by the factor
$\frac{d}{(4\pi)^{d/2}\Gamma(d/2)}\tilde{F}_{\nu}$. Hence, the coefficient in
front of the second exponential equals

$\frac{d}{(4\pi)^{d/2}\Gamma(d/2)}\tilde{F}_{\nu}\frac{c_{-}}{m_{-}^{2}}=\frac{\lambda^{2}d^{3}(N+2)\Gamma^{3}(d/2)}{4608N^{2}\pi^{\frac{3}{2}d+3}m^{10}}.$
(154)

Rescaling the coupling constant and substituting $N=1$ and $d=3$ into Eq.
(154), we see that it is in agreement with the corresponding coefficient in
our Eqs. (90)–(92).

The coefficient in front of the first exponential in Eq. (133) equals

$\displaystyle\frac{d}{(4\pi)^{d/2}\Gamma(d/2)}\tilde{F}_{\nu}\frac{c_{+}}{m_{+}^{2}}$
$\displaystyle=$
$\displaystyle\frac{d\Gamma(d/2)}{4\pi^{\frac{d}{2}+1}m^{2}}\bigg{(}1-\frac{\lambda
d(N+2)\Gamma(d/2)}{24N\pi^{\frac{d}{2}+1}m^{4}}$ $\displaystyle+$
$\displaystyle\frac{\lambda^{2}d^{2}(N+2)(4N+11)\Gamma^{2}(d/2)}{1152N^{2}\pi^{d+2}m^{8}}\bigg{)}.$

After suitable substitutions, the expression (LABEL:GS28) coincides with the
coefficient of the first term of Eqs. (90)–(92).

## References

  * (1) N. A. Chernikov, E. A. Tagirov, Ann. Inst. H. Poincare Phys. Theor. A 9, 109 (1968)
  * (2) E. A. Tagirov, Ann. Phys. 76, 561 (1973)
  * (3) T. S. Bunch, P. C. W. Davies, Proc. Roy. Soc. Lond. A 360, 117 (1978)
  * (4) A. D. Linde, Phys. Lett. 116B, 335 (1982)
  * (5) A. A. Starobinsky, Phys. Lett. 117B, 175 (1982)
  * (6) A. Vilenkin, L. H. Ford, Phys. Rev. D 26, 1231 (1982)
  * (7) B. Allen, Phys. Rev. D 32, 3136 (1985)
  * (8) B. Allen, A. Folacci, Phys. Rev. D 35, 3771 (1987)
  * (9) A. A. Starobinsky, JETP Lett. 30, 682 (1979)
  * (10) N. C. Tsamis, R. P. Woodard, Nucl. Phys. B 724, 295 (2005)
  * (11) T. Prokopec, N. Tsamis, R. Woodard, Ann. Phys. 323, 1324-1360 (2008)
  * (12) A. Y. Kamenshchik, T. Vardanyan, Phys. Rev. D 102, 065010 (2020)
  * (13) A. A. Starobinsky, Lect. Notes Phys. 246, 107 (1986)
  * (14) A. A. Starobinsky, J. Yokoyama, Phys. Rev. D 50, 6357 (1994)
  * (15) V. Vennin, A. A. Starobinsky, Eur. Phys. J. C 75, 413 (2015)
  * (16) H. Collins, R. Holman, T. Vardanyan, JHEP 1711, 065 (2017)
  * (17) C. P. Burgess, R. Holman, G. Tasinato, M. Williams, JHEP 1503, 090 (2015)
  * (18) C. P. Burgess, R. Holman, G. Tasinato, JHEP 1601, 153 (2016)
  * (19) J. Tokuda and T. Tanaka, JCAP 02, 014 (2018)
  * (20) L. Pinol, S. Renaux-Petel, Y. Tada, JCAP 04, 048 (2021)
  * (21) M. Baumgart and R. Sundrum, JHEP 2020, 119 (2020)
  * (22) V. Gorbenko and L. Senatore, arXiv: 1911.00022 [hep-th]
  * (23) F. Finelli, G. Marozzi, A. A. Starobinsky, G. P. Vacca, G. Venturi, Phys. Rev. D 79, 044007 (2009)
  * (24) F. Finelli, G. Marozzi, A. A. Starobinsky, G. P. Vacca, G. Venturi, Phys. Rev. D 82, 064020 (2010)
  * (25) D. Marolf, I. A. Morrison, Phys. Rev. D 84, 044040 (2011)
  * (26) F. Gautier and J. Serreau, Phys. Lett. B 727, 541 (2013)
  * (27) S. Hollands, Commun. Math. Phys. 319, 1 (2013)
  * (28) E. T. Akhmedov, F. K. Popov, V. M. Slepukhin, Phys. Rev. D 88, 024021 (2013)
  * (29) B. Garbrecht, G. Rigopoulos and Y. Zhu, Phys. Rev. D 89, 063506 (2014)
  * (30) F. Gautier and J. Serreau, Phys. Rev. D 92, 105035 (2015)
  * (31) V.K. Onemli, Phys. Rev. D 91, 103537 (2015)
  * (32) G. Karakaya and V. Onemli, Phys. Rev. D 97, 123531 (2018)
  * (33) T. Markkanen, A. Rajantie, S. Stopyra and T. Tenkanen, JCAP 08, 001 (2019)
  * (34) H. Kitamoto, Phys. Rev. D 100, 025020 (2019)
  * (35) G. Moreau and J. Serreau, Phys. Rev. D 101, 045015 (2020)
  * (36) T. Markkanen and A. Rajantie, JCAP 03, 049 (2020)
  * (37) J. S. Schwinger, J. Math. Phys. 2, 407 (1961)
  * (38) L. V. Keldysh, Zh. Eksp. Teor. Fiz. 47, 1515 (1964) [Sov. Phys. JETP 20, 1018 (1965)]
  * (39) P. M. Bakshi, K. T. Mahanthappa, J. Math. Phys. 4, 1 (1963); P. M. Bakshi, K. T. Mahanthappa, J. Math. Phys. 4,12 (1963)
  * (40) R. D. Jordan, Phys. Rev. D 33, 444 (1986)
  * (41) S. Weinberg, Phys. Rev. D 72, 043514 (2005)
  * (42) H. Collins, arXiv:1101.1308 [astro-ph.CO]
  * (43) A. M. Polyakov, arXiv:1209.4135 [hep-th]
  * (44) M. Beneke and P. Moch, Phys. Rev. D 87, 064018 (2013)
  * (45) A.Y. Kamenshchik, A.A. Starobinsky and T. Vardanyan, work in progress
  * (46) I.S. Gradstein and I.M. Ryzhik, Table of Integrals, Series, and Products, Academic Press, 2015
  * (47) H. Risken, The Fokker-Planck Equation, Springer, New York, 1989
  * (48) R. Parentani and J. Serreau, Phys. Rev. D 87, 045020 (2013)