Spatial correlations of dark energy

from quantum fluctuations in inflation

Enis Belgacem${}^{a}~{}$111e-mail<EMAIL_ADDRESS>Tomislav Prokopeca
222e-mail<EMAIL_ADDRESS>
aInstitute for Theoretical Physics and EMME$\Phi$, Faculty of Science, Utrecht
University,

Princetonplein 5, 3584 CC Utrecht, The Netherlands

This paper contains a detailed study of the properties of a simple model
attempting to explain dark energy as originated from quantum fluctuations of a
light spectator scalar field in inflation. In [1] we recently outlined how
Starobinsky’s stochastic formalism can be used to study the spatial
correlations imprinted on dark energy by its quantum origin in this model and
we studied their possible role in relieving the Hubble tension. Here we
provide a more comprehensive derivation of the results in [1] and we refine
some of our estimates, comparing to the approximate results obtained
previously. Among the main results, we analyze the non-coincident correlators
predicted by a full field theoretical treatment and their relation with those
computed within the stochastic formalism. We find that in the region where
stochastic theory predicts significant sub-Hubble correlators it is in
disagreement with field theoretical predictions. However, agreement can be
restored by introducing a reduced speed of sound for the scalar field. We also
discuss an alternative approach to the problem of studying correlators within
the stochastic formalism based directly on the evolution of probability
distributions. We find that the two approaches give the same answer for
2-point functions of the field, but not for 4-point functions relevant to
density correlators and we discuss the behaviour of the two methods with
respect to Wick’s theorem.

###### Contents

  1. 1 Introduction
  2. 2 Quantum dark energy model
    1. 2.1 Stochastic formalism
    2. 2.2 Non-coincident 2-point IR correlators
      1. 2.2.1 de Sitter inflation
      2. 2.2.2 Radiation epoch
      3. 2.2.3 Matter epoch
    3. 2.3 Non-coincident 4-point IR correlators
      1. 2.3.1 de Sitter inflation
      2. 2.3.2 Radiation epoch
      3. 2.3.3 Matter epoch
    4. 2.4 Matching at current time
  3. 3 Comparing stochastic formalism and quantum field theory
    1. 3.1 Mode functions
    2. 3.2 The field 2-point function
      1. 3.2.1 Coincident correlator
      2. 3.2.2 Super-Hubble regime
      3. 3.2.3 Full comparison
      4. 3.2.4 Effect of a reduced speed of sound
    3. 3.3 Free fields and Wick’s theorem
  4. 4 Stochastic solution from Fokker-Planck equation
    1. 4.1 One-point probability distribution
    2. 4.2 Two-point probability distribution
      1. 4.2.1 Solution for the joint two-point PDF
      2. 4.2.2 Summary on joint PDF and non-coincident correlators
  5. 5 Quantum backreaction
  6. 6 Theoretical aspects of the Hubble tension
    1. 6.1 Distribution of the global squared Hubble rate
    2. 6.2 Joint distribution of the local and global squared Hubble rates
  7. 7 Conclusions and outlook
  8. A Sources of 4-point function stochastic equations
  9. B Correlators in matter + cosmological constant epoch
    1. B.1 Alternative solution
  10. C Effect of non-zero initial conditions in inflation

##  1 Introduction

The discovery of the accelerated expansion of the Universe in 1998 by the
Supernova Cosmology Project [2] and the High-Z Supernova Search Team [3] has
triggered the quest for the origin of cosmic acceleration. The simplest
explanation, namely a cosmological constant, equivalent to a perfect fluid
with equation of state parameter $-1$, has proven to be able to explain a
diversified number of observations. Therefore it has been fully enrolled in
the standard cosmological model, which also requires cold dark matter, to give
the current $\Lambda$CDM model. However, both the theoretical shortcomings
(the cosmological constant problem [4, 5, 6, 7]) and the tensions between the
values of cosmological parameters predicted by probing different ages and
scales of the Universe, most notably in $H_{0}$ (Hubble tension, e.g. between
CMB data [8] and supernovae luminosity distance measurements [9]) and in the
parameters ($\Omega_{M}$, $\sigma_{8}$) affecting galaxy number counts (see
[10]), have stimulated the search for alternatives to $\Lambda$CDM. In such
models, the accelerated expansion at late times is not merely due to a
cosmological constant, but to a dynamical dark energy. A typical investigation
line is given by quintessence models (see [11] for a review), where a field
(typically scalar) evolves in such a way that its energy density becomes
relevant, and eventually lead the expansion, at late times. This desired
behavior can be achieved by the existence of a tracker solution [12, 13, 14],
in which the energy density of the scalar field decreases slightly less
rapidly than radiation during the radiation-dominated epoch and, after matter-
radiation equality, starts decreasing less rapidly than matter density and
eventually behaves like dark energy. In quintessence models, the potential
energy of the scalar field is the essential ingredient for the late-time
acceleration, but it is also possible to obtain dark energy from a scalar
field with non-quadratic kinetic terms, in the so-called ${\it k}$-essence
models [15, 16]. Quintessence and ${\it k}$-essence models are examples of
scalar-tensor theories of gravity with self-accelerated solutions. A larger
set of scalar-tensor theories used in dark energy studies is the Horndeski
class [17, 18, 19], which is made by the most general covariant scalar-tensor
theories with second-order equations of motion. Examples of theories in the
Horndeski class are Brans-Dicke theory [20], galileons [21] and $f(R)$ gravity
[22]. It has also been shown that a wider class of theories (beyond Horndeski
[23, 24], in turn included in the class of degenerate higher-order scalar-
tensor (DHOST) theories [25], see also [26] for a review) can evade
Ostrogradsky instabilities [27], despite leading to equations of motion of
higher order. Quintessence models (and the more general class of theories
mentioned before) attempting to explain dark energy without a cosmological
constant do not properly address the initial conditions issue about the value
the quintessence field needs to start from. Furthemore, they usually assume
perfect correlations of dark energy on large scales, typically without any
justification.

An alternative line of research is based on the possibility that dark energy
could emerge at late times from the amplification of inflationary quantum
fluctuations inherited in the subsequent phases of the Universe evolution [28,
29, 30]. This was shown to be possible with a non-minimally coupled light
scalar field. In such a model, inflation itself, by amplifying infrared
quantum fluctuations, provides the natural initial conditions for the
evolution of the field which will later backreact on the cosmological
expansion, manifesting as dark energy. In [31], based on the result of [32,
29, 30] that the quantum backreaction can be largely ascribed to infrared (IR)
modes, a suitable application of Starobinsky’s stochastic formalism was used
to study the time evolution of the energy density and pressure of a light non-
minimally coupled spectator scalar field. The conclusions agree with the field
theoretic treatment in [32], where the renormalized stress-energy tensor was
evaluated. The approach in [31] required evaluation of coincident 2-point IR
correlators. Later, it was also shown in [1] that the same quantum
fluctuations also imply unperfect spatial correlations of dark energy over a
length scale determined by the comoving Hubble horizon at the beginning of
inflation and that they are significant today for a duration of inflation of
the order of $60$ e-foldings. In [1], energy density correlations were
expressed in terms of non-coincident 4-point functions, whose evolution was
again studied with stochastic formalism. An interesting application of dark
energy spatial correlations is the study of the Hubble tension and the same
work [1] showed a remarkable reduction of the tension from $4.4~{}\sigma$ in
$\Lambda$CDM down to possibly $1~{}\sigma$ in the simple dark energy model
considered. Working out the observational predictions of the model is clearly
of fundamental importance to test in and compare it to $\Lambda$CDM. It was
found in [33] that the model is slightly favored with respect to $\Lambda$CDM,
although not at a statistically significant level. More recently, the work in
[34] has studied a class of phenomenological models of dark energy which
exhibit spatial correlations and their implications on the observed luminosity
distances of supernovae, presenting a discussion on the detectability of such
effects. The result is that the ongoing Dark Energy Survey (DES) [35] is not
able to detect the expected signal, but the upcoming LSST survey [36] can
succeed, thanks to its large sky coverage which reduces the effect of cosmic
variance.

In Section 2 we summarize the techniques and slightly extend the results of
[1] and [31] to compute non-coincident 2-point and 4-point functions within
the stochastic approach. Starobinsky’s stochastic formalism [37, 38] manages
to describe the dynamics of the infrared (long) modes in terms of classical
stochastic equations, where the coupling between long and short modes appears
as a noise source. The stochastic formalism has been used in the long-
wavelength approach of [39]. Its assumptions have been better understood and
the theory was used in the separate universe picture [40, 41] and developed
into the stochastic $\Delta N$ formalism, see e.g. [42, 60, 44]. One of its
prominent applications is in understanding formation of primordial back-holes
[45, 46, 47].

It is remarkable that stochastic theory maps the quantum problem into a
classical stochastic description. However, this is only possible in certain
regimes. In particular, as we will see in Section 3, when considering non-
coincident field correlators, only the behavior for super-Hubble separations
is captured by the stochastic formalism, while shorter scales require a full
field theoretic treatment. Stochastic theory has been tested against quantum
field theory for a scalar field in an exact de Sitter background, without
taking into account the backreaction of the scalar field on the metric, which
breaks the exact de Sitter assumption [48, 49].

One of the main goals of the present paper is to test the stochastic formalism
predictions for non-coincident correlators of a spectator field, not only in
inflation, but also at later times (radiation and matter epochs). In Section
3, we prove that, even though there is agreement at the end of inflation,
stochastic formalism and quantum field theory disagree at near-Hubble and sub-
Hubble scales in matter-dominated epoch, with stochastic theory predicting a
steeper decrease of correlators with distance compared to the full field
theoretic answer. This disagreement can be attributed to the neglect of
spatial gradients by the stochastic formalism, which play an important role in
the post-inflationary evolution of near-Hubble modes. We propose and show that
agreement can be restored by introducing a reduced speed of sound $c_{s}\leq
1$. Such a feature is quite common in Horndeski scalar-tensor theories and it
generally appears in the effective field theory of dark energy (cf. eq. (39)
of [50] or, similarly in the context of inflation, eq. (38) of [51]).
Likewise, we can assume the speed of sound $c_{s}$ as an extra free parameter
in the quantum dark energy model considered here.

Section 3 also shows that 4-point functions predicted by the full field theory
treatment are related to 2-point functions by Wick’s theorem333Wick’s theorem
applies to perturbative quantum field theory and is a consequence of the
Gaussian nature of the free field appearing in the unperturbed Hamiltonian.
This is the case for eqs. (2.1)–(2.2) as long as the quantum backreaction
effects, which make the metric a quantum field, are negligible and
$\hat{\Phi}$ is a free scalar field with a time-dependent mass determined by
the classical background metric. At late times, due to the backreaction of the
field on cosmological expansion, non-Gaussianities are unavoidably generated
because the problem cannot be seen anymore as the evolution of a free field in
a given FLRW background. Indeed, the quantum nature of the field necessarily
produces local fluctuations of the Hubble rate, thus the field and Hubble rate
dynamics have to be solved consistently. Dark energy is contained in the
expectation value of these quantum effects. This idea will be made more
precise in Section 5, where a first approximation is applied to refine the
matching between the model and the current cosmological parameters.. On the
contrary, the 4-point functions evaluated according to the formulation of
stochastic formalism in Section 2 violate Wick’s theorem (more precisely the
formulation there only agrees with Wick’s theorem in the coincident regime and
for large spatial separations, but it fails to reproduce the correct behavior
at intermediate scales). Starting from this disagreement, in Section 4 we
develop an alternative way to apply the stochastic formalism to extract non-
coincident correlations functions, applying the techniques in [56] to the non-
equilibrium regime. This approach is based on the time evolution of classical
probability distributions via Fokker-Planck equations, rather than the coupled
equations for IR correlators in Section 2. We find that the formulation of
stochastic theory in Section 4 agrees with Wick’s theorem and, since this was
also the outcome of the quantum field theory treatment, we expect that the
stochastic theory answer for the 4-point functions in Section 4 is more
trustworthy than the corresponding results in Section 2.

In Section 5, we extend the matching of the model parameters to the usual
cosmological parameters $H_{0}$, $\Omega_{M}$ and $\Omega_{\Lambda}$ presented
in [1] and [31], by taking into account the quantum backreaction effect of the
scalar field on the Hubble rate, which is relevant in the most recent stages
of the Universe evolution.

Then, in Section 6, we discuss how to apply the results from previous sections
to the Hubble tension problem. This part of the paper is meant to extend the
discussion presented in Section III of [1] and to go beyond the approximations
used there. It also provides a possible path for a more consistent numerical
calculation of the Hubble tension probability and the key ideas can be applied
to any model predicting spatial fluctuations of dark energy with known
probability distribution.

In Section 7 we conclude by discussing the main results of the manuscript.

##  2 Quantum dark energy model

The simple dark energy model that we want to consider is based on the
following action for a light non-minimally coupled scalar field $\Phi$ in
$D=4$ spacetime dimensions.

$\boxed{S[\Phi]=\int\\!d^{4}\\!x\,\sqrt{-g}\,\biggl{\\{}-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\Phi\partial_{\nu}\Phi-\frac{1}{2}m^{2}\Phi^{2}-\frac{1}{2}\xi
R\Phi^{2}\biggr{\\}}}\,.$ (2.1)

For a given background metric $g_{\mu\nu}$, the non-minimal coupling term
between the field $\Phi$ and the metric with Ricci scalar $R$ contributes to
the effective mass of the field $M^{2}\equiv m^{2}+\xi R$.

Following [1] and [31], we specialize the action (2.1) to a FLRW background
with metric element $ds^{2}=-dt^{2}+a^{2}(t)d\vec{x}^{2}$, where $t$ is
cosmological time, $\vec{x}$ are the comoving coordinates and $a(t)$ is the
scale factor. Then the field depends on cosmological time and comoving
coordinates, which we denote altogether as $x$, so that
$\Phi(x)=\Phi(t,\vec{x})$. If $\dot{\Phi}$ is the time derivative of $\Phi$
and $\vec{\nabla}\Phi$ its spatial gradient, then the action (2.1) specializes
to

$S[\Phi]\equiv\int\\!dt~{}d^{3}x\,\mathcal{L}_{\Phi}=\int\\!dt~{}d^{3}x\,a^{3}\left\\{\frac{1}{2}\dot{\Phi}^{2}-\frac{\bigl{(}\vec{\nabla}\Phi\bigr{)}^{2}}{2a^{2}}-\frac{1}{2}\left[m^{2}+6\xi(2-\epsilon)H^{2}\right]\Phi^{2}\right\\}\,,$
(2.2)

where $H(t)\equiv\dot{a}(t)/a(t)$ is the (global) Hubble expansion rate and
$\epsilon(t)\equiv-\dot{H}(t)/H^{2}(t)$ is the principal slow-roll parameter.
The (squared) effective mass is

$M^{2}(t)=m^{2}+6\xi\left(2-\epsilon(t)\right)H^{2}(t)\,.$ (2.3)

The assumption of light field means that $m/H(t)<1$ throughout the history of
the Universe. In inflation and radiation epochs we can actually assume
$m/H(t)\ll 1$, but in matter-dominated epoch the mass $m$ plays an important
role because, as we will see, it is necessary for generating a contribution to
the energy-momentum tensor scaling as dark energy and eventually leading the
expansion, which is the ultimate effect we are looking for.

Canonical momentum $\Pi(x)$ is defined as

$\Pi(x)\equiv\frac{\partial\mathcal{L}_{\Phi}}{\partial\dot{\Phi}(x)}=a^{3}\dot{\Phi}(x)\,.$
(2.4)

We can then quantize the scalar field in the classical FLRW background metric
by canonical quantization, promoting $\Phi(x)$ and $\Pi(x)$ to quantum
operators $\hat{\Phi}(x)$ and $\hat{\Pi}(x)$, satisfying equal-time canonical
commutation rules

$\bigl{[}\hat{\Phi}(t,\vec{x}),\hat{\Pi}(t,\vec{x}^{\,\prime})\bigr{]}=i\delta^{3}(\vec{x}\\!-\\!\vec{x}^{\,\prime})\,,\qquad\bigl{[}\hat{\Phi}(t,\vec{x}),\hat{\Phi}(t,\vec{x}^{\,\prime})\bigr{]}=0=\bigl{[}\hat{\Pi}(t,\vec{x}),\hat{\Pi}(t,\vec{x}^{\,\prime})\bigr{]}\,.$
(2.5)

The classical Hamilton equations are mapped into the Heisenberg equations

$\displaystyle\frac{d}{dt}\hat{\Phi}(t,\vec{x})-a^{-3}(t)\,\hat{\Pi}(t,\vec{x})=0\,,$
(2.6) $\displaystyle
a^{-3}(t)\,\frac{d}{dt}\hat{\Pi}(t,\vec{x})-\frac{\nabla^{2}}{a^{2}(t)}\hat{\Phi}(t,\vec{x})+M^{2}(t)\,\hat{\Phi}(t,\vec{x})=0\,,$
(2.7)

which can be combined to get the second order equation of motion444Here we
assume that the backreaction of $\hat{\Phi}$ on the expansion of the Universe
is negligible, and therefore $H(t)$ is the global (classical) expansion rate,
defined as
$H^{2}(t)\equiv\langle\Omega|\hat{H}^{2}\left(t,\vec{x}\right)|\Omega\rangle$,
where the vacuum state $|\Omega\rangle$ obeys eq. (2.16). The global Hubble
rate $H(t)$ does not depend on $\vec{x}$ because of the spatial homogeneity of
the background metric and vacuum state. The local (squared) Hubble rate
operator $\hat{H}^{2}\left(t,\vec{x}\right)$ is discussed in Section 5. for
the quantum field $\hat{\Phi}(t,\vec{x})$,

$\frac{d^{2}}{dt^{2}}\hat{\Phi}(t,\vec{x})+3H(t)\frac{d}{dt}\hat{\Phi}(t,\vec{x})-\frac{\nabla^{2}}{a^{2}(t)}\hat{\Phi}(t,\vec{x})+M^{2}(t)\,\hat{\Phi}(t,\vec{x})=0\,.$
(2.8)

In Section 3 we will be interested in a slight generalization of the equation
of motion, with the introduction of a speed of sound $c_{s}\leq 1$, assumed
constant for simplicity, namely

$\frac{d^{2}}{dt^{2}}\hat{\Phi}(t,\vec{x})+3H(t)\frac{d}{dt}\hat{\Phi}(t,\vec{x})-c_{s}^{2}\frac{\nabla^{2}}{a^{2}}\hat{\Phi}(t,\vec{x})+M^{2}(t)\,\hat{\Phi}(t,\vec{x})=0\,.$
(2.9)

The original model (2.1) has $c_{s}=1$ but, as already discussed in the
Introduction 1, several scalar-tensor theories of gravity predict a reduced
speed of sound, notably galileons and more generally theories in the Horndeski
class.

Unless explicitly specified, we always assume $c_{s}=1$ in this Section. The
introduction of creation/annihilation operators and the construction of the
Fock space is well known and the relevant parts have also been reviewed in [1]
and [31]. Therefore we just limit ourselves to summarize some of the steps,
starting by moving to momentum space and expanding the field and canonical
momentum operators in creation/annihilation operators $\hat{b}(\vec{k})$ and
$\hat{b}^{{\dagger}}(\vec{k})$ for the Fourier mode $\vec{k}$:

$\displaystyle\hat{\Phi}(t,\vec{x})$ $\displaystyle=$
$\displaystyle{\int}\\!\\!\frac{d^{3}k}{(2\pi)^{\frac{3}{2}}}\biggl{\\{}e^{i\vec{k}\cdot\vec{x}}\varphi(t,k)\,\hat{b}(\vec{k})+e^{-i\vec{k}\cdot\vec{x}}\varphi^{*}(t,k)\,\hat{b}^{{\dagger}}(\vec{k})\biggr{\\}}\,,$
(2.10a) $\displaystyle\hat{\Pi}(t,\vec{x})$ $\displaystyle=$ $\displaystyle
a^{3}(t){\int}\\!\\!\frac{d^{3}k}{(2\pi)^{\frac{3}{2}}}\biggl{\\{}e^{i\vec{k}\cdot\vec{x}}\dot{\varphi}(t,k)\,\hat{b}(\vec{k})+e^{-i\vec{k}\cdot\vec{x}}\dot{\varphi}^{*}(t,k)\,\hat{b}^{{\dagger}}(\vec{k})\biggr{\\}}\,.$
(2.10b)

