\section{Introduction}

Cosmic structures in the universe are understood to be seeded by some
pre-existing super-Hubble cosmological fluctuations. Their
gravitational collapse starts when their size becomes smaller than the Hubble radius,
an inevitable outcome in any decelerating Friedmann-Lema\^itre
spacetime. Observational evidence of this mechanism is present in
the Cosmic Microwave Background (CMB) data by the correlation patterns
associated with the polarization and temperature angular power
spectra~\cite{Hu:1996yt, Planck:2018nkj}, as well as in the statistics
of the large-scale structures observed at lower redshifts~\cite{Fonseca:2015laa,
  EUCLID:2020jwq}.

Cosmic Inflation, an early era of accelerated cosmic expansion, is the
prime candidate to explain the origin of the super-Hubble
fluctuations. They are of quantum origin, stretched to length scales
much larger than the Hubble radius during
inflation~\cite{Starobinsky:1979ty, Starobinsky:1980te, Guth:1980zm,
  Linde:1981mu, Albrecht:1982wi, Linde:1983gd, Mukhanov:1981xt,
  Mukhanov:1982nu, Starobinsky:1982ee, Guth:1982ec, Hawking:1982cz,
  Bardeen:1983qw}. At the same time, inflation smooths out any
pre-existing inhomogeneity and one of the historical motivations for
Cosmic Inflation is that the spatial curvature of spacetime,
$\Omega_{\uK}$, should be exponentially small at the end of inflation (at
most $e^{-60}$). This is required to account for the current bound
$|\Omega_{\uK_{0}}|<3 \times 10^{-3}$ today, coming from the Planck CMB
data and Baryon Acoustic Oscillations (BAO) measurements.

Intuitively, the existence, today, of Hubble-sized curvature
fluctuations suggests that these could be confused with a small
non-vanishing spatial curvature of the local background metric. In
particular, these modes are expected to induce a limitation on our
ability to measure very small values of the curvature density
parameter~\cite{Waterhouse:2008vb, Buchert:2008zm, Vardanyan:2009ft,
  Leonard:2016evk}. More than being a nuisance, we will show that
super-Hubble (hence ``conserved'') fluctuations do create spatial
curvature.
  

In order to deal with fluctuations over a background metric when both are
intertwined, we can start from the inhomogeneous metric proposed in
Refs.~\cite{Salopek:1990jq, Creminelli:2004yq, Kolb:2004jg, Lyth:2004gb}
\begin{equation}
  \dd s^2 = -\dd \tau^2 + a^2(\tau) e^{2 \zeta(\tau,\vb*{x})}
    \delta_{ij}\dd x^i \dd x^j.
\label{eq:metric}
\end{equation}
This metric is not fully general as inhomogeneities are all contained
in one scalar function $\zeta$. However, as discussed in
Refs.~\cite{Salopek:1990jq, Creminelli:2004yq, Kolb:2004jg,
  Lyth:2004gb}, this is the most generic metric in absence of vector-
and tensor-type inhomogeneities, and, in the gauge where fixed time
slices have uniform energy density and fixed spatial worldlines are
comoving with matter. At super-Hubble scales this reduces to
the synchronous gauge supplemented by some additional conditions that
fix it uniquely. The quantity $\zeta(\tau,\vb*{x})$ can be shown to
be ``conserved'' at large distances. As such, it provides a non-linear
generalization of the constant-energy-density curvature
perturbation~\cite{Rigopoulos:2004gr, Langlois:2005ii}.

Historically, this metric has been intensively discussed in the
attempts to explain the acceleration of the Universe by the
backreaction of super-Hubble
inhomogeneities~\cite{Kolb:2005me,Barausse:2005nf}. But, as realized
soon after~\cite{Hirata:2005ei, Kolb:2005da, Geshnizjani:2005ce}, the
sole observable effects of super-Hubble fluctuations are to modify the
spatial curvature. Up to our knowledge, the only works having
addressed how super-Hubble modes affect the spatial curvature are
Refs.~\cite{Geshnizjani:2005ce, Kleban:2012ph}, based, however, on
perturbative gradient expansions or linear perturbation theory
only. When the non-perturbative terms of our derivation can be
neglected, we recover some of their results.

The paper is organized as follows. In \cref{sec:curvdef}, we derive an
exact expression for the curvature density parameter $\Omega_{\uK}$ in
terms of the non-linear curvature perturbation $\zeta$. This promotes
$\Omega_{\uK}$ to a stochastic quantity, and in \cref{sec:stats} we
calculate its moments as well as its probability density function,
assuming Gaussian statistics for $\zeta$. Finally, we conclude by
discussing how the statistics of the curvature density parameter is
modified in the presence of non-Gaussian super-Hubble fluctuations, or
if inflation lasted for a very long time.

\section{Curvature density parameter}
\label{sec:curvdef}

When spatial curvature is included, the Friedmann-Lema\^{\i}tre-Robertson-Walker
(FLRW) line-element reads
\begin{equation}
  \dd s^2 = -\dd \tau^2 + a^2(\tau) 
  \frac{  \delta_{ij}\dd x^i \dd x^j}{\left(1+\frac{K}{4}\delta_{mn}x^m x^n\right)^2}\,,
