\section{Introduction}
\begin{figure*}[!t]
    \centering
    \includegraphics[scale=0.56]{figures/pipeline-white.drawio.pdf}
    \caption{SAINE Workflow and Pipeline.}
    \label{fig:pipeline}
\end{figure*}
A precise classification of publications across and within disciplines is key not only for a fast and comprehensive search to guide researchers to relevant material but also to identify the novelty of research, the standing and significance of scholars as well as their home institutions, and of the relative growth of fields of work.

Machine learning develops into being not only \textit{a} but \textit{the} customary approach to establish such a classification. Clearly, one would expect a search that is geared to identifying a high-quality corpus of keywords to benefit crucially from supervision. Existing classifications of academic output are based on a blend of (supervised) author-chosen and (unsupervised) machine-chosen keyword lists, where the composition of the blend is unknown to the researcher.

Prevailing systems of keywords for academic publications are lists based on abstracts in a discipline, field, and subfield distilled from
\begin{itemize}
\item unsupervised machine learning (from word or phrase frequencies);
\item supervised learning (mostly from keyword self-reporting by authors);
\item semi-supervised learning (a mixture of the two; e.g., as done by Microsoft Academic Graph (MAG) described in~\citet{sinha2015overview, 10.3389/fdata.2019.00045, wang2020microsoft}).
\end{itemize}

To design an annotation and inference engine that helps us understand the publication space better, \textbf{we should cater to the following needs:} (1) a simple user interface with clear annotation instructions; (2) a reproducible pipeline across various disciplines; (3) good support for inference tailored to downstream tasks (e.g., model retraining) in meta-science studies.

Among the existing open-source annotation tools, Label Studio \cite{LabelStudio} suits our needs. Note that \citet{gayoso2019review} have reviewed extensively on the annotation tools that facilitate classification tasks. The reviewed tools are either not open-sourced or are domain-specific and hence do not suit our purposes.