The mode function $\varphi(t,k)$ depends only on $t$ and on the norm of the
comoving momentum $k\\!=\\!\|\vec{k}\|$, due to spatial homogeneity and
isotropy of the FLRW background. The canonical commutation rules (2.5) are
imposed by requiring that

$\bigl{[}\hat{b}(\vec{k}),\hat{b}^{{\dagger}}(\vec{k}^{\,\prime})\bigr{]}=\delta^{3}(\vec{k}\\!-\\!\vec{k}^{\,\prime})\,,\qquad\bigl{[}\hat{b}(\vec{k}),\hat{b}(\vec{k}^{\prime})\bigr{]}=0=\bigl{[}\hat{b}^{{\dagger}}(\vec{k}),\hat{b}^{{\dagger}}(\vec{k}^{\prime})\bigr{]}\,,$
(2.11)

with the Wronskian normalization condition for the mode function

$\varphi(t,k)\dot{\varphi}^{*}(t,k)-\dot{\varphi}(t,k)\varphi^{*}(t,k)=i~{}a^{-3}(t)\,.$
(2.12)

The equation of moton (2.8) implies for the mode function $\varphi(t,k)$,

$\ddot{\varphi}(t,k)+3H(t)\dot{\varphi}(t,k)+\left[\frac{k^{2}}{a^{2}(t)}+M^{2}(t)\right]\,\varphi(t,k)=0\,,$
(2.13)

or, for a generic speed of sound $c_{s}$ in eq. (2.9),

$\ddot{\varphi}(t,k)+3H(t)\dot{\varphi}(t,k)+\left[c_{s}^{2}\frac{k^{2}}{a^{2}(t)}+M^{2}(t)\right]\,\varphi(t,k)=0\,.$
(2.14)

The energy-momentum tensor is evaluated from the variation of the action (2.1)
with respect to the metric. The corresponding quantum operator
$\hat{T}_{\mu\nu}$ is

$\hat{T}_{\mu\nu}=\partial_{\mu}\hat{\Phi}\,\partial_{\nu}\hat{\Phi}-\frac{1}{2}g_{\mu\nu}g^{\alpha\beta}\partial_{\alpha}\hat{\Phi}\,\partial_{\beta}\hat{\Phi}-\frac{m^{2}}{2}g_{\mu\nu}\hat{\Phi}^{2}+\xi\bigl{[}G_{\mu\nu}-\nabla_{\mu}\nabla_{\nu}+g_{\mu\nu}\square\bigr{]}\hat{\Phi}^{2}\,,$
(2.15)

where $G_{\mu\nu}$ is the Einstein tensor, $\nabla_{\mu}$ denotes covariant
derivative and $\square\equiv g^{\alpha\beta}\nabla_{\alpha}\nabla_{\beta}$ is
the covariant d’Alembert operator. On the FLRW background metric, one can
determine the quantum energy density
$\rho_{Q}\equiv\langle\hat{\rho}_{Q}\rangle\equiv-\langle\hat{T}^{0}{}_{0}\rangle$
and pressure $p_{Q}\delta^{i}_{j}\equiv\langle\hat{T}^{i}{}_{j}\rangle$ by
taking expectation values on the homogeneous and isotropic vacuum state
$|\Omega\rangle$ annihilated by all operators $\hat{b}(\vec{k})$, which means

$\hat{b}(\vec{k})|\Omega\rangle=0\,.$ (2.16)

When taking into account the aforementioned homogeneity and isotropy, the
final results are

$\rho_{Q}=\frac{H^{2}}{2}\left\\{\left[\left(\frac{m}{H}\right)^{2}+6\xi\right]\Bigl{\langle}\hat{\Phi}^{2}\Bigr{\rangle}+6\xi\frac{\Bigl{\langle}\left\\{\hat{\Phi},\hat{\Pi}\right\\}\Bigr{\rangle}}{a^{3}H}+\frac{\Bigl{\langle}\hat{\Pi}^{2}\Bigr{\rangle}}{a^{6}H^{2}}+\frac{\Bigl{\langle}\left(\vec{\nabla}\hat{\Phi}\right)^{2}\Bigr{\rangle}}{a^{2}H^{2}}\right\\}\,,$
(2.17) $\displaystyle
p_{Q}=\frac{H^{2}}{2}\left\\{\left[-2\xi(3-2\epsilon)-\left(\frac{m}{H}\right)^{2}(1-4\xi)+24\xi^{2}(2-\epsilon)\right]\Bigl{\langle}\hat{\Phi}^{2}\Bigr{\rangle}+2\xi\frac{\Bigl{\langle}\left\\{\hat{\Phi},\hat{\Pi}\right\\}\Bigr{\rangle}}{a^{3}H}\right.\,$
$\displaystyle\left.+(1-4\xi)\frac{\Bigl{\langle}\hat{\Pi}^{2}\Bigr{\rangle}}{a^{6}H^{2}}-\frac{1-12\xi}{3}\frac{\Bigl{\langle}\left(\vec{\nabla}\hat{\Phi}\right)^{2}\Bigr{\rangle}}{a^{2}H^{2}}\right\\}\,.$
(2.18)

The correlators in eqs. (2.17–2.18) are the coincident
$\Bigl{\langle}\hat{\Phi}^{2}(t,\vec{x})\Bigr{\rangle}$,
$\Bigl{\langle}\hat{\Pi}^{2}(t,\vec{x})\Bigr{\rangle}$,
$\Bigl{\langle}\bigl{\\{}\hat{\Phi}(t,\vec{x}),\hat{\Pi}(t,\vec{x})\bigr{\\}}\Bigr{\rangle}$
(with the usual definition of anticommutator
$\\{\hat{A},\hat{B}\\}\equiv\hat{A}\hat{B}+\hat{B}\hat{A}$) and
$\Bigl{\langle}\left(\vec{\nabla}\hat{\Phi}(t,\vec{x})\right)^{2}\Bigr{\rangle}$.
They only depend on time $t$ and not on the comoving position $\vec{x}$, again
because of the assumed spatial homogeneity and isotropy of the background
metric $g_{\mu\nu}(t)$ and vacuum state $|\Omega\rangle$. When one replaces
the fields/canonical momenta in eqs. (2.17–2.18) by their free field
expansions in eqs. (2.10a–2.10b) (or eqs. (2.20a–2.20b) in the stochastic
approximation that we will discuss), one obtains one-loop results for the
energy density and pressure. In preparation for subsection 2.3, we also write
explicitly the energy density operator (which can be obtained from eq. (2.17)
by removing the expectation values):

$\displaystyle\hat{\rho}_{Q}\left(t,\vec{x}\right)=\frac{H^{2}(t)}{2}\left\\{\left[\left(\frac{m}{H(t)}\right)^{2}+6\xi\right]\hat{\Phi}^{2}\left(t,\vec{x}\right)+6\xi\frac{\left\\{\hat{\Phi}\left(t,\vec{x}\right),\hat{\Pi}\left(t,\vec{x}\right)\right\\}}{a^{3}(t)H(t)}\right.\,$
$\displaystyle\left.+\frac{\hat{\Pi}^{2}\left(t,\vec{x}\right)}{a^{6}(t)H^{2}(t)}+\frac{\left(\vec{\nabla}\hat{\Phi}\left(t,\vec{x}\right)\right)^{2}}{a^{2}(t)H^{2}(t)}\right\\}\,.$
(2.19)

In the next subsection 2.1, when studying the dynamics of the infrared (IR)
modes which are super-Hubble, following [1] and [31] we will neglect the
contribution of spatial gradients
$\left(\vec{\nabla}\hat{\Phi}(t,\vec{x})\right)^{2}$. This is justified
because, for a given mode with comoving wavenumber $k$, the expectation values
of spatial gradient terms entering in eqs. (2.17–2.18) are suppressed by a
factor $k^{2}/(aH)^{2}\ll 1$. We will comment again on the role of spatial
gradients in subsection 3.2.3 when discussing the comparison between
stochastic formalism (briefly revised in subsection 2.1) and the full quantum
field theory (QFT) results. We will see that they are actually relevant after
matter-radiation equality, but their contribution can be made small by a
reduced speed of sound $c_{s}<1$.

###  2.1 Stochastic formalism

As anticipated in the Introduction 1, Starobinsky’s stochastic formalism has
been applied in [31] to determine the time evolution of the coincident
correlators (which are 2-pt functions of field/canonical momentum) appearing
in eqs. (2.17–2.18). The reason lies upon the findings in [32, 29, 30] that
the energy-momentum tensor of a very light non-minimally coupled scalar field,
backreacting on the cosmological expansion, is dominated by the infrared (IR)
modes. The same idea has been applied in [1] to study non-coincident 4-pt
functions entering the density-density correlator. Here we review those
results and include in our treatment the non-coincident 2-pt functions. This
will also allow us to comment on the validity or violation of Wick’s theorem
expected for Gaussian fields and it will partly serve as a basis for the
investigations in the next Sections.

At cosmological time $t$, we set the separation of modes at a scale $\mu
a(t)H(t)$ where $0<\mu<1$ is a dimensionless constant. Then modes with
comoving wavenumber $k<\mu a(t)H(t)$ are said to be long-wavelength modes
(they are super-Hubble) while those with $k>\mu a(t)H(t)$ are short (sub-
Hubble).

The long-wavelength parts of the field and canonical momentum operators,
denoted by $\hat{\phi}(t,\vec{x})$ and $\hat{\pi}(t,\vec{x})$ in contrast to
the full field $\hat{\Phi}(t,\vec{x})$ and canonical momentum
$\hat{\Pi}(t,\vec{x})$ operators in eqs. (2.10a–2.10b), are defined as

$\displaystyle\hat{\phi}(t,\vec{x})$ $\displaystyle=$
$\displaystyle{\int}\\!\\!\frac{d^{3}k}{(2\pi)^{\frac{3}{2}}}\,\theta\left(\mu
a(t)H(t)\\!-\\!\|\vec{k}\|\right)\biggl{\\{}e^{i\vec{k}\cdot\vec{x}}\varphi(t,k)\,\hat{b}(\vec{k})+e^{-i\vec{k}\cdot\vec{x}}\varphi^{*}(t,k)\,\hat{b}^{{\dagger}}(\vec{k})\biggr{\\}}\,,$
(2.20a) $\displaystyle\hat{\pi}(t,\vec{x})$ $\displaystyle=$ $\displaystyle
a^{3}(t){\int}\\!\\!\frac{d^{3}k}{(2\pi)^{\frac{3}{2}}}\,\theta\left(\mu
a(t)H(t)\\!-\\!\|\vec{k}\|\right)\biggl{\\{}e^{i\vec{k}\cdot\vec{x}}\dot{\varphi}(t,k)\,\hat{b}(\vec{k})+e^{-i\vec{k}\cdot\vec{x}}\dot{\varphi}^{*}(t,k)\,\hat{b}^{{\dagger}}(\vec{k})\biggr{\\}}\,.\qquad$
(2.20b)

In eqs. (2.20a–2.20b), the separation between short and long modes has been
set by a top-hat function (Heaviside step function) with transition at $k=\mu
a(t)H(t)$ where $0<\mu\ll 1$ is a dimensionless factor that selects the super-
Hubble UV cutoff of the stochastic theory. The introduction of $\mu\neq 1$
allows for control of the dependence of physical quantities on the UV cutoff
$\mu H(t)$ of the stochastic theory, which limits physical wavenumbers to
$\frac{k}{a(t)}\leq\mu H(t)$.

The dynamics of the long modes is described by a modified version of the
Heisenberg equations of motion,

$\displaystyle\frac{d}{dt}\hat{\phi}(t,\vec{x})-a^{-3}(t)\,\hat{\pi}(t,\vec{x})=\hat{f}_{\phi}(t,\vec{x})\,,$
(2.21) $\displaystyle
a^{-3}(t)\,\frac{d}{dt}\hat{\pi}(t,\vec{x})-\frac{\nabla^{2}}{a^{2}(t)}\hat{\phi}(t,\vec{x})+M^{2}(t)\,\hat{\phi}(t,\vec{x})=a^{-3}(t)\hat{f}_{\pi}(t,\vec{x})\,,$
(2.22)

where the “stochastic forces” on the right-hand sides are due to modes
crossing the separation scale $\mu aH$ and they are given by (see [31]):

$\displaystyle\hat{f}_{\phi}(t,\vec{x})$ $\displaystyle=$ $\displaystyle\mu
aH^{2}\left(1-\epsilon\right){\int}\\!\\!\frac{d^{3}k}{(2\pi)^{\frac{3}{2}}}\,\delta\left(\|\vec{k}\|-\mu
aH\right)\biggl{\\{}e^{i\vec{k}\cdot\vec{x}}\varphi(t,k)\,\hat{b}(\vec{k})+e^{-i\vec{k}\cdot\vec{x}}\varphi^{*}(t,k)\,\hat{b}^{{\dagger}}(\vec{k})\biggr{\\}}\,,$
(2.23) $\displaystyle\hat{f}_{\pi}(t,\vec{x})$ $\displaystyle=$
$\displaystyle\mu
a^{4}H^{2}\left(1-\epsilon\right){\int}\\!\\!\frac{d^{3}k}{(2\pi)^{\frac{3}{2}}}\,\delta\left(\|\vec{k}\|-\mu
aH\right)\biggl{\\{}e^{i\vec{k}\cdot\vec{x}}\dot{\varphi}(t,k)\,\hat{b}(\vec{k})+e^{-i\vec{k}\cdot\vec{x}}\dot{\varphi}^{*}(t,k)\,\hat{b}^{{\dagger}}(\vec{k})\biggr{\\}}\,.\qquad$
(2.24)

###  2.2 Non-coincident 2-point IR correlators

When trading the expectation values in eqs. (2.17–2.18) for their respective
long-wavelength parts and neglecting spatial gradients, the problem of
determining the time evolution of energy density and pressure reduces to the
study of IR coincident correlators (rescaled to have the same dimensions)
$\Delta_{\phi,\phi}(t)$, $\Delta_{\phi,\pi}(t)$ and $\Delta_{\pi,\pi}(t)$
defined in eqs. (30)-(32) of [31].

Similarly to [1], we can generalize them to the following non-coincident
equal-time 2-point correlators, which only depends on time $t$ and relative
comoving separation $r=\|\vec{x}_{2}-\vec{x}_{1}\|$ between the points
considered555The correlator $\Delta_{\phi,\pi}(t,r)$ in eq. (2.26) can be
equivalently defined replacing
$\Bigl{\langle}\bigl{\\{}\hat{\phi}(t,\vec{x}_{1}),\hat{\pi}(t,\vec{x}_{2})\bigr{\\}}\Bigr{\rangle}$
by
$\Bigl{\langle}\hat{\phi}(t,\vec{x}_{1})\hat{\pi}(t,\vec{x}_{2})+\hat{\pi}(t,\vec{x}_{1})\hat{\phi}(t,\vec{x}_{2})\Bigr{\rangle}$,
because it only depends on $\vec{x}_{1}$ and $\vec{x}_{2}$ through their
relative distance $r$. These properties follow from the homogeneity and
isotropy of the metric $g_{\mu\nu}(t)={\rm
diag}\left(-1,a^{2}(t),a^{2}(t),a^{2}(t)\right)$ and the state
$|\Omega\rangle$.:

$\displaystyle\Delta_{\phi,\phi}(t,r)$ $\displaystyle\equiv$
$\displaystyle\Bigl{\langle}\hat{\phi}(t,\vec{x}_{1})\hat{\phi}(t,\vec{x}_{2})\Bigr{\rangle}\,,$
(2.25) $\displaystyle\Delta_{\phi,\pi}(t,r)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{a^{3}(t)H(t)}\Bigl{\langle}\bigl{\\{}\hat{\phi}(t,\vec{x}_{1}),\hat{\pi}(t,\vec{x}_{2})\bigr{\\}}\Bigr{\rangle}\,,$
(2.26) $\displaystyle\Delta_{\pi,\pi}(t,r)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{a^{6}(t)H^{2}(t)}\Bigl{\langle}\hat{\pi}(t,\vec{x}_{1})\hat{\pi}(t,\vec{x}_{2})\Bigr{\rangle}\,.$
(2.27)

The coincident correlators in eqs. (30)-(32) of [31] are obtained by setting
$r=0$ in eqs. (2.25–2.27) and then energy density and pressure are
approximated starting from (2.17–2.18) as

$\rho_{Q}(t)\approx\frac{H^{2}}{2}\left\\{\left[\left(\frac{m}{H}\right)^{2}+6\xi\right]\Delta_{\phi,\phi}(t,0)+6\xi\Delta_{\phi,\pi}(t,0)+\Delta_{\pi,\pi}(t,0)\right\\}\,,$
(2.28) $\displaystyle
p_{Q}(t)\approx\frac{H^{2}}{2}\left\\{\left[-2\xi(3-2\epsilon)-\left(\frac{m}{H}\right)^{2}(1-4\xi)+24\xi^{2}\left(2-\epsilon\right)\right]\Delta_{\phi,\phi}(t,0)\right.\,$
$\displaystyle\left.+2\xi\Delta_{\phi,\pi}(t,0)+(1-4\xi)\Delta_{\pi,\pi}(t,0)\right\\}\,.$
(2.29)

The time evolution of the general non-coincident 2-pt correlators in eqs.
(2.25–2.27) is conveniently described by switching from cosmological time $t$
to the number of e-foldings

$N(t)\equiv\ln{\left(\frac{a(t)}{a_{\rm in}}\right)}\,$ (2.30)

measured from the beginning of inflation, where $a_{\rm in}$ is the scale
factor at the beginning of inflation (corresponding to $N=0$).

Formally the time evolution equations are the same as in eqs. (33)-(35) given
in [31] for coincident correlators,

$\displaystyle\frac{\partial}{\partial
N}\Delta_{\phi,\phi}-\Delta_{\phi,\pi}=n_{\phi,\phi}\,,$ (2.31)
$\displaystyle\frac{\partial}{\partial
N}\Delta_{\phi,\pi}+(3-\epsilon)\Delta_{\phi,\pi}+2\Bigl{(}\frac{M}{H}\Bigr{)}^{\\!2}\Delta_{\phi,\phi}-2\Delta_{\pi,\pi}=n_{\phi,\pi}\,,$
(2.32) $\displaystyle\frac{\partial}{\partial
N}\Delta_{\pi,\pi}+2(3-\epsilon)\Delta_{\pi,\pi}+\Bigl{(}\frac{M}{H}\Bigr{)}^{\\!2}\Delta_{\phi,\pi}=n_{\pi,\pi}\,,$
(2.33)

but now the noise sources $n$’s on the right-hand sides are functions of $t$
and $r$, so eqs. (36)-(38) of [31], connecting the noise sources to the
stochastic forces (2.23–2.24), are generalized by

$\displaystyle n_{\phi,\phi}(t,r)$ $\displaystyle=$
$\displaystyle\frac{1}{H(t)}\Bigl{\langle}\left\\{\hat{f}_{\phi}(t,\vec{x}_{1}),\hat{\phi}(t,\vec{x}_{2})\right\\}\Bigr{\rangle}\,,$
$\displaystyle n_{\phi,\pi}(t,r)$ $\displaystyle=$
$\displaystyle\frac{1}{a^{3}(t)H^{2}(t)}\Bigl{\langle}\left\\{\hat{f}_{\phi}(t,\vec{x}_{1}),\hat{\pi}(t,\vec{x}_{2})\right\\}+\left\\{\hat{f}_{\pi}(t,\vec{x}_{1}),\hat{\phi}(t,\vec{x}_{2})\right\\}\Bigr{\rangle}\,,$
$\displaystyle n_{\pi,\pi}(t,r)$ $\displaystyle=$
$\displaystyle\frac{1}{a^{6}(t)H^{3}(t)}\Bigl{\langle}\left\\{\hat{f}_{\pi}(t,\vec{x}_{1}),\hat{\pi}(t,\vec{x}_{2})\right\\}\Bigr{\rangle}\,.$
(2.34)

Their expressions in terms of the mode function $\varphi(t,k)$, which
generalize eqs. (39)-(42) of [31], are:

$\displaystyle n_{\phi,\phi}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi^{2}}(\mu
aH)^{3}(1-\epsilon)\left[|\varphi(t,k)|^{2}\right]_{k=\mu aH}j_{0}(\mu
aHr)\,,$ $\displaystyle n_{\phi,\pi}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi^{2}}\mu^{3}a^{3}H^{2}(1-\epsilon)\left[\frac{\partial}{\partial
t}|\varphi(t,k)|^{2}\right]_{k=\mu aH}j_{0}(\mu aHr)\,,$ $\displaystyle
n_{\pi,\pi}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi^{2}}\mu^{3}a^{3}H(1-\epsilon)\left[|\dot{\varphi}(t,k)|^{2}\right]_{k=\mu
aH}j_{0}(\mu aHr)\,,$ (2.35)

