complete treatment. For simplicity, here we opt against the inclusion of such
effects. energy density contribution which dominates the expansion until the
backreaction of the scalar field is large, $\hat{1\\!\\!1}$ is the identity on
the Hilbert space of quantum states and $\hat{\rho}_{Q}(t,\vec{x})$ is the
energy density operator of the scalar field given by eq. (2.19). The quantum
version (3.47) of the Friedmann equation will be the starting point of Section
5.

Contrarily to the quantum field theory answer (3.45), even when fixing the
classical background metric (as FLRW) so that the spectator scalar field is
Gaussian, the result of stochastic formalism in eqs. (2.89)–(2.90) does not
respect Wick’s theorem. This is the main reason for studying an alternative
stochastic approach in the next Section 4.

##  4 Stochastic solution from Fokker-Planck equation

In this Section, we discuss an alternative solution for the time evolution of
correlators during inflation within stochastic theory. Instead of deriving and
solving a set of differential equations obeyed by the correlators as in eqs.
(2.31)–(2.33), we will work with probability distributions for the field,
extending the formalism developed by Starobinsky and Yokoyama in [56]. The
result for the 2-point field correlator
$\langle\hat{\phi}(t,\vec{x})\hat{\phi}(t,\vec{y})\rangle$ will be the same as
the one derived from the system of equations in subsection 2.2, but remarkably
the two methods give a different prediction for the dependence of the 4-pt
function $\langle\hat{\phi}^{2}(t,\vec{x})\hat{\phi}^{2}(t,\vec{y})\rangle$ on
the relative distance between the two points $\vec{x}$ and $\vec{y}$. Indeed,
the method discussed in this Section will provides a result obeying Wick’s
theorem while the outcome obtained from the method described in subsection 2.3
did not satisfy it.

For simplicity we focus on the evolution during de Sitter inflation, but the
treatment could be extended to post-inflationary phases taking the
inflationary results as initial conditions. Similarly to what we discussed in
the previous approach of subsections 2.2 and 2.3, in post-inflationary epochs
the contribution of stochastic sources is negligible with respect to the
initial conditions inherited from inflation.

If $H_{I}$ is the constant Hubble rate in de Sitter inflation, in the slow-
roll regime the long-wavelength part (coarse-grained) of the scalar field
obeys a Langevin equation with stochastic white noise, see eqs. (9)-(11) of
[56]. This observation is the basis for simplifying the quantum problem into a
classical stochastic process.

###  4.1 One-point probability distribution

Following [56], let us start with the evolution of the one-point probability
distribution function (PDF)
$\rho_{1}\left[\varphi(t,\vec{x})\right]=\rho_{1}(\varphi,t)$ expressing the
probability density function for the value of the coarse-grained field at
comoving position $\vec{x}$ at the time $t$. Due to space homogeneity this
cannot depend on the position $\vec{x}$, but only on time $t$ and the field
value $\varphi$. Coincident correlators and, more generally, expectation
values of a function $F_{1}$ of the field at a single spacetime point
$(t,\vec{x})$ can be computed as

$\langle F_{1}\left(\varphi(t,\vec{x})\right)\rangle=\int{\rm
d}\varphi~{}\rho_{1}(\varphi,t)F_{1}(\varphi)$ (4.1)

where we normalize $\rho_{1}(\varphi,t)$ to unity:

$\int{\rm d}\varphi~{}\rho_{1}(\varphi,t)=1\,,$

the integral being evaluated on the support of $\rho_{1}(\varphi,t)$.

Following [56], just like for Brownian motion, the Langevin equation for the
coarse-grained field implies a Fokker-Planck equation for the evolution of the
probability distribution (see eq. (12) of [56]):

$\frac{\partial}{\partial
t}\rho_{1}(\varphi,t)=\hat{\Gamma}_{\varphi}\rho_{1}(\varphi,t)\equiv\frac{1}{3H_{I}}\frac{\partial}{\partial\varphi}\left(V^{\prime}(\varphi)\rho_{1}(\varphi,t)\right)+\frac{H_{I}^{3}}{8\pi^{2}}\frac{\partial^{2}}{\partial\varphi^{2}}\rho_{1}(\varphi,t)$
(4.2)

where the second equality defines the action of the derivative operator
$\hat{\Gamma}_{\varphi}$ on $\rho_{1}(\varphi,t)$ and $V(\varphi)$ is the
potential of the field $\varphi$. The subscript “1” in $\rho_{1}(\varphi,t)$
just reminds that this is a one-point PDF. Eq. (4.2) has been used in [56] to
study dynamical mass generation in de Sitter from a self-interacting potential
$V(\varphi)$, e.g. quartic self-interaction. In the present work what is
relevant for us is just the case of (2.1), which corresponds to a free field
and therefore a non-interacting potential
$V(\varphi)=\frac{1}{2}M^{2}\varphi^{2}$. Then
$V^{\prime}(\varphi)=M^{2}\varphi$ where $M^{2}=m^{2}+\xi R=m^{2}+12\xi
H_{I}^{2}$ is the effective squared mass, which is constant in de Sitter. Thus
we can reduce (4.2) to

$\frac{1}{H_{I}}\frac{\partial}{\partial
t}\rho_{1}(\varphi,t)=\beta\frac{\partial}{\partial\varphi}\left(\varphi\rho_{1}(\varphi,t)\right)+\frac{H_{I}^{2}}{8\pi^{2}}\frac{\partial^{2}}{\partial\varphi^{2}}\rho_{1}(\varphi,t)\,,$
(4.3)

where we introduced the dimensionless parameter $\beta$ as

$\beta\equiv\frac{1}{3}\frac{M^{2}}{H_{I}^{2}}=\frac{1}{3}\frac{m^{2}}{H_{I}^{2}}+4\xi\,.$
(4.4)

At the beginning of inflation $t=t_{\rm in}$ we take the initial PDF to be

$\rho_{1}(\varphi,t_{\rm in})=\delta(\varphi)\,,$ (4.5)

so that the field has no condensate and all higher order self-correlations are
zero, in agreement with the initial conditions that we used in subsections 2.2
and 2.3. The solution of the Fokker-Planck equation in the non-interacting
case with Dirac delta initial condition is known to be a Gaussian
distribution. It is straightforward to check the validity of the following
solution (properly normalized to one):

$\rho_{1}(\varphi,t)=\frac{1}{\sigma(t)\sqrt{2\pi}}\exp{\left[-\frac{\varphi^{2}}{2\sigma^{2}(t)}\right]}\,,\qquad\sigma^{2}(t)=\frac{H_{I}^{2}}{8\pi^{2}}\frac{1-e^{-2\beta
H_{I}(t-t_{\rm in})}}{\beta}\,.$ (4.6)

The property $\sigma^{2}(t_{\rm in})=0$ ensures compliance with the initial
condition $\rho_{1}(\varphi,t_{\rm in})=\delta(\varphi)$. We are ultimately
interested in $\beta\approx-4|\xi|<0$ corresponding to a negative non-minimal
coupling dominating over the very light mass, but the solution (4.6) is valid
for any sign of $\beta$. It is also interesting to observe that the free
massless minimally coupled case $\beta=0$ can be obtained by taking the limit
$\beta\to 0$ of (4.6), which then gives the random walk behavior with variance
growing linearly in time $\sigma^{2}(t)=H_{I}^{3}(t-t_{\rm in})/(4\pi^{2})$.
The presence of a non-zero mass (i.e. $\beta$) is responsible for deviations
from this pure random walk evolution. When $\beta>0$ the variance
$\sigma^{2}(t)$ asymptotically approaches the finite value
$H_{I}^{2}/(8\pi^{2}\beta)$ and the existence of a stationary solution is a
consequence of the stability of the potential $V(\varphi)=(3/2)\beta
H_{I}^{2}\varphi^{2}$, which has $V^{\prime\prime}(\varphi)=M^{2}=3\beta
H_{I}^{2}$ positive for $\beta>0$ (upward concavity of the parabola). However,
in the most relevant case for us $\beta<0$ and the variance grows
exponentially, which is in agreement with what one would expect from an
unstable potential (downward concave parabola).

We can immediately read the coincident 2-point function
$\Delta_{\phi\phi}(t)=\langle\hat{\phi}^{2}(t,\vec{x})\rangle$ from (4.6)
because it is simply given by the variance $\sigma^{2}(t)$. In terms of the
number of e-foldings $N=H_{I}(t-t_{\rm in})$ measured from the beginning of
inflation one has

$\Delta_{\phi\phi}(N)=\frac{H_{I}^{2}}{8\pi^{2}}\frac{1-e^{-2\beta
N}}{\beta}\,,$ (4.7)

which, when $\beta\simeq-4|\xi|$, shows the same exponential growth with $N$
predicted before in [31] and written in eq. (2.52).

###  4.2 Two-point probability distribution

We now turn our attention to non-coincident correlators. Their study requires
knowledge of the joint probability distribution at two different spacetime
points. Of course, for our purposes it is sufficient to consider fields at the
same time $t$ but different positions $\vec{x}_{1}$ and $\vec{x}_{2}$. Our
starting point is eq. (73) in [56], which is the Fokker-Planck equation for
the joint two-point PDF $\rho_{2}(\varphi_{1},t;\varphi_{2},t)$:

$\frac{\partial}{\partial
t}\rho_{2}(\varphi_{1},t;\varphi_{2},t)=\hat{\Gamma}_{\varphi_{1}}\rho_{2}+\hat{\Gamma}_{\varphi_{2}}\rho_{2}+j_{0}\left(\mu
a(t)H_{I}r\right)\frac{H_{I}^{3}}{4\pi^{2}}\frac{\partial^{2}}{\partial\varphi_{1}\partial\varphi_{2}}\rho_{2}\,.$
(4.8)

The subscript “2” in $\rho_{2}$ signals that this is a two-point joint PDF
while the arguments $\varphi_{1},t$ and $\varphi_{2},t$ mean that we are
considering the probability density for the field at $\vec{x}_{1}$ and time
$t$ to have value $\varphi_{1}$ and for the field at $\vec{x}_{2}$ and same
time $t$ to have value $\varphi_{2}$. The distance $r$ appearing in the
argument of the $j_{0}$ spherical Bessel function is the comoving distance
between $\vec{x}_{1}$ and $\vec{x}_{2}$ given by
$r=\|\vec{x}_{2}-\vec{x}_{1}\|$. The derivative operators
$\hat{\Gamma}_{\varphi_{1}}$ and $\hat{\Gamma}_{\varphi_{2}}$ have the same
form as $\hat{\Gamma}_{\varphi}$ defined in (4.2), but now with respect to
${\varphi_{1}}$ and ${\varphi_{2}}$ respectively. Just like for (4.1),
expectation values of a function $F_{2}$ of fields at the two spacetime points
$(t,\vec{x}_{1})$ and $(t,\vec{x}_{2})$ can be computed using
$\rho_{2}(\varphi_{1},t;\varphi_{2},t)$ as

$\langle
F_{2}\left(\varphi(t,\vec{x}_{1});\varphi(t,\vec{x}_{2})\right)\rangle=\int{\rm
d}\varphi_{1}\int{\rm
d}\varphi_{2}~{}\rho_{2}(\varphi_{1},t;\varphi_{2},t)F_{2}(\varphi_{1};\varphi_{2})\,,$
(4.9)

with the normalization condition

$\int{\rm d}\varphi_{1}\int{\rm
d}\varphi_{2}~{}\rho_{2}(\varphi_{1},t;\varphi_{2},t)=1\,.$

We now explain how to obtain $\rho_{2}(\varphi_{1},t;\varphi_{2},t)$ solving
eq. (4.8). The reader who is mainly interested in the results for correlators
might want to skip this part and jump directly to the summary in subsection
4.2.2.

####  4.2.1 Solution for the joint two-point PDF

Following [56], since we are dealing with the coarse-grained field (long-
wavelength part), the spherical Bessel function $j_{0}\left(\mu
a(t)H_{I}r\right)$ should be approximated by a Heaviside step function:

$j_{0}\left(\mu a(t)H_{I}r\right)\simeq\theta\left(1-\mu a(t)H_{I}r\right)\,.$
(4.10)

This corresponds exactly to the approximation that we use for noise sources in
eq. (2.48). Given the two comoving positions $\vec{x}_{1}$ and $\vec{x}_{2}$
separated by a comoving distance $r$, we can distinguish two periods in the
time evolution of $\rho_{2}(\varphi_{1},t;\varphi_{2},t)$ depending on whether
the Heaviside step function (4.10) evaluates to unity or zero. The first phase
starts from the beginning of inflation $t_{\rm in}$ and lasts until the time
$t_{r}\equiv t_{\rm in}-\ln{(\mu a_{\rm in}H_{I}r)}/H_{I}$. In this time
interval the Heaviside function evaluates to unity. The second phase is
$t>t_{r}$, and the Heaviside function becomes zero. In the first phase $t_{\rm
in}<t<t_{r}$ one can check that the 2-point Fokker-Planck equation (4.8)
admits a solution of the form

$\rho_{2}(\varphi_{1},t;\varphi_{2},t)=\rho_{1}(\varphi_{1},t)\delta(\varphi_{1}-\varphi_{2})\,,$
(4.11)

where $\rho_{1}(\varphi_{1},t)$ is the one-point PDF studied in the previous
subsection 4.1 evaluated at time $t$ for the field value $\varphi_{1}$. The
existence of such a solution was pointed out in eq. (75) of [56] in the
stationary late-times regime, but it is actually valid in full generality with
time dependence, so that is applicable to study the time evolution of
correlators. As a check of the initial conditions that we are using, at the
beginning of inflation $t_{\rm in}$ the solution (4.11) gives
$\rho_{2}(\varphi_{1},t_{\rm in};\varphi_{2},t_{\rm
in})=\delta(\varphi_{1})\delta(\varphi_{2})$, because of the initial condition
used for the one-point PDF $\rho_{1}(\varphi_{1},t_{\rm
in})=\delta(\varphi_{1})$. This initial value of $\rho_{2}(\varphi_{1},t_{\rm
in};\varphi_{2},t_{\rm in})$ implies $\langle\varphi(t_{\rm
in},\vec{x}_{1})\varphi(t_{\rm in},\vec{x}_{2})\rangle=0$ which is equivalent
to the zero initial conditions for correlators used in Section 2.

It is important to notice that the first phase $t_{\rm in}<t<t_{r}$ really
exists only as long as $t_{r}>t_{\rm in}$, i.e. when $r<r_{0}\equiv 1/(\mu
a_{\rm in}H_{I})$. Indeed when $r>r_{0}$ the Heaviside function evaluates to
zero already at the beginning of inflation and it keeps being null at later
times. Therefore when $r>r_{0}$ we are immediately in the “second” phase,
already from the beginning of inflation $t_{\rm in}$ and we will get back to
this case later.

For now let us assume $r<r_{0}$ and discuss the solution for the two-point
joint PDF in the second phase $t>t_{r}$, when the last term in (4.8) drops to
zero. The initial conditions are those inherited from the first phase at time
$t_{r}$, i.e.

$\rho_{2}(\varphi_{1},t_{r};\varphi_{2},t_{r})=\rho_{1}(\varphi_{1},t_{r})\delta(\varphi_{1}-\varphi_{2})\,.$
(4.12)

The relevant information for the solution at $t>t_{r}$ can be extracted from
eq. (79) of [56] paying attention to the fact that we are interested in the
full time evolution and not in the stationary regime. The solution at
$t>t_{r}$ can be expressed as

$\rho_{2}(\varphi_{1},t;\varphi_{2},t)=\int{\rm
d}\varphi_{r}~{}\rho_{1}(\varphi_{r},t_{r})\Pi(\varphi_{1},t|\varphi_{r},t_{r})\Pi(\varphi_{2},t|\varphi_{r},t_{r})\,,$
(4.13)

where $\Pi(\varphi_{1},t|\varphi_{r},t_{r})$ is the conditional probability
density that the field at position $\vec{x}_{1}$ and time $t$ has value
$\varphi_{1}$ given that at the same position $\vec{x}_{1}$ the field at time
$t_{r}$ had value $\varphi_{r}$. A completely analogous meaning is attributed
to $\Pi(\varphi_{2},t|\varphi_{r},t_{r})$ but now the fields refer to position
$\vec{x}_{2}$. The integral in (4.13) goes over all possible values of
$\varphi_{r}$. The conditional probability
$\Pi(\varphi_{1},t|\varphi_{r},t_{r})$ satisfies a Fokker-Planck equation (see
eq. (26) of [56]):

$\frac{\partial}{\partial
t}\Pi(\varphi_{1},t|\varphi_{r},t_{r})=\hat{\Gamma}_{\varphi_{1}}\Pi(\varphi_{1},t|\varphi_{r},t_{r})\,,$
(4.14)

with initial conditions

$\Pi(\varphi_{1},t_{r}|\varphi_{r},t_{r})=\delta(\varphi_{1}-\varphi_{r})\,.$
(4.15)

A similar equation holds with respect to the second argument
$\varphi_{r},t_{r}$ (see eq. (27) of [56]), but it will not be needed. In the
same way $\Pi(\varphi_{2},t|\varphi_{r},t_{r})$ follows the equation

$\frac{\partial}{\partial
t}\Pi(\varphi_{2},t|\varphi_{r},t_{r})=\hat{\Gamma}_{\varphi_{2}}\Pi(\varphi_{2},t|\varphi_{r},t_{r})\,,$
(4.16)

with initial conditions

$\Pi(\varphi_{2},t_{r}|\varphi_{r},t_{r})=\delta(\varphi_{2}-\varphi_{r})\,.$
(4.17)

It is straightforward to check that if the conditional probabilities satisfy
eqs. (4.14) and (4.16) with their initial conditions (4.15) and (4.17)
respectively, then (4.13) is a solution of (4.8) (with the spherical Bessel
function $j_{0}$ approximated by zero as it is appropriate for $t>t_{r}$) with
initial conditions (4.12). Thus the problem completely reduces to solving eqs.
(4.14) and (4.16).

The solution for conditional probabilities can be expressed in terms of the
Mehler kernel282828Defining
$\Psi=\Pi\exp{\left[-\frac{\beta}{2}\left(H_{I}t-4\pi^{2}\frac{\varphi_{1}^{2}}{H_{I}^{2}}\right)\right]}$
and changing variables from $t$ and $\varphi_{1}$ to
$\tau\equiv\frac{1}{2}|\beta|H_{I}t$ and $\chi\equiv
2\pi\sqrt{|\beta|}\frac{\varphi_{1}}{H_{I}}$, eq. (4.14) is equivalent to
$\frac{\partial}{\partial\tau}\Psi=\left(\frac{\partial^{2}}{\partial\chi^{2}}-\chi^{2}\right)\Psi$.
The Green’s function of this equation is known as the Mehler kernel, see also
Section V.A of [60]. Its importance in physics stems from being the propagator
of the quantum harmonic oscillator (upon changing $\tau$ into the imaginary
$i\tau$), see eq. (2.5.18) of Sakurai’s textbook [61] for its expression in
the quantum mechanics context., which can be obtained by making a Gaussian
ansatz, checking its validity and solving for the mean and variance as
functions of time, thus reducing the problem to ordinary differential
equations. The result is that eq. (4.14) with initial condition (4.15) is
solved by a Gaussian distribution with mean and variance evolving with time as
follows:

$\displaystyle\Pi(\varphi_{1},t|\varphi_{r},t_{r})=\frac{1}{\sigma_{r}(t)\sqrt{2\pi}}\exp{\left[-\frac{\left(\varphi_{1}-\mu_{r}(t)\right)^{2}}{2\sigma_{r}^{2}(t)}\right]}\,,$
$\displaystyle\mu_{r}(t)=\varphi_{r}~{}e^{-\beta
H_{I}(t-t_{r})}\,,\qquad\,\sigma_{r}^{2}(t)=\frac{H_{I}^{2}}{8\pi^{2}}\frac{1-e^{-2\beta
H_{I}(t-t_{r})}}{\beta}\,.$ (4.18)