In this system demonstration, we utilize a set of standard open-source software, mainly Label Studio \cite{LabelStudio}, \href{https://mlflow.org/}{MLflow} and \href{https://fastapi.tiangolo.com/}{FastAPI} to configure an annotation and inference engine for scientific publication annotations. In this demonstration, we illustrate the benefit of using supervised learning based on pre-established keyword lists and abstracts, and how annotators can help us better understanding the importance of \textit{supervised} learning in establishing a classification of academic publications.

This system is built on top of the hitherto largest scale of \textit{multi-class} hierarchical classification studies across all disciplines in both \textit{single-label} and \textit{multi-label} settings (cf.~\citet{rao2023hierarchical}). There, we have built a supervised hierarchical classification system that associates every publication with at least one and potentially several disciplines, fields, and subfields.

With the annotations above, we conduct a small user study with domain experts using our \textbf{annotation} engine. We then invoke our \textbf{inference} engine to fine-tune the base models in \citet{rao2023hierarchical}. The comparison between the base and fine-tuned models shows that the proposed annotation and inference system is able to benefit the further development of a more accurate classification.

To summarize, the paper presents a scientific annotation and inference engine called SAINE, which is based on open-source software like Label Studio and MLflow. \textbf{The main contributions of the paper are:}
(1) The demonstration of using SAINE in understanding the space for scholarly publications, particularly in hierarchical discipline classifications.
(2) The result of a user study, which shows that user input collected with the help of SAINE can help better understand the classification process.
(3) The ability of SAINE to benefit the further development of a more accurate classification, demonstrated through the comparison between the base and fine-tuned models.
(4) The potential of SAINE to support downstream meta-science projects and foster greater transparency and understanding of scientific research.

Overall, the paper presents the benefits of supervised learning and the importance of having a simple user interface with clear annotation instructions, reproducible pipelines, and good support for inference in scientific publication annotations. The \href{https://app.heartex.com/user/signup/?token=e2435a2f97449fa1}{live demo website} and \href{https://youtu.be/yToO-G9YQK4}{demonstration video} are also available for those interested in further exploring SAINE. The codebase for development is publically available \href{https://gitlab.ethz.ch/raox/science-clf/-/tree/labelstudio}{under this link} and collocates with the codebase of \citet{rao2023hierarchical}.

In Figure~\ref{fig:pipeline} we illustrate the workflow in SAINE by assigning the roles of "Administrator", "Annotators", "Label Studio", and "MLflow" to each task in the pipeline. The sections are organized as follows. Section~\ref{sec:labelstudio} introduces the functionality of Label Studio and its fit to our annotation needs, as well as our annotation guidelines for experts. Section~\ref{sec:implementation} specifies the annotation design for the field of \textit{Economics} and discusses the annotation results. Section~\ref{sec:inf} discusses the integration of annotation results into the pre-trained base models and fine-tuned ones with MLflow. We then conclude this system demonstration with a discussion of system limitations, ethics, and broader impact statements.

\section{Annotating Scientific Articles with Label Studio}
\label{sec:labelstudio}
We briefly introduce functionalities of Label Studio and illustrate why Label Studio is a suitable tool to annotate scientific articles. \textbf{The goals of annotations in the present context are three-fold:}
\textbf{(1) [To judge the appropriateness of an assigned category.]} Experts judge whether a scientific article has a correctly assigned category. If a category does not suit the abstract's content, the annotator will pick a new category from a pre-defined list.
\textbf{(2) [To evaluate keywords.]} Experts evaluate keywords assigned by MAG and mark missing keywords in the abstract.
\textbf{(3) [To calculate inter-annotator agreement.]} The annotation engine should permit an efficient way of calculating inter-annotator agreement (IAA) scores among annotators.

\subsection{Label Studio}
Label Studio is a powerful and versatile annotation tool that can handle various types of annotation tasks. Here are some of the features of Label Studio that make it a suitable tool for the annotation tasks.

\begin{itemize}
    \item \textbf{[Customizable interface.]} Label Studio allows to design a customized annotation interface. One can create a pre-defined list of categories for the experts to choose from and provide them with clear instructions on how to evaluate the assigned category.
    \item \textbf{[Multiple annotation types.]} Label Studio supports various types of annotation, including text classification, entity recognition, and sequence labeling. Therefore, different types of annotation can be used to evaluate keywords, mark missing keywords, and judge predicted categories.
    \item \textbf{[Collaboration and version control.]} Label Studio enables multiple experts to work on the same project simultaneously, allowing efficient and collaborative annotation. It also includes a version control system that tracks changes to the annotations, facilitating easy comparison and IAA evaluation.
    \item \textbf{[Inter-annotator agreement (IAA).]} Label Studio has built-in tools to calculate IAA scores. These tools can help evaluate the consistency and reliability of chosen annotations.
    \item \textbf{[Integrations with machine-learning models.]} Label Studio also provides integration with various machine learning (ML) models. Although we do not use the integrated ML functions, Label Studio allows us to export the annotation results in JSON, with which we improve the classification models using the annotated data in the inference engine.
\end{itemize}

Overall, Label Studio offers a powerful and customizable annotation platform that can handle relevant annotation tasks, facilitate efficient collaboration among experts, and efficiently compute IAA.

\begin{figure*}[!t]
    \centering
    \includegraphics[scale=0.25]{./figures/labelstudio_admin.png}
    \caption{Administrative Panel of Annotation Tasks in Label Studio.}
    \label{fig:labelstudio_admin}
\end{figure*}

In Figure~\ref{fig:labelstudio_admin} we demonstrate the administrative panel of the project manager. The project manager uses this panel to assign annotation tasks to each registered annotator and can monitor the annotation progress. The manager can also adjust the assigned annotations based on individual progress. The "Filters" and "Order (Annotation results)" tabs make it easy to inspect tasks by annotation progress (e.g., "Annotators", "Agreement", "Completed", "Total annotations per task").

\subsection{Annotation Guidelines}
When a publication is annotated, each annotator is provided with the abstract, the keywords offered by MAG, and the assigned category based on the keywords provided by MAG. The categories of a discipline classification (like JEL) are assigned to MAG publications on the basis of the keywords. Therefore, MAG's keywords help us identify potential misalignments and better understand the classifiers we built.

The annotation samples provided in the annotation engine are stratified sampled (ratio: 2e-5) across all classes of the training set introduced by \citet{rao2023hierarchical} for one discipline. Each annotator is required to judge whether a category is correctly assigned to an abstract. If not, the annotator is required to select the suitable one from a predefined list. The annotator is also required to evaluate MAG-generated keywords and make corrections (by removing unqualified keywords/marking suitable keywords from the abstract). Figure~\ref{fig:labelstudio_multipleannotator} shows two annotations of one publication. Label Studio makes it easy to navigate among the annotations generated by various annotators on an identical instance. Note that, as we discussed in \citet{rao2023hierarchical}, our multi-class hierarchical classification system is modularized in both \textit{single-label} and \textit{multi-label} settings. The current annotation engine is equipped with both annotation functionalities. For the sake of system demonstration and user study in Section~\ref{sec:implementation}, we discuss the \textit{single-label} setting. More details on the multi-label setting are provided in Appendix~\ref{app:multi-label}.

\begin{figure*}[!t]
    \centering
    \includegraphics[scale=0.4]{./figures/labelstudio_multipleannotator.png}
    \caption{Publication Annotations by Multiple Annotators.}
    \label{fig:labelstudio_multipleannotator}
\end{figure*}

\section{Implementation: User Study in \textit{Economics}}
\label{sec:implementation}

We now use \textit{Economics} as a discipline to show how we utilize the annotation engine to collect expert annotations.

\subsection{Annotation Design}
We invited three economist experts from the Chair of Applied Economics at ETH Zurich to join the annotation project by accessing \href{https://app.heartex.com/user/signup/?token=e2435a2f97449fa1}{this link}. The annotation guidelines are stated at \href{https://gitlab.ethz.ch/raox/science-clf/-/blob/labelstudio/labelstudio/annotation_guideline.md}{here}. Of the three experts, one has annotated all provided instances (Annotator 1), one has annotated 10\% of the instances (Annotator 2), and one has annotated a subset of instances with an ex-ante denomination in Urban and Spatial Economics only (Annotator 3). Each annotator was provided a user panel shown in Figure~\ref{fig:labelstudio_annotator} in Appendix~\ref{app:more-function}.

\subsection{Annotation Results in Label Studio}

Altogether, 788 instances of abstracts and keywords from MAG had to be annotated for a \textit{single-label} classification. In Economics, a standardized field and subfield system with keywords exists, and it is called \href{https://www.aeaweb.org/econlit/jelCodes.php?view=jel}{the Journal of Economic Literature (JEL) classification system}. This system is known to all academic economists and serves as a guiding principle to associate an article or a topic with a specific subfield in \textit{Economics}. The subfields in the JEL categories are associated with keywords.

We report the annotation time and IAA scores that are automatically calculated by Label Studio (see \href{https://labelstud.io/guide/stats.html}{the official documentation} for the steps). The final task agreement score is calculated by averaging all IAA scores for each annotation pair. Table~\ref{tab:iaa} illustrates the IAA scores amongst three experts. Annotators 1, 2, and 3 have annotated 788, 181, and 99 instances, respectively. The annotations overlap between the pairs is 4 or 7\% of the overlapping instances (Annotators 2 and 3), 99 or 100\% of the overlapping instances (Annotators 1 and 3), and 181 or 100\% of the overlapping instances (Annotators 1 and 2). The median annotation time of Annotators 1-3 per instance was 17.7s, 29.8s, and 40.9s, respectively. The annotators were entitled to disapprove of the assigned category based on MAG upon suggesting an alternative category. Marking and filling in missing keywords is time consuming, reading the MAG-generated keywords can help to some extent the annotation speed. However, all annotators reported that the MAG-provided keywords could be a source of error for wrongly assigned categories. As discussed among the annotators after they underwent the annotations separately, the category they found the best was for \textit{Mathematical \& Quantitative Methods}, and it was worst for \textit{Macroeconomics} and \textit{Public Economics}.

\begin{table}[!t]
\centering
\resizebox{\linewidth}{!}{
\begin{tabular}{cccc}
\toprule
                     & \textbf{Annotator 1} & \textbf{Annotator 2} & \textbf{Annotator 3} \\
\midrule
\textbf{Annotator 1} &                      & 55\%                 & 55\%                 \\
\textbf{Annotator 2} & 55\%                 &                      & 27\%                 \\
\textbf{Annotator 3} & 55\%                 & 27\%                 &    \\
\bottomrule
\end{tabular}}
\caption{Annotator Agreement Matrix Among Three Expert Annotators.}
\label{tab:iaa}
\end{table}

\begin{figure*}[!t]
    \centering
    \includegraphics[scale=0.3]{figures/inference_3_crop.png}
    \caption{Inference Engine with \textit{MLflow} Integration. API: Batch Inference by Model.}
    \label{fig:inference}
\end{figure*}

\section{Inference Engine: Incorporating Annotation Results into the Existing Classification Pipeline}
\label{sec:inf}
We illustrate the pipeline using the discipline \textit{Economics} as discussed in Section \ref{sec:implementation}.

\subsection{Post-processing of Annotation Results}
We downloaded the annotation results in JSON (\href{https://gitlab.ethz.ch/raox/science-clf/-/blob/labelstudio/labelstudio/labeled_data/labelstudio_econ.json}{here}) of all experts and post-processed them following the protocols below, before feeding them into the pre-trained base models of various neural networks as discussed in \citet{rao2023hierarchical}. In total, we obtained 1,068 partly overlapping annotations (incl.~"Skip", "(Dis)agree", keywords, added categories). The basic statistics on the number of instances of "Agree", "Disagree" and "Not ECON" are 498, 297, and 268, respectively.

Here is the post-processing procedure. (1) We removed abstracts that were inadequately classified as belonging in Economics from the sample (206 of 788 instances). Additionally, we deleted 5 instances due to bad annotations. For example, no one labeled this sample ("Skip"), or an annotator chose "Disagree" but did not choose a new category.
(2) For each remaining instance, we counted the percentages of \textit{"Agree"} and \textit{"Disagree"} verdicts relative to the label generated on the basis of MAG keywords. If strictly more experts agreed than disagreed with MAG, the original label was preserved (for 351 of the 577 valid instances).  Otherwise, we took the label suggested by the majority of annotating experts (for 226 of the 577 valid instances).
(3) In the case of ties, we randomly pick a label from the suggested annotations (for 22 of the 226 category-renewed instances). Following this protocol, we obtained 561 instances with expert-curated labels to fine-tune the base models.
\begin{table*}[!t]
\centering
\resizebox{0.8\linewidth}{!}{
\begin{tabular}{cccccc}
\toprule
\textbf{Model} & \textbf{Model\_{FT}} & $\Delta$ (\textbf{Model\_{FT}} - \textbf{Model}) & \textbf{Model} = \textbf{Model\_{FT}} & \multicolumn{1}{l}{\textbf{Total}} & \textbf{Model Type}  \\
(1) & (2) & (3) & (4) & (5) & (6) \\
\midrule
54    & 58        & 4                         & 108               & \multirow{4}{*}{190}      & CNN         \\
59    & 69        & 10                        & 148               &                           & \textbf{RNN}         \\
39    & 39        & 0                         & 190               &                           & DNN         \\
31    & 33        & 2                         &  37               &                           & Transformer \\
\bottomrule
\end{tabular}
}
\caption{Results of the User Study. \textit{FT}: Fine-tuned.}
\label{tab:m-vs-mft}
\end{table*}

\subsection{Fine-tuning Pre-trained Base Models}
We used the 561 labels generated by the experts as a fine-tuning set on the base models reported in \citet{rao2023hierarchical} on the discipline of \textit{Economics} (model-1). We compared the inference performances of the base model (\textbf{Model} in Table~\ref{tab:m-vs-mft}) with those of the fine-tuned model (\textbf{Model\_FT} in Table~\ref{tab:m-vs-mft}) on various neural network architectures, Deep Neural Network (DNN), Recurrent Neural Network (RNN), Convolutional Neural Network (CNN), and Transformers. To benchmark the differences in performances between \textbf{Model} and \textbf{Model\_FT}, we created a small test set from the \href{https://www.ssrn.com/index.cfm/en/}{Social Science Research Network (SSRN)}, which is a website that provides a platform for researchers to share and distribute their research papers and other scholarly work in the social sciences and other related fields. We decided to use the \textit{Economics} SSRN publications because they come with human-currated JEL categories, keywords, and abstracts.

Concretely, we built \href{https://gitlab.ethz.ch/raox/science-clf/-/tree/labelstudio/labelstudio/ssrn_crawler}{a crawler} to download the publication space in \href{https://papers.ssrn.com/sol3/displayjel.cfm}{Economics publications in SSRN}, where all contained research articles in Economics are multi-category-indexed. This means, each publication there is indexed by at least one JEL code and it allows multiple JEL codes per publication. We could easily validate with our multi-label engine in principal, but we focus on \textit{single-label} classifications for this user study.

To create this test set, we randomly sampled 10 instances from each of 19 JEL field classes, which resulted in a sample of 190 test instances. In the implementation of hierarchical classifications reported in \citet{rao2023hierarchical}, we have used \href{https://mlflow.org/}{MLflow} to track and manage ML experiments, with which we have saved all pre-trained base models. Now, based on them, we could seamlessly integrate model fine-tuning and inference with various models.
The inference engine API has been implemented using \href{https://fastapi.tiangolo.com/}{FastAPI} with help from \href{https://pydantic-docs.helpmanual.io/}{Pydantic}. We illustrate the batch inference API in Figure~\ref{fig:inference}, with which users can feed the test set into various models (base or fine-tuned) and obtain predictions. In Appendix~\ref{app:inferece} we provide more details about the inference engine.

\subsection{Benefits of Expert Annotations}
We present the results of user studies in Table~\ref{tab:m-vs-mft}. Specifically, we inspect two types of statistics, the correct predictions of the base and fine-tuned models in Columns (1)-(2), and the identical predictions of the base and fine-tuned models in Column (4). Since each publication is multi-JEL-category-indexed, we count the prediction as "correct" if the indices include the predicted category. Column (1) is the base model trained with the model type specified in Column (6). Column (2) presents the results of the fine-tuned (supervised) model. Column (4) shows that out of a total of 190 test instances, identical predictions were generated by the base and fine-tuned models. We see that fine-tuning with user-generated results has brought benefits to all models except DNN because DNN predicts for all test examples only one class. RNN is the best performer when considering the benefits resulting from expert supervision, because the $\Delta$ in correct predictions has increased the most according to Column (3). Interestingly, fine-tuning a pre-trained Transformer model may not always result in a significant improvement in performance, as we see from a comparison with other base models. However, the current fine-tuning set is too small to draw firm conclusions in this regard.

\section{Conclusions}
In this system demonstration, we utilize a set of standard open-source software (mainly Label Studio \cite{LabelStudio}, \href{https://mlflow.org/}{MLflow} and \href{https://fastapi.tiangolo.com/}{FastAPI}) to configure an annotation and inference engine for scientific publications (SAINE). This system is built on top of hitherto largest multi-class hierarchical classification study across all disciplines in both single-label and multi-label settings (cf.~\citet{rao2023hierarchical}). We illustrate the functionality of the system with a user study in \textit{Economics} and show that the expert inputs into our system can help better understanding the classification process, which benefits the development of a stronger model in the next iteration. We plan to open-source the data and codebase and invite collaborative work in the direction of meta-science.

\section*{Limitations}

Label Studio has some limitations in incorporating existing ML pipelines into the annotation engine, especially, when using customary code. We will discuss this with the developers at Label Studio and see how we can bring the annotation engine and the ML pipeline closer to each other.

In terms of annotator selection, at the moment we have to select the experts for each discipline. However, we would wish to rank the annotators by their field expertise and weigh their annotations, accordingly. One future idea is to automatically compute an associative score between a third-party academic product such as Google Scholar and the publication space. The project \href{https://github.com/poloclub/people-map}{PeopleMap} provides interesting techniques to generate researcher profiles based on their research interests and publications taking as input the Google Scholar profile URLs of researchers. At this stage, Label Studio developers suggest that we add a self-declarative questionnaire to each annotator, which can be used as meta-data on annotators when quantifying the annotation confidence score. Due to time constraints, we have not yet added this questionnaire, as the experts in the current user study are selected by our project PI and have strong expertise in \textit{Economics}.

In terms of annotation efforts, given enough time and annotators, we can better benchmark annotation quality in terms of JEL categories and annotator expertise. Annotator 1 has compared the user experience and performance of \href{https://openai.com/blog/chatgpt/}{ChatGPT} in conducting annotations on a random sample of 100 instances. It turns out that ChatGPT captures high-quality keywords and achieves more than 90\% of the IAA score with that annotator. Considering our annotators' feedback that it is time-consuming to extract keywords for humans, it makes sense to use ChatGPT as an annotation-assisting engine for keyword extraction and label assignment.

\section*{Ethics Statement}
We acknowledge that our system may involve processing potentially sensitive data, and we take data privacy and ethical considerations very seriously. In accordance with ethical guidelines of "ACM Code of Ethics", we have taken steps to protect the privacy of experts who may be affected by our research. We have also made efforts to ensure that our system and its annotations are unbiased and fair. We believe that our work will help foster greater transparency and understanding in scientific research, and we welcome collaboration and feedback from the scientific community to further advance the ethical and responsible use of AI in research.

\section*{Broader Impact Statement}
Our annotation engine and inference engine can further support downstream meta-science projects. We list a few interesting questions we can answer using our pipeline.
\begin{enumerate}
    \item \textbf{[For students.] }Which fields of research are more impactful/growing?
    \item \textbf{[For policy makers.]} How to design education for cross-/inter-/pluridisciplinary studies?
    \item \textbf{[For department and tenure committees.]} How to benchmark output and impact levels of an untenured scholar?
    \item \textbf{[For funding institutions.]} How to measure/quantify inter-/pluri-disciplinary standards for institutions such as \href{https://snis.ch/}{SNIS} and \href{https://www.snf.ch/en}{SNSF} which emphasize the interdisciplinarity of research?
    \item \textbf{[For librarians.]} How can one effectively organize bibliographical resources across disciplines and departments in one university?
\end{enumerate}

We plan to add other disciplines covered by \citet{rao2023hierarchical} to our annotation engine. We would also like to incorporate subjective (self-declaration) and objective measurements (e.g., Google Scholar profile integration) into the annotation pipeline. This may help develop confidence scores of one annotation/annotator.

\section*{Acknowledgements}
We thank the colleagues at \href{https://ds3lab.inf.ethz.ch/}{DS3Lab} for providing valuable feedback when prototyping the system design. Without the strong support of our expert group at the \href{https://cae.ethz.ch/}{Chair of Applied Economics}, the user study would not have been possible. We thank Ms.~Piriyakorn Piriyatamwong for her technical support to our project. We appreciate that Label Studio has offered us an academic license for the project, which allows us to invite more experts to contribute in the long run. The user agreement and terms of an academic license are listed \href{https://labelstud.io/academic/}{here}.


\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}.