where $j_{0}(z)\equiv\frac{\sin(z)}{z}$ is the $0$-th order spherical Bessel
function. Comparing with eqs. (39)-(42) of [31], one immediately realizes that
the only modification with respect to the coincident case is the factor
$j_{0}(\mu aHr)$ in the stochastic sources (2.35), which correctly reduces to
$1$ when $r\to 0$.

Similarly to [31], we solve the system of equations (2.31–2.33) in the three
epochs of de Sitter inflation ($\epsilon=0$), radiation-domination
($\epsilon=2$) and matter-domination ($\epsilon=3/2$). The most important
modification is that here we are considering non-coincident 2-pt correlators
instead of the coincident ones of [31], but this is just a simple version of
the technically more complicated problem of non-coincident 4-pt functions
already solved in [1]. Therefore the strategy is already set up and we can
rely on it to construct the solution to eqs. (2.31–2.33). In Appendix B, we
also discuss a refinement for the last stages of evolution by studying
correlators in a Universe containing both matter and a cosmological constant.
This is still an approximation, but it is supposed to (partially) take into
account in analytical form the backreaction of the scalar field on the FLRW
background in the most recent e-foldings of evolution before the current time.
As discussed in Section IV of [31], evolving the correlators in a given
classical FLRW background metric is legitimate as long as the quantum
backreaction on cosmological expansion is small with respect to the classical
sources of energy density guiding the expansion. This is the case throughout
all epochs except for recent times when quantum backreaction manifests as dark
energy and later takes the lead of cosmological expansion. We begin with the
evolution in a de Sitter inflationary epoch.

####  2.2.1 de Sitter inflation

Following [31] and [1], in de Sitter inflation we use the Chernikov-Tagirov-
Bunch-Davies (CTBD) mode function. We also include the possibility of a
reduced speed of sound $c_{s}$: looking at the difference between eqs. (2.13)
and (2.14), it just amounts to replacing $k\to c_{s}k$ in eq. (45) of [31].
Therefore the de Sitter mode function is

$\varphi(t,k)=\sqrt{\frac{\pi}{4a^{3}(t)H_{I}}}~{}H_{\nu_{I}}^{(1)}\left(\frac{c_{s}k}{a(t)H_{I}}\right)\,,\qquad\nu_{I}\equiv\sqrt{\tfrac{9}{4}-\bigl{(}\tfrac{M}{H_{I}}\bigr{)}^{\\!2}}\,,$
(2.36)

where we denoted by $H_{I}$ is the constant Hubble rate in inflation and by
$M^{2}=m^{2}+12\xi H_{I}^{2}$ the constant effective squared mass, while
$H_{\nu_{I}}^{(1)}$ the Hankel function of the first kind. Eq. (2.36)
satisfies the Wronskian condition (2.12).

As in [31], since we are dealing with long modes $k<\mu aH_{I}\ll aH_{I}$, we
can simplify the mode function to

$\varphi(t,k)\approx-\frac{i}{\sqrt{\pi}}2^{\nu_{I}-1}\Gamma(\nu_{I})a(t)^{\nu_{I}-3/2}H_{I}^{\nu_{I}-1/2}c_{s}^{-\nu_{I}}k^{-\nu_{I}}\,.$
(2.37)

For $(m/H_{I})^{2}\ll 1$ and $|\xi|\ll 1$, working at leading order in

$X\equiv M^{2}/H_{I}^{2}=(m/H_{I})^{2}+12\xi\,,$ (2.38)

one gets similarly to [31] (but with the extra factor $j_{0}(\mu aH_{I}r)$)
the expression for the noise sources in de Sitter inflation, which we write in
(row) vector form, to leading order in $X$, as666We used
$\nu_{I}=\sqrt{\frac{9}{4}-X}\simeq\frac{3}{2}-\frac{X}{3}$ and just kept
$\nu_{I}\simeq\frac{3}{2}$ in the power of $c_{s}$.

$\left(n_{\phi,\phi},n_{\phi,\pi},n_{\pi,\pi}\right)\approx
c_{s}^{-3}\frac{H_{I}^{2}}{4\pi^{2}}j_{0}(\mu
aH_{I}r)~{}\left(1,-\frac{2}{3}X,\frac{1}{9}X^{2}\right)\,.$ (2.39)

We discuss in detail the solution of eqs. (2.31–2.33) for the non-coincident
2-pt correlators. First, eqs. (2.31–2.33) can be written in matrix form by
introducing:

1) the column vector $\Delta_{(\rm 2)}$ of correlators whose corresponding row
form (transpose vector $\Delta_{(\rm 2)}^{\rm T}$) is

$\Delta_{(\rm 2)}^{\rm
T}=\left(\Delta_{\phi,\phi},\Delta_{\phi,\pi},\Delta_{\pi,\pi}\right)\,,$
(2.40)

2) the column vector $n_{(\rm 2)}$ of noise sources whose corresponding row
form (transpose vector $n_{(\rm 2)}^{\rm T}$) is

$n_{(\rm 2)}^{\rm T}=\left(n_{\phi,\phi},n_{\phi,\pi},n_{\pi,\pi}\right)\,,$
(2.41)

3) the inflationary (de Sitter) constant matrix evolution for 2-pt functions
(using $\epsilon=0$)

$A_{(\rm 2),I}=\begin{pmatrix}0&-1&0\\\ 2X&3&-2\\\ 0&X&6\\\ \end{pmatrix}\,,$
(2.42)

where the subscript “$(\rm 2),I$” is meant to remind that this matrix refers
to the evolution of 2-pt correlators in inflation. These definitions lead to
the following equivalent form of eqs. (2.31–2.33) in de Sitter,

$\frac{\partial}{\partial N}\Delta_{(\rm 2)}(N,r)+A_{(\rm 2),I}~{}\Delta_{(\rm
2)}(N,r)=n_{(\rm 2)}(N,r)\,.$ (2.43)

The general solution is simple because of the constancy of the matrix $A_{(\rm
2)}$ in de Sitter. It is

$\Delta_{({\rm 2})}(N,r)=\exp{\left[-A_{(\rm 2),I}N\right]}\Delta_{({\rm
2})}(0,r)+\int_{0}^{N}dN^{\prime}~{}\exp{\left[A_{(\rm
2),I}(N^{\prime}-N)\right]}~{}n_{({\rm 2})}(N^{\prime},r)\,.$ (2.44)

As in [31] and [1], we assume zero initial conditions for the correlators
because we are interested in their growth purely generated by the inflationary
expansion, thus $\Delta_{({\rm 2})}(0,r)=0$. We briefly discuss the effect of
non-zero initial conditions inherited from a pre-inflationary epoch in
Appendix C.

The computation of (2.44) is conveniently done by diagonalizing the matrix
$A_{(\rm 2),I}$. Let us call $B$ the change-of-basis matrix (i.e. the matrix
whose columns are the eigenvectors of $A_{(\rm 2),I}$) and by $B^{-1}$ its
inverse. If $\lambda_{j},~{}j=1,2,3$ are the eigenvalues of $A_{(\rm 2)}$,
then the $i$-th 2-pt correlator $\Delta_{(\rm 2),i}(N,r)$ is given by

$\Delta_{({\rm
2}),i}(N,r)=\sum_{j=1}^{3}B_{ij}e^{-\lambda_{j}N}\int_{0}^{N}dN^{\prime}~{}e^{\lambda_{j}N^{\prime}}\sum_{k=1}^{3}(B^{-1})_{jk}~{}n_{({\rm
2}),k}(N^{\prime},r)\,.$ (2.45)

Starting from eq. (2.39), we write

$n_{({\rm 2}),k}(N^{\prime},r)\approx
c_{s}^{-3}\frac{H_{I}^{2}}{4\pi^{2}}~{}j_{0}\left(\mu a_{\rm in}H_{I}{\rm
e}^{N^{\prime}}r\right)~{}\alpha_{k}$ (2.46)

where $\alpha$ is the constant vector with components
$\left(1,-\frac{2}{3}X,\frac{1}{9}X^{2}\right)$ and then eq. (2.45) gives

$\Delta_{({\rm 2}),i}(N,r)\approx
c_{s}^{-3}\frac{H_{I}^{2}}{4\pi^{2}}~{}\sum_{j=1}^{3}\sum_{k=1}^{3}B_{ij}~{}(B^{-1})_{jk}~{}\alpha_{k}~{}F_{j}(N,r)\,,$
(2.47)

where

$F_{j}(N,r)=e^{-\lambda_{j}N}\int_{0}^{N}dN^{\prime}~{}e^{\lambda_{j}N^{\prime}}~{}j_{0}\left(\mu
a_{\rm in}H_{I}r~{}{\rm e}^{N^{\prime}}\right)$ (2.48)

Using the approximation777We have checked numerically that it is a good
approximation. $j_{0}(z)\approx\theta\left(1-z\right)$ (see also [1]) and
introducing the definition

$r_{0}\equiv\left(\mu a_{\rm in}H_{I}\right)^{-1}\,,$ (2.49)

it is straightforward to show that

$F_{j}(N,r)\approx\begin{cases}\frac{1-e^{-\lambda_{j}N}}{\lambda_{j}}&\text{if
$r<r_{0}~{}e^{-N}$}\\\
\frac{e^{-\lambda_{j}N}}{\lambda_{j}}~{}\left[\left(\frac{r}{r_{0}}\right)^{-\lambda_{j}}-1\right]&\text{if
$r_{0}e^{-N}<r<r_{0}$}\\\ 0&\text{if $r>r_{0}$}\,.\end{cases}$ (2.50)