The time evolution of the variance $\sigma_{r}^{2}(t)$ is perfectly analogous
to the solution for the one-point PDF (4.6). The novel element is the drift of
the distribution, namely the translation of its mean value $\mu_{r}(t)$. This
effect was not present in (4.6) because of the zero initial conditions (4.5),
which kept the mean to zero also at later times. On the contrary, a generic
$\varphi_{r}$ appears in the initial conditions (4.15). In the same way, the
solution of (4.16) with initial conditions (4.17) is

$\Pi(\varphi_{2},t|\varphi_{r},t_{r})=\frac{1}{\sigma_{r}(t)\sqrt{2\pi}}\exp{\left[-\frac{\left(\varphi_{2}-\mu_{r}(t)\right)^{2}}{2\sigma_{r}^{2}(t)}\right]}\,,$
(4.19)

with the same $\mu_{r}(t)$ and $\sigma_{r}(t)$ as in (4.18).

The only other result needed for evaluating the joint two-point PDF according
to eq. (4.13) is the one-point PDF $\rho_{1}(\varphi_{r},t_{r})$ which is
readily obtained by substituting its arguments in the general formula (4.6):

$\rho_{1}(\varphi_{r},t_{r})=\frac{1}{\sigma(t_{r})\sqrt{2\pi}}\exp{\left[-\frac{\varphi_{r}^{2}}{2\sigma^{2}(t_{r})}\right]}\,,\qquad\sigma^{2}(t_{r})=\frac{H_{I}^{2}}{8\pi^{2}}\frac{1-e^{-2\beta
H_{I}(t_{r}-t_{\rm in})}}{\beta}\,.$ (4.20)

One can then perform the integration (4.13) paying attention to the dependence
on $\varphi_{r}$ in $\mu_{r}(t)$ given by (4.18). This leads after a few steps
to the following two-point joint PDF for $t>t_{r}$:

$\displaystyle\rho_{2}(\varphi_{1},t;\varphi_{2},t)=\frac{1}{2\pi\sqrt{{\rm
det}(\Sigma(t,r))}}\exp\left[-\frac{1}{2}(\varphi_{1},\varphi_{2})\Sigma^{-1}(t,r)\begin{pmatrix}\varphi_{1}\\\
\varphi_{2}\\\ \end{pmatrix}\,\right]\,,$
$\displaystyle\Sigma(t,r)=\begin{pmatrix}\sigma^{2}(t)&\sigma^{2}(t)-\sigma_{r}^{2}(t)\\\
\sigma^{2}(t)-\sigma_{r}^{2}(t)&\sigma^{2}(t)\end{pmatrix}\,,$ (4.21)

which is a two-dimensional Gaussian distribution with zero mean for both
$\varphi_{1}$ and $\varphi_{2}$ and with covariance matrix $\Sigma(t,r)$
specified above. Notice that the diagonal terms of $\Sigma(t,r)$ are both
equal to $\sigma^{2}(t)$ consistently with the one-point PDF in eq. (4.6) that
we are supposed to obtain when marginalizing over one of the two variables. It
is worth reminding that, for a given comoving spatial separation $r$ between
the two points, the result (4.2.1) is valid when $t>t_{r}>t_{\rm in}$, which
using $t_{r}=t_{\rm in}-\frac{1}{H_{I}}\ln\left(\frac{r}{r_{0}}\right)$
translates into $N>\ln\left(\frac{r_{0}}{r}\right)>0$. The last inequality
makes sense only when $r<r_{0}$.

We also remind that when $t_{\rm in}<t<t_{r}$, which translates into
$0<N<\ln\left(\frac{r_{0}}{r}\right)$, then the solution for
$\rho_{2}(\varphi_{1},t;\varphi_{2},t)$ has been written in eq. (4.11). Again
this requires $r<r_{0}$.

Now we study $\rho_{2}(\varphi_{1},t;\varphi_{2},t)$ in the case $r>r_{0}$.
For such distances the Bessel function in (4.8) (and approximated by a
Heaviside step function) is always zero from the very beginning of inflation.
Then the solution is given by

$\rho_{2}(\varphi_{1},t;\varphi_{2},t)=\int{\rm d}\varphi_{\rm
in}~{}\rho_{1}(\varphi_{\rm in},t_{\rm in})\Pi(\varphi_{1},t|\varphi_{\rm
in},t_{\rm in})\Pi(\varphi_{2},t|\varphi_{\rm in},t_{\rm in})\,,$ (4.22)

with the same meaning of terms as in (4.13), except for the replacement of
$\varphi_{r},t_{r}$ by $\varphi_{\rm in},t_{\rm in}$ because when $r>r_{0}$
there are no longer two distinct phases in the evolution of
$\rho_{2}(\varphi_{1},t;\varphi_{2},t)$ and the $t_{r}$ makes no sense.
Therefore one has to integrate over values $\varphi_{\rm in}$ of the field at
time $t_{\rm in}$, instead of $\varphi_{r}$ at $t_{r}$. The conditional
probabilities are now given by

$\displaystyle\Pi(\varphi_{1},t|\varphi_{\rm in},t_{\rm
in})=\frac{1}{\sigma(t)\sqrt{2\pi}}\exp{\left[-\frac{\left(\varphi_{1}-\mu_{\rm
in}(t)\right)^{2}}{2\sigma^{2}(t)}\right]}\,,$ $\displaystyle\mu_{\rm
in}(t)=\varphi_{\rm in}~{}e^{-\beta H_{I}(t-t_{\rm
in})}\,,\qquad\,\sigma^{2}(t)=\frac{H_{I}^{2}}{8\pi^{2}}\frac{1-e^{-2\beta
H_{I}(t-t_{\rm in})}}{\beta}\,,$ (4.23)

and similarly

$\Pi(\varphi_{2},t|\varphi_{\rm in},t_{\rm
in})=\frac{1}{\sigma(t)\sqrt{2\pi}}\exp{\left[-\frac{\left(\varphi_{2}-\mu_{\rm
in}(t)\right)^{2}}{2\sigma^{2}(t)}\right]}\,.$ (4.24)

Notice that $\sigma(t)$ in eqs. (4.23–4.24) is exactly the same defined for
the one-point PDF (4.6). This is again a consequence of considering $r>r_{0}$,
in which case one has initial conditions set at $t_{\rm in}$, instead of the
other time $t_{r}$ introduced when $r<r_{0}$. Finally, due to the initial
conditions (4.5), one has

$\rho_{1}(\varphi_{\rm in},t_{\rm in})=\delta(\varphi_{\rm in})\,.$ (4.25)

The integration in (4.22) immediately gives the solution for $r>r_{0}$:

$\displaystyle\rho_{2}(\varphi_{1},t;\varphi_{2},t)=\frac{1}{2\pi\sigma^{2}(t)}\exp{\left[-\frac{\varphi_{1}^{2}+\varphi_{2}^{2}}{2\sigma^{2}(t)}\right]}=\rho_{1}(\varphi_{1},t)\rho_{1}(\varphi_{2},t)\,,$
(4.26)

which is a Gaussian with diagonal covariance matrix, from which we infer the
absence of correlations between $\varphi_{1}$ and $\varphi_{2}$ when
$r>r_{0}$.

####  4.2.2 Summary on joint PDF and non-coincident correlators

Before moving to the calculation of non-coincident 2-point and 4-point
functions, let us summarize the results on the joint two-point probability
distribution $\rho_{2}(\varphi_{1},t;\varphi_{2},t)$ resulting from the
solution of the Fokker-Planck equation (4.8). We introduced the scale

$r_{0}\equiv\frac{1}{\mu a_{\rm in}H_{I}}\,,$ (4.27)

which sets the separation between long and short modes of the stochastic
theory at the beginning of inflation. In (4.4) we also defined for convenience
a dimensionless parameter related to the effective mass of the scalar field as
$\beta\equiv\frac{1}{3}\frac{M^{2}}{H_{I}^{2}}=\frac{1}{3}\frac{m^{2}}{H_{I}^{2}}+4\xi$.
We computed $\rho_{2}(\varphi_{1},t;\varphi_{2},t)$ depending on the number of
e-foldings since the beginning of inflation $N=H_{I}(t-t_{\rm in})$ and on the
separation $r=\|\vec{x}_{2}-\vec{x}_{1}\|$ between the comoving positions
$\vec{x}_{1}$ and $\vec{x}_{2}$. The results from (4.11), (4.2.1) and (4.26)
can be effectively summarized as follows:

$\boxed{\rho_{2}(\varphi_{1},N;\varphi_{2},N)=\begin{cases}\frac{1}{\sigma(N)\sqrt{2\pi}}\exp{\left[-\frac{\varphi_{1}^{2}}{2\sigma^{2}(N)}\right]}\delta(\varphi_{1}-\varphi_{2})&\text{if
$r<r_{0}~{}e^{-N}$}\\\ \frac{1}{2\pi\sqrt{{\rm
det}(\Sigma(N,r))}}\exp\left[-\frac{1}{2}(\varphi_{1},\varphi_{2})\Sigma^{-1}(N,r)\begin{pmatrix}\varphi_{1}\\\
\varphi_{2}\\\ \end{pmatrix}\,\right]&\text{if $r_{0}~{}e^{-N}<r<r_{0}$}\\\
\frac{1}{2\pi\sigma^{2}(N)}\exp{\left[-\frac{\varphi_{1}^{2}+\varphi_{2}^{2}}{2\sigma^{2}(N)}\right]}&\text{if
$r>r_{0}$}\,,\end{cases}}$ (4.28)

where $\sigma^{2}(N)$ is the variance of the one-point PDF defined in (4.6),
which we rewrite in terms of $N$ as

$\sigma^{2}(N)=\frac{H_{I}^{2}}{8\pi^{2}}\frac{1-e^{-2\beta N}}{\beta}\,,$
(4.29)

and $\Sigma(N)$ is the matrix defined in (4.2.1), which, upon using
$t_{r}=t_{\rm in}-\frac{1}{H_{I}}\ln\left(\frac{r}{r_{0}}\right)$, can be
equivalently expressed as

$\Sigma(N,r)=\frac{H_{I}^{2}}{8\pi^{2}\beta}\begin{pmatrix}1-e^{-2\beta
N}&e^{-2\beta N}\left[\left(\frac{r}{r_{0}}\right)^{-2\beta}-1\right]\\\
e^{-2\beta
N}\left[\left(\frac{r}{r_{0}}\right)^{-2\beta}-1\right]&1-e^{-2\beta
N}\end{pmatrix}\,.$ (4.30)

One can immediately recognize from eq. (4.28) that
$\rho_{2}(\varphi_{1},N;\varphi_{2},N)$ exhibits perfect correlation between
$\varphi_{1}$ and $\varphi_{2}$ for separations $r<r_{0}~{}e^{-N}$ (due to the
Dirac delta function $\delta(\varphi_{1}-\varphi_{2})$) and it is completely
uncorrelated for separations $r>r_{0}$. For intermediate separations between
$\vec{x}_{1}$ and $\vec{x}_{2}$, correlations go from their maximum value to
the minimum one, interpolating between the two extremes according to the off-
diagonal element(s) of $\Sigma(N,r)$ in eq. (4.30).

Quantitatively, the computation of the 2-point function
$\langle\varphi_{1}\varphi_{2}\rangle$ follows from (4.9) giving, as a
function of the separation $r$ and the number of e-foldings $N$,

$\langle\varphi_{1}\varphi_{2}\rangle(N,r)=\begin{cases}\frac{H_{I}^{2}}{8\pi^{2}}\frac{1-e^{-2\beta
N}}{\beta}&\text{if $r<r_{0}~{}e^{-N}$}\\\
\frac{H_{I}^{2}}{8\pi^{2}}e^{-2\beta
N}\frac{\left(\frac{r}{r_{0}}\right)^{-2\beta}-1}{\beta}&\text{if
$r_{0}~{}e^{-N}<r<r_{0}$}\\\ 0&\text{if $r>r_{0}$}\,.\end{cases}$ (4.31)

The result can also be expressed in terms of the correlation coefficient of
the joint distribution (4.28), which we will denote by $C(N,r)$. For a
distribution like $\rho_{2}(\varphi_{1},N;\varphi_{2},N)$ having
$\langle\varphi_{1}\rangle=\langle\varphi_{2}\rangle=0$, the correlation
coefficient is defined as customary as

$C(N,r)\equiv\frac{\langle\varphi_{1}\varphi_{2}\rangle(N,r)}{\sqrt{\langle\varphi_{1}\rangle^{2}(N)}\sqrt{\langle\varphi_{2}\rangle^{2}(N)}}=\frac{\langle\varphi_{1}\varphi_{2}\rangle(N,r)}{\sigma^{2}(N)}\,,$
(4.32)

that is by rescaling the 2-point function
$\langle\varphi_{1}\varphi_{2}\rangle$ with respect to its value at
coincidence. From (4.31) it follows that

$C(N,r)=\begin{cases}1&\text{if $r<r_{0}~{}e^{-N}$}\\\
\frac{\left(\frac{r}{r_{0}}\right)^{-2\beta}-1}{e^{2\beta N}-1}&\text{if
$r_{0}~{}e^{-N}<r<r_{0}$}\\\ 0&\text{if $r>r_{0}$}\,,\end{cases}$ (4.33)

which interpolates continuously from $1$ to $0$ with a power law. For
$\beta\simeq-4|\xi|$ and $N$ large enough so that $e^{-8|\xi|N}\ll 1$, this
simplifies into

$\boxed{C(N,r)\simeq\begin{cases}1&\text{if $r<r_{0}~{}e^{-N}$}\\\
1-\left(\frac{r}{r_{0}}\right)^{8|\xi|}&\text{if $r_{0}~{}e^{-N}<r<r_{0}$}\\\
0&\text{if $r>r_{0}$}\,.\end{cases}}$ (4.34)

It agrees with the result found with another method in eq. (2.54), thus
providing further confirmation for its validity. The function $C(N,r)$ is
plotted against $r/r_{0}$ in Fig. 3. The precise choice of $N$ is irrelevant
as long as it satisfies $e^{-8|\xi|N}\ll 1$ . We use $\xi=-0.06$ for the
negative non-minimal coupling.

The correlation coefficient $C(N,r)$ determines the shape of the contour
levels of $\rho_{2}(\varphi_{1},N;\varphi_{2},N)$. Indeed it is
straightforward to show from (4.28) that the $1\sigma$ contour of
$\rho_{2}(\varphi_{1},N;\varphi_{2},N)$ can be written as

$\frac{\left(\frac{\varphi_{1}}{\sigma(N)}\right)^{2}-2~{}C(N,r)\frac{\varphi_{1}}{\sigma(N)}\frac{\varphi_{2}}{\sigma(N)}+\left(\frac{\varphi_{2}}{\sigma(N)}\right)^{2}}{1-C^{2}(N,r)}=1\,,$
(4.35)

where $\sigma(N)$ is the standard deviation of the one-point PDF evolving with
time as in eq. (4.29). The $1\sigma$ regions corresponding to (4.35) are
plotted in Fig. 4 at late times for different separations $r$.

Figure 3: The function $C(N,r)$ in (4.34), equal to the non-coincident 2-point
function $\Delta_{\phi\phi}(N,r)$ rescaled by its value at coincidence,
plotted as a function of the separation $r$ between the points (rescaled by
$r_{0}=1/(\mu a_{\rm in}H_{I})$) at late times in inflation, $e^{8|\xi|N}\gg
1$. We used a negative non-minimal coupling $\xi=-0.06$. The function stays
constant at zero for $r/r_{0}>1$. The other scale $r/r_{0}=e^{-N}$ appearing
in eq. (4.34) and setting the upper limit for the coincident regime is too
small to be visible in the plot, $e^{-N}\ll 1$.  Figure 4: The $1\sigma$
contour plot of the joint distribution $\rho_{2}(\varphi_{1},N;\varphi_{2},N)$
in (4.28) at a given large N (with $e^{8|\xi|N}\gg 1$), for different values
of the comoving separation $r=\|\vec{x}_{2}-\vec{x}_{1}\|$ between the points
$\vec{x}_{1}$ and $\vec{x}_{2}$, which the field values $\varphi_{1}$ and
$\varphi_{2}$ refer to. The variables $\varphi_{1}$ and $\varphi_{2}$ have
been rescaled by their (equal) standard deviation $\sigma(N)$ written in
(4.29). The $1\sigma$ contour is an ellipse and it goes from being degenerate
as the red segment for $r<r_{0}e^{-N}$ to becoming a circle (with the blue
contour in the figure) when $r>r_{0}$. At intermediate distances,
$r_{0}e^{-N}<r<r_{0}$, the shape of the ellipse is determined by the
corresponding correlation coefficient $C(N,r)$ according to eq. (4.35).

Let us now study the correlator $\langle\varphi_{1}^{2}\varphi_{2}^{2}\rangle$
(which is the 4-point function $\Delta_{\phi^{2},\phi^{2}}$ in the notation of
Section 2) predicted by the Fokker-Planck solution (4.28). Again the
evaluation is done by using (4.9) with the PDF in (4.28) to get from direct
integration or by applying Isserlis’ theorem [62] (valid for Gaussian
distributions)

$\langle\varphi_{1}^{2}\varphi_{2}^{2}\rangle(N,r)=\begin{cases}3\left(\frac{H_{I}^{2}}{8\pi^{2}}\frac{1-e^{-2\beta
N}}{\beta}\right)^{2}&\text{if $r<r_{0}~{}e^{-N}$}\\\
\left(\frac{H_{I}^{2}}{8\pi^{2}}\frac{1-e^{-2\beta
N}}{\beta}\right)^{2}+2\left[\frac{H_{I}^{2}}{8\pi^{2}}e^{-2\beta
N}\frac{\left(\frac{r}{r_{0}}\right)^{-2\beta}-1}{\beta}\right]^{2}&\text{if
$r_{0}~{}e^{-N}<r<r_{0}$}\\\
\left(\frac{H_{I}^{2}}{8\pi^{2}}\frac{1-e^{-2\beta
N}}{\beta}\right)^{2}&\text{if $r>r_{0}$}\,,\end{cases}$ (4.36)

or, dividing by
$\sigma^{4}(N)=\left(\frac{H_{I}^{2}}{8\pi^{2}}\frac{1-e^{-2\beta
N}}{\beta}\right)^{2}$,

$\boxed{\frac{\langle\varphi_{1}^{2}\varphi_{2}^{2}\rangle(N,r)}{\sigma^{4}(N)}=\begin{cases}3&\text{if
$r<r_{0}~{}e^{-N}$}\\\ 1+2~{}C^{2}(N,r)&\text{if $r_{0}~{}e^{-N}<r<r_{0}$}\\\
1&\text{if $r>r_{0}$}\,.\end{cases}}$ (4.37)

As one can check this result obeys Wick’s theorem. Indeed this is just a
consequence of the Gaussianity of $\rho_{2}(\varphi_{1},N;\varphi_{2},N)$ in
(4.28). This is in agreement with the QFT calculation in Section 3 for super-
Hubble scales. For this reason we think that it correctly captures the
classical stochastic limit. On the contrary the result (2.90) produced by
solving the coupled system of equations (2.78)–(2.83) violates Wick’s theorem
and is in disagreement with the QFT treatment. We do not understand yet the
origin of the disagreement, considering that in the slow-roll regime, where
the contribution of correlators involving the canonical momentum
$\hat{\pi}(t,\vec{x})$ is subdominant, the two approaches should produce the
same results based on the same approximation used for the source
$n_{\phi^{2},\phi^{2}}$ and for the last term in the Fokker-Planck equation
(4.8), where $j_{0}(z)\rightarrow\theta(1-z)$. A possibility is that the
approach of Section 2 does not fully catch the interaction between long and
short modes for higher-order correlators, but only for 2-point correlators, so
that the sources appearing on the right-hand sides of eqs. (2.78)–(2.83) would
need to be modified with respect to those used, given in eq. (A.1).