\label{eq:metric:curved:FLRW}
\end{equation}
where $K$ is a constant, and its Ricci scalar is given by
\begin{equation}
\label{eq:Ricci:curved:FLRW}
R=6\frac{\dot{a}^2}{a^2}+6\frac{\ddot{a}}{a}+\frac{6}{a^2}K\, .
\end{equation}

The metric \eqref{eq:metric} can be viewed as an inhomogeneous
generalization of a flat, i.e., $K=0$ FLRW spacetime having a
space-dependent scale factor
\begin{equation}
b(\tau,\vb*{x}) \equiv a(\tau) e^{\zeta(\tau,\vb*{x})},
\end{equation}
from which one can derive the Ricci scalar
\begin{equation}
R = 6 \dfrac{\dot{b}^2}{b^2} + 6
\dfrac{\ddot{b}}{b} + 2 \dfrac{(\vnabla
  b)^2}{b^4} - 4 \dfrac{\Delta b}{b^3}\,.
\label{eq:ricci}
\end{equation}
We now split $\zeta(\tau,\vb*{x})=\xi(\vb*{x}) + \zeta_{\us}(\tau,\vb*{x})$ into a
conserved part $\xi$ (super-Hubble), and time-dependent
fluctuations $\zeta_{\us}$ (sub-Hubble). Expanding in the (presumably
small) short-length part, one has
\begin{equation}
b(\tau,\vb*{x}) = a(\tau) e^{\xi(\vb*{x})} \left[1 + \zeta_{\us}(\tau,\vb*{x}) +
\cdots \right],
\label{eq:pertexpand}
\end{equation}
and upon defining
\begin{equation}
\tilde{a}(\tau,\vb*{x}) = a(\tau) e^{\xi(\vb*{x})}
\end{equation}
one is led to
\begin{equation}
R = 6 \dfrac{\dot{\tilde{a}}^2}{\tilde{a}^2} +
6\dfrac{\ddot{\tilde{a}}}{\tilde{a}} +
\dfrac{6}{\tilde{a}^2}\left[-\dfrac{2}{3} \Delta \xi -
    \dfrac{1}{3} \left(\vnabla \xi\right)^2 \right] + \cdots.
\label{eq:ricciexpand}
\end{equation}
The omitted terms in this expression are the ones appearing in the
linear theory of cosmological perturbations, in the synchronous gauge,
completed by all possible non-linear corrections involving powers of
$\zeta_{\us}(\tau,\vb*{x})$ and products with $\tilde{a}(\tau,\vb*{x})$. The mixed
terms involving both $\tilde{a}(\tau,\vb*{x})$ and powers of
$\zeta_{\us}(\tau,\vb*{x})$ were precisely the ones discussed in the early
works on backreaction and are non-observable~\cite{Geshnizjani:2005ce,
  Hirata:2005ei, Kolb:2005da}. As can be checked in
\cref{eq:ricciexpand}, the terms we have kept are invariant by a
constant shift of $\xi(\vb*{x})$, up to a redefinition of $a(\tau)$.

Since $\xi(\vb*{x})$ varies on super-Hubble length scales only, so does
$\tilde{a}(\tau,\vb*{x})$, hence any observer will identify
$\tilde{a}(\tau,\vb*{x})$ as the FLRW scale factor of their local Hubble
patch. Let us notice that in the gauge we work in, the Hubble radius
is the same for all observers
since~\cite{Matarrese:2003ye,Kolb:2004jg}
\begin{equation}
\tilde{H} \equiv \dfrac{\dot{\tilde{a}}}{\tilde{a}} = \dfrac{\dot{a}}{a} =
H,
\end{equation}
which does not depend on $\vb*{x}$. An important remark is that
\cref{eq:Ricci:curved:FLRW,eq:ricciexpand} coincide upon identifying
\begin{equation}
K = -\dfrac{2}{3} \Delta \xi - \dfrac{1}{3} \left(\vnabla \xi \right)^2,
\label{eq:K}
\end{equation}
which is indeed constant since $\xi$ is conserved, 
and whose measurable curvature density parameter reads
\begin{equation}
\Omega_{\uK} = -\dfrac{K}{\tilde{a}^2 \tilde{H}^2} = -\dfrac{K e^{-2\xi}}{a^2 H^2}\,.
\label{eq:OmegaK}
\end{equation}
Let us stress that Eq.~\eqref{eq:K} is exact in the sense that all the
terms omitted involve $\zeta_{\us}(\tau,\vb*{x})$, hence they are time-dependent
and cannot be absorbed in $K$.
Equation~\eqref{eq:K} makes also explicit that only gradients
of super-Hubble inhomogeneities have a non-trivial effect.

\section{Statistics}
\label{sec:stats}