The eigenvalues $\lambda_{j}$ of the matrix $A_{(\rm 2)}$ are
$\left\\{\lambda_{1},\lambda_{2},\lambda_{3}\right\\}=\left\\{3-\sqrt{9-4X},3,3+\sqrt{9-4X}\right\\}$
while the change-of-basis matrix can be chosen (up to any constant factor,
which cancels in eq. (2.47) as

$B=\begin{pmatrix}1&1&1\\\ -\lambda_{1}&-3&-\lambda_{3}\\\
\frac{1}{4}\lambda_{1}^{2}&X&\frac{1}{4}\lambda_{3}^{2}\\\ \end{pmatrix}\,.$
(2.51)

At leading order in $X$, the approximated set of egenvalues is
$\left\\{\frac{2}{3}X,3,6-\frac{2}{3}X\right\\}$. Due to the result for
$F_{j}(N,r)$ in (2.50), one can realize that, after a few e-foldings, the sum
in eq. (2.47) is dominated by the smallest eigenvalue and, with some algebra,
the result for the correlators in eq. (2.47), at leading order in $X$, is

$\begin{pmatrix}\Delta_{\phi,\phi}(N,r)\\\ \Delta_{\phi,\pi}(N,r)\\\
\Delta_{\pi,\pi}(N,r)\\\ \end{pmatrix}\approx
c_{s}^{-3}\frac{H_{I}^{2}}{4\pi^{2}}\begin{pmatrix}1\\\ -\frac{2}{3}X\\\
\frac{1}{9}X^{2}\\\
\end{pmatrix}\times\begin{cases}\frac{3}{2X}\left(1-e^{-\frac{2}{3}XN}\right)&\text{if
$r<r_{0}~{}e^{-N}$}\\\
\frac{3}{2X}~{}e^{-\frac{2}{3}XN}~{}\left[\left(\frac{r}{r_{0}}\right)^{-\frac{2}{3}X}-1\right]&\text{if
$r_{0}e^{-N}<r<r_{0}$}\\\ 0&\text{if $r>r_{0}$}\,.\end{cases}$ (2.52)

It is clear from eq. (2.52) that the amplification of inflationary quantum
fluctuations works at its best for $X<0$ so that the factor
$e^{-\frac{2}{3}XN}$ grows exponentially with $N$. This was already pointed
out in [31] and [1]; the best conditions for a negative $X$ are a negative
non-minimal coupling $\xi<0$ and a very light field so that
$(m/H_{I})^{2}\ll|\xi|\ll 1$, where the last condition ensures consistency of
the approximation $|X|\ll 1$ used so far (because $X\simeq-12|\xi|$). The
requirement $X<0$ for the best enhancement of fluctuations (and therefore a
smaller number of inflationary e-foldings required to match the dark energy
content of the Universe today) can be understood from the effective potential
of the scalar field. We comment again on this in Section 4.1 after eq. (4.6).

With $X\simeq-12|\xi|$ the amplitude of correlators in (2.52) grows
exponentially. At the end of inflation, lasting $N_{I}$ e-foldings, and
assuming $e^{8|\xi|N_{I}}\gg 1$, the 2-pt correlators are

$\boxed{\begin{pmatrix}\Delta_{\phi,\phi}(N_{I},r)\\\
\Delta_{\phi,\pi}(N_{I},r)\\\ \Delta_{\pi,\pi}(N_{I},r)\\\ \end{pmatrix}\simeq
c_{s}^{-3}\frac{H_{I}^{2}}{32\pi^{2}|\xi|}~{}e^{8|\xi|N_{I}}\begin{pmatrix}1\\\
8|\xi|\\\ 16\xi^{2}\\\ \end{pmatrix}s_{(\rm 2)}(r)\,}\,,$ (2.53)

where

$\boxed{s_{(\rm 2)}(r)\simeq\begin{cases}1&\text{if $r<r_{0}~{}e^{-N_{I}}$}\\\
1-\left(\frac{r}{r_{0}}\right)^{8|\xi|}&\text{if
$r_{0}~{}e^{-N_{I}}<r<r_{0}$}\\\ 0&\text{if $r>r_{0}$}\end{cases}}\,.$ (2.54)

The spatial profile at the end of inflation is described by the function
$s_{(\rm 2)}(r)$ (the subscript “$(\rm 2)$” stands for the 2-pt IR
correlators). It depends on the scale $r_{0}=\left(\mu a_{\rm
in}H_{I}\right)^{-1}$ defined in eq. (2.49), which is (up to $\mu^{-1}$) the
comoving Hubble length at the beginning of inflation.

The results (2.53–2.54) generalize to the non-coincident case those in eq.
(58) of [31]. The one-loop energy density $\rho_{Q}(N_{I})$ and pressure
$p_{Q}(N_{I})$ at the end of inflation can be found from the coincident
correlators $\Delta_{\phi,\phi}(N_{I},0)$, $\Delta_{\phi,\pi}(N_{I},0)$ and
$\Delta_{\pi,\pi}(N_{I},0)$ using eqs. (2.28)–(2.29) with $\epsilon=0$ and
$\xi=-|\xi|<0$. The correlators $\Delta_{\phi,\pi}(N_{M},0)$ and
$\Delta_{\pi,\pi}(N_{M},0)$ can been neglected in eqs. (2.69)–(2.70) because
their contribution to energy density and pressure is suppressed by another
factor of $|\xi|$ with respect to $\Delta_{\phi,\phi}(N_{M},0)$. At leading
order in $|\xi|$ and $(m/H)^{2}$,

$\rho_{Q}(N_{I})\approx\left(\frac{m^{2}}{2}-3|\xi|H_{I}^{2}\right)\Delta_{\phi,\phi}(N_{I},0)\,,$
(2.55) $\displaystyle
p_{Q}(N_{M})\approx\left(-\frac{m^{2}}{2}+3|\xi|H_{I}^{2}\right)\Delta_{\phi,\phi}(N_{I},0)\,.$
(2.56)

Eqs. (2.69)–(2.70) show a cosmological constant type contribution
$p_{Q}\approx-\rho_{Q}$ and they agree with eq. (59) of [31]. Of course, the
non-coincident correlators in eq. (2.53) contain more information than simply
$\rho_{Q}$ and $p_{Q}$, because they describe spatial dependence.

####  2.2.2 Radiation epoch

The 2-pt IR correlators at the end of inflation are then inherited in
radiation epoch. In the subsequent evolution, like in [31] and [1], the
contribution of noise sources appearing in the right-hand side of eqs.
(2.31)–(2.33) is negligible with respect to the initial conditions which have
been amplified by inflation with a factor $\sim e^{8|\xi|N_{I}}/|\xi|$. In
post-inflationary epochs the stochastic problem reduces to the classical
evolution of the stochastic initial conditions for the field inherited from
inflation. When neglecting noise sources, the set of eqs. (2.31)–(2.33)
becomes a linear homogeneous system of equations. This implies that the
problem can be reduced to the evolution of coincident correlators, because the
spatial dependence $s_{(\rm 2)}(r)$ in eq. (2.54) inherited from inflation
cannot be further modified by the noise sources. The evolution of coincident
correlators in radiation and matter epochs was solved in [31]. An alternative
way to find them can be obtained by adapting the matrix notation used in eq.
(2.43) to radiation and matter epochs.

In radiation epoch, let us start measuring the number of e-foldings from its
beginning (i.e. from the end of inflation, assuming an instantaneous
reheating). In the limit888Corrections due to the mass $m$ are evaluated in
eqs. (76)-(78) of [31]. They are not necessary to obtain the results in eqs.
(2.62)–(2.63), which agree with [31]. of very light field
$(m/H)^{2}\ll|\xi|\ll 1$, we can study the evolution of correlators neglecting
the mass $m$. Then, in radiation epoch, the 2-pt correlators $\Delta_{(\rm
2)}(N,r)$ evolve from their initial value $\Delta_{(\rm 2)}(N,r)$ as

$\Delta_{({\rm 2})}(N,r)=\exp{\left[-A_{(\rm 2),R}N\right]}\Delta_{({\rm
2})}(0,r)\,,$ (2.57)

where $A_{(\rm 2),R}$ is the constant (in the massless limit) evolution matrix
in radiation epoch ($\epsilon=2$)

$A_{(\rm 2),R}=\begin{pmatrix}0&-1&0\\\ 0&1&-2\\\ 0&0&2\\\ \end{pmatrix}\,.$
(2.58)

From the upper triangular form of $A_{(\rm 2),R}$ one can immediately read the
eigenvalues $\\{0,1,2\\}$, then diagonalize the matrix and compute

$\exp{\left[-A_{(\rm
2),R}N\right]}=\begin{pmatrix}1&1-e^{-N}&\left(1-e^{-N}\right)^{2}\\\
0&e^{-N}&2e^{-N}\left(1-e^{-N}\right)\\\ 0&0&e^{-2N}\\\ \end{pmatrix}\,.$
(2.59)

At the end of radiation epoch (matter-radiation equality), lasting
$N_{R}\approx 50$ e-foldings, this gives

$\exp{\left[-A_{(\rm 2),R}N_{R}\right]}\approx\begin{pmatrix}1&1&1\\\
0&e^{-N_{R}}&2e^{-N_{R}}\\\ 0&0&e^{-2N_{R}}\\\ \end{pmatrix}\,.$ (2.60)

Taking into account the initial conditions inherited from inflation in eq.
(2.53) and their hierarchy
$\Delta_{\phi,\phi}(N_{I},r)\gg\Delta_{\phi,\pi}(N_{I},r)\gg\Delta_{\pi,\pi}(N_{I},r)$
for $|\xi|\ll 1$, the result of eq. (2.57) at the end of radiation epoch
(matter-radiation equality) is

$\boxed{\begin{pmatrix}\Delta_{\phi,\phi}(N_{R},r)\\\
\Delta_{\phi,\pi}(N_{R},r)\\\ \Delta_{\pi,\pi}(N_{R},r)\\\
\end{pmatrix}\simeq\begin{pmatrix}\Delta_{\phi,\phi}(N_{I},r)\\\
\Delta_{\phi,\pi}(N_{I},r)~{}e^{-N_{R}}\\\
\Delta_{\pi,\pi}(N_{I},r)~{}e^{-2N_{R}}\\\
\end{pmatrix}\simeq\begin{pmatrix}\Delta_{\phi,\phi}(N_{I},r)\\\ 0\\\ 0\\\
\end{pmatrix}}\,.$ (2.61)

In the last equality we highlighted that $\Delta_{\phi,\pi}(N_{R},r)$ and
$\Delta_{\pi,\pi}(N_{R},r)$ are completely irrelevant at matter-radiation
equality (and even much earlier). We call $H_{\rm eq}\equiv
H\left(N_{R}\right)$ the Hubble rate at matter-radiation equality.

The corresponding energy density and pressure are obtained from the coincident
correlators by eqs. (2.28)–(2.29) with $\epsilon=2$ and they can be
approximated as

$\rho_{Q}(N_{R})\approx\left(\frac{m^{2}}{2}-3|\xi|H_{\rm
eq}^{2}\right)\Delta_{\phi,\phi}(N_{R},0)\,,$ (2.62) $\displaystyle
p_{Q}(N_{R})\approx\left(-\frac{m^{2}}{2}-|\xi|H_{\rm
eq}^{2}\right)\Delta_{\phi,\phi}(N_{R},0)\,.$ (2.63)

We identify a term behaving like a cosmological constant (CC) and another like
a negative radiation contribution in agreement with eqs. (79), (80) of [31].
Note that the mass $m$ is fundamental to obtain the CC-like contribution with
opposite energy density and pressure. This is not in contradiction with the
fact that we neglected it to compute the result (2.57), because this
assumption was only used to further simplify the time evolution of
correlators, which then enter the energy-momentum tensor via eqs.
(2.28)–(2.29).

####  2.2.3 Matter epoch

The evolution in matter epoch can be determined with the same strategy as for
radiation epoch, namely by neglecting the stochastic sources in the right-hand
sides of eqs. (2.31)–(2.33) and using the left-hand sides of the same
equations to evolve the 2-pt IR correlators from their initial conditions at
matter-radiation equality given in eq. (2.61). Again, we neglect the mass $m$
contribution (see Section IV.C of [31] for the corrections due to it) when
evolving the correlators. However the mass $m$ will still play a fundamental
role for the energy density and pressure, just like we already discussed after
eqs. (2.62)–(2.63).

Let us start measuring the number of e-foldings in matter era from its
beginning, namely from matter-radiation equality. When neglecting the mass $m$
in the time evolution of correlators, the structure of the solution is
analogous to eq. (2.57) obtained in radiation epoch. It takes the form

$\Delta_{({\rm 2})}(N,r)=\exp{\left[-A_{(\rm 2),M}N\right]}\Delta_{({\rm
2})}(0,r)\,,$ (2.64)

where $A_{(\rm 2),M}$ is now the constant matrix (using $\epsilon=3/2$ in
matter era)

$A_{(\rm 2),M}=\begin{pmatrix}0&-1&0\\\ -6|\xi|&3/2&-2\\\ 0&-3|\xi|&3\\\
\end{pmatrix}\,.$ (2.65)

The eigenvalues of the matrix $A_{(\rm 2),M}$ are
$\left\\{\frac{3-\sqrt{9+48|\xi|}}{2},\frac{3}{2},\frac{3-\sqrt{9+48|\xi|}}{2}\right\\}$
or, at leading order in $|\xi|\ll 1$, the set
$\left\\{-4|\xi|,\frac{3}{2},3+4|\xi|\right\\}$. Introducing the diagonal
matrix $\Lambda_{M}\equiv{\rm
diag}\left(\lambda_{1},\lambda_{2},\lambda_{3}\right)\equiv{\rm
diag}\left(-4|\xi|,\frac{3}{2},3+4|\xi|\right)$, we write $A_{(\rm
2),M}=B_{M}\Lambda_{M}B_{M}^{-1}$, where $B_{M}$ is the change-of-basis matrix
given, at leading order in $|\xi|$, by

$B_{M}=\begin{pmatrix}0&-1&0\\\ -6|\xi|&3/2&-2\\\ 0&-3|\xi|&3\\\
\end{pmatrix}\,.$ (2.66)

Then the entries $(ij)$ of the matrix exponential appearing in eq. (2.64) are

$\left(\exp{\left[-A_{(\rm
2),M}N\right]}\right)_{ij}=\sum_{k=1}^{3}\left(B_{M}\right)_{ik}~{}e^{-\lambda_{k}N}~{}\left(B_{M}^{-1}\right)_{kj}\,.$
(2.67)

The sum above is soon dominated by the smallest eigenvalue
$\lambda_{1}\approx-4|\xi|<0$, which gives an exponential growth, while the
other eigenvalues give an exponential decay. Using this observation and the
initial condition (2.61) with negligible $\Delta_{\phi,\pi}$ and
$\Delta_{\pi,\pi}$ components, one easily finds that after $N_{M}$ e-foldings
in matter era the 2-pt correlators given by eq. (2.64) are

$\boxed{\begin{pmatrix}\Delta_{\phi,\phi}(N_{M},r)\\\
\Delta_{\phi,\pi}(N_{M},r)\\\ \Delta_{\pi,\pi}(N_{M},r)\\\
\end{pmatrix}\simeq\begin{pmatrix}1\\\ 4|\xi|\\\ 4\xi^{2}\\\
\end{pmatrix}e^{4|\xi|N_{M}}~{}\Delta_{\phi,\phi}(N_{R},r)}\,.$ (2.68)

The exponential growth as $e^{4|\xi|N_{M}}$ of correlators in matter era is
confirmed by the quantum field theory calculation in eq. (3.37) of this paper,
as well as by [32], but it was missed in eqs. (87)-(89) of [31]. The result
(2.68) serves to reaffirm this growth in matter era also in the context of
stochastic formalism. As we will see in Section 3 within the field theoretic
treatment, the fact that both inflation and matter era give rise to an
exponential growth of correlators for $\xi<0$ and
$\left(m/H\right)^{2}\ll|\xi|\ll 1$ is due to the similarity of their mode
functions in the massless limit: in inflation and in matter era the CTBD mode
functions are Hankel functions of the same order.

The coincident correlators in matter era $\Delta_{\phi,\phi}(N_{M},0)$,
$\Delta_{\phi,\pi}(N_{M},0)$ and $\Delta_{\pi,\pi}(N_{M},0)$ determine the
energy density and pressure via eqs. (2.28)–(2.29), using $\epsilon=3/2$. At
leading order in $|\xi|$ and $\left(m/H(N_{M})\right)^{2}$, they are
approximated999Similarly to inflation, the correlators
$\Delta_{\phi,\pi}(N_{M},0)$ and $\Delta_{\pi,\pi}(N_{M},0)$ have been
neglected in eqs. (2.69)–(2.70) because their contribution to energy density
and pressure is suppressed by another factor of $|\xi|$ with respect to
$\Delta_{\phi,\phi}(N_{M},0)$. by

$\rho_{Q}(N_{M})\approx\left(\frac{m^{2}}{2}-3|\xi|H^{2}(N_{M})\right)\Delta_{\phi,\phi}(N_{M},0)\,,$
(2.69) $\displaystyle
p_{Q}(N_{M})\approx-\frac{m^{2}}{2}\Delta_{\phi,\phi}(N_{M},0)\,.$ (2.70)

By inspection of eqs. (2.69)–(2.70), we can see a CC-like (cosmological
constant) term and a negative matter-like contribution, in agreement with eq.
(90) of [31]. This result can be used to match the model to the current
cosmological parameters and reproduce the right amount of dark energy. The
consequences of this step are revised in the subsection 2.4. Before doing
that, we discuss the other important feature of the dark energy model of
quantum origin treated in this paper, namely spatial correlations. For this
purpose, a study of non-coincident 4-pt functions is needed.

###  2.3 Non-coincident 4-point IR correlators

Since in the model described in this paper dark energy is a consequence of
quantum fluctuations amplified by inflation, one can expect that important
physical information can be extracted not only from the expectation value of
the energy-momentum tensor at a single point (as in eqs. (2.28)–(2.29)), but
also from correlations between the energy-momentum tensor at two different
space points. Thus, a source of information on the features of dark energy
predicted by the simple action (2.1) is in the correlator
$\langle\hat{\rho}_{Q}(t,\vec{x}_{1})\hat{\rho}_{Q}(t,\vec{x}_{2})\rangle$
between energy densities101010The correlator
$\langle\hat{\rho}_{Q}(x)\hat{\rho}_{Q}(x^{\prime})\rangle$ is an essential
ingredient in building the
$\langle\hat{T}_{\mu\nu}(x)\hat{T}_{\rho\sigma}(x^{\prime})\rangle$ correlator
which, together with the local contribution from the graviton 4-pt vertex,
builds the one-loop graviton self-energy. Therefore, studying the
$\langle\hat{\rho}_{Q}(x)\hat{\rho}_{Q}(x^{\prime})\rangle$ correlator teaches
us something on the off-coincident one-loop graviton self-energy. at the two
comoving positions $\vec{x}_{1}$ and $\vec{x}_{2}$, evaluated at the same time
$t$. This was one of the main subjects of study in [1]. Here we want to
complete the presentation given in that paper and show all the
definitions/equations needed, which were omitted in [1] for brevity. In
stochastic formalism one approximates the occurrences of field/canonical
momentum in correlators with their long-wavelength parts. Starting from eq.
(2.19) and neglecting spatial gradients, the density-density correlator (at
equal times) is expressed in terms of six 4-point IR correlators as follows:

$\displaystyle\langle\hat{\rho}_{Q}(t,\vec{x}_{1})\hat{\rho}_{Q}(t,\vec{x}_{2})\rangle\approx\frac{H^{4}}{4}\biggl{\\{}\Bigl{[}\Bigl{(}\frac{m}{H}\Bigr{)}^{\\!2}+6\xi\Bigr{]}^{2}\Delta_{\phi^{2},\phi^{2}}(t,r)+6\xi\Bigl{[}\Bigl{(}\frac{m}{H}\Bigr{)}^{\\!2}+6\xi\Bigr{]}\Delta_{\phi^{2},\phi\pi}(t,r)\,$
$\displaystyle+36\xi^{2}\Delta_{\phi\pi,\phi\pi}(t,r)+\Bigl{[}\Bigl{(}\frac{m}{H}\Bigr{)}^{\\!2}+6\xi\Bigr{]}\Delta_{\phi^{2},\pi^{2}}(t,r)+6\xi\Delta_{\phi\pi,\pi^{2}}(t,r)+\Delta_{\pi^{2},\pi^{2}}(t,r)\biggr{\\}}\,.$
(2.71)

The complete definitions of the 4-point functions needed111111The definition
of the first correlator $\Delta_{\phi^{2},\phi^{2}}(t,r)$ was already given in
eq. (9) of [1]. are:

$\displaystyle\Delta_{\phi^{2},\phi^{2}}(t,r)$ $\displaystyle\equiv$
$\displaystyle\Bigl{\langle}\hat{\phi}^{2}(t,\vec{x}_{1})\hat{\phi}^{2}(t,\vec{x}_{2})\Bigr{\rangle}\,,$
(2.72) $\displaystyle\Delta_{\phi^{2},\phi\pi}(t,r)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{a^{3}(t)H(t)}\Bigl{\langle}\hat{\phi}^{2}(t,\vec{x}_{1})\left\\{\hat{\phi}(t,\vec{x}_{2}),\hat{\pi}(t,\vec{x}_{2})\right\\}+\left\\{\hat{\phi}(t,\vec{x}_{1}),\hat{\pi}(t,\vec{x}_{1})\right\\}\hat{\phi}^{2}(t,\vec{x}_{2})\Bigr{\rangle}\,,$
(2.73) $\displaystyle\Delta_{\phi\pi,\phi\pi}(t,r)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{a^{6}(t)H^{2}(t)}\Bigl{\langle}\left\\{\hat{\phi}(t,\vec{x}_{1}),\hat{\pi}(t,\vec{x}_{1})\right\\}\left\\{\hat{\phi}(t,\vec{x}_{2}),\hat{\pi}(t,\vec{x}_{2})\right\\}\Bigr{\rangle}\,,$
(2.74) $\displaystyle\Delta_{\phi^{2},\pi^{2}}(t,r)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{a^{6}(t)H^{2}(t)}\Bigl{\langle}\hat{\phi}^{2}(t,\vec{x}_{1})\hat{\pi}^{2}(t,\vec{x}_{2})+\hat{\pi}^{2}(t,\vec{x}_{1})\hat{\phi}^{2}(t,\vec{x}_{2})\Bigr{\rangle}\,,$
(2.75) $\displaystyle\Delta_{\phi\pi,\pi^{2}}(t,r)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{a^{9}(t)H^{3}(t)}\Bigl{\langle}\hat{\pi}^{2}(t,\vec{x}_{1})\left\\{\hat{\phi}(t,\vec{x}_{2}),\hat{\pi}(t,\vec{x}_{2})\right\\}+\left\\{\hat{\phi}(t,\vec{x}_{1}),\hat{\pi}(t,\vec{x}_{1})\right\\}\hat{\pi}^{2}(t,\vec{x}_{2})\Bigr{\rangle}\,,$
(2.76) $\displaystyle\Delta_{\pi^{2},\pi^{2}}(t,r)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{a^{12}(t)H^{4}(t)}\Bigl{\langle}\hat{\pi}^{2}(t,\vec{x}_{1})\hat{\pi}^{2}(t,\vec{x}_{2})\Bigr{\rangle}\,.$
(2.77)

where, as always, $r\equiv\|\vec{x}_{2}-\vec{x}_{1}\|$ is the relative
comoving distance between the two points. Similarly to the 2-pt IR correlator,
also these non-coincident 4-pt IR correlators evolve in time under the effect
of noise sources. In terms of the number of e-foldings (2.30), the full system
of equations describing the process (including eqs. (10), (11) in [1]) is

$\displaystyle\frac{\partial}{\partial
N}\Delta_{\phi^{2},\phi^{2}}-\Delta_{\phi^{2},\phi\pi}=n_{\phi^{2},\phi^{2}}\,,$
(2.78) $\displaystyle\frac{\partial}{\partial
N}\Delta_{\phi^{2},\phi\pi}+(3-\epsilon)\Delta_{\phi^{2},\phi\pi}-2\Delta_{\phi\pi,\phi\pi}-2\Delta_{\phi^{2},\pi^{2}}+4\Bigl{(}\frac{M}{H}\Bigr{)}^{\\!2}\Delta_{\phi^{2},\phi^{2}}=n_{\phi^{2},\phi\pi}\,,$
(2.79) $\displaystyle\frac{\partial}{\partial
N}\Delta_{\phi\pi,\phi\pi}+2(3-\epsilon)\Delta_{\phi\pi,\phi\pi}-2\Delta_{\phi\pi,\pi^{2}}+2\Bigl{(}\frac{M}{H}\Bigr{)}^{\\!2}\Delta_{\phi^{2},\phi\pi}=n_{\phi\pi,\phi\pi}\,,$
(2.80) $\displaystyle\frac{\partial}{\partial
N}\Delta_{\phi^{2},\pi^{2}}+2(3-\epsilon)\Delta_{\phi^{2},\pi^{2}}-\Delta_{\phi\pi,\pi^{2}}+\Bigl{(}\frac{M}{H}\Bigr{)}^{\\!2}\Delta_{\phi^{2},\phi\pi}=n_{\phi^{2},\pi^{2}}\,,$
(2.81) $\displaystyle\frac{\partial}{\partial
N}\Delta_{\phi\pi,\pi^{2}}+3(3-\epsilon)\Delta_{\phi\pi,\pi^{2}}-4\Delta_{\pi^{2},\pi^{2}}+2\Bigl{(}\frac{M}{H}\Bigr{)}^{\\!2}\Delta_{\phi^{2},\pi^{2}}+2\Bigl{(}\frac{M}{H}\Bigr{)}^{\\!2}\Delta_{\phi\pi,\phi\pi}=n_{\phi\pi,\pi^{2}}\,,$
(2.82) $\displaystyle\frac{\partial}{\partial
N}\Delta_{\pi^{2},\pi^{2}}+4(3-\epsilon)\Delta_{\pi^{2},\pi^{2}}+\Bigl{(}\frac{M}{H}\Bigr{)}^{\\!2}\Delta_{\phi\pi,\pi^{2}}=n_{\pi^{2},\pi^{2}}\,.$
(2.83)

The sources appearing at the right-hand sides of the equations above can be
expressed in terms of the field operator $\hat{\phi}$, canonical momentum
operator $\hat{\pi}$ and stochastic forces $\hat{f}_{\phi}$, $\hat{f}_{\pi}$.
They are given by eq. (A.1) in Appendix A. Starting from those definitions,
they can be written in momentum space in terms of integrals involving the mode
function $\varphi(t,k)$. The resulting expressions are listed in eq. (A.2).

One can then study how the 4-pt correlators evolve during different
cosmological epochs (de Sitter, radiation domination and matter domination)
via their coupled equations (2.78)–(2.83), following essentially the same
steps that we already explained for the evolution of 2-pt correlators in eqs.
(2.31)–(2.33). Since this was already done in [1]) (see in particular its
Appendix A), we recall the results from there and make additional comments.

####  2.3.1 de Sitter inflation

Similarly to the definition (2.40), let us collect the 4-pt functions
(2.72)–(2.77) in a column vector $\Delta_{(\rm 4)}$, whose corresponding row
form (transpose vector $\Delta_{(\rm 4)}^{\rm T}$) is

$\Delta_{(\rm 4)}^{\rm
T}\equiv\left(\Delta_{\phi^{2},\phi^{2}},\Delta_{\phi^{2},\phi\pi},\Delta_{\phi\pi,\phi\pi},\Delta_{\phi^{2},\pi^{2}},\Delta_{\phi\pi,\pi^{2}},\Delta_{\pi^{2},\pi^{2}}\right)\,,$
(2.84)

and we can do the same thing for the noise sources defining a column vector
$n_{(\rm 4)}$, such that

$n_{(\rm 4)}^{\rm
T}\equiv\left(n_{\phi^{2},\phi^{2}},n_{\phi^{2},\phi\pi},n_{\phi\pi,\phi\pi},n_{\phi^{2},\pi^{2}},n_{\phi\pi,\pi^{2}},n_{\pi^{2},\pi^{2}}\right)\,.$
(2.85)

The system (2.78)–(2.83) in de Sitter inflation takes the compact form

$\frac{\partial}{\partial N}\Delta_{(\rm 4)}(N,r)+A_{(\rm 4),I}~{}\Delta_{(\rm
4)}(N,r)=n_{(\rm 4)}(N,r)\,,$ (2.86)

where the constant matrix $A_{(\rm 4),I}$ and the noise functions $n_{(\rm
4)}(N,r)$ were evaluated in eqs. (A4) and (A5) of [1]). The matrix $A_{(\rm
4),I}$, in terms of $X\equiv(M/H_{I})^{2}=(m/H_{I})^{2}+12\xi$, is

$A_{(\rm 4),I}=\begin{pmatrix}0&-1&0&0&0&0\\\ 4X&3&-2&-2&0&0\\\
0&2X&6&0&-2&0\\\ 0&X&0&6&-1&0\\\ 0&0&2X&2X&9&-4\\\ 0&0&0&0&X&12\\\
\end{pmatrix}\,.$ (2.87)

We also rewrite here the expression for the noise sources (as a row vector)
including the effect of a general speed of sound, which was not studied in
[1]:

$n_{(\rm 4)}^{\rm
T}(N,r)=c_{s}^{-6}\left(\frac{H_{I}}{2\pi}\right)^{4}\left[1+2j_{0}(\mu a_{\rm
in}{\mathrm{e}}^{N}H_{I}r)\right]\left(\frac{1}{4|\xi|},4,16|\xi|,8|\xi|,64\xi^{2},64|\xi|^{3}\right)\,.$
(2.88)

Using a method perfectly analogous to that in subsection 2.2.1, one can prove
(see eq. (A7) of [1]) that, assuming zero initial conditions at the beginning
of inflation and working with $\xi<0$ and $(m/H_{I})^{2}\ll|\xi|\ll 1$, then
at the end of inflation the non-coincident 4-pt correlators evolve into

$\boxed{\Delta_{(\rm 4)}^{\rm T}(N_{I},r)\simeq
c_{s}^{-6}\frac{H_{I}^{4}}{1024\pi^{4}\xi^{2}}e^{16|\xi|N_{I}}s_{(\rm
4)}(r)~{}(1,16|\xi|,64\xi^{2},32\xi^{2},256|\xi|^{3},256\xi^{4})}\,,$ (2.89)

where the function $s_{(\rm 4)}(r)$ containing the spatial dependence is:

$\boxed{s_{(\rm 4)}(r)\simeq\begin{cases}3&\text{if $r<r_{0}~{}e^{-N_{I}}$}\\\
3-2\left(\frac{r}{r_{0}}\right)^{16|\xi|}&\text{if
$r_{0}~{}e^{-N_{I}}<r<r_{0}$}\\\ 1&\text{if $r>r_{0}$}\,.\end{cases}}$ (2.90)