Figure 5: Comparison between the spatial profile predicted by eq. (4.37) (blue
solid curve) and by eq. (2.90) (red dashed curve). The blue solid curve,
computed with the approach to stochastic formalism based on Fokker-Planck
equations discussed in this Section 4, obeys Wick’s theorem. On the contrary,
the red dashed curve, computed with the other approach to stochastic formalism
adopted in Section 2, based on linear systems of equations for the time
evolution of correlators, does not respect Wick’s theorem. We used $\xi=-0.06$
for the non-minimal coupling, as in Fig. 3.

Due to the Gaussian nature of the probability distribution (4.28), Wick’s
theorem (or direct integration via (4.9)) can be used to compute higher order
correlators, like $\langle\varphi_{1}^{n}\varphi_{2}^{n}\rangle$ with $n$ a
natural number, in terms of the correlation function $C(N,r)$ given in (4.33)
and the one-point variance $\sigma^{2}(N)$ in (4.29). For simplicity we will
omit the space ($r$) and time ($N$) dependence of $C(N,r)$ and $\sigma(N)$ in
the following expressions. By inspection of the combinatorics, it can be shown
that

$\frac{\langle\varphi_{1}^{n}\varphi_{2}^{n}\rangle}{\sigma^{2n}}=\sum_{k=0\atop
n-k={\rm even}}^{n}k!\left[{n\choose
k}(n-k-1)!!\right]^{2}C^{k}=\sum_{k=0\atop n-k={\rm
even}}^{n}\left[\frac{\Gamma(n+1)}{\Gamma\left(\frac{n-k+2}{2}\right)}\right]^{2}\frac{C^{k}}{2^{n-k}\Gamma(k+1)}\,,$
(4.38)

where the sum goes only on $k$ with the same parity as $n$, as specified above
by the condition that $n-k$ must be an even number. The last equality uses the
double factorial identity
$(2m-1)!!=\frac{(2m)!}{2^{m}m!}=\frac{\Gamma(2m+1)}{2^{m}\Gamma(m+1)}$,
applied to the natural number $m=(n-k)/2$.

Eq. (4.38) is a polynomial of degree $n$ in the variable $C$ with parity
determined by $n$ (even function in $C$ if $n$ is even, or odd if $n$ is odd).
The result can also be expressed in an equivalent (but less transparent) way
in terms of Gauss’ hypergeometric function $~{}_{2}F_{1}$ as

$\frac{\langle\varphi_{1}^{n}\varphi_{2}^{n}\rangle}{\sigma^{2n}}=\Gamma(n+1)~{}C^{n}~{}_{2}F_{1}\left(\frac{1-n}{2},-\frac{n}{2};1;\frac{1}{C^{2}}\right)\,.$
(4.39)

The same calculation in the framework of Section 2 would be much more
difficult because it would require solving a coupled system of equations
involving both the field and the canonical momentum. Furthermore, we have
already seen for $n=2$ (corresponding to the 4-pt function
$\Delta_{\phi^{2},\phi^{2}}$) that the procedure in Section 2 fails to
resproduce the correct form of non-coincident correlators expected from Wick’s
theorem, but it only succeeds in giving the right values in the coincident
(and very short distance) regime ($r<r_{0}~{}e^{-N}$) and in the large
distance regime ($r>r_{0}$) corresponding respectively to $C=1$ (perfectly
correlated) and $C=0$ (completely uncorrelated).

The result for $C=1$ (coincidence) is readily obtained from Wick’s theorem to
be $\frac{\langle\varphi^{2n}\rangle}{\sigma^{2n}}=(2n-1)!!$. In Section 5 we
will see an application of this coincident value when evaluating the quantum
backreaction of the field on the expansion rate at late times, and how that
affects the energy balance between classical and quantum contributions to the
total energy density of the Universe.

At large distances $r>r_{0}$ (where $C=0$), eq. (4.38) saturates to $0$ if $n$
is odd and $\left[(n-1)!!\right]^{2}$ if $n$ is even, because only the term
$k=0$ contributes to the sum.

The most general 2-point correlator
$\langle\varphi_{1}^{n_{1}}\varphi_{2}^{n_{2}}\rangle$, with $n_{1}$ and
$n_{2}$ arbitrary natural numbers obeys

$\displaystyle\frac{\langle\varphi_{1}^{n_{1}}\varphi_{2}^{n_{2}}\rangle}{\sigma^{n_{1}+n_{2}}}$
$\displaystyle=$ $\displaystyle\sum_{\begin{subarray}{c}k=0\\\ n_{1}-k={\rm
even}\\\ n_{2}-k={\rm even}\end{subarray}}^{{\rm
min}(n_{1},n_{2})}k!{n_{1}\choose k}{n_{2}\choose
k}(n_{1}-k-1)!!(n_{2}-k-1)!!~{}C^{k}\,$ (4.40) $\displaystyle=$
$\displaystyle\sum_{\begin{subarray}{c}k=0\\\ n_{1}-k={\rm even}\\\
n_{2}-k={\rm even}\end{subarray}}^{{\rm
min}(n_{1},n_{2})}\frac{\Gamma(n_{1}+1)\Gamma(n_{2}+1)}{\Gamma\left(\frac{n_{1}-k+2}{2}\right)\Gamma\left(\frac{n_{2}-k+2}{2}\right)}\frac{C^{k}}{2^{\frac{n_{1}+n_{2}}{2}-k}\Gamma(k+1)}\,.$

Note that the sums in eq. (4.40) are non-empty only if $n_{1}$ and $n_{2}$
have the same parity, because $k$ runs over numbers with the same parity as
both $n_{1}$ and $n_{2}$. This is in agreement with the fact that only when
$n_{1}+n_{2}$ is even it is possible to fully decompose the correlator in a
number of contracted Wick pairs. When $n_{1}$ and $n_{2}$ do not have the same
parity the result (4.40) must be simply read as $0$. Of course, due to the
symmetry under exchange of the positions where the fields $\varphi_{1}$ and
$\varphi_{2}$ are evaluated, one could take $n_{1}\leq n_{2}$, but that is not
necessary for the validity of (4.40).

Just like for eq. (4.39), again an expression in terms of Gauss’
hypergeometric function is possible:

$\displaystyle\frac{\langle\varphi_{1}^{n_{1}}\varphi_{2}^{n_{2}}\rangle}{\sigma^{n_{1}+n_{2}}}$
$\displaystyle=$
$\displaystyle\frac{1+(-1)^{n_{1}+n_{2}}}{2}~{}2^{-\frac{|n_{1}-n_{2}|}{2}}~{}\frac{\Gamma\left(\max(n_{1},n_{2})+1\right)}{\Gamma\left(\frac{|n_{1}-n_{2}|}{2}+1\right)}\,$
(4.41) $\displaystyle\times$ $\displaystyle
C^{\min(n_{1},n_{2})}~{}_{2}F_{1}\left(\frac{1-\min(n_{1},n_{2})}{2},-\frac{\min(n_{1},n_{2})}{2};\frac{|n_{1}-n_{2}|}{2}+1;\frac{1}{C^{2}}\right)\,,$

where the coefficient $\frac{1+(-1)^{n_{1}+n_{2}}}{2}$ is basically a
Kronecker delta on the relative parity of $n_{1}$ and $n_{2}$, giving $1$ if
$n_{1}$ and $n_{2}$ have the same parity, and $0$ otherwise. Equivalence
between the functions of $C$ in eqs. (4.40) and (4.41) has been checked with a
symbolic mathematical software up to large values of $n_{1}$ and $n_{2}$,
confirming its correctness. When $n_{1}=n_{2}\equiv n$, then eqs. (4.40–4.41)
reduce to eqs. (4.38–4.39).

##  5 Quantum backreaction

In the most recent stages of evolution, it is not possible to describe the
evolution of the scalar field $\hat{\Phi}$ as a free field in a given
classical background metric. The latter is the approximation used throughout
the previous sections and, as we already mentioned in the quantum field theory
treatment (Section 3), it corresponds to the one-loop approximation for the
energy momentum tensor and the density-density correlator. In recent
cosmological epochs, one cannot neglect the quantum backreaction of the scalar
field on the expansion of the Universe, because the energy-momentum tensor of
the scalar field becomes more and more important and eventually takes the lead
of the expansion. As a consequence, the quantum nature of the gravitational
field has to be taken into account and the metric unavoidably exhibits quantum
fluctuations. In particular, the Hubble expansion rate is a fluctuating
quantity and it is properly represented as a quantum operator at each
spacetime point $\hat{H}\left(t,\vec{x}\right)$, obeying the energy
conservation constraint. Assuming zero spatial curvature, the energy
constraint takes the form of a quantum Friedmann equation:

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

where $\rho_{C}(t)$ is the contribution to energy density coming from
classical matter (baryons and cold dark matter292929More precisely, $\rho_{C}$
only accounts for a fraction of matter because, as we have seen,
$\hat{\rho}_{Q}$ also contains a small part which scales like matter. and
$\hat{\rho}_{Q}(t,\vec{x})$ is the quantum energy density due to the scalar
field (part of which yields dark energy).

The quantum scalar field in turn obeys an equation motion which can be read
from eq. (2.8), but the Hubble rate is now a quantum operator $\hat{H}$, as
well as the Ricci scalar $\hat{R}$. Neglecting explicit spatial gradients, the
equation of motion for $\hat{\Phi}$ reads

$\ddot{\hat{\Phi}}+3\hat{H}\dot{\hat{\Phi}}+m^{2}\hat{\Phi}+\xi\hat{R}\hat{\Phi}=0\,.$
(5.2)

Eqs. (5.1) and (5.2) constitute the leading order approximation in spatial
gradients to the dynamical field equations and the relevant gravitational
constraint equation. This is a standard approximation scheme used for studying
inflationary dynamics, known as the separate Universe approach [39, 40].

Solving the full quantum equation of motion (5.2) goes beyond the scope of
this paper and we will limit ourselves to take into account the quantum
character of the metric in eq. (5.1) via $\hat{H}^{2}$, but we neglect it for
the dynamics of the field $\hat{\Phi}$ in eq. (5.2), which we simply
approximate as

$\ddot{\hat{\Phi}}+3H\dot{\hat{\Phi}}+m^{2}\hat{\Phi}+\xi R\hat{\Phi}=0\,.$
(5.3)

In other words, we only include the quantum backreaction effect in the gravity
(energy) constraint equation (5.1), but not in the full dynamics of the field.
Note that eq. (5.3), which is the leading order approximation of eq. (5.2), is
linear in $\hat{\Phi}$ and therefore it preserves Gaussianity of the initial
state, while eq. (5.2) necessarily generates non-Gaussianities.

Analogously to the expansion that we will develop in this Section, one could
also write down an expansion for the operators $\hat{H}$ and $\hat{R}$
appearing in eq. (5.2). By defining
$\bar{H}\equiv\sqrt{\langle\hat{H}^{2}\rangle}$, then, at linear order in the
(squared) Hubble rate fluctuation
$\hat{\delta}_{H^{2}}\equiv\frac{\hat{H}^{2}-\bar{H}^{2}}{\bar{H}^{2}}$, one
gets
$\hat{H}\simeq\bar{H}\left(\hat{1\\!\\!1}+\frac{1}{2}\hat{\delta}_{H^{2}}\right)$.
Similarly the Ricci operator $\hat{R}=12\hat{H}^{2}+6\dot{\hat{H}}$ can also
be expanded using the second Friedmann equation for $\dot{\hat{H}}$.

As we will see, the approximation made when limiting ourselves to eqs. (5.1)
and (5.3) will still allow us to find some interesting results from the
resummation of scalar loops, arising from the energy conservation equation
(5.1).

As a first step, let us recall that, in the model with action (2.1), the
quantum energy density operator receives contributions quadratic in the
field/canonical momentum. Neglecting the subdominant
$\\{\hat{\Phi},\hat{\Pi}\\}$ and $\hat{\Pi}^{2}$, neglecting spatial gradients
and considering a negative non-minimal coupling $\xi<0$, the quantum energy
density, which has to be used in eq. (5.1), can be approximated from eq.
(2.19) as

$\hat{\rho}_{Q}=\left(\frac{m^{2}}{2}\hat{1\\!\\!1}-3|\xi|\hat{H}^{2}\right)\hat{\Phi}^{2}\,,$
(5.4)

where $\hat{1\\!\\!1}$ is the identity operator acting on the Hilbert space of
quantum states. At leading order in $|\xi|$, the pressure in matter
era303030We are using eq. (2.18) without taking expectation values and with
$\epsilon=3/2$ (matter-dominated epoch). Spatial gradients have been
neglected. is

$\hat{p}_{Q}=-\frac{m^{2}}{2}\hat{\Phi}^{2}\,.$ (5.5)

As we see by comparing energy density and pressure, there are two
contributions to the energy density: one ($m^{2}/2~{}\hat{\Phi}^{2}$) behaves
like a cosmological constant (CC) (i.e. equation of state parameter $-1$) and
the other ($-3|\xi|\hat{H}^{2}\hat{\Phi}^{2}$) mimics a non-relativistic
matter-like contribution (pressureless).

An important detail that should not be underestimated is that, in the right-
hand side of (5.4), the (squared) Hubble parameter itself appears as an
operator (as denoted by its hat). This is the proper way to take into account
the fact that, due to the intrinsic quantum nature of the field $\hat{\Phi}$,
the local expansion rate necessarily inherits a quantum behavior, making it a
fluctuating quantity. Indeed, the squared Hubble parameter is properly
represented by an operator $\hat{H}^{2}(t,\vec{x})$, which obeys eq. (5.1).
But then, since $\hat{H}^{2}$ is a quantum operator in (5.1), it should also
be treated in the same way when it appears in (5.4).

The quantum nature of $\hat{H}^{2}$ in (5.4) is the quantum backreaction
effect that we want to discuss in this Section. A consequence of it, which is
readily obtained by combining (5.4) and (5.1), is the local relation between
expansion rate and field

$\hat{H}^{2}=\frac{1}{3M_{P}^{2}}\left(\rho_{C}\hat{1\\!\\!1}+\frac{m^{2}}{2}\hat{\Phi}^{2}\right)\left(\hat{1\\!\\!1}+|\xi|\frac{\hat{\Phi}^{2}}{M_{P}^{2}}\right)^{-1}\,,$
(5.6)

where the last (inverse) operator is exactly due to the quantum backreaction,
which therefore renormalizes the Planck mass $M_{P}$. For simplicity of
notation we omitted space-time dependence.

We would like to match our quantum DE model to the Universe as it appears
today $t=t_{0}$, which is (mostly) made up by non-relativistic matter (cold
dark matter and baryons) and dark energy well described by a cosmological
constant. Firstly, the expectation value of $\hat{H}^{2}(t_{0},\vec{x})$
should be a proxy for the global cosmological squared Hubble parameter today
$H_{0}^{2}$,

$\langle\hat{H}^{2}\rangle=H_{0}^{2}\,.$ (5.7)

This condition basically defines what we mean by the (squared) Hubble
parameter today within the quantum model. For simplicity we do not write
explicitly that quantities and expectation values are considered at the
current time $t=t_{0}$. Furthermore spatial homogeneity implies that
conditions do not depend on the comoving position $\vec{x}$.

As a second step, we require that the CC part of $\hat{\rho}_{Q}$ should
account for the cosmological constant observed, therefore

$\frac{m^{2}}{2}\langle\hat{\Phi}^{2}\rangle=3M_{P}^{2}H_{0}^{2}\Omega_{\Lambda}\,,$
(5.8)

where $M_{P}=(8\pi G)^{-1/2}$ is the reduced Planck mass and
$\Omega_{\Lambda}$ is the fraction of energy density today due to a
cosmological constant. Recalling the definition of the positive dimensionless
quantity $\alpha$ (see eq. (2.103 and use $H_{\rm
DE}=H_{0}\sqrt{\Omega_{\Lambda}}$)

$\alpha\equiv\frac{1}{6|\xi|\Omega_{\Lambda}}\left(\frac{m}{H_{0}}\right)^{2}\,,$
(5.9)

then eq. (5.8) gives

$\langle\hat{\Phi}^{2}\rangle=\frac{M_{P}^{2}}{\alpha|\xi|}\equiv\sigma^{2}\,,$
(5.10)

which is also the variance of the field
$\sigma^{2}\equiv\langle\hat{\Phi}^{2}\rangle-\langle\hat{\Phi}\rangle^{2}$
because we always assume $\langle\hat{\Phi}\rangle=0$.

The classical energy density today $\rho_{C}$ is then determined by eq. (5.7)
where $\hat{H}^{2}$ obeys the Friedmann equation (5.6).

We then use the formal geometric series expansion

$\left(\hat{1\\!\\!1}+\frac{|\xi|}{M_{P}^{2}}\hat{\Phi}^{2}\right)^{-1}=\sum_{n=0}^{\infty}\left(-|\xi|\frac{\hat{\Phi}^{2}}{M_{P}^{2}}\right)^{n}$
(5.11)

and the results from Wick’s theorem at coincidence in the previous Section,
$\langle\hat{\Phi}^{2n}\rangle=\sigma^{2n}(2n-1)!!$ and similarly
$\langle\hat{\Phi}^{2n+2}\rangle=\sigma^{2n+2}(2n+1)!!$. Introducing the
classical energy density fraction today

$\Omega_{C}\equiv\frac{\rho_{C}}{3M_{P}^{2}H_{0}^{2}}\,,$ (5.12)

one gets after a few steps

$\langle\hat{H}^{2}\rangle=H_{0}^{2}\sum_{n=0}^{\infty}\left(-\frac{1}{\alpha}\right)^{n}\left[\Omega_{C}(2n-1)!!+\Omega_{\Lambda}(2n+1)!!\right]\,,$
(5.13)

and reindexing the second part of the series, the result can be written as

$\langle\hat{H}^{2}\rangle=H_{0}^{2}\left[\alpha\Omega_{\Lambda}+\left(\Omega_{C}-\alpha\Omega_{\Lambda}\right)\sum_{n=0}^{\infty}\left(-\frac{1}{\alpha}\right)^{n}(2n-1)!!\right]\,.$
(5.14)

This contains a formal power series in $1/\alpha$, or equivalently an
asymptotic expansion at infinity in the variable $\alpha$. The convergence
radius in the variable $1/\alpha$ is actually zero, meaning that, written like
that, the series diverges for every finite $\alpha$. Nevertheless, if we can
find a resummation for it, we would then be able to fix $\Omega_{C}$ (and
therefore $\rho_{C}$) in terms of the other parameters of the problem by
imposing the condition (5.7). A clever way to do this resummation consists in
taking an alternative route which does not rely on any series expansion. Of
course, at the end the two methods should agree, in the sense that the
asymptotic expansion at $\alpha\to\infty$ of the result that we are going to
find, despite not convergent, should match exactly (5.14).

Note that Wick’s theorem is nothing but a way to say that $\hat{\Phi}$ follows
a Gaussian distribution. More precisely we can consider a classical stochastic
variable with the same statistical properties as the quantum field
$\hat{\Phi}$. Indeed that is the picture that one should have in mind for what
we called $\varphi$ throughout the previous Section 4, where we used
Starobinsky’s stochastic formalism. We remind that we are assuming zero VEV
for $\hat{\Phi}$ and therefore, like in eq. (4.6), we can say that the
equivalent classical stochastic variable $\varphi$ follows a distribution
$P(\varphi)$ which is Gaussian (see the discussion after eq. (5.3) about
preservation of Gaussianity in the approximation adopted here), given by

$P(\varphi)=\frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{\varphi^{2}}{2\sigma^{2}}\right)\,,$
(5.15)

where $\sigma^{2}$ is given in eq. (5.10) and $\varphi$ can assume any real
value. This could be applied at any time as long as $\hat{\Phi}$ is a free
field, so that $\varphi$ is Gaussian. We proved explicitly the validity of
this requirement for inflation in the previous Section and it keeps being true
also in radiation and matter era if one treats $\hat{\Phi}$ as a free field,
in the spirit of perturbation theory.

We are interested in the current cosmological time $t_{0}$, as we want to
determine the value $\rho_{C}$ (or $\Omega_{C}$) today in terms of the model
parameters $\xi$ and $m$ (or equivalently $\xi$ and $\alpha$) and the
cosmological parameters $H_{0}$ and $\Omega_{\Lambda}=1-\Omega_{M}$.

Expectation values can then be evaluated using $P(\varphi)$ and mapping
quantum operators to their corresponding classical stochastic quantities,
which means that $3M_{P}^{2}\langle H^{2}\rangle$ is given by

$3M_{P}^{2}\langle\hat{H}^{2}\rangle=\int_{-\infty}^{+\infty}{\rm
d}\varphi~{}P(\varphi)\frac{\rho_{C}+\frac{m^{2}}{2}\varphi^{2}}{1+|\xi|\frac{\varphi^{2}}{M_{P}^{2}}}\,.$
(5.16)

Using eq. (5.15) and changing the integration variable to the dimensionless
$z\equiv\varphi/\sigma$ the expectation value reduces to

$\langle\hat{H}^{2}\rangle=H_{0}^{2}\int_{-\infty}^{+\infty}{\rm
d}z~{}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{z^{2}}{2}\right)\frac{\Omega_{C}+\Omega_{\Lambda}z^{2}}{1+\frac{z^{2}}{\alpha}}\,,$
(5.17)