Current cosmological measurements~\cite{Planck:2019kim} imply that
$\zeta$ has Gaussian statistics and can thus be treated as a random
Gaussian field, with vanishing mean and higher-point correlation
functions entirely determined by the power spectrum
\begin{equation}
\ev{\zeta(\vb*{k}) \zeta(\vb*{k}')} = \left(2 \pi\right)^3 \delta(\vb*{k}+\vb*{k}')P_\zeta(k).
\label{eq:Pk}
\end{equation}
This is also in agreement with the most favored inflationary scenarios, where the
mean values are identified with vacuum expectation values of quantum operators in
the Bunch-Davis vacuum.
Later on, we will also use the spherical power spectrum $\calP_\zeta(k)$
defined by
\begin{equation}
\calP_\zeta(k) = \dfrac{k^3}{2\pi^2} P_\zeta(k) \simeq \calP_*,
\end{equation}
where the last approximation holds for a scale-invariant power
spectrum. 

From \cref{eq:K,eq:OmegaK}, $\Omega_{\uK}$ can therefore also be seen as a
stochastic quantity, though its non-linear dependence on $\xi$, and
thus on $\zeta$, implies that it does not feature Gaussian
statistics. In particular, its expectation value does not necessarily
vanish.

Let us make the decomposition $\zeta(\tau,\vb*{x})=\xi(\vb*{x}) + \zeta_{\us}(\tau,\vb*{x})$ explicit in
Fourier space,
\begin{equation}
\begin{aligned}
  \zeta(\tau,\vb*{x})  & = \dfrac{1}{\left(2\pi \right)^3} \int\dd^3 \vb*{k} \,
  \heaviside{k_{\sigma}-k} \zeta(\vb*{k}) e^{i\vb*{k}\vdot\vb*{x}} \\ & +
  \dfrac{1}{\left(2 \pi \right)^3} \int \dd^3 \vb*{k} \heaviside{k-k_{\sigma}}
    \zeta(\tau,\vb*{k}) e^{i \vb*{k} \vdot \vb*{x}},
  \end{aligned}
\label{eq:fouriersplit}
\end{equation}
where we have introduced a wavenumber $k_{\sigma}$ below which all
Fourier modes $\zeta(\tau,k<k_{\sigma}) = \zeta(k)$ can be approximated
as time-independent. Based on the theory of cosmological
perturbations, and its generalizations~\cite{Rigopoulos:2004gr,
  Langlois:2005ii}, this wavenumber is at most of the order of the
conformal Hubble parameter at the observer's time, say $\tau_\zero$,
namely $k_{\sigma} \lesssim \tilde{a}(\tau_\zero) \tilde{H}(\tau_\zero)$. Let us
remark the presence of $\tilde{a}(\tau_\zero,\vb*{x})$, instead of
$a(\tau_\zero)$, in this expression. A priori, this would induce an
extra-dependence on $\vb*{x}$ in \cref{eq:fouriersplit}, where one should
write $k_{\sigma}(\vb*{x})$. In order to circumvent this issue, we can, for
now, simply choose the cutoff $k_{\sigma}$ to be sufficiently small such
that it encompasses all possible spatial modulations of
$\tilde{a}(\tau_\zero,\vb*{x})$. In other words, we define
\begin{equation}
k_{\sigma} \equiv \sigma a_{\zero} H_{\zero}\,,
\label{eq:ksigma}
\end{equation}
where in principle $\sigma<e^{{\min}_{\vb*{x}}(\xi)}$. As such, we can
identify the conserved quantity with
\begin{equation}
\xi(\vb*{x})=\dfrac{1}{\left(2\pi\right)^3} \int \dd^3\vb*{k} \,
\heaviside{k_{\sigma}-k} \zeta(\vb*{k})
e^{i \vb*{k} \vdot \vb*{x}}\,.
\label{eq:zetacdef}
\end{equation}

\subsection{Mean value}

The mean value of the curvature density parameter reads
\begin{equation}
\ev{\Omega_{\uK}} = -\dfrac{\ev{K e^{-2 \xi}}}{a^2 H^2}\,,
\label{eq:vevdefOmegaK}
\end{equation}
where $\xi$ is given by \cref{eq:zetacdef}.
The curvature scalar $K$, given in \cref{eq:K}, can be split into
two terms $K = K_1 + K_2$ with
\begin{equation}
K_1 \equiv -\dfrac{2}{3} \Delta \xi, \qquad K_2 \equiv
-\dfrac{1}{3} \left(\vnabla \xi \right)^2.
\label{eq:K12def}
\end{equation}
Therefore, one needs the Laplacian and the squared gradient 
of $\xi$. They read
\begin{equation}
\Delta \xi = - \int \dfrac{\dd^3 \vb*{k}}{\left(2\pi\right)^3}
\heaviside{k_{\sigma}-k} k^2 \zeta(\vb*{k}) e^{i \vb*{k} \vdot \vb*{x}},
\label{eq:laplacian}
\end{equation}
and
\begin{equation}
  \begin{aligned}
    \left(\vnabla \xi \right)^2 & = -\int \dfrac{\dd^3
      \vb*{p} \dd^3\vb*{q}}{\left(2\pi\right)^6} \, \heaviside{k_{\sigma}-p}
    \heaviside{k_{\sigma}-q} \\ & \times  \vb*{p} \vdot
    \vb*{q} \, \zeta(\vb*{p}) \zeta(\vb*{q})
    e^{i\left(\vb*{p} + \vb*{q}\right)\vdot\vb*{x}},
  \end{aligned}
  \label{eq:gradmag2}
\end{equation}
from which one can immediately calculate
\begin{equation}
\begin{aligned}
  \ev{K} & = \ev{K_2} = -\dfrac{1}{3} \int \dfrac{\dd^3
    \vb*{k}}{\left(2\pi\right)^3} \heaviside{k_{\sigma}-k} k^2 P_\zeta(k) \\
  & = -\dfrac{1}{3} \int_0^{k_{\sigma}} \dd k k \calP_\zeta(k) \simeq
  -\dfrac{1}{6} k_{\sigma}^2 \calP_*,
\end{aligned}
\label{eq:evK}
\end{equation}
the rightmost equality holding only for a scale-invariant power-spectrum.    

The term $e^{-2\xi}$ appearing in \cref{eq:vevdefOmegaK} can be
expressed in terms of $\zeta(\vb*{k})$ by using the series representation
\begin{equation}
e^{-2\xi} = \sum_{n=0}^{+\infty} \dfrac{(-2)^n}{n!} \xi^n,
\label{eq:expm2zetac}
\end{equation}
with
\begin{equation}
\xi^n = \int \dfrac{\dd^3 \vb*{k}_1 \dots \dd^3 \vb*{k}_n}{\left(2
  \pi\right)^{3n}} \left[\prod_{j=1}^n \heaviside{k_{\sigma} - k_j}
\zeta(\vb*{k}_j) \right] e^{i \vb*{x} \vdot \sum_j \vb*{k}_j}.
\label{eq:zetacn}
\end{equation}
As can be seen in \cref{eq:vevdefOmegaK}, the mean value of the
curvature density parameter requires the explicit determination of an
infinite number of terms, the non-vanishing ones being of the form
$\ev{K_1 \xi^{2p + 1}}$ and $\ev{K_2 \xi^{2p}}$. From
\cref{eq:Pk,eq:K12def,eq:zetacn}, one can make extensive use of the
Wick's theorem to reduce all the expectation values to a few
two-point functions with the following diagrammatic rules:
\begin{equation}
  \begin{aligned}
\ev{K_1 K_1} & \equiv  \diagKoneKone  = \ev{K_1^2},\\
\ev{\xi \xi} & \equiv  \diagZetaZeta  = \ev{\xi^2},\\
\ev{K_1\xi} & \equiv  \diagKoneZeta  = - 2 \ev{K},\\
 \ev{K_2} & \equiv \diagKtwoSelf  = \ev{K}.
\end{aligned}
\label{eq:diagrules}
\end{equation}
Let us notice that, due to the inner product structure of
\cref{eq:gradmag2}, the $K_2$ vertices have two ``legs'' that can only
connect to other $K_2$ vertices.  From \cref{eq:laplacian}, one has
\begin{equation}
  \ev{K_1^2}  = \dfrac{4}{9} \int_0^{k_{\sigma}} \dd k k^3 \calP_\zeta(k) \simeq
  \dfrac{1}{9} k_{\sigma}^4 \calP_*,
\label{eq:m2K1}
\end{equation}
which allows us to express the second moment of the curvature scalar
as
\begin{equation}
\ev{K^2} = \ev{K_1^2} + \dfrac{5}{3} \ev{K}^2 \simeq \dfrac{1}{9}
k_{\sigma}^4 \calP_* \left(1+ \dfrac{5}{12} \calP_*\right).
\label{eq:m2K}
\end{equation}
In \cref{eq:diagrules}, we also need the variance of the conserved
quantity $\xi$. It can be determined from \cref{eq:zetacdef} and reads
\begin{equation}
\ev{\xi^2} = \int_{k_{\varepsilon}}^{k_{\sigma}}\dd k \dfrac{\calP_\zeta(k)}{k} \simeq
\calP_* \ln\left(\dfrac{k_{\sigma}}{k_{\varepsilon}} \right) \simeq \calP_* N_\uinf,
\label{eq:varzetac}
\end{equation}
where we have introduced an expected infrared cutoff $k_{\varepsilon}$. Indeed,
in the context of Cosmic Inflation, the ratio between the largest and
shortest lengths being amplified is precisely given by the total amount of
stretching generated by the accelerated expansion, the so-called total
number of {$e$-folds} $N_\uinf$. For the measured value of $\calP_* =
2.1\times 10^{-9}$~\cite{Planck:2018jri}, and a not too long
inflationary era $N_\uinf \ll 10^{9}$, $\ev{\xi^2}$ is a small
quantity.

Denoting by $W_{2p}=(2p)!/(p!2^p)$ the number of Wick's contractions
between $p$ pairs, one obtains
\begin{equation}
\begin{aligned}
  \ev{K_1 \xi^{2p+1}} & = (2p+1) \left(\diagKoneZeta\right) \times
  W_{2p} \left(\diagZetaZeta\right)^{p}\\
  & =  -\dfrac{(2p+1)!}{p! \, 2^{p-1}} \ev{K}
  \ev{\xi^2}^p,
\end{aligned}
\label{eq:evK1zetap}
\end{equation}
and
\begin{equation}
\begin{aligned}
  \ev{K_2 \xi^{2p}} & = \diagKtwoSelf \times W_{2p}
  \left(\diagZetaZeta\right)^{p} \\
  &= \dfrac{(2p)!}{p! \, 2^p} \ev{K} \ev{\xi^2}^p.
\label{eq:evK2zetap}
\end{aligned}
\end{equation}
The infinite series obtained by combining
\cref{eq:evK1zetap,eq:evK2zetap,eq:expm2zetac} can be resumed and one gets the
exact expression
\begin{equation}
\ev{\Omega_{\uK}} = - \dfrac{5}{a^2 H^2} \ev{K} e^{2 \ev{\xi^2}}.
\label{eq:evOmegaK}
\end{equation}
Making use of \cref{eq:evK,eq:varzetac}, for a scale-invariant power
spectrum, \cref{eq:evOmegaK} simplifies to
\begin{equation}
\ev{\Omega_{\uK_{0}}} \simeq \dfrac{5}{6} \dfrac{k_{\sigma}^2}{a_{\zero}^2 H_{\zero}^2} \calP_*
e^{2\calP_* N_\uinf} \simeq \dfrac{5}{6} \sigma^2 \calP_*,
\label{eq:evOmegaKnum}
\end{equation}
which saturates for $\sigma=1$ at $\ev{\Omega_{\uK_{0}}} \simeq 1.7 \times
10^{-9}$, a barely open universe.

\subsection{Variance}

There is little hope to measure such a small value of $\ev{\Omega_{\uK}}$,
but $\Omega_{\uK}$ being a stochastic variable, its realizations are also
dictated by the higher moments, the second one being given by
\begin{equation}
\ev{\Omega_{\uK}^2} = \dfrac{\ev{K^2 e^{-4 \xi}}}{a^4 H^4} =
\dfrac{\ev{\left(K_1^2 + 2 K_1 K_2 + K_2^2\right)e^{-4 \xi}}}{a^4 H^4}\,.
\label{eq:m2defOmegaK}
\end{equation}
Using again a series representation for the exponential,
\cref{eq:m2defOmegaK} can be expanded in an infinite sum requiring the
calculation of the non-vanishing terms $\ev{K_1^2 \xi^{2p+2}}$,
$\ev{K_1 K_2 \xi^{2p+1}}$, $\ev{K_2^2 \xi^{2p}}$, with $p\ge0$.
Using the diagrammatic rules of \cref{eq:diagrules}, one gets
\begin{equation}
  \begin{aligned}
   &  \ev{K_1^2 \xi^{2p+2}}  = \diagKoneKone \times W_{2p+2}
    \left(\diagZetaZeta\right)^{p+1} \\ & +
    (2p+1)\left(\diagKoneZeta\right) \times 2p
    \left(\diagKoneZeta\right) \\
    & \times W_{2p} \left(\diagZetaZeta\right)^p \\
    & = \dfrac{(2p+1)!}{p!\, 2^{p}}
    \ev{K_1^2}\ev{\xi^2}^{p+1} + 4 \dfrac{(2p+2)!}{p!\,2^p}\ev{K}^2\ev{\xi^2}^p,
  \end{aligned}
  \label{eq:m2K1zetap}
\end{equation}
together with
\begin{equation}
\begin{aligned}
  \ev{K_1K_2 \xi^{2p+1}} & = \diagKtwoSelf \times (2p+1)
  \left(\diagKoneZeta\right) \\ & \times W_{2p}
  \left(\diagZetaZeta\right)^{p} \\
  & = -2 \dfrac{(2p+1)!}{p!\,2^{p}}
  \ev{K}^2 \ev{\xi^2}^p,
\end{aligned}
\label{eq:K1K2zetap}
\end{equation}
and
\begin{equation}
\begin{aligned}
  \ev{K_2^2 \xi^{2p}} & = \left[\left(\diagKtwoSelf\right)^2  + 2
  \left(\diagKtwoCycleII\right) \right] \\ & \times
  W_{2p}\left(\diagZetaZeta\right)^{p}  = \dfrac{5}{3} \dfrac{(2p)!}{p!\,2^p}
  \ev{K}^2 \ev{\xi^2}^p.
\end{aligned}
\label{eq:m2K2zetap}
\end{equation}
Summing all the terms coming from the expansion of
\cref{eq:m2defOmegaK} gives the exact expression
\begin{equation}
\ev{\Omega_{\uK}^2} = \dfrac{1}{a^4 H^4} \left(\ev{K^2} + 80 \ev{K}^2\right)e^{8\ev{\xi^2}}.
\label{eq:m2OmegaK}
\end{equation}
For a scale-invariant power spectrum, using
\cref{eq:evK,eq:varzetac,eq:m2K}, one obtains
\begin{equation}
  \ev{\Omega_{\uK_{0}}^2}  \simeq \dfrac{1}{9} \dfrac{k_{\sigma}^4}{a_{\zero}^4 H_{\zero}^4}
  \calP_*\left(1+\dfrac{245}{12} \calP_*\right)e^{8\calP_*
    N_\uinf}
   \simeq \dfrac{1}{9} \sigma^4 \calP_*.
\label{eq:m2OmegaKnum}
\end{equation}
Using \cref{eq:evOmegaKnum} for $\sigma=1$, the standard deviation of
$\Omega_{\uK_{0}}$ is given by
\begin{equation}
\sqrt{\ev{\Omega_{\uK_{0}}^2}-\ev{\Omega_{\uK_{0}}}^2} \simeq \dfrac{\sigma^2}{3}
\sqrt{\calP_*} \simeq 1.5 \times 10^{-5}.
\end{equation}

To summarise, \cref{eq:evOmegaK,eq:m2OmegaK} show that, in a Universe
filled with cosmological fluctuations stretched over super-Hubble
scales, the curvature density parameter is not vanishingly small but
is promoted to a stochastic variable. At any time in the cosmic
history, we therefore expect an observer to measure a realization of
$\Omega_{\uK}$ dominated by its standard deviation, i.e., at about $1.5
\times 10^{-5}$. However, because the probability distribution is, a
priori, non-Gaussian, the rarity of extreme values of $\Omega_{\uK}$ could
be affected by the higher moments and we now turn to their
calculation.

\subsection{Higher moments}

All the higher moments $\ev{\Omega_{\uK}^n}$ with $n>2$ can be explicitly
calculated with the same method as the one employed for the mean value and
the variance. Expanding the exponential in series and using the binomial
expansion of $(K_1+K_2)^n$ shows that one has to determine the mean
value of combinations of the form $\ev{K_1^p K_2^q
  \xi^m}=\ev{K_2^q}\ev{K_1^p \xi^m}$. Those can all be
expressed in terms of powers of $\ev{\xi^2}$, $\ev{K}$ and
$\ev{K^2}$ by using the diagrammatic rules of \cref{eq:diagrules}.

The only new subtlety consists in evaluating the terms in $\ev{K_2^q}$
that need to be decomposed into ``self-cycles''. For instance, the
third moment requires to evaluate
\begin{equation}
\begin{aligned}
  \ev{K_2^3} & = 2^2 2!\left(\diagKtwoCycleIII\right) + 2 \left(\diagKtwoCycleII\right) \times 3
\left(\diagKtwoSelf\right) \\ & + \left(\diagKtwoSelf\right)^3,
\end{aligned}
\label{eq:meanK23}
\end{equation}
and one obtains
\begin{equation}
\begin{aligned}
  \ev{\Omega_{\uK}^3} & = -\dfrac{\ev{K}}{a^6 H^6}\left(39 \ev{K^2} +
  \dfrac{19430}{9} \ev{K}^2\right) e^{18 \ev{\xi^2}}.
\end{aligned}
\end{equation}
Similarly, the fourth moment is given by
\begin{equation}
\begin{aligned}
\ev{\Omega_{\uK}^4} & = \dfrac{1}{a^8 H^8}\left(3 \ev{K^2}^2 + 1728
\ev{K^2}\ev{K}^2 \right. \\ & \left. + \dfrac{736682}{9}
\ev{K}^4\right) e^{32\ev{\xi^2}},
\end{aligned}
\end{equation}
and so on and so forth. These expressions are not particularly
illuminating, but the leading order terms of all the moments are
diagrammatically tractable and one can show that, for a
scale-invariant power spectrum, the standardized moments $ \tilde{\mu}_n$ (the moments
divided by the $n^\textrm{th}$ power of the standard deviation) verify
\begin{equation}
\begin{aligned}
  \tilde{\mu}_{n=2p} &\simeq W_{n} \,e^{ \left( 2n^2-4n
    \right)\ev{\xi^2}},\\ \tilde{\mu}_{n=2p+1} & \simeq n W_{n-1}
  \left(1+4n\right) \dfrac{\sqrt{\calP_*}}{2} \, e^{ \left( 2n^2-4n
    \right)\ev{\xi^2}}.
\end{aligned}
\label{eq:mutilde}
\end{equation}
All odd standardized moments are suppressed by the factor
$\sqrt{\calP_*}$ with respect to the even ones. Moreover, provided
the exponential terms in \cref{eq:mutilde} are close to unity, i.e.,
for $n^2 \ev{\xi^2} \ll 1$, the even moments exactly match the
ones associated with a Gaussian probability distribution. As such,
$\ev{\Omega_{\uK}^n}$ shows significant deviations compared to the
Gaussian expectations only for large values of $n \gtrsim
1/\ev{\xi^2}$. To better assess the effect of these higher moments, we next turn our attention to 
 the functional form of the $\Omega_{\uK}$'s probability
distribution.

\subsection{Probability distribution}

\begin{figure}
  \begin{center}
    \includegraphics[width=\figw]{probomega}
    \caption{Probability distribution function for $\bar{\Omega}_{\uK} =(a
      H/k_{\sigma})^2 \Omega_{\uK}$ (red curve) for unrealistically large
      values of $\calP_* = 10^{-3}$ (and $N_\uinf=100$), compared to a
      Gaussian of same mean and variance (black curve). Notice that
      the most probable value of $\bar{\Omega}_{\uK}$ is slightly negative
      (closed universe) whereas the mean value remains slightly
      positive (open universe).}
    \label{fig:proba}
  \end{center}
\end{figure}

The probability density function of $\Omega_{\uK}$ can be determined by
noticing that \cref{eq:K,eq:OmegaK} imply that $\Omega_{\uK}$ can be seen
as a non-linear functional over five
stochastic Gaussian variables, $\vb*{\Xi} \equiv (\xi, \Delta \xi,
\vnabla\xi)$. As such, defining $\bar{\Omega}_{\uK} \equiv (a^2 H^2/k_{\sigma}^2)
\Omega_{\uK}$ and marginalizing over the five-dimensional space associated
with $\vb*{\Xi}$, one has
\begin{equation}
P(\bar{\Omega}_{\uK}) = \int \dfrac{\dd[5]\vb*{\Xi}}{\left(2\pi\right)^{5/2}}
\, \delta\left(\bar{\Omega}_{\uK} + \dfrac{K}{k_{\sigma}^2} e^{-2 \xi}\right)
\dfrac{e^{-\frac{1}{2} \vb*{\Xi}^{\mathrm{T}} \vb*{\Sigma}^{-1} \vb*{\Xi}}}{\sqrt{\det \vb*{\Sigma}}},
\label{eq:Pmulti}
\end{equation}
where the five-dimensional covariance matrix $\vb*{\Sigma}$ is completely
determined by the diagrammatic rules of \cref{eq:diagrules}. All but
one integral appearing in \cref{eq:Pmulti} can be analytically reduced
and, after some algebra, one obtains
\begin{equation}
\begin{aligned}
  P(\bar{\Omega}_{\uK}) & = \dfrac{k_{\sigma}^2}{4
    \pi}\sqrt{\dfrac{27\sqrt{2}\Sigma_\omega}{\ev{\xi^2}\abs{\ev{K}}^{3}}} \int_{-\infty}^{+\infty} \dd x \,  e^{-\frac{x^2}{2\ev{\xi^2}} + 2x}
   \\ & \times  e^{-\frac{1}{2}\bar{\omega}^2(\bar{\Omega}_{\uK},x)}\hermiteb{-\frac{3}{2}}{\dfrac{3
      \Sigma_\omega}{\sqrt{8} \abs{\ev{K}}} - \dfrac{\bar{\omega}(\bar{\Omega}_{\uK},x)}{\sqrt{2}}},
\end{aligned}
\label{eq:Pone}
\end{equation}
where we have defined
\begin{equation}
\begin{aligned}
  \Sigma_\omega^2 & \equiv \ev{K_1^2} - 4 \dfrac{\ev{K}^2}{\ev{\xi^2}},\\
  \bar{\omega}(\bar{\Omega}_{\uK},x) &\equiv e^{2x}\dfrac{\bar{\Omega}_{\uK}}{(\Sigma_\omega/k_{\sigma}^2)} +
  \dfrac{2 \abs{\ev{K}}}{\ev{\xi^2} \Sigma_\omega} x.
\end{aligned}
\end{equation}
In \cref{eq:Pone}, $\hermite{\nu}{x}$ stands for the generalized Hermite
polynomial of fractional order, defined from the parabolic cylinder
functions~\cite{1980tisp.book.G} as
$\hermite{\nu}{x}\equiv2^{\nu/2}e^{x^2/2} D_\nu(\sqrt{2}x)$. This
distribution shows that, for $\ev{\xi^2} \simeq \calP_* N_\uinf \ll
1$, one can use the approximation
\begin{equation}
e^{-\frac{x^2}{2\ev{\xi^2}} + 2x} \simeq \sqrt{2\pi \ev{\xi^2}}
e^{2\ev{\xi^2}} \delta\left(x-2\ev{\xi^2}\right),
\label{eq:deltax}
\end{equation}
to simplify the integral over $x$ in \cref{eq:Pone}. Remarking that,
in this limit, the argument of the Hermite function is dominated by
the first term, which is a constant scaling as $1/\sqrt{\calP_*}$,
$P(\bar{\Omega}_{\uK})$ is therefore close to a Gaussian distribution over the
quantity $\bar{\omega}(\bar{\Omega}_{\uK},2\ev{\xi^2})$. In other words, for
$\ev{\xi^2}\ll 1$, the distribution of $\bar{\Omega}_{\uK}$ is almost
Gaussian, with a width given by $\Sigma_\omega/k_{\sigma}^2
\simeq\sqrt{\calP_*}/3$ and a peak located at a very small
\emph{negative} value
\begin{equation}
\left. \bar{\Omega}_{\uK} \right|_{\max} \simeq \dfrac{4}{k_{\sigma}^2} \ev{K}
e^{-4\ev{\xi^2}} \simeq -\dfrac{2}{3} \calP_*.
\label{eq:OmegaKbar:max}
\end{equation}
For the curvature parameter today, one would get the most probable
value at $\left.\Omega_{\uK_{0}}\right|_{\max} \simeq -1.4\times 10^{-9}$,
a barely closed universe. Let us notice the different sign than the
mean value of \cref{eq:evOmegaKnum}, the distribution is indeed
slightly skewed by the Hermite function. This can be seen in
\cref{fig:proba}, where we have plotted $P(\bar{\Omega}_{\uK})$ for an
irrealistically large value of $\calP_*=10^{-3}$. These distortions
are also apparent in the odd moments of \cref{eq:mutilde} which are,
as already noted, all proportional to $\sqrt{\calP_*}$.

\begin{figure}
  \begin{center}
    \includegraphics[width=\figw]{probomega_tails}
    \caption{Probability distribution function for $\bar{\Omega}_{\uK} =(a
      H/k_{\sigma})^2 \Omega_{\uK}$ (red curve) for the currently favored
      value of $\calP_* = 2.1 \times 10^{-9}$ and for a large number
      of {$e$-folds} $N_\uinf=10^8$. The variance $\ev{\xi^2}$ is no
      longer a small quantity and the distribution acquires heavy
      tails. Even though the width at half-maximum is
      $\order{\sqrt{\calP_*}}$, substantial values of $|\bar{\Omega}_{\uK}|$
        are not rare anymore. For comparison, the black curve shows a
        Gaussian of same mean and variance.}
    \label{fig:tails}
  \end{center}
\end{figure}

When $\calP_* N_\uinf = \ev{\xi^2}$ increases, \cref{eq:deltax} is
no longer accurate and all the terms of \cref{eq:Pone} are
relevant. The distribution now acquires heavy tails, kicking in at
increasingly smaller values of $\left|\bar{\Omega}_{\uK}\right|$ and erasing
the Gaussian profile in the neighboorhood of
$\left.\bar{\Omega}_{\uK}\right|_{\max}$. In \cref{fig:tails}, we have plotted
$P(\bar{\Omega}_{\uK})$, in logarithmic scales, for
$\calP_*=2.1\times10^{-9}$ and for a large number of {$e$-folds}
$N_\uinf=10^8$.  These heavy tails imply that large values of $\vert
\Omega_{\uK_{0}}\vert$ are (much) more likely than what a Gaussian profile
would imply. Their existence is also manifest in the moments of
\cref{eq:mutilde} through the exponential coefficients involving
$\ev{\xi^2}$. Such an effet is reminiscent of the non-linear
mapping of vacuum quantum fluctuations encountered in the context of
stochastic inflation~\cite{Pattison:2017mbe, Ezquiaga:2019ftu}.

Finally, let us mention that numerical computations of $\ev{\Omega_{\uK}}$
and $\ev{\Omega_{\uK}^2}$ based on using the distribution of \cref{eq:Pone}
do match the values we can get from \cref{eq:evOmegaK,eq:m2OmegaK}.

\section{Discussion}
\label{sec:discuss}

If inflation lasts for a long period, then substantial values of
$\Omega_{\uK_{0}}$ might be produced. Indeed, letting $\sigma^2\simeq
e^{-\ev{\xi^2}}$ to implement the condition stated below
\cref{eq:ksigma}, \cref{eq:m2OmegaKnum} becomes
$\ev{\Omega_{\uK_{0}}^2}^{1/2} \simeq \sqrt{\calP_*}
e^{3\calP_*N_\uinf}/3$. For this value not to exceed the observational
bound $\vert\Omega_{\uK_{0}}\vert<3\times 10^{-3}$, with
$\calP_*=2.1\times 10^{-9}$ this leads to $N_\uinf<7\times 10^8$. This
suggests that scenarios leading to phases of inflation lasting for
more than a billion {$e$-folds} might be disfavored by current
cosmological data. Notice that this bound becomes more stringent if
one accounts for the slightly red observed spectral index.

Let us note, however, that when the above bound on $N_\uinf$ is saturated,
$\ev{\xi^2}\simeq 1.5$.  A priori, our non-linear formulas do not
require $\ev{\xi^2}$ to be small, hence they can still be used in
that case. In particular, although one can see that all the moments
are becoming exponentially large with $\ev{\xi^2}$, \cref{eq:Pone}
shows that $P(\bar{\Omega}_{\uK})$ remains well-defined. Nonetheless, the fact
that the scale $k_{\sigma}$ must be set in a way that accommodates
potentially large values of $\zeta$ suggests that our formalism may
not be best suited in that case, and the upper bound we have obtained
on $N_\uinf$ must be taken with care.

If inflation lasts even longer, $\Omega_{\uK}$ gets even larger and our
formalism needs to be extended in at least two ways. First, when
$|\Omega_{\uK}|$ becomes of order unity, or more, the metric associated
with \cref{eq:metric} is not acceptable anymore. For instance, a large
negative curvature density parameter would imply a compact manifold
and this demands another coordinate system than the one of
\cref{eq:metric}. Second, when $|\Omega_{\uK}|$ becomes sizeable, it opens
up a channel of backreaction of the curvature perturbation onto the
background dynamics, which in turn alters the inflationary
amplification of the curvature perturbations
themselves~\cite{Handley:2019anl, Letey:2022hdp}. This mechanism might
be tractable in an extended stochastic-inflation
formalism~\cite{Starobinsky:1986fx, Goncharov:1987ir,
  Starobinsky:1994bd, Vennin:2015hra, Grain:2020wro}, which we plan to
develop in a future work.
  
Finally, our results would be modified if curvature perturbations are
non-Gaussian at non-observably large scales. This is, strictly
speaking, not excluded, although it would require very specific
early-universe models for which curvature perturbations are Gaussian
at observable scales today (in order to satisfy the tight constraints
on non-Gaussianities~\cite{Planck:2019kim}), and non-Gaussian at
larger scales.  Another hypothesis that could be broken is that
$\xi$ is conserved by adiabaticity. The presence of entropic modes
today could invalidate this assumption, but, as for non-Gaussianities,
their presence during inflation is also disfavored by current
data. Some exotic effects would then be required to find them only
now~\cite{Gaztanaga:2022bdh}.