The scale $r_{0}\equiv\left(\mu a_{\rm in}H_{I}\right)^{-1}$ was already
introduced in eq. (2.49). It is interesting to compare the spatial dependence
of the 4-pt functions with the one of 2-pt functions. For definiteness, let us
focus on $\Delta_{\phi^{2},\phi^{2}}(N_{I},r)$ and
$\Delta_{\phi,\phi}(N_{I},r)$ given by the first component of eqs. (2.89) and
(2.53), respectively.

One would expect from Wick’s theorem (which holds at one-loop level) that

$\Delta_{\phi^{2},\phi^{2}}(N_{I},r)\stackrel{{\scriptstyle?}}{{=}}\left[\Delta_{\phi,\phi}(N_{I},0)\right]^{2}+2\left[\Delta_{\phi,\phi}(N_{I},r)\right]^{2}\,.$
(2.91)

Comparing the results

$\Delta_{\phi^{2},\phi^{2}}(N_{I},r)\simeq
c_{s}^{-6}\frac{H_{I}^{4}}{1024\pi^{4}\xi^{2}}e^{16|\xi|N_{I}}s_{(\rm
4)}(r)\,$ (2.92)

and

$\Delta_{\phi,\phi}(N_{I},r)\simeq
c_{s}^{-3}\frac{H_{I}^{2}}{32\pi^{2}|\xi|}e^{8|\xi|N_{I}}s_{(\rm 2)}(r)\,,$
(2.93)

eq. (2.91) is equivalent to check whether

$s_{(\rm 4)}(r)\stackrel{{\scriptstyle?}}{{=}}\left[s_{(\rm
2)}(0)\right]^{2}+2\left[s_{(\rm 2)}(r)\right]^{2}\,.$ (2.94)

But using eq. (2.54), the right-hand side of eq. (2.94) is

$\left[s_{(\rm 2)}(0)\right]^{2}+2\left[s_{(\rm
2)}(r)\right]^{2}\simeq\begin{cases}3&\text{if $r<r_{0}~{}e^{-N_{I}}$}\\\
3-4\left(\frac{r}{r_{0}}\right)^{8|\xi|}+2\left(\frac{r}{r_{0}}\right)^{16|\xi|}&\text{if
$r_{0}~{}e^{-N_{I}}<r<r_{0}$}\\\ 1&\text{if $r>r_{0}$}\,.\end{cases}$ (2.95)

This is clearly different from the function $s_{(\rm 4)}(r)$ in eq. (2.90).
More precisely they are equal only for $r<r_{0}~{}e^{-N_{I}}$ (basically the
coincident regime) and $r>r_{0}$ (large distance regime), but they disagree at
intermediate scales $r_{0}~{}e^{-N_{I}}<r<r_{0}$ when interpolating between
the extreme values of $3$ and $1$. We infer that Wick’s theorem is not
satisfied by the result for 4-pt correlators (2.89)–(2.90). This is a problem
because the spectator scalar field $\hat{\Phi}$, described by the action
(2.1), is free121212This is true when quantum gravitational effects are
neglected and when the quantum backreaction of the scalar is negligibly small.
when considering a given classical background metric, which is legitimate
throughout all the evolution of the Universe, except for the very few most
recent e-foldings of evolution in matter era (see also the discussion after
eq. (3.46)). We do not fully understand the origin of this problem, but we
think it lies in the form of the 4-pt stochastic sources $n_{(\rm 4)}(N,r)$.
It is possible that they do not fully catch the interaction between short and
long modes. Section 3.3 shows that, as expected for free fields, quantum field
theory (QFT) predicts that Wick’s theorem is obeyed. QFT says that the spatial
dependence $s_{(\rm 2)}(r)$ of 2-pt functions in eq. (2.54) is correct, while
$s_{(\rm 4)}(r)$ for 4-pt functions in eq. (2.90) is not. An alternative way
to apply (and rescue) stochastic formalism is proposed in Section 4 and it
works with classical probability distributions instead of systems of equations
for IR correlators. One of the results therein is that the joint probability
distribution of fields at two points is Gaussian and therefore correlators
evaluated from it obey Wick’s theorem, see Fig. 5.

####  2.3.2 Radiation epoch

Once inflation ends, the Unverse enters an epoch where expansion is dominated
by radiation. In this case $\epsilon=2$ and the quantity $(M/H)^{2}$ appearing
in eqs. (2.78)–(2.83) is small since it reduces to the mass contribution
$(m/H)^{2}$. For our purposes it can be safely neglected, like we did in
subsection 2.2.2. Furthermore, similarly to the case of 2-pt functions, the
stochastic sources $n_{(\rm 4)}$ are much less relevant than the initial
conditions inherited from inflation, because only initial conditions contain
an enhancement factor that is exponential in $N_{I}$. Then, with the same
technique used for 2-pt functions in subsection 2.2.2, the initial conditions
(2.89) are easily evolved through radiation domination and, at the end of it
(lasting $N_{R}\approx 50$ e-foldings until matter-radiation equality) the
4-point functions are (see also Appendix A.2 of [1])

$\boxed{\begin{pmatrix}\Delta_{\phi^{2},\phi^{2}}(N_{R},r)\\\
\Delta_{\phi^{2},\phi\pi}(N_{R},r)\\\ \Delta_{\phi\pi,\phi\pi}(N_{R},r)\\\
\Delta_{\phi^{2},\pi^{2}}(N_{R},r)\\\ \Delta_{\phi\pi,\pi^{2}}(N_{R},r)\\\
\Delta_{\pi^{2},\pi^{2}}(N_{R},r)\\\
\end{pmatrix}\simeq\begin{pmatrix}\Delta_{\phi^{2},\phi^{2}}(N_{I},r)\\\
\Delta_{\phi^{2},\phi\pi}(N_{I},r)~{}e^{-N_{R}}\\\
\Delta_{\phi\pi,\phi\pi}(N_{I},r)~{}e^{-2N_{R}}\\\
\Delta_{\phi^{2},\pi^{2}}(N_{I},r)~{}e^{-2N_{R}}\\\
\Delta_{\phi\pi,\pi^{2}}(N_{I},r)~{}e^{-3N_{R}}\\\
\Delta_{\pi^{2},\pi^{2}}(N_{I},r)~{}e^{-4N_{R}}\\\
\end{pmatrix}\simeq\begin{pmatrix}\Delta_{\phi^{2},\phi^{2}}(N_{I},r)\\\ 0\\\
0\\\ 0\\\ 0\\\ 0\\\ \end{pmatrix}}\,.$ (2.96)

Note that, within the approximations mentioned, $\Delta_{\phi^{2},\phi^{2}}$
is constant during the radiation period, while the other correlators are
suppressed by a factor of $e^{-N_{R}}$ for each occurrence of a canonical
momentum in their definition. This is perfectly analogous to the result for
2-point functions in eq. (2.61). At the end of radiation domination
($N_{R}\approx 50$), all the correlators involving canonical momentum are
negligible, which is the content of the last equality in eq. (2.96). Only
$\Delta_{\phi^{2},\phi^{2}}$ survives (roughly unchanged) until matter-
radiation equality. In matter-dominated epoch it will “leak” again into the
other correlators.

####  2.3.3 Matter epoch

At matter-radiation equality, non-relativistic matter takes the lead of
cosmological expansion. Eqs. (2.69), (2.70) show that the scalar field
contributes to the energy-momentum tensor with a cosmological constant (CC)
type portion and a negative matter-like portion. The CC-like part is what
manifests as dark energy at recent cosmological times. In Section II.C and
Appendix A.3 of [1], the evolution of the scalar field was studied in a
Universe containing both matter and a cosmological constant, as a way to take
into account the backreaction of the scalar field on the expansion. Appendix B
can be used to review (and even refine) those results. Here, for a simple and
more fair comparison with the discussion on 2-pt functions in subsection
2.2.3, we give the results for 4-pt functions in pure matter epoch. Using the
same approximations of subsection 2.2.3, one arrives with similar calculations
to the following 4-pt IR correlators in matter era, after $N_{M}$ e-foldings
from matter-radiation equality:

$\boxed{\Delta_{(\rm 4)}^{\rm
T}(N_{M},r)\simeq\Delta_{\phi^{2},\phi^{2}}(N_{R},r)~{}e^{8|\xi|N_{M}}~{}\left(1,8|\xi|,16\xi^{2},8\xi^{2},32|\xi|^{3},16\xi^{4}\right)}\,.$
(2.97)

At leading order in $|\xi|$ and $\left(m/H(N_{M})\right)^{2}$, the
corresponding density-density correlator computed from eq. (2.71) is dominated
by $\Delta_{\phi^{2},\phi^{2}}(N_{M},r)$ and gives

$\boxed{\langle\hat{\rho}_{Q}(N_{M},\vec{x}_{1})\hat{\rho}_{Q}(N_{M},\vec{x}_{2})\rangle\approx\left[\rho_{Q}(N_{M})\right]^{2}~{}s_{(\rm
4)}(r)\qquad\text{with~{}~{}$r\equiv\|\vec{x}_{2}-\vec{x}_{1}\|$}}\,,$ (2.98)

where we also made use of eqs. (2.68), (2.69) and the function $s_{(\rm
4)}(r)$ was defined in (2.90). Once again, the scale of spatial correlations
is determined by $r_{0}$ appearing in $s_{(\rm 4)}(r)$ and therefore by the
comoving Hubble length at the beginning of inflation (from eq. (2.49)).

###  2.4 Matching at current time

The model parameters have to be matched to the current cosmological parameters
in order to study its actual predictions. In the last few e-foldings the
backreaction of the scalar field on the background metric due to its energy-
momentum tensor becomes very important, to the point that eventually it
becomes a fundamental contribution to the total energy density of the Universe
and in the future it will lead the expansion. A full solution of the problem
would unavoidably be numerical. However some simplified understanding is
possible in analytical form. For the moment we approximate the FLRW background
simply as a matter-dominated Universe, but in Section 5 we will refine the
matching by a more consistent account for the quantum backreaction at late
times. As usual in the literature, let us call $\Omega_{M}$ the energy density
fraction today due to non-relativistic matter (baryons and dark matter),
$\Omega_{R}$ the little fraction due to energy density in radiation and,
assuming zero spatial curvature131313CMB data from Planck combined with baryon
acoustic oscillations (BAO) give the constraint $\Omega_{K}=0.0007\pm 0.0019$
on the spatial curvature parameter, see [8]., $\Omega_{\Lambda}\simeq
1-\Omega_{M}$ is the dark energy contribution (in the form of a cosmological
constant). The Hubble parameter today is $H_{0}$.

In the dark energy model studied in this paper, the cosmological constant type
contribution to the energy density (due to the scalar field) can be read from
the first term of eq. (2.69) and it is

$\rho_{\Lambda}(N_{M})\simeq\frac{m^{2}}{2}\Delta_{\phi,\phi}(N_{M},0)\simeq\frac{m^{2}}{2}c_{s}^{-3}\frac{H_{I}^{2}}{32\pi^{2}|\xi|}e^{8|\xi|N_{I}+4|\xi|N_{M}}$
(2.99)

Today, after $N_{M}=\ln\left(\frac{\Omega_{M}}{\Omega_{R}}\right)$ e-foldings
from matter-radiation equality, this quantity has to match the amount of dark
energy density
$\rho_{\Lambda}(N_{M})=3M_{P}^{2}H_{0}^{2}\Omega_{\Lambda}=\frac{3}{8\pi}m_{P}^{2}H_{0}^{2}\Omega_{\Lambda}$,
where $m_{P}\equiv G^{-1/2}\simeq 1.2\times 10^{19}~{}{\rm GeV}$ is the Planck
mass in natural units and $M_{P}\equiv\left(8\pi G\right)^{-1/2}\simeq
2.4\times 10^{18}~{}{\rm GeV}$ is the reduced Planck mass.

This gives the required length of inflation in terms of the other parameters:

$N_{I}=\frac{1}{8|\xi|}\ln\biggl{[}24\pi|\xi|\Bigl{(}\frac{m_{P}}{H_{I}}\Bigr{)}^{2}\Bigl{(}\frac{H_{DE}}{m}\Bigr{)}^{2}c_{s}^{3}\biggr{]}-\frac{N_{M}}{2}\,,$
(2.100)

where $H_{\rm DE}\equiv H_{0}\sqrt{\Omega_{\Lambda}}$ is the Hubble parameter
when dark energy will completely dominate the cosmological expansion. This
situation will be realized provided that the CC-like term dominates over the
matter-like term in $\rho_{Q}$ at late times. From eq. (2.69), this
requirement imposes

$|\xi|<\frac{1}{6}\left(\frac{m}{H_{\rm DE}}\right)^{2}$ (2.101)

with the assumption $\xi<0$. The hypothesis that the scalar field is light,
which was used in all derivations and allows the best enhancement of quantum
fluctuations, is true at all stages of cosmological evolution (until late
times) if

$\frac{m}{H_{\rm DE}}<1\,.$ (2.102)

Eq. (2.100) refines the estimates in [31] and [1], which did not take into
account the growth of correlators in matter era and includes the effect of the
speed of sound $c_{s}$.

For $\Omega_{M}=0.3$ and $\Omega_{R}=9.1\times 10^{-5}$, then the number of
e-foldings of matter era (from matter-radiation equality until today) is
$N_{M}\simeq 8.1$. Therefore, looking at eq. (2.100) the growth in matter era
only decreases the required number of e-foldings by about $4$.

It is convenient to introduce the dimensionless parameter141414The condition
(2.102) also imposes $\alpha<\frac{1}{6|\xi|}$. Together with the inequality
in eq. (2.103), it means that, for a given value of $\xi<0$, the parameter
$\alpha$ can take values $1<\alpha<\frac{1}{6|\xi|}$. This requires
$|\xi|<\frac{1}{6}$. $\alpha$ defined as

$\boxed{\alpha\equiv\frac{1}{6|\xi|}\left(\frac{m}{H_{\rm
DE}}\right)^{2}>1}\,$ (2.103)

where the inequality follows from eq. (2.101) and we rewrite the result
(2.100) as

$\boxed{N_{I}=\frac{1}{8|\xi|}\ln\biggl{[}\frac{4\pi}{\alpha}\Bigl{(}\frac{m_{P}}{H_{I}}\Bigr{)}^{2}c_{s}^{3}\biggr{]}-\frac{1}{2}\ln\left(\frac{\Omega_{M}}{\Omega_{R}}\right)}\,.$
(2.104)

Then, using (2.99), the energy density today (cosmological time $t_{0}$) due
to quantum fluctuations of the scalar field is found from (2.69) to be

$\rho_{\rm Q}(t_{0})\simeq
3M_{P}^{2}H_{0}^{2}\left(\Omega_{\Lambda}-\frac{1}{\alpha}\right)\,.$ (2.105)

The rest of the energy density in the current Universe must be due to the the
classical energy density of matter $\rho_{C}(t_{0})$ that was used as a
background for the evolution of the scalar field, such that
$3M_{P}^{2}H_{0}^{2}=\rho_{C}(t_{0})+\rho_{\rm Q}(t_{0})$, giving

$\rho_{C}(t_{0})\simeq
3M_{P}^{2}H_{0}^{2}\left(\Omega_{M}+\frac{1}{\alpha}\right)\,.$ (2.106)

In Section 5, we study how a simple, but physically more consistent,
modelization of quantum backreaction, modifies eqs. (2.105)–(2.106), see in
particular eq. (5.24).

We conclude this section by discussing how the scale of spatial correlations
$r_{0}=\left(\mu a_{\rm in}H_{I}\right)^{-1}$, governing e.g. density
correlations in eq. (2.98), is related to the other cosmological/model
parameters. Assuming for simplicity an instantaneous reheating after inflation
and matching the inflationary Hubble parameter $H_{I}$ with that at the
beginning of the radiation era, we find the length scale of spatial variations
in units of the current Hubble horizon $H_{0}^{-1}$ (see also eq. (20) of
[1]):

$\boxed{r_{0}H_{0}\simeq\mu^{-1}e^{N_{I}}\left(\frac{H_{I}}{H_{0}}\right)^{-\frac{1}{2}}\Omega_{R}^{-\frac{1}{4}}=\mu^{-1}\left(\frac{H_{0}\sqrt{\Omega_{R}}}{H_{I}\Omega_{M}}\right)^{\frac{1}{2}}\biggl{[}\frac{4\pi}{\alpha}\Bigl{(}\frac{m_{P}}{H_{I}}\Bigr{)}^{2}c_{s}^{3}\biggr{]}^{\frac{1}{8|\xi|}}}\,.$
(2.107)

As an example, let us consider $\Omega_{M}=0.3$, $\Omega_{R}=9.1\times
10^{-5}$, $\Omega_{\Lambda}\simeq 0.7$, $H_{0}\simeq 10^{-33}~{}{\rm eV}$,
$H_{I}\simeq 10^{13}~{}{\rm GeV}$, $\xi=-0.06$ (which is in the ballpark of
the values used for relieving the Hubble tension in [1]), $\alpha\simeq 1$
(its minimum allowed value) and a speed of sound equal to speed of light
$c_{s}=1$. Then (using $m_{P}\simeq 1.22\cdot 10^{19}~{}{\rm GeV}$) the number
of inflationary e-foldings required from eq. (2.104) is $N_{I}\simeq 60$,
while (taking for definiteness $\mu=1$) the scale $r_{0}$ from eq. (2.107) is
$r_{0}\simeq 0.25H_{0}^{-1}$.

It is also interesting to comment on the effect of a reduced speed of sound
$c_{s}<1$, with respect to the base model with $c_{s}=1$. According to eq.
(2.104), this produces a reduction in the required number of e-foldings to
match the desired amount of dark energy today:

$N_{I}(c_{s})\approx
N_{I}(1)-\frac{3}{8|\xi|}\ln{\left(\frac{1}{c_{s}}\right)}\,.$ (2.108)

The origin of such a reduction is in the increased amplitude of quantum
fluctuations in inflation when $c_{s}<1$, ultimately determined by the mode
function (2.37), proportional to $c_{s}^{-\nu_{I}}\approx c_{s}^{-3/2}$, which
enters the noise sources $n_{(\rm 2)}$ for 2-pt correlators. A larger
amplitude when $c_{s}<1$ implies that a smaller number of inflationary
e-foldings is sufficient to reach the desired amount of dark energy, as eq.
(2.100) correctly predicts.

As a consequence, from eq. (2.107), the scale $r_{0}$ ruling the shape of
spatial correlations is affected by $c_{s}$ as

$r_{0}(c_{s})\simeq c_{s}^{\frac{3}{8|\xi|}}r_{0}(1)\,.$ (2.109)

This effect could be used to increase the detectability of the model by
studying fluctuations of luminosity distances from supernovae with surveys
like LSST. In [34] it has been shown that the amplitude of the angular power
spectrum for luminosity distances fluctuations is proportional to
$\left(r_{0}H_{0}\right)^{-16|\xi|}$. Using eq. (2.109), the effect of the
speed of sound on this amplitude is therefore

$\left(r_{0}(c_{s})H_{0}\right)^{-16|\xi|}\simeq
c_{s}^{-6}\left(r_{0}(1)H_{0}\right)^{-16|\xi|}\,,$ (2.110)

amounting to an amplification by $c_{s}^{-6}$, which would augment the chances
of detection.

##  3 Comparing stochastic formalism and quantum field theory

Stochastic formalism has been mainly studied in inflation, and its predictions
have been tested against a full quantum field theory (QFT) treatment only in
the coincident limit. Here we want to check the stochastic results against QFT
more generally, namely for non-coincident correlators and up to recent epochs.
Stochastic theory describes the evolution of super-Hubble modes and, as
expected, it agrees with QFT for super-Hubble separations. However, for
intermediate non-coincident separations there is no reason a priori to trust
the stochastic predictions and we find (in general) disagreement with QFT. As
we are dealing with a free spectator field (in inflation and in subsequent
epochs, except for very recent cosmological times, quantum gravitational
effects can be neglected and a classical treatment of gravity applies), the
stochastic formulation can be tested by comparing with one-loop QFT results.

As the starting point of the field theoretic analysis, let us consider the
non-coincident 2-point functions at equal times $\Delta_{\Phi\Phi}(t,r)$,
$\Delta_{\Phi\Pi}(t,r)$ and $\Delta_{\Pi\Pi}(t,r)$ defined from the full free
fields/canonical momenta (not their long-wavelength part) as