whose result can be written in terms of the complementary error
function313131The complementary error function erfc(z) can be defined as ${\rm
erfc}(z)\equiv\frac{2}{\sqrt{\pi}}\int_{z}^{+\infty}{\rm d}t~{}{\rm
e}^{-t^{2}}$. as

$\langle\hat{H}^{2}\rangle=H_{0}^{2}\left[\alpha\Omega_{\Lambda}+\left(\Omega_{C}-\alpha\Omega_{\Lambda}\right){\rm
e}^{\alpha/2}\sqrt{\frac{\pi\alpha}{2}}~{}{\rm
erfc}\left(\sqrt{\frac{\alpha}{2}}\right)\right]\,.$ (5.18)

This is the resummation that we are looking for. In order to check agreement
of eq. (5.18) with the result from application of Wick’s theorem in eq.
(5.14), the following asymptotic expansion for the complementary error
function ${\rm erfc}(x)$ is needed:

${\rm erfc}(x)=\frac{{\rm
e}^{-x^{2}}}{x\sqrt{\pi}}\sum_{n=0}^{\infty}\left(-\frac{1}{2x^{2}}\right)^{n}(2n-1)!!\,.$
(5.19)

Using $x=\sqrt{\alpha/2}$, one immediately finds agreement with eq. (5.14),
confirming the correctness of the derivation. We are now ready to impose the
condition (5.7), which enables us to express $\Omega_{C}$ today as

$\Omega_{C}(t_{0})=\alpha\Omega_{\Lambda}+\left(1-\alpha\Omega_{\Lambda}\right)\left[{\rm
e}^{-\alpha/2}\sqrt{\frac{2}{\alpha\pi}}\frac{1}{{\rm
erfc}\left(\sqrt{\frac{\alpha}{2}}\right)}\right]\,,$ (5.20)

where we explicitly wrote that it refers to a value today.

As a consequence, the quantum energy density fraction today
$\Omega_{Q}(t_{0})\equiv\frac{\langle\hat{\rho}_{Q}\rangle(t_{0})}{3M_{P}^{2}H_{0}^{2}}=1-\Omega_{C}(t_{0})$
is

$\Omega_{Q}^{(\rm
br)}(t_{0})=\left(\Omega_{\Lambda}-\frac{1}{\alpha}\right)\left[{\rm
e}^{-\alpha/2}\sqrt{\frac{2\alpha}{\pi}}\frac{1}{{\rm
erfc}\left(\sqrt{\frac{\alpha}{2}}\right)}-\alpha\right]\,,$ (5.21)

where the superscript “(br)”, standing for “backreaction”, is there to stress
that this is the correct expression taking into account the backreaction due
to the quantum nature of the field $\hat{\Phi}$ which makes the Hubble
parameter itself a stochastic quantity. This result should be compared with
what one gets when neglecting the quantum backreaction effect. Without quantum
backreaction one would get the simplified formula (see eq. (2.105))

$\Omega_{Q}^{(\rm no-br)}(t_{0})=\Omega_{\Lambda}-\frac{1}{\alpha}\,.$ (5.22)

This is true because, when neglecting the quantum nature of $\hat{H}^{2}$ in
(5.4), taking expectation values today gives

$\langle\hat{\rho}_{Q}^{(\rm no-
br)}\rangle=\left(\frac{m^{2}}{2}\hat{1\\!\\!1}-3|\xi|H_{0}^{2}\right)\langle\hat{\Phi}^{2}\rangle\,,$
(5.23)

which, upon substituting (5.10), leads to (5.22).

Comparing eqs. (5.21) and (5.22) we see that the effect of quantum
backreaction is encoded in a multiplicative function of the dimensionless
parameter $\alpha$ that we remind to be defined as
$\alpha\equiv\frac{m^{2}}{6|\xi|\Omega_{\Lambda}H_{0}^{2}}$. The ratio between
$\Omega_{Q}^{(\rm br)}(t_{0})$ and $\Omega_{Q}^{(\rm no-br)}(t_{0})$ is

$\boxed{\frac{\Omega_{Q}^{(\rm br)}(t_{0})}{\Omega_{Q}^{(\rm no-
br)}(t_{0})}={\rm e}^{-\alpha/2}\sqrt{\frac{2\alpha}{\pi}}\frac{1}{{\rm
erfc}\left(\sqrt{\frac{\alpha}{2}}\right)}-\alpha}\,.$ (5.24)

In the limit of large $\alpha$ the ratio goes to unity323232The first few
terms of the expansion of (5.24) in $\frac{1}{\alpha}$ are
$1-\frac{2}{\alpha}+\frac{10}{\alpha^{2}}+{\mathcal{O}}\left(\frac{1}{\alpha^{3}}\right)$.,
which makes sense if one thinks about the limit $|\xi|\to 0$ with $m/H_{0}$
finite. In this case the full quantum energy density behaves like dark energy
(in the form of a cosmological constant), so both (5.21) and (5.22) reduce to
$\Omega_{\Lambda}$. For finite $\alpha$, the ratio between the quantum energy
densities today with or without accounting for the backreaction can be seen in
Fig. 6, where we plot and discuss as a function of $\alpha$. There is a
physical lower limit for $\alpha$, which is obtained by requiring that the
dark energy contribution to the quantum energy density in eq. (5.4) eventually
dominates over its matter-like part. If this is true then at late times in the
future $t\to+\infty$, one should get $\langle\hat{\rho}_{Q}(t)\rangle>0$, or
equivalently,

$\lim_{t\to+\infty}\Omega_{Q}^{(\rm br)}(t)>0\,.$ (5.25)

It is quite straightforward333333Let us imagine to live in another epoch, say
at time $t$, and repeat the matching procedure that we discussed starting from
eq. (5.8) at that time $t$, instead of $t_{0}$. What changes is just that we
should refer to the energy density fraction made by the cosmological constant
at the new time $\Omega_{\Lambda}(t)$ and to the Hubble parameter at that same
time $H(t)$. Note that, for a cosmological constant,
$\Omega_{\Lambda}(t)H^{2}(t)=\rho_{\Lambda}/(3M_{P}^{2})$ is constant because
the energy density $\rho_{\Lambda}$ for a cosmological constant is indeed
constant. Thus there is no time dependence left in the dimensionless
$\alpha(t)\equiv m^{2}/(6|\xi|H^{2}(t)\Omega_{\Lambda}(t))$, which is then the
same as in eq. (5.9),
$\alpha(t)=\alpha=(mM_{P})^{2}/(2|\xi|\rho_{\Lambda})={\rm constant}$. This
explains why the only modification needed in eq. (5.21) to make it valid at
any time $t$ is the replacement of $\Omega_{\Lambda}$ (which referred to
$t_{0}$) with $\Omega_{\Lambda}(t)$. Incidentally, the constancy of $\alpha$
also implies that the result for the ratio in eq. (5.24) is valid at any time
$t$ and not only today (time $t_{0}$). to extend eq. (5.21) to any time $t$
instead of today ($t_{0}$) just by replacing $\Omega_{\Lambda}$ measured today
with the same quantity at time $t$, which is the fraction of energy density
made up by the cosmological constant at time $t$. When $t\to+\infty$ that
fraction approaches $1$, because the entire energy budget of the Universe will
be completely made by the cosmological constant, therefore

$\lim_{t\to+\infty}\Omega_{Q}^{(\rm
br)}(t)=\left(1-\frac{1}{\alpha}\right)\left[{\rm
e}^{-\alpha/2}\sqrt{\frac{2\alpha}{\pi}}\frac{1}{{\rm
erfc}\left(\sqrt{\frac{\alpha}{2}}\right)}-\alpha\right]\,.$ (5.26)

Since the factor in square brackets in (5.26) is always positive (it is
between $0$ and $1$, as plotted in Fig. 6), then the condition (5.25) imposes

$\alpha>1\,.$ (5.27)

Figure 6: The ratio between the average quantum energy densities evaluated
with ($\Omega_{Q}^{(\rm br)}$) or without ($\Omega_{Q}^{(\rm no-br)}$)
accounting for the backreaction due to the quantum nature of the field
$\hat{\Phi}$ (which makes the Hubble parameter itself a stochastic quantity),
as a function of the dimensionless combination of parameters
$\alpha\equiv\frac{(m/H_{0})^{2}}{6|\xi|\Omega_{\Lambda}}$. We denoted by
$H_{0}$ and $\Omega_{\Lambda}$ the Hubble parameter and the cosmological
constant energy density fraction today, while $m$ and $\xi<0$ are the mass and
negative non-minimal coupling of the quantum spectator field $\hat{\Phi}$, as
they appear in eq. (2.1). Only values of $\alpha>1$, starting from the red
vertical line, are physically relevant so that dark energy eventually
completely dominates over matter in the expansion of the Universe. The plotted
ratio, given by the right-hand side of eq. (5.24), saturates to $1$ for large
$\alpha$ values and it is roughly $0.53$ when $\alpha=1$, which is its minimum
value. This means that neglecting the quantum backreaction leads to
overestimating the fraction of energy density due to the quantum field
$\hat{\Phi}$ by a factor of $1.9$ for $\alpha=1$ (and less for larger
$\alpha$).

##  6 Theoretical aspects of the Hubble tension

This section contains the main theoretical steps needed to properly address
the Hubble tension problem in the quantum dark energy model studied throughout
the paper or other generalizations of it. In particular, we will refine some
of the assumptions about probability distributions used in the previous work
[1], where it was shown that the same model considered here can relieve the
Hubble tension. It is known that the $\Lambda$CDM standard cosmological model
the value of the Hubble parameter inferred from early-Universe probes is in
tension (at the level of more than $4\sigma$) with the expansion rate measured
in the local Universe. We refer to them as global and local measurements of
the Hubble parameter, respectively. A brief overview on the Hubble tension
problem and on different measures of the Hubble rate is contained in [1] and
we refer the reader to that work (and references therein) for a more complete
picture. Here we limit ourselves to recall only those concepts strictly
necessary to put the Hubble tension in the context of quantum dark energy
models and turn it into a probability statement. In that respect, the
fundamental idea is that in a Universe where the expansion rate is a
fluctuating quantity like in eq. (5.6), due to the quantum nature of the
underlying fundamental field(s), it is natural that different probes of the
expansion produce a different outcome depending on which length/time scales
they are sensible to. The real matter is understanding how much their results
can be apart and whether their difference is sufficient to explain the Hubble
tension. Following [1] and due to its appearance in the Friedmann equation
(5.6), we focus on the square of the Hubble parameter. Let us suppose that
local measurements of the Hubble parameter (squared) can be seen as a spatial
average over some volume $V_{1}$ at the current cosmological time $t_{0}$, and
similarly we consider the global measurements as a spatial average are over a
larger volume $V_{2}$ (including $V_{1}$) at the same time $t_{0}$343434Real
measurements are always done on the past light cone. However here we do not
aim at describing a precise global or local measurement procedure, but we just
want to encompass the general idea that some kind of averaging is necessarily
part of the answer since a probe explores different regions of the sky. This
physical ingredient is already contained in the simple procedure of spatial
averaging adopted here. For a study of the effect of quantum dark energy
fluctuations (and therefore a fluctuating Hubble rate) on specific
cosmological observables, e. g. supernovae luminosity distances, see [34].. As
a proxy to quantify the Hubble tension, we consider the following conditional
probability (see [1]), with “P” standing for “probability”:

${\rm
P}\left(\left[H^{2}\right]_{V_{1}}>H_{1}^{2}\bigg{|}\left[H^{2}\right]_{V_{2}}<H_{2}^{2}\right)\,,$
(6.1)

where $\left[H^{2}\right]_{V_{i}},~{}(i=1,2)$ denotes the spatial average of
the squared Hubble rate today over the local ($i=1$) or global ($i=2$)
volumes, i.e.

$\left[H^{2}\right]_{V_{i}}\equiv(1/V_{i})\int_{V_{i}}d^{3}x\,H^{2}(t_{0},\vec{x})\,,\qquad(i=1,2)\,.$
(6.2)

The quantity $H^{2}(t_{0},\vec{x})$ is meant as a classical stochastic
variable with the same statistical properties as the quantum operator
$\hat{H}^{2}(t_{0},\vec{x})$ in eq. (5.6), just like we already introduced the
classical stochastic variable $\varphi$ following the same statistics as the
quantum field $\hat{\Phi}$.

As for $H_{1}$ and $H_{2}$, for definiteness we choose the local value $H_{1}$
as the one measured by the SH0ES collaboration [9] from luminosity distances
of supernovae Ia, and the global value $H_{2}$ as the one inferred by the
Planck mission [8] from CMB temperature fluctuations assuming the $\Lambda$CDM
cosmology:

$\displaystyle H_{1}$ $\displaystyle=$ $\displaystyle 74.03~{}{\rm
km/s~{}Mpc^{-1}}\,,$ $\displaystyle H_{2}$ $\displaystyle=$ $\displaystyle
67.4~{}{\rm km/s~{}Mpc^{-1}}\,.$ (6.3)

Keeping in mind that the local value $H_{1}$ is larger than the global value
$H_{2}$ (see eq. (6.3)), we can explain why the probability in eq. (6.1) is a
measure of the Hubble tension: it expresses how likely it is that local
measurements give a result at least as large as the one actually measured,
given that global measurements give a value at least as low as the one they
really provide. If the probability (6.1) is large enough, then the Hubble
tension is relieved by the quantum nature of dark energy. Indeed, showing that
this actually happens was the main result of the letter [1]. The volumes
$V_{1}$ and $V_{2}$ are taken as spheres centered around the observer with
radii $R_{1}$ and $R_{2}$. In [1], we considered a global radius given by the
Hubble horizon $R_{2}=H_{2}^{-1}\simeq 4400$ Mpc and a local radius
$R_{1}=100$ Mpc (we also studied another choice $R_{1}=1$ Gpc for the local
radius to mimic other kinds of measurements from the H0LiCOW collaboration
instead of SH0ES). The precise numerical values are not relevant for what we
will discuss, which only deals with the theoretical approach for evaluating
the probability (6.1).

By making use of Bayes’ theorem, we can rewrite the conditional probability in
eq. (6.1) in terms of a joint probability as

${\rm
P}\left(\left[H^{2}\right]_{V_{1}}>H_{1}^{2}\bigg{|}\left[H^{2}\right]_{V_{2}}<H_{2}^{2}\right)=\frac{{\rm
P}\left(\left[H^{2}\right]_{V_{1}}>H_{1}^{2}\bigcap\left[H^{2}\right]_{V_{2}}<H_{2}^{2}\right)}{{\rm
P}\left(\left[H^{2}\right]_{V_{2}}<H_{2}^{2}\right)}\,.$ (6.4)

Evaluating the denominator requires knowing the distribution of the global
variable $\left[H^{2}\right]_{V_{2}}$, while the numerator requires the joint
probability distribution of the two variables $\left[H^{2}\right]_{V_{1}}$
(local) and $\left[H^{2}\right]_{V_{2}}$ (global). We shall study aspects of
both of them, starting with the denominator, focusing on the methodology to
approach the problem and not on numerical predictions. The latter were already
obtained in ref.[1] thanks to considerable simplifications stemming from the
assumption that the squares of the fields $\varphi^{2}(t_{0},\vec{x})$ can be
approximately treated as Gaussian variables (with correlation between fields
at different positions). Such a simplified approach takes into account the
main physical ingredient for a possible relief of the Hubble tension in the
context of a quantum dark energy model, namely the spatial correlations, and
the main lesson from [1] is that spatially correlated dark energy can actually
help with the Hubble tension issue. Nevertheless we know, e.g. from the
results of Section 4, that the field $\varphi(t_{0},\vec{x})$ it is which
follows a Gaussian distribution, not its square. The distribution of
$\varphi^{2}$ at a given position is a chi-squared distribution, and the
spatially averaged variables have the complication of spatial correlations,
which create dependence between fields at different positions. The major
purpose of the rest of this Section is to explore the possible path to a more
rigorous statistical treatment of the problem. As a general rule for notation,
capital letter “P” will be reserved for finite probabilities (i.e. integrated
over a region), while we will use the small “$p$” for probability distribution
functions (PDF).

###  6.1 Distribution of the global squared Hubble rate

As a first step towards calculating the denominator ${\rm
P}\left(\left[H^{2}\right]_{V_{2}}<H_{2}^{2}\right)$ of eq. (6.4), we
discretize space in the global volume $V_{2}$ by a large set of $N_{2}$ points
$\left\\{\vec{x}_{i}\right\\}_{i=1,2,\dots,N_{2}}$ evenly distributed in the
volume $V_{2}$. This will make some manipulations more easy to understand, in
the same spirit in which, for example, in the path integral approach to
quantum mechanics it is useful to deal with an integral over a finite number
of variables and only later take the continuum limit to get a functional
integration. We will discuss some formal aspects of the continuum limit,
however in practice a discretized approach is more suitable for numerical
computations. Since we always refer to quantities at the current cosmological
time $t_{0}$ we will sometimes omit it from notation, for simplicity.

Due to discretization, the definition of the global squared Hubble parameter
is now,

$\left[H^{2}\right]_{2}\equiv\frac{1}{N_{2}}\sum_{i=1}^{N_{2}}H^{2}_{i}\,,$
(6.5)

where $H^{2}_{i}\equiv H^{2}\left(t_{0},\vec{x}_{i}\right)$ is the Hubble rate
at the point $\vec{x}_{i}$. $H_{i}^{2}$ can be considered as a local expansion
parameter (squared), defined by the geodesic deviation equation or the
Raychaudhuri equation. Using eq. (5.6), one can relate the squared Hubble rate
$H_{i}^{2}$ and the classical field variable $\varphi_{i}$ at $\vec{x}_{i}$.
It is convenient to work with dimensionless variables to simplify as much as
possible the notation and the following calculations. Thus we introduce the
dimensionless field variable

$z_{i}\equiv\varphi_{i}/\sigma\,,$ (6.6)

where $\sigma^{2}=\langle\varphi_{i}^{2}\rangle$ is the variance of the field
already introduced in eq. (5.10)(independent from the position because of
statistical spatial homogeneity). Note that a similar variable $z$ was already
introduced for the integration in eq. (5.17). The Friedmann equation (5.6) at
$\vec{x}_{i}$ becomes

$\frac{H_{i}^{2}}{H_{0}^{2}}=\frac{\Omega_{C}+\Omega_{\Lambda}z_{i}^{2}}{1+\frac{z_{i}^{2}}{\alpha}}\,,$
(6.7)