$\displaystyle\Delta_{\Phi\Phi}(t,r)$ $\displaystyle\equiv$
$\displaystyle\Bigl{\langle}\hat{\Phi}(t,\vec{x}_{1})\hat{\Phi}(t,\vec{x}_{2})\Bigr{\rangle}\,,$
$\displaystyle\Delta_{\Phi\Pi}(t,r)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{a^{3}(t)H(t)}\Bigl{\langle}\bigl{\\{}\hat{\Phi}(t,\vec{x}_{1}),\hat{\Pi}(t,\vec{x}_{2})\bigr{\\}}\Bigr{\rangle}\,,$
$\displaystyle\Delta_{\Pi\Pi}(t,r)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{a^{6}(t)H^{2}(t)}\Bigl{\langle}\hat{\Pi}(t,\vec{x}_{1})\hat{\Pi}(t,\vec{x}_{2})\Bigr{\rangle}\,,$
(3.1)

where $r\equiv\|\vec{x}_{2}-\vec{x}_{1}\|$ as usual. The correlators in (3.1)
can be obtained from Wightman functions defined as

$i\Delta^{(+)}\left(x;x^{\prime}\right)\equiv\langle\hat{\Phi}(x)\hat{\Phi}(x^{\prime})\rangle\equiv\left[i\Delta^{(-)}\left(x;x^{\prime}\right)\right]^{*}\,,$
(3.2)

by the following relations151515Both the positive- and negative-frequency
Wightman functions can be used for the equal-time correlators in (3.3) because
they are equal, at equal times, due to spatial homogeneity and isotropy of the
FLRW metric and the vacuum state $|\Omega\rangle$. From eq. (3.2) it follows
that
$\left[i\Delta^{(+)}\left(x;x^{\prime}\right)\right]_{t^{\prime}=t}=\left[i\Delta^{(-)}\left(x;x^{\prime}\right)\right]_{t^{\prime}=t}$
is real, which is consistent with eqs. (3.4).:

$\displaystyle\Delta_{\Phi\Phi}(t,r)$ $\displaystyle=$
$\displaystyle\left[i\Delta^{(\pm)}\left(x;x^{\prime}\right)\right]_{t^{\prime}=t}\,,$
$\displaystyle\Delta_{\Phi\Pi}(t,r)$ $\displaystyle=$
$\displaystyle\frac{1}{H(t)}\frac{\partial}{\partial
t}\left[i\Delta^{(\pm)}\left(x;x^{\prime}\right)\right]_{t^{\prime}=t}\,,$
$\displaystyle\Delta_{\Pi\Pi}(t,r)$ $\displaystyle=$
$\displaystyle\left[\frac{1}{H(t)}\frac{1}{H(t^{\prime})}\frac{\partial}{\partial
t}\frac{\partial}{\partial
t^{\prime}}i\Delta^{(\pm)}\left(x;x^{\prime}\right)\right]_{t^{\prime}=t}\,,$
(3.3)

where $x=\left(t,\vec{x}_{1}\right)$ and
$x^{\prime}=\left(t^{\prime},\vec{x}_{2}\right)$ are spacetime points.

The goal of this Section is to compare off-coincident correlators ($r\neq 0$)
computed with the stochastic formalism to those evaluated in QFT. Coincident
correlators ($r=0$) in stochastic formalism have already been compared with
the QFT prediction in [31, 32]. For the off-coincident correlators, it is
sufficient to work in $D=4$ spacetime dimensions, while the coincident ones
require to work in general spacetime dimensions $D$ to allow for dimensional
regularization of divergences. Therefore, as long as we focus on off-
coincident correlators, we work in $D=4$ spacetime dimensions as in previous
Sections of this paper. Only in subsection 3.2.1, when discussing the results
in the coincident limit based on [31, 32], we will use general $D$ spacetime
dimensions.

It is straightforward to express the correlators in (3.1) in terms of the mode
function $\varphi(t,k)$ using the decomposition in creation/annihilation
operators given in eqs. (2.10a–2.10b), to get

$\displaystyle\Delta_{\Phi\Phi}\left(t,r\right)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi^{2}}\int_{k_{0}}^{\infty}{\rm
d}k~{}k^{2}j_{0}(kr)|\varphi(t,k)|^{2}\,,$
$\displaystyle\Delta_{\Phi\Pi}\left(t,r\right)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi^{2}}\int_{k_{0}}^{\infty}{\rm
d}k~{}k^{2}j_{0}(kr)\frac{1}{H(t)}\frac{\partial}{\partial
t}|\varphi(t,k)|^{2}\,,$ $\displaystyle\Delta_{\Pi\Pi}\left(t,r\right)$
$\displaystyle=$ $\displaystyle\frac{1}{2\pi^{2}}\int_{k_{0}}^{\infty}{\rm
d}k~{}k^{2}j_{0}(kr)\Big{|}\frac{1}{H(t)}\dot{\varphi}(t,k)\Big{|}^{2}\,,$
(3.4)

where $j_{0}$ is the spherical Bessel functions of order $0$,
$j_{0}(z)=\sin(z)/z$. The infrared (IR) sector of the integral has been
regularized with a time-independent comoving cutoff $k_{0}$. We can anticipate
$k_{0}$ to be of the order of the inverse comoving Hubble horizon at the
beginning of inflation $k_{0}={\mathcal{O}}(1)a_{\rm in}H_{I}$. Smaller values
of the comoving wavenumber $k<k_{0}$ (i.e. modes with longer comoving
wavelength $\lambda=2\pi/k$) correspond to scales which are already super-
Hubble when inflation starts (during inflation the comoving Hubble sphere
shrinks, so longer modes are the first to exit the Hubble sphere).

For a comparison between between stochastic theory and quantum field theory
predictions, we will focus on the correlator
$\Delta_{\Phi,\Phi}\left(t,r\right)$. Following ref. [32], where the authors
studied the renormalized energy-momentum tensor for the same model (2.1), we
work in conformal time $\eta$, defined from cosmological time $t$ by the
relation $d\eta=dt/a(t)$ and rescale the mode function defining

$U(\eta,k)\equiv a(\eta)~{}\varphi\left(t(\eta),k\right)\,.$ (3.5)

The non-coincident correlator $\Delta_{\Phi,\Phi}$ written in the first of
eqs. (3.4) is then

$\Delta_{\Phi\Phi}\left(\eta,r\right)=\frac{1}{2\pi^{2}a^{2}(\eta)}\int_{k_{0}}^{\infty}{\rm
d}k~{}k^{2}j_{0}(kr)|U(\eta,k)|^{2}\,.$ (3.6)

Let us focus for the moment on a trivial speed of sound $c_{s}=1$. We will
discuss the effect of a reduced $c_{s}$ in subsection 3.2.4.

Starting from eq. (2.13) and denoting by “prime” (′) a derivative with respect
to conformal time, then $U(\eta,k)$ obeys the equation,

$U^{\prime\prime}+\left[k^{2}+m^{2}a^{2}-(1-6\xi)(2-\epsilon){\mathcal{H}}^{2}\right]U=0\,,$
(3.7)

and the Wronskian normalization (2.12) implies

$U~{}U^{\prime*}-U^{*}~{}U^{\prime}=i\,.$ (3.8)

In eqs. (3.7–3.8), the scale factor $a$, the reduced Hubble parameter
${\mathcal{H}}\equiv a^{\prime}/a=aH$ and the parameter
$\epsilon\equiv-\dot{H}/H^{2}=1-{\mathcal{H}}^{\prime}/{\mathcal{H}}^{2}$ are
functions of the conformal time $\eta$, while $U$ is function of $\eta$ and
$k$, as already said.

###  3.1 Mode functions

In the three epochs of de Sitter inflation, radiation domination and matter
domination, $\epsilon$ takes the constant values $0$, $2$ and $3/2$,
respectively161616When $\epsilon$ is constant, a useful relation for the
background evolution is ${\mathcal{H}}a^{\epsilon-1}={\rm constant}$, which
follows from integrating
$\epsilon=1-{\mathcal{H}}^{\prime}/{\mathcal{H}}^{2}$. In particular, in de
Sitter inflation ${\mathcal{H}}/a=H_{I}={\rm constant}$, in radiation epoch
${\mathcal{H}}a={\rm constant}$ and in matter epoch ${\mathcal{H}}^{2}a={\rm
constant}$.. In a period of constant $\epsilon$ the solutions for the mode
functions $U(\eta,k)$ are known exactly in the massless limit. The massless
approximation is definitely good for the very light non-minimally coupled
scalar we are dealing with, except for the very few most recent e-foldings in
matter era when the field backreacts on the FLRW background (leading to a dark
energy component, as we have seen). Furthermore, working in the massless limit
is sufficient for an interesting comparison between the stochastic formalism
and field theory, which is the goal of this Section. In this case, the mode
function $U_{\epsilon}(\eta,k)$ which solves eqs. (3.7–3.8) can be expressed
as a linear combination of the Chernikov-Tagirov-Bunch-Davies (CTBD) mode
functions in the form (see [32])

$U_{\epsilon}(\eta,k)=\alpha_{\epsilon}(k)u_{\epsilon}(\eta,k)+\beta_{\epsilon}(k)u_{\epsilon}^{*}(\eta,k)\,,$
(3.9)

where the CTBD mode function $u_{\epsilon}(\eta,k)$ is

$\displaystyle
u_{\epsilon}(\eta,k)=\sqrt{\frac{\pi}{4|1-\epsilon|\mathcal{H}(\eta)}}H_{\nu_{\epsilon}}^{(1)}\left(\frac{k}{|1-\epsilon|\mathcal{H}(\eta)}\right)\,,$
$\displaystyle{\rm
with}~{}~{}~{}\nu_{\epsilon}\equiv\sqrt{\frac{1}{4}+(1-6\xi)\frac{2-\epsilon}{(1-\epsilon)^{2}}}\,.$
(3.10)

$H_{\nu_{\epsilon}}^{(1)}$ is the Hankel function of the first kind with order
$\nu_{\epsilon}$.

The coefficients $\alpha_{\epsilon}(k)$ and $\beta_{\epsilon}(k)$ are known as
Bogolyubov coefficients. Since the CTBD mode function $u_{\epsilon}(\eta,k)$
in (3.10) satisfies the Wronskian condition
$u_{\epsilon}~{}u_{\epsilon}^{\prime*}-u_{\epsilon}^{*}~{}u_{\epsilon}^{\prime}=i$,
then the same condition (3.8) for the more general mode functions
$U_{\epsilon}(\eta,k)$ requires

$|\alpha_{\epsilon}(k)|^{2}-|\beta_{\epsilon}(k)|^{2}=1\,.$ (3.11)

The form of the Bogolyubov coefficients in different epochs depends on the
initial conditions at the the beginning of inflation, as well as on the
precise way one matches periods with different $\epsilon$, i.e. reheating and
matter-radiation transition. Following [1],[31], [32], we assume that during
inflation the full mode function is simply the CTBD mode function specialized
to de Sitter171717See also subsection III.C of [32] for interesting comments
on the initial state. ($\epsilon=0$),

$U_{I}(\eta,k)=\sqrt{\frac{\pi}{4{\mathcal{H}}(\eta)}}H_{\nu_{I}}^{(1)}\left(\frac{k}{{\mathcal{H}}(\eta)}\right)\,,$
(3.12)

where ${\mathcal{H}}(\eta)=a(\eta)H_{I}$, with $H_{I}$ the constant Hubble
parameter in de Sitter. We remind that, in the case of de Sitter inflation
($\epsilon=0$), the order of the Hankel function (3.10) is
$\nu_{0}=\sqrt{\frac{9}{4}-12\xi}$ and it can be shown that a mass $m$ can be
simply included by promoting this index to

$\nu_{I}\equiv\sqrt{\frac{9}{4}-12\xi-\left(\frac{m}{H_{I}}\right)^{2}}\,.$
(3.13)

As for the subsequent phases of radiation and matter domination, the
Bogolyubov coefficients were determined in [32] in the approximation of a
sudden transition between different $\epsilon$ epochs181818Limitations of this
approximation were discussed in refs. [29, 30].. Adopting a similar procedure
as in [32] (where the focus was on the computation of energy density and
pressure), in the massless limit191919The authors of ref. [32] also discussed
the exact massive CTBD mode function in radiation era (and an approximation
for it) as well as the approximation at
$\mathcal{O}\Bigl{(}\frac{ma}{\mathcal{H}}\Bigr{)}^{\\!2}$ of the mode
function in matter era. Here we only work in the massless limit. we arrive, by
imposing continuity of the mode function and its first derivative with respect
to conformal time, at the following expressions for the full mode functions in
radiation and matter epochs202020The result in radiation era (3.14) looks
simpler than in matter era (3.15) because in radiation era the non-minimal
coupling $\xi$ plays no role. Indeed, the order $\nu_{\epsilon}$ of the
massless CTBD mode function (3.10) reduces to 1/2 when $\epsilon=2$ and the
Hankel function of order 1/2 is $H_{1/2}^{(1)}(z)=-i\sqrt{\frac{2}{\pi
z}}~{}{\rm e}^{iz}$.:

$\displaystyle U_{R}(\eta,k)$ $\displaystyle\approx$ $\displaystyle
U_{I}(\eta_{\rm
end},k)\left[\left(\nu_{I}-\frac{1}{2}\right)\frac{{\mathcal{H}}_{\rm
end}}{k}\sin{\left(\frac{k}{{\mathcal{H}}(\eta)}\right)}-\left(\nu_{I}-\frac{3}{2}\right)\cos{\left(\frac{k}{{\mathcal{H}}(\eta)}\right)}\right]\,,$
(3.14) $\displaystyle U_{M}(\eta,k)$ $\displaystyle\approx$ $\displaystyle
U_{I}(\eta_{\rm
end},k)\left(\nu_{I}-\frac{1}{2}\right)\frac{{\mathcal{H}}_{\rm
end}}{{\mathcal{H}}_{\rm eq}}\left(\frac{{\mathcal{H}}_{\rm
eq}}{{\mathcal{H}}(\eta)}\right)^{\frac{1}{2}}\,$ (3.15)
$\displaystyle\times\frac{1}{2}\left[\left(\frac{3}{2}+\nu_{M}\right)\Gamma(\nu_{M})\left(\frac{{\mathcal{H}}_{\rm
eq}}{k}\right)^{\nu_{M}}J_{\nu_{M}}\left(\frac{2k}{{\mathcal{H}}(\eta)}\right)+\left(\nu_{M}\rightarrow-\nu_{M}\right)\right]\,,$

where $\nu_{M}$ in eq. (3.15) is the Hankel index of eq. (3.10) specialized to
matter era212121In the massless case, where eq. (3.10) is valid, matter era
and de Sitter inflation give the same order $\nu_{\epsilon}$ of the Hankel
function. This is true because $\epsilon$ affects $\nu_{\epsilon}$ only
through the combination $\frac{2-\epsilon}{(1-\epsilon)^{2}}$, which evaluates
to $2$ for both $\epsilon=0$ and $\epsilon=\frac{3}{2}$. More generally
$\epsilon$ and $\frac{3-2\epsilon}{2-\epsilon}$ give the same
$\nu_{\epsilon}$. An important difference is, however, that while in inflation
comoving modes exit the Hubble sphere, in matter era they enter it., namely

$\nu_{M}=\sqrt{\frac{9}{4}-12\xi}\,$ (3.16)

and $J_{\nu_{M}}$ is the Bessel function of order $\nu_{M}$. In the last line
we denoted by $\left(\nu_{M}\rightarrow-\nu_{M}\right)$ a contribution equal
to previous expression, but with $\nu_{M}$ replaced by $-\nu_{M}$.

In eqs. (3.14–3.15), ${\mathcal{H}}_{\rm end}$ and ${\mathcal{H}}_{\rm eq}$
are the reduced Hubble parameter at the end of inflation and at matter-
radiation equality, respectively. The matching with the inflationary mode
function at the end of inflation originates the factor $U_{I}(\eta_{\rm
end},k)$ in eqs. (3.14–3.15). For $k\ll{\mathcal{H}}_{\rm end}$, it can be
approximated as

$U_{I}(\eta_{\rm
end},k)\approx-\frac{i}{\sqrt{\pi}}2^{\nu_{I}-1}\Gamma(\nu_{I})\mathcal{H}_{\rm
end}^{\nu_{I}-1/2}k^{-\nu_{I}}\,.$ (3.17)

It is also useful to have at hand the mode functions approximations for long
modes (low $k$) because they will be useful for the infrared (IR) part of the
integrals in the next subsection. At leading order in $k/{\mathcal{H}}(\eta)$,
the three mode functions (3.12, 3.14, 3.15) reduce to

$\displaystyle U_{I}^{\rm(IR)}(\eta,k)$ $\displaystyle\approx$
$\displaystyle-\frac{i}{\sqrt{\pi}}2^{\nu_{I}-1}\Gamma(\nu_{I})\mathcal{H}(\eta)^{\nu_{I}-1/2}k^{-\nu_{I}}\,,$
(3.18) $\displaystyle U_{R}^{\rm(IR)}(\eta,k)$ $\displaystyle\approx$
$\displaystyle U_{I}(\eta_{\rm
end},k)\left[\left(\nu_{I}-\frac{1}{2}\right)\frac{{\mathcal{H}}_{\rm
end}}{{\mathcal{H}}(\eta)}-\left(\nu_{I}-\frac{3}{2}\right)\right]\,,$
$\displaystyle U_{M}^{\rm(IR)}(\eta,k)$ $\displaystyle\approx$ $\displaystyle
U_{I}(\eta_{\rm
end},k)\left(\nu_{I}-\frac{1}{2}\right)\frac{{\mathcal{H}}_{\rm
end}}{{\mathcal{H}}_{\rm
eq}}\frac{1}{2}\left[\left(1+\frac{3}{2\nu_{M}}\right)\left(\frac{{\mathcal{H}}_{\rm
eq}}{{\mathcal{H}}(\eta)}\right)^{\frac{1}{2}+\nu_{M}}+\left(\nu_{M}\rightarrow-\nu_{M}\right)\right]\,.$

From eq. (3.18), it is straightforward to check the continuity of the mode
function and of its first derivative with respect to conformal time in the
transition between two epochs.

###  3.2 The field 2-point function

Once the mode functions have been determined, we can move to the determination
of $\Delta_{\Phi\Phi}$ using eq. (3.6). Here we outline a strategy to compute
the required regularized integral. First one can formally write it as

$\Delta_{\Phi\Phi}\left(\eta,r\right)=\Delta_{\Phi\Phi}^{(\infty)}-\Delta_{\Phi\Phi}^{(k_{0})}\,,$
(3.19)

where we defined

$\displaystyle\Delta_{\Phi\Phi}^{(\infty)}\equiv\frac{1}{2\pi^{2}a^{2}(\eta)}\int_{k_{\rm
IR}}^{\infty}{\rm d}k~{}k^{2}j_{0}(kr)|U(\eta,k)|^{2}\,,$
$\displaystyle\Delta_{\Phi\Phi}^{(k_{0})}\equiv\frac{1}{2\pi^{2}a^{2}(\eta)}\int_{k_{\rm
IR}}^{k_{0}}{\rm d}k~{}k^{2}j_{0}(kr)|U(\eta,k)|^{2}\,.$ (3.20)

The quantity $k_{\rm IR}$, used as the lower extreme of integration in eq.
(3.20), is a deep IR cutoff, $k_{\rm IR}\ll k_{0}$. It regulates the IR
divergences, which would otherwise affect $\Delta_{\Phi\Phi}^{(\infty)}$ and
$\Delta_{\Phi\Phi}^{(k_{0})}$ if one had used $k=0$ as a lower extreme. The
divergence is there222222This can be seen from eq. (3.25). Since $j_{0}(0)=1$,
the integral (3.25) would diverge for $\nu_{I}\geq\frac{3}{2}$ if one had used
$k=0$ as a lower extreme instead of $k=k_{\rm IR}$. The condition
$\nu_{I}\geq\frac{3}{2}$ applies to our case $(m/H_{I})^{2}\ll|\xi|\ll 1$ with
$\xi<0$, giving from eq. (3.13), $\nu_{I}\simeq\frac{3}{2}+4|\xi|$. in the
case of a massless (or very light) scalar field with negative non-minimal
coupling $\xi<0$, which we are interested in, reducing eq. (3.13) to

$\nu_{I}\simeq\frac{3}{2}+4|\xi|\,.$ (3.21)

One can compute expressions for the integrals $\Delta_{\Phi\Phi}^{(\infty)}$
and $\Delta_{\Phi\Phi}^{(k_{0})}$, and subtract them to obtain the result of
(3.19), where the contributions from $k_{\rm IR}$ cancel exactly.