where $\alpha$ was defined in eq. (5.9) and $H_{0}^{2}$ is the expectation
value of the squared Hubble rate today (again space-independent). As we
discussed in Section 5, just before eq. (5.20), the classical energy density
fraction $\Omega_{C}$ is determined by requiring that $\langle
H_{i}^{2}\rangle=H_{0}^{2}$ (where the ensemble average $\langle\dots\rangle$
coincides with the quantum state average), resulting in the expression (5.20).
As for the value of $H_{0}$ (which is again marginal for our mostly
theoretical discussion), we can suppose that $H_{2}$ in (6.3), being an
average result over the large volume $V_{2}$, is a good approximation for
$H_{0}$, so that

$H_{0}\simeq H_{2}\,.$ (6.8)

Introducing the dimensionless variable

$h_{i}\equiv\frac{H_{i}^{2}}{H_{0}^{2}}\,,$ (6.9)

eq. (6.7) is equivalent to

$h_{i}=\frac{\Omega_{C}+\Omega_{\Lambda}z_{i}^{2}}{1+\frac{z_{i}^{2}}{\alpha}}\,.$
(6.10)

It follows from $\langle H_{i}^{2}\rangle=H_{0}^{2}$ that $\langle
h_{i}\rangle=1$, which is guaranteed for $\Omega_{C}$ given by eq. (5.20).
After introducing the average variable

$\left[h\right]_{2}\equiv\frac{1}{N_{2}}\sum_{i=1}^{N_{2}}h_{i}\,,$ (6.11)

where the subscript “2” is meant to remind that it refers to an average over
$N_{2}$ points (i.e. the global average), the probability ${\rm
P}\left(\left[H^{2}\right]_{V_{2}}<H_{2}^{2}\right)$ is, in the discretized
picture and with $H_{0}\simeq H_{2}$, equal to ${\rm
P}\left(\left[h\right]_{2}<1\right)$.

Evaluating such a probability is possible if one manages to find how
$\left[h\right]_{2}$ is distributed. Note that $\langle h_{i}\rangle=1$
implies

$\langle\left[h\right]_{2}\rangle=1\,.$ (6.12)

However higher momenta (determining variance, skewness, kurtosis, etc…) are
not so trivial to determine due to the correlations between different $h_{i}$,
whose origin is in the non-zero off-coincident quantum correlations. We know
that the variables $\left\\{z_{i}\right\\}_{i=1,2,\dots,N_{2}}$ are simply
obtained by rescaling the fields
$\left\\{\varphi_{i}\right\\}_{i=1,2,\dots,N_{2}}$ by the constant $\sigma$.
Thus their distribution is readily determined by the the distribution of
fields. In Section 4 we confirmed, by explicitly solving the Fokker-Planck
equation (4.8), that the fields at a given couple of points separated by
comoving distance $r$ (and at equal times) follow a Gaussian distribution with
correlation coefficient $C(r)$ imprinted by inflation as in eq. (4.34). This
is a consequence of the fields being free. Due to the rescaling (6.6), it is
clear from eq. (4.32) that the correlator $\langle z_{i}z_{j}\rangle$ is given
by

$C_{ij}\equiv\langle z_{i}z_{j}\rangle=C(\|\vec{x}_{i}-\vec{x}_{j}\|)\,.$
(6.13)

In particular, when $i=j$ it correctly gives a variance $C_{ii}\equiv\langle
z_{i}^{2}\rangle=C(0)=1$ in agreement with $z_{i}=\varphi_{i}/\sigma$ and
$\langle\varphi_{i}^{2}\rangle=\sigma^{2}$. Since we are now dealing with
fields at $N_{2}$ points it is straightforward to write down a joint
multivariate Gaussian distribution such that each pair of points
$\\{z_{i},z_{j}\\}$ has a correlator $C_{ij}$ given by eq. (6.13). This is
reached by

$p\left(z_{1},z_{2},\dots,z_{N_{2}}\right)=\frac{1}{\sqrt{(2\pi)^{N_{2}}{\rm
det}(C)}}\exp{\left(-\frac{1}{2}z^{T}\cdot C^{-1}\cdot z\right)}\,,$ (6.14)

where $z^{T}$ collects all the variables $z_{i}$ in a row vector as
$z^{T}\equiv\left(z_{1},z_{2},\dots,z_{N_{2}}\right)$ and $z$ is the
corresponding column vector. We also introduced the symmetric and positive
definite $N_{2}\times N_{2}$ matrix $C$ with elements $C_{ij}$, its inverse
$C^{-1}$ and used its determinant ${\rm det}(C)$ in the normalization
coefficient to ensure that

$\int_{-\infty}^{+\infty}{\rm d}z_{1}\int_{-\infty}^{+\infty}{\rm
d}z_{2}~{}\dots\int_{-\infty}^{+\infty}{\rm
d}z_{N_{2}}~{}p\left(z_{1},z_{2},\dots,z_{N_{2}}\right)=1\,.$ (6.15)

We can now formulate in simple terms the temporary goal problem as follows,
making reference to eqs. (6.10), (6.11) and (6.14): we are interested in the
probability distribution function of the average variable $\left[h\right]_{2}$
defined in eq. (6.11), where each $h_{i}$ is given in terms of $z_{i}$ by eq.
(6.10), knowing that the variables
$\left\\{z_{i}\right\\}_{i=1,2,\dots,N_{2}}$ follow the Gaussian distribution
in eq. (6.14).

A useful tool to find the distribution of $\left[h\right]_{2}$ is its moment
generating function (MGF), defined as

$M_{\left[h\right]_{2}}(s)\equiv\langle{\rm
e}^{s\left[h\right]_{2}}\rangle\,,$ (6.16)

where $s$ is an auxiliary variable. As a reminder of the main features of the
MGF, all the moments of the distribution for $\left[h\right]_{2}$ are encoded
in $M_{\left[h\right]_{2}}(s)$ because

$M_{\left[h\right]_{2}}(s)=\sum_{n=0}^{\infty}\langle\left[h\right]_{2}^{n}\rangle\frac{s^{n}}{n!}\,,$
(6.17)

which means that

$\langle\left[h\right]_{2}^{n}\rangle=\frac{{\rm d}^{n}}{{\rm
d}s^{n}}M_{\left[h\right]_{2}}(s)\bigg{|}_{s=0}\,.$ (6.18)

A good strategy to determine the PDF of $\left[h\right]_{2}$ is to first
compute its MGF as

$\displaystyle M_{\left[h\right]_{2}}(s)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{+\infty}{\rm d}z_{1}\int_{-\infty}^{+\infty}{\rm
d}z_{2}~{}\dots\int_{-\infty}^{+\infty}{\rm
d}z_{N_{2}}~{}p\left(z_{1},z_{2},\dots,z_{N_{2}}\right){\rm
e}^{s\left[h\right]_{2}\left(z_{1},z_{2},\dots,z_{N_{2}}\right)}\,$ (6.19)
$\displaystyle=$ $\displaystyle\int\frac{{\rm d}z_{1}{\rm d}z_{2}\dots\rm
dz_{N_{2}}}{\sqrt{(2\pi)^{N_{2}}{\rm
det}(C)}}~{}\exp{\left[-\frac{1}{2}z^{T}\cdot C^{-1}\cdot
z+s\frac{1}{N_{2}}\sum_{i=1}^{N_{2}}\frac{\Omega_{C}+\Omega_{\Lambda}z_{i}^{2}}{1+\frac{z_{i}^{2}}{\alpha}}\right]}\,,$

where in the second line a shorthand notation for the integration over all
variables $z_{i}$ was used. However, the integral above is too complicated to
be computed exactly due to the denominators of the form
$1+\frac{z_{i}^{2}}{\alpha}$ inside the exponential, which are in turn
originated by the quantum backreaction discussed in the previous Section. As
we have seen explicitly in the computation of $\Omega_{Q}=1-\Omega_{C}$
leading to Fig. 6, the backreaction effect decreases with larger $\alpha$ and
gives a correction by a factor of roughly $0.5$ when $\alpha$ is close to its
minimum allowed value of $1$. Since we have seen a case where the quantum
backreaction effect does not change the results by orders of magnitude, but
only by a factor between 0.5 and 1, it is worth to use also here some
approximations to reduce the complication of integration, so that one can
extract some analytical understanding. Therefore let us approximate353535Note
that $z_{i}$ follows a Gaussian distribution centered at zero and with unitary
variance while $\alpha$ is at least 1. Therefore it is not so likely that
$\frac{z_{i}^{2}}{\alpha}$ is large, which justifies using the truncated
expansion in eq. (6.20). Clearly the approximation works at its best for large
$\alpha$ values.,

$h_{i}=\frac{\Omega_{C}+\Omega_{\Lambda}z_{i}^{2}}{1+\frac{z_{i}^{2}}{\alpha}}=\Omega_{C}+\left(\Omega_{\Lambda}-\frac{\Omega_{C}}{\alpha}\right)z_{i}^{2}+{\mathcal{O}}\left(z_{i}^{4}\right)\,,$
(6.20)

which gives

$\left[h\right]_{2}\simeq\Omega_{C}+\left(\Omega_{\Lambda}-\frac{\Omega_{C}}{\alpha}\right)\frac{1}{N_{2}}\sum_{i=1}^{N_{2}}z_{i}^{2}\,.$
(6.21)

We remark that this simplification does not obey exactly the result
$\langle\left[h\right]_{2}\rangle$ in eq.(6.12), indeed because of its
approximate nature. From eq. (6.21) and $\langle z_{i}^{2}\rangle=1$ one gets
$\langle\left[h\right]_{2}\rangle\simeq\Omega_{\Lambda}+\Omega_{C}\left(1-\frac{1}{\alpha}\right)$,
which only goes to $1$ in the limit of large $\alpha$, since then $\Omega_{C}$
approaches $\Omega_{M}=1-\Omega_{\Lambda}$.

Using the approximation (6.21), the MGF in eq. (6.19) simplifies to

$\displaystyle M_{\left[h\right]_{2}}(s)$ $\displaystyle\simeq$
$\displaystyle\int\frac{{\rm d}z_{1}{\rm d}z_{2}\dots\rm
dz_{N_{2}}}{\sqrt{(2\pi)^{N_{2}}{\rm
det}(C)}}~{}\exp{\left\\{-\frac{1}{2}z^{T}\cdot C^{-1}\cdot
z+s\left[\Omega_{C}+\left(\Omega_{\Lambda}-\frac{\Omega_{C}}{\alpha}\right)\frac{1}{N_{2}}\sum_{i=1}^{N_{2}}z_{i}^{2}\right]\right\\}}\,$
(6.22) $\displaystyle=$ $\displaystyle\frac{{\rm e}^{s\Omega_{C}}}{\sqrt{{\rm
det}\left[{\rm
Id}_{2}-2s\left(\Omega_{\Lambda}-\frac{\Omega_{C}}{\alpha}\right)\frac{1}{N_{2}}C\right]}}\,,$

where the integral is solved in the last line and ${\rm Id}_{2}$ is the
$N_{2}\times N_{2}$ identity matrix.

This approximate result, which basically depends on the characteristic
polynomial of the covariance matrix $C$, can be used to extract some insight
on the continuum limit. An expansion of eq. (6.22) is in order for this
purpose363636It comes from the matrix identities ${\rm det}(A)=\exp{\left[{\rm
Tr}\log(A)\right]}$ and $-\log\left({\rm
Id}_{2}-X\right)=\sum_{n=1}^{\infty}\frac{1}{n}X^{n}$ applied to $A\equiv{\rm
Id}_{2}-X$ and $X\equiv
2s\left(\Omega_{\Lambda}-\frac{\Omega_{C}}{\alpha}\right)\frac{1}{N_{2}}C$.:

$M_{\left[h\right]_{2}}(s)\simeq\exp{\left\\{\Omega_{C}s+\frac{1}{2}\sum_{n=1}^{\infty}\frac{1}{n}\left[2\left(\Omega_{\Lambda}-\frac{\Omega_{C}}{\alpha}\right)s\right]^{n}\frac{1}{N_{2}^{n}}{\rm
Tr}\left(C^{n}\right)\right\\}}\,,$ (6.23)

where ${\rm Tr}\left(C^{n}\right)$ is the trace of the n-th power of the
matrix $C$. The continuum limit is obtained from the following observation on
the coefficients $f_{n}\equiv\frac{1}{N_{2}^{n}}{\rm Tr}\left(C^{n}\right)$:

$\displaystyle f_{n}\equiv\frac{1}{N_{2}^{n}}{\rm
Tr}\left(C^{n}\right)=\frac{1}{N_{2}^{n}}\sum_{i_{1}=1}^{N_{2}}\sum_{i_{2}=1}^{N_{2}}\dots\sum_{i_{n}=1}^{N_{2}}C_{i_{1}i_{2}}C_{i_{2}i_{3}}\dots
C_{i_{n}i_{1}}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\rightarrow\frac{1}{V_{2}^{n}}\int_{V_{2}}{\rm
d}^{3}\vec{x}_{1}\int_{V_{2}}{\rm d}^{3}\vec{x}_{2}~{}\dots\int_{V_{2}}{\rm
d}^{3}\vec{x}_{n}~{}C(\|\vec{x}_{1}-\vec{x}_{2}\|)C(\|\vec{x}_{2}-\vec{x}_{3}\|)\dots
C(\|\vec{x}_{n}-\vec{x}_{1}\|)$ $\displaystyle\,.$ (6.24)

Note that $f_{n}\leq 1$ (with equality holding for $n=1$) because the function
$C$ lies between $0$ and $1$, see eq. (4.34) and Fig. 3. Despite the elegant
obtention of the continuum limit via the procedure outlined in eq. (6.24), for
numerical calculations it is still preferable to use the discretized version
(6.22) and check convergence for large $N_{2}$, instead of trying to resum the
series (6.23). Furthermore, the analytical computation of the integrals in the
last line of eq. (6.24) is trivial when $n=1$ and still doable for $n=2$, but
it becomes already too complicated for $n=3$, so that Monte Carlo integration
has to be used.

For $n=1$ one gets $f_{1}=1$, while for $n=2$, the coefficient
$f_{2}=\frac{1}{V_{2}^{2}}\int_{V_{2}}{\rm d}^{3}\vec{x}_{1}\int_{V_{2}}{\rm
d}^{3}\vec{x}_{2}~{}C^{2}(\|\vec{x}_{1}-\vec{x}_{2}\|)$, where $V_{2}$ is the
spherical global volume $V_{2}=\frac{4}{3}\pi R_{2}^{3}$, can be reduced to a
single integration on the relative distance
$r\equiv\|\vec{x}_{2}-\vec{x}_{1}\|$. More precisely, the computation can be
done373737As a summary of the main steps involved, one goes from the
coordinates $\vec{x}_{1}$ and $\vec{x}_{2}$, measured with respect to the
origin located at the center of the sphere, to the relative and center-of-mass
coordinates, given by $\vec{r}\equiv\vec{x}_{2}-\vec{x}_{1}$ and
$\vec{R}_{CM}\equiv\left(\vec{x}_{1}+\vec{x}_{2}\right)/2$, respectively. Then
the integration over the center-of-mass coordinate $\vec{R}_{CM}$ is done for
a given relative distance $\vec{r}$. This must take into account the
geometrical limits $\|\vec{R}_{CM}+\frac{1}{2}\vec{r}\|\leq R_{2}$ and
$\|\vec{R}_{CM}-\frac{1}{2}\vec{r}\|\leq R_{2}$, which can be seen as the
intersection of two spheres, so the problem is essentially reduced to
Euclidean geometry. The relative distance $r=\|\vec{r}\|$ between two points
inside a sphere of radius $R_{2}$ can take values between $0$ and $2R_{2}$ and
its distribution function gives eq. (6.26). This is properly normalized as a
one-dimensional PDF with $\int_{0}^{2R_{2}}{\rm d}r~{}p(r)=1$ and we also
checked its expression numerically by randomly generating pairs of points
inside the sphere. in terms of the probability distribution $p(r)$ for the
relative distance between two points randomly chosen with uniform probability
inside the sphere of radius $R_{2}$:

$f_{2}=\int_{0}^{2R_{2}}{\rm d}r~{}p(r)C^{2}(r)\,,$ (6.25)

where

$p(r)=\frac{12}{R_{2}}\left[\left(\frac{r}{2R_{2}}\right)^{2}-\frac{3}{2}\left(\frac{r}{2R_{2}}\right)^{3}+\frac{1}{2}\left(\frac{r}{2R_{2}}\right)^{5}\right]\,,\qquad{\rm
defined~{}for}~{}r\in[0,2R_{2}]\,.$ (6.26)

The moments of the distribution for $\left[h\right]_{2}$ can be extracted from
the series (6.23) and computed in terms of $f_{n}$’s. It is actually more
convenient to work with cumulants so that the PDF of $\left[h\right]_{2}$ can
be approximated by series expansions like the Gram-Charlier type A series
[63]–[64], which however is not always guaranteed to converge. The cumulant
generating function (CGF) $K_{\left[h\right]_{2}}(s)$ is defined as the
logarithm of the MGF and the cumulants $\kappa_{n}$ are defined from the
coefficients of the expansion in $s$, as383838In the language of quantum field
theory (and up to factors of $i$), the MGF corresponds to the partition
function $Z[J]$ (generating functional of Green’s functions), while the CGF
corresponds to generating functional of connected Green’s functions $W[J]$
which obeys $Z[J]={\rm e}^{iW[J]}$, with the external current $J(x)$ playing
the role of the source variable $s$. In this sense the cumulants are the
connected parts of correlators.

$\log{\left[M_{\left[h\right]_{2}}(s)\right]}\equiv
K_{\left[h\right]_{2}}(s)\equiv\sum_{n=1}^{\infty}\kappa_{n}\frac{s^{n}}{n!}\,.$
(6.27)

Therefore, the cumulants can be read immediately from the argument of the
exponential in (6.23). The first one $\kappa_{1}$ is the mean value393939See
also the discussion on $\langle\left[h\right]_{2}\rangle$ right after eq.
(6.21).,
$\kappa_{1}=\langle\left[h\right]_{2}\rangle\simeq\Omega_{\Lambda}+\Omega_{C}\left(1-\frac{1}{\alpha}\right)$.
The second one $\kappa_{2}$ is the variance,
$\kappa_{2}=\langle\left[h\right]_{2}^{2}\rangle-\langle\left[h\right]_{2}\rangle^{2}=2\left(\Omega_{\Lambda}-\frac{\Omega_{C}}{\alpha}\right)^{2}f_{2}$.
For $n\geq 2$ the cumulants are

$\kappa_{n}=\frac{\Gamma(n)}{2}\left[2\left(\Omega_{\Lambda}-\frac{\Omega_{C}}{\alpha}\right)\right]^{n}f_{n}\,,\qquad
n\geq 2\,.$ (6.28)

The Gram-Charlier series gives the distribution
$p\left(\left[h\right]_{2}\right)$ as an expansion in cumulants around a
reference distribution, usually taken to be Gaussian. Truncating for
simplicity to the $4^{\rm th}$ cumulant, the PDF of $\left[h\right]_{2}$ is

$p\left(\left[h\right]_{2}\right)\simeq\frac{1}{\sqrt{2\pi\kappa_{2}}}\exp{\left[-\frac{1}{2}\frac{\left(\left[h\right]_{2}-\kappa_{1}\right)^{2}}{\kappa_{2}}\right]}\left[1+\frac{\kappa_{3}}{6\kappa_{2}^{3/2}}H_{3}\left(\frac{\left[h\right]_{2}-\kappa_{1}}{\sqrt{\kappa_{2}}}\right)+\frac{\kappa_{4}}{24\kappa_{2}^{2}}H_{4}\left(\frac{\left[h\right]_{2}-\kappa_{1}}{\sqrt{\kappa_{2}}}\right)+\dots\right]\,,$
(6.29)

where $H_{3}(x)=x^{3}-3x$ and $H_{4}(x)=x^{4}-6x^{2}+3$ are Hermite
polynomials and the dots denote terms dependent on higher order cumulants with
perfectly known structure [63]–[64], but that we omit for brevity. The
physical information about dark energy spatial correlations is encoded in the
cumulants $\kappa_{n}$ of eq. (6.28) via the integrals $f_{n}$ in eq. (6.24).