Let us see explicitly how this works during inflation. The expression of
$\Delta_{\Phi\Phi}^{(\infty)}$, using the mode function $U_{I}(\eta,k)$ in
(3.12), is

$\displaystyle\Delta_{\Phi\Phi,I}^{(\infty)}\left(\eta,r\right)$
$\displaystyle=$ $\displaystyle\frac{1}{2\pi^{2}a^{2}(\eta)}\int_{k_{\rm
IR}}^{\infty}{\rm d}k~{}k^{2}j_{0}(kr)|U_{I}(\eta,k)|^{2}\,$ (3.22)
$\displaystyle=$ $\displaystyle\frac{1}{8\pi
a^{2}(\eta){\mathcal{H}}(\eta)}\sqrt{\frac{\pi}{2r}}\int_{k_{\rm
IR}}^{\infty}{\rm
d}k~{}k^{\frac{3}{2}}J_{\frac{1}{2}}(kr)\bigg{|}H_{\nu_{I}}\left(\frac{k}{{\mathcal{H}}(\eta)}\right)\bigg{|}^{2}\,,$

where the subscript “$I$” in the first line stands for quantities computed
during inflation (thus ${\mathcal{H}}(\eta)=a(\eta)H_{I}$) and we used the
relation $j_{0}(z)=\sqrt{\frac{\pi}{2z}}J_{\frac{1}{2}}(z)$ between spherical
and ordinary Bessel functions.

By analytic continuation of eq. (6.578.10) in [58] or directly from eq. (27)
of [57], one finds for the integral in (3.22)

$\Delta_{\Phi\Phi,I}^{(\infty)}\left(\eta,r\right)=\frac{H_{I}^{2}}{16\pi^{2}}~{}\Gamma\left(\frac{3}{2}+\nu_{I}\right)\Gamma\left(\frac{3}{2}-\nu_{I}\right)~{}_{2}F_{1}\left(\frac{3}{2}+\nu_{I},\frac{3}{2}-\nu_{I};2;1-\frac{{\mathcal{H}}^{2}(\eta)r^{2}}{4}\right)-\Delta_{\rm
IR}\left(\eta,r\right)\,,$ (3.23)

where $\Delta_{\rm IR}\left(\eta,r\right)$ is the contribution due to the
lower extreme of integration $k_{\rm IR}$; it will be evaluated in eq. (3.28).

As for the infrared integral $\Delta_{\Phi\Phi,I}^{(k_{0})}$ in inflation, one
can use the low-$k$ approximation for the mode function (3.18) instead of the
full result (3.12) because $k_{0}\ll{\mathcal{H}(\eta)}$, to get232323 More
precisely, since $k_{0}={\mathcal{O}}(1){\mathcal{H}}_{\rm in}$ is of the
order of the inverse comoving Hubble radius at the beginning of inflation, the
relation $k_{0}\ll{\mathcal{H}(\eta)}$ is not valid in the very early stages
of inflation, but after a few e-foldings
$N(\eta)=\ln{\left(\frac{a(\eta)}{a_{\rm in}}\right)}$, one can say that
$\frac{k_{0}}{{\mathcal{H}}(\eta)}={\mathcal{O}}(1)~{}{\rm e}^{-N(\eta)}\ll
1$, thus decaying exponentially with $N(\eta)$, where we used that
${\mathcal{H}(\eta)}=a(\eta)H_{I}$ in inflation.

$\Delta_{\Phi\Phi,I}^{(k_{0})}\left(\eta,r\right)=\frac{2^{2\nu_{I}-3}\Gamma^{2}(\nu_{I}){\mathcal{H}}(\eta)^{2\nu_{I}-1}}{\pi^{3}a^{2}(\eta)}\int_{k_{\rm
IR}}^{k_{0}}{\rm d}k~{}k^{2-2\nu_{I}}j_{0}(kr)\,.$ (3.24)

Let us rewrite the integral appearing in (3.24) in terms of a dimensionless
integration variable $z=kr$ as

$\int_{k_{\rm IR}}^{k_{0}}{\rm
d}k~{}k^{2-2\nu_{I}}j_{0}(kr)=r^{2\nu_{I}-3}\int_{k_{\rm IR}r}^{k_{0}r}{\rm
d}z~{}z^{2-2\nu_{I}}j_{0}(z)\,.$ (3.25)

Now we use the formula242424Eq. (3.26) can be checked by the series expansions
$j_{0}(z)={\sum}_{n=0}^{\infty}\frac{\left(-z^{2}\right)^{n}}{(2n+1)!}$ and
$~{}_{1}F_{2}\left(\frac{3}{2}-\nu_{I};\frac{3}{2},\frac{5}{2}-\nu_{I};-\frac{z^{2}}{4}\right)={\sum}_{n=0}^{\infty}\frac{\left(\frac{3}{2}-\nu_{I}\right)_{n}}{\left(\frac{3}{2}\right)_{n}\left(\frac{5}{2}-\nu_{I}\right)_{n}}\frac{1}{n!}\left(-\frac{z^{2}}{4}\right)^{n}$,
where $(c)_{n}=\frac{\Gamma(c+n)}{\Gamma(c)}$ is the Pochhammer symbol.

$\int_{z_{\rm IR}}^{z_{0}}{\rm
d}z~{}z^{2-2\nu_{I}}j_{0}(z)=\frac{z_{0}^{3-2\nu_{I}}}{3-2\nu_{I}}~{}_{1}F_{2}\left(\frac{3}{2}-\nu_{I};\frac{3}{2},\frac{5}{2}-\nu_{I};-\frac{z_{0}^{2}}{4}\right)-\frac{z_{\rm
IR}^{3-2\nu_{I}}}{3-2\nu_{I}}~{}_{1}F_{2}\left(\frac{3}{2}-\nu_{I};\frac{3}{2},\frac{5}{2}-\nu_{I};-\frac{z_{\rm
IR}^{2}}{4}\right)\,,$ (3.26)

and substitute $z_{\rm IR}=k_{\rm IR}r$ and $z_{0}=k_{0}r$, to obtain

$\displaystyle\Delta_{\Phi\Phi,I}^{(k_{0})}\left(\eta,r\right)=\frac{H_{I}^{2}}{\pi^{3}}\frac{\Gamma^{2}(\nu_{I})}{3-2\nu_{I}}\\!\\!\\!\\!\\!\\!$
$\displaystyle\left[\left(\frac{2{\mathcal{H}}(\eta)}{k_{0}}\right)^{2\nu_{I}-3}~{}_{1}F_{2}\left(\frac{3}{2}-\nu_{I};\frac{3}{2},\frac{5}{2}-\nu_{I};-\frac{k_{0}^{2}r^{2}}{4}\right)\right.\,$
(3.27) $\displaystyle~{}\left.-\left(\frac{2{\mathcal{H}}(\eta)}{k_{\rm
IR}}\right)^{2\nu_{I}-3}~{}_{1}F_{2}\left(\frac{3}{2}-\nu_{I};\frac{3}{2},\frac{5}{2}-\nu_{I};-\frac{k_{\rm
IR}^{2}r^{2}}{4}\right)\right]\,.$

The second line of eq. (3.27) contains the deep IR cutoff $k_{\rm IR}$ and it
corresponds exactly to the quantity $\Delta_{\rm IR}\left(\eta,r\right)$
appearing in eq. (3.23), in the sense that

$\Delta_{\rm
IR}\left(\eta,r\right)=\frac{H_{I}^{2}}{\pi^{3}}\frac{\Gamma^{2}(\nu_{I})}{3-2\nu_{I}}\left(\frac{2{\mathcal{H}}(\eta)}{k_{\rm
IR}}\right)^{2\nu_{I}-3}~{}_{1}F_{2}\left(\frac{3}{2}-\nu_{I};\frac{3}{2},\frac{5}{2}-\nu_{I};-\frac{k_{\rm
IR}^{2}r^{2}}{4}\right)\,.$ (3.28)

Subtracting eq. (3.27) from eq. (3.23), as specified in (3.19), the deep IR
term $\Delta_{\rm IR}\left(\eta,r\right)$ cancels out and one finally finds

$\displaystyle\Delta_{\Phi\Phi,I}\left(\eta,r\right)$ $\displaystyle=$
$\displaystyle-\Delta_{\Phi\Phi,I}^{(k_{0})}\left(\eta,r\right)+\Delta_{\Phi\Phi,I}^{(\infty)}\left(\eta,r\right)\,$
$\displaystyle=$
$\displaystyle\frac{H_{I}^{2}}{\pi^{2}}\left[\frac{\Gamma^{2}(\nu_{I})}{(2\nu_{I}-3)\pi}\left(\frac{2{\mathcal{H}}(\eta)}{k_{0}}\right)^{2\nu_{I}-3}~{}_{1}F_{2}\left(\frac{3}{2}-\nu_{I};\frac{3}{2},\frac{5}{2}-\nu_{I};-\frac{k_{0}^{2}r^{2}}{4}\right)\right.\,$
$\displaystyle\left.+\frac{1}{16}~{}\Gamma\left(\frac{3}{2}+\nu_{I}\right)\Gamma\left(\frac{3}{2}-\nu_{I}\right)~{}_{2}F_{1}\left(\frac{3}{2}+\nu_{I},\frac{3}{2}-\nu_{I};2;1-\frac{{\mathcal{H}}^{2}(\eta)r^{2}}{4}\right)\right]\,.\,$

This is the result during inflation for the field 2-point function (3.6).

####  3.2.1 Coincident correlator

Let us consider the coincident limit $r=0$ in inflation. The coincident
correlator has been computed within the stochastic approximation in [31] and
compared to the field theory result of [32]. Here we summarize the QFT result,
referring the reader to [32] for more details. In the coincident limit, one
cannot work in $D=4$ spacetime dimensions from the beginning, but it is
required to work with general $D$ spacetime dimensions in order to regularize
divergences (via dimensional regularization). As a starting point for
understanding this, let us consider eq. (3.2), which was obtained in $D=4$,
and specialize it to $r=0$. Using ${\mathcal{H}}(\eta)=a(\eta)H_{I}$, eq.
(3.2) reduces to

$\displaystyle\Delta_{\Phi\Phi,I}\left(\eta,0\right)$ $\displaystyle=$
$\displaystyle\frac{H_{I}^{2}}{\pi^{2}}\left[\frac{\Gamma^{2}(\nu_{I})}{(2\nu_{I}-3)\pi}\left(\frac{2a_{\rm
in}H_{I}}{k_{0}}\right)^{2\nu_{I}-3}\left(\frac{a(\eta)}{a_{\rm
in}}\right)^{2\nu_{I}-3}\right.\,$ (3.30)
$\displaystyle\left.+\frac{1}{16}~{}\Gamma\left(\frac{3}{2}+\nu_{I}\right)\Gamma\left(\frac{3}{2}-\nu_{I}\right)~{}_{2}F_{1}\left(\frac{3}{2}+\nu_{I},\frac{3}{2}-\nu_{I};2;1\right)\right]\,,$

where, in the first line of (3.30), we multiplied and divided by the same
power of the scale factor at the beginning of inflation $a_{\rm in}$. The
second line of eq. (3.30) originates from
$\Delta_{\Phi\Phi,I}\left(\eta,0\right)$ and it is a divergent constant, as
one can see from Gauss’ identity
$~{}_{2}F_{1}(a,b;c;1)=\frac{\Gamma(c)\Gamma(c-a-b)}{\Gamma(c-a)\Gamma(c-b)}$
because
$c-a-b=2-\left(\frac{3}{2}+\nu_{I}\right)-\left(\frac{3}{2}-\nu_{I}\right)=-1$,
which is a simple pole of the gamma function. This divergence can be cured by
dimensional regularization, working in $D$ spacetime dimensions instead of
$D=4$ and promoting $D$ to a complex variable. In general $D$ spacetime
dimensions, one can prove[59] (see also eq. (30) of [57]) that the result
(3.23) for $\Delta_{\Phi\Phi,I}^{(\infty)}\left(\eta,r\right)$ (omitting to
write the deep IR contribution $\Delta_{\rm IR}(\eta,r)$, which cancels out in
the final expression (3.2) for $\Delta_{\Phi\Phi,I}\left(\eta,r\right)$)
generalizes to

$\displaystyle\Delta_{\Phi\Phi,I}^{(\infty)}\left(\eta,r\right)$
$\displaystyle=$
$\displaystyle\frac{H_{I}^{D-2}}{(4\pi)^{\frac{D}{2}}}\frac{\Gamma\left(\frac{D-1}{2}+\nu_{I,D}\right)\Gamma\left(\frac{D-1}{2}-\nu_{I,D}\right)}{\Gamma\left(\frac{D}{2}\right)}\,$
(3.31)
$\displaystyle\times~{}_{2}F_{1}\left(\frac{D-1}{2}+\nu_{I,D},\frac{D-1}{2}-\nu_{I,D};\frac{D}{2};1-\frac{{\mathcal{H}}^{2}(\eta)r^{2}}{4}\right)\,,$

with $\nu_{I,D}\equiv\sqrt{\left(\frac{D-1}{2}\right)^{2}-D(D-1)\xi}$ for a
massless non-minimally coupled field. In the coincident limit $r=0$, by
applying Gauss’ identity, one gets

$\Delta_{\Phi\Phi,I}^{(\infty)}\left(\eta,0\right)=\frac{H_{I}^{D-2}}{(4\pi)^{\frac{D}{2}}}\frac{\Gamma\left(\frac{D-1}{2}+\nu_{I,D}\right)\Gamma\left(\frac{D-1}{2}-\nu_{I,D}\right)\Gamma\left(1-\frac{D}{2}\right)}{\Gamma\left(\frac{1}{2}+\nu_{I,D}\right)\Gamma\left(\frac{1}{2}-\nu_{I,D}\right)}\,,$
(3.32)

and the factor $\Gamma\left(1-\frac{D}{2}\right)$ is responsible for the
divergence at $D=4$, having a simple pole there. Expanding around $D=4$ and
introducing a renormalization energy scale $\bar{\mu}$, one gets

$\displaystyle\Delta_{\Phi\Phi,I}^{(\infty)}\left(\eta,0\right)=-\frac{H_{I}^{2}}{8\pi^{2}}\left(1-6\xi\right)\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
$\displaystyle\left[2\frac{\bar{\mu}^{D-4}}{D-4}+\ln\left(\frac{H_{I}^{2}}{4\pi\bar{\mu}^{2}}\right)+\gamma_{E}-\frac{\xi}{1-6\xi}+\psi\left(\frac{1}{2}+\nu_{I}\right)+\psi\left(\frac{1}{2}-\nu_{I}\right)\right.\,$
(3.33) $\displaystyle\left.+{\mathcal{O}}\left(D-4\right)\right]\,,$

where $\nu_{I}$ in eq. (3.33) refers to the value in $D=4$ spacetime
dimensions, $\nu_{I}\equiv\nu_{I,4}=\sqrt{\frac{9}{4}-12\xi}$. The function
$\psi(z)\equiv\Gamma^{\prime}(z)/\Gamma(z)$ is the digamma function and
$\gamma_{E}=-\psi(1)\approx 0.577$ is the Euler-Mascheroni constant. The
divergence at $D=4$ is removed by the addition of suitable counterterms so
that the coincident correlator
$\Delta_{\Phi\Phi,I}^{(\infty)}\left(\eta,0\right)$ in eq. (3.33) reduces to a
finite constant, which is then added to the contribution in the first line of
(3.30), coming from $-\Delta_{\Phi\Phi,I}^{(k_{0})}\left(\eta,0\right)$, to
get a final result for the full renormalized252525In fact one never
renormalizes the propagator, but rather physical quantities. In general, the
renormalized propagator cannot be used for full loop calculations, except for
some specific cases. One such example is the expectation value of the energy-
momentum tensor, where the coincident Feynman propagator appears:
$\langle\hat{T}_{\mu\nu}\rangle\supset-\frac{m^{2}}{2}g_{\mu\nu}\left[i\Delta_{F}(x;x)\right]^{(\rm
ren)}$. This quantity can be renormalized by a cosmological constant
counterterm. Here, the renormalized coincident correlator is considered for
the purpose of making comparison with the stochastic 2-point functions of
Section 2. coincident correlator $\Delta_{\Phi\Phi,I}^{(\rm
ren)}\left(\eta,0\right)$. For comparison with the prediction of stochastic
formalism, we choose the (non-minimal) subtraction scheme in such a way that
$\Delta_{\Phi\Phi,I}^{(\rm ren)}\left(\eta_{\rm in},0\right)$ is zero at the
beginning of inflation $(\eta=\eta_{\rm in})$. Then the time evolution of the
renormalized correlator $\Delta_{\Phi\Phi,I}^{(\rm ren)}\left(\eta,0\right)$
is given by

$\Delta_{\Phi\Phi,I}^{(\rm
ren)}\left(\eta,0\right)=H_{I}^{2}\frac{\Gamma^{2}(\nu_{I})}{(2\nu_{I}-3)\pi^{3}}\left(\frac{2a_{\rm
in}H_{I}}{k_{0}}\right)^{2\nu_{I}-3}\left[\left(\frac{a(\eta)}{a_{\rm
in}}\right)^{2\nu_{I}-3}-1\right]\,.$ (3.34)

For $\nu_{I}>\frac{3}{2}$, as in the case $\xi<0$ where (3.21) holds, the
contribution $\left(\frac{a(\eta)}{a_{\rm
in}}\right)^{2\nu_{I}-3}\approx\left(\frac{a(\eta)}{a_{\rm
in}}\right)^{8|\xi|}$ originated from
$-\Delta_{\Phi\Phi,I}^{(k_{0})}\left(\eta,0\right)$ is amplified by the growth
of the scale factor. Therefore, after a few e-foldings
$N_{I}(\eta)=\ln{\left(\frac{a(\eta)}{a_{\rm in}}\right)}$ of inflation, a
good approximation for eq. (3.34) is

$\Delta_{\Phi\Phi,I}^{(\rm
ren)}\left(\eta,0\right)\approx\frac{H_{I}^{2}}{32\pi^{2}|\xi|}\left(\frac{2a_{\rm
in}H_{I}}{k_{0}}\right)^{8|\xi|}{\rm e}^{8|\xi|N_{I}(\eta)}\,,$ (3.35)

where we used $\nu_{I}\simeq\frac{3}{2}+4|\xi|$ and worked at leading order in
$|\xi|\ll 1$. Comparing the quantum field theory result (3.35) with the
stochastic prediction (2.53) at coincidence (and with $c_{s}=1$), we find the
same growth with $N_{I}$. We also find the same amplitude for $k_{0}=2a_{\rm
in}H_{I}$, which is compatible with the expectation that the IR cutoff $k_{0}$
should be of the order of the inverse comoving Hubble radius at the beginning
of inflation $k_{0}={\mathcal{O}}(1)a_{\rm in}H_{I}$ (as discussed below eq.
(3.4)).

By applying the same logic, the coincident correlators in radiation epoch,
$\Delta_{\Phi\Phi,R}\left(\eta,0\right)$, and in matter epoch,
$\Delta_{\Phi\Phi,M}\left(\eta,0\right)$, are also dominated by their infrared
parts, namely $-\Delta_{\Phi\Phi,R}^{(k_{0})}\left(\eta,0\right)$ (and
$-\Delta_{\Phi\Phi,M}^{(k_{0})}\left(\eta,0\right)$). They can be evaluated
using the low-$k$ mode functions in the second and third line of eq. (3.18).
The result is that, in radiation epoch, the coincident correlator stays
roughly constant, as we find for the ratio with respect to its initial value
(i. e. at the end of inflation, at conformal time $\eta_{\rm end}$),

$\frac{\Delta_{\Phi\Phi,R}\left(\eta,0\right)}{\Delta_{\Phi\Phi,I}\left(\eta_{\rm
end},0\right)}=\left[\nu_{I}-\frac{1}{2}-\left(\nu_{I}-\frac{3}{2}\right){\rm
e}^{-N_{R}(\eta)}\right]^{2}\,,$ (3.36)

where $N_{R}(\eta)$ is the number of e-foldings in radiation epoch measured
with respect to the end of inflation, so that $N_{R}(\eta_{\rm end})=0$. For
$\nu_{I}\simeq\frac{3}{2}+4|\xi|$ the ratio in eq. (3.36) only grows from $1$
at the beginning of radiation era ($N_{R}=0$) to
$\left(\nu_{I}-\frac{1}{2}\right)^{2}\simeq 1+8|\xi|$ for large value of
$N_{R}$ (radiation epoch lasts $\approx 50$ e-foldings), which means that
$\Delta_{\Phi\Phi,R}\left(\eta,0\right)$ is approximately constant, in
agreement with the stochastic result (2.61).