Alternatively, a technique to recover the PDF,
$p\left(\left[h\right]_{2}\right)$, which does not rely on the cumulant
expansion, is via numerical evaluation of the inverse Laplace transform of the
MGF $M_{\left[h\right]_{2}}(s)$. The starting point is the relation
implied404040The rigorous limits of integration, i. e. the possible values of
$\left[h\right]_{2}$, can be deduced from eqs. (6.10), (6.11) and are given by
$\left[h\right]_{2}\in[{\rm
min}\left(\Omega_{C},\alpha\Omega_{\Lambda}\right),{\rm
max}\left(\Omega_{C},\alpha\Omega_{\Lambda}\right)]$ depending on the sign of
$\Omega_{\Lambda}-\frac{\Omega_{C}}{\alpha}$. When $\alpha$ is not too close
to $1$, so that the approximations (6.20) and (6.21) work well (better for
large $\alpha$), $\Omega_{\Lambda}-\frac{\Omega_{C}}{\alpha}>0$ and the MGF in
eq. (6.22) can be fully trusted, the range of values for $\left[h\right]_{2}$
can be taken as $\left[h\right]_{2}\in[\Omega_{C},+\infty)$. by eq. (6.16),

$M_{\left[h\right]_{2}}(s)\simeq\int_{\Omega_{C}}^{+\infty}{\rm
d}\left[h\right]_{2}~{}{\rm
e}^{s\left[h\right]_{2}}~{}p\left(\left[h\right]_{2}\right)\,.$ (6.30)

After a shift of variable and recalling the definition of Laplace transform

${\mathcal{L}}\left\\{f(x)\right\\}(s)\equiv\int_{0}^{+\infty}{\rm d}x~{}{\rm
e}^{-sx}f(x)\,,$ (6.31)

eq. (6.30) can be recast into

$M_{\left[h\right]_{2}}(s)\simeq{\rm
e}^{s\Omega_{C}}{\mathcal{L}}\left\\{p(x+\Omega_{C})\right\\}(-s)\,,$ (6.32)

where $x$ is a dummy variable, and then eq. (6.32) is inverted as

$p\left(\left[h\right]_{2}\right)\simeq{\mathcal{L}}^{-1}\left\\{{\rm
e}^{s\Omega_{C}}M_{\left[h\right]_{2}}(-s)\right\\}\left(\left[h\right]_{2}-\Omega_{C}\right)\,,$
(6.33)

where ${\mathcal{L}}^{-1}$ is the inverse Laplace transform, which can be
evaluated numerically or expressed as a line integral (called Bromwich or
Fourier-Mellin integral), see e.g. Section 15.12 of [65]. Once the PDF
$p\left(\left[h\right]_{2}\right)$ is known, either from eq. (6.29) or via eq.
(6.33), the goal probability ${\rm
P}\left(\left[H^{2}\right]_{V_{2}}<H_{2}^{2}\right)={\rm
P}\left(\left[h\right]_{2}<1\right)$ can be found by integration.

###  6.2 Joint distribution of the local and global squared Hubble rates

Let us now discuss briefly the numerator ${\rm
P}\left(\left[H^{2}\right]_{V_{1}}>H_{1}^{2}\bigcap\left[H^{2}\right]_{V_{2}}<H_{2}^{2}\right)$
of eq. (6.4). As for the denominator, we discretize space and consider the
same $N_{2}$ points in the global volume $V_{2}$ as before, with $N_{1}<N_{2}$
of them belonging to the local volume $V_{1}$ contained in $V_{2}$. Without
loss of generality, we can suppose that points $\vec{x}_{i}$ with
$i=1,2,\dots,N_{1}$ belong to $V_{1}$ and points $\vec{x}_{i}$ with
$i=N_{1}+1,N_{1}+2,\dots,N_{2}$ lie in $V_{2}$ but outside of $V_{1}$.
Following eq. (6.11) we consider the variables

$\left[h\right]_{1}\equiv\frac{1}{N_{1}}\sum_{i=1}^{N_{1}}h_{i}\,,\qquad\qquad\left[h\right]_{2}\equiv\frac{1}{N_{2}}\sum_{i=1}^{N_{2}}h_{i}\,,$
(6.34)

and, using the identification in eq. (6.8), the probability ${\rm
P}\left(\left[H^{2}\right]_{V_{1}}>H_{1}^{2}\bigcap\left[H^{2}\right]_{V_{2}}<H_{2}^{2}\right)$
is equivalent in the discretized picture to ${\rm
P}\left(\left[h\right]_{1}>\frac{H_{1}^{2}}{H_{2}^{2}}\bigcap\left[h\right]_{2}<1\right)$,
where $H_{1}$ and $H_{2}$ are the local and global measures of the Hubble
parameter in (6.3). Computing such a probability requires the joint PDF,
$p\left(\left[h\right]_{1},\left[h\right]_{2}\right)$, of the local and global
variables $\left[h\right]_{1}$ and $\left[h\right]_{2}$.

Following eq. (6.16) we define the joint MGF for $\left[h\right]_{1}$ and
$\left[h\right]_{2}$ as

$M_{\left[h\right]_{1},\left[h\right]_{2}}(s_{1},s_{2})\equiv\langle{\rm
e}^{s_{1}\left[h\right]_{1}+s_{2}\left[h\right]_{2}}\rangle\,,$ (6.35)

where we introduced two auxiliary variables $s_{1}$ and $s_{2}$.

Eq. (6.17) generalizes to

$M_{\left[h\right]_{1},\left[h\right]_{2}}(s_{1},s_{2})=\sum_{n_{1}=0}^{\infty}\sum_{n_{2}=0}^{\infty}\langle\left[h\right]_{1}^{n_{1}}\left[h\right]_{2}^{n_{2}}\rangle\frac{s_{1}^{n_{1}}}{n_{1}!}\frac{s_{2}^{n_{2}}}{n_{2}!}\,,$
(6.36)

and eq. (6.18) generalizes to

$\langle\left[h\right]_{1}^{n_{1}}\left[h\right]_{2}^{n_{2}}\rangle=\frac{{\rm
d}^{n_{1}+n_{2}}}{{\rm d}s_{1}^{n_{1}}{\rm
d}s_{2}^{n_{2}}}M_{\left[h\right]_{1},\left[h\right]_{2}}(s_{1},s_{2})\Big{|}_{s_{1}=0\atop
s_{2}=0}\,.$ (6.37)

In the same approximations leading to eq. (6.22), it can be shown that the
jont MGF is

$M_{\left[h\right]_{1},\left[h\right]_{2}}(s_{1},s_{2})\simeq\frac{{\rm
e}^{(s_{1}+s_{2})\Omega_{C}}}{\sqrt{{\rm det}\left[{\rm
Id}_{2}-2\left(\Omega_{\Lambda}-\frac{\Omega_{C}}{\alpha}\right)\left(\frac{s_{2}}{N_{2}}C+\frac{s_{1}}{N_{1}}{\rm
Id}_{1}C\right)\right]}}\,,$ (6.38)

where ${\rm Id}_{1}$ is the $N_{2}\times N_{2}$ matrix given by the identity
on the first $N_{1}$ rows and columns (the local subspace) and with zero
entries otherwise. Hence all the three matrices ${\rm Id}_{1}$, ${\rm Id}_{2}$
and $C$ are $N_{2}\times N_{2}$ in size and their combination in eq. (6.38) is
a legitimate operation.

A similar expansion to eq. (6.23) applied to the result (6.38) and the
introduction of the cumulant generating function
$K_{\left[h\right]_{1},\left[h\right]_{2}}(s_{1},s_{2})\equiv\log{\left[M_{\left[h\right]_{1},\left[h\right]_{2}}(s_{1},s_{2})\right]}$,
gives

$K_{\left[h\right]_{1},\left[h\right]_{2}}(s_{1},s_{2})\simeq\Omega_{C}(s_{1}+s_{2})+\frac{1}{2}\sum_{n=1}^{\infty}\frac{1}{n}\left[2\left(\Omega_{\Lambda}-\frac{\Omega_{C}}{\alpha}\right)\right]^{n}{\rm
Tr}\left[\left(\frac{s_{2}}{N_{2}}C+\frac{s_{1}}{N_{1}}{\rm
Id}_{1}C\right)^{n}\right]\,.$ (6.39)

The joint cumulants $\kappa_{n_{1},n_{2}}$ are retrieved from the expansion

$K_{\left[h\right]_{1},\left[h\right]_{2}}(s_{1},s_{2})\equiv\sum_{n_{1}=0}^{\infty}\sum_{n_{2}=0}^{\infty}\kappa_{n_{1},n_{2}}\frac{{s_{1}}^{n_{1}}}{n_{1}!}\frac{{s_{2}}^{n_{2}}}{n_{2}!}\,,$
(6.40)

with $\kappa_{0,0}=0$. Since the matrices $C$ and ${\rm Id}_{1}C$ appearing in
eq. (6.39) do not commute, the general expression for $\kappa_{n_{1},n_{2}}$
is not straightforward to obtain. A simple interesting case with a non-trivial
result, which will also enable us to take its continuum limit, is
$n_{1}=1,n_{2}=1$. Using the cyclic property of trace and the symmetry of $C$,
one gets

$\displaystyle\kappa_{1,1}$ $\displaystyle=$ $\displaystyle
2\left(\Omega_{\Lambda}-\frac{\Omega_{C}}{\alpha}\right)^{2}\frac{1}{N_{1}N_{2}}{\rm
Tr}\left({\rm
Id}_{1}C^{2}\right)=2\left(\Omega_{\Lambda}-\frac{\Omega_{C}}{\alpha}\right)^{2}\frac{1}{N_{1}N_{2}}\sum_{i=1}^{N_{1}}\sum_{j=1}^{N_{2}}C_{ij}^{2}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\,$
(6.41) $\displaystyle\rightarrow$ $\displaystyle
2\left(\Omega_{\Lambda}-\frac{\Omega_{C}}{\alpha}\right)^{2}\frac{1}{V_{1}V_{2}}\int_{V_{1}}{\rm
d}^{3}\vec{x}_{1}\int_{V_{2}}{\rm
d}^{3}\vec{x}_{2}~{}C^{2}(\|\vec{x}_{1}-\vec{x}_{2}\|)\,,$

where in the last line the continuum limit involves integration of correlators
between pairs of points, with the first point in the local volume
$V_{1}=\frac{4}{3}\pi R_{1}^{3}$ and the second one in the global volume
$V_{2}=\frac{4}{3}\pi R_{2}^{3}$. Similarly to eq. (6.25), also the
integration in the last line of eq. (6.41) can be recast as a single integral
over the relative distance $r\equiv\|\vec{x}_{1}-\vec{x}_{2}\|$,

$\frac{1}{V_{1}V_{2}}\int_{V_{1}}{\rm d}^{3}\vec{x}_{1}\int_{V_{2}}{\rm
d}^{3}\vec{x}_{2}~{}C^{2}(\|\vec{x}_{1}-\vec{x}_{2}\|)=\int_{0}^{R_{2}+R_{1}}{\rm
d}r~{}p_{12}(r)C^{2}(r)\,,$ (6.42)

where now $p_{12}(r)$ is the probability distribution for the relative
distance between a point randomly chosen with uniform probability from the
sphere of radius $R_{1}$ and a point randomly chosen with uniform probability
from the sphere of radius $R_{2}>R_{1}$ with the same center. Its exact
expression is414141The strategy is the same as for $p(r)$ in eq. (6.26), but
with a slightly more complicated geometric setup. Note that when
$R_{1}=R_{2}$, then $p_{12}(r)$ reduces to the aforementioned $p(r)$. Again we
checked the correctness of the formula (6.43) numerically (for different
values of the radii ratio $R_{1}/R_{2}<1$) by randomly generating pairs of
points inside the spheres. It is also correctly normalized as a one-
dimensional distribution, namely $\int_{0}^{R_{2}+R_{1}}{\rm
d}r~{}p_{12}(r)=1$.

$p_{12}(r)=\begin{cases}\frac{3r^{2}}{R_{2}^{3}}\,,&\text{if $0\leq r\leq
R_{2}-R_{1}$}\,,\\\
\frac{3\left[r^{5}-6\left(R_{1}^{2}+R_{2}^{2}\right)r^{3}+8\left(R_{1}^{3}+R_{2}^{3}\right)r^{2}-3\left(R_{2}^{2}-R_{1}^{2}\right)^{2}r\right]}{16R_{1}^{3}R_{2}^{3}}\,,&\text{if
$R_{2}-R_{1}\leq r\leq R_{2}+R_{1}$}\,.\end{cases}$ (6.43)

Other cumulants and their continuum limit can obtained case by case from eqs.
(6.39–6.40) and then one might attempt to approximate the joint PDF,
$p\left(\left[h\right]_{1},\left[h\right]_{2}\right)$, by existing
generalizations to two variables of the Gram-Charlier series (6.29), see e.g.
[66]. Another possibility is to compute numerically the joint MGF in eq.
(6.38) and, also numerically, take its inverse Laplace transform in two
variables, again a generalization of the approach used for eq. (6.33). We did
not dig into the practical feasibility and challenges of such a program, but
we expect the numerical problem to be computationally expensive, unless a fast
routine is used. Of course, if one succeeds in evaluating
$p\left(\left[h\right]_{1},\left[h\right]_{2}\right)$, then the goal
probability ${\rm
P}\left(\left[H^{2}\right]_{V_{1}}>H_{1}^{2}\bigcap\left[H^{2}\right]_{V_{2}}<H_{2}^{2}\right)={\rm
P}\left(\left[h\right]_{1}>\frac{H_{1}^{2}}{H_{2}^{2}}\bigcap\left[h\right]_{2}<1\right)$
is given by integration of
$p\left(\left[h\right]_{1},\left[h\right]_{2}\right)$ over the proper two-
dimensional region.

##  7 Conclusions and outlook

In this paper we study the features of dark energy predicted by a simple
model, namely a light scalar field non-minimally coupled to the metric, as in
eq. (2.1). Building on previous work in [1, 29, 30, 31, 32], we see how this
model is able to produce dark energy from the amplification of quantum
fluctuations in inflation (and, less relevantly, in matter-dominated epoch).
Due to its quantum nature, dark energy emerging at later times from the
backreaction of the scalar field on the Universe expansion, is predicted to
exhibit spatial correlations of a precise form.

A first technique used to get these results starts from the observation in
[32, 29, 30] that the field backreaction is mostly made up of the energy-
momentum tensor of infrared modes. For this reason, in Section 2, we apply
Starobinsky’s stochastic formalism to study how correlators of long-wavelength
field/canonical momentum evolve from inflation until current times. The
discussion in Section 2 completes the treatment in [1] and [31] and it
presents the non-coincident 2-point and 4-point correlators at equal times.
Even starting from zero initial conditions, the accelerated expansion in
inflation amplifies quantum fluctuations and creates the spatial profile of
correlators. We see how stochastic theory predicts that, up to the overall
amplitude, the spatial shape of field correlators generated in inflation is
transmitted, roughly unaltered, throughout the subsequent epochs of cosmic
history. The scale entering the spatial profile of correlators is of the order
of the comoving Hubble horizon at the beginning of inflation, see eqs. (2.54),
(2.90). The quantum fields at two points separated by a larger comoving
distance lose any information about each other and fluctuate independently.

Interestingly, 2-point and 4-point correlators have profiles well described by
simple power laws, but unexpectedly their relation fails to reproduce Wick’s
theorem. Apart from the few most recent e-foldings when the quantum field
backreacts on the metric generating unavoidably non-Gaussianities, during all
previous stages of evolution the backreaction is negligible and the field is
Gaussian because it is well approximated by a free field in a given classical
background metric. Therefore Wick’s theorem should be respected in these
conditions. We checked numerically for possible contributions lost in
approximations when solving the system of equations for the 4-point
correlators, finding no evidence for that. For this reason, we think that the
only explanation of the problem lies in the noise sources of 4-point
functions, which would actually include some other contribution from short and
long wavelength modes interaction, which is not taken into account by eq.
(A.1). Further research is needed to understand what exactly goes wrong.

The issue with Wick’s theorem in the stochastic approach of Section 2 is one
of the main reasons for looking into another approach by which the dynamics of
infrared modes could be studied. This is done in Section 4, following the
ideas introduced by Starobinsky and Yokoyama in [56]. The strategy consists in
studying the evolution of the probability distribution of a classical field
configuration subjected to stochastic noise reproducing the coupling to short-
wavelength modes. This allows us to show, in an elegant way, that the joint
distribution for the fields at two points is Gaussian and correlators computed
from it obey Wick’s theorem, differently from the previous results of Section
2, thus rescuing stochastic formalism. We also provide a concise form for the
most general higher-order correlator between fields at two points.

Independently of its implementation, the only way to really assess whether
(and in which conditions) stochastic formalism provides a good approximation
for the evolution of the field and for dark energy generation, is to compare
its predictions to the results from quantum field theory. This is the subject
of Section 3, which presents a study for the non-coincident 2-point field
correlator and compares its full result to the stochastic theory
approximation. Even in its relative simplicity, this is to our knowledge the
first test of stochastic formalism in the non-coincident regime beyond
inflation, namely in radiation and matter epochs. It confirms that stochastic
theory gives the right amplitude of correlators and it also reproduces the
right behavior for super-Hubble separations between the two points, but at
intermediate scales the situation is more intriguing. Stochastic formalism,
focusing on long modes, is not designed to catch the shape of correlations at
sub-Hubble scales and we showed that in matter-dominated epoch, a full quantum
field theory treatment predicts that correlations persist at deeper distances
than the stochastic approximation suggests. The discrepancy is mostly
attributable to the energy content of spatial gradients, which are ignored in
stochastic formalism. We support quantitatively this argument by studying the
effect of a reduced speed of sound, which suppresses the contribution of
spatial gradients. The consequence is a better agreement between QFT and
stochastic formalism when the speed of sound is smaller.

The exploration of the observational consequences of the model is crucial to
test it using cosmological data and compare its performance to $\Lambda$CDM.
At the background level, it was found in [31] that Euclid and LSST could be
able to test the redshift-dependence of the dark energy equation of state
predicted. Encouraging results in model comparison with $\Lambda$CDM also come
from [33].

It is remarkable that the simple dark energy model considered in this paper
shows a rich phenomenology of dark energy, which is not just a time-dependent
and spatially homogeneous dark energy, but rather a spatially-correlated
quantum field, whose fluctuations can be used to compute interesting effects
on physical observables and ultimately test the model. An example in this
direction is [34], where the consequences of spatial correlations on the
luminosity distance of supernovae is quantified and compared to the signal
expected from a perturbed $\Lambda$CDM Universe (null hypothesis).
Furthermore, a reduced speed of sound would increase the signal of dark energy
fluctuations imprinted on luminosity distance correlations, to the point that
the LSST survey could detect it.

In [1] it was investigated the possibility that a dark energy of quantum
origin could relieve the Hubble tension. The result is affirmative because,
within the approximations used there, the tension goes from more than
$4\sigma$ in $\Lambda$CDM down to $1\sigma$ in the quantum dark energy model.
Independently of the details related to the model specifics, this shows that
proposals of this kind have a great potential for the Hubble tension issue. In
Section 6 of the present paper we review and further expand the theoretical
aspects of the Hubble tension within our dark energy model, highlighting the
role of spatial correlations in providing different answers for the average
Hubble rate, depending on the scale probed. The path from spatial correlations
of dark energy to the Hubble tension drawn in Section 6 can be used to extract
a few general techniques. It could serve as a basis for numerical refinements
and also be applied to other fluctuating dark energy models.

Finally, the techniques developed here can be used for a more accurate
modeling of not just dark energy, but also of dark matter, given that the
ultimate origin of these mysterious components is quantum [52, 53, 54, 55].

Acknowledgments. This work is part of the Delta ITP consortium, a program of
the Netherlands Organisation for Scientific Research (NWO) that is funded by
the Dutch Ministry of Education, Culture and Science (OCW) - NWO projectnumber
24.001.027.