With the same strategy we find that in matter era, the coincident correlator
grows with respect to its initial value at matter-radiation equality (which
corresponds to conformal time $\eta_{\rm eq}$) by a factor

$\frac{\Delta_{\Phi\Phi,M}\left(\eta,0\right)}{\Delta_{\Phi\Phi,R}\left(\eta_{\rm
eq},0\right)}=\frac{1}{4}\left[1+\frac{3}{2\nu_{M}}+\left(1-\frac{3}{2\nu_{M}}\right){\rm
e}^{-\nu_{M}N_{M}(\eta)}\right]^{2}~{}{\rm
e}^{\left(\nu_{M}-\frac{3}{2}\right)N_{M}(\eta)}\,,$ (3.37)

where $N_{M}(\eta)$ is the number of e-foldings in matter epoch measured
starting from matter-radiation equality, meaning that $N_{M}(\eta_{\rm
eq})=0$. Using $\nu_{I}\simeq\frac{3}{2}+4|\xi|$ (leading order in $|\xi|$ of
(3.16) for $\xi<0$), then eq. (3.37) is well approximated by simply ${\rm
e}^{4|\xi|N_{M}(\eta)}$, in agreement with the stochastic formalism prediction
in eq. (2.68).

Thus quantum field theory confirms the one-loop prediction of the stochastic
formalism in the coincident limit.

####  3.2.2 Super-Hubble regime

Stochastic theory, following the dynamics of super-Hubble modes, is expected
to be trustworthy for super-Hubble separations $r$, i.e. when
${\mathcal{H}}(\eta)r\gg 1$. Here we want to check its validity in this regime
against the quantum field theory prediction
$\Delta_{\Phi\Phi,I}\left(\eta,r\right)=-\Delta_{\Phi\Phi,I}^{(k_{0})}\left(\eta,r\right)+\Delta_{\Phi\Phi,I}^{(\infty)}\left(\eta,r\right)$
that we worked out in eq. (3.2).

For this purpose, in inflation it is useful to rewrite the contribution from
$\Delta_{\Phi\Phi,I}^{(\infty)}\left(\eta,r\right)$ in eq. (3.23), using the
hypergeometric transformation formula (9.132.1) in [58] (and omitting to write
the deep IR contribution $\Delta_{\rm IR}(\eta,r)$, which cancels out in the
final expression for $\Delta_{\Phi\Phi,I}\left(\eta,r\right)$ in eq. (3.2)),
to get

$\displaystyle\Delta_{\Phi\Phi,I}^{(\infty)}\left(\eta,r\right)=\frac{H_{I}^{2}}{16\pi^{2}}\left[\frac{\Gamma\left(\frac{3}{2}-\nu_{I}\right)\Gamma\left(2\nu_{I}\right)}{\Gamma\left(\frac{1}{2}+\nu_{I}\right)}\left(\frac{{\mathcal{H}}(\eta)r}{2}\right)^{2\nu_{I}-3}~{}_{2}F_{1}\left(\frac{3}{2}-\nu_{I},\frac{1}{2}-\nu_{I};1-2\nu_{I};\frac{4}{{\mathcal{H}}^{2}(\eta)r^{2}}\right)\right.\,$
$\displaystyle\left.+\frac{\Gamma\left(\frac{3}{2}+\nu_{I}\right)\Gamma\left(-2\nu_{I}\right)}{\Gamma\left(\frac{1}{2}-\nu_{I}\right)}\left(\frac{{\mathcal{H}}(\eta)r}{2}\right)^{-2\nu_{I}-3}~{}_{2}F_{1}\left(\frac{3}{2}+\nu_{I},\frac{1}{2}+\nu_{I};1+2\nu_{I};\frac{4}{{\mathcal{H}}^{2}(\eta)r^{2}}\right)\right]\,.$
(3.38)

For super-Hubble distances with ${\mathcal{H}}(\eta)r\gg 1$ and for
$\nu_{I}\simeq\frac{3}{2}+4|\xi|$, the leading contribution to
$\Delta_{\Phi,\Phi}^{(\infty)}$ comes from the power
$\left(\frac{{\mathcal{H}}(\eta)r}{2}\right)^{2\nu_{I}-3}$ in the first line
of (3.38) (the hypergeometric functions are close to $1$ for small arguments
$\frac{4}{{\mathcal{H}}^{2}(\eta)r^{2}}$), resulting in

$\displaystyle\Delta_{\Phi\Phi,I}^{(\infty)}\left(\eta,r\right)\simeq-\frac{H_{I}^{2}}{32\pi^{2}|\xi|}\left(\frac{{\mathcal{H}}(\eta)r}{2}\right)^{8|\xi|}=-\frac{H_{I}^{2}}{32\pi^{2}|\xi|}\left(\frac{a_{\rm
in}H_{I}r}{2}\right)^{8|\xi|}{\rm e}^{8|\xi|N_{I}(\eta)}\,.$ (3.39)

Adding the contribution from
$-\Delta_{\Phi\Phi,I}^{(k_{0})}\left(\eta,r\right)$ (see eq. (3.27)) and using
$k_{0}=2a_{\rm in}H_{I}$ which we required from matching the coincident
correlators from quantum field theory and stochastic formalism (see the
discussion below eq. (3.35)), one finds262626The function
$~{}_{1}F_{2}\left(\frac{3}{2}-\nu_{I};\frac{3}{2},\frac{5}{2}-\nu_{I};-\frac{k_{0}^{2}r^{2}}{4}\right)$
in eq. (3.27) can be approximated as $1$ even at super-Hubble scales
(${\mathcal{H}}(\eta)r\gg 1$) and not only in the coincident regime, because
its argument $-\frac{k_{0}^{2}r^{2}}{4}$ is still small. The IR cutoff $k_{0}$
sets an upper limit for comoving scales $r\lesssim 1/k_{0}$. Another way to
understand the smallness of $-\frac{k_{0}^{2}r^{2}}{4}$ is to use
$k_{0}=2a_{\rm in}H_{I}=2a(\eta)H_{I}{\rm
e}^{-N(\eta)}=2{\mathcal{H}}(\eta){\rm e}^{-N(\eta)}$ and therefore
$-\frac{k_{0}^{2}r^{2}}{4}=-\left({\mathcal{H}}(\eta)r\right)^{2}{\rm
e}^{-2N(\eta)}$. This is suppressed by ${\rm e}^{-2N(\eta)}$, so that for
$N(\eta)$ large enough, $-\frac{k_{0}^{2}r^{2}}{4}$ is small even for large
values $1\ll{\mathcal{H}}(\eta)r\ll{\rm e}^{2N(\eta)}$.

$\Delta_{\Phi\Phi,I}\left(\eta,r\right)=-\Delta_{\Phi\Phi,I}^{(k_{0})}\left(\eta,r\right)+\Delta_{\Phi\Phi,I}^{(\infty)}\left(\eta,r\right)\simeq\frac{H_{I}^{2}}{32\pi^{2}|\xi|}~{}{\rm
e}^{8|\xi|N_{I}(\eta)}\left[1-\left(\frac{a_{\rm
in}H_{I}r}{2}\right)^{8|\xi|}\right]\,.$ (3.40)

Eq. (3.40) agrees with the stochastic prediction in eqs. (2.53)–(2.54), whose
spatial profile is proportional to $1-\left(\mu a_{\rm
in}H_{I}r\right)^{8|\xi|}$ where $\mu$ sets the UV cutoff of the stochastic
formalism. The two results have the same amplitude (for $k_{0}=2a_{\rm
in}H_{I}$) and a profile described by a power law, proportional to
$1-\left(\frac{r}{r_{0}}\right)^{8|\xi|}$. In both cases $r_{0}$ is of the
order of $\left(a_{\rm in}H_{I}\right)^{-1}$, which is the comoving Hubble
horizon at the beginning of inflation.

####  3.2.3 Full comparison

We have seen that stochastic formalism is able to reproduce QFT results in the
coincident limit (subsection 3.2.1) and for super-Hubble separations
(subsection 3.2.2). At intermediate non-coincident separations (sub-Hubble and
near-Hubble scales) there is no reason to trust stochastic formalism. Hence it
is interesting to compare its predictions and the field theory results. The
stochastic formalism predicts that the spatial dependence (normalized with
respect to the value at coincidence) of the field 2-point function inherited
from inflation is not substantially altered in radiation and matter epochs and
therefore it is the same today as it was at the end of inflation. This was a
consequence of the negligibility of stochastic sources after the end of
inflation, so that the time evolution could only modify the amplitude (in
stochastic formalism, $\Delta_{\phi\phi}$ grows in matter era, while it is
roughly constant in radiation era), but not the spatial dependence of
correlators. In Fig. 1 we compare this stochastic prediction to the full
quantum field theory result, highlighting how, according to field theory, the
spatial profile gets modified today (assuming the matter domination
approximation) with respect to the end of inflation. The stochastic prediction
for the spatial profile works much better at the end of inflation and in
radiation epoch than in matter epoch. In other words the stochastic formalism
fails to predict the right spatial dependence of field correlators in matter
era and today. This can be largely understood as an effect of the non-
negligibility of spatial gradients at late times in the life of the Universe,
which were dropped by the stochastic formalism calculation in Section 2. For a
given Fourier mode, with comoving momentum $k$, the effect of spatial
gradients on the equation of motion of the field is determined at conformal
time $\eta$ by $\left(\frac{k}{{\mathcal{H}}(\eta)}\right)^{2}$. Since the
comoving Hubble radius ${\mathcal{H}}^{-1}(\eta)$ was much smaller at the end
of inflation than today, the gradients play a largely more relevant role today
than in the early Universe. This helps us to understand the problem with
stochastic formalism today (or in general at late times), because a steep
decrease in the field correlator like the one predicted by stochastic theory
(green curve in Fig. 1) would imply a large energy of the field configuration
associated to the spatial gradients and it is therefore energetically
disfavored. On the contrary, quantum field theory predicts today (orange curve
in Fig. 1) a softer decrease of correlations with distance, thus diminishing
the energy content of the field configuration. We also note that the coherent
oscillations on super-Hubble scales predicted by QFT are absent in the
stochastic approximation.

####  3.2.4 Effect of a reduced speed of sound

Stochastic theory fails to reproduce the correct spatial dependence of the
field correlator at sub-Hubble scales in matter epoch, as shown by Fig. 2.
However, the discussion in subsection 3.2.3 contains a hint on how agreement
could be restored even at late times (today). Based on the role of spatial
gradients, a simple strategy consists in decreasing their effect via the
introduction of a reduced speed of sound. Such a situation is common in many
studied and motivated modified gravity theories, e. g. scalar-tensor theories
in the Horndeski class.

Let us introduce a constant speed of sound $c_{s}\leq 1$. The equation of
motion for the quantum field $\hat{\Phi}(t,\vec{x})$ is given in (2.9) and the
equation of motion (3.7) for the mode function $U\left(\eta,k\right)$ gets
simply modified into

$U^{\prime\prime}+\left[c_{s}^{2}k^{2}+m^{2}a^{2}-(1-6\xi)(2-\epsilon){\mathcal{H}}^{2}\right]U=0\,,$
(3.41)

i.e. $k$ is replaced by $c_{s}k$. One can then repeat the calculations of this
Section 3 to predict the spatial profile of the field 2-point function in the
presence of a constant reduced speed of sound. The result is plotted in Fig.
(2) and it confirms the intuition that a $c_{s}<1$ enhances the stochastic
formalism performance, reconciling it with quantum field theory.

Figure 1:  The field correlator $\Delta_{\Phi\Phi}$ computed with quantum
field theory, as a function of distance. Distances are measured in units of
the current Hubble horizon $H_{0}^{-1}$, thus a value of $1$ on the horizontal
axis stands for a separation whose physical size today is the Hubble horizon
$H_{0}^{-1}$. Correlators are normalized with respect to their respective
value at coincidence. For simplicity and better readability of the plot, we
considered three key moments in the history of the Universe: end of inflation
(red), matter-radiation equality (blue) and today (orange). The stochastic
formalism prediction (green) does not change relevantly from the end of
inflation until today. We see that, at sub-Hubble scales, the stochastic
theory only agrees with quantum field theory at early times (from the end of
inflation until matter-radiation equality), but it fails at later times in
matter-dominated era, whose result has been extrapolated until today. The plot
assumes a sudden transition between different epochs and the values chosen for
the cosmological parameters are: $\Omega_{M}=0.3$, $\Omega_{R}=9.1\times
10^{-5}$ (matter and radiation energy density fractions today),
$H_{0}=10^{-33}~{}{\rm eV}$ (Hubble parameter today) and $H_{I}=10^{13}~{}{\rm
GeV}$ (Hubble parameter in inflation). They imply a length of radiation epoch
of $N_{R}\simeq 57$ e-foldings and the matter era started $N_{M}\simeq 8.1$
e-foldings ago from now. We worked in the massless limit for the field, with a
negative non-minimal coupling $\xi=-0.06$. For simplicity, the duration of
inflation assumed is $N_{I}=61$ e-foldings, such that, given the other
parameters listed before, the comoving Hubble length at the beginning of
inflation is equal to the Hubble horizon today. Note that this is not
necessarily the realistic number of e-foldings needed to reproduce the right
amount of dark energy today: the correct value also depends on the mass $m$ of
the scalar field and the speed of sound $c_{s}$, as in eq. (2.100). In this
plot, we assumed a speed of sound $c_{s}=1$ for the field. The effect of
$c_{s}<1$ is depicted in Fig. 2.

Figure 2:  Effect of a reduced speed of sound $c_{s}<1$ on the spatial
dependence of the field correlator $\Delta_{\Phi,\Phi}$, computed with quantum
field theory: $c_{s}=0.2$ (top left), $c_{s}=0.1$ (top right), $c_{s}=0.05$
(bottom left) and $c_{s}=0.01$ (bottom right). The plots should be compared to
Fig. 1 (which uses $c_{s}=1$). For $c_{s}=0.2$ the spatial dependence at the
end of inflation and at matter radiation-equality (already well-fitted by the
stochastic formalism for $c_{s}=1$) are indistinguishable and correlations at
sub-Hubble scales today decrease more rapidly than in the case with $c_{s}=1$.
By diminishing $c_{s}$, the spatial profiles at the three different times
considered (end of inflation, matter-radiation equality and today) are barely
distinguishable, in incrementally better agreement with the prediction from
stochastic theory (constant spatial dependence after the end of inflation, up
to the evolution in time of the overall amplitude). The cosmological
parameters and the number of inflationary e-foldings have been chosen as for
Fig. 1.

###  3.3 Free fields and Wick’s theorem

We conclude this Section by discussing the QFT prediction for the higher-order
correlators (computed on the vacuum state $|\Omega\rangle$ defined by the
property (2.16)). We focus on

$\Delta_{\Phi^{2},\Phi^{2}}(t,r)\equiv\Bigl{\langle}\hat{\Phi}^{2}(t,\vec{x}_{1})\hat{\Phi}^{2}(t,\vec{x}_{2})\Bigr{\rangle}\,,$
(3.42)

which is the first of the 4-point functions considered in subsection 2.3, but
with the full field instead of just its long-wavelength part. Using the
decomposition (2.10a–2.10b) and the identities

$\displaystyle\Bigl{\langle}\hat{b}(\vec{k}_{1})\hat{b}(\vec{k}^{\prime}_{1})\hat{b}^{{\dagger}}(\vec{k}_{2})\hat{b}^{{\dagger}}(\vec{k}^{\prime}_{2})\Bigr{\rangle}$
$\displaystyle=$
$\displaystyle\delta^{3}\left(\vec{k}_{1}-\vec{k}_{2}\right)\delta^{3}\left(\vec{k}^{\prime}_{1}-\vec{k}^{\prime}_{2}\right)+\delta^{3}\left(\vec{k}_{1}-\vec{k}^{\prime}_{2}\right)\delta^{3}\left(\vec{k}^{\prime}_{1}-\vec{k}_{2}\right)\,,$
$\displaystyle\Bigl{\langle}\hat{b}(\vec{k}_{1})\hat{b}^{{\dagger}}(\vec{k}^{\prime}_{1})\hat{b}(\vec{k}_{2})\hat{b}^{{\dagger}}(\vec{k}^{\prime}_{2})\Bigr{\rangle}$
$\displaystyle=$
$\displaystyle\delta^{3}\left(\vec{k}_{1}-\vec{k}^{\prime}_{1}\right)\delta^{3}\left(\vec{k}_{2}-\vec{k}^{\prime}_{2}\right)\,,$
(3.43)

which follow from eqs. (2.11) and (2.16), one gets

$\Delta_{\Phi^{2},\Phi^{2}}(t,r)=\left[\frac{1}{2\pi^{2}}\int_{k_{0}}^{\infty}{\rm
d}k~{}k^{2}|\varphi(t,k)|^{2}\right]^{2}+2\left[\frac{1}{2\pi^{2}}\int_{k_{0}}^{\infty}{\rm
d}k~{}k^{2}j_{0}(kr)|\varphi(t,k)|^{2}\right]^{2}\,,$ (3.44)

where, again, we regularized the infrared sector with the IR cutoff $k_{0}$.
Comparing to the first of eqs. (3.4), we immediately recognize that

$\Delta_{\Phi^{2},\Phi^{2}}(t,r)=\left[\Delta_{\Phi\Phi}(t,0)\right]^{2}+2\left[\Delta_{\Phi\Phi}(t,r)\right]^{2}\,.$
(3.45)

This result agrees with Wick’s theorem, i.e. the statement that higher-order
correlators of free fields can be expressed as the sum of all possible pair
contraction products:

$\displaystyle\Delta_{\Phi^{2},\Phi^{2}}(t,r)=\Bigl{\langle}\hat{\Phi}(t,\vec{x}_{1})\hat{\Phi}(t,\vec{x}_{1})\hat{\Phi}(t,\vec{x}_{2})\hat{\Phi}(t,\vec{x}_{2})\Bigr{\rangle}$
$\displaystyle=$
$\displaystyle\Bigl{\langle}\hat{\Phi}(t,\vec{x}_{1})\hat{\Phi}(t,\vec{x}_{1})\Bigr{\rangle}\Bigl{\langle}\hat{\Phi}(t,\vec{x}_{2})\hat{\Phi}(t,\vec{x}_{2})\Bigr{\rangle}\,$
(3.46) $\displaystyle+$ $\displaystyle
2\Bigl{\langle}\hat{\Phi}(t,\vec{x}_{1})\hat{\Phi}(t,\vec{x}_{2})\Bigr{\rangle}\Bigl{\langle}\hat{\Phi}(t,\vec{x}_{1})\hat{\Phi}(t,\vec{x}_{2})\Bigr{\rangle}\,$
$\displaystyle=$
$\displaystyle\left[\Delta_{\Phi,\Phi}(t,0)\right]^{2}+2\left[\Delta_{\Phi,\Phi}(t,r)\right]^{2}\,.$

This is a consequence of the one-loop approximation to QFT and therefore the
use of free quantum fields in the correlators. Wick’s theorem also applies
more generally in QFT than just eq. (3.46), where some $n$-point function is
given in terms of products of 2-point functions plus irreducible $n$-pt
functions. Canonical quantization via eqs. (2.10a), (2.10b), (2.11) applies to
free fields, which is the case for the action (2.1) as long as one considers a
fixed classical background metric. In the most recent cosmological times, the
backreaction of the scalar field becomes large (eventually dominating the
expansion as dark energy) and therefore one should treat the metric as a
quantum field to take into account the strong fluctuations induced by the
energy-momentum tensor of the scalar field. A simple way to approximate these
last stages of non-Gaussian evolution is to assume that the Friedmann equation
is still valid, but the Hubble rate is a local quantum operator
$\hat{H}(t,\vec{x})$ (therefore fluctuating):

$3M_{P}^{2}\hat{H}^{2}(t,\vec{x})=\rho_{C}(t)\hat{1\\!\\!1}+\hat{\rho}_{Q}(t,\vec{x})\,,$
(3.47)

where $M_{P}$ is the reduced Planck mass $M_{P}=\left(8\pi
G\right)^{-\frac{1}{2}}$, $\rho_{C}(t)$ is the classical272727$\rho_{C}$ is
dominated by dark matter, which may also be of quantum origin and therefore
with a fluctuating character [52, 53, 54, 55], which should be included in a