##  Appendix A Sources of 4-point function stochastic equations

The stochastic sources needed in eqs. (2.78)–(2.83) are

$\displaystyle n_{\phi^{2},\phi^{2}}(t,r)$ $\displaystyle=$
$\displaystyle\frac{1}{H(t)}\Bigl{\langle}\left\\{\hat{f}_{\phi}(t,\vec{x}_{1}),\hat{\phi}(t,\vec{x}_{1})\right\\}\hat{\phi}^{2}(t,\vec{x}_{2})+\hat{\phi}^{2}(t,\vec{x}_{1})\left\\{\hat{f}_{\phi}(t,\vec{x}_{2}),\hat{\phi}(t,\vec{x}_{2})\right\\}\Bigr{\rangle}\,,$
$\displaystyle n_{\phi^{2},\phi\pi}(t,r)$ $\displaystyle=$
$\displaystyle\frac{1}{a^{3}(t)H^{2}(t)}\Bigl{\langle}\hat{\phi}^{2}(t,\vec{x}_{1})\left\\{\hat{f}_{\phi}(t,\vec{x}_{2}),\hat{\pi}(t,\vec{x}_{2})\right\\}+\left\\{\hat{f}_{\phi}(t,\vec{x}_{1}),\hat{\pi}(t,\vec{x}_{1})\right\\}\hat{\phi}^{2}(t,\vec{x}_{2})\Bigr{.}$
$\displaystyle\Bigl{.}+\hat{\phi}^{2}(t,\vec{x}_{1})\left\\{\hat{\phi}(t,\vec{x}_{2}),\hat{f}_{\pi}(t,\vec{x}_{2})\right\\}+\left\\{\hat{\phi}(t,\vec{x}_{1}),\hat{f}_{\pi}(t,\vec{x}_{1})\right\\}\hat{\phi}^{2}(t,\vec{x}_{2})\Bigr{.}$
$\displaystyle\Bigl{.}+\left\\{\hat{f}_{\phi}(t,\vec{x}_{1}),\hat{\phi}(t,\vec{x}_{1})\right\\}\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\\}\left\\{\hat{f}_{\phi}(t,\vec{x}_{2}),\hat{\phi}(t,\vec{x}_{2})\right\\}\Bigr{\rangle}\,,$
$\displaystyle n_{\phi\pi,\phi\pi}(t,r)$ $\displaystyle=$
$\displaystyle\frac{1}{a^{6}(t)H^{3}(t)}\Bigl{\langle}\left\\{\hat{\phi}(t,\vec{x}_{1}),\hat{\pi}(t,\vec{x}_{1})\right\\}\left\\{\hat{f}_{\phi}(t,\vec{x}_{2}),\hat{\pi}(t,\vec{x}_{2})\right\\}\Bigr{.}$
$\displaystyle\Bigl{.}+\left\\{\hat{f}_{\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{.}$
$\displaystyle\Bigl{.}+\left\\{\hat{\phi}(t,\vec{x}_{1}),\hat{\pi}(t,\vec{x}_{1})\right\\}\left\\{\hat{\phi}(t,\vec{x}_{2}),\hat{f}_{\pi}(t,\vec{x}_{2})\right\\}+\left\\{\hat{\phi}(t,\vec{x}_{1}),\hat{f}_{\pi}(t,\vec{x}_{1})\right\\}\left\\{\hat{\phi}(t,\vec{x}_{2}),\hat{\pi}(t,\vec{x}_{2})\right\\}\Bigr{\rangle}\,,$
$\displaystyle n_{\phi^{2},\pi^{2}}(t,r)$ $\displaystyle=$
$\displaystyle\frac{1}{a^{6}(t)H^{3}(t)}\Bigl{\langle}\left\\{\hat{f}_{\phi}(t,\vec{x}_{1}),\hat{\phi}(t,\vec{x}_{1})\right\\}\hat{\pi}^{2}(t,\vec{x}_{2})+\hat{\pi}^{2}(t,\vec{x}_{1})\left\\{\hat{f}_{\phi}(t,\vec{x}_{2}),\hat{\phi}(t,\vec{x}_{2})\right\\}\Bigr{.}$
$\displaystyle\Bigl{.}+\left\\{\hat{f}_{\pi}(t,\vec{x}_{1}),\hat{\pi}(t,\vec{x}_{1})\right\\}\hat{\phi}^{2}(t,\vec{x}_{2})+\hat{\phi}^{2}(t,\vec{x}_{1})\left\\{\hat{f}_{\pi}(t,\vec{x}_{2}),\hat{\pi}(t,\vec{x}_{2})\right\\}\Bigr{\rangle}\,,$
$\displaystyle n_{\phi\pi,\pi^{2}}(t,r)$ $\displaystyle=$
$\displaystyle\frac{1}{a^{9}(t)H^{4}(t)}\Bigl{\langle}\hat{\pi}^{2}(t,\vec{x}_{1})\left\\{\hat{f}_{\pi}(t,\vec{x}_{2}),\hat{\phi}(t,\vec{x}_{2})\right\\}+\left\\{\hat{f}_{\pi}(t,\vec{x}_{1}),\hat{\phi}(t,\vec{x}_{1})\right\\}\hat{\pi}^{2}(t,\vec{x}_{2})\Bigr{.}$
$\displaystyle\Bigl{.}+\hat{\pi}^{2}(t,\vec{x}_{1})\left\\{\hat{\pi}(t,\vec{x}_{2}),\hat{f}_{\phi}(t,\vec{x}_{2})\right\\}+\left\\{\hat{\pi}(t,\vec{x}_{1}),\hat{f}_{\phi}(t,\vec{x}_{1})\right\\}\hat{\pi}^{2}(t,\vec{x}_{2})\Bigr{.}$
$\displaystyle\Bigl{.}+\left\\{\hat{f}_{\pi}(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\\}+\left\\{\hat{\phi}(t,\vec{x}_{1}),\hat{\pi}(t,\vec{x}_{1})\right\\}\left\\{\hat{f}_{\pi}(t,\vec{x}_{2}),\hat{\pi}(t,\vec{x}_{2})\right\\}\Bigr{\rangle}\,,$
$\displaystyle n_{\pi^{2},\pi^{2}}(t,r)$ $\displaystyle=$
$\displaystyle\frac{1}{a^{12}(t)H^{5}(t)}\Bigl{\langle}\left\\{\hat{f}_{\pi}(t,\vec{x}_{1}),\hat{\pi}(t,\vec{x}_{1})\right\\}\hat{\pi}^{2}(t,\vec{x}_{2})+\hat{\pi}^{2}(t,\vec{x}_{1})\left\\{\hat{f}_{\pi}(t,\vec{x}_{2}),\hat{\pi}(t,\vec{x}_{2})\right\\}\Bigr{\rangle}\,.$
(A.1)

The stochastic sources can be computed in terms of the mode function
$\varphi(t,k)$. The results depend on integrals over Fourier modes, which we
regulate at low momenta with an IR cutoff $k_{0}$ (while $\mu aH$ plays the
role of an UV cutoff). Their expression is

$\displaystyle n_{\phi^{2},\phi^{2}}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi^{4}}(\mu
aH)^{3}(1-\epsilon)\left[|\varphi(t,k)|^{2}\right]_{k=\mu aH}\int_{k_{0}}^{\mu
aH}dk~{}k^{2}|\varphi(t,k)|^{2}\left[1+2j_{0}(\mu aHr)j_{0}(kr)\right]\,,$
$\displaystyle n_{\phi^{2},\phi\pi}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi^{4}}\mu^{3}a^{3}H^{2}(1-\epsilon)\left\\{\left[\frac{\partial}{\partial
t}|\varphi(t,k)|^{2}\right]_{k=\mu aH}\int_{k_{0}}^{\mu
aH}dk~{}k^{2}|\varphi(t,k)|^{2}\left[1+2j_{0}(\mu aHr)j_{0}(kr)\right]\right.$
$\displaystyle\left.+\left[|\varphi(t,k)|^{2}\right]_{k=\mu
aH}\int_{k_{0}}^{\mu aH}dk~{}k^{2}\frac{\partial}{\partial
t}|\varphi(t,k)|^{2}\left[1+2j_{0}(\mu aHr)j_{0}(kr)\right]\right\\}\,,$
$\displaystyle n_{\phi\pi,\phi\pi}$ $\displaystyle=$
$\displaystyle\frac{1}{\pi^{4}}\mu^{3}a^{3}H(1-\epsilon)\left\\{\frac{1}{2}\left[\frac{\partial}{\partial
t}|\varphi(t,k)|^{2}\right]_{k=\mu aH}\int_{k_{0}}^{\mu
aH}dk~{}k^{2}\frac{\partial}{\partial t}|\varphi(t,k)|^{2}\left[1+j_{0}(\mu
aHr)j_{0}(kr)\right]\right.$ $\displaystyle\left.+j_{0}(\mu
aHr)\left[|\varphi(t,k)|^{2}\right]_{k=\mu aH}\int_{k_{0}}^{\mu
aH}dk~{}k^{2}|\dot{\varphi}(t,k)|^{2}\right.$ $\displaystyle\left.+j_{0}(\mu
aHr)\left[|\dot{\varphi}(t,k)|^{2}\right]_{k=\mu aH}\int_{k_{0}}^{\mu
aH}dk~{}k^{2}|\varphi(t,k)|^{2}+\frac{1}{2}\left(\frac{\mu
H}{a}\right)^{3}j_{0}(\mu aHr)\frac{j_{1}(\mu aHr)}{\mu aHr}\right\\}\,,$
$\displaystyle n_{\phi^{2},\pi^{2}}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi^{4}}\mu^{3}a^{3}H(1-\epsilon)\left\\{\left[|\varphi(t,k)|^{2}\right]_{k=\mu
aH}\int_{k_{0}}^{\mu aH}dk~{}k^{2}|\dot{\varphi}(t,k)|^{2}\right.$
$\displaystyle\left.+\left[|\dot{\varphi}(t,k)|^{2}\right]_{k=\mu
aH}\int_{k_{0}}^{\mu aH}dk~{}k^{2}|\varphi(t,k)|^{2}\right.$
$\displaystyle\left.+j_{0}(\mu aHr)\left[\frac{\partial}{\partial
t}|\varphi(t,k)|^{2}\right]_{k=\mu aH}\int_{k_{0}}^{\mu
aH}dk~{}k^{2}\frac{\partial}{\partial t}|\varphi(t,k)|^{2}-\left(\frac{\mu
H}{a}\right)^{3}j_{0}(\mu aHr)\frac{j_{1}(\mu aHr)}{\mu aHr}\right\\}\,,$
$\displaystyle n_{\phi\pi,\pi^{2}}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi^{4}}\mu^{3}a^{3}(1-\epsilon)\left\\{\left[\frac{\partial}{\partial
t}|\varphi(t,k)|^{2}\right]_{k=\mu aH}\int_{k_{0}}^{\mu
aH}dk~{}k^{2}|\dot{\varphi}(t,k)|^{2}\left[1+2j_{0}(\mu
aHr)j_{0}(kr)\right]\right.$
$\displaystyle\left.+\left[|\dot{\varphi}(t,k)|^{2}\right]_{k=\mu
aH}\int_{k_{0}}^{\mu aH}dk~{}k^{2}\frac{\partial}{\partial
t}|\varphi(t,k)|^{2}\left[1+2j_{0}(\mu aHr)j_{0}(kr)\right]\right\\}\,,$
$\displaystyle n_{\pi^{2},\pi^{2}}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi^{4}}\frac{\mu^{3}a^{3}}{H}(1-\epsilon)\left[|\dot{\varphi}(t,k)|^{2}\right]_{k=\mu
aH}\int_{k_{0}}^{\mu aH}dk~{}k^{2}|\dot{\varphi}(t,k)|^{2}\left[1+2j_{0}(\mu
aHr)j_{0}(kr)\right]\,.$ (A.2)

where $j_{0}(z)\equiv\frac{\sin z}{z}$ and
$j_{1}(z)\equiv-\frac{d}{dz}j_{0}(z)$ are the spherical Bessel functions of
order $0$ and $1$ respectively, which are evaluated at $z=kr$.

In de Sitter inflation with constant Hubble parameter $H_{I}$, assuming
$\xi<0$ and $(m/H_{I})^{2}\ll|\xi|\ll 1$, their leading behavior is given by
eq. (2.88).

##  Appendix B Correlators in matter + cosmological constant epoch

This Appendix contains the evolution of the scalar field correlators after
matter-radiation equality, including both matter and a cosmological constant
(CC) as the fluids leading the expansion. This could serve as a possible
refinement of the evolution in a pure matter-dominated Universe given in
subsections 2.2.3 and 2.3.3, because the energy-momentum tensor of the scalar
field contains a part behaving like a cosmological constant that contributes
to the expansion. For simplicity, we focus on 2-pt functions, but it is
straightforward to extend the treatment to 4-pt functions. Part of the work
was already done in Appendix A.3 of [1]. Our goal is to solve the system of
equations (2.31–2.33) in a matter+CC Universe to predict the 2-pt correlators
today.

We set the number of e-foldings $N$ to $N=0$ at matter-radiation equality and
denote by $N_{0}=\ln{\left(\frac{\Omega_{M}}{\Omega_{R}}\right)}\simeq 8.1$
the current time (today). As usual, $\Omega_{R}\simeq 9.1\times 10^{-5}$,
$\Omega_{M}\simeq 0.3$ and $\Omega_{\Lambda}\simeq 1-\Omega_{M}\simeq 0.7$ are
the radiation, matter and cosmological constant fractions of energy density
today. The parameter $\epsilon=-\dot{H}/H^{2}$ as a function of $N$ is

$\epsilon(N)=\frac{3}{2}\frac{1}{1+\frac{\Omega_{\Lambda}}{\Omega_{M}}e^{3(N-N_{0})}}\,.$
(B.1)

Following Appendix A.3 of [1], let us change variable from $N$ to

$n(N)\equiv\frac{1}{1+\frac{\Omega_{\Lambda}}{\Omega_{M}}e^{3(N-N_{0})}}\,.$
(B.2)

Then $\epsilon(n)=(3/2)n$ and the ratio between the effective mass and the
Hubble rate $M(n)/H(n)$ is

$\left(\frac{M(n)}{H(n)}\right)^{2}=\left(\frac{m}{H_{\rm
DE}}\right)^{2}(1-n)-12|\xi|\left(1-\frac{3}{4}n\right)\,,$ (B.3)

where $H_{\rm DE}=H_{0}\sqrt{\Omega_{\Lambda}}$. Matter-radiation equality
$N=0$ corresponds to $n_{\rm
eq}\equiv\left(1+\frac{\Omega_{\Lambda}}{\Omega_{M}}e^{-3N_{0}}\right)^{-1}$
while the current time $N_{0}$ gives $n_{0}\equiv\Omega_{M}$. The derivatives
with respect to $N$ are easily transformed to derivatives with respect to $n$
by using $dn/dN=3n(n-1)$. Neglecting stochastic noise sources, as legitimate
after inflation, eqs. (2.31–2.33) are equivalent, in vector/matrix notation,
to

$\frac{{\mathrm{d}}}{{\mathrm{d}}n}\Delta_{(\rm 2)}(n,r)+A_{(\rm
2)}(n)~{}\Delta_{(\rm 2)}(n,r)=0\,,$ (B.4)

where the matrix $A_{(\rm 2)}(n)$ can be decomposed as $A_{(\rm
2)}(n)=b_{1}(n)B_{1}+b_{2}(n)B_{2}+b_{3}(n)B_{3}$, with

$b_{1}(n)\equiv\frac{2-n}{2n(n-1)}\,,\qquad
b_{2}(n)\equiv-\frac{\left(\frac{m}{H_{\rm
DE}}\right)^{2}}{3n}-|\xi|\frac{4-3n}{n(n-1)}\,,\qquad
b_{3}(n)\equiv\frac{1}{3n(n-1)}\,,$ (B.5)

and

$B_{1}=\begin{pmatrix}0&0&0\\\ 0&1&0\\\ 0&0&2\\\ \end{pmatrix}\,,\qquad
B_{2}=\begin{pmatrix}0&0&0\\\ 2&0&0\\\ 0&1&0\\\ \end{pmatrix}\,,\qquad
B_{3}=\begin{pmatrix}0&-1&0\\\ 0&0&-2\\\ 0&0&0\\\ \end{pmatrix}\,.$ (B.6)

The exact solution of (B.4) is

$\Delta_{(\rm 2)}(n,r)=T\exp\left[-\int_{n_{\rm eq}}^{n}dn^{\prime}~{}A_{(\rm
2)}(n^{\prime})\right]\Delta_{(\rm 2)}(n_{\rm eq},r)\,,$ (B.7)

where $T\exp$ denotes the time-ordered exponential, accounting for the non-
commutativity of $B(n)$ matrices for different $n$.

As we will see at the end of this Appendix, in particular in Fig. 7,
neglecting the time ordering only induces errors at a few percent level.
Therefore, before embarking in a refined estimate using another method, let us
mention what the result is when time-ordering is not considered. In this case,
as already done for 4-pt functions in Appendix A.3 of [1], one finds that at
leading order in $|\xi|$ and $\left(\frac{m}{H_{\rm DE}}\right)^{2}$, the 2-pt
correlators today have a relatively simple approximate expression424242This
result refers strictly to the current time
$N_{0}=\ln{\left(\frac{\Omega_{M}}{\Omega_{R}}\right)}$. At intermediate times
$0<N<N_{0}$ the correlators have more complicated expressions.

$\begin{pmatrix}\Delta_{\phi,\phi}(N_{0},r)\\\ \Delta_{\phi,\pi}(N_{0},r)\\\
\Delta_{\pi,\pi}(N_{0},r)\\\ \end{pmatrix}\simeq\begin{pmatrix}1\\\ 4\zeta\\\
4\zeta^{2}\\\ \end{pmatrix}e^{4\zeta N_{0}}~{}\Delta_{\phi,\phi}(0,r)\,,$
(B.8)

where $\Delta_{\phi,\phi}(0,r)$ is the only 2-pt correlators inherited from
radiation epoch (see eq. (2.61)) and $\zeta$ is a combination of parameters
(the same as in eq. (15) of [1]) defined as

$\zeta\equiv|\xi|\frac{\frac{3}{2}\frac{\ln\left(\Omega_{R}\right)}{\ln\left(\Omega_{M}\right)}-\frac{1}{6|\xi|}\left(\frac{m}{H_{\rm
DE}}\right)^{2}}{\frac{3}{2}\frac{\ln\left(\Omega_{R}\right)}{\ln\left(\Omega_{M}\right)}-1}\,.$
(B.9)

We now propose another strategy to study the evolution of correlators, which
circumvents the computation of time-ordered exponentials. We checked
numerically that it reproduces the full time-ordered solution (B.7). Finally,
Fig. 7 will show that neglecting time-ordering in the exponential (B.7)
actually gives a surprisingly accurate result.

###  B.1 Alternative solution

The full quantum field $\hat{\Phi}$ obeys the second order equation of motion
(2.8). Trading cosmological time $t$ for the number of e-folding $N$ and
omitting for simplicity of notation the spacetime dependence of $\hat{\Phi}$,

$\partial_{N}^{2}\hat{\Phi}+\left(3-\epsilon(N)\right)\partial_{N}\hat{\Phi}-\frac{\nabla^{2}}{a^{2}(N)H^{2}(N)}\hat{\Phi}+\frac{M^{2}(N)}{H^{2}(N)}\hat{\Phi}=0\,.$
(B.10)

As we already did before, for super-Hubble (infrared) modes $k\ll aH$ we can
neglect spatial gradients, which here means dropping the Laplacian term. In
matter+CC epoch we switch from $N$ to the variable $n=\frac{2}{3}\epsilon(n)$
as in eqs. (B.1)–(B.2), and then eq. (B.10) becomes

$\left[n^{2}(1-n)^{2}\partial_{n}^{2}-\frac{3}{2}n^{2}(1-n)\partial_{n}+\left(\frac{m}{3H_{DE}}\right)^{2}(1-n)-\frac{|\xi|}{3}(4-3n)\right]\hat{\Phi}=0\,.$
(B.11)

Notice that, when $N$ increases, $n$ decreases, therefore the goal is to solve
the equation (B.11) from the initial $n_{\rm
eq}=\left(1+\frac{\Omega_{\Lambda}}{\Omega_{M}}e^{-3N_{0}}\right)^{-1}$ close
to 1 (but slightly smaller), down to smaller $n$ (our current time
corresponding to $n_{0}=\Omega_{M}$ and the infinitely far future being $n\to
0$).

One can show that, with a suitable (time-dependent) rescaling of the field
$\hat{\Phi}$, eq. (B.11) reduces to a hypergeometric equation. The procedure
consists in introducing a new variable $\hat{\tilde{\Phi}}$ related to the old
$\hat{\Phi}$ by,

$\hat{\Phi}=n^{\alpha}(1-n)^{\beta}\hat{\tilde{\Phi}}\,,$ (B.12)

where $\alpha$ and $\beta$ are constants which are chosen in such a way that
$\hat{\tilde{\Phi}}$ satisfies a hypergeometric equation. One can then check
that, taking $\alpha$ and $\beta$ as solutions of the quadratic equations,

$\displaystyle\alpha^{2}-\alpha+\left(\frac{m}{3H_{\rm
DE}}\right)^{2}-\frac{4}{3}|\xi|$ $\displaystyle=$ $\displaystyle 0\,,$
$\displaystyle\beta^{2}+\frac{1}{2}\beta-\frac{1}{3}|\xi|$ $\displaystyle=$
$\displaystyle 0\,,$ (B.13)

then (B.11) simplifies to

$n(1-n)\partial_{n}^{2}\hat{\tilde{\Phi}}+\left[\alpha-\left(\alpha+\beta+\frac{3}{2}\right)n\right]\partial_{n}\hat{\tilde{\Phi}}-(\alpha+\beta)\left(\alpha+\beta+\frac{1}{2}\right)\hat{\tilde{\Phi}}=0\,,$

which is a hypergeometric equation. Note that, under the assumption of light
field $m/H_{\rm DE}<1$ that we are using in this work, we are guaranteed that
the solutions $\alpha$ and $\beta$ of (B.1) are real numbers. However we still
have to make a choice for the solutions since each of those quadratic
equations has two solutions. This will be done in a moment after requiring
identification of (B.1) with the standard form of hypergeometric equation
(with variable $n$)

$n(1-n)\partial_{n}^{2}\hat{\tilde{\Phi}}+\left[c-\left(a+b+1\right)n\right]\partial_{n}\hat{\tilde{\Phi}}-ab~{}\hat{\tilde{\Phi}}=0\,,$
(B.14)

so that $c=\alpha$, $a+b+1=\alpha+\beta+\frac{3}{2}$ and
$ab=(\alpha+\beta)\left(\alpha+\beta+\frac{1}{2}\right)$. The solution (up to
the irrelevant symmetry between $a$ and $b$ built in (B.14) itself) is,
without loss of generality,

$\displaystyle a$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[\frac{1}{2}+\alpha+\beta+\sqrt{\left(\frac{1}{2}+\alpha+\beta\right)\left(\frac{1}{2}-3\alpha-3\beta\right)}\right]\,,$
$\displaystyle b$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[\frac{1}{2}+\alpha+\beta-\sqrt{\left(\frac{1}{2}+\alpha+\beta\right)\left(\frac{1}{2}-3\alpha-3\beta\right)}\right]\,,$
$\displaystyle c$ $\displaystyle=$ $\displaystyle\alpha\,.$ (B.15)

Note that $a$ and $b$ are real if the quantity under square root appearing in
(B.1) is non-negative, i.e. when
$-\frac{1}{2}\leq\alpha+\beta\leq\frac{1}{6}$. This gives a way to select the
solutions for $\alpha$ and $\beta$ in (B.1). For our purposes, it is
sufficient to work at linear order in the small quantities $(m/H_{\rm
DE})^{2}$ and $|\xi|$. It is then straightforward to check that the only
choice which complies with the condition
$-\frac{1}{2}\leq\alpha+\beta\leq\frac{1}{6}$ for small but otherwise
arbitrary values of $(m/H_{\rm DE})^{2}$ and $|\xi|$ is

$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(1-\sqrt{1-\left(\frac{2m}{3H_{\rm
DE}}\right)^{2}+\frac{16}{3}|\xi|}\right)\,,$ $\displaystyle\beta$
$\displaystyle=$
$\displaystyle\frac{1}{4}\left(-1+\sqrt{1+\frac{16}{3}|\xi|}\right)\,.$ (B.16)

Due to the singular points at $n=0$, $n=1$ and $n\to\infty$ of the
hypergeometric differential equation (B.14), the choice of two linearly
independent solutions for it depends on the interval where one has to solve
it. For our problem, as already anticipated, the interval of $n$ values needed
is
$\left[\Omega_{M},\left(1+\frac{\Omega_{\Lambda}}{\Omega_{M}}e^{-3N_{0}}\right)^{-1}\right]$
(if we want to evolve from matter-radiation equality until today) or
$\left(0,\left(1+\frac{\Omega_{\Lambda}}{\Omega_{M}}e^{-3N_{0}}\right)^{-1}\right]$
if we also want to extrapolate our predictions to the future (of course the
singular point $n=0$ is never reached, since it corresponds to an infinite
scale factor). In these intervals (and for non-integer values of $c-a-b$) two
independent solutions are known to be
$~{}_{2}F_{1}\left(a,b;1+a+b-c;1-n\right)$ and
$(1-n)^{c-a-b}~{}_{2}F_{1}\left(c-a,c-b;1+c-a-b;1-n\right)$, where
$~{}_{2}F_{1}$ is the Gauss’ hypergeometric function. We can simplify some of
the coefficients by using (B.1) and then write the solution of (B.14) as

$\hat{\tilde{\Phi}}(n,\vec{x})=\hat{A}(\vec{x})~{}_{2}F_{1}\left(a,b;\frac{3}{2}+\beta;1-n\right)+\hat{B}(\vec{x})(1-n)^{-\frac{1}{2}-\beta}~{}_{2}F_{1}\left(\alpha-a,\alpha-b;\frac{1}{2}-\beta;1-n\right)\,.$

In the equation above we wrote explicitly the dependence of the field on the
spatial (comoving) position $\vec{x}$, which is reflected on the right-hand
side by the quantum operators $\hat{A}(\vec{x})$ and $\hat{B}(\vec{x})$ (time-
independent). Thus, taking aso into account the factors in (B.12) and the
relation (B.1), we find

$\displaystyle\hat{\Phi}(n,\vec{x})=F_{A}(n)\hat{A}(\vec{x})+F_{B}(n)\hat{B}(\vec{x})\,,$
(B.17)

where we introduced a compact notation for the time-dependent functions

$\displaystyle F_{A}(n)$ $\displaystyle=$ $\displaystyle
n^{\alpha}(1-n)^{\beta}~{}_{2}F_{1}\left(a,b;\frac{3}{2}+\beta;1-n\right)\,,$
$\displaystyle F_{B}(n)$ $\displaystyle=$ $\displaystyle
n^{\alpha}(1-n)^{-\frac{1}{2}}~{}_{2}F_{1}\left(\alpha-a,\alpha-b;\frac{1}{2}-\beta;1-n\right)\,,$
(B.18)

with $a$ and $b$ determined by $\alpha$ and $\beta$ according to eq. (B.1). We
can also write the corresponding solution for the canonical momentum operator
$\hat{\Pi}(n,\vec{x})$ by looking at eq. (2.6) and transforming from the
cosmological time $t$ variable to the $x$ variable to get for the rescaled
momentum $\hat{\Pi}/(a^{3}H)$ relevant for our correlators

$\frac{\hat{\Pi}(n,\vec{x})}{a^{3}(n)H(n)}=-3n(1-n)~{}\partial_{n}\hat{\Phi}(n,\vec{x})\,.$
(B.19)

Upon inserting (B.17) into (B.19) and after taking derivative of the
hypergeometric functions, we obtain

$\frac{\hat{\Pi}(n,\vec{x})}{a^{3}(n)H(n)}=G_{A}(n)\hat{A}(\vec{x})+G_{B}(n)\hat{B}(\vec{x})\,,$
(B.20)

where the functions $G_{A}(n)$ and $G_{B}(n)$ are defined as

$\displaystyle G_{A}(n)$ $\displaystyle=$
$\displaystyle-3n^{\alpha+1}(1-n)^{\beta+1}\left[\left(\frac{\alpha}{n}-\frac{\beta}{1-n}\right)~{}_{2}F_{1}\left(a,b;\frac{3}{2}+\beta;1-n\right)\right.\,$
$\displaystyle\left.-\frac{(\alpha+\beta)(\alpha+\beta+\frac{1}{2})}{\frac{3}{2}+\beta}~{}_{2}F_{1}\left(a+1,b+1;\frac{5}{2}+\beta;1-n\right)\right]\,,$
$\displaystyle G_{B}(n)$ $\displaystyle=$
$\displaystyle-3n^{\alpha+1}(1-n)^{\frac{1}{2}}\left[\left(\frac{\alpha}{n}+\frac{1/2}{1-n}\right)~{}_{2}F_{1}\left(\alpha-a,\alpha-b;\frac{1}{2}-\beta;1-n\right)\right.\,$
(B.21)
$\displaystyle\left.-\frac{\alpha^{2}+\beta(\alpha+\beta+\frac{1}{2})}{\frac{1}{2}-\beta}~{}_{2}F_{1}\left(\alpha-a+1,\alpha-b+1;\frac{3}{2}-\beta;1-n\right)\right]\,.$

We want to use these results to evaluate the coincident 2-pt functions
$\Delta_{\phi\phi}$, $\Delta_{\phi\pi}$ and $\Delta_{\pi\pi}$ defined in
(2.25–2.27). We begin by observing that eqs. (B.17) and (B.20) imply

$\begin{pmatrix}\Delta_{\phi\phi}(n)\\\ \Delta_{\phi\pi}(n)\\\
\Delta_{\pi\pi}(n)\\\
\end{pmatrix}=\mathcal{M}(n)\begin{pmatrix}\langle\hat{A}^{2}\rangle\\\
\langle\\{\hat{A},\hat{B}\\}\rangle\\\ \langle\hat{B}^{2}\rangle\\\
\end{pmatrix}\,,$ (B.22)

where the matrix $\mathcal{M}(n)$ is defined as

$\mathcal{M}(n)=\begin{pmatrix}F_{A}^{2}(n)&F_{A}(n)F_{B}(n)&F_{B}^{2}(n)\\\
2F_{A}(n)G_{A}(n)&F_{A}(n)G_{B}(n)+G_{A}(n)F_{B}(n)&2F_{B}(n)G_{B}(n)\\\
G_{A}^{2}(n)&G_{A}(n)G_{B}(n)&G_{B}^{2}(n)\\\ \end{pmatrix}\,.$ (B.23)

The constant correlators $\langle\hat{A}^{2}\rangle$,
$\langle\\{\hat{A},\hat{B}\\}\rangle$ and $\langle\hat{B}^{2}\rangle$ (which
are also space-independent due to the homogeneity of the FLRW background) are
determined by the values of $\Delta_{\phi\phi}(n_{\rm eq})$,
$\Delta_{\phi\pi}(n_{\rm eq})$ and $\Delta_{\pi\pi}(n_{\rm eq})$ at matter-
radiation equality inherited from radiation epoch, so that we conclude

$\begin{pmatrix}\Delta_{\phi\phi}(n)\\\ \Delta_{\phi\pi}(n)\\\
\Delta_{\pi\pi}(n)\\\ \end{pmatrix}=\mathcal{M}(n)\mathcal{M}^{-1}(n_{\rm
eq})\begin{pmatrix}\Delta_{\phi\phi}(n_{\rm eq})\\\ \Delta_{\phi\pi}(n_{\rm
eq})\\\ \Delta_{\pi\pi}(n_{\rm eq})\\\ \end{pmatrix}\,.$ (B.24)

The same equation also holds for non-coincident correlators as one can see by
considering fields at two comoving positions $\vec{x}$ and $\vec{y}$. As we
know, homogeneity and isotropy of the FLRW background imply that the spatial
dependence of these correlators only appears through the relative distance
$r=\|\vec{x}-\vec{y}\|$. Then from eqs. (B.17) and (B.20) it follows that

$\begin{pmatrix}\Delta_{\phi\phi}(n,r)\\\ \Delta_{\phi\pi}(n,r)\\\
\Delta_{\pi\pi}(n,r)\\\
\end{pmatrix}=\mathcal{M}(n)\begin{pmatrix}\langle\hat{A}(\vec{x})\hat{A}(\vec{y})\rangle\\\
\langle\hat{A}(\vec{x})\hat{B}(\vec{y})+\hat{B}(\vec{x})\hat{A}(\vec{y})\rangle\\\
\langle\hat{B}(\vec{x})\hat{B}(\vec{y})\rangle\\\ \end{pmatrix}\,,$ (B.25)

where $\mathcal{M}(n)$ is the same matrix defined in (B.23). In terms of
initial conditions at matter-radiation equality

$\begin{pmatrix}\Delta_{\phi\phi}(n,r)\\\ \Delta_{\phi\pi}(n,r)\\\
\Delta_{\pi\pi}(n,r)\\\ \end{pmatrix}=\mathcal{M}(n)\mathcal{M}^{-1}(n_{\rm
eq})\begin{pmatrix}\Delta_{\phi\phi}(n_{\rm eq},r)\\\ \Delta_{\phi\pi}(n_{\rm
eq},r)\\\ \Delta_{\pi\pi}(n_{\rm eq},r)\\\ \end{pmatrix}\,,$ (B.26)

which shows the same time evolution as in (B.24). Fig. 7 shows the comparison
between the full result for the coincident correlators (B.24) computed using
the technique explained here against the approximation of neglecting time-
ordering in the matrix exponentials of (B.7). The values of the cosmological
parameters chosen are given in the caption of the figure. For the most
relevant correlator $\Delta_{\phi\phi}(N)$ the difference between the two
methods at the current epoch ($N_{0}\simeq 8.1$ e-folds since matter-radiation
equality) is less than $2\%$ level. Thus the approximation without time-
ordering works quite well, at the few percent level. The deviation grows to
$6\%$ in the future at $N=12$.

Figure 7: The coincident 2-point functions $\Delta_{\phi\phi}$,
$\Delta_{\phi\pi}$ and $\Delta_{\pi\pi}$ in the cosmological era dominated by
non-relativistic matter and cosmological constant as functions of the number
of e-foldings starting from matter-radiation equality. The full lines refer to
the results in terms of hypergeometric functions given in (B.24), while the
dashed curves are computed by neglecting time-ordering of matrix exponentials
in (B.7). The black vertical line marks the current epoch ($N_{0}\simeq 8.1$).
The curves are normalized with respect to the initial value of
$\Delta_{\phi\phi}$ at matter-radiation equality. The correlators
$\Delta_{\phi\pi}$ and $\Delta_{\pi\pi}$ have utterly negligible initial
values because of their exponential suppression through radiation epoch. The
values used for the cosmological density fractions today are
$\Omega_{R}=9.1\times 10^{-5}$, $\Omega_{M}=0.3$, $\Omega_{\Lambda}\simeq
1-\Omega_{M}\simeq 0.7$, while for the model we adopt $\xi=-0.06$ and
$m/H_{\rm DE}=0.6$. For the largest correlator $\Delta_{\phi\phi}$, the
relative difference between the full result and the approximated is less than
$2\%$ today and $6\%$ in the future when $N=12$. The difference is even
smaller for the less relevant correlators $\Delta_{\phi\pi}$ and
$\Delta_{\pi\pi}$.

##  Appendix C Effect of non-zero initial conditions in inflation

Here we investigate the consequences of initial conditions for the
inflationary era inherited from a pre-inflationary epoch. If $P_{\phi}(t_{\rm
in},k)$ is the power spectrum of the field at the beginning of inflation
inherited from the pre-inflationary epoch (where $t_{\rm in}$ is cosmological
time at the beginning of inflation and $k$ is a comoving wavenumber), then its
2-point function is

$\Delta_{\phi\phi}(t_{\rm in},r)=\frac{1}{2\pi^{2}}\int_{k_{0}}^{\mu a_{\rm
in}H_{I}}\frac{dk}{k}j_{0}(kr)P_{\phi}(t,k)\,.$ (C.1)

Let us suppose that the field power spectrum follows a distribution of the
form

$P_{\phi}(t_{\rm in},k)=\left(\frac{k}{a_{\rm
in}}\right)^{2}\left[n_{0}\left(\frac{a_{\rm in}T_{\rm
in}}{k}\right)^{\sigma}+\frac{1}{2}\right]\,.$ (C.2)

The equation above is obtained naturally by supposing that the average
occupation number for a mode with momentum $k$ follows a thermal-like
distribution $\bar{n}_{k}=n_{0}\left(\frac{a_{\rm in}T_{\rm
in}}{k}\right)^{\sigma}$, where we allow for a generic power $\sigma$ and
$T_{\rm in}$ is a generic mass scale parameter that plays the role of a
temperature (but it is not necessarily the physical temperature).

The term $1/2$ in eq. (C.2) is the zero point energy. Notice that, apart from
this $1/2$ addend, the case $\sigma=2$ corresponds to a scale-invariant
spectrum.

Using (C.2) into (C.1) and approximating the Bessel function as a Heaviside
$\theta$ function, we get for $\sigma<2$ that

$\Delta_{\phi,\phi}(t_{\rm in},r)=\frac{(\mu
H_{I})^{2}}{2\pi^{2}}d^{2}(r)\left[\frac{n_{0}}{2-\sigma}\left(\frac{T_{\rm
in}}{\mu H_{I}d(r)}\right)^{\sigma}+\frac{1}{4}\right]\,.$ (C.3)

The spatial dependence has been encoded in the function $d(r)$ given by

$d(r)\simeq\begin{cases}1&\text{if $r<\frac{1}{\mu a_{\rm in}H_{I}}$}\\\
\frac{1}{\mu a_{\rm in}H_{I}r}&\text{if $\frac{1}{\mu a_{\rm
in}H_{I}}<r\lesssim\frac{1}{k_{0}}$}\,,\end{cases}$ (C.4)

where we have truncated distances up to $1/k_{0}$ since $k_{0}$ is the IR
cutoff. We want to follow the evolution of 2-pt correlators during inflation.

As for the initial conditions of $\Delta_{\phi,\pi}$ and $\Delta_{\pi,\pi}$ we
suppose $\Delta_{\phi,\pi}(t_{\rm in},r)=0$ and $\Delta_{\pi,\pi}(t_{\rm
in},r)=R\Delta_{\phi,\phi}(t_{\rm in},r)$, where we introduced a parameter $R$
as the ratio between the amplitude of the initial power spectrum of the
canonical momentum $\hat{\pi}$, compared to that of the field $\hat{\phi}$.

Evolving these initial conditions, at the end of a phase of de Sitter
inflation, the (coincident for simplicity) 2-point functions $\Delta_{(\rm
2)}=\left(\Delta_{\phi,\phi},\Delta_{\phi,\pi},\Delta_{\pi,\pi}\right)$ are in
row vector form

$\Delta_{(\rm 2)}^{\rm T}(N_{I})=\left(1+\frac{R}{9}\right)\frac{(\mu
H_{I})^{2}}{2\pi^{2}}e^{8|\xi|N_{I}}\left[n_{0}\left(\frac{T_{in}}{\mu
H_{I}}\right)^{\sigma}\left(\frac{1}{2-\sigma}+\frac{1}{\tau+8|\xi|}\right)+\frac{1}{4}+\frac{1}{16|\xi|}\right]\left(1,8|\xi|,16\xi^{2}\right)\,,$
(C.5)

where the parameter $\tau$ includes both the possibilities